Ingredients for 21 cm Intensity Mapping

The ΛCDM model describes how the initial quantum perturbations in the primordial universe grow and give rise to the cosmic web: large accumulations of matter in the form of dark matter halos accrete material through filaments and sheets that surround enormous diffuse regions in space. This model has been successful in explaining a very diverse set of cosmological observables, including, among others, the anisotropies in the cosmic microwave background (CMB), the spatial distribution of galaxies, the statistical properties of the Lyα forest, the abundance of galaxy clusters, and correlations in the shapes of galaxies induced by gravitational lensing.

The ΛCDM model has free parameters that describe physical quantities such as the geometry of the universe, the amount of cold dark matter (CDM) and baryons, the sum of the neutrino masses, the expansion rate of the universe, the nature of dark energy, and the initial conditions of the universe. The current quest in cosmology is to determine the values of these parameters as precisely as possible, by exploiting the fact that they influence the spatial distribution of matter. Thus, by examining the statistical properties of matter tracers such as galaxies and cosmic neutral hydrogen, the spatial distribution of matter can be inferred and the value of the cosmological parameters can be constrained.

The amount of information that can be extracted from cosmological surveys depends on several factors, such as the volume being covered or the range in scales where theoretical predictions are reliable. For example, in the case of the CMB, theoretical predictions are extremely precise, because the radiation we observe was produced when the fluctuations were in the linear regime. Tracing the large-scale structure of the universe at low redshifts, through spectroscopic galaxy surveys, represents a complementary approach to extracting cosmological information, where much larger volumes can be surveyed but theoretical predictions are more uncertain. For galaxy surveys, the volume that can be probed also limits the method, because at high redshifts, galaxies are fainter and their spectroscopic detection is challenging.

A different way to trace the matter field is through 21 cm intensity mapping (Bharadwaj & Sethi 2001; Bharadwaj et al. 2001; Battye et al. 2004; McQuinn et al. 2006; Chang et al. 2008; Loeb & Wyithe 2008; Bull et al. 2015; Santos et al. 2015; Villaescusa-Navarro et al. 2015). The method consists of carrying out a low angular resolution survey where the total 21 cm flux from unresolved sources is measured on large areas of the sky at different frequencies. The emission arises from the hyperfine splitting of the ground state of neutral hydrogen into two levels because of the spin–spin interaction between the electron and proton. An electron located in the upper energy level can decay into the lower energy state by emitting a photon with a rest wavelength of 21 cm. This method has several advantages over traditional approaches. First, given that the observable is the 21 cm line, the method is spectroscopic in nature. Second, very large cosmological volumes can be surveyed in a fast and efficient manner. Third, the amplitude of the signal depends only on the abundance and clustering of neutral hydrogen, and so cosmic H i can be traced from z = 0 to z ≃ 50.¹⁴

Current, upcoming, and future surveys such as the Giant Meterwave Radio Telescope (GMRT),¹⁵ the Ooty Radio Telescope (ORT),¹⁶ the Canadian Hydrogen Intensity Mapping Experiment (CHIME),¹⁷ the Five-hundred-meter Aperture Spherical Telescope (FAST),¹⁸ Tianlai,¹⁹ BINGO (Baryon acoustic oscillations In Neutral Gas Observations),²⁰ ASKAP (The Australian Square Kilometer Array Pathfinder),²¹ MeerKAT (The South African Square Kilometer Array Pathfinder),²² HIRAX (The Hydrogen Intensity and Real-time Analysis eXperiment),²³ and the SKA (The Square Kilometer Array)²⁴ will sample the large-scale structure of the universe in the post-reionization era by detecting 21 cm emission from cosmic neutral hydrogen (Choudhuri et al. 2016; Carucci et al. 2017b; Marthi et al. 2017; Sarkar et al. 2018a, 2018b, 2016b).

In order to extract information from those surveys, the observational data has to be compared with theoretical predictions. To linear order, the amplitude and shape of the 21 cm power spectrum is given by

$$\begin{eqnarray}&&{P}_{21\mathrm{cm}}(k,\mu )={\bar{T}}_{b}^{2}[{({b}_{{\rm{H}}{\rm{I}}}+f{\mu }^{2})}^{2}{P}_{{\rm{m}}}(k)+{P}_{\mathrm{SN}}],\end{eqnarray} \tag{ 1 }$$

where ${\bar{T}}_{b}$ is the mean brightness temperature, b_{H i} is the H i bias, f is the linear growth rate, μ = k_z/k, ${P}_{{\rm{m}}}(k)$ is the linear matter power spectrum, and P_SN is the H i shot noise. Here, k_z is the projection of k along the line of sight, which we take to be the z-axis.

At redshifts z ∈ [0, 5], we have relatively good knowledge of the abundance of cosmic neutral hydrogen, and therefore, of ${\bar{T}}_{b}$. On the other hand, little is known about the value of the H i bias and H i shot noise in that redshift interval. It is important to determine their values since the signal-to-noise ratio and range of scales where information can be extracted critically depend on them. One of the purposes of our work is to measure these quantities at different redshifts. Moreover, it is important to determine the regime where linear theory is accurate. In this work, we investigate in detail at which redshifts and scales the clustering of H i in real space (i.e., the H i bias) and in redshift space (Kaiser factor) becomes nonlinear.

In order to optimize what can be learned from the surveys mentioned above, theoretical predictions in the mildly and fully nonlinear regimes are needed. The halo model provides a reasonably accurate framework for predicting the abundance and clustering of H i from linear to fully nonlinear scales. To apply this method, several ingredients are needed for a given cosmological model: (1) the linear matter power spectrum, P_lin(k, z), (2) the halo mass function, n(M, z), (3) the halo bias, b(M, z), (4) the average H i mass that a halo of mass M hosts at redshift z, M_{H i}(M, z), which we refer to as the halo H i mass function,²⁵ and (5) the mean density profile of neutral hydrogen within halos of mass M at redshift z, ${\rho }_{{\rm{H}}{\rm{I}}}(r| M,z)$. In addition, the halo model is formulated under the assumption that all H i are confined to dark matter halos. With the above ingredients in hand, one can write the fully nonlinear H i power spectrum as the sum of one-halo and two-halo terms:

$$\begin{eqnarray}&&{P}_{{\rm{H}}{\rm{I}}}(k,z)={P}_{{\rm{H}}{\rm{I}},1{\rm{h}}}(k)+{P}_{{\rm{H}}{\rm{I}},2{\rm{h}}}(k),\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}\begin{array}{rcl}{P}_{{\rm{H}}{\rm{I}},1{\rm{h}}}(k,z) & = & \displaystyle \frac{1}{{({\rho }_{{\rm{c}}}^{0}{{\rm{\Omega }}}_{{\rm{H}}{\rm{I}}}(z))}^{2}}{\displaystyle \int }_{0}^{\infty }{dMn}(M,z){M}_{{\rm{H}}\,{\rm{I}}}^{2}\\ & & \times \,(M,z){| {u}_{{\rm{H}}{\rm{I}}}(k| M,z)| }^{2},\end{array}\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\begin{array}{l}{P}_{{\rm{H}}{\rm{I}},2{\rm{h}}}(k,z)=\displaystyle \frac{{P}_{\mathrm{lin}}(k,z)}{{({\rho }_{{\rm{c}}}^{0}{{\rm{\Omega }}}_{{\rm{H}}{\rm{I}}}(z))}^{2}}\\ \qquad \times \,{\left[{\displaystyle \int }_{0}^{\infty }{dMn}(M,z)b(M,z){M}_{{\rm{H}}{\rm{I}}}(M,z)| {u}_{{\rm{H}}{\rm{I}}}(k| M,z)| \right]}^{2},\end{array}\end{eqnarray} \tag{ 4 }$$

where ${\rho }_{{\rm{c}}}^{0}$ is the critical density of the universe today and ${{\rm{\Omega }}}_{{\rm{H}}{\rm{I}}}(z)={\bar{\rho }}_{{\rm{H}}{\rm{I}}}(z)/{\rho }_{{\rm{c}}}^{0}$, with ${\bar{\rho }}_{{\rm{H}}{\rm{I}}}(z)$ being the mean H i density at redshift z. $u(k| M,z)$ is the Fourier transform of the normalized H i density profile: $u(x| M,z)={\rho }_{{\rm{H}}{\rm{I}}}(x| M,z)/{M}_{{\rm{H}}{\rm{I}}}(M,z)$.

Some of the goals of our work are to quantify (1) the amount of H i outside of halos, (2) the form of the halo H i mass function, and (3) the density profiles of H i within halos.

While the halo model is a powerful analytic framework, it does not model accurately a number of things, e.g., the transition between the one-halo and two-halo terms (Massara et al. 2014). Thus, its accuracy can be severely limited by that. A more precise modeling can be achieved by painting H i on top of dark matter halos according to the H i halo model ingredients (Villaescusa-Navarro et al. 2014), i.e., more like H i halo occupation distribution (HOD) modeling.

Since 21 cm intensity mapping observations are carried out in redshift space, modeling the abundance and spatial distribution of H i in halos is not enough. A complete description also requires knowing the distribution of H i velocities. In this work, we investigate the H i bulk velocities, the velocity dispersion of H i inside halos, and the amplitude of the Fingers of God in the power spectrum.

The standard halo model does not account for the various complexities expected in the real universe, e.g., whether the H i density profiles depend not only on mass but also on the galaxy population (blue/red) of the halo, whether the clustering of halos depends not only on mass but also on environment, and so forth. These questions can, however, be addressed with hydrodynamic simulations, and in this paper we investigate them in detail.

We also study some quantities that can help us to improve our knowledge of the spatial distribution of H i: the probability distribution function of H i, the relation between the overdensities of matter and H i, the contribution of central and satellite galaxies to the total H i mass content in halos, the H i column density distribution function, and the cross-section of damped Lyα systems (DLAs).

We carry out our analysis using the IllustrisTNG Project, state-of-the-art hydrodynamic simulations that follow the evolution of billions of resolution elements representing CDM, gas, black holes, and stars in the largest volumes ever explored at such mass and spatial resolution. Given the realism of our hydrodynamic simulations, we always aim to connect our results to the underlying physical processes. We note that previous works have studied the H i content of simulated galaxies in detail (Duffy et al. 2012; Davé et al. 2013; Bird et al. 2014; Lagos et al. 2014; Bahé et al. 2016; Faucher-Giguère et al. 2016; Crain et al. 2017; Marinacci et al. 2018; Xie et al. 2017; Zoldan et al. 2017).

We also show that once the most important ingredients for modeling the abundance and clustering properties of H i have been calibrated using full hydrodynamic simulations, less costly dark-matter-only simulations or approximate methods such as COLA (Tassev et al. 2013), peak-patch, or Pinocchio (Monaco et al. 2002; Munari et al. 2017) can be used to generate accurate 21 cm maps. Those maps can then be used to study other properties of H i in the fully nonlinear regime, such as the 21 cm bispectrum or the properties of H i voids. In this work, we investigate the accuracy achieved by creating 21 cm maps from N-body simulations with respect to full hydrodynamic simulations.

This paper is organized as follows. In Section 2, we describe the characteristics of the IllustrisTNG simulations and the method we use to estimate the mass of neutral hydrogen associated with each gas cell. We consider different properties of the abundance of H i in Sections 3–12:

• 1.  
  In Section 3, we compare the overall H i abundance in our simulations to observations.
• 2.  
  In Section 4, we quantify the fraction of H i within halos and galaxies.
• 3.  
  In Section 5, we study the halo H i mass function.
• 4.  
  In Section 6, we investigate the density profiles of H i inside halos.
• 5.  
  In Section 7, we quantify the fraction of the H i mass in halos in the central and satellite galaxies.
• 6.  
  In Section 8, we examine the probability distribution function (pdf) of the H i density and compare it with the total matter density pdf.
• 7.  
  In Section 9, we compute the H i column density distribution function for the absorbers with high column density and quantify the cross-section and bias of DLAs.
• 8.  
  In Section 10, we consider the bulk velocities of H i inside halos.
• 9.  
  In Section 11, we investigate the velocity dispersion of H i inside halos and compare it against that of matter.
• 10.  
  In Section 12, we quantify the relation between the overdensity of matter and H i.

We investigate H i clustering in Sections 13–16:

• 1.  
  In Section 13, we present the amplitude and shape of the H i bias and investigate how well perturbation theory can reproduce H i clustering in real space.
• 2.  
  In Section 14, we show that the clustering of dark matter halos in general depends on their H i masses for fixed halo mass.
• 3.  
  In Section 15, we quantify the amplitude of the H i shot noise.
• 4.  
  In Section 16, we study the clustering of H i in redshift space.

In Section 17, we estimate the accuracy that can be achieved by simulating H i through a combination of N-body simulations with the results derived in the previous sections rather than through full hydrodynamic simulations. Finally, we provide the main conclusions of our work in Section 18. During the course of the discussion, we provide fitting formulae that can be used to reproduce our results.

We emphasize that our results are not converged against resolution (see Appendix A). This should, however, not be seen as a limitation of our claims, as we present results for TNG100, i.e., at the resolution to which the model parameters have been tuned to reproduce a set of galaxy properties.

2.1. The IllustrisTNG Simulations

2.2. Modeling the Hydrogen Phases

We now describe the method we use to quantify the fraction of hydrogen that is in each phase (neutral, ionized, or molecular) for each Voronoi cell in the simulation.

For non-star-forming gas, we use the division between neutral and ionized mass fractions that is calculated in the IllustrisTNG runs on the fly and is included in the simulation outputs. This breakdown assumes primordial chemistry in photoionization equilibrium with the cosmic background radiation (Faucher-Giguère et al. 2009), including a density-dependent attenuation thereof to account for self-shielding following Rahmati et al. (2013a).

For star-forming gas, we post-process the outputs of the simulations to account for the multiphase interstellar medium, including the presence of molecular hydrogen, H₂. The values stored in the simulation output are based on the mass-weighted temperature between the cold and hot phases according to the Springel & Hernquist (2003) model and are thus expected to underestimate the neutral hydrogen fraction. Instead, we set the temperature of star-forming cells to T = 10⁴ K and recalculate the equilibrium neutral hydrogen fraction, also including the self-shielding correction.

The above procedure gives the fraction of hydrogen that is neutral: M_NH/M_H, with ${M}_{\mathrm{NH}}={M}_{{\rm{H}}{\rm{I}}}+{M}_{{{\rm{H}}}_{2}}$. We then compute the H₂ fraction, ${f}_{{{\rm{H}}}_{2}}={M}_{{{\rm{H}}}_{2}}/{M}_{\mathrm{NH}}$, employing the Krumholz, McKee, & Tumlinson (KMT) model (Krumholz et al. 2008, 2009; McKee & Krumholz 2010; see also Sternberg et al. 2014), which we briefly review here.

The molecular hydrogen fraction, ${f}_{{{\rm{H}}}_{2}}$, which we assume to be non-zero only for star-forming gas, is estimated through

