Does the HCN/CO Ratio Trace the Star-forming Fraction of Gas? I. A Comparison with Analytical Models of Star Formation

A current challenge for understanding star formation in molecular clouds is determining the fraction of gas that is converted into stars over the lifetime of a cloud. Observations show that star formation sites primarily lie in regions of dense molecular gas in Milky Way clouds (Lada et al. 1991a, 1991b; Helfer & Blitz 1997a, 1997b), and these regions are confined to ∼0.1 pc scales within larger molecular structures in the form of clumps or filaments (André et al. 2016). Within these gas structures, the fraction of gravitationally bound gas continues to form stars. Analytical models of star formation rely on estimates of the self-gravitating gas fraction, f_grav, to then predict the star formation rate (SFR; e.g., Krumholz & McKee 2005; Hennebelle & Chabrier 2011; Padoan & Nordlund 2011; Federrath & Klessen 2012; Burkhart & Mocz 2019), which makes f_grav an important parameter to constrain observationally.

Extragalactic observations rely on molecular transitions with high critical densities, n_crit ≳ 10⁴ cm⁻³, to gain information on the dense gas in other galaxies. The most commonly used dense molecular gas tracer in extragalactic studies is the HCN J = 1–0 transition (Gao & Solomon 2004a, 2004b). Under the common assumption that the total emissivity of HCN traces the dense gas mass, ${I}_{\mathrm{HCN}}\propto {{\rm{\Sigma }}}_{\mathrm{dense}}$, the ratio of the HCN and CO emissivities, ${I}_{\mathrm{HCN}}/{I}_{\mathrm{CO}}$, is proportional to the fraction of molecular gas in the dense phase, f_dense. If the dense gas mass traced by HCN is also self-gravitating, then this line ratio is a simple observational method for estimating f_grav. However, the interstellar medium (ISM) of galaxies resides at a range of densities ≲1 − ≫ 10⁸ cm⁻³, and molecular transitions are sensitive to a continuum of these densities, including some fraction below their critical density (Shirley 2015; Leroy et al. 2017a). Recent studies within the Milky Way have also shown that HCN may predominantly trace moderate gas densities (Kauffmann et al. 2017), rather than denser gas associated with star formation. The fraction of CO-dark or CO-faint gas may also contribute to variations in ${I}_{\mathrm{HCN}}/{I}_{\mathrm{CO}}$ (Grenier et al. 2005; Wolfire et al. 2010; Bolatto et al. 2013). On average in the Milky Way, ∼30% of molecular hydrogen appears to be CO-faint (Pineda et al. 2013; Langer et al. 2014). This fraction is higher in lower-metallicity gas (Pineda et al. 2013; Jameson et al. 2018; Chevance et al. 2020) and appears to be more significant for clouds at lower masses (Grenier et al. 2005; Bolatto et al. 2013) and at higher radii in the Milky Way (Pineda et al. 2013; Langer et al. 2014; Chevance et al. 2020). This would likely result in overestimates of f_dense in lower-mass clouds, assuming HCN emission is not affected in the same way. This is likely less significant in more extreme systems (e.g., luminous IR galaxies, LIRGs, or ultraluminous IR galaxies, ULIRGs) or galaxy centers, such as those studied in this paper.

The dynamics of the gas in the ISM are set by a combination of gravity, turbulence, and magnetic fields, and these processes act together to set the spatial structure of gas in ISM clouds. These processes also likely play a role in setting the star formation properties of the ISM. One type of analytical model, gravoturbulent models of star formation, aims to predict both the observed structure of gas clouds and their star formation properties. In particular, gravoturbulent models of star formation predict the shape of the gas volume density probability distribution function (n − PDF) and the star formation efficiency (_ff) over a freefall time (t_ff) ³ (Krumholz & McKee 2005; Hennebelle & Chabrier 2011; Padoan & Nordlund 2011; Federrath & Klessen 2012; Imara & Burkhart 2016; Burkhart 2018). ⁴ Studies have found that the column density PDF of the diffuse (n < 1 cm⁻³) component of gas in the Milky Way and M33 is consistent with a lognormal PDF (see Berkhuijsen & Fletcher 2008; Hill et al. 2008; Tabatabaei et al. 2008; Burkhart et al. 2015). The seminal analytical work by Vazquez-Semadeni (1994) showed that if the turbulent ISM develops a series of isothermal and interacting supersonic shocks, the gas would naturally follow a lognormal PDF (see Vazquez-Semadeni 1994; Padoan et al. 1997; Scalo et al. 1998; Nordlund & Padoan 1999.) In this picture, the shocks amplify each other via a turbulent cascade of energy, and this multiplicative process causes the gas density PDF to assume a lognormal shape (see Vazquez-Semadeni 1994; Padoan et al. 1997; Scalo et al. 1998; Nordlund & Padoan 1999; Chen et al. 2018).

Observations of molecular regions of the ISM reveal that the gas column density PDF takes on a different form at high densities. The highest-density regions within more evolved molecular clouds contribute a power-law tail to the gas column density PDF (see Chen et al. 2018), with some cloud PDFs being almost entirely power law (e.g., Kainulainen et al. 2009; Schneider et al. 2013; Lombardi et al. 2015; Schneider et al. 2015, 2016; Alves et al. 2017). This power law has also been observed in simulations that develop self-gravitating gas (see Ballesteros-Paredes et al. 2011; Collins et al. 2012; Schneider et al. 2015; Burkhart et al. 2017; Padoan et al. 2017.) These results strongly suggest that the gas density PDF in a star-forming molecular cloud is likely a combination of a lognormal and power-law shape, and that the power-law tail is potentially the result of gas becoming self-gravitating. The fraction of gas within this power-law tail would then be the self-gravitating gas fraction, f_grav. The density at which the n − PDF transitions from a lognormal to power-law shape would then represent the moment at which gas becomes available to star formation, n_SF. The fraction of dense gas mass above n_SF (f_grav) relative to the total mass of a star-forming cloud is then

$$\begin{eqnarray}&&{f}_{\mathrm{grav}}=\displaystyle \frac{M(n\gt {n}_{\mathrm{SF}})}{M}.\end{eqnarray} \tag{ 1 }$$

Turbulent models of star formation estimate _ff by integrating over a purely lognormal gas density PDF, also above a gas threshold density (see Krumholz & McKee 2005; Hennebelle & Chabrier 2011; Padoan & Nordlund 2011; Federrath & Klessen 2012). This threshold density is also meant to capture the point at which gas becomes self-gravitating in the ISM, so that the fraction of gas above this threshold is f_grav in these models. However, lognormal-only models fail to explain observed variations in _ff and Mach number, ${ \mathcal M }$, seen in some galaxies (e.g., Leroy et al. 2017b). For example, giant molecular clouds (GMCs) in M51 show a weak anticorrelation between velocity dispersion (which is proportional to ${ \mathcal M }$; see Equation (11)) and the star formation efficiency of gas per freefall time, _ff. ⁵ Leroy et al. (2017a) argue that this anticorrelation may reflect differences in the dynamical state of their clouds with galactocentric radius.

Burkhart & Mocz (2019) show that including a power-law tail in the gas volume density PDF reproduces the observed variations in _ff in M51, without requiring changes in the dynamical state of the clouds and without explicitly setting n_SF. They find a slight anticorrelation between _ff and ${ \mathcal M }$ for virialized clouds (α_vir ≈ 1). This anticorrelation coincides with an increasing depletion time with ${ \mathcal M }$, in agreement with the findings from PAWS (Leroy et al. 2017b). _ff ∝ f_grav implies f_grav may also anticorrelate with ${ \mathcal M }$. Thus, without needing to explicitly set n_SF, it is ultimately a decrease in f_grav with respect to increasing ${ \mathcal M }$ that leads to the anticorrelation between σ_rv and _ff (Burkhart & Mocz 2019). However, in starbursts where ${ \mathcal M }$ is higher, we see higher _ff and shorter t_dep on average (Wilson et al. 2019). Starbursts typically have enhanced HCN/CO ratios in addition to a shorter depletion time t_dep (Kennicutt & De Los Reyes 2021). If ${I}_{\mathrm{HCN}}/{I}_{\mathrm{CO}}\propto {f}_{\mathrm{dense}}$ and f_dense ∝ f_grav, then the result from BM19 appears to differ from what is observed in starbursts.

A potential explanation for these differences may be differences in the timescale for star formation, which may be set by the environment in which a gas cloud is immersed. The star formation law can be written as (Krumholz et al. 2012)

