Resolved Measurements of the CO-to-H2 Conversion Factor in 37 Nearby Galaxies

We measure the CO-to-H2 conversion factor (α CO) in 37 galaxies at 2 kpc resolution, using the dust surface density inferred from far-infrared emission as a tracer of the gas surface density and assuming a constant dust-to-metal ratio. In total, we have ∼790 and ∼610 independent measurements of α CO for CO (2–1) and (1–0), respectively. The mean values for α CO (2–1) and α CO (1–0) are 9.3−5.4+4.6 and 4.2−2.0+1.9M⊙pc−2(Kkms−1)−1 , respectively. The CO-intensity-weighted mean is 5.69 for α CO (2–1) and 3.33 for α CO (1–0). We examine how α CO scales with several physical quantities, e.g., the star formation rate (SFR), stellar mass, and dust-mass-weighted average interstellar radiation field strength ( U¯ ). Among them, U¯ , ΣSFR, and the integrated CO intensity (W CO) have the strongest anticorrelation with spatially resolved α CO. We provide linear regression results to α CO for all quantities tested. At galaxy-integrated scales, we observe significant correlations between α CO and W CO, metallicity, U¯ , and ΣSFR. We also find that α CO in each galaxy decreases with the stellar mass surface density (Σ⋆) in high-surface-density regions (Σ⋆ ≥ 100 M ⊙ pc−2), following the power-law relations αCO(2–1)∝Σ⋆−0.5 and αCO(1–0)∝Σ⋆−0.2 . The power-law index is insensitive to the assumed dust-to-metal ratio. We interpret the decrease in α CO with increasing Σ⋆ as a result of higher velocity dispersion compared to isolated, self-gravitating clouds due to the additional gravitational force from stellar sources, which leads to the reduction in α CO. The decrease in α CO at high Σ⋆ is important for accurately assessing molecular gas content and star formation efficiency in the centers of galaxies, which bridge “Milky Way–like” to “starburst-like” conversion factors.

not directly observable in many cases in the the cold molecular ISM due to its high transition energies (hν/k B ∼ 510 K) for the lowest rotational levels.As a result, low-J CO emission lines are the most widely used tracer for the molecular ISM, as the second most abundant molecule with strong millimeter rotational lines that can be excited at typical temperatures in molecular clouds.The standard practice is to use a CO-to-H 2 conversion factor α CO , as follows where Σ mol [M ⊙ pc −2 ] is the surface density of molecular ISM (including mass of Helium), α CO [M ⊙ pc −2 (K km s −1 ) −1 ] is the "CO-to-H 2 conversion factor", and W CO [K km s −1 ] is the integrated intensity of the CO emission at rest frame frequency.The conventional α CO in the Milky Way (MW) is α MW CO = 4.35 M ⊙ pc −2 (K km s −1 ) −1 for CO (1−0) 1 .In mass surface density units with a factor of 1.36 for Helium included, this is equivalent to α CO (1−0) = 4.35 M ⊙ pc −2 (K km s −1 ) −1 .(Solomon et al. 1987;Strong & Mattox 1996;Abdo et al. 2010).The "(1−0)" and "(2−1)" symbols represent the CO rotational transition which α CO is derived for and W CO is measured, i.e.CO (1−0) stands for the CO J = 1 → 0 rotational transition at ∼ 115 GHz (λ ∼ 2.6 mm) and CO (2−1) stands for the CO J = 2 → 1 rotational transition at ∼ 230 GHz (λ ∼ 1.3 mm), respectively.In this work, we focus on the 12 C 16 O isotopologue and use CO for 12 C 16 O for simplicity.
CO (1−0) had been the most frequently measured CO transition, and α CO (1−0) is thus the fiducial case for α CO in the literature.Meanwhile, CO emission from the next highest rotational level, i.e.CO (2−1), has become more and more common with modern instruments, e.g.ALMA, throughout the last decade.Directly deriving α CO (2−1) has attained its own importance.Thus we will present both α CO (1−0) and α CO (2−1) in this work.
With the precise measurements of CO emission from modern instruments, our understanding of α CO has become the factor that limits our ability to precisely study the molecular ISM and star formation in nearby galaxies.Observations have shown at least two main trends in the variation of α CO .The first one is that α CO tends to increase at lower metallicity or lower dust-to-gas ratios (e.g. up to ∼ 1 dex higher than the MW value at ∼ 0.2 solar metallicity, Israel 1997b; Leroy et al. 2011).This enhanced α CO is often explained by the decrease in CO emission relative to the cloud mass defined by H 2 as shielding for CO weakened at lower metallicity (Papadopoulos et al. 2002;Grenier et al. 2005;Wolfire et al. 2010;Genzel et al. 2012;Planck Collaboration et al. 2011;Accurso et al. 2017;Gong et al. 2020;Madden et al. 2020;Hirashita 2023).This phenomena is often known as the "CO-dark gas." The second trend is that α CO appears to be lower in the central ∼ kpc of some galaxies (it can be a factor 5-10 times lower, Israel 2009aIsrael ,b, 2020;;Sandstrom et al. 2013;Teng et al. 2022Teng et al. , 2023b)).It is also observed to be lower in (ultra-)luminous infrared galaxies (U/LIRGs, Downes et al. 1993;Downes & Solomon 1998;Papadopoulos et al. 2012;Herrero-Illana et al. 2019).This trend towards lower α CO in galaxy centers and starbursts likely results from a combination of higher gas temperature, larger line width, and lower CO optical depth, which breaks the relationship between molecular cloud mass and line width that one would expect from isolated, virialized clouds (Shetty et al. 2011;Bolatto et al. 2013).This phenomena is often referred to as the "starburst conversion factor."Because galaxy centers and U/LIRGs tend to be bright in CO and thus easily observed, understanding this starburst conversion factor is important to making best use of a wide range of extragalactic observations in characterizing the star formation efficiency, gas dynamics and H I-to-H 2 transition conditions.Bolatto et al. (2013) proposed a formula for conversion factor treating the CO-dark gas and the starburst trend independently and simultaneously.This formula aimed to predict both the spatially resolved measurements from Sandstrom et al. (2013, including galaxy centers) and the (U)LIRGs measurements in Downes & Solomon (1998).The formula reads: Here Z ′ is the metallicity relative to the solar value, which traces the CO-dark gas effect; Σ 100 Total is the total surface density (Σ Total = Σ gas + Σ ⋆ ) in unit of 100 M ⊙ pc −2 , which is the proposed observational tracer and threshold for regions where the decrease of α CO occurs, i.e. galaxy centers and (U)LIRGs.The authors found that with a threshold at Σ 100 Total ≥ 1 (a threshold related to self-gravitating giant molecular clouds), the relation α CO ∝ Σ 100 Total −0.5 reproduces the trend found in galaxy centers and ULIRG samples.A similar formula was also suggested by Ostriker & Shetty (2011), where the authors suggested a power-law relationship between α CO and Σ gas to describe the decrease in α CO needed for their simulations to match observations.They showed that a relation of α CO ∝ Σ −0.5 gas over the surface density range 10 2 to 10 3 M ⊙ pc −2 brings the observations and simulations into agreement.
To better characterize how α CO depend on local environments, spatial resolved measurements of α CO are required.To measure α CO , one needs to measure Σ mol independent of (single-line) CO intensity, and then divide it by the measured CO intensity.This could be achieved by several methodologies, e.g. using virial mass estimates (see the review in Mc-Kee & Ostriker 2007), modelling multiple spectral lines (e.g.Cormier et al. 2018;Teng et al. 2022Teng et al. , 2023b)), converting γray emission (e.g.Abdo et al. 2010;Ackermann et al. 2012), and tracing gas mass with dust mass (Israel 1997a,b;Leroy et al. 2007Leroy et al. , 2009bLeroy et al. , 2011;;Bolatto et al. 2011;Planck Collaboration et al. 2011;Sandstrom et al. 2013;Schruba et al. 2017;den Brok et al. 2023).Based on existing resources, most of the methods are not practical for a survey in tens of nearby galaxies due to the requirement in target brightness or total observing time.The most feasible methodology is to use dust as a tracer for gas mass, where dust mass can be derived from infrared (IR) data observed with Herschel (Pilbratt et al. 2010), WISE (Wright et al. 2010), and Spitzer.
In this work, we measure α CO with dust as the tracer for total gas mass.We use dust masses derived from modeling the far-IR spectral energy distribution (SED) to trace the total gas mass.The key assumption we make is a constant fraction of heavy elements locked in the solid phase, i.e. a dustto-metals ratio (D/M), which allows us to convert measurements of dust surface density, H I surface density and metallicity into molecular gas mass.The assumption of approximately constant D/M is supported by dust evolution models (Dwek 1998;Hirashita & Kuo 2011;Asano et al. 2013;Feldmann 2015) and kpc-scale measurements (Issa et al. 1990;Leroy et al. 2011;Draine et al. 2014;Vílchez et al. 2019;Chiang et al. 2021) in high-metallicity (12 + log(O/H) ≳ 8.2) galaxy disks, matching the region of interest in this work.In simulations (e.g.Aoyama et al. 2020;Choban et al. 2022), an approximately constant D/M results from efficient dust growth in the ISM, i.e. the majority of the refractory elements are locked in solid grains quickly.Although there are also studies that found variations in D/M with both depletion measurements (Jenkins 2009;Jenkins & Wallerstein 2017;Roman-Duval et al. 2019) and emission measurements (Roman-Duval et al. 2014, 2017;Chiang et al. 2018;De Vis et al. 2019), no widely agreed-upon prescription for the environmental dependence of D/M has been found thus far.The other reason we assume a constant D/M is that we anticipate the variation of D/M (≤ 2 times) to be smaller in comparison to that of α CO (up to ∼ 10 times) in normal galaxy disks.
Given the challenges in measuring spatial resolved α CO , there have been few studies with large samples of resolved measurements in galaxies.Sandstrom et al. (2013) looked at the overlap of IR from KINGFISH (Kennicutt et al. 2011), CO from HERACLES (Leroy et al. 2009a), and H I from THINGS (Walter et al. 2008) and made measurements of α CO in 26 galaxies with ∼40 ′′ resolution elements, which is the resolved study with one of the largest sample size.Most other studies in the field, e.g.Hunt et al. (2015), Accurso et al. (2017) and COMING (Yasuda et al. 2023), focused on galaxy-scale α CO .Moreover, a survey of α CO at fixed physical scale, which allows us to evaluate the environmental dependence of α CO fairly, is also missing as previous spatially resolved studies (e.g.Leroy et al. 2011;Schruba et al. 2012;Sandstrom et al. 2013) tend to perform their analysis at fixed angular resolution.In this work, we will measure α CO across 37 nearby galaxies at a fixed ∼2 kpc resolution.This study is made possible by several surveys of resolved CO intensities in the past two decades: the Nobeyama CO Atlas of nearby galaxies (CO Atlas, Kuno et al. 2007), the HERA CO Line Extragalactic Survey (HERACLES, Leroy et al. 2009a), the CO Multi-line Imaging of Nearby Galaxies project (COMING, Sorai et al. 2019), the Physics at High Angular resolution in Nearby Galaxies project (PHANGS-ALMA, Leroy et al. 2021), and recent IRAM 30m observations (PI: A. Schruba; see Leroy et al. 2022).
This paper is presented as follows.In Section 2, we explain our methodology for deriving α CO from the data.We describe the data sets necessary for this work and how we constrain other physical quantities from observations in Section 3. We present the α CO measurements and their correlations with local and galaxy-integrated conditions in Section 4. In Section 5, we investigate an observed power-law relation between α CO and Σ ⋆ in the high-surface-density regions, and provide a prescription for α CO based on our findings.We discuss the physical interpretations of our results and how it compares to the literature findings in Section 6.Finally, we summarize our findings in Section 7.

