A VIRTUAL SKY WITH EXTRAGALACTIC H i AND CO LINES FOR THE SQUARE KILOMETRE ARRAY AND THE ATACAMA LARGE MILLIMETER/SUBMILLIMETER ARRAY^*

Here, we recapitulate the galaxy simulation presented in earlier studies. This simulation relies on three consecutive layers: (1) a simulation of the cosmic evolution of dark matter (Springel et al. 2005); (2) a semi-analytic simulation of the evolution of galaxies on the dark matter skeleton (Croton et al. 2006; De Lucia & Blaizot 2007); and (3) a post-processing to split the cold hydrogen masses associated with each galaxy into H i and H₂ (Obreschkow et al. 2009a).

For the dark matter simulation, we adopted the Millennium Simulation (Springel et al. 2005), an N-body dark matter simulation within the standard Λ cold dark matter (ΛCDM) cosmology. This simulation uses a cubic simulation box with periodic boundary conditions and a comoving volume of (500 h⁻¹ Mpc)³. The Hubble constant was fixed to H₀ = 100 h km s^-1 Mpc⁻¹ with h = 0.73. The other cosmological parameters were chosen as Ω_matter = 0.25, Ω_baryon = 0.045, Ω_Λ = 0.75, and σ₈ = 0.9. The simulation box contains ∼10¹⁰ particles with individual masses of 8.6 × 10⁸ M_☉. This mass resolution allows the identification of structures as low in mass as the Small Magellanic Cloud (see Springel et al. 2005).

For the second simulation layer, i.e., the cosmic evolution of the galaxies distributed on the dark matter skeleton, we adopt the semi-analytic model of De Lucia & Blaizot (2007) (see also Croton et al. 2006). In this macroscopic model, all galaxies are represented by a list of global properties, such as position, velocity, and total masses of gas, stars, and black holes. These properties are evolved using empirically or theoretically motivated formulae for mechanisms, such as gas cooling, re-ionization, star formation, gas heating by supernovae, starbursts, black hole accretion, black hole coalescence, and the formation of stellar bulges via disk instabilities. The resulting virtual galaxy catalog (hereafter the "DeLucia-catalog") contains the positions, velocities, merger histories, and intrinsic properties of ∼3 × 10⁷ galaxies at 64 cosmic time steps. The free parameters in the semi-analytic model were tuned to various observations in the local universe (see Croton et al. 2006). Therefore, despite the simplistic implementation and the possible incompleteness of this model, the DeLucia-catalog nonetheless provides a good fit to the joint luminosity/color/morphology distribution of observed low-redshift galaxies (Cole et al. 2001; Huang et al. 2003; Norberg et al. 2002), the bulge-to-black hole mass relation (Häring & Rix 2004), the Tully–Fisher relation (Giovanelli et al. 1997), and the cold gas metallicity as a function of stellar mass (Tremonti et al. 2004).

In this paper, we are particularly interested in the cold gas masses of the galaxies in the DeLucia-catalog. These cold gas masses are the net result of (1) gas accretion by cooling from a hot halo (dominant mode) and galaxy mergers, (2) gas losses by star formation and feedback from supernovae, and (3) cooling flow suppression by feedback from accreting black holes. The DeLucia-catalog does not distinguish between molecular and atomic cold gas, but simplistically treats all cold gas as a single phase. The atomic and molecular phases are therefore dealt with in the third simulation layer.

The third simulation layer, i.e., the subdivision of the cold hydrogen mass of each galaxy into H i and H₂ distributions (Obreschkow et al. 2009a), relies on an analytic model for the mass distributions of H i and H₂ within regular galaxies. In this model, the column densities of H i and H₂, Σ_H i and $\Sigma _{\rm H_2}$, respectively, are given by

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

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

where r denotes the galactocentric radius, r_disk is a scale length, R^c_mol is the H₂/H i mass ratio at the galaxy center, and $\tilde{\Sigma }_{\rm H}$ is a normalization factor. We derived Equations (1) and (2) based on a list of empirically supported assumptions, the most important of which are: (1) the cold gas of regular galaxies resides in a flat disk (see Leroy et al. 2008 for local spiral galaxies, Young 2002 for local elliptical galaxies, Tacconi et al. 2006 for galaxies at higher redshifts); (2) the surface density of the total hydrogen component (H i+H₂) is well described by an axially symmetric exponential profile (Leroy et al. 2008); (3) the local H₂/H i mass ratio scales as a power of the gas pressure of the ISM outside molecular clouds (Blitz & Rosolowsky 2006).

Using Equations (1) and (2), we can characterize the H i- and H₂-content of every simulated galaxy in the DeLucia-catalog. The resulting hydrogen simulation successfully reproduces many local observations of H i and H₂, such as mass functions (MFs), mass–diameter relations, and mass–velocity relations (Obreschkow et al. 2009a). This success is quite surprising, since our model for H i and H₂ only introduced one additional free parameter to match the observed average space density of cold gas in the local universe (Obreschkow et al. 2009a). A key prediction of this simulation is that the H₂/H i ratio of most regular galaxies increases dramatically with redshift, hence causing a clear signature of cosmic downsizing in the H₂-MF (Obreschkow & Rawlings 2009b).