$$\begin{eqnarray}{f}_{{{\rm{H}}}_{2}}=\left\{\begin{array}{ll}1-\displaystyle \frac{0.75s}{1+0.25s} & \mathrm{if}\ s\,\lt \,2\\ 0 & \mathrm{if}\ s\geqslant 2\end{array}\right.,\end{eqnarray} \tag{ 5 }$$

where s is given by

$$\begin{eqnarray}&&s=\displaystyle \frac{\mathrm{log}(1+0.6\chi +0.01{\chi }^{2})}{0.6{\tau }_{c}},\end{eqnarray} \tag{ 6 }$$

and

$$\begin{eqnarray}&&\chi =0.756(1+3.1{Z}^{0.365}),\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{\tau }_{c}={\rm{\Sigma }}{\sigma }_{d}/{\mu }_{{\rm{H}}}.\end{eqnarray} \tag{ 8 }$$

In the above equations, Z represents the gas metallicity in units of solar metallicity (Allende Prieto et al. 2001); σ_d is the cross-section of dust, which we estimate from σ_d = Z × 10⁻²¹ cm²; μ_H is the mean mass per hydrogen nucleus, ${\mu }_{{\rm{H}}}\,=2.3\times {10}^{-24}\,{\rm{g}}$; and Σ is the surface density of the gas, which we compute as Σ = ρR, where ρ is the gas density and $R={(3V/4\pi )}^{1/3}$ with V being the volume of the Voronoi cell.

It is possible that our treatment may underestimate the H₂ fractions since (1) the molecular hydrogen fractions go to zero at low densities in the KMT model, (2) we assign molecular hydrogen only to star-forming cells, and (3) it is pessimistic to estimate the surface density from the cell radii. However, we believe that a more precise treatment of H₂ will not affect our results, as its overall abundance is small and therefore not much H i will be transformed into H₂. In order to test this more explicitly, we have considered two extreme cases in which (1) no H₂ is modeled, and (2) all hydrogen in star-forming cells is in molecular form. We have computed the value Ω_{H i}(z) (see Section 3) and did not find significant changes for the first case. On the other hand, at z ∈ [1–3], the prescription where we put all hydrogen in star-forming cells in molecular form can change results by a factor of ≃2–3. We note that Diemer et al. (2018) have recently investigated the impact of different H₂ models on the H i mass function, finding that the results are rather insensitive to the method used to simulate molecular hydrogen at z = 0, although larger differences show up at higher redshifts. We thus believe that our conclusions are robust against our H₂ treatment at low redshift. Quantifying the uncertainty in our results arising from the H₂ model at higher redshift is beyond the scope of this paper.

We note that in our approach we have considered ionization only from the UV background. In other words, we are ignoring the contribution of ionizing photons from, e.g., local sources (Miralda-Escudé 2005; Schaye 2006; Rahmati et al. 2013b) or X-ray heating from the intracluster medium (Kannan et al. 2016).

Figure 1 shows the z = 1 spatial distribution of H i and gas in slices of 5 h⁻¹ Mpc depth throughout the entire TNG100 simulation box as well as in zoomed-in regions thereof. We see that the Lyα forest dominates the abundance of H i in terms of volume, but the H i inside galaxies dominates in terms of mass.

Zoom In Zoom Out Reset image size

Here, we study the overall abundance of neutral hydrogen in the IllustrisTNG simulations. In Figure 2, we show the value of Ω_{H i}(z) from TNG300 (solid black) and TNG100 (dashed black). In this plot, we also indicate measurements from different observations (Rao et al. 2006; Zwaan & Prochaska 2006; Lah et al. 2007; Martin et al. 2010; Songaila & Cowie 2010; Braun 2012; Noterdaeme et al. 2012; Delhaize et al. 2013; Rhee et al. 2013; Crighton et al. 2015).

Zoom In Zoom Out Reset image size

The agreement between the results from our simulations and observations is good, although the simulations tend to overpredict the amount of H i at redshifts z < 0.5 and underpredict the H i abundance at 2 < z < 3.5. Compared to earlier studies with hydrodynamic simulations (e.g., Davé et al. 2013; Bird et al. 2014) and semi-analytic models (Lagos et al. 2014), however, our results agree better with observations at z > 1 and are comparable to the agreement found by Rahmati et al. (2015).

We have investigated whether different H₂ models can improve the agreement between our results and observations at low redshift, finding that they play a very minor role and therefore our conclusions are robust against it. We however believe that the inclusion of the ionization from sources (not incorporated here due to its highly demanding computational needs) can alleviate that tension. Whether this disagreement between observations and simulations is due to the subgrid physics of our simulations or to our approximate modeling of H i is beyond the scope of this work.

The overall H i mass in the simulations depends on resolution, such that the simulation with higher resolution, TNG100, contains between ≃2.5% and ≃40% more H i in the redshift interval z = 0–5 than TNG300. This is mainly a consequence of the fact that TNG300 cannot resolve the smallest halos where H i resides but we believe that this is also due to our H i/H₂ postprocessing (see Appendix A for resolution tests). This can be seen in Figure 4, where we show the H i mass within halos versus the total halo mass. In TNG300, halos with masses only above 2 × 10⁹ h⁻¹ M_⊙ can be resolved, assuming a minimum of 50 CDM particles in a halo. From Figure 4, we see that the amount of H i in halos below that mass is not negligible at high redshift, thus we would expect that Ω_{H i} will be lower at high redshift in TNG300 in comparison with TNG100, as we find. We note that at low redshift, the results of TNG100 and TNG300 are very similar. Since our results for the H i abundance in halos are not converged (see Appendix A), this agreement is just a coincidence. It is not clear to us whether a higher resolution simulation will output a larger or smaller value of Ω_{H i} at z ≃ 0.

In this paper, we examine the most important properties of cosmic neutral hydrogen over a wide range of redshifts. Not being able to resolve the H i that is contained within the smallest halos impacts our results in several ways. For example, the values of the H i bias, H i shot noise, the H i halo mass function, or the amplitude of the H i Fingers-of-God effect will be affected by this. For this reason, from now on we focus our analysis on the TNG100 simulation.

The fraction of the total H i mass that resides within halos is an important ingredient for theoretical frameworks that aim to model the abundance and clustering properties of cosmic neutral hydrogen, such as the halo model (Cooray & Sheth 2002; Barnes & Haehnelt 2014; Villaescusa-Navarro et al. 2014; Padmanabhan & Refregier 2017; Padmanabhan et al. 2016, 2017). In particular, these methods make the assumption that all H i are confined within halos, whose properties, such as spatial distribution or abundance, are well-described by analytic models and/or numerical simulations.

In contrast to the gas in halos, the properties of the gas in the intergalactic medium (IGM) are more difficult to model analytically (see, however, Iršič & McQuinn 2018). The standard approach has been to characterize the gas in the IGM using numerical simulations. If a significant amount of H i is found outside halos, any standard H i halo model will need to be complemented with either simulations or further analytic ingredients. Below, we determine the amount of H i that is outside of halos to quantify the limitations of standard H i halo models.

We have computed the H i mass inside each FoF, FoF-SO, and galaxy in the simulation at several redshifts. In Figure 3, we show, for each of these object types, the fraction of the total H i mass in the simulation that resides inside all objects combined.

Zoom In Zoom Out Reset image size

We find that at redshifts z ≤ 2, more than 99% of all H i is contained within FoF halos. While a significant fraction of the baryons lie outside these regions, the IGM is highly ionized at these times. At these redshifts, the fraction of H i within galaxies is larger than 95%. We note that subfind may not identify any subhalo/galaxy within an FoF halo. This could happen for several reasons, like low density or virialization not having been reached. Thus, we conclude that at these redshifts, ≃5% of the cosmic neutral hydrogen is outside galaxies²⁷ but inside halos.

The fraction of H i within FoF halos and galaxies decreases monotonically with redshift. At redshift z = 5, the H i inside FoF halos only accounts for 88% of the total H i, while the mass within galaxies is 80%. These results are in qualitative agreement with Villaescusa-Navarro et al. (2014), who studied the H i outside halos using a different set of hydrodynamic simulations. Moreover, we find that our results are not significantly affected by mass resolution, since the same analysis carried out for the TNG300 simulation gives similar results.

We consider our finding that the fraction of H i outside halos increases with redshift to be reasonable. At high redshift, the gas in the IGM is denser and the amplitude of the UV background is lower, and so it is easier for that gas to host higher fractions of neutral hydrogen (see Appendix B for further details).

The fraction of H i outside FoF-SO halos is not negligible, varying from 10% at z = 0 to 33% at z = 5. On average, the ratio of H i mass in FoF-SO halos to that in FoF halos is similar to the ratio between their total masses. Thus, FoF's host more H i than FoF-SO halos simply because they are larger and more massive. This also tells us that the regions beyond the virial radius of typical halos are neither H i-poor nor H i-rich, while when this is examined specifically in massive halos we find these regions to be H i rich (see Appendix D for further details).

We thus conclude that while the standard assumption that all H i lies within halos is reasonable at z ≤ 2, at high redshift it begins to break down since a small fraction is located outside halos (∼10% at z = 5). The numbers derived here can be used to quantify the limitations of H i halo models that target the distribution of H i at high redshift.

In the previous section, we have shown that most of the H i is inside halos, justifying the use of H i halo models to characterize the spatial distribution of H i. As discussed in the introduction, besides the linear matter power spectrum, halo mass function, and halo bias, we need to know the halo H i mass function (i.e. the average H i mass hosted by a halo of mass M at redshift z) and the spatial distribution of H i inside halos. Below, we investigate the former: M_{H i}(M, z).

We emphasize the paramount importance of this function by noting that knowing it is sufficient for predicting the amplitude and shape of the 21 cm power spectrum to linear order (see Equation (1)):

$$\begin{eqnarray*}&&{P}_{21\mathrm{cm}}(k,\mu ,z)={\bar{T}}_{b}{(z)}^{2}[{({b}_{{\rm{H}}{\rm{I}}}(z)+f(z){\mu }^{2})}^{2}{P}_{{\rm{m}}}(k,z)+{P}_{\mathrm{SN}}(z)],\end{eqnarray*}$$

where ${\bar{T}}_{b}(z)$, b_{H i}(z), and P_SN(z) can all be derived from M_{H i}(M, z) as

$$\begin{eqnarray}&&{\bar{T}}_{b}(z)=189h\left(\displaystyle \frac{{H}_{0}{(1+z)}^{2}}{H(z)}\right){{\rm{\Omega }}}_{{\rm{H}}{\rm{I}}}(z)\,\mathrm{mK},\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{{\rm{H}}{\rm{I}}}(z)=\displaystyle \frac{1}{{\rho }_{{\rm{c}}}^{0}}{\int }_{0}^{\infty }n(M,z){M}_{{\rm{H}}{\rm{I}}}(M,z){dM},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{b}_{{\rm{H}}{\rm{I}}}(z)=\displaystyle \frac{1}{{\rho }_{{\rm{c}}}^{0}{{\rm{\Omega }}}_{{\rm{H}}{\rm{I}}}(z)}{\int }_{0}^{\infty }n(M,z)b(M,z){M}_{{\rm{H}}{\rm{I}}}(M,z){dM},\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{P}_{\mathrm{SN}}(z)=\displaystyle \frac{1}{{({\rho }_{{\rm{c}}}^{0}{{\rm{\Omega }}}_{{\rm{H}}{\rm{I}}}(z))}^{2}}{\int }_{0}^{\infty }n(M,z){M}_{{\rm{H}}\,{\rm{I}}}^{2}(M,z){dM},\end{eqnarray} \tag{ 12 }$$

where n(M, z) and b(M, z) are the halo mass function and halo bias, respectively. Knowledge of this function can be used to understand the impact of different phenomena on the amplitude and shape of the 21 cm power spectrum such as neutrino masses (Villaescusa-Navarro et al. 2015), warm dark matter (Carucci et al. 2015), or modified gravity (Carucci et al. 2017a).

For each dark matter halo in the simulation, we have computed its enclosed H i mass. In Figure 4, we show the H i mass versus halo mass for each single FoF halo in the simulation at redshifts 0, 1, 2, 3, 4, and 5. The color indicates the number of halos in each bin. We show this map rather than the M_{H i}(M, z) function since the former contains more information, such as the scatter in M_{H i}(M, z). The vertical lines represent a halo mass comprising 50 CDM particles. Its purpose is just to provide the reader with a halo mass for barely resolved halos. We, however, emphasize that 50 CDM particles will not be enough to have internal halo properties, such as concentration, converged.

Zoom In Zoom Out Reset image size

The halo H i mass function increases monotonically with halo mass. Two trends can be identified: (1) in the high-mass end, M_{H i}(M, z) can be approximated by a power law, and (2) in the low-mass end, it has a sharp cutoff. A good fit to our results is given by

$$\begin{eqnarray}&&{M}_{{\rm{H}}{\rm{I}}}(M,z)={M}_{0}{\left(\displaystyle \frac{M}{{M}_{\min }}\right)}^{\alpha }\exp (-{({M}_{\min }/M)}^{0.35}).\end{eqnarray} \tag{ 13 }$$

The free parameters are M_min, which sets the cutoff mass in M_{H i}(M, z); α, which controls the slope of the function at the high-mass end; and M₀, which determines the overall normalization and represents ≃40% of the H i mass of a halo of mass ${M}_{\min }$. We have fitted our results to this function and give the best fitting values for both FoF and FoF-SO halos in Table 1. The green lines in Figure 4 indicate the best fits at each redshift.

At redshifts z ≥ 3, for FoF halos, α is ≃0.75, while it decreases at lower redshifts: α = 0.60 at z = 2, α = 0.53 at z = 1, and α = 0.24 at z = 0. We interpret this as a result of several physical processes, such as AGN feedback and ram-pressure and tidal stripping, being more efficient at removing gas from galaxies at lower redshift than at higher redshifts.

The value of M_min decreases monotonically with redshift, from ≃2 × 10¹² h⁻¹ M_⊙ at z = 0 to 2 × 10¹⁰ h⁻¹ M_⊙ at z = 5. This indicates that as the redshift increases, lower mass halos host H i. In Appendix C, we discuss the physical origin of the cutoff in the halo H i mass function and the relative importance of supernova feedback and gas stripping, and the UV background for this.

Similar conclusions can be reached when computing ${M}_{{\rm{H}}{\rm{I}}}(M,z)$ using FoF-SO halos (see Table 1). The derived values of α, M₀, and M_min are roughly compatible, within the errors, between FoF and FoF-SO halos. There are, however, systematic differences between the best-fit values from FoF and FoF-SO, which is because the total amount of H i in a given halo is sensitive to its definition, as we will see below.

For FoF-SO halos, we can relate halo masses to circular velocities through ${V}_{\mathrm{circ}}=\sqrt{{GM}/R}$, where G is the gravitational constant and M and R are the halo mass and radius. Expressing ${M}_{\min }$ in terms of circular velocities, we obtain ${V}_{\mathrm{circ}}({M}_{\min })\simeq 180\pm 20\,\mathrm{km}\,{{\rm{s}}}^{-1}$ for z ∈ [0, 2] and ${V}_{\mathrm{circ}}({M}_{\min })\,\simeq 120\pm 20\,\mathrm{km}\,{{\rm{s}}}^{-1}$ for z ∈ [3, 5]. This suggests that the minimum halo mass that can host H i depends primarily on the depth of its gravitational potential and that at lower redshifts, the potential has to be deeper since astrophysical processes such as AGN feedback, tidal stripping, and so forth are more effective at removing gas from small halos.

Since our parameterization of the halo H i mass function does not have a "hard" cutoff, the value of M_min only represents a mass scale at which the halo H i mass function changes its trend. In other words, halos with masses around M_min host a significant amount of H i. It is also very interesting to quantify the cutoff in the halo H i mass function more rigidly, i.e., so that halos below a certain mass contain a negligible amount of H i.

We have calculated the halo mass at which 98% of all H i in halos is above that mass and only 2% is in smaller halos. We term this halo mass as a "hard cutoff mass," M_hard,

$$\begin{eqnarray}&&\displaystyle \frac{{\int }_{0}^{{M}_{\mathrm{hard}}}n(M,z){M}_{{\rm{H}}{\rm{I}}}(M,z){dM}}{{\int }_{0}^{\infty }n(M,z){M}_{{\rm{H}}{\rm{I}}}(M,z){dM}}=0.02,\end{eqnarray} \tag{ 14 }$$

and we show the corresponding values in Table 1. For FoF-SO halos, we obtain 1.3 × 10¹⁰, 6.1 × 10⁹, 2.5 × 10⁹, 7.6 × 10⁸, 2.3 × 10⁸, and 1.7 × 10⁸ h⁻¹ M_⊙ at redshifts 0, 1, 2, 3, 4, and 5, respectively. These can be transformed to circular velocities, giving 34, 35, 31, 24, 19, and 18 km s⁻¹ at redshifts 0, 1, 2, 3, 4, and 5, correspondingly. The values we infer do not change much if we use a threshold equal to 99%. We thus conclude that at redshifts z ≤ 2, only halos with circular velocities above about 30 km s⁻¹ host H i, while at redshifts z ≥ 3, the H i is only in halos with circular velocities above ∼20 km s⁻¹.

A more conventional parameterization of the halo H i mass function

$$\begin{eqnarray}&&{M}_{{\rm{H}}{\rm{I}}}(M,z)={M}_{0}{\left(\displaystyle \frac{M}{{M}_{\min }}\right)}^{\alpha }\exp (-{M}_{\min }/M)\end{eqnarray} \tag{ 15 }$$

also reproduces our results. The fit in Equation (13) is, however, preferred for our results since at high redshift, ${M}_{{\rm{H}}{\rm{I}}}(M,z)$ falls more slowly for very low halo masses than the standard profile. In order to facilitate the comparison with works using the above parameterization (e.g., Bagla et al. 2010; Castorina & Villaescusa-Navarro 2017; Obuljen et al. 2018a; Padmanabhan et al. 2017; Villaescusa-Navarro et al. 2017; Pénin et al. 2018), we also provide the best-fit values for the more conventional fit in Appendix E.

Another important ingredient in describing the spatial distribution of cosmic neutral hydrogen using H i halo models is the density profile of H i inside halos (see Equations (2)–(4)). In this section, we investigate the spatial distribution of H i inside simulated dark matter halos.

Since FoF halos can have very irregular shapes, we have computed the H i profiles inside FoF-SO halos. For each FoF-SO halo, we have computed the H i mass within narrow spherical shells up to the virial radius, and from them the H i profile. Figure 5 shows the individual H i profiles for halos in a narrow mass bin at different redshifts with gray lines. The large halo-to-halo scatter is surprising and highlights the fact that individual H i profiles, as opposed to dark matter ones, are far from universal.

Zoom In Zoom Out Reset image size

The scatter is particularly large toward the centers of massive halos, as discussed below, as well as for halos with masses around or below the cutoff we observe in Figure 4. This is expected, as the halo H i mass function also exhibits large scatter in that range. As we will see later in Section 14, the clustering of halos in that mass range depends significantly on their H i mass. Thus, it is likely that the H i content of these halos is influenced by their environment, so small halos around more massive ones may have lost or gained a significant fraction of their H i mass due to related effects.

The scatter generally tends to be lower at higher redshifts and, in particular, is small in halos with masses above 10¹⁰ h⁻¹ M_⊙ at redshift z = 5. This is related to the lower scatter we find at high redshift in the halo H i mass function, M_{H i}(M, z) (see Figure 4). We speculate that this originates from a reduced role that AGN feedback and environmental gas stripping play at earlier times.

The blue lines in Figure 5 show the mean and the standard deviation of the H i profiles from all halos that lie in each mass bin and redshift, while the red lines display the median. Clearly, in some cases they differ substantially. This behavior can be partially attributed to the H i profiles arising from two distinct populations: i.e., H i-rich blue galaxies versus H i-poor red ones (Nelson et al. 2018). This can clearly be seen in the panel in Figure 5 corresponding to halos in the mass range M ∈ [1–2] × 10¹² h⁻¹ M_⊙ at z = 0. In this range, some halos have a core in their H i profiles while others do not. The reason is that the central galaxy of some halos is experiencing AGN feedback (those with holes in the profile) and are therefore becoming red, while the galaxies in the other halos are not yet being affected by AGN feedback, remaining blue (Nelson et al. 2018).

We find that the H i density profiles of small halos (M ≲ 10¹² h⁻¹ M_⊙) increase toward their halo center. We note, however, that the amplitude of the H i profile tends to saturate, i.e., the slope of the profiles declines significantly toward the halo center. For example, at z = 0 and z = 1 and for halos with masses larger than ${10}^{11}\,{h}^{-1}\,{M}_{\odot }$, the mean H i profiles change slope around ∼20 h⁻¹ kpc. This is expected since neutral hydrogen at high densities will turn into molecular hydrogen and stars on short timescales. For higher halo masses (M ≃ 10¹³ h⁻¹ M_⊙), the H i density profile exhibits a hole in the center. This is caused by AGN feedback in the central galaxy of those halos. We note that higher densities in the center of halos can give rise to the formation of molecular hydrogen, which can produce a similar effect (Marinacci et al. 2018). Holes, which extend even further than in groups, are also found in the H i profiles of galaxy clusters, which we, however, do not show here since there are only a few of them and only at low redshift.

We illustrate these features of the H i profiles in Figure 6, where we show the spatial distribution of H i in and around four individual halos with masses 10¹¹, 10¹², 10¹³, and 10¹⁴ h⁻¹ M_⊙ at redshifts 0 or 1. We have selected these halos by requiring that their H i density profiles are close to the mean. It can be seen that H i is localized in the inner regions of small halos, while for groups, the central galaxy exhibits a hole produced by AGN feedback. For galaxy clusters, the central regions have little H i. This happens because the central galaxy is an H i-poor elliptical, and ram-pressure and tidal stripping are very efficient in removing the gas content of galaxies passing near the center. The analysis in this section suggests that analytical approaches to the distribution of DLAs employing a universal H i profile, e.g., (Padmanabhan et al. 2017) will not be able to reproduce observations.

Zoom In Zoom Out Reset image size

In order to quantitatively investigate what the effective average H i density profile is across different halo masses and redshift, we use the mean measured H i density profile and test two models of H i density that both include an exponential cutoff on small scales.

First we consider a simple power law with an exponential cutoff on small scales—Model 1:

$$\begin{eqnarray}&&{\rho }_{{\rm{H}}{\rm{I}}}(r)=\displaystyle \frac{{\rho }_{0}}{{r}^{{\alpha }_{\star }}}\exp (-{r}_{0}/r),\end{eqnarray} \tag{ 16 }$$

where ρ₀ is the overall normalization, α_⋆ is the slope parameter, and r₀ is the inner radius at which the density drops and the profile changes its slope.

Second, we consider an altered Navarro–Frenk–White (NFW) profile (Maller & Bullock 2004; Barnes & Haehnelt 2014), found to be a good fit to the multiphase gas distribution at high redshifts in hydrodynamical simulations, with an exponential cutoff on small scales—Model 2:

$$\begin{eqnarray}&&{\rho }_{{\rm{H}}{\rm{I}}}(r)=\displaystyle \frac{{\rho }_{0}{r}_{s}^{3}}{(r+3/4{r}_{s}){(r+{r}_{s})}^{2}}\exp (-{r}_{0}/r),\end{eqnarray} \tag{ 17 }$$

where ρ₀ is the overall normalization and r_s is the scale radius of the H i cloud. In both cases, the overall normalization, ρ₀, is fixed such that the volume integral of the model density profile integrated up to the virial radius of a given halo matches the mean total H i mass obtained from the density profile found in simulations (blue lines in Figure 5). We are then left with two free parameters for each model: {α_⋆, r₀} and {r_s, r₀}. We fit these models to the measured mean H i density profiles, limiting our analysis only to the scales above r ≥ 2h⁻¹ kpc. For the uncertainties in the density profiles, we use the scatter among different galaxies (blue error bars in Figure 5) and assume that these uncertainties are uncorrelated between different scales.

The best-fit values along with the 68% confidence intervals are presented in Table 2, while in Figure 31 in Appendix F we show the best-fit results for Model 1. Based on the resulting best-fit χ², we find that both Model 1 and 2 are good fits for all the considered redshifts and halo masses, except for the most massive halo bin M_h = 10¹⁴ [h⁻¹ M_⊙] at z = 0. We find the difference in the best-fit χ² between the two models to be negligible. This is to be expected since the models are rather similar and have the same slope on large scales. In the case of Model 1, we find the H i density profile slope to be consistent with a value of α_⋆ = 3 for all halo masses and redshifts. The inner radius r₀ depends on the halo mass and is larger for larger halo masses at a fixed redshift, while at a fixed halo mass, it increases with increasing redshift. For example, for halos with M_h ≤ 10¹¹ h⁻¹ Mpc and z ≤ 2, r₀ is below the minimum scale considered, and the uncertainties are rather large. In the case of Model 2, we find a similar behavior. The inferred values of r₀ are consistent between the two models, with Model 2 having larger uncertainties due to the degeneracy between parameters r₀ and r_s.

Table 2.  Best-fit Values of the Parameters Determining the H i Density Profiles

Model 1—Power Law + Exponential Cutoff: α⋆, log10 r0 [h−1 Mpc]
z	Mh = 109 [h−1 M⊙]		Mh = 1010 [h−1 M⊙]		Mh = 1011 [h−1 M⊙]		Mh = 1012 [h−1 M⊙]		Mh = 1013 [h−1 M⊙]		Mh = 1014 [h−1 M⊙]
0	⋯	⋯	${3.04}_{-0.03}^{+0.04},$	$-{3.59}_{-0.92}^{+0.85}$	${3.03}_{-0.02}^{+0.03},$	$-{2.8}_{-1.2}^{+0.5}$	${3.02}_{-0.03}^{+0.03},$	$-{2.32}_{-1.15}^{+0.33}$	${3.00}_{-0.04}^{+0.04},$	$-{1.71}_{-0.12}^{+0.09}$	${2.92}_{-0.03}^{+0.03},$	$-{1.91}_{-0.14}^{+0.11}$
1	${3.3}_{-0.7}^{+1.3},$	$-{2.5}_{-1.6}^{+1.1}$	${3.05}_{-0.02}^{+0.02},$	$-{3.72}_{-0.84}^{+0.77}$	${3.02}_{-0.02}^{+0.02},$	$-{3.3}_{-1.1}^{+0.7}$	${3.00}_{-0.02}^{+0.03},$	$-{2.32}_{-0.28}^{+0.16}$	${2.99}_{-0.03}^{+0.03},$	$-{1.77}_{-0.11}^{+0.09}$	⋯	⋯
2	${3.07}_{-0.08}^{+0.10},$	$-{3.2}_{-1.2}^{+0.9}$	${3.03}_{-0.02}^{+0.01},$	$-{3.64}_{-0.89}^{+0.78}$	${3.01}_{-0.01}^{+0.01},$	$-{2.75}_{-0.68}^{+0.26}$	${3.00}_{-0.02}^{+0.02},$	$-{2.18}_{-0.12}^{+0.09}$	${2.98}_{-0.01}^{+0.02},$	$-{1.74}_{-0.05}^{+0.04}$	⋯	⋯
3	${3.05}_{-0.02}^{+0.02},$	$-{3.63}_{-0.93}^{+0.85}$	${3.02}_{-0.02}^{+0.02},$	$-{3.1}_{-1.1}^{+0.5}$	${3.00}_{-0.01}^{+0.01},$	$-{2.52}_{-0.20}^{+0.13}$	${3.00}_{-0.02}^{+0.02},$	$-{2.09}_{-0.07}^{+0.06}$	⋯	⋯	⋯	⋯
4	${3.04}_{-0.02}^{+0.02},$	$-{3.3}_{-1.0}^{+0.7}$	${3.00}_{-0.01}^{+0.01},$	$-{2.46}_{-0.24}^{+0.15}$	${3.00}_{-0.01}^{+0.01},$	$-{2.32}_{-0.08}^{+0.07}$	${2.99}_{-0.01}^{+0.01},$	$-{2.04}_{-0.04}^{+0.03}$	⋯	⋯	⋯	⋯
5	${3.03}_{-0.02}^{+0.02},$	$-{2.9}_{-1.2}^{+0.5}$	${3.00}_{-0.01}^{+0.01},$	$-{2.28}_{-0.12}^{+0.09}$	${3.00}_{-0.01}^{+0.01},$	$-{2.18}_{-0.05}^{+0.04}$	${3.00}_{-0.01}^{+0.01},$	$-{2.02}_{-0.03}^{+0.03}$	⋯	⋯	⋯	⋯
Model 2—Altered NFW + Exponential Cutoff: log10rs [h−1 Mpc], log10r0 [h−1 Mpc]
0	⋯	⋯	$-{4.0}_{-0.7}^{+0.7},$	$-{3.8}_{-0.8}^{+0.8}$	$-{3.7}_{-0.9}^{+0.6},$	$-{3.4}_{-1.0}^{+0.7}$	$-{3.2}_{-1.1}^{+0.7},$	$-{3.1}_{-1.3}^{+0.8}$	$-{3.0}_{-1.3}^{+1.0},$	$-{1.8}_{-1.0}^{+0.1}$	$-{2.3}_{-1.7}^{+0.5},$	$-{2.6}_{-1.6}^{+0.7}$
1	$-{2.8}_{-1.5}^{+1.8},$	$-{3.3}_{-1.1}^{+1.2}$	$-{4.0}_{-0.6}^{+0.5},$	$-{3.7}_{-0.8}^{+0.6}$	$-{3.9}_{-0.7}^{+0.6},$	$-{3.6}_{-0.9}^{+0.7}$	$-{3.0}_{-1.2}^{+0.4},$	$-{2.9}_{-1.4}^{+0.6}$	$-{2.2}_{-1.6}^{+0.3},$	$-{2.4}_{-1.7}^{+0.6}$	⋯	⋯
2	$-{3.7}_{-0.8}^{+0.9},$	$-{3.7}_{-0.9}^{+0.9}$	$-{3.8}_{-0.7}^{+0.5},$	$-{3.5}_{-0.9}^{+0.5}$	$-{3.6}_{-0.9}^{+0.5},$	$-{3.3}_{-1.1}^{+0.6}$	$-{2.8}_{-1.3}^{+0.3},$	$-{2.6}_{-1.5}^{+0.4}$	$-{1.8}_{-0.2}^{+0.1},$	$-{3.0}_{-1.3}^{+0.8}$	⋯	⋯
3	$-{3.8}_{-0.8}^{+0.6},$	$-{3.5}_{-0.9}^{+0.6}$	$-{3.6}_{-0.8}^{+0.5},$	$-{3.3}_{-1.1}^{+0.5}$	$-{3.3}_{-1.0}^{+0.4},$	$-{3.0}_{-1.3}^{+0.5}$	$-{2.8}_{-1.3}^{+0.4},$	$-{2.4}_{-1.4}^{+0.3}$	⋯	⋯	⋯	⋯
4	$-{3.6}_{-0.9}^{+0.5},$	$-{3.2}_{-1.1}^{+0.5}$	$-{3.2}_{-1.1}^{+0.4},$	$-{3.0}_{-1.3}^{+0.5}$	$-{2.9}_{-1.3}^{+0.3},$	$-{2.7}_{-1.4}^{+0.4}$	$-{2.6}_{-1.4}^{+0.3},$	$-{2.4}_{-1.2}^{+0.3}$	⋯	⋯	⋯	⋯
5	$-{3.4}_{-1.0}^{+0.5},$	$-{3.1}_{-1.1}^{+0.5}$	$-{2.7}_{-1.1}^{+0.2},$	$-{3.1}_{-1.3}^{+0.7}$	$-{2.5}_{-1.1}^{+0.1},$	$-{3.3}_{-1.2}^{+1.0}$	$-{2.2}_{-0.4}^{+0.1},$	$-{3.1}_{-1.2}^{+0.8}$	⋯	⋯	⋯	⋯

Note. We show the resulting parameters for the two different models considered (see text): an altered NFW profile with an exponential cutoff on small scales (top) and a simple power law with an exponential cutoff on small scales (bottom), as a function of the halo mass (columns) and redshift (rows).

Download table as:  ASCIITypeset image

We note that other observational and simulation studies have found that the H i surface density profile of galaxies can be reproduced by an exponential profile (Obreschkow et al. 2009; Wang et al. 2014). Based on these studies, other spherically averaged density models have been used in the literature, e.g., an exponential profile (Padmanabhan et al. 2017). We find that using an exponential profile for the spherically averaged profile does not reproduce our mean data very well.

It is interesting to quantify what fraction of H i mass inside halos comes from their central and satellite galaxies. This will help us to better understand the H i density profiles (see Section 6) and improves our intuition for the amplitude of the H i Fingers-of-God effect (see Section 16).

For each FoF halo, we have computed its total H i mass, the H i mass within its central galaxy, and the H i mass inside its satellites We emphasize that our definition of central galaxy departs significantly from that used in observations, as we consider the central galaxy to be the most massive subhalo. In general, this subhalo hosts the particle at the minimum of the gravitational potential and is therefore the classical central galaxy, but it also has a significant spatial extent and particles far away from the halo center can be associated with this subhalo, unless they are bound to a satellite.

We take narrow bins in halo mass and compute the average H i mass for each of the above quantities. The outcome is shown in Figure 7. The black lines in the upper panels show the halo H i mass function, while the red and blue ones display the average H i mass inside the central and satellite galaxies as a function of halo mass. The bottom panels show the fraction of H i mass within halos that comes from the central and the satellite galaxies.

Zoom In Zoom Out Reset image size

Aside from very low-mass halos (M ≲ 10⁹ h⁻¹ M_⊙), the fraction of H i in the central galaxy decreases with halo mass, while the fraction of H i in satellites increases, independent of redshift. For small halos, nearly all of the H i are located in the central galaxy, as expected. For halos of masses ∼5 × 10¹² h⁻¹ M_⊙, the fraction of H i in the central galaxy and in satellites is roughly the same, almost independent of redshift. For more massive halos, the total H i mass is dominated by the H i in satellite galaxies.

At high redshift, the contribution of satellites to the total H i mass in small (M ∈ [10¹⁰–10¹¹] h⁻¹ M_⊙) halos is non-negligible: ≃20%. At z = 0 and for galaxy clusters M ≥ 10¹⁴ h⁻¹ M_⊙, the contribution of the central galaxy to the total H i mass is negligible, as expected.

The H i mass in very low-mass halos ($M\lesssim {10}^{9}\,{h}^{-1}\,{M}_{\odot }$) is small, in particular at low redshift. For some of these halos, the sum of the H i mass in the central and satellite subhalos is much less than the total H i mass. In such cases, some of the H i mass was determined by SUBFIND to be unbound. In a fraction of these halos, no bound structure was identified by SUBFIND altogether, rendering the combined H i masses of the central and satellites, which by definition do not exist, to be zero, even if the FoF group contains some H i.

We now study other quantities that, although not ingredients for H i halo models, will help us better understand the spatial distribution of neutral hydrogen. One of those quantities is the density pdf, which we investigate in detail in this section and compare to that of matter.

We compute the density fields of neutral hydrogen and total matter in the whole simulation volume using cloud-in-cell (CIC) interpolation on a grid with 2048³ cells in real space, namely ≈36.6 h⁻¹ kpc across each grid cell. We then smooth those fields with top-hat filters of radii 1 and $5\,{h}^{-1}\mathrm{Mpc}$. We have chosen those values for the smoothing scale, R, as a compromise between large and small scales. On one hand, the volumes of our simulations do not allow us to explore values much larger than ≃5 h⁻¹ Mpc, while on the other hand we take 1 h⁻¹ Mpc as a representative estimate of the nonlinear regime. In Figure 8, we show the pdf's, computed as the number of cells in a given interval in overdensity, over the total number of cells, divided by the width of the overdensity interval.

Zoom In Zoom Out Reset image size

While the density pdf of the matter field is highly non-Gaussian at low redshift, for either of the two smoothing scales we consider here, at high redshift it becomes more nearly Gaussian, as expected. At redshifts z ≥ 4 and for R = 5 h⁻¹ Mpc, the matter pdf can be well approximated by a log-normal distribution

$$\begin{eqnarray}&&\mathrm{pdf}(\delta )=\displaystyle \frac{1}{\sqrt{2\pi {\sigma }^{2}}}\exp \left[\displaystyle \frac{{(\mathrm{log}(1+\delta )+{\sigma }^{2}/2)}^{2}}{2{\sigma }^{2}}\right],\end{eqnarray} \tag{ 18 }$$

where we consider σ to be a free parameter. We find σ = 0.21 at z = 5 and σ = 0.28 at z = 4. At lower redshifts and for smaller values of the smoothing scale, the log-normal function does not provide a good fit to our results, as expected (see e.g., Uhlemann et al. 2016).

The H i density pdf exhibits a different behavior compared to the matter pdf. First, the abundance of large H i overdensities remains roughly constant with redshift, independent of the smoothing scale considered. Second, for a smoothing scale of $1\,{h}^{-1}\,\mathrm{Mpc}$, the H i pdf hardly changes with redshift.

For redshifts z ≥ 3, a log-normal distribution characterizes our results relatively well: for R = 1 h⁻¹ Mpc, σ ≃ 1.9, while for for R = 5 h⁻¹ Mpc. σ ≃ 1 at z = 3, σ ≃ 0.9 at z = 4, and σ ≃ 0.8 at z = 5. At lower redshifts, a log-normal distribution does not provide a good match to the simulations.

To understand the physical origin of the differences between the pdf's of H i and matter density, it is useful to relate the width of the pdf to the amplitude of the H i power spectrum. This is possible since the amplitude of the H i power spectrum represents a measurement of the variance of the field at a given scale. Low values of the H i power spectrum indicate that H i is distributed homogeneously, while higher values mean that spatial variations in H i density can be large.

One of the reasons that the H i density pdf is roughly similar across redshifts—while this is less true of the matter density pdf—is that the amplitude of the H i power spectrum depends more weakly on redshift than does the matter power spectrum (see Section 13). Thus, the variance of the H i pdf is necessarily smaller than that of the matter pdf. Since the amplitude of the H i power spectrum is larger than that of the matter power spectrum at high redshifts, the variance of the H i pdf is larger than that of the matter pdf at those redshifts, as we find. Finally, in the central region of a void, the matter density will be low, but the H i density will be even lower.²⁸ Thus, we find that there will be more cells with low H i overdensity than with low matter overdensity.

It can be seen for both H i and matter that the distributions are broader for a smoothing scale of 1 h⁻¹ Mpc than for 5 h⁻¹ Mpc. This is expected since when smoothing over larger scales, any field will become more homogeneous and therefore the width of its pdf will become smaller. This can also be quantified through the amplitude of the power spectrum, using the same reasoning as above.

Another quantity commonly employed to study the abundance of neutral hydrogen in the post-reionization era is the H i column density distribution function (H i CDDF), defined as

$$\begin{eqnarray}&&{f}_{{\rm{H}}{\rm{I}}}({N}_{{\rm{H}}{\rm{I}}})=\displaystyle \frac{{d}^{2}n({N}_{{\rm{H}}{\rm{I}}})}{{{dN}}_{{\rm{H}}{\rm{I}}}{dX}},\end{eqnarray} \tag{ 19 }$$

where n is the number of lines of sight with column densities between N_{H i} and N_{H i} + dN_{H i}, and ${dX}={H}_{0}{(1+z)}^{2}/H(z){dz}$ is the absorption distance. This quantity can be inferred directly from observations of the Lyα forest.

Here, we investigate the H i CDDF, focusing on absorbers with high column densities: DLAs, ${N}_{{\rm{H}}{\rm{I}}}\geqslant {10}^{20.3}\,{\mathrm{cm}}^{-2}$. We also examine the DLA cross-sections, which are required both observationally (Font-Ribera et al. 2012; Alonso et al. 2018; Pérez-Ràfols et al. 2018) and theoretically (Castorina & Villaescusa-Navarro 2017).

In Figure 9, we show an example of the spatial distribution of gas and H i around a massive halo at redshift z = 3. As in this case, DLAs correspond to gas in galaxies, gas recently stripped from galaxies, and gas in streams.

Zoom In Zoom Out Reset image size

The H i CDDF at redshifts 0, 1, 2, 3, 4, and 5 is computed using the following procedure (we refer the reader to Appendix B of Villaescusa-Navarro et al. 2014 for further details). We approximate each Voronoi cell by a uniform sphere with radius equal to $R={(3V/4\pi )}^{1/3}$, where V is the volume of the cell, and determine the H i column density of a line through it from ${N}_{{\rm{H}}{\rm{I}}}={\rho }_{{\rm{H}}{\rm{I}}}d$, where ρ_{H i} = M_{H i,cell}/V and d is the length of the segment intersecting the sphere. The simulation volume is projected along the z-axis and a grid with 20,000 × 20,000 points is overlaid. Each point is considered to be a line of sight, and the column density along it is estimated as the sum of the column densities of all Voronoi cells contributing to it. Since our box size is relatively small, the probability of encountering more than a single absorber with a large column density along the line of sight is negligible. Thus, if the column density of a given line of sight is larger than ∼10¹⁹ cm⁻², it can be attributed to a single absorber. We repeated the tests carried out in Villaescusa-Navarro et al. (2014) to verify that (1) the grid is fine enough to achieve convergence in the CDDF, and (2) the results do not change if the CDDF is computed by slicing the box into slabs of different widths.

We show the results in Figure 10. We find excellent agreement with the observations, which are shown as black points with error bars, at redshifts [1.8–3.5] (Péroux et al. 2005), [2.0–3.5] (Noterdaeme et al. 2012), [3.5–5.4] (Crighton et al. 2015), and [1.5–5.0] (Zafar et al. 2013). The differences between the observed and simulated CDDFs, e.g., the amplitude of the H i CDDF around 10²⁰–10²¹ cm⁻², are related to the mismatch between ${{\rm{\Omega }}}_{{\rm{H}}{\rm{I}}}$ from observations and TNG100 (see Figure 2), since

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{{\rm{H}}{\rm{I}}}(z)=\displaystyle \frac{{m}_{{\rm{H}}}{H}_{0}}{c{\rho }_{c}^{0}}{\int }_{0}^{\infty }{f}_{{\rm{H}}{\rm{I}}}({N}_{{\rm{H}}{\rm{I}}},z){N}_{{\rm{H}}{\rm{I}}}{{dN}}_{{\rm{H}}{\rm{I}}},\end{eqnarray} \tag{ 20 }$$

where m_H is the mass of the hydrogen atom and c is the speed of light. In agreement with previous works, we find that the H i CDDF exhibits a weak dependence on redshift (see, e.g., Rahmati et al. 2013a). This self-similarity can be associated with the weak redshift dependence that we observe in the high-overdensity tail of the H i density pdf for small smoothing scales (see Figure 8).

Zoom In Zoom Out Reset image size

Next, we examine the cross-section of DLAs. For each dark matter halo of the simulation, the area covered by DLAs with different column densities is computed. Then, all halos within mass bins are selected, and the mean and standard deviation of their DLA cross-sections are determined. As shown in Figure 11, we find that, for fixed column density, the DLA cross-section increases with halo mass, while the cross-section decreases with column density for halos of fixed mass.

Zoom In Zoom Out Reset image size

The cross-section of the DLAs is well fitted by the following function

$$\begin{eqnarray}&&\sigma (M| {N}_{{\rm{H}}{\rm{I}}},z)=A{\left(\displaystyle \frac{M}{{h}^{-1}{M}_{\odot }}\right)}^{\alpha }(1-{e}^{-{(M/{M}_{0})}^{\beta }}).\end{eqnarray} \tag{ 21 }$$

Here, A is a parameter that controls the overall normalization of function, while α sets the slope of the cross-section for large halo masses, and M₀ determines the characteristic halo mass where the DLA cross-section exponentially decreases at a rate controlled by β.

We fit our results at redshifts z ∈ [2, 4] using the above form and find α = 0.82 in the large majority of the cases, while β is well approximated by $\beta =0.85\cdot {\mathrm{log}}_{10}({N}_{{\rm{H}}{\rm{I}}}/{\mathrm{cm}}^{-2})-16.35$. There is also a strong correlation between A and M₀, given by $A\cdot {M}_{0}=0.0141\,{h}^{-2}\,\mathrm{kpc}\,{M}_{\odot }$. The only dependence on redshift enters through M₀, the value of which is given in Table 3. We find that M₀ decreases with redshift, in agreement with the halo H i mass function, which implies that less massive halos can host H i at higher redshifts.

The fits to the simulation results, shown as dashed lines in Figure 11, are a good approximation for column densities below 10²² cm⁻², but apparently less so at higher column densities; e.g., the fit for column densities above 10^22.5 cm⁻² and low halo masses is several orders of magnitude below the mean. This is mainly an illusion of the fact that some error bars are larger than the value. The reduced χ² obtained from the fits in all cases is below 0.35. The preferred value for α is slightly larger at higher redshifts, but the redshift dependence is so weak that for simplicity we did not use it in our fitting. The largest discrepancy between the fit and our results occurs at z = 2 for the DLAs with column densities larger than 10²² cm⁻². In Castorina & Villaescusa-Navarro (2017), it was suggested that a very good fit to the column density distribution and the DLA bias can be obtained assuming the differential cross-section, dσ/dN_{H i}, is roughly independent of column density. This implies that the linear bias of different absorbers will be very similar, and the measurements in the BOSS survey of Pérez-Ràfols et al. (2018) confirm this simple picture, although with large error bars. As discussed above, our fit to Equation (21) indicates that the slope α to a very good approximation is not a function of column density, and the N_{H i} dependence of β can be moved to the normalization constant A using their tight correlation. If we then look at the expression for the linear bias of a given absorber,

$$\begin{eqnarray}&&{b}_{{N}_{{\rm{H}}{\rm{I}}}}(z)=\displaystyle \frac{{\int }_{0}^{\infty }b(M,z)n(M,z)\sigma (M| {N}_{{\rm{H}}{\rm{I}}},z){dM}}{{\int }_{0}^{\infty }n(M,z)\sigma (M| {N}_{{\rm{H}}{\rm{I}}},z){dM}},\end{eqnarray} \tag{ 22 }$$

we note that A cancels between the numerator and the denominator, and the only column density dependence is left in β, and it is rather small. The analytical calculation in Castorina & Villaescusa-Navarro (2017) therefore agrees with the measurements in IllustrisTNG, and future observations will tell us whether or not our current understanding of the cross-section is correct.

We have used the above expression to estimate the bias of the DLAs. We take the DLAs cross-section (for absorbers with N_{H i} > 10²⁰ cm⁻²) and halo mass function from our simulations and use the formula in Sheth et al. (2001) to compute the halo bias. We obtain values of the DLA bias equal to 1.7 at z = 2 and 2 and z = 3. Considering that the DLA bias follows a linear relation between z = 2 and z = 3, we obtain b_DLA(z = 2.3) = 1.8, in agreement with the latest observations by Pérez-Ràfols et al. (2018): b_DLAs = 1.99 ± 0.11. We have repeated the above calculations using our fit for the DLA cross-section, taking the halo mass function from Sheth & Tormen (2002) or Crocce et al. (2010), and find that our results barely change.

We believe that the above calculation should be considered as a lower bound. In other words, the halo bias may be underestimated when calculated using Sheth et al. (2001). The reason is that we obtain a value of the H i bias (see Section 13), computed without any assumption, of 2 at z = 2 and 2.56 at z = 3, i.e., b_{H i}(z = 2.3) = 2.17. Following the theoretical arguments in Castorina & Villaescusa-Navarro (2017), it is reasonable to expect that b_{H i} ≃ b_DLAs. We thus conclude that both estimations of the DLA bias are in agreement with observations.

In Section 4, we showed that nearly all of the H i at redshifts z ≤ 5 reside within halos. Thus, the elements needed to describe the abundance and spatial distribution of H i in real space through H i halo models are the halo H i mass function and the H i density profiles. To model the distribution of H i in redshift space, an additional ingredient is required: the velocity distribution of H i inside halos. This quantity can be used in both H i halo models and H i HOD models. The accuracy that can be achieved with the former may not be high, due to the limitations of the formalism itself. On the other hand, H i HOD, i.e., painting H i on top of dark matter halos from either N-body or fast numerical simulations like COLA (Tassev et al. 2013), can produce highly accurate results. Hence, we examine the velocity distribution of H i inside halos, beginning with the H i bulk velocity in this section and continuing with the H i velocity dispersion in the next section, and, in both cases, comparing with the results for all matter.

For each dark matter halo in the simulation, we have computed the H i bulk velocity as

$$\begin{eqnarray}&&{{\boldsymbol{V}}}_{{\rm{H}}{\rm{I}}}=\displaystyle \frac{{\displaystyle \sum }_{i}{M}_{{\rm{H}}{\rm{I}},{\rm{i}}}{{\boldsymbol{V}}}_{{\rm{H}}{\rm{I}},{\rm{i}}}}{{\displaystyle \sum }_{i}{M}_{{\rm{H}}{\rm{I}},{\rm{i}}}},\end{eqnarray} \tag{ 23 }$$

where the sum runs over all gas cells belonging to the halo, and M_{H i,i} and ${{\boldsymbol{V}}}_{{\rm{H}}{\rm{I}},{\rm{i}}}$ are the H i mass and peculiar velocity of cell i, respectively. The peculiar velocity of halos, ${{\boldsymbol{V}}}_{h}$, is computed in a similar manner, but summing over all resolution elements in the halo (gas, CDM, stars, and black holes) and weighting their velocities by their corresponding masses.

Here, we examine (1) whether the peculiar velocity of the H i points in the same direction as the halo peculiar velocity, and (2) whether the modulus of the H i peculiar velocity is the same as that of the halo peculiar velocity. The first point is addressed by computing the angle between the peculiar velocities of H i and the halo from

$$\begin{eqnarray}&&\cos (\alpha )=\displaystyle \frac{{{\boldsymbol{V}}}_{{\rm{H}}{\rm{I}}}\cdot {{\boldsymbol{V}}}_{h}}{| {{\boldsymbol{V}}}_{{\rm{H}}{\rm{I}}}| | {{\boldsymbol{V}}}_{h}| }\end{eqnarray} \tag{ 24 }$$

for each halo in the simulation. We do not consider halos with total H i masses below 10⁵ h⁻¹ M_⊙, since we expect the H i peculiar velocities of those halos to be uncorrelated with the halo peculiar velocity. For example, the H i in such halos mass can be from a single cell that is partially self-shielded and not bound to the halo. Moreover, in the limit where the H i mass is close to zero, the H i velocity dispersion is not well defined. Thus, in order to avoid such circumstances, we adopt the above threshold, which corresponds to the mass of ≃1/5 of a completely self-shielded gas cell. However, we find that this threshold does not have a significant impact on our results. Choosing a different value hardly changes our results, with the only consequence being that the scatter of very small halos is affected. We then take narrow bins in halo mass and compute the mean value of $\cos (\alpha )$ and its standard deviation. The results are shown in the upper panels of Figure 12.

Zoom In Zoom Out Reset image size

For small halos, M ≲ 10¹² h⁻¹ M_⊙, $\cos (\alpha )\simeq 1$, indicating that the H i and halo peculiar velocities are aligned. This is expected because the H i is mainly located in the inner regions in low-mass halos, which usually traces well the peculiar velocity of the halo. For smaller halos, the value of $\cos (\alpha )$ deviates from 1, with increased scatter. This happens for halo masses below the cutoff scale, $\sim {M}_{\min }$. In at least some cases, this is likely due to halos acquiring H i through an unusual mechanism, e.g., by passing through an H i-rich filament, so that the H i bulk velocity will not be correlated with the halo peculiar velocity. On the other hand, we find significant misalignments between the H i and halo peculiar velocities for the most massive halos at any redshift. This is because the H i content of these halos is largely contributed by satellites, whose peculiar motions do not necessarily trace those of the halo. We return to this point below.

Further, in the bottom panels of Figure 12, we show the average and standard deviation of the ratio between the moduli of the H i and halo peculiar velocities, $| {{\boldsymbol{V}}}_{{\rm{H}}{\rm{I}}}| /| {{\boldsymbol{V}}}_{h}|$. This quantity is again calculated in narrow bins in halo mass, for all halos with H i mass larger than 10⁵ h⁻¹ M_⊙.

For small halos, the moduli of the H i and halo peculiar velocities are essentially the same. For halos with masses below ∼M_min, the modulus ratio can be larger than 1, and its scatter increases. This is for the same reason as above: the H i content of some of those halos may not be bound to the halos and are instead part of a filament. For massive halos, the modulus of the H i peculiar velocity can be much larger than that of the halo peculiar velocity. As earlier, this is because the H i peculiar velocity is dominated by the H i in satellites, whose peculiar velocities do not perfectly trace the halo peculiar velocity.

To corroborate the assertion that the peculiar velocities of satellites do not trace the halo peculiar velocity in either modulus or direction, we have performed the following test. We compute the peculiar velocities of halo satellites and compared their mean, weighted by the total mass of each satellite, against the peculiar velocity of their host halo. The velocities of the satellites do not have the same modulus or direction as those of the host halo, showing similar trends to those for H i, with differences increasing with halo mass.

Thus, for small halos, where most of the H i is in the central galaxy, the H i bulk velocity traces the halo peculiar velocity well, in both modulus and direction. On the other hand, the contribution of satellites to the total H i mass in halos increases with mass, and since the bulk velocities of satellites do not trace the halo peculiar velocity, the H i bulk velocities will diverge, in modulus and direction, from the halo peculiar velocity, with differences increasing with halo mass.

For each halo in the simulation, we have computed the 3D velocity dispersion of its H i from

$$\begin{eqnarray}&&{\sigma }_{{\rm{H}}\,{\rm{I}}}^{2}=\displaystyle \frac{{\displaystyle \sum }_{i}{M}_{{\rm{H}}{\rm{I}},{\rm{i}}}| {{\boldsymbol{V}}}_{{\rm{H}}{\rm{I}},{\rm{i}}}-{{\boldsymbol{V}}}_{{\rm{H}}{\rm{I}}}{| }^{2}}{{\displaystyle \sum }_{i}{M}_{{\rm{H}}{\rm{I}},{\rm{i}}}},\end{eqnarray} \tag{ 25 }$$

where the sums run over all gas cells belonging to the halo. ${M}_{{\rm{H}}{\rm{I}},{\rm{i}}}$ and ${{\boldsymbol{V}}}_{{\rm{H}}{\rm{I}},{\rm{i}}}$ are the H i mass and velocity of gas cell i, and ${{\boldsymbol{V}}}_{{\rm{H}}{\rm{I}}}$ is the H i bulk velocity, computed as in Equation (24).

We then take narrow bins in halo mass and compute the mean and standard deviation of the H i velocity dispersion. Figure 13 shows the results, as well as a comparison with matter, whose properties are calculated analogously, but considering all mass elements within halos, i.e., gas, CDM, stars, and black holes.

Zoom In Zoom Out Reset image size

As expected, the velocity dispersion of both H i and CDM increases with halo mass, independent of redshift. The results can be represented by a simple power law,

$$\begin{eqnarray}&&{\sigma }_{v}(M)={\sigma }_{10}{\left(\displaystyle \frac{M}{{10}^{10}{h}^{-1}{M}_{\odot }}\right)}^{\alpha },\end{eqnarray} \tag{ 26 }$$

where σ₁₀ and α are free parameters with best-fit values provided in Table 4.

The mean H i velocity dispersions are always equal to or smaller than the matter velocity dispersions. For large halo masses, both exhibit the same amplitude, but for low-mass halos, the H i velocity dispersion is less than that of matter. The typical halo masses where the velocity dispersions diverge is around 10¹² h⁻¹ M_⊙, with higher redshifts exhibiting departures at larger masses. This behavior is embedded in the slope of the σ_{H i}(M) relation, the value of which, α, is larger than that of σ_m(M) at all redshifts, and more so at higher redshifts.

For very small halos, and particularly at low redshift, the velocity dispersion of H i is much smaller than that of matter. This is a consequence of several factors. First, σ_m(M) is artificially high, as the relation flattens out toward low masses. By comparing to a version of Figure 13 (not shown) generated from a lower resolution analog of the same cosmological volume (TNG100-2; see Nelson et al. 2018), we conclude that this is due to finite numerical resolution, driven by particles in the outskirts of the halos, which does not apply to the H i, which is centrally concentrated. In addition, we have examined a few individual low-mass halos and found that in some cases, the H i arises from just a few cells, or even just a single one. In those cases, the H i bulk velocity will be set by these few cells and the H i velocity dispersion will be artificially suppressed due to sampling, again from a finite resolution. As we move to more massive halos, the contribution of H i from satellites increases, and those satellites trace the underlying matter distribution more closely.

The scatter in the H i velocity dispersion for very low-mass halos is typically much smaller than the scatter in the matter velocity dispersion. One reason for this is that the H i velocity dispersion has been computed only for halos with total H i masses above 10⁵ h⁻¹ M_⊙. Without such a threshold, the scatter in the H i velocity dispersion would be much larger. That is because if the H i mass in a halo is very low, it will often not be bound to the halo, e.g., the halo is crossing a filament that hosts a small amount of H i. In that case, H i velocity dispersions can be large. For example, several highly ionized unbound gas cells can produce a large, unphysical, velocity dispersion.

The relation between H i and matter is given, to linear order, by δ_{H i} = b_{H i}δ_m + , where is the stochasticity. Below, we examine whether or not this relation reproduces our results in real space and the amplitude of the stochasticity.

As for the density pdf's (see Section 8), we compute the density fields of H i and matter on a grid with 2048³ cells employing the CIC mass-assignment scheme. We then compute the overdensity of each field and smoothed them with a top-hat filter of radius 1 or 5 h⁻¹ Mpc. Next, we randomly select a subset of cells and make a scatter plot between the overdensities of H i and matter for each chosen cell. The results are shown in Figures 14 (R = 1 h⁻¹ Mpc) and 15 (R = 5 h⁻¹ Mpc).

Zoom In Zoom Out Reset image size

For R = 1 h⁻¹ Mpc, two trends can be distinguished. The Lyα forest shows up as cells with matter overdensities below the mean and very low H i overdensities because the gas there is mostly ionized. For large matter overdensities, the H i within halos is self-shielded. The density ${\rho }_{{\rm{m}}}/{\bar{\rho }}_{{\rm{m}}}\,\simeq \,0.4$ marks a transition from one regime to the other, indicating that H i self-shielding does not take place in lower matter overdensities. In all cases, the H i overdensity increases with matter overdensity.

At higher redshifts, the range occupied by matter and H i overdensities is smaller. As the universe becomes more homogeneous, fluctuations are smaller. This behavior can also be seen in the pdf's of Figure 8. The scatter in the overdensity relations also decreases toward higher redshift.

The dashed black lines show the predictions from linear theory, δ_{H i} = b_{H i}δ_m, where b_{H i} is the linear H i bias measured from the simulation (see Section 13 and Table 5). As expected, linear theory is not accurate in this regime because the bias is not linear on the smoothing scale considered. An exception is for z = 1, where linear H i bias reproduces the results reasonably well. This is because the H i bias is relatively flat at that redshift (see Figure 16).

Table 5.  Values of the H i Bias and H i Shot Noise from the Simulation at Different Redshifts

z	0	1	2	3	4	5
bH i	0.84	1.49	2.03	2.56	2.82	3.18
${P}_{{\rm{H}}\,{\rm{I}}}^{\mathrm{SN}}$	104	124	65	39	14	7
[(h−1 Mpc)3]

Download table as:  ASCIITypeset image

As we move to a larger smoothing radius, the morphology of the results changes, as can be seen in Figure 15, where R = 5 h⁻¹ Mpc. Now, the H i and matter overdensities extend over a smaller range, because the smoothing is over a larger scale, making the field more homogeneous. In addition, the Lyα forest is no longer visible because in the neighborhood of the highly ionized H i in filaments, i.e., the Lyα forest, there will always be some halo within 5 h⁻¹ Mpc that contains self-shielded H i gas.

Zoom In Zoom Out Reset image size

As above, we find that H i overdensities increase with matter overdensities at all redshifts. However, at high redshift, the slope of the relation becomes more pronounced. This behavior can be partly explained by linear H i bias, shown as dashed black lines. As with R = 1 h⁻¹ Mpc, linear bias can explain the results relatively well at z = 1. At other redshifts, linear bias is more accurate than for smaller smoothing scales, as expected, but the agreement is not good for both large and small matter overdensities. Again, the scatter reduces with redshift.

We now examine different aspects of H i clustering in detail. In this section, we focus on the amplitude and shape of the H i bias.

The relation between the clustering of H i and that of dark matter involves the H i bias through ${P}_{{\rm{H}}{\rm{I}}}(k)={b}_{{\rm{H}}\,{\rm{I}}}^{2}(k){P}_{{\rm{m}}}(k)$. The matter power spectrum, the quantity that contains information on the values of the cosmological parameters, can thus be inferred only if bias is understood (see Pénin et al. 2018 for a detailed discussion on the H i scale-dependence bias). On linear scales, the H i bias is constant, but on small scales, we expect to see scale dependence. It is important to determine the scales on which the H i bias is dependent and whether or not analytic models can reproduce that behavior.

We have computed the H i and matter auto-power spectrum and the H i−matter cross-power spectrum of the simulation at redshifts 0, 1, 2, 3, 4, and 5. The H i bias is then obtained using two different definitions: ${b}_{{\rm{H}}{\rm{I}}}(k)=\sqrt{{P}_{{\rm{H}}{\rm{I}}}(k)/{P}_{{\rm{m}}}(k)}$ and b_{H i}(k) = P_{H i−m}(k)/P_m(k). While the latter is "preferred," as it does not suffer from stochasticity, the former is closer to observations. The results are shown in Figure 16.

Zoom In Zoom Out Reset image size

The amplitude of the H i bias on large scales increases with redshift, from ≃0.85 at z = 0 to ≃3.20 at z = 5. On the largest scales that can be probed with TNG100, the amplitude of the H i bias is independent of the method used to estimate it. We assume that the values on large scales are equal to the linear H i bias, but note that there could be small corrections to those because of box size. The linear H i bias at different redshifts is given in Table 5. These values can be reproduced from the halo H i mass function as

$$\begin{eqnarray}&&{b}_{{\rm{H}}{\rm{I}}}(z)=\displaystyle \frac{{\int }_{0}^{\infty }n(M,z)b(M,z){M}_{{\rm{H}}{\rm{I}}}(M,z){dM}}{{\int }_{0}^{\infty }n(M,z){M}_{{\rm{H}}{\rm{I}}}(M,z){dM}},\end{eqnarray} \tag{ 27 }$$

and therefore, the amplitude of the H i bias is sensitive to the astrophysical parameters α and M_min (see Section 5). Note, however, that the agreement between the above expression and the simulation results is not perfect because, among other things, our models for the halo mass function and halo bias do not include corrections for baryonic effects.

At z = 0, the H i bias exhibits a scale dependence even on the largest scales we can probe, due to the fact that the matter power spectrum at the scales probed by TNG100 is not in the linear regime at such low redshift. The dip in the H i bias at k ≃ 1 h Mpc⁻¹, which has also been found in observations (Anderson et al. 2018), is interesting. At z = 1, the bias remains almost constant down to rather small scales, k ≃ 1 h Mpc⁻¹. These trends agree with the findings of Springel et al. (2018), who studied galaxy bias for different galaxy populations at different redshifts. At high redshifts, z ≥ 2, the H i bias already exhibits a dependence on scale at k = 0.3 h Mpc⁻¹, even though these scales are close to linear at those redshifts. Our results are also in qualitative agreement with those of Sarkar et al. (2016a), who studied the H i bias by painting H i on top of dark matter halos. The scale dependence of the bias is not necessarily a bad thing, as long we can use perturbative methods to predict the shape of the H i power spectrum. For this purpose, we have compared the measurements of the H i power spectrum in TNG100 to analytical calculations using Lagrangian Perturbation Theory (LPT). The first-order LPT solution is the well-known Zel'dovich approximation (ZA) (Zel'dovich 1970; White 2014), for which we can simply write

