# Low-J CO Line Ratios from Single-dish CO Mapping Surveys and PHANGS-ALMA

Adam K. Leroy, Erik Rosolowsky, Antonio Usero, Karin Sandstrom, Eva Schinnerer, Andreas Schruba, Alberto D. Bolatto, Jiayi Sun, Ashley. T. Barnes, Francesco Belfiore, et al.

► To cite this version:

Adam K. Leroy, Erik Rosolowsky, Antonio Usero, Karin Sandstrom, Eva Schinnerer, et al. Low-J CO Line Ratios from Single-dish CO Mapping Surveys and PHANGS-ALMA. The Astrophysical Journal, 2022, 927, 10.3847/1538-4357/ac3490. insu-03667457

HAL Id: insu-03667457

https://insu.hal.science/insu-03667457

Submitted on 14 May 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Low-J CO Line Ratios from Single-dish CO Mapping Surveys and PHANGS-ALMA

Adam K. Leroy<sup>1,2</sup>, Erik Rosolowsky<sup>3</sup>, Antonio Usero<sup>4</sup>, Karin Sandstrom<sup>5</sup>, Eva Schinnerer<sup>6</sup>, Andreas Schruba<sup>7</sup>,
Alberto D. Bolatto<sup>8</sup>, Jiayi Sun (<u>孙嘉懿</u>)<sup>1</sup>, Ashley. T. Barnes<sup>9</sup>, Francesco Belfiore<sup>10</sup>, Frank Bigiel<sup>9</sup>,
Jakob S. den Brok<sup>9</sup>, Yixian Cao<sup>11</sup>, I-Da Chiang (<u>江宜達</u>)<sup>5</sup>, Mélanie Chevance<sup>12</sup>, Daniel A. Dale<sup>13</sup>,
Cosima Eibensteiner<sup>9</sup>, Christopher M. Faesi<sup>14</sup>, Simon C. O. Glover<sup>15</sup>, Annie Hughes<sup>16</sup>, María J. Jiménez Donaire<sup>4,17</sup>,
Ralf S. Klessen<sup>15,18</sup>, Eric W. Koch<sup>19</sup>, J. M. Diederik Kruijssen<sup>12</sup>, Daizhong Liu<sup>7</sup>, Sharon E. Meidt<sup>20</sup>, Hsi-An Pan<sup>6,21</sup>,
Jérôme Pety<sup>22,23</sup>, Johannes Puschnig<sup>9</sup>, Miguel Querejeta<sup>4</sup>, Toshiki Saito<sup>6</sup>, Amy Sardone<sup>1,2</sup>, Elizabeth J. Watkins<sup>12</sup>,
Axel Weiss<sup>24</sup>, and Thomas G. Williams<sup>6</sup>

<sup>1</sup> Department of Astronomy, The Ohio State University, 140 West 18th Avenue, Columbus, OH 43210, USA; leroy.42@osu.edu

<sup>2</sup> Center for Cosmology and Astroparticle Physics, 191 West Woodruff Avenue, Columbus, OH 43210, USA

<sup>3</sup> Department of Physics, University of Alberta, Edmonton, AB, T6G 2E1, Canada

<sup>4</sup> Observatorio Astronómico Nacional (IGN), C/Alfonso XII, 3, E-28014 Madrid, Spain

<sup>5</sup> Center for Astrophysics and Space Sciences, Department of Physics, University of California, San Diego, 9500 Gilman Drive, La Jolla, CA 92093, USA

<sup>6</sup> Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117, Heidelberg, Germany

<sup>7</sup> Max-Planck-Institut für extraterrestrische Physik, Giessenbachstraße 1, D-85748 Garching, Germany

<sup>8</sup> Department of Astronomy, University of Maryland, College Park, MD 20742, USA

<sup>9</sup> Argelander-Institut für Astronomie, Universität Bonn, Auf dem Hügel 71, D-53121 Bonn, Germany

<sup>10</sup> INAF—Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, I-50157, Firenze, Italy