$$\begin{eqnarray}&&{t}_{\mathrm{ff}}{{\rm{\Sigma }}}_{\mathrm{SFR}}={\epsilon }_{\mathrm{ff}}{{\rm{\Sigma }}}_{\mathrm{gas}},\end{eqnarray} \tag{ 2 }$$

where Σ_gas is the gas surface density, Σ_SFR is the SFR surface density, and t_ff is the freefall time and is set by the self-gravity of a cloud.

Burkhart & Mocz (2019) demonstrate the connection between f_grav and the instantaneous efficiency of the gas, _inst ≈ ₀ f_grav, which reflects both the local efficiency, ₀ (set, e.g., by stellar feedback), and f_grav. This local efficiency may correlate with its observational analog, the star formation efficiency per freefall time, _ff (Krumholz & McKee 2005; Lee et al. 2016). Furthermore, if turbulence plays a significant role in setting n_SF, then we may find a correlation between _ff, observed velocity dispersions of gas, and turbulent pressure. Turbulent models of the ISM predict a dependence of ${\sigma }_{{\rm{n}}/{{\rm{n}}}_{0}}$, the density variance of the volume density PDF (n − PDF) on the sonic Mach number, ${ \mathcal M }={\sigma }_{v,3D}/{c}_{s}$, ⁶ within individual star-forming clouds. Lada et al. (1994) found a correlation between the extinction, A_V, in the dark cloud IC 5126 and the standard deviation of A_V, with extinction increasing with dispersion (Goodman et al. 2009). Kainulainen & Tan (2013) find a correlation between measurements of the velocity dispersion from ¹²CO and¹³CO and the density contrast (N/N₀) of the column density PDFs (N − PDFs) derived from IR data in several Milky Way clouds. Combined, these correlations imply that the CO velocity dispersion may be sensitive to the density variance of N − PDF, ${\sigma }_{{\rm{N}}/{{\rm{N}}}_{0}}$ (where ${\sigma }_{N/{N}_{0}}\propto \mathrm{ln}(N/{N}_{0})$), and therefore is a probe of the ISM physics.

In external galaxies where resolution is limited, molecular line ratios are an additional tool for assessing the n − PDF shape. Shirley (2015) showed that molecular transitions have an emissivity ⁷ that extends over a range of gas densities, including a significant amount of emission at densities below the critical density associated with that transition. ⁸ To determine whether molecular line ratios remain sensitive to n − PDF shape, Leroy et al. (2017a) modeled molecular line emissivities and explored a range of n − PDF shapes. They reported that dense gas tracers (e.g., HCN and HCO⁺) are more sensitive to changes in the shape of the n − PDF than lower-density tracers such as CO. Combining molecular line ratios of dense gas tracers with information on kinematics is therefore a promising tool for assessing the n − PDF information of clouds in external galaxies.

To assess the relation between the HCN/CO ratio and f_grav, we examine more extreme star-forming environments in which turbulence is stronger (e.g., mergers, starbursts, (U)LIRGs, or barred galaxies). We study a sample of 10 (U)LIRGs and disk galaxy centers that have archival CO, HCN, CN, and HCO⁺ J = 1–0 data, in addition to 93 GHz radio continuum. In this paper, we focus on general trends of the HCN/CO ratio, the SFR surface density, and _ff, and we compare these modeled trends with the observed trends in our sample. We use the EMPIRE sample of galaxies (Jiménez-Donaire et al. 2019) as a comparison, which predominantly targets normal regions of star formation within galaxy disks.

Our sample consists of 10 nearby (z < 0.03) galaxies, including the dense centers of five disk galaxies and five mergers and (U)LIRGs. We list these galaxies and their basic properties in Table 1. Four of the five disk galaxies in our sample are also barred. For each galaxy, we image archival ALMA data of the HCN, CN, CO, and HCO⁺ J = 1–0 transitions, in addition to the radio continuum emission at 93 GHz. The data are uv-matched and tapered to a common beam for each individual galaxy. The sample selection and data reduction process are presented in Wilson et al. (2023) in detail, except for NGC 4038/9, NGC 1808, and NGC 3351, which have been reduced and imaged separately at higher velocity resolutions to be included in this analysis.

2.1. Moment Maps

2.2. Molecular Gas Surface Densities

One of the main goals of this study is to assess variations in the fraction of gas in the dense phase as traced by the HCN and CO molecular line luminosities. Our analysis allows us to estimate trends in the HCN and CO luminosity-to-mass conversion factors (see Bolatto et al. 2013), which we define as

$$\begin{eqnarray}&&{\alpha }_{\mathrm{mol}}=\displaystyle \frac{{{\rm{\Sigma }}}_{\mathrm{mol}}}{{I}_{\mathrm{mol}}}\qquad [{M}_{\odot }\,{\mathrm{pc}}^{2}\,{\left({\rm{K}}\ \mathrm{km}\ {{\rm{s}}}^{-1}\right)}^{-1}],\end{eqnarray} \tag{ 3 }$$

where Σ_mol is the mass surface density of the molecular gas, including helium, I_mol = L_mol/A_pix is the intensity in units of K km s⁻¹, and L_mol is total luminosity over the physical area of a pixel, A_pix. To calculate the molecular gas mass surface density, Σ_mol, we use

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\mathrm{mol}}={\alpha }_{\mathrm{mol}}\ {I}_{\mathrm{mol}}\ \cos (i).\end{eqnarray} \tag{ 4 }$$

We apply inclination angles only to the disk galaxies in this sample, as inclination angles are typically uncertain in mergers and (U)LIRGs. Thus, the uncorrected measurements of galaxies with nonzero inclinations will result in overestimates of Σ_mol and other surface densities.

To facilitate comparison with other studies, our fiducial value is α_CO = 1.1 [M_⊙ (K km s⁻¹ pc²)⁻¹] for our sample of galaxies, the (U)LIRG value including helium (Downes et al. 1993), which is about four times lower than the Milky Way value. This lower value is motivated by evidence that gas in these systems is subject to more extreme excitation mechanisms, e.g., higher temperatures and densities (see Downes et al. 1993; Bolatto et al. 2013 and references therein). Additionally, the gas traced by CO in these systems often shows broad line widths, potentially reducing the opacity of the CO transition (Downes et al. 1993; Bolatto et al. 2013). Downes et al. (1993) also suggest that CO may be subthermally excited (T_ex < T_kin) in starbursts and (U)LIRGs. We choose a fixed value of α_CO = 4.35 [M_⊙ (K km s⁻¹ pc²)⁻¹] for the EMPIRE sample of galaxies because these are mostly disk galaxies and are likely more similar to the Milky Way than to starbursts. Because we already have NGC 3627 in our sample, we drop NGC 3627 from the EMPIRE data.

The HCN conversion factor is less certain. Historically, ${\alpha }_{\mathrm{HCN}}\approx 13.6$ [M_⊙ (K km s⁻¹ pc²)⁻¹] has been used (Gao & Solomon 2004a, 2004b), which is appropriate for a virialized cloud core with a mean density ${n}_{{{\rm{H}}}_{2}}\sim 3\times {10}^{4}$ cm⁻³ and brightness temperature T_B ∼ 35 K (e.g., Radford et al. 1991, including helium). This HCN conversion factor assumes that this molecular transition is optically thick and that the gas it traces is in local thermodynamic equilibrium (LTE). In Bemis (2020) and Bemis & Wilson (in prep.), we explore variations in the relative values of ${\alpha }_{\mathrm{HCN}}$ and α_CO using an nonLTE radiative transfer analysis. If all of the above assumptions are true, then the fraction of dense gas traced by HCN/CO is given by

$$\begin{eqnarray}&&{f}_{\mathrm{dense}}=\displaystyle \frac{{\alpha }_{\mathrm{HCN}}}{{\alpha }_{\mathrm{CO}}}\displaystyle \frac{{I}_{\mathrm{HCN}}}{{I}_{\mathrm{CO}}},\end{eqnarray} \tag{ 5 }$$