CALCULATING α CO
To calculate α CO , we first estimate Σ mol without using CO emission and an adopted conversion factor.In this study, we use a dust-based strategy to estimate Σ mol by assuming a value for the fraction of metals locked in solid phase, i.e. the dust-to-metals ratio (D/M).The D/M is defined as: where Σ gas = Σ atom + Σ mol is the total neutral gas surface density, Σ mol is the molecular gas surface density, and the metallicity Z converts Σ gas to a "metal mass surface density".By replacing Σ mol in the above equation with the definition of α CO in Eq. 1, we have: where Σ dust , Z, Σ atom , and W CO are measurable quantities, thus, by assigning a value of D/M, we can then calculate α CO with our data set.The uncertainty of α CO is propagated from the uncertainties of Σ dust , metallicity, Σ atom , and W CO .The typical uncertainty of the pixel-by-pixel α CO measurements in this work is in the range of 0.2-0.5 dex2 .Dust-based α CO measurements in the literature usually have formulae similar with Eq. 3 but with different assumptions.For example, Israel (1997b) and Leroy et al. (2009b) assumed a fixed dust-to-gas ratio (D/G) in their sample galaxies to derive α CO .On the other hand, Sandstrom et al. (2013, also see Leroy et al. (2011), den Brok et al. (2023) and Yasuda et al. (2023)) assumed that the D/G remains approximately constant in a certain spatial region, e.g.kpc scale or entire galaxy.With this assumption, the authors are able to derive α CO by minimizing the scatter in D/G in a group of nearby pixels, treating α CO and D/G as free parameters.This method has the advantage in not forcing the value of D/G and the disadvantage in sacrificing the spatial resolution.
Generally speaking, we could assume different D/M values in each pixel when we apply Eq. 4.However, despite the previous efforts on studying the evolution of D/M (De Vis et al. 2019;Aoyama et al. 2020;Péroux & Howk 2020;Chiang et al. 2021;Choban et al. 2022;Roman-Duval et al. 2022), we do not have a well established prescription of how D/M depends on local environments, or the prescription is not more accurate than simply assuming a constant D/M.Studies have shown that at 12 + log(O/H) > 8.2, the D/M falls in the range3 between 0.4-0.7,e.g.0.5 (Jenkins 2009, taking F ⋆ = 1), 0.72 (Leroy et al. 2011), 0.46 (Rémy-Ruyer et al. 2014), 0.68 (Draine et al. 2014), 0.7 (Feldmann 2015), 0.56 (Chiang et al. 2018), and 0.40-0.58(Chiang et al. 2021).Meanwhile, studies have shown that in galaxy centers, α CO can vary by a factor of 10 (Sandstrom et al. 2013;Bolatto et al. 2013;Israel 2020;Teng et al. 2022Teng et al. , 2023b)), which is a significantly larger dynamic range than D/M.In the following, we will take a constant D/M = 0.55 (mean of 0.4-0.7)as the fiducial case for deriving α CO .A possible drawback by assuming a constant D/M is that D/M has been shown to decrease toward lower metallicity (e.g.Rémy-Ruyer et al. 2014;Hirashita & Kuo 2011;Chiang et al. 2018;De Vis et al. 2019).Thus, the assumed D/M value that is appropriate for galaxy centers is likely too high for low-metallicity regions (usually the outer disks), resulting in an underestimation of α CO (Eq.4) in the outer disk.We will discuss the case of a varying D/M with a toy model in Appendix A. We will also discuss the possible uncertainties in Σ dust derivation in Section 3.

CO (1−0) and CO (2−1) cases
In this multi-wavelength, multi-galaxy study, we do not always have both the lowest-J CO emission lines for all the target galaxies.Studies have been using CO line ratios to convert the intensity between CO emission lines.For an indepth discussion on low-J CO line ratios, we refer the readers to Leroy et al. (2022).In the literature, perhaps the most frequently used method to treat different CO line coverage is converting everything to W CO (1−0) with a constant CO line ratio (e.g.Sandstrom et al. 2013;Sun et al. 2020;Chiang et al. 2021).This method allows us to uniformly use α CO (1−0) for calculating Σ mol in the study.Theoretically, the line ratio can vary with excitation conditions like gas temperature and linewidth.Thus, we expect R 21 to trace the local environmental conditions.Taking the line ratio for W CO (1−0) and W CO (2−1) as an example: where R 21 is usually falls in the range of 0.3 to 0.9 with mean value ∼ 0.65 for normal star-forming galaxies, and it is expected to be higher in galaxy centers (Leroy et al. 2009a(Leroy et al. , 2022(Leroy et al. , 2023;;den Brok et al. 2021;Yajima et al. 2021).
In this work, however, we will not adopt the simple strategy as our fiducial case because most of our target galaxies have CO (2−1) data, that α CO (2−1) has attained its own importance due to modern observations, and that the variation in R 21 is non-negligible.We will present four solutions of α CO in parallel, two without any conversions between CO (2−1) and CO (1−0), and two with different prescriptions of the line ratio: for galaxies without W CO (1−0) data, converted with a constant R 21 .
For the third method, we adopt the constant R 21 = 0.65 from Leroy et al. (2022).For the last method, we adopt the MIRdependent formula suggested by Leroy et al. (2023), namely We follow the suggestion in Leroy et al. (2023) Leroy et al. 2022Leroy et al. , 2023)).

