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

# Millimeter-wave Line Ratios and Sub-beam Volume Density Distributions

Adam K. Leroy¹, Antonio Usero², Andreas Schruba³, Frank Bigiel⁴, J. M. Diederik Kruijssen⁵,⁶, Amanda Kepley⁷,
Guillermo A. Blanc⁸,⁹,¹⁰, Alberto D. Bolatto¹¹, Diane Cormier⁶, Molly Gallagher¹, Annie Hughes¹²,¹³, Maria J. Jiménez-Donaire⁶,
Erik Rosolowsky¹⁴, and Eva Schinnerer⁵

¹Department of Astronomy, The Ohio State University, 140 West 18th Avenue, Columbus, OH 43210, USA

²Observatorio Astronmico Nacional (IGN), C/Alfonso XII, 3, E-28014 Madrid, Spain

³Max-Planck-Institut für extraterrestrische Physik, Giessenbachstraße 1, D-85748 Garching, Germany

⁴Institute für theoretische Astrophysik, Zentrum für Astronomie der Universität Heidelberg, Albert-UEberle Str. 2, D-69120 Heidelberg, Germany

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

⁶Max Planck Institute für Astronomie, Königstuhl 17, D-69117, Heidelberg, Germany

⁷National Radio Astronomy Observatory, 520 Edgemont Road, Charlottesville, VA 22903, USA

⁸Departamento de Astronomía, Universidad de Chile, Casilla 36-D, Santiago, Chile

⁹Centro de Astrofísica y Tecnologías Afines (CATA), Camino del Observatorio 1515, Las Condes, Santiago, Chile

¹⁰Visiting Astronomer, Observatories of the Carnegie Institution for Science, 813 Santa Barbara St, Pasadena, CA 91101, USA

¹¹Department of Astronomy, Laboratory for Millimeter-wave Astronomy, and Joint Space Institute, University of Maryland, College Park, MD 20742, USA

¹²CNRS, IRAP, 9 av. du Colonel Roche, BP 44346, F-31028 Toulouse cedex 4, France

¹³Université de Toulouse, UPS-OMP, IRAP, F-31028 Toulouse cedex 4, France

¹⁴Department of Physics, University of Alberta, Edmonton, AB, Canada

Received 2016 September 18; revised 2016 November 10; accepted 2016 November 11; published 2017 January 31

## Abstract

We explore the use of mm-wave emission line ratios to trace molecular gas density when observations integrate over a wide range of volume densities within a single telescope beam. For observations targeting external galaxies, this case is unavoidable. Using a framework similar to that of Krumholz & Thompson, we model emission for a set of common extragalactic lines from lognormal and power law density distributions. We consider the median density of gas that produces emission and the ability to predict density variations from observed line ratios. We emphasize line ratio variations because these do not require us to know the absolute abundance of our tracers. Patterns of line ratio variations have the potential to illuminate the high-end shape of the density distribution, and to capture changes in the dense gas fraction and median volume density. Our results with and without a high-density power law tail differ appreciably; we highlight better knowledge of the probability density function (PDF) shape as an important area. We also show the implications of sub-beam density distributions for isotopologue studies targeting dense gas tracers. Differential excitation often implies a significant correction to the naive case. We provide tabulated versions of many of our results, which can be used to interpret changes in mm-wave line ratios in terms of adjustments to the underlying density distributions.

**Key words:** galaxies: ISM – galaxies: star formation – ISM: clouds – ISM: molecules – ISM: structure – radio lines: ISM

*Supporting material:* machine-readable tables

## 1. Introduction

Gas volume density, $\rho$, plays a central role in most theories explaining star formation in molecular clouds (e.g., see Gao & Solomon 2004b; McKee & Ostriker 2007; Lada et al. 2010; Padoan & Nordlund 2011; Federrath & Klessen 2012; Krumholz et al. 2012, among many others). Many models, especially those based on turbulent clouds (e.g., Padoan & Nordlund 2002; Krumholz & McKee 2005; Feldmann et al. 2011; Hennebelle & Chabrier 2011), consider the free fall time $\tau_{ff} \propto \rho^{-0.5}$ as the governing timescale for star formation, with an impact on the output cluster population (e.g., Kruijssen 2012). Another class of observationally motivated “threshold models” (Gao & Solomon 2004b; Wu et al. 2005, 2010; Heiderman et al. 2010; Lada et al. 2010, 2012; Evans et al. 2014), posit that stars form only in the densest parts of clouds, with the star formation rate driven by the gas mass above a gas volume density threshold.

This important theoretical role for gas density agrees with observations of local clouds. These show star formation confined to the highest density regions (Lada & Lada 2003; Heiderman et al. 2010; Lada et al. 2010; André et al. 2014),

and that the amount of high column density material in a cloud correlates with its ability to form stars (Kainulainen et al. 2009).

Testing these theories requires estimating gas density across galaxies. Such estimates also offer the prospect to understand how galaxies set the density of their gas. The nearby galaxy population offers access to a range of conditions and a clean external perspective not available in the Milky Way, and so offers many advantages as a laboratory to explore the role of gas density. The challenge to making such measurements is that the concentrations of dense gas seen in Milky Way clouds are very compact, ~0.1–1 pc (e.g., Lada & Lada 2003; André et al. 2014) and lie within larger, lower density superstructures. As a result, even within the ~40 pc beam of the highest resolution gas maps targeting nearby galaxies (e.g., Schinnerer et al. 2013), a vast range of gas volume densities are convolved together within an individual extragalactic beam.

Given this limitation, multi-line spectroscopy has become the standard approach to investigate the density of gas across galaxies (e.g., Gao & Solomon 2004b; García-Burillo et al. 2012; Kepley et al. 2014; Chen et al. 2015; Usero et al. 2015; Bigiel et al. 2016). By contrasting a line excited at high
gas densities, e.g., HCN (1-0), with a line excited at nearly all gas densities, e.g., CO (1-0), one can estimate the fraction of dense gas. In a pioneering study, Gao & Solomon (2004a) carried out a large survey of HCN in bright galaxies and used this approach to show that, to first order, the ratio of IR-to-HCN luminosity appears constant across the Galaxy population, while the ratio to IR-to-CO does not. Gao & Solomon and following authors (e.g., Wu et al. 2005, 2010; Lada et al. 2012) interpreted this as evidence for a universal gas density threshold for star formation.

More recent studies of normal star-forming galaxies by García-Burillo et al. (2012), Usero et al. (2015), and Bigiel et al. (2016) reveal a more complex relationship between high critical density lines, CO emission, and star formation. Usero et al. (2015) and Bigiel et al. (2016) demonstrated that both HCN (1-0)/CO (1-0), tracing the dense gas fraction, and IR/HCN (1-0), tracing the efficiency with which dense gas forms stars, depend on environment. These results have the same sense as those obtained contrasting local clouds with those in the Milky Way center (Longmore et al. 2013; Kruijssen et al. 2014).

So far, these extragalactic studies have mainly considered a two-phase molecular medium in which gas is either “dense,” and so emits HCN or HCO$^{+}$, or “not dense,” and so emits CO. This is a poor representation of the cold interstellar medium (ISM). Turbulent theories predict a wide range of densities, distributed approximately lognormally within any given cloud (e.g., Padoan & Nordlund 2002). Milky Way observations that resolve individual clouds show a large range of column densities, implying a large range of volume densities (e.g., Kainulainen et al. 2009; Abreu-Vicente et al. 2015). The functional form of the density distribution is still debated (e.g., Lombardi et al. 2015; Schneider et al. 2015), but all observations and models agree that a wide range of densities coexist within any individual region.

Further complicating the issue, emission from high density tracers like HCN (1-0) or HCO$^{+}$ (1-0) does not come exclusively from gas above some threshold collider density. The emissivity (line emission per mass) of a molecule does peak at some collider density, and that peak depends on the critical density, optical depth, and temperature of the molecule in question. But regions with lower densities can still emit; they merely do so with lower efficiency. If low density regions outnumber higher density regions, then they may contribute appreciably to, or even dominate, emission from that molecule. This ability of gas to emit below the nominal critical density (even after modification for line trapping) has been emphasized in the Galactic literature (e.g., Mangum & Shirley 2015; Shirley 2015, and references therein) but less explored in other galaxies.

In this paper, we consider the impact of these two effects on extragalactic observations that integrate over a wide range of distributions. We model line emission from regions that contain a wide distribution of densities, and account for a realistic dependence of emissivity on collider density. To do this, we treat a cloud as an ensemble of one-zone non-equilibrium models (calculated using RADEX; van der Tak et al. 2007) that share an optical depth, and so an escape probability. This approach follows Krumholz & Thompson (2007), who laid out this cloud model but focused on how density distributions influence observed star formation scaling relations (see also the closely related paper by Narayanan et al. 2008). Here, we

investigate the implications of sub-beam density distributions and realistic emissivities on observed millimeter line emission.

The paper proceeds as follows. We describe our model and calculations in Section 2. Then, in Section 3.1 we note several results from one-zone models that are important to interpret emission from density distributions. In Section 3.2 we show results integrated over realistic density distributions. We consider the median density producing emission, the pattern of line ratio changes induced by changing sub-grid volume density distributions, and the ability of line ratio variations to gauge changes in the dense gas mass fraction and median density. We also consider the implications of sub-beam density distributions for isotopologue ratios and dense gas conversion factors. We summarize our conclusions and discuss implications and future directions in Section 4.

### 1.1. Definitions of Densities

We discuss a number of densities throughout the paper. For clarity, we summarize these here. All densities refer to molecular hydrogen, H$_2$. We neither account for helium nor consider collisions with electrons or atomic hydrogen. Given that we focus on relative statements and cold, dark gas, these approximations should have minimal impact.

1. *Collider Density ($n_{H2}$):* The volume density of H$_2$ molecules. This is an input to our one zone models. When we refer to density distributions, we mean collider density distributions.

2. *Critical Density:* For a given transition and kinetic temperature, the collider density at which collisional de-excitation balance radiative de-excitation, taking no account of line trapping. See Shirley (2015) for a recent review.

3. *Effective Critical Density:* For a given transition and kinetic temperature, the collider density at which collisional de-excitation balance radiative de-excitation, taking into account radiative line trapping. See Shirley (2015) for a recent review.

4. *Most Effective Density for Emission ($n_{eff}$):* For a given transition, kinetic temperature, and optical depth, the minimum density at which the emissivity of the line (defined below) reaches 95% of its peak value. This is a new quantity defined in this paper. It is similar, but not identical, to the effective critical density. Often, effective critical density is used as a short-hand for the density at which gas is best at emitting. This quantity measures this directly. We present tabulated results $n_{eff}$ below.

5. *Median Density for Emission ($n_{med}^{emis}$):* For a given transition, kinetic temperature, optical depth, and density distribution, the collider density below which half of the line emission is generated. A closely related quantity is the main focus of Krumholz & Thompson (2007).

6. *Median Density by Mass ($n_{med}^{mass}$):* For a given density distribution, the collider density below which half the mass lies.

We also refer to the dense gas fraction, $f_{dense}$. For a given density distribution, this is the fraction of mass in the distribution that lies above a density threshold, $n_{thresh}$. We adopt $n_{thresh} = 10^{4.5} \text{ cm}^{-3}$ by default, but also consider a range of possible threshold values.
<table><caption>**Table 1**
Molecular Transitions Considered</caption><thead><tr><th>Transition</th><th>Frequency (GHz)</th><th>τ<sup>a</sup></th><th>n<sub>H2</sub> (ε<sub>max</sub>)<sup>b</sup> (cm<sup>-3</sup>)</th><th>X<sub>mol</sub><sup>c</sup></th></tr></thead><tbody><tr><td>HCN<sup>d</sup> J = 1 → 0</td><td>88.630</td><td>1</td><td>2 × 10<sup>5</sup></td><td>10<sup>-8</sup></td></tr><tr><td>HCO<sup>+</sup> J = 1 → 0</td><td>89.189</td><td>1</td><td>4 × 10<sup>4</sup></td><td>10<sup>-8</sup></td></tr><tr><td>HNC J = 1 → 0</td><td>90.663</td><td>1</td><td>1 × 10<sup>5</sup></td><td>10<sup>-8</sup></td></tr><tr><td>CS J = 2 → 1</td><td>97.980</td><td>1</td><td>7 × 10<sup>4</sup></td><td>10<sup>-8</sup></td></tr><tr><td><sup>13</sup>CO<sup>e</sup> J = 1 → 0</td><td>110.201</td><td>0.1</td><td>8 × 10<sup>2</sup></td><td>2 × 10<sup>-6</sup></td></tr><tr><td><sup>12</sup>CO J = 1 → 0</td><td>115.271</td><td>10</td><td>1 × 10<sup>2</sup></td><td>10<sup>-4</sup></td></tr><tr><td><sup>12</sup>CO J = 2 → 1</td><td>230.538</td><td><sup>*f</sup></td><td>1 × 10<sup>3</sup></td><td>10<sup>-4</sup></td></tr><tr><td><sup>12</sup>CO J = 3 → 2</td><td>345.796</td><td><sup>*f</sup></td><td>9 × 10<sup>3</sup></td><td>10<sup>-4</sup></td></tr></tbody></table>

**Notes.** Data for models for these species taken from Leiden Atomic Molecular Database (LAMBDA; Schöier et al. 2005).

<sup>a</sup> Representative optical depth assumed for this line. When not otherwise noted, this τ is used throughout the paper.

<sup>b</sup> Minimum collider density, $n_{H2}$, at which the emissivity reaches 95% of its peak value for T = 25 K and the representative τ.

<sup>c</sup> Fiducial abundance of the molecule adopted in this paper. This divides out of many aspects of the analysis.

<sup>d</sup> We ignore the splitting of the HCN 1-0 line because the component spacing is small compared to typical extragalactic line widths.

<sup>e</sup> Results for C<sup>18</sup>O are essentially equivalent to those for <sup>13</sup>CO if both lines are taken to be optically thin. We plot only <sup>13</sup>CO for clarity.

<sup>f</sup> For internal consistency, we calculate emission from CO (2-1) and CO (3-2) using the N/Δv calculated for CO (1-0). Thus, these lines are pinned to an assumed τ for CO (1-0).

## 2. Method

We consider a suite of commonly observed rotational transitions of molecules in the 80–115 GHz range. These lines, listed in Table 1, are bright enough that they can be surveyed in other galaxies. They are excited at a range of densities, so that ratios among them offer the prospect to constrain changing density distributions. Focusing on the 80–115 GHz range, we aim to de-emphasize excitation concerns due to temperature, highlighting the effect of changing density. We include the $J = 2 \to 1$ and $J = 3 \to 2$ transitions of CO, because these are observed across many nearby galaxies and often serve as bulk gas tracers, for practical reasons (e.g., Leroy et al. 2009; Wilson et al. 2012). This suite of lines closely resembles that mapped by the IRAM Large Program EMPIRE (Bigiel et al. 2016), the survey of Usero et al. (2015), and ALMA mapping by M. J. Gallagher et al. (2016, in preparation). A goal of this paper is to help inform interpretation of those and similar observations.

For each molecule, we create a grid of one-zone models (Section 2.1), each describing a single set of physical conditions. From each, we calculate the emissivity, $\epsilon$, which will be our figure of merit for much of the paper (Section 2.2). We combine these one-zone models to simulate emission from a realistic distribution of densities. To do this, we follow Krumholz & Thompson (2007) and make the simplifying assumption that a single escape probability, and so a single optical depth, describes each transition throughout the beam (Section 2.3). We combine this assumption with an assumed density distribution (Section 2.4) to calculate the beam-averaged emissivity, $\langle \epsilon \rangle$. Ratios of $\langle \epsilon \rangle$ for different lines but otherwise matched conditions can be observed as line ratios. Thus, this approach allows us to explore the impact of changing density distributions on observed line ratios. While useful and a large improvement over a one zone treatment, this