Despite its consistency with existing observations, we emphasize that the presented model for the cosmic evolution of H i and H₂ is simplistic and uncertain. In particular, at high z, the model assumptions may significantly differ from the reality. For example, high-z galaxies are likely to be more disturbed due to higher merger rates and long dynamic timescales compared to their age. There is also evidence that cold gas disks become more turbulent with redshift (e.g., Förster Schreiber et al. 2006; Genzel et al. 2008). Uncertainties from these and other model limitations are discussed briefly in Section 4 and in depth in Section 6 of Obreschkow et al. (2009a).

We shall now describe how the cubic simulation box is transformed into a virtual sky field. This procedure can be regarded as a fourth simulation layer on top of the hierarchical simulation described in Section 2.1.

The method adopted here closely follows the one described by Blaizot et al. (2005), namely the building of a chain of replicated simulation boxes along the line of sight, as shown in Figure 1. At any position in this chain, the galaxies are drawn from the cosmic time in the simulation, which corresponds to the look-back time, at which the galaxy is seen by the observer $\mathcal {O}$. Since our galaxy simulation uses 64 discrete time steps, we describe each galaxy in the cone by its properties at the closest available time step,³ in terms of redshift. This defines the spherical shells of identical cosmic time, which are separated by dashed lines in Figure 1. The relatively narrow redshift separation of these shells ensures that the assigned galaxy properties cannot differ significantly from the properties at their exact look-back time.

Zoom In Zoom Out Reset image size

The same galaxy can appear once in every box in Figure 1, but with different intrinsic properties due to the cosmic evolution. However, the position of the repeated galaxy in comoving coordinates will be very similar, which can result in spurious radial features for the observer $\mathcal {O}$ (see Figure 1 in Blaizot et al. 2005). To suppress this effect, we randomize the galaxy positions by applying random symmetry operations to each box in the chain. These operations consist of 90° rotations, inversions, and continuous translations.⁴ Applying these symmetry operations also removes the non-physical periodicity of 500 h⁻¹ Mpc associated with the side length of the periodic simulation box. But we emphasize that applying the symmetry operations does not provide information on scales larger than the simulation box. Symmetry operations can, however, introduce unwanted small-scale density variations at the interface of two neighboring boxes. These and other limitations of this method are discussed by Blaizot et al. (2005).

From the randomized chain of replicated simulation boxes, an observing cone can be extracted (shaded region in Figure 1). Each galaxy in this cone is projected onto the celestial sphere centered about the vernal point (R.A. = 0, decl. = 0). The Euclidian projection formulas for arbitrary large angles are

$$\begin{eqnarray} {\rm R.A.} & = & \arctan \left(\frac{r_{x}}{r_{z}}\right), \end{eqnarray}\vspace*{-12pt} \tag{ 3 }$$

$$\begin{eqnarray} {\rm decl.} & = & \arctan \left(\frac{r_{y}}{\sqrt{r_{x}^2+r_{z}^2}}\right), \end{eqnarray} \tag{ 4 }$$

where r_x, r_y, and r_z are the comoving coordinates of the galaxy relative to the observer (see Figure 1). The "cosmological redshift" z of each galaxy is computed directly from its comoving distance D_C = (r²_x + r²_y + r²_z)^1/2, while the Doppler-shift-corrected "apparent redshift" is computed as $\tilde{z}=z\,{+}V_{r}/c$, where V_r is the peculiar recession velocity of the galaxy relative to the Hubble flow.

Figure 1 shows that the opening angle φ of the virtual sky field is set by the maximal comoving distance D_C,max via

$$\begin{equation} \varphi = 2\,\arcsin \frac{s_{\rm box}}{2\,D_{\rm C,max}}, \end{equation} \tag{ 5 }$$

where s_box is the comoving side length of the simulation box. Given a value of s_box and a choice of cosmological parameters, Equation (5) implies a one-to-one relation between φ and the maximal redshift z_max.

Figure 2 displays the relation between φ, D_C,max, and z_max for the cosmological parameters of the Millennium Simulation (Section 2.1) and three different side lengths s_box. The choice s_box = 500 h⁻¹ Mpc (solid line) corresponds to the box of the Millennium Simulation. s_box = 62.5 h⁻¹ Mpc (dashed line) corresponds to the "Milli-Millennium" Simulation, a small test version of the Millennium Simulation. s_box = 2 h⁻¹ Gpc (dash-dotted line) corresponds to the giant simulation box of the Horizon-4π Simulation, a dark matter stimulation with ∼10 times less mass resolution than the Millennium Simulation (Prunet et al. 2008; Teyssier et al. 2009; see also Section 4.4).

Zoom In Zoom Out Reset image size

We shall now assign apparent line fluxes to each galaxy in the mock observing cone constructed in Section 2.2. In general, the frequency-integrated line flux S of any emission line can be computed from the frequency-integrated luminosity (= intrinsic power) via

$$\begin{equation} S= \frac{L}{4\pi D_{\rm L}^2}, \end{equation} \tag{ 6 }$$

where D_L = (1 + z) D_C is the luminosity distance to the source. We note that Equation (6) takes a different form for velocity-integrated fluxes and brightness temperatures. For reference, a summary of all relations between frequency-integrated and velocity-integrated fluxes, luminosities, and brightness temperatures has been compiled in Appendix A of Obreschkow et al. (2009b).

In Sections 2.3.1 and 2.3.2, we shall summarize the models used to estimate the frequency-integrated luminosities L (hereafter simply "luminosities") of different emission lines.

2.3.1. Conversion of H i Mass into H i-line Luminosities

We evaluate the H i luminosities L_H i from the H i masses M_H i of the simulated galaxies via the standard conversion (e.g., Meyer et al. 2004)