DATA
We measure α CO in 37 nearby galaxies in this study.To measure α CO with our dust-based methodology (Section 2), the data sets required are dust surface density (Σ dust , from IR SED modelling), CO low-J rotational line integrated intensity (W CO ), atomic gas surface density (Σ atom , from H I 21 cm line emission), and metallicity (Z, from gas-phase oxygen abundance in H II regions).We first select our sample galaxies from the dust catalog of z=0 Multiwavelength Galaxy Synthesis (z0MGS) (Leroy et al. 2019;Chastenet et al. 2021, J. Chastenet et al. in preparation).From this large sample, we pick the 49 galaxies with both low-J CO rotational line and H I data available, including our own new H I data sets.We design our study with a common resolution of 2 kpc, which draws a limit in sample selection at distance ∼ 20 Mpc since the worst resolution data in our sample usually have angular resolution around 20 ′′ .High-inclination (> 80 • ) targets are also excluded.We further exclude 4 galaxies 4 that satisfy all above conditions but do not satisfy the signal-to-noise ratio conditions that will be described in Section 3.3.The selection yields 37 galaxies in our sample.We list the properties of these galaxies in Table 1.(2021) and J. Chastenet et al. (in preparation).We briefly summarize the process below.We obtain the dust emission SED in the IR observed by two space telescopes: the 3.4, 4.6, 12 and 22 µm bands with the Wide-field Infrared Survey Explorer (WISE, Wright et al. 2010), and the 70, 100, 160 and 250 µm bands with the Herschel Space Observatory (Pilbratt et al. 2010).The Herschel and WISE maps are first convolved to SPIRE 250 PSF and then to a 21 ′′ Gaussian PSF using the SPIRE 250-to-Gauss-21 ′′ kernel from Aniano et al. (2011).The 21 ′′ PSF is the "moderate" Gaussian from Aniano et al. (2011) that provides relatively high angular resolution without amplifying image artifacts.Finally, these maps are convolved to the desired resolution: a Gaussian PSF with spatial resolution corresponding to FWHM of 2 kpc.
After convolving the IR maps, we fit the dust SED with the Draine & Li (2007) physical dust model with the dust opacity calibration derived in Chastenet et al. (2021).This calibration is based on metallicity measured with "direct" electrontemperature-based methods, which is consistent with the strong lien calibration adopted in this work (S-calibration in Pilyugin & Grebel 2016), and yields reasonable D/M.Thus, the calibration ties dust opacity, D/M and metallicity into one framework.The complete set of data products from the fitting includes the maps of dust mass surface density (Σ dust ), interstellar radiation field (the minimum radiation field U min and the fraction of dust heated by the power-law radiation field γ) and the fractional dust mass in the form of polycyclic aromatic hydrocarbons (q PAH ).The maximum radiation field is fixed at U max = 10 7 and the power-law index for radiation field distribution is fixed at α = 2. From the fitted U min and γ, we can derive the dust-mass-averaged radiation field U, which is the fiducial tracer for radiation field in this work.
We also note that we assume fixed dust properties in our dust SED fitting throughout this study, which is the most frequently adopted strategy in the literature.The accuracy of Σ dust estimates, derived by fitting the IR SED with dust emission models, can be affected by variations in the dust opacity.Interstellar dust grains are not uniform in their chemical composition, size distribution, and shape, leading to variations in their opacity (e.g.Draine & Li 2007;Hirashita & Voshchin-nikov 2014;Draine & Hensley 2021).In the MW, Stepnik et al. (2003) found that the dust opacity increases by ∼ 3 times from the diffuse ISM to the dense clouds.The authors argued that the increase in dust opacity is resulted from the deficit of small grains due to grain-grain coagulation.It is challenging to measure the variation in opacity of interstellar dust since it is degenerate with the environmental dependence of α CO and D/M.Moreover, many of the mechanisms that affects dust opacity, e.g.grain-grain coagulation and ice mantles, are smoothed out in kpc-scale extra-galactic studies (Galliano et al. 2018), meaning that its variation is likely less observable than the other degenerate factors like α CO and D/M.We note that there are extra-galactic studies that attribute all the variations in dust and gas properties to dust opacity to evaluate its variation, e.g.Clark et al. (2019) found that dust opacity change by a factor ∼ 2 within M74 and ∼ 5 within M83.
Atomic gas surface density.We trace the atomic gas surface density (Σ atom ) with the H I 21 cm integrated intensity (I H I ), assuming the opacity is negligible (e.g., Walter et al. 2008): where the 1.36 factor accounts for the mass of helium.
We obtain W H I from both literature and new data, as listed in Table 1.The two new H I surveys are the EveryTHINGS survey (P.I.K. M. Sandstrom; I. Chiang et al. in preparation) and the PHANGS-VLA survey (P.I.D. Utomo).The EveryTHINGS survey targets nearby galaxies with Herschel photometric data but without high-resolution H I observations, while the PHANGS-VLA survey focuses on galaxies in the Physics at High Angular resolution in Nearby Galaxies (PHANGS) project5 .Both surveys have their data observed with the C-and D-configurations of Karl G. Jansky Very Large Array (VLA)6 , which yield an angular resolutions in the range of 20 ′′ to 30 ′′ .Both surveys provide new highsensitivity 21 cm H I observations in tens of nearby galaxies.We did not include WHISP (Swaters et al. 2002) data because the galaxies that only have WHISP data have resolution coarser than 2 kpc after convolving to a circular PSF.
CO low-J rotational lines.The integrated intensity of CO low-J rotational lines traces the molecular gas surface density (Eq.1), and is key to this study.We use the compilation of CO mapping assembled by Leroy et al. (2022Leroy et al. ( , 2023) ) from several publicly available CO (1−0) and CO (2−1) data: • CO (1−0) data from the COMING survey (Sorai et al. 2019) and the CO Atlas (Kuno et al. 2007).
The source of CO data for each galaxy is listed in Table 1, where CO (1−0) and CO (2−1) are listed separately.All these literature measurements focus on the 12 C 16 O isotopologue, hereafter CO for simplicity.
Surface densities of stellar mass and SFR.We trace the surface densities of stellar mass and SFR (Σ ⋆ and Σ SFR , respectively) using the data and conversion formulae presented in the z0MGS survey (Leroy et al. 2019).We utilize the z0MGS compilation of the background-subtracted WISE (Wright et al. 2010) λ ∼ 3.4 and 22 µm (hereafter WISE1 and WISE4, respectively) data and the Galaxy Evolution Explorer (GALEX, Martin et al. 2005) λ ∼ 154 nm (hereafter FUV) data.
We use WISE1 data to trace stellar mass surface density (Σ ⋆ ) with: ⋆ is calculated from the galaxy-by-galaxy SFRto-WISE1 ratio, a "specific SFR-like" quantity, with the prescription calibrated in Appendix A of Leroy et al. (2019).
We use FUV and WISE4 data to trace the SFR surface density (Σ SFR ) also with the conversion formula suggested by z0MGS (Leroy et al. 2019;Belfiore et al. 2023).For galaxies with both FUV and WISE4 available, we use: For NGC3953 and NGC4689, where only WISE4 is available, we use: For both WISE and GALEX maps, we blank the regions with foreground stars identified in the z0MGS database.We interpolate the intensities in the blanked region with a Gaussian kernel FWHM = 22.5 ′′ (the adopted WISE and GALEX maps have FWHM = 15 ′′ ) with the function interpolate_replace_nans in astropy.convolution.This interpolation is done on the maps before any convolution, reprojection, or unit conversion.Regarding the WISE maps, this treatment is only done for the maps used for calculating Σ ⋆ and Σ SFR .For the WISE maps used in dust SED fitting, we refer the readers to J. Chastenet et al. (in preparation).
Specific SFR.With the measurements of Σ SFR and Σ ⋆ , we calculate the specific SFR (sSFR) as: where i is the summation over pixels in a galaxy.Meanwhile, we calculate the spatially resolved sSFR (rsSFR) as: Metallicity.We use the abundance of oxygen, 12 + log(O/H), to trace the metallicity (Z).We assume a fixed abundance pattern, i.e. a constant oxygen-to-total-metal mass ratio.The conversion from 12 + log(O/H) to Z then becomes: where 0.0134 and 8.69 are the adopted Z ⊙ and 12 + log(O/H) ⊙ , respectively (Asplund et al. 2009).
We calculate 12 + log(O/H) for each pixel as a function of the galactocentric distance by adopting the radial gradient of 12 + log(O/H) derived from measurements in H II regions.We use data from two surveys: the PHANGS-MUSE survey (Emsellem et al. 2022;Groves et al. 2023) and the Zurita et al. ( 2021) compilation.We use the Pilyugin & Grebel (2016) S-calibration7 (hereafter PG16S) as the fiducial calibration for 12 + log(O/H).PG16S is a calibration method that shows good agreement with direct metallicity measurements (Croxall et al. 2016;Kreckel et al. 2019).Since PG16S only relies on lines covered by MUSE, the 12 + log(O/H) measurement can be expanded to the full PHANGS-MUSE dataset in our future works.
For galaxies in the PHANGS-MUSE survey, we adopt radial 12 + log(O/H) gradients presented in Santoro et al. (2022), which are calculated with the PG16S calibration.For galaxies that only appear in the Zurita et al. (2021) emission data compilation, we use the Zurita et al. (2021) data to calculate the PG16S 12 + log(O/H) and then fit the radial 12 + log(O/H) gradient in these galaxies.We only consider galaxies that have at least 5 measurements spanning at least 0.5R 25 in the Zurita et al. (2021) data table.The uncertainties of the 12 + log(O/H) gradient is not explicitly provided in either work.We will assume a 0.1 dex scatter for 12+log(O/H) derived from gradients as suggested in Berg et al. (2020).
For galaxies without measurements of 12 + log(O/H) in either Zurita et al. (2021) or Santoro et al. (2022), we use the two-step strategy proposed in Sun et al. (2020) to estimate their 12 + log(O/H).First, we use a mass-metallicity relation to predict 12 + log(O/H) at one effective radius (R e ) in a given galaxy.Second, we extend the prediction with a radial gradient of −0.1 dex/R e suggested by Sánchez et al. (2014).We characterize the mass-metallicity with a function of the form: 12 where between the best fit and data is 0.13 dex, meaning that there is still intrinsic scatter in 12 + log(O/H) that is not explained by the mass-metallicity relation and the adopted radial gradients, e.g. the azimuthal variations (Williams et al. 2022).
We use this fitted relation to derive the radial metallicity gradient of galaxies without metallicity measurements with the M ⋆ and R e listed in Table 1.We will assume an uncertainty of 0.15 dex (rounding up 0.013 + 0.13 dex) for galaxies with this type of 12 + log(O/H) measurements.
Studies have reported that the PG16S calibration could result in 12 + log(O/H) value lower than other calibrations (e.g.De Vis et al. 2019).Aligning with that, there are also studies reporting that the 12 + log(O/H) calibrated with PG16S in Orion Nebula and other H II regions in the solar neighborhood have values ∼ 0.2 dex lower than the solar reference value (e.g.Esteban et al. 2022).This effect could lead to an underestimate of Z' and thus an overestimate in α CO with our methodology.For consistency with the direct metallicity calibration used in Chastenet et al. (2021) for calibrating dust opacity, we will use PG16S in this work.