model has several shortcomings. We describe these and note directions for future improvement in Section 2.6.

### 2.1. Grid of One-zone Models

For each molecule, we use the RADEX code (van der Tak et al. 2007) to predict the emission from a series of one zone models. Each model assumes a single collider (H<sub>2</sub>) volume density, $n_{H2}$, kinetic temperature, $T_{kin}$, and column-per-line width, $N_{mol}/\Delta v$, of the molecule.<sup>15</sup> Given these inputs, RADEX solves for the level populations and predicted emission without assuming local thermodynamic equilibrium (LTE).

Our model grid covers a range of kinetic temperatures, $T_{kin} = 10–300$ K, and a range of H<sub>2</sub> collider densities, $n_{H2} = 10^1–10^8$ cm<sup>-3</sup>. This range of collider (H<sub>2</sub>) densities spans from well below to well above the critical density of each transition of interest. The range of column-per-line width, N/Δv, is $10^{13}–10^{19}$ cm<sup>-2</sup> (km s<sup>-1</sup>)<sup>-1</sup> for the CO molecules and $10^{10}–10^{16}$ cm<sup>-2</sup> (km s<sup>-1</sup>)<sup>-1</sup> for the other molecules. We chose these ranges to span from optically thin ($\tau \lesssim 0.1$) to optically thick ($\tau \gtrsim 10$) for each molecule.

We space grid points logarithmically along each axis of the grid, sampling 60 points along the $T_{kin}$ and $N_{mol}/\Delta v$ axes and 70 points along the $n_{H2}$ axis. At each ($n_{H2}$, $T_{kin}$, $N_{mol}/\Delta v$) point, we record the intensity, I, of each low J line in K km s<sup>-1</sup>, and its optical depth, τ. The result is a grid of $\tau(n_{H2}, T_{kin}, N_{mol}/\Delta v)$ and $I(n_{H2}, T_{kin}, N_{mol}/\Delta v)$ that serves as the backbone of our calculations.

All models assume the expanding sphere (classic large velocity gradient, LVG) geometry (Sobolev 1960; Goldreich & Kwan 1974). We use the version of RADEX released on 2011 November 30 and adopt molecular data from the Leiden Atomic Molecular Database (Schöier et al. 2005) for HCN, HNC, HCO<sup>+</sup>, CS, <sup>13</sup>CO, <sup>12</sup>CO. We also ran model grids for C<sup>18</sup>O. The results are very similar to those for <sup>13</sup>CO and we only present the <sup>13</sup>CO results here. To use our model for C<sup>18</sup>O, we suggest to divide the emissivity relative to H<sub>2</sub> by the <sup>13</sup>CO/C<sup>18</sup>O abundance ratio.

### 2.2. Emissivity

For each one zone model, we calculate the emissivity, $\epsilon$, of each transition. We define emissivity as the intensity, I, divided by the column density, N,

$$ \epsilon = \frac{I}{N}. \qquad (1) $$

Thus, the emissivity, $\epsilon$, measures how effectively gas emits in the specified transition for the conditions in the model. It has units of K km s<sup>-1</sup> (cm<sup>-2</sup>)<sup>-1</sup>. Along with a mean molecular mass, the emissivity can be recast to have units of K km s<sup>-1</sup> (M<sub>☉</sub> pc<sup>-2</sup>)<sup>-1</sup>. Therefore, $\epsilon$ is the mass-to-light ratio of gas with the given ($n_{H2}$, $T_{kin}$, $N_{mol}/\Delta v$).

In practice, we are interested in the emitted intensity per unit total gas (H<sub>2</sub>) rather than emitted intensity per unit mass of HCN or CO. Therefore, our calculations use the column of the

<footnote n="15"><sup>15</sup> We treat the column density (of the observed molecule) per unit line width, $N_{mol}/\Delta v$, rather than $N_{mol}$ and $\Delta v$ separately. In the escape probability treatment, it is this quantity that maps to optical depth rather than the two quantities separately.</footnote>
emitting molecule, $N_{mol}$. Thus, when we write $\epsilon$ below we mean

$$ \epsilon = \frac{I}{N_{H2}} = \frac{I}{N_{mol} X (mol)^{-1}}. \qquad (2) $$

Here, $X$ (mol) refers to the abundance of the molecule, with $X$ (mol) $\equiv N_{mol}/N_{H2}$. Our models yield $I/N_{mol}$. We note our fiducial values of $X$ (mol) for each molecule in Table 1. These are $10^{-8}$ for CS, HCN, HCO$^+$, HNC, $10^{-4}$ for CO, and $2 \times 10^{-6}$ for $^{13}$CO (i.e., 50 times lower than CO). These are motivated by existing abundance estimates (e.g., Martín et al. 2006), but for many of the results in this paper, the exact values of $X$ (mol) are not important. Instead, we will be interested in how the ability of molecular gas to emit in these lines varies as a function of changing density distributions, optical depth, and so on. As discussed in the next section, our model treatment implies modest changes in at least some values of $X$ (mol) within the model framework.

## 2.3. Distributions of Densities With Shared $\tau$ and $T_{kin}$

Observations of a whole cloud or a large part of a galaxy will blend emission from gas at many densities. To simulate this, we combine one-zone models, which each have a single associated $n_{H2}$, to model a distribution of densities. Our key simplifying assumption is that the emission for each line comes from gas that shares a single $T_{kin}$ and $\tau$. That is, we consider a distribution of densities, but take the gas to be isothermal and adopt a fixed optical depth within each beam.

The condition of fixed $\tau$, and thus, fixed escape probability, follows Krumholz & Thompson (2007), albeit our implementation differs in detail. For a given transition, $T_{kin}$, and $n_{H2}$, $\tau$ depends on $N_{mol}/\Delta v$. Although column densities at the scales relevant to radiative transfer can be difficult to gauge at extragalactic distances, $\tau$ can be accessed by observations, including via comparison of rarer isotopologues to the main gas tracers (e.g., Jiminez-Donaire et al. MNRAS submitted).

For the rest of the paper, we will not focus on the actual values of the column density, $N_{mol}$, or column-per-line width, $N_{mol}/\Delta v$. Instead, for each transition we will consider $\epsilon$ ($n_{H2}$, $T_{kin}$, $\tau$), where $\tau$ is the optical depth of that transition. We derive this quantity from our grid of one-zone models by selecting the appropriate $n_{H2}$ and $T_{kin}$ and then interpolating along the curve of $N_{mol}/\Delta v$ versus $\tau$ for our model grid to select $N_{mol}/\Delta v$ that gives the desired value of $\tau$. Then, we interpolate within our grid to calculate $\epsilon$ ($n_{H2}$, $T_{kin}$, $\tau$). Because of the fine spacing of the model grid, the interpolation has only a minor influence on the results.

Fixing $T_{kin}$ and $\tau$, we consider a distribution of volume densities, $P(n_{H2})$. We sum over all densities to derive the beam averaged emissivity:

$$ \langle \epsilon \rangle = \frac{\int n_{H2} P(n_{H2}) \epsilon (n_{H2}, T, \tau) dn_{H2}}{\int n_{H2} P(n_{H2}) dn_{H2}}. \qquad (3) $$

The beam-averaged emissivity, $\langle \epsilon \rangle$, resembles the one-zone emissivity (Equation (2)) in that it describes the intensity of the line per column of H$_2$. However, $\langle \epsilon \rangle$ now measures the effective emissivity of a whole distribution of densities convolved together.

We weight by $n_{H2}$ in Equation (3) because $P(n_{H2})$ represents the probability of finding a gas volume density $n_{H2}$ in a given volumetric cell within the cloud. The amount of emission from

gas at $n_{H2}$ will trace the total amount of mass at that density, rather than the amount of volume. The amount of mass at $n_{H2}$ is given by $n_{H2}P(n_{H2})$.

**Multiple Transitions of the Same Molecule:** Picking $\tau$ for one transition while fixing $n_{H2}$ and $T_{kin}$ implies $N_{mol}/\Delta v$. With $n_{H2}$, $T_{kin}$, and $N_{mol}/\Delta v$ specified, the emission from other transitions of the same molecule is also determined. Thus, once we specify $\tau$ for CO (1-0), the optical depths of CO (2-1) and CO (3-2) are no longer free parameters. In this paper, this issue arises mainly for CO. In general, our approach is to specify an observationally motivated $\tau$ for a single transition for each molecule. We then calculate emission from the other transitions self-consistently (see Rosolowsky et al. 2008, for a similar approach treating Galactic observations of ammonia). This means that the optical depth of these other transitions may not be fixed across the model distribution. In exchange, the calculated line ratios among multiple transitions of a single molecule will make physical sense.

**Implicit Changes in Relative Abundances:** A corollary of this point is that our model implies modest changes in the relative abundance of the molecules that we treat. For each molecule at each density, we pick $N_{mol}/\Delta v$ to yield our fiducial $\tau$ for the specified transition. Because we fix the distribution of $n_{H2}$, the ratio $N_{mol}/n_{H2}$ will change from density to density to keep $\tau$ fixed. As a result, the relative values of, e.g., $N_{HCN}$ and $N_{CO}$ will vary within a model distribution. This does not represent an inherent inconsistency, but does assume that mild abundance variations will conspire to yield fixed $\tau$. We discuss alternative formulations below (Section 2.6). The Krumholz & Thompson (2007) model suffers from the same internal inconsistency. They formulate this as a cloud radius that varies from transition to transition, so that in their model $N_{mol}/\Delta v$ can vary because $N_{mol} \approx n_{H2} X (mol)/R_{mol}$ varies with $R_{mol}$, where $R_{mol}$ is the radius of the emitting region. Because the radii of the HCN-emitting region, $R_{HCN}$, and the CO-emitting region, $R_{CO}$, can differ, the total abundance will also vary in this formulation.

## 2.4. Realistic Density Distributions

We explore emission from lognormal density distributions, sometimes modified to include a power law tail at high density. Lognormal distributions are expected, based on simulations (e.g., Vazquez-Semadeni 1994; Elmegreen & Scalo 2004; Mac Low & Klessen 2004), and we follow the description outlined for turbulent gas in Padoan & Nordlund (2002). Observations do support the existence of lognormal column density distributions within a cloud (e.g., Kainulainen et al. 2009; Rathborne et al. 2014; Abreu-Vicente et al. 2015). However, power law descriptions of the column density distribution may be equally valid given observational uncertainties (e.g., Lombardi et al. 2015). Even when a lognormal describes the bulk of the gas in a cloud, self-gravity may lead the densest gas to exhibit power-law density distributions (e.g., Kainulainen et al. 2009, 2013; Stutz & Kainulainen 2015).

Note that the appropriate form of the density distribution on scales larger than a cloud, including the $\gtrsim$kpc regions observed by many mm-wave spectral studies, remains a subject of investigation. On larger (60 pc) scales, molecular gas surface densities do appear to exhibit an approximately lognormal distribution (Leroy et al. 2016), although processes other than turbulence may create the lognormal surface density distributions observed at these large scales (Berkhuijsen &
Fletcher 2015). Simulations show density distributions that can depend strongly on environment and the nature and strength of stellar feedback (e.g., Hopkins et al. 2013). At the least, one might expect the superposition of a diverse population of lognormal and power law distributions. Here, we begin the simplest assumption and note further explanation as an area for future work.

**Lognormal Distributions:** In the Padoan & Nordlund (2002) formulation, the volume elements in a cloud show a distribution of gas volume densities described by:

