A scale dependent bias on linear scales: the case for HI intensity mapping at z=1

TL;DR

This paper demonstrates that higher order non-linear corrections significantly affect the HI power spectrum at z=1, impacting bias, BAO features, and cosmological inferences in future intensity mapping surveys.

Contribution

It provides a full one-loop derivation of the HI power spectrum, revealing scale-dependent biases and non-linear effects crucial for precision cosmology.

Findings

HI bias drops by up to 15% on linear scales at z=1

Scale dependence in the power spectrum affects BAO peak ratios

Full non-linear corrections cause a fivefold difference in power spectrum amplitude

Abstract

Neutral hydrogen (HI) will soon be the dark matter tracer observed over the largest volumes of Universe thanks to the 21 cm intensity mapping technique. To unveil cosmological information it is indispensable to understand the HI distribution with respect to dark matter. Using a full one-loop derivation of the power spectrum of HI, we show that higher order corrections change the amplitude and shape of the power spectrum on typical cosmological (linear) scales. These effects go beyond the expected dark matter non-linear corrections and include non-linearities in the way the HI signal traces dark matter. We show that, on linear scales at z = 1, the HI bias drops by up to 15% in both real and redshift space, which results in underpredicting the mass of the halos in which HI lies. Non-linear corrections give rise to a significant scale dependence when redshift space distortions arise, in particular on the scale range of the baryonic acoustic oscillations (BAO). There is a factor of 5 difference between the linear and full HI power spectra over the full BAO scale range, which will modify the ratios between the peaks. This effect will also be seen in other types of survey and it will be essential to take it into account in future experiments in order to match the expectations of precision cosmology.

Figures

Tables

Table 1. Table 1: Free parameters of the MHIMh prescriptions along with the associated HI biases and mean temperatures at z = 1 𝑧 1 z=1 .

Model	Parameters	$b_1$	$b_2$	$b_3$	$b_{\mathrm{eff}}$	$T_{\mathrm{HI}}\times {10}^4$ K
Bagla10	None	0.93	-0.41	0.62	0.80	1.65
AGN	$\alpha =0.73$, $\gamma =2$	0.91	-0.27	0.41	0.82	12.14
21cm	$\alpha$ = 0.15	0.96	-0.42	0.60	0.81	2.43
HOD A	$\log v_{c,0}$ = 1.58 , $\log v_{c,1}$ = 3.14, $\alpha =0.17$ , $\beta =-0.5$	1.00	-0.35	0.38	0.82	4.38
HOD B	${\log }_{10}{v}_{c,0}$ = 1.56 , $\alpha =0.09$ , $\beta =-0.58$	0.96	-0.37	0.49	0.85	2.55
DLA50	$\alpha =0.13$	1.64	0.56	-1.27	1.74	49.94

Equations

$$\overline{T}(z)=566h\left(\frac{H_0}{H(z)}\right)\left(\frac{{\Omega }_{\mathrm{HI}}(z)}{0.003}\right){\left(1+z\right)}^2\mu \mathrm{K}$$

$$T(z,\mathbf{x})=\overline{T}(z)\left(1+{\delta }_{\mathrm{HI}}(\mathbf{x})\right)$$

$$⟨T(z,\mathbf{k})T^{\star }(z,{\mathbf{k}}^{\prime })⟩=(2\pi {)}^3{P}_{\mathrm{HI}}(k,z){\delta }^3(\mathbf{k}-{\mathbf{k}}^{\prime })$$

$$P_{\mathrm{HI}}(z,k)=P_{\mathrm{HI}}^{11}(z,k)+P_{\mathrm{HI}}^{22}(z,k)+P_{\mathrm{HI}}^{13}(z,k)$$

$$P_{\mathrm{HI}}^{11}(k)$$

$$P_{\mathrm{HI}}^{22}(k)$$

$$P_{\mathrm{HI}}^{13}(k)$$

$${\sigma }_{\Lambda }^2={\int }_{k_{\min }}^{k_{\max }}\frac{{\mathrm{d}}^{3}k}{(2\pi {)}^3}{P}_{\mathrm{m}}(k)$$

$$P_{\mathrm{HI}}(k,\mu )=P_{\mathrm{HI}}^{11}(k,\mu )+P_{\mathrm{HI}}^{22}(k,\mu )+P_{\mathrm{HI}}^{13}(k,\mu )$$

$$P_{\mathrm{HI}}^{11}(k,\mu )={\overline{T}}^2{\left[b_1+f{\mu }^2\right]}^2{P}_{\mathrm{m}}^{11}(k)$$

$$P_{\mathrm{HI}}^{22}(k,\mu )$$

$$P_{\mathrm{HI}}^{13}(k,\mu )$$

$${\rho }_{\mathrm{HI}}=\int \mathrm{d}M\frac{\mathrm{d}n}{\mathrm{d}M}{M}_{\mathrm{HI}}(M)$$

$$b_n^{\mathrm{HI}}=\frac{1}{{\rho }_{\mathrm{HI}}}\int \mathrm{d}M\frac{\mathrm{d}n}{\mathrm{d}M}{b}_n^{\mathrm{h}}(M)M_{\mathrm{HI}}(M)$$