<sup>11</sup> Aix Marseille Université, CNRS, CNES, LAM (Laboratoire d'Astrophysique de Marseille), F-13388 Marseille, France

<sup>12</sup> Astronomisches Rechen-Institut, Zentrum für Astronomie der Universität Heidelberg, Mönchhofstraße 12-14, D-69120 Heidelberg, Germany

<sup>13</sup> Department of Physics and Astronomy, University of Wyoming, Laramie, WY 82071, USA

<sup>14</sup> University of Massachusetts–Amherst, 710 North Pleasant Street, Amherst, MA 01003, USA

<sup>15</sup> Universität Heidelberg, Zentrum für Astronomie, Institut für Theoretische Astrophysik, Albert-Ueberle-Str 2, D-69120 Heidelberg, Germany

<sup>16</sup> Université de Toulouse, UPS-OMP, IRAP, F-31028 Toulouse cedex 4, France

<sup>17</sup> Centro de Desarrollos Technológicos, Observatorio de Yebes (IGN), E-19141 Yebes, Guadalajara, Spain

<sup>18</sup> Universität Heidelberg, Interdisziplinäres Zentrum für Wissenschaftliches Rechnen, Im Neuenheimer Feld 205, D-69120 Heidelberg, Germany

<sup>19</sup> Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA

<sup>20</sup> Sterrenkundig Observatorium, Universiteit Gent, Krijgslaan 281 S9, B-9000 Gent, Belgium

<sup>21</sup> Department of Physics, Tamkang University, No. 151, Yingzhuan Road, Tamsui District, New Taipei City 251301, Taiwan

<sup>22</sup> Institut de Radioastronomie Millimétrique (IRAM), 300 Rue de la Piscine, F-38406 Saint Martin d'Hères, France

<sup>23</sup> Sorbonne Université, Observatoire de Paris, Université PSL, CNRS, LERMA, F-75014, Paris, France

<sup>24</sup> Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, D-53121 Bonn, Germany

Received 2021 July 22; revised 2021 September 22; accepted 2021 September 23; published 2022 March 10

Abstract

We measure the low-J CO line ratios $R_{21} \equiv \text{CO}(2-1)/\text{CO}(1-0)$, $R_{32} \equiv \text{CO}(3-2)/\text{CO}(2-1)$, and $R_{31} \equiv \text{CO}(3-2)/\text{CO}(1-0)$ using whole-disk CO maps of nearby galaxies. We draw $\text{CO}(2-1)$ from PHANGS-ALMA, HERACLES, and follow-up IRAM surveys; $\text{CO}(1-0)$ from COMING and the Nobeyama CO Atlas of Nearby Spiral Galaxies; and $\text{CO}(3-2)$ from the James Clerk Maxwell Telescope Nearby Galaxy Legacy Survey and Atacama Pathfinder Experiment Large APEX Sub-Millimetre Array mapping. All together, this yields 76, 47, and 29 maps of $R_{21}$, $R_{32}$, and $R_{31}$ at $20'' \sim 1.3$ kpc resolution, covering 43, 34, and 20 galaxies. Disk galaxies with high stellar mass, $\log(M_*/M_{\odot}) = 10.25-11$, and star formation rate (SFR) = $1-5 M_{\odot} \text{ yr}^{-1}$, dominate the sample. We find galaxy-integrated mean values and a 16%–84% range of $R_{21} = 0.65$ (0.50–0.83), $R_{32} = 0.50$ (0.23–0.59), and $R_{31} = 0.31$ (0.20–0.42). We identify weak trends relating galaxy-integrated line ratios to properties expected to correlate with excitation, including SFR/$M_*$ and SFR/$L_{\text{CO}}$. Within galaxies, we measure central enhancements with respect to the galaxy-averaged value of $\sim 0.18_{-0.14}^{+0.09}$ dex for $R_{21}$, $0.27_{-0.15}^{+0.13}$ dex for $R_{31}$, and $0.08_{-0.09}^{+0.11}$ dex for $R_{32}$. All three line ratios anticorrelate with galactocentric radius and positively correlate with the local SFR surface density and specific SFR, and we provide approximate fits to these relations. The observed ratios can be reasonably reproduced by models with low temperature, moderate opacity, and moderate densities, in good agreement with expectations for the cold interstellar medium. Because the line ratios are expected to anticorrelate with the $\text{CO}(1-0)$-to-$H_2$ conversion factor, $\alpha_{\text{CO}}^{1-0}$, these results have general implications for the interpretation of CO emission from galaxies.

Unified Astronomy Thesaurus concepts: Galaxies (573); Interstellar molecules (849); Interstellar medium (847); Molecular gas (1073)

Supporting material: machine-readable tables

1. Introduction

Rotational line emission from carbon monoxide (CO) represents the main way to trace the distribution, kinematics, and physical conditions in the molecular interstellar medium (ISM) in external galaxies (e.g., see reviews by Bolatto et al.

<figure><img src="image_2.png" /></figure>

Original content from this work may be used under the terms
of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further
distribution of this work must maintain attribution to the author(s) and the title
of the work, journal citation and DOI.
2013a; Klessen & Glover 2016). After several decades focused primarily on the fundamental ¹²CO (1–0) transition at ν≈115 GHz (e.g., Young & Scoville 1991; Young et al. 1995; Helfer et al. 2003), improvements in (sub)millimeter facilities over the last 15 yr have enabled extensive mapping of nearby galaxies in ¹²CO (2–1) and ¹²CO (3–2) at ν≈231 and 345 GHz (e.g., Leroy et al. 2009; Wilson et al. 2012). In the last decade, the Atacama Large Millimeter/submillimeter Array (ALMA) has come online and revolutionized our view of molecular line emission from galaxies while also accelerating the trend toward observing multiple CO lines. Thanks to its excellent site and submillimeter sensitivity, ALMA can often map CO (2–1) and CO (3–2) several times faster than CO (1–0) at matched resolution and sensitivity. As a result, ALMA surveys of nearby galaxies have targeted all three low-$J$ CO lines: CO (3–2), CO (2–1), and CO (1–0) (e.g., García-Burillo et al. 2014; Hirota et al. 2018; Leroy et al. 2021a).

Meanwhile, studies of redshifted CO emission have also become common, tracing the molecular gas at earlier cosmic epochs. Driven by similar technical considerations, these studies currently focus on CO (3–2), CO (2–1), or even higher $J$ transitions (e.g., see reviews by Carilli & Walter 2013; Hodge & da Cunha 2020; Tacconi et al. 2020). In the near future, observations at high redshift may become even more diverse, as the proposed next-generation Very Large Array (Murphy et al. 2018) will vastly improve our ability to observe CO (1–0) emission at intermediate and high redshift.

This increased diversity of CO line observations at low and high $z$ makes the ability to translate between results obtained using these different CO lines crucial. Despite the proliferation of CO (2–1) and CO (3–2) studies, many surveys still target CO (1–0), including xCOLD GASS (Saintonge et al. 2017) and CARMA EDGE (Bolatto et al. 2017), the largest low-$z$ single-dish and interferometric CO surveys to date. Critical work informing our interpretation of CO emission has also built on observations of a single transition; e.g., Donovan Meyer et al. (2013) focused on CO (1–0) emission, Sandstrom et al. (2013) studied CO (2–1), and Wilson et al. (2008) employed CO (3–2). Well-understood, observationally tested translations between the various low-$J$ CO lines are required to link these efforts.

Indeed, translations between the different transitions are not straightforward because the ratios among CO (3–2), CO (2–1), and CO (1–0) also reflect physical conditions in the molecular gas. The observed ratios emerge from an interplay among the distributions of collider density; kinetic temperature, $T_{kin}$; and column density per line width (see Section 2). These, in turn, depend on the structure, kinematics, and heating mechanisms at play in the cold ISM. The ratios of low-$J$ CO lines thus represent a potentially powerful observational probe of the local physical conditions in the molecular ISM. This potential is complicated by degeneracies in their interpretation and the modest dynamic range in their observed values. This limited dynamic range places relatively strict requirements on observations aiming to measure these line ratios.

In contrast to commonly used “dense gas tracers” like HCN (1–0) and HCO$^+$ (1–0), the CO lines are bright and can be studied across a range of environments (see Usero et al. 2015, regarding relative line strengths). Numerical simulations can now resolve CO chemistry and predict CO line emission over whole molecular clouds, large parts of a spiral galaxy, or even entire dwarf galaxies (e.g., Glover & Clark 2012; Peñaloza

et al. 2017, 2018; Gong et al. 2020; Hu et al. 2021), but such calculations remain extremely challenging for tracers of higher-density gas (e.g., Onus et al. 2018). A combined observational, numerical, and analytical approach that leverages ratios among the low-$J$ CO lines and their isotopologues represents a promising path forward to diagnose physical conditions in the molecular gas of galaxies. This approach can become even more powerful when paired with high-resolution imaging of the CO emission (e.g., see Gallagher et al. 2018b), which places constraints on the mean density and kinematics of the cold gas (e.g., see Sun et al. 2018, 2020; Rosolowsky et al. 2021). Of course, spectroscopy targeting multiple CO transitions or isotopologues has a long history (e.g., Papadopoulos & Seaquist 1999; Israel & Baas 2001, 2003; Bayet et al. 2004, 2006; Kamenetzky et al. 2014, Israel et al. 2015; 2017). However, most previous work has focused on single-pointing or galaxy-integrated measurements, with a heavy emphasis on galaxy centers and starburst galaxies, including many ultraluminous and luminous infrared galaxies. Resolved studies that measure the ratios among multiple CO lines over the full area of “normal” star-forming main-sequence galaxies remain relatively scarce.

This paper presents new measurements of the CO (2–1)/CO (1–0), CO (3–2)/CO (2–1), and CO (3–2)/CO (1–0) line ratios for nearby galaxies ($D < 40$ Mpc, median ~14 Mpc) based on maps of CO emission from Kuno et al. (2007), HERACLES (Leroy et al. 2009), the James Clerk Maxwell Telescope (JCMT) Nearby Galaxy Legacy Survey (NGLS; Wilson et al. 2012), the CO Multiline Imaging of Nearby Galaxies (COMING) survey (Sorai et al. 2019), PHANGS-ALMA (Leroy et al. 2021a), new IRAM 30 m CO (2–1) observations (PI: A. Schruba), and new Atacama Pathfinder Experiment (APEX) Large APEX Sub-Millimetre Array (LASMA) CO (3–2) observations (PI: A. Weiss). We measure both resolved and integrated CO line ratios using mapping surveys. All of these surveys except PHANGS-ALMA use receiver arrays on single-dish telescopes to cover large areas quickly (e.g., see Schuster et al. 2007). Restricting our focus to mapping data allows us to construct identical matched apertures when measuring integrated ratios. This avoids the common issue of mismatched beams, which plagued some earlier studies that relied on pointed observations. These surveys also target many of the largest, closest, best-studied galaxies, so the ratios for individual targets are of particular interest. Finally, because we analyze maps, we can measure line ratios associated with distinct regions to, e.g., test for a dependence of excitation on galactocentric radius or star formation rate (SFR) surface density (e.g., following den Brok et al. 2021; Yajima et al. 2021).

Leroy et al. (2009), Wilson et al. (2012), and Leroy et al. (2013b) calculated these ratios based on first versions of HERACLES and the JCMT NGLS, and Yajima et al. (2021) recently combined HERACLES and COMING. But the number and quality of CO maps of galaxies have grown significantly compared to any study currently in the literature, particularly with the release of PHANGS-ALMA. Quite a few studies have examined these ratios in individual galaxies and noted local variations in individual ratios (e.g., Crosthwaite & Turner 2007; Koda et al. 2012; Ueda et al. 2012; Vlahakis et al. 2013; Druard et al. 2014; Law et al. 2018; Koda et al. 2020), but so far, there has been relatively little attempt to synthesize these mapping measurements (though see the beam-matched,
single-pointing measurements by Saintonge et al. 2017; Lamperti et al. 2020).

As a practical matter, we present our study in the context of the PHANGS-ALMA CO (2–1) survey. PHANGS-ALMA mapped CO (2–1) across 90 nearby galaxies at ≲150 pc resolution. To aid in the interpretation of these data, we also aim to improve our empirical understanding of the CO (2–1)/CO (1–0) and CO (3–2)/CO (2–1) ratios. Ultimately, we expect this to improve our ability to estimate the molecular mass and infer an appropriate CO-to-H$_2$ conversion factor for these data. This work complements three other recent or forthcoming studies. Den Brok et al. (2021) used nine new, high-quality 22'' resolution CO (1–0) maps from the IRAM 30 m telescope to investigate the resolved CO (2–1)/CO (1–0) ratio. This work also complements the study by T. Saito et al. (in preparation), which investigates the behavior of the CO (2–1)/CO (1–0) ratio at a higher 4''–8'' resolution in four PHANGS-ALMA targets. In scope, our study resembles the recent thorough investigation by Yajima et al. (2021), but we take advantage of a larger database of CO (2–1) maps and include CO (3–2) in our analysis.

After framing some theoretical and observational expectations (Section 2), we describe the data that we use and our measurements (Section 3). Then we measure galaxy-integrated line ratios (Section 4.1) and compare them to the galaxy’s integrated properties (Section 4.2). Then we examine the resolved behavior of the ratio as a function of galactocentric radius, local SFR, and stellar mass surface density (Section 4.3). Finally, we discuss the implications of our measurements and next steps (Section 5) and then summarize our results (Section 6).

## 2. Expectations

Throughout this paper, we refer to the line ratios as

$$ 
\begin{aligned}
R_{21} &\equiv I_{2-1}/I_{1-0}, \\
R_{32} &\equiv I_{3-2}/I_{2-1}, \\
R_{31} &\equiv I_{3-2}/I_{1-0}.
\end{aligned} 
\qquad (1) 
$$

Here $I_{2-1}$, for example, refers to the velocity-integrated specific intensity of the CO (2–1) line, with analogous definitions for the other lines. All intensities, luminosities, and line ratios in this paper are calculated in Rayleigh–Jeans brightness temperature units. Line-integrated intensities are presented in K km s⁻¹, and luminosities are given in K km s⁻¹ pc⁻².

In these kelvin units, we expect a line ratio of 1 for all ratios for an optically thick source in local thermodynamic equilibrium (LTE) when both transitions sit securely on the Rayleigh–Jeans tail given the source temperature, $T_{kin}$. Note, however, that under real conditions, the Rayleigh–Jeans criterion, $h\nu \ll kT_{kin}$, may not be satisfied. This will happen, for example, considering emission in high-frequency transitions from low-temperature sources. When the Rayleigh–Jeans criterion is not met because $T_{kin}$ has low enough values relative to the frequencies of the observed transitions, this will drive the “thermal” value of the line ratio observed from an optically thick source to values below 1. In Appendix B, we illustrate the expected opaque LTE value for the relevant ratios and also show the effects of the cosmic microwave background (CMB). For purposes of reading the observational results in this paper, the key point is that at relevant temperatures, $T_{kin} \sim 10–20$ K,

the expected ratio for opaque, thermalized gas can be as low as ~0.7.

*Expectations from models.* Theoretically, the observed line ratios depend on the distributions of temperature, $T_{kin}$; collider density, $n_{H_2}$; and column density of CO per line width in the gas, $N/\Delta\nu$. A full discussion of the interplay of these quantities with $R_{32}$, $R_{21}$, and $R_{31}$ lies beyond the scope of this work, and we refer the reader to Bolatto et al. (2013b), Shirley (2015), Leroy et al. (2017), and Peñaloza et al. (2017), each of which touches on some aspects of the topic.

As a brief summary, we illustrate the behavior of $R_{32}$ and $R_{21}$ in Figure 1, which plots results from a set of model calculations following Leroy et al. (2017). In the figure, each point shows the line ratios predicted from a model that has a lognormal distribution of collider densities described by a mean density, $n_0$, and a width, $\sigma$. Each model also has a single fixed $T_{kin}$ and a single value of $N_{CO}/\Delta\nu$, which we adopt for each individual density layer. We use RADEX (van der Tak et al. 2007) with data from the Leiden Atomic and Molecular Database (LAMDA; Schöier et al. 2005) to calculate predicted emission. The calculations generally follow Leroy et al. (2017) with the distinction that here we fix $N_{CO}/\Delta\nu$ rather than the optical depth, $\tau$, of a particular line, as in that paper. Because these exact calculations may be of general use and are not fully reported in Leroy et al. (2017), we report the model grid as a machine-readable table in Appendix A.

Figure 1 illustrates the combined effects of temperature, density, and optical depth on the line ratios. In the left panel, each line shows fixed $T_{kin}$ and a fixed density distribution, while we vary $N_{CO}/\Delta\nu$, the total column density of CO molecules normalized to the line width. Here $N_{CO}/\Delta\nu$ affects the optical depth and escape probability and, through these, the critical density and level populations. The figure shows how the low opacities yielded by low $N_{CO}/\Delta\nu$ can lead to high, “superthermal” line ratios with values >1 in the case of low- opacity gas in LTE. Alternatively, for low-density gas, low $N_{CO}/\Delta\nu$ can yield very low line ratios, indicating subcritically excited gas. Meanwhile, higher $N_{CO}/\Delta\nu$ tends to drive gas closer to optically thick LTE and toward line ratios of ~1.

The right panel shows how at fixed $N_{CO}/\Delta\nu$, the density distribution and temperature also play important roles. Their exact impact depends on the $N_{CO}/\Delta\nu$. In general, a higher density or temperature at fixed $N_{CO}/\Delta\nu$ generally drive both ratios toward higher values. The variations for optically thin gas are more extreme, even allowing line ratios above 1, while optically thicker gas shows more dynamic range in $R_{32}$ at the densities illustrated because of the higher excitation requirements of those transitions.

*Expectations from previous observations.* Previous observations established some basic expectations for low-J CO line ratios in nearby galaxies:

1. Normal star-forming galaxies show $R_{21}$ in the range ~0.4–0.9 (e.g., Leroy et al. 2013b; den Brok et al. 2021; Yajima et al. 2021). Here $R_{31}$ likely shows lower values, sometimes as low as ~0.2, in normal galaxies (e.g., Mao et al. 2010; Wilson et al. 2012) but also a larger range of reported values (e.g., Mauersberger et al. 1999; Mao et al. 2010; Lamperti et al. 2020).

2. Starburst and active galaxies show higher, closer-to-thermal (i.e., ~1) ratios (e.g., Mauersberger et al. 1999; Weiß et al. 2005; Mao et al. 2010; Lamperti et al. 2020;
<figure><img src="image_3.png" /><figcaption>**Figure 1.** The CO line ratios produced by model calculations (gray points) with key trends illustrated. Shown are the predicted line ratios $R_{32}$ vs. $R_{21}$ for model calculations using RADEX (van der Tak et al. 2007) and lognormal density distributions following Leroy et al. (2017). Each gray point shows an individual model. As described in Appendix A, each model is characterized by a kinetic temperature ($T_{kin}$), a total CO column density per line width ($N_{CO}/\Delta v$), a lognormal density distribution width ($\sigma$), and a mean H$_2$ collider density ($n_0$). To illustrate how the line ratios change when varying density ($n_0$), temperature ($T_{kin}$), and opacity (set by $N_{CO}/\Delta v$) while holding other parameters fixed, a subset of the data is shown by colored lines. Both panels show the same model grid using gray points while changing the reference. In the left panel, the colored lines show the effect of varying $N_{CO}/\Delta v$, which sets opacity and escape probability, while holding temperature and density constant at a few representative values. In the right panel, dashed lines show the result of varying temperature ($T_{kin}$), and solid lines show the results of varying the mean density ($n_0$) while holding $N_{CO}/\Delta v$ fixed. For all reference lines, we fix $\sigma = 0.6$ dex. Note that these lines are intended to be illustrative, not to span the full phase space of the modeling. The full model grid spanning a plausible range of temperatures and densities is available as a machine-readable table in Appendix A.</figcaption></figure>

Yajima et al. 2021), consistent with having both higher densities and hotter gas.

3. The central parts of normal star-forming galaxies show systematically higher $R_{21}$ (e.g., Braine & Combes 1992; Braine et al. 1993; Leroy et al. 2009, 2013b; Israel 2020; den Brok et al. 2021; Yajima et al. 2021), consistent with higher densities and hotter gas in the central parts of these galaxies (e.g., Mangum et al. 2013; Sun et al. 2020, among many others) and observations showing high temperatures and densities in the center of our own Milky Way (e.g., Ao et al. 2013; Ginsburg et al. 2016; Krieger et al. 2017).

4. In addition to the contrast between normal galaxies and starbursts and between disks and galaxy centers, there is statistical evidence that regions with hotter dust, higher SFR surface density, or shorter depletion times show higher line ratios within normal star-forming galaxies (e.g., Lamperti et al. 2020; den Brok et al. 2021; Yajima et al. 2021).

5. Given the modest dynamic range of the observed ratios and the need to combine multiple instruments, calibration uncertainties can imply significant scatter in line ratio measurements (see excellent discussions in den Brok et al. 2021; Yajima et al. 2021). For single-pointing observations with single-dish telescopes, uncertain aperture corrections also represent a significant source of uncertainty. These systematics, in addition to the limited sensitivity of millimeter-wave telescopes before ALMA, may help explain why many results in the literature show substantial scatter.

These general trends are largely born out by detailed studies of individual galaxies (e.g., Koda et al. 2012, 2020), though there remains disagreement in the literature about the behavior of the ratios, e.g., relative to spiral arms or within individual targets. Some of this may reflect the fact that at high resolution, line ratios can show detailed variations that track the location of individual heating sources or vary across spiral arms and bars (e.g., Ueda et al. 2012; Law et al. 2018; T. Saito et al., in preparation).

### 3. Measurements

Table 1 summarizes the survey combinations and targets for each line ratio. We use new CO (2–1) maps from the PHANGS-ALMA survey (Leroy et al. 2021a), CO (2–1) maps from HERACLES on the IRAM 30 m telescope (Leroy et al. 2009), and another set of IRAM 30 m CO (2–1) maps that cover mostly Virgo Cluster targets (A. Schruba et al., in preparation). We draw CO (1–0) maps from two Nobeyama 45 m surveys, the COMING survey (Sorai et al. 2019) and the Nobeyama CO Atlas of Nearby Spiral Galaxies (hereafter the NRO Atlas; Kuno et al. 2007). We compare these to CO (3–2) maps from the JCMT NGLS (Wilson et al. 2012) and a new APEX LASMA mapping project (J. Puschnig et al., in preparation). In total, as summarized in Table 1, this leads to 152 map pairs with 43 unique galaxies mapped in both CO (2–1) and CO (1–0), 34 mapped in CO (3–2) and CO (2–1), and 20 mapped in CO (3–2) and CO (1–0). A total of 16 unique galaxies have a measurement, not a limit, for all three lines.

To consider a line ratio measurement, we require CO emission to be securely detected in at least one transition, so
<table><caption>**Table 1**
Summary of Map Pairs in This Paper</caption><thead><tr><th>Survey Pair</th><th>Sample Size<br>(Meas./LL/UL)<sup>a</sup></th></tr></thead><tbody><tr><td colspan="2"><i>R</i><sub>21</sub> ≡ CO (2–1)/CO (1–0)</td></tr><tr><td>PHANGS-ALMA+COMING</td><td>10/1/0</td></tr><tr><td>PHANGS-ALMA+NRO Atlas</td><td>18/0/0</td></tr><tr><td>HERA+COMING</td><td>23/1/0</td></tr><tr><td>HERA+NRO Atlas</td><td>23/0/0</td></tr><tr><td>Total map pairs</td><td>76</td></tr><tr><td>Unique galaxies</td><td>43</td></tr><tr><td colspan="2"><i>R</i><sub>32</sub> ≡ CO (3–2)/CO (2–1)</td></tr><tr><td>NGLS+PHANGS-ALMA</td><td>11/0/2</td></tr><tr><td>NGLS+HERA</td><td>22/0/5</td></tr><tr><td>APEX+PHANGS-ALMA</td><td>5/0/0</td></tr><tr><td>APEX+HERA</td><td>2/0/0</td></tr><tr><td>Total map pairs</td><td>47</td></tr><tr><td>Unique galaxies</td><td>34</td></tr><tr><td colspan="2"><i>R</i><sub>31</sub> ≡ CO (3–2)/CO (1–0)</td></tr><tr><td>NGLS+COMING</td><td>13/0/1</td></tr><tr><td>NGLS+NRO Atlas</td><td>11/0/0</td></tr><tr><td>APEX+COMING</td><td>2/0/0</td></tr><tr><td>APEX+NRO Atlas</td><td>2/0/0</td></tr><tr><td>Total map pairs</td><td>29</td></tr><tr><td>Unique galaxies</td><td>20</td></tr></tbody></table>

**Notes.** Surveys: COMING is described by Sorai et al. (2019), NRO Atlas refers to the survey presented by Kuno et al. (2007), PHANGS-ALMA is described by Leroy et al. (2021a), HERA refers to HERACLES (Leroy et al. 2009, 2013b) and a follow-up Virgo Cluster survey (PI: A. Schruba), NGLS refers to the JCMT survey by Wilson et al. (2012) supplemented by a few follow-up or archival JCMT observations, and APEX refers to APEX LASMA mapping by J. Puschnig et al. (in preparation).

that we can at least obtain a limit on the line ratio. For each target that meets this criterion, we estimate the integrated line ratios and compare these to the integrated properties of the galaxy. For targets with a high enough signal-to-noise ratio (S/N), we also measure the line ratio in individual 20" regions, with this 20" scale picked because it represents the common angular resolution achievable by all of the mapping surveys used in our analysis, with APEX LASMA being the limiting data set. The median distance to the target across all of our measurements is ~14 Mpc, where this scale corresponds to ~1.3 kpc, and the 16th–84th percentile range of distance is 9–18 Mpc, implying physical beam sizes of ~0.9–1.8 kpc. This resolution is typically sufficient to resolve the disk of the galaxy but not isolate individual molecular clouds or resolve features like spiral arms or bars. Using these measurements, we correlate the line ratios with local conditions in the galaxy disk, including galactocentric radius, $r_{gal}$; local SFR surface density, $\Sigma_{SFR}$; and stellar mass surface density, $\Sigma_{*}$.

We make separate line ratio measurements for each galaxy and specific survey pair. This can lead to the case where we measure the same line ratio multiple times for a single galaxy; e.g., NGC 0628 appears in both HERACLES and PHANGS-ALMA. We use these duplicated observations to help assess the systematic uncertainty, confirming that key uncertainties in the field still often relate to calibration differences among

telescopes (see below, Figure 2, and more discussion in den Brok et al. 2021). We report all line ratio pairs, but when searching for possible correlations and fitting scaling relations, we adopt only a single value of a line ratio per galaxy using the following priority: PHANGS-ALMA over HERACLES, COMING over the NRO Atlas, and APEX over the JCMT NGLS.

*Conventions.* We correct all quoted surface densities for the effects of inclination. Our stellar mass and SFR maps assume a Chabrier (2003) initial mass function (IMF) and are calibrated to be on the same scale as the GALEX-WISE-SDSS Legacy Survey (Salim et al. 2016, 2018).

### 3.1. CO Data

Table 1 summarizes the sources of our line ratios, which come from combining data from ALMA, the IRAM 30 m telescope, the Nobeyama Radio Observatory (NRO) 45 m telescope, and the JCMT. Specifically, we use the following individual surveys.

**PHANGS-ALMA CO(2–1) data.** Leroy et al. (2021a) described the selection and observations of PHANGS-ALMA, and Leroy et al. (2021b) described the data processing, imaging, and data product creation. Here we use the combined interferometric and total power CO (2–1) cubes convolved to our common resolution of 20". These PHANGS-ALMA CO (2–1) data have a median native resolution of 1"/3, much higher than our working resolution. Because they include total power data, as well as short-spacing 7 m array data, we expect them to have the correct global flux scale and be sensitive to extended emission, and thus to be well suited to this analysis after convolution. Leroy et al. (2021a) confirmed an overall good agreement between the PHANGS-ALMA CO (2–1) and lower-resolution single-dish measurements, which we also show below.

ALMA provides a total power calibration based on regular monitoring of quasars, and the gain uncertainty associated with ALMA at these frequencies is nominally 5%–10%. We take 10% as a conservative estimate, though we note that this likely overestimates the true uncertainty. In Leroy et al. (2021b), we verified that the internal stability of the PHANGS-ALMA total power data appears very good, with fluxes repeatable at the ~3% level from day to day. A few cubes do suffer from 2%–7% gain uncertainties due to issues described in Leroy et al. (2021b). The PHANGS-ALMA data have extremely good sensitivity compared to the other data in this paper, but their field of view tends to be more limited than the other maps, with PHANGS-ALMA typically covering 70% of the total mid-IR (MIR) emission from its target galaxy. As described in Section 3.4, we account for this issue in our analysis.

**IRAM 30 m HERA CO(2–1) data.** We also analyze CO (2–1) maps from HERACLES (Leroy et al. 2009) and another 21 galaxies observed by the IRAM 30 m telescope as part of a survey focused on the Virgo Cluster (PI: Schruba; A. Schruba et al., in preparation). These new data were observed in a manner identical to HERACLES and reduced following the same procedures. Both data sets have a native FWHM beam size of 13"/3 and large extent, usually covering beyond the optical radius of the galaxy. The calibration uncertainty associated with the HERA data is <20% based on a detailed gain analysis presented in den Brok et al. (2021) and a bootstrapping analysis in Leroy et al. (2009) and comparing to PHANGS-ALMA (Leroy et al. 2021a). This calibration

<footnote><sup>a</sup> Entries report number of map pairs yielding a measured line ratio (Meas.), a lower limit (LL), or an upper limit (UL).</footnote>
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Consistency among galaxy-integrated CO fluxes in repeat observations. Each panel shows the ratio among the galaxy-integrated CO luminosity, $L_{CO}$, estimated using different surveys targeting the same galaxy. Each point indicates one galaxy, with the error bar showing only statistical uncertainties. The x-axis shows the CO flux associated with the denominator (i.e., COMING or PHANGS-ALMA). A solid gray line in both panels indicates the expectation for perfect agreement, and the shaded region and colored line show a ±1 standard deviation about the median ratio. The NRO atlas shows 0.91 times the CO luminosity of COMING, with ±12% robustly estimated scatter. Where the two samples overlap in this study, IRAM 30 m HERA mapping shows median 1.03 times higher luminosity than PHANGS with ±12% robustly estimated scatter.</figcaption></figure>

uncertainty can also apply within maps, reflecting uneven gains among the receiver array. Den Brok et al. (2021) noted that the two most extreme cases, NGC 3627 and NGC 5194, were both observed early in the life of the HERA instrument (Schuster et al. 2007), when observing procedures had not yet been optimized. To be conservative, we adopt a nominal uncertainty of 20%, near the upper limit of the plausible calibration uncertainty. This adopted calibration uncertainty agrees with the consistency check shown in Figure 2 and the results of the checks in Leroy et al. (2021a).

*NRO CO(1–0) data.* We utilize CO (1–0) maps obtained by the NRO from the COMING survey (Sorai et al. 2019). These data have a native resolution of 17" and cover large areas in each target. We also compare to CO (1–0) maps from the NRO Atlas presented by Kuno et al. (2007). These have 15" resolution and higher sensitivity than COMING, but the Kuno et al. (2007) data suffer from more visible mapping artifacts and poor baselines compared to Sorai et al. (2019). Following Yajima et al. (2021), we take the gain uncertainty for both data sets to be 25%, reflecting a combination of true calibration uncertainties and pointing errors. Similar to the HERA case noted above, these systematic uncertainties can apply within a galaxy and do not only reflect an overall scaling from galaxy to galaxy.

*JCMT CO(3–2) data.* We also compare to maps of CO (3–2) emission from the JCMT NGLS (Wilson et al. 2012) and individual galaxy follow-up programs (PI: E. Rosolowsky) observed under projects M09BC15, M10BC06, and M12AC03. The follow-up program data were calibrated using the STARLINK software package (Currie et al. 2014) using the observatory-recommended pipelines. The JCMT data initially had 14.′′5 resolution. We convolved them to a resolution of 20′′ for further processing and translated them into the main beam temperature scale using an efficiency of 0.6 at 345 GHz (following Wilson et al. 2012).

After this convolution, we inspected the data and found that they could be improved by fitting and subtracting low-order baselines. For each line of sight, we defined a velocity range of interest based on the velocity range of CO emission seen in the other CO data for the galaxy. Specifically, we defined a reference cube used to define the baseline region, with the priority given to PHANGS-ALMA, then IRAM HERA data, then COMING data, then NRO Atlas data. Because the CO (3–2) data tend to have a lower S/N than the lower-J cubes, using the other cubes as a template to fit the baseline should be well defined and impose little or no bias. Then we considered each line of sight in the JCMT cube. Along that line of sight, we excluded velocities detected above S/N = 3 in either the reference CO cube or the JCMT data themselves from the fit and focused the baseline fit on regions near the detected line emission. We defaulted to a ±50 km s⁻¹ fitting window but adjusted this slightly from galaxy to galaxy. Then we fit and subtracted a baseline from each line of sight using iterative outlier rejection of the JCMT data to refine the fit. We used an order-zero fit; i.e., we subtracted the mean for all galaxies except NGC 2403, where we used a linear fit. In a few galaxies, we also blanked regions of the JCMT cube that were clearly dominated by artifacts. Finally, after inspecting the data, we dropped a few potential targets where the data remained clearly dominated by artifacts despite our baseline fitting. By virtue of using the other CO data as a prior for baseline fits, we also effectively required all JCMT targets to be detected in another CO transition.

We lack a detailed characterization of the calibration and pointing uncertainty for the JCMT data, but according to the JCMT website,<fnref n="25" /> the nominal calibration accuracy of the JCMT is 10% before accounting for uncertainties in pointing or

<footnote n="25">²⁵ https://www.eaobservatory.org/jcmt/instrumentation/heterodyne/calibration/</footnote>
suboptimal conditions. Empirically, the integrated fluxes that we measure frequently vary by ≳10%–20% when we compare results before and after our rebaselining. We adopt an overall calibration uncertainty 20% as a conservative estimate, consistent with, e.g., Yajima et al. (2021) and Sorai et al. (2019) for NRO and den Brok et al. (2021) and Leroy et al. (2009) for IRAM HERA maps. As with the IRAM HERA maps, we consider this to represent the upper envelope of plausible uncertainties.

**APEX CO(3–2) data.** We also compare to five maps of CO (3–2) emission obtained using the LASMA receiver array on the APEX telescope (Güsten et al. 2008).<fnref n="26" /> LASMA is a 7 pixel, single polarization array receiver that can observe from $ν=268$ to 375 GHz. At the $ν≈345$ GHz of CO (3–2), the array has a beam size of ~18"/5. After reduction and convolution during gridding, these maps end up having a beam size of ~20", which sets the common resolution of our data. These maps were obtained as part of APEX projects m-0103.f-9520a-2019 and m-0104.f-9516a-2019 (PI: A. Weiss) and will be presented in detail in J. Puschnig et al. (in preparation). They target galaxies that have PHANGS-ALMA imaging, and the areal extent is designed to match that of the PHANGS-ALMA maps almost exactly. The final data cubes have been gridded to our common velocity resolution of 10 km s⁻¹. At 20" resolution and with 10 km s⁻¹ channels, the cubes have rms noise of ~9–11 mK. The baselines and data quality appear excellent, with few visible artifacts and emission visible in most individual channels. Based on advice from the APEX team, we scale the maps assuming that LASMA has 0.9× the nominal efficiency of APEX at these frequencies in order to account for a partial shadowing of two of the LASMA receivers. In most respects, the data resemble the other array receiver data. We expect the telescope to recover all the flux from the source, and the data should be well suited to stacking. The overall calibration of the data represents the main source of uncertainty for many of our calculations, and we assume LASMA+APEX to have an rms gain uncertainty of ±20%, in line with the other facilities. We defer more details of these data to J. Puschnig et al. (in preparation).

## 3.2. CO Processing

We rebin and downsample all of the CO data to have a velocity resolution of ~10 km s⁻¹. Then we convolve all data cubes to a common angular resolution of 20" and reproject them onto the astrometric grid of the stellar mass maps described below.

For each CO (1–0) and CO (2–1) cube, we produce a three-dimensional high-completeness “signal mask” that includes the volume of the cube where CO emission is detected in either CO (2–1) or CO (1–0), as well as some surrounding volume. As discussed above, the JCMT cubes tend to have a lower S/N and more artifacts than the other data, so at this stage, we rebaseline them as described above. Given this situation, we apply the mask constructed based on the CO (2–1) and CO (1–0) to CO (3–2). Because we have relatively few APEX maps, we treat them the same as the JCMT data. On visual inspection, we do not see any evidence that this approach causes us to miss CO (3–2) emission in our targets. Moreover, as described below, our checks show the masks to have very high completeness for CO (2–1) and CO (1–0); given that the

CO (2–1) data are more sensitive than the CO (3–2) data, we do not expect this choice to bias our results in any important way.

We construct these signal masks following a variation of the masking scheme in Rosolowsky & Leroy (2006) and the “broad masking” approach in Leroy et al. (2021a, 2021b). First, we estimate the noise from signal-free regions of each cube and then calculate the S/N for each pixel in each cube. We then construct individual “signal” masks for each CO (1–0) and CO (2–1) data cube. Each mask began with a high-significance core identified based on a threshold S/N value. We expanded this initial mask to include all adjacent regions of the cube with lower but still significant signal. Then the masks are further dilated by 20 km s⁻¹ in velocity and several beam sizes in each spatial direction. Most of the masking was done at 30" resolution to improve the S/N, but we also included any bright emission seen only at 20" resolution in the mask. We adjusted the exact S/N thresholds used in the masking for each data set based on visual inspection until they yielded a mask that encompassed all visible CO emission in all cubes with a comfortable margin in both velocity and spatial extent. Based on this visual inspection, we also slightly lowered the core S/N threshold in the case of a few compact galaxies with faint CO emission.

For each galaxy, we created a final mask by combining all signal masks from all individual CO (2–1) and CO (1–0) surveys. Any pixel included in any mask is included in the final mask. We adopt this approach aiming at high completeness and minimal bias; i.e., we try to include all likely CO emission in the mask, even if this increases the noise (i.e., this is a “broad” mask following Leroy et al. 2021b). In addition to visual inspection, we verified the completeness of the maps by comparing the integrated CO flux to a direct integral of the cube over the velocity width of the galaxy. For COMING, the masks include a median 97% of the CO emission with <0.03 dex scatter. For PHANGS-ALMA, the masks include 100% of the CO emission, on average, with <0.01 dex scatter.

Finally, we applied this combined mask to all cubes for that galaxy. We collapse this masked cube to construct an integrated intensity (“moment zero”) map for each data cube. We calculate the associated statistical uncertainty from error propagation using the rms noise estimated from the signal-free parts of the cube.

**CO luminosities.** For comparing with the integrated properties of galaxies, we also calculate the CO luminosity, $L_{CO}$, implied by each map. To do this, we adopt the distances compiled in Anand et al. (2021) for PHANGS-ALMA and follow Leroy et al. (2019) for other targets.

In some cases, this calculation is complicated by the fact that the CO line maps do not cover the entire area of the galaxy. In particular, this often affects PHANGS-ALMA (see above and Leroy et al. 2021b). In these cases, we apply an aperture correction that uses WISE3 emission as the template for CO emission. This approach is discussed in more detail in Leroy et al. (2021a), who showed that WISE3 offers the best available template to construct such aperture corrections (consistent with the findings by Chown et al. 2021 that WISE3 correlates very strongly with CO emission).

Note that when we report CO luminosities in Table 2, we give only a single value of $L_{CO}$ (1–0), $L_{CO}$ (2–1), and $L_{CO}$ (3–2) for each galaxy. In choosing which CO luminosity to report, we prefer COMING values over NRO Atlas values, PHANGS-ALMA values over IRAM 30 m HERA values, and APEX

<footnote n="26">²⁶ And see http://www.mpfifr-bonn.mpg.de/5278286/lasma.</footnote>
<table><caption>**Table 2**
Integrated CO Line Ratios and Galaxy Properties</caption><thead><tr><th>Galaxy</th><th>Line Ratio</th><th>Survey Pair</th><th>log<sub>10</sub> R</th><th>Dist. (Mpc)</th><th>log<sub>10</sub> M<sub>☉</sub> (M<sub>☉</sub>)</th><th>log<sub>10</sub> SFR (M<sub>☉</sub> yr<sup>-1</sup>)</th><th>log<sub>10</sub> L<sub>CO</sub><sup>1-0</sup> (K km s<sup>-1</sup> pc<sup>2</sup>)</th><th>log<sub>10</sub> L<sub>CO</sub><sup>2-1</sup></th></tr></thead><tbody><tr><td>ic0750</td><td>R31</td><td>JCMTCOMING</td><td>-0.44 ± 0.04</td><td>17.10</td><td>10.18</td><td>0.20</td><td>8.72</td><td>...</td><td>8.27</td></tr><tr><td>ngc0253</td><td>R21</td><td>PHANGSNROATLAS</td><td>-0.28 ± 0.00</td><td>3.70</td><td>10.64</td><td>0.70</td><td>9.26</td><td>8.96</td><td>...</td></tr><tr><td>ngc0337</td><td>R21</td><td>HERACOMING</td><td>-0.18 ± 0.08</td><td>19.50</td><td>9.80</td><td>0.11</td><td>8.19</td><td>7.98</td><td>...</td></tr><tr><td>ngc0628</td><td>R21</td><td>HERACOMING</td><td>-0.23 ± 0.02</td><td>9.84</td><td>10.34</td><td>0.24</td><td>8.93</td><td>8.66</td><td>8.14</td></tr><tr><td>ngc0628</td><td>R21</td><td>PHANGSCOMING</td><td>-0.30 ± 0.01</td><td>9.84</td><td>10.34</td><td>0.24</td><td>8.93</td><td>8.66</td><td>8.14</td></tr><tr><td>ngc0628</td><td>R31</td><td>JCMTCOMING</td><td>-0.78 ± 0.03</td><td>9.84</td><td>10.34</td><td>0.24</td><td>8.93</td><td>8.66</td><td>8.14</td></tr><tr><td>ngc0628</td><td>R32</td><td>JCMTHERA</td><td>-0.56 ± 0.03</td><td>9.84</td><td>10.34</td><td>0.24</td><td>8.93</td><td>8.66</td><td>8.14</td></tr><tr><td>ngc0628</td><td>R32</td><td>JCMTPHANGS</td><td>-0.43 ± 0.02</td><td>9.84</td><td>10.34</td><td>0.24</td><td>8.93</td><td>8.66</td><td>8.14</td></tr><tr><td>ngc0925</td><td>R32</td><td>JCMTHERA</td><td>-0.06 ± 0.07</td><td>9.16</td><td>9.79</td><td>-0.17</td><td>...</td><td>7.52</td><td>7.52</td></tr><tr><td>ngc1068</td><td>R21</td><td>PHANGSNROATLAS</td><td>-0.15 ± 0.01</td><td>13.97</td><td>10.91</td><td>1.64</td><td>9.47</td><td>9.34</td><td>...</td></tr></tbody></table>

**Note.** This table is a stub. The full version of the table appears as a machine-readable table in the online version of the paper. Columns are as follows. Galaxy: the name of the galaxy. Line Ratio: the reported line ratio. Survey Pair: shorthand for the pair of surveys used to make the measurement. log<sub>10</sub> R: the log<sub>10</sub> of the measured ratio, with uncertainty; in the case of limits, we report the 4σ upper or lower limit as the value. D: the adopted distance in megaparsecs, following Anand et al. (2021). log<sub>10</sub> M<sub>☉</sub>: log of the stellar mass. log<sub>10</sub> SFR: log of the SFR. log<sub>10</sub> L<sub>CO1-0</sub>, L<sub>CO2-1</sub>, and L<sub>CO3-2</sub>: log<sub>10</sub> of the best-estimate CO luminosity in the noted transition with aperture corrections applied. For the luminosity, we report only one best-estimate value per galaxy. That is, this the single best estimate of L<sub>CO</sub>. We give preference to COMING over NRO Atlas and PHANGS-ALMA over IRAM HERA data. Note that because the ratios R are measured over matched apertures inside the galaxies, they do not match the ratios of luminosities by construction.

(This table is available in its entirety in machine-readable form.)

values over JCMT values. Because we provide only a single luminosity, and because the luminosities include aperture corrections while the reported ratios use exactly matched apertures, we note that dividing our quoted CO luminosities will not yield exactly the same value as the line ratio reported in the table. That is, our reported R<sub>21</sub> uses a matched field of view and is calculated for each survey pair (see Section 3.4), while the CO luminosities are aperture-corrected, and we report only one value for each transition.

### 3.3. Stellar Masses and SFRs

We estimate SFRs and stellar masses based on Galaxy Evolution Explorer (GALEX; Martin et al. 2005) far-ultraviolet (FUV) and near-ultraviolet (NUV) images, Spitzer IRAC near-infrared (NIR) maps (Fazio et al. 2004), and Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010) NIR and MIR imaging. The GALEX and WISE maps were created as part of the z = 0 Multiwavelength Galaxy Synthesis (Leroy et al. 2019). The IRAC maps were obtained mostly by the S<sup>4</sup>G survey (Sheth et al. 2010). Leroy et al. (2019, 2021a) give details of the conversion from these bands to SFR and M<sub>☉</sub>. We use the same calculations described in Leroy et al. (2021a), which we carried out for PHANGS-ALMA, the targets of the HERA surveys, and the targets of the Nobeyama surveys in a self-consistent way. Figure 3 shows the resulting estimated SFR and M<sub>☉</sub> for our targets plotted over a large set of local galaxies (from Leroy et al. 2019).