Uniform Data Processing
The analyses in this work are done at a common resolution of 2 kpc for all data.For H I, CO, Σ ⋆ , and Σ SFR , we convolve them to a circular Gaussian PSF with a FWHM corresponding to 2 kpc, using the astropy.convolutionpackage (Astropy Collaboration et al. 2013, 2018, 2022).The images are then reprojected onto a grid with a pixel size of one third of the FWHM (that is, we oversample at roughly the Nyquist sampling rate) with the astropy affiliated package reproject.The convolution and reprojection of the dust maps are done in J. Chastenet et al. (in preparation).They convolve the IR maps into the final resolution using kernels from Aniano et al. (2011).Note that the convolution is done before SED fitting for dust properties.The galactocentric radius and metallicity are directly calculated on the final pixel grids.All the surface density and surface brightness quantities presented in this work have been corrected for inclination, listed in Table 1.

S/N Mask and Completeness
S/N Mask.For statistical quantities that only involve α CO , e.g.mean values and percentiles, we masked out the low signal-to-noise ratio (S/N) pixels.Specifically, we blank pixels with S/N < 1 in W CO and Σ mol .Note that Σ mol here is derived from Σ dust (with IR photometry), metallicity and Σ atom (Eq.4), thus the uncertainty of Σ mol is propagated from the uncertainties of these three quantities and the IR photometry.
Completeness.For statistical quantities that involve α CO and another quantity (X), e.g.correlations and linear regression, in addition to the S/N mask, we only calculate with data with high (> 50%) completeness in X as the trend confidence criteria.The completeness for data with X i ≤ X < X f , or [X i , X f ), is defined as: where the definition of "S/N > 1" is the same as the S/N mask earlier this subsection.We show the completeness and the 50% threshold in our data set in Figure 1.For most quantities, the CO (2−1) data has a better completeness coverage than the CO (1−0) data.Note that at the high-U end, the completeness fluctuates around 50%.We treat all data with logU > 0.75 as incomplete for simplicity.

RESULTS
In total, we measure resolved α CO values across 37 galaxies, including ∼ 790 and ∼ 610 independent measurements8 from CO (2−1) and CO (1−0) data, respectively.We summarize the measurements in Table 2 and the distribution of α CO in Figure 2.For each galaxy, we report the mean, COintensity-weighted mean, 16 th -84 th percentiles, and number of pixel-by-pixel measurements of α CO .
Besides the simple mean, we also calculate the W COweighted mean, which better reflects the α CO value to be adopted for data at coarser resolution.The W CO -weighted mean for α CO (2−1) and α CO (1−0) are 5.69 and 3.33, respectively.Both values are lower than the simple mean, which indicates that the α CO values are lower in brighter regions.Unless specified otherwise, we use W CO -weighted mean α CO for galaxy-integrated analysis in the following content.
When we include CO (2−1) data for galaxies without CO (1−0) measurements with a variable R 21 (values with † in Table 2), the mean and W CO -weighted mean of α CO (1−0) increase to 4.72 and 3.48, respectively, indicating that galaxies without CO (1−0) measurements in this data set tend to have larger α CO .This is also visualized in the bottom panel of Figure 2. Also, we find that the distribution of α CO (1−0) does not differ much between the two R 21 prescriptions adopted in this work: the fixed R 21 = 0.65 and the I W4 -dependent R 21 .
In Table 2, we also list the W CO recovery fraction (W CO %), which is the percentage of W CO recovered (above the S/N mask) in the pixel-by-pixel analysis in each galaxy.We notice a few galaxies with low W CO recovery fraction, meaning that there are significant numbers of pixels with W CO detection removed from analysis.The main reason for NGC3631 and most galaxies with recovery fraction above 50% is low sensitivity in dust/IR data.In NGC3631, > 85% of pixels removed have S/N < 1 in the IR bands.In NGC3198 and NGC4625, the sensitivity in dust/IR only explains < 60% of pixels removed.The rest of the pixels were removed due to S/N < 1 in Σ mol , a combined effect of Σ dust , Σ atom and 12 + log(O/H).This type of pixels have S/N > 1 in W CO and S/N < 1 in Σ mol , likely indicating a small α CO .

Correlations with Local Conditions
We measure the pixel-by-pixel correlations between α CO and several parameters describing local physical conditions.These results are summarized in Table 3 and visualized in Figure 3-Figure 4. Note that the correlations and linear regressions9 are calculated with data in the complete zone only (with data completeness > 50%; see Section 3.3 for details).The errorbars in Figure 3-Figure 4 include both the scatter within a bin and the uncertainties of pixel-by-pixel measurements.We first bootstrap the measurements by 1000 times with uncertainties, and then sample the 16th-and 84percentiles in each bin from the bootstrapped sample as the errorbars.We apply the same method for visualizing the other binned data in this work.
logU has the strongest negative correlation with α CO (2−1) , meaning that α CO (2−1) decreases toward regions with stronger interstellar radiation field strength.This is consistent with the picture that α CO decreases with higher gas temperature and larger line width (Bolatto et al. 2013).It is also the case that a higher logU might correspond to a lower Σ dust as a caveat of our fitting methodology (equivalent to "fixing β" in modified blackbody models, see Shetty et al. 2009a,b).However, since we do not see a strong Σ dust -to-α CO (2−1) correlation, the above effect should be minor.Several other quantities also show moderate (negative) correlations with α CO (2−1) , i.e.Σ SFR and rsSFR.Studies have shown strong correlations between logU and Σ SFR (e.g.Hirashita & Chiang 2022;Chiang et al. 2023, J. Chastenet et al. in preparation).Another quantity that shows moderate correlation is W CO .α CO is expected to anti-correlate with W CO due to either external pressure or other dynamical effects (e.g.Bolatto et al. 2013;Sun et al. 2018).The power-law index for W CO is within the previously reported range of −0.32 to −0.54 (Narayanan et al. 2012;Gong et al. 2020;Hunt et al. 2023).Compared to α CO (2−1) , the correlations between α CO (1−0) and local conditions are overall weaker.log rsSFR has the strongest correlation with α CO (1−0) in all three R 21 cases, and logU is the second strongest.
For all combinations of α CO and local conditions, we perform linear regression with the functional form: where x10 is the local condition and both m and b are free parameters, all of which are listed in Table 3.In the same table, we also report the root-mean-square error (∆ rms ) between the measured and fitted log α CO .Most formulas have ∆ rms ∼ 0.2 dex.W CO , U and Σ SFR have the strongest correlations and smallest ∆ rms in general.One quantity that is often used for parametrizing α CO is 12 + log(O/H).In our measurements, α CO has moderate to weak correlation with 12 + log(O/H) in all cases.The slope from linear regression (m) is −2.6 for α CO (2−1) and −1.8 to −2.3 for α CO (1−0) .These values are mildly steeper than most literature values (∼ −1.6 to −2.0, Bolatto et al. 2013;Hunt et al. 2015;Accurso et al. 2017;Sun et al. 2020), but are still within previous reported range, e.g.−2.0 to −2.8 in Schruba et al. (2012).Note that our data set is less suitable for an in-depth study on 12 + log(O/H) since > 80% of our data is concentrated in a small dynamic range of 8.4 ≤ 12 + log(O/H) ≤ 8.6.
We also calculate the correlation between physical quantities and normalized α CO .In this calculation, α CO is normalized by the W CO -weighted mean in each galaxy.For most quantities, the correlation become weaker after normalization.Meanwhile, Σ ⋆ has a stronger correlation with normalized α CO in all cases in Table 3, indicating that Σ ⋆ traces the intra-galaxy α CO variations after normalization of galaxy-togalaxy differences.
Overall, we have shown that α CO (2−1) has stronger correlations with local conditions than α CO (1−0) .Among the local quantities, W CO , U, and Σ SFR usually have stronger correlations with α CO .We do not see a strong correlation coefficient between α CO and 12 + log(O/H), one of the most frequently used quantity to model α CO .