$$M_{\mathrm{HI}}(M)=\frac{f_{3}M}{1+\frac{M}{M_{\max }}}\text{mbox}forM\ge M_{\min }$$

$$v_{\mathrm{circ}}=30\sqrt{1+z}{\left(\frac{M}{10^{10}{M}_{\odot }}\right)}^{1/3}\text{mbox}km/s$$

$$M_{\mathrm{HI}}(M)=\alpha f_{\mathrm{H},\mathrm{c}}M\exp \left[-{\left(\frac{v_{c,0}}{v_c(M)}\right)}^3\right]\exp \left[-{\left(\frac{v_{c,1}}{v_c(M)}\right)}^3\right]$$

$$v_c(M)=96.6\text{mbox}km/s{\left(\frac{{\Delta }_v{\Omega }_m{h}^2}{24.4}\right)}^{1/6}{\left(\frac{1+z}{3.3}\right)}^{1/2}{\left(\frac{M}{10^{11}{M}_{\odot }}\right)}^{1/3}$$

$$M_{\mathrm{HI}}(M)$$

$$M_{\mathrm{HI}}(M)=\alpha f_{\mathrm{H},\mathrm{c}}M{\left(\frac{M}{10^{11}{h}^{-1}{M}_{\odot }}\right)}^{\beta }\exp \left[-{\left(\frac{v_{c,0}}{v_c(M)}\right)}^3\right]$$

$$b_{\mathrm{HI}}^{\mathrm{eff}}(k)=\frac{1}{{\overline{T}}_{\mathrm{HI}}}\sqrt{\frac{P_{\mathrm{HI}}(k)}{P_{\mathrm{m}}^{11}(k)}}$$

$$b_1^{\mathrm{HI}}\to b_{\mathrm{HI}}^{\mathrm{eff}}\approx b_1^{\mathrm{HI}}+\frac{1}{2}\left(b_3^{\mathrm{HI}}+\frac{68}{21}{b}_2^{\mathrm{HI}}\right){\sigma }_{\Lambda }^2$$

$$b_1^{\mathrm{HI}}\to b_{\mathrm{HI}}^{\mathrm{eff}}(k,\mu )\approx b_1^{\mathrm{HI}}+\frac{1}{2}\left[\left(b_3^{\mathrm{HI}}+\frac{68}{21}{b}_2^{\mathrm{HI}}\right){\sigma }_{\Lambda }^2+I_R(k,\mu )\right]$$

$$P_{\mathrm{m}}^{13}(k)$$

$$P_{\theta }^{13}(k)$$

$$I_R(k,\mu )$$

$$B_1(r)$$

$$B_2(r)$$

$$B_3(r)$$

$$B_4(r)$$

$$B_5(r)$$

$$P_{\mathrm{m}}^{\mathrm{NL}}(k)=P_{\mathrm{m}}^{11}(k)+P_{\mathrm{m}}^{22}(k)+P_{\mathrm{m}}^{13}(k)$$

$$P_{\mathrm{m}}^{22}(k)=\frac{1}{2}\int \frac{{\mathrm{d}}^3{k}_1}{(2\pi {)}^3}{F}_2^2({\mathbf{k}}_1,{\mathbf{k}}_2)P_{\mathrm{m}}^{11}(k_2)P_{\mathrm{m}}^{11}(k_1)$$

$$P_{\mathrm{m}}^{\mathrm{NL}}(k,\mu )$$

Full text

††pagerange: A scale dependent bias on linear scales: the case for HI intensity mapping at $z=1$–References††pubyear: 2016

1 Introduction

Future neutral hydrogen (HI) experiments such as the Square Kilometer Array (SKA, Santos et al., 2015), its pathfinder MeerKAT, the Canadian Hydrogen Intensity Mapping Experiment (CHIME, Newburgh et al., 2014), the Hydrogen Intensity and Real-time Analysis eXperiment (HIRAX, Newburgh et al., 2016), and the Baryon acoustic oscillations In Neutral Gas Observations (BINGO, Battye et al., 2013) will map the cosmological neutral hydrogen within unprecedented volumes of the Universe thanks to the line intensity mapping (IM) technique. This technique relies on the measurement of the HI integrated intensity from hundreds of galaxies in one single large voxel (3D pixel) instead of detecting individual HI galaxies. The observed volumes will allow unrivalled constraints on cosmology (e.g. Bull et al., 2015). Nevertheless, to achieve the expected levels of accuracy one needs to understand how HI relates to the underlying dark matter distribution. Current HI observations are quite sparse. At $z\sim 0$, HI is observed in emission at 21 cm but current detections weaken quickly and vanish at $z>0.1$. At intermediate and higher redshifts, the main tracer of HI is Damped Ly-$\alpha$ systems (DLAs), objects with $N_{\mathrm{HI}}>10^{20.3}$ cm*-2* displaying a 21 cm line in absorption in the spectrum of a distant quasar. As they are optically thick, hydrogen into their midst is self-shielded and remains neutral. DLAs are actually thought to host most of the neutral gas within $0<z<5$ (Wolfe et al., 1986; Lanzetta et al., 1995; Wolfe et al., 2005) and hence to contain a significant reservoir of neutral gas for star formation at high redshift. Combining emission and absorption measurements, the redshift evolution of the fraction density of HI $\Omega_{\mathrm{HI}}$ has been shown to decrease slightly from high to low redshift (Crighton et al., 2015; Sánchez-Ramírez et al., 2016). Such a mild but somewhat steady evolution leads to the picture of a balance between consumption and replenishment of the gas reservoir.

