© ESO, 2017

1. Introduction

Star formation within galactic disks was proceeding much faster in the first half of the history of the Universe: the cosmic star-formation rate density declined by a factor of approximately ten since z = 1 (Madau et al. 1998; Hopkins & Beacom 2006). The mean star-formation rate with respect to the total stellar mass also decreases with decreasing redshift: star-forming galaxies form a “main sequence” in the star-formation/stellar-mass space (e.g., Brinchmann et al. 2004; Daddi et al. 2007; Elbaz et al. 2007; Noeske et al. 2007a,b; Salim et al. 2007; Whitaker et al. 2012). In addition, galactic starbursts (e.g., ultraluminous infrared galaxies (ULIRGs) or smm-galaxies) represent outliers from the main sequence. The slope and offset of the Ṁ_∗–M_∗ powerlaw main sequence relation change with redshift (Speagle et al. 2014). The most prominent change is an increasing offset with increasing redshift, that is, the specific star-formation rate (Ṁ_∗/M_∗) increases significantly with increasing z (by a factor of ~6 for galaxies with masses of M_∗ ~ 3 × 10¹⁰M_⊙; Pannella et al. 2015). The slope and scatter of this correlation, the evolution of its normalization with cosmic time, contain crucial and still poorly known information on galaxy evolution (e.g., Karim et al. 2011; Rodighiero et al. 2011; Wuyts et al. 2011; Sargent et al. 2012).

Two factors can be invoked to explain the higher star-formation efficiency: (i) higher gas fractions (see, e.g., Combes et al. 2013); and (ii) dynamical trigger of interactions, whose frequency increases with redshift (e.g., Conselice et al. 2009; Kartaltepe et al. 2012).

In this work, we take a closer look at the gas content or fraction and the associated star-formation rate in main sequence and starburst galaxies at z = 0 and z ~ 1–2. We look preferentially at local starburst galaxies, ULIRGs, and high-redshift starbursts smm-galaxies.

Genzel et al. (2010) and Tacconi et al. (2012) showed that star-forming galaxies at z = 1–2 have higher gas fractions (~33%) and higher star-formation efficiencies with respect to the molecular gas (SFE = SFR/M_H2 ~ 1/0.7 Gyr^-1) than local spiral galaxies (~10% and SFE ~ 1/2.0 Gyr^-1; e.g., Bigiel et al. 2008; Leroy et al. 2008). Local ULIRGs and high-redshift smm-galaxies have the highest star-formation efficiencies (e.g., Pope et al. 2013).

The uncertainties of the determined star-formation rates are typically ~50% (e.g., Leroy et al. 2008). Since molecular hydrogen at temperatures below 100 K is not directly detectable, one has to rely on CO, HCN, or HCO⁺ observations to derive H₂ gas masses. Unfortunately, the associated conversion factors have high uncertainties (approx. a factor of two, e.g., Bolatto et al. 2013). Our understanding of the gas content and thus the star-formation efficiency is limited by these uncertainties.

A complementary way to determine galactic gas masses is the direct modeling of molecular emission. Narayanan & Krumholz (2014) combined numerical simulations of disc galaxies and galaxy mergers with molecular line radiative transfer calculations to develop a model for the physical parameters that drive variations in CO spectral line energy distributions (SLEDs) in galaxies in terms of the star-formation-rate density. Their model was able to reproduce the SLEDs of galaxies over a dynamic range of approximately 200 in star-formation-rate surface density. However, the CO high-J transitions (J> 8) of ULIRGs are difficult to reproduce within the model (Fig. 2 of Kamenetzky et al. 2016).

Bournaud et al. (2015) modeled the intensity of CO emission lines, based on hydrodynamic simulations of spirals, mergers, and high-redshift galaxies with very high resolutions (3 pc and 10³M_⊙) and detailed models for the phase-space structure of the interstellar gas, including shock heating, stellar feedback processes, and galactic winds. The simulations were analyzed with a large velocity gradient (LVG) model to compute the local emission in various molecular lines in each resolution element, radiation transfer, opacity effect, and the intensity emerging from galaxies to generate synthetic spectra for various CO transitions. This model reproduced the known properties of CO spectra and CO-to-H₂ conversion factors in nearby spirals and starbursting major mergers.

Alternatively, galactic gas disks can be modeled analytically, assuming axis-symmetry. Krumholz & Thompson (2007) provided a simple model for understanding how Kennicutt-Schmidt laws, which relate the star-formation rate to the mass or surface density of gas as inferred from some particular line, depend on the line chosen to define the correlation. They assume a probability distribution for the mass fraction of gas at a given density and calculate the molecular emission with an escape probability formalism. The model gas clouds have constant temperature, Mach number, and optical depth. Their results showed that for a turbulent medium, the luminosity per unit volume in a given line, provided that this line can be excited at temperatures lower than the mean temperature in a galaxy’s molecular clouds, increased faster than linearly with the density for molecules with critical densities larger than the median gas density. The star-formation rate also rose superlinearly with the gas density, and the combination of these two effects produced a close to linear correlation between star-formation rate and line luminosity.

Kazandjian et al. (2015) investigated the effect of mechanical heating on atomic fine-structure and molecular lines and on their ratios. They tried to use those ratios as a diagnostic to constrain the amount of mechanical heating in an object and also study its significance on estimating the H₂ mass. Equilibrium photodissociation models (PDRs) were used to compute the thermal and chemical balance for the clouds. The equilibria were solved for numerically using the optimized version of the Leiden PDR-XDR code. Large velocity-gradient calculations were done as post-processing on the output of the PDR models using RADEX (van der Tak et al. 2007). They showed that high-J CO line ratios and ratios involving HCN are very sensitive to mechanical heating.

In this work, we investigate a different analytical approach, where we extend the model of galactic clumpy gas disks presented in Vollmer & Leroy (2011). The model has a large-scale and a small-scale part. The large-scale part gives the surface density, turbulent velocity, disk height, and gas viscosity. The small-scale part begins at densities where gas clouds become self-gravitating. The non-self-gravitating and self-gravitating clouds obey different scaling relations, which are set by observations. For clouds of a certain density, the area filling factor is calculated. The gas and dust temperatures are calculated by the heating and cooling equilibrium. Dense clouds are heated by turbulent, mechanical, and cosmic-ray heating. For all model clouds, the size, density, temperature, and velocity dispersion are known. The molecular abundances of individual gas clouds are determined by a detailed chemical network involving the cloud lifetime, density, and temperature. Molecular line emission is calculated with an escape probability formalism. The model is applied to samples of local spiral galaxies, ULIRGs, high-z star-forming galaxies, and smm-galaxies. The model simultaneously calculates the total gas mass, H₂ mass, the gas velocity dispersion, Hi mass, IR luminosity, IR SED, CO SLED, HCN(1–0), and HCO⁺(1–0) emission of a galaxy given its size, integrated star-formation rate, stellar mass radial profile, rotation curve, and Toomre Q parameter. In addition, the temperature, density, velocity dispersion, and molecular abundance of a gas cloud at a given density can be retrieved.

This article presents a sophisticated model and justifies the physics behind it, shows the results we obtain for different choices of input parameters, and compares with observations. The physical processes included in the model are described in detail in Sect. 2. The steps involved in the calculations and in determining the uncertainties are described respectively in Sects. 3 and 4. Not all spirals are the same so rather than define a “typical” member of each of the four classes (local spirals, local ULIRGs, submillimeter (smm) galaxies, and high-z star-forming galaxies), we use true samples of real objects for which we think we can estimate the appropriate values of the input parameters. The samples and their origin are described in Sect. 5 and the results of the calculations for the samples are shown in detail in Sect. 6. Hence, the reader interested in how well the model reproduces the observations can go directly to Sect. 5 (short) or even Sect. 6. Section 7 evaluates the influence of the choice of the chemical network, the Toomre Q parameter, and the length scale parameter δ in terms of the effect on line and continuum emission. The importance of each of the heating and cooling processes is described in Sect. 8 and that of the assumed gas properties (mass, velocity dispersion, and free-fall time) in Sect. 9. Finally, we give our conclusion in Sect. 10.

2. The analytical model

Compared to the model described in Sect. 2 of Vollmer & Leroy (2011) that is based on Vollmer & Beckert (2003; VB03), the present model does not include a break radius, where the star-formation timescale changes from the free fall timescale to the molecular formation timescale. In addition, we included in this more advanced model (i) the determination of the dense gas fraction; (ii) ISM scaling relations; (iii) the determination of dust and gas temperatures; (iv) a chemical network for the determination of molecular abundances; (v) a formalism for the photodissociation of molecules by the interstellar radiation field; and (vi) the determination of dust and molecular line emission.

The model considers the warm, cool neutral, and molecular phases of the ISM as a single, turbulent gas. We assume this gas to be in vertical hydrostatic equilibrium with the midplane pressure balancing the weight of the gas and stellar disk. The gas is assumed to be clumpy, so that the local density is enhanced relative to the average density of the disk. Using this local density, we calculate two timescales relevant to star formation: the free-fall timescale of an individual clump and the characteristic timescale for H₂ to form on grains. The free-fall timescale is taken as the governing timescale for star formation. The star-formation rate is used to calculate the rate of energy injection by supernovae. This rate is related to the turbulent velocity dispersion and the driving scale of turbulence. These quantities, in turn, provide estimates of the clumpiness of gas in the disk (i.e., the contrast between local and average density) and the rate at which viscosity moves matter inward.

The model relies on several empirical calibrations: e.g., the relationship between star-formation rate and energy injected into the ISM by supernovae, the H₂ formation timescale (and its dependence on metallicity), and the turbulent dimension of the ISM (used to relate the driving length scale to the characteristic cloud size modulo a free parameter δ, which is constrained by observations). As far as possible, these are drawn from observations of the Milky Way.

The model only contains two free parameters. First, there is an unknown scaling factor relating the driving length of turbulence to the size of gravitationally bound clumps, which we call δ. Second, the mass accretion rate, Ṁ, which is related to the driving length and turbulent velocity, is a free parameter. From a detailed comparison of local spiral galaxies from the THINGS survey, Vollmer & Leroy (2011) found δ = 5 ± 3. For simplicity, we assume δ = 5 in this work. Moreover, the Toomre Q parameter of the gas is set to the observed values for the local spirals (2–8, e.g., Leroy et al. 2008), and to Q = 1.5 for the ULIRGs, high-z star-forming and submillimeter galaxies.

In the remainder of this section, we discuss our assumptions in slightly more detail, justify them via comparison to observation and theory, and note the physics that we neglect.

2.1. The interstellar medium (ISM)

Following, for example, Mac Low & Klessen (2004), the warm, cool neutral, and molecular phases of the ISM are viewed as a single entity. Locally, the exact phase of the gas depends on the local pressure, metallicity, stellar radiation field, stellar winds, and shocks. Here, we view these factors as secondary, making a few simplifying assumptions. The equilibrium between the different phases of the ISM and the equilibrium between turbulence and star formation depends on three local timescales: the turbulent crossing time , the molecule formation timescale , and the local free-fall timescale of a cloud. In addition, photodissociation of molecules is taken into account.

2.1.1. The fraction of dense gas

To calculate the mass fraction between two gas densities, we use the density probability distribution function of Padoan et al. (1997) for overdensities x: (1)where the standard deviation, σ, is given by (2)and ℳ = v_turb/c_s is the Mach number with the sound speed c_s. The mass fraction of gas with overdensities exceeding x is then (3)The overdensity for a given density ρ₁ is calculated with respect to the midplane gas density x = ρ₁/ρ.

For the calculation of the molecular line emission we divide the ISM into two density bins: (i) densities ρ₁ equal to or higher than that of the self-gravitating clouds (see Sect. 2.7): ρ₁ ≥ ρ/ Φ_V and (ii) densities ρ₂ equal to or higher than that of the cool neutral medium: ρ_CNM ≤ ρ₂ ≤ ρ/ Φ_V, where Φ_V is the volume-filling factor of the largest self-gravitating structures in the disk, computed via a procedure described in Sect. 2.7. Following Wolfire et al. (2003), we set the minimum density of the cool neutral medium to (4)With respect to Eq. (35) of Wolfire et al. (2003), we set the normalized FUV radiation field and the normalized dust abundances and gas metallicities (Eq. (28)), where α is the constant of the molecule-formation timescale (Eq. (26)). Moreover, we set the normalized total ionization rate by cosmic rays and EUV/X-rays . If the CNM density exceeds the midplane density, we set ρ_CNM = ρ.

The mass fraction of the self-gravitating clouds with respect to the diffuse clouds is a major unknown. Using the lognormal pdf of Padoan et al. (1997) neglects self-gravitation, which can change the shape of the pdf significantly (e.g., Schneider et al. 2015). Based on the findings of the latter authors, we adopt the following recipe:

• density bin (ii): (5)

• density bin (i): (6)where x_CNM and x_sg are the overdensities of the cool neutral medium and self-gravitating clouds.

The normalization factor is (7)Within the density bin, the mass fraction of clouds of overdensity between x₁ and x₂ is calculated as the difference between the mass fractions. For the determination of the Mach number, we adopt the temperature of the cool neutral medium (~100 K) to calculate the sound speed. For the self-gravitating, clouds we adopt the temperature of the molecular cloud (10–30 K). This prescription conserves mass, that is, .

2.1.2. ISM scaling relations

We assume different scaling relations for the two density regimes: (i) for non-selfgravitating clouds, we adopt the scaling relations found for galactic Hi by Quiroga (1983): ρ_cl ∝ l^-2, v_turb,cl ∝ l^1/3, and thus v_turb,cl = v_turb(ρ_cl/ρ)^{− 1/6}, where v_turb and ρ are the turbulent velocity and the density of the disk, respectively. Since the minimum density considered in this work is 100 cm^-2, the maximum turbulent velocity of diffuse clouds is ~v_turb/ 2 ~ 5 km s^-1; (ii) for self-gravitating clouds, we adopt the scaling relations of Lombardi et al. (2010): ρ_cl ∝ l^-1.4, v_turb,cl ∝ l^1/2. As described in Sect. 2.7, the scale of the largest self-gravitating clouds l_cl is smaller than the turbulent driving length scale l_driv by a factor δ = l_driv/l_cl. We assume that the turbulent velocity dispersion of the largest self-gravitating clouds of density ρ_sg is , where v_turb is the velocity dispersion of the disk. For the assumed value of δ = 5 (see Sect. 2.7), this yields a velocity dispersion of the largest self-gravitating clouds of v_turb,cl ~ 5 km s^-1, which is consistent with observations (e.g., Solomon et al. 1987). We thus obtain for clouds with ρ_cl>ρ_sg. Alternatively, we assume ρ_cl ∝ l^-1 and (Solomon et al. 1987).

It turned out that the different scaling relations for self-gravitating clouds result in consistent molecular line luminosities within ~10%, except for the HCN emission of local spirals where the difference is ~20%. The models including the Solomon et al. (1987) scaling always reproduce observations slightly better. Therefore, the following results are based on the Solomon et al. (1987) scaling.