Correlation with Galaxy-Averaged Quantities
Besides kpc-scale variations, we also test how galaxyaveraged α CO (2−1) and α CO (1−0) vary between galaxies, and how their variations correlate with galaxy-averaged properties.The results are visualized in Figure 5-Figure 6.We use the symbol <X> to represent the simple averaged value of quantity X in each galaxy, i.e.Σ gal i X i /Σ gal i , where Σ gal i is the summation over all pixels in a galaxy.There are three quantities for which we do not apply simple averages: (1) as mentioned in Section 4, we will present the W CO -weighted mean of α CO ; (2) logU gal is calculated as log Σ gal i Σ dust,i U i /Σ gal i Σ dust,i to reflect the dustmass-weighted averaged ISRF; (3) sSFR is calculated as The errorbars in Figure 5-Figure 6 show the 16th-84th percentiles of the corresponding quantity.In Figure 6, we also include α CO (1−0) calculated with CO (2−1) data with a I W4 -dependent R 21 .We only include one R 21 prescription here for clarity.
We report the correlation coefficients and the corresponding p-values in each panel of Figure 5-Figure 6.Compared to the results in Section 4.1, we note that whether a quantity has significant correlation with α CO and the strength of the correlation often differ between the spatially resolved case and the galaxy-averaged case.Several quantities show significant correlations with <α CO >. 12 + log(O/H) and <Σ SFR > seem to show stronger correlation with <α CO > for both CO (2−1) and CO (1−0) than in the spatially resolved case.The insignificant correlation between <Σ ⋆ > and α CO (1−0) is consistent with the findings in Carleton et al. (2017) and Dunne et al. (2022) assuming that Σ ⋆ dominates the total mass surface density and that CO (1−0) dominates the CO measurements.

POWER-LAW DEPENDENCE OF THE CONVERSION FACTOR ON STELLAR MASS SURFACE DENSITY
In the α CO prescription proposed in Bolatto et al. (2013), the authors use a power law with Σ Total (=Σ ⋆ +Σ gas ) to trace the changes in α CO due to CO emissivity variations (related to gas temperature and opacity) and a threshold in Σ Total to trace where the effects become important.Inspired by their work and motivated by the necessity of improving α CO prescriptions in galaxy centers, we examine whether a similar functional form applies to our α CO measurements.Furthermore, as shown in the previous section, the correlation between Σ ⋆ and α CO improves after normalizing α CO the their W CO -weighted mean, which could fit into the picture of separating CO-dark and starburst components in Bolatto et al. (2013).With the WISE full-sky observations, the resolved Σ ⋆ for all nearby galaxies is widely available, which makes this kind of prescription easy to apply.In this study, we will focus on the α CO -to-Σ ⋆ relation instead of Σ Total because our data set is mostly Σ ⋆ -dominated (50% with Σ gas /Σ Total < 0.2; 99.5% with Σ gas /Σ Total < 0.5).
We present the correlations between α CO and Σ ⋆ in Figure 7 for both α CO (2−1) (left panels) and α CO (1−0) (right panels).In the top panels, we present the profile of measured α CO versus Σ ⋆ in each galaxy.For α CO (2−1) , most galaxies have their α CO anti-correlate with Σ ⋆ at Σ ⋆ > 100 M ⊙ pc −2 aside from a few exceptions.It is similar for α CO (1−0) , but The filled markers show the galaxies with WCO recovery fraction ≥ 50%, while the empty markers show the < 50% ones (Table 2).The correlation coefficients and p-values are labeled at the lower left in each panel, highlighting the significant (p-value< 0.05) ones.We use weighted averaged for αCO, U and sSFR, and simple averages for the other quantities.
Please see Section 4.2 for details.
with a flatter α CO -to-Σ ⋆ slope.At the low-Σ ⋆ end, some galaxies still have negative α CO -to-Σ ⋆ correlations while the others have strong positive correlations.In the middle panels, we show the collective behavior across galaxies using a binned average as a function of Σ ⋆ , with each bin spanning ∼ 0.1 dex in Σ ⋆ .We find that in regions with high-Σ ⋆ , α CO generally decreases with Σ ⋆ , which is consistent with the negative powerlaw index in the Bolatto et al. (2013) formula.There appears to be galaxy-to-galaxy variation in the value of α CO , but good agreement in the rate of how fast α CO decreases with Σ ⋆ .To better illustrate this phenomena, we normalize α CO in each galaxy at a threshold Σ ⋆,T ≡ 100 M ⊙ pc −2 (a threshold inspired by Bolatto et al. 2013) and show the normalized α CO in the bottom panel of Figure 7.The normalization in each galaxy (α CO,gal,T ) is defined as the median α CO of pixels with their Σ ⋆ within Σ ⋆,T ± 0.05 dex.
In the remainder of this section, we will focus on analyzing the scaling relation between α CO and Σ ⋆ in a sub-sample of galaxies with at least 5 measurements with Σ ⋆ > Σ ⋆,T (29 galaxies for CO (2−1) and 25 galaxies for CO (1−0), see the N pix,100 column in Table 2).We use a power law to characterize this scaling relation: where a (the power-law index) and b (the offset) are free parameters.Since both α CO and Σ ⋆ are normalized in the formula, we expect b ∼ 0 (and b gal ∼ 0) if α CO monotonically decreases with Σ ⋆ .By default, we fit Eq. 17 with all data.When fitting data in individual galaxies only, we will describe the parameters as a gal and b gal .We exclude data from galaxies that do not satisfy the following criteria: (1) at least 5 measurements at Σ ⋆ > Σ ⋆,T ; (2) spanning at least 0.1 dex in Σ ⋆ at Σ ⋆ > Σ ⋆,T .Since all criteria are Σ ⋆,T -dependent, we expect the size of sub-sample space to vary with Σ ⋆,T .We will visualize the galaxies not satisfying the last two criteria in Section 5. of a simple power law.The difference between the a values for CO (2−1) and CO (1−0) data is consistent with the finding that R 21 ∝ I −0.2 MIR (Leroy et al. 2023).The uncertainties of a and b are estimated from 1000 rounds of bootstrap resampling.In each round, we select 29 (25) galaxies for CO (2−1) (CO (1−0)) data with replacement from our sample galaxies to fit the power law.We then take the difference between the best-fit parameter and the 84th-(16th-) percentile from the 1000 bootstraps as the upper (lower) uncertainty.
We will measure how the α CO -to-Σ ⋆ relation varies according to the adopted D/M and Σ ⋆,T in the remainder of this section.In Appendix A, we also test how the results would change with an internal variations of D/M.Our toy model in Appendix A shows that a could be up to ∼ 0.2 steeper than the constant D/M case.