Notwithstanding, even though measurements are used altogether, 21 cm emission and 21 cm absorption line surveys might not target the same population of objects. This is crucial in order to clarify what is the HI bias. DLAs at high redshift might not belong to the same population than HI galaxies at low redshift. Properties of DLA hosts remain largely unknown, either because the background quasar is several magnitudes brighter or because they are too faint to be detected by current spectrographs ($m>25$, Rauch et al., 2008; Cen, 2012). When it comes to the mass of their host dark matter halos, there seems to be a tension between 21 cm low redshift galaxies and DLAs. There are only a handful of measurements of HI and DLA biases. Martin et al. (2012) measured $b_{\mathrm{HI}}\sim 0.8$ at $z\sim 0$ in the ALFALFA survey while Chang et al. (2010); Masui et al. (2013), and Switzer et al. (2013) measured the product \Omega_{\mathrm{HI}}$$b_{\mathrm{HI}} in IM data taken with the Green Bank Telescope at $z\sim 0.8$. The latter used the IM data in auto-correlation while the former cross-correlated them with galaxy surveys to circumvent the contamination of foregrounds residuals. Font-Ribera et al. (2012) measured $b_{\mathrm{DLA}}=2.17\pm 0.2$ at $z\sim 2.3$ in the Baryon Oscillation Spectroscopic Survey. Such a value leads to host dark matter halos of $10^{11.5}$M⊙ as compared to the $10^{9-11}$M⊙ found with 21 cm measurements as well as in simulations (Pontzen et al., 2008; Rahmati & Schaye, 2014). To reconcile bias measurements, Padmanabhan et al. (2016) argued that there must be a significant change in the properties of HI-bearing systems.

The knowledge of the HI bias requires to understand how HI populates dark matter halos. Even though it is widely accepted that HI is within galaxies at $z<5$, today a simple relation between dark matter halo and HI masses (MHIMh) is used, with a certain gas profile when necessary (which is not for the bias). The MHIMh relation is measured in hydrodynamical simulations often assuming a simple power law (Davé et al., 2013; Kim et al., 2016; Villaescusa-Navarro et al., 2016), inspired from observations (Bagla et al., 2010) or parametrised and fitted on data (Barnes & Haehnelt, 2010, 2014; Padmanabhan et al., 2016; Padmanabhan & Refregier, 2017). Lately, Padmanabhan & Kulkarni (2017) derived the MHIMh relation using the abundance matching technique where the halo mass function is matched to the HI mass function. Even if these different schemes can strongly differ, they all lead to similar values of the linear HI bias.

On non-linear scales, the HI bias has been barely investigated yet, while it contains a wealth of information on cosmology and, above all, on the MHIMh relation. To date, the scale dependence of the HI bias has been measured in two ways : in hydrodynamical simulations (Villaescusa-Navarro et al., 2014) and in N-body simulations where halos are populated with HI through an empirical relation (Sarkar et al., 2016; Seehars et al., 2016). However these methods suffer from a few limitations. In hydrodynamical simulations, the bias is sensitive to the physical processes that are included in the simulation, and, first and foremost, to the resolution. For instance, several zoom-in simulations will not lead to the same value of the bias at the same common scale. In addition, hydrodynamic simulations can hardly access linear scales and both approaches are computationally heavy. The investigation of the influence of the MHIMh scheme on the HI bias on the full scale range requires a more flexible approach.

We use the full one-loop calculation of Umeh et al. (2016) and Umeh (2016) to compute the non-linear power spectrum of HI. It relies on high order HI biases that are computed with the halo model for each MHIMh prescription. We compute the power spectrum and the bias of HI in both real and redshift space and show that non-linear terms have a significant contribution on linear scales. We limit our analysis to $z=1$ which is one of the most targeted redshifts for BAO measurements.

This paper is organised as follows. We begin with reviewing the theoretical framework of the HI power spectrum and listing the MHIMh relations we use in Sect. 2. Second, we compute the power spectrum and bias in real space with which we examine the mass of halos in which HI lie in Sect. 3. Third, we carry a similar analysis in redshift space and discuss its cosmological implications in Sect. 4. We conclude in Sect. 5

Throughout the article, we use the Planck 2014 Cosmology (Planck Collaboration et al., 2014).

2 Modelling the HI power spectrum

2.1 The power spectrum

The average HI brightness temperature is given by (Battye et al., 2013)

[TABLE]

where the HI density fraction is defined as $\Omega_{\mathrm{HI}}=\rho_{\mathrm{HI}}/\rho_{c,0}$ with $\rho_{c,0}$ is the critical density of the Universe today. The fluctuating part is