Briefly, to estimate the SFR, we use the best available combination of ultraviolet and MIR data, preferring more stable combinations of tracers whenever available. In order of most preferred to least preferred, we use FUV+WISE4, NUV+WISE4, FUV+WISE3, NUV+WISE3, WISE4 only, and WISE3 only. We adopt the conversions reported in Table 7 of Leroy et al. (2019) and apply them as detailed in Section 3 of that paper. These conversions are calibrated to reproduce galaxy-integrated SFR values calculated for the Sloan Digital Sky Survey (SDSS) main galaxy sample based on UV-to-IR CIGALE spectral energy distribution modeling by Salim et al. (2016, 2018). As discussed in that paper, these agree well with

previous calibrations using similar bands (e.g., Salim et al. 2007; Kennicutt & Evans 2012; Leroy et al. 2012; Janowiecki et al. 2017), usually within ~0.1 dex. In Leroy et al. (2021a), we showed that the resolved Σ<sub>SFR</sub> estimates agree with high-quality Balmer decrement-based Σ<sub>SFR</sub> measurements from PHANGS-MUSE (E. Emsellem et al. submitted) within ~20%, on average, but that the UV+IR maps likely overestimate Σ<sub>SFR</sub> in regions of low SFR, with the most likely explanation being contamination by IR cirrus (see Groves et al. 2012; Leroy et al. 2012), but other effects including stochastic sampling of the IMF or issues with extinction correction also remain possible. The magnitude of the effect may reach up to a factor of 2 for Σ<sub>SFR</sub> ≲ 10<sup>-3</sup> M<sub>☉</sub> yr<sup>-1</sup> kpc<sup>-2</sup>.

We base our stellar mass estimates on NIR emission at 3.6 µm (IRAC1) or 3.4 µm (WISE1). After subtracting a background, we flag stars and replace them with interpolated values from similar galactocentric radii. Then we convert from NIR intensity to stellar mass surface density using a mass-to-light ratio that depends on the ratio of SFR to WISE1. This quantity serves as a proxy for the specific SFR, SFR/M<sub>☉</sub>, which is a strong predictor of the WISE1 mass-to-light ratio in Salim et al. (2016, 2018). Leroy et al. (2021a) described the detailed calculations and presented comparisons to results from resolved stellar mass estimates from optical spectral mapping by PHANGS-MUSE (E. Emsellem et al. 1999 A&A submitted). Leroy et al. (2019) presented the motivation for the approach based on matching the Salim et al. (2016, 2018) estimates.

We measure integrated M<sub>☉</sub> and integrated SFR by directly integrating all pixels within 2r<sub>25</sub>. Based on comparisons among different methods and bands, we adopt uncertainties of 0.1 dex for both M<sub>☉</sub> and SFR estimates. When relevant, we calculate offsets from the star-forming main sequence exactly as described by Leroy et al. (2021a).

### 3.4. Line Ratio Measurements

Before proceeding, we reproject all data, which have already been convolved to 20", onto a grid with a pixel size equal to 20"; i.e., we work with pixels equal to the FWHM beam size.
This leads to a moderate undersampling of the maps in exchange for rendering the individual measurements mostly independent. We consider that the convolution to 20" has removed most sampling effects present in the on-the-fly single-dish maps (e.g., see Mangum et al. 2007).

### 3.4.1. Integrated Line Ratios

We calculate each integrated line ratio over the area where both surveys involved have coverage and where the combined mask described in Section 3.2 indicates the presence of CO emission. For a given galaxy, we denote this matched area $A$, and we calculate the line ratio, $R_{ul}$, as the ratio of the sum of emission from each line:

$$R_{ul} = \frac{\sum_{x,y \in A} I_u(x, y)}{\sum_{x,y \in A} I_l(x, y)} \qquad (2)$$

This ratio-of-sums approach weights the calculated $R_{ul}$ by intensity, and when the maps and mask cover the whole galaxy, $R_{ul}$ will match the result expected from an unresolved, single-pointing measurement.

We follow standard error propagation to estimate the statistical uncertainty on the measurement. The uncertainty in each measurement is the sum in quadrature of this statistical uncertainty with the calibration uncertainties for both telescopes: $\sigma^2 = \sigma_{stat}^2 + \sigma_{cal,u}^2 + \sigma_{cal,l}^2$. The calibration uncertainty frequently dominates the total error budget.

We require both lines to be detected at a statistical S/N > 4 to report a ratio, i.e., S/N > 4 before accounting for the calibration uncertainties. For cases where only one line is detected at the required significance, we estimate and report an upper or lower limit using the $4\sigma$ statistical uncertainty in the undetected line to define the limit.

**Literature data.** We compare our galaxy-integrated measurements to recent measurements combining IRAM 30 m CO (1–0) maps from the EMPIRE survey with PHANGS-ALMA, HERACLES, and a new M51 CO (2–1) map (den Brok et al. 2021). In that case, we use the same procedure to calculate $M_*$ and SFR described above.

We also compare to the 28" single-dish APEX+IRAM 30 m line ratio measurements presented by Saintonge et al. (2017). These have closely matched beams, and their stellar masses and SFR values are calculated on a system similar to our own. They do not report their exact aperture correction for the IRAM 30 m data but note it to be between 2% and 10%. We apply a 5% upward correction to all IRAM 30 m luminosities and include a 15% overall calibration uncertainty in addition to their reported statistical error.

**Effect of the CMB.** The observed brightness temperature reflects only the contrast against the CMB, such that the measured brightness temperature, $T_b$, will be $T_b = (1 - e^{-\tau})(T_{ex} - T_{CMB})$ for each transition (e.g., see Eckart et al. 1990; Bolatto et al. 2013a; Zschaechner et al. 2018, among many other discussions), with $T_{ex}$ the relevant excitation temperature. This can imply modest corrections to the line ratios, especially for cold clouds. However, this radiative transfer proceeds only at the scale of the molecular clouds themselves. The brightness temperatures in our current work are heavily affected by beam dilution. Without measuring the clumping of CO emission at subresolution, we cannot calculate an appropriate correction for the CMB. These values can, in principle, be measured for PHANGS-ALMA (e.g., following Leroy et al. 2013a), but we lack similar high-resolution

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Galaxies with beam-matched line ratio measurements in SFR-$M_*$ space. Shown is the location of our targets in SFR-$M_*$ space, plotted over the larger sample of galaxies within $d < 50$ Mpc, with SFR and $M_*$ calculated in a consistent way (Leroy et al. 2019). The black line shows an estimate of the SFR-$M_*$ relationship for local star-forming galaxies (Leroy et al. 2019, in excellent agreement with Catinella et al. 2018). Red circles, blue hexagons, and green squares show galaxies with measurements (not limits) in this paper, with the top panel showing galaxies with $R_{21}$ measured and the bottom panel showing $R_{32}$ and $R_{31}$. For comparison, we also plot the properties of EMPIRE targets from den Brok et al. (2021) and the $R_{21}$ properties of the subset of xCOLD GASS with beam-matched $R_{21}$ measurements from Saintonge et al. (2017). Despite the relatively large number of measurements, our $R_{21}$ measurements remain highly clustered in this space, biased toward high-mass and high-SFR galaxies. The measurements involving the CO (3–2) line show slightly more homogeneous sampling of parameter space despite their overall lower number and higher typical uncertainty.</figcaption></figure>

templates for the other data, and the measurements for PHANGS-ALMA represent future work. We note the effect, do not apply any CMB correction, and leave an improved treatment for future work. See Appendix B for more details.

**Consistency among integrated measurements.** Because surveys targeting the same line sometimes share targets, we make repeated measurements for several galaxy ratio pairs. Figure 2 checks for internal consistency within our measurements. The left panel shows the ratio of CO (1–0) luminosity estimated from the NRO Atlas to that from COMING, and the right panel shows the ratio of CO (2–1) luminosity estimated using HERA to that from PHANGS. In both panels, we use the integrated galaxy luminosity and thus trust the aperture corrections described above to account for any differences in area covered. We do not show a panel for CO (3–2). Only one
galaxy is detected in both JCMT and APEX, NGC 3627, and there the luminosity inferred from the APEX data is ~1.4 times that calculated from the JCMT.

Overall, Figure 2 illustrates that the CO (2–1) measurements are mostly consistent across the two surveys, with a median ratio only a few percent different from 1.0. We do observe significant scatter, with an rms variation of about ±12%, much larger than the statistical noise. This agrees with Leroy et al. (2021a) and mostly validates the calibration uncertainties adopted above. The situation for CO (1–0) is similar, with the NRO Atlas ~8% lower than COMING, on average, and a scatter of about ±12% from a relatively low sample size. Based on Sorai et al. (2019), we expect that COMING has better overall calibration compared to the NRO Atlas.

Overall, Figure 2 shows that systematic uncertainties related to calibration, pointing, etc. impose an uncertainty that has an rms of order 10%–20% on individual CO line measurements from galaxies. We will see in the rest of the paper that this uncertainty is comparable to the range of variation in the line ratios across the galaxy population.

### 3.4.2. Resolved, Binned, Normalized Line Ratios

In Section 4.3, we compare line ratios to location within a galaxy. The challenges here are the limited S/N of individual measurements and the need to account for the substantial galaxy-to-galaxy calibration uncertainties.

We focus on three quantities: galactocentric radius, $r_{gal}$; $\Sigma_{SFR}$, the SFR per unit area; and $\Sigma_{SFR}/\Sigma_*$, the local specific SFR. We consider the area covered by both surveys and extract measurements of both relevant CO lines, $r_{gal}$, $\Sigma_{SFR}$, and $\Sigma_{SFR}/\Sigma_*$ for all 20" pixels in this overlap region. We calculate $r_{gal}$ using the orientations and distances in Leroy et al. (2021a), drawing on Lang et al. (2020) and Anand et al. (2021). For cases outside PHANGS-ALMA, we prefer orientation parameters from $S^4G$ (Sheth et al. 2010; Muñoz-Mateos et al. 2015), where available, and follow the compilation in Leroy et al. (2019) otherwise. We calculate bins for both physical $r_{gal}$, expressed in units of kiloparsecs, and $r_{gal}$ normalized to the effective half-mass radius, $r_{eff}$, calculated in Leroy et al. (2021a). We also calculate $\Sigma_{SFR}$ and $\Sigma_{SFR}/\Sigma_*$ as described above and in Leroy et al. (2021a).

To account for the limited S/N of individual measurements, we define a set of bins in each quantity of interest. Then, within each galaxy, we identify the pixels belonging to each bin and then sum all data for each line. As with the global line ratios, we divide the summed, binned values by one another to estimate the line ratio in that bin. As above, we propagate statistical uncertainties following standard error propagation, and we use an S/N threshold of 4 to determine whether a bin is a detection (both numerator and denominator have S/N > 4), an upper limit (only denominator has S/N > 4), or a lower limit (only numerator has S/N > 4). After calculating the line ratio, we account for uncertainty associated with the calibration by normalizing each measured line ratio by the galaxy-average value. That is, in Section 4.3, we measure the enhancement or depression of the ratio relative to its mean value in any given galaxy. This should remove any global gain calibration uncertainty term, though not local calibration uncertainties, e.g., due to pointing uncertainties or pixel-to-pixel gain variations. This also removes any real galaxy-to-galaxy scatter in the mean line ratio, so that this analysis focuses on how these

