SIMULATION OF THE COSMIC EVOLUTION OF ATOMIC AND MOLECULAR HYDROGEN IN GALAXIES

Before addressing the subcomposition of cold hydrogen, we note that the total cold hydrogen mass M_H can be inferred from the total cold gas mass M_gas by a constant factor M_H = 0.74 M_gas, which corresponds to the universal abundance of hydrogen (e.g., Arnett 1996) that changes insignificantly with cosmic time. The remaining gas is composed of helium (He) and a minor fraction of heavier elements, collectively referred to as metals (Z). The DeLucia catalog gives an estimate for the metal mass in cold gas M_Z, and hence we shall compute the masses of cold hydrogen and He as

$$\begin{equation} \begin{array}{rcl} M_{\rm H}&\!\! =\!\! & (M_{\rm gas}-M_{Z})\times \beta, \\ M_{\rm He}&\!\! =\!\! & (M_{\rm gas}-M_{Z})\times (1-\beta), \end{array} \end{equation} \tag{ 2 }$$

where the hydrogen fraction β = 0.75 is chosen slightly above 0.74 to account for the subtraction of the 1%–2% metals in Equation (2).

The subdivision of the cold hydrogen mass $M_{\rm H}=M_{\rm H\,\hbox{\scriptsize {\sc i}}}+M_{{\rm H}_2}$ depends on the galaxy and evolves with cosmic time. We shall tackle this complexity using the variable H₂/H i ratio $R_{\rm mol}^{\rm galaxy}\equiv M_{{\rm H}_2}/M_{\rm H\,\hbox{\scriptsize {\sc i}}}$, hence

$$\begin{equation} \begin{array}{rcl} M_{\rm H\,\hbox{\scriptsize {\sc i}}}&\!\! =\!\! & M_{\rm H}\times \big(1+R_{\rm mol}^{\rm galaxy}\big)^{-1}, \\ M_{{\rm H}_2}&\!\! =\!\! & M_{\rm H}\times \big(1+{R_{\rm mol}^{\rm galaxy}}^{-1}\big)^{-1}. \end{array} \end{equation} \tag{ 3 }$$

Detailed observations of H i and CO in nearby regular spiral galaxies revealed that virtually all cold gas of these galaxies resides in flat, often approximately axially symmetric, disks (e.g., Walter et al. 2008; Leroy et al. 2008). CO maps recently obtained for five nearby elliptical galaxies (Young 2002) show that even these galaxies, who carry most of their stars are in a spheroid, have most of their cold gas in a disk. There is also empirical evidence that most cold gas in high-redshift galaxies resides in disks (e.g., Tacconi et al. 2006). Based on these findings, we assume that galaxies generally carry their cold atomic and molecular cold gas in flat disks with axially symmetric surface density profiles $\Sigma _{\rm H\,\hbox{\scriptsize {\sc i}}}(r)$ and $\Sigma _{\rm H_2}(r)$, where r denotes the galactocentric radius in the plane of the disk. Using these functions, R^galaxy_mol can be expressed as

$$\begin{equation} R_{\rm mol}^{\rm galaxy}= \frac{2\pi \int _0^\infty { d}r\,r\,\Sigma _{\rm H_2}(r)}{2\pi \int _0^\infty { d}r\,r\,\Sigma _{\rm H\,\hbox{\scriptsize {\sc i}}}(r)}. \end{equation} \tag{ 4 }$$

To solve Equation (4), we shall now derive an analytic model for $\Sigma _{\rm H\,\hbox{\scriptsize {\sc i}}}(r)$ and $\Sigma _{\rm H_2}(r)$. To this end, we analyzed the observed density profiles $\Sigma _{\rm H\,\hbox{\scriptsize {\sc i}}}(r)$ and $\Sigma _{\rm H_2}(r)$ presented by Leroy et al. (2008) for 12 nearby spiral galaxies of The H i Nearby Galaxy Survey (THINGS).⁵ In general, the surface density of the total hydrogen component (H i + H₂) is well fitted by a single exponential profile,

$$\begin{equation} \Sigma _{\rm H\,\hbox{\scriptsize {\sc i}}}(r)+\Sigma _{\rm H_2}(r) = \tilde{\Sigma }_{\rm H}\exp (-r/r_{\rm disk}), \end{equation} \tag{ 5 }$$

where r_disk is a scale length and $\tilde{\Sigma }_{\rm H}\equiv M_{\rm H}/(2\pi r_{\rm disk}^2)$ is a normalization factor, which can be interpreted the maximal surface density of the cold hydrogen disk.

Equation (5) can be solved for $\Sigma _{\rm H\,\hbox{\scriptsize {\sc i}}}(r)$ and $\Sigma _{\rm H_2}(r)$, if we know the local H₂/H i ratio in the disk, i.e., the radial function $R_{\rm mol}(r)\equiv \Sigma _{\rm H_2}(r)/\Sigma _{\rm H\,\hbox{\scriptsize {\sc i}}}(r)$. Following the theoretical prediction that R_mol(r) scales as some power of the gas pressure (Elmegreen 1993), Blitz & Rosolowsky (2006) presented compelling observational evidence for this power law based on 14 nearby spiral galaxies of various types. Perhaps the most complete empirical study of R_mol(r) today has recently been published by Leroy et al. (2008), who analyzed the correlations between R_mol(r) and various disk properties in 23 galaxies of the THINGS catalog. This study confirmed the power-law relation between R_mol(r) and pressure. On theoretical grounds, Krumholz et al. (2009) argued that R_mol(r) is most fundamentally driven by density rather than pressure. However, by virtue of the thermodynamic relation between pressure and density, it is, in the context of this paper, irrelevant which quantity is considered, and the density law for R_mol(r) by Krumholz et al. (2009) is indeed consistent with the pressure laws by Blitz & Rosolowsky (2006) and Leroy et al. (2008). Here, we shall apply the pressure law

$$\begin{equation} R_{\rm mol}(r) = [P(r)/P_\ast ]^\alpha, \end{equation} \tag{ 6 }$$

where P(r) is the kinematic midplane pressure outside molecular clouds, and P_* = 2.35 × 10⁻¹³ Pa and α = 0.8 are empirical values adopted from Leroy et al. (2008).

Elmegreen (1989) showed that the equations of hydrostatic equilibrium for an infinite thin disk with gas and stars exhibit a simple approximate solution for the macroscopic kinematic midplane pressure P(r) of the ISM,

$$\begin{equation} P(r) = \frac{\pi }{2}\,G\,\Sigma _{\rm gas}(r)\big (\Sigma _{\rm gas}(r)+f_{\sigma }(r)\,\Sigma _{\rm stars}^{\rm disk}(r)\big), \end{equation} \tag{ 7 }$$

where G is the gravitational constant, Σ_gas(r) is the surface density of the total cold gas component (H i + H₂ + He+ metals), Σ^disk_stars(r) is the surface density of stars in the disk (thus excluding the bulge stars of early-type spiral galaxies and elliptical galaxies), and f_σ(r) ≡ σ_gas/σ_stars,z is the ratio between the vertical velocity dispersions of gas and stars. The impact of supernovae and other small-scale effects on the gas pressure are implicitly included in Equation (7) via the velocity dispersion σ_gas. For Σ^disk_stars = 0, Equation (7) reduces to P(r) = 0.5 π G Σ_gas(r)², which is sometimes used as an approximation for the ISM pressure in gas-rich galaxies (e.g., Crosthwaite & Turner 2007).

To simplify Equation (7), we note that Σ_gas(r) can be expressed as Σ_gas(r) = M_gas/(2π r²_disk)exp(−r/r_disk), which is identical to Equation (5) up to the constant factor correcting for helium and metals. To find a similar expression for Σ^disk_stars(r), we analyzed the stellar surface densities Σ_stars(r) of the 12 THINGS galaxies mentioned before. In agreement with many other studies (e.g., Courteau et al. 1996), we found that Σ_stars(r) is generally well approximated by a double-exponential profile, i.e., the sum of an exponential profile Σ^bulge_stars for the bulge and an exponential profile Σ^disk_stars for the disk. On average, the scale length of the stellar disk $\tilde{r}_{\rm disk}$ is 30%–50% smaller than the gas scale length r_disk, which traces the fact that stars form in the more central H₂-dominated parts of galaxies. Indeed, several observational studies revealed that the stellar scale length is nearly identical to that of molecular gas (e.g., Young et al. 1995; Regan et al. 2001; Leroy et al. 2008), and hence smaller than the scale length of H i or the scale length r_disk of the total cold gas component. For simplicity, we shall here assume $r_{\rm disk}=2\,\tilde{r}_{\rm disk}$ for all galaxies, such that Σ^disk_stars(r) = 4 M_stars/(2π r²_disk)exp(−2 r/r_disk). Finally, we approximate the dispersion ratio f_σ(r) as f_σ(r) = f⁰_σ exp(r/r_disk), where f⁰_σ is a constant. This approximation is motivated by empirical evidence that the gas dispersion σ_gas remains approximately constant across galactic discs (e.g., Dickey et al. 1990; Boulanger & Viallefond 1992; Leroy et al. 2008), combined with theoretical and observational studies showing that the stellar velocity dispersion σ_stars,z decreases approximately exponentially with a scale length twice that of the stellar surface density (e.g., Bottema 1993). Within those approximations, Equation (7) reduces to