$$dP(\ln n') \propto \exp\left(-\frac{(\ln n' - \overline{\ln n'})^2}{2\sigma^2}\right) d\ln n', \quad (4)$$

where $dP$ is the fraction of cells with volume densities in a logarithmic step $d\ln n'$ centered on $n'$, $n' = n_{H2}/n_0$ is the volume density normalized by the mean volume density, $n_0$, and $\sigma$ is the width of the distribution. The distribution does not peak at the mean volume density, $n' \equiv 1$, but at the mean of the logarithm, $\overline{\ln n'}$. The two are related via $\overline{\ln n'} = -0.5\sigma^2$ (see Padoan & Nordlund 2002, for details). In practice, we work in $\log_{10}$ rather than $\ln$, so the $\sigma$ that we quote are in dex and differ from those in the Padoan & Nordlund (2002) formalism by a factor of $\ln 10 \approx 2.3$.

Equation (4) describes the distribution of volume densities in the cloud. This differs from the distribution of mass because higher-density cells hold more mass. The fraction of mass in a given logarithmic density step is $dm \propto ndP$. Thus, as discussed by Padoan & Nordlund (2002), most mass resides in dense substructures that occupy a relatively small fraction of the volume.

In practice, we implement Equation (4) numerically, also creating a mass distribution from $dm = ndP$. Each lognormal distribution is specified by two numbers: the mean density, $n_0$, and the rms logarithmic width, bookkeeping in dex (i.e., $\log_{10} \sigma$). Following Padoan & Nordlund (2002), for a distribution resulting from isothermal supersonic turbulence, $\sigma$ is closely related to the three-dimensional turbulent Mach number, $M$. In that case,

$$\sigma \approx 0.43 \sqrt{\ln(1 + 0.25M^2)} \text{ dex.} \quad (5)$$

Typical Mach numbers in spiral and starburst galaxies span the range $M \approx 5-100$, so we will mainly be concerned with $\sigma \approx 0.6-1.2$ dex (e.g., see Leroy et al. 2016). Unless otherwise noted, we adopt a fiducial $\sigma = 0.8$ dex, corresponding to $M \sim 10$.

**Power Law Tails:** Observations and simulations suggest that self-gravity in the densest parts of a cloud leads to the formation of a power law tail in $P(n)$ at high $n_{H2}$. This effect is observed to vary from cloud-to-cloud, with the amount of mass in the tail correlating with the star-forming activity (Kainulainen et al. 2009; Abreu-Vicente et al. 2015).

We explore this effect by incorporating a power law tail into some of our density distributions. In these cases, we follow Federrath & Klessen (2013) and take $P(n_{H2})$ to be a power law of form

$$dP(\ln n') \propto \exp(\alpha \ln n') \text{ where } n' > n'_{\text{thresh}}. \quad (6)$$

Here, $n'_{thresh}$ is the threshold for the onset of the power law tail in units of the density normalized by the mean, $n_0$, and $\alpha$ is the

<figure><img src="image_2.png" /><figcaption>Figure 1. Illustration of model density distributions by volume (blue) and mass (green), with (dashed) and without (solid) power law tails at high densities. Vertical lines indicate the mean density (gray) and the adopted threshold for the onset of a power law tail (pink), when one is present. Although a large part of the volume holds low-density gas, most of the mass resides in a small part of the volume with high densities (see Padoan & Nordlund 2002). Similarly, although the power law tail adds only a small amount of additional volume at high densities, it contains ~20% of the mass for this model.</figcaption></figure>

slope of the power law in logarithmic units. Note that this is offset by one power from the slope in non-logarithmic units, so that $dP(n) \propto n^{\alpha-1}$ (see Federrath & Klessen 2013). As a result, for any $\alpha < -1$ the mass distribution ($\propto \int n \times n^{\alpha-1} dn$) will converge. In these models, a lognormal distribution still describes gas at $n' < n'_{thresh}$.

Following Federrath & Klessen (2013) and Vallini et al. (2016), we take $\ln n'_{thresh} \approx 3.8$. That is, we set the power law to begin at ~45 times the mean density of the cloud. We take $\alpha = -1.5$, intermediate in the range of -1 to -2.5 found by Federrath & Klessen (2013), albeit toward the steep end of the values observed for column density distribution in Orion (Stutz & Kainulainen 2015). In principle, the appropriate threshold for the onset of self-gravity depends on the Mach Number and other quantities (see Federrath & Klessen 2013), with our adopted value appropriate for $M \approx 7$ (Krumholz & McKee 2005; Padoan & Nordlund 2011). Parameter studies of $n'_{thresh}$ and $\alpha$ will be useful, but lie beyond the scope of this paper.

Figure 1 illustrates our model distributions of volume densities for an example case with mean density $n_0$, $H_2 = 10^3$ cm$^{-3}$ and width $\sigma = 0.8$ dex. We show both the amount of volume and the amount of mass at each density, and illustrate the distribution with and without a power law tail. In this case, the power law tail integrated from $n'_{thresh}$ through $n_{H2} = 10^7$ cm$^{-3}$ holds ~20% of the total mass. The same density range in the pure lognormal case holds only ~7% of the mass. Thus, despite the modest fraction of volume in the power law tail, it has an appreciable impact on the mass distribution.
## 2.5. Emission From Distributions of Densities

We create density distributions for a range of $\sigma$ and $n_0$ appropriate to normal spiral galaxies and modest starbursts. In each case, we make two versions of the distribution: one with a power law tail and one without. From these density distributions, we calculate the corresponding mass distributions. We note the volume density below which 50% of the mass lies, and the fraction of the mass in the power law tail, if one is implemented. Finally, we calculate the fraction of the mass that lies at densities above $n_{H2} = 10^{4.5} \text{ cm}^{-3}$, a value commonly conflated with “dense” gas.

We then combine our mass distributions with our grid of one zone models, assuming a fixed $T_{\text{kin}}$ and $\tau$. For each model and each line in Table 1, we calculate $\langle \epsilon \rangle$, as defined in Equation (3). This is the emission per unit mass of the whole observed region. We also record the volume density below which 50% of the light is emitted.

The tables present results from these calculations, which serve as the basis for the plots and discussion in the rest of this paper.

## 2.6. Limitations of the Model and Directions for Improvement

Our major simplifying assumption is that all zones within a beam share the same $\tau$. That is, our model has a vastly simplified geometry. As mentioned above, fixing $\tau$ imposes modest variations in the relative abundances of different species. These also exist in the Krumholz & Thompson (2007) model, manifesting as different cloud radii for each molecule.

An alternative approach would be to specify $\tau$ for a single transition and a single species (e.g., CO) and then adopt a fixed abundance pattern to derive $N_{\text{mol}}/\Delta v$ for all other species, allowing their optical depth to vary. We consider assuming fixed $\tau$ to impose moderate optical depth on the key transitions (something apparently required by observations, e.g., Meier et al. 2015, Jimenez-Donaire et al. 2017) more reliable than adopting an uncertain abundance pattern.

One can imagine alternative configurations, with coupling between adjacent zones or built-in correlations, e.g., between $\tau$ and $n_{H2}$. A natural alternative to fixed $\tau$ would be to impose some dynamical state on the individual zones and so link $n_{H2}$ to $N/\Delta v$. Such models have been frequently invoked in the context of the CO-to-H2 conversion factor (e.g., Maloney & Black 1988; Bolatto et al. 2013) and used to treat ensembles of clouds (e.g., Aalto et al. 2015). The validity of imposing such a constraint on each small range of gas densities within a cloud is unclear, but this represents a direction for exploration. In particular, assuming a dynamical state for the power law tail, when present, represents a natural addition to the model.

We assume that a single, constant temperature, $T_{\text{kin}}$, describes our clouds, and consider a restricted range, $T_{\text{kin}} \approx 15-35 \text{ K}$ for most of the paper, taking $T_{\text{kin}} = 25 \text{ K}$ by default. The model grid includes a larger range of temperatures. Although these are not the focus of the paper, accounting for the effects of $T_{\text{kin}}$ variations will be key to explore with the contrast between central molecular zones and galaxy disks. This contrast between these two common environments is dramatic and readily observed in both the Milky Way (e.g., Longmore et al. 2013; Kruijssen et al. 2014) and other galaxies (e.g., Usero et al. 2015). Density is thought to play a key role in this contrast (e.g., Krumholz & Kruijssen 2015; Krumholz

et al. 2016), but temperature variations will also need to be accounted for. More, resolved Galactic clouds do show significant temperature substructure, and a density–temperature or optical depth-temperature correlation could be added to the model in place of fixed $T_{\text{kin}}$. However, at present we have no general prescription to implement.

Finally, beyond the variations needed to fix $\tau$, we assume fixed abundances, $X(\text{mol})$, throughout our models. One could modify Equations (2) and (3) to account for varying abundances. These could change within the distribution as a function of density or appear as overall changes in the normalization of $\langle \epsilon \rangle$ between measurements. Implementing such chemical variations is beyond the scope of this paper, but even simple PDR models show significant chemical substructure within clouds, so this is a productive direction for future investigation. Here, we begin using the simplest approach, and trust that the strongest abundance variations occur at the lowest and highest densities, whereas most of the mass, and the focus of this study, rests at intermediate scales.

To address these concerns, comparison to PDR models (e.g., Le Bourlot et al. 2012), numerical simulations that can capture complex cloud geometry (e.g., Glover & Clark 2012), and detailed observations of Galactic clouds offer a clear way forward. The DESPOTIC code by Krumholz (2014) offers the prospect to make several next steps; it allows multi-zone clouds and implements some basic PDR structure.

For the remainder of the paper, we emphasize that insight can be gained from our simple approach, which already offers a large improvement over one-zone models.

## 3. Results

Our models predict how mm-wave line emission depends on changes in the sub-beam density distribution, and so give a framework to explore density variations across galaxies. Before discussing these results, we show the dependence of emissivity, $\epsilon$, on $n_{H2}$ and $\tau$ in our grid of one-zone models (Section 3.1). We then consider emission from distributions of densities (Section 3.2). We examine the median density from which emission emerges (Section 3.3), the emergent line ratio pattern (Section 3.4), and the ability to estimate changes in the density distribution based on line ratio variations (Section 3.5). We also use our models to explore the generality of dense gas conversion factors (Section 3.6) and complications in the use of optically thin isotopologues to gauge optical depths (Section 3.7).

### 3.1. Emissivity, Density, and Optical Depth in One-zone Models

The emissivity, $\epsilon$, varies across our grid of one-zone models as a function of $n_{H2}$, $T_{\text{kin}}$, and $\tau$. The variation of $\epsilon$ with $n_{H2}$ will be key to understand our results for density distributions. Figure 2 shows $\epsilon(n_{H2})$ for fixed $T_{\text{kin}} = 25 \text{ K}$ and the $\tau$ in Table 1. The absolute emissivity of the molecules varies from line to line. To highlight the dependence on $n_{H2}$, we normalize results for each line by the maximum $\epsilon$ for any density in the calculation.

Figure 2 illustrates the sensitivity of our line suite to a wide range of collider densities. It also shows that, although $\epsilon$ does peak at some $n_{H2}$ value, it can be substantial even at much lower densities. In this calculation, $\epsilon$ remains $\gtrsim 0.1$ times its peak value even at $n_{H2}$ ten times lower than the density for peak
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Visualization of our one-zone model grid. Normalized emissivity as a function of collider volume density, $n_{H2}$, for our considered transitions at their nominal optical depths (Table 1; the curves in the figure match the order of the legend from left to right). Here, normalized emissivity is defined as emission per unit column density divided by the maximum for that line at any density. The suite of lines that we consider is excited at a range of volume densities, so that observing all of these lines in a beam gives sensitivity to the distribution of volume densities within the beam. The dependence of emissivity on volume density is not a perfect step function. Instead, material at densities below the nominal effective density still emits, only somewhat less efficiently.</figcaption></figure>

emissivity. That is, molecules still emit effectively well below the density that maximizes their emissivity. This point has been recognized and widely emphasized in the Galactic literature (see Evans 1999; Shirley 2015). In extragalactic observations, where observations necessarily average over a wide range of densities in a single beam, it is potentially even more important than in the Milky Way.

For the optically thick lines, the drop in $\epsilon$ at high $n_{H2}$ results from our assumption of fixed $\tau$. The optical depth for a given transition depends on the column density of molecules in the lower state for that transition ($J=0$ for most of our molecules, $J=1$ for CS). At densities near the peak, the molecules do not obey LTE. Sub-thermal excitation leads to a higher fraction of molecules in the lower, absorbing state. At higher densities, as LTE sets in, the column needed to reach our adopted $\tau$ will be higher because more molecules will be in the higher $J$ states. If we remake Figure 2, fixing $N_{mol}/\Delta v$ instead of $\tau$, then we do not see a similar drop in $\epsilon$ at high $n_{H2}$ for optically thick lines. For optically thin lines, $\epsilon$ will still peak near the critical density for fixed $N_{mol}/\Delta v$, because at higher $n_{H2}$, near LTE, some additional potentially emitting molecules will be excited out of the emitting state.

In addition to collider density, optical depth sets $\epsilon$ at fixed $T_{kin}$. At high $\tau$, emitted photons are re-absorbed rather than escaping, allowing more time for collisional de-excitation. This line trapping (Goldreich & Kwan 1974; Scoville & Solomon 1974) lowers the effective critical density of a line, with the well-known result that, at high $\tau$, the effective critical density is adjusted by a factor $\tau^{-1}$ relative to the critical density in the absence of line trapping (for recent treatments, see Scoville et al. 2015; Shirley 2015).

Thus, high $\tau$ lowers the $n_{H2}$ for peak $\epsilon$ at a given $T_{kin}$. The left panel of Figure 3 illustrates the effect of changing $\tau$ on the normalized emissivity versus $n_{H2}$ plot seen in Figure 2. The right panel shows how changing $\tau$ affects the peak of the $\epsilon$ versus density curve over the range where $\tau \approx 1$. For convenience, we include a machine-readable table (Table 2) reporting the collider density, $n_{H2}$, at which the emissivity reaches 95% of its peak value given some $T_{kin}$ and $\tau$.

At fixed $T_{kin}$, the absolute value of $\epsilon$ depends on both $n_{H2}$ and $\tau$. Figure 4 shows the interplay of these two factors for CO and

HCN. We plot $\epsilon$, normalized to the peak of the whole plane, over a wide range of $\tau$ and $n_{H2}$. The figure shows the decrease in $n_{H2}$ for maximum $\epsilon$ with increasing $\tau$ (the peaks move left as one moves up). It also illustrates that at or around the critical density, $\epsilon$ becomes higher for lower $\tau$. As one would expect, the most emissive gas has low optical depth and $n_{H2}$ near the critical density of the transition.

## 3.2. Emission for Density Distributions

We combine our one-zone calculations to simulate emission from a distribution of densities within the beam. As described in Section 2.3, we calculate the distribution of mass as a function of volume density, $n_{H2}P(n_{H2})$ for realistic distributions of densities. Then, we multiply the mass distribution by the one-zone emissivity shown in the previous section. From this product, we calculate the beam-averaged emissivity, $\langle\epsilon\rangle$, and the density below which 50% of the light is emitted for each spectral line, $n_{med}^{emis}$ (see Section 1.1).

Figure 5 illustrates the interplay of emissivity and density distribution at the core of the model. The red curve shows $\epsilon(n_{H2})$ of HCN for $T_{kin} = 25$ K and $\tau = 1$. The blue curves show a series of lognormal distributions with varying mean density, $n_0$, and width, $\sigma$. The product of the emissivity and the mass distribution, shown as a purple curve, indicates how much emission emerges from each density. The plots take $\epsilon$ and the mass distribution to be normalized at their peak values, but we do not normalize the emission, so that the purple curves can be compared among panels to see how the integrated brightness of the line varies.

## 3.3. What Densities Produce Emission?

Figure 5 illustrates how the density of gas producing HCN emission depends on the sub-beam density distribution. Perhaps most strikingly, when the gas at densities far below the critical density significantly outmasses gas at high density (as in the top left panels), most HCN emission comes from modest-density gas. When much of the mass lies at high $n_{H2}$, as in the bottom row, the density for emission tracks the density of the gas and HCN traces the total gas reservoir, as CO does for
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Visualization of our one-zone model grid. (left) Normalized emissivity (y-axis) as a function of collider volume density (x-axis) for ¹²CO (1-0) considering a range of optical depths, τ (color). Lines with lower optical depths (blue) require higher collider densities to reach their maximum emissivity. The result is that thin CO transitions, like ¹³CO (1-0) or C¹⁸O (1-0), emit maximally well at higher densities than thick transitions. This is the well-known effect of “line trapping” lowering the effective critical density of a transition (e.g., see Goldreich & Kwan 1974; Shirley 2015). (right) Most effective density for emission from one zone non-LTE RADEX models as a function of assumed optical depth for our considered transitions over the range τ = 0.1 to τ = 1 (the lines proceed in order of the legend from bottom to top of the plot). At the low end, this approaches the classical critical density. At the high end, the standard factor of τ⁻¹ to account for line trapping applies.</figcaption></figure>

<table><caption>**Table 2**
Density for Maximum Emissivity in One-zone Models</caption><thead><tr><th>Line</th><th>T<sub>kin</sub><br>(K)</th><th>log<sub>10</sub> τ</th><th>log<sub>10</sub> n<sub>eff</sub><sup>0.95</sup><br>(cm<sup>-3</sup>)</th></tr></thead><tbody><tr><td>12CO10</td><td>10.0</td><td>-1.00</td><td>3.38</td></tr><tr><td>12CO10</td><td>10.0</td><td>-0.90</td><td>3.37</td></tr><tr><td>12CO10</td><td>10.0</td><td>-0.80</td><td>3.36</td></tr><tr><td>12CO10</td><td>10.0</td><td>-0.70</td><td>3.35</td></tr><tr><td>12CO10</td><td>10.0</td><td>-0.60</td><td>3.33</td></tr><tr><td>12CO10</td><td>10.0</td><td>-0.50</td><td>3.30</td></tr><tr><td>12CO10</td><td>10.0</td><td>-0.40</td><td>3.28</td></tr><tr><td>12CO10</td><td>10.0</td><td>-0.30</td><td>3.25</td></tr><tr><td>12CO10</td><td>10.0</td><td>-0.20</td><td>3.21</td></tr><tr><td>12CO10</td><td>10.0</td><td>-0.10</td><td>3.17</td></tr></tbody></table>

**Note.** The table reports the minimum collider density at which the emissivity, ϵ, reaches 95% of its peak value for a given T<sub>kin</sub> and τ.

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

lower density clouds (see Krumholz & Thompson 2007; Narayanan et al. 2008).

We identify the characteristic density for emission, $n_{\text{med}}^{\text{emis}}$, as the density below which half the emission emerges. For HCN, we plot $n_{\text{med}}^{\text{emis}}$ as a function of $n_0$, $\sigma$, and $\tau$ in Figures 6 and 7. Figure 6 shows the effect of varying the mean density, $n_0$, and distribution width, $\sigma$ on $n_{\text{med}}^{\text{emis}}$ for models with and without a power law tail. The range of $\sigma$ spans from 0.4–1.2, corresponding to $M \approx 2.5–100$. Although we illustrate the case for HCN, other dense gas tracers show similar behaviors with shifted density scales.

For a lognormal distribution (left panel in Figure 6), the median density for HCN emission depends on both the mean,

$n_0$, and the width, $\sigma$. A high $n_0$ shifts the distribution of masses to higher $n_{\text{H2}}$. A high $\sigma$ does the same, both by shifting the peak of the mass distribution and by extending the wing of the distribution to higher densities. Both effects yield higher $n_{\text{med}}^{\text{emis}}$ (see Section 2.4 or compare across a row in Figure 5), and for pure lognormal distributions, $\sigma$ and $n_0$ exert similar effects. This interplay leads to the curved contours in the left panel of Figure 6. To achieve a given $n_{\text{med}}^{\text{emis}}$, a distribution can either have high $\sigma$ or high $n_0$ or both.

The absolute value of $n_{\text{med}}^{\text{emis}}$ varies over a wide range in the left panel of Figure 6. Lognormal distributions with low $n_0$ and low $\sigma$ will generate most of their HCN emission from low $n_{\text{med}}^{\text{emis}} \lesssim 10^3–10^4$ cm⁻³. Such distributions do not have high beam-averaged emissivity (the purple curve is low in Figure 5), but they could be common. Cloud-averaged mean densities $\lesssim 10^{2.5}$ cm⁻³ and Mach numbers $M \lesssim 5–10$ are widespread in nearby galaxies (e.g., Bolatto et al. 2008; Hughes et al. 2013b; Leroy et al. 2016). If the sub-beam density distributions in such clouds were pure lognormals, then much of their HCN emission could arise from gas at moderate densities, in contrast to the normal extragalactic view of HCN as a high-density tracer.

The right panel of Figure 6 shows that the inclusion of a power tail in the density distribution (see Section 2.4) changes this picture somewhat. Such a tail can add high $n_{\text{H2}}$, high $\epsilon$ material to clouds with otherwise low mean densities. When we include such a tail in the calculation, the range of $n_{\text{med}}^{\text{emis}}$ narrows considerably across the parameter space in Figure 6. As long as $\sigma \gtrsim 0.6 \text{ dex}$ ($M \gtrsim 5$), the median density for emission is in the range $n_{\text{H2}} \gtrsim 10^{4.5}$ cm⁻³. For narrower density distributions ($\sigma \lesssim 0.6 \text{ dex}$) with modest $n_0$, our formulation of the power law tail does not add as much high-density material and $n_{\text{med}}^{\text{emis}}$ still reaches low values. This may only reflect an issue with our
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Visualization of our one-zone model grid. Emissivity, $\epsilon$, of CO (left) and HCN (right) at fixed $T_{\text{kin}}$ as a function of optical depth, $\tau$, and collider density, $n_{\text{H2}}$. Contours show 10%, 20%, 30%, ..., 90% of the peak $\epsilon$ over the whole plane. $T_{\text{kin}} = 25$ K for all models. The plots show decrease of $n_{\text{H2}}$ for maximum $\epsilon$ with increasing $\tau$ and the decrease in $\epsilon$ with increasing $\tau$ at fixed $n_{\text{H2}}$.</figcaption></figure>

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

adopted threshold for the onset of the power law; we fix this at a multiple of $n_0$, whereas a more general physical prescription might reflect $\sigma$ and other physical properties of the gas (e.g., see Krumholz & McKee 2005; Federrath & Klessen 2012).

The inclusion of a power law tail diminishes the influence of $\sigma$ on $n_{\text{med}}^{\text{emis}}$ compared to the pure lognormal case. Although this results partially from our model formulation, the physical point is general: in the case where self-gravitating structures dominate the high end of the density distribution, the turbulent Mach number may have diminished effect on emission from high-density tracers. This can be tested, for example, by correlation of density sensitive ratios like HCN/CO or CS/CO with measurements of the turbulent velocity dispersion, e.g., from high-resolution CO imaging.

Figure 6 adopts a fixed optical depth, $\tau = 1$, for HCN. Variations in $\tau$ will also affect $n_{\text{med}}^{\text{emis}}$. In Figure 7, we plot $n_{\text{med}}^{\text{emis}}$ for HCN, holding $\sigma$ fixed at 0.8 dex, but allowing $\tau$ to vary from 0.1 to 10. Again, the left panel shows results for a pure lognormal, and the right shows results for a distribution with a power law tail at high densities. The overall sense of both panels is that a high $\tau$ leads to a lower $n_{\text{med}}^{\text{emis}}$, but the effects are modest over the plausible range of $\tau$ for HCN.

For the lognormal case, $\tau$ exerts a weaker influence than variations in $\sigma$. Variations across the plausible range of $\tau$ (optically thin up to $\tau \sim 10$) change $n_{\text{med}}^{\text{emis}}$ by only a factor of a few at fixed $n_0$. Optical depth has a stronger influence when a power law tail is present and the mean density, $n_0$, is low. In this case, a higher $\tau$ shifts the density for maximum emission (the effective critical density) to lower $n_{\text{emis}}$. This allows lower-density material on the tail to emit effectively. The result is that $\tau$ plays a main role in setting $n_{\text{med}}^{\text{emis}}$ for modest $n_0 \lesssim 10^3 \text{ cm}^{-3}$ when a power law tail is present.

The plots here illustrate $n_{\text{med}}^{\text{emis}}$ for HCN. Similar trends hold for our other lines. In general, the interplay of a realistic density distribution and $\epsilon$ ($n_{\text{H2}}$) can lead to a wide range of emitting densities for lognormal distributions. Including a power law

tail, however, mitigates the variation for many cases. Optical depth exerts a more modest influence on $n_{\text{med}}^{\text{emis}}$ than the density distribution in the lognormal case, but can play a significant role when a power law tail is present and $n_0$ is modest.

We include tabulated results for $n_{\text{med}}^{\text{emis}}$ as a function of $\tau$, $n_0$, $\sigma$, and $T_{\text{kin}}$ for each line in Table 3. These can be used to reproduce the results of this and the following sections.

### 3.4. Beam-averaged Emissivity and Line Ratio Patterns

The convolution of emissivity and the density distribution determines the overall ability of the gas to emit. For each model, we record the beam-averaged emissivity, $\langle\epsilon\rangle$, defined in Equation (3). Here, $\langle\epsilon\rangle$ measures emission per unit mass, averaged over the whole density distribution, and so quantifies how well a given density distribution emits in a given transition. In Figure 5, $\langle\epsilon\rangle$ corresponds to the integral of the purple curves, and so varies across the figure.

Figure 8 shows $\langle\epsilon\rangle$ as a function of $n_0$ and $\sigma$ for HCN (1-0) and $^{12}$CO (1-0) at their fiducial $\tau$. For HCN (1-0), shown in the left panel, $\langle\epsilon\rangle$ changes by more than two orders of magnitude across the plotted parameter space. The sense of the variation follows Figure 5: high $n_0$ and high $\sigma$ lead to higher $\langle\epsilon\rangle$ because they yield more dense gas in the distribution. The right panel shows $\langle\epsilon\rangle$ for CO (1-0) at $\tau = 10$. In contrast to HCN (1-0), CO (1-0) shows little variation in $\langle\epsilon\rangle$, reflecting that the low densities needed to excite emission are present across the whole model grid. These weak variations in $\langle\epsilon\rangle$ help motivate the use of low J CO as a tracer of the total H$_2$ column.

The top left panel in Figure 9 shows how $\langle\epsilon\rangle$ changes with mean density, $n_0$, for each of our target lines. This plot assumes a pure lognormal distribution with $\sigma = 0.8$ dex ($M \sim 11$). The density-sensitive transitions (HCN (1-0), HNC (1-0), HCO$^+$ (1-0), CS (2-1)), show strong variations in $\langle\epsilon\rangle$ as $n_0$ changes. These transitions are much brighter in denser gas. As a result, $\langle\epsilon\rangle$ increases by more than an order of magnitude as $n_0$
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Illustration of our model for HCN emission. The purple curves show emission as a function of collider density, $n_{\text{H2}}$, for lognormal distributions with a range of mean densities, $n_0$ (in cm$^{-3}$), and widths ($\sigma$ in dex). The blue curves show the distribution of mass, the red curve—which is the same in all panels—shows $\epsilon(n_{\text{H2}})$ for HCN at $T_{\text{kin}} = 25$ K, $\tau = 1$. Their product, the purple curve, shows the light emitted as a function of density. Even as $\epsilon(n_{\text{H2}})$ remains fixed, the density of gas emitting HCN (1-0) changes dramatically with the sub-beam density distribution.</figcaption></figure>

goes from $\sim 10^2$ to $10^5$ cm$^{-3}$. Meanwhile, the CO transitions show weak variations in $\langle \epsilon \rangle$ across the same range of $n_0$.

By observing ratios between the dense gas tracers and the CO lines, one can leverage these differences in $\langle \epsilon \rangle$ to constrain the sub-beam density distribution. The top right panel of Figure 9 shows $\langle \epsilon \rangle$ for each line divided by $\langle \epsilon \rangle$ for CO (1-0). These $\langle \epsilon \rangle$ ratios can be directly observed as line ratios within a matched beam. That is, the top right panel of Figure 9 shows predicted line ratio patterns for lognormal density distributions with a range of mean densities, $n_0$.

The pattern in Figure 9 changes with $n_0$. Line ratios between high-density tracers and CO (1-0) change by more than an order of magnitude as $n_0$ varies. Variations among transitions that trace similar densities are much smaller. For example, ratios among the CO lines change by a factor of $\lesssim 2$ as $n_0$ changes more than two orders of magnitude.

The exact ratios in Figure 9 depend on our adopted abundances, which are often poorly known. Even if we do

not know the true abundances, they may remain approximately fixed among targets. In this case, comparing line ratios offers a way to diagnose changes in physical conditions. That is, translating an HCN/CO ratio to a dense gas fraction requires knowledge of the abundance of HCN and CO. However, comparing two HCN/CO ratios to one another may allow a comparative statement about the mean density or the dense gas fraction in which factors like abundance divide out.

The bottom panels in Figure 9 illustrate this approach for a pure lognormal (left) and a lognormal distribution with a power law tail (right). Again, we plot ratios relative to CO (1-0). Now, however, we normalize each line ratio pattern by the pattern calculated for $n_0 = 10^3$ cm$^{-3}$. Thus, the figure shows how we expect line ratios to change as $n_0$ changes. For example, if $n_0$ changes from $10^3$ cm$^{-3}$ to $10^4$ cm$^{-3}$, then we expect HCN/CO to increase by 0.5 dex, about a factor of 3. These changes are independent of the absolute value of the adopted abundances,
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Median density for emission of HCN (1-0), $n_{\text{med}}^{\text{emis}}$, as a function of the mean density, $n_0$, and width, $\sigma$, for lognormal density distributions without (left) and with (right) a power law tail (Section 2.4). Contours are spaced by 0.1 dex. The range of $\sigma$ corresponds to turbulent Mach numbers $M \approx 2.5-100$. For pure lognormal distributions, the median density for HCN emission varies strongly with the combination of $n_0$ and $\sigma$. It can reach low values, $n_{\text{H2}} \lesssim 10^3-10^4 \text{ cm}^{-3}$, for values of $n_0$ and $\sigma$ observed for extragalactic molecular clouds. A power law tail, if present, changes the picture (right panel). Such a tail ensures the presence of some high density, high emissivity material even for low $n_0$. In this case, $n_{\text{H2}}^{\text{emis}}$ remains $\gtrsim 10^{4.5} \text{ cm}^{-3}$ except for very narrow, low $n_0$ distributions. Because the power law tail now sets the high density mass distribution, the influence of $\sigma$ also decreases when such a tail is present.</figcaption></figure>

they require only that the abundance remains fixed between locations.

In the next section, we explore the quantitative relationship between changing line ratios and changes in the dense gas fraction or median density. Before proceeding, we note some qualitative behavior from Figure 9. First, large changes in $n_0$ induce large changes in the ratios between dense gas tracers and total gas tracers. This is the strongest behavior in the figures. Internal ratios among CO lines or dense gas tracers show relatively weak contrasts, but a change in $n_0$ leads to significant changes between any dense gas tracer and CO.

Second, the magnitude of this variation depends in detail on the sub-grid density distribution adopted. That is, the pure lognormal distributions (left) and those that include a power law tail at high densities (right) show distinct results. The model with a power law tail shows weaker variations in line ratios across the same range of $n_0$. It also shows less variation among the different high-density tracers, especially at low mean density (the bluer colors). For the lower end of $n_0$, $\sim 10^2-10^3 \text{ cm}^{-3}$, the four high-density tracers all show about the same strength of variations despite their different effective critical densities.

Meanwhile, the pure lognormal shows strong variations in line ratios that “flare” as a function of effective critical density. That is, variations in line ratio as a function of $n_0$ are strongest for our highest critical density line, HCN. They appear weakest for $\text{HCO}^+$, which has the lowest effective critical density among our dense gas tracers. This flaring pattern is also evident in the lognormal with a power law, but the behavior is much weaker and appears more prominently at high $n_0$ than at low $n_0$, where it is almost absent.

The differences between the pure lognormal and the case with a power law tail can be understood from the shapes of the two distributions. To first order, each tracer picks up material at or above some effective critical density. In the power law case,

as long as two densities both lie on the power law tail, the ratio of material above these two densities will be fixed. That is:

$$ \frac{M(n > n_1)}{M(n > n_2)} = \frac{\int_{n_1}^{\infty} n^{\alpha} dn}{\int_{n_2}^{\infty} n^{\alpha} dn} \propto \frac{n_1^{1+\alpha}}{n_2^{1+\alpha}}, \quad (7) $$

where we have used the definitions above<fnref n="16" /> and take $\alpha < -1$. Equation (7) illustrates that once $\alpha$ and two threshold densities on the power law tail are determined, then the ratio of mass above these two densities on the power law tail will be fixed. In other words, the power law itself is scale-free: $n_0$ only affects the lognormal part and the interface between the two.

No such fixed ratio holds for the lognormal. The integral above some threshold density falls rapidly with increasing threshold ($\propto 0.5(1 - \text{erf} n)$). As a result, for any pair of densities, the ratio will change as the mean of the distribution changes. The sense of the change is that the mass above a threshold density close to (but still above) the mean density changes less than the mass above a threshold density at a high value. Thus, for the lognormal distribution, the lines with highest effective critical densities (HCN, HNC) show the strongest variations in their ratios. Those with lower effective critical densities show weaker variations. This leads to the flaring seen in Figure 9. It will also lead to a nonlinear mapping between any given line ratio and dense gas fraction, assuming that dense gas fraction is defined relative to some fixed threshold density.

Temperature variations create weaker changes in the line ratio pattern than density variations. Figure 10 shows the effect of varying $T_{\text{kin}}$ while other inputs to the model remain fixed. Over the range $T_{\text{kin}} = 10-50 \text{ K}$, line ratio variations induced by

<footnote n="16">¹⁶ Recall that $\alpha$ is the slope in logarithmic units, so that the integral in normal space goes as $\int n \times n^{\alpha-1} \propto \int n^{\alpha}$.</footnote>
<figure><img src="image_10.png" /><figcaption>Figure 7. As in Figure 6, but now plotting optical depth, $\tau$, on the y-axis while holding the distribution width fixed at $\sigma = 0.8$ dex. Optical depth plays a modest role for the pure lognormal case (left panel) at low $n_0$. It has a stronger impact on $n_{\text{med}}^{\text{emis}}$ when the density distribution includes a power law tail (right panel). In this case, a high optical depth allows gas with densities on the tail to emit effectively at lower $n_{\text{H2}}$ and so lowers $n_{\text{med}}^{\text{emis}}$.</figcaption></figure>

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

changing temperature show lower magnitude and less flaring than those that we calculate for changing $n_0$. This reflects that a higher $T_{\text{kin}}$ also renders the CO line brighter. A higher $T_{\text{kin}}$ does have some effect on the line ratio pattern, however, because high $T_{\text{kin}}$ lowers the critical density and thus makes it easier to excite emission from dense gas tracer at intermediate densities.

*Implication for Observations:* Figures 9 and 10 show how one can approach observations of an ensemble of lines (e.g., as in the EMPIRE survey; Bigiel et al. 2016). Density variations should induce large changes in ratios between the CO lines and high-density tracers. When the lines are sorted by effective critical density, changes in the mean density for a lognormal or similar steep and curving distribution should produce flaring line ratio variations. That is, the lines with highest effective critical density respond most strongly to density changes. When a power law tail is present, such flaring will be weaker, and the lines that sample the power law tail should change in lock step. The line ratio variations themselves will also be weaker in the presence of a power law tail, because such a tail ensures the presence of some high-density gas.

Thus, ensembles of lines can in principle be used to constrain the characteristic shape of the sub-resolution density distribution. Another application is to use a suite of lines to help control for abundance variations. In the bottom row of Figure 9, we have assumed the abundance of each molecule to remain fixed between beams, although we have not assumed any particular abundance values. In the case that one line exhibits a distinct behavior while the others indicate density variations, abundance variations (which do happen) represent a logical hypothesis.

Figure 9 does illustrate an issue with our target line suite. Our lines sort into two groups, which broadly exhibit similar behaviors. All four CO lines show weak line ratio changes with density, whereas all four dense gas tracers show strong variations. We have chosen the $\lambda \approx 3$ mm lines most readily observable in other galaxies, but this suite leaves us with limited sensitivity to intermediate densities.

### 3.5. Quantitative Estimates of Changing Gas Density

Beyond the qualitative analysis in the previous section, we would like to use mm-wave line ratios to infer quantitative changes in the underlying density distribution. We expect changes in the ratio, e.g., of HCN-to-CO to reflect differences in the fraction of dense gas, $f_{\text{dense}}$. Assuming that some characteristic density distribution holds from beam-to-beam, these ratios will also reflect changes in the median density of gas by mass, $n_{\text{med}}^{\text{mass}}$.

We test our ability to infer changes in $f_{\text{dense}}$ and $n_{\text{med}}^{\text{mass}}$ using a grid of models. Across this grid, we vary the mean density from $n_0 = 10^2$ to $10^4$ cm$^{-3}$, the distribution width from $\sigma = 0.6$ to 1.2 dex ($M \approx 5$ to 100), and the temperature from $T_{\text{kin}} = 15$ to 35 K. These values of $n_0$ and $\sigma$ span most of the values observed for whole clouds in normal and starburst galaxies (e.g., Leroy et al. 2015, 2016). We calculate line ratio patterns for each model in this grid (Table 3). Then, for all possible model pairs, we measure the logarithmic change in each line ratio,

$$ \Delta \log_{10} \text{HCN/CO} = \frac{\log_{10} I_{\text{HCN},1}/I_{\text{CO},1}}{\log_{10} I_{\text{HCN},2}/I_{\text{CO},2}}, \quad (8) $$

where the subscripts 1 and 2 refer to some pair of models in the grid and $\Delta \log_{10} \text{HCN/CO}$ is expressed in dex.

For each model pair, we also measure the logarithmic change in the median density by mass, $n_{\text{med}}^{\text{mass}}$, and the dense gas mass fraction. Here, $n_{\text{med}}^{\text{mass}}$ is defined as the density below which half the mass in the distribution lies. The dense gas fraction, $f_{\text{dense}}$, is defined as the fraction of mass above some threshold density, $n_{\text{thresh}}$. We take $n_{\text{thresh}} = 10^{4.5} \text{ cm}^{-3}$ as our default threshold and explore the appropriate definition of $n_{\text{thresh}}$ below.

#### 3.5.1. Estimating the Dense Gas Fraction

A primary use for ratios like HCN/CO is to estimate the fraction of dense gas mass. Figure 11 shows how well $\Delta \log_{10} \text{HCN/CO}$ predicts changes in the dense gas mass
<table><caption>**Table 3**
Beam Averaged Emissivity and Density Distributions</caption><thead><tr><th>Line</th><th>Distribution</th><th>log<sub>10</sub> n<sub>0</sub><br>(cm<sup>-3</sup>)</th><th>σ<br>(dex)</th><th>log<sub>10</sub> X<sub>mol</sub><br>(cm<sup>-3</sup>)</th><th>log<sub>10</sub> n<sub>med</sub><sup>mass</sup><br>(cm<sup>-3</sup>)</th><th>log<sub>10</sub> f<sub>dense</sub><br>(7)</th><th>T<sub>kin</sub><br>(K)</th><th>log<sub>10</sub> τ<br>(9)</th><th>log<sub>10</sub> &lt;ε&gt;<br>(K km s<sup>-1</sup> / cm<sup>-2</sup> / 10)</th><th>log<sub>10</sub> n<sub>med</sub><sup>emis</sup><br>(cm<sup>-3</sup>)</th></tr><tr><th>(1)</th><th>(2)</th><th>(3)</th><th>(4)</th><th>(5)</th><th>(6)</th><th>(7)</th><th>(8)</th><th>(9)</th><th>(10)</th><th>(11)</th></tr></thead><tbody><tr><td>12CO10</td><td>LOGNORMAL+TAIL</td><td>2.00</td><td>0.6</td><td>-4.00</td><td>2.41</td><td>-1.63</td><td>15.0</td><td>-1.00</td><td>-19.20</td><td>2.78</td></tr><tr><td>12CO10</td><td>LOGNORMAL+TAIL</td><td>2.00</td><td>0.7</td><td>-4.00</td><td>2.60</td><td>-1.31</td><td>15.0</td><td>-1.00</td><td>-19.14</td><td>2.98</td></tr><tr><td>12CO10</td><td>LOGNORMAL+TAIL</td><td>2.00</td><td>0.8</td><td>-4.00</td><td>2.80</td><td>-1.14</td><td>15.0</td><td>-1.00</td><td>-19.10</td><td>3.15</td></tr><tr><td>12CO10</td><td>LOGNORMAL+TAIL</td><td>2.00</td><td>0.9</td><td>-4.00</td><td>2.98</td><td>-1.03</td><td>15.0</td><td>-1.00</td><td>-19.07</td><td>3.27</td></tr><tr><td>12CO10</td><td>LOGNORMAL+TAIL</td><td>2.00</td><td>1.0</td><td>-4.00</td><td>3.11</td><td>-0.96</td><td>15.0</td><td>-1.00</td><td>-19.06</td><td>3.36</td></tr><tr><td>12CO10</td><td>LOGNORMAL+TAIL</td><td>2.00</td><td>1.1</td><td>-4.00</td><td>3.22</td><td>-0.92</td><td>15.0</td><td>-1.00</td><td>-19.04</td><td>3.42</td></tr><tr><td>12CO10</td><td>LOGNORMAL+TAIL</td><td>2.00</td><td>1.2</td><td>-4.00</td><td>3.30</td><td>-0.88</td><td>15.0</td><td>-1.00</td><td>-19.04</td><td>3.47</td></tr><tr><td>12CO10</td><td>LOGNORMAL+TAIL</td><td>2.25</td><td>0.6</td><td>-4.00</td><td>2.66</td><td>-1.50</td><td>15.0</td><td>-1.00</td><td>-19.11</td><td>2.93</td></tr><tr><td>12CO10</td><td>LOGNORMAL+TAIL</td><td>2.25</td><td>0.7</td><td>-4.00</td><td>2.85</td><td>-1.19</td><td>15.0</td><td>-1.00</td><td>-19.08</td><td>3.12</td></tr><tr><td>12CO10</td><td>LOGNORMAL+TAIL</td><td>2.25</td><td>0.8</td><td>-4.00</td><td>3.04</td><td>-1.01</td><td>15.0</td><td>-1.00</td><td>-19.05</td><td>3.28</td></tr></tbody></table>

**Note.** The table reports results for line emission for distributions of densities observed within a beam. Columns report: (1) emission line, (2) distribution shape, (3) mean density for lognormal part of density distribution, (4) rms dispersion in lognormal part of density distribution, (5) adopted abundance of species relative to H<sub>2</sub>, (6) median density by mass for the adopted distribution, (7) dense gas mass fraction, defined as gas above $n_{H2} = 10^{4.5} \text{ cm}^{-3}$, (8) adopted kinetic temperature, (9) adopted optical depth, (10) calculated beam-averaged emissivity of H<sub>2</sub>, (11) median density for emission.

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

fraction across our model grid. We plot results for pure lognormal distributions (left) and distributions with power law tails (right) separately. Figure 11 adopts a threshold density $n_{thresh} = 10^{4.5} \text{ cm}^{-3}$.

In both cases, $\Delta \log_{10} \text{HCN/CO}$ can predict changes in $f_{dense}$, but the two distributions yield different results. For the case with a power law tail, $\Delta \log_{10} \text{HCN/CO}$ relates to $f_{dense}$ via a roughly linear relation. We expect such behavior based on the arguments in Section 3.5, with the ratio in the mass above any two densities on the power law fixed for fixed $\alpha$. For most of our models, both HCN and $n_{thresh}$ used to define $f_{dense}$ sample power law part of the distribution. In this case, we expect a simple, linear relation between the dense gas tracer and $f_{dense}$. The scatter in our model reflects temperature variations and cases where the lognormal part of the distribution contributes to the calculation.

The relationship between $f_{dense}$ and $\Delta \log_{10} \text{HCN/CO}$ appears steeper for the pure lognormal (left panel of Figure 11). Following our argument in Section 3.5, the steep, curving form of the lognormal renders the ratio of area above any two densities sensitive to the mean of the distribution. This, in turn, leads to a steeper than linear relation between $f_{dense}$ and $\Delta \log_{10} \text{HCN/CO}$.

We plot results for HCN/CO, but report results for all of our dense gas tracers in Table 4. Following the logic above, for the power law tail case, all tracers show slopes near unity with modest ($\approx 50\%$) scatter across the model grid. When a fixed-slope power law tail is present, dense gas tracers behave as expected and do a good job of tracing the fraction of mass in dense gas. Tracers with lower effective critical densities show mildly nonlinear slopes and shallow scatter as the lognormal portion of the distribution begins to contribute to the emission in some cases. For a pure lognormal, the relationships between $\Delta \log_{10} \text{HCN/CO}$ and $\Delta \log_{10} f_{dense}$ are steeper, with power law index from 1.4 to 2.3. In this case, the line ratio variations predict $f_{dense}$ with larger scatter, a factor of $\sim 2$ across the grid.

Figure 11 shows results for dense gas defined by a threshold density $n_{thresh} = 10^{4.5} \approx 3 \times 10^4 \text{ cm}^{-3}$. This is also the definition adopted by Gao & Solomon (2004a). Figure 12 shows how the results depend on our definition of dense gas.

We vary $n_{thresh}$ and calculate (1) the best-fit power law slope relating $\Delta f_{dense}$ to $\Delta \text{HCN/CO}$, (2) the scatter among all model pairs about this best-fit relation, and (3) the scatter in $\Delta \log_{10} f_{dense}$ relative to the total scatter in the model grid. This latter is relevant because, for low values of $n_{thresh}$, the grid itself shows little variation; most models have a large fraction of their mass above $n_{thresh} \approx 10^{3.5} \text{ cm}^{-3}$. The contrast between the scatter about the fit and the scatter in the grid highlight the accuracy with which HCN/CO picks out $f_{dense}$.

Figure 12 again shows distinct behavior for the pure lognormal (left) and the case with a power law tail at high density (right). In the case of a power law tail, any threshold density above $\log_{10} n_{H2} \approx 4.5 \text{ cm}^{-3}$ yields a nearly linear relation between $\Delta \text{HCN/CO}$ and $\Delta f_{dense}$. In the case of the pure lognormal, the exact scaling depends strongly on the threshold density chosen. In the pure lognormal case, high threshold densities yield a steep slope and larger scatter about the fit.

Thus, when a fixed-slope power law tail is present in the density distribution, $\Delta \log_{10} \text{HCN/CO}$ and other dense gas tracers appear to be stable tracers of the dense gas mass fraction for most reasonable definitions. They capture $\Delta \log_{10} f_{dense}$ with modest scatter over a range of plausible temperature, density, and Mach number variations. For the case of a pure lognormal distribution, the relationship between line ratios and $f_{dense}$ is more unstable, and depends sensitively on adopted definitions. We return to these different results for different distributions below. Briefly, we prefer the case with a power law tail as more physical, given the presence of ongoing star formation in most regions targeted by extragalactic surveys. But addressing this question represents a natural next topic for observational and theoretical work.

### 3.5.2. Estimating the Median Density

The CO lines constrain the integral under the distribution, whereas dense gas tracers access the integral only above some high density. If the shape of the density distribution is known, combining these two pieces of information constrains $n_{med}^{mass}$.

Figure 13 shows results from our model grids using pairs of models to relate the logarithmic change in $n_{med}^{mass}$ to
<figure><img src="image_12.png" /><figcaption>**Figure 8.** Beam-averaged emissivity, $\langle \epsilon \rangle$, of HCN (1-0) (left) and CO (1-0) (right) for z lognormal density distribution as a function of the mean, $n_0$, and width, $\sigma$, of the distribution. As expected, the emissivity of HCN varies according to the fraction of high density material, which is higher for high $n_0$ and high $\sigma$. Optically thick CO exhibits only weak emissivity variations, consistent with its use as a tracer of the total H$_2$ column. The ratio of the two lines provides a diagnostic of the mass of high-density material in the beam.</figcaption></figure>

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

$\Delta \log_{10}$ HCN/CO. Again, gray contours enclose 68% and 95% of the model pairs at a given $\Delta \log_{10}$ HCN/CO. The colored points show the models binned by change in line ratio. Black lines show the best fit slope for change in line ratios <1.5 dex. Table 4 reports the best fit slopes and scatter for each dense gas tracer.

Changing line ratios can predict changes in the median $n_{H2}$ by mass across our model grid. However, they do so with a factor of ~2–3 rms accuracy and the slope of the best-fit scaling depends on the adopted sub-beam distribution. To first order, this calculation works by taking the integral above the effective critical density of some dense gas tracer, normalizing by the integral under the whole distribution (traced by CO), and relating this to the median of the distribution. The mapping is not perfectly one to one across the model grid, but does work.

This calculation works about equally well for the two distributions, but it seems better posed for the pure lognormal. In this case, the whole distribution of densities obeys a single simple functional form. Thus, it makes sense to constrain the whole distribution by contrasting the high-density wing and the integral under the whole curve. In the case of the power law tail, the results rely sensitively on the adopted threshold joining the two distributions and slope of the power law tail.

### 3.5.3. Using Ensembles of Lines

So far, we have considered how changes in a single line ratio map to changes in the underlying density distribution. We saw results for different lines in Table 4 and discussed the qualitative use of multiple lines in Section 3.4. Modern surveys will often capture many of these transitions (e.g., see Bigiel et al. 2016), so that a pattern of ratios analogous to what we saw in Section 3.4 are available for each source.

In Figure 14, we explore the use of ensembles of lines to fit the underlying density distribution. We use the model grid from the previous sections, varying $n_0 = 10^2-10^4 \text{ cm}^{-3}$, $\sigma = 0.6-1.2$ dex, and $T_{\text{kin}} = 15-35 \text{ K}$. Again, we consider

pairs of models. A pair of model implies a change in line ratios, dense gas fraction, $n_{\text{med}}^{\text{mass}}$, $n_0$, and $\sigma$. We draw 1000 random pairs of models from the larger model grid. For each drawn model pair, A, we imagine that we have observed the line ratio variations implied by this model pair. Then, we search the whole grid for all model pairs that produce the same line ratio variations within a tolerance of 0.05 dex (~12%). This is a typical calibration uncertainty for mm-wave observations. These are “good fit” model pairs, that would match the line ratio pattern for A within the uncertainties.

We compare the variation in $f_{\text{dense}}$, $n_0$, $\sigma$, and the median $n_{H2}$ for good fit model pairs to those known for A. Based on the scatter in these parameters for good-fit model pairs about the true values for A, we gauge the rms accuracy with which a given line suite can access the true variations in each quantity.

We repeat this exercise using only the HCN/CO ratio, then using both the HCN/CO and HCO$^+/\text{CO}$ ratio, then using CO along with all four dense gas tracers, and finally using all lines in our model. For comparison, we measure the scatter across the model grid itself. We show the results for a lognormal distribution with (left) and without (right) a tail in Figure 14.

Figure 14 illustrates a few key points. First, constraints on the distribution width, $\sigma$, and mean density, $n_0$, are individually weak. The combination of the two produces $n_{\text{med}}^{\text{mass}}$, and this quantity is better constrained. Fitting using only HCN/CO already improves the accuracy dramatically compared to randomly drawing from the model grid. Adding additional lines further improves the ability to recover the $n_{\text{med}}^{\text{mass}}$ for lognormal distributions with and without tails. Each additional line helps better constrain this quantity, including folding in the full suite of CO transitions (“All Lines”).

Unsurprisingly, the dense gas fraction is the best-constrained quantity in this exercise. Here, too, additional lines add accuracy, albeit the gain from additional dense gas tracers is relatively modest, improving the accuracy of the fit by ~0.1 dex. Again, folding in $^{13}$CO and several low-$J$ CO lines
<figure><img src="image_14.png" /></figure>

**Figure 9.** Beam-averaged emissivity, $\langle\epsilon\rangle$, and line ratio patterns as a function of changing mean density, $n_0$. We adopt our fiducial $\tau$ and take $\sigma = 0.8$ dex and $T_{\text{kin}} = 25$ K throughout. (top left) $\langle\epsilon\rangle$ as a function of $n_0$ for each line, normalized to $\langle\epsilon\rangle$ for that line at $n_0 = 10^3 \text{ cm}^{-3}$. Here, we plot results for pure lognormal distributions. The emissivity of a distribution in the CO lines varies only weakly with its mean density. However, $\langle\epsilon\rangle$ for the high-density tracers depends strongly on the sub-grid density distribution. (top right) Implied line ratio pattern relative to CO for our fiducial abundances and different $n_0$, again for pure lognormal distributions. Line ratios between dense gas tracers and CO lines vary strongly with density. Internal ratios among CO lines and dense gas tracers vary more weakly. The exact pattern depends on our adopted abundances. (bottom row) Line ratio patterns, now normalized to the pattern at $n_0 = 10^3 \text{ cm}^{-3}$, for various $n_0$. We calculate results for both a pure lognormal distribution (bottom left) and one with a power law tail (bottom right). These variations in line ratios are robust to the absolute abundance pattern, but will be affected by abundance variations. For both distributions, the most density-sensitive lines show the strongest variations (e.g., HCN (1-0), HNC (1-0), CS (2-1)) relative to CO. In the lognormal case, a changing density distribution creates a “flaring” pattern of line ratio variations when the lines are sorted by effective critical density. This signature is far less pronounced when the power law tail dominates the density distribution (e.g., this is the case at low $n_0$ in the bottom right plot).

improves our recovery of $f_{\text{dense}}$, dropping uncertainties to ~0.1 dex.

The exact numerical results in Figure 14 depend on our adopted model grid, knowing the optical depths of our lines, and other details. However, the qualitative conclusions should be robust. The mean density and distribution width are individually less well-constrained than $n_{\text{med}}^{\text{mass}}$, which they combine to set. In the lognormal case, one does some constraint on $\sigma$, but the low dynamic range expected is an issue. The dense gas mass fraction is already constrained well

by a single line ratio. This accuracy can be improved by the combination of many lines. Perhaps surprisingly, adding a suite of low-$J$ CO lines can substantially improve the fit.

Note that, in practice, the model grid does not represent a suite of equally likely conclusions. Using external knowledge of the temperature or the distribution width (e.g., from the Mach number gauged via measured velocity dispersions) will improve the ability to infer the full density distribution. Moreover, as discussed above, a key role of multiple dense gas tracers is to provide some robustness to abundance and
<figure><img src="image_15.png" /><figcaption>**Figure 10.** As in the bottom row of Figure 9, but now varying $T_{\text{kin}}$ while other quantities remain fixed. We plot line ratios relative to those expected for $T_{\text{kin}} = 35$ K, so that the figure shows the expected changes to the observed line ratio pattern due to temperature variations. These are much weaker than the density variations seen in Figure 9.</figcaption></figure>

optical depth variations. That role is not reflected in Figure 14.

### 3.6. Dense Gas Conversion Factors

We have emphasized tracking changes in the gas density via changing line ratios. A more direct approach is to simply translate the flux of HCN, HCO$^{+}$, or another dense gas tracer into a mass of dense gas. Gao & Solomon (2004a, 2004b) suggest $\alpha_{\text{HCN}} \approx 10 M_{\odot} \text{pc}^{-2} (\text{K km s}^{-1})^{-1}$ to translate HCN luminosity to the mass of gas at $n_{\text{H2}} \gtrsim 3 \times 10^4 \text{ cm}^{-3}$. Our calculations yield conversion factors for each line and each grid point.

How good of an approximation is a fixed $\alpha_{\text{HCN}}$ or a fixed $\alpha_{\text{HCO}^+}$? Considering the same range of densities, distribution widths, and temperatures as above, we calculate the conversion factor relating each dense gas tracer to the dense gas mass for each model. Following the previous section and Gao & Solomon (2004a), we adopt a threshold of $n_{\text{H2}} \gtrsim 3 \times 10^4 \text{ cm}^{-3}$ to define dense gas. In this exercise, we consider distributions with and without a power law tail together.

For each line, we consider models with these physical parameters described by a range of optical depths. At each $\tau$, we find the mean and the scatter in $\log_{10} \alpha$, which we plot in Figure 15. As a control, we show the scatter in the dense gas conversion factor for CO (1-0) as a gray region. We expect (and find) CO on its own to do a poor job of tracing the dense gas. When the scatter in a dense gas conversion factor approaches this value, the utility of that conversion factor is also limited. All conversion factors produced by the model will scale with our adopted abundance, taken as $X(\text{mol}) = 10^{-8}$ for all lines here. We summarize the results in Table 5.

We calculate $\alpha_{\text{dense,HCN}}$ between about 3 and $30 M_{\odot} \text{pc}^{-2} (\text{K km s}^{-1})^{-1}$. For matched abundances, HCN and HNC yield

similar conversion factors. HCO$^{+}$ produces a lower conversion factor, reflecting its lower effective critical density, and CS (2-1) has a moderately high conversion factor. In each case, the scatter at fixed optical depth is modest in the optically thin case, but becomes larger when the line becomes optically thick. Most simply, this reflects temperature variations, which will linearly affect $\alpha$ when the line is thick.

Accounting for the two times higher HCN abundance assumed by Gao & Solomon (2004b), these calculations agree well with their $\alpha_{\text{HCN}} \approx 10 M_{\odot} \text{pc}^{-2} (\text{K km s}^{-1})^{-1}$ for densities above $\sim 3 \times 10^4/\tau \text{ cm}^{-3}$. This general agreement belies significant uncertainty. Both the optical depth and the abundance (and thus chemistry) must be known to calculate the conversion factor using this approach. In the optically thick regime, the temperature also plays a large role. To first order, all of these quantities act linearly, so that even a factor of 2 uncertainty in each (which seems optimistic) implies a factor of $\sim 3.5$ uncertainty.

To understand the uncertainty in dense gas conversion factors, it may be helpful to contrast these with the CO-to-H$_2$ conversion factor (see review in Bolatto et al. 2013). As Figure 8 shows, $\langle \epsilon \rangle$ (and so $\alpha_{\text{CO}}$) of CO (1-0) varies only weakly across a wide range of density distributions, so that $\alpha_{\text{CO}}$ does vary strongly as a function of the sub-beam density distribution, all other things being equal. Meanwhile, the assumption of a fixed dynamical state, either marginal boundedness or virialization, has been widely argued and observed for whole clouds (e.g., see Bolatto et al. 2013; Heyer & Dame 2015). More, the $^{12}\text{CO}$ lines usually are optically thick, if only because CO represents the dominant carbon reservoir in the well-shielded parts of molecular clouds. In the CO case, the line width of a cloud is set by the cloud dynamical state. The high opacity combines with modest temperature variations to fix the specific intensity in the line at a weakly varying brightness temperature. Then, at fixed metallicity, with $\alpha_{\text{CO}} \propto \rho^{0.5} T_{\text{kin}}^{-1}$ (e.g., Maloney & Black 1988; Narayanan et al. 2012; Bolatto et al. 2013, and many others).

Neither high opacity nor a fixed dynamical state can be taken as a given for the gas traced by HCN, HCO$^{+}$, or similar lines. Unlike CO, these do not represent a main carbon reservoir and their absolute abundances are often quite low. These lines are observed to have some opacity, but not the pervasive optical thickness of $^{12}\text{CO}$. Meanwhile, the densities traced by these lines represents a subset of the cloud. For clouds of changing mean density, virial parameter, and with different self-gravitating substructures, the dense gas tracers may or may not trace exactly the emission from this tail. Thus, although Gao & Solomon (2004b) argue for a fixed dynamical state for HCN-emitting regions, there does not appear to be a strong a priori reason to take specifically the HCN-emitting gas to be in a particular, universal dynamical state.

Instead of assuming a dynamical state for a particular line, the natural place to introduce such a criteria to our model is to couple the column and line width integrated over the power law tail, which could be taken at least to be gravitationally bound. The issue is that this couples the column density of H$_2$ to the line width, again requiring one to adopt an abundance for $N_{\text{H2}}/\Delta v$ to imply an optical depth in any given line.

Thus, Figure 15 and Table 5 give approximate conversion factors consistent with previous work. These highlight the systematic variation of $\alpha$ with optical depth. They also give a sense of the relative conversion factors for different lines with
<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 11.** Mapping between changes in HCN-to-CO line ratio (x-axis) to changes in the dense gas mass fraction $f_{\text{dense}}$ (defined relative to $n_{\text{thresh}} > 10^{4.5} \text{ cm}^{-3}$). We calculate cloud-integrated emissivities for a range of densities ($\log_{10} n_0 = 2-4 \text{ cm}^{-3}$), distribution widths ($\sigma = 0.6-1.2 \text{ dex}$), and temperatures ($T_{\text{kin}} = 15-35 \text{ K}$). Then, we pair calculations and measure the logarithmic contrast in line ratios and dense gas fractions. Contours show the 68% and 95% range of model pairs at for each change in line ratio. Colored points show binned results. Lines in both panels indicate the best fit scaling for $\Delta \log_{10} \text{HCN/CO} < 1.5$.</figcaption></figure>

<table><caption>**Table 4**
Line Ratio Predictors of Gas Density</caption><thead><tr><th rowspan="2">Line Relative<br>to CO (1-0)</th><th colspan="2">Lognormal + Tail</th><th colspan="2">Pure Lognormal</th></tr><tr><th>Slope<sup>a</sup></th><th>Scatter<sup>b</sup><br>(dex)</th><th>Slope<sup>a</sup></th><th>Scatter<sup>b</sup><br>(dex)</th></tr></thead><tbody><tr><td colspan="5" style="text-align:center; font-weight:bold;">Predicting f<sub>dense</sub></td></tr><tr><td>HCO<sup>+</sup></td><td>1.21</td><td>0.19</td><td>2.31</td><td>0.35</td></tr><tr><td>CS</td><td>1.08</td><td>0.18</td><td>1.86</td><td>0.33</td></tr><tr><td>HNC</td><td>1.05</td><td>0.15</td><td>1.70</td><td>0.27</td></tr><tr><td>HCN</td><td>0.98</td><td>0.16</td><td>1.41</td><td>0.26</td></tr><tr><td colspan="5" style="text-align:center; font-weight:bold;">Predicting Median n<sub>H2</sub> by Mass</td></tr><tr><td>HCO<sup>+</sup></td><td>2.03</td><td>0.44</td><td>1.68</td><td>0.51</td></tr><tr><td>CS</td><td>1.85</td><td>0.41</td><td>1.54</td><td>0.47</td></tr><tr><td>HNC</td><td>1.83</td><td>0.35</td><td>1.47</td><td>0.42</td></tr><tr><td>HCN</td><td>1.73</td><td>0.36</td><td>1.34</td><td>0.42</td></tr></tbody></table>

**Notes.** Power law index relating the change in line ratio to the change in median H₂ density by mass, $n_{\text{med}}^{\text{mass}}$, and dense gas fraction, $f_{\text{dense}}$. Models run over mean densities from $\log_{10} n_0$ [cm⁻³] = 2 to 4, $T_{\text{kin}}$ = 15 to 35 K, and $\sigma$ = 0.6 to 1.2 dex. We report results with and without a power law tail separately. Models assume our fiducial optical depths.

<sup>a</sup> Best fit slope, $\alpha$, of a scaling $y = \alpha x$, that goes through the origin, with $x$ the logarithmic change in line ratios, e.g., $\log_{10} (\text{HCN/CO})_1 - \log_{10} (\text{HCN/CO})_2$, and $y$ the logarithmic change in $n_{\text{med}}^{\text{mass}}$.

<sup>b</sup> Indicative logarithmic scatter about the power law scaling. Reflective of ability of ratio to predict changes in gas density. For comparison, the whole model grid (with no fit) has scatter ~1 dex in $n_{\text{med}}^{\text{mass}}$ and ~0.5–0.7 dex in $f_{\text{dense}}$.

<sup>c</sup> Fraction of mass at densities above $n_{\text{H2}} = 10^{4.5} \text{ cm}^{-3}$.

different abundances. But these calculations remain subject to large systematic uncertainties related to the sub-beam density distribution, opacity, abundance (chemistry), and temperature.

Many of these factors can divide out of a relative approach. Our view remains that, because $\tau$ is observationally accessible, fixing $\tau$ and comparing line ratios offers a more robust path than invoking an absolute conversion factors.

### 3.7. Implications for Isotopologue Studies

Our model assumes that the optical depth is known, or can be known. The main driver for this assumption is that rare, optically thin isotopologues offer the chance to gauge the optical depth of the main dense gas tracing lines. Alternatively, such observations may constrain the isotopic abundance, which can give clues to recent nucleosynthesis or enrichment patterns as well as cloud chemistry.

The optical depth of a line changes its $\epsilon$ ($n_{\text{H2}}$), and thus its sensitivity to gas density. When a distribution of densities exists within the beam, the thin transition will arise from a denser subset of the gas than the optically thick main line. In distributions like those that we consider, this can have the effect of depressing emission from optically thin isotopologues of optically thick dense gas tracers. The effect will be strongest when the mean density of the cloud is lower than the effective critical density of the lines in question, such that density strongly affects emissivity. Thus, we expect this to be a larger effect for HCN or HCO⁺ than for CO.

Figure 16 shows this effect. The key physics appear in the ratio of the beam averaged emissivity at some optical depth, $\langle\epsilon\rangle(\tau)$, to the beam-averaged emissivity in the optically thin case, $\langle\epsilon\rangle(\tau) \approx 0.1)$. In the limit of $\tau \gg 1$, for LTE and otherwise fixed physical conditions, we expect this ratio to approach $\tau^{-1}$. If the relative abundance of the isotopologues being studied is known, one is optically thin, and one is optically thick, then the abundance ratio is multiplied by this $\tau^{-1}$ to produce the expected line ratio (e.g., CO typically has
<figure><img src="image_18.png" /><figcaption>**Figure 12.** Effect of choice of density threshold to define $f_{\text{dense}}$ on the ability of HCN/CO to predict $f_{\text{dense}}$. Best fit slope, scatter, and improvement over scatter in the model grid for various definitions of $f_{\text{dense}}$ for a pure lognormal (left) and a lognormal with a power law tail (right). When a power law tail is present, our results are robust to choice of threshold above ~$10^{4.5}$ cm$^{-3}$, whereas the appropriate fit for the pure lognormal depends sensitively on the adopted threshold.</figcaption></figure>

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

$\tau \approx 10$ and a $^{12}\text{C}/^{13}\text{C} \approx 50$; therefore typical $^{12}\text{CO}/^{13}\text{CO}$ ratios are $\sim 50 \times 10^{-1} \approx 5$).

Figure 16 plots deviations from this simple case due to differential density sensitivities, which also incorporate excitation. In the left panel, we hold all other things fixed, and we compare the emissivity of HCN at $\tau = 1, 2, 5,$ and 10 to HCN at $\tau = 0.1$. For distributions centered at high densities, the expectation of a ratio approaching $\tau^{-1}$ roughly holds. However, as we change the mean density of the distribution, $n_0$, the ratio changes. The sense is that, for lower mean densities, the beam-averaged emissivity of the optically thin gas becomes much worse relative to the optically thick gas. This is most extreme for the lognormal case at low mean densities. However, even when a power law tail is present, the effect can be large.

The right panel shows the same plot, taking $\tau = 5$ for all four dense gas tracers and $^{12}\text{CO}$ (1-0). Although $^{12}\text{CO}$ shows some differential excitation effects at low densities, the effect here is far more modest. All of the dense gas tracers show significant deviations from the simple expectation for mean densities $n_0 \lesssim 10^3$ cm$^{-3}$, with $\text{HCO}^+$ requiring lower densities to show these variations. These differential excitation effects become important when an appreciable amount of gas lies below the critical density of the line in question. This makes isotopologue studies using CO somewhat robust, and $\text{HCO}^+$— which has an intermediate critical density—more robust than HCN or HNC, all other things being equal.

Taking the lognormal with a tail, Figure 16 shows that, for $\tau \approx 5$, the magnitude of the correction from the simple case is a factor of ~2–3. That is, the ratio of an optically thin tracer like $\text{H}^{13}\text{CN}$ to HCN will be ~2–3 times lower than the expectation for the simple case if $n_0 \lesssim 10^3$ cm$^{-3}$. The difference, due to differential excitation, should not be confused with optical depth or abundance effects.

To apply these results generally, one can take our tabulated calculations and, for matched conditions, compare $\langle \epsilon \rangle(\tau)/\langle \epsilon \rangle(\tau \approx 0.1)$ to either the simple expectation ($\sim \tau^{-1}$)

or the value found at high densities. This factor represents the part of the line ratio due only to differential excitation and not abundance or optical depth effects.

Optically thin isotopologues remain a crucial way to probe the optical depth of mm-wave lines and isotopic abundance, especially in other galaxies. Indeed, given the importance of optical depth to understand the density-sensitivity of the lines that we discuss, more such observations are crucial. However, Figure 16 cautions that sub-resolution density distributions and differential excitation should be borne in mind when interpreting observed ratios. Without taking these into account, it would be easy to solve for an optical depth lower than the true value when observing low-density gas. Considering isotopologue ratios for HCN, HNC, and $\text{HCO}^+$, M. J. Jimenez-Donaire et al. (2017, in preparation) demonstrate the need for such corrections in real measurements. Their synthesis of the literature and new measurements and limits show high HCN/$\text{H}^{13}\text{CN}$ and $\text{HCO}^+/\text{H}^{13}\text{CO}^+$ ratios in regions where $n_0$ should be $\lesssim 10^3$ cm$^{-3}$. These ratios yield different interpretations in the LTE case and the case described here.

## 4. Discussion and Summary

Integrating over a wide range of densities within a single telescope beam is unavoidable when studying molecular gas at extragalactic distances. To date, there has been limited effort to interpret spectral line observations in the context of realistic sub-beam distributions, with modeling mostly considering one or two-phase media. To move forward, mm-wave line emission must be modeled taking into account realistic sub-resolution distributions, similar to the treatment of stellar populations in population synthesis or multiple dust populations in IR SED modeling.

In this paper, we consider the interaction of realistic density distributions with the emissivity of common mm-wave transitions. We focus on a suite of transitions in the 85–115 GHz range that are commonly used to trace density.
<figure><img src="image_20.png" /><figcaption>**Figure 13.** Using changes in the HCN/CO line ratio to infer changes in the median $n_{H2}$ by mass. As in Figure 11, but now both panels relate HCN-to-CO line ratio (x-axis) to changes in $n_{\text{med}}^{\text{mass}}$ (y-axis).</figcaption></figure>

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

Building on work by Krumholz & Thompson (2007), our model treats emission within a beam as the sum of emission from a collection of one-zone models that share a characteristic optical depth, $\tau$, and temperature, $T_{\text{kin}}$, but have a realistic distribution of collider densities, $n_{\text{H2}}$.

To implement this, we use the one-zone non-LTE code RADEX (van der Tak et al. 2007) to calculate the emissivity, $\epsilon$, of gas in each line across a range of $n_{\text{H2}}$, $T_{\text{kin}}$, and $\tau$ giving $\epsilon$ ($n_{\text{H2}}$, $T_{\text{kin}}$, $\tau$). We calculate emission from a density distribution by combining one-zone models in which the escape probability (via $\tau$) is fixed, while varying the collider density, $n_{\text{H2}}$. We implement two density distributions: a pure lognormal and a lognormal distribution that exhibits a power law tail at high densities. Theoretical, numerical, and observational work all commonly invoke this combination of density distributions as a reasonable description of the cold, turbulent gas that produces low-$J$ mm-wave line emission (e.g., Vazquez-Semadeni 1994; Padoan & Nordlund 2002; Krumholz & McKee 2005; Kainulainen et al. 2009; Federrath & Klessen 2013, among many others), albeit there are some caveats from recent Milky Way work (Lombardi et al. 2015).

Whereas Krumholz & Thompson (2007) and the closely related paper by Narayanan et al. (2008) mainly focused on star formation scaling relations, we are interested in the use of mm-wave line emission to trace gas density within a telescope beam. Because of the small physical scales involved, spectroscopy offers almost the only way to access the small-scale density distribution across a wide range of extragalactic environments. Despite the faintness of the high-effective-density lines like HCN (1-0), HCO$^+$ (1-0), and CS (2-1), the Atacama Large Millimeter/submillimeter Array, the Green Bank Telescope (GBT), and the IRAM telescopes are now able to regularly observe these transition across the disks of normal galaxies (e.g., Kepley et al. 2014; Usero et al. 2015; Bigiel et al. 2016, M. J. Gallagher et al. 2016 in preparation, among others).

The paper includes a large set of tabulated results available as online-only material (Tables 2 and 3). These report:

1. The minimum collider density, $n_{\text{H2}}$, at which each line achieves 95% of its maximum emissivity for a given $T_{\text{kin}}$ and $\tau$ in a one-zone model (Table 2). This is closely related to the effective critical density (e.g., Shirley 2015, and references therein). Having such values tabulated for the non-LTE case, cast directly in terms of emissivity and calculated across the transition from optically thin to optically thick lines, may be useful to readers.

2. The beam-averaged emissivity, $\langle\epsilon\rangle$, and median density for emission in each studied transition for a range of $T_{\text{kin}}$, $\tau$, and density distributions comparable to those found across the local galaxy population (Table 3). From these, line ratio patterns can be constructed for any adopted set of optical depth and temperature.

### 4.1. Using Line Ratios to Trace Density

We examine the ability of mm-wave line ratios to trace the sub-beam density distribution. For these purposes, the interaction between the emissivity, $\epsilon$ ($n_{\text{H2}}$, $T_{\text{kin}}$, $\tau$), and the density distribution, $P(n_{\text{H2}})$, is crucial. From a mixture of one-zone models and combinations that model a density distribution, our conclusions include:

1. Gas can still emit effectively at densities well below its effective critical density. That is, the emissivity, $\epsilon$ ($n_{\text{H2}}$), of a mm-wave transition still has high values for at least a decade in $n_{\text{H2}}$ below the density of peak emissivity (see Figure 2).

2. The flux that emerges from each density in the cloud is the product of the mass at that density and the emissivity at that density. For realistic distributions, $P(n_{\text{H2}})$ drops with increasing $n_{\text{H2}}$ near the critical density of the dense
<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /><figcaption>**Figure 14.** Accuracy of recovered model inputs requiring all line ratios to match the model within a tolerance of 0.05 dex. If fitting observations, these would also be good fits to the data. The two panels show results for a pure lognormal (left) and a lognormal distribution with a powerlaw tail (right). From left to right in each panel, we show the scatter in the model grid itself (not a fit), a fit using only the HCN/CO ratio, and then fits adding in HCO$^{+}$, all dense gas tracers, and all lines. The figure shows that $f_{\text{dense}}$ can be recovered most accurately. The median density by mass, $n_{\text{med}}^{\text{mass}}$, can also be fit if the distribution is known. Adding more lines improves the fit. The most significant constraints come from using the first line and folding in all lines, including those with critical densities $\lesssim$ the mean of the distribution.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Figure 15.** Conversion factors relating line emission to dense gas mass for a dense gas threshold of $n_{\text{thresh}} = 10^{4.5} \text{ cm}^{-3}$. The gray region shows the range of *dense* gas conversion factors for CO, and so represents a control (this much scatter offers little predictive power). The colors show results for HCO$^{+}$ (purple), HCN (red), HNC (blue), and CS (green), all using our fiducial abundance of $10^{-8}$.</figcaption></figure>

gas tracers. This often leads to the case where the density distribution and emissivity exert competing effects: there is more mass at lower densities, but gas emits better at high densities (see Figure 5).

<table><caption>**Table 5**
Dense Gas Conversion Factors</caption><thead><tr><th>Line</th><th>τ &le; 1</th><th>1 &le; τ &le; 3</th><th>τ &ge; 3</th></tr><tr><td></td><td colspan="3" style="text-align:center;">α<sub>dense</sub> ± scatter<br>(M<sub>☉</sub> pc<sup>-2</sup> (K km s<sup>-1</sup>)<sup>-1</sup>) ± (dex)</td></tr></thead><tbody><tr><td>HCO<sup>+</sup></td><td>1.5 ± 0.16 dex</td><td>2.3 ± 0.23 dex</td><td>5.3 ± 0.34 dex</td></tr><tr><td>CS</td><td>8.2 ± 0.11 dex</td><td>11.4 ± 0.16 dex</td><td>23.7 ± 0.28 dex</td></tr><tr><td>HNC</td><td>3.4 ± 0.08 dex</td><td>4.8 ± 0.14 dex</td><td>9.8 ± 0.25 dex</td></tr><tr><td>HCN</td><td>4.5 ± 0.06 dex</td><td>6.3 ± 0.10 dex</td><td>12.2 ± 0.21 dex</td></tr></tbody></table>

**Note.** Median and scatter in $\alpha_{\text{dense}}$ across the model grid where $n_0 = 10^2-10^4 \text{ cm}^{-3}$, $\sigma = 0.6-1.2$ dex, and $T_{\text{kin}} = 15-35$ K. Here, dense gas is defined using a threshold density $n_{\text{thresh}} = 10^{4.5} \text{ cm}^{-3}$. All calculations assume an abundance $X(\text{mol}) = 10^{-8}$. The quoted scatter does not account for uncertainties in this quantity, which affects the answer linearly. Note that the table quotes linear $\alpha_{\text{dense}}$ but logarithmic scatter (in dex).

3. For a pure lognormal distribution, the interaction between $\epsilon$ ($n_{\text{H2}}$) and $P(n_{\text{H2}})$ leads to large variations in the median density producing emission from high-density tracers. We illustrate this for HCN (1-0), and a similar case holds for other lines. In cases of low mean density, $n_0$, and narrow distribution width, $\sigma$ (expected for low Mach number), the median density of gas producing HCN emission can fall well below the critical density of HCN. When this happens, the overall emissivity of the gas in HCN also drops. Thus, for regions in which modest density gas produces most HCN emission, HCN/CO and similar line ratios will also be low (see Figure 6).

4. The presence of a power law tail in the density distribution tends to suppress variations in the median density for emission. Such tails are expected for self-gravitating gas and lack a preferred scale. As a result,
<figure><img src="image_25.png" /></figure>

**Figure 16.** Effect of differential excitation on measured isotopologue ratios. Beam averaged emissivity of an optically thick line divided by that of the same line when optically thin. The simple LTE case with $\tau \gg 1$ and LTE predicts a ratio of $\tau^{-1}$. However, taking into account a distribution of densities, emission from the optically thin line is often weaker than one would expect. This results in apparent suppression of the line ratios often used to gauge optical depth, one that must be accounted for to avoid highly biased results in many common cases. An analogous effect should be included to estimate $T_{\text{kin}}$ from multi-J observations that integrate over density distributions.

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

they produce much weaker variations in the median density for emission. For example, HCN emission from a distribution with our fiducial power law tail almost always arises from $\sim 10^4-10^5 \text{ cm}^{-3}$ gas (see Figure 6).

The sub-beam density distributions also affects the beam-averaged emissivity, $\langle \epsilon \rangle$, so that some lines emit better for certain density distributions.

5. Variations in $\langle \epsilon \rangle$ are strongest for the dense gas tracers, which have high $\langle \epsilon \rangle$ when more dense gas is present. $\langle \epsilon \rangle$ remains more nearly constant for the low-J CO lines, reflecting their common use as total gas tracers. As a result of these differential variations in $\langle \epsilon \rangle$, line ratios have the power to capture changes in the density distribution within the beam (see Figures 8 and 9).

The sense of these variations is that lines with effective critical densities low compared to the mean density, like the low-J CO transitions, vary little in their emissivity. Transitions with higher effective critical densities, compared to the mean density, see their emissivity vary strongly as the density distribution changes. Thus, ratios like HCN-to-CO have the power to probe the fraction of dense gas or the mean gas density within a beam.

The exact value of $\langle \epsilon \rangle$, or any individual line ratio, depends on our adopted abundances. These are often poorly known. To help circumvent this concern, we highlight the power of line ratio variations.

6. In the case that two regions have different density distributions and the abundances remain fixed (but also unknown), the measured variation in line ratios has the power to trace the change in density distribution within the beam, independent of the absolute adopted abundance (see Figure 9).

7. These line ratio variations will be strongest for lines with the highest effective critical densities. Thus, the signature of a changing density distribution is a “flaring” pattern of line ratios when the lines are sorted as a function of increasing critical density (see Figure 9).

8. The shape of this flare reflects the high-end form of the density distribution. For a self-similar distribution like a power law, the ratios of all lines on the power law tail relative to CO will vary in lock-step. For a downward curving, steep distribution like a lognormal, line ratio variations will increase in strength with increasing critical density (see Figure 9).

Beyond accessing the shape of the high density tail, line ratio patterns allow some prospect to control for abundance variations. If several high-density tracers are observed, then consistency (or lack thereof) among their variations can help distinguish changing abundances from changing density.

9. We also examine the impact of temperature variations on changing line ratio patterns. In general, $T_{\text{kin}}$ variations exerts a weaker influence on line ratios than changing density because they affects all lines, including those with with low critical densities. However, $T_{\text{kin}}$ does still have a some impact on the line ratios because it lowers the effective critical density of the dense gas tracers (see Figure 10).

We quantify how changes in the line ratio map to variations in the dense gas mass fractions, $f_{\text{dense}}$, and $n_{\text{med}}^{\text{mass}}$. Here, the shape of the high end of the density distribution plays a large role.

10. When a power law tail describes the distribution at high densities, HCN/CO, HCO$^+/CO$, and similar ratios trace variations in $f_{\text{dense}}$ with approximately linear slope. This remains true for most reasonable definitions of “dense”
gas. That is, if a power law tail is present and approximately universal, then the interpretation of dense gas tracers may be relatively straightforward (see Figure 11).

11. On the other hand, if the density obeys a steep, curving distribution like a lognormal, the interpretation of line ratio variations becomes more complex. We calculate scaling relations relating variations in HCN/CO and other line ratios to variations in $f_{\text{dense}}$, but caution that these are nonlinear and depend on the adopted threshold density. Still within a factor of ~2–3 scatter, these ratios can capture variations in the dense gas mass fraction and $n_{\text{med}}^{\text{mass}}$ (see Figures 11 and 13).

Given strong priors on the Mach number or mean density, the grid of models that we consider could be narrowed considerably. In this case, the inference of changes in $f_{\text{dense}}$ and $n_{\text{med}}^{\text{mass}}$ from line ratio variations can be even more precise.

Current surveys often capture a suite of line ratios. We explore the use of an ensemble of ratios to estimate changes in the sub-resolution density distribution. We do this by directly fitting line ratio variations to our grid of models. This approach can be implemented generally using the online tables included with the paper.

12. We show that including multiple lines improves the accuracy with which variations in the sub-beam density distribution are recovered. Variations in the mean density, $n_0$, and distribution width, $\sigma$, are, in general, less well-recovered than their combination, $n_{\text{med}}^{\text{mass}}$. Changes in the dense gas fraction, $f_{\text{dense}}$, tend to be better recovered than $n_{\text{med}}^{\text{mass}}$ (see Figure 14).

13. Including more lines in the fit steadily improves the accuracy with which the model variations are recovered. The best fit comes from including all dense gas tracers and multiple tracers of lower-density gas, here meaning the low-$J$ and optically thin CO lines (see Figure 14).

## 4.2. Differential Excitation and Isotopologue Studies

Our model specifies the optical depth rather than the abundance. One motivation for this is that the optical depth of mm-wave transition can be constrained by observations. The most direct way to do this is to pair observations of the main tracers with optically thin isotopologues.

14. We show that, in the presence of the sub-beam density distributions considered by our analysis, differential excitation between the main line and the optically thin isotopologue complicates this analysis. Even treating the two molecules as identical, the optically thin line has different $\epsilon$ ($n_{\text{H2}}$). This, convolved with the sub-beam density distribution, can affect the line ratio. If one interprets observations of the optically thin tracer without taking these effects into account, the optical depth inferred can be biased by a factor of $\gtrsim 2-3$ (see Figure 16).

15. The key physics behind this effect are captured in the ratio $\langle\epsilon\rangle(\tau)/\langle\epsilon\rangle(\tau=0.1)$ for fixed $T_{\text{kin}}$ and density distribution. This is the ratio of beam-averaged emissivity of a molecule at some optical depth to the beam-averaged emissivity for that molecule when optically thin. We show this factor for some representative cases, and it can be calculated from our tabulated model grid. When a

power law tail is present and the mean density is low compared to the critical density of the thin line, we arrive at the factor of 2–3 connection mentioned above. For a pure lognormal centered at low mean density, the effects can be even more extreme (see Figure 16).

16. Although we do not focus on temperature in this paper, a related approach to differential excitation should be used to infer temperatures from line ratios. $\epsilon$ ($n_{\text{H2}}$) varies with $T_{\text{kin}}$ and with transition; e.g., see the case for CO in Figure 2. When using multi-$J$ HCN or HCO$^+$ observations to constrain $T_{\text{kin}}$, it will be important into account the convolution of $\epsilon$ ($n_{\text{H2}}$, $T_{\text{kin}}$) and $P(n_{\text{H2}})$ in order to constrain $T_{\text{kin}}$.

## 4.3. Dense Gas Conversion Factors

We emphasize differential measurements, because these may help control for an unknown absolute abundance scale. However, it will still often be of interest to discuss absolute masses or mass fractions of dense gas. To this end, we examine the dense gas “conversion factors” implied by our models (see Figure 15). These are formally the ratio $\langle\epsilon\rangle$ to $f_{\text{dense}}$, and give the factor by which line luminosity should be multiplied to calculate the dense gas mass.

Our derived conversion factors depend on the optical depth, abundance, $T_{\text{kin}}$, and the density distribution. This renders them substantially uncertain. The overall normalization that we find agrees with previous work for matched assumptions. However, we emphasize that arguments based on a fixed dynamical state, at a particular density inside a cloud, remain aggressive. Our calculations assume a fixed abundance for all four dense gas tracers. In this case at fixed $\tau$, $\alpha_{\text{HCO}^+} < \alpha_{\text{HNC}} \approx \alpha_{\text{HCN}} < \alpha_{\text{CS}}$.

## 4.4. Next Steps

We have presented a framework for interpreting mm-wave line ratio observations. A number of natural next directions recommend themselves.

**Observations:** First and foremost, we need a large set of observations targeting individual regions in other galaxies. These should include transitions that sample a range of densities, similar to the line suite considered in this paper. With such a large database, it will be possible to place empirical constraints on the shape of the high-density probability density function (PDF) as well as identify which tracers show coherent variations and which show signatures of strong abundance variations.

Ideally, these observations can be paired with high resolution imaging of ISM structure. Such observations have the prospect to constrain the mean density and Mach number at the scale of whole clouds (e.g., Leroy et al. 2016). Such external constraints on $n_0$ and $\sigma$ offer the prospect to validate the models proposed here. Then, if the models are validated, they offer the prospect to place priors on the model grid when it is applied to individual regions.

We have taken the optical depth to be known, or at least observable. However, the isotopologues of the dense gas tracers discussed here are very faint. Their emission is further suppressed by the differential excitation effects discussed here. As a result, constraints on the optical depth of the dense gas tracers are actually quite weak, to date. Systematic observation of H$^{13}$CN, H$^{13}$CO$^+$, H$^{13}$CS, etc. in the disks of normal galaxies will be essential to break a key degeneracy in the model.
Fortunately, current mm-wave facilities can make many of these key observations. ALMA, the GBT, the IRAM telescopes, and other modern facilities have the ability to survey dense gas tracers and cloud properties across nearby galaxies. Several recent, ongoing, and planned surveys promise to expand our knowledge of dense gas tracers across nearby galaxies (e.g., Usero et al. 2015; Bigiel et al. 2016; M. J. Gallagher et al. 2016, in preparation; and upcoming efforts using ARGUS on the GBT). Meanwhile, ALMA and NOEMA are rapidly expanding our knowledge of cloud-scale ISM structure, often in overlapping targets (e.g., Hughes et al. 2013a; Leroy et al. 2016).

Despite good prospects for the next few years, many key lines remain very faint. Systematic surveys of optically thin isotopologues (e.g., M. J. Jimenez Donaire et al., in preparation) or imaging of dense gas tracers in normal galaxies at high physical resolution (e.g., Rosolowsky et al. 2011) will remain very challenging in normal star-forming galaxies beyond the Local Group. Similarly, sensitivity considerations limit the suite of available lines for extragalactic studies compared to those used in the Galaxy. The long-term prospects for such studies may require an increase even beyond the capabilities of ALMA. One natural facility to enable such studies would be a Next Generation Very Large Array optimized for mm-wave studies. Large receiver arrays on large single-dish telescopes offer another promising direction.

**Power Law Tails:** A recurring result of our analysis is that the shape of the density distribution at high $n_{H2}$ exerts a large influence. The results for a power law tail and a pure lognormal distribution differ substantially. The main way to improve the interpretation of high-density tracers is to resolve the questions of how common a power law tail is, whether it has a universal slope, and the threshold density for its onset.

Most current extragalactic surveys target large areas of star-forming galaxies, often selected for local or global activity. We expect such locales to harbor self-gravitating gas that serves as the immediate fuel for star formation. In that sense, we do have a basic expectation that a power law tail should be present at some level in most star-forming regions. We do not, however, have a strong expectation for the universality of the slope of this tail or the density at which it becomes the dominant contributor to the PDF. Our power law models show an optimistic case, in which these factors do not vary much from beam to beam.

Galactic studies have made great progress in identifying power law tails in column density distributions and linking these to environment in the cloud and Galaxy (e.g., Abreu-Vicente et al. 2015; Stutz & Kainulainen 2015). From some combination of these studies, spectral surveys and highly resolved observations of the nearest galaxies, we can hope to understand a characteristic slope, typical range of variation, and condition for onset of a power law tail. For modeling distributions in other galaxies, an approach like that of Draine et al. (2007) may be most appropriate, in which the strength of a power law tail relative to the lognormal acts as a free parameter in the model.

We thank the anonymous referee for a constructive report that improved this paper. This work was carried out as part of the “Star Formation and Feedback in Nearby Galaxies” (SFNG) collaboration. The work of A.K.L. is partially supported by the National Science Foundation under grant

Nos. 1615105 and 1615109. A.U. acknowledges support from Spanish MINECO grants AYA2012-32295 and FIS2012-32096. J.M.D.K. gratefully acknowledges financial support in the form of an Emmy Noether Research Group from the Deutsche Forschungsgemeinschaft (DFG), grant number KR4801/1-1. A.H. acknowledges support from the Centre National d'Etudes Spatiales (CNES). The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. E.R. is supported by a Discovery Grant from NSERC of Canada. G.B. is supported by CONICYT/FONDECYT, Programa de Iniciacion, Folio 11150220. E.S. acknowledges financial support by the DAGAL network from the People Programme (Marie Curie Actions) of the European Union's Seventh Framework Programme FP7/2007-2013/under REA grant agreement number PITN-GA-2011-289313. A.U. acknowledges support from Spanish MINECO grants FIS2012-32096 and ESP2015-68964.

## References

Aalto, S., Martín, S., Costagliola, F., et al. 2015, A&A, 584, A42

Abreu-Vicente, J., Kainulainen, J., Stutz, A., Henning, T., & Beuther, H. 2015, A&A, 581, A74

André, P., Di Francesco, J., Ward-Thompson, D., et al. 2014, in Protostars and Planets VI, ed. H. Beuther et al. (Tucson, AZ: Univ. Arizona Press), 27

Berkhuijsen, E. M., & Fletcher, A. 2015, MNRAS, 448, 2469

Bigiel, F., Leroy, A. K., Jiménez-Donaire, M. J., et al. 2016, ApJ, 822, L26

Bolatto, A. D., Leroy, A. K., Rosolowsky, E., Walter, F., & Blitz, L. 2008, ApJ, 686, 948

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

Chen, H., Gao, Y., Braine, J., & Gu, Q. 2015, ApJ, 810, 140

Draine, B. T., Dale, D. A., Bendo, G., et al. 2007, ApJ, 663, 866

Elmegreen, B. G., & Scalo, J. 2004, ARA&A, 42, 211

Evans, N. J., II 1999, ARA&A, 37, 311

Evans, N. J., II, Heiderman, A., & Vutisalchavakul, N. 2014, ApJ, 782, 114

Federrath, C., & Klessen, R. S. 2012, ApJ, 761, 156

Federrath, C., & Klessen, R. S. 2013, ApJ, 763, 51

Feldmann, R., Gnedin, N. Y., & Kravtsov, A. V. 2011, ApJ, 732, 115

Gao, Y., & Solomon, P. M. 2004a, ApJS, 152, 63

Gao, Y., & Solomon, P. M. 2004b, ApJ, 606, 271

García-Burillo, S., Usero, A., Alonso-Herrero, A., et al. 2012, A&A, 539, A8

Glover, S. C. O., & Clark, P. C. 2012, MNRAS, 426, 377

Goldreich, P., & Kwan, J. 1974, ApJ, 189, 441

Heiderman, A., Evans, N. J., II, Allen, L. E., Huard, T., & Heyer, M. 2010, ApJ, 723, 1019

Hennebelle, P., & Chabrier, G. 2011, ApJ, 743, L29

Heyer, M., & Dame, T. M. 2015, ARA&A, 53, 583

Hopkins, P. F., Narayanan, D., Murray, N., & Quataert, E. 2013, MNRAS, 433, 69

Hughes, A., Meidt, S. E., Colombo, D., et al. 2013a, ApJ, 779, 46

Hughes, A., Meidt, S. E., Schinnerer, E., et al. 2013b, ApJ, 779, 44

Jiménez-Donaire, M. J., Cormier, D., Bigiel, F., et al. 2017, MNRAS, 466, 49

Kainulainen, J., Beuther, H., Henning, T., & Plume, R. 2009, A&A, 508, L35

Kainulainen, J., Federrath, C., & Henning, T. 2013, A&A, 553, L8

Kepley, A. A., Leroy, A. K., Frayer, D., et al. 2014, ApJ, 780, L13

Kruijssen, J. M. D. 2012, MNRAS, 426, 3008

Kruijssen, J. M. D., Longmore, S. N., Elmegreen, B. G., et al. 2014, MNRAS, 440, 3370

Krumholz, M. R. 2014, MNRAS, 437, 1662

Krumholz, M. R., Dekel, A., & McKee, C. F. 2012, ApJ, 745, 69

Krumholz, M. R., & Kruijssen, J. M. D. 2015, MNRAS, 453, 739

Krumholz, M. R., Kruijssen, J. M. D., & Crocker, R. M. 2016, arXiv:1605.02850

Krumholz, M. R., & McKee, C. F. 2005, ApJ, 630, 250

Krumholz, M. R., & Thompson, T. A. 2007, ApJ, 669, 289

Lada, C. J., Forbrich, J., Lombardi, M., & Alves, J. F. 2012, ApJ, 745, 190

Lada, C. J., & Lada, E. A. 2003, ARA&A, 41, 57

Lada, C. J., Lombardi, M., & Alves, J. F. 2010, ApJ, 724, 687

Le Bourlot, J., Le Petit, F., Pinto, C., Roueff, E., & Roy, F. 2012, A&A, 541, A76

Leroy, A. K., Bolatto, A., Bot, C., et al. 2009, ApJ, 702, 352

Leroy, A. K., Bolatto, A. D., Ostriker, E. C., et al. 2015, ApJ, 801, 25
Leroy, A. K., Hughes, A., Schruba, A., et al. 2016, ApJ, 831, 16

Lombardi, M., Alves, J., & Lada, C. J. 2015, A&A, 576, L1

Longmore, S. N., Bally, J., Testi, L., et al. 2013, MNRAS, 429, 987

Mac Low, M.-M., & Klessen, R. S. 2004, RvMP, 76, 125

Maloney, P., & Black, J. H. 1988, ApJ, 325, 389

Mangum, J. G., & Shirley, Y. L. 2015, PASP, 127, 266

Martín, S., Mauersberger, R., Martín-Pintado, J., Henkel, C., & García-Burillo, S. 2006, ApJS, 164, 450

McKee, C. F., & Ostriker, E. C. 2007, ARA&A, 45, 565

Meier, D. S., Walter, F., Bolatto, A. D., et al. 2015, ApJ, 801, 63

Narayanan, D., Cox, T. J., Shirley, Y., et al. 2008, ApJ, 684, 996

Narayanan, D., Krumholz, M. R., Ostriker, E. C., & Hernquist, L. 2012, MNRAS, 421, 3127

Padoan, P., & Nordlund, Å. 2002, ApJ, 576, 870

Padoan, P., & Nordlund, Å. 2011, ApJ, 730, 40

Rathborne, J. M., Longmore, S. N., Jackson, J. M., et al. 2014, ApJ, 795, L25

Rosolowsky, E., Pineda, J. E., & Gao, Y. 2011, MNRAS, 415, 1977

Rosolowsky, E. W., Pineda, J. E., Foster, J. B., et al. 2008, ApJS, 175, 509

Schinnerer, E., Meidt, S. E., Pety, J., et al. 2013, ApJ, 779, 42

Schneider, N., Bontemps, S., Girichidis, P., et al. 2015, MNRAS, 453, 41

Schöier, F. L., van der Tak, F. F. S., van Dishoeck, E. F., & Black, J. H. 2005, A&A, 432, 369

Scoville, N., Sheth, K., Walter, F., et al. 2015, ApJ, 800, 70

Scoville, N. Z., & Solomon, P. M. 1974, ApJ, 187, L67

Shirley, Y. L. 2015, PASP, 127, 299

Sobolev, V. V. 1960, Moving Envelopes of Stars (Cambridge: Harvard Univ. Press)

Stutz, A. M., & Kainulainen, J. 2015, A&A, 577, L6

Usero, A., Leroy, A. K., Walter, F., et al. 2015, AJ, 150, 115

Vallini, L., Ferrara, A., Pallottini, A., & Gallerani, S. 2016, arXiv:1606.08464

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

Vazquez-Semadeni, E. 1994, ApJ, 423, 681

Wilson, C. D., Warren, B. E., Israel, F. P., et al. 2012, MNRAS, 424, 3050

Wu, J., Evans, N. J., II, Gao, Y., et al. 2005, ApJ, 635, L173

Wu, J., Evans, N. J., II, Shirley, Y. L., & Knez, C. 2010, ApJS, 188, 313