variables drive relative changes in each line ratio within a galaxy.

We note the following details regarding bin construction.

1. When considering physical galactocentric radius, in units of kiloparsecs, we use linearly spaced bins 1 kpc in width with the first bin centered at $r_{gal} = 0$ kpc and the last one centered at 10 kpc. Note that, as discussed above, this implies some slight over- or undersampling of the data because the range of distances to the galaxies means that our 20" resolution corresponds to different physical resolution across the sample.

2. When considering normalized galactocentric radius, we normalize by the half-mass radius, $r_{eff}$, calculated following Leroy et al. (2021a). Our bins have a width of 0.5 times $r_{eff}$ with the first bin centered at $r_{gal} = 0$ and the outermost bin centered at $r_{gal} = 3r_{eff}$.

3. When considering $\Sigma_{SFR}$, we bin the data by $\log_{10} \Sigma_{SFR} / \langle \Sigma_{SFR} \rangle$. Here $\langle \Sigma_{SFR} \rangle$ is the galaxy-averaged SFR surface density. We calculate via $\langle \Sigma_{SFR} \rangle = 0.5 \text{ SFR}/(\pi r_{eff}^2)$, i.e., the surface density implied by placing half of the galaxy-integrated star formation within the effective radius, $r_{eff}$, measured for the mass by Leroy et al. (2021a). Normalizing in this way allows us to focus on how the internal structure of the line ratio tracks the local SFR with fewer concerns about how the overall amplitude of the SFR or the calibration of our SFR tracer varies from galaxy to galaxy. This makes sense given our similar galaxy-by-galaxy normalization of the CO line ratio for this analysis.

4. For specific SFR, we calculate $\Sigma_{SFR}/\Sigma_*$, normalize by the integrated galaxy-averaged $\text{sSFR} = \text{SFR}/M_*$, and then bin $\log_{10}(\Sigma_{SFR}/\Sigma_*)/\text{sSFR}$ in bins of 0.25 dex from -1.25 to +1.25 dex about the galaxy average.

This binning procedure is functionally equivalent to a stacking procedure within each bin similar to that used by, e.g., Cormier et al. (2018), Jiménez-Donaire et al. (2019), or den Brok et al. (2021). It has the advantage of retaining information from individual pixels with a modest S/N and so avoids some biases present in direct pixel-by-pixel analysis. We record bins in which both lines are detected at S/N > 4 as measurements and record upper and lower limits using the $4\sigma$ value for the limiting line. Typically, for $R_{21}$, any limits are lower limits because the CO (2–1) maps are more sensitive than the CO (1–0) maps. For $R_{32}$ and $R_{31}$, our limits are mostly upper limits because the CO (3–2) maps lack sensitivity compared to the CO (2–1) and CO (1–0) maps.

# 4. Results

## 4.1. Global Line Ratios

In total, as reported in Table 1, we study 152 pairs of overlapping surveys. For each measured $R_{21}$, $R_{32}$, and $R_{31}$, Table 2 gives the name, survey pair, adopted distance to the galaxy, estimated SFR, $M_*$, and CO luminosity, $L_{CO}$, in each line. Note that, as discussed above, we only quote one CO luminosity per line for each galaxy. This is the $L_{CO}$ that we suggest using as characteristic of the galaxy, not the value used in the calculation of the line ratio, because our quoted $L_{CO}$ includes an aperture correction. When a galaxy was covered in the same line by multiple surveys, we chose which to use for
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Galaxy-integrated CO line ratios. Colored data points show galaxy-integrated line ratios measured from maps, including limits. The error bars on the individual points reflect both statistical and calibration errors, as described in Section 3. For each survey pair or literature result, the black circle shows the median, and the error bar shows the ±1σ reported range for that survey pair or literature study. The top row shows results for $R_{21}$, the middle row shows results for $R_{32}$, and the bottom row shows results for $R_{31}$. In each panel, the thick gray horizontal line and shaded gray region show the median and 16th–84th percentile range from all mapping data (see Table 3) in this study.</figcaption></figure>

$L_{CO}$ following the same prioritization of surveys as noted in Section 3.4.

Figure 4 shows the distribution of the three measured line ratios, along with ranges reported for nearby galaxies in the literature. The first part of each plot shows results for all mapping surveys, then we separate the results according to the survey pair used for the measurement. Table 3 reports basic results for the distributions combining all mapping surveys.

$R_{21} \equiv \text{CO}(2-1)/\text{CO}(1-0)$. Treating all data equally, we find a median $R_{21} = 0.61$ and a mean $R_{21} = 0.65$ with a 16th–84th percentile range from 0.50 to 0.83. This reflects 76 mapping measurements from 43 galaxies in this work and nine measurements from den Brok et al. (2021). Note that this distribution allows repeat measurements when the same galaxy was targeted by multiple surveys and will therefore weigh “popular” targets more heavily. Since our goal here is to
<table><caption>**Table 3**
Galaxy-integrated CO Line Ratios from Mapping Surveys</caption><thead><tr><th>Ratio</th><th>Mean</th><th>Median</th><th>16th Percentile</th><th>84th Percentile</th></tr></thead><tbody><tr><td>R<sub>21</sub><sup>a</sup></td><td>0.65</td><td>0.61</td><td>0.50</td><td>0.83</td></tr><tr><td>R<sub>32</sub></td><td>0.50</td><td>0.46</td><td>0.23</td><td>0.59</td></tr><tr><td>R<sub>31</sub></td><td>0.31</td><td>0.29</td><td>0.20</td><td>0.42</td></tr></tbody></table>

Notes. See Figure 4.

a Includes EMPIRE measurements from den Brok et al. (2021).

synthesize the current literature, we consider this approach reasonable, and we use only a single best value for each galaxy when fitting the scaling relations below. As Figure 3 shows, the targets of the surveys we consider emphasize high-mass, high-SFR galaxies (see also Section 5.3).

As Figure 4 shows, our $R_{21}$ values agree well with previous results for normal star-forming galaxies. We find an almost identical median value to the 0.64 measure for EMPIRE galaxies in den Brok et al. (2021), though our data show higher scatter than theirs. This likely reflects both the high quality of the EMPIRE CO (1–0) maps presented in den Brok et al. (2021) and the narrower range of galaxy properties sampled by EMPIRE, which we illustrate in Figure 3. We also find almost perfect agreement with Yajima et al. (2021), who also derived a median of $R_{21} = 0.61$ with a scatter of $\pm 0.19$. The sample in Yajima et al. (2021) represents a subset of our own, so we expect this close match. Finally, our measurements also agree reasonably well with previous HERACLES results by Leroy et al. (2013b), who found a median $R_{21}$ of 0.67$^{27}$ with a scatter of $\pm 0.16$ dex or $\pm 44\%$, corresponding to a range of $\sim 0.46–0.97$. Finally, we also agree well with the 0.58 median and $\sim 0.49–0.77$ $1\sigma$ range for literature single-pointing measurements compiled by den Brok et al. (2021).

Our measurements appear slightly lower than the xCOLD GASS measurements by Saintonge et al. (2017). We attribute this mostly to selection effects, though given that the line ratios drop with radius (Section 4.3), there could also be some mild impact from the limited xCOLD GASS beam size. As Figure 3 shows, the xCOLD GASS $R_{21}$ measurements target a wider range of stellar mass than our current $R_{21}$ sample. The lower-mass galaxies and high-SFR/$M_\star$ galaxies in the xCOLD GASS sample likely shift the median to the higher average value of $R_{21} = 0.8$ that they report. For reference, if we include the Saintonge et al. (2017) measurements in our sample, the combined data set has a median $R_{21} = 0.64$, mean $R_{21} = 0.74$, and 16th–84th percentile range of 0.52–0.91.

This paper does not focus on individual targets, but we briefly note that the three high $R_{21} > 1.0$ values seen in Figure 4 are all consistent with $R_{21} \lesssim 1$ within the uncertainties. These are NGC 1087 in PHANGS-ALMA+COMING, NGC 2976 in HERA+NRO Atlas, and NGC 4536 in HERA+NRO Atlas. Given the sample size and magnitude of the uncertainties, we expect a few such outliers.

Summarizing, Figure 4 shows overall good convergence among recent studies of $R_{21}$. Adopting a typical value of $R_{21} = 0.65$ with a 30% uncertainty that reflects scatter across the galaxy population represents a good assumption for high-mass, $\log_{10} M_\star [M_\odot] \gtrsim 10.25$ (see Figure 3), $z=0$ galaxies on

the main sequence of star-forming galaxies (Leroy et al. 2013b; den Brok et al. 2021; Yajima et al. 2021).

$R_{32} \equiv \text{CO}(3-2)/\text{CO}(2-1)$. We find a median $R_{32} = 0.46$ with a mean $R_{32} = 0.50$ and a 16th–84th percentile range of 0.23–0.59. Because CO (3–2) is comparatively faint and the CO (2–1) data have higher S/N, upper limits affect our $R_{32}$ distribution more than the other two ratios. We have treated the upper limits as equal to our minimum measured ratio for quantifying the distribution. This choice mainly affects our 16th percentile estimate. With seven of 40 $R_{32}$ measurements being upper limits, the 16th percentile quoted is set to the lowest measured ratio.

Compared to $R_{21}$ and $R_{31}$, $R_{32}$ has the least extensive sample of previous beam-matched or mapping-based studies of whole nearby galaxies; we are only aware of the work by Wilson et al. (2012), who found a mean $R_{32}$ of 0.36 with $\pm 0.13$ scatter using earlier versions of the same data that we use here. This is moderately lower than our calculated mean value. We attribute part of the difference to revisions to the adopted IRAM main beam efficiency (see above and Leroy et al. 2013a) after Wilson et al. (2012) made their measurements and to our ability to match the areas used for the calculations in this paper, which was not possible in Wilson et al. (2012). Taking these factors into account, the measurements appear roughly consistent. We also note that in Figure 4, the APEX and JCMT NGLS data show hints of an offset. As far as we can tell, this reflects a mixture of small number statistics and perhaps the choice to focus the initial APEX LASMA mapping on bright, actively star-forming targets, which may have more excited molecular gas. Only one target overlaps between the two surveys, NGC 3627, and for that case, we do find an $\sim 30\%$ higher CO (3–2) luminosity from APEX than JCMT, but this is only a single target.

It does appear that $R_{32}$ has a larger dynamic range than $R_{21}$. Because lower limits confuse the 16th percentile estimate for $R_{32}$, we compare the interquartile (25th–75th percentile) ranges for the two ratios and find an $\sim 0.14$ dex range for $R_{21}$ and $\sim 0.27$ dex for $R_{32}$. Though caveats related to a small sample size and the effect of lower limits still apply, this agrees with the expectation (see Section 2) that $R_{32}$ shows significant excitation variations across the range of real conditions found in molecular gas and indicates that the ratio has the potential to act as a strong diagnostic of the local excitation of molecular gas.

$R_{31} \equiv \text{CO}(3-2)/\text{CO}(1-0)$. We find a median $R_{31} = 0.29$ and a mean $R_{31} = 0.30$ with a 16th–84th percentile range of 0.20–0.41. Though the samples used to calculate the ratios vary, our $R_{31}$, $R_{32}$, and $R_{21}$ values are approximately “closed,” as expected, with $R_{31} \sim R_{32} \times R_{21}$, implying reasonable self-consistency within our measurements.

The literature reports a wide range of values. Our measurement is high compared to the $0.18 \pm 0.06$ reported by Wilson et al. (2012) comparing NRO Atlas and JCMT NGLS data. Note, however, that Wilson et al. (2012) cautioned that they did not match areas for the comparison, so their lower value could simply reflect a mismatch in apertures. Our measurements do agree well with the results for a smaller set of galaxies from Wilson et al. (2009). They used the NGLS and NRO Atlas to study regions in individual galaxies and found ratios in the $\sim 0.3–0.4$ range.

Our measurements lie on the low end of the range $R_{31} = 0.2–1.0$ found by Mauersberger et al. (1999), but note

<footnote>$^{27}$ Note that the earlier 0.8 recommended by Leroy et al. (2009) based on HERACLES was revised down by Leroy et al. (2013b) based on updated estimates of the IRAM 30 m main beam efficiency.</footnote>
that their sample also includes many starburst and active galaxies and only ~seven normal spiral galaxies. A similar case holds for Mao et al. (2010), who obtained R<sub>31</sub> = 0.7 ± 0.5. Partially based on those studies, we expect much higher, approaching thermal, R<sub>31</sub> in active and dense galaxies, so the contrast with our more quiescent targets seems reasonable. This appears to be mostly born out by our the results in Section 4.2. Our values also appear low compared to the mean R<sub>31</sub> = 0.55 and ±0.23 scatter found for 28 star-forming galaxies by Lamperti et al. (2020) using single-pointing measurements. The R<sub>31</sub> ratio clearly shows a wide dynamic range and so good promise as a diagnostic. However, the state of R<sub>31</sub> observations remains fairly limited, though not quite so much as for R<sub>32</sub>.

With R<sub>31</sub> < R<sub>32</sub> < R<sub>21</sub> < 1, on average, these values appear consistent with the standard picture that most low-J CO emission from nearby star-forming galaxies comes from optically thick clouds with the CO (2–1) and CO (3–2) transitions moderately subthermally excited (e.g., see Section 2, Weiß et al. 2005; Bolatto et al. 2013a). We discuss this more in Section 5.

## 4.2. Comparison to Integrated Galaxy Properties

As discussed in Section 2, these ratios are both expected and observed to vary between galaxies. We test for correlations between all three line ratios and global galaxy properties in Figures 5 and 6. Before doing so, we again highlight the relatively narrow range of galaxy properties covered by current mapping surveys, visualized in Figure 3 and discussed in Section 5.3. Especially for R<sub>21</sub>, the current measurements focus on high-mass galaxies, mostly in the range log<sub>10</sub>M<sub>⋆</sub>[M<sub>☉</sub>] ~ 10.25–10.75, with SFRs just above the star-forming main sequence. Though the data are sparser for the ratios involving CO (3–2), these measurements do span a larger range of stellar mass and SFR at fixed stellar mass.

Figure 5 shows the rank correlation coefficients relating each line ratio to M<sub>⋆</sub>, SFR, specific SFR (SFR/M<sub>⋆</sub>), offset from the main sequence of star-forming galaxies (ΔMS), CO luminosity, and SFR per unit L<sub>CO</sub> (SFR/L<sub>CO</sub>). The colored bars show the expected correlation and 1σ scatter for the null hypothesis. In the cases of M<sub>⋆</sub>, SFR, SFR/M<sub>⋆</sub>, and ΔMS, we adopt the simple null hypothesis that the line ratio and the other variable are not correlated. Then, the expected correlation will be zero, and the 1σ scatter reflects the range from randomly re-pairing the variables. Recall that here, unlike in Figure 4, we use only a single best estimate of each line ratio for each galaxy.

In the case of L<sub>CO</sub> and SFR/L<sub>CO</sub>, the two variables used in the correlation will be correlated by construction. In this case, we construct the null hypothesis as follows. First, we measure the logarithmic scatter in the real line ratio across our data set. Then, as the null hypothesis, we assume a fixed underlying line ratio and that this scatter is evenly distributed across the variables involved. We generate the expectations for the null hypothesis by randomly applying this scatter repeatedly to each variable and calculating the rank correlation coefficient. Here we consider the individual CO luminosities, not the line ratios, as the underlying variables for this exercise. This ensures that the null hypothesis captures the built-in correlation between the axes. For example, when we correlate the line ratio with L<sub>CO,upp</sub>, i.e., the luminosity in the numerator, we measure the scatter in the log<sub>10</sub> of the line ratio, σ, and then we apply √2 times σ in the model noise to each of L<sub>CO,upp</sub> and L<sub>CO,low</sub>. Then we construct the model line ratios R = L<sub>CO,upp</sub>/L<sub>CO,low</sub> before

<figure><img src="image_8.png" /></figure>

**Figure 5.** Rank correlation of line ratios with galaxy properties. Shown are the Spearman rank correlation coefficients relating global galaxy properties to measured line ratios, neglecting the limits. The bars show the 16%–84% range of coefficients expected under the null hypothesis of no real correlation (note that in the cases with correlated axes, the rank correlation coefficient for no true correlation may not be centered on zero; see the text). Here L<sub>CO,low</sub> refers to the luminosity of the lower CO transition in the ratio, e.g., CO (1–0) in R<sub>21</sub>, while L<sub>CO,upp</sub> refers to the upper transition, e.g., CO (3–2) in R<sub>31</sub>. As the contrast between results for L<sub>CO,upp</sub> and L<sub>CO,low</sub> demonstrates, the correlated axes play an important but hard-to-avoid role in determining correlations with L<sub>CO</sub> and SFR/L<sub>CO</sub>. Most correlations show marginal significance but are consistent with low-mass, high-SFR/M<sub>⋆</sub>, and high-SFR/L<sub>CO</sub> galaxies showing higher line ratios, especially among the measurements involving CO (3–2). See also Figure 6.

calculating the expected correlation. That is, we define the null hypothesis to be the case where the line ratio is fixed and the variance matches the observed variance and is randomly distributed among the relevant variables.

Overall, Figures 5 and 6 show a consistent sense of variations. Lower-mass galaxies, which also have lower L<sub>CO</sub>, higher SFR/M<sub>⋆</sub>, and higher SFR/L<sub>CO</sub>, tend to show higher line ratios. The absent or weak correlation with SFR can be understood as competing effects; low-M<sub>⋆</sub> galaxies have higher SFR/M<sub>⋆</sub> but also a lower overall SFR. The radiation field in low-mass galaxies may be more intense due to a higher local Σ<sub>SFR</sub>, but the integrated SFR will still be lower. In general, correlations with integrated galaxy properties appear stronger for R<sub>32</sub> and R<sub>31</sub> compared to R<sub>21</sub>. This partially reflects the broader range of galaxy properties covered by those
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Correlations between line ratios and galaxy-integrated properties. Shown are the galaxy-integrated $R_{21}$, $R_{32}$, and $R_{31}$ as a function of quantities expected to relate to excitation, specifically SFR/$M_*$ and SFR/$L_{CO,low}$. Each point represents one map pair for one galaxy, with the black points showing limits. Error bars reflect both statistical and calibration uncertainties. The numbers in the top left corners report the Spearman $\rho$ relating the quantities on the x- and y-axes. The gray lines show highly approximate scaling relations that go through the data and then saturate at the thermal $log_{10}R = 0$; see Section 4.2.1, Equation (3), and Table 4. As discussed in Section 4.2, the axes in the right column are correlated, and this may exert a significant effect on the observed relation.</figcaption></figure>

measurements (Figure 3) and may also reflect that $R_{32}$ and $R_{31}$ have more sensitivity than $R_{21}$ to the range of conditions found in normal galaxies (see Section 2).

These trends make physical sense and agree with the limited previous measurements. Physically, elevated SFR/$M_*$ may trace more intense radiation fields and stronger heating of the
gas, suggesting higher temperatures. The anticorrelation with $M_*$ may reflect the impact of dust shielding. Based on the existence of the mass–metallicity relation (e.g., Tremonti et al. 2004; Kewley & Ellison 2008), we expect the low-mass members of our sample to also have lower dust-to-gas ratios (e.g., Leroy et al. 2011; Rémy-Ruyer et al. 2014; Casasola et al. 2020) and more intense radiation fields. In literature studies, CO line ratios do appear enhanced in low-metallicity regions or galaxies (e.g., Lequey et al. 1994; Bolatto et al. 2003; Druard et al. 2014; Kepley et al. 2016; Cicone et al. 2017, among many others). Higher SFR/$L_{CO}$ may indicate poorly shielded, low-metallicity gas in which the CO persists only in the core of a molecular cloud (e.g., see discussion in Glover & Clark 2012; Schruba et al. 2012; Bolatto et al. 2013b; Rubio et al. 2015). Alternatively, higher SFR/$L_{CO}$ can indicate more efficiently star-forming gas, which will often be denser gas with more nearby heating sources. These are both factors that can lead to higher line ratios, especially $R_{32}$ and $R_{31}$ (see Section 2). Given that our sample skews toward relatively massive and thus nearly solar metallicity targets, we expect that these density and heating effects likely represent the main drivers of the observed correlations. The correlations that we see agree with the results of Lamperti et al. (2020), who showed a correlation between SFR/$L_{CO,low}$ and $R_{31}$, and Yajima et al. (2021), who used a subset of the data we consider here and showed a correlation between $R_{21}$ and SFR/$L_{CO,low}$. Qualitatively, Figure 5 echoes other results at low and high redshift that show a correlation between normalized star formation activity and excitation (e.g., Weiß et al. 2005; Bolatto et al. 2013b; Liu et al. 2021).