[TABLE]

with $\delta_{\mathrm{HI}}(\mathbf{x})$ the HI density fluctuation at position $\mathbf{x}$, hence, in Fourier space

[TABLE]

Carrying a full one-loop derivation of the HI brightness temperature in Perturbation Theory (Bernardeau et al., 2002), the power spectrum of HI in real space at redshift $z$ is

[TABLE]

where $P_{\mathrm{HI}}^{11}(z,k)$ is the linear power spectrum (tree level) while $P_{\mathrm{HI}}^{22}(z,k)$ and $P_{\mathrm{HI}}^{13}(z,k)$ are the non-linear corrections. For clarity purposes we will not specify the redshift dependence in the following. Following Umeh et al. (2016) and Umeh (2016) the three terms of $P_{\mathrm{HI}}(k)$ are

[TABLE]

where $k_{2}=|\mathbf{k}_{1}-\mathbf{k}|$. $b_{1},b_{2}$, and $b_{3}$ are the linear, second and third order HI biases, respectively. The latter are the higher terms of the bias expanded in Taylor series, which means assuming that the HI bias is local. $F_{2}$ is the non-linear density kernel defined in Appendix A. Finally $\sigma_{\Lambda}$, the variance of the dark matter field, is

[TABLE]

For simplicity, we set ${k_{\mathrm{max}}}$ to the non-linear dispersion scale, $k_{\mathrm{NL}}=0.2\,h(1+z)^{2/(2+n_{s})}$ Mpc*-1* with $n_{s}$ the spectral index.

In redshift space, the 3D power spectrum of HI on linear and quasi-linear scales at scale $k$, and $\mu$, the cosine of the angle between the line of sight and the separation vector $\mathbf{k}$, writes

[TABLE]

Following Kaiser (1987), the linear term in redshift space is

[TABLE]

with $\mu=k_{\parallel}/k$, $P_{\mathrm{m}}(k)$ the linear power spectrum of matter, and $f$ the linear growth rate. We compute the former using the transfer function of Eisenstein & Hu (1998) and assume $f(z)=\Omega_{\mathrm{m}}(z)^{\gamma}$ with $\gamma=0.55$ for $\Lambda$CDM (Peebles, 1980; Linder, 2005).

Following Umeh et al. (2016) and Umeh (2016) the one-loop corrections are

[TABLE]

where $P^{13}_{\mathrm{m}}(k)$ and $P^{13}_{\theta}(k)$ are the third order matter power spectrum and velocity field power spectrum, respectively. Their expressions along with that of $I_{R}(k,\mu)$ are listed in Appendix A. Finally, several kernels are involved in the computation of the $P_{\mathrm{HI}}^{22}(k,\mu)$ term : $G_{2}$ induced by peculiar velocities at second order and $K_{R}$ arises from non-linear mode coupling (velocity-velocity and velocity-density, Bernardeau et al., 2002). Their expressions are also given in Appendix A.

2.2 HI quantities

HI-related quantities such as densities and biases are computed using the halo model that provides a description of the clustering of dark matter halos at both linear and non-linear scales (Cooray & Sheth, 2002). It relies on the halo mass function $\mathrm{d}n/\mathrm{d}M$ and the associated $n-$th order halo biases $b_{n}^{h}(M)$ measured in N-body simulations. We use the prescriptions of Sheth & Tormen (1999).

The comoving density of HI writes

[TABLE]

The $n$-th order HI biases are

[TABLE]

where $M_{\mathrm{HI}}(M)$ is the relation between the HI mass and the halo mass (see Sect. 2.3). Fig. 1 shows an example of a set of biases. Note that only the first order bias is always positive while the two others change sign. All of them increase for high halo masses. We will use the terms linear and first order bias interchangeably.

2.3 The HI mass - halo mass relation

The distribution of HI within the Large Scale Structure is rather unclear today. It is believed that in the post-reionization era most of HI lies within galaxies while only a negligible fraction is diffuse (Seehars et al., 2016). It is often simply parametrised by relating the mass of HI to the mass of its host dark matter halo through a simple power law including, or not, a cut-off at small and high halo masses. We compile here several MHIMh relations that have been used or estimated using both hydrodynamical simulations and parametrised models fitted on data measurements. We also consider a DLA model.

Bagla10 : One relation that has been widely used is that of Bagla et al. (2010). It has been inspired from quasar observations and assumes that there is no HI in high mass halos :

[TABLE]

where $f_{3}$ comes from the normalisation to $\Omega_{\mathrm{HI}}$. This prescription is commonly used for studies of 21 cm intensity mapping (amongst others, Villaescusa-Navarro et al., 2014; Sarkar et al., 2016; Seehars et al., 2016). $M_{\mathrm{min}}$ and $M_{\mathrm{max}}$ are the limits for a dark matter halo to host HI. They assume that only halos with 30 km/s $<v_{\mathrm{circ}}<$200 km/s host HI, which translates to lower and upper bounds, $M_{\mathrm{min}}$ and $M_{\mathrm{max}}$, through

[TABLE] 2. 2.