2.1.3. Gas and dust temperatures

Neufeld & Kaufman (1993) and Neufeld et al. (1995) considered the radiative cooling of fully shielded molecular astrophysical gas over a wide range of temperatures (10 K ≤ T_g ≤ 2500 K) and H₂ densities (10³ cm^-3 ≤ n(H₂) ≤ 10¹⁰ cm^-3). Their model for the radiative cooling of molecular gas includes a detailed treatment of the interstellar chemistry that determines the abundances of important coolant molecules, and a detailed treatment of the excitation of the species H₂, CO, H₂O, HCl, O₂, C, O, and their isotopic variants where important.

For simplicity, we only take the main cooling agents, CO, H₂, and H₂O, into account. We assume CO and H₂O abundances of x_CO = 10^-4(Z/Z_⊙) and x_H2O = 10^-6(Z/Z_⊙). According to Fig. 2 of Neufeld et al. (1995), we may underestimate, in this way, the cooling rates by approximately a factor of 2. However, for densities n(H₂) > 10⁵ cm^-3 and low temperatures (T ~ 20 K), the discrepancy increases up to a factor 3–4.

Neufeld & Kaufman (1993) and Neufeld et al. (1995) defined the molecular cooling rate as Λ_g = Ln(H₂)n(M), where n(H₂) and n(M) are the H₂ and coolant particle densities. The rate coefficient L depends on n(H₂), the gas temperature T, and (8)where g = 1 is a dimensionless geometrical factor and dv_turb/ dl is the turbulent velocity gradient. Neufeld & Kaufman (1993) and Neufeld et al. (1995) provided an analytical expression for L (Eq. (5) of Neufeld & Kaufman 1993) as a function of a set of parameters which depend on n(H₂), T, and . Since for each model cloud, n(H₂), T, and are known, we calculated the molecular line cooling Λ_g by interpolating the tabulated values of this parameter set.

To investigate the differences between our gas cooling and that proposed by Goldsmith (2001), we calculated these quantities for light and massive self-gravitating molecular clouds of different sizes, densities, column densities, and velocity dispersions: (i) light clouds with M_cl = 10⁴M_⊙: l_cl = ζ^-110 pc, n_cl = ζ³ 380 cm^-3, N_cl = ζ² 10²² cm^-2, and km s^-1 (plus signs in Fig. 1); and (ii) massive clouds with M_cl = 6 × 10⁴M_⊙: l_cl = ζ^-110 pc, n_cl = ζ³ 2350 cm^-3, N_cl = ζ² 7 × 10²² cm^-2, and km s^-1 (triangles in Fig. 1) with 1 ≤ ζ ≤ 10.

Our simplified cooling prescription is in good agreement with that of Goldsmith (2001) (~0.2 dex) for the less massive clouds (Fig. 1). For the more massive clouds, our cooling prescription gives values up to a factor 4 (0.6 dex) higher than those of Goldsmith (2001) for the highest cloud densities. Overall, the ratio between our cooling and that of Goldsmith (2001) is approximately 0.3 dex. Since the dependence of cooling on temperature is approximately Γ ∝ T^2.5−3.0 (Goldsmith 2001), the corresponding uncertainty on the gas temperature is 0.24 dex (a factor of 1.7) at most and 0.12 dex (a factor of 1.3) overall.

We thus conclude that our cooling prescription and that of Goldsmith (2001) are the same within a factor of 2. We prefer to use the reduced Neufeld line cooling instead of the cooling function proposed by Goldsmith (2001), because it takes into account the cloud column density and velocity dispersion.

For the calculation of the thermal balance within molecular clouds, one needs to consider processes affecting the gas and the dust in addition to the radiative gas cooling discussed above.

We assume gas heating via turbulence and cosmic rays: (9)Photoelectric heating by UV photons within photodissociation regions is neglected because the local FUV field plays a minor role for the CO luminosity of a giant molecular cloud (Wolfire et al. 1993; see Sect. 2.11). The factor of in the expression for the turbulent heating is somewhat lower than the factor of 0.42 advocated by Mac Low (1999). Following Nelson & Langer (1997), the constant η is chosen such that for Γ_CR = 6.4 × 10^-28(n(H₂)/cm^-3) erg cm^-3 s^-1. Furthermore, η includes the attenuation factor (Σ/(0.9 M_⊙ pc^-2))^0.021exp(−Σ/(9 × 10⁴M_⊙ pc^-2)) described by Padovani & Galli (2013) and the CR advection by a galactic wind (1 for local spiral and high-z star-forming galaxies; 1 for non-self-gravitating clouds and 140 for the self-gravitating clouds in ULIRGs; 140 for smm-galaxies; see Sect. 2.12).

The dust is heated by the interstellar UV and optical radiation field: (10)where n_d is the density of dust grains and σ_d the absorption cross section of a grain. Following Goldsmith (2001), we set n_dσ_d = 7.4 × 10^-22(n(H₂)/cm^-3) cm^-1. The ratio between the interstellar UV/optical and total radiation field is assumed to be (11)where F₀ = 5.3 × 10^-3 erg cm^-2 s^-1 (Goldsmith 2001). We assume that the UV radiation is emitted by young massive stars whose surface density is proportional to the star-formation rate per unit area . The optical light stems from the majority of disks stars (Mathis et al. 1983; Draine 2011) whose surface density is Σ_∗. The normalizations of and Σ_∗ are set by observations of the ISRF at the solar radius: 6.7% of the total stellar light is emitted in the UV (Mathis et al. 1983; Draine 2011). This implies that the local Galactic star-formation rate is , which is (i) approximately a factor of two lower than the value given by Kennicutt & Evans (2012); (ii) consistent with the local star-formation rate at ~0.75 × R₂₅ in the sample of nearby spiral galaxies of Leroy et al. (2008); and (iii) a reasonable value for a gas disk at R = 8 kpc with v_rot = 200 km s^-1, Q ~ 3, and Ṁ = 0.2M_⊙ yr^-1. Furthermore, we allow for an additional factor k, which plays the role of U_min in the Draine & Li (2007) models. We set k = 1 for all galaxies except the local spirals where k = 2. This additional factor is (i) needed to reproduce the observed infrared spectra energy distributions and (ii) consistent with the distribution found for nearby spiral galaxies by Dale et al. (2012).

In the presence of dust and gas, the interstellar radiation field is attenuated. For this attenuation we adopted the mean extinction of a sphere of constant density (12)with τ = (Z/Z_⊙) Σ/(15 M_⊙ pc^-2). For high optical depths, Eq. (12) becomes , that is, the attenuation decreases to very low values. However, when the molecular clouds become optically thick in the near-infrared (at τ_V ~ 10), radiative transfer effects become important; a significant infrared radiation field builds up which heats the dust until deep into the molecular clouds. To take this additional heating term into account, we set I/I(0) = 0.246 if 3 (τ^-1−2 τ^-2 + 2 τ^-3)−6 exp(−τ) τ^-3< 0.246 (see Appendix A).

The model dust temperature in the absence of collisional dust heating for the local Galactic interstellar radiation field with I/I(0) is T_dust = 18.4 K. This temperature lies between the equilibrium temperature of silicate (16.4 K) and graphite (22.3 K) for the local Galactic interstellar radiation field (Draine 2011).

The expression for the radiative heating of dust grains yields (13)We assume the dust mass absorption coefficient of the following form: (14)with λ₀ = 250 μm, κ₀ = 0.48 m² kg^-1 (Dale et al. 2012), and a gas-to-dust ratio of (including helium; Rémy-Ruyer et al. 2014). Our gas-to-dust ratio is a factor of 1.5 lower than the solar gas-to-dust ratio (Sofia & Meyer 2001a,b). We set the slope β = 1.5 for the local spiral (Dale et al. 2012) and high-z star-forming galaxies, and β = 2.0 for the ULIRGs (Klaas et al. 2001) and smm-galaxies. Adapting the dust cooling rate of Goldsmith (2001) yields (15)Following Goldsmith (2001), the dust cooling energy transfer between dust and gas due to collisions is (16)where ΔT = T_g−T_d K.

To determine the thermal balance of gas and dust, coupled together by the gas-dust collisions, we solve the following equations simultaneously:(17)and (18)

2.1.4. CO, HCN, and HCO⁺ abundances from chemical network

Chemical modeling is carried out using the Nautilus gas-grain code presented in detail in Hersant et al. (2009), Semenov et al. (2010), and Ruaud et al. (2015). This code computes the abundances of chemical species (atoms and molecules) as a function of time by solving the rate equations for a network of reactions. For gas-phase reactions, we use the kida.uva.2014 network (Wakelam et al. 2015)¹ comprising 489 species and 6992 reactions. For grain surface reactions, we use the desorption, diffusion, activation barrier energies along with a set of grain surface reactions, all from the KIDA database. Both, thermal and non-thermal desorption processes are taken into account, the latter consisting mainly of CR-induced desorption following the formalism presented by Hasegawa et al. (1993).

The model parameters are time, density, gas temperature, grain temperature, UV flux, cosmic ray ionization rate, and the elemental abundances of the elements C, O, and N (C/H = 1.7 × 10^-4, O/H = 2.4 × 10^-4, N/H = 6.2 × 10^-5).

Grids of models were obtained by varying the cloud lifetime (20 log spaced steps between 10³ and 10⁸ yr), the cloud density (20 log spaced steps between 10³ and 10⁸ cm^-3), and cloud gas temperatures (20 log spaced steps between 10 and 300 K for local spirals and 10 and 800 K for ULIRGs). For each type of cloud, the CO, HCN, and HCO+ abundances were interpolated on the grid given the lifetime, density, and temperature of the cloud.

The ISM chemistry also depends on the dust temperature and the cosmic ray ionization rate, which were assumed to be constant for each galaxy sample:

• local spiral galaxies: T_dust = 15 K, ζ_CR = 3 × 10^-17 s^-1,

• ULIRGs, smm, high-z star-forming galaxies: T_dust = 30 K, ζ_CR = 1.3 × 10^-15 s^-1.

Testing of the influence of these parameters on the molecular line emission showed that the dust temperature plays a minor role. For the choice of the CR ionization rates, we refer to Sect. 2.12. To take the gas metallicity into account, all abundances are multiplied by (Z/Z_⊙).

2.2. Supernova-driven turbulence

First, we assume that the gas is turbulent, so that the turbulent velocity is the relevant one throughout the disk (making the exact temperature of the gas largely irrelevant; for simplicity, we assume a constant sound speed of c_s = 6 km s^-1 for the warm neutral medium). We assume that this turbulence is driven by SNe and that they input their energy in turbulent eddies that have a characteristic length scale, l_driv, and a characteristic velocity, v_turb. This driving length scale may be the characteristic length scale of a SN bubble, but it does not have to be so. It may be set by the interaction of multiple SN bubbles or of a SN with the surrounding ISM. We note that based on simulations, the assumption of a single driving scale may be a simplification (Joung & Mac Low 2006). The VB03 model does not address the spatial inhomogeneity of the turbulent driving nor the mechanics of turbulent driving and dissipation. It is assumed that the energy input rate into the ISM due to SNe, Ė_SN, is cascaded to smaller scales without losses by turbulence. At scales smaller than the size of the largest self-gravitating clouds, the energy is dissipated via radiation from cloud contraction and star formation. We refer to Mac Low & Klessen (2004) for a review of these topics. We limit our analytical model to the first energy sink, which is the scale where the clouds become self-gravitating.

We can connect the energy input into the ISM by SNe directly to the star-formation rate. With the assumption of a constant initial mass function (IMF) independent of environment one can write (19)where ΔA is the unit surface element of the disk and the CO disk thickness is assumed to be l_driv. The gas disk viscosity is defined as ν = l_drivv_turb (VB03 and Sect. 2.6). Following Vollmer & Beckert (2001), the turbulent energy dissipation rate is . The turbulent dissipation timescale is (20)This result is in agreement with numerical simulations of turbulence that show a decay of turbulence on an approximately crossing timescale (e.g., Stone et al. 1998; Mac Low et al. 1998).

The factor of proportionality ξ relates the local SN energy input to the local star-formation rate and is assumed to be independent of local conditions. ξ is normalized using Galactic observations by integrating over the Galactic disk and results in ξ = 4.6 × 10^-8 (pc/yr)² (see VB03). The adopted energy that is injected into the ISM is erg based on numerical studies by Thornton et al. (1998). The final two parts of Eq. (19) assume that stars form over a characteristic scale equal to the driving length and equate energy output from SNe with the energy transported by turbulence (see VB03).

In the outer galactic disk, where the star-formation activity is very low, turbulence can be maintained by the energy gained via accretion within the gravitational potential of the galaxy (Vollmer & Beckert 2002). In this case, the energy injection rate is (21)This energy source represents an addition to the model described in Vollmer & Leroy (2011). The total energy injection rate is (22)

2.3. Star formation in molecular clouds

Second, we assume that stars form out of gravitationally bound clouds. We take the local gravitational free-fall time, given by (23)to be the relevant timescale for star formation. Here, G is the gravitational constant and ρ_cl the density of a single cloud.

Cloud collapse, and thus star formation can only proceed if enough molecules form during the cloud collapse to allow the gas to continue cooling². Vollmer & Leroy (2011) assumed that, where the timescale for H₂ formation is long (compared to the free-fall time), the relevant timescale for star formation is the H₂ formation timescale. Since we aim at reproducing the ISM properties within the optical disks, we decided not to include this complication in our present model.

2.4. Molecular fraction

We follow two approaches to calculate the fraction of molecular hydrogen: (i) based on the lifetimes of the molecular clouds and (ii) based on photodissociation of molecules by the interstellar radiation field. In Sect. 6.3, we show that both formalisms lead to consistent results for the molecular fraction.

2.4.1. Molecular fraction based on the lifetimes of the molecular clouds

This approach assumes that molecular clouds are relatively short-lived, appearing and disappearing over approximately a free-fall time (equivalently, by our construction, a turbulent crossing time); otherwise they might reach chemical equilibrium even when the H₂ formation time is long compared to the free-fall time. Accordingly, we estimate the molecular ratio in the disk from the ratio of a cloud lifetime (i.e., the crossing or free-fall time) to the H₂ formation time scale: (24)The molecular fraction is (25)We take the characteristic time to form H₂ out of H to be approximately (26)where α is a coefficient that depends on the gas phase metallicity and temperature (Draine & Bertoldi 1996) and ρ_cl is the density of a single cloud.