with ${\alpha }_{\mathrm{HCN}}/{\alpha }_{\mathrm{CO}}=3.2$ implied from the above discussion. We do this assuming that the physical effects impacting α_CO will impact ${\alpha }_{\mathrm{HCN}}$ in a similar way. Simulations show that a slight reduction in ${\alpha }_{\mathrm{HCN}}$ does happens in presence of higher gas temperatures (T_kin = 20 K; see Onus et al. 2018), similar to our expectation for α_CO. Furthermore, by using ${\alpha }_{\mathrm{HCN}}/{\alpha }_{\mathrm{CO}}=3.2$, our plots of f_dense are a direct reflection of the HCN/CO ratio. We distinguish this calculation of f_dense from analytical estimates of the fraction of dense star-forming gas, which is generally defined as the mass fraction above the density at which gas becomes self-gravitating, f_grav.

2.3. The Radio Continuum SFR

We detect radio continuum emission at 93 GHz in all of our sources and use this as our SFR tracer. The radio continuum is a combination of thermal (T) free–free emission and nonthermal (NT) synchrotron emission from regions with massive star formation, spanning ∼1–100 GHz. At 93 GHz, we are in the regime where thermal free–free emission from young star-forming regions will likely dominate the radio continuum emission. Nonthermal emission is expected to contribute ∼25% to the radio continuum luminosity at this frequency (see Murphy et al. 2011; Wilson et al. 2019), assuming an electron temperature T_e ∼ 10⁴ K and nonthermal spectral index α_NT ∼ 0.84. This fraction will change if there are variations in either T_e or α_NT. We adopt fixed values for α_NT = 0.84 and T_e = 10⁴ K, and we use the composite calibration from Murphy et al. (2011), which accounts for both thermal and nonthermal contributions to the SFR,

$$\begin{eqnarray}&&\begin{array}{l}\left(\displaystyle \frac{{\mathrm{SFR}}_{\nu }}{{M}_{\odot }{\mathrm{yr}}^{-1}}\right)\\ \quad =\,{10}^{-27}\left[2.18{\left(\displaystyle \frac{{T}_{{\rm{e}}}}{{10}^{4}\ {\rm{K}}}\right)}^{0.45}{\left(\displaystyle \frac{\nu }{\mathrm{GHz}}\right)}^{-0.1}\right.\\ \qquad {\left.+15.1{\left(\displaystyle \frac{\nu }{\mathrm{GHz}}\right)}^{-{\alpha }^{\mathrm{NT}}}\right]}^{-1}\left(\displaystyle \frac{{L}_{\nu }}{\mathrm{erg}\,{{\rm{s}}}^{-1}\ {\mathrm{Hz}}^{-1}}\right).\end{array}\end{eqnarray} \tag{ 6 }$$

At 93 GHz, the radio continuum emission from star-forming regions may overlap with the lower-frequency tail of the dust SED. Wilson et al. (2019) estimate a ∼10% contribution from dust at 93 GHz for IRAS 13120-5453, NGC 3256, and NGC 7469, the three most IR-luminous galaxies in our sample. We adopt this correction factor of ∼10% for all sources, but we acknowledge that the emissivity of dust at 93 GHz may vary between our sources. We further mask pixels that may be contaminated with emission from an active galactic nucleus (AGN; see Table 1).

Variations in T_e and α_NT will impact our SFR estimates from the radio continuum at 93 GHz. Electron temperatures typical of HII regions are T_e ∼ 5 × 10³ − 10⁴ K, which will produce a change of ∼30% (Murphy et al. 2011). In contrast to this, α_NT has been observed as low as ∼0.5, which would also give a change in luminosity ∼30% at 93 GHz. Wilson et al. (2019) also find evidence of a significant fraction (up to 50%) of nonthermal emission in NGC 7469 at 93 GHz by comparing with an archival radio continuum map at 8 GHz. We adopt an uncertainty of 30% for SFR estimates derived from the radio continuum.

We compare the results of our sample with those of the EMPIRE survey (Jiménez-Donaire et al. 2019), for which there are publicly available single-dish (IRAM 30m) observations of HCN and CO. We estimate Σ_SFR in EMPIRE galaxies using 24 μm IR maps from the Spitzer Space Telescope. These IR data are convolved to a 15'' Gaussian beam using Aniano et al. (2011) Gaussian kernels and further smoothed to a 33'' Gaussian beam using CASA (McMullin et al. 2007). Backgrounds are subtracted and SFRs are derived using the Rieke et al. (2009) calibration.

2.3.1. Star Formation Timescales and Efficiency

We use the SFR and molecular gas surface densities to estimate the depletion time of the total (mol = CO) and dense (mol = HCN) molecular gas content,

$$\begin{eqnarray}&&{t}_{\mathrm{dep}}=\displaystyle \frac{{{\rm{\Sigma }}}_{\mathrm{mol}}}{{{\rm{\Sigma }}}_{\mathrm{SFR}}},\end{eqnarray} \tag{ 7 }$$

where Σ_mol is estimated using Equation (4). To estimate the dimensionless star formation efficiency, we compare this depletion timescale with the freefall timescale,

$$\begin{eqnarray}&&{t}_{\mathrm{ff}}=\sqrt{\displaystyle \frac{3\pi }{32G\rho }},\end{eqnarray} \tag{ 8 }$$

where ρ is a characteristic density of the gas associated with star formation. We calculate gas density assuming a fixed line-of-sight (LOS) depth via ρ ≈ Σ_mol/R, where R is equivalent to half of the cloud or gas depth along the LOS. We do not know R, so we assume a fixed value of 100 pc for all galaxies when estimating ρ. We also explored estimating the freefall time using the Krumholz et al. (2012) prescriptions for t_ff in the GMC and Toomre regimes. We find little qualitative difference in the results between using t_ff as described above and the Krumholz et al. (2012) prescription.

The efficiency of the star formation process is then estimated by comparing the observed depletion timescales with estimates of the freefall time,

$$\begin{eqnarray}&&{\epsilon }_{\mathrm{ff}}=\displaystyle \frac{{t}_{\mathrm{ff}}}{{t}_{\mathrm{dep}}}\mathrm{,}\end{eqnarray} \tag{ 9 }$$

which is just the star formation law (Equation (2)) rewritten in terms of depletion time.

2.4. Velocity Dispersion

We measure the one-dimensional velocity dispersion of the molecular gas in our sources, σ_v , using the CO J = 1–0 transition. The velocity dispersions, σ_v,meas, are measured directly from moment 2 maps. We correct for broadening of the line due to the finite spectral resolution of our data using Rosolowsky & Leroy (2006),

$$\begin{eqnarray}&&{\sigma }_{v}=\sqrt{{\sigma }_{v\mathrm{meas}}^{2}-\displaystyle \frac{\delta {v}^{2}}{2\pi }},\end{eqnarray} \tag{ 10 }$$

where δ v is the channel resolution at which we image the data. This value is then converted into the three-dimensional velocity dispersion via ${\sigma }_{v,3D}=\sqrt{3}{\sigma }_{v}$.

2.4.1. Velocity Dispersion as a Tracer of Mach Number

As previous studies have done for Milky Way clouds (Kainulainen & Federrath 2017), we used CO line widths as an indicator of the Mach number of the gas in our galaxies,

$$\begin{eqnarray}&&{ \mathcal M }=\displaystyle \frac{\sqrt{3}{\sigma }_{{\rm{v}},1{\rm{D}}}}{{c}_{s}},\end{eqnarray} \tag{ 11 }$$

where c_s is the thermal sound speed, and σ_v,1D is the observed one-dimensional velocity dispersion. There are limits to this approach, and the ability of the CO line to trace cloud turbulence may be limited by its optical depth (Goodman et al. 2009; Burkhart et al. 2013), and the observed line width can include disk rotation or other large-scale motions, but alternative measures of gas kinematics in extragalactic clouds are lacking.

To estimate c_s , we consider an ideal gas equation of state (P_th = nk_B T/μ m_H), such that the sound speed is ${c}_{s}\,=\sqrt{{k}_{{\rm{B}}}T/\mu {m}_{{\rm{H}}}}$, where k_B, T, μ, and m_H are the Boltzmann constant, gas kinetic temperature, the mass of a hydrogen atom, and the mean particle weight. For a gas in which hydrogen is primarily in molecular form, μ = 2.33 (assuming cosmic abundances; Kauffmann et al. 2008). For μ = 2.33 and a temperature range of T = 10–100 K, c_s ≈ 0.2–0.6 km s⁻¹. The choice of c_s can have a significant impact on the estimate of ${ \mathcal M }$. We choose an intermediate value c_s = 0.4 km s⁻¹, corresponding to T ∼ 45 K, as our fiducial value.

2.5. Uncertainties from a Multiscale Sample