$$\begin{equation} P(r)\approx \frac{GM_{\rm gas}}{8\pi \,r_{\rm disk}^4}\big (M_{\rm gas}+\langle f_{\sigma }\rangle M_{\rm stars}^{\rm disk}\big)\exp (-2\,r/r_{\rm disk}), \end{equation} \tag{ 8 }$$

where 〈f_σ〉 ≡ 4f⁰_σ is a constant, which can be interpreted as the average value of f_σ(r) weighted by the stellar surface density, since ∫2π r Σ^disk_stars(r)f_σ(r)/M^disk_stars = 4f⁰_σ.

Substituting P(r) in Equation (6) for Equation (8), we obtain

$$\begin{equation} \frac{\Sigma _{\rm H_2}(r)}{\Sigma _{\rm H\,\hbox{\scriptsize {\sc i}}}(r)} \equiv R_{\rm mol}(r) = R_{\rm mol}^{c}\,\exp (-1.6\,r/r_{\rm disk}) \end{equation} \tag{ 9 }$$

with

$$\begin{equation} R_{\rm mol}^{c}= \big[K\,r_{\rm disk}^{-4}\,M_{\rm gas}\,\big (M_{\rm gas}+\langle f_{\sigma }\rangle \,M_{\rm stars}^{\rm disk}\big)\big]^{0.8}, \end{equation} \tag{ 10 }$$

where $K\equiv G/(8\pi \,P_\ast)=11.3\rm \;m^4\;kg^{-2}$. Equation (9) reveals that the H₂/H i ratio R_mol(r) is described by an exponential profile with scale length r_disk/1.6. It should be emphasized that the central value R^c_mol does not necessarily correspond to the H₂/H i ratio measured at the center of real galaxies due to an additional H₂ enrichment caused by the central stellar bulge. However, R^c_mol represents the extrapolated cental H₂/H i ratio of the exponential profile, which approximates R_mol(r) in the outer, disk-dominated galaxy parts.

We can now solve Equations (5) and (9) for the atomic and molecular surface density profiles, i.e.,

$$\begin{eqnarray} \Sigma _{\rm H\,\hbox{\scriptsize {\sc i}}}(r) & = & \frac{\tilde{\Sigma }_{\rm H}\,\exp (-r/r_{\rm disk})}{1+R_{\rm mol}^{c}\exp (-1.6\,r/r_{\rm disk})},\, \end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray} \Sigma _{\rm H_2}(r) & = & \frac{\tilde{\Sigma }_{\rm H}\,R_{\rm mol}^{c}\,\exp (-2.6\,r/r_{\rm disk})}{1+R_{\rm mol}^{c}\exp (-1.6\,r/r_{\rm disk})}.\, \end{eqnarray} \tag{ 12 }$$

These model profiles can be checked against the observed H i - and H₂-density profiles of the nearby galaxies analyzed by Leroy et al. (2008). In particular, we can test the limitations implied by our assumption that $r_{\rm disk}=2\,\tilde{r}_{\rm disk}$. To this end, we selected two observed regular spiral galaxies with $r_{\rm disk}\approx 2\,\tilde{r}_{\rm disk}$ (NGC 3184) and $r_{\rm disk}\approx \tilde{r}_{\rm disk}$ (NGC 5505). To evaluate Equations (11) and (12) for those galaxies, we require the quantities M^disk_stars, M_gas, r_disk, and 〈f_σ〉. M^disk_stars, M_gas, and r_disk were determined by fitting a single exponential profile to the total cold gas component and a double-exponential profile (bulge and disk) to the stellar component, and 〈f_σ〉 was chosen as 〈f_σ〉 = 0.4, i.e., the value given by Elmegreen (1993) for nearby galaxies. As shown in Figure 1, the resulting model profiles $\Sigma _{\rm H\,\hbox{\scriptsize {\sc i}}}(r)$ and $\Sigma _{\rm H_2}(r)$ approximately match the empirical data. The fact that the fit is rather good for both galaxies demonstrates that the quality of the model predictions does not sensibly depend on the goodness of the model assumption $r_{\rm disk}\approx 2\,\tilde{r}_{\rm disk}$. Similarly good fits are indeed found for most of the 12 THINGS galaxies, for which Leroy et al. (2008) published radial H i - and H₂-density profiles.

Zoom In Zoom Out Reset image size

Equations (11) and (12) can be solved for the maximal surface densities of H i and H₂. $\Sigma _{\rm H\,\hbox{\scriptsize {\sc i}}}(r)$ exhibits its maximum at the radius $r_{\rm H\,\hbox{\scriptsize {\sc i}}}^{\max }=0.625\,r_{\rm disk}\,\ln (3/5\times R_{\rm mol}^{c})$, as long as R^c_mol>5/3. Galaxies in this category show an H i drop toward their center, such as observed in most galaxies in the THINGS catalog (Walter et al. 2008). By contrast, disk galaxies with R^c_mol ⩽ 5/3 have H i-density profiles peaking at the center, $r_{\rm H\,\hbox{\scriptsize {\sc i}}}^{\max }=0$. Galaxies with such small values of R^c_mol have low gas densities by virtue of Equation (10), such as the irregular galaxies NGC 4214 and NGC 3077 (see profiles in Leroy et al. 2008). $\Sigma _{\rm H_2}(r)$ given in Equation (12) always peaks at the disk center, $r_{\rm H_2}^{\max }=0$.

The maximal values of $\Sigma _{\rm H\,\hbox{\scriptsize {\sc i}}}(r)$ and $\Sigma _{\rm H_2}(r)$, called $\Sigma _{\rm H\,\hbox{\scriptsize {\sc i}}}^{\rm max}\equiv \Sigma _{\rm H\,\hbox{\scriptsize {\sc i}}}(r_{\rm H\,\hbox{\scriptsize {\sc i}}}^{\max })$ and $\Sigma _{\rm H_2}^{\rm max}\equiv \Sigma _{\rm H_2}(r_{\rm H_2}^{\max })$, can be computed as

$$\begin{eqnarray} \Sigma _{\rm H\,\hbox{\scriptsize {\sc i}}}^{\rm max}\big /\tilde{\Sigma }_{\rm H}&=& \left\lbrace \begin{array}{@{}ll@{}} 1/(1+R_{\rm mol}^{c}) & \rm { if }\quad R_{\rm mol}^{c}\le 5/3\quad \\ 0.516\,{R_{\rm mol}^{c}}^{-5/8} & \rm { if }\quad R_{\rm mol}^{c}>5/3, \end{array}\right. \end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray} \Sigma _{\rm H_2}^{\rm max}\big /\tilde{\Sigma }_{\rm H}&=& R_{\rm mol}^{c}\big /\big(1+R_{\rm mol}^{c}\big). \end{eqnarray} \tag{ 14 }$$

The density profiles of Equations (11) and (12) can be substituted into Equation (4). The exact solution of Equation (4) is quite unhandy, but R^galaxy_mol only depends on R^c_mol and an excellent approximation, accurate to better than 5% over the 9 orders of magnitude R^c_mol = 10⁻³to10⁶ (covering the most extreme values at all redshifts), is given by

$$\begin{equation} R_{\rm mol}^{\rm galaxy}= \big (3.44\,{R_{\rm mol}^{c}}^{-0.506}+4.82\,{R_{\rm mol}^{c}}^{-1.054}\big)^{-1}. \end{equation} \tag{ 15 }$$

Equations (10) and (15) constitute a physical prescription to estimate the H₂/H i ratio of any regular galaxy based on four global quantities: the disk stellar mass M^disk_stars, the cold gas mass M_gas, the scale radius of the cold gas disk r_disk, and the dispersion parameter 〈f_σ〉. In Obreschkow & Rawlings (2009), we showed that the H₂/H i ratios inferred from this model are consistent with observations of nearby galaxies. Moreover, this model presumably extends to high redshifts, since it essentially relies on the fundamental relation between pressure and molecular fraction and on a few other physical assumptions with weak or absent dependence on cosmic epoch. However, a critical discussion of the limitations of this model is presented Sections 6.2 and 6.3.

We applied the model given in Equations (10) and (15) together with Equations (2) and (3) to the DeLucia catalog in order to assign H i, H₂, and He masses to the simulated galaxies. The quantities M^disk_stars and M_Z used in these Equations are directly contained in the DeLucia catalog, and M_gas was inferred from the given cold gas masses via the correction of Equation (1). The dispersion parameter 〈f_σ〉 is approximated by the constant 〈f_σ〉 = 0.4, consistent with the local observational data used by (Elmegreen 1989, but see discussion in Section 6.3).

The remaining and most subtle ingredient for our prescription of Equations (10) and (15) is the scale radius r_disk. This radius can be estimated from the virial radius r_vir of the parent halo, but their relation is intricate. Even modern N-body plus smoothed particle hydrodynamics (SPH) simulations of galaxy formation still struggle at reproducing observed disk diameters (Kaufmann et al. 2007), and hence the more simplistic semianalytic approaches are likely to require some empirical adjustment. Mo et al. (1998) studied the case of a flat exponential disk in an isothermal singular halo. When assuming that the disk's mass can be neglected for its rotation curve, they find

$$\begin{equation} r_{\rm disk}= \frac{\lambda \times \xi }{\sqrt{2}}\,r_{\rm vir}, \end{equation} \tag{ 16 }$$