The coefficient of the molecular formation timescale α₀ is assumed to be metallicity dependent (Tielens & Hollenbach 1985). Because we admit external gas accretion, the metallicity of the star-forming ISM mainly depends on the ratio of accretion to star formation rate a. Small a< 1 lead to a metallicity derived from a closed box model, whereas in the case of a> 1, the metallicity is equal to the true stellar yield y_true (Köppen & Edmunds 1999). For gas fractions higher than 0.04, the difference between the two solutions is less than a factor of two. Moreover, Dalcanton (2007) showed that the effective yield y_eff = Z_gas/ ln(1 /f_gas), where Z_gas is the gas metallicity and f_gas the gas fraction, for disk galaxies with a rotation velocity higher than 100 km s^-1 is approximately constant, that is, for these galaxies, a closed box model can be applied. We thus feel confident to estimate the gas phase metallicity based on a closed box model using the gas fraction: (27)where Σ_∗ is the stellar surface density and α₀ = 3.6 × 10⁷ yr M_⊙ pc^-3. Adopting a stellar and gas surface density of Σ_∗ = 40M_⊙ pc^-2 and Σ_gas = 10 M_⊙ pc^-2 at the solar radius of the Galaxy yields α_⊙ = 2.2 × 10⁷ yr M_⊙ pc^-3, which corresponds to the value used by Hollenbach & Tielens (1997). Within this framework, the metallicity is (28)It turned out that the gas phase metallicities of the ULIRG and high-z star-forming galaxy samples are underestimated by up to a factor 10 with our simple closed-box model (see Sect. 6.1). To remedy the situation, we adopted the following heuristic recipe for all galaxies: (29)A possible explanation for this recipe is accretion of pre-enriched gas onto or into the galactic disks in which case the closed-box model underestimates the metallicity. Less gas depletion in starburst galaxies might also play a role. Since we expect starburst galaxies to host galactic winds (see, e.g., Veilleux et al. 2005), the ejection of metals due to these outflows (leaky box model) is assumed to be much less significant than the addition of metals from accretion.

2.4.2. Molecular fraction based on photodissociation

For the determination of the H₂ column density of a gas cloud, we take into account (i) photo-dissociation of H₂ molecules and (ii) the influence of the finite cloud lifetime on the H₂ formation.

For the photo-dissociation of H₂ molecules, we follow the approach of Krumholz et al. (2008, 2009). These authors solved the idealized problem of determining the location of the atomic-to-molecular transition in a uniform spherical gas cloud bathed in a uniform isotropic dissociating radiation field. It is assumed that the transition from atomic to molecular gas occurs in an infinitely thin shell. The cloud has a constant inner molecular and outer atomic gas density. The inner molecular core and the outer atomic shell are assumed to be in thermal pressure equilibrium. The atomic gas density is taken to be the density of the cool neutral medium (Eq. (4)). The H₂ to Hi ratio is (30)with s = (Σ_cl/ 1 M_⊙)(Z/Z_⊙)/(4 τ_Hi). The Hi optical depth is (31)with the dimensionless radiation field strength χ, which we set to . Here, we assume a constant ratio between the inner molecular and outer atomic gas density, which is of the order of 10. For and solar metallicity, the transition between a molecular- and atomic-dominated cloud occurs at Σ_cl ≃ 20M_⊙ pc^-2. The H₂ fraction of the cloud is f_H2 = R_H2/ (1 + R_H2). This treatment insures a proper separation of Hi and H₂ in spiral galaxies, that is, clouds of low density (~100 cm^-3) and low column density (~10²¹ cm^-2) are fully atomic, whereas clouds of high density, that is, GMCs, (≥1000 cm^-3) and high column density (≥10²² cm^-2) are fully molecular. In starburst regions (e.g., in ULIRGs), where gas densities and surface densities are much higher, this treatment has no effect, since the gas will be fully molecular.

In a second step, we take into account the molecular fraction due to the finite lifetime of the gas cloud . The total molecular fraction of a cloud is . The molecular fraction due to the finite lifetime has the highest influence on f_mol at large galactic radii.

We now go from the H₂ mass fraction to the CO mass fraction. In an externally irradiated gas cloud, a significant H₂ mass may lie outside the CO region, that is, it is dark in the outer regions of the cloud where the gas phase carbon resides in C or C⁺. In this region, H₂ self-shields or is shielded by dust from UV photodissociation, whereas CO is photodissociated. Following Wolfire et al. (2010), the dark gas mass fraction for a cloud of constant density is (32)with (33)and A_V = 2 (Z/Z_⊙)N_cl/ (1.9 × 10²¹ cm^-2) where N_cl is the H₂ column density. The CO mass fraction is then .

Since the attenuation of the UV radiation field leading to Eq. (32) is mainly caused by dust, we expect HCN to survive everywhere where the ISRF is attenuated enough to permit a high CO abundance. Thus, the HCN abundance should approximately follow the CO abundance, unless there is a strong X-ray/cosmic ray flux that is not attenuated by dust. In the absence of a proper theoretical model for the HCN dissociation, we thus assume the same dissociation rate for HCN as for CO.

2.5. Vertical disk structure

In the model, the disk scale height is determined unambiguously by the assumption of hydrostatic equilibrium and the turbulent pressure (Elmegreen 1989): (34)where ρ is the average density, v_turb the gas turbulent velocity in the disk, the stellar vertical velocity dispersion, and Σ the surface density of gas and stars. The stellar velocity dispersion is calculated by , where the stellar vertical height is taken to be H_∗ = l_∗/ 7.3 with l_∗ being the stellar radial scale length (Kregel et al. 2002). We neglect thermal, cosmic ray, and magnetic pressure.

2.6. Treatment as an accretion disk

The turbulent motion of clouds is expected to redistribute angular momentum in the gas disk, like an effective viscosity would do. This allows accretion of gas towards the center and makes it possible to treat the disk as an accretion disk (e.g., Pringle 1981). This gaseous turbulent accretion disk rotates in a given gravitational potential Φ with an angular velocity , where R is the disk radius. The disk has an effective turbulent viscosity that is responsible for mass accretion and outward angular momentum transport. In this case, the turbulent velocity is driven by SN explosions, which stir the disk and lead to viscous transport of angular momentum. In addition, star formation removes gas from the viscous evolution. Following Lin & Pringle (1987), the evolution of the gas surface density is given by (35)where ν is the gas disk viscosity, Ω the angular velocity, and is the external mass accretion rate. In contrast to Lin & Pringle (1987), we assume a continuous and non-zero external gas mass accretion rate.

Forbes et al. (2014) presented an analytical approach based on Eq. (35). They showed that galaxies tend to be in a slowly evolving equilibrium state wherein new accretion is balanced by star formation, galactic winds, and radial transport of gas through the disc by gravitational instability-driven torques. For a stationary gas disk in such an equilibrium, where star formation is balanced by external accretion, the local mass and momentum conservation together with yield: (36)where Ṁ is the mass-accretion rate within the disk. In the absence of external mass accretion, the gas disk can be assumed to be stationary as long as the star-formation timescale t_∗ exceeds the viscous timescale t_ν = R²/ν. For and t_∗<t_ν, Eq. (36) is not valid. In this case the gas disk is rapidly turned into stars within the gas consumption time (2 Gyr, Evans 2008). Since most spiral galaxies still have a significant amount of gas, we think that spiral galaxies are generally not in this state. Solving the time dependent Eq. (35) is beyond the scope of this work and we apply Eq. (36).

The viscosity is related to the driving length scale and characteristic velocity of the SN-driven turbulence by ν = v_turbl_driv (VB03). Because the lifetime of a collapsing and star-forming cloud () is smaller than the turnover time of the large-scale eddy (l_driv/v_turb), the turbulent and clumpy ISM can be treated as one entity for the viscosity description.

2.7. Clumpiness

A critical factor in the model is the relationship between the density of individual clouds, ρ_cl, and the average density of the disk, ρ. It is the density of individual clouds that is relevant to the timescale for star formation. In this model, the two are related by the volume filling factor, Φ_V, so that .

Here, ρ_cl refers to the density of the largest self-gravitating structures in the disk, so that for these structures, the turbulent crossing time and gravitational free-fall time are equal. The scale of such a cloud, l_cl, is smaller than the driving length scale, l_driv, by a factor δ, which we do not know a priori. Following Vollmer & Leroy (2011), we set δ = 5.

Shear, due to differential galactic rotation, could stabilize clouds, modifying the timescale for collapse. However, this effect is mainly important when the ratio of the cloud to disk surface density is lower than the ratio of cloud to disk velocity dispersion, which is not the case over most of the disk in a typical spiral. Typical GMC surface densities are ~200 M_⊙ pc^-2 (Solomon et al. 1987), whereas disk surface densities only exceed 100M_⊙ pc^-2 in the very center of spiral galaxies (Leroy et al. 2008).

We can calculate the turbulent timescale for the cloud, , for a fractal ISM: (37)where D is the fractal dimension (see, e.g., Frisch 1995) of the ISM. We assume D = 2 for a compressible, self-gravitating fluid, which is close to the findings of Elmegreen & Falgarone (1996). Once δ and thus are specified, we can solve for the density of the corresponding scale by setting . The volume filling factor is then defined by comparing ρ_cl and ρ. Once the volume filling factor is known (from δ or l_cl), we can calculate the local star-formation rate, , via (38)where is the local free-fall timescale determined by corresponds to the contraction timescale (Ostriker et al. 1999) of clouds of constant density in Virial equilibrium. Since, in our model, the lifetime of a cloud is the free-fall as suggested by Ballesteros-Paredes & Hartmann (2007), this implies that during the cloud lifetime, approximately % of the cloud mass turns into stars.

The vertically integrated star-formation rate in the inner disk where is (39)that is, it is the mass flux density of the turbulent ISM into the regions of star formation.

2.8. Thermal dust emission

The dust temperature T_d of a gas cloud of given density and size illuminated by a local mean radiation field is calculated by solving Eq. (18). With the dust mass absorption coefficient of Eq. (14), the dust optical depth is (40)where Σ_cl is the cloud surface density in g/cm². The infrared emission at a given wavelength at a given galactic radius R is calculated in the following way: (41)where B(λ,T_d) is the Planck function and Φ_A = 1.5 (ΔM/M) (Σ/Σ_cl) the area filling factor. The factor 1.5 takes into account that the mean cloud surface density is 1.5 times lower than the surface density in the cloud center Σ_cl = ρ_cll_cl. The integration of Eq. (41) yields the total infrared emission at a given galactic radius R: (42)At a given wavelength, λ, the effective background temperature of the thermal dust emission T_{eff dust} is determined by I_dust(λ) = B(λ,T_{eff dust}). We do not subtract the cosmic infrared background, because it is always much smaller than the dust emission within the region of interest. The total infrared luminosity is given by (43)To reproduce the observed total infrared luminosity, it is essential to take into account the diffuse warm neutral medium, which is not taken into account for the molecular line emission, because the gas is in atomic form. We do so by explicitly calculating the dust infrared emission based on the proper dust temperature (Eq. (18)), density of ρ/ 2, area filling factor , and gas mass fraction of (ΔM/M)_WNM = (1−(ΔM/M)_{CNM + mol}).

2.9. Molecular line emission

A molecular line source is usually observed by chopping the telescope’s beam between on- and off-source positions and measuring the difference in antenna temperatures. In general, the difference in brightness temperatures is (44)where τ is the optical depth of the line, ν the frequency of the observations, h and k the Planck and Boltzmann constants, and T_ex and T_bg the excitation and background brightness temperatures, respectively.

Considering only a single collider (H₂) for simplicity, the excitation temperature is (45)where T_∗ = hν_ul/k, n is the gas density, nq_ul the collisional de-excitation rate, and A_ul the Einstein coefficients of the transition ul. The background brightness temperature T_bg is the sum of the effective emission temperatures of the galaxy’s dust T_{eff dust} and the cosmic background at the galaxy redshift T_CMB (see Eq. (17) of da Cunha et al. 2013). For optically thin transitions, the ratio of the radiative and collisional rates is just the ratio of the density to the critical density for the transition (46)We use T_CMB = 2.73 K for the local spiral galaxies and ULIRGs, T_CMB = 6.0 K for the high-z star-forming galaxies, and T_CMB = 8.19 K for the submillimeter galaxies. This corresponds to a cosmic microwave background of T_CMB = 2.73 (1 + z) K (see, e.g., Carilli & Walter 2013) and mean redshifts of ⟨ z ⟩ = 0, 1.2, 2, respectively.

For optically thick transitions, the upper-level population can be enhanced due to absorption of line photons, leading to excitation temperatures higher than those expected simply due to H₂ collisions, since the line photons emitted upon spontaneous decay cannot easily escape the cloud. This so-called radiative trapping of the line photons builds up the radiation field at the frequency of the line, leading to enhanced excitation of the upper state via photon absorption. The escape probability formalism can be used to treat this optically thick situation (see, e.g., Scoville 2013).

This formalism is applicable to situations in which systematic velocity gradients are large compared to the small-scale thermal motions. The line photons from one region of the cloud are then incoherent with other regions due to the Doppler shift; they can then only interact with molecules in the local region near where they were emitted. In the photon trapping regime, the spontaneous decay rates (A_ul) and thus the critical density (n_crit) used in analyzing the equilibrium molecular excitation are reduced by a factor β equal to the effective probability for escape of line photons from the emission region (Scoville & Solomon 1974; Goldreich & Kwan 1974). For a spherical cloud of uniform density, Draine (2011) gives (47)The critical density in the optically thick case is then (48)For our analytic analysis, we follow Scoville et al. (2015) and use the sum of the collision rate coefficients out of the upper level J to any other rotational level (both below and above J) since all of these transitions couple the level to the gas kinetic temperature.

For the determination of the optical depth of a molecular emission line, we follow Draine (2011). The line-center optical depth, from cloud center to edge, for a transition from level J + 1 to level J is (49)where r_cl = l_cl/ 2 is the cloud radius, λ the wavelength of the observations, and g_J = 2J + 1 the transition weights. Following Draine (2011), we adopt the following expression for the CO line optical depth: (50)where n₃ = n/ (10³ cm^-3) and R₁₉ = r_cl/ (10¹⁹ cm). The fraction of molecules of species X in a given rotational level is (51)where B₀ is the rotation constant of a molecule of species X.

In summary, we consider two-level molecular systems in which the level populations are determined by a balance of collisions with H₂, spontaneous decay and line photon absorption, and stimulated emission with τ> 1. Our final expression for the CO line optical depth reads (52)where the rotation temperature is B₀/k = 2.77 K. We use a normalization, which is different from that of Draine (2011), because we assume the canonical x(CO) = 10^-4. In a second step, the HCN abundances are calculated using a chemical network (see Sect. 2.1.4).

For simplicity, we neglected the hyperfine structure of HCN. In this simplified treatment, we can write (53)where we assumed a HCN abundance of x(HCN) = 2 × 10^-8. In a second step, the HCN abundances are calculated using a chemical network (see Sect. 2.1.4). The rotation constant (B₀/k = 2.13 K) and Einstein coefficients are those of HCN. In the present work, we only investigate the HCN(1–0) transition.

The rotation constants, Einstein coefficients, and collision rates were taken from the Leiden Atomic and Molecular Database (LAMDA; Schöier et al. 2005). The CO collision rates were provided by Yang et al. (2010). The HCN collision rates were taken from the He–HCN rate coefficients calculated by Dumouchel et al. (2010), scaled by a factor of 1.36 to go to HCN–H₂ (see Green & Thaddeus 1976).

The HCO⁺ emission was calculated in the same way.