Although the pattern in Figure 5 makes physical sense, the trends are not particularly significant. The p-values relating $M_*$ to $R_{31}$ and $R_{32}$ are only 0.24–0.29. For SFR/$M_*$, the p-values for $R_{31}$ and $R_{32}$ are ~0.06–0.13, more significant but still indicating only weak correlations. The other significant correlations involve $L_{CO,low}$, the luminosity of the lower CO transition in the line ratio. We report these because the results make physical sense and are interesting, but the line ratio (~$L_{CO,upp}/L_{CO,low}$) and a quantity involving $L_{CO,low}$ are correlated by construction. To see this, contrast the significant correlations seen in Figure 5 for all line ratios and $L_{CO,low}$ and the lack of similar significant correlations for $L_{CO,upp}$. For the moment, we only caution that these results include the effects of correlated measurements and should be taken as indicative but likely overstate the significance of the correlation in the data.

With these caveats in mind, Figures 6 and 7 visualize the correlations between each line ratio and global quantities: SFR/$M_*$, SFR/$L_{CO,low}$, $M_*$, and $L_{CO,low}$. The data show large scatter but also evidence for overall correlation with the sense that higher line ratios emerge from galaxies with high SFR/$M_*$ and/or SFR/$L_{CO}$. The correlations of line ratios with $M_*$ and $L_{CO,low}$ appear weaker.

### 4.2.1. Approximate Scaling Relations

Given the current state of the data, especially the large systematic uncertainties and uneven sampling of the SFR–$M_*$ space, rigorous scaling relations to predict the line ratios are probably not feasible. Nonetheless, we find it useful to illustrate current best fits and sketch a plausible approach to line ratio scaling relations.

We adopt the following principles.

1. The line ratio should not go above 1. Although there are a few measurements with $R > 1$ in the literature, such ratios likely reflect optically thin gas, which is not expected or observed to be the dominant emitting component over large parts of galaxies. One should cap any fit at $R \sim 1$ (see Section 2 and Appendix B).

2. The line ratio also should not drop to arbitrarily low values, even when a galaxy has a low SFR, specific SFR, or stellar mass. The appropriate lower limit for the ratios is less immediately clear than the upper limit, but based on Section 2, we suggest that $R_{21} \gtrsim 0.5$, $R_{32} \gtrsim 0.2$, and $R_{31} \gtrsim 0.1$ represent reasonable lower bounds.

3. Lacking any strong physical motivation for another functional form, we assume a power law intermediate between these two limits.

Thus, we suggest that the following broken power law represents a reasonable way to predict $y = \log_{10} R$ from some variable $x$ that shows a positive correlation with $R$:

$$y = \begin{cases} y_{\text{low}}, & \text{if } x < x_{\text{low}} \\ y_{\text{low}} + m(x - x_{\text{low}}), & \text{if } x_{\text{low}} < x < x_{\text{high}}, \\ 0.0, & \text{if } x > x_{\text{high}} \end{cases} \quad (3)$$

where

$$y \equiv \log_{10} R \quad \text{and} \quad m = \frac{0.0 - y_{\text{low}}}{x_{\text{high}} - x_{\text{low}}}.$$

Each relation has three parameters: $y_{low}$, the $\log_{10}$ value of the line ratio at low values of the x-axis; $x_{low}$, the x-axis value where the ratio begins to increase; and $x_{high}$, the x-axis value where the line ratio reaches the thermal value $y \equiv \log_{10} R = 0.0$, i.e., $R = 1$. We adopt the values of $y_{low}$ described above based on physical expectations, though we note that, particularly for $R_{21}$, a lower value could be plausible. Then, for both $x = \log_{10} \text{SFR}/M_*$ and $x = \log_{10} \text{SFR}/L_{CO,low}$, we fit $x_{low}$ and $x_{high}$ by minimizing $\chi^2$ while varying $x_{low}$ and $x_{high}$. We report the best-fit values in Table 4 and illustrate them in Figure 6, but we emphasize again that in our current data, the measured correlations have only marginal significance, especially for $R_{21}$. Improvements to both the data quality and the sample studied are needed before we will have a strong predictive relation to estimate $R$ based on other galaxy properties.

### 4.3. Correlations with Local Conditions

As described in Section 3.4.2, we also separate galaxies into individual regions to measure how changes in local conditions inside a galaxy relate to variations in the line ratios. Figures 8 and 9 and Tables 3 and 6 summarize the results of these measurements.

We conduct this analysis by binning the data within each galaxy according to several properties of interest: galactocentric radius, SFR surface density ($\Sigma_{SFR}$), and local specific SFR ($\Sigma_{SFR}/\Sigma_*$). The binning allows us to improve the S/N of individual measurements. Despite the averaging, many bins still lack a detection, and Figures 8 and 9 also show upper and lower limits. Here a detection has S/N > 4 in both lines, while a limit has S/N > 4 in only one line. We do not analyze bins with neither line detected.
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Correlations between line ratios and galaxy-integrated properties. Shown are the galaxy-integrated $R_{21}$, $R_{32}$, and $R_{31}$ as a function of $M_*$ and $L_{CO,low}$. Each point represents one map pair for one galaxy, with the black points showing limits. Error bars reflect both statistical and calibration uncertainties. The numbers in the top left corners report the Spearman $\rho$ relating the quantities on the x- and y-axes. As discussed in Section 4.2, the axes in the right column are correlated, and this may exert a significant effect on the observed relation. The correlations here appear even weaker than those in Figure 6, and we do not show any possible functional relations.</figcaption></figure>

We focus this analysis on the relative behavior of the line
ratio and these local quantities. In order to do this, we
normalize all line ratio measurements to their average value for

the galaxy. We perform a similar normalization of $\Sigma_{SFR}$ and $\Sigma_{SFR}/\Sigma_*$, normalizing the measurements for each individual galaxy by the galaxy-averaged values, $\langle\Sigma_{SFR}\rangle$ and $\langle\text{SFR}/M_\ast\rangle$.
<table><caption>**Table 4**
Approximate Scaling Relations</caption><thead><tr><th>Ratio</th><th>Quantity</th><th>y<sub>low</sub></th><th>x<sub>low</sub></th><th>x<sub>high</sub></th><th>Scatter<sup>a</sup></th></tr></thead><tbody><tr><td>R<sub>21</sub></td><td>SFR/M<sub>*</sub></td><td>-0.3</td><td>-10.78</td><td>-9.02</td><td>0.07</td></tr><tr><td>R<sub>21</sub></td><td>SFR/L<sub>CO,low</sub><sup>b</sup></td><td>-0.3</td><td>-9.12</td><td>-7.52</td><td>0.07</td></tr><tr><td>R<sub>32</sub></td><td>SFR/M<sub>*</sub></td><td>-0.7</td><td>-11.25</td><td>-9.36</td><td>0.10</td></tr><tr><td>R<sub>32</sub></td><td>SFR/L<sub>CO,low</sub><sup>b</sup></td><td>-0.7</td><td>-9.50</td><td>-7.54</td><td>0.08</td></tr><tr><td>R<sub>31</sub></td><td>SFR/M<sub>*</sub></td><td>-1.0</td><td>-11.25</td><td>-9.02</td><td>0.08</td></tr><tr><td>R<sub>31</sub></td><td>SFR/L<sub>CO,low</sub><sup>b</sup></td><td>-1.0</td><td>-9.50</td><td>-7.60</td><td>0.08</td></tr></tbody></table>

**Notes.** Coefficients for the indicative scaling relations following Equation (3) and illustrated in Figure 6. The fits are derived from assuming $y_{low}$ and then fitting $x_{low}$ and $x_{high}$ by minimizing $\chi^2$. We report all quantities, including the scatter, in dex. As shown in the figure, these should be taken as approximate, and we expect future work to revise them considerably.

<sup>a</sup> “Scatter” reports the median absolute value of the residuals about the fit in dex.

<sup>b</sup> Here $L_{CO,low}$ refers to the CO luminosity of the lower transition in the line ratio, e.g., CO (1–0) in $R_{21}$.

We present the trends with galactocentric radius in both physical units of kiloparsecs and radius normalized by the half-mass radius, $r_{eff}$. This normalization means that our analysis mostly controls for calibration uncertainties and galaxy-to-galaxy scatter. More details, including the exact bin definitions, can be found in Section 3.4.2.

**CO excitation and galactocentric radius.** All three line ratios show significant anticorrelation with radius and enhanced values relative to the galaxy mean in the inner parts of galaxies. This has been seen before for samples of galaxies by Leroy et al. (2009), Wilson et al. (2009), Leroy et al. (2013b), den Brok et al. (2021), and Yajima et al. (2021), as well as for many individual galaxies, though note that only Wilson et al. (2009) used CO (3–2) among these studies.

In Figure 8, we see that radial gradients in all three line ratios appear to be a general feature. This is true whether we express the galactocentric radius in physical units (left panels) or normalize to the effective radius (right panels). To characterize the gradient, we fit a linear function of the form

$$ \log_{10} \frac{R}{\langle R \rangle} = a \times Q + b, \quad (4) $$

where $R/\langle R \rangle$ is the line ratio relative to the galaxy mean and $Q$ is $r_{gal}$, either expressed in kiloparsecs or normalized to $r_{eff}$, the half-light radius. We conduct a simple $\chi^2$ minimization fit to all of the median values (colored points in Figure 8) in each bin. We calculate these medians accounting for lower limits in $R_{21}$ and upper limits in $R_{32}$ and $R_{31}$ and use bins that have at least five nonlimit measurements and for which the 16%, 50%, and 84% values are not limits. We report $a$ and $b$ for the fits for each line in Table 5.

We find gradients, given by the $a$ values in Table 5, of -0.011, -0.009, and -0.028 dex kpc$^{-1}$ or -0.043, -0.090, and -0.162 dex $r_{eff}$. This implies, on average, a mild but significant change in each of the line ratios across the disk. The gradients appear weaker for $R_{21}$ than for $R_{32}$ or $R_{31}$, but they are statistically significant in all cases based on the $p$-value associated with the Spearman rank correlation coefficient reported in the table.

**Central enhancements.** For $R_{21}$ and $R_{31}$ in Figure 8, much of the apparent gradient is driven by high $R$ in the innermost parts

of galaxies. In Table 6, we quantify the central enhancement in each line ratio relative to the disk-averaged value. We report the median and 16th–84th percentile values in the innermost 1 kpc wide bin. We measure median enhancements of 0.18 dex for $R_{21}$ and 0.27 dex for $R_{31}$ and a weaker median enhancement of 0.08 dex for $R_{32}$.

Recall that the samples are not matched among the three lines, but the qualitative point seems clear: both $R_{21}$ and $R_{31}$ appear enhanced in the inner parts of normal star-forming galaxies, and $R_{32}$ appears at least mildly enhanced. In addition to agreeing with previous measurements, this trend makes physical sense given that both gas density and star formation activity tend to rise toward the inner parts of galaxies (e.g., Usero et al. 2015; Gallagher et al. 2018a; Jiménez-Donaire et al. 2019, for studies showing this via HCN observations). Recall that any galaxy without central gas detected at $S/N > 4$ in both lines will not appear in this analysis.

Are the centers special or simply the end points of a steady radial gradient? To test this, Table 5 also includes fits that exclude the central, $r_{gal} = 0$, bin. Other than excluding these $r_{gal} = 0$ data, the fits are identical. These “no center” fits still all show negative gradients. In other words, all of the line ratios decline with galactocentric radius outside of the galaxy center, and the radial trends do not appear to be purely products of central enhancements. However, for $R_{21}$ and $R_{31}$, the radial gradients excluding the center appear much weaker than those including the center. In fact, the $R_{21}$ trends appear almost flat outside the inner two data points. This demonstrates that much of the overall decline in the ratio with radius is driven by the behavior of the inner galaxy, rather than a smooth gradient. This also agrees with the visual impression from Figure 8. In short, bright galaxy centers do appear special, showing evidence of enhanced excitation, and there also appear to be weak radial gradients outside just the central region.

**CO excitation and star formation activity.** Figure 9 shows two closely related trends. We plot each line ratio as a function of normalized $\Sigma_{SFR}$ and $\Sigma_{SFR}/\Sigma_{*}$. Both variables trace the local star formation activity, which should relate directly to the heating of the gas and indirectly to the gas density. Again, we report rank correlation coefficients and $\chi^2$ minimization fits to the median values in Table 9 (now $Q = \log_{10} \Sigma_{SFR}/\langle\Sigma_{SFR}\rangle$ or $\log_{10} \Sigma_{SFR}/\Sigma_{*}/\langle\Sigma_{SFR}/\Sigma_{*}\rangle$).

The left column of Figure 9 shows more or less continuous trends relating enhancement in each line ratio to enhancement in $\Sigma_{SFR}$. The figure does give some hint of “flattening” in $R_{21}$ and $R_{31}$, at least at low $\Sigma_{SFR}/\langle\Sigma_{SFR}\rangle$. This might relate to the structure seen in the radial profiles, in which the central enhancements appear strong while the ratio, at least for $R_{21}$, appears flatter outside the center. More high-S/N observations of the extended, low-$\Sigma_{SFR}$ parts of galaxy disks are needed to pursue this further.

Table 5 reports the gradient in each line ratio per dex change in $\Sigma_{SFR}/\langle\Sigma_{SFR}\rangle$. For $R_{21}$, this is an ~0.13 dex change in $R_{21}/\langle R_{21} \rangle$ per dex change in $\Sigma_{SFR}/\langle\Sigma_{SFR}\rangle$. This agrees well with the slope relating $R_{21}$ to total infrared surface brightness, $\Sigma_{TIR}$, measured by den Brok et al. (2021) for nine EMPIRE targets.

The line ratios all appear at least moderately enhanced in regions with higher star formation activity, and again, these trends make physical sense. More intense activity correlates with higher-density gas (e.g., Usero et al. 2015; Gallagher et al. 2018a; Jiménez-Donaire et al. 2019), stronger radiation fields,
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Correlations between line ratios and galactocentric radius. Each plot shows line ratio enhancement or deficit relative to the galaxy average, $log_{10} R/(\langle R \rangle)$, as a function of galactocentric radius. Individual points showing binned results for individual galaxy plus survey pair combinations are calculated as described in Section 3. Arrows show $4\sigma$ upper or lower limits. Colored points indicate the median in each radial bin considering all galaxies and survey pairs with a measurement in that bin, including limits. The black vertical error bars show the 16th–84th percentile of measurement in that bin. We only plot median values, where a nonlimit 16%, 50%, and 84% value could be inferred from our measurements. The solid colored lines show a least-squares fit to all bins with at least five data points (see Table 5). The left panels show bins constructed considering physical radius, while the right panels normalize the radius by the half-mass radius, $r_{\text{eff}}$.</figcaption></figure>

and higher cosmic-ray densities, all leading to warmer gas. The results here show that after accounting for galaxy-to-galaxy variations, internal correlations between star formation activity and galactocentric radius are evident and widespread inside nearby, normal, star-forming galaxies.

The right column of Figure 9 shows how changes in $R$ relate to variations in the specific SFR, $\Sigma_{\text{SFR}}/\Sigma_{\star}$. In principle, this quantity has appealing properties as a predictor of line ratio variations; its behavior should be more scale-independent than

$\Sigma_{\text{SFR}}$, and the normalization by $\Sigma_{\star}$ removes some overall scaling effects to isolate the intensity of star formation.

Figure 9 and Table 5 show that we do observe significant correlations with the expected sense between $R_{21}$ and $\Sigma_{\text{SFR}}/\Sigma_{\star}$ and $R_{31}$ and $\Sigma_{\text{SFR}}/\Sigma_{\star}$, but $R_{32}$ does not show a formally significant correlation. The correlations between $R_{21}$ and $R_{32}$ and $\Sigma_{\text{SFR}}/\Sigma_{\star}$ are also weaker than those with $\Sigma_{\text{SFR}}$ (see Table 5), which is reflected in the noisier appearance of the correlations in Figure 9, especially at extreme values.
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Correlations between line ratios and star formation activity. Each plot shows line ratio enhancement or deficit relative to the galaxy average, $log_{10} R/(R_*)$, as a function of local SFR surface density and specific SFR. Before binning, we normalize $\Sigma_{SFR}$ and $\Sigma_{SFR}/\Sigma_*$ by their galaxy averages, so that the plots show how enhancements in the line ratio relative to the galaxy average correlate with enhancements in SFR or SFR/$M_*$ relative to the galaxy average. As in Figure 8, gray points show results for individual map pairs calculated as described in Section 3, and arrows show $4\sigma$ upper or lower limits. Colored points indicate the median result for all map pairs with a measurement in that bin, including limits, and so distill the general trend for the whole data set. The black vertical error bars show the 16th–84th percentile of individual measurements. The solid colored lines show a least-squares fit to all median values with at least five data points (see Table 5).</figcaption></figure>

Contrasting the two columns suggests that $\Sigma_{SFR}$ may be a somewhat better predictor of line ratio variations within normal galaxies than $\Sigma_{SFR}/\Sigma_*$. This could make physical sense if gas density and the overall gas reservoir play more important roles than the strength of the interstellar radiation field and the intensity of star formation. Alternatively, it may simply reflect that $\Sigma_{SFR}$ shows a larger dynamic range and more regular structure within our target galaxies than $\Sigma_{SFR}/\Sigma_*$.

**Overall.** Our resolved analysis offers a consistent first-order picture. The inner parts of galaxies often host the most intense star formation activity, and these regions appear enhanced in all line ratios in all plots. These enhancements also correlate with increased $\Sigma_{SFR}$ and, perhaps to a lesser extent, higher $\Sigma_{SFR}/\Sigma_*$. The trends become weaker outside the inner parts of galaxies. Tables 5 and 6 give quantitative estimates of gradients and central enhancements.
<table><caption>**Table 5**
Correlations between Line Ratio Variations and Local Conditions</caption><thead><tr><th>Ratio</th><th>Quantity</th><th>a</th><th>b</th><th>ρ<sup>a</sup></th></tr></thead><tbody><tr><td>R<sub>21</sub>/⟨R<sub>21</sub>⟩</td><td>r<sub>gal</sub> [kpc]</td><td>-0.011</td><td>0.038</td><td>-0.37</td></tr><tr><td>R<sub>21</sub>/⟨R<sub>21</sub>⟩</td><td>r<sub>gal</sub>—no center<sup>b</sup></td><td>-0.003</td><td>-0.003</td><td>...</td></tr><tr><td>R<sub>32</sub>/⟨R<sub>32</sub>⟩</td><td>r<sub>gal</sub> [kpc]</td><td>-0.009</td><td>0.051</td><td>-0.35</td></tr><tr><td>R<sub>32</sub>/⟨R<sub>32</sub>⟩</td><td>r<sub>gal</sub>—no center<sup>b</sup></td><td>-0.007</td><td>0.043</td><td>...</td></tr><tr><td>R<sub>31</sub>/⟨R<sub>31</sub>⟩</td><td>r<sub>gal</sub> [kpc]</td><td>-0.028</td><td>0.125</td><td>-0.48</td></tr><tr><td>R<sub>31</sub>/⟨R<sub>31</sub>⟩</td><td>r<sub>gal</sub>—no center<sup>b</sup></td><td>-0.019</td><td>0.077</td><td>...</td></tr><tr><td>R<sub>21</sub>/⟨R<sub>21</sub>⟩</td><td>r<sub>gal</sub>/r<sub>eff</sub></td><td>-0.043</td><td>0.045</td><td>-0.48</td></tr><tr><td>R<sub>21</sub>/⟨R<sub>21</sub>⟩</td><td>r<sub>gal</sub>/r<sub>eff</sub>—no center<sup>b</sup></td><td>-0.021</td><td>0.006</td><td>...</td></tr><tr><td>R<sub>32</sub>/⟨R<sub>32</sub>⟩</td><td>r<sub>gal</sub>/r<sub>eff</sub></td><td>-0.090</td><td>0.097</td><td>-0.51</td></tr><tr><td>R<sub>32</sub>/⟨R<sub>32</sub>⟩</td><td>r<sub>gal</sub>/r<sub>eff</sub>—no center<sup>b</sup></td><td>-0.099</td><td>0.111</td><td>...</td></tr><tr><td>R<sub>31</sub>/⟨R<sub>31</sub>⟩</td><td>r<sub>gal</sub>/r<sub>eff</sub></td><td>-0.162</td><td>0.177</td><td>-0.60</td></tr><tr><td>R<sub>31</sub>/⟨R<sub>31</sub>⟩</td><td>r<sub>gal</sub>/r<sub>eff</sub>—no center<sup>b</sup></td><td>-0.131</td><td>0.135</td><td>...</td></tr><tr><td>R<sub>21</sub>/⟨R<sub>21</sub>⟩</td><td>Σ<sub>SFR</sub>/⟨Σ<sub>SFR</sub>⟩</td><td>0.129</td><td>0.019</td><td>0.55</td></tr><tr><td>R<sub>32</sub>/⟨R<sub>32</sub>⟩</td><td>Σ<sub>SFR</sub>/⟨Σ<sub>SFR</sub>⟩</td><td>0.075</td><td>0.023</td><td>0.38</td></tr><tr><td>R<sub>31</sub>/⟨R<sub>31</sub>⟩</td><td>Σ<sub>SFR</sub>/⟨Σ<sub>SFR</sub>⟩</td><td>0.261</td><td>0.057</td><td>0.63</td></tr><tr><td>R<sub>21</sub>/⟨R<sub>21</sub>⟩</td><td>Σ<sub>SFR</sub>/Σ<sub>*</sub>/⟨Σ<sub>SFR</sub>/Σ<sub>*</sub>⟩</td><td>0.112</td><td>0.025</td><td>0.43</td></tr><tr><td>R<sub>32</sub>/⟨R<sub>32</sub>⟩</td><td>Σ<sub>SFR</sub>/Σ<sub>*</sub>/⟨Σ<sub>SFR</sub>/Σ<sub>*</sub>⟩</td><td>0.031</td><td>0.024</td><td>0.27</td></tr><tr><td>R<sub>31</sub>/⟨R<sub>31</sub>⟩</td><td>Σ<sub>SFR</sub>/Σ<sub>*</sub>/⟨Σ<sub>SFR</sub>/Σ<sub>*</sub>⟩</td><td>0.131</td><td>0.077</td><td>0.35</td></tr></tbody></table>