$$\begin{equation} \frac{L_{{\rm H}\,\hbox{\scriptsize {\sc i}}}}{{L}_{\odot }}=6.27\times 10^{-9}\cdot \frac{M_{{\rm H}\,\hbox{\scriptsize {\sc i}}}}{{M}_{\odot }}. \end{equation} \tag{ 7 }$$

The H i line or "21 cm line" at a rest-frame frequency of 1.420 GHz originates from the photon-mediated transition between the two spin states of the proton–electron system in the electronic ground state. The upper spin state has a low spontaneous decay rate of f = 2.9 × 10⁻¹⁵ Hz. This frequency is about five orders of magnitude smaller than that of H i–H i collisions (Binney & Merrifield 1998). Hence, the two spin states are in thermal equilibrium with the kinetic state of the gas, which implies a spin-temperature far above the spin excitation temperature T_ex ≈ 0.07 K. Therefore, the spin systems are in the high-temperature limit, where 3/4 of all systems are in the upper (three-fold degenerate) state. The radiative power emitted per atom can therefore be calculated as L_H i = 0.75 f h_p 1.4 GHz (h_p is the Planck constant), which readily reduces to Equation (7).

We have neglected H i self-absorption, since this seems to affect only massive spiral galaxies when observed almost edge-on (Ferrière 2001; Wall 2006). This assumption should also be valid for high-redshift galaxies, since their H i masses were not much larger than today, as can be inferred from Lyα-absorption measurements against distant quasars (Lah et al. 2007) and as is predicted by our simulation (Obreschkow & Rawlings 2009b).

H i in collapsed structures, i.e., in galaxies, is generally warm enough (≳50 K) that the Cosmic Microwave Background (CMB) can be safely neglected as an observing background for all galaxies at z < 10. Only in the intergalactic medium (IGM) during the cosmic Epoch of Re-ionization (EoR) can H i appear in 21 cm absorption against the CMB (e.g., Iliev et al. 2002).

2.3.2. Conversion of H₂ Mass into CO-line Luminosities

Having assigned an integrated line flux to each galaxy in the simulation, we can now refine their attributes by characterizing each line with a profile. To this end, we depart from the edge-on line profiles evaluated in Obreschkow et al. (2009a) for each galaxy in the simulation. We represented those profiles by normalized flux densities Ψ(V), where V is the velocity measured in the rest frame of the center of the observed galaxy. The normalization condition, ∫_VΨ(V)dV = 1, implies that Ψ(V) only needs to be multiplied by the velocity-integrated flux (in units of Jy km s⁻¹) in order to obtain actual flux densities (in units of Jy). For each galaxy, we calculated two profiles Ψ(V), one for the H i component and one for the H₂ component (associated with CO).⁵ This calculation relied on a detailed mass model based on the halo, disk, and bulge of the galaxies, combined with our model for the H i- and H₂-surface densities given in Equations (1) and (2). For practical reasons, the resulting line profiles were reduced to five parameters (see Figure 8): the normalized flux density at the line center Ψ₀; the normalized peak flux density Ψ_max, usually corresponding to the two peaks of a double-horn profile; the line width w_peak between the two peaks of the double-horn profile; the line width w₅₀ at the 50 percentile level of the peak flux density; and the line width w₂₀ at the 20 percentile level. The original normalized line profile can be approximately recovered from these parameters using the formulas in Appendix A.

The remaining task consists of correcting the line profiles for the inclination i of each galaxy.⁶ i is defined as the angle between the line of sight and the galaxy's rotation axis; hence i = 0° corresponds to face-on galaxies and i = 90° corresponds to edge-on galaxies. In the absence of a random gas velocity dispersion, apparent line widths w^obs could be computed from the edge-on line widths w, via w^obs = w · sin i, while apparent normalized flux densities would scale as Ψ^obs = Ψ/sin i.

Here, we assume that the cold gas has a random, isotropic velocity dispersion characterized by a Gaussian velocity distribution in each space dimension. The observed line profile of a face-on galaxy (i = 0°) then takes the shape of a Gaussian function. Under no inclination can the line profile become more narrow than this Gaussian function, nor can the normalized line flux densities become higher than the peak of this Gaussian function. Let σ be the standard deviation of the Gaussian velocity dispersion. Then the minimum line widths are given by $w_{20}^{\rm obs}=2\sqrt{-2\,\ln (0.2)}\approx 3.6\,\sigma$, $w_{50}^{\rm obs}=2\sqrt{-2\,\ln (0.5)}\approx 2.4\,\sigma$, and w^obs_peak = 0, and the maximum normalized flux densities are $\Psi _0=(\sigma \sqrt{2\pi })^{-1}$ and $\Psi _{\rm max}=(\sigma \sqrt{2\pi })^{-1}$.

In addition, the maximal normalized flux density Ψ_max cannot differ from the central normalized flux density Ψ₀ by an arbitrarily large amount due to the line profile smoothing. Explicitly, the slope in the emission line between the points Ψ₀ and Ψ_max cannot exceed the maximal slope of the Gaussian velocity function, which is equal to 0.24 σ⁻². This requirement translates into an upper bound for Ψ_max equal to Ψ^obs₀ + 0.12 σ⁻² w^obs_peak. A set of equations respecting all of these conditions is given by