Dependence of the power-law index on adopted D/M and Σ ⋆,T
To test the robustness of our results for potentially different dust properties, we expand the assumed D/M from the single value (0.55) assumed in the previous section to the possible range of D/M, i.e. 0.1 ≤ D/M ≤ 1.We do not go to even lower D/M values because our methodology relies on the existence of certain amount of dust.We derive α CO and fit the α CO -to-Σ ⋆ power-law relation at each assumed D/M.The results are shown in Figure 8.We highlight the fit results with 0.4 ≤ D/M ≤ 0.7, which is the D/M value inferred from literature introduced in Section 2. Same as in the previous calculations, the uncertainties of parameters a and b in the fitting parameters are estimated from 1000 rounds of bootstrap resampling.
As shown in the top panel of Figure 8, we find that the power-law index (a) is invariant with assumed D/M throughout the range we examine for both α CO (2−1) and α CO (1−0) .The average a in the range of 0.4 ≤ D/M ≤ 0.7 is −0.48 +0.08  −0.09 and −0.22 +0.08  −0.09 for α CO (2−1) and α CO (1−0) , respectively.The statistical uncertainties in a for both CO transitions are around ±0.1 dex.The result implies that as long as the D/M stays roughly constant within each galaxy, we can recover similar behavior in the Σ ⋆ dependence of α CO with a power law Due to the nature of the definition of b in Eq. 17, we expect b ∼ 0. This is seen in most D/M values we examine as |b| stay below 0.03, shown in the middle panel of Figure 8.This indicates that the power law parameterization reasonably fits the observed data regardless of the assumed D/M value.However, the b values for α CO (2−1) seems biased toward the positive end.This could result from a steeper logα CO -to-logΣ ⋆ slope toward higher Σ ⋆ , which results in a positive offset in the power law at relatively lower Σ ⋆ .In the bottom panel, we show the ∆ rms value of each fit as an indicator of goodness of fit.All fits have ∆ rms below 0.2 dex, and the fits around 0.4 ≤ D/M ≤ 0.7 have ∆ rms ∼ 0.14 dex.
We further test if the chosen threshold in stellar mass surface density, Σ ⋆,T , will affect the fitting results.We fix Figure 7. Measured αCO plotted against stellar mass surface density.The left panels show the αCO(2−1) measurements while the right panels show the αCO(1−0) measurements.Top panels: the measured αCO in each galaxy, where the solid lines show the median in each Σ⋆ bin.Middle panels: the overall binned median of αCO, where the errorbars show the 16th-and 84th-percentiles.Bottom panels: similar to the middle panels, but the y-axis is normalized to the αCO value at Σ⋆,T = 100 M⊙ pc −2 (shown in vertical dashed line) in each galaxy.The galaxies with αCO increasing with Σ⋆ are removed in the bottom panels.The 2d histogram in all panels shows the overall data distribution.D/M = 0.55 and fit the power-law relation at Σ ⋆,T ranges from 30 to 300 M ⊙ pc −2 .Note that the number of galaxies included in the subsample changes in each case due to the threshold in Σ ⋆ .The results are shown in Figure 9.
In the top panel of Figure 9, we notice that the power-law index (a) has a larger dynamic range than the case where we alter D/M, but the index stays negative throughout the Σ ⋆,T range we examine.The power-law index for α CO (2−1) stays within ±0.1 of the fiducial case, and the indices for α CO (1−0) is consistent with the fiducial value at Σ ⋆,T > 60 M ⊙ pc −2 .There is a weak trend that |a| becomes larger toward larger Σ ⋆,T .The small b values indicate that the power-law function form applies in general.We also show the ∆ rms values in the bottom panel.We have ∆ rms ≤ 0.2, with a weak trend of smaller ∆ rms toward higher Σ ⋆,T .

Galaxy-to-galaxy variations
In this section, we examine the possible variation of the α CO -to-Σ ⋆ relation between individual galaxies, mainly how the variation in α CO,gal,T (normalization of α CO at Σ ⋆,T , see Section 5) and a gal (power-law index, Eq. 17) correlate with galaxy-averaged properties.By understanding what sets a gal and α CO,gal,T , we can build a prescription of α CO considering the α CO -to-Σ ⋆ relation and galaxy-to-galaxy variations.The results are visualized in Figure 10 and Figure 11.
In the upper panels of Figure 10 and Figure 11, we show how the power-law index (a gal ) varies with 7 selected galaxyaveraged properties and whether the galaxy is barred or not.The set of properties is the same as the ones in Figure 5 and Figure 6.None of the properties show significant correlation with a gal .Meanwhile, the standard deviation of a gal is 0.30 for both CO (2−1) and CO (1−0).
In the lower panels of Figure 10 and Figure 11, we show how the normalization in each galaxy (α CO,gal,T ) varies with galaxy-averaged properties.The standard deviation of α CO,gal,T is 0.2 dex for both CO (2−1) and CO (1−0).For CO (2−1), <W CO >, logU gal , <Σ SFR >, and sSFR show significant correlations with α CO,gal,T with similar strength.For CO (1−0), <12 + log(O/H)>, logU gal , <Σ ⋆ >, <Σ SFR > and sSFR show significant correlations with α CO,gal,T with similar strength.We use these significant correlations to fit empirical relations for α CO,gal,T and summarize the results in Table 4.The fitted empirical relations do not differ significantly in terms of ∆ rms .Meanwhile, the fit with U has the smallest 12 the statistical uncertainties in the fitted parameters among the parameters for both CO (2−1) and CO (1−0).Besides U, <Σ SFR > also has small statistical uncertainties and is available for more galaxies.

DISCUSSION
12 Considering the product of log x and the uncertainties in m.

General suggestions for α CO prescriptions
In Section 4.1, we present how the measured α CO correlates with local physical quantities and provide linear regression for each quantity at 2 kpc scale in Table 2.These measurements consider the statistical behavior of the overall sample.Among the quantities, W CO , U, and Σ SFR usually have the strongest correlations with α CO and smallest ∆ rms from linear regression.We would suggest the readers go with these parameters if the parameter space of their sample overlaps with this study (see Figure 1 for the completeness of each quantity).
Meanwhile, we explicitly explore the relation between α CO and Σ ⋆ in Section 5 as a possible tracer for starburst α CO at 2 kpc scale, and find strong correlation between α CO and Σ ⋆ after normalization at some Σ ⋆ threshold.There are two ways to adopt these results.The first one is a stand-alone prescription combining the indices from Section 5.1 and normalization from Table 4, using galaxy-averaged U gal as an example: where α CO is in unit of M ⊙ pc −2 (K km s −1 ) −1 , and the U galdependent normalization could be replaced with other quantities listed in Table 4, e.g.<Σ SFR >.Please refer to Section 5 for relevant uncertainties.We note that the possible variation of the power-law index could be up to ∼ 0.2 due to internal variations of D/M (Appendix A).
On the other hand, one key mechanism that sets α CO , the CO-dark gas, is likely not parameterized by our formula (see Section 6.2).This is because the CO-dark gas effect is relatively weak in the metallicity span of our sample.However, both the CO-dark gas and the "starburst α CO " effects should be considered for an α CO prescription to be applied through all environments.Thus, another suggestion we have is to make a Bolatto et al. (2013)-style combination (also see E. Schinnerer & A. K. Leroy ARA&A submitted) of our Σ ⋆ power-law term with existing α CO prescriptions tracing the CO-dark gas effect, e.g.Wolfire et al. (2010), Narayanan et al. (2012), Schruba et al. (2012), Hunt et al. (2015), Accurso et al. (2017), or Sun et al. (2020).That is, assuming the adopted existing CO-dark prescription is α CO-dark CO , we suggest: Under this functional form, we expect the normalization (α CO,gal,T ) to be taken into account by the α CO-dark CO term.Since this formula does not include our own normalization, the Σ ⋆ power-law can be extended to lower Σ ⋆ threshold.For  α CO (1−0) case, one can simply replace the power-law index with −0.22.
One of the future directions we will take is to study how α CO correlates with physical quantities at cloud scales instead of kpc scales and build α CO prescriptions accordingly.The advantage in this direction is that the physical quantities at cloud scale are more strongly linked to fundamental physics of CO emission and dynamics of molecular clouds.For instance, based on the α CO measurements in this work, Teng et al. (2023a) have reported α CO dependence with cloud-scale velocity dispersion which likely traces CO opacity change.We also refer the readers to Teng et al. (2022Teng et al. ( , 2023b) ) for more details.
The other possible future direction for dust-based α CO is a new strategy that simultaneously allows variations in α CO , dust properties (e.g.D/M in Appendix A and dust opacity in Section 3.1) and metallicity at best-possible resolution.The Leroy et al. (2011) and Sandstrom et al. (2013, also see their Appendix A) strategy is a good demonstration of the concept for most items on the list except the resolution.A more sophisticated strategy would help identify the next step forward on α CO prescriptions.We are also interested in investigating whether the Σ ⋆ -dependence still applies to Σ gas -dominated environments.