The spatial resolution of the data in our sample spans ∼50–900 pc. We consider how this range may impact our measurements from the context of a turbulent ISM, and how we can best interpret the measured σ_v and line ratios in our sources. The beam-filling fraction, ϕ_ff,beam is less of an uncertainty in this study because the sources with the lowest resolution are (U)LIRGs and likely have filling fractions approaching unity, while the spiral galaxy centers have resolutions approaching cloud scales. This reduces the issue of variations in the beam-filling fractions from galaxy to galaxy. A significant source of uncertainty instead comes from variations in the relative filling fractions of the HCN/CO transitions,

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{HCN}/\mathrm{CO}}\approx \displaystyle \frac{{r}_{\mathrm{HCN}}^{2}}{{r}_{\mathrm{CO}}^{2}},\end{eqnarray} \tag{ 12 }$$

where r_HCN and r_CO are the average radial extents of HCN and CO, respectively. We explore variations in the filling fraction in Bemis (2020) and Bemis & Wilson (in prep.) from a radiative transfer perspective.

In addition to the uncertainty in ${{\rm{\Phi }}}_{\mathrm{HCN}/\mathrm{CO}}$, the velocity dispersion in a turbulent medium is scale dependent (ℓ) such that (Larson 1981; Heyer et al. 2009)

$$\begin{eqnarray}&&{\sigma }_{v}({\ell })={\sigma }_{v,L}{\left(\displaystyle \frac{{\ell }}{L}\right)}^{p}\mathrm{,}\end{eqnarray} \tag{ 13 }$$

where p ∼ 0.33–0.5 depends on the type of turbulence (e.g., p ≈ 0.33 for Kolmogorov, Larson 1981), L can be defined as the cloud diameter (or turbulent injection scale), and σ_v,L is the cloud dispersion at that scale. A sample of galaxies at different physical scales will therefore be affected by this scale dependence, such that velocity dispersions will be lower at smaller scales and vice versa. However, this relation should saturate at the turbulent injection scale, which may mitigate some of this uncertainty. If larger-scale cloud-cloud motions do not affect the velocity dispersions of CO, then this saturation should be detected, although there is still an uncertainty in L. This means that lower-resolution observations may still be useful for assessing σ_v.

We consider several prescriptions of f_grav in the context of analytical models of star formation and compare the predictions of these models to measurements of the ${I}_{\mathrm{HCN}}/{I}_{\mathrm{CO}}$ ratio and the star formation properties of our galaxies. We use Equation (2), which is relevant to gravoturbulent models of star formation (KMD12).

_ff is calculated by integrating over the star-forming portion of the density-weighted gas density PDF (e.g., Krumholz & McKee 2005; Hennebelle & Chabrier 2011; Padoan & Nordlund 2011; Federrath & Klessen 2012),

$$\begin{eqnarray}&&{\epsilon }_{\mathrm{ff}}={\epsilon }_{0}{\int }_{{n}_{\mathrm{SF}}}^{\infty }\displaystyle \frac{{t}_{\mathrm{ff}}({n}_{0})}{{t}_{\mathrm{ff}}(n)}\displaystyle \frac{n}{{n}_{0}}p(n){\rm{d}}n.\end{eqnarray} \tag{ 14 }$$

Here p(n) is the volumetric density PDF of the cloud, n − PDF, and n₀ is its mean density. The n − PDF characterizes the distribution of densities of the entire gas within the cloud, including the diffuse atomic and molecular components. We briefly review specific terms within this equation individually.

p(n): Isothermal gas in the presence of supersonic turbulence naturally becomes distributed such that its n − PDF is roughly lognormal (Vazquez-Semadeni 1994; Nordlund & Padoan 1999; Wada & Norman 2001). Earlier formalisms of gravoturbulent models adopt a purely lognormal form of p(n) (e.g., Krumholz & McKee 2005; Padoan & Nordlund 2011; Federrath & Klessen 2012). In terms of the logarithmic volume density $s=\mathrm{ln}(n/{n}_{0})$, this is given by

$$\begin{eqnarray}&&{{\rm{p}}}_{s}=N\displaystyle \frac{1}{\sqrt{2\pi {\sigma }_{s}^{2}}}\exp \left(-\displaystyle \frac{{\left(s-{s}_{0}\right)}^{2}}{2{\sigma }_{s}^{2}}\right),\end{eqnarray} \tag{ 15 }$$

where N is a normalization constant, ⁹ ${\sigma }_{s}^{2}$ is the logarithmic density variance, which depends on the underlying physics of the gas and sets the width of p(n), and ${s}_{0}=-0.5\,{\sigma }_{s}^{2}$. Newer formalisms (i.e., Burkhart 2018) of gravoturbulent models of star formation suggest that the n − PDF may then evolve to include a high-density power-law tail that contains gas that is becoming gravitationally unstable to collapse. Observations of clouds in the Milky Way support a composite form of p(n), where a power-law tail is clearly present in addition to a more diffuse component of gas. Burkhart (2018) presents this as a piecewise function,