$$\begin{eqnarray} w_{20}^{\rm obs} & = & (w_{20}-3.6\,\sigma)\cdot \sin {i}+3.6\,\sigma, \end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray} w_{50}^{\rm obs} & = & (w_{50}-2.4\,\sigma)\cdot \sin {i}+2.4\,\sigma, \end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray} w_{\rm peak}^{\rm obs} & = & w_{\rm peak}\cdot \sin {i}, \end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray} \Psi _{0}^{\rm obs} & = & \min \!\left(\!\frac{\Psi _{0}}{\sin {i}},\frac{1}{\sigma \sqrt{2\pi }}\right), \end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray} \Psi _{\rm max}^{\rm obs} & = & \min \!\left(\frac{\Psi _{\rm max}}{\sin {i}},\frac{1}{\sigma \sqrt{2\pi }},\Psi _{0}^{\rm obs}\!+\!\frac{0.12w_{\rm peak}^{\rm obs}}{\sigma ^2}\!\right)\!. \end{eqnarray} \tag{ 12 }$$

For all the line profiles, we here adopt σ = 8 km s⁻¹ to remain consistent with Obreschkow et al. (2009a). We note, however, that high-redshift galaxies may have higher velocity dispersions (Förster Schreiber et al. 2006) perhaps due to an intense ongoing accretion of gas. Another limitation of the presented line model is that all CO lines have by definition the same line shape. This assumption nevertheless approximately agrees with simultaneous observations of different emission lines (e.g., Weiss et al. 2007; Greve et al. 2009).

We shall finalize our simulation of line-emitting galaxies by ascribing an angular distribution (flux per unit solid angle) in the sky to each emission line. To this end, we assume that, for a fixed galaxy, the angular distributions of the H i flux and the different CO fluxes are proportional to the angular distributions of H i and H₂, respectively. We emphasize, however, that the conversion factor between H₂ masses and CO fluxes depend on the J level of the CO transition and on the galaxy, such as outlined in Section 2.3.2.

Given this assumption, we can infer the line flux densities per unit solid angle from the surface densities Σ_H i(r) and $\Sigma _{\rm H_2}(r)$ given in Equations (1) and (2). We only need to normalize these densities to the respective line fluxes and to replace the scale radius r_disk by the apparent scale radius r_disk/D_A, where D_A = (1 + z)⁻¹ D_C is the angular diameter distance.

The shapes of the line-emitting regions need to be corrected for the inclinations i. If q_HI,0 and $q_{\rm H_2,0}$ respectively denote the intrinsic axis ratios of the atomic and molecular gas in galaxies, then the apparent axis ratios are given by (Equation (1) in Kannappan et al. 2002)

$$\begin{eqnarray} q_{{\rm H}\,\hbox{\scriptsize {\sc i}}}^2 & = & \cos ^2 i+q_{\rm HI,0}^2\,\sin ^2 i, \end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray} q_{\rm H_2}^2 & = & \cos ^2 i+q_{\rm H_2,0}^2\,\sin ^2 i. \end{eqnarray} \tag{ 14 }$$

These relations satisfy $q_{{\rm H}\,\hbox{\scriptsize {\sc i}}}=q_{\rm H_2}=1$, if i = 0° (face-on), and q_H i = q_HI,0 and $q_{\rm H_2}=q_{\rm H_2,0}$, if i = 90° (edge-on).

We assume that all galaxies have the same values for respectively q_HI,0 and $q_{\rm H_2,0}$. In local spiral galaxies, we typically find q_HI,0 = 0.1, as can be seen from high-resolution maps of the edge-on spiral galaxies NGC 891 and NGC 4565 (Rupen 1991). To our knowledge, no reliable estimate of $q_{\rm H_2,0}$ for disk galaxies is available. However, simultaneous CO(1–0) and optical observations revealed that the density of stars in nearby galaxies correlates with the density of molecular gas, when averaged over sufficiently large (> kcp) scales (Richmond & Knapp 1986; Leroy et al. 2008), probably as a natural consequence of the formation of stars from molecular gas. Therefore, we shall assume that the intrinsic aspect ratio of molecular gas $q_{\rm H_2,0}$ is identical to that of stellar disks. The latter is ∼0.1, as can be seen in the sample of 34 nearby edge-on spiral galaxies studied by Kregel et al. (2002). We therefore adopt $q_{\rm H_2,0}=0.1$. We stress that $q_{\rm HI,0}=q_{\rm H_2,0}$ does not contradict the fact that characteristic scale radii and scale heights of H i distributions are generally larger than those of H₂ distributions (Leroy et al. 2008).

The assumption of constant values for q_HI,0 and $q_{\rm H_2,0}$ for all galaxies at all redshifts may not be verified at high redshifts due to timescale arguments. In fact, a significant fraction of the galaxies at z > 5 may have an age comparable to their dynamical timescale. Their cold gas distribution might therefore be thicker and less circularly symmetric than the flat gas disks seen today. However, no reliable estimates of q_HI,0 and $q_{\rm H_2,0}$ beyond the local universe are available today.

By successively applying the simulation steps described in Section 2, we have constructed an observing cone of line-emitting galaxies. Figure 3 shows a longitudinal slice of this cone with a thickness of 10 Mpc. This slice corresponds to a diagonal cut, as illustrated in Figure 3 and has an opening angle of 58. Each pixel inside this slice corresponds to a galaxy. The structure of the cosmic web appears clearly, as well as the increasing filamentarity of this structure with decreasing redshift. The color scales represent the H₂/H i mass ratios of the galaxies. We can clearly recognize the pressure-driven cosmic decline of this ratio (see Obreschkow & Rawlings 2009b).

Zoom In Zoom Out Reset image size