$$\begin{eqnarray}&&{P}_{{\rm{H}}{\rm{I}}}={b}_{{\rm{H}}\,{\rm{I}}}^{2}{P}_{\mathrm{ZA}}(k)+N,\end{eqnarray} \tag{ 28 }$$

and then fit for the two free parameters in the above equation. The constant piece takes care of the shot noise and any other term that is scale independent and uncorrelated with the H i field and therefore can be treated as noise in a cosmological analysis (Seljak & Vlah 2015). Given the small volume of TNG100, a perturbative analysis makes sense only at high redshifts, where linear and mildly nonlinear modes are contained in the box, thus we restrict the comparison of Equation (28) with the measurements in the simulation to z ≥ 2. The upper panel in Figure 17 shows the measurements of the H i power spectrum at different redshifts, using the same color scheme as in the previous figures. The points with error bars have been shifted horizontally to avoid overlap and facilitate the visual comparison with the theoretical models. The dashed lines display the fit to Equation (28) including all the modes up to k_max = 1 h Mpc⁻¹. The fit is quite accurate, despite its simple functional form, and it confirms the H i distribution as an ideal tracer for cosmological studies. The continuous lines show the next to leading order, i.e., one-loop, calculation in LPT (Vlah et al. 2016; Modi et al. 2017), which includes an improved treatment of nonlinearities in the matter fields as well as several nonlinear bias parameters. Up to the scale we include in the fit, there is no difference between the two approaches, with the one-loop calculation also working on smaller scales not included in the analysis. The fact that the ZA works so well in describing the simulation measurements could vastly simplify the cosmological analysis and interpretation of 21 cm surveys observing at high redshift. For instance, interferometric surveys with a large instantaneous field of view like CHIME will be forced to include all of the complications arising from the curved sky, which are very easily handled in the ZA (Castorina & White 2018a, 2018b).