**Notes.** Correlations between each line ratio, normalized to the galaxy value, and various local conditions. These are calculated by binning the data within individual galaxies to increase the S/N as described in Section 3. Columns *a* and *b* report the slope and intercept for a linear fit of the form in Equation (4) carried out on the median binned data. Column *ρ* reports the Spearman rank correlation coefficient relating enhancements in the line ratio to the local quantity considering the individual bins (gray points in Figures 8 and 9).

<sup>a</sup> The Spearman rank coefficient only for bins with both lines detected at S/N > 4.

<sup>b</sup> These fits exclude the innermost radial bin.

<table><caption>**Table 6**
Central Enhancements in Line Ratios</caption><thead><tr><th>Ratio</th><th>Central Enhancement [dex]</th></tr></thead><tbody><tr><td>R<sub>21</sub>/⟨R<sub>21</sub>⟩</td><td>0.18 (0.04–0.27)</td></tr><tr><td>R<sub>32</sub>/⟨R<sub>32</sub>⟩</td><td>0.08 (–0.01–0.19)</td></tr><tr><td>R<sub>31</sub>/⟨R<sub>31</sub>⟩</td><td>0.27 (0.14–0.42)</td></tr></tbody></table>

**Note.** Enhancement in line ratio at galaxy center relative to the galaxy-averaged value. The table reports the median value of log<sub>10</sub>R/⟨R⟩ calculated in a 1 kpc bin centered at r<sub>gal</sub> = 0 kpc (see Figure 8). The quoted range gives the 16th–84th percentile range of measurements in that bin.

A note on correlations with Σ<sub>mol</sub> and SFR/CO. Note that we deliberately avoid stacking by SFR/CO because at the modest S/N present in the individual pixels, the effect of the correlated axes becomes overwhelming (see den Brok et al. 2021, for more discussion). We did explore this direction. As expected based on previous results, the resolved trends relating R/⟨R⟩ to Σ<sub>SFR</sub>/Σ<sub>mol</sub> or Σ<sub>mol</sub> appear highly significant (see also Yajima et al. 2021). However, as we also saw in the previous section, interpreting these trends is difficult due to the correlated axes, and in these cases, we lack an independent, high-S/N quantity to stack the data (i.e., we would have to use the CO itself as the x-axis) or stack by a third quantity (e.g., radius) and carefully handle upper limits in the stacks. For now, we note only that these correlations appear significant but driven by correlated axes and that this will be a crucial trend to carefully analyze in

data with higher S/Ns, excellent interline calibration, and well-understood uncertainties.

# 5. Discussion

## 5.1. Comparison to Simple Models Assuming Lognormal Density Distributions

What do these measured ratios imply for physical conditions in the cold gas? Figure 10 compares our measurements to the models described in Section 2 and Figure 1. In the top left panel, we replot the models from Figure 1 but using a logarithmic grid. In that same panel, we also indicate the 16th–84th percentile range of all measured R<sub>21</sub> and R<sub>32</sub> (i.e., the values from Figure 4 and Table 3) as shaded colored regions, and we plot points for all 16 galaxies that have measured CO (1–0), CO (2–1), and CO (3–2) luminosities. Note that here we plot only the ratios implied by the luminosities in Table 2, because this allows us to ensure a matched area among all measurements.

The other three panels of Figure 10 visualize the set of conditions in the model grids that produce both R<sub>21</sub> and R<sub>32</sub> in the measured 16th–84th percentile range. That is, these are the model grid entries that lie inside the square overlap region in the top left panel and appear as darker shaded points in that panel. The most common conditions in the model grid appear comparable to conditions expected in the cold ISM. We find log<sub>10</sub>(N<sub>CO</sub>/Δν) [cm<sup>-2</sup> (km s<sup>-1</sup>)<sup>-1</sup>] ≈ 16.5–17. For a standard CO abundance N(CO)/N(H<sub>2</sub>) ~ 10<sup>-4</sup> (e.g., van Dishoeck & Black 1988) and a typical molecular cloud line full width of ~10 km s<sup>-1</sup> (e.g., Sun et al. 2018, 2020; Rosolowsky et al. 2021), this N<sub>CO</sub>/Δν implies a total log<sub>10</sub>N(H<sub>2</sub>) [cm<sup>-2</sup>] ≈ 21.5–22, in reasonable agreement with column densities for molecular clouds or resolved surface densities in maps of cold gas in galaxies.

The models that produce the observed line ratios generally have low temperatures (see bottom right panel in Figure 10), preferring T<sub>kin</sub> < 20 K and, frequently, T<sub>kin</sub> = 10–15 K, the lowest temperatures that we modeled.

Meanwhile, the bottom left panel in Figure 10 shows the mean, n<sub>0</sub>, and width, σ, of the lognormal distribution of densities used in the models. Each individual model combines an ensemble of single-density models in a way that mimics a lognormal distribution of densities. Because a wider distribution with higher σ includes higher-density gas, it can produce higher excitation and line ratios (e.g., see Leroy et al. 2017). As a result, the mean density, n<sub>0</sub>, and the width, σ, somewhat trade off.<sup>28</sup> The overall sense of the panel is that the viable models tend to include gas with densities in the range of n<sub>H2</sub> ~ 300–1000 cm<sup>-3</sup>. We expect that the addition of a power-law tail of densities (e.g., as expected for self-gravitating gas; see Krumholz & Thompson 2007; Federrath & Klessen 2013; Burkhart 2018) could affect the CO (3–2) emission but would likely yield mostly similar mean densities.

In short, our observed line ratios can be broadly produced by cold, intermediate-density, intermediate column density gas. Though our models include density distributions and so exhibit a range of optical depth, τ, we find it useful to note the implied optical depth for a single zone model with the same average properties that we consider. A one-zone RADEX model with T<sub>kin</sub> = 15 K, n<sub>H2</sub> = 1000 cm<sup>-3</sup>, and N<sub>CO</sub>/Δν = 3 × 10<sup>16</sup> cm<sup>-2</sup> (km s<sup>-1</sup>)<sup>-1</sup> yields τ ~ 4 for CO (1–0), τ ~ 9 for

<footnote n="28"><sup>28</sup> The mean ln n depends on both n<sub>0</sub> and σ<sup>2</sup>, following Padoan & Nordlund (2002) for a lognormal distribution of densities ln n = ln n<sub>0</sub> + σ<sup>2</sup>/2.</footnote>
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Comparison of measurements and model calculations. The top left panel shows the same model grid visualized in Figure 1 but with logarithmic axes. Individual points show $L_{CO}$-based estimates for $R_{21}$ and $R_{32}$ for galaxies with all three lines measured. Note that we use the $L_{CO}$ values and now show only one point per galaxy to homogenize the area considered and use the best available data. Shaded regions show the 16th–84th percentile range for each ratio from Table 3. Darker points highlight entries in the model grid that lie within the measured 16th–84th percentile range for both ratios. In the remainder of the panels, we visualize the properties of the model grid entries that lie in the overlap of the average $R_{32}$ and $R_{21}$ constraints. The top right panel shows the distribution of viable models in $N_{CO}/\Delta v - n_0$ space, the bottom left panel shows $\sigma - n_0$ space, and the bottom right panel shows $n_0 - T_{kin}$ space. The line ratios can be most readily produced by cold, moderate-opacity (inferred from the column per line width), moderate-density gas (from $n_0$ and $\sigma$).</figcaption></figure>

CO (2–1), and $\tau \sim 5$ for CO (3–2) (suggesting heavy excitation to the $J=2$ state). For CO (1–0), these are also in good agreement with the opacities implied by contrasting $^{13}$CO (1–0) and $^{12}$CO (1–0) measurements (e.g., see summary in Roman-Duval et al. 2010; Cao et al. 2017; Cormier et al. 2018).

Though these conditions appear reasonable, we note that not all of our measurements can be readily explained by the models that we consider. This appears entirely reasonable, given that the simple models we use are most appropriate for individual molecular clouds or parts of galaxies. In the top left panel of

Figure 10, we see observations with higher $R_{32}$ at lower $R_{21}$ than the model readily produces, i.e., up and to the left of the model points, though note that most of the data still lie within $\sim 1\sigma$ of the model grid.

More sophisticated models may be able to explain the observed line ratios by mixing different models considering, e.g., multimodal density distributions or distributions of $T_{kin}$ or column density. This would naturally reflect the blending of clouds and regions in different physical states that we expect to be averaged together in our whole-galaxy measurements.
Alternatively, abundance variations, the impact of cosmic rays, geometry, and coupling between zones can all broaden the range of parameter space covered by the models (e.g., Bemis & Wilson 2019; Bisbas et al. 2019). Comparisons to numerical models of galaxy disks that include chemistry and radiation transfer (e.g., Gong et al. 2020; Tress et al. 2020) will help illuminate line ratios that we might expect from a realistic blend of conditions.

## 5.2. Line Ratios, the CO-to-H₂ Conversion Factor, and Scaling Relations

As discussed by both den Brok et al. (2021) and Yajima et al. (2021), an immediate implication of a systematic dependence of the line ratios on local conditions is that the slopes of the scaling relations, e.g., measured between $\Sigma_{SFR}$ and $\Sigma_{mol}$, will differ depending on the observed transition. This is in addition to offsets in normalization that reflect $R_{21}$, $R_{32}$, and $R_{31}$.

At the simplest level, our mean line ratios represent factors that can be used to place relations measured using different transitions on the same scale. That is, our median $R_{21}$ and $R_{31}$ from Table 3 can and should be used to renormalize relations derived using CO (2–1) or CO (3–2) onto a consistent scale with those derived using CO (1–0). Equivalently, our line ratios can be combined with the standard Milky Way $\alpha_{CO}^{1-0} = 4.35 M_{\odot}pc^{-2}(K \text{ km s}^{-1})^{-1}$ to yield

$$ 
\begin{aligned}
\alpha_{\text{CO}}^{2-1} &= 6.7 M_{\odot}pc^{-2}(K \text{ km s}^{-1})^{-1}, \\
\alpha_{\text{CO}}^{3-2} &= 14.0 M_{\odot}pc^{-2}(K \text{ km s}^{-1})^{-1}.
\end{aligned} 
\quad (5) $$

Equation (5) helps set the normalization of any relation involving $\Sigma_{mol}$ or $M_{mol}$. Beyond this, Table 5 and Figure 9 suggest that, when considering the slope of any power-law-style scaling relation involving $\Sigma_{SFR}$, one should expect differences of ~0.1 between CO (3–2) and CO (2–1), ~0.15 between CO (2–1) and CO (1–0), and ~0.25–0.3 between CO (3–2) and CO (1–0). These values come from the slope of the gradient fit between each line ratio and $\Sigma_{SFR}$. The most aggressive application of these trends would be to use them to adjust the power-law index of scaling relations measured using different lines to place them on a common scale. In this case, we expect that

$$ 
\begin{aligned}
\text{For } \Sigma_{\text{SFR}} &\propto I_{1-0}^a \\
\text{If } R_{21} \propto \Sigma_{\text{SFR}}^{\beta} \text{ then } \Sigma_{\text{SFR}} &= I_{2-1}^{a/(1+a\beta)} \\
\text{If } R_{32} \propto \Sigma_{\text{SFR}}^{\gamma} \text{ then } \Sigma_{\text{SFR}} &= I_{3-2}^{a/(1+a\gamma)}
\end{aligned} 
\quad (6) $$

so that, e.g., for a linear $\Sigma_{SFR} \propto I_{1-0}^{1.0}$ then if we adopt the coefficients in Table 5, one expects $\Sigma_{SFR} \propto I_{2-1}^{0.88}$ or $\Sigma_{SFR} \propto I_{3-2}^{0.79}$. Alternatively, these can be inverted to yield the index associated with the fundamental CO (1–0) transition via

$$ 
\begin{aligned}
\text{For } \Sigma_{\text{SFR}} &\propto I_{2-1}^b \text{ and } R_{21} \propto \Sigma_{\text{SFR}}^{\beta} \\
\text{then } \Sigma_{\text{SFR}} &= I_{1-0}^{b/(1-b\beta)} \\
\text{And for } \Sigma_{\text{SFR}} &\propto I_{3-2}^c \text{ and } R_{32} \propto \Sigma_{\text{SFR}}^{\gamma} \\
\text{then } \Sigma_{\text{SFR}} &= I_{1-0}^{c/(1-c\gamma)}.
\end{aligned} 
\quad (7) $$

Two issues complicate this approach. First, recall that, as discussed in Sections 2 and 4.2.1 and Appendix B, the line ratios likely saturate at or near 1. They also seem unlikely to

drop to arbitrarily low values. As a result, the scaling of the line ratios with $\Sigma_{SFR}$ likely occurs only over a bounded range. Because the scaling likely only occurs over a limited range of $\Sigma_{SFR}$, simply adjusting the index of a power-law fit by these values will likely overcorrect when the scaling relation is measured over a high dynamic range. For example, $R_{21}$ variations seem unlikely to strongly affect the slope of scaling relations measured for highly excited major mergers. A more conservative approach would be to only apply such a correction to more quiescent, “normal” galaxies and add a one-sided component to the systematic uncertainty estimate to reflect the uncertainty in the magnitude of the term. For example, if one measures a slope $b$ for $\Sigma_{SFR}$ versus $\Sigma_{mol}$ using CO$_{21}$, it would be reasonable and conservative to report a best estimate $b^{\text{corr}}=b/(1-0.15b)$ with an additional component $\Delta b \sim b^{\text{corr}} - b$ added to the systematic uncertainty, e.g., as $^{+0.0}_{-\Delta b}$.

Second, and perhaps even more important, the changes in physical conditions associated with variations in $R_{21}$, $R_{32}$, or $R_{31}$ also imply changes in the CO-to-H$_2$ conversion factor for CO (1–0). As density, temperature, and opacity vary, so will the emissivity of the gas in CO (1–0), generally with the sense that gas with higher $R_{21}$ and $R_{31}$ will show higher emissivity.

Modulo changes in the CO abundance, $N_{CO}/N_{H_2}$, the emissivity of the gas relates directly to $\alpha_{CO}^{1-0}$. Of course, abundance variations do occur, and the amount of “CO-faint” gas represents a major consideration in the behavior of $\alpha_{CO}$ across the galaxy population (e.g., Wolfire et al. 2010; Glover & Mac Low 2011; Leroy et al. 2011; Sandstrom et al. 2013; Bolatto et al. 2013a). In fact, one of the findings in Section 4.2 is that CO excitation appears to increase in low-mass galaxies, where CO-faint gas will be more prevalent. The same diminished shielding that leads to the CO abundance variations may relate to the gas heating (e.g., Peñaloza et al. 2018), so emissivity can change nonlinearly with abundance variations. But for many purposes, it is useful to think about $\alpha_{CO}^{1-0}$ as being a separable problem, with the emissivity of CO and the abundance of CO-faint gas representing distinct factors that, combined, set $\alpha_{CO}^{1-0}$ (e.g., see Bolatto et al. 2013a). We proceed by discussing only the emissivity portion of the problem.

Figure 11 illustrates the behavior of CO (1–0) emissivity $\propto \alpha_{CO}^{-1}$ in the model grid entries that match our measured ratios in Figure 10. Simply knowing the $^{12}$CO line ratios does not determine $\alpha_{CO}$, but the figure shows that variations in $R_{21}$ (or the other ratios) will correlate with variations in $\alpha_{CO}^{1-0}$ as long as other factors, like $N_{CO}/N_{H_2}$ and $N_{CO}/\Delta v$, remain approximately fixed. As discussed in Section 2, $N_{CO}/\Delta v$ relates closely to optical depth and is expected to vary across galaxies and especially to vary some between galaxies and galaxy centers, e.g., in response to a changing dynamical state of the molecular gas or the emergence of a diffuse CO component (e.g., Downes & Solomon 1998; Bolatto et al. 2013a). Still, observations of the $^{13}$CO/$^{12}$CO ratio in normal galaxies show a relatively narrow range (e.g., Cormier et al. 2018, and references therein), implying a relatively stable optical depth, and it may be reasonable to imagine that within any given galaxy, one can focus on a single “color” in Figure 11.

This simple illustration agrees qualitatively with the more detailed discussion and numerical results in Gong et al. (2020). Gong et al. (2020) simulated portions of galaxies with realistic chemistry and radiative transfer to predict scaling relations relating $\alpha_{CO}^{2-1}$ and $\alpha_{CO}^{1-0}$, finding $\alpha_{CO}^{1-0} \propto R_{21}^{-0.9}$, similar to the scaling at fixed $N_{CO}/\Delta v$ in Figure 11.
<figure><img src="image_14.png" /><figcaption>**Figure 11.** The $\alpha_{CO}^{1-0}$ tracks line ratio variations. Shown is the inverse of CO (1–0) emissivity in the model grid, $N_{CO}/I_{CO}$ (1–0), as a function of $\log_{10} R_{21}$. Points show model grid entries that satisfy the mean $R_{32}$ and $R_{21}$ constraints in Figure 10, colored by the column density per line width, $N_{CO}/\Delta v$, of that model. For a fixed CO abundance, $N_{CO}/N_{H_2}$, the inverse of the emissivity is proportional to $\alpha_{CO}^{1-0}$, the CO-to-H$_2$ conversion factor for CO (1–0). The offsets among model families with different $N_{CO}/\Delta v$ reflect that the column per line width directly relates to the opacity of the line. The plot shows that for approximately fixed $N_{CO}/\Delta v$, increasing $\log_{10} R_{21}$ corresponds to decreasing $\alpha_{CO}^{1-0}$, with the proportionality roughly linear (dashed lines) or slightly steeper. As emphasized by Gong et al. (2020), the low-$J$ CO line ratios offer one of our best prospects to empirically map out likely changes in $\alpha_{CO}^{1-0}$.</figcaption></figure>