The mock observing cone translates into a virtual sky field when projected onto a sphere using Equations (3) and (4). Figure 4 displays the H i- and CO-flux densities of galaxies between z = 1 and z = 1.1 in a small extract of this virtual sky. The surface densities of the galaxies have been modeled using Equations (1) and (2) in the manner described in Section 2.5. The more massive galaxies in the field reveal ring-like H i distributions with CO-rich central regions. By contrast, some of the smaller galaxies with low surface brightness (LSB), have most of their H i in the center, with nearly no detectable CO. In general, CO (and hence H₂) is more compact than H i. All these simulated features compare well to observed radial distributions (averaged on suitably large scales) of H i and CO in nearby galaxies (Leroy et al. 2008) as demonstrated in Obreschkow et al. (2009a) and Obreschkow & Rawlings (2009a). Nevertheless, we emphasize that the idealized nature (e.g., axial symmetry) of the simulated galaxies misses out a list of structural features found in many real galaxies, such as spiral arms, warps, and lopsidedness.

Zoom In Zoom Out Reset image size

The coloring of the CO-surface densities in Figure 4 (right) represents the CO(5–4)/CO(1–0) line ratio in terms of brightness temperatures. For normal galaxies, without a particular source of heating, this ratio is much smaller than unity (yellow coloring); however, for some galaxies with strong heating by an ongoing SB or AGN, the higher order lines can get excited (red coloring). These mechanisms and our model to assess them are discussed in Obreschkow et al. (2009b).

Figure 9 in Appendix C shows a 3 times larger sky field than Figure 4 at the three redshifts z ≈ 1, z ≈ 3, and z ≈ 6. The progression from z ≈ 1 to z ≈ 6 in Figure 9 reveals two notable features. Firstly, galaxy sizes decrease with redshift. In fact, the angular diameter distances at z ≈ 1 and z ≈ 3 are virtually identical according to the cosmology adopted in this paper (Section 2.1). Therefore, the galaxy sizes of these two virtual sky maps can be compared directly. The angular diameter distance at z ≈ 6 is about 25% smaller, hence the same physical scales appear slightly oversized. The size evolution of the galaxies reflects the cosmic evolution of the volume/mass ratio of the dark matter halos (Gunn & Gott 1972). We discussed the impact of this evolution on the surface densities of H i and H₂ in Obreschkow & Rawlings (2009a).

Secondly, the CO(5–4)/CO(1–0) line ratios of CO-detectable galaxies increases with redshift. In fact, at z ≈ 6 no galaxy with a line ratio significantly below unity (i.e., with yellow coloring) can be seen. This feature relies partially on the compactness of the galaxies, which, according to our model for CO lines (Obreschkow et al. 2009b), allows an efficient heating by star formation. An additional reason for the absence of low CO(5–4)/CO(1–0) line ratios at z ≳ 6 is that molecular gas in galaxies with no significant star formation and no AGN will be hardly detectable in CO due to its near thermal equilibrium with the CMB.

In the simulated observing cone, we can readily count the number of galaxies per redshift interval with line fluxes above a certain threshold. This dN/dz analysis is a key step toward the prediction of the number of sources detectable with any particular telescope and survey strategy. In this section, we focus on the number of sources detected in a peak-flux-density-limited survey, and we restrict the presented results to the H i-, CO(1–0)-, CO(5–4)-, and CO(10–9)-emission lines. Results for other CO-emission lines and/or for integrated flux limited surveys are presented in Appendix D.

Figure 5 shows the dN/dz functions for six different peak flux density limits, logarithmically spaced between 1 mJy and 10 nJy. Peak flux densities for each source and emission line are calculated in the way described in Section 2.4. This method accounts for the different gas distributions, rotation curves, and inclinations of the galaxies. Every source with a peak flux density above the peak flux density limit is considered as detected, while all other sources are considered as non-detected. Different aspects of Figure 5 will be discussed in detail over the following paragraphs.

Zoom In Zoom Out Reset image size

3.2.1. Cosmic Variance

The simulated dN/dz functions shown in Figure 5 (solid and dotted lines) correspond to a bin size of Δz = 0.1 and a sky field of 4 × 4 deg², extracted from one particular realization of the mock observing cone, that is one random choice of symmetry operations for the replicated simulation boxes (see Section 2.2). The wiggles visible in the simulated dN/dz functions are physical. Similar wiggles can be expected for a real sky survey of a sky field of 4 × 4 deg² with a redshift bin size of Δz = 0.1. The fact that the amplitude of those wiggles does not decrease as $1/\sqrt{{d}N/{d}z}$ clearly uncovers the presence of the cosmic large-scale structure, also visible in Figure 3.

To quantify the effects of cosmic variance, we now consider the dN/dz functions extracted from five different random realizations of the mock observing cone. Figure 6 shows the corresponding dN/dz functions for a peak flux limited H i survey with a peak flux limit of 1 μ Jy. Each function uses a bin size of Δz = 0.1 and a small sky field of 1 × 1 deg² in order to make the effects of cosmic structure obvious. As a rough estimate the log-scatter between the different dN/dz functions is about 0.1 dex. From this small scatter, we conclude that cosmic variance is, in most cases, negligible compared to the uncertainties of the semi-analytic galaxy model.

Zoom In Zoom Out Reset image size