AGN : Nevertheless, Villaescusa-Navarro et al. (2016) measured the MHIMh relation in hydrodynamical simulations including AGN feedback and show that there is HI in halos that have $v_{\mathrm{circ}}>$200 km/s. They measured $M_{\mathrm{HI}}(M)=e^{\alpha}\,M^{\gamma}$ and fit $\alpha$ and $\gamma$ up to redshift 2. 3. 3.

DLA50 : A prescription adapted from DLA studies (Barnes & Haehnelt, 2010, 2014) by Padmanabhan et al. (2016)

[TABLE]

where $\alpha$ is the ratio of HI within halos and cosmic HI, $f_{\mathrm{H,c}}=(1-Y_{p})\Omega_{b}/\Omega_{m}$ is the cosmic hydrogen fraction with $Y_{p}$ the cosmological helium fraction by mass, and $v_{c}(M)$ is the virial velocity of a halo (Bullock et al., 2001) :

[TABLE]

with $\Delta_{v}$ the mean overdensity of the halo that we take to be 200. For DLAs, Padmanabhan et al. (2016) considered $v_{c,0}=50$ km/s and an infinite $v_{c,1}$. They fitted $\alpha$ to measurements between redshift 0 and 4 (column density distributions, biases, $\Omega_{\mathrm{HI}}$and the incidence rate). 4. 4.

21cm : Padmanabhan et al. (2016) adapted Eq. 17 to 21 cm IM observations using ad-hoc velocity cuts $v_{c,0}=30$ km/s and $v_{c,1}=200$ km/s. Similarly to the DLA50 model, Padmanabhan et al. (2016) fitted $\alpha$ on the same measurements. Note that in both latter cases the slope is fixed and equal to unity which is higher than what is measured in hydro-simulations. 5. 5.

HOD A : Padmanabhan & Refregier (2017) improved Eq. 17 by introducing a flexible slope, $\beta$, as well as the velocity cut-offs:

[TABLE]

where $\alpha$, $\beta$, $v_{c,0}$, and $v_{c,1}$ are free parameters and fitted on data measurements. 6. 6.

HOD B : Lastly, Padmanabhan et al. (2017) fitted an updated version of Eq. 19

[TABLE]

on all the available measurements including galaxy clustering. Their free parameters are $\beta$ and $\alpha$.

All these prescriptions are shown in Fig. 2 at $z=1$. They vary in shape, amplitude, and slope. Clearly the DLA50 scheme favours high halo masses as compared to the other models. We limit our analysis to $z=1$, the values of the free parameters are given in Table 1.

3 The HI power spectrum in real space

In this section, we compute the non-linear HI power spectrum in real space using all the above MHIMh models. We first describe the non-linear contributions to the power spectrum and show that the bias is neither constant nor linear on so-called linear scales. We discuss the implications of that effective bias on our understanding of the distribution of HI and compare this modeling approach to others.

3.1 A non linear bias on linear scales

The total HI power spectrum along with the non-linear contributions are shown in Fig. 3 for model HOD A. Contrary to our expectations, both $P_{\mathrm{HI}}^{22}$ and $P_{\mathrm{HI}}^{13}$ terms have significant contributions on linear scales. These contributions arise from the coupling of short and long wavelength modes. The $P_{\mathrm{HI}}^{13}$ term is negative and proportional to the matter power spectrum. Therefore, on linear scales, it lowers the amplitude of the HI power spectrum by $\sim 25$ % as compared to a standard biased power spectrum. Hence, the actual HI bias is lower than the linear HI bias. The $P_{\mathrm{HI}}^{22}$ term is constant on linear scales which induces a scale dependence of the HI bias on the largest scales. The flat contribution to the $P_{\mathrm{HI}}^{22}$ term (the dot-dashed line) is simply proportional to $b_{2}^{\mathrm{HI}}$ while, on those scales, the $P_{\mathrm{HI}}^{13}$ term is a function of $b_{1}^{\mathrm{HI}},b_{2}^{\mathrm{HI}}$, and $b_{3}^{\mathrm{HI}}$. The latter depend on the MHIMh prescriptions and are listed in Table 1. While the $b_{1}^{\mathrm{HI}}$s vary by $\sim 13\%$ amongst the different prescriptions, the variation strongly increases at higher orders. Indeed $b_{2}^{\mathrm{HI}}$s and $b_{3}^{\mathrm{HI}}$s differ by $35\%$ and $74\%$, respectively. Hence, the shape of the HI bias depends on the MHIMh prescription as shown in Fig. 4. Hereafter, we will callHI effective bias the following :

[TABLE]