Zoom In Zoom Out Reset image size

It is well known that the clustering of halos depends primarily on mass. However, mass is not the only variable that determines halo clustering; there is also a dependence on halo age (Sheth & Tormen 2004; Gao et al. 2005), concentration (Wechsler et al. 2006), subhalo abundance (Wechsler et al. 2006; Gao & White 2007), halo shape (Faltenbacher & White 2010), spin (Gao & White 2007), and environment (Xu & Zheng 2018; Han et al. 2018; Salcedo et al. 2018). Here, we identify a new secondary bias, originating from the H i content of halos.

Differently from the previous quantities, which are properties of the dark matter halos, the H i content is more related to the properties of the galaxies inside a halo, e.g., whether galaxies are red or blue. Thus, a study of this kind can only be carried out using hydrodynamic simulations.

To study this issue, we apply the following procedure. First, all halos whose total mass is within a relatively narrow mass bin are selected. The H i mass inside each of those halos, and the median value are determined. Next, the halos are split into two categories: H i rich and H i-poor, depending on whether the H i content of a particular halo is above or below the median, respectively. Finally, we compute the power spectrum of 1) all halos, 2) H i-rich halos, and 3) H i-poor halos.

The results are shown in Figure 18 at redshifts 0, 1, 2, 3, and 4 and for three different mass bins, M ∈ [10⁸–10⁹] h⁻¹ M_⊙, M ∈ [10⁹–10¹⁰] h⁻¹ M_⊙, and M ∈ [10¹⁰–10¹¹] h⁻¹ M_⊙. The black lines indicate the power spectrum of all halos, while the blue and red lines represent the power spectra of the H i-rich and H i-poor halos, respectively.