We verified that the model brightness temperature is consistent (within 10-20%) with the brightness temperature calculated by RADEX (van der Tak et al. 2007) for densities, gas kinetic temperatures, column densities, and linewidths typical for giant molecular clouds.

2.10. HCN infrared pumping

HCN has a large dipole moment and therefore does not trace dense gas if there is another excitation mechanism that is faster than the H₂ collisions and independent of gas density. One such excitation path is through a vibrationally excited state, to which molecules can be pumped by infrared radiation (Carroll & Goldsmith 1981). The first vibrationally excited state of HCN is its bending state (v₂ = 1) 1024 K above the ground with an emitting wavelength of λ = 14 μm (Sakamoto et al. 2010). Following Sakamoto et al. (2012), we define an equivalent gas density (54)where T₀ = 1024 K corresponds to the energy gap between the two vibrational levels v = 0 and 1, A_vib = 3.7 s^-1 is the Einstein coefficient for the vibrational transition, and γ_J,J−1 is the collisional rate coefficient. As already stated in Sect. 2.9, we follow Scoville et al. (2015) and use the sum of the collision rate coefficients out of the upper level J to any other rotational level (both below and above J) since all of these transitions couple the level to the gas kinetic temperature. T_vib is the equivalent blackbody temperature at λ = 14 μm of the local and global background radiation. HCN IR-pumping is implemented in the model by replacing the cloud density n_cl by n_equiv if n_equiv > n_cl in the HCN emission calculations (Sect. 2.9).

HCO⁺ has a similar vibrational bending state at λ = 12 μm. The associated excitation through radiative pumping is less than half as significant as that provided by the HCN molecule (Imanishi et al. 2016). For this work, we did not take IR-pumping of HCO⁺ into account.

2.11. Photon-dominated regions

Photodissociation regions (PDRs) are regions of a gas cloud where the physics and chemistry is dominated by penetrating FUV photons. The structure of PDRs can be described by a plane-parallel, semi-infinite slab illuminated by an intense FUV field G₀, measured in units of the Habing (1968) interstellar radiation field (=1.6 × 10^-3 erg cm^-2 s^-1). At the surface of the cloud, an atomic surface layer is created by the incoming FUV photons. The transition from atomic to molecular hydrogen occurs at the depth approximately A_V ~ 2. As the FUV photons are attenuated by the dust, the phase of carbon shifts from C⁺ to C and CO at A_V ~ 4. Deep inside the cloud (A_V > 10), HCN and HCO⁺ are formed.

Wolfire et al. (1993) found that the integrated CO(1–0) luminosity of giant molecular clouds increases by only ~10% between G₀ = 1 and G₀ = 100 models. The luminosities are similar because the higher incident FUV field forces clumps to become optically thick in the CO(1–0) transition deeper into the cloud where dust extinction lowers the radiation field to a value near G₀ ~ 1. This results in similar gas temperatures in the optically thick clumps for both models. Thus, the local FUV radiation field plays a minor role as long as G₀> 1.

Loenen et al. (2008) subdivided PDRs into two types according to the cycle of star formation; UV-dominated high-density (n ≥ 10⁵ cm^-3) PDRs from deeply embedded young stars, and lower-density (n = 10^4.5 cm^-3) PDRs that are dominated by mechanical feedback from supernova shocks. Due to the short duty cycle of the evolutionary stage involving young stars compared to the second stage involving supernovae, most of the luminous infrared galaxies of their sample are observed to be in the later stage of their evolution. Since we can reproduce the multi-transition CO, HCN(1–0), and HCO⁺ luminosities without the inclusion of PDRs (Sect. 6), we suggest that the conclusion of Loenen et al. (2008) is also valid for the ULIRGs, smm, and high-z star-forming galaxies.

Within our model framework, PDRs are taken into account through the dissociation of molecules (Sect. 2.4.2), but not through their molecular line emission.

2.12. Galactic winds

Galactic winds are driven by multiple supernova explosions during a starburst or by AGN activity. The following statements are based on the review by Veilleux et al. (2005). The minimum star-formation rate that creates a galactic wind is Ṁ_∗ ~ 5M_⊙ yr^-1 or M_⊙ kpc^-2 yr^-1. Galactic winds remove mass and, if they are magnetized, angular momentum from the galactic gas disk. Mass outflow rates range between 0.1M_⊙ yr^-1 and 10M_⊙ yr^-1, with a trend for the rate to increase with increasing star-formation rate. The mass-to-light conversion factors are highly uncertain for galactic winds. The ratio between the mass outflow and star-formation rate varies between 0.01 and 10.

By writing νΣ = Ṁ/ (2 π) (Eq. (58)), we ignored galactic external accretion and outflows (wind). Including galactic winds, Eq. (35) becomes (55)Within the framework of our equilibrium model, we thus assumed . This might not be true for the compact starbursts as the ULIRGs. In this case, the constant Ṁ/ (2 π) links the turbulent dissipation timescale t_turb to the star-formation rate per unit surface (combining the mass, momentum, and energy conservation equations of Eq. (58)): (56)Galactic winds also remove cosmic ray particles from the galactic disk. Following Suchkov et al. (1993; see also Papadopoulos 2010), the CR energy densities in starburst galaxies driving a galactic wind scale with respect to that of the Galaxy as (57)where v_diff is the diffusion velocity at which CRs escape from quiescent disks such as the Milky Way while v_wind is the velocity of a star-formation-induced wind at which CRs are advected out of the star-forming regions. The diffusion velocity is set by the Alfvén velocity , where B is the magnetic field. Here, we assumed that the SN explosion rate is proportional to the star-formation rate. For the chemical network (Sect. 2.1.4), we assumed for the ULIRGs, smm, and high-z star-forming galaxies. Since our model yields for the ULIRGs and smm-galaxies, the mean ratio is . Assuming v_diff = 10 km s^-1 gives a mean wind velocity of v_wind ~ 1500 km s^-1. This represents approximately four times the escape velocity, two times the observed mean velocity of a molecular wind (Feruglio et al. 2015; Sakamoto et al. 2014), and is close to terminal wind velocities (e.g., Veilleux et al. 2005).

For the ULIRG model, we decided to apply the full CR heating to the non-self-gravitating clouds. This modification did not change the CO flux luminosities, but increased the HCN(1–0) luminosities by a factor of ~1.5 putting them closer to the observed luminosities. This treatment implies that (i) the initially spherical wind escapes mostly vertically from the star-forming regions without touching many of the non-self-gravitating clouds (“champagne effect”); (ii) the SN wind that blows in the direction of the galactic disk is absorbed by the non-self-gravitating clouds; and (iii) the solid angle around the star-forming region occupied by non-self-gravitating clouds is substantial (i.e., a wind opening angle ≳90^⊙). Further testing of this hypothesis on local starburst galaxies is needed to clarify the situation.

In the present model, the assumed CR ionization rate of the high-z star-forming galaxies is consistent with an absence of a galactic wind. In a future project, we will investigate the influence of a lower CR energy density due to a wind on the chemical network and thus the molecular line emission of these galaxies.

3. Method

The outline of our model/method is presented in Fig. 2.

The VB03 model yields the following system of equations to describe a turbulent clumpy galactic accretion disk: (58)The meaning of the variables is given in Table 1. For the global comparison between the observed and the model radial profiles, we solve the set of equations given above numerically.

The free parameters of the analytical model are the Toomre parameter Q and the disk mass accretion rate Ṁ. For the local spirals, we set Q to values derived in Vollmer & Leroy (2011) (see Table B.1). For the ULIRGs, high-z star-forming galaxies, and submillimeter galaxies, we assume a constant Q = 1.5 (see Tables B.2, B.4, B.3). The mass accretion rate Ṁ is determined by the total star-formation rate of the galactic disks (see Eqs. (58)).

For the calculation of the infrared dust emission and molecular line emission from the galactic disk, we divide the ISM into two density bins (see Sect. 2.1.2): (i) non-self-gravitating clouds with densities ρ₂ equal or higher than that of the cool neutral medium: max(n_CNM,100 cm^-3) ≤ ρ₂ ≤ ρ/ Φ_V and (ii) self-gravitating clouds with densities ρ₁ ≥ ρ/ Φ_V. The calculation of the molecular emission is done in the following steps:

• 1.

  Calculation of the temperatures of the dust and cool neutralmedium T_CNM and the self-gravitating clouds T_cl according to Sect. 2.1.3. Even if the gas is not molecular at densities of 100 cm^-3, we use the CO, H₂, H₂O cooling function for model consistency, because we want to apply the same code for local spirals and ULIRGs. For clouds of lowest densities in local galaxies, this yields temperatures between 80 and 150 K, well in the range of observed CNM temperatures (Kulkarni & Heiles 1987; Dickey & Lockman 1990). In addition, in these clouds, almost all H₂ and CO is dissociated (Sect. 2.4.2) and the molecular line emission is very weak (Sect. 2.9) due to the low gas density in these low-density clouds. Therefore, the high uncertainty of the calculated temperature of these clouds does not affect the total molecular line emission of the galactic disk.

• 2.

  We assume different scaling relations for the two density regimes: (i) non-self-gravitating clouds and (ii) self-gravitating clouds (Sect. 2.1.2).

• 3.

  Within the two density regimes, the properties (size, column density, turbulent velocity, temperature) of clouds with over-densities 2^N with N = 1,2,... are calculated according to the scaling relations. The mass fraction of each density bin is calculated via the pdfs described in Sect. 2.1.1. This procedure insures mass conservation, that is, .

• 4.

  The fraction of molecular mass f_H2 and mass emitting in the molecular line f_CO is determined according to Sect. 2.4.2. In the absence of a detailed theory of HCN dissociation, we use the same dark gas fraction for HCN as for CO.

• 5.

  The molecular line emission () is calculated according to Sect. 2.9 using the density, size, and temperature of the clouds in each density bin.

• 6.

  The total dust infrared emission at a given galactic radius R is determined by Eq. (41), the TIR luminosity is calculated via Eqs. (42) and (43), and the total CO or HCN molecular line emissions are calculated in the following way: (59)where the factor 2.35 links the turbulent velocity to the linewidth. The cosmic and local dust infrared backgrounds are taken into account. The total flux density is calculated by integrating the radial profile (60)The molecular fraction is (61)The total molecular gas mass is (62)

The Hi surface density is Σ_HI = Σ−Σ_H2 and the total Hi mass is (63)The upper limit of the galactic radius is R₀ = 4.5 × R_∗ for the local spiral galaxies and high-z star-forming galaxies (van der Kruit & Searle 1982) and R₀ = 3.4 × l_∗ for the ULIRGs and submillimeter galaxies, where l_∗ is the exponential scalelength of the stellar disk. The latter normalization was chosen to reproduce the spatial extents of ULIRGs (Downes & Solomon 1998).

4. Uncertainties

All the different steps of our model calculations have associated uncertainties. In the following, we discuss the main uncertainties for each step.

• Analytical model (Sect. 2,Eq. (58)): the constant ξ that links the star-formation rate to the energy injection rate is calibrated to the Galaxy. We do not know if this calibration also holds for the densest starburst regions; for example, in ULIRGs. In the presence of outbreaking supernovae giving rise to a galactic wind, one might expect a lower value of ξ. Since the total and molecular gas masses are only weakly dependent on ξ, we expect an uncertainty of the order of a few 10% on the determination of the gas mass associated to the uncertainty of ξ. For our calculations, we regard ξ as constant. The mass accretion rate Ṁ is determined by the total star-formation rate which is typically uncertain to a factor of 1.5 (Leroy et al. 2008, Genzel et al. 2015). Due to an uncertain AGN contribution to the IR luminosity, the uncertainty is expected to be higher for ULIRGs. Since Ṁ_∗ ∝ Ṁ^{0.4 to 0.6} (Vollmer & Leroy 2011), we estimate the uncertainty associated to the determination of Ṁ from the integrated star-formation rate to be of the order of a few tens of percent. The uncertainty on the derived total gas mass is of the same order. Another uncertainty comes from the unknown velocity dispersion of the disk stars or the thickness of the galactic stellar disk. The applied relation H_∗/l_∗ = 7.3 for local spiral galaxies (Kregel et al. 2002) has a dispersion of ~30%. We estimate the uncertainty on the gas pressure and thus the gas surface density and total gas mass to be small for local spirals. However, it is expected to be higher for ULRIGs, high-z star-forming galaxies, and submillimeter galaxies that might deviate from this relation. The mass and momentum conservation equation Ṁ/ (2 π) = νΣ represents another source of uncertainty. This equation implies that the mass loss through star formation is balanced by radial or external gas accretion. In most cases, radial accretion within the disk is negligible. In the absence of this external accretion, we expect that the gas surface density in the central part of the galaxy decreases because of star formation. We take this phenomenon into account by increasing Q towards the galaxy center in some of the local spiral galaxies. For ULIRGs, high-z star-forming galaxies and submillimeter galaxies we assume a constant Q. This might be justified by a large gas mass and star-formation rate of these galaxies. Vollmer & Leroy (2011) have shown that the system of equations including the mass and momentum conservation equation describes the radial H₂, star formation, turbulent velocity, and molecular fraction profiles of local spiral galaxies in a satisfying way. For the metallicity, we use a simple closed-box model. As long as the star-formation timescale is smaller than the gas accretion timescale, this assumption is justified. In ULIRGs and submillimeter galaxies, this might not be the case, leading to an overestimation of the metallicity and thus the molecular abundances X. This will increase the molecular line emission ∝X^0.4 (Scoville & Solomon 1974). We are thus confident that the model equations describe the physics of a star-forming galactic disk within a factor of 1.5–2.

• 2.

  Density probability distribution functions (pdfs) (Sect. 2.1.1): the role of gas self-gravitation is not included in the pdf of Padoan et al. (1997). In the presence of self-gravitation, the shape of the density pdf is altered at high densities, that is, the mass fraction of high-density gas is higher in a pdf including self-gravitation than in a pdf without self-gravitation (see, e.g., Schneider et al. 2015). We estimate the uncertainty of the mass calculation due to different pdfs to be a factor of 2.

• 3.

  Gas and dust temperature calculations (Sect. 2.1.3): the main uncertainty comes here from the gas cooling function, which is assumed to be dominated by CO, H₂, and H₂O line cooling. Based on Neufeld et al. (1995), this leads to a possible underestimation of the total cooling by up to a factor of 3 to 4. On the other hand, the comparison with the cooling function proposed by Goldreich (2001) shows an overall discrepancy of a factor of 2. We believe that, on average, our cooling function is uncertain to a factor of ~2. The resulting gas temperature is uncertain by approximately a factor of ~1.3.

• 4.

  Molecular line emission (Sect. 2.9): sources of uncertainties are manyfold: (i) we only calculate a single transition at a time and not the full system of transitions of a given molecule; (ii) we assume spherical geometry of all clouds, whereas a sheet- or filament-like geometry is not excluded for non-self-gravitating clouds. This directly affects the area-filling factor Φ_A; (iii) the area-filling factor is small enough so that line emission of a smaller cloud is not absorbed by a larger cloud. This effect might only play a role in ULIRGs; (iv) we assume abundances that are proportional to the closed-box metallicity; we therefore neglect gas depletion and chemistry effects. We had to modify our simple closed-box model for starburst galaxies to take these effects into account (Eq. (29)). Of all these sources of uncertainty we estimate (iv) to be most important, leading to an uncertainty associated to the molecular line emission of a factor of a 2 for CO and of a few for HCN.