Interpreting of the environmental dependence of α CO
In this section, we will discuss the physical interpretations of the correlations between α CO and the physical quantities we present in the previous sections.As we mentioned in Section 1, we expect two main trends in the variation of α CO : (1) the "CO-dark gas" trend, where α CO increases toward lower metallicity as shielding for CO weakens; (2) the "starburst conversion factor" trend, where α CO decreases toward galaxy centers and U/LIRGs with the decrease in CO optical depth or increase in CO excitation.
Regarding the CO-dark gas trend, we observe moderate to weak anti-correlation between α CO and 12 + log(O/H) at kpc scale (Section 4.1), and moderate anti-correlation at galaxy scale (Section 4.2).One possible explanation for the weak correlation is that the statistical significance becomes weaker with the small dynamic range of our 12 + log(O/H) data: 80% of the 12 + log(O/H) measurements fall within a 0.2 dex range from 8.4 to 8.6.The α CO -Σ ⋆ relation we present in Section 5 is unlikely caused by this CO-dark gas effect since the dynamic range in 12 + log(O/H) is even smaller for data above the Σ ⋆,T threshold.Another explanation is that the CO-dark gas effect is weaker at nearly solar metallicity (e.g.Wolfire et al. 2010;Glover & Mac Low 2011;Hunt et al. 2015).However, we note that recent simulations show that there is significant fraction of CO-dark gas ( f dark ) up to solar metallicity, e.g.Gong et al. (2018) found f dark ranges from 26%-79%.These studies found that f dark correlates with ex-tinction and/or W CO (Smith et al. 2014;Gong et al. 2018Gong et al. , 2020;;Hu et al. 2022).
We interpret U and (r)sSFR as empirical tracers for regions with high SFR, where the "starburst conversion factor" trend matters and lowers down α CO .Some studies have also argued that α CO could decrease with increased radiation field due to CO dissociation (Israel 1997b;Wolfire et al. 2010;Accurso et al. 2017).However, we did not observe this trend, and one possible explanation is that the CO dissociation effect should be weak as long as the CO emission is optically thick (Wolfire et al. 2010;Bolatto et al. 2013).Σ SFR , which simultaneously traces the UV radiation and starburst regions, also have moderate anti-correlation with α CO across all cases.In general, we observe a stronger correlation between Σ SFR and α CO (2−1) than α CO (1−0) .We also observe moderate anticorrelation between W CO and α CO .This is consistent with the theoretical assumption under optically thick assumption Bolatto et al. (2013) and recent observations (Hunt et al. 2023).
We interpret the α CO -to-Σ ⋆ anti-correlation as the increase of velocity dispersion of molecular gas from additional gravity.Bolatto et al. (2013, also see Hirashita (2023)) suggested that in high-Σ ⋆ environments, the molecular gas experiences gravitational potential from stellar sources, ending up with a total pressure larger than isolated, virialized clouds.This larger pressure results in gas line width larger than one would expect from self-gravitating clouds.This increase in gas line width scales with total mass (stars and gas) in the system, or α CO ∝ M mol /(M mol +M ⋆ ) 0.5 .The above functional form approximates a α CO -to-Σ ⋆ power law in M ⋆ -dominated regions.
For the argument to hold, the CO emission must be optically thick.Bolatto et al. (2013) mentioned that the only possible structure for molecular gas that satisfies this scenario is an extended molecular medium.

Comparing to previous α CO surveys
In this section, we compare our measurements to α CO values obtained in previous dust-based α CO surveys.Since we will cover studies with both α CO (2−1) and α CO (1−0) measurements, we will convert α CO (2−1) values in literature and this work to α CO (1−0) with R 21 = 0.65 for simplicity and uniformity.First, we compare our α CO maps with Sandstrom et al. (2013).They measured spatially resolved α CO in 26 nearby, star-forming galaxies using a dust-based methodology (also see Leroy et al. 2011;den Brok et al. 2023;Yasuda et al. 2023).They assume that the variation of D/G is minimal in a few-kpc-scale "solution pixel" consisting of 37 samples in a hexagonal region, and fit D/G and α CO simultaneously from data by minimizing the variation in D/G.Sandstrom et al. (2013) used CO (2−1) data from HERACLES (Leroy et al. 2009a) and R 21 = 0.7.We rescale their results with R 21 = 0.65 for uniformity.Compared to our work, the Sandstrom et al. (2013) metholodogy has larger degrees of freedom for the spatial variation of D/G and D/M.Meanwhile, it is more difficult to push their methodology to a larger sample size at fixed physical resolution.
There are 13 galaxies that are studied in both Sandstrom et al. (2013, see their Figure 7) and this work.We show the α CO measurements from both works as a function of galactocentric radius (R g in terms of R 25 .) in Figure 12.We adopt R 25 values in this work instead of the Sandstrom et al. (2013) values.The simple mean α CO of all measurements in Sandstrom et al. ( 2013) is ∼ 2.3 M ⊙ pc −2 (K km s −1 ) −1 .Similar with Sandstrom et al. (2013), we find a weak to moderate (but significant) positive correlation between α CO and R g .When we normalize the α CO in each galaxy by the mean α CO in each galaxy (<α CO, gal >), both works show a flat trend with radius in mid-to outer-disk.Sandstrom et al. (2013) data shows a more significant decrease in α CO in the galaxy center.Several factors could contribute to the difference in the inner most radial bins.If we calculate the W CO -weighted mean instead of the median in each bin, the difference in the bin with smallest radii will decrease by 0.1 dex, which partially explains the discrepancy.Another possible explanation is that some of the measurements with small Σ mol (and possibly small α CO ) are removed due to small S/N; however, they are taken into account in Sandstrom et al. (2013).It is not clear whether the difference in resolution is a cause.When we calculate α CO at 1 kpc resolution with the galaxies with distance within 10 Mpc, there is no clear trend of the resulting α CO with resolution.The fixed D/M is unlikely to be a major cause since the Sandstrom et al. (2013) results are consistent with a D/G-metallicity power law, see their Figure 13.
For comparing galaxy-averaged α CO measurements, we include another previous study: the COMING survey (Sorai et al. 2019;Yasuda et al. 2023).The COMING survey solves D/G and α CO (1−0) simultaneously by minimizing a χ 2 value defined by the difference between D/G × (Σ atom + α CO (1−0) W CO (1−0) ) and Σ dust derived from dust SED fitting.Here, we quote their "global" result, where the authors fit all data within one galaxy to retrieve one set of D/G and α CO values.
We compare our measured α CO in each galaxy with Sandstrom et al. ( 2013) and the COMING survey (Sorai et al. 2019;Yasuda et al. 2023) in Figure 13.For Sandstrom et al. (2013) and this work, we adopt the W CO -weighted mean.For the COMING survey, we adopt their "global" result.The α CO measured in the three works are in general consistent with each other within uncertainties.Our measurements made with CO (2−1) and CO (1−0) agree with each other.When there is a difference, it is more often that the one derived with CO (2−1) has a slightly larger value.We also note that in several galaxies with signatures of active galactic nucleus (AGN; see classification in Kennicutt et al. 2011), there is larger offset between our measurements and literature, e.g.(2013).For both works, we display αCO = 0.65αCO(2−1) for uniformity.We only include measurements from the galaxies that are included in both works.The vertical dashed line shows the completeness threshold in Rg/R25 in this work.In the top panel, the horizontal cyan dashed-dotted line shows the mean αCO in Sandstrom et al. (2013), and the orange dashed line shows the mean αCO from this work.The two lines are close to each other.In the bottom panel, the horizontal line shows αCO = <αCO, gal>, where <αCO, gal> is the mean αCO in each galaxy.
NGC3627 and NGC4725; however, there are also galaxies with AGN show consistent results, e.g.NGC4536, NGC4736 and NGC5055.Thus, having AGN is not the only cause for the mismatch, and it is likely that the type of nuclei activities do not dominate the kpc-scale α CO values (e.g.Sandstrom et al. 2013).The adopted dust SED fitting method is also unlikely the cause for the difference in NGC4725 since we have a lower estimate of Σ dust , which should yield smaller α CO .
Our measurements made with CO (1−0) generally agree with the COMING survey.
We examine how α CO scales with several physical quantities, i.e.W CO , metallicity, Σ dust , ISRF, Σ ⋆ , Σ SFR and (r)sSFR.At 2 kpc scale, all quantities have significant local correlation  2), and half-filled symbols show literature values.We only include galaxies that are measured in at least one of the literature survey.The errorbar for this work shows the 16th-and 84th-percentiles (Table 2).The mean and errorbar of previous works are quoted from Table 4 of Sandstrom et al. (2013, with rescaling for R21) and Table 3 of Yasuda et al. (2023).The mean values are WCO-weighted mean for this work and Sandstrom et al. (2013), and the "global" result for Yasuda et al. (2023).with α CO .Among them, the strength of the ISRF (U), Σ SFR and W CO have the strongest anti-correlation with spatially resolved α CO .We provide linear regression of α CO with all the quantities tested, along with the corresponding performance and uncertainties in Table 3.
At galaxy-integrated scale, most quantities have significant correlation with W CO -weighted mean α CO .U, Σ SFR , W CO and 12 + log(O/H) have significant correlations with α CO for both CO (1−0) and CO (2−1) cases.
When we normalized resolved α CO measurements by the W CO -weighted mean in each galaxy, we find an increased correlation strength between normalized α CO and Σ ⋆ .After examining through Σ ⋆ bins, we find that in regions with high stellar mass surface densities (Σ ⋆ ≥ 100 M ⊙ pc −2 ), the α CO decreases with Σ ⋆ .Specifically, we find: It also has little dependence on the adopted ratio between CO rotational lines.
When fitting the power-law relation within individual galaxies, we find significant dependence of the normalization of α CO in each galaxy on several quantities.Among them, the linear regression to logU gal has the minimal statistical uncertainties.Thus, we recommend using Σ ⋆ and logU gal to predict α CO at high-Σ ⋆ environments.
This decrease of α CO in the high-Σ ⋆ region is likely due to the increased CO brightness with increased line width.The line width is larger than self-gravitating clouds due to the additional gravity from stellar sources, and the structure satisfying this scenario is likely an extended molecular medium.Understanding the decrease in α CO at high Σ ⋆ is important for accurately assessing molecular gas content and star-formation efficiency in the centers of galaxies and bridges the "MW-like" to "starburst" conversion factor.tract NAS5-26555.ADB acknowledges support from the National Science Foundation through grants AST-2108140 and AST-2307441.EWK acknowledges support from the Smithsonian Institution as a Submillimeter Array (SMA) Fellow and the Natural Sciences and Engineering Research Council of Canada.JC acknowledges support from ERC starting grant #851622 DustOrigin This work uses observations made with ESA Herschel Space Observatory.Herschel is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA.The Herschel spacecraft was designed, built, tested, and launched under a contract to ESA managed by the Herschel/Planck Project team by an industrial consortium under the overall responsibility of the prime contractor Thales Alenia Space (Cannes), and including Astrium (Friedrichshafen) responsible for the payload module and for system testing at spacecraft level, Thales Alenia Space (Turin) responsible for the service module, and Astrium (Toulouse) responsible for the telescope, with in excess of a hundred subcontractors.
This paper makes use of the VLA data with project codes 14A-468, 14B-396, 16A-275 and 17A-073, which has been processed as part of the EveryTHINGS survey.This paper makes use of the VLA data with legacy ID AU157, which has been processed in the PHANGS-VLA survey.The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.This publication makes use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration.This research made use of Astropy,13 a communitydeveloped core Python package for Astronomy (Astropy Collaboration et al. 2013, 2018, 2022).This research has made CO-intensity-weighted mean.(b)  The percentiles are calculated with non-weighted data.(c) Fraction of WCO recovered (above the S/N mask) in each galaxy.Galaxies with CO recovery fraction < 50% will be visualized differently in figures showing galaxy-averaged values.(d) Number of pixel-by-pixel measurements.(e) Number of pixel-by-pixel with valid αCO measurements at Σ⋆ ≥ 100 M⊙ pc −2 .This will be discussed later in Section 5. ( †) αCO(1−0) calculated with R21(IW4) (Eq. 6) and WCO(2−1) due to no WCO(1−0) data.