Zoom In Zoom Out Reset image size

Going toward smaller scales, the amplitudes of the different power spectra first flatten and then rise back up. This happens (1) because we approach the shot-noise limit and (2) due to aliasing. The shot-noise level of the H i-rich/-poor halos is different from that of all halos, as the latter contain, by definition, twice as many halos than the former. Thus, on small scales, the amplitude of the H i-rich and H i-poor halos is expected to be the same but higher by a factor of 2 than that of all halos, as is indeed seen.

For all redshifts and mass intervals considered, the clustering of H i-rich galaxies is different from that of H i-poor ones, showing that halo clustering depends not only on mass but on H i content as well. The difference in the clustering of H i-poor and H i-rich halos decreases, in general, with halo mass. At z = 0, the amplitude of the power spectrum of the H i-rich and H i-poor halos can be almost one order of magnitude different for halos in the [10⁹–10¹⁰] h⁻¹ M_⊙ mass bin. The largest differences are seen for halos with masses around or below M_min at that particular redshift, namely, around the mass scale where the H i content starts being exponentially suppressed (see Section 5).

At z = 0, and for the mass bin intervals considered here, the H i-poor halos are more strongly clustered than the H i rich halos. On the other hand, at high redshift, the situation is the opposite, and H i-rich halos are more strongly clustered than H i-poor halos. At z = 1, we find that depending on the halo mass considered, H i-rich halos can be more or less clustered than H i-poor halos. We note that Guo et al. (2017) found that the clustering of H i galaxies from the ALFALFA survey at z ∼ 0 shows a bias systematically lower than mass-selected halos, in agreement with our results for H i-rich halos.