• 5.

  H₂ and CO dissociation (Sect. 2.4.2): we treat the photo-dissociated region in molecular gas clouds in a crude way. It is assumed that the outer Hi envelope has the density of the CNM, which has not to be the case. The dark gas fraction (H₂ without CO) also depends on the density of the cloud envelope or the density profile of the cloud. We estimate the uncertainty of the H₂ and CO-emitting mass calculation to be approximately a factor of 2.

• 6.

  Scaling relations (Sect. 2.1.2): whereas the scaling relations for self-gravitating GMCs are relatively robust, those for diffuse (non-self-gravitating) clouds are less well established. This will lead to higher uncertainties for galaxies with very dense gas disks, especially for ULIRGs, where the molecular line emission from the diffuse clouds dominates the total emission. We estimate the uncertainty on the calculation of the mass fraction at a given density due to the adopted scaling relations to be a factor of 2 for the ULIRGs and 30–50% for the high-z star-forming and submillimeter galaxies.

In summary, the most important uncertainties are due to the analytical model, the choice of the pdf, and the adopted molecular abundances. All uncertainties are of the order of a factor of 2. In addition, the uncertain scaling relations of ULIRGs might add an additional factor to the molecular line emission calculation of ULIRGs.

5. Galaxy samples

The input parameters of our model are: the exponential scale-length of the stellar disk l_∗, the stellar mass M_∗, the rotation velocity v_rot, and the star-formation rate Ṁ_∗. We apply our model to four galaxy samples for which these parameters are determined observationally in a uniform way: local spiral galaxies, local ULIRGs, high-z star-forming galaxies, and submillimeter galaxies. Following Boissier et al. (2003), for all galaxies, we assume a rotation curve of the form (64)where v_flat and l_flat represent the velocity at which the rotation curve is flat and the length scale over which it approaches this velocity, respectively. Moreover, we assumed δ = 5 for all galaxies. The stellar surface density is assumed to be of the form: (65)with .

5.1. Local spiral galaxies

The sample of local spiral galaxies (Table B.1) is taken from Leroy et al. (2008). For this work, we did not consider dwarf galaxies, which will be the subject of a subsequent article. The spiral galaxies have total stellar masses in excess of 10¹⁰M_⊙. The gas masses are derived from IRAM 30 m CO(2–1) HERACLES (Leroy et al. 2009) and VLA Hi THINGS (Walter et al. 2008) data. The star-formation rate was derived from Spitzer MIR and GALEX UV data (Leroy et al. 2008). The total infrared luminosities are taken from Dale et al. (2012). Following Vollmer & Leroy (2011), the Toomre parameter of NGC 628, NGC 3198, NGC 5194, and NGC 7331 was set to Q(R) = Q + 3 exp(−2 R/l_∗), that of NGC 3351 by Q(R) = Q−4 exp(−(2 R/l_∗)²). The Toomre parameter Q was assumed to be constant for all other galaxies (see Table B.1).

5.2. Local ULIRGs

The ULIRG sample (Table B.2) is taken from Downes & Solomon (1998). These authors derived the spatial extent, rotation velocity, gas mass, and dynamical mass M_dyn for local ULIRGs from PdB interferometric CO-line observations. The total infrared luminosities are taken from Graciá-Carpio et al. (2008). We adopted the star-formation rates based on FIR data from Graciá-Carpio et al. (2008). We calculated the stellar mass as M_∗ = M_dyn−M_gas and assumed that the stellar scale length is approximately equal to the observed CO scale length.

5.3. High-z star-forming galaxies

The high-z star-forming sample (Table B.4) is taken from PHIBSS (Tacconi et al. 2013), the IRAM PdB high-z blue sequence CO(3–2) survey of the molecular gas properties in massive, main-sequence star-forming galaxies at z = 1–1.5. For our purpose, we only took the disk galaxies from PHIBSS. The stellar masses given by Tacconi et al. (2013) were derived from SED fitting, assuming a Chabrier IMF. Following Genzel et al. (2010), we calculated the star-formation rate from the total infrared luminosity using . Their star-formation rates are based on the sum of the observed UV- and IR-luminosities, or an extinction-corrected Hα luminosity. Their half-light radii were derived from Sersic fits to the HST ACS and/or WFC3 CANDELS data (Grogin et al. 2011). To estimate the characteristic circular velocities, Tacconi et al. (2013) took the isotropic virial estimate , where Δv_FWHM is the CO(3–2) linewidth for unresolved galaxies without a velocity gradient, and v_circ = 1.3 (Δv_blue−red/ (2sin(i))) if a velocity gradient (Δv_blue−red) indicative of rotation is detected in a galaxy with an inclination i. Since the inclination angle is difficult to determine in these high-z star-forming disk galaxies we adopted the following strategy for the determination of the rotation velocity: if the assumed rotation velocity is ; otherwise v_rot = v_circ. In this way, the rotation velocity of 15 out of 45 galaxies was increased by more than 50%.

5.4. Submillimeter galaxies

The smm-galaxy sample (Table B.3) was drawn from Genzel et al. (2010). The total infrared luminosities are based on the 850μm flux densities (Genzel et al. 2010). We calculated the star-formation rate from the total infrared luminosity using . Stellar masses are from the SED fits in Erb et al. (2006) and Förster Schreiber et al. (2009) and rotation velocities and half-light radii are from the data in the same references with the methods discussed in Förster Schreiber et al. (2009). When the stellar mass was not available, it was set to 10¹¹M_⊙. The CO-line observations were made in the CO(2–1), CO(3–2), and CO(4–3) lines.

6. Results

The aim of the present work is to directly compare infrared and molecular line luminosities. Before doing so, the distribution of the galaxy metallicities, dust SEDs, TIR luminosities, and dust temperatures are presented.

6.1. Metallicity

As described in Sect. 2, we use a simple close-box model for the metallicity (Eq. (28)). The oxygen abundance is then calculated by 12 + log (O/H) = log (Z/Z_⊙) + 8.7, with 12 + log (O/H) = 8.7 being the solar oxygen abundance (Asplund et al. 2005).

For the global mean metallicity, the average is weighted by optical emission from a disk where a uniform distribution of stars is mixed homogeneously with dust (“mixed” model of McLeod et al. 1993): (66)with τ = Σ/(7.5 M_⊙ pc^-2). The mean oxygen abundances resulting from our modeling ⟨ 12 + log (O/H) ⟩ are presented in Fig. 3.

The oxygen abundances of our local spiral galaxy sample range is 8.7 ≤ 12 + log (O/H) ≤ 9.0. This is close to the findings of Moustakas & Kennicutt (2006) who studied 14 nearby disk galaxies with integrated spectrophotometry and observations of more than 250 individual Hii regions; their oxygen abundances based on the McGaugh (1991) calibration also fall in the range 8.6 ≤ 12 + log (O/H) ≤ 9.0, with most of the galaxies having an oxygen abundance of 8.8–9.0.

The oxygen abundances of our ULIRG sample fall in the range 8.6 ≤ 12 + log (O/H) ≤ 9.4. The mean oxygen abundance of the sample is 12 + log (O/H) = 9.0 ± 0.3. Compared to the oxygen abundances determined by Rupke et al. (2008) and Kilerci Eser et al. (2014), with lower and upper limits of 8.4 and 9.0, respectively, four out of nine ULIRGs show model oxygen abundances in excess of 9.0. The model metallicities are thus up to a factor of 2.5 higher than the observed metallicities derived from optical emission line diagnostics.

The oxygen abundances or metallicities of the smm-galaxy sample are in the range 8.7 ≤ 12 + log (O/H) ≤ 9.4, that is, the metallicities are mostly supersolar. Three out of ten galaxies show metallicities in excess of 12 + log (O/H) = 9.4. We decided not to arbitrarily modify these obviously too high model metallicities. The mean oxygen abundance of the smm-galaxy sample is 12 + log (O/H) = 9.3 ± 0.4. Swinbank et al. (2004) found slightly subsolar metallicities in their sample of 30 smm-galaxies at a median redshift of z ~ 2.4. Tecza et al. (2004) found a supersolar oxygen abundance of 12 + log (O/H) = 9.0 for SMM J14011+0252 at z = 2.57. Nagao et al. (2012) found a solar metallicity of the submillimeter galaxy LESS J033229.4–275619 at z = 4.76. The model metallicities of our smm-galaxy sample are thus a factor 2–4 higher than the observed metallicities.

The oxygen abundances of our high-z star-forming galaxy sample range between 8.0 ≤ 12 + log (O/H) ≤ 9.0. The mean oxygen abundance of the sample is 12 + log (O/H) = 8.6 ± 0.2, thus very close to solar metallicity. This is consistent with the results of Shapley et al. (2004) who found solar and possibly supersolar metallicities in high-z star-forming galaxies. It is also consistent with the mean metallicity of 12 + log (O/H) = 8.7 ± 0.2 of the sample of 50 galaxies at z ~ 1.2 in the MASSIV survey (Queyrel et al. 2012).

We conclude that the integrated model metallicities of the local spiral and high-z star-forming galaxies are in good agreement with observations. On the other hand, the model metallicities of half of the galaxies of the ULIRG and smm-galaxy sample show metallicities that are 2–4 times higher than the observed metallicities. This leads to a potential overestimation of the molecular line emission (∝X^0.4; Scoville & Solomon 1974) by a factor of 1.3–1.7.

6.2. Dust SED, TIR luminosity, and dust temperature

For the direct comparison between the model and observed dust IR spectral energy distributions (SED), we extracted all available photometric data points for our galaxy samples from the CDS VizieR database³. Since the flux densities in the different catalogs are determined within different apertures, we only take the highest flux densities for a given wavelength range around a central wavelength λ₀ (0.75 ≤ λ/λ₀ ≤ 1.25). In this way, only the outer envelope of the flux density distribution is selected. We did not attempt to remove flux densities below this outer envelope, which most probably corresponds to apertures that do not include the whole object or are erroneous. Since our dust model does not include stochastically heated small grains and PAHs, the observed IR flux densities for λ ≲ 50 μm cannot be reproduced by the model. Flux densities at wavelengths ≥70μm are available for all local spirals/ULIRGs, 9 out of 10 smm-galaxies, and 30 out of 44 high-z star-forming galaxies.

The model dust IR SEDs of the local spiral galaxies reproduce the observed SEDs very well (Fig. C.1). Only the flux densities at λ> 200 μm NGC 4736, NGC 4535, NGC 6946, and NGC 3627 are somewhat overestimated by the model.

The models of the ULIRGs reproduce the existing observations very well (Fig. C.2). The comparison between the model and observed SED is difficult for Arp 220, because the observed SED contains the whole system, whereas the model SEDs are made separately for the Disk, Western, and Eastern nuclei. For the comparison, we assumed that 30%, 20%, and 30% of the observed total flux densities are emitted by the Disk, Western, and Eastern nuclei, respectively.

The models of the smm-galaxies reproduce the existing observations in a satisfactory way (Fig. C.3). Only for SMM J123549+6215 are the infrared flux densities underestimated by approximately a factor of 2.

As for the local spiral galaxy sample, the model dust IR SEDs of the high-z star-forming galaxies reproduce the observed SEDs well (Fig. C.4) at almost all wavelengths, especially those of the z ~ 1.5 (EGS) sample. The model SEDs of the z ~ 2.5 sample underestimate the observed SEDs by up to a factor of 2 for BzK 4171 and BzK 16000 and overestimate the observed SEDs by approximately a factor of 2 for BzK 17999.

The comparison of the model total IR luminosities (from 10 μm to 1000 μm) to the observed total IR luminosities is presented in Fig. 4. The model reproduces the total IR luminosities within a factor of 2. Only for two smm-galaxies and two high-z star-forming galaxies does the model underestimate the total IR luminosities by more than a factor of ~2.5.

We fitted modified Planck functions with β = 1.5 (see Sect. 2.8) to the model dust IR SEDs to derive dust temperatures. The modified Planck functions are shown as red dashed lines in Figs. C.1 to C.4. The resulting distribution of dust temperatures for the different galaxy samples are presented in Fig. 5.

The dust temperatures of the local spiral galaxies range between 19 and 24 K. This is in excellent agreement with Dale et al. (2012; Fig. 10), which is not surprising given the good fit to the observed IR dust SEDs.

The dust temperatures of the ULIRGs range between 39 and 72 K with a mean of 50 ± 11, in reasonable agreement with the results of Symeonidis et al. (2013) who found that the majority of (U)LIRGs at all redshifts have mean dust temperatures between 25 and 45 K using IRAS- and Herschel-selected samples, and Hwang et al. (2012) who found a dust temperature range between 35 and 43 K based on Herschel IR SEDs. Our dust temperatures are somewhat lower than the temperature distribution (61 ± 9) found by Klaas et al. (2001) for local ULIRGs.

The dust temperatures of the smm-galaxies range between 31 and 64 K. The smm-galaxy dust temperatures cover approximately the same range as the ULIRG dust temperatures, but their mean dust temperature is somewhat smaller 43 ± 10 K compared to 50 ± 11 K for the ULIRG sample. This is in good agreement with the dust temperatures (30–45 K) of smm-galaxies with IR luminosities <10¹³L_⊙ of Hwang et al. (2010; Fig. 3).

The dust temperatures of the high-z star-forming galaxies range between 27 and 48 K with a mean of 33 ± 4 K. This is in good agreement with (i) the mean temperature of the stacked z ~ 1 sample (32 ± 2) of Magdis et al. (2012; Table 2) and (ii) the mean temperature of 30 K found for z ~ 1 sources by Magnelli et al. (2014; Eq. (4)) based on Herschel observations.

We conclude that our model reproduces the dust IR SEDs of all galaxy samples. The IR flux densities of one smm-galaxy and two high-z star-forming galaxies at z ~ 2.5 underestimated by a factor of ~2. The derived dust temperatures of all galaxy samples are consistent with MIR and FIR observations.

6.3. Molecular fraction and HI mass

In the present model, the molecular fraction is calculated for each cloud of scale l_cl. It is determined by the photodissocation of H₂ molecules or the finite lifetime of a cloud (see Sect. 3). As a first step, we compare this molecular fraction to the molecular fraction defined by Vollmer & Leroy (2011) as f_mol = t_ff/t_mol/ (1 + t_ff/t_mol), where the free-fall and molecule formation times are those of self-gravitating clouds at a galactic radius R (Fig. 6).

The mean deviation between the two molecular fractions is approximately 30%, with a maximum deviation of approximately 50% at large galactic radii. At small galactic radii, the molecular fraction of the present model is ~30% smaller, whereas at large galactic radii, it is up to ~50% higher than that of Vollmer & Leroy (2011). We thus conclude that both prescriptions are comparable. This is quite surprising, because the Vollmer & Leroy (2011) prescription is only based on the properties of the self-gravitating clouds. We interpret this result as evidence for the dominant role of self-gravitating clouds for the formation of molecular hydrogen in local spiral galaxies.