where λ is the spin parameter of the halo and ξ is the ratio between the specific angular momentum of the disk (angular momentum per unit mass) and the specific angular momentum of the halo. The model behind Equation (16) does not distinguish between different scale radii for stars and cold gas, but it assumes a single exponential disk in hydrostatic equilibrium without including the effects of star formation. It is therefore natural to identify r_disk in Equation (16) with the cold gas scale radius r_disk of Equations (10)–(12).

In the Millennium Simulation, r_vir was calculated from the virial mass M_vir using the relation

$$\begin{equation} M_{\rm vir}= \frac{4\pi }{3}r_{\rm vir}^3 \times 200 \rho _{c}(z), \end{equation} \tag{ 17 }$$

where ρ_c(z) is the critical density for closure ρ_c(z) = 3H²(z)/(8πG) and M_vir is the virial mass of the halo. For central halos, M_vir was approximated as M₂₀₀, i.e., the mass in the region with an average density equal to 200ρ_c(z); and for subhalos, M_vir was approximated as the total mass of the gravitationally bound simulation particles.

The spin parameter λ was calculated directly from the N-body Millennium Simulation according to the definition λ ≡ J_haloE_halo^1/2G⁻¹M_halo^−5/2, where J_halo denotes the angular momentum of the halo, E_halo its energy, and M_halo its total mass. For the satellite galaxies in the DeLucia catalog, i.e., the ones without halo (see Section 2), the value of λ × r_vir was approximated as the respective value of the original galaxy halo before its disappearance.

The only missing parameter for the calculation of r_disk via Equation (16) is the angular momentum ratio ξ. It is of order unity for evolved disk galaxies (e.g., Fall & Efstathiou 1980; Zavala et al. 2008), but its exact value is uncertain because of the difficulty of measuring the spin of dark matter halos and convergence issues in numerical simulations. For example, Kaufmann et al. (2007) showed that N-body plus SPH simulations with as many as 10⁶ particles per galaxy do not reach convergence in angular momentum, because of the difficulties to model the transport of angular momentum.

Here, we shall chose ξ, such that our simulation reproduces the empirical relation between the galaxy baryon mass M_bary and the stellar scale radius $\tilde{r}_{\rm disk}\approx r_{\rm disk}/2$, measured in the local universe. To ensure consistency with Section 3.1, we use again the data from the THINGS galaxies analyzed by Leroy et al. (2008). Of the 23 galaxies in this sample, we reject the six irregular objects, since our model for H i and H₂ assumes regular galaxies. The remaining 17 galaxies cover all spiral types and include H i-rich and H₂-rich galaxies. For each galaxy, we adopted the stellar scale radii $\tilde{r}_{\rm disk}$ and baryon masses $M_{\rm bary}=M_{\rm stars}+M_{\rm H\,\hbox{\scriptsize {\sc i}}}+M_{{\rm H}_2}$ directly from the data presented in Leroy et al. (2008).⁶ The Hubble types T, i.e., the numerical stage indexes along the revised Hubble sequence of the RC2 system (de Vaucouleurs et al. 1976), were drawn from the HyperLeda database (Paturel et al. 2003).

The resulting empirical relation between M_bary and $\tilde{r}_{\rm disk}$ is displayed in Figure 2. The scatter of these data probably underestimates the true scatter caused by all galaxies, since we exclusively considered non-interacting regular spiral galaxies. The zero point of the mean relation between M_bary and $\tilde{r}_{\rm disk}$ depends on the morphological galaxy type, as is revealed by the distinction of early- and late-type spiral galaxies in Figure 2. Given identical baryon masses, the stellar scale radii of early-type galaxies tend to be smaller than the radii of late-type galaxies. This trend was also detected in other data samples (e.g., data from Kregel et al. 2002 shown in Obreschkow & Rawlings 2009). Several reasons could explain this finding: (1) early-type galaxies have more massive stellar bulges, which present an additional central potential that contracts the disk; (2) bulges often form from disk instabilities, occurring preferably in systems with relatively low angular momentum, and hence early-type galaxies are biased toward smaller angular momenta and smaller scale radii; (3) larger bulges like those of lenticular and elliptical galaxies, often arise from galaxy mergers, which tend to reduce the specific angular momenta and scale radii.

Zoom In Zoom Out Reset image size

To parameterize the dependence of the scale radius on the morphological galaxy type, the latter shall be quantified using the stellar mass fraction of the bulge, $\mathcal {\bm B}\equiv M_{\rm stars}^{\rm bulge}/M_{\rm stars}$. In the observed sample, the values of $\mathcal {\bm B}$ can be approximately inferred from the Hubble type T. Here, we shall use the relation

$$\begin{equation} \mathcal {\bm B}= (10-T)^2/256, \end{equation} \tag{ 18 }$$

which approximately parameterizes the mean behavior of 146 moderately inclined barred and unbarred local spiral galaxies analyzed by Weinzirl et al. (2009). Equation (18) satisfies the boundary conditions $\mathcal {\bm B}=1$ for T = −6 (i.e., pure spheriods) and $\mathcal {\bm B}=0$ for T = 10 (pure disks).

We shall approximate the relation between M_bary and $\tilde{r}_{\rm disk}$ as a power law with an additional term for the observed secondary dependence on morphological type,

$$\begin{equation} \log \left(\frac{\tilde{r}_{\rm disk}}{\rm kpc}\right) = a_0+a_1\log \left(\frac{M_{\rm bary}}{10^{10}\;{M}_{\odot }}\right)+a_2\,\mathcal {\bm B}, \end{equation} \tag{ 19 }$$

where a₀, a₁, and a₂ are free parameters. The best fit to the empirical data in terms of a maximum likelihood approach is given by the choice a₀ = 0.3, a₁ = 0.4, and a₂ = −0.6.

The fit of Equation (19) is displayed in Figure 2 for early-type spiral galaxies ($T=2\leftrightarrow \mathcal {\bm B}=0.25$, solid line) and for late-type spiral galaxies ($T=10\leftrightarrow \mathcal {\bm B}=0$, dashed line). For comparison, the mean power-law relations for the early- and late-type spiral galaxies in the DeLucia catalog at z = 0 are displayed as dotted and dash-dotted lines for the choice ξ = 1. The simulated scale radii using ξ = 1 are consistent with the observed ones for very massive galaxies (M_bary ≈ 10¹¹ M_☉), but less massive galaxies in the simulation turn out slightly too large if ξ = 1. Furthermore, the morphological dependence of the simulation is too small compared to the observations. This can be corrected ad hoc by introducing a variable ξ that depends on both M_bary and $\mathcal {\bm B}$, i.e.,

$$\begin{equation} \log (\xi) = b_0+b_1\log \left(\frac{M_{\rm bary}}{10^{10}\;{M}_{\odot }}\right)+b_2\,\mathcal {\bm B}, \end{equation} \tag{ 20 }$$

where b₀, b₁, and b₂ are free parameters. The parameters minimizing the rms deviation between the simulated galaxies and the empirical model of Equation (19) are b₀ = −0.1, b₁ = 0.3, and b₂ = −0.6. In addition, we chose a lower limit for ξ equal to 0.5, in order to prevent unrealistically small-scale radii.

We emphasize that Equation (20) is merely an empirical correction; this choice of ξ should not be considered as an estimate of the true ratio between the specific angular momenta of the disk and the halo, but it also accounts for the imperfection of the simplistic halo model by Mo et al. (1998), for missing physics in the semianalytic modeling, and for possible systematic errors in the spin parameters λ of the Millennium Simulation. The average value of Equation (20) over all galaxies in the DeLucia catalog is 〈ξ〉 = 0.7 (with σ = 0.2), which is approximately consistent with ξ ≈ 1 of modern high-resolution simulations of galaxy formation (e.g., Zavala et al. 2008), even though the latter still suffer from issues with the transport of angular momentum as mentioned above.

Using Equations (16) and (20), we estimated a scale radius r_disk for each galaxy in the DeLucia catalog. A sample of 10² simulated early- and late-type spiral galaxies at z = 0 is shown in Figure 2. Given r_disk as well as M_stars, M_gas, and 〈f_σ〉 = 0.4, we then applied Equations (2), (3), (10), and (15) in order to subdivide the non-metallic cold gas mass (M_gas − M_Z) of each galaxy into H i, H₂, and He.

We shall now compare the H i and H₂ masses predicted by our model of Sections 3.1 and 3.2 to recent observations in the local universe. From the viewpoint of the simulation, a fundamental output are the MFs of H i and H₂, while the available observational counterparts are the LFs of the H i-emission line (Zwaan et al. 2005) and the CO(1–0)-emission line (Keres et al. 2003). Therefore, either the simulated data or the observed data need a luminosity-to-mass (or vise versa) conversion to compare the two. Section 3.3 focuses on the MFs, adopting the standard luminosity-to-mass conversion for H i used by Zwaan et al. (2005) and the CO-luminosity-to-H₂-mass conversion of Obreschkow & Rawlings (2009). As a complementary approach, Section 3.4 will focus on the LFs, which will require a model for the conversion of simulated H₂ masses into CO luminosities.