Although the halo mass bins are fairly narrow, is not unreasonable to suspect that the most massive halos will have larger H i masses, and therefore, this could introduce some natural splitting that arises just from halo mass and not from H i secondary bias. In order to test this, we split the halos according to their median total halo mass and repeated the above analysis. We find that the clustering of the two samples is almost indistinguishable, ruling out the possibility that halo mass is affecting our results.

At low redshift, small halos near big ones are more likely to be stripped of their gas content. Thus, H i-poor halos should be more strongly clustered than H i-rich halos, as we find. This possibility has been recently suggested to explain the secondary bias that arises from several halo properties in Salcedo et al. (2018) and Han et al. (2018).

On the other hand, at high redshift, gas stripping by nearby halo neighbors should be less effective, as the largest halos are not yet very massive and there has been less physical time for these processes to operate. We speculate that at high redshift, regions around massive halos are richer in H i than other regions. For example, in regions with higher density, we would expect the filaments to be slightly more dense and therefore will host more H i. Halos connected by those filaments may thus become H i rich. This naive picture can be seen in Figure 27, where at high redshift, the filaments in the denser regions host more H i than in less dense regions.

An important consideration in any cosmological survey is shot noise, as its amplitude determines the maximum scale where cosmological information can be extracted from (see Equation (1)). However, it can also be used to learn about the galaxy population hosting the H i (Wolz et al. 2017, 2018). The purpose of this section is to quantify the amplitude of the H i shot noise from our simulations.

We now illustrate why computing the H i shot noise is slightly more complicated than determining the shot noise for other tracers, such as halos, where the value of the shot noise is simply given by the amplitude of the power spectrum on small scales.

The solid lines of Figure 19 show the H i power spectrum at redshifts 0 and 4. It can be seen that those power spectra receive contributions from both the one- and two-halo terms, i.e., the power spectrum on very small scales does not become constant, in contrast with the halo power spectrum. This happens simply because there is structure in H i inside halos and galaxies (see Figure 6).

Zoom In Zoom Out Reset image size

In order to isolate the contribution of the H i shot noise, i.e., to avoid the one-halo term contribution, we do the following. We compute the total amount of H i inside every halo in the simulation and place that H i mass in the halo center. We then compute the H i power spectrum of that configuration. Since in that case there is no H i structure inside halos, there is no one-halo term, and the amplitude of the H i power spectrum on small scales is just the H i shot noise. The dashed lines in Figure 19 display the results.

It can be seen that on large scales, the amplitude of the H i power spectrum of the two configurations is essentially identical. The small differences at z = 4 arise from the H i that is outside of halos (see Figure 3), whose contribution is not accounted for with our procedure. That contribution should, however, have a negligible impact on the amplitude of the H i shot noise. On small scales, the lack of the one-halo term when we artificially place the H i at the halo center makes it possible to isolate the value of the H i shot noise. We determine the value of the H i shot noise by averaging the amplitude of the H i power spectrum on scales k ∈ [20–30] h Mpc⁻¹ and show the results in Table 5.

The H i shot noise can also be computed using the halo model framework (assuming all H i are in halos) as (Castorina & Villaescusa-Navarro 2017)

$$\begin{eqnarray}\begin{array}{rcl}{P}_{{\rm{H}}\,{\rm{I}}}^{\mathrm{SN}}(z) & = & \mathop{\mathrm{lim}}\limits_{k\to 0}{P}_{1{\rm{h}},{\rm{H}}{\rm{}}{\rm{i}}}(k,z)=\displaystyle \frac{1}{{[{\rho }_{{\rm{c}}}^{0}{{\rm{\Omega }}}_{{\rm{H}}{\rm{I}}}(z)]}^{2}}\\ & & \times \,{\displaystyle \int }_{0}^{\infty }n(M,z){M}_{{\rm{H}}\,{\rm{I}}}^{2}(M,z){dM}\\ & = & \displaystyle \frac{{\displaystyle \int }_{0}^{\infty }n(M,z){M}_{{\rm{H}}\,{\rm{I}}}^{2}(M,z){dM}}{{\left[{\displaystyle \int }_{0}^{\infty }n(M,z){M}_{{\rm{H}}{\rm{I}}}(M,z){dM}\right]}^{2}}.\end{array}\end{eqnarray} \tag{ 29 }$$

If we use the Sheth & Tormen (ST) formula (Sheth & Tormen 2002) for the halo mass function in this expression²⁹ and our fitting formula for M_{H i}(M, z) (see Table 1), we obtain values for the H i shot noise in agreement with those measured directly in the simulation.

In order to validate our results, we have estimated the H i shot noise by measuring the stochasticity between the H i and the matter fields (Seljak et al. 2009)

$$\begin{eqnarray}&&{P}_{{\rm{H}}{\rm{I}}}(k)-\displaystyle \frac{{P}_{{\rm{H}}\,{\rm{I}}-{\rm{m}}}^{2}(k)}{{P}_{{\rm{m}}}(k)},\end{eqnarray} \tag{ 30 }$$

which at low enough k should correspond to the uncorrelated part of the H i power spectrum, i.e., P_SN. However, the shot-noise amplitudes we obtain in this case are much larger than those found using our fiducial procedure. By repeating this approach for halos, we also find that we cannot recover the standard $1/\bar{n}$ result. Our findings are in agreement with those of other studies, see Hand et al. (2017) and reference therein, and point toward some lack of understanding of the noise properties in low-mass halos even in gravity-only simulations. For better comparison with previous work, which always considered Poissonian shot noise, we assume P_SN is given by the fiducial halo model procedure, and defer a study of this unexpected disagreement to future work.

We find that the H i shot noise is low at all redshifts. This is in broad agreement with the analytical work of Castorina & Villaescusa-Navarro (2017), whose results at high redshift had, however, a large theoretical error. Our shot-noise values are slightly lower than those in that paper, mostly because the values of M_min at high z that we find here are lower than those considered in Castorina & Villaescusa-Navarro (2017).

The low values of shot noise we find have important implications. First, 21 cm intensity mapping experiments that aim at measuring the BAO peak (Obuljen et al. 2018a) position such as CHIME, OOTY, BINGO, HIRAX, or SKA will barely be affected by shot noise. We illustrate this in Figure 20, where we plot, in blue,

$$\begin{eqnarray}&&{{nP}}_{0.2}(z)=\displaystyle \frac{{P}_{{\rm{H}}{\rm{I}}}(k=0.2\,h\,{\mathrm{Mpc}}^{-1},z)}{{P}_{\mathrm{SN}}(z)}\end{eqnarray} \tag{ 31 }$$

as a function of redshift. We obtain values of nP_0.2 above ≃15 at all redshifts, showing that shot-noise contamination on BAO scales is minimal. This is a great advantage over traditional methods such as galaxy surveys.

Zoom In Zoom Out Reset image size

Second, the shot-noise levels are very low at high redshift. This means that shot noise does not erase the cosmological information on small scales. The red line shows the value of nP at k = 0.5 h Mpc⁻¹, which is still much larger than one and implies that these modes can be in principle measured in the cosmic variance limit. With accurate theory predictions such as the one described in the previous section, this information can be retrieved and used to reduce errors in the values of the cosmological parameters.

Twenty-one centimeter intensity mapping observations probe the spatial distribution of cosmic neutral hydrogen in redshift space, not in real space. Thus, it is of utmost importance to understand the impact of peculiar velocities on the clustering of H i. In this section, we study the clustering of neutral hydrogen in redshift space.

Here, we make use of the plane-parallel approximation to displace the positions of particles and Voronoi cells from real (${\boldsymbol{x}}$) to redshift space (${\boldsymbol{s}}$) through

$$\begin{eqnarray}&&{\boldsymbol{s}}={\boldsymbol{x}}+\displaystyle \frac{1+z}{H(z)}{{\boldsymbol{v}}}_{| | }({\boldsymbol{r}}),\end{eqnarray} \tag{ 32 }$$

where ${{\boldsymbol{v}}}_{| | }({\boldsymbol{r}})$ is the peculiar velocity of the particle/cell along the line of sight. We use the three Cartesian axes as different lines of sight. Our results represent the mean over the three axes.

We have computed the clustering of H i in redshift space and show the results in Figure 21. Although we have computed the H i monopoles, quadrupoles, and hexadecapoles, the latter two are too noisy so we restrict our analysis to the monopoles.

Zoom In Zoom Out Reset image size

We show with red/blue lines the monopoles in real/redshift space in this figure. The two main physical processes governing the effect of peculiar velocities can clearly be seen. On large scales, the clustering of H i in redshift space is enhanced due to the Kaiser effect (Kaiser 1987). On small scales, the peculiar velocities of H i inside halos give rise to the Fingers of God, suppressing the amplitude of the H i power spectrum.

The bottom part of each panel of Figure 21 shows the ratio between the monopoles in redshift and real space as a solid black line. The black dashed line displays the prediction of linear theory for that ratio, i.e.,

$$\begin{eqnarray}&&\displaystyle \frac{{P}_{0}^{s}(k)}{{P}_{0}^{r}(k)}=1+\displaystyle \frac{2}{3}\beta +\displaystyle \frac{1}{5}{\beta }^{2},\end{eqnarray} \tag{ 33 }$$

where β = f/b_{H i}, with $f\simeq {{\rm{\Omega }}}_{{\rm{m}}}^{0.545}(z)$ being the linear growth rate. We have estimated β using the values of the linear H i bias from Table 5.

As expected, linear theory is not able to describe our results at low redshift. This is because the scales we probe in TNG100 are too small for linear theory to hold. On the other hand, we find that the Kaiser factor can reasonably well explain the monopole ratio down to k ≃ 0.3, 0.5, and 1.0 h Mpc⁻¹ at redshifts 3, 4, and 5, respectively.

Given the results in Section 13, our findings are not surprising, and we plan to compare the measurements in TNG100 to redshift-space LPT analytical predictions in upcoming work.

As an alternate way to visualize the consequences of redshift-space distortions, we show the two-dimensional power spectrum of cosmic H i in the top two rows of Figure 22, at redshifts from 0 to 5. The Kaiser effect manifests itself as a squeezing of isopower contours along the perpendicular direction on large scales. It is apparent down to relatively small scales at high redshift, as expected.

Zoom In Zoom Out Reset image size

On small scales, the Fingers of God arise as isopower contours propagating farther in the perpendicular than the radial direction. At low redshift and on small scales, we find that isopower contours exhibit a very weak dependence on the perpendicular direction, i.e., ${P}_{{\rm{H}}\,{\rm{I}}}^{s}({k}_{| | },{k}_{\perp })\simeq {P}_{{\rm{H}}\,{\rm{I}}}^{s}({k}_{| | })$.

We show the results for the 2D matter power spectrum in the bottom rows of Figure 22. It can be seen that (1) the Kaiser effect is visible down to smaller scales than in the H i field at high redshift and (2) the magnitude of the Fingers of God is lower in the matter field than in the H i field.

The halo model can be used to understand why the magnitude of the Fingers of God is higher in H i than in matter. Each dark matter halo will show up in redshift space not as a sphere but as an ellipsoid, due to internal peculiar velocities. The eccentricity of those ellipsoids will depend primarily on halo mass: the velocity dispersion in large halos will be higher than that in small halos, so their eccentricity will be larger. The amplitude of the matter power spectrum on small scales will be dominated by small halos, whose velocity dispersions are not large.

On the other hand, we know for H i that since there is a cutoff in the halo H i mass function, halos below a certain mass will not contribute to the H i power spectrum. Thus, it is expected that the amplitude of the Fingers of God will be larger in the H i field than in the matter field because (1) the velocity dispersions of H i and matter/CDM are similar (see Section 11), and (2) in the H i field, we do not have the contribution of small halos that dominates the amplitude of the power on small scales.

In order to corroborate this hypothesis, we have taken all halos with masses above 10¹⁰ h⁻¹ M_⊙, and we have computed the power spectrum of the matter inside them. In that case, we observe very similar results to those from the H i, i.e., the amplitude of the Fingers of God of that particle distribution is much larger than the one for all matter.

The impact of redshift-space distortions on the H i power spectrum has recently been studied in Sarkar & Bharadwaj (2018) using N-body simulations. Following their work, we have tried to model the 2D H i power spectrum using the following phenomenological expression

$$\begin{eqnarray}&&{P}_{{\rm{H}}\,{\rm{I}}}^{s}(k,\mu )={(1+\beta {\mu }^{2})}^{2}{P}_{{\rm{H}}\,{\rm{I}}}^{r}(k)\displaystyle \frac{1}{{\left(1+\tfrac{1}{2}{k}^{2}{\mu }^{2}{\sigma }_{p}^{2}\right)}^{2}},\end{eqnarray} \tag{ 34 }$$

where β = f/b_{H i}, ${P}_{{\rm{H}}\,{\rm{I}}}^{r}(k)$ is the fully nonlinear H i power spectrum, and σ_p is a phenomenological parameter that accounts for the Fingers-of-God effect. The reason for using this expression is that even if the H i bias is nonlinear, we will recover the Kaiser factor when computing the monopole ratio (see Figure 21) at high redshift, where the amplitude of the Fingers of God is small.

By fitting the above expression to our 2D H i power spectrum down to k = 1 h Mpc⁻¹ and assuming Gaussian errors, we find that this approach works reasonably well, χ²/d.o.f ≃ [1.5–2.2], with σ_p = 1.73, 2.09, 1.37, 0.93, 0.34, and 0 at redshifts 0, 1, 2, 3, 4, and 5, respectively. However, by computing the monopole from the above expression

$$\begin{eqnarray}&&{P}_{{\rm{H}}\,{\rm{I}}}^{0}(k)=\displaystyle \frac{1}{2}{\int }_{-1}^{1}{P}_{{\rm{H}}\,{\rm{I}}}^{s}(k,\mu )d\mu \end{eqnarray} \tag{ 35 }$$

and comparing with our measurements, we do not find good agreement.

We have also repeated the above exercise for the models considered in Sarkar & Bharadwaj (2018), concluding that none of those represent a good description of our results. We leave a more detailed analysis of this issue to future work.