Since the model yields the molecular fraction of the ISM, we can calculate the Hi surface density radial profile and the Hi mass. The comparison between the observed and model Hi masses is presented in Fig. 7. The model reproduces the observed Hi masses within a factor of 2. Especially the Hi masses of the galaxies with M_HI> 10¹⁰M_⊙ are underestimated by a factor of ~2. This is intrinsic to the model and cannot be compensated by a modification of the free model parameters (Q and δ).

The associated radial Hi surface density profiles together with the observed Hi profiles from Leroy et al. (2008) are presented in the upper and lower panels of Fig. 8. Whereas the model profiles are all monotonically declining, most of the observed profiles are constant in the range 1 ≤ R/l_∗ ≤ 4. In addition, most of the observed profiles show a depression in the central part of the galactic disks, whereas the models often show a maximum toward the galaxy center. The latter difference can be explained by (i) an underestimation of the model molecular fraction or (ii) by the ionization of atomic hydrogen in the inner part of the galactic disks, that is, the inclusion of a warm ionized medium into the model.

In order to investigate the effect of ionization by the UV radiation of massive stars in the galactic disk, we use the ionization-recombination equilibrium to calculate the number surface density of ionized gas following Maloney (1993): (67)where ϕ is the flux of ionizing photons and n the gas number density in cm^-3. The surface density of ionized gas is then Σ_ion = 1.36 × m_p × N_ion. Newborn massive stars are preferentially located in high-density regions, which they illuminate. As long as the area filling factor of the high-density gas surrounding these stars is not too large, the UV photons can escape the Hii regions and ionize the warm neutral medium. This scenario is supported by Hα observations of Thilker et al. (2002) and Oey et al. (2007) who found a fraction of diffuse Hα emission of 0.5–0.6 for local spiral galaxies. Based on Galactic Hi observations (Dickey & Lockman 1990), we used a constant gas number density of n = 5 cm^-3 for the warm neutral medium, which corresponds to the observed midplane density of 0.6 cm^-3 (Dickey & Lockman 1990) and a volume-filling factor of 0.12. For the relation between the ionizing photon flux and the star-formation rate we use (68)based on the star-formation rate calibration of Kennicutt (1998). We then calculated the Hi surface density as (69)The resulting Hi surface density is presented in the middle panel of Fig. 8. As the observed Hi surface density, the model Hi surface density is approximately constant Σ_HI ~ 10M_⊙ pc^-2 between R = l_∗ and R = 4 × l_∗. We conclude that the inclusion of the warm ionized medium in the model leads to radial Hi surface density profiles that are well comparable to observations. However, for the galaxies with the highest Hi masses (NGC 2841, NGC 3198, NGC 3521, and NGC 5055), the model Hi surface density profiles and total Hi masses are underestimated by approximately a factor of two.

6.4. CO(1–0) and HCN(1–0) radial profiles for local spiral galaxies

The CO and HCN emission can be spatially resolved in local spiral galaxies only (at D = 10 Mpc, the beam sizes are ~20′′ or ~1 kpc). Whereas CO maps are frequently found in the literature (e.g., Wong & Blitz 2002; Leroy et al. 2008), HCN maps are rare (e.g. Chen et al. 2015; Bigiel et al. 2016). The H₂ and dense gas surface densities are usually calculated with a CO-H₂ conversion factor of α_CO = 4.36M_⊙ pc^-2/(K km s^-1) (e.g., Bolatto et al. 2013) and α_HCN = 10M_⊙ pc^-2/(K km s^-1) (e.g., Gao & Solomon 2004). The model SFR-Σ_H2 relation of the local spiral galaxies is presented in the top panel of Fig. 9 together with the observed relation (e.g., Bigiel et al. 2008; Leroy et al. 2008), which corresponds to a constant star-formation rate timescale of 2 × 10⁹ yr. The model SFR-Σ_H2(CO) relations are well consistent with the observed relation within a scatter of approximately 0.2 dex. Overall, they show a somewhat flatter slope than the observed relation.

The situation is more complex for the dense gas (HCN) than for H₂. The observed –Σ_dense relations show different slopes; a linear slope was found by Gao & Solomon (2004) and Garcia-Burillo et al. (2012) for integrated SFR and gas masses, whereas a sub-linear slope of 0.7 was found by Usero et al. (2015) using single pointings in nearby spiral galaxies. For a direct comparison between the model and observed –Σ_dense relations, the conversion factor between the total infrared luminosity and the star-formation rate c_TIR must be taken into account, SFR(M_⊙ yr^-1) = c_TIRL_TIR(L_⊙). The different conversion factors are c_TIR = 2.0, 1.7, 1.5 × 10^-10 used by Gao & Solomon (2004), Garcia-Burillo et al. (2012), and Usero et al. (2015), respectively. The model conversion factor for the local spiral galaxy sample is c_TIR = 0.9 × 10^-10, that is, a factor of ~2 lower than the values used in the literature. For a consistent comparison between the model and observations, we multiplied the star-formation rates of the observed relation by . The resulting model –Σ_dense relation (lower panel of Fig. 9) has a slope consistent with observations, but shows a negative offset of 0.2–0.3 dex with respect to the observed relations.

We conclude that the resolved –Σ_H2 and –Σ_dense relations are consistent with available IR, CO, and HCN observations within the uncertainties (see Sect. 4).

6.5. Integrated CO, HCN(1–0), and HCO⁺(1–0) flux densities

For the ULIRG, smm, and high-z star-forming galaxies, we can only compare the integrated model HCN, CO, and HCO⁺ luminosities to observations. The comparison between the model and the observed CO luminosities is shown in Fig. 10, where the observed transitions were used: CO(2–1) for the local spiral galaxies, CO(1–0) for the ULIRGs, CO(3–2) for the smm-galaxies, and CO(3–2) for the high-z star-forming galaxies. The corresponding mean ratio L_model/L_obs and its uncertainty are presented in Table 2 (preferred model). We observe approximately linear correlations between the model and observed CO luminosities. The model CO luminosities of the smm and high-z star-forming galaxies are ⟨ log (L_CO,obs/L_CO,model) ⟩ ~ 0.2 dex smaller than observed; the ratio is only ~0.1 dex for the local spirals and ULIRGs. The corresponding model and observed L_TIR– relations for the CO(1–0) line are shown in Fig. 11. Line ratios from the literature are applied to determine the CO(1–0) fluxes (CO(2–1)/CO(1–0) = 0.7, CO(3–2)/CO(1–0) = 0.77, Genzel et al. 2010; and CO(3–2)/CO(1–0) = 0.5, Tacconi et al. 2013). The differences between the model and observed relation are again due to the underestimation (0.2 dex) of the model CO luminosities of the smm and high-z star-forming galaxies with respect to observations. Overall, the model L_TIR– relations are consistent with the observed relations.

For the HCN(1–0) luminosities, the situation is more complex. We only compare L_TIR– relations (upper panel of Fig. 12) for the local spirals, smm, and high-z star-forming galaxies, because the overlap between the model and observed samples is very small. For the ULIRG sample, there are six out of nine galaxies in common between the Downes & Solomon (1998) and Graciá-Carpio et al. (2008) samples, which permits a direct comparison of the HCN(1–0) luminosities (lower panel of Fig. 12). Since the observed relations have different slopes, we compare the model relation to the relations observed by Gao & Solomon (2004) and Graciá-Carpio et al. (2008). Based on the L_TIR– relation, the model overestimates the HCN(1–0) luminosities of the local spiral galaxies by ~0.2 dex and underestimates those of ULIRGs by ~0.3 dex. However, the direct comparison of model and observed HCN(1–0) luminosities for ULIRGs yields an underestimation of only 0.13 dex. This means that the Downes & Solomon (1998) ULIRG sample contains mainly HCN-bright galaxies. Concerning the HCN emission of smm-galaxies, Gao & Solomon (2007) claimed that the FIR/HCN ratios in these high-redshift sources lie systematically above the FIR/HCN correlation established for nearby galaxies by approximately a factor of 2. This behavior is well reproduced by the model. Since there are no HCN detections of high-z star-forming galaxies in the literature, we can only suggest that their HCN emission might be a factor of 3 lower than expected from observed L_TIR– relations.

The model HCN/CO ratio is compared to observations in Fig. 13 for the four galaxy samples. As expected, the model points for the local spiral galaxies lie below whereas those of the ULIRGs lie above the observed correlation. Within the model, the smm-galaxies follow the same correlation as the ULIRGs, whereas the HCN emission of the high-z star-forming galaxies is significantly lower than expected by the observed correlation.

The model integrated HCO⁺ luminosity can be compared to observations via the L_TIR– relation (upper panel of Fig. 14). The observed relations (Juneau et al. 2009; Garcia-Burillo et al. 2012) do not differ significantly. Within the model, the ULIRGs and high-z star-forming galaxies follow the observed relations. The smm-galaxies lie somewhat below the observed relation. The model HCO⁺ luminosities are a factor of approximately three higher than the expected HCO⁺ luminosities assuming (e.g., Nguyen et al. 1992; Brouillet et al. 2005; Knudsen et al. 2007). The direct comparison of the model and observed HCO⁺ luminosities for the ULIRG sample shows good agreement (lower panel of Fig. 14). The model thus seems to overestimate the HCO⁺ luminosities of the local spiral galaxies. We note that the HCO⁺ emission strongly depends on the cosmic ray ionization rate used for the chemical network.

As an additional consistency check, the HCN/HCO⁺ ratio as a function of the total infrared luminosity is shown in Fig. 15. Observations of local spiral galaxies (e.g., Nguyen et al. 1992; Brouillet et al. 2005; Knudsen et al. 2007) show . As expected, the HCN/HCO⁺ ratio of the local spiral galaxies is smaller by ~0.2 dex than the observed ratio. The model HCN/HCO⁺ ratio for the ULIRGs are well within the observed range (Graciá-Carpio et al. 2008).

We conclude that, overall, the CO luminosities are underestimated by up to a factor of 1.5. The model HCN luminosities for local spirals are overestimated by a factor of 1.5 and those for ULIRGs are underestimated by a factor of 1.5 with respect to observations. The model HCN luminosities are consistent with observations for the smm-galaxies. The model HCO⁺ luminosities are consistent observations for all galaxy samples except the local spiral galaxies, where they are significantly overestimated.

6.6. CO SLEDs

The lowest three rotational transitions of CO, which trace the cooler gas component, are relatively easily accessible with ground-based radio and submillimeter telescopes, and have been observed in many local galaxies. It was not until the launch of the Herschel Space Observatory (Pilbratt et al. 2010) that the CO ladder up to J = 13 was generally available for the ISM within our Galaxy and in nearby galaxies. Early SPIRE observations showed much brighter high-J CO emission than would be predicted by cool (T_kin< 50 K) molecular gas in giant molecular clouds, the type of gas responsible for the CO(1–0) emission. A warmer, denser (higher pressure) component of molecular gas is responsible for the emission of mid- (J = 4–3 to J = 6–5) and high-J (J = 7–6 and above) CO lines. Far-IR CO rotational lines, with J_upper ≥ 13, arise from states 500–7000 K above the ground state and have critical densities of 10⁶ to 10⁸ cm^-3 (Hailey-Dunsheath et al. 2012).

It is well known locally that empirically and observationally, different CO excitation properties characterize spiral galaxies as opposed to merging-driven ULIRGs (e.g., Daddi et al. 2015). The latter are much more highly excited in their high-J CO transitions (Weiss et al. 2007; Papadopoulos et al. 2012). The model CO SLEDs of all galaxies of the four samples and the mean CO SLEDs are presented in Fig. 16. The latter can be directly compared to the observed mean SLEDs in Fig. 8 of Daddi et al. (2015). We can only compare the model CO SLEDs of the local spiral galaxies to that of the inner Milky Way, which shows, as expected, higher intensities for transition with J_upper ≥ 3. The model CO SLEDs of the ULIRGs are consistent in shape and absolute values with observations (Daddi et al. 2015). The shape of the model CO SLEDs of the smm and high-z star-forming galaxies is different from that of the observed SLEDs; whereas the shape of the model CO SLEDs is concave, that of the observed CO SLED is convex. This difference is mainly due to the CO(3–2) flux, which in the model is a factor of two higher than observed.

In a sample of ULIRGs, higher CO transitions (J_upper> 6) were observed by Kamenetzky et al. (2016). Their mean CO SLEDs for galaxies with total infrared luminosities between 3 × 10¹¹ and 10¹²L_⊙ and higher than 10¹²L_⊙ are shown together with the model ULIRG CO SLEDs in Fig. 17. The two compact residual disks of Arp 220 (E and W) cannot be directly compared to this sample because Arp 220 would be seen as one entity (Disk + East + West) by the Herschel satellite. The shape and the absolute values of the observed CO SLEDs are well reproduced by the model.

6.7. Integrated CO and HCN conversion factors

Since the integrated molecular fraction and the line emission are calculated within the model (see Sect. 6.3), the integrated mass-to-light conversion factors can be determined. For the CO line emission we use two approaches: (i) the model CO(1–0) is adopted for all galaxies (upper panel of Fig. 18) and (ii) the model flux of the observed CO line (CO(2–1) for the local spiral galaxies, CO(1–0) for the ULIRGs, CO(3–2) for the smm-galaxies, and CO(3–2) for the high-z star-forming galaxies) is calculated and a line ratio given in the literature is applied to determine the CO(1–0) flux (CO(2–1)/CO(1–0)= 0.7, CO(3–2)/CO(1–0) = 0.77, Genzel et al. 2010; and CO(3–2)/CO(1–0)= 0.5, Tacconi et al. 2013; lower panel of Fig. 18). The mass-to-light conversion factor is then calculated by . The mean conversion factors are given in Table 3.

The model CO(1–0) conversion factor varies between 2M_⊙(K km s^-1 pc²)^-1 and 6M_⊙(K km s^-1 pc²)^-1 with one galaxy showing α_CO = 10M_⊙(K km s^-1 pc²)^-1. The mean conversion factor of local spirals (⟨ α_CO ⟩ = 4.7 ± 1.8M_⊙(K km s^-1 pc²)^-1) is close to the observed value of α_CO = 4.3M_⊙(K km s^-1 pc²)^-1 (Bolatto et al. 2013). The mean conversion factor of the ULIRGs (⟨ α_CO ⟩ = 1.7 ± 0.4M_⊙(K km s^-1 pc²)^-1) is twice the usually assumed conversion factor for dense starburst galaxies (α_CO = 0.8M_⊙(K km s^-1 pc²)^-1; Downes & Solomon 1998) while that of smm-galaxies is similar to it. The model CO(1–0) conversion factor of high-z star-forming galaxies is intermediate between those of local spiral galaxies and ULIRGs/smm-galaxies (⟨ α_CO ⟩ = 2.6 ± 0.9M_⊙(K km s^-1 pc²)^-1).