We define the MFs as ϕ_x(M_x) ≡ dρ_x/dlog M_x, where ρ_x(M_x) is the space density (number per comoving volume) of galaxies containing a mass M_x of the constituent x (H i, H₂, He, etc.). Given a mass, such as $M_{\rm H\,\hbox{\scriptsize {\sc i}}}$, for each galaxy in the DeLucia catalog, the derivation of the corresponding MF only requires the counting of the number of sources per mass interval. We chose 60 mass intervals, logarithmically spaced between 10⁸ M_☉ and 10¹¹ M_☉, giving about ∼10⁶ galaxies per mass interval in the central mass range, while keeping the mass error relatively small (Δlog(M) < 0.05). Since MFs combine units of mass and length, they generally depend on the Hubble constant H₀, or the dimensionless Hubble parameter h, defined as H₀ = 100 h km s⁻¹ Mpc⁻¹. Although MFs are often plotted in units making no assumption on h (e.g., M_☉ h⁻² for the mass scale), this is impossible when observations are compared to cosmological simulations. The reason is that simulated masses in the Millennium Simulation scale to first order as h⁻¹, whereas empirical masses, when determined from electromagnetic fluxes, are proportional to the square of the distance and hence scale as h⁻². For all plots in this paper, we shall therefore use h = 0.73, which corresponds to the value adopted by the Millennium Simulation (see Section 2).

Figure 3 displays the H i and H₂ MF of our simulation (solid lines), as well as the corresponding empirical MFs for the local universe (points with error bars). The empirical H i MF was obtained by Zwaan et al. (2005) based on 4315 galaxies of the HI-Parkes All Sky Survey (HIPASS) and the empirical H₂ MF was derived in Obreschkow & Rawlings (2009) from the CO LF presented by Keres et al. (2003). Both empirical MFs approximately match the simulated data. We note, however, that the consistency between observation and simulation decreases if we skip the overall correction of the cold gas masses in the DeLucia catalog by the constant factor ζ (dotted lines), which, as argued in Section 2 can be justified by a fraction of the disk gas being electronically excited or ionized.

Zoom In Zoom Out Reset image size

Our simulation slightly overpredicts the observed number of the largest H i and H₂ masses, i.e., the ones in the exponential tail of the MFs in Figure 3. These tails contain the most massive systems, whose emergent luminosities are most likely to be biased by opacity and thermal effects. Including these effects would probably correct the space density of massive systems toward the simulated MFs. Additionally, we note that the presented empirical MFs neglect mass measurement errors, which might have an important effect on the slope of the exponential tails. Another difference between observations and simulation are the spurious bumps in the low-mass range of the simulated MFs, i.e., $\log (M_{\rm H\,\hbox{\scriptsize {\sc i}}})\approx 8.5$ and $\log (M_{{\rm H}_2})\approx 8.0$, where the number density is about doubled compared to observations. This feature can also be seen in the optical b_J-band LF shown by Croton et al. (2006) and stems from an imprecision in the number density of the smallest galaxies in the DeLucia catalog, where the mass resolution of the Millennium Simulation implies a poorly resolved merger history. The overdensity of sources around this resolution limit roughly balances the mass of even smaller galaxies, i.e., $\log (M_{\rm H\,\hbox{\scriptsize {\sc i}}}/{M}_{\odot })\ll 8.0$, that are missing in the simulation.

The universal gas densities of the simulation, expressed relative to the critical density for closure, are $\Omega _{\rm H\,\hbox{\scriptsize {\sc i}}}^{\rm sim}=3.4\times 10^{-4}$ and $\Omega _{{\rm H}_2}^{\rm sim}=1.1\times 10^{-4}$, in good agreement with the observations $\Omega _{\rm H\,\hbox{\scriptsize {\sc i}}}^{\rm obs}=(3.6\pm 0.4)\times 10^{-4}$ (Zwaan et al. 2005) and $\Omega _{{\rm H}_2}^{\rm obs}=(0.95\pm 0.37)\times 10^{-4}$ (Obreschkow & Rawlings 2009).

Figure 4 shows our simulated H i MF and H₂ MF together with the MF for the cold gas metals given in the original DeLucia catalog and the MF for He as trivially derived using Equation (2). This picture reveals that in the cold gas of the local universe He is probably more abundant than H₂, but less abundant than H i.

Zoom In Zoom Out Reset image size

We shall now consider the H i masses in elliptical and spiral galaxies separately. This division is based on the Hubble type T, where we consider galaxies with T < 0 as "ellipticals" and galaxies with T ⩾ 0 as "spirals"—a rough separation that neglects other types like irregular galaxies as well as various subclassifications. In the simulation, T is computed from the bulge mass fraction according to Equation (18).

The simulated H i MFs of both elliptical and spiral galaxies are shown in Figure 5. In order to determine the observational counterparts, we split the HIPASS galaxy sample into elliptical and spiral galaxies according to the Hubble types provided in the HyperLeda reference database (Paturel et al. 2003). For both subsamples, the H i MF was evaluated using the 1/V_max method (Schmidt 1968), where V_max was estimated from the analytic completeness function for HIPASS, which characterizes the completeness of each source given its H i-peak flux density S_p and integrated H i-line flux S_int (Zwaan et al. 2004). In order to estimate the uncertainties of the MFs, we derived them for 10⁴ random half-sized subsets of the HIPASS sample—a bootstrapping approach. The standard deviation of the 10⁴ values of $\log (\phi _{\rm H\,\hbox{\scriptsize {\sc i}}})$ for each mass bin was divided by $\sqrt{2}$ to estimate the 1σ errors of $\log (\phi _{\rm H\,\hbox{\scriptsize {\sc i}}})$ for the full sample.

Zoom In Zoom Out Reset image size

Figure 5 demonstrates that our simulation successfully reproduces the H i masses of both spiral and elliptical galaxies for H i masses greater than ∼10⁹ M_☉, although the nearly perfect match between simulation and observation may be somewhat coincidental due to the uncertainties of the Hubble types T calculated via Equation (18). For H i masses smaller than 10⁹ M_☉, the morphological separation seems to breakdown (shaded zone in Figure 5). Indeed the H i-mass range $M_{\rm H\,\hbox{\scriptsize {\sc i}}}\lesssim 10^9\;{M}_{\odot }$ approximately corresponds to the stellar mass range M_stars ≲ 4 × 10⁹ M_☉, for which morphology properties are poorly resolved (see Section 2).

The characteristic radio line of H i stems from the hyperfine energy level splitting of the hydrogen atom and lies at 1.42 GHz rest-frame frequency. The velocity-integrated luminosity of this line $L_{\rm H\,\hbox{\scriptsize {\sc i}}}$ can be calculated from the H i mass via

$$\begin{equation} \frac{M_{\rm H\,\hbox{\scriptsize {\sc i}}}}{{M}_{\odot }} = 1.88\times 10^4\times \frac{L_{\rm H\,\hbox{\scriptsize {\sc i}}}}{\rm Jy\;km\;s^{-1}\;Mpc^2}. \end{equation} \tag{ 21 }$$

Equation (21) neglects H i-self-absorption effects, but this is likely to be a problem only for the largest disk galaxies observed edge-on (Rao et al. 1995). The strict proportionality between $L_{\rm H\,\hbox{\scriptsize {\sc i}}}$ and $M_{\rm H\,\hbox{\scriptsize {\sc i}}}$ assumed in Equation (21) means that the H i LF is geometrically identical to the H i MF.

By contrast, the H₂ masses used for the empirical H₂ MFs in Figures 3 and 4 rely on measurements of the CO(1–0) line, i.e., the 115 GHz radio line stemming from the fundamental rotational relaxation of the most abundant CO isotopomer $\rm ^{12}C^{16}O$. Here, we only consider this line, but luminosities of other CO lines can be estimated using approximate empirical line ratios (e.g., Braine et al. 1993; Righi et al. 2008).

The CO(1–0)-to-H₂ conversion generally depends on the galaxy and the cosmic epoch, and it is often represented by the dimensionless factor

$$\begin{equation} X \equiv \frac{N_{{\rm H}_2}/{\rm cm}^{-2}}{I_{\rm CO}/({\rm K\;km\;s}^{-1})}\times 10^{-20}, \end{equation} \tag{ 22 }$$

where $N_{{\rm H}_2}$ is the column density of H₂ molecules and I_CO is the integrated CO(1–0)-line intensity per unit surface area defined via the surface brightness temperature in the Rayleigh–Jeans approximation. The definition of Equation (22) implies the mass–luminosity relation (e.g., review by Young & Scoville 1991)

$$\begin{equation} \frac{M_{{\rm H}_2}}{{M}_{\odot }} = 313\times X\times \frac{L_{\rm CO}}{\rm Jy\;km\;s^{-1}\;Mpc^2}, \end{equation} \tag{ 23 }$$

where L_CO is the velocity-integrated luminosity of the CO(1–0) line.

As discussed in Obreschkow & Rawlings (2009), the theoretical and observational determination of the X-factor is a subtle task with a long history. Most present-day studies assume a constant X-factor X_c, such as

$$\begin{equation} X_{c}= 2, \end{equation} \tag{ 24 }$$

which is typical for spiral galaxies in the local universe (Leroy et al. 2008). By contrast, Arimoto et al. (1996) and Boselli et al. (2002) suggested that X is variable, X_v, and approximately inversely proportional to the metallicity O/H, i.e., the ratio between the number of oxygen ions and hydrogen ions in the hot ISM. Using their data, we found that (Obreschkow & Rawlings 2009)

$$\begin{equation} \log (X_{v})=(-2.9\pm 0.2)-(1.02\pm 0.05)\log (\rm O/H). \end{equation} \tag{ 25 }$$