The right panel of Fig. 4 shows the ratio $1/\overline{T}_{\mathrm{HI}}\sqrt{P_{\mathrm{HI}}(k)/P_{\mathrm{m}}^{\mathrm{NL}}(k)}$ where $P_{\mathrm{m}}^{\mathrm{NL}}(k)$ is the non-linear matter power spectrum computed using the same perturbation theory framework (See Appendix A). The normalisation by the mean HI temperature is to focus only on the bias and avoid additional amplitude variations. Indeed, the HI temperature is a function of $\Omega_{\mathrm{HI}}$, therefore of the MHIMh relation (see Eq. 1 and Table 1). On the largest scales, there is only a few percent difference between the different prescriptions. We won’t discuss them here as General Relativity corrections must be taken into account on ultra large scales. On non-linear scales, regardless of the MHIMh relation, we recover a bias well below 1, meaning that HI galaxies are highly anti-biased while they are only slightly on linear scales (Marín et al., 2010; Martin et al., 2012). At $k>0.2\,h\,$Mpc*-1* our biases are of the same order of magnitude than that of Sarkar et al. (2016) as shown in the right panel of Fig. 4. Their dip is deeper because their bias on linear scales is higher than ours. Nevertheless, they are consistent as explained in the following and in Sect. 3.3.

On linear scales, there is, at most, a difference of 10% in the amplitude of the models and 15% at $k=1\,h$ Mpc*-1*. Regardless of the MHIMh prescription, the effective bias $b_{\mathrm{HI}}^{\mathrm{eff}}$, is lower than the linear one. Their values at $k=0.01\,h$ Mpc*-1* are listed in Table 1 along with the associated linear biases. On linear scales, effective biases are always lower by 10-15% than their linear counterpart. They can be approximated, in real space, by

[TABLE]

It has consequences on our understanding of HI within the Large Scale Structure. The assumption that the measured HI bias on large scales is linear leads to an underestimation of that linear bias and therefore of the halo mass hosting HI. Hence HI lies in slightly more massive halos than thought.

3.2 In which halos does HI lie?

Currently, there is a tension between halo masses of HI-bearing systems, in particular between observations of HI galaxies at low redshift and those of DLAs at higher redshift as highlighted by Padmanabhan et al. (2016). They fitted all HI available measurements (DLA incidence rates, column densities, biases, HI fractional densities and biases) at several redshifts with both DLA- and 21 cm-based models (our schemes DLA50 and 21cm, amongst others). They showed that the 21 cm based model fitting all measurements systematically underpredicts the DLA bias. Similarly, DLA models, that are tuned to reproduce the DLA bias, always overpredict \Omega_{\mathrm{HI}}$$b_{\mathrm{HI}}. Their two DLA models have low velocity cut-offs of 50 and 90 km/s, which implies that there is no or only a low amount of neutral hydrogen in low mass halos. While Barnes & Haehnelt (2014) suggested that it could be caused by a strong stellar feedback, it remains inconsistent with HI observations at low redshift. In addition, the discrepancy holds when varying the HI concentration. It is important to note that the discrepancy is only at the level of the biases, hence, the tension is between the host halos of low and high redshift HI-bearing systems. Padmanabhan et al. (2016) argued that there must be a dramatic change in the properties of these systems over $0<z<3$ to have these halo masses on the same evolution path. This idea is strengthened by Padmanabhan & Refregier (2017) who introduced a halo occupation model inspired by both DLA and 21 cm emission framework (our HOD A model), using the same dataset as the former together with the HI mass function at $z\sim 0$. Again, most of the observables are relatively well fitted but the DLA bias is, again, underpredicted while the high mass part of the HI mass function is overpredicted. Lately, Padmanabhan et al. (2017) carried a similar analysis with an updated version of the MHIMh relation (our HOD B model) adding the 2-point correlation function (2PCF) of HI galaxies at small scales. By levering some degrees of freedom in the HI concentration, they did improve the overall quality of the fit but with an overpredicted 2PCF on large scales, a high mass tail of the HI mass function, and a too low DLA bias. It is clear that both models predict too many objects in high mass halos and a HI bias that is too high. Indeed, fits are driven towards high halo masses by the DLA bias measurement, which might be flawed as it is inconsistent with most observations and simulations. The latter statement seems inconsistent with our previous argument which is that HI lies in more massive halos than we think. The non-linear corrections to the bias are of the order of 15% at most while the discrepencies between the predicted and measured HI biases are, at least, of 50%. Therefore the systematic error due to the assumption of a linear bias on large scales is concealed by the error induced by the DLA bias.

We also consider a MHIMh relation adapted to DLAs, the DLA50 model. It is obvious from Fig. 2 that it favours higher mass halos as compared to any other prescription. This translates to a higher linear bias and to a change of sign of $b_{2}^{\mathrm{HI}}$ and $b_{3}^{\mathrm{HI}}$(see Table 1). $P_{\mathrm{HI}}^{13}$ becomes positive on large scales as shown in Fig. 3 which adds power to the HI power spectrum. Hence, the effective bias is higher than its linear counterpart in the case of DLAs. The amplitude of the power spectrum rises by 13% which translates in 7% on the effective bias for the DLA50 model. Of course, the additional power increases when going towards even higher mass halos. For instance, using a $v_{c,0}$ = 90 km/s instead of $v_{c,0}$ = 50 km/s, which translates to minimum halo masses of $10^{10.04}$ and $10^{10.80}\,h^{-1}\,M_{\odot}$ at $z=1$, leads to an increase of power of 24% and 11% at the power spectrum and bias levels, respectively. Thus, DLA models overpredict the HI bias even more than previously thought which enhances the tension between DLAs and 21 cm biases preventing any reconciliation.