However, the comoving volume per unit solid angle and unit redshift varies as a function of redshift. Therefore, the scatter due to cosmic variance varies with redshift. It is largest at the lowest redshifts (z < 0.5), where the comoving surface per unit solid angle is small, and at the highest redshifts (z > 5), where the radial comoving distance per unit of redshift is small. In these redshift regimes, effects of cosmic variance should therefore be estimated, when comparing simulated data to observations.

3.2.2. Completeness

3.2.3. Parameterization of the dN/dz plots

The simulated dN/dz functions can easily be recovered from the on-line database of the sky simulation (see Appendix B). Alternatively, we also approximated the dN/dz functions by analytic fits of the form

$$\begin{equation} \frac{{d}N/{d}z}{\rm deg^2} = 10^{c_1}\cdot z^{c_2}\cdot \exp (-c_3\cdot z), \end{equation} \tag{ 15 }$$

where c₁, c₂, and c₃ are free parameters. The best parameters in terms of an rms minimization are shown in Table 1 for various emission lines detected with different limits for the peak flux densities and integrated fluxes. Analytic dN/dz functions for intermediate flux limits can be approximately inferred by linearly interpolating the parameters c₁, c₂, and c₃.

Table 1. Parameters for the Analytic Fit Formula of Equation (15) for dN/dz Peak-Flux-Density-Limited and Integrated-Flux-Limited Surveys

		Limiting Peak Flux Density (Jy)						Limiting Integrated Flux (Jy km s−1)
Parameters		10−8	10−7	10−6	10−5	10−4	10−3	10−6	10−5	10−4	10−3	10−2	10−1
	c1	6.55	6.87	6.73	5.75	4.56	6.62	6.61	6.94	6.53	5.79	5.21	6.33
H i	c2	2.54	2.85	2.32	1.14	0.43	2.64	2.66	2.91	2.09	1.33	1.23	2.60
	c3	1.42	2.17	3.09	3.95	6.86	35.49	1.51	2.34	2.86	3.92	8.38	30.52
	zc	2.3	1.1	0.5	0.2	0.1	0.1	2.1	0.9	0.4	0.1	0.1	0.1
	c1	6.32	6.40	6.37	5.92	5.06	4.15	6.31	6.37	6.30	5.91	5.27	5.17
CO(1–0)	c2	2.16	2.29	2.05	1.33	0.54	0.31	2.14	2.23	1.96	1.43	0.95	1.37
	c3	1.11	1.27	1.46	1.55	1.72	3.46	1.11	1.23	1.38	1.47	1.69	5.19
	zc	6.0	5.1	2.2	0.9	0.3	0.1	6.0	5.1	2.2	0.9	0.4	0.1
	c1	6.29	6.34	6.42	6.25	5.64	4.64	6.29	6.33	6.38	6.21	5.67	4.87
CO(2–1)	c2	2.10	2.20	2.27	1.81	1.07	0.35	2.09	2.18	2.18	1.84	1.23	0.70
	c3	1.07	1.15	1.34	1.51	1.58	1.94	1.06	1.14	1.28	1.46	1.49	1.82
	zc	5.8	6.0	3.8	1.5	0.6	0.2	5.8	6.0	3.8	1.5	0.6	0.2
	c1	6.28	6.32	6.40	6.32	5.79	4.86	6.28	6.31	6.37	6.25	5.81	5.09
CO(3–2)	c2	2.08	2.16	2.27	1.95	1.22	0.47	2.08	2.14	2.21	1.89	1.36	0.87
	c3	1.06	1.12	1.27	1.44	1.48	1.58	1.06	1.11	1.23	1.36	1.43	1.56
	zc	5.8	6.0	4.9	2.0	0.7	0.3	5.8	6.0	4.7	2.0	0.8	0.3
	c1	6.28	6.31	6.37	6.26	5.74	4.84	6.28	6.30	6.34	6.19	5.77	5.03
CO(4–3)	c2	2.08	2.13	2.22	1.87	1.18	0.52	2.08	2.12	2.17	1.79	1.34	0.86
	c3	1.06	1.10	1.21	1.32	1.36	1.38	1.06	1.09	1.18	1.25	1.33	1.36
	zc	5.8	6.0	5.2	2.1	0.7	0.3	5.8	6.0	5.2	2.0	0.8	0.3
	c1	6.28	6.30	6.32	6.05	5.45	4.53	6.28	6.29	6.30	6.00	5.53	4.74
CO(5–4)	c2	2.07	2.11	2.13	1.56	0.96	0.54	2.07	2.11	2.09	1.52	1.20	0.83
	c3	1.05	1.09	1.17	1.15	1.14	1.14	1.05	1.08	1.15	1.09	1.17	1.15
	zc	5.8	6.0	5.1	1.7	0.6	0.2	5.8	6.0	4.7	1.6	0.6	0.2
	c1	6.28	6.29	6.17	5.62	4.95	4.17	6.28	6.28	6.13	5.62	5.08	4.35
CO(6–5)	c2	2.07	2.09	1.81	1.04	0.78	0.93	2.07	2.07	1.76	1.11	1.02	1.13
	c3	1.05	1.09	1.07	0.89	0.92	1.08	1.05	1.08	1.03	0.87	0.97	1.06
	zc	5.8	6.0	6.2	0.9	0.3	0.1	5.8	6.0	6.0	0.9	0.4	0.1
	c1	6.27	6.22	5.72	5.01	4.45	3.91	6.27	6.20	5.70	5.08	4.55	4.06