At first sight, the empirical negative dependence of X on the metallicity seems to contradict the fact that $\rm ^{12}C^{16}O$ is optically thick for 115 GHz radiation. Indeed, the radiated luminosity should not depend on the density of metals, as long as the latter is high enough for the radiation to remain optically thick (Kutner & Leung 1985). However, detailed theoretical investigations (e.g., Maloney & Black 1988) of the sizes and temperatures of molecular clumps were indeed able to explain, and in fact predict, the negative dependence of X on metallicity.

Equation (25) links the metallicity of the hot ISM to the X-factor of cold molecular clouds, and it is likely a consequence of a more fundamental relation between cold gas metallicity and X. To uncover such a relation, we assume that the O/H metallicity of the cold ISM in local galaxies is approximated by O/H of the hot ISM. Given an atomic mass of 16 for oxygen, the fact that hydrogen makes up a fraction 0.74 of the total baryon mass, and assuming that oxygen accounts for a fraction of 0.4 of the mass of all metals (based on Arnett 1996; Kobulnicky & Zaritsky 1999), Equation (25) translates to

$$\begin{equation} X_{v}\approx 0.04\,M_{\rm gas}/M_{Z}, \end{equation} \tag{ 26 }$$

where M_gas is the total cold gas mass and M_Z is the mass of metals in cold gas. Equation (26) only relates cold gas properties to each other and therefore is more fundamental than Equation (25).

To evaluate Equation (26) for each galaxy in the simulation, we used the cold gas metal masses M_Z given in the DeLucia catalog. Those masses are reasonably accurate as demonstrated by De Lucia et al. (2004) and Croton et al. (2006) through a comparison of the simulated stellar mass–metallicity relation to the empirical mass–metallicity relation obtained from 53,000 star-forming galaxies in the SDSS (Tremonti et al. 2004). For most galaxies at z = 0, the simulation yields metal fractions M_Z/M_gas ≈ 0.01–0.04 in the local universe, thus implying X_v ≈ 1–4 in agreement with observed values (e.g., Boselli et al. 2002).

Figure 6 displays the simulated CO LF for the variable X-factor X_v (solid line) and the constant X-factor X_c (dashed line) together with the empirical CO LF (Keres et al. 2003), adjusted to h = 0.73. The comparison supports the variable X-factor of Equation (26) against X_c = 2 (and the same conclusion is found for other constant values of X_c). Using Equation (26) also has the advantage that the cosmic evolution of the X-factor due to the evolution of metallicity is implicitly accounted for. Nevertheless, Equation (26) may not be appropriate at high redshift as discussed in Section 6.3.

Zoom In Zoom Out Reset image size

To account for the narrowness of the emission lines observed in the central gas regions of many galaxies (e.g., Sauty et al. 2003), we require a halo model with vanishing velocity at the center, as opposed to, for example, the commonly adopted singular isothermal sphere with a density ρ_halo(r) ∼ r⁻² and a constant velocity profile. We chose the Navarro–Frenk–White (NFW; Navarro et al. 1995, 1996) model, which relies on high-resolution numerical simulations of dark matter halos in equilibrium. These simulations revealed that halos of all masses in a variety of dissipationless hierarchical clustering models are well described by the spherical density profile

$$\begin{equation} \rho _{\rm halo}(r) = \rho _0\,{ [(r/r_{s})(1+r/r_{s})^2 ]}^{-1}, \end{equation} \tag{ 31 }$$

where ρ₀ is a normalization factor and r_s is the characteristic scale radius of the halo. This profile is also supported by the Hubble Space Telescope analysis of the weak lensing induced by the galaxy cluster MS 2053-04 at redshift z = 0.58 (Hoekstra et al. 2002). ρ_halo(r) varies as r⁻¹ at the halo center and continuously steepens to r⁻³ for r → ∞. It passes through the equilibrium profile of the self-gravitating isothermal sphere, i.e., ρ_halo(r) ∼ r⁻², at r = r_s.

The definition of r_vir in the Millennium Simulation given in Equation (17) implies that

$$\begin{equation} \rho _0 = \frac{200}{3}\,\frac{\rho _{c}\,c_{\rm halo}^3}{\ln (1+c_{\rm halo})-c_{\rm halo}/(1+c_{\rm halo})}, \end{equation} \tag{ 32 }$$

where c_halo ≡ r_vir/r_s is referred to as the halo concentration parameter. Most numerical models predict that c_halo scales with the virial mass M_vir, defined as the mass inside the radius r_vir, according to a power law (e.g., Navarro et al. 1997; Bullock et al. 2001; Dolag et al. 2004; Hennawi et al. 2007). Here, we shall use the result of Hennawi et al. (2007),

$$\begin{equation} c_{\rm halo}=\frac{12.3}{1+z}\left(\frac{M_{\rm vir}}{1.3\times 10^{13}\,h^{-1}\,{M}_{\odot }}\right)^{-0.13}, \end{equation} \tag{ 33 }$$

which is consistent with recent empirical values of the matter concentration in galaxy clusters derived from X-ray measurements and strong lensing data (Comerford & Natarajan 2007).

For a spherical halo, the circular velocity profile is given by V^halo_c²(r) = GM_halo(r)/r with $M_{\rm halo}(r)=4\pi \int _0^rd\tilde{r}\,\tilde{r}^2\,\rho _{\rm halo}(\tilde{r})$. Using Equations (31) and (32), this implies that

$$\begin{equation} {V_{c}^{\rm halo}}^2(x) = \frac{G\!M_{\rm vir}}{r_{\rm vir}}\times \frac{\ln (1+c_{\rm halo}x)-\frac{c_{\rm halo}x}{1+c_{\rm halo}x}}{x[\ln (1+c_{\rm halo})-\frac{c_{\rm halo}}{1+c_{\rm halo}}]}, \end{equation} \tag{ 34 }$$

where x ≡ r/r_vir (thus c_halox = r/r_s). This velocity vanishes at the halo center, then climbs to a maximal value $V_{\max }=1.65\,r_{s}\sqrt{G\rho _0}$ at r = 2.16 r_s, from where it decreases monotonically with r, typically reaching 0.65–0.95V_max at r = r_vir (x = 1) with the extremes corresponding, respectively, to c_halo = 25 and c_halo = 5. For larger radii, the velocity asymptotically approaches the point-mass velocity profile V^halo_c²(r) = GM_vir/r.

For simplicity, we assume in this section that the galactic disk is described by a single exponential surface density for stars and cold gas,

$$\begin{equation} \Sigma ^{\rm disk}(r) = \frac{M^{\rm disk}}{2\pi \,r_{\rm disk}^2}\,\exp \left(-\frac{r}{r_{\rm disk}}\right), \end{equation} \tag{ 35 }$$

where M^disk is the total disk mass, taken as the sum of the cold gas mass and the stellar mass in the disk. In most real galaxies, the stellar surface densities are slightly more compact (see Section 3.1), but we found that including this effect does not significantly modify the shape of the atomic and molecular emission lines. In fact, the radius, which maximally contributes to the disk mass, i.e., the maximum of r Σ^disk(r), is r = r_disk. Therefore, we expect the gravitational potential to differ significantly from the point-mass potential only for r of order r_disk or smaller. Applying Poisson's equation to the surface density of Equation (35), the gravitational potential in the plane of the disk becomes

$$\begin{equation} \varphi _{\rm disk}(r) = -\frac{G\!M^{\rm disk}}{2\pi r_{\rm disk}^2}\int \!\!\!\int _D\frac{\exp (-\tilde{r}/r_{\rm disk})\,\tilde{r}\,d\tilde{r}\,d\theta }{(r^2{+\tilde{r}^2}{-2}r\tilde{r}\cos \,\theta)^{1/2}}, \end{equation} \tag{ 36 }$$

where the integration surface D is given by $\tilde{r}\in [0,\infty)$, θ = [0, 2π).

The velocity profile for circular orbits in the plane of the disk can be calculated as V^disk_c² = r dφ_disk/dr. The integral in Equation (36) is elliptic, and hence there are no exact closed-from expressions for φ_disk and V^disk_c. However, in this study we numerically found that an excellent approximation is given by

$$\begin{eqnarray} {V_{c}^{\rm disk}}^2(x) &\approx & \frac{G\!M^{\rm disk}}{r_{\rm vir}} \\ &\times &\frac{c_{\rm disk}+4.8c_{\rm disk}\exp [-0.35c_{\rm disk}x-3.5/(c_{\rm disk}x)]}{c_{\rm disk}x+(c_{\rm disk}x)^{-2}+2(c_{\rm disk}x)^{-1/2}},\quad \nonumber \end{eqnarray} \tag{ 37 }$$

where c_disk ≡ r_vir/r_disk is the disk concentration parameter in analogy to the halo concentration parameter c_halo of Section 5.1. Equation (37) is accurate to less than 1% over the whole range r = 0–10 r_disk and it correctly converges toward the circular velocity of a point-mass potential, V^halo_c²(r) = GM^disk/r, for r → ∞.