The situation changes only slightly if higher J transitions are used instead of CO(1–0). The only remarkable change is the decrease of the CO conversion factor of high-z star-forming galaxies by 30% to (⟨ α_CO ⟩ = 1.6 ± 0.6M_⊙(K km s^-1 pc²)^-1). However, the distribution of the CO conversion factors for high-z star-forming galaxies shows a tail with conversion factors comparable to that of the Galaxy.

The HCN conversion factor is usually given with respect to the dense gas, that is, gas with densities exceeding n = 3 × 10⁴ cm^-3 (e.g., Gao & Solomon 2004). For completeness, we give the HCN(1–0)–M_H2 and the HCN(1–0)–M_dense conversion factors in Fig. 19. Only the latter conversion factor can be compared to the literature (lower panel of Fig. 19), where α_HCN = 10M_⊙(K km s^-1 pc²)^-1 is assumed (e.g., Gao & Solomon 2004). The mean HCN conversion factors are α_HCN = 21 ± 6, 33 ± 17, and 59 ± 21M_⊙(K km s^-1 pc²)^-1 for the local spiral galaxies/ULIRGs, smm-galaxies, and high-z star-forming galaxies, respectively (Table 3).

The mean HCO⁺ conversion factors are α_{HCO +} = 11 ± 2, 17 ± 5, 19 ± 11, and 25 ± 7M_⊙(K km s^-1 pc²)^-1 for the local spiral galaxies, ULIRGs, smm-galaxies, and high-z star-forming galaxies, respectively (Table 3). Since the HCO⁺(1–0) luminosity of the local spiral galaxies is significantly overestimated by the model, the associated conversion factor is a lower limit. We thus find a relatively uniform HCO⁺ conversion factor α_{HCO +} ~ 20M_⊙(K km s^-1 pc²)^-1 for all galaxy samples.

We conclude that all model mass-to-light conversion factors are consistent with the values used in the literature within a factor of two. Both, the HCN and HCO⁺ emission trace the dense molecular gas to a factor of approximately two for the local spiral galaxies, ULIRGs and smm-galaxies. For the high-z star-forming galaxies, HCO⁺ might be the better tracer, but this needs to be confirmed.

7. Variation of model parameters

Only the preferred model is presented in Sect. 6. A natural question to ask is how our model results depend on the different model assumption and free parameters. To answer this question, we replace the chemical network by constant molecular abundances and vary the free parameters Q and δ. In addition, in Sect. 8, we remove the cloud substructure, cosmic ray heating, and HCN IR-pumping.

7.1. Importance of the chemical network

Initially, we adopted constant molecular abundances for the CO, HCN, and HCO⁺ molecules, which are then scaled with metallicity: x_CO = 10^-4, x_HCN = x_HCO⁺ = 2 × 10^-8. Whereas the CO abundance corresponds to the canonical value (e.g., Draine 2011), the assumption of a constant HCN and HCO⁺ is not well justified. The resulting ratios between the model and observed line luminosities are shown in Table 2.

The mean of the ratios between the model and observed CO luminosities are for the local spiral, −0.10 ± 0.16 for the ULIRG, −0.19 ± 0.12 for the smm-galaxy, and −0.20 ± 0.22 for the high-z star-forming galaxy samples, respectively. Thus, the observed CO luminosities are reproduced for all galaxy samples by the model within ~0.2 dex or a factor of ~1.6.

The resulting HCN(1–0) luminosities are compared to the observed CO luminosities in Fig. 20. In order to compare all model HCN luminosities to observations even in the absence of HCN measurements, we assumed (Gao & Solomon 2004). The mean of the ratios between the model and observed HCN luminosities are for the local spiral, −0.47 ± 0.15 for the ULIRG, −0.44 ± 0.07 for the smm-galaxy, and −0.75 ± 0.08 for the high-z star-forming galaxy samples, respectively. Whereas the model HCN luminosities agree well with the observed HCN luminosities for the local spiral galaxies, the model underestimates the HCN luminosities by up to a factor of four for the ULIRG, smm-galaxy, and high-z star-forming galaxies.

Observations of HCN(1–0)/HCO⁺(1–0) ratios of large galaxy samples are rare. Graciá-Carpio et al. (2006) and Juneau et al. (2009) found in local ULIRGs. The model with constant abundances systematically underestimates the HCN/HCO⁺ ratio by approximately a factor of three for all galaxy samples (Fig. 21). We conclude that, as expected, our model with constant HCN(1–0) and HCO⁺(1–0) abundances does not lead to HCN and HCO⁺ luminosities that are comparable to observations.

We conclude that models with the canonical CO abundances reproduce observations as well as models using a detailed chemical network. However, constant HCN and HCO⁺ abundances yield HCN and HCO⁺ line luminosities that are at least a factor of two smaller than the observed line luminosities.

7.2. Toomre Q

To investigate the influence of the Toomre Q parameter on the model CO and HCN line emission, we tested Q = 1 instead of Q = 1.5 for the ULIRGs, smm, and high-z star-forming galaxies (see Table 2). The comparison of the model CO and HCN luminosities and CO SLED with observations is presented in Fig. 22. The resulting ratios between the model and observed line luminosities are shown in Table 2. Whereas the model CO and HCO⁺ line luminosities are barely affected, the model HCN line luminosity decreases by ~15% when Q = 1 instead of Q = 1.5. The CO SLEDs of the smm and high-z star-forming galaxies do not change significantly with respect to Q = 1.5, but the ULIRG CO SLED increases by ~30% (lower panel of Fig. 22).

We conclude that the low J CO, HCN, and HCO⁺ line emissions of all galaxies are not significantly affected when Q is decreased from 1.5 to 1. However, a variation of Q changes the gas velocity dispersion significantly (Sect. 9.2).

7.3. The scale parameter δ

As described in Sect. 2.7, the scale of the largest self-gravitating clouds l_cl is smaller than the turbulent driving length scale l_driv by a factor δ = l_driv/l_cl. For a typical cloud size of 20 pc and a driving length scale of 100 pc, one obtains δ = 5, which is the mean value determined for local spiral galaxies by Vollmer & Leroy (2011). To investigate the influence of the scale parameter δ on the model CO and HCN line emission, we used δ = 15 instead of δ = 5 for all galaxies. The comparison of the model CO and HCN luminosities and CO SLED with observations is presented in Fig. 23. The resulting ratios between the model and observed line luminosities are shown in Table 2.

As for the Q = 1, the model CO, HCN, and HCO⁺ line luminosities are barely affected. The effect on the CO SLEDs is also comparable to the Q = 1 models: whereas the CO SLEDs of the local spiral galaxies, smm, and high-z star-forming galaxies do not change significantly with respect to Q = 1.5, the ULIRG CO SLED increases by ~30%.

We conclude that the low J CO, HCN, and HCO⁺ line emissions of all galaxies are not significantly affected, but the CO SLED of the ULIRGs is increased by ~30% when δ is increased by a factor of three. In contrast to the decrease of Q, the increase of δ does not result in a significant increase of the gas velocity dispersion (Sect. 9.2).

8. Galactic physics

In this section we examine the role of non-self-gravitating clouds for molecular line emission and investigate how different recipes for galactic physics influence the model results.

8.1. The role of non-self-gravitating clouds for molecular line emission

Within the framework of the analytical model (Sect. 2), the scale parameter δ determines the density of the self-gravitating clouds ρ_cl = ρ/φ_V via . At this characteristic spatial length scale, the scaling relations for the density and velocity dispersion change (see Sect. 2.1.2). We would like to know which fractions of molecular line emission originate in self-gravitating and non-self-gravitating gas clouds in our galaxy samples. To answer this question, the ratio between the molecular line emission of non-self-gravitating clouds and the total line emission is shown as a function of the total infrared luminosity in Fig. 24.

In the local spiral galaxies, approximately half of the total CO emission originates in non-self-gravitating clouds (f_nsg,CO ~ 0.5). This is consistent with the ratio between CO emission from diffuse clouds and the total emission I_diff/I_tot = 0.6 ± 0.1 determined by Polk et al. (1988)⁴.

In the ULIRGs, the fraction is f_nsg,CO ~ 0.8. The smm-galaxies also show high fractions of emission from non-self-gravitating clouds, however there are also two smm-galaxies with f_nsg,CO ~ 0.4. The high-z star-forming galaxies show an approximately flat distribution with 0.3 ≤ f_nsg,CO ≤ 0.9.

For the HCN and HCO⁺ emission, the situation is different. Since the critical densities for these transitions are much higher than for the CO emission, it is expected that the denser self-gravitating clouds dominate the HCN and HCO⁺ emission. This is indeed the case for the local spiral galaxies and the high-z star-forming galaxies. However, in the ULIRGs and smm galaxies, the HCN and HCO⁺ emission mostly originates in non-self-gravitating clouds (0.7 ≤ f_{nsg,HCN/HCO⁺} ≤ 0.9). Even in local spiral galaxies, the model predicts that ~30% of the HCN(1–0) emission is emitted by non-self-gravitating clouds. In these clouds, the effective critical density might be as low as cm^-3 due to radiative trapping, approximately 30 times lower than the critical density in the optically thin limit (Shirley 2015). Future combined interferometric and single dish HCN(1–0) observations of local spiral galaxies will be able to test our prediction.

We conclude that compact starburst galaxies (ULIRGs and smm-galaxies) have, on average, significantly higher fractions of molecular line emission that originates in non-self-gravitating clouds. This effect is more pronounced for HCN(1–0) and HCO⁺(1–0) emission. In compact starburst galaxies, the molecular line emission is most frequently dominated by emission from non-self-gravitating clouds/filaments. This is due to the highly compact nature of ULIRG and smm-galaxie centers, such that the density of the intercloud medium rivals or even exceeds the density of Galactic giant molecular clouds.

8.2. The role of cloud substructure

The molecular line emission of giant molecular clouds shows substructure (e.g., Heyer & Dame 2015). This substructure is taken into account by the model (Sect. 2). To investigate the role of GMC substructure, we tested a model where self-gravitating clouds with uniform densities are assumed. The resulting CO and HCN luminosities are presented in Fig. 25.

It might be expected that the higher densities in cloud substructures, that is, cores, increase the molecular line emission. To our surprise, the changes, with respect to the model including substructure of self-gravitating clouds, are marginal (Table 2). As expected, these small changes are most visible at high-J CO transitions.

This effect can be explained by a balance between the absence of high-density gas and the increase of the area-filling factor of emission from self-gravitating clouds with a uniform gas density.

8.3. Cosmic ray heating

The molecular line emission depends on the gas density, column density, and temperature. Within the model framework, the gas temperature within the dense gas clouds depends on the turbulent and cosmic ray heating (see Sect. 2.1.3). The cosmic ray ionization rate is a factor of 40 higher in ULIRGs and smm-galaxies than in local spiral galaxies. To investigate the influence of the cosmic ray heating, we made a model calculation without this heating term. The resulting CO and HCN luminosities are presented in Fig. 26. The CO(1–0) luminosities of all galaxies are not altered by the absence or presence of cosmic ray heating (Table 2). This implies that, at the relevant density regime of 300–1000 cm^-3, turbulent heating dominates over cosmic ray heating.

One would expect that the situation changes for molecular line transitions with higher critical densities because the cosmic ray heating is proportional to the gas density, whereas the turbulent heating is proportional to the square root of the gas density (see Sect. 2.1.3). Surprisingly, the model ULIRG CO emission of transitions with upper J of 2 ≥ J ≥ 11 only changes by at most ±0.1 dex. Four ULIRGs even show a slight increase of CO emission around J = 7 when the CR ray heating is suppressed. We suspect that the effect is due to a changing chemistry that requires thorough investigation, something beyond the scope of this work.

There is no effect of CR heating on HCN(1–0) and HCO⁺(1–0) emission for the local spiral and smm-galaxies. In the absence of CR heating, the HCN(1–0) line luminosity decreases by a factor of approximately 1.5 for ULIRGs and high-z star-forming galaxies. The HCO⁺ (1–0) line emission shows the same behavior as the HCN(1–0) line emission for local spiral galaxies, smm-galaxies, and high-z star-forming galaxies. For ULIRGs, the situation is more complicated; whereas the overall comparison with the observed L_IR–L_HCN relation shows a decreasing line emission in the absence of CR heating, the direct comparison of the model and observed HCN luminosities is consistent with no change.

We conclude that within the framework of this model, the inclusion of CR heating does not significantly affect the CO emission, but increases the HCN(1–0), and HCO⁺(1–0) emission by at most a factor of two. This factor is needed to reproduce the observed HCN(1–0) emission of the ULIRGs.

8.4. Infrared-pumping

Infrared pumping of HCN (Sect. 2.10) via the 14 μm bending modes is suggested to play an important role for the HCN(1–0) emission of ULIRGs (e.g., Aalto et al. 2015). To investigate the influence of infrared pumping on the HCN(1–0) line emission, we calculated a model without infrared pumping. The results are presented in Fig. 27 and Table 2. The HCN(1–0) emission of local spiral galaxies and smm-galaxies is not significantly affected by HCN infrared pumping. The HCN(1–0) emission of ULIRGs and high-z star-forming galaxies decreases by 20–30% when the infrared pumping is suppressed. Thus, within the framework of our model, HCN infrared pumping has a measurable but relatively small effect on the HCN(1–0) of ULIRGs and high-star-forming galaxies.

9. Physical properties of the galaxies

9.1. Gas fraction

The observed high gas fractions (M_H2/ (M_H2 + M_∗)) of high-z star-forming galaxies (Tacconi et al. 2013) strongly depend on the assumed CO conversion factor α_CO. Tacconi et al. (2013) used a Galactic conversion factor. Our CO conversion factor is a factor of approximately two smaller (Sect. 6.7). Genzel et al. (2015) claimed that the H₂ mass estimates of the high-z star-forming galaxies with a Galactic conversion factor agree to better than 50% with the H₂ mass estimates based on the infrared SEDs. We suppose that our molecular mass estimate is consistent with that of Tacconi et al. (2013) given the uncertainties of the observed and model (see Sect. 4) CO and dust conversion factors of approximately a factor of two. The infrared luminosities of high-z star-forming galaxies (~10¹¹L_⊙) are between those of local spiral galaxies (~10¹⁰L_⊙) and ULIRGs (~10¹²L_⊙). It seems thus reasonable to find conversion factors for high-z star-forming galaxies that range between the Galactic conversion factor and that for ULIRGs.

The model H₂ gas fractions are presented in Fig. 28.

Whereas the H₂ gas fraction of local spirals is approximately 5%, that of the high-z star-forming galaxies range between 10 and 30% with a mean of ~20%. Tacconi et al. (2013) found a H₂ gas fraction of 0.33. If we use a Galactic conversion factor, we recover these high gas fractions. However, we determined a conversion factor for the high-z star-forming galaxies which is half of the Galactic conversion factor. So, at least in our model counterparts of the high-z star-forming galaxies, we know that the H₂ gas fraction is lower.