This relationship between emissivity and the line ratio implies that the trends in Figures 8 and 9 and Tables 5 and 6 can also be interpreted as closely relating to variations in $\alpha_{CO}^{1-0}$, albeit with plenty of caveats. In this case, $\alpha_{CO}^{1-0}$ drops by ~0.2 dex, or a factor of ~1.6, on average, toward the centers of galaxies, and $\alpha_{CO}^{1-0}$ shows gradients of order 0.1–0.2 dex decade$^{-1}$ in $\Sigma_{SFR}$. Note that these variations all refer to $\alpha_{CO}^{1-0}$. Translating, e.g., CO (2–1) intensity into $\Sigma_{mol}$ further requires multiplying by $R_{21}^{-1}$, so that $\Sigma_{mol} \propto \alpha_{CO}^{1-0} R_{21}^{-1} I_{CO}^{2-1} \sim R_{21}^{-2} I_{CO}^{2-1}$. In practice, e.g., as discussed by Gong et al. (2020), there will be distinct relationships between $\alpha_{CO}^{2-1}$ and $R_{21}$ or $R_{31}$. The key point here is not the exact prescription but that we would expect the relationship between $\alpha_{CO}^{2-1}$ or $\alpha_{CO}^{3-2}$ and the line ratios to be steeper than that for $\alpha_{CO}^{1-0}$. The conservative recommendations above would account for this additional steepness by normalizing the scaling relations to effectively operate in terms of CO (1–0), but they would not account for the implied variations of $\alpha_{CO}^{1-0}$.

Qualitatively, this highlights that mapping the low-$J$ CO line ratios offers one of our best options to trace out $\alpha_{CO}^{1-0}$ variations in detail. Despite the ambiguities, these line ratios represent a handle on excitation that can be surveyed across many galaxies. Because of the brightness of CO emission, ALMA or NOEMA can map these lines with high angular resolution; e.g., the PHANGS-ALMA survey prompted this analysis. Though ambiguities in the interpretation of these ratios exist, progress on numerical simulations (e.g., Szűcs et al. 2016; Peñaloza et al. 2018; Gong et al. 2020; Tress et al. 2020) is cause for optimism, and many of the alternative approaches to map out $\alpha_{CO}$ variations, e.g., isotopologue mapping or dust observations, involve a priori unknown additional free parameters, like the isotopologue abundance or coupling between dust opacity, dust-to-gas ratio, and environment (see Bolatto et al. 2013b, for extensive discussion).

### 5.3. Biases, Selection Effects, and Next Directions

To the best of our knowledge, our sample reflects the current state of the literature. However, our analysis also highlights that the field of extragalactic line ratios remains in development. We highlight three important issues, already visible from the analysis above.

Most map pairs target massive galaxies on the star-forming main sequence. As shown in Figures 3, 6, and 7, our measurements still span only a limited range of galaxy properties. To be concrete, galaxies with $R_{21}$ measurements (not limits) span a 16%–84% range of $\log_{10} M_*$ [$M_\odot$] = 10.3–10.8 (median 10.5) and $\log_{10} \text{SFR}/M_*$ [yr$^{-1}$] = −10.5 to −9.8. Those for $R_{32}$ span a 16%–84% range of $\log_{10} M_*$ [$M_\odot$] = 10.1–10.8 (median 10.4) and $\log_{10} \text{SFR}/M_*$ [yr$^{-1}$] = −10.5 to −9.9. And galaxies with measured $R_{31}$ have a 16%–84% range of $\log_{10} M_*$ [$M_\odot$] = 10.2–10.8 (median 10.4) and $\log_{10} \text{SFR}/M_*$ [yr$^{-1}$] = −10.5 to −9.9.

These values represent only a narrow range concentrated near the high-mass end of the star-forming main sequence. Combined with the already low number of measurements, this limits the ability to fit robust scaling relations and contributes to the weakness of the correlations measured between integrated galaxy properties and the line ratios (Section 4.2). An obvious path to make progress will be to expand the set of well-measured, beam-matched line ratios (either from mapping or carefully constructed galaxy-integrated experiments) to include lower-mass, more actively star-forming, and more quiescent galaxies. In the near future, populating the regime of $\log_{10} M_*$ [$M_\odot$] ~ 9.5–10.2 with high-quality measurements would already dramatically improve our understanding of how line ratios vary in the $z=0$ galaxy population (see the sparse coverage in Figure 7).

Many maps that do exist have limited sensitivity, especially CO (3–2). Figures 8 and 9 show many limits at large radii in the stacked profiles of the maps that do exist. For $R_{32}$, there are a large number of unconstraining upper limits at large radii, reflecting the poor sensitivity of the CO (3–2) compared to the CO (2–1) maps. For $R_{21}$, the situation is a bit better but more extended, sensitive CO (1–0) maps are needed. Meanwhile, the radial extent of $R_{31}$ coverage remains very poor; clearly, more sensitive mapping is needed in both transitions.

Calibration issues induce scatter of the same order as the dynamic range of the physical variations in the ratios. In Table 3, the 16%–84% range of $R_{21}$ is 0.22 dex, with $R_{31}$ and $R_{32}$ showing slightly wider ranges of 0.32 and 0.41 dex. As
discussed in Section 3.1 and shown in Figure 2, current calibration uncertainties are imperfectly known but likely 20%–25% for the non-ALMA data. This uncertainty in the calibration alone introduces ≈0.1–0.15 dex rms scatter in the line ratios. Expressed as a 16%–84% range, this is 0.2–0.3 dex, of the same order as the range in line ratios themselves. The best path forward appears to be larger internally consistent or carefully cross-calibrated data sets.

## 6. Summary

We combine a large set of publicly available maps of CO (1–0), CO (2–1), and CO (3–2) emission from nearby galaxies with the new PHANGS-ALMA CO (2–1) survey to measure low-*J* CO line ratios for 152 CO map pairs (see Table 1). The full sample spans $M_* \sim 10^9-10^{11} M_\odot$ but consists mostly of relatively massive ($\log_{10} M_*$ [$M_\odot$] $\approx 10.25-11.0$) star-forming (SFR $\sim 1-5 M_\odot$ yr$^{-1}$) galaxies that lie near the star-forming main sequence (see Figure 3).

These maps of low-*J* CO line ratios across normal galaxies complement earlier detailed multitransition studies (e.g., Israel & Baas 2001, 2003; Papadopoulos et al. 2002; Bayet et al. 2004, 2006; Papadopoulos et al. 2012; Kamenetzky et al. 2014, 2017; Israel et al. 2015), which tended to focus on galaxy-integrated measurements, galaxy centers, and starburst galaxies but often employed a wider set of transitions than we have available here, including CO isotopologues and higher *J* lines. We also extend earlier mapping work on samples by Wilson et al. (2012), Leroy et al. (2013b), Yajima et al. (2021), and den Brok et al. (2021), which used smaller samples and earlier versions of some of the same data and focused mostly on a single line ratio. For this study, a key addition and the motivating data set is the new PHANGS-ALMA CO (2–1) mapping survey (Leroy et al. 2021a).

Tables 2–4 report our results for whole galaxies.

1. Integrating over whole map pairs, we find $R_{21} = 0.50-0.83$ with mean 0.65, $R_{32} = 0.22-0.58$ with mean 0.47, and $R_{31} = 0.16-0.44$ with mean 0.32. We compare these to literature measurements of the same ratios by studying nearby galaxies (Section 4.1, Figure 4) and find overall consistency, though our $R_{31}$ values appear somewhat higher than previously reported for nearby galaxies. Area-matched $R_{32}$ and $R_{31}$ measurements for nearby galaxies remain relatively scarce, dominated by the JCMT NGLS (Wilson et al. 2012), and we note this as a productive area for new observations.

We search for correlations between the low-*J* CO line ratios and integrated galaxy properties (Section 4.2, Figures 5 and 6, Table 4). Such correlations remain hard to discern, partially due to the limited diversity of galaxies with measured ratios (Figure 3) and partially because the calibration uncertainties associated with the data are of the same order as the dynamic range of the ratio in the nearby galaxy population (Section 3, Figure 2). A secondary issue is that some of the most physically meaningful comparisons (e.g., between the line ratios and SFR per CO) involve correlated axes that can strongly influence the inferred trends (see Section 4.1). Despite this:

2. We identify a consistent set of marginally significant correlations between all three line ratios, $R_{21}$, $R_{32}$, and

$R_{31}$, and quantities that trace normalized star formation activity. The line ratios anticorrelate with stellar mass ($M_*$) and CO luminosity ($L_{CO}$) and positively correlate with specific SFR ($SFR/M_*$) and SFR per CO ($SFR/L_{CO}$). These correlations have the sense that both low-mass dwarf galaxies and starburst galaxies should show high line ratios. This agrees with the physical expectations and previous observations that dwarf galaxies have high-excitation, poorly shielded molecular gas and that starburst galaxies have high CO excitation.

We measure local variations of each line ratio within the galaxies (Section 4.3, Figures 8 and 9, Tables 5 and 6). In this analysis, we control for galaxy-to-galaxy scatter and global calibration uncertainties by focusing on the line ratio normalized to the galaxy average. We examine how line ratio variations within galaxies correlate with galactocentric radius, local star formation surface density, and local specific SFR, also normalized to the galaxy average. We find:

3. Most galaxies with measurements show enhanced values of all ratios in the central 1 kpc or $0.5r_{eff}$ wide bin in our analysis (Table 6). Here $R_{31}$ shows the strongest central enhancements, 0.27 dex, on average, followed by $R_{21}$ with median 0.18 dex and $R_{32}$ with 0.08 dex. These central enhancements have been noted before, especially in $R_{21}$, but this study represents the largest systematic measurement for all three lines to date. The sense of these variations agrees with the expectation that, compared to the gas in galaxy disks, the gas in galaxy centers is denser and more actively star-forming, can be heated by active galactic nuclei, and perhaps includes an optically thinner diffuse component due to high turbulence.

4. Within galaxies, all three line ratios also show significant internal gradients as a function of radius and $\Sigma_{SFR}$. The $R_{21}$ and $R_{31}$ also show significant gradients as a function of $\Sigma_{SFR}/\Sigma_*$. We report fits to the magnitude of these gradients in Table 5. This behavior agrees with the expectation that the more active parts of galaxies host hotter, denser gas and reflects the same underlying trend as the central enhancements.

Finally, we note some implications of our measurements:

6. We consider a set of simple non-LTE models that treat density distributions (Section 2 and Appendix A) and note the physical conditions in the models that satisfy our measurements of all three lines (Section 5.1). Our observed ratios can be broadly reproduced by cold gas with moderate density and intermediate column density per line width, $N_{CO}/\Delta v$.

7. Following Gong et al. (2020) and illustrated using our own model grid, we highlight that these line ratio variations also imply corresponding variations in $\alpha_{CO}^{1-0}$ because the changes in the physical conditions tracked by the CO line ratio imply changes in the CO emissivity (Section 5.2). These $\alpha_{CO}$ variations will compound with those caused by “CO-faint” gas to produce the overall variations of $\alpha_{CO}$ in galaxies.

Together, these results paint a basic picture of how the low-*J* CO line ratios vary across the local galaxy population, with dwarf and starburst galaxies showing enhanced excitation and galaxy centers and high-$\Sigma_{SFR}$ regions also being more excited.
We close by emphasizing that much more work is needed. The major limitations of the present study, and the field in general, are S/N in the CO (1–0) and CO (3–2) data, uncertainties in the calibration, and the limited sampling of the full galaxy population. The PHANGS-ALMA CO (2–1) maps offer a high-quality, high-S/N, well-calibrated starting point that spans the local galaxy population. One next major step will be to measure a large, diverse sample of galaxies in all three lines with well-controlled calibration, excellent S/N, and good resolution. This will sharpen our knowledge of integrated galaxy scaling relations and allow us to investigate resolved line excitation variations across galaxies without the need for stacking or the aggressive normalizations performed here.

We thank the anonymous referee for a constructive and thorough report that significantly improved the quality of this work.

This work was carried out as part of the PHANGS collaboration.

The work of A.K.L. and J.S. was partially supported by the National Science Foundation (NSF) under grant Nos. 1615105 and 1653300, as well as the National Aeronautics and Space Administration (NASA) under ADAP grants NNX16AF48G and NNX17AF39G.

E.R. acknowledges the support of the Natural Sciences and Engineering Research Council of Canada (NSERC), funding reference No. RGPIN-2017-03987, and computational support from Compute Canada.

A.U. acknowledges support from Spanish funding grants PGC2018-094671-B-I00 (MCIU/AEI/FEDER) and PID2019-108765GB-I00 (MICINN).

K.S. and I.D.C. acknowledge funding support from National Science Foundation grant No. 1615728 and NASA ADAP grants NNX16AF48G and NNX17AF39G.

E.S., H.A.P., T.S., and T.G.W. acknowledge funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No. 694343).

J.Pe. acknowledges support by the Programme National "Physique et Chimie du Milieu Interstellaire" (PCMI) of CNRS/INSU with INC/INP, cofunded by CEA and CNES.

M.C. and J.M.D.K. gratefully acknowledge funding from the Deutsche Forschungsgemeinschaft (DFG) in the form of the Emmy Noether Research Group (grant No. KR4801/1-1) and the DFG Sachbeihilfe (grant No. LR4801/2-1), as well as the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program via the ERC Starting Grant MUSTANG (grant agreement No. 714907).

R.S.K. and S.C.O.G. acknowledge funding from the European Research Council via the ERC Synergy Grant "ECOGAL—Understanding our Galactic ecosystem: From the disk of the Milky Way to the formation sites of stars and planets" (project ID 855130). They also acknowledge support from the DFG via the Collaborative Research Center (SFB 881, Project ID 138713538) "The Milky Way System" (subprojects A1, B1, B2, and B8) and the Heidelberg cluster of excellence (EXC 2181-390900948) "STRUCTURES: A unifying approach to emergent phenomena in the physical world, mathematics, and complex data," funded by the German Excellence Strategy.

A.T.B., F.B., and J.d.B. would like to acknowledge funding from the European Research Council (ERC) under the

European Union's Horizon 2020 research and innovation program (grant agreement No.726384/Empire).

C.E. acknowledges funding from the Deutsche Forschungsgemeinschaft (DFG) Sachbeihilfe, grant No. BI1546/3-1.

M.Q. acknowledges support from the research project PID2019-106027GA-C44 from the Spanish Ministerio de Ciencia e Innovación.

A.S. is supported by an NSF Astronomy and Astrophysics Postdoctoral Fellowship under award AST-1903834.

This work is based on observations carried out under project Nos. 169-06, 053-07, 122-07, 160-06, 218-05, 058-08, 212-08, 196-13, 078-14, and 190-14 with the IRAM 30 m telescope. IRAM is supported by INSU/CNRS (France), MPG (Germany), and IGN (Spain).

This work is based on COMING and the Nobeyama Atlas of Nearby Spiral Galaxies, which are both legacy programs of the Nobeyama 45 m radio telescope, which is operated by Nobeyama Radio Observatory, a branch of the National Astronomical Observatory of Japan. We gratefully acknowledge both teams for making their data public and acknowledge the hard work of the teams and staff in obtaining the data.

We also thank the JCMT NGLS team and the JCMT staff for their hard work and making their data public, and we acknowledge helpful correspondence with the JCMT staff regarding calibration uncertainties. The James Clerk Maxwell Telescope is operated by the Joint Astronomy Centre on behalf of the Science and Technology Facilities Council of the United Kingdom, the Netherlands Organisation for Scientific Research, and the National Research Council of Canada.

This paper makes use of the following ALMA data, which have been processed as part of the PHANGS-ALMA CO (2–1) survey:

ADS/JAO.ALMA#2012.1.00650.S
ADS/JAO.ALMA#2013.1.00803.S
ADS/JAO.ALMA#2013.1.01161.S
ADS/JAO.ALMA#2015.1.00121.S
ADS/JAO.ALMA#2015.1.00782.S
ADS/JAO.ALMA#2015.1.00925.S
ADS/JAO.ALMA#2015.1.00956.S
ADS/JAO.ALMA#2016.1.00386.S
ADS/JAO.ALMA#2017.1.00392.S
ADS/JAO.ALMA#2017.1.00766.S
ADS/JAO.ALMA#2017.1.00886.L
ADS/JAO.ALMA#2018.1.01321.S
ADS/JAO.ALMA#2018.1.01651.S
ADS/JAO.ALMA#2018.A.00062.S
ADS/JAO.ALMA#2019.1.01235.S
ADS/JAO.ALMA#2019.2.00129.S

ALMA is a partnership of ESO (representing its member states), NSF (USA), and NINS (Japan), together with NRC (Canada), NSC and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO, and NAOJ. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.

## Appendix A

### CO Line Emission from Density Distributions

Table 7 presents predicted line ratios, $R_{21}$, $R_{32}$, and $R_{31}$, and CO (1–0) emissivity, $\epsilon_{1-0}$, for models with distributions of densities and a fixed temperature, $T_{\text{kin}}$, and column of CO per
<table><caption>**Table 7**
CO Emission from Multidensity Models</caption><thead><tr><th>T<sub>kin</sub><br>(K)</th><th>N<sub>CO</sub>/Δν<br>(cm<sup>-2</sup> (km s<sup>-1</sup>)<sup>-1</sup>)</th><th>n<sub>0,H2</sub><br>(cm<sup>-1</sup>)</th><th>σ<br>(dex)</th><th>ε<sub>1-0</sub><br>(K km s<sup>-1</sup> (cm<sup>-2</sup>)<sup>-1</sup>)</th><th>R<sub>21</sub></th><th>R<sub>32</sub></th><th>R<sub>31</sub></th></tr></thead><tbody><tr><td>10</td><td>1.000E+15</td><td>3.160E+01</td><td>0.2</td><td>7.779E-21</td><td>0.17</td><td>0.05</td><td>0.01</td></tr><tr><td>10</td><td>1.000E+15</td><td>3.160E+01</td><td>0.3</td><td>9.791E-21</td><td>0.18</td><td>0.05</td><td>0.01</td></tr><tr><td>10</td><td>1.000E+15</td><td>3.160E+01</td><td>0.4</td><td>1.317E-20</td><td>0.19</td><td>0.05</td><td>0.01</td></tr><tr><td>10</td><td>1.000E+15</td><td>3.160E+01</td><td>0.5</td><td>1.818E-20</td><td>0.22</td><td>0.06</td><td>0.01</td></tr><tr><td>10</td><td>1.000E+15</td><td>3.160E+01</td><td>0.6</td><td>2.490E-20</td><td>0.26</td><td>0.07</td><td>0.02</td></tr><tr><td>10</td><td>1.000E+15</td><td>3.160E+01</td><td>0.7</td><td>3.306E-20</td><td>0.32</td><td>0.10</td><td>0.03</td></tr><tr><td>10</td><td>1.000E+15</td><td>3.160E+01</td><td>0.8</td><td>4.208E-20</td><td>0.40</td><td>0.12</td><td>0.05</td></tr><tr><td>10</td><td>1.000E+15</td><td>3.160E+01</td><td>0.9</td><td>5.130E-20</td><td>0.48</td><td>0.16</td><td>0.07</td></tr><tr><td>10</td><td>1.000E+15</td><td>3.160E+01</td><td>1.0</td><td>6.012E-20</td><td>0.56</td><td>0.19</td><td>0.11</td></tr><tr><td>10</td><td>1.000E+15</td><td>3.160E+01</td><td>1.1</td><td>6.811E-20</td><td>0.66</td><td>0.23</td><td>0.15</td></tr></tbody></table>

Note. This table is a stub. The full version of the table appears as a machine-readable table in the online version of the paper. Columns are as follows: T<sub>kin</sub>—kinetic temperature of all zones in the model; N<sub>CO</sub>/Δν—column of CO per line width in all zones in the model; n<sub>0,H2</sub>—mean collider density for the lognormal distribution of densities in the model; σ—rms width σ of the lognormal distribution of densities in the model; ε<sub>1-0</sub>—emissivity of H<sub>2</sub> in the CO (1–0) transition assuming a fixed CO/H<sub>2</sub> abundance of 1E-4; R<sub>21</sub>—CO (2–1)/CO (1–0) line ratio for the model; R<sub>32</sub>—CO (3–2)/CO (2–1) line ratio for the model; R<sub>31</sub>—CO (3–2)/CO (1–0) line ratio for the model, and see Leroy et al. (2017) for more details regarding the calculation.

(This table is available in its entirety in machine-readable form.)