Like in Section 3.2, we shall use Equation (16) to compute r_disk and c_disk in the DeLucia catalog. This approach is slightly inconsistent because Equation (16) was derived by Mo et al. (1998) under the assumption that the disk is supported by an isothermal halo with ρ_halo(r) ∼ r⁻², while in Section 5.1 we have assumed the more complex NFW profile. We argue, however, that Equation (16) with the empirical correction of Equation (20) is sufficiently accurate, as it successfully reproduces the observed relation between M_bary and r_disk (see Figure 2) as well as the relation between $M_{\rm H\,\hbox{\scriptsize {\sc i}}}$ and $r_{\rm H\,\hbox{\scriptsize {\sc i}}}$ (see Figure 8). It can also be shown that, for realistic values of the halo concentration c_halo (10–20 for one-galaxy systems) and the spin parameter λ (0.05–0.1, e.g., Mo et al. 1998), the scale radius of the halo r_s and the disk radius r_disk are similar. Hence, the main mass contribution of the disk comes indeed from galactocentric radii, where the halo profile is approximately isothermal, thus justifying the assumption made by Mo et al. (1998) to derive Equation (16).

Many models for the surface brightness or surface density profiles of bulges have been proposed (e.g., overview by Balcells et al. 2001). A rough consensus seems established that no single surface density profile can describe a majority of the observed bulges, but that they are generally well matched by the class of Sersic functions (Sersic 1968), Σ^bulge(r) ∼ exp [ − (r/r_bulge)^1/n], where the exponent n depends on the morphological type (Andredakis et al. 1995), such that n ≈ 4 for lenticular/early-type galaxies (de Vaucouleurs profile) and n ≈ 1 for the bulges of late-type galaxies (exponential profile). Courteau et al. (1996) find slightly steeper profiles with n = 1–2 for nearly all spirals in a sample of 326 spiral galaxies using deep optical and IR photometry, and they show that by imposing n = 1 for all late-type galaxies, the ratio between the exponential scale radius of the bulge r_bulge and the scale radius of the disk r_disk is roughly constant, $r_{\rm bulge}\approx 0.1\,\tilde{r}_{\rm disk}\approx 0.05\,r_{\rm disk}$. We shall therefore assume that all bulges have an exponential projected surface density,

$$\begin{equation} \Sigma ^{\rm bulge}(r)=\frac{M^{\rm bulge}}{2\pi \,r_{\rm bulge}^2}\,\exp \left(-\frac{r}{r_{\rm bulge}}\right) \end{equation} \tag{ 38 }$$

with r_bulge = 0.05 r_disk.

For simplicity, we assume that the bulges of all galaxies are spherical and thus described by a radial space density function ρ_bulge(r). This function is linked to the projected surface density via Σ^bulge(r) = ∫^∞_−∞dz ρ_bulge[(r² + z²)^1/2]. Numerically, we find that this model for ρ_bulge(r) is closely approximated by the Plummer model (Plummer 1911), more often used in the context of clusters,

$$\begin{equation} \rho _{\rm bulge}(r) \approx \frac{3\,M^{\rm bulge}}{4\pi \,r_{\rm Plummer}^3}\left[1+\left(\frac{r}{r_{\rm Plummer}}\right)^2\right]^{-5/2}, \end{equation} \tag{ 39 }$$

with a characteristic Plummer radius r_Plummer ≈ 1.7 r_bulge. The circular velocity profile V^bulge_c corresponding to Equation (39) is given by V^bulge_c²(r) = GM^bulge(r)/r with $M^{\rm bulge}(r)=4\pi \int _0^r{ d}\tilde{r}\,\tilde{r}^2\,\rho _{\rm bulge}(\tilde{r})$. This solves to

$$\begin{equation} {V_{c}^{\rm bulge}}^2(x) = \frac{G\!M^{\rm bulge}}{r_{\rm vir}}\times \frac{(c_{\rm bulge}x)^2c_{\rm bulge}}{[1+(c_{\rm bulge}x)^2]^{3/2}}, \end{equation} \tag{ 40 }$$

where c_bulge ≡ r_vir/r_Plummer ≈ 12 c_disk is the bulge concentration parameter.

The addition rule for gravitational potentials implies that the circular velocity profile of the combined halo–disk–bulge system in the plane of the disk is given by

$$\begin{equation} V_{c}^2(x)={V_{c}^{\rm halo}}^2(x)+{V_{c}^{\rm disk}}^2(x)+{V_{c}^{\rm bulge}}^2(x), \end{equation} \tag{ 41 }$$

where x ≡ r/r_vir as in Sections 5.1–5.3. According to Equations (34), (37), and (40), this profile is determined by six parameters: the three form parameters c_halo, c_disk, and c_bulge and the three mass scales M_vir/r_vir, M^disk/r_vir, and M^bulge/r_vir. The form parameters were calculated as explained in Sections 5.1–5.3, while the mass scales were directly adopted from the DeLucia catalog. For the satellite galaxies with no resolved halo (see Section 2), M_vir and r_vir were approximated as the corresponding quantities of the original galaxy halo just before its disappearance. An exemplar circular velocity profile for a galaxy in the DeLucia catalog at redshift z = 0 is shown in Figure 9.

Zoom In Zoom Out Reset image size

In order to evaluate the profile of a radio emission line associated with any velocity profile V_c(r), we shall first consider the line profile of a homogeneous flat ring with constant circular velocity V_c and a total luminosity of unity. If a point of the ring is labeled by the angle γ it forms with the line of sight (see Figure 10), the apparent projected velocity of that point is given by V_obs = V_c sin γ. The ensemble of all angles γ ∈ [0, 2π) therefore spans a continuum of apparent velocities V_obs ∈ (−V_c, V_c) with a luminosity density distribution $\tilde{\psi }(V_{\rm obs},V_{c})\sim { d}\gamma /{ d}V_{\rm obs}$. Imposing the normalization condition $\int { d}V_{\rm obs}\tilde{\psi }(V_{\rm obs})=1$, we find that the edge-on line profile of the ring is given by

$$\begin{equation} \tilde{\psi }(V_{\rm obs},V_{c}) = \left\lbrace\!\! \begin{array}{ll} \frac{1}{\pi \sqrt{V_{c}^2-V_{\rm obs}^2}} & {\rm if }\quad |V_{\rm obs}|<V_{c}\\ 0 & \rm otherwise. \end{array}\right. \end{equation} \tag{ 42 }$$

This profile exhibits spurious divergent singularities at |V_obs| → V_c, which, in reality, are smoothed by the random, e.g., turbulent, motion of the gas. We assume that this velocity dispersion is given by the constant $\sigma _{\rm gas}=8\rm \;km\;s^{-1}$, which is consistent with the velocity dispersions observed across the disks of several nearby galaxies (e.g., Shostak & van der Kruit 1984; Dickey et al. 1990; Burton 1971). The smoothed velocity profile is then given by

$$\begin{equation} \psi (V_{\rm obs},V_{c})\!=\!\frac{\sigma _{\rm gas}^{-1}}{\!\sqrt{2\pi }}\!\int \!\!dV\!\exp \bigg [\!\frac{(V_{\rm obs}\!\!-\!V)^2}{-2\sigma _{\rm gas}^2}\bigg ]\tilde{\psi }(V,V_{c}), \end{equation} \tag{ 43 }$$

which conserves the normalization ∫dV_obsψ(V_obs) = 1. Some examples of the functions $\tilde{\psi }(V_{\rm obs},V_{c})$ and ψ(V_obs, V_c) are plotted in Figure 11.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

From the edge-on line profile ψ(V_obs, V_c) of a single ring and the face-on surface densities of atomic and molecular gas, $\Sigma _{\rm H\,\hbox{\scriptsize {\sc i}}}(r)$ and $\Sigma _{\rm H_2}(r)$, we can now evaluate the edge-on profiles of emission lines associated with the entire H i and H₂ disks, respectively. Since H₂ densities are most commonly inferred from CO detections, we shall hereafter refer to all molecular emission lines as "the CO line." The edge-on line profiles (or "normalized luminosity densities") $\Psi _{\rm H\,{\mathsc i}}(V_{\rm obs})$ and Ψ_CO(V_obs) are given by

$$\begin{eqnarray} \Psi _{\rm H\,{\mathsc i}}(V_{\rm obs}) & = & \frac{2\pi }{M_{\rm H\,\hbox{\scriptsize {\sc i}}}}\int _0^\infty \!\!\!\!dr\,r\,\Sigma _{\rm H\,\hbox{\scriptsize {\sc i}}}(r)\,\psi (V_{\rm obs},V_{c}(r)), \end{eqnarray} \tag{ 44 }$$

$$\begin{eqnarray} \Psi _{\rm CO}(V_{\rm obs}) & = & \frac{2\pi }{M_{{\rm H}_2}}\int _0^\infty \!\!\!\!dr\,r\,\Sigma _{\rm H_2}(r)\,\psi (V_{\rm obs},V_{c}(r)). \end{eqnarray} \tag{ 45 }$$

These two functions satisfy the normalization conditions $\int { d}V_{\rm obs}\Psi _{\rm H\,{\mathsc i}}(V_{\rm obs})=1$ and ∫dV_obsΨ_CO(V_obs) = 1. To obtain intrinsic luminosity densities, $\Psi _{\rm H\,{\mathsc i}}(V_{\rm obs})$ must be multiplied by the integrated luminosity of the H i line (see Equation (21)) and Ψ_CO(V_obs) must be multiplied by the integrated luminosity of the considered molecular emission line, e.g., the integrated luminosity of the CO(1–0) line given in Equation (23).

Figure 12 displays the line profiles $\Psi _{\rm H\,{\mathsc i}}(V_{\rm obs})$ and Ψ_CO(V_obs) for the exemplar galaxy with the velocity profile shown in Figure 9. All line profiles produced by our model are mirror symmetric, but they can, in principle, differ significantly from the basic double-horned function ψ(V_obs). Especially for CO, where the emission from the bulge can play an important role, several local maxima can sometimes be found in the line profile, in qualitative agreement with various observations (e.g., Lavezzi & Dickey 1998).