3.3 Consistency with other modeling approaches and clustering analysis

It is the coupling between small and large scale modes that gives rise to an effective bias different from the linear one. Therefore, the mismatch exists for any tracer of dark matter. Hence one can wonder why it is not predicted by any other modeling approaches and why it has never been noticed in any clustering analysis. The answer to the first question is straightforward : models are constructed to predict a linear bias on linear scales. The procedure for modeling the clustering of any tracer is a distribution of halos coming either from the halo model or dark matter simulations which are filled in with the tracer. Therefore, only non-linearities coming from the evolution of the distribution of dark matter are present and not the ones coming from the distribution of the tracer, which is not the case in our approach. Indeed, it is the distribution of the tracer, the HI brightness temperature precisely, that has been perturbed. Therefore, our biases in Fig. 4 are consistent with that of Sarkar et al. (2016). They used a dark matter simulation in which they defined halos that are assigned an HI mass through the Bagla10 model. On the largest scales, they measure a HI bias of 0.92 fully consistent with the linear bias of 0.94 computed through Eq. 14.

Lastly, why a mismatch has not been noticed in clustering analysis yet as linearities are missing in current models? On linear scales, it is an offset of 15% on the bias for HI at most and the scale dependency is only of a few percent so it can be well concealed in the error bars. For instance, when fitting the parameters of a halo occupation distribution on a correlation function on both linear and non-linear scales, the halo mass thresholds hosting a central galaxy and one satellite galaxy would be found to be lower than it is in reality. Notwithstanding, this systematic error is lower than the statistical error on the fitted parameters in current clustering analysis. It won’t be the case with stage IV experiments such as Euclid and SKA. In addition, in the era of precision cosmology, ignoring these corrections will lead to flawed estimations of cosmological parameters.

4 The HI power spectrum in redshift space

In this section, we extend the previous analysis to the anisotropic power spectrum of HI in redshift space. We first adopt a theoretical point of view, investigating the power spectrum as a function of $k$ and $\mu$ to understand the effects of RSDs and second, we compare the expected linear power spectrum to the full one in the transverse and radial directions.

4.1 The HI effective bias on linear scales

We begin by investigating the different contributions to the HI power spectrum in the two extreme directions: in the transverse one where $\mu=0$, meaning that RSDs are null, and in the radial direction where $\mu=1$ and, hence, RSDs are maximal. Fig. 5 shows the different terms contributing to the HI power spectrum for the prescription HOD A. We recover similar behaviours to those in real space. The $P_{\mathrm{HI}}^{22}$ term is constant on large scales ($k<0.02\,h\,$Mpc*-1*) and rises towards small scales. The amplitude of the rise increases with $\mu$ as RSDs come in, they are contained in the $G_{2}$ and $K_{R}$ terms as shown in the lower panel. Therefore, on small scales ($k>0.2\,h\,$Mpc*-1*) non-linear contributions are maximal for $\mu=1$ where fingers of God are recovered. Again, the $P_{\mathrm{HI}}^{13}$ term is negative and thus removes power to the HI power spectrum on linear scales. The amplitude of the removal decreases with $\mu$: at the power spectrum level it lowers from $\sim 25\%$ at $\mu=0$ to $\sim 13\%$ at $\mu=1$, respectively. The former is similar to the real space case. The effective bias in redshift space is

[TABLE]

The expression of $I_{R}$ is given in Appendix A and it is worth noticing that the effective bias is also a function of the growth factor. We will explore this in more details in future work.

For $\mu\neq 0$ the effective bias cannot be computed directly because of RSD effects so we compare the full HI power spectrum to the linear Kaiser prediction $P_{\mathrm{HI}}^{\mathrm{11}}(k,\mu)=\overline{T}^{2}\left[b_{1}+f\,\mu^{2}\right]^{2}\,P_{\mathrm{m}}^{11}(k)$ in Fig. 6 for all MHIMh models. Regardless of the scale, they lead to ratios that are within 10% and those differences lowers with $\mu$. On linear scales, the effective bias gets closer to the linear one as $\mu$ increases. We can also notice in the lower panel that RSD effects impact the power spectrum only at $\mu>0.2$. On smaller scales, the rise is due to non-linear effects, only, at $\mu=0$ and also to RSDs for $\mu>0$. The slope of the rise scales with $\mu$ and it is exactly over the BAO scale range.

4.2 A scale dependent HI bias on BAO scales