$$\begin{eqnarray}{{\rm{p}}}_{s}=\left\{\begin{array}{ll}N\tfrac{1}{\sqrt{2\pi {\sigma }_{s}^{2}}}\exp \left(-\tfrac{{\left(s-{s}_{0}\right)}^{2}}{2{\sigma }_{s}^{2}}\right), & s\lt {s}_{t}\\ {{NCe}}^{-{\alpha }_{\mathrm{PL}}s}, & s\gt {s}_{t}\end{array}\right.\end{eqnarray} \tag{ 16 }$$

where s_t is the logarithmic form of the transition density between the two components of the n − PDF, and α_PL is the slope of the power-law tail. The factors N and C are normalization constants given by Burkhart (2018).

σ_s : Numerical studies have shown that if the turbulence is super-Alfvénic and the magnetic field (B) follows a power-law relation with gas density, B ∝ n^1/2, then the logarithmic density variance is given by (Molina et al. 2012; Federrath & Klessen 2012)

$$\begin{eqnarray}&&{\sigma }_{s}^{2}=\mathrm{ln}\left(1+{b}^{2}{{ \mathcal M }}^{2}\displaystyle \frac{\beta }{\beta +1}\right)\mathrm{,}\end{eqnarray} \tag{ 17 }$$

where ${ \mathcal M }$ is the Mach number, β = P_th/P_mag is the plasma beta (which characterizes the ratio of thermal pressure to magnetic pressure), and b is the turbulent forcing parameter (which characterizes the relative amount of solenoidal b = 1 or compressive b = 0.33 turbulence within the gas). Without a clear prescription for β and b in our sample, we simply take β → ∞ , which assumes that B = 0 G. We take an intermediate value b = 0.4, which assumes that turbulence is a mixture of compressive and solenoidal forcing. Nonisothermal gas will add intermittency, resulting in deviations from a lognormal shape of p(n) (Federrath & Banerjee 2015). For the purpose of this work, we assume that an underlying lognormal is a reasonable intermediate approximation to part of the n − PDF (Federrath & Banerjee 2015).

n_SF: Gravity overcomes supportive processes and begins the process of collapse above the density threshold, n_SF. The processes competing with gravity include some combination of magnetic support, internal turbulent motions, and external ISM pressures. Any gas with densities higher than n_SF is then potentially gravitationally unstable and star forming, so that n_SF serves as the lower limit of the integration that determines _ff.

t_ff(n₀)/t_ff(n): Within the integral in Equation (14), the integral is a freefall time factor, which converts the integrand into a dimensionless mass rate equivalent to the mass per freefall time (Federrath & Klessen 2012). This time factor is treated differently depending on the analytical model being considered, e.g., single-freefall (SFF) time (Krumholz & McKee 2005; Padoan & Nordlund 2011) versus multi-freefall (MFF) time models (Hennebelle & Chabrier 2011; Federrath & Klessen 2012), and these differences are summarized in Federrath & Klessen (2012) and Burkhart (2018). MFF models keep this factor in the integral, which predicts different rates of collapse for different densities, while SFF models take this factor out of the integral. Due to this difference, MFF models predict higher _ff.

₀: The prefactor in Equation (14), ₀, is the local dimensionless efficiency of star formation. ₀ depends on additional processes, such as the level of stellar feedback dispelling some of the gas that is already above n_SF (Burkhart & Mocz 2019). The mass that is converted into stars is then M_* = ₀ M(n > n_SF).

3.1. Different Formalisms of Gravoturbulent Models of Star Formation

3.2. Predictions of Analytical Models

3.3. Connecting Observations to Theory

To compare models with data from our sample, we create a sets of models for each of the three unique prescriptions. We produce LN models over the ${\sigma }_{{\rm{v}}}^{2}/R-{{\rm{\Sigma }}}_{\mathrm{mol}}$ parameter space that encompasses our data in Figure 1. For comparison, we also show the Sun et al. (2018) measurements of cloud-scale observations in nearby galaxies in Figure 1. The mean trend of this parameter space is the empirical trend between velocity dispersion and molecular gas surface density found by Sun et al. (2018) for a subset of PHANGS galaxies at cloud scales,

$$\begin{eqnarray}&&{\mathrm{log}}_{10}\ \left(\displaystyle \frac{{\sigma }_{{\rm{v}}}}{\mathrm{km}\,{{\rm{s}}}^{-1}}\right)={\mathrm{log}}_{10}\ \left(\displaystyle \frac{{{\rm{\Sigma }}}_{\mathrm{mol}}}{{10}^{2}\ {M}_{\odot }{\mathrm{pc}}^{-2}}\right)+0.85.\end{eqnarray} \tag{ 24 }$$

This is shown in Figure 1 as the dotted line. Sun et al. (2018) assume a radius of 40 pc for their sources. Because we assume a fixed LOS depth of 100 pc, we scale the right side of Equation (24) by a factor of $\sqrt{100\ \mathrm{pc}/40\ \mathrm{pc}}$.

Zoom In Zoom Out Reset image size

The vertical spread in Figure 1 is from variations in σ_v that we consider for the LN models, and it roughly reproduces the spread in ${\sigma }_{{\rm{v}}}^{2}/R-{{\rm{\Sigma }}}_{\mathrm{mol}}$ seen in the Sun et al. (2018) sample as well as in our sample of galaxies. For the LN+PL B18 models, we impose Equation (24) and consider a range of α_PL = 1.0–2.3, instead of the range in σ_v considered for the LN models.

We note that the observational efficiency is estimated by taking the ratio of the freefall time to the depletion time, _ff = t_ff/t_dep. As Burkhart & Mocz (2019) point out, a more physically meaningful efficiency for observations of individual star-forming clouds may be the local instantaneous efficiency of star formation at time t, _inst ≈ ₀(t)f_grav(t), where _inst < _ff. This in particular applies to clouds that have evolved sufficiently to have a power-law tail, where f_grav is the fraction of mass that is gravitationally bound (and therefore in the power-law tail). _ff itself is an average of the star formation efficiency of a cloud over the freefall time, which is not directly observable. For simplicity, we assume that our observational estimates of _ff are indeed sensitive to the theoretical definition of _ff. We also compare observational trends of _ff with model trends of f_grav ∝ _inst.

We adopt a constant local efficiency of ₀ = 1% for all models. This value of ₀ is consistent with estimates of _ff from observations, and simulations imply low efficiencies, which range from 0.01 − ≲ 20% (Evans et al. 2009; Lada et al. 2010; Ostriker & Shetty 2011; Krumholz 2014; Zamora-Avilés & Vázquez-Semadeni 2014; Lee et al. 2016; Semenov et al. 2017; Grudić et al. 2018). We summarize the observational measurements and the quantities that we estimate from them in the top half of Table 2. In the bottom half of Table 2, we list the equations used to estimate the same quantities using model outputs.

Table 2. Observational Quantities (Top Panel) Compared to Model Inputs and Outputs (Bottom Panel)

	Observational Estimate	Equation	Unit
	Σmol,obs	αCO(ICO/K km s−1)	M⊙ pc−2
	Σmol,dense,obs	$3.2\,{\alpha }_{\mathrm{CO}}({I}_{\mathrm{HCN}}/$K km s−1)	M⊙ pc−2
(1)	ΣSFR,obs	1.14 × 10−29(L93 GHz/erg s−1 Hz)	M⊙ yr−1 kpc−2
(2)	Pturb,obs	$(3/2){{\rm{\Sigma }}}_{\mathrm{mol}}{\sigma }_{{\rm{v}}}^{2}/R/{k}_{{\rm{B}}}$	K cm−3
(3)	tdep,obs	Σmol/ΣSFR	yr
	n0,obs	Σmol/R	kg cm−3
	tff,obs	$\sqrt{3\pi /32\,G\,{n}_{0,\mathrm{obs}}}$	yr
(4)	ff,obs	tff,obs/tdep,obs	⋯
(5)	fdense	$3.2\,{I}_{\mathrm{HCN}}/{I}_{\mathrm{CO}}$	⋯
(6)	tdep,dense,obs	Σmol,dense/ΣSFR	yr
	Model Estimates	Equation	Unit
(1)	ΣSFR,model	ff,model × Σmol,grid/tff,grid	M⊙ yr−1 kpc−2
(2)	Pturb,model	$(3/2){n}_{0}{\sigma }_{{\rm{v}}}^{2}/{k}_{{\rm{B}}}$	K cm−3
(3)	tdep,model	1/ff,model × tff,grid	yr
(4)	ff,model	Equations (13) or (27) in Burkhart (2018)	⋯
(5)	fgrav	Equations (1), (20) in Burkhart & Mocz (2019)	⋯
(6)	tdep,grav	fgrav/ff,model × tff,grid	yr
(6)	tdep,dense,model	f(n > 104.5 cm−3)/ff,model × tff,grid	yr

Note. We number matching observational and model analogs in the leftmost column. We fix T_kin = 45 K, b = 0.4, and β → ∞ in this analysis, and we model over the observed ranges of Σ_mol and σ_v. Quantities with the subscript "grid" are input from the grid shown in Figure 1.

Download table as:  ASCIITypeset image

In this section, we review general trends between observation-based quantities and the predictions of those from analytical models of star formation. We present a series of plots showing the observational quantities listed in the top half of Table 2 that we use to compare against the model predictions described in the bottom half of Table 2. We focus on variations in star formation timescales (depletion time and freefall time), star formation efficiency, and the dense gas fraction. We do not directly assess variations in emissivity in this paper and defer an analysis of emissivity (and therefore conversion factors) to Bemis (2020) and Bemis & Wilson (in prep.), where we present a modeling of the HCN and CO emissivities for our sample of galaxies on a pixel-by-pixel basis.

The results from analytical models are compared to the data in Figures 2–11. We consider fiducial models that correspond to the parameters that best fit the observed Kennicutt–Schmidt relation (KS; see Figure 3). For each model, we calculate the χ² of the data relative to this fiducial model, and this value is shown in each panel. We note that the best-fit value for α_PL is 1.9 for the LN+PL models. The best-fit relation between Σ_mol − σ_v for the LN-only models is in agreement with the Sun et al. (2018) relation (Equation (24)) within the uncertainty.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We show separate plots with the data and their measurement uncertainties alone, followed by plots comparing the data with the models. We refer to Spearman rank coefficients as a measure of the strength and direction of correlations between two parameters. For reference, a coefficient of r_s = 0 indicates that there is no monotonic relation between the two parameters, negative values indicate negative monotonic relations, and positive values indicate positive monotonic relations. In the following discussion, we use the following definitions: 0.1 ≤ ∣r_s ∣ < 0.4 is considered a weak correlation, 0.4 ≤ ∣r_s ∣ < 0.7 is considered a moderate correlation, and ∣r_s ∣ ≥ 0.7 is considered a strong correlation.

4.1. The Kennicutt–Schmidt Relation

4.2. Does HCN/CO Trace the Star-forming Gas Fraction?

We compare the observed dense gas fraction (${f}_{\mathrm{dense}}={\alpha }_{\mathrm{HCN}}/{\alpha }_{\mathrm{CO}}\,{I}_{\mathrm{HCN}}/{I}_{\mathrm{CO}}$) with estimates of turbulent pressure, P_turb, in Figure 4. Our data show a weak positive trend between f_dense and P_turb on average, and this relation also holds within most of the individual galaxies in our sample. We plot model predictions of f_grav versus P_turb in Figure 5 and show the outline of the data relation between f_dense and P_turb in the background. For comparison, in the bottom rows of Figure 5, we plot the fraction of gas above a fixed density n = 10^4.5 cm⁻³ for all of the three model n − PDFs prescriptions.

• 1.  
  LN PN11: The LN varying-threshold models also predict a negative trend between f_grav and P_turb. On average, f_grav is lower than f_dense predicted by the data. Similar to the LN+PL B18 models, we see a positive trend between f(n > 10^4.5 cm⁻³) and P_turb. However, there is very little spread in the model f(n > 10^4.5 cm⁻³) versus P_turb relation. This indicates that the spread observed in the f_dense versus P_turb relation in our data is poorly reproduced by variations in σ_v alone.
• 2.  
  LN+PL B18: For constant α_PL, these models predict a negative trend between f_grav and P_turb, and this becomes steeper for larger α_PL. The data primarily overlap with f_grav for models with shallower values of α_PL. In contrast, we see a positive trend between f(n > 10^4.5 cm⁻³) predicted by the models and P_turb. Lower values of α_PL produce a flatter relation between f(n > 10^4.5 cm⁻³) at low P_turb. This results in a broader spread in this relation at lower P_turb. In general, the spread observed in the f_dense versus P_turb relation in our data is poorly reproduced by variations in α_PL alone.
• 3.  
  LN Fixed: f_grav predicted by the fixed density-threshold models is identical to f(n > 10^4.5 cm⁻³) because we set the threshold for these models to n = 10^4.5 cm⁻³). As with the other models, we see a positive relation between f(n > 10^4.5 cm⁻³) and P_turb. The spread in the data is poorly reproduced by these models.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The predictions of f_grav from the varying-threshold models show negative trends with P_turb for constant α_PL (LN+PL B18 models), in contrast to the positive trends observed with f(n > 10^4.5 cm⁻³) for all models. f_grav is equivalent to f(n > 10^4.5 cm⁻³) in the fixed-threshold models. From these results, the weak positive trend of ${f}_{\mathrm{dense}}\propto {I}_{\mathrm{HCN}}/{I}_{\mathrm{CO}}$ with P_turb seen in our data is more consistent with the positive trend between f(n > 10^4.5 cm⁻³) and P_turb seen in the models, and it supports the conclusion that ${I}_{\mathrm{HCN}}/{I}_{\mathrm{CO}}$ traces gas above a roughly constant density, but not necessarily the f_grav predicted by gravoturbulent models of star formation.