Zoom In Zoom Out Reset image size

For every galaxy in the DeLucia catalog, we computed the edge-on line profiles $\Psi _{\rm H\,\hbox{\scriptsize {\sc i}}}(V_{\rm obs})$ and Ψ_CO(V_obs), from which we extracted the line parameters indicated in Figure 12. $\Psi _{\rm H\,\hbox{\scriptsize {\sc i}}}^0\equiv \Psi _{\rm H\,\hbox{\scriptsize {\sc i}}}(0)$ and Ψ⁰_CO ≡ Ψ_CO(0) are the luminosity densities at the line center, and $\Psi _{\rm H\,\hbox{\scriptsize {\sc i}}}^{\rm max}$ and Ψ^max_CO are the peak luminosity densities, i.e., the absolute maxima of $\Psi _{\rm H\,\hbox{\scriptsize {\sc i}}}(V_{\rm obs})$ and Ψ_CO(V_obs). $w_{\rm H\,\hbox{\scriptsize {\sc i}}}^{\rm peak}$ and w^peak_CO are the line widths measured between the left and the right maximum. These values vanish if the line maxima are at the line center, such as found, for example, in slowly rotating systems. $w_{\rm H\,\hbox{\scriptsize {\sc i}}}^{50}$, w⁵⁰_CO, $w_{\rm H\,\hbox{\scriptsize {\sc i}}}^{20}$, and w²⁰_CO are the line widths measured at, respectively, the 50 percentile level or the 20 percentile level of the peak luminosity densities—the two most common definitions of observed line widths.

We shall now check the simulated line widths against observations by analyzing their relation to the mass of the galaxies. Here, we shall refer to all line width versus mass relations as TFRs, since they are generalized versions of the original relation between line widths and optical magnitudes of spiral galaxies (Tully & Fisher 1977). A variety of empirical TFRs have been published, such as the stellar mass TFR and the baryonic TFR (McGaugh et al. 2000). The latter relates the baryon mass (stars + gas) of spiral disks to their line widths (or circular velocities) and is probably the most fundamental TFR detected so far, obeying a single power law over 5 orders of magnitude in mass. We have also investigated the less fundamental empirical TFR between $M_{\rm H\,\hbox{\scriptsize {\sc i}}}$ and $w_{\rm H\,\hbox{\scriptsize {\sc i}}}^{50}$—hereafter the H i-TFR—using the spiral galaxies of the HIPASS catalog. Assuming no prior knowledge on the inclinations of the HIPASS galaxies, but taking an isotropic distribution of their axes as given, we found that the most-likely relation is

$$\begin{equation} \log \left(\frac{M_{\rm H\,\hbox{\scriptsize {\sc i}}}}{{M}_{\odot }}\right) = 2.86+2.808\times \log \left(\frac{w_{\rm H\,\hbox{\scriptsize {\sc i}}}^{50}}{\rm km\;s^{-1}}\right) \end{equation} \tag{ 46 }$$

for the Hubble parameter h = 0.73. Relative to Equation (46), the HIPASS data exhibit a Gaussian scatter with σ = 0.38 in $\log (M_{\rm H\,\hbox{\scriptsize {\sc i}}})$. Our method to find Equation (46) will be detailed in a forthcoming paper (D. O. Obreschkow et al. 2009, in preparation), especially dedicated to the H i-TFR.

Figures 13(a)–(d) show four TFRs at redshift z = 0. Each figure represents 10³ simulated galaxies (black), randomly drawn from the simulation with a probability proportional to their cold gas mass in order to include the rare objects in the high end of the MF. Spiral and elliptical galaxies are distinguished as black dots and open circles.

Zoom In Zoom Out Reset image size

Figure 13(a) shows the simulated H i-TFR together with the empirical counterpart given in Equation (46). This comparison reveals good consistency between observation and simulation for spiral galaxies. However, the elliptical galaxies lie far off the H i-TFR. In fact, simulated elliptical galaxies generally have a significant fraction of their cold hydrogen in the molecular phase, consistent with the galaxy-type dependence of the H₂/H i ratio first identified by Young & Knezek (1989). Therefore, H i is a poor mass tracer for elliptical galaxies, both in simulations and observations, leading to their offset from the TFR when only H i masses are considered. There seems to be no direct analog to the H i-TFR for elliptical galaxies.

Figures 13(b) and (c), respectively, display the simulated stellar mass TFR and the baryonic TFR, together with the observed data of McGaugh et al. (2000) corrected for h = 0.73. These data include various galaxies from dwarfs to giant spirals, whose edge-on line widths were estimated from the observed ones using the inclinations determined from the optical axis ratios. Figures 13(b) and (c) reveal a surprising consistency between simulation and observation. In Figure 13(b), both the simulated and observed data show a systematic offset from the power-law relation for all galaxies with $w_{\rm H\,\hbox{\scriptsize {\sc i}}}^{20}\lesssim 200\rm \;km\;s^{-1}$. Yet, the power-law relation is restored as soon as the cold gas mass is added to the stellar mass (Figure 13(c)), thus confirming that the TFR is indeed fundamentally a relation between mass and circular velocity.

It is interesting to consider the prediction of the simulation for the most fundamental TFR, i.e., the one between the total dynamical mass, taken as the virial mass M_vir, and the circular velocity, represented by the line width $w_{\rm H\,\hbox{\scriptsize {\sc i}}}^{20}$. This relation is shown in Figure 13(d) and reveals indeed a two to three times smaller scatter in $\log (\rm mass)$ than the baryonic TFR, hence confirming its fundamental character.

Although the simulated elliptical galaxies shown in Figures 13(b)–(d) roughly align with the respective TFRs for spiral galaxies, their scatter is larger. This is caused by the mass domination of the bulge, which leads to steep circular velocity profiles V_c(r) with a poorly defined terminal velocity. Therefore, line widths obtained by averaging over the whole elliptical galaxy are weak tracers of its spin. This picture seems to correspond to observed elliptical galaxies, where the central line widths, corresponding to the velocity dispersion in the bulge-dominated parts, are more correlated to the stellar mass than the line widths of the whole galaxy (see Faber–Jackson relation, e.g., Faber & Jackson 1976).

Kassin et al. (2007) noted that S_0.5 ≡ (0.5 V²_c + σ²_gas)^1/2 is a better kinematic estimator than the circular velocity V_c, in the sense that it markedly reduces the scatter in the stellar mass TFR. However, since our model assumes a constant gas velocity dispersion σ_gas for all galaxies, it is not possible to investigate this estimator.

The predicted evolution of the four TFRs in Figures 13(a)–(d) is shown in Figures 14(a)–(d). In all four cases, the simulation predicts two important features: (1) galaxies of identical mass have broader lines (and larger circular velocities) at higher redshift, and (2) the scatter of the TFRs generally increases with redshift. The first feature is mainly a consequence of the mass–radius–velocity relation of the dark matter halos assumed in the Millennium Simulation (see Croton et al. 2006; Mo et al. 1998). This relation predicts that, given a constant halo mass, V_c scales as (1 + z)^1/2 for large z. Furthermore, the ratios $M_{\rm H\,\hbox{\scriptsize {\sc i}}}/M_{\rm vir}$ and M_stars/M_vir on average decrease with increasing redshift, explaining the stronger evolution found in Figures 14(a) and (b) relative to Figures 14(c) and (d).

Zoom In Zoom Out Reset image size

The increase of scatter in the TFRs with redshift is a consequence of the lower degree of virialization at higher redshifts, which, in the model, is accounted for via the spin parameter λ of the halos. λ is more scattered at high redshift due to the young age of the halos and the higher merger rates. More scatter in λ leads to more scatter in the radius r_disk via Equation (16) and thus to more scatter in the circular velocity V_c via Equations (37) and (40).

Current observational databases of resolved CO-line profiles are much smaller than H i databases and their signal/noise characteristics are inferior. Nevertheless efforts to check the use of CO-line widths for probing TFRs (e.g., Lavezzi & Dickey 1998) have led to the conclusion that in most spiral galaxies the CO-line widths are very similar to H i-line widths, even though the actual line profiles may radically differ. Figure 15 shows our simulated relation between $w_{\rm H\,\hbox{\scriptsize {\sc i}}}^{20}$ and w²⁰_CO, as well as the linear fit to observations of 44 galaxies in different clusters (Lavezzi & Dickey 1998). These observations are indeed consistent with the simulation. The simulated elliptical galaxies tend to have higher $w_{\rm CO}^{20}/w_{\rm H\,\hbox{\scriptsize {\sc i}}}^{20}$ ratios than the spiral ones, due to fast circular velocity of the bulge component implied by its own mass.

Zoom In Zoom Out Reset image size