Figure 1 .
Figure1.The completeness of our data set at each αCO-quantity pair.The 50% completeness is marked with a horizontal black line.In the statistical calculations, e.g.linear regression and correlation, we only use data in the parameter range with completeness > 50%.

Figure 2 .
Figure 2. The distribution of measured αCO(2−1) and αCO(1−0).The mean value of each type of measurement is marked in vertical lines with the corresponding color.

Figure 3 .
Figure 3.Our measured αCO(2−1), measured with WCO(2−1) data only, as a function of environmental parameters.The orange circles show the median of binned data.The errorbars include both the scatter within a bin and uncertainties of pixel-by-pixel measurements.The orange line shows the linear regression between αCO(2−1) and the quantity on the horizontal axis.The empty circles indicate bins where the quantity in the horizontal axis is incomplete, i.e. less than 50% of the pixels have S/N > 1 in WCO and Herschel bands (see Section 3.3).The dashed vertical line shows the 50% boundary.The gray shaded region shows the hexagonal-binned pixel-by-pixel measurements.The gray horizontal line shows the value of α MW CO , converted to αCO(2−1) assuming R21 = 0.65.

Figure 5 .
Figure5.Scaling relations between galaxy-averaged αCO(2−1) and environmental parameters.The errorbars show the 16th-84th percentiles in each galaxy.All data points are calculated with WCO(2−1) data only.The filled markers show the galaxies with WCO recovery fraction ≥ 50%, while the empty markers show the < 50% ones (Table2).The correlation coefficients and p-values are labeled at the lower left in each panel, highlighting the significant (p-value< 0.05) ones.We use weighted averaged for αCO, U and sSFR, and simple averages for the other quantities.Please see Section 4.2 for details.

Figure 6 .
Figure 6.Similar with Figure 5 but for αCO(1−0) measurements.The orange points show the measurements with WCO(1−0) data only, and the purple points show measurements with WCO(2−1) data with a IW4-dependent R21.The fixed R21 results are similar to the ones from IW4-dependent R21.They are not displayed for clarity.

Figure 8 .
Figure 8.The dependence of fitted power-law index and offset on assumed D/M.The cyan points show results for αCO(2−1), and the orange points show results for αCO(1−0).For αCO(1−0), we only show results from WCO(1−0) data only since fitting results including WCO(2−1) data have mininal difference.We highlight the region with D/M inferred from literature, i.e. 0.4 ≤ D/M ≤ 0.7, while the empty circles show results outside that range.Top: the power-law index (a).The dashed lines show the mean value of a in the range of 0.4 ≤ D/M ≤ 0.7.Middle: the offset (b).The dotted line shows b = 0, where expect the fitting result to be if αCO monotonically decreases with Σ⋆.Bottom: the ∆rms of each fit.

Figure 9 .
Figure 9.The dependence of fitted power-law index and offset on adopted threshold Σ⋆,T.The cyan points show results for αCO(2−1), and the orange points show results for αCO(1−0).For αCO(1−0), we only show results from WCO(1−0) data only since fitting results including WCO(2−1) data have mininal difference.Top: the powerlaw index (a).The dashed lines show the average a values with Σ⋆,T = 100 M⊙ pc −2 in the range of 0.4 ≤ D/M ≤ 0.7.These are the same lines as Figure 8. Middle: the offset (b).The dotted line shows b = 0, where expect the fitting result to be if αCO monotonically decreases with Σ⋆.Bottom: the ∆rms of each fit.All the calculation here are done with D/M = 0.55.

Figure 10 .
Figure10.Scaling relations between agal (the power-law index), αCO,gal,T and galaxy-averaged environmental parameters for CO (2−1) data.Galaxies that are barred and not barred are colored in orange and cyan, respectively.The black dashed line in the agal show the a value calculated with overall data.The correlation coefficients and p-values are labeled at the lower left in each panel, highlighting the significant (p-value< 0.05) ones.

Figure 12 .
Figure12.The relation between measured αCO and galacto-centric radius in this work andSandstrom et al. (2013).For both works, we display αCO = 0.65αCO(2−1) for uniformity.We only include measurements from the galaxies that are included in both works.The vertical dashed line shows the completeness threshold in Rg/R25 in this work.In the top panel, the horizontal cyan dashed-dotted line shows the mean αCO inSandstrom et al. (2013), and the orange dashed line shows the mean αCO from this work.The two lines are close to each other.In the bottom panel, the horizontal line shows αCO = <αCO, gal>, where <αCO, gal> is the mean αCO in each galaxy.

Figure 13 .
Figure13.The mean αCO values in each galaxy from this work (both CO (2−1) data with R21 = 0.7 and CO (1−0) data),Sandstrom et al. (2013,  S13)  and the COMING survey(Sorai et al. 2019; Yasuda et al. 2023, Y23).Circles show αCO values derived with CO (2−1) data, and triangles show αCO values derived with CO (1−0) data.Filled symbols show the results from this work, empty symbols show the ones with low CO recovery fraction (Table2), and half-filled symbols show literature values.We only include galaxies that are measured in at least one of the literature survey.The errorbar for this work shows the 16th-and 84th-percentiles (Table2).The mean and errorbar of previous works are quoted from Table4ofSandstrom et al. (2013, with rescaling for R21) and Table3ofYasuda et al. (2023).The mean values are WCO-weighted mean for this work andSandstrom et al. (2013), and the "global" result forYasuda et al. (2023).

Table 1
(continued) 5 (seeSánchez et al. 2019, and  references therein).aandb are free parameters.We fit the function with 12 + log(O/H) at one R e from galaxies with the available measurements listed in Table1.The best-fit parameters are a = 8.56 ± 0.02 and b = 0.010 ± 0.002, which are robust under resampling.Meanwhile, the typical statistical uncertainty in the 12 + log(O/H) data used for fitting is ∼ 0.013 dex.However, the root mean square error (∆ rms )

Table 2 .
Statistics of αCO measurements.All the αCO values are in unit of [M⊙ pc −2

Table 3 .
The correlation and linear regression between pixel-bypixel αCO measurements and local physical quantities.
(1) The linear regression formula is log αCO = mx + b.An uncertainty of ±0.0 represents that the rounded uncertainty of the parameter is smaller than 0.01.(2) Correlation coefficients calculated with αCO normalized by WCO-weighted mean value in each galaxy.We underline the cases where the correlation of normalized αCO is stronger than the one without normalization.

Table 4 .
Dependence of αCO,gal,T on galaxy-integrated quantities.The linear regression formula is log αCO,gal,T = m log x + d for most quantities and log αCO,gal,T