4.2.1. The Star Formation Efficiency per Freefall Time

As a check on the results above, we also consider _ff. In Figure 6 we plot our data as a function of P_turb (left panel) and f_dense (right panel). In Figure 7 we compare the model predictions of _ff with P_turb (top row), f(n > 10^4.5 cm⁻³) (middle row), and f_grav (bottom row). _ff shows a weak negative correlation with P_turb and a weak positive correlation with f_dense (Figure 6) when considering our sample as a whole. These results agree with the qualitative predictions of the varying density-threshold models (discussed in Section 3), where higher ${ \mathcal M }\sim {\sigma }_{{\rm{v}}}$ (and therefore ${P}_{\mathrm{turb}}\sim {\sigma }_{v}^{2}$) yields a lower _ff.

• 1.  
  LN PN11:The negative trend between _ff and P_turb seen in the data is not reproduced by the LN PN11 models (middle column, top row of Figure 7). However, the LN PN11 models do reproduce the positive trend between _ff and f_dense seen in the data when we consider f(n > 10^4.5 cm⁻³) (middle row of Figure 7). There is some overlap between _ff and model f_grav (bottom row of Figure 7), but this trend does not track the positive relation between _ff and f_dense seen in the data.
• 2.  
  LN+PL B18:Similar to the LN PN11 models, the negative trend between _ff and P_turb seen in the data is not reproduced by any single value of α_PL for these models (left column, top row of Figure 7). Rather, the data overlap with these models for a range of α_PL, with higher values of _ff corresponding to lower values of α_PL and lower P_turb.Moreover, similar to the LN PN11 models, the LN+PL B18 models do reproduce the positive trend between _ff and f_dense seen in the data when we consider f(n > 10^4.5 cm⁻³) (middle row of Figure 7). In this case, the vertical spread in the data corresponds to variations in α_PL. Models with lower α_PL also show an upturn in _ff at higher f(n > 10^4.5 cm⁻³) consistent with the trend seen in the data.At first look, the positive trend between _ff and f_dense seen in the data also appears consistent with the trend between _ff and model f_grav (bottom row of Figure 7). However, this requires that f_dense ∝ f_grav, for which we find evidence to the contrary in the previous section.
• 3.  
  LN Fixed: The fixed-threshold models predict an increase in _ff with P_turb up to P_turb ≈ 10⁷ K cm⁻³, after which this relation flattens and _ff slightly turns over. This is not seen in our data. These models predict a positive relation between _ff and f_grav up to $\mathrm{log}\ {f}_{\mathrm{grav}}\approx -0.5$, above which _ff turns over. There is no evidence for this turnover in our data.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The varying-threshold models are able to reproduce the positive trend seen between _ff and f_dense when considering f(n > 10^4.5 cm⁻³). Some of the scatter is also reproduced by these models when considering a range in α_PL or a range in σ_v (see Figure 7).

4.2.2. Total Gas Depletion Times

Now we consider whether estimates of the total gas depletion times may give insight into the discrepancy between f_grav and observational estimates of f_dense. We plot t_dep as a function of P_turb and f_dense for our sample of galaxies in Figure 8. We find moderate negative correlations between t_dep and the dense gas fraction traced by ${I}_{\mathrm{HCN}}/{I}_{\mathrm{CO}}$. t_dep estimated from the data appears relatively constant with P_turb, but the Spearman rank coefficient of the combined data indicates a moderate negative correlation (see Figure 8).

• 1.  
  LN PN11: These models are able to reproduce the average trends seen in our data (i.e., negative trend between t_dep and P_turb and negative trend between t_dep and f(n > 10^4.5 cm⁻³)). Variations in P_turb alone do not reproduce the spread seen in our data. t_dep versus f_grav does not track the relation between t_dep and f_dense seen in the data.
• 2.  
  LN+PL B18: We find good agreement between these models and the data when comparing t_dep, P_turb, and f(n > 10^4.5 cm⁻³) (Figure 9). Although the data overlap with t_dep versus f_grav, the models show a much larger spread than that observed in the data.
• 3.  
  LN Fixed: These models predict that t_dep decreases until P_turb ∼ 10⁷ K cm⁻³ until it reaches a roughly constant value. This behavior is not seen in our data. t_dep monotonically decreases with f(n > 10^4.5 cm⁻³) and f_grav for these models.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The varying-threshold models best reproduce the observed trend between t_dep and P_turb. The varying-threshold models also perform better than the fixed-threshold models at reproducing the observed trend between t_dep and f_dense. Yet again, the trends in the data are best reproduced by the models when considering f(n > 10^4.5 cm⁻³) rather than f_grav.

4.2.3. Dense Gas Depletion Times

We plot the depletion time of the dense gas as traced by HCN, t_dep,dense, as a function of P_turb (left) and dense gas fraction (right) in Figure 10. We see a moderate positive correlation between t_dep,dense and f_dense in ∼3 of our sources and the EMPIRE galaxies. The Spearman rank coefficient shows a weak positive correlation for the combined data set. At first, this appears to contradict the expectation of star formation models, but all models except for those with a fixed threshold are able to roughly reproduce this increase in t_dep,dense with gas fraction (see Figure 11).

• 1.  
  LN PN11: These models predict a turnover in t_dep(n >10^4.5 cm⁻³) with P_turb and in t_dep(n > 10^4.5 cm⁻³) with f(n > 10^4.5 cm⁻³) (see the middle column, top and middle rows in Figure 11). This qualitatively agrees with the average trends seen in our data. We find that t_dep,grav versus f_grav does not track the relation between t_dep and f_dense seen in the data. Variations in P_turb alone do not reproduce the spread seen in our data.
• 2.  
  LN+PL B18: For α_PL > 1.7, these models predict a turnover in t_dep(n > 10^4.5 cm⁻³) around P_turb ∼ 10⁷ K cm⁻³ and for $\mathrm{log}\ f(n\gt {10}^{4.5}\ {\mathrm{cm}}^{-3})\sim -0.5$. For α_PL ≤ 1.7, t_dep(n > 10^4.5 cm⁻³) consistently decreases with P_turb and f(n > 10^4.5 cm⁻³). The observed trends in the data are in qualitative agreement with the trends seen for α_PL > 1.7. We include t_dep, grav versus f_grav for completeness, but again find that the models show a much larger spread than is observed in the data.