The line profiles and widths studied in this section correspond to galaxies observed edge-on. First-order corrections for spiral galaxies seen at an inclination i ≠ 90° can be obtained by dividing the normalized luminosity densities $\Psi _{\rm H\,\hbox{\scriptsize {\sc i}}}^0$, $\Psi _{\rm H\,\hbox{\scriptsize {\sc i}}}^{\rm max}$, Ψ⁰_CO, and Ψ^max_CO by sin  i, while multiplying the line widths $w_{\rm H\,\hbox{\scriptsize {\sc i}}}^{50}$, $w_{\rm H\,\hbox{\scriptsize {\sc i}}}^{20}$, $w_{\rm H\,\hbox{\scriptsize {\sc i}}}^{\rm peak}$, w⁵⁰_CO, w²⁰_CO, and w^peak_CO by sin  i. More elaborate corrections, accounting for the isotropy of the velocity dispersion σ_gas, were given by Makarov et al. (1997) and Kannappan et al. (2002).

The consistency question arises, because the DeLucia catalog relies on a simplified version of a Schmidt–Kennicutt law (Schmidt 1959; Kennicutt 1998), i.e., a prescription where the SFR scales as some power of the surface density of the ISM. However, in a smaller-scaled picture, new stars are bred inside molecular clouds, and hence it must be verified whether our prescription to assign H₂ masses to galaxies is compatible with the macroscopic Schmidt–Kennicutt law. Our prescription exploited the empirical power law between the pressure of the ISM and its molecular content, as first presented by Blitz & Rosolowsky (2004, 2006). Based on this power-law relation, Blitz & Rosolowsky (2006) themselves formulated an alternative model for the computation of SFRs in galaxies, which seems more fundamental than the Schmidt–Kennicutt law. Applying both models for star formation to six molecule-rich galaxies in the local universe, they showed that their new pressure-based law predicts SFRs nearly identical to the ones predicted by the Schmidt–Kennicutt law. Therefore, our choice to divide cold hydrogen in H i and H₂ according to pressure is indeed consistent with the prescription for SFRs used by De Lucia & Blaizot (2007) and Croton et al. (2006).

A first limitation of our simulation comes from the assumption that the surface densities of H i and H₂ are axially symmetric (no spiral structures, no central bars, no warps, no satellite structures). In general, our models describe all galaxies as regular galaxies—as do all semianalytic models for the Millennium Simulation. Hence, the simulation results cannot be used to predict the H i and H₂ properties of irregular galaxies.

While our models allowed us to reproduce the observed relation between $M_{\rm H\,\hbox{\scriptsize {\sc i}}}$ and $r_{\rm H\,\hbox{\scriptsize {\sc i}}}$ remarkably well for various spiral galaxies (e.g., Figure 8), it tends to underestimate the size of H i distributions in elliptical galaxies. For example, observations by Morganti et al. (1997) show that seven nearby E- and S0-type galaxies all have very complex H i distributions, often reaching far beyond the corresponding radius of a mass-equivalent disk galaxies. The patchy H i distributions found around elliptical galaxies are probably due to mergers and tidal interactions, which could not be modeled in any of the semianalytic schemes for the Millennium Simulation.

Another limitation arises from neglecting stellar bulges as an additional source of disk pressure in Equation (7) (Elmegreen 1989). Especially, the heavier bulges of early-type spiral galaxies could introduce a positive correction of the central pressure and hence increase the molecular fraction, thus leading to very sharp H₂ peaks in the galaxy centers, such as observed, for example, in the SBb-type spiral galaxy NGC 3351 (Leroy et al. 2008). Our model for the H₂ surface density of Equation (12) fails at predicting such sharp peaks, although the predicted total H i and H₂ mass and the corresponding radii and line profiles are not significantly affected by this effect.

Additional limitations are likely to occur at higher redshifts, where our models make a number of assumptions based on low-redshift observations. Furthermore, the underlying DeLucia catalog itself may suffer from inaccuracies at high redshift, but we shall restrict this discussion to possible issues associated with the models in this paper.

Regarding the subdivision of hydrogen into atomic and molecular material (Section 3), our most critical assumption is the treatment of all galactic disks as regular exponential structures in hydro-gravitational equilibrium. This model is very likely to deviate more from the reality at high redshift, where galaxies were generally less virialized and mergers were much more frequent (de Ravel et al. 2008). Less virialized disks are thicker, which would decrease the average pressure and fraction of molecules compared to our model. Yet, galaxy mergers counteract this effect by creating complex shapes with locally increased pressures, where H₂ forms more efficiently, giving rise to merger-driven starbursts. Therefore, it is unclear whether the assumption of regular disks tends to underestimate or overestimate the H₂/H i ratios.

Another critical assumption is the high-redshift validity of the local relation between the H₂/H i ratio R_mol and the external gas pressure P (Equation 6). This relation is not a fundamental thermodynamic relation, but represents the effective relation between the average H₂/H i ratio and P, resulting from complex physical processes such as cloud formation, H₂ formation on metallic grains, and H₂ destruction by the photodissociative radiation field of stars and supernovae. Therefore, the R_mol–P relation could be subjected to a cosmic evolution resulting from the cosmic evolution of the cold gas metallicity or the cosmic evolution of the photodissociative radiation field. However, the metallicity evolution is likely to be problem only at the highest redshifts (z ≳ 10). Observations in the local universe show that spiral galaxies with metallicities differing by a factor 5 fall on the same R_mol–P relation (Blitz & Rosolowsky 2006). Yet, the average cold gas metallicity of the galaxies in the DeLucia catalog is only a factor 1.9 (3.6) smaller at z = 5 (z = 10) than in the local universe. These predictions are consistent with observational evidence from the SDSS that stellar metallicities were at most a factor 1.5–2 smaller at z ≈ 3 than today (Panter et al. 2008). The effect of the cosmic evolution of the photodissociative radiation field on the R_mol–P relation is difficult to assess. Blitz & Rosolowsky (2006) argued that the ISM pressure and the radiation field both correlate with the surface density of stars and gas, and therefore the radiation field is correlated to pressure. This is supported by observations in the local universe showing that the R_mol–P relation holds true for dwarf galaxies and spiral galaxies spanning almost 3 orders of magnitude in SFR. For those reasons, the R_mol–P relation could indeed extend surprisingly well to high redshifts.

In the expression for the disk pressure in Equation (7) (Elmegreen 1989), we assumed a constant average velocity dispersion ratio 〈f_σ〉. Observations suggest that V_c/σ_gas decreases significantly with redshift (Förster Schreiber et al. 2006; Genzel et al. 2008), and therefore 〈f_σ〉 perhaps increases. This would lead to even higher H₂/H i ratios than predicted by our model. However, according to Equation (10) this is likely to be a problem only for galaxies with M_stars>M_gas, while most galaxies in the simulation at z>2 are indeed gas dominated.

Regarding cold gas geometries and velocity profiles, the most important limitation of our model again arises from the simplistic treatment of galactic disks as virialized exponential structures. Very young galaxies (≲$10^8\,\rm yr$) or galaxies undergoing a merger do not conform with this model, and therefore the predicted velocity profiles may be unreal and the disk radii may be meaningless. This is not just a limitation of the simulation, but it reveals a principal difficulty to describe galaxy populations dominated by very young or merging objects with quantities such as $r_{\rm H\,\hbox{\scriptsize {\sc i}}}$ or $w_{\rm H\,\hbox{\scriptsize {\sc i}}}^{50}$, which are common and useful for isolated systems in the local universe.

The radio line widths predicted by our model (Section 5.4) may be underestimated at high redshift, due to the assumption of a constant random velocity dispersion σ_gas. Förster Schreiber et al. (2006) found V_c/σ_gas ≈ 2–4 for 14 UV-selected galaxies at redshift z ≈ 2. This result suggests that radio lines at z ≈ 2 should be about 20%–30% broader than predicted by our model, and therefore the evolution of the TFRs could be slightly stronger than shown in Figure 14.

In summary, the H i and H₂ properties predicted for galaxies at high redshift are generally uncertain, even though no significant, systematic trend of the model errors could be identified. Perhaps the largest deviations from the real universe occur for very young galaxies or merging objects, while isolated field galaxies, typically late-type spiral systems, might be well described by the model at all redshifts.

In Section 3.4, we ascribed CO(1–0) luminosities to the H₂ masses using the metallicity-dependent X-factor of Equation (26). This model neglects several important aspects: (1) the projected overlap of molecular clouds, which is negligible in the local universe, may become significant at high redshifts, where galaxies are denser and richer in molecules; (2) the temperature of the cosmic microwave background (CMB) increases with redshift, hence changing the level population of the CO molecule (Silk & Spaans 1997; Combes et al. 1999); (3) the CMB presents a background against which CO sources are detected; (4) the molecular material in the very dense galaxies, such as ultraluminous infrared galaxies, may be distributed smoothly rather than in clouds and clumps (Downes et al. 1993); (5) the higher SFRs in early galaxies probably led to higher gas temperatures, hence changing the CO-level population.⁷ Combes et al. (1999) presented a simplistic simulation of the cosmic evolution of X, taking the cosmic evolution of metallicity and points (1) and (2) into account. They found that for an H₂-rich disk galaxy 〈X〉 increases by a factor 1.8 from redshift z = 0 to z = 5. This value approximately matches the average increase of X by a factor 2 predicted by our simulation using the purely metallicity-based model of Equation (26). This indicates that the effects of (1) and (2) approximately balance each other. If a more elaborate model for X becomes available, the latter can be directly applied to correct our CO predictions. In fact, the X-factor only affects the CO(1–0) luminosity L_CO calculated via Equation (23), but has no effect on the line properties considered in this paper, namely, the line widths w^peak_CO, w⁵⁰_CO, and w²⁰_CO and the normalized luminosity densities Ψ⁰_CO and Ψ^max_CO.