We will carry on the analysis using only the model HOD A. The bottom panel of Fig. 6 shows a scale dependence of the HI bias that is enhanced by RSD effects over the BAO scale range : the bias rises from 10 % at $\mu=0$ to a factor 2 at $\mu=1$. To adopt an observational point of view, we change the coordinates to transverse and radial directions in Fig. 7. The two top panels show the linear and total HI power spectra. At first glance, non-linear terms shift the turnover of the power spectrum towards higher $k_{\perp}$ ($>0.1\,h\,$Mpc*-1*). They slightly enhance the signal along the $k_{\perp}$ direction while it is boosted along $k_{\parallel}$ by RSD effects. The lower panel of Fig. 7 shows the ratio between the linear and total HI power spectra. On large scales, $k_{\parallel},k_{\perp}<0.01\,h\,$Mpc*-1*, we recover a maximum ratio of 25%. At small $k_{\parallel}$ and towards large $k_{\perp}$, non-linearities increase the amplitude of the HI power spectrum by a factor 2 while at large $k_{\parallel}$ RSD effects dominate non-linear ones and make any $k_{\perp}$ dependence vanish. Over the BAO scale range, the HI power spectrum increases by a factor of 5 and 2 in the radial and transverse directions, respectively. Therefore, both non-linearities and RSD effects modify the ratios between the BAO peaks. It is therefore necessary to take non-linearities into account when estimating cosmological parameters. To circumvent the contamination by non linear effects, one would preferentially measure the BAO peaks in the transverse direction but Villaescusa-Navarro et al. (2017) showed that, in single dish mode, beyond a certain size, the beam of the instrument smears the wiggles out in the transverse direction and that BAOs can only be detected in the radial direction. This is a limitation for the SKA and Meerkat but not for BINGO, CHIME or HIRAX as they will have a higher angular resolution.

5 Conclusion

Radio telescopes are about to open a new window of observation on the Universe, in particular, using 21 cm intensity mapping. We investigate the non-linear power spectrum of HI in both real and redshift space and in the light of the relation between the halo mass and the HI mass.

Our main result is that on linear scales, the HI bias is not constant but scale-dependent. Using a full 1-loop development in perturbation theory of the power spectrum of HI we show that non-linear contributions remove power to the HI power spectrum on linear scales in both real and redshift space at $z=1$. This result is contrary to our expectations and is not found in other modeling approaches. Commonly, a distribution of dark matter halos is ‘painted’ with a baryonic tracer, so that only non-linearities coming from the distribution of dark matter are taken into account and not the ones coming from the evolution of the distribution of the tracer.

In real space, the effective bias of HI is 10-15% lower than the linear one, depending on the MHIMh relation. The assumption that the observed HI bias is linear underpredicts the actual linear bias, and hence, the mass of halos hosting HI.

In redshift space, the effective bias is also lower than its linear counterpart up to 15% and its scale dependence is highly sensitive to RSD effects. Over the BAO scale range, the HI bias rises with a slope that steepens with $\mu$. Regardless of the MHIMh prescription, the difference between the linear and the full HI power spectra reaches a factor 5 which can lead to a modification of the ratios between BAO peaks. Therefore, it will be crucial to take non-linearities into account when estimating cosmological parameters.

The different MHIMh relations lead to variations of 15% at most on the HI bias. It is within the error bars on any of the current HI bias measurements so it is not an issue at the moment. Nevertheless, it will be indispensable for the upcoming HI surveys. It is worth noting that the observable is the product $T_{\mathrm{HI}}\,b_{\mathrm{HI}}$ where the HI temperature is also a function of the MHIMh relation through $\Omega_{\mathrm{HI}}$. This product differs up to a factor of 7 between the different prescriptions.

Thorough forecasts of the effect of non-linearities on the estimation of BAO peaks and in a broader way, cosmological parameters, are required, including the redshift evolution as HI is positively biased at higher redshift, therefore non-linearities add power to the power spectrum of HI. Lastly, this effect is not only present in HI intensity mapping surveys but in any galaxy surveys.

Acknowledgments

The authors would like to thank Roy Marteens, Chris Clarkson, Jose Fonseca, and Vincent Desjacques for useful discussions. They also would like to thank Debanjan Sarkar for providing them the data points of their paper. AP, OU and MGS acknowledge support from the South African Square Kilometre Array Project as well as from the National Research Foundation.

Appendix A Perturbation Theory formulas

The necessary kernels are : $F_{2}$ the non-linear density kernel, $G_{2}$ induced by peculiar velocities at second order, and $K_{R}$ arises from non-linear mode coupling (velocity-velocity and velocity-density, Bernardeau et al., 2002) {IEEEeqnarray}rCl F_2(k_1,k_2) &= 57 + 12 k1⋅k2k1k2 [ k1k2 + k2k1] + 27 [ k1⋅k2k1k2 ]^2

G_2(k_1,k_2) = 37 + 12 k1⋅k2k1k2 [ k1k2 + k2k1] + 47 [ k1⋅k2k1k2 ]^2

K_R(k_1,k_2) = f  b_1 μ_1^2 + f b_1  μ_2^2

+μ_1μ_2 [ f  b_1k1k2 +f  b_1 k2k1 ]

+f^2 [ 2μ_1^2 μ_2^2 + μ_1μ_2( μ_1^2k1k2 + μ_2^2k2k1 ) ]

The matter and velocity power spectra at third order write

[TABLE]

The last component of the $P_{\mathrm{HI}}^{13}(k,\mu)$ is

[TABLE]

with

[TABLE]

Lastly, in this framework, the full 1-loop matter power spectrum in real space is

[TABLE]

where the second order term writes

[TABLE]

In redshift space, the full 1-loop matter power spectrum is

[TABLE]