The Breakdown Scale of HI Bias Linearity
Zhenyuan Wang, Yangyao Chen, Yi Mao, Houjun Mo, Huiyuan Wang, Hong Guo, Cheng Li, Jian Fu, Yipeng Jing, Jing Wang, Xiaohu Yang, Zheng Zheng
DDraft version November 23, 2020Typeset using L A TEX twocolumn style in AASTeX63
The Breakdown Scale of H I Bias Linearity
Zhenyuan Wang ,
1, 2
Yangyao Chen ,
1, 3
Yi Mao , Houjun Mo ,
3, 1
Huiyuan Wang , Hong Guo , Cheng Li , Jian Fu, Yipeng Jing , Jing Wang , Xiaohu Yang , and Zheng Zheng Department of Astronomy, Tsinghua University, Beijing 100084, China Department of Astronomy and Astrophysics, The Pennsylvania State University, University Park, PA 16802, USA Department of Astronomy, University of Massachusetts, Amherst MA 01003-9305, USA Key Laboratory for Research in Galaxies and Cosmology, Department of Astronomy, University of Science and Technology of China, Hefei, Anhui 230026, China Key Laboratory for Research in Galaxies and Cosmology, Shanghai Astronomical Observatory, Shanghai 200030, China Department of Astronomy, and Tsung-Dao Lee Institute, Shanghai Jiao Tong University, Shanghai 200240, China Kavli Institute for Astronomy and Astrophysics, Peking University, Beijing 100871, China Department of Physics and Astronomy, University of Utah, 115 South 1400 East, Salt Lake City, UT 84112, USA (Received 27-Apr-2020; Revised 08-Nov-2020; Accepted 16-Nov-2020)
Submitted to ApJABSTRACTThe 21 cm intensity mapping experiments promise to obtain the large-scale distribution of H I gas at thepost-reionization epoch. In order to reveal the underlying matter density fluctuations from the H I mapping, it isimportant to understand how H I gas traces the matter density distribution. Both nonlinear halo clustering andnonlinear effects modulating H I gas in halos may determine the scale below which the H I bias deviates fromlinearity. We employ three approaches to generate the mock H I density from a large-scale N-body simulationat low redshifts, and demonstrate that the assumption of H I linearity is valid at the scale corresponding to thefirst peak of baryon acoustic oscillations, but breaks down at 𝑘 (cid:38) . ℎ Mpc − . The nonlinear effects of haloclustering and H I content modulation counteract each other at small scales, and their competition results in amodel-dependent “sweet-spot” redshift near 𝑧 =1 where the H I bias is scale-independent down to small scales.We also find that the linear H I bias scales approximately linearly with redshift for 𝑧 ≤ Keywords:
H I line emission (690), Line intensities (2084), Galaxy dark matter halos (1880), Large-scalestructure of the universe (902) INTRODUCTIONNeutral hydrogen (H I ) atoms, which are expected to be con-tained in halos at low redshifts (0 . (cid:46) 𝑧 (cid:46) (Chen2012), CHIME (Bandura et al. 2014), HIRAX (Newburghet al. 2016), BINGO (Battye 2013), and SKA (Pritchard et al.2015), which will survey the H I mass distribution in very large Corresponding author: Yi [email protected] http://tianlai.bao.ac.cn https://chime-experiment.ca https://hirax.ukzn.ac.za volumes, provide a promising way to constrain the expansionhistory and structure formation in the Universe, thereby un-veiling the nature of dark energy.These 21 cm intensity mapping experiments, despite lowangular resolutions, can be used to detect large-scale featuresin the cosmological density field (Chang et al. 2008; Loeb &Wyithe 2008). For this purpose, it is important to understandhow accurately H I gas traces the matter density fluctuations.In general, the power spectrum of the H I gas distribution isrelated to that of the underlying matter through a bias relation, 𝑃 HI ( 𝑘 ) = 𝑏 𝑃 m ( 𝑘 ) , where 𝑏 HI is the bias factor. It is, there-fore, necessary to understand the bias factor, 𝑏 HI , in particularits scale dependence, in order to use 𝑃 HI ( 𝑘 ) to infer the dis-tribution of mass in the universe. Note that the measurementof the baryon acoustic oscillations (BAO) can be obtained byusing a template of wiggles in the power spectrum, which isleast sensitive to the nonlinear bias. But the nonlinear bias a r X i v : . [ a s t r o - ph . C O ] N ov Z. Wang et al.can affect the broadband shape of power spectrum which alsocontains a wealth of cosmological information. In particular,it is important to determine the breakdown scale below whichthe H I bias deviates from linearity, which is the focus of thispaper. At quasi-linear scales, large-scale structure perturba-tion theory (see Desjacques et al. 2018; d’Amico et al. 2020and references therein), which incorporates the higher-orderbias parameters, may be developed to model the nonlinear H I clustering (e.g. Modi et al. 2019).After cosmic reionization, most H I gas is expected to bein galaxies, thanks to their high density and low temperature,while the neutral fraction in the intergalactic medium is verylow, about 10 − . Furthermore, fluctuations in the ionizationfield are not expected to affect the H I power spectrum onlarge scales (Wyithe & Loeb 2009). Thus, the distributionof the H I gas may be understood in terms of its relationwith galaxies, or with dark matter halos in which galaxiesreside (Cai et al. 2016, 2017; Cui et al. 2017). Gas andstar-formation processes can, in principle, change the H I gasdistribution in dark matter halos, and potentially introducenonlinear bias in the relationship between H I gas and darkmatter (Guo et al. 2020). In addition, it is well-known fromN-body simulations that the distribution of dark matter halostraces the underlying matter distribution nonlinearly at smallscales (Jeong & Komatsu 2009; Nishizawa et al. 2013). Thesenonlinearities, albeit at small scales (i.e., the size of halos),might spoil the H I linearity assumption even on large scales,because of mode coupling on different scales.Previous studies of H I bias either employed oversimpli-fied H I -halo mass relation (similar to the fitting formula inKhandai et al. 2011) applied to N-body simulations (Baglaet al. 2010; Guha Sarkar et al. 2012; Sarkar et al. 2016;Padmanabhan et al. 2016; Padmanabhan & Refregier 2017;Padmanabhan et al. 2017; Sarkar & Bharadwaj 2018), ormodelled the H I gas using hydrodynamic simulations, suchas IllustrisTNG (Villaescusa-Navarro et al. 2018), Illustrisand Osaka (Ando et al. 2019). However, the volumes ofgas simulations, typically (cid:46) ( ℎ − Mpc ) , are usually toosmall to be valid on BAO scales ( ∼ ℎ − Mpc).Given its importance, in this paper, we study the relation-ship between H I gas and dark matter on large scales, us-ing three – empirically, numerically, and observationally ori-ented, respectively – approaches to model H I gas in halos ofdifferent masses, and using halos in a large N-body simulationto construct the H I gas distribution on large scales. Our simu-lation volume, ( ℎ − Mpc ) , is sufficiently large so that thefinite box effect on the power spectrum and bias is negligibleon BAO scales (Klypin & Prada 2019). The use of differentmodels for HI gas in halos also allows us to draw genericconclusions that are independent of our ignorance about thedetails of galaxy formation in dark matter halos. The rest of this paper is organized as follows. In Section 2,we describe the methodology of modelling the H I gas. Weshow the results and discussions in Section 3, and make con-cluding remarks in Section 4. MOCKING THE H I GAS DISTRIBUTIONOur H I mock data is constructed from the results of a large-scale, high-resolution N-body simulation, ELUCID (Wanget al. 2016), of the Λ CDM universe, performed with theL-Gadget code, a memory-optimized version of Gadget-2(Springel 2005), in a comoving volume of 500 ℎ − Mpc oneach side using 3072 particles. We refer the readers toWang et al. (2016) for details of this simulation. To findhalos, we use the FoF algorithm with a linking length of 0.2times the mean particle separation. The SUBFIND algorithm(Springel et al. 2001) is employed to resolve the sub-structures(i.e. subhalos) in each FoF halo and to build the merger trees.We adopt an empirical model (Lu et al. 2014) to construct thestar formation histories of galaxies in those halos with massesabove 10 ℎ − 𝑀 (cid:12) (about 30 N-body particles). To fullytrace the star formation history, we develop a Monte Carlomethod to append unresolved progenitors to the leaf-halos ofthe merger tree (Chen et al. 2019). The H I gas is then assignedto halos with masses above 10 ℎ − 𝑀 (cid:12) using a star formationmodel (Krumholz et al. 2008, 2009a,b; Krumholz 2013) thatprovides the full information about the star formation history.Finally, the H I gas is smoothed onto grids to compute theH I power spectrum. The key ingredients of our method aredetailed below. The background cosmology is consistent withthat given by the WMAP five-year data (Dunkley et al. 2009): Ω 𝑚 = . Ω Λ = . Ω 𝑏 = . ℎ = . 𝑛 𝑠 = . 𝜎 = .
8. 2.1.
Star formation history
For resolved halos with 𝑀 h ≥ ℎ − 𝑀 (cid:12) , we follow theempirical model for star formation rate (SFR) as describedin Lu et al. (2014) (their “Model III”). The SFR of a central galaxy is assumed to depend only on the mass of its host halo, 𝑀 h , and redshift 𝑧 ,SFR ( 𝑀 h , 𝑧 ) = 𝜀 𝑓 𝑏 𝑀 h 𝜏 ( 𝑋 + ) 𝛼 (cid:18) 𝑋 + 𝑅𝑋 + (cid:19) 𝛽 (cid:18) 𝑋𝑋 + 𝑅 (cid:19) 𝛾 . (1)Here 𝜀 is the overall efficiency, 𝑓 𝑏 = Ω 𝑏 / Ω 𝑚 is the cosmicbaryon fraction, 𝜏 = [ /( 𝐻 )] ( + 𝑧 ) − / describes thedynamical timescale of halos at a redshift 𝑧 , the variable 𝑋 ≡ 𝑀 h / 𝑀 𝑐 where 𝑀 𝑐 is a characteristic mass scale. Othervariables are parametrized as 𝛼 = 𝛼 ( + 𝑧 ) 𝛼 (cid:48) , and 𝛾 = 𝛾 a if 𝑧 < 𝑧 𝑐 , or, otherwise, 𝛾 = ( 𝛾 a − 𝛾 b ) [( 𝑧 + )/( 𝑧 𝑐 + )] 𝛾 (cid:48) + 𝛾 b .The free parameters ( 𝜀 , 𝑅 , 𝑀 𝑐 , 𝛼 , 𝛼 (cid:48) , 𝛽 , 𝛾 a , 𝛾 b , 𝛾 (cid:48) , 𝑧 𝑐 ) canbe found by fitting the observed galaxy stellar mass functionsand a composite local cluster conditional galaxy luminosityfunction at the 𝑧 -band, as shown in Lu et al. (2014) (their Table I Bias Linearity 33). For unresolved halos with 𝑀 h < ℎ − 𝑀 (cid:12) , MonteCarlo trees are adopted to extend their assembly historiesdown to 10 ℎ − 𝑀 (cid:12) (Chen et al. 2019).This model (Lu et al. 2014) assumes that, during galaxymergers, the SFR is under exponential decay in satellite galax-ies where the gas can be stripped. As such, the H I gas isdominated by the contributions from central galaxies. Whilethis may not be true for big halos (Villaescusa-Navarro et al.2018), we neglect the H I gas from satellite galaxies, for sim-plicity.With empirical star formation and merger models, we cantrace the mass growth of each central galaxy from its mergertree, and obtain its stellar mass 𝑀 ∗ . For a given halo mass,the stellar mass 𝑀 ∗ may not be the same in different halosbecause of their different merger histories.2.2. Star formation model
To connect the surface density of SFR (cid:164) Σ ∗ and that of gasmass Σ 𝑔 , we follow the star formation model developed inKrumholz et al. (2008, 2009a,b); Krumholz (2013), (cid:164) Σ ∗ = 𝑓 H 𝜖 ff Σ g t ff , (2)where 𝜖 ff = 0.01, t ff = [ Σ 𝑔 /( 𝑀 (cid:12) pc − )] − . Myr. Assum-ing that the gas is cold and comprised of H and H I , the H fraction is given by 𝑓 H = (cid:40) − ( 𝑠 + . 𝑠 ) , if 𝑠 (cid:54) , otherwise (3)The variable 𝑠 = ln ( + . 𝜒 + . 𝜒 )/( . 𝜏 c ) , where 𝜏 𝑐 = 𝑐 𝑍 𝑜 Σ 𝑔 /( g cm − ) , and the clumping factor 𝑐 = . 𝑍 𝑜 , we adopt the average metallicity-stellar mass rela-tion from the FIRE simulation (Ma et al. 2015), log 𝑍 𝑜 = . [ log ( 𝑀 ∗ / 𝑀 (cid:12) ) − ] + .
93 exp (− . 𝑧 ) − .
74. Theradiation field parameter 𝜒 is estimated (Krumholz 2013)as 𝜒 = 𝐺 (cid:48) / 𝑛 CNM , where 𝐺 (cid:48) = (cid:164) Σ ∗ / (cid:164) Σ ∗ , , (cid:164) Σ ∗ , = . × − 𝑀 (cid:12) pc − Myr − , and 𝑛 CNM is the density of cold neu-tral medium (CNM) in units of cm − . In molecular-poorregions, the CNM density is 𝑛 CNM , hydro ≈ Σ g /( 𝑀 (cid:12) pc − ) ,while in molecular-rich regions, the CNM density is 𝑛 CNM , = 𝐺 (cid:48) /[( . / . )( + 𝑍 . 𝑜 )] . In general, 𝑛 CNM = max { n CNM , , n CNM , hydro } .2.3. Disk size
To connect the surface density and the total density, weassume that the gas surface density follows an exponen-tial profile, Σ 𝑔 ( 𝑟 ) = Σ e − 𝑟 / 𝑅 𝑔 . We assume the gas diskto stellar disk size ratio 𝑅 𝑔 / 𝑅 ∗ = . 𝑅 𝑔 / 𝑅 ∗ = . Figure 1.
The HI-halo mass relation derived from different modelsat 𝑧 =
0. We show the results using the LK model (blue), the TKmodel (magenta), and the AH model (green). Here we also includethe scatter points (gray dots) and 1 𝜎 envelope (blue dashed lines)for the LK model, and the error bars for the AH model. 𝑧 ≈ . 𝑅 ∗ ( 𝑀 ∗ ) = 𝑅 ( 𝑀 ∗ / 𝑀 ) . (cid:2) ( / ) + ( / )( 𝑀 ∗ / 𝑀 ) . (cid:3) ( . / . ) , where 𝑅 = . kpc, 𝑀 = . 𝑀 (cid:12) . The disk size evolveswith redshift as 𝑅 ∗ ( 𝑧, 𝑀 ∗ ) = 𝑅 ∗ ( 𝑀 ∗ ) [( + 𝑧 )/ . ] − . .2.4. H I -halo mass relation In our above modelling, for a fixed stellar mass 𝑀 ∗ , a givenvalue of disk central density Σ determines Σ 𝑔 ( 𝑟 ) at someradius in the disk. The aforementioned star formation modelis employed to solve for (cid:164) Σ ∗ ( 𝑟 ) numerically from Σ 𝑔 ( 𝑟 ) , whichgives the H I surface density Σ HI ( 𝑟 ) . By integrating over thedisk, we can find a correlation between the SFR and theH I mass for a central galaxy, given 𝑀 ∗ . For each halo, wecompute the SFR using the aforementioned empirical model,and 𝑀 ∗ from halo merger history. Finally, the H I mass iscomputed by interpolation using its correlation with SFR. OurH I gas model, which incorporates the empirical SFR model(Lu et al. 2014) and the star formation model (Krumholzet al. 2008, 2009a,b; Krumholz 2013), is dubbed “LK model”,which stands for “ L u et al. + K rumholz et al. model”.To test the model dependence of H I bias, we also as-sign the H I mass inside a halo by using the average H I -halo mass relation obtained from two other approaches.One approach uses the IllustrisTNG simulation (their gasdata)(Villaescusa-Navarro et al. 2018) and the same starformation model(Krumholz et al. 2008, 2009a,b; Krumholz2013). This model is dubbed “TK model” herein, whichstands for “Illustris T NG + K rumholz et al. model”. Theaverage H I -halo mass relation in the other approach was ob-tained by using the updated measurements of ALFALFA sur-vey and HOD model (Guo et al. 2017) (only available at 𝑧 = A LFALFA data + H OD model”. Following the fitting for-mula of average H I -halo mass relation in Villaescusa-Navarro Z. Wang et al. Table 1.
Parameter values used for the TK and AH models.Model 𝑧 𝛼 𝑀 [ 𝑀 (cid:12) / ℎ ] 𝑀 min [ 𝑀 (cid:12) / ℎ ] TK 0 0.24 4 . × . × . × . × . × . × . × . × AH 0 0.12 2 . × . × et al. (2018), we use the following expression for both TK andAH models, 𝑀 HI ( 𝑀 h , 𝑧 ) = 𝑀 (cid:18) 𝑀 h 𝑀 min (cid:19) 𝛼 exp (cid:34) − (cid:18) 𝑀 min 𝑀 h (cid:19) . (cid:35) . (4)The bestfit parameter values, as listed in Table 1, are takenfrom Villaescusa-Navarro et al. (2018) (their Table 1 for FoFhalos) for the TK model, and obtained by 𝜒 -fitting the 𝑀 HI - 𝑀 ℎ data at 𝑧 = I -halo mass relation for centralgalaxies at 𝑧 =
0. Our results (LK model) are compared withpredictions from the IllustrisTNG simulation (TK model),and the results from updated ALFALFA observations (AHmodel). All results agree well for low-mass halos ( 𝑀 h < ℎ − 𝑀 (cid:12) ). We checked that this agreement holds well athigher redshifts (0 < 𝑧 <
2) between LK and TK models. Formassive halos, nevertheless, our model underestimates the H I mass, for two possible reasons. First, the H I mass in theTK model includes the contributions from both central andsatellite galaxies, while both our model and the AH modelonly consider those from the central galaxies. Secondly, ourempirical model might underestimate the SFR for massivehalos. However, the contribution of H I gas from massivehalos is generally not important due to the sharp decrease ofthe halo mass function towards the massive end. In addition,the slope of H I -halo mass curve declines at the high mass end,which further suppresses the contribution of H I gas insidethe massive halos. We will further discuss the impact ofH I modelling in the high-mass end on the linear H I bias inSection 3.2 below.2.5. H I Power spectrum
The H I mass in each halo is smoothed onto a uniform gridwith 1024 cells, and we compute the H I power spectrum fromthe FFT. We only keep the power spectrum for wavenumberless than a quarter of Nyquist number ( 𝑘 < . ℎ Mpc − )to avoid the alias effect. In Fourier space, we can define ascale-dependent effective bias, 𝑏 HI ( 𝑘 ) , 𝛿 HI ( 𝒌 ) = 𝑏 HI ( 𝑘 ) 𝛿 m ( 𝒌 ) + 𝜖 ( 𝒌 ) , (5)where 𝜖 ( 𝒌 ) is a stochastic component which does not correlatewith the density field, 𝛿 m . On large scales, we expect 𝑏 HI is ascale-independent linear bias. The H I bias can be estimated using the auto-power spectrumof H I gas, 𝑏 uncorrHI , auto ( 𝑘 ) = [ 𝑃 HI ( 𝑘 )/ 𝑃 m ( 𝑘 )] / , if the shot noiseis uncorrected. The leading-order mass-weighted H I shotnoise is estimated by shuffling H I gas randomly, i.e. 𝑃 SN = 𝑉 − (cid:104) 𝜖 ( 𝒌 ) 𝜖 (− 𝒌 )(cid:105) , and then subtracted from the raw powerspectrum. After correcting for shot noise, we have 𝑏 HI , auto ( 𝑘 ) = √︄ 𝑃 HI ( 𝑘 ) − 𝑃 SN 𝑃 m ( 𝑘 ) . (6)The assumption of H I linearity can be tested by checking if theH I bias, 𝑏 HI , auto ( 𝑘 ) , is equal to the scale-independent linearbias at large scales. Of course, the H I bias is expected to bescale-dependent at small scales due to nonlinear evolution.The H I bias may also be estimated using the cross-powerspectrum between H I density and total matter density, 𝑏 HI , cross ( 𝑘 ) = 𝑃 HI , m ( 𝑘 ) 𝑃 m ( 𝑘 ) . (7)This estimator avoids the shot noise automatically. However,in this paper, we choose to estimate the H I bias based on theauto-power spectrum of H I gas, because the 21 cm intensitymapping measures the auto-power spectrum of the 21 cmbrightness temperature. As shown below in Section 3.3, theresults from these two estimators are in good agreement. Thuswe neglect the subscript “auto” throughout this paper exceptin Section 3.3. RESULTS AND DISCUSSION3.1.
Generic behavior
In Figure 2, we show the H I bias from different H I -halo massrelations at different redshifts (except that the AH model isonly at 𝑧 =
0) as well as the halo bias. In all three models, theH I bias remains a constant at large scales for 𝑘 (cid:46) . ℎ Mpc − ,i.e. we confirm that, generically, H I gas is indeed a linearbiased tracer at the first BAO peak. However, the linearityassumption begins to break down at the second BAO peak. Totest whether this break-down scale relies on the halo resolutionin our simulation, we vary the minimum halo mass from10 ℎ − 𝑀 (cid:12) to 10 ℎ − 𝑀 (cid:12) , and find that while the amplitudeof H I bias depends on the halo mass cutoff similar to that of thehalo bias, the linearity break-down scale is almost unchanged.Also, to test the effect of satellite galaxies, we estimate the H I masses from satellite galaxies and assign them to the centersof subhalos, using the LK model at 𝑧 =
0. We find thatincluding satellites does not change the shape of the H I powerspectrum significantly on scales 𝑘 (cid:46) ℎ Mpc − .The behaviors at small scales are more interesting, as mostof the H I gas resides only inside halos after cosmic reioniza-tion. Figure 2 shows that nonlinear halo clustering alwaysenhances the halo power spectrum at small scales (before cor-rected for shot noise). However, Figure 1 shows that H I mass I Bias Linearity 5 − − . . . . b ( k ) z = 0 halo, M HI (LK model)HI (TK model) HI (AH model)1st BAO peak2nd BAO peak − − . . b ( k ) z = 1 halo, M HI (LK model)HI (TK model) − − k [ h /Mpc] . . b ( k ) z = 2 halo, M HI (LK model)HI (TK model) − − k [ h /Mpc] . . . b ( k ) z = 3 halo, M HI (LK model)HI (TK model)
Figure 2.
The bias of halo mass density fluctuations (red) and HI mass density fluctuations derived from the LK(blue), TK(magenta), andAH(green) models at 𝑧 =
0, 1, 2, and 3, respectively, with respect to the matter density fluctuations, with shot-noised corrected (thick solid lines)and uncorrected (thin solid lines). The dashed lines indicate the constant linear bias which is estimated by averaging over 𝑘 = .
025 — 0.075 ℎ Mpc − (we neglect the smallest 𝑘 -mode due to its relatively large cosmic variance). The dot-dashed vertical lines mark the wavenumbers ofthe first (black) and second (grey) BAO peaks. is suppressed in large halos. This suppression decreases theH I density fluctuations at small scales relative to the level offluctuations caused by halos (see Fig. 2). The H I suppressioneffect is stronger at lower redshifts as more massive halosform. The competition between these two opposite effects,namely the nonlinear effects in halo clustering and those mod-ulating the H I gas in halos, determines the evolution of theH I bias at small scales. As shown in Figure 3, for both LKand TK models, the H I bias at small scales is enhanced withrespect to the linear bias at high redshifts, just like the non-linear halo bias, while the H I bias is actually suppressed atsmall scales at 𝑧 = 𝑘 (cid:38) . ℎ Mpc − (Jeong & Komatsu 2009; Nishizawa et al.2013) from N-body simulations. Naively, this sets the genericscale for the breakdown of linearity in H I bias, since mostof the H I gas resides inside halos after cosmic reionization.Nevertheless, the nonlinearity of the H I content significantlyaffects the level of H I fluctuations with respect to halo bias-ing, thereby modulating the breakdown scale and making itredshift-dependent, as shown in Figure 3. In particular, thesetwo nonlinear effects appear to balance each other at a tran-sition time where the H I bias is linear down to small scales.In the LK model, this “sweet-spot” redshift is near 𝑧 = . 𝑘 (cid:39) . ℎ Mpc − .In the TK model, the transition takes place at 𝑧 ≈
1, withthe linearity extending down to 𝑘 (cid:39) . ℎ Mpc − . Thus, the“sweet-spot” redshift is likely to be near 𝑧 =
1, although theexact value is model-dependent. − − . . . . . . . b H I ( k ) LK model z =3.0 z =2.0 z =1.2 z =1.0 z =0.5 z =0.0 − − k [ h /Mpc] . . . . . . . b H I ( k ) TK model
Figure 3.
The redshift evolution of the HI bias from 𝑧 = 𝑘 = .
025 — 0.075 ℎ Mpc − . The arrow marksthe scale at which the HI bias deviates from the linear bias at the1.5% level. The dot-dashed vertical lines mark the scales of the first(black) and second (grey) BAO peaks. Linear H I bias Z. Wang et al. . . . . . . . z . . . . b H I , li n e a r T K m o d e l L K m o d e l − − k [ h /Mpc] . . . . . P H I ( k , z ) / P H I ( k , z = ) TK z =0 z =1 z =2 z =3 10 − − k [ h /Mpc] . . . . . P H I ( k , z ) / P H I ( k , z = ) LK z =0 z =1 z =2 z =3 Figure 4.
The redshift evolution of the HI linear bias in the LKmodel (blue dots) and the TK model (magenta dots). We fit thedata linearly between 𝑧 = 𝑧 = 𝑃 HI ( 𝑘, 𝑧 )/ 𝑃 HI ( 𝑘, 𝑧 = ) for both models ininsets.
10 11 12 13 14 15 log M h [ M (cid:12) /h ] . . . . . ∆ b H I / ∆ l og M h LKTKAH
Figure 5.
The contribution to the linear HI bias from differentlogarithmic halo mass bin for the LK (blue), TK (magenta), and AH(green) model at 𝑧 = The linear H I bias (i.e. the constant H I bias averaged overlarge scales) increases with redshift, as shown in Figure 4.We find an interesting feature in both LK and TK models. Ingeneral, the H I bias varies approximately linearly with red-shift. This linear relation is almost exact between 𝑧 = <
10% for 𝑧 < <
15% for 2 < 𝑧 < I bias canbe written as 𝑏 HI , linear ( 𝑧 ) = (cid:2) 𝐷 HI ( 𝑧 )/ 𝐷 m ( 𝑧 ) (cid:3) 𝑏 HI , linear ( ) ,where 𝐷 HI and 𝐷 m are the linear growth functions of the H I and matter density fluctuations, respectively, i.e. 𝐷 HI ( 𝑧 ) = (cid:2) 𝑃 HI ( 𝑧 )/ 𝑃 HI ( ) (cid:3) / and 𝐷 m ( 𝑧 ) = (cid:2) 𝑃 m ( 𝑧 )/ 𝑃 m ( ) (cid:3) / . Asshown in the insets of Figure 4, the H I density power spec-trum varies only slightly with redshift, i.e. 𝐷 HI ( 𝑧 ) ≈
1. Thesimilar result was also found in Villaescusa-Navarro et al.(2018). The reason that H I clustering only weakly variesat 0 < 𝑧 < 𝐷 m ( 𝑧 ) ∝ ( + 𝑧 ) − (Cooray & Sheth 2002).These two effects combined lead to the linear scaling rela-tion, 𝑏 HI , linear ( 𝑧 ) ∝ ( + 𝑧 ) , which we find to be generic. There are two reasons why this relation is not exactly linear.First, 𝐷 m ( 𝑧 ) is suppressed at 𝑧 < I power spectrum has small, non-monotonous,evolution with redshift. As an illustration, consider a case inwhich 𝐷 HI ( 𝑧 ) =
1, but 𝐷 m ( 𝑧 ) takes the value from the linearperturbation theory (including the effect of dark energy). Wefind that the prediction of the linear H I bias in this case agreeswith the actual results in both models, with <
15% error. Thisis consistent with the fact that the H I power spectrum reachesits maximum at 𝑧 (cid:39) −
2, with the values at 𝑧 = I bias can be model-dependent. Fig-ures 2–4 show that in general the TK model predicts a highervalue of linear H I bias than the LK and AH models. This dif-ference might be attributed to the contributions of the H I gasin massive halos. We can understand this with halo model,in which the linear bias can be written as the integration ofcontributions from halos with different mass, 𝑏 HI , linear ( 𝑧 ) = ∫ 𝑀 max 𝑀 min 𝑛 ( 𝑀 h , 𝑧 ) 𝑏 ( 𝑀 h , 𝑧 ) 𝑀 HI ( 𝑀 h , 𝑧 ) 𝑑𝑀 h ∫ 𝑀 max 𝑀 min 𝑛 ( 𝑀 h , 𝑧 ) 𝑀 HI ( 𝑀 h , 𝑧 ) 𝑑𝑀 h . We calculate the prediction of linear H I bias in halo modelusing the fitting formula of the halo bias 𝑏 ( 𝑀 h , 𝑧 ) and haloabundance 𝑛 ( 𝑀 h , 𝑧 ) in Tinker et al. (2008, 2010), and theaverage H I -halo mass relation for all three models at 𝑧 =
0, andfind the results agree quite well with the bias directly measuredfrom the simulation. In Fig. 5, we show the contribution tothe H I bias from each logarithmic halo mass bin of finitestepsize, Δ 𝑏 HI ( 𝑀 h ) Δ log 𝑀 h = 𝑛 ( 𝑀 h , 𝑧 ) 𝑏 ( 𝑀 h , 𝑧 ) 𝑀 HI ( 𝑀 h , 𝑧 ) ( Δ 𝑀 h / Δ log 𝑀 h ) ∫ 𝑀 max 𝑀 min 𝑛 ( 𝑀 h , 𝑧 ) 𝑀 HI ( 𝑀 h , 𝑧 ) d 𝑀 h . Coincidentally, the linear galaxy bias also typically scales linearly with1 + 𝑧 , because for a passively evolving population, 𝑏 gal ( 𝑧 ) − = [ 𝑏 ( 𝑧 ) − ] 𝐷 ( 𝑧 ) / 𝐷 ( 𝑧 ) (Fry 1996; Skibba et al. 2014), and in a matter-dominateduniverse, 𝐷 ( 𝑧 ) ∝ ( + 𝑧 ) − . However, this cannot explain the nearly linearscaling of H I bias evolution we find herein, because the above relation onlyholds for a tracer with conservative total number, i.e. a passively evolvingpopulation, and therefore the bias is predicted to be either always greateror always smaller than unity. But Figure 4 shows that the linear H I biascrosses the unity between 𝑧 = 𝑧 = I Bias Linearity 7 − − k [ h /Mpc] . . . . . b H I ( k ) Figure 6.
The HI bias defined by auto-power spectrum after correct-ing for shot noise, 𝑏 HI , auto (dots) (see equation 6), and cross-powerspectrum, 𝑏 HI , cross (dot-dashed lines) (see equation 7), at variousredshifts 𝑧 = .
0, 2.0, 1.2, 1.0, 0.5, and 0 (for dots and lines fromtop to bottom, respectively) in the LK model. The upward anddownward arrows mark the scale at which the HI bias 𝑏 HI , cross and 𝑏 HI , auto deviate from the linear bias at the 1.5% level, respectively. (For the 𝑛 th -bin, Δ 𝑀 = 𝑀 𝑛 + − 𝑀 𝑛 , Δ log 𝑀 = log 𝑀 𝑛 + − log 𝑀 𝑛 .) For all three models, Fig. 5 shows that the peakcontribution appears at 𝑀 h = − 𝑀 (cid:12) / ℎ , i.e. theintermediate-mass halos contribute most to the linear H I bias.If we add up the contributions from different halo mass bins,we find that the massive halos of 𝑀 h = − 𝑀 (cid:12) / ℎ only contribute to 3 .
7% of the linear H I bias in TK model,but contribute to about 30% of the linear halo bias. This in-dicates that the decreasing slope of the H I -halo mass relationat the high mass end indeed further suppresses the linear H I bias. Since the H I mass is more suppressed in the massivehalos in the LK and AH model than in the TK model, thisexplains why the linear H I bias is smaller in the former. Wealso point out that since the contribution at our lower masslimit 𝑀 h = 𝑀 (cid:12) / ℎ does not vanish in Fig. 5, especiallyfor the LK and AH model, our results of linear H I bias may beoverestimated due to the neglect of unresolved smaller-masshalos which smooth out the fluctuations.3.3. Auto- vs. cross-power spectrum
In Figure 6, we compare the H I bias obtained from theauto-power spectrum with that obtained from the cross-powerspectrum. The two results agree very well with each otherdown to very small scales. We note that the small differencebetween them does not affect any of the conclusions reachedabove. 3.4. Comparison with previous work
Pénin et al. (2018) performed an analytical calculation thataccounts for the contribution from nonlinear matter fluctu-ations and nonlinear H I modulation, using the combinationof a perturbation theory and a halo model. They employedsix different fitting formulae for the H I -halo mass relation at 𝑧 =
1. Their results show that the H I bias is scale-independentin the range 𝑘 = . − . ℎ Mpc − . However, their re- sults indicate that the H I bias is weakly scale-dependentat 𝑘 ∼ . ℎ Mpc − and significantly scale-dependent at 𝑘 > . ℎ Mpc − (see their figure 4), which is different fromour results. Because of the limitation of our simulation vol-ume, our results are reliable only for 𝑘 > . ℎ Mpc − , whichmakes it difficult to test the presence of scale-dependence onultra-large scales. The redshift 𝑧 = I bias is scale-independentdown to scales smaller than 𝑘 ∼ . ℎ Mpc − . The differenceon small scales may be due to the different methodologiesadopted in the two investigations. While perturbative calcu-lations can provide important insights, numerical simulationscan account for nonlinear effects more accurately.Umeh et al. (2016) and Pénin et al. (2018) investigatednonlinear effects in observations , such as nonlinear redshift-space distortion and nonlinear lensing, based on perturbativecalculations, and found that these nonlinearities can also pro-duce scale-dependent bias on large scales. We will explorethese effects in numerical simulations in the future.Spinelli et al. (2020) investigated the H I content in ha-los using N-body simulations and semi-analytical model ofgalaxy evolution and gas content. This approach is similar toours but the two differ in details. Their simulation box is thesame as ours, both in a comoving volume of 500 ℎ − Mpc oneach side. They found that the H I bias is scale-independent onlarge scales, which is similar to ours, but their results are muchnoisier (see their Figure 12). They also found that the H I biasis enhanced at high redshifts. However, the H I bias in theirresults is roughly scale-independent down to 𝑘 ∼ ℎ Mpc − up to 𝑧 =
2. In contrast, our results show that the small-scalebias evolves from being enhanced at higher redshift to beingsuppressed at 𝑧 <
1, and that there is a sweet-spot redshiftnear 𝑧 =
1. This difference is likely due to different treat-ments in the H I contents of halos, which can lead to differentsuppression of the nonlinear H I modulation.Villaescusa-Navarro et al. (2018) made a comprehensiveanalysis of H I gas distribution based on the IllustrisTNGmegneto-hydrodynamic simulation. They employed the sim-ilar gas model (Krumholz et al. 2008, 2009a,b; Krumholz2013) to ours , to divide the cold hydrogen gas withineach cell into atomic and molecular components. Dueto the limited simulation volume ( ∼ ℎ − Mpc on eachside), Villaescusa-Navarro et al. (2018) cannot test the scale-dependence of the H I bias at the first BAO peak scale. In While Spinelli et al. (2020) published earlier than us, our original preprintwas posted on arXiv eight months earlier. Nevertheless, since the atomic hydrogen is the dominant component in thecold gas, the H I content inside halos is mostly determined by the gas physicsin hydro simulations, and less affected by the gas model during the post-processing. To see this, Diemer et al. 2018 applied different cold gas modelto the same hydro simulation, but found the similar H I gas mass in galaxiesin their Figure 4. Z. Wang et al.comparison, our work applied the average H I -halo mass rela-tion from Villaescusa-Navarro et al. (2018), nevertheless, tohalos resolved from our N-body simulation with large enoughvolume (500 ℎ − Mpc on each side) — i.e. our TK model —to avoid the finite box effect on the bias at the BAO scales.Therefore, we can directly confirm from simulation that theH I bias is linear at the scale corresponding to the first BAOpeak. On the other hand, both Villaescusa-Navarro et al.(2018) and our work find that the H I bias linearity becomesto break down at 𝑘 (cid:38) . ℎ Mpc − generically, the smallestwavenumber presented in Villaescusa-Navarro et al. (2018).Modi et al. (2019) investigated the clustering of H I gas us-ing the Hidden Valley
N-body simulation which has a largecomoving volume of 1 ℎ − Gpc on each side at 𝑧 = −
6, sotheir work is complementary to ours regarding the focusedregime of redshift. The Hidden Valley simulation can re-solve halos down to 10 𝑀 (cid:12) / ℎ , so they can incorporate theH I gas inside smaller-mass halos at high redshifts than ourwork. They adopted fitting formulae of average H I -halo massrelation similar to Eq. (4) herein, in order to assign H I massto halos, and explore both two-point correlation function andpower spectrum. Similar to our results, they also found thatthe H I bias becomes scale-dependent at 𝑘 (cid:38) . ℎ Mpc − .Their results indicate a strong scale-dependence of H I biasat the large 𝑘 at the high redshift, which is consistent withthe behavior in our results for the redshifts higher than thesweet-spot 𝑧 ≈ SUMMARYIn this paper, we use a large N-body simulation to explorethe H I bias for 21 cm intensity mapping experiments at lowredshifts. We adopt three models, LK, TK and AH, represent-ing empirically, numerically, and observationally oriented ap-proaches, respectively, to assign H I mass to dark matter halosand to account for uncertainties in the H I -halo mass relation.We confirm that the H I gas distribution is a linearly bi-ased tracer of the total dark matter density field on the scalescorresponding to the first BAO peak. However, the H I linear-ity assumption breaks down at 𝑘 > . ℎ Mpc − . The exactbreakdown scale is redshift-dependent, because the nonlin-ear effects that modulate the H I gas in halos evolve withtime. This H I nonlinearity, which is caused by the nonlinearhalo clustering and nonlinear H I content modulation, is in-trinsic and not related to the instrumental and observationaleffects. This imposes a challenge to the upcoming 21 cm intensity mapping experiments in their capabilities to extractcosmological information from the broadband shape of the21 cm power spectrum in this 𝑘 -range where a large numberof modes are located. The result is particularly important forforecasting cosmological constraints with upcoming 21 cmintensity mapping experiments. It is, therefore, necessary tobetter model the H I power spectrum beyond the linear regime,e.g. applying the large-scale structure perturbation theory atthe quasi-linear scales. We note, however, that cosmologicalconstraints from the BAO measurement of the 21 cm powerspectrum is not affected by the nonlinear bias.We find the existence of a characteristic redshift above andbelow which the small scale H I bias is enhanced and sup-pressed relative to the linear bias, respectively. For redshiftsclose to this “sweet spot”, the H I bias is linear down to smallscales. For example, for the LK model, the characteristicredshift is (cid:39) .
2, at which the linearity of the bias extendsfrom large scales all the way down to 𝑘 (cid:39) . ℎ Mpc − . How-ever, the exact value of this “sweet spot" redshift dependsboth on the H I -halo mass relation and on nonlinear clusteringof halos. Determining the “sweet-spot" redshift observation-ally can, therefore, also provide valuable information on starformation and clustering of dark matter halos.Finally, we also find that the linear H I bias is an approxi-mately linear function of redshift for 𝑧 ≤
3. This may makecross-checks between different redshifts more powerful forinterpreting observational data.ACKNOWLEDGEMENTSThis work is supported by the National Key R&D Programof China (Grant No.2018YFA0404502, 2018YFA0404503,2017YFB0203302), and the National Natural Science Foun-dation of China (NSFC Grant No.11673014, 11761141012,11821303, 11543006, 11833005, 11828302, 11922305,11733004, 11773049, 11761131004, 11673015, 11421303,11721303, U1531123). YM and JW were also supported inpart by the Chinese National Thousand Youth Talents Pro-gram. JF acknowledges the support by the Youth innovationPromotion Association CAS and Shanghai Committee of Sci-ence and Technology (Grant No.19ZR1466700). HJM wasalso supported in part by the NSF (Grant No. AST-1517528).We are grateful to Xuelei Chen, Kai Hoffmann, Adam Lidz,Matt McQuinn and Francisco Villaescusa-Navarro for use-ful discussions, and the anonymous referee for constructivecomments.REFERENCES