CO(7–6)	c2	2.06	1.96	1.14	0.69	1.13	1.58	2.05	1.93	1.18	0.80	1.18	1.72
	c3	1.05	1.06	0.80	0.65	0.88	1.24	1.05	1.04	0.79	0.66	0.86	1.14
	zc	6.0	5.2	1.3	0.4	0.2	0.1	6.0	5.2	1.2	0.4	0.1	0.1
	c1	6.22	5.76	5.00	4.42	4.13	3.68	6.20	5.74	5.08	4.48	4.19	3.78
CO(8–7)	c2	1.97	1.25	0.66	1.04	1.85	2.22	1.95	1.27	0.80	0.99	1.83	2.06
	c3	1.04	0.78	0.51	0.62	1.03	1.45	1.02	0.78	0.55	0.58	0.97	1.19
	zc	6.0	6.6	0.5	0.2	0.1	0.1	6.0	1.3	0.5	0.2	0.1	0.1
	c1	5.70	4.95	4.35	4.10	3.84	3.33	5.69	5.03	4.43	4.12	3.93	3.43
CO(9–8)	c2	1.24	0.74	1.05	1.81	2.31	2.17	1.27	0.87	1.04	1.70	2.20	2.27
	c3	0.73	0.46	0.48	0.81	1.14	1.43	0.73	0.51	0.48	0.73	1.04	1.21
	zc	6.2	0.4	0.1	0.1	0.1	0.1	1.2	0.4	0.1	0.1	0.1	0.1
	c1	4.81	4.29	4.03	3.88	3.52	2.95	4.89	4.35	4.07	3.91	3.65	3.00
CO(10–9)	c2	0.86	1.29	1.80	2.29	2.27	2.23	0.95	1.26	1.76	2.21	2.35	2.03
	c3	0.42	0.48	0.66	0.97	1.11	1.47	0.47	0.48	0.64	0.89	1.06	1.06
	zc	0.3	0.1	0.1	0.1	0.1	0.1	0.3	0.1	0.1	0.1	0.1	0.1

Download table as:  ASCIITypeset image

3.2.4. Basic Conclusions

An important conclusion from Figure 5 is that H i surveys at high redshift (z ≳ 2) will be difficult compared to CO surveys. For example, in order to detect the same number of sources at z ≈ 4, an H i survey will need to be approximately 10 times more sensitive than a CO(1–0) survey and over 100 times more sensitive than a CO(5–4) survey. However, we emphasize that with single dish receivers this effect is partially offset by the fact that the instantaneous field of view increases as λ², where λ is the observed wavelength. This scaling implies that single dish H i surveys are likely to yield much more integration time per unit solid angle, which effectively increases the sensitivity. Furthermore, new technologies, such as aperture arrays (Carilli & Rawlings 2004), are currently being developed to provide an enormous instantaneous field of view and hence a very high effective sensitivity for the detection of H i.

Figure 7 (left) shows a comparison of the simulated dN/dz functions for different emission lines observed with an identical peak flux density limit of 1 μJy. The flat slope of the dN/dz function for CO(10–9) reflects that this line is boosted by SBs, which were more abundant and effective (more compact galaxies) at high z (see Obreschkow et al. 2009b). CO(1–0) reveals the steepest slope of all the CO lines in the dN/dz-plot. On one hand, this feature indicates that local galaxies are dominated by low-order excitations of the CO molecule, consistent with empirical data (Braine et al. 1993). On the other hand, CO(1–0) becomes nearly invisible in normal (i.e., no AGN, no SB) galaxies at high redshift (z > 7) due to a near thermal equilibrium between the molecular gas and the CMB (Obreschkow et al. 2009b; see also Combes et al. 1999; Papadopoulos et al. 2000). The even steeper slope of the dN/dz function for H i originates from the cosmic decline of the H₂/H i ratio in galaxies described in Obreschkow & Rawlings (2009b).

Zoom In Zoom Out Reset image size

Figure 7 (right) shows the same dN/dz functions as Figure 7 (left), but for the case of no galaxy evolution. These functions were obtained by constructing a mock observing cone using only the simulation box at z = 0. The comparison of Figure 7 (left) to Figure 7 (right) reveals that H i at high redshifts will be much harder to detect than predicted by a no-evolution model. Qualitatively, the same conclusion applies to low-order CO-emission lines, but the effect is less significant. In contrast, our simulation predicts that the high-order CO-emission lines will be easier to detect than suggested by a no-evolution model, since these lines will be strongly boosted by SBs at high redshift.

Our simulation is inevitably bound to the ΛCDM cosmology with the cosmological parameters given in Section 2. The empirical uncertainty of these parameters may be a source of systematic errors in our predictions. To analyze the errors associated with the uncertainty of the Hubble constant, we can study the change of our predictions in the linear expansion⁷ of h. This analysis shows that varying h between 0.6 and 0.8 does not significantly affect the dN/dz functions, i.e., not more than a factor 2. Additionally, Wang et al. (2008) showed that the lower value for the fluctuation amplitude σ₈ found by WMAP-3 compared to the value used in the Millennium Simulation is almost entirely compensated by an increase in halo bias. Caution should nevertheless be applied when relying on predictions from a single cosmological model.

An additional limitation of the Millennium Simulation is the mass resolution of 8.6 × 10⁸ M_☉ per particle. This mass scale sets the completeness limit in our hydrogen simulation to $M_{{\rm H}\,\hbox{\scriptsize{\sc i}}}+M_{{\rm H}_2}\approx 10^8\;{M}_{\odot }$ (Section 3.2; Obreschkow et al. 2009a). Moreover, galaxies with $M_{{\rm H}\,\hbox{\scriptsize{\sc i}}}+M_{{\rm H}_2}\lesssim 10^9\;{M}_{\odot }$ normally sit a the centers of dark matter halos with poorly resolved merger histories. Therefore, their properties may not have converged in the semi-analytic simulation (Obreschkow et al. 2009a; Croton et al. 2006).