• 3.  
  LN Fixed: As predicted in Section 3.2, these models return a fixed t_dep(n > 10^4.5 cm⁻³) regardless of the variations in P_turb or f(n > 10^4.5 cm⁻³).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The LN+PL B18 models are able to qualitatively reproduce the trends observed in t_dep,dense, P_turb, and f_dense for α_PL > 1.7. The LN PN11 models are also able to reproduce the average trends, but variations in P_turb do not reproduce the observed scatter. These results again refute the interpretation of the ${I}_{\mathrm{HCN}}/{I}_{\mathrm{CO}}$ ratio as a a tracer of f_grav. The observed trends are best reproduced by the varying-threshold models when considering f(n > 10^4.5) cm⁻³, but not with f_grav, and our results therefore refute this interpretation.

In this work, we explore how well the HCN-to-CO ratio traces the star-forming fraction of gas in molecular clouds across the more extreme environments of U/LIRGs, mergers, and galaxy centers using gravoturbulent models of star formation (Krumholz & McKee 2005; Hennebelle & Chabrier 2011; Padoan & Nordlund 2011; Federrath & Klessen 2012; Burkhart 2018; Burkhart & Mocz 2019). Previous studies find that HCN may be tracing gas above a fairly constant density and that the mass traced by HCN appears to scale linearly with the SFR, indicating a constant average depletion time of the dense gas that spans clouds (see Wu et al. 2005) to entire galaxies, including disk galaxies and U/LIRGs (see Gao & Solomon 2004a, 2004b). This result is consistent with the predictions of the LN-only models with a fixed threshold that we consider in our analysis (see the right columns of Figures 3, 5, 7, 9, and 11). These models are able to qualitatively reproduce some of the average trends involving the dense molecular gas traced by HCN, such as the decrease in t_dep with increasing ${f}_{\mathrm{dense}}\propto {I}_{\mathrm{HCN}}/{I}_{\mathrm{CO}}$ (Figure 9). However, these models fail to capture the spread observed in many of the trends of our data and also the observed trends with P_turb.

The behavior of star formation laws observed in disk galaxies may give some insight into the discrepancy between the fixed-threshold models and observations. There is evidence of nonlinear scalings between dense gas mass and the SFR within disk galaxies that appears to be correlated with galactic radius in disk galaxies (e.g., Chen et al. 2015; Usero et al. 2015; Gallagher et al. 2018). Relative to their disks, the centers of late-type galaxies show longer dense gas depletion times despite larger dense gas fractions traced by ${I}_{\mathrm{HCN}}/{I}_{\mathrm{CO}}$ (e.g., Chen et al. 2015; Usero et al. 2015; Gallagher et al. 2018). Similarly, the nuclei of the merging Antennae galaxies also display longer dense gas depletion times despite higher dense gas fractions compared to the overlap region of this merging system (Bigiel et al. 2015; Bemis & Wilson 2019). Ambient pressure of the ISM is higher in the centers of disk galaxies compared to larger radii, which is apparent from observational estimates of P_turb from CO across the disks of nearby disk galaxies (see Gallagher et al. 2018). If higher P_turb is indicative of stronger turbulent support, then gravoturbulent models of star formation with varying star formation thresholds may offer an explanation for these observed trends (i.e., Krumholz & McKee 2005; Padoan & Nordlund 2011; Burkhart 2018).

The varying-threshold models considered in this paper (i.e., Padoan & Nordlund 2011; Burkhart 2018) predict that the gas density required for the onset of collapse (and therefore star formation) increases in the presence of supportive processes such as solenoidal turbulence. This is apparent in Figure 5 (left two columns, top row) where the LN+PL B18 and LN P11 models show a decrease in f_grav with increasing P_turb. Turbulence can also increase the mean gas density of the molecular cloud as a whole by widening the n − PDF. This is also apparent in Figure 5 (left two columns, bottom row) where the LN+PL B18 and LN P11 models show an increase in f(n > 10^4.5 cm⁻³) with increasing P_turb. The difference in behavior between f_grav and f(n > 10^4.5 cm ⁻³) with P_turb reflects an important prediction of these varying-threshold models: the supportive effect of turbulence can have a more significant impact on star formation than the increase in dense gas mass due to the widening of the n − PDF.

A wider n − PDF would naturally result in enhanced ${I}_{\mathrm{HCN}}/{I}_{\mathrm{CO}}$ if this ratio were tracing gas above a roughly constant density (see Leroy et al. 2017a; Bemis 2020). At the same time, the SFR itself can appear suppressed relative to the total dense gas mass because the onset of star formation occurs at a higher density than the mean density of the gas traced by ${I}_{\mathrm{HCN}}/{I}_{\mathrm{CO}}$. The results of our work support this picture. In particular, the LN+PL B18 and LN P11 models are able to qualitatively reproduce the average trends observed in our data when considering f(n > 10^4.5 cm⁻³) as a proxy for the dense gas fraction traced by ${I}_{\mathrm{HCN}}/{I}_{\mathrm{CO}}$ (see Figures 5, 7, 9, and 11).

One apparent discrepancy between the data and the varying-threshold models is the behavior of _ff with P_turb (see the top row in Figure 7). Our sample shows a weak negative trend between _ff and P_turb and the EMPIRE sample has a moderate negative trend. In contrast, the varying-threshold models predict an increase in _ff with P_turb (see Figure 7). The fixed-threshold LN model predicts that _ff initially increases with P_turb and then decreases toward higher pressures. The discrepancy between the models and data is likely due in part to inaccurate estimates of t_ff from our data. When estimating the mean density of the molecular gas traced by CO, we assume a fixed LOS depth of 100 pc for all measurements, including those from the galaxies in the EMPIRE sample. We can also see from Figure 9 that the models predict a steeper decline in t_dep with P_turb than what is observed in the data. These two things combined likely contribute to the discrepancy between the data and the models when considering _ff as a function of P_turb. Due to the weakness of the observed trend in our data and the uncertainty in t_ff, it is difficult to determine whether _ff truly increases or decreases with P_turb. The apparent agreement between _ff and f(n > 10^4.5 cm⁻³) predicted by the varying-threshold models and our data may then be a reflection of the stronger connection between star formation and denser gas that is better traced by HCN relative to the bulk molecular gas traced by CO (see Gao & Solomon 2004a, 2004b).

We also find that the scatter of the data is only poorly reproduced by variations in P_turb when considering the LN-only models (e.g., Figures 5, 7, 9, and 11). One parameter that we did not vary in this analysis is α_vir. We consider here how this would impact our results. For the Padoan & Nordlund (2011) models, higher α_vir increases n_SF such that more unbound clouds will have lower _ff. Variations in α_vir would introduce some scatter into the results, primarily in the relations that depend on _ff and t_dep. However, variations in ${ \mathcal M }$ have an impact on both the width of the n − PDF (Equation (17)) and the shift in n_SF (Equation (19)) for the LN PN11 models, while changes in α_vir primarily have an impact on n_SF. Changes in the dense gas fraction indicated by f(n > 10^4.5 cm⁻³) are not analytically dependent on α_vir and therefore would need an additional explanation.

The LN+PL B18 models offer an alternative explanation for the scatter of the data: variations in power-law slope (α_PL). Simulations show that a power-law tail develops as a consequence of cloud evolution as gas becomes self-gravitating, and the slope of the power-law tail becomes more shallow as the fraction of bound gas (i.e., f_grav) increases (see Ballesteros-Paredes et al. 2011; Collins et al. 2012; Schneider et al. 2015; Burkhart et al. 2017; Padoan et al. 2017). This is supported by observations of some Milky Way clouds, which show a power-law tail at high densities (e.g., Kainulainen et al. 2009; Schneider et al. 2013; Lombardi et al. 2015; Schneider et al. 2015, 2016; Alves et al. 2017). Additionally, there is observational evidence for a connection between power-law slope and _ff that matches the predictions of the LN+PL B18 models (Federrath & Klessen 2013; Burkhart & Mocz 2019).

The scatter in our data is well reproduced for α_PL = 1–2.3 (see the left column of Figures 3, 5, 7, 9, and 11). These models are also able to reproduce the average trends of our data. Out of the three models we consider, the LN+PL B18 models perform best at reproducing the observational trends and scatter of our data and the EMPIRE sample. In the context of the LN+PL B18 models, the scatter in our data is therefore driven in part due to cloud evolution (Ballesteros-Paredes et al. 2011), such that the n − PDF of clouds evolves from a lognormal shape to a composite lognormal and power law over time as more gas becomes gravitationally bound. The underlying trends are still driven largely by variations in turbulence, which we impose in our model parameter space by using the observational relation between gas surface density and velocity dispersion found by Sun et al. (2018; see Equation (24)).