unit line width, N<sub>CO</sub>/Δν. Following Leroy et al. (2017), we consider lognormal density distributions with a mean density, n<sub>0</sub>, and width in log density space, σ. The predicted ratios represent the sum over an ensemble of zones that share T<sub>kin</sub> and N<sub>CO</sub>/Δν but each have a distinct density.

The methods mostly follow Leroy et al. (2017) with atomic data from LAMDA (Schöier et al. 2005) and calculations via RADEX (van der Tak et al. 2007). Distinct from Leroy et al. (2017), we assume a single fixed N<sub>CO</sub>/Δν across all zones. That paper assumed a fixed optical depth, τ, and selected N/Δν to match that opacity. The current approach is better suited to treat multitransition measurements like CO line ratios, while adopting a fixed τ as in Leroy et al. (2017) may be more appropriate to single-transition studies, e.g., HCN (1–0) only, where τ may be at least roughly constrained from isotopologue studies (e.g., Jiménez-Donaire et al. 2017). We note that for a fixed abundance, N<sub>CO</sub>/N<sub>H2</sub>, N<sub>CO</sub>/Δν will relate to the collider density n<sub>H2</sub> within a layer via a combination of size and line width that can vary from zone to zone. Since we model only a single species here, there is no implied inconsistency other than the zones having variable size or structure. This makes no difference to our modeling, which does not consider zone size as an important variable. However, as discussed by Leroy et al. (2017) when modeling multiple species, fixing the optical depth can lead to implied zone-to-zone abundance variations.

The models have the advantage of incorporating a realistic mixture of densities, which is certainly present in any low-resolution observation of galaxies. They have a suite of caveats, discussed at length in Leroy et al. (2017). Here we only note that we have assumed that the zones all share a fixed T<sub>kin</sub> and N<sub>CO</sub>/Δν, and they do not “shadow” one another; i.e., we observe the linear combination of emission from all zones. We also only consider a lognormal density distribution in this paper because the modeling is not a central focus of this work. Beyond this, all of the usual caveats related to RADEX

modeling apply, and we refer the reader to van der Tak et al. (2007) for more details.

## Appendix B
### Illustration of the Effect of the CMB and the Rayleigh-Jeans Approximation

Figure 12 illustrates the impact of the CMB and use of the Rayleigh-Jeans approximation on measured CO line ratios. These effects are well known, but we are not aware of a clean illustration of the impact of both on this set of ratios, and we found these plots useful to interpret our measurements, so we include them here (as above, we note discussions in Eckart et al. 1990; da Cunha et al. 2013; Bolatto et al. 2013a; Zschaechner et al. 2018).

The right panel shows the effect of measuring line ratios in contrast against the 2.73 K CMB (Fixsen et al. 1996). The effect of the CMB is to selectively suppress low-J emission, leading the measured ratio to be higher than the true value expected for the source without any CMB. By contrast, the effect of using the Rayleigh-Jeans approximation, shown in the left panel, is that for sources with temperatures in the range of real molecular clouds, we expect “thermal” line ratios (i.e., the value for an opaque source in LTE) <1, with lower values for colder objects. The two effects somewhat cancel out.

As noted above, the magnitude of the CMB effect on real observations of galaxies can be difficult to gauge without information on the small-scale structure of the emission. The radiative transfer involving the CMB plays out on the scale of individual clouds, and the intensity is then subject to a large beam dilution effect before entering the sort of measurements presented in this paper. Both effects are accounted for by models (e.g., RADEX; van der Tak et al. 2007), but coupling those models to observations can require estimates (or the addition of a free parameter) of beam filling.
<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

**Figure 12.** Effect of the CMB and Rayleigh-Jeans (RJ) approximation on measured ratios. For reference, we illustrate the effect of the RJ approximation (left) and CMB (right) for line ratios measured from perfectly opaque sources ($\tau \gg 1$). The left panel shows the expected $R_{\text{true}}$ for a thermalized, opaque source when R is expressed in RJ brightness temperature. The increased inaccuracy of the RJ approximation at low temperature leads to “thermal” line ratios with values of $R_{\text{true}} < 1$. The right panel shows how the measured $R_{\text{meas}}$ deviates from $R_{\text{true}}$ due to the effects of the CMB, plotting $R_{\text{true}}/R_{\text{meas}}$ as a function of the excitation temperature for a CMB temperature of 2.73 K. The CMB tends to increase the observed line ratio because it affects low-$J$ lines more. The two effects act in opposite directions, with $R_{\text{true}} < 1$ but then $R_{\text{meas}}$ increased by the effects of the CMB. Both are well handled by modeling but important to bear in mind when interpreting measured line ratios. Correcting observations with low physical resolution for the effects of the CMB can be particularly challenging.

## ORCID iDs

Adam K. Leroy https://orcid.org/0000-0002-2545-1700

Erik Rosolowsky https://orcid.org/0000-0002-5204-2259

Antonio Usero https://orcid.org/0000-0003-1242-505X

Karin Sandstrom https://orcid.org/0000-0002-4378-8534

Eva Schinnerer https://orcid.org/0000-0002-3933-7677

Alberto D. Bolatto https://orcid.org/0000-0002-5480-5686

Jiayi Sun (孙嘉懿) https://orcid.org/0000-0003-0378-4667

Ashley. T. Barnes https://orcid.org/0000-0003-0410-4504

Francesco Belfiore https://orcid.org/0000-0002-2545-5752

Frank Bigiel https://orcid.org/0000-0003-0166-9745

Jakob S. den Brok https://orcid.org/0000-0002-8760-6157

Yixian Cao https://orcid.org/0000-0001-5301-1326

I-Da Chiang (江宜達) https://orcid.org/0000-0003-2551-7148

Mélanie Chevance https://orcid.org/0000-0002-5635-5180

Daniel A. Dale https://orcid.org/0000-0002-5782-9093

Cosima Eibensteiner https://orcid.org/0000-0002-1185-2810

Christopher M. Faesi https://orcid.org/0000-0001-5310-467X

Simon C. O. Glover https://orcid.org/0000-0001-6708-1317

Annie Hughes https://orcid.org/0000-0002-9181-1161

María J. Jiménez Donaire https://orcid.org/0000-0002-9165-8080

Ralf S. Klessen https://orcid.org/0000-0002-0560-3172

Eric W. Koch https://orcid.org/0000-0001-9605-780X

J. M. Diederik Kruijssen https://orcid.org/0000-0002-8804-0212

Daizhong Liu https://orcid.org/0000-0001-9773-7479

Sharon E. Meidt https://orcid.org/0000-0002-6118-4048

Hsi-An Pan https://orcid.org/0000-0002-1370-6964

Jérôme Pety https://orcid.org/0000-0003-3061-6546

Johannes Puschnig https://orcid.org/0000-0003-1111-3951

Miguel Querejeta https://orcid.org/0000-0002-0472-1011

Toshiki Saito https://orcid.org/0000-0002-2501-9328

Amy Sardone https://orcid.org/0000-0002-5783-145X

Elizabeth J. Watkins https://orcid.org/0000-0002-7365-5791

Axel Weiss https://orcid.org/0000-0003-4678-3939

Thomas G. Williams https://orcid.org/0000-0002-0012-2142

## References

Anand, G. S., Lee, J. C., Van Dyk, S. D., et al. 2021, MNRAS, 501, 3621

Ao, Y., Henkel, C., Menten, K. M., et al. 2013, A&A, 550, A135

Bayet, E., Gerin, M., Phillips, T. G., & Contursi, A. 2004, A&A, 427, 45

Bayet, E., Gerin, M., Phillips, T. G., & Contursi, A. 2006, A&A, 460, 467

Bemis, A., & Wilson, C. D. 2019, AJ, 157, 131

Bisbas, T. G., Schruba, A., & van Dishoeck, E. F. 2019, MNRAS, 485, 3097

Bolatto, A. D., Leroy, A., Israel, F. P., & Jackson, J. M. 2003, ApJ, 595, 167

Bolatto, A. D., Warren, S. R., Leroy, A. K., et al. 2013b, Natur, 499, 450

Bolatto, A. D., Wolfire, M., & Leroy, A. K. 2013a, ARA&A, 51, 207

Bolatto, A. D., Wong, T., Utomo, D., et al. 2017, ApJ, 846, 159

Braine, J., & Combes, F. 1992, A&A, 264, 433

Braine, J., Combes, F., Casoli, F., et al. 1993, A&AS, 97, 887

Burkhart, B. 2018, ApJ, 863, 118

Cao, Y., Wong, T., Xue, R., et al. 2017, ApJ, 847, 33

Carilli, C. L., & Walter, F. 2013, ARA&A, 51, 105

Casasola, V., Bianchi, S., De Vis, P., et al. 2020, A&A, 633, A100

Catinella, B., Saintonge, A., Janowiecki, S., et al. 2018, MNRAS, 476, 875

Chabrier, G. 2003, PASP, 115, 763

Chown, R., Li, C., Parker, L., et al. 2021, MNRAS, 500, 1261

Cicone, C., Bothwell, M., Wagg, J., et al. 2017, A&A, 604, A53

Cormier, D., Bigiel, F., Jiménez-Donaire, M. J., et al. 2018, MNRAS, 475, 3909

Crosthwaite, L. P., & Turner, J. L. 2007, AJ, 134, 1827

Currie, M. J., Berry, D. S., Jenness, T., et al. 2014, in ASP Conf. Ser. 485, Astronomical Data Analysis Software and Systems XXIII, ed. N. Manset & P. Forshay (San Francisco, CA: ASP), 391

da Cunha, E., Groves, B., Walter, F., et al. 2013, ApJ, 766, 13

den Brok, J. S., Chatzigiannakis, D., Bigiel, F., et al. 2021, MNRAS, 504, 3221

Donovan Meyer, J., Koda, J., Momose, R., et al. 2013, ApJ, 772, 107

Downes, D., & Solomon, P. M. 1998, ApJ, 507, 615

Druard, C., Braine, J., Schuster, K. F., et al. 2014, A&A, 567, A118

Eckart, A., Downes, D., Genzel, R., et al. 1990, ApJ, 348, 434
Emsellem, E., Wolk, S. J., Walter, F. M., et al. 1999, AJ, submitted
Fazio, G. G., Hora, J. L., Allen, L. E., et al. 2004, ApJS, 154, 10
Federrath, C., & Klessen, R. S. 2013, ApJ, 763, 51
Fixsen, D. J., Cheng, E. S., Gales, J. M., et al. 1996, ApJ, 473, 576
Gallagher, M. J., Leroy, A. K., Bigiel, F., et al. 2018a, ApJ, 858, 90
Gallagher, M. J., Leroy, A. K., Bigiel, F., et al. 2018b, ApJL, 868, L38
García-Burillo, S., Combes, F., Usero, A., et al. 2014, A&A, 567, A125
Ginsburg, A., Henkel, C., Ao, Y., et al. 2016, A&A, 586, A50
Glover, S. C. O., & Clark, P. C. 2012, MNRAS, 426, 377
Glover, S. C. O., & Mac Low, M. 2011, MNRAS, 412, 337
Gong, M., Ostriker, E. C., Kim, C.-G., & Kim, J.-G. 2020, ApJ, 903, 142
Groves, B., Krause, O., Sandstrom, K., et al. 2012, MNRAS, 426, 892
Güsten, R., Baryshev, A., Bell, A., et al. 2008, Proc. SPIE, 7020, 702010
Helfer, T. T., Thornley, M. D., Regan, M. W., et al. 2003, ApJS, 145, 259
Hirota, A., Egusa, F., Baba, J., et al. 2018, PASJ, 70, 73
Hodge, J. A., & da Cunha, E. 2020, RSOS, 7, 200556
Hu, C.-Y., Sternberg, A., & van Dishoeck, E. F. 2021, ApJ, 920, 44
Israel, F. P. 2020, A&A, 635, A131
Israel, F. P., & Baas, F. 2001, A&A, 371, 433
Israel, F. P., & Baas, F. 2003, A&A, 404, 495
Israel, F. P., Rosenberg, M. J. F., & van der Werf, P. 2015, A&A, 578, A95
Janowiecki, S., Catinella, B., Cortese, L., et al. 2017, MNRAS, 466, 4795
Jiménez-Donaire, M. J., Bigiel, F., Leroy, A. K., et al. 2017, MNRAS, 466, 49
Jiménez-Donaire, M. J., Bigiel, F., Leroy, A. K., et al. 2019, ApJ, 880, 127
Kamenetzky, J., Rangwala, N., & Glenn, J. 2017, MNRAS, 471, 2917
Kamenetzky, J., Rangwala, N., Glenn, J., Maloney, P. R., & Conley, A. 2014,
ApJ, 795, 174
Kennicutt, R. C., & Evans, N. J. 2012, ARA&A, 50, 531
Kepley, A. A., Leroy, A. K., Johnson, K. E., Sandstrom, K., & Chen, C. H. R.
2016, ApJ, 828, 50
Kewley, L. J., & Ellison, S. L. 2008, ApJ, 681, 1183
Klessen, R. S., & Glover, S. C. O. 2016, SAAS, 43, 85
Koda, J., Sawada, T., Sakamoto, K., et al. 2020, ApJL, 890, L10
Koda, J., Scoville, N., Hasegawa, T., et al. 2012, ApJ, 761, 41
Krieger, N., Ott, J., Beuther, H., et al. 2017, ApJ, 850, 77
Krumholz, M. R., & Thompson, T. A. 2007, ApJ, 669, 289
Kuno, N., Sato, N., Nakanishi, H., et al. 2007, PASJ, 59, 117
Lamperti, I., Saintonge, A., Koss, M., et al. 2020, ApJ, 889, 103
Lang, P., Meidt, S. E., Rosolowsky, E., et al. 2020, ApJ, 897, 122
Law, C. J., Zhang, Q., Ricci, L., et al. 2018, ApJ, 865, 17
Lequeux, J., Le Bourlot, J., Des Forets, G. P., et al. 1994, A&A, 292, 371
Leroy, A. K., Bigiel, F., de Blok, W. J. G., et al. 2012, AJ, 144, 3
Leroy, A. K., Bolatto, A., Bot, C., et al. 2009, ApJ, 702, 352
Leroy, A. K., Bolatto, A., Gordon, K., et al. 2011, ApJ, 737, 12
Leroy, A. K., Hughes, A., Liu, D., et al. 2021b, ApJS, 255, 19
Leroy, A. K., Lee, C., Schruba, A., et al. 2013a, ApJL, 769, L12
Leroy, A. K., Sandstrom, K. M., Lang, D., et al. 2019, ApJS, 244, 24
Leroy, A. K., Schinnerer, E., Hughes, A., et al. 2021a, ApJS, 257, 43
Leroy, A. K., Usero, A., Schruba, A., et al. 2017, ApJ, 835, 217
Leroy, A. K., Walter, F., Sandstrom, K., et al. 2013b, AJ, 146, 19
Liu, D., Daddi, E., Schinnerer, E., et al. 2021, ApJ, 909, 56
Mangum, J. G., Darling, J., Henkel, C., et al. 2013, ApJ, 779, 33
Mangum, J. G., Emerson, D. T., & Greisen, E. W. 2007, A&A, 474, 679
Mao, R.-Q., Schulz, A., Henkel, C., et al. 2010, ApJ, 724, 1336
Martin, D. C., Fanson, J., Schiminovich, D., et al. 2005, ApJL, 619, L1
Mauersberger, R., Henkel, C., Walsh, W., & Schulz, A. 1999, A&A, 341, 256

Muñoz-Mateos, J. C., Sheth, K., Regan, M., et al. 2015, ApJS, 219, 3
Murphy, E. J., Bolatto, A., Chatterjee, S., et al. 2018, in ASP Conf. Ser. 517, Science with a Next Generation Very Large Array, ed. E. Murphy (San Francisco, CA: ASP), 3
Onus, A., Krumholz, M. R., & Federrath, C. 2018, MNRAS, 479, 1702
Padoan, P., & Nordlund, Å. 2002, ApJ, 576, 870
Papadopoulos, P. P., & Seaquist, E. R. 1999, ApJ, 516, 114
Papadopoulos, P. P., Thi, W.-F., & Viti, S. 2002, ApJ, 579, 270
Papadopoulos, P. P., van der Werf, P., Xilouris, E., Isaak, K. G., & Gao, Y. 2012, ApJ, 751, 10
Peñaloza, C. H., Clark, P. C., Glover, S. C. O., & Klessen, R. S. 2018, MNRAS, 475, 1508
Peñaloza, C. H., Clark, P. C., Glover, S. C. O., Shetty, R., & Klessen, R. S. 2017, MNRAS, 465, 2277
Rémy-Ruyer, A., Madden, S. C., Galliano, F., et al. 2014, A&A, 563, A31
Roman-Duval, J., Israel, F. P., Bolatto, A., et al. 2010, A&A, 518, L74
Rosolowsky, E., Hughes, A., Leroy, A. K., et al. 2021, MNRAS, 502, 1218
Rosolowsky, E., & Leroy, A. 2006, PASP, 118, 590
Rubio, M., Elmegreen, B. G., Hunter, D. A., et al. 2015, Natur, 525, 218
Saintonge, A., Catinella, B., Tacconi, L. J., et al. 2017, ApJS, 233, 22
Salim, S., Boquien, M., & Lee, J. C. 2018, ApJ, 859, 11
Salim, S., Lee, J. C., Janowiecki, S., et al. 2016, ApJS, 227, 2
Salim, S., Rich, R. M., Charlot, S., et al. 2007, ApJS, 173, 267
Sandstrom, K. M., Leroy, A. K., Walter, F., et al. 2013, ApJ, 777, 5
Schöier, F. L., van der Tak, F. F. S., van Dishoeck, E. F., & Black, J. H. 2005, A&A, 432, 369
Schruba, A., Leroy, A. K., Walter, F., et al. 2012, AJ, 143, 138
Schuster, K. F., Kramer, C., Hitschfeld, M., García-Burillo, S., & Mookerjea, B. 2007, A&A, 461, 143
Sheth, K., Regan, M., Hinz, J. L., et al. 2010, PASP, 122, 1397
Shirley, Y. L. 2015, PASP, 127, 299
Sorai, K., Kuno, N., Muraoka, K., et al. 2019, PASJ, 71, S14
Sun, J., Leroy, A. K., Schinnerer, E., et al. 2020, ApJL, 901, L8
Sun, J., Leroy, A. K., Schruba, A., et al. 2018, ApJ, 860, 172
Szűcs, L., Glover, S. C. O., & Klessen, R. S. 2016, MNRAS, 460, 82
Tacconi, L. J., Genzel, R., & Sternberg, A. 2020, ARA&A, 58, 157
Tremonti, C. A., Heckman, T. M., Kauffmann, G., et al. 2004, ApJ, 613, 898
Tress, R. G., Smith, R. J., Sormani, M. C., et al. 2020, MNRAS, 492, 2973
Ueda, J., Iono, D., Petitpas, G., et al. 2012, ApJ, 745, 65
Usero, A., Leroy, A. K., Walter, F., et al. 2015, AJ, 150, 115
van der Tak, F. F. S., Black, J. H., Schöier, F. L., Jansen, D. J., & van Dishoeck, E. F. 2007, A&A, 468, 627
van Dishoeck, E. F., & Black, J. H. 1988, ApJ, 334, 771
Vlahakis, C., van der Werf, P., Israel, F. P., & Tilanus, R. P. J. 2013, MNRAS, 433, 1837
Weiß, A., Walter, F., & Scoville, N. Z. 2005, A&A, 438, 533
Wilson, C. D., Petitpas, G. R., Iono, D., et al. 2008, ApJS, 178, 189
Wilson, C. D., Warren, B. E., Israel, F. P., et al. 2009, ApJ, 693, 1736
Wilson, C. D., Warren, B. E., Israel, F. P., et al. 2012, MNRAS, 424, 3050
Wolfire, M. G., Hollenbach, D., & McKee, C. F. 2010, ApJ, 716, 1191
Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868
Yajima, Y., Sorai, K., Miyamoto, Y., et al. 2021, PASJ, 73, 257
Young, J. S., & Scoville, N. Z. 1991, ARA&A, 29, 581
Young, J. S., Xie, S., Tacconi, L., et al. 1995, ApJS, 98, 219
Zschaechner, L. K., Bolatto, A. D., Walter, F., et al. 2018, ApJ, 867, 111