In our model, the highest H₂ gas fractions, exceeding 30%, are found in ULIRGs. The widely varying H₂ gas fractions of the smm-galaxies greatly depend on the determination of the stellar masses from observations, which have large uncertainties.

9.2. Gas velocity dispersion

An important finding in many high-redshift star-forming disc galaxies is the high gas velocity dispersion. Values of 50–100 km s^-1 (Genzel et al. 2006; Law et al. 2009; Förster Schreiber et al. 2009; Vergani et al. 2012; Tacconi et al. 2013) are frequently observed. Whether or not these high-velocity dispersions are mandatory for the high-z star-forming galaxies remains to be elucidated. To find an initial answer to this question, we compare our preferred model, Q ~ 1.5 (Fig. 29), to the Q ~ 1 model (Fig. 30). The upper panels show the comparison between the observed (Hα and CO) and modelled gas velocity dispersions. The galaxies within the different samples (local spiral galaxies, ULIRGs, smm, high-z star-forming galaxies) are not the same. Clearly, the decrease of Q at fixed star-formation rate from Q ~ 1.5 to Q ~ 1 has a big effect on the gas velocity dispersions, which decreases by a factor of 2 for the ULIRGs and smm-galaxies and by a factor of 1.6 for the high-z star-forming galaxies. A lower Q decreases the velocity dispersion, but also increases the star-formation rate at fixed gas surface density. Thus, at fixed star-formation rate, a lower Q decreases both the gas surface density and the velocity dispersion. Our preferred model (Q ~ 1.5) better reproduces the observed velocity dispersions.

A direct comparison of observed and modeled gas velocity dispersions for single galaxies in presented in the lower panels of Figs. 29 and 30. Within the preferred Q ~ 1.5 model the gas velocity dispersions of two ULIRGs and three high-z star-forming galaxies exceed the observed values by approximately a factor of two. On the other hand, the gas velocity dispersions of four ULIRGs is approximately two times lower than the observed values. Given that the measured velocity dispersion of ULIRGs and high-z star-forming galaxies that are barely spatially resolved can be easily dominated by non-circular gas motions, one expects the model gas velocity dispersions to be systematically smaller than the observed velocity dispersions. This is only the case for the Q ~ 1 model.

We conclude that the preferred Q ~ 1.5 model better reproduces the high gas velocity dispersions observed in ULIRGs and high-star-forming galaxies. However, if the observed linewidths are dominated by non-circular gas motions, the Q ~ 1 model is consistent with available observations. We recall that the HCN(1–0) emission of the Q ~ 1 model is somewhat smaller than that of the preferred model leading to a less good reproduction of the available HCN(1–0) observations (see Sect. 7.2).

9.3. Star-formation laws

Genzel et al. (2010) and Daddi et al. (2010) found a long-lasting star-formation mode for disk galaxies and a more rapid mode for starbursts. The two modes can be unified to a single star-formation law, if the star-formation timescale with respect to the molecular gas (t_SF = M_H2/SFR) observed in CO emission is assumed to be proportional to the dynamical timescale, that is, the angular velocity of the galaxy: . Our model directly yields the local star-formation rate, the H₂ gas surface density, and the angular velocity.

The Kennicutt-Schmidt relation for the integrated star-formation rates and H₂ gas masses of our model is presented in Fig. 31.

The slopes of the relation for the local spiral galaxies and ULIRGs/smm galaxies are approximately unity (1.0 and 1.2). We observe a factor of approximately 50 between the star-formation efficiencies SFE_H2 of the local spiral galaxies on the one hand, and the ULIRGs/smm-galaxies on the other. The SFR–M_H2 relation of the high-z star-forming galaxies is intermediate between those of the local spiral galaxies and ULIRGs/smm-galaxies with a slope of 1.7.

The Kennicutt-Schmidt relation with respect to the molecular gas surface density of our model galaxies is presented in the upper panel of Fig. 32. For the area calculation, the stellar scalelength was adopted. The slope of the relation for the local spiral galaxies derived by the IDL routine robust_linefit is 1.5. The fitted slopes of the relations for ULIRGs and high-z star-forming galaxies are 1.5 and 1.4, respectively. The slope of the relation for smm-galaxies is 2, that of the combined ULIRG and smm-galaxy sample is 1.6. We observe an offset of a factor of approximately 7 between the relation of the high-z star-forming and the local spiral galaxies.

The model star formation per unit area as a function of the model molecular gas surface density multiplied by the angular velocity (or divided by the dynamical timescale) Σ_H2Ω is presented in the lower panel of Fig. 32. For the calculation of the angular velocity Ω = v_rot/R, we followed Daddi et al. (2010) and took optical radius R = R₂₅ = 4.5 × l_∗ for the local spiral galaxies (see also Kennicutt 1998) and the half-light radius for the ULIRGs, smm-galaxies, and high-z star-forming galaxies. The slope of this relation is 1.0 for the entire sample. Compared to that of the local spiral galaxies, the molecular star-formation efficiency of the ULIRGs, smm, and high-z star-forming galaxies is approximately twice as high. With for the local spiral galaxies, this ratio would increase to a factor of 5.

We conclude that the model Kennicutt-Schmidt laws for the integrated H₂ masses and surface densities do not show the same slopes. The integrated Kennicutt-Schmidt law has a slope ofapproximately 1 for the local spirals, ULIRGs, and smm-galaxies, whereas the slope is 1.7 for high-z star-forming galaxies. The model shows Kennicutt-Schmidt laws with respect to the molecular gas surface density with slopes of approximately 1.5 for local spiral galaxies, ULIRGs, and high-z star-forming galaxies. The slope for the smm-galaxies is approximately 2. The model star-formation rate per unit area is, as observed, approximately proportional to the molecular gas surface density divided by the dynamical timescale.

10. Conclusions

The theory of clumpy gas disks (Vollmer & Beckert 2003) provides the large-scale and small-scale properties of galactic gas disks. Large-scale properties considered are the gas surface density, density, disk height, turbulent driving length scale, velocity dispersion, gas viscosity, volume filling factor, and molecular fraction. Small-scale properties are the mass, size, density, turbulent, free-fall, and molecular formation timescales of the most massive self-gravitating gas clouds. These quantities depend on the stellar surface density, the angular velocity Ω, the disk radius R, and three additional parameters, which are the Toomre parameter Q of the gas, the mass accretion rate Ṁ, and the ratio δ between the driving length scale of turbulence and the cloud size. The large-scale part of the model disk is governed by vertical pressure equilibrium, a constant Toomre Q parameter, conservation of the turbulent energy flux, a relation between the gas viscosity and the gas surface density, a star-formation recipe (Sect. 3), and a simple closed-box model for the gas metallicity. The small-scale part is divided into two parts according to gas density: non-self-gravitating and self-gravitating gas clouds. The mass fraction at a given density is determined by a density probability distribution involving the overdensity and the Mach number (Sect. 2.1.1). Both density regimes are governed by different observed scale relations (Sect. 2.1.2). The dense gas clouds are mechanically heated by turbulence. In addition, they are heated by cosmic rays. The gas temperature is calculated through the equilibrium between gas heating and cooling via molecular line emission (CO, H₂, H₂O; Sect. 2.1.3). The dust temperature is determined by the equilibrium between radiative heating and cooling and the heat transfer between gas and dust (Sect. 2.8). The molecular line emission calculation is based on the escape probability formalism (Sect. 2.9). An important ingredient for the line emission is the area-filling factor of the gas clouds, which is a result of the small-scale part of the analytic disk model. The molecular abundances of individual gas clouds are determined by a detailed chemical network involving the cloud lifetime, density, and temperature (Sect. 2.1.4). H₂ and CO dissociation in photodissociation regions are taken into account (Sect. 2.4.2). Moreover, a simple formalism for HCN infrared pumping is applied to the HCN line emission (Sect. 2.10).

The stellar radiation field is constrained by the observed IR luminosity and SED. The normalization of the cosmic ray ionization rate is constrained by the observed HCO⁺(1–0) emission. The density and temperature structure of the clumpy gas disk is constrained by the observed HCN(1–0) and multi-transition CO emission. This model is applied to samples of local spiral galaxies, ULIRGs, smm, and high-z star-forming galaxies (Sect. 5).

Based on the comparison between the model results and observations available in the literature (see Table 2) we conclude that

• 1.

  The following observed quantities are consistent withobservations:

  • global metallicities (Fig. 3);

  • total infrared luminosities and dust SEDs (Fig. 4, Appendix C);

  • dust temperatures (Fig. 5);

  • Hi masses and radial profiles of the local spiral galaxies (Figs. 7, 8);

  • CO luminosities (Fig. 10);

  • HCO⁺ luminosities of the ULIRGs, smm, and high-z star-forming galaxies (Fig. 14);

  • CO SLEDs up to J = 6 (Fig. 16);

  • CO SLEDs up to J = 12 for the ULIRGs (Fig. 17).

• 2.

  The model HCN radial profiles are a factor 1.5–2 higher than the observed profiles (Fig. 9).

• 3.

  The model HCN luminosities are a factor of 1.5 higher/lower than the observed luminosities for the local spiral galaxies/ULIRGs (Fig. 10).

• 4.

  The model HCO⁺ luminosities of the local spiral galaxies are a factor of ~3 higher than the observed HCO⁺ luminosities (Fig. 14); the HCO⁺ emission mainly depends on the CR ionzation rate used in the chemical network.

• 5.

  All model conversion factors (mass-to-light and SFR-to-light) have uncertainties of a factor of two.

• 6.

  The model CO conversion factors deduced when including CO-dark H₂ are α_CO = 4.7 ± 1.8, 1.7 ± 0.4, 1.4 ± 0.7, 2.6 ± 0.9M_⊙(K km s^-1 pc²)^-1 for the local spirals, ULIRGs, smm, and high-z star-forming galaxies (Fig. 18); the model CO conversion factor of the ULIRGs is a factor of two higher than the value derived by Downes & Solomon (1998), the CO conversion factor of the high-star-forming galaxies is a factor of two lower than assumed by Genzel et al. (2010) and Tacconi et al. (2010).

• 7.

  The model HCN-dense gas conversion factor is α_HCN = 21 ± 6M_⊙(K km s^-1 pc²)^-1 for the local spiral galaxies and ULIRGs; this is a factor of two higher than the value used in the literature (e.g., Gao & Solomon 2004); the model HCN-dense gas conversion factor is 33 ± 17M_⊙(K km s^-1 pc²)^-1 for the smm-galaxies, and 59 ± 21M_⊙(K km s^-1 pc²)^-1 for the high-z star-forming galaxies.

• 8.

  Both, the HCN and HCO⁺ emission trace the dense molecular gas to a factor of approximately 2 for the local spiral galaxies, ULIRGs and smm-galaxies.

We tested the influence of constant abundances (Sect. 7.1), Toomre Q (Q = 1 instead of Q = 1.5; Sect. 7.2), and of the scale parameter δ (δ = 15 instead of δ = 5; Sect. 7.3). The changes in molecular emission are minor (<0.2 dex). The Q = 1 and δ = 15 model overestimate the CO SLEDs of the ULIRGs. The Q = 1 model yields a lower HCN(1–0) emission than the Q = 1.5 model. Since the Q = 1.5 model already underestimates the HCN(1–0) luminosity with respect to observations, the Q = 1.5 model is our preferred model.

Whereas the CO emission is robust against the variation of model parameters and chemistry, the HCN and HCO⁺ emission is most sensitive to the chemistry of the ISM.

Within the model framework ~60% of the CO(1–0) emission of local spiral galaxies and high-z star-forming galaxies is emitted in non-self-gravitating clouds. This fraction increases to ~80% for ULIRGs. Whereas ~80% of the HCN(1–0) and HCO⁺ emission originates in non-self-gravitating clouds, this fraction decreases to ~30% for local spirals and high-z star-forming galaxies (Sect. 8.1).

The resulting CO, HCN, and HCO⁺ line emission does not change significantly if cloud-substructure is not taken into account (Sect. 8.2). Ignoring the cosmic ray heating (Sect. 8.3) leads to low CO(1–0) emission from high-z star-forming galaxies and low HCN(1–0) emission from ULIRGs.

The effect of HCN infrared pumping is small but measurable (20–30%; Sect. 8.4).

The gas velocity dispersion varies significantly with the Toomre Q parameter. The Q = 1.5 model yields high-velocity dispersions (v_disp ≫ 10 km s^-1) consistent with available observations of high-z star-forming galaxies and ULIRGs (Fig. 29). However, we note that these high-velocity dispersions may not be mandatory for starburst galaxies (Fig. 30).

The model yields molecular star-formation laws ( –Σ_H2) with slopes of 1.5/2 for the local spiral galaxies, ULIRGs, and high-z star-forming galaxies/smm galaxies. The model star-formation rate per unit area is, as observed, proportional to the molecular gas surface density divided by the dynamical timescale (lower panel of Fig. 32). There is a pronounced offset between the –Σ_H2Ω relations of the local spirals on the one hand, and ULIRGs, smm, and high-z star-forming galaxies on the other.

We conclude that our relatively simple analytic model (Sect. 2), together with the recipes for the molecular line emission (Fig. 2), captures the essential physics of galactic clumpy gas disks.

---

Acknowledgments

We would like to thank the anonymous referee for his/her suggestion to include the IR radiative transfer into our model and for improving the manuscript significantly. P.G. acknowledges funding by the European Research Council (Starting Grant 3DICE, grant agreement 336474, PI: V. Wakelam). This research has made use of the SIMBAD database and the VizieR catalogue access tool, operated at CDS, Strasbourg, France. This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.

References

Appendix A: Mean dust temperatures of molecular clouds

As described in Sect. 2.1.3, we assume that the main radiative heating agent is the UV emission component of the interstellar radiation field. For this emission component, we use an attenuation factor based on the mean extinction of a sphere of constant density (Eq. (12)). For high optical depths, Eq. (12) becomes , that is, the attenuation decreases to very low values. However, when the molecular clouds become optically thick in the near-infrared (at τ_V ~ 10), radiative transfer effects become important: a significant infrared radiation field builds up, which heats the dust until deep into the molecular clouds. To take this effect into account, we calculated 1D radiation transfer models for spherical clouds of constant density that are illuminated by an interstellar radiation field. For this purpose, we used the code TRANSPHERE (Dullemond et al. 2002). The clouds have a size of 20 pc and densities of n = (10¹,10²,10³,10⁴,10⁵) cm^-3. The clouds are illuminated by an interstellar radiation field of (10,100,1000) times the Galactic ISRF at the solar radius. For each cloud model, the resulting infrared SED was fitted by a modified Planck function to give the effective dust temperature of the cloud (Fig. A.1). According to Eqs. (13) and (15), the relation between the dust temperature and I/I(0) is T_dust ∝ (I/I(0))^1/5.5 (solid line for τ_V ≤ 10 and dashed line for τ_V> 10 in Fig. A.1). It turned out that the effective

dust temperatures obtained from the radiative transfer models for τ_V> 10 are well-fitted by a constant attenuation factor I/I(0) = 0.246.

Appendix B: Galaxy samples

Appendix C: Infrared SEDs

All Tables

All Figures