The H i properties studied in this paper can be used to generate mock 21 cm maps. In this section, we study (1) whether less computationally expensive simulations can be used to create 21 cm maps and (2) the importance of accounting for the one-halo term when making mocks.

The right column of Figure 23 shows 21 cm maps created from the spatial distribution of H i in the TNG100 simulation. From top to bottom, these maps show the 21 cm map in real space with 3' angular resolution (top), 21 cm map in redshift space with 3' angular resolution (top middle), 21 cm map in real space with 03 angular resolution (bottom middle), and 21 cm map in redshift space with 03 angular resolution (bottom). All of those maps are centered at a frequency of 710 MHz (z = 1) and have a bandwidth of 1 MHz (≃5 h⁻¹ Mpc).

Zoom In Zoom Out Reset image size

The procedure used to create these maps is as follows. First, the H i density field is computed by assigning H i masses of gas cells (either in real or redshift space) to a grid of 2048³ cells using the nearest-grid-point (NGP) mass-assignment scheme. Then, we select a slice of the H i density field grid whose width is taken to reproduce the desired frequency bandwidth. Next, that slice is projected onto a two-dimensional grid, and H i densities are transformed to brightness temperatures through

$$\begin{eqnarray}&&{T}_{b}({\boldsymbol{x}})=189\,h\left(\displaystyle \frac{{H}_{0}{(1.0+z)}^{2}}{H(z)}\right)\displaystyle \frac{{\rho }_{{\rm{H}}{\rm{I}}}({\boldsymbol{x}})}{{\rho }_{{\rm{c}}}}\,\mathrm{mK}.\end{eqnarray} \tag{ 36 }$$

Finally, we convolve that grid with a Gaussian filter of radius R = θχ, where θ is the desired angular resolution and χ is the comoving distance to redshift z (or frequency ν).

The 21 cm maps on the left column of Figure 23 have been created from an N-body simulation that shares the same initial conditions as TNG100, i.e., TNG100-1-DM, but whose computational cost is over an order of magnitude lower. In that simulation, we have placed an H i mass given by M_{H i}(M, z) from Table 1 on the center of each dark matter halo. Then, we have followed the procedure outlined above to create the 21 cm maps. To make 21 cm maps in redshift space, we displace the H i mass that we put in the center of each halo according to the peculiar velocity of that halo.

It can be seen that the qualitative agreement between the maps from the full hydrodynamic simulation and the cheaper N-body one is very good. In Figure 24, we quantify this visual agreement by computing the H i power spectrum of the hydrodynamic and N-body simulations in real space at different redshifts.

Zoom In Zoom Out Reset image size

We find that the amplitude of the two H i power spectra can be significantly different, between 10% and 40%. We attribute these differences to the inaccuracies of our M_{H i}(M, z) function (e.g., we do not require our fit to reproduce Ω_{H i}(z)), to the effects of baryons on the halo mass function and halo clustering, which are ignored in this exercise, and to the omission of the one-halo term at very high k. For the latter, we emphasize that in this section we do not consider the spatial distribution of H i inside halos, which we leave to future work. However, once the difference in amplitude is taken into account, we find that shapes differ by only ≃5% from the largest scales down to k = 1 h Mpc⁻¹ at all redshifts. This is rather remarkable considering that the one-halo term was not accounted for. This exercise demonstrates that populating dark matter halos from a computationally cheap N-body simulation with H i can yield results that are reasonable, at least in shape, at a few percent level in the fully nonlinear regime at all relevant redshifts.

In Figure 25, we show the comparison between the H i power spectra from IllustrisTNG and the N-body simulation in redshift space at several redshifts. As for the comparison in real space, the amplitude of the two power spectra can be quite different for the same reasons outlined above. However, even when the normalization offset is taken into account, there are larger shape differences than in real space. For instance, at z = 1, the discrepancies up to k = 1 h Mpc⁻¹ can be as large as 45%. That difference declines with redshift, so that at z = 5 it is only 6%, very similar to the differences we find in real space.

Zoom In Zoom Out Reset image size

The larger differences in redshift space can be attributed to the Fingers of God. Here, we did not attempt to model the one-halo term and, therefore, the H i Fingers of God are not present in our mock maps, while they are in those created from IllustrisTNG. As we saw in Section 16, the H i Fingers of God can propagate to relatively large scales and affect the amplitude and shape of the H i power spectrum. At higher redshift, the magnitude of the Fingers of God is lower, so the agreement between IllustrisTNG and the N-body mocks is expected to improve, as we observe. Finally, the agreement between the different maps also depends on bandwidth. For larger bandwidths, the Fingers-of-God effects will have a smaller impact. We could also argue that beam smoothing in IM surveys will further reduce the effect of the Fingers of God in the final measured power spectrum.

We emphasize that our results are not converged against resolution (see Appendix A). This should, however, not be seen as a limitation of our claims, as we present results for TNG100, i.e., at the resolution to which the model parameters have been tuned to reproduce a set of galaxy properties. Besides, the purpose of this section is not to claim an absolute calibration of 21 cm maps via N-body simulations, but instead that given the spatial distribution of H i (at a given resolution in our case), most of the information can be condensed into a set of simple parameters and reproduced through N-body simulations.

We conclude that 21 cm intensity maps can be created via less computationally expensive simulations, like N-body, or via fast simulations, e.g., COLA (Tassev et al. 2013) or Pinocchio (Monaco et al. 2002), instead of expensive hydrodynamic simulations. It is, however, very important to account for the one-halo term, i.e., the Fingers of God, as expected when modeling the distribution of H i in redshift space.

A goal of current and upcoming radio telescopes is to map the spatial distribution of matter by detecting 21 cm emission from cosmic neutral hydrogen. The very large volumes that can be sampled through 21 cm intensity mapping observations will place tight constraints on the values of the cosmological parameters (Bull et al. 2015; Villaescusa-Navarro et al. 2015; Obuljen et al. 2017; Sprenger et al. 2018). In order to extract maximum information from these surveys, accurate theory predictions are needed.

Theory predictions to linear order are well known. For instance, the amplitude and shape of the 21 cm power spectrum are given by

$$\begin{eqnarray}&&{P}_{21\mathrm{cm}}(k,\mu )={\bar{T}}_{b}^{2}{({b}_{{\rm{H}}{\rm{I}}}+f{\mu }^{2})}^{2}{P}_{{\rm{m}}}(k)+{P}_{\mathrm{SN}},\end{eqnarray} \tag{ 37 }$$

where ${\bar{T}}_{b}\propto {{\rm{\Omega }}}_{{\rm{H}}{\rm{I}}}$ is the mean brightness temperature, b_{H i} is the H i bias, f is the linear growth rate, μ = k_z/k, P_m(k) is the linear matter power spectrum, and P_SN is the H i shot noise. While the value of ${\bar{T}}_{b}$ is relatively well known from several observations across the redshift range z ∈ [0, 5], little is known about b_{H i} and ${P}_{\mathrm{SN}}$. In this work, we have quantified them and studied the scales where linear theory holds, e.g., when the H i bias becomes nonlinear.

Accurate theory predictions in the mildly or fully nonlinear regimes will allow us to recover the large amount of information embedded in the spatial distribution of H i on small scales. There are several techniques for accomplishing this, such as perturbation theory, H i halo models, or numerical simulations. The purpose of our work has been to study the ingredients that these techniques employ. For example, H i halo models require the halo H i mass function and the H i density profile as inputs.

Our purpose here is not limited to analytic approaches, but also to understand how H i is distributed across the universe and how it evolves with time. We have shown that even with a subset of the ingredients studied in this work, one can model the spatial distribution of H i in the fully nonlinear regime without the use of computationally expensive hydrodynamic simulations, but using H i HOD models.

We have carried out our analysis using IllustrisTNG (Marinacci et al. 2018; Naiman et al. 2018; Pillepich et al. 2018a; Springel et al. 2018; Nelson et al. 2018), a sophisticated series of cosmological hydrodynamical simulations that have been shown to be in broad agreement with many basic observables, and run at an unprecedented combination of large volume and high resolution, therefore providing an excellent testbed for accurately investigating the distribution of H i from the disks of spiral galaxies to cosmological scales.

We outline the main conclusions of our work below:

• 1.  
  We find that almost all H i in the universe are inside halos: from more than 99% at z = 0 to around 88% at z = 5. The fraction of H i outside halos increases with redshift because the gas in the IGM is denser and the amplitude of the UV background decreases with redshift (at z > 2). This justifies the use of halo models to model the distribution of H i in the universe, but quantifies their limitations at high redshifts. The fraction of H i inside galaxies is slightly lower than that in halos. At z = 0, ≃ 97% of all H i are inside galaxies, while at z = 5 this number decreases to ≃80%.
• 2.  
  We find that the halo H i mass function, i.e., the average H i mass hosted by a halo of mass M at redshift z, is well reproduced by a function like

  $$\begin{eqnarray*}&&{M}_{{\rm{H}}{\rm{I}}}(M,z)={M}_{0}{\left(\displaystyle \frac{M}{{M}_{\min }}\right)}^{\alpha }\exp (-{({M}_{\min }/M)}^{0.35}),\end{eqnarray*}$$

  where M₀, M_min, and α are free parameters. The best-fit values are given in Table 1 for both FoF and FoF-SO halos. The value of α increases with redshift, likely indicating that at low redshift, processes such as ram-pressure and tidal stripping, and AGN feedback make galaxies in clusters H i poor. We find that M_min decreases with redshift. On the other hand, only halos with circular velocities above around 30 km s⁻¹ host a significant H i mass fraction. Although the fit is slightly worse, our halo H i mass function can also be well reproduced by the function

  $$\begin{eqnarray*}&&{M}_{{\rm{H}}{\rm{I}}}(M,z)={M}_{0}{\left(\displaystyle \frac{M}{{M}_{\min }}\right)}^{\alpha }\exp (-{M}_{\min }/M),\end{eqnarray*}$$

  the best-fit values of which are given in Table 6 to facilitate comparison with previous works.
• 3.  
  We find that the H i density profiles inside halos exhibit a large halo-to-halo variation. H i profiles are sensitive to the physical processes that occur in and around halos, such as AGN feedback and tidal stripping. The H i profile of small halos (M ≲ 10¹² h⁻¹ M_⊙) is not cuspy, but its amplitude saturates. This is expected as H i at high densities will turn into molecular hydrogen and then stars over short time periods. More massive halos exhibit holes in their centers. For galaxy groups, this is mostly due to AGN feedback, while in galaxy clusters the holes are large and generated by a combination of AGN feedback, and ram-pressure and tidal stripping. We find that the average H i density profiles are universal and can be reproduced by an expression like

  $$\begin{eqnarray}&&{\rho }_{{\rm{H}}{\rm{I}}}(r)=\displaystyle \frac{{\rho }_{0}}{{r}^{{\alpha }_{\star }}}\exp (-{r}_{0}/r)\end{eqnarray} \tag{ 38 }$$

  or

  $$\begin{eqnarray}&&{\rho }_{{\rm{H}}{\rm{I}}}(r)=\displaystyle \frac{{\rho }_{0}{r}_{s}^{3}}{(r+3/4{r}_{s}){(r+{r}_{s})}^{2}}\exp (-{r}_{0}/r),\end{eqnarray} \tag{ 39 }$$

  where ρ₀, α_*, r₀, and r_s are free parameters. We fix the value of ρ₀ by requiring that ${M}_{{\rm{H}}{\rm{I}}}(M)={\int }_{0}^{{R}_{v}}4\pi {r}^{2}{\rho }_{{\rm{H}}{\rm{I}}}(r){dr}$, where R_v is the halo virial radius. The best-fit values for α_*, r₀, and r_s can be found in Table 2.
• 4.  
  We find that the H i mass in small/big halos is mostly located in its central/satellite galaxies. The fraction of the total H i mass in halos that is within the central galaxy decreases with halo mass, while the opposite trend occurs for the satellites. For halos of masses ∼5 × 10¹² h⁻¹ M_⊙, the H i mass in the central galaxy is similar to that of the satellites, almost independent of redshift. The H i mass fraction in the central galaxy of clusters is negligible at z = 0. At high redshift, $z\geqslant 2$, the fraction of the halo H i mass in satellites is roughly 20% for small halos M ∈ [10¹⁰–10¹¹] h⁻¹ M_⊙.
• 5.  
  We find that the pdf of the H i density field is quite different from that of the matter field. In general, the H i pdf is broader, indicating that the H i is more clustered than matter. The amplitude of the pdf for low overdensities is higher for H i than for matter, indicating that H i voids are emptier than matter voids. At high redshift, the H i and matter density pdf can be well reproduced by a log-normal, while at low redshift, the log-normal is not a good description of our results.
• 6.  
  We find that the H i column density distribution function is nearly constant across redshifts, in agreement with previous studies and with observations. In the redshift range z ∈ [2, 4], we find that the DLA cross-section depends on both halo mass and H i column density, and its mean value can be well reproduced by

  $$\begin{eqnarray}&&{\sigma }_{\mathrm{DLAs}}(M| {N}_{{\rm{H}}{\rm{I}}},z)={{AM}}^{\alpha }(1-{e}^{-{(M/{M}_{0})}^{\beta }}),\end{eqnarray} \tag{ 40 }$$

  where α = 0.82, $\beta =0.85{\mathrm{log}}_{10}({N}_{{\rm{H}}{\rm{I}}}/{\mathrm{cm}}^{-2})-16.35$, A · M₀ = 0.0141 h⁻² kpc M_⊙, and the best-fit values of M₀ are given in Table 3. We argued that the small dependence of the above relation on column density implies that the bias of different absorbers will be very similar. We estimate the DLA bias using two methods and find agreement with observations in both.
• 7.  
  We find that for small halos, M ≲ 10¹² h⁻¹ M_⊙, the bulk velocities of H i inside halos trace very well, in modulus and direction, the peculiar velocity of the halos they reside in. On the other hand, for bigger halos, we observe differences, in modulus and direction, between the H i and halo peculiar velocities. This happens because while for small halos most of the H i are in the central galaxy, for larger halos a significant H i mass is in satellites, whose peculiar velocities do not trace that of the halo.
• 8.  
  We find that the velocity dispersion of H i inside halos can be well reproduced by a simple power law,

  $$\begin{eqnarray}&&\sigma (M)={\sigma }_{10}{\left(\displaystyle \frac{M}{{10}^{10}{h}^{-1}{M}_{\odot }}\right)}^{\alpha },\end{eqnarray} \tag{ 41 }$$

  where σ₁₀ and α are free parameters whose best-fit values are given in Table 4. While at z = 0 the mean velocity dispersions of CDM and H i are similar (for halos above ≃5 × 10¹⁰ h⁻¹ M_⊙), at higher redshifts they diverge for small halos, with H i having a lower amplitude than CDM. The mass where they diverge increases with redshift, but is typically around ${10}^{12}\,{h}^{-1}\,{M}_{\odot }$. In general, for fixed mass and redshift, the variance in the velocity dispersion of H i is larger than that of CDM, reflecting the larger variation in H i profiles than CDM profiles inside halos.
• 9.  
  We find that the values of the H i bias on the largest scales we can probe is equal to 0.84, 1.49, 2.03, 2.56, 2.82, and 3.18 at redshifts 0, 1, 2, 3, 4, and 5, respectively. While the H i bias is relatively flat down to k ≃ 1 h Mpc⁻¹ at z = 1, it is already nonlinear at k ≃ 0.3 h Mpc⁻¹ at z ≥ 3. Our results suggest that the H i bias becomes more nonlinear with redshift. This is expected as the value of the linear bias increases with redshift. We have shown that the perturbative approaches based on LPT are able to reproduce the clustering measurement up to k = 1 h/Mpc⁻¹, therefore making possible, at least in principle, to extract cosmological information from such small scales.
• 10.  
  We identify a new secondary halo bias. Halos of the same mass are clustered differently depending on their H i mass. At low redshift, H i-poor halos are more clustered than H i-rich halos. However, at high redshift, the situation is reversed, and H i rich halos cluster more strongly than H i-poor halos. We believe that this is mainly driven by environment. At low redshift, small halos may lose their gas due to stripping by a larger neighboring halo, so H i-poor halos will be more clustered than field H i-rich halos. On the other hand, at high redshift, gas stripping is likely less effective, so H i-rich halos will be found around larger halos, and therefore, their clustering will be higher.