A long list of limitations associated with the semi-analytic galaxy simulation and our post-processing to assign extended H i and H₂ properties has been considered by Obreschkow et al. (2009a). The bottom line of this discussion is that, at z ≳ 5, the simulation becomes very uncertain because the geometries and matter content of regular galaxies are virtually unconstrained from an empirical viewpoint. The young age and short merger intervals of these galaxies may, in fact, have caused them to deviate substantially from the simplistic disk-gas model. At z ≲ 5, the predictions of our H i and H₂ properties are more certain, as they are consistent with available observations. For example, two measurements of CO(2–1)-line emission in regular galaxies at z ≈ 1.5 (Daddi et al. 2008) are consistent with the H₂-MF at this redshift (Obreschkow & Rawlings 2009b). Furthermore, the predicted comoving space density of H₂ evolves proportionally to the observed space density of star formation (e.g., Hopkins 2007) within a factor 2 out to at least z = 3. At z = 0, the simulated H i-MF and CO(1–0)-luminosity function are consistent with the observations of Zwaan et al. (2005) and Keres et al. (2003). Additionally, the local sizes and line widths of H i and CO match the local observations (Obreschkow et al. 2009a, and references therein).

We shall now highlight some specific limitations associated with the emission lines considered in this paper.

We emphasize that at high redshift, the simulated cosmic H i-space density Ω_H i falls below the inferences from Lyα absorption against distant QSOs by a factor ∼2. As mentioned in Obreschkow & Rawlings (2009b), this could reflect a serious limitations of the semi-analytic models implied by the treatment of all cold hydrogen (H i+H₂) as a single phase. Consequently, our dN/dz predictions for H i could be slightly pessimistic. If we believe the empirical estimations of Ω_H i, the offset of our H i masses by a factor ∼2 can be readily accounted for by artificially decreasing the flux limit of the simulated survey by a factor 2. For typical H i surveys in the redshift range z = 0.5–10, this would increase the number of detectable sources by a factor 2–4.

We have limited our predictions for H i to H i emission from galaxies. However, in the EoR, the IGM was not completely ionized and therefore acted as an additional source of H i emission or absorption (Iliev et al. 2002). It may therefore be necessary to analyze the implications of intergalactic H i on the detectability of galactic H i at z ≳ 6 (Becker et al. 2001). On a theoretical level, such an analysis could result from combining the simulation presented in this paper with a simulation of the EoR (e.g., Baek et al. 2009; Santos et al. 2008).

The main discussion of these limitations is given in Obreschkow & Rawlings (2009b), where we introduced our model for the conversion between H₂ and CO. The most serious sources of uncertainty appear to be the heating of molecular gas by SBs and AGNs, the overlap of molecular clouds at high redshift, and the possible presence of nuclear molecular disks in high-redshift galaxies. By contrast, the often-discussed effects of the CMB and the cosmic evolution of the metallicity seem relatively well understood today. Overall, the uncertainty in the predicted CO luminosities increases with redshift and with the J level of the CO transition.

The highest uncertainties, i.e., those for the higher order CO lines at high redshift, can be close to a factor 10. The dN/dz functions in this regime are therefore expected to deviate significantly from our predictions. Such deviations will uncover much of the physics of CO-line emission. In fact, in Obreschkow & Rawlings (2009b), we have explained in detail how different deviations of the CO-luminosity functions from our predictions can be translated into physical interpretations.

The largest coherence scale of our sky simulation is defined by the size of the periodic simulation box of the underlying dark matter simulation (Millennium Simulation, Springel et al. 2005). The side length of this box is s_box = 500 h⁻¹ Mpc, which sets the smallest extractable wave number to k = 2π/s_box ≈ 0.013 h. This value is comparable to the wave number of the first peak in the CDM power spectrum (e.g., Springel et al. 2005). Therefore, the presented simulation allows us to study the power spectrum of H i- and CO-lines and to extract the baryon acoustic oscillations (BAOs); however the position and the amplitude of the first peak of the BAOs will be very poorly constrained.

By contrast, the SKA will have the potential to improve on present measurements of the baryonic power spectrum by at least an order of magnitude in amplitude, and it will detect power in spatial frequencies far below the first acoustic peak. Such a detection could set a primordial constraint on cosmological parameters, especially on the equation of state of dark energy (Blake et al. 2004; Abdalla et al. 2009). Therefore, a simulation of such a detection is regarded as a necessary step in designing the SKA. Yet, this requirement represents a major challenge since no current simulation of cosmic structure is large enough to accurately follow the largest acoustic oscillations, while simultaneously resolving structures small enough to allow the assembly of typical galaxies.

A circumvention of this numerical predicament could result from merging two simulations with different length-scales (see, e.g., Angulo et al. 2008). For example, we could adopt the Horizon-4π dark matter simulation (Prunet et al. 2008; Teyssier et al. 2009), which has a giant box side length of s_box = 2 h⁻¹ Gpc, yet 10 times less mass resolution than the Millennium Simulation. Each dark matter halo of the Horizon-4π Simulation could then be populated with the resolved dark matter substructure and the galaxies contained in comparable halos of the Millennium Simulation.