We find the following conclusions from our comparison between observation and the predictions of analytical models of star formation:

• 1.  
  KS relation: For the model parameter space considered, none of the models reproduce the multislope KS relation as found in Wilson et al. (2019). Variations in the slope of the Σ_mol − σ_v relation can produce variations in the KS slope, such that a shallower (steeper) Σ_mol − σ_v relation would produce a shallower (steeper) KS relation. Furthermore, systematic changes in ₀ or power-law slope with Σ_mol could also affect the KS slope. An increasing (decreasing) power-law slope with Σ_mol would produce a shallower (steeper) KS relation. Similarly, a decreasing (increasing) ₀ with Σ_mol produces a shallower (steeper) KS slope.
• 2.  
  f_dense and f_grav: For the PN11 and B18 models, f_grav decreases with P_turb and increases for the fixed-threshold models. f(n > 10^4.5 cm⁻³) increases with P_turb for all models. The behavior of f(n > 10^4.5 cm⁻³) with P_turb is in qualitative agreement with our data, which show an increase in ${f}_{\mathrm{dense}}=({\alpha }_{\mathrm{HCN}}/{\alpha }_{\mathrm{CO}}){I}_{\mathrm{HCN}}/{I}_{\mathrm{CO}}$ with observational estimates of P_turb. We therefore conclude that the ${I}_{\mathrm{HCN}}/{I}_{\mathrm{CO}}$ ratio likely does not track f_grav as predicted by the varying-threshold models, but rather the fraction of gas above some fixed density, such as f(n > 10^4.5 cm⁻³).
• 3.  
  _ff: For the PN11 and B18 models, _ff increases with P_turb on average. For the fixed-threshold models, _ff increases with P_turb until P_turb/k_B ≈ 10^7.5 K cm⁻³ where it turns over. All models predict an increase in _ff with f(n > 10^4.5 cm⁻³), on average, with only the fixed-threshold models showing a turnover in this relation at high fractions. _ff decreases with f_grav for the the PN11 models and the B18 models for fixed α_PL. For the fixed-threshold models, _ff increases with f_grav until it turns over at high fractions. Under the assumption of a fixed LOS depth, our data show a weak decrease in _ff with P_turb and increase with f_dense. This is in qualitative agreement with the PN11 and B18 models, assuming ${I}_{\mathrm{HCN}}/{I}_{\mathrm{CO}}$ is tracking f(n > 10^4.5 cm⁻³), but not f_grav.
• 4.  
  Total gas depletion time: All models predict a decrease in t_dep with both P_turb and f(n > 10^4.5 cm⁻³). The fixed-threshold models predict a flattening of t_dep with higher P_turb, but a constant decrease with f(n > 10^4.5 cm⁻³). In contrast, the PN11 models and the B18 models show a constant decrease in t_dep with P_turb, but a steepening decline in t_dep with higher f(n > 10^4.5 cm⁻³). t_dep increases with f_grav for the PN11 models and the B18 models for fixed α_PL. The trends in data are in better qualitative agreement with the PN11 and B18 models, again assuming ${I}_{\mathrm{HCN}}/{I}_{\mathrm{CO}}$ is tracking f(n > 10^4.5 cm⁻³), but not f_grav.
• 5.  
  Dense gas depletion time: The PN11 and B18 models (for α_PL > 1.7) predict a turnover in the depletion time of dense gas, t_dep(n > 10^4.5 cm⁻³), with P_turb and f(n > 10^4.5 cm⁻³), while the fixed-threshold models predict a constant dense gas depletion time (Gao & Solomon 2004a, 2004b). We note that α_PL < 1.7 produces a nearly constant decrease in t_dep(n >10^4.5 cm⁻³) with P_turb and f(n > 10^4.5 cm⁻³). t_dep,grav increases with f_grav for the PN11 models and the B18 models for fixed α_PL. The trends in data are in better qualitative agreement with the PN11 and B18 models, again assuming ${I}_{\mathrm{HCN}}/{I}_{\mathrm{CO}}$ is tracking f(n > 10^4.5 cm⁻³), but not f_grav.
• 6.  
  Variations in PL slope (α_PL) (lognormal+power-law B18 models) are able to explain the scatter in the data. The scatter in our data is only poorly reproduced by variations in P_turb (or ${ \mathcal M }$) alone.
• 7.  
  The varying-threshold models (lognormal+power-law B18 and lognormal-only PN11) are able to better reproduce the average trends of our data and the EMPIRE sample compared to the fixed-threshold lognormal-only models.
• 8.  
  The average trends in the data with P_turb and f_dense are likely set by variations in turbulence determined by local environment within each galaxy.

In summary, the HCN/CO ratio likely is a reliable tracer of gas above a constant density, such as n_SF ≈ 10^4.5 cm⁻³, but not necessarily f_grav. If varying thresholds for star-forming gas exist in nature, then the HCN/CO ratio is a poor tracer of f_grav.

The authors thank the referee for their constructive report and suggestions that have improved the paper. A.B. acknowledges support from an Ontario Trillium Scholarship. The research of CDW is supported by grants from the Natural Sciences and Engineering Research Council of Canada and the Canada Research Chairs program. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2011.0.00467.S, ADS/JAO.ALMA#2011.0.00525.S, ADS/JAO.ALMA#2011.0.00772.S, ADS/JAO.ALMA#2012.1.00165.S, ADS/JAO.ALMA#2012 .1.00185.S, ADS/JAO.ALMA#2012.1.01004.S, ADS/JAO.ALMA#2013.1.00218.S, ADS/JAO.ALMA#2013.1.00247.S, ADS/JAO.ALMA#2013.1.00634.S, ADS/JAO.ALMA#2013.1.00885.S, ADS/JAO.ALMA#2013.1.00911.S, ADS/JAO.ALMA#2013.1.01057.S, ADS/JAO.ALMA#2015.1.00993.S, ADS/JAO.ALMA#2015.1.01177.S, ADS/JAO.ALMA#2015.1.01286.S, ADS/JAO.ALMA#2015.1.01538.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ.

Software: This work made use of the following software packages: astropy (Astropy Collaboration et al. 2013, 2018,2022), casa (McMullin et al. 2007), cmasher (van der Velden 2020), matplotlib (Hunter 2007), numpy (Harris et al. 2020), scipy (Virtanen et al. 2020), and spectral-cube (https://spectral-cube.readthedocs.io/en/latest/).

We use the following expressions to calculate moments and their corresponding uncertainties:

$$\begin{eqnarray}&&{\rm\small{MOMENT}}0:\quad I={{\rm{\Sigma }}}_{i}{T}_{i}{\rm{\Delta }}v\end{eqnarray} \tag{ A1 }$$

$$\begin{eqnarray}&&\delta I=\delta \bar{T}{\rm{\Delta }}v\sqrt{{N}_{\mathrm{chan}}}\end{eqnarray} \tag{ A2 }$$

$$\begin{eqnarray}&&{\rm\small{MOMENT}}1:\quad \bar{v}=\displaystyle \frac{{{\rm{\Sigma }}}_{i}{T}_{i}{v}_{i}}{{{\rm{\Sigma }}}_{i}{T}_{i}}\end{eqnarray} \tag{ A3 }$$

$$\begin{eqnarray}&&\delta \bar{v}=\displaystyle \frac{{N}_{\mathrm{chan}}{\rm{\Delta }}v}{2\sqrt{3}}\displaystyle \frac{\delta I}{I}\end{eqnarray} \tag{ A4 }$$

$$\begin{eqnarray}&&{\rm\small{MOMENT}}2:\qquad {\sigma }_{v}=\sqrt{\displaystyle \frac{{{\rm{\Sigma }}}_{i}{T}_{i}{\left({v}_{i}-\bar{v}\right)}^{2}}{{{\rm{\Sigma }}}_{i}{T}_{i}}}\end{eqnarray} \tag{ A5 }$$

$$\begin{eqnarray}&&\delta {\sigma }_{v}=\displaystyle \frac{{\delta }_{I}}{I}\displaystyle \frac{{\left({N}_{\mathrm{chan}}{\rm{\Delta }}v\right)}^{2}}{8\sqrt{5}}\displaystyle \frac{1}{{\sigma }_{v}}.\end{eqnarray} \tag{ A6 }$$

A full derivation of the uncertainties on Equations (A4) and (A6) is given in Wilson et al. (2023).