• 11.  
  We quantify the amplitude of the H i shot noise to be 104, 124, 65, 39, 14, and 7 (h⁻¹ Mpc)³ at redshifts 0, 1, 2, 3, 4, and 5, respectively. These low levels imply that BAO measurements through 21 cm intensity mapping are hardly affected by shot noise. Furthermore, the very low shot-noise levels at high redshift suggest that a large amount of cosmological information can be extracted from the clustering of H i on small scales.
• 12.  
  We find that the relation between ρ_m and ${\rho }_{{\rm{H}}{\rm{I}}}$ cannot be explained with linear theory for smoothing scales ≤5 h⁻¹ Mpc at any redshift. The scatter in that relation decreases with redshift, and much larger H i overdensities can be found for the same matter overdensities.
• 13.  
  We find that the Kaiser factor alone cannot explain the clustering of H i in redshift space at low redshift, as expected, given the small volume of our simulations. But, at high redshift, the ratio between the monopoles in redshift and real space can be explained with linear theory down to 0.3, 0.5, and 1 h Mpc⁻¹ at redshifts 3, 4, and 5, respectively. This is rather surprising, taking into account that the H i bias becomes nonlinear already at k = 0.3 h Mpc⁻¹ at those redshifts.
• 14.  
  We find that the two-dimensional H i power spectrum in redshift space exhibits large differences with respect to the matter field. Those differences arise mainly because the amplitude of the Fingers of God is higher for H i than for matter. This can be understood by taking into account that H i resides only in relatively massive halos. While the amplitude of the matter power spectrum on small scales is dominated by small halos with low velocity dispersion, only H i halos above ≃M_min, i.e., with larger velocity dispersion, can contribute. We find that standard phenomenological models to describe the clustering in 2D in redshift space are not adequate for reproducing our results.
• 15.  
  We show that accurate 21 cm maps can be created from N-body simulations, rather than full hydrodynamic simulations, by using the ingredients studied in our work. In real space and without modeling the one-halo term, the agreement in the shape of the 21 cm power spectrum from N-body and IllustrisTNG is around 5% down to 1 h Mpc⁻¹ at all redshifts. In redshift space, however, the lack of the one-halo term, i.e., the H i Fingers of God, induces much larger errors in the 21 cm power spectrum from N-body versus hydrodynamic simulations at low redshift, e.g., 45% at z = 1. Modeling the one-halo term is thus crucial for creating mock 21 cm maps.

We emphasize that our results are not converged against resolution (see Appendix A). This should, however, not be seen as a limitation of our claims, as we present results for TNG100, i.e., at the resolution to which the model parameters have been tuned to reproduce a set of galaxy properties.

The H i properties investigated in this work will help to improve our knowledge of the way neutral hydrogen is distributed across the universe. The different quantities we have studied can be used as input to analytic approaches like H i halo models or to create very accurate mock 21 cm maps.

The Python/Cython scripts written to carry out the analysis performed in our work can be found in https://github.com/franciscovillaescusa/Pylians/tree/master/HI_Illustris. Our scripts made use of the Pylians Python routines, publicly available at https://github.com/franciscovillaescusa/Pylians.

We thank the anonymous referee for his/her report, which helped us to improve the quality of this work. We thank Uros Seljak for detailed comments on shot noise in discrete tracers. We also thank Sandrine Codis, Philippe Berger, Andreu Font-Ribera, Ariyeh Maller, Anze Slosar, Amiel Sternberg, Martin White, and Cora Uhlemann for useful discussions. The work of F.V.N., S.G., and D.N.S. is supported by the Simons Foundation. The analysis of the simulations has been carried out on the Odyssey cluster at Harvard University. The IllustrisTNG simulations were run on the HazelHen Cray XC40 supercomputer at the High-Performance Computing Center Stuttgart (HLRS) as part of project GCS-ILLU of the Gauss Centre for Supercomputing (GCS).

In this section, we study whether our results are converged or are at least converging against changes in particle mass and spatial resolution. In Figure 2, we found that the values of Ω_{H i}(z) for TNG100 and TNG300 are in agreement to better than 10% below z = 1 and to better than 25% below z ∼ 3, but can disagree more substantially at higher redshifts. We interpreted that result considering the fact that TNG300 cannot resolve halos that host H i as small as those that TNG100 resolves. Besides, we also stated that for well-resolved halos, the total H i mass inside them is resolution-dependent.

We now explicitly show these results by comparing the M_{H i}(M, z) function from TNG100-1, TNG100-2, and TNG100-3 in the left panel of Figure 26. TNG100-1 corresponds to what we have called TNG100 in this work, while TNG100-2 and TNG100-3 only differs from TNG100-1 in mass resolution: TNG100-1/-2/-3 follow initially 1820³/910³/455³ CDM particles and 1820³/910³/455³ gas cells.

Zoom In Zoom Out Reset image size

As can be seen, for halo masses above ∼10¹² h⁻¹ M_⊙, our results (dis)agree on average by up to a factor of 1.6–1.7 (1.3–1.4) within a factor of 64 (8) in particle mass resolution. At the low-mass end, the results of the different simulations can (dis)agree by up to factors of 3 or more. This is not surprising: while, e.g., halos of 10⁹ h⁻¹ M_⊙ are relatively better resolved for TNG100-1, they only host a few hundred resolution elements in TNG100-2, and they are below the mass resolution for TNG100-3. On the other hand, at the high-mass end, we find that the average H i mass still depends on resolution, with larger H i masses at progressively better resolutions. This is interesting and possibly related to the H i/H₂ postprocessing, as in fact the total gas fraction within the virial radius of halos is consistent to better than 10%–20% across the TNG100-1, TNG100-2, and TNG100-3 resolution levels (see Figure A2 of Pillepich et al. 2018a).

A lack of convergence can be seen in the H i power spectrum, the results for which we show in the middle panel of Figure 26. The three different resolution runs are consistent with one another to better than 30% across spatial scales, but within this agreement (or discrepancy, according to the point of view and needs of the reader), the amplitude of the power spectrum at different spatial separations does not necessarily change monotonically at increasingly better resolutions.

The right panel of Figure 26 shows the results for the H i density profile within halos of masses M ∈ [1–2] × 10¹² h⁻¹ M_⊙. Within the studied resolution range, the TNG results (dis)agree within a factor of 1.4–1.6 outside radii of  10 kpc, but can disagree more substantially in the innermost regions of the halos. In fact, at radii near the virial radius, results are almost converged to better than 20%, suggesting that differences among resolutions are not due to small differences in the halo masses or radii.

All previous results are at z = 0. We find similar results at higher redshifts for the considered quantities above. We have repeated the above exercise with a simulation with 25 h⁻¹ Mpc and 2 × 1024³ CDM particles and gas cells, i.e., at a higher mass resolution than TNG100-1, reaching similar conclusions.

We emphasize that our H i mass results, although they are not fully converged at the TNG100 or TNG300 resolution or are not converging against resolution, are still valuable, as they come from a physical model that is consistent (to different levels of agreement) with a large number of observed galaxy properties and statistics.³⁰

In Section 4, we found that while at z ≤ 2 most of the universe H i mass is inside halos, at higher redshift, an increasing fraction of it is located outside halos. In order to visualize this effect, we show in Figure 27 the spatial distribution of H i in a slice of 5 h⁻¹ Mpc width across 10 × 10 (h⁻¹ Mpc)². We show in that figure the H i column density in comoving units to facilitate the comparison across redshifts. We note that our color palette may produce the incorrect impression that at z = 5 there is much more H i than at z = 0. We have explicitly checked that the sum of all column densities across all pixels in our figures gives a similar value across redshifts, indicating that Ω_{H i} is very similar in all panels.

Zoom In Zoom Out Reset image size

It can be seen that at low redshift, most of the H i mass is inside galaxies, while the hydrogen in the filaments is highly ionized. At higher redshifts, on the other hand, the gas in the intergalactic medium becomes denser and the filaments contain a larger amount of H i. The lower amplitude of the UV background at those redshifts facilitates gas self-shielding. Given these effects, it is thus natural that the fraction of H i outside dark matter halos increases with redshift.

We saw in Section 5 that the halo H i mass function exhibits a cutoff at low masses. In this appendix, we shed light on the physical origin of that feature.

The lack of H i gas in small halos may be due to (1) a deficit in the abundance of gas in those halos, (2) gas being present but highly ionized, or a mixture of (1) and (2). We quantify the gas content of dark matter halos in TNG100 by computing their gas fraction, i.e., the ratio between the gas mass to the total mass. In the left panel of Figure 28, we show the average gas fraction as a function of halo mass at different redshifts. In this plot, we show results only for halos whose masses are above the mass of 50 CDM particles.

Zoom In Zoom Out Reset image size

We find that the gas fraction of the smallest halos shown at high redshifts is around ∼0.12. As we go to higher halo masses, some stars form, and supernova feedback expels the gas of these halos. In even more massive halos, the gravitational potential is deeper, and hence the mass-loading factors of galactic winds in the TNG model are lower, rendering it more difficult for supernova feedback to eject gas from halos. This explains the dip we observe around 5 × 10⁹ h⁻¹ M_⊙ halos at high redshift. In halos with masses above ∼3 × 10¹¹ h⁻¹ M_⊙, AGN feedback becomes effective at expelling their gas, explaining the peak at around that mass. Finally, as we go to even higher mass halos, a smaller fraction of the gas can be expelled by AGN feedback since the gravitational potential becomes deeper, explaining the dip around 3 × 10¹² h⁻¹ M_⊙ halos.

It is interesting to note that the gas fraction of small halos decreases quickly with redshift. We will see below that this is caused by the UV background. At low redshift, the gas fraction of small halos is small, so it is reasonable to expect that very little H i is found in those halos. However, as we go to higher redshifts, the gas fraction of halos at the hard cutoff mass ${M}_{\mathrm{hard}}$ (defined such that halos with masses below M_hard host only 2% of all H i that are in halos; see Equation (14)) is rather large. We show this with colored stars in that figure. Thus, the cutoff in the halo H i mass function cannot be explained, at high redshift, by the lack of gas in halos. A better explanation is that the gas in these halos is highly ionized.

The average H i mass to gas mass ratio is shown in the right panel of Figure 28. We find that halos at the hard cutoff mass M_hard exhibit similarly low H i to gas mass fractions across redshifts: ≃5%. Thus, while there is a significant amount of gas in these halos at high redshift, H i formation is impeded, likely due to the gas being at high temperature and diffuse. The low H i to gas mass ratio shows up at all redshifts for low-mass halos. However, this is the case particularly at low redshift, where the H i to gas mass fraction of small halos is practically zero, showing how difficult it is to form H i in these halos. We speculate that this is related to the fact that for fixed halo mass, the density of gas and CDM increases with redshift. Thus, for example, while for halos with masses ≃10⁹ h⁻¹ M_⊙ the gas fraction increases by only ≃60% from z = 3 to z = 5, the H i to gas mass ratio changes by more than 400%, showing how denser gas enables the formation of H i.

We thus conclude that the reason why there is almost no H i in small halos is because the gas in these halos is highly ionized, presumably because its low density and high temperature prevent H i formation.

It is interesting to understand why the gas fraction of small halos decreases with redshift. We argue that this effect is due to the UV background. In order to demonstrate our claim, we have run three hydrodynamic simulations with radiative cooling and star formation but without feedback (neither galactic winds nor AGN). In two of them, no heating by the UV background is included, one with 2 × 256³ CDM+baryon resolution elements and another with 2 × 512³ CDM+baryon resolution elements. In the third one, we have heating by the UV background with 2 × 512³ CDM+baryon resolution elements. In all cases, the simulation box is 25 h⁻¹ Mpc across.

We have computed the average baryon fraction, i.e., the mass in gas and stars over the total mass, in these simulations and show the results in Figure 29. The top panel shows the baryon fraction for these simulations and TNG100 at z = 0. We find that in simulations without feedback, the gas content of small halos exhibits a cutoff that occurs at higher masses when the UVB is present.

Zoom In Zoom Out Reset image size

The bottom panels of Figure 29 show the time evolution of the average baryon fraction as a function of halo mass. It can be seen that in both types of simulations, with and without UV background, the gas fraction of small halos decreases with redshift. We speculate that in the case of no UV background, many small halos lie in the vicinity of massive halos, whose presence can strip the gas from these small halos, and further, that the same mechanism may explain the dependence of halo clustering on H i mass we described in Section 14. Although not shown here, we find that our results at z > 5 for the case with no UV background are in agreement with those by Qin et al. (2017). For the case with UV background and no feedback, our results are also in agreement with Qin et al. (2017) for redshifts below the time where the UV background is turned on.

The effect of the UV background on the baryon fraction of small halos is more pronounced, as can be seen in the bottom-right panel of Figure 29. We believe that the reason for this behavior is that the hot intergalactic gas cannot cluster in small halos since their gravitational potential is not deep enough (Okamoto et al. 2008; Bose et al. 2018).

We thus conclude that at low redshift, small halos have a low gas fraction. There are several mechanisms that can remove the gas from such halos, such as tidal stripping by neighbors, the heating of the intergalactic medium by the UV background, and supernova feedback. Our results suggest the most effective one for the lowest masses is the presence of the UV background. The little gas inside those halos is highly ionized, so no H i is found within them.

We found in Section 4 that the H i mass inside FoF and FoF-SO halos is quite different. Here, we determine the reason for this difference.

We have computed the total mass inside FoF and FoF-SO halos, and we find that the former host ∼9% more mass than the latter at z = 0. This difference is almost equal to the deficit in H i mass we find between FoF-SO and FoF halos at that redshift (see Figure 3). Similar results hold at higher redshift, where the differences in total mass and H i are slightly larger (≃25% at z = 5). This indicates that the deficit in H i mass we find in FoF-SO halos with respect to FoF is simply due to the fact that the latter host a larger total mass, and therefore more H i.

We note, however, that there are some situations where the difference in H i mass can be much larger than the difference in total mass. We illustrate one of these situations in Figure 30. We have selected the gas cells belonging to an FoF halo of mass M = 5 × 10¹³ h⁻¹ M_⊙ at redshift z = 0. The total mass inside the FoF and FoF-SO halos are nearly the same, while the H i masses vary by a factor of 2. In Figure 30, we plot the column density of gas and H i from the gas cells belonging to that halo. In the same plot, we mark the radius of the corresponding FoF-SO halo with a white line. It can be seen that while most of the gas in the FoF halo is inside the virial radius of the SO halo, the situation is quite different for H i. The H i mass outside the SO radius is almost equal to the one inside. The reason is that the FoF algorithm links the external galaxies to the main halo, and these galaxies are rich in H i. We have found that this situation is usual for the most massive halos at each redshift.

Zoom In Zoom Out Reset image size

In Section 5, we found that the halo H i mass function, M_{H i}(M, z), from IllustrisTNG can be well reproduced by a fitting formula like ${M}_{0}{x}^{\alpha }\exp (-1/{x}^{0.35})$, where $x=M/{M}_{\min }$. In the literature, however, an expression of the type ${M}_{0}{x}^{\alpha }\exp (-1/x)$ has been widely used (e.g., Bagla et al. 2010; Castorina & Villaescusa-Navarro 2017; Obuljen et al. 2017; Padmanabhan et al. 2017; Villaescusa-Navarro et al. 2017; Pénin et al. 2018). We have fit our results to that function, and while the reduced χ² is larger than that with our fiducial function, it is still a good fit to the underlying data. We thus provide in Table 6 the best-fit values of the parameters α, M₀, and M_min of the latter expression in order to help in the comparison with previous works. We note that the value of M_hard, i.e., the hard cutoff mass (see Section 5), is not affected by using a different parameterization for the halo H i mass function. Thus, the values we quote in Table 1 are valid also here.

The points with error bars in Figure 31 show the mean H i density profiles within halos of different halo masses at different redshifts (they are the same as the blue lines in Figure 6). The solid lines represent the best fit obtained when we fit those results with a H i density profile as

$$\begin{eqnarray}&&{\rho }_{{\rm{H}}{\rm{I}}}(r)=\displaystyle \frac{{\rho }_{0}{r}_{s}^{3}}{(r+3/4{r}_{s}){(r+{r}_{s})}^{2}}\exp (-{r}_{0}/r).\end{eqnarray} \tag{ 42 }$$

In each panel of the plot, we show the best-fit values of r_s and r₀, and the value of the reduced χ². The value of ρ₀ is fixed by requiring that ${M}_{{\rm{H}}{\rm{I}}}=4\pi {\int }_{0}^{{R}_{v}}{r}^{2}{\rho }_{{\rm{H}}{\rm{I}}}(r){dr}$, where R_v is the halo virial radius.

Zoom In Zoom Out Reset image size