{'fraction_non_alphanumeric': 0.03122879724203693, 'fraction_numerical': 0.019327746502117982, 'mean_word_length': 4.314168778015981, 'pattern_counts': {'":': 0, '<': 17, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The BINGO Project VI: HI Halo Occupation Distribution and Mock Building

BINGO (Baryon Acoustic Oscillations from Integrated Neutral Gas Observations.) is a radio telescope designed to survey from 980 MHz to 1260 MHz, observe the neutral Hydrogen (HI) 21-cm line and detect BAO (Baryon Acoustic Oscillation) signal with Intensity Mapping technique. Here we present our method to generate mock maps of the 21-cm Intensity Mapping signal covering the BINGO frequency range and related test results. (Abridged)


Introduction
The spin-flip of the electrons in neutral hydrogen (Hi) emit or absorb photons in the wavelength of 21-cm and therefore indicate their location in the Universe. As a tracer of the underlying matter density field, 21-cm intensity mapping can reveal the to-mographic Baryon Acoustic Oscillation signal in a cheap and fast fashion (Wyithe et al. 2008). Therefore, several different experiments have been proposed to measure the 21-cm intensity map, such as CHIME (Vanderlinde & Chime Collaboration 2014), Tianlai (Chen 2012 (Santos et al. 2015), etc. By analyzing the large scale structure of the Universe using 21-cm intensity mapping, we can have better understanding about dark matter and dark energy (Kovetz et al. 2017), which are crucial to understand the Universe. More detailed overview of BINGO project is described in the our companion papers I and II (Abdalla et al. 2021a;Wuensche et al. 2021).
At low redshift Hi mostly contained inside galaxies and dark matter halos (Villaescusa-Navarro et al. 2018). Since z ∼ 10, the UV photons from the first generation stars started reionizing the Universe . At z < 5, the Universe has been mostly ionized (Becker et al. 2001;Fan et al. 2006b,a). The Hi can be "protected" in high density regions Wolfe et al. 2005;Zafar et al. 2013), due to sufficient column depth for self-shielding (Pritchard & Loeb 2012). The state-of-art hydrodynamic simulation IllustrisTNG (Villaescusa-Navarro et al. 2018) has shown that more than 95% Hi gas is inside dark matter halo virial radius.
Therefore, it is valid to assume that Hi gas traces the distribution of galaxies and to use the distribution of galaxies and dark matter halos to estimate the distribution of Hi gas with high resolution N-body simulations. Zhang et al. (2020) has used the ELUCID simulation (Wang et al. 2016) together with a semianalytical model (Luo et al. 2016) to construct the Hi mock. They have carefully studied the Redshift Space Distortion (RSD) effect and shot noise for intensity mapping.
By observing the 21-cm emission line of Hi, we can determine the redshift of the emission source. Then, we can estimate the distance of the sources by applying Hubble's law. However, they are not the real distances of the sources. The redshift of the source is mainly contributed from the recessional velocity and the peculiar velocity. The estimation about the position of the source is distorted by the peculiar velocity. Such distortion is known as RSD effect. Since the peculiar velocity of the Hi gas is mainly affected by gravity from the clustering of matter, the RSD effect also contains rich information about the large scale structure. The RSD effect has been used to measure the growth factor of the Universe (Icaza-Lizaola et al. 2019), test General Relativity (Jullo et al. 2019;Anagnostopoulos et al. 2019) and test different cosmological models (Costa et al. 2017;Yang et al. 2019;Cheng et al. 2019). The forecast of constraints from BINGO is described in our companion paper VII (Costa et al. 2021).
The RSD effect can be classified into two parts, Kaiser effect (Kaiser 1987) and Finger-of-God (FoG) effect (Jackson 1972). The moving of matter towards high density regions under gravity leads to the Kaiser effect. The Kaiser effect is dominant at large scales, which squeezes the distribution of galaxies in the line-ofsight direction. On the other hand, the random motion of matter inside high density regions introduces the FoG effect, which is dominant at small scales. It elongates the distribution of galaxies, making them look like fingers pointing at the observer. The RSD effect in 21-cm intensity mapping is well studied in Sarkar & Bharadwaj (2019), and later improved in Zhang et al. (2020).
There are various factors that contribute to the 21-cm intensity map, and the leading terms are the distribution of Hi gas in real space and the RSD effect (Hall et al. 2013). The contributions from gravitational redshift, the integrated Sachs-Wolfe effect, gravitational lensing etc., can be neglected under current consideration. However, more detailed study of these effects might be necessary in the future experiments with higher sensitivity and resolution. Using the information of positions and velocity of halos and galaxies from simulations is sufficient to create a mock 21-cm intensity map for BINGO.
In Sec. 2, we give the overview about the goals and methods of building mock map. We introduce the Hi halo occupation distribution (HOD) model built from ELUCID simulation in Sec. 3. To validate the HOD model and the halo mass resolution, we test our method using ELUCID simulation in Sec. 4. In Sec. 5, we introduce the method of building a full-sky 21-cm intensity mock map. The measurements of the mock are presented in Sec. 6. Finally we summarize the key points we have got in this study and discuss the future works in Sec. 7. This is the sixth of a series of companion papers presenting the BINGO project. The companion paper I is the project paper (Abdalla et al. 2021a), paper II describes the instrument , paper III the optics (Abdalla et al. 2021b), paper IV simulations (Liccardo et al. 2021), paper V data analysis, correlations and component separation (Fornazier et al. 2021), and paper VII forecasts (Costa et al. 2021).

Overview
It is very popular in cosmological surveys to use simulated mocks to mimic real observations before the observations are taken. There are two major purposes for building the mock: 1 Provide a mock data challenge, test the entire data analysis pipeline and see whether we can have the correct constraints of cosmological parameters. 2 Construct the covariance matrix, which is essential in understanding the source of error and uncertainty range of the parameter constraints.
For 21-cm intensity mapping, the data product of the mock is the brightness temperature distribution in 2D spherical surface at different redshifts. Such data product consists of the following four components: 1 signal, the 21-cm emission from extra-galactic sources, which mainly comes from Hi gas in the galaxies and dark matter halos. 2 foreground, coming from radio wave emission in the Milky Way and the extra-galactic sources, which mainly consists of synchrotron radiation and free-free emission. 3 noise, including thermal noise, shot noise and other sources of noise. 4 mask, cover part of the sky for removing too strong foreground, usually around very bright stars or directly toward the Milky Way, or to account for non observed regions.
In this paper, we will only discuss about the 21-cm signal. The other three components are discussed in our companion paper IV and V (Liccardo et al. 2021;Fornazier et al. 2021). In order to generate the mock 21-cm intensity map, the key part is determining the distribution of Hi. There are 6 ways to generate the Hi gas distribution: 1 Hydro: Directly read out the gas particle properties from hydrodynamic simulations, such as IllustrisTNG (Villaescusa-Navarro et al. 2018). 2 SAM: Use the Hi gas mass information from semi-analytical galaxy catalog based on N-body simulations, such as ELU-CID (Zhang et al. 2020). 3 N-body: Use empirical Hi mass-halo mass relation to populate Hi gas in dark matter halos, such as HIR4 (Asorey et al. 2020). 4 Fast Halo: Again, relate Hi mass with halo mass using empirical relations, but the dark matter halo catalog is generated by fast simulation method, such as COLA-HALO (Koda et al. 2016).
5 Model: Populate Hi mass on lognormal dark matter density distribution (Alonso et al. 2014;Xavier et al. 2016) or even Gaussian realization (Asorey et al. 2020). 6 Poisson Halo: Use empirical Hi mass-halo mass relation to populate Hi gas in dark matter halos, but the sub-resolution dark matter halos are generated by Poisson sampling on dark matter density fields measured in N-body simulations (Seehars et al. 2016).
Clearly, Hydro is the most expensive and accurate method. However, Hydro is too expensive to construct a full sky mock. Even the state-of-art hydrodynamic simulation cannot fulfill our requirement. Semi-analytical galaxy catalog is built from merger-tree traced by N-body simulations, which requires subhalo identification. SAM method is less accurate than hydrodynamic simulation, but more efficient. Populating Hi mass in dark matter halos identified from N-body simulation or fast simulations is even more efficient, but less accurate. Using dark matter halos Poisson sampled from dark matter density field can largely reduce the requirement of simulation resolution, which is very useful for large scale application. The lognormal (Gaussian) density generation method is the most efficient among these six methods, but it is the least accurate one. We need to find the balance between accuracy and efficiency for different purpose of mock building. We also need to notice that, in order to construct a full-sky light cone from z = 0.12 to z = 0.45 without using periodic boundary, the box size of such simulation will be larger than 2.4Gpc/h. Currently, we do not have such large scale hydrodynamic simulation. Therefore, in this paper, we will introduce our empirical relations from IllustrisTNG and ELUCID semi-analytical model, and apply that onto Horizon Run 4 simulation (Kim et al. 2015) light cone to generate the mock map. This method can also be applied on fast halo generation catalogs for calculating covariance matrix in the future.

ELUCID SAM catalog
We use the semi-analytical galaxy catalog from ELUCID Nbody simulation (Wang et al. 2016;Luo et al. 2016). The ELU-CID simulation is an N-body simulation with 3072 3 particles in a periodic cubic box of 500h −1 Mpc on a side. WMAP5 cosmology (Dunkley et al. 2009) was assumed in the ELUCID simulation. By following the the merger tree of dark matter halos in the simulations, Luo et al. (2016) constructed the galaxy catalog using their semi-analytical model, which considered the physical process like gas cooling, star formation, supernova feedback, etc. In this work for studying the Hi gas halo occupation distribution, we use the position, velocity, Hi mass and dark matter halo mass of the galaxies from the galaxy catalog, the lower limit of the dark matter halo mass is about 1.85 × 10 9 M /h. We assume that all the Hi mass in the Universe is inside galaxies and their hosting halos, concentrated in the center of the galaxies. As described in Villaescusa-Navarro et al. (2018), from IllustrisTNG simulation (Naiman et al. 2018;Nelson et al. 2018;Pillepich et al. 2018;Springel et al. 2018;Nelson et al. 2019) results, it was shown that at z < 3.0, more than 90% of the Hi gas is inside the galaxies and more than 95% of the Hi mass is inside the halos. Therefore, this assumption is reasonable.
In this galaxy catalog, we also have the information about whether the galaxy is the central galaxy of the halo or a satellite galaxy. This allows us to further separate the galaxies into central galaxies and satellite galaxies for further discussion. Although IllustrisTNG is the state-of-art hydrodynamic simulation, its box size is still quite small (at most 205 Mpc/h for IllustrisTNG-300 simulation and smaller for the other IllustrisTNG simulations). Therefore, a larger box size simulation like ELUCID with high particle resolution can provide better statistics about the Hi halo occupation distribution in very massive clusters. The 21-cm intensity mapping using a semi-analytical model galaxy catalog from MDPL2 simulation has been considered (Cunnington et al. 2020). However the MDPL2 simulation can only resolve halos with mass > 3 × 10 10 M /h (Klypin et al. 2016), one order of magnitude larger than the ELUCID simulation. Without sufficient resolution, it will be hard to tell the Hi distribution in low mass halos like 10 9 M /h. The ELUCID SAM (semi-analytical model) catalog is a good choice for studying the Hi halo occupation distribution in both the high mass end and low mass end. We will use this catalog to study the Hi distribution in halos throughout this section.

Central and Satellite Galaxies
In general, gas cannot efficiently cool and form galaxies in both too massive dark matter halos or too light dark matter halos. However, due to different formation history, central galaxies and satellite galaxies can host very different fraction of Hi gas. Padmanabhan & Refregier (2017) provided a halo model and an Hi mass-halo mass relation, but the subhalos inside main halos was not taken into account separately. An Hi mass-halo mass relation was provided by Guo et al. (2017), which was fitted to ∼ 16000 galaxies in the range of 0.0025 < z < 0.05 from ALFALFA. The importance of satellite galaxies has been illustrated in an updated halo model by Paul et al. (2018). By studying the different Hi gas distribution in satellite galaxies and central galaxies using the SAM galaxy catalog at different redshift, we hope to provide a different Hi mass-halo mass relation, which is more suitable for our redshift range. We show the number of galaxies as a function of their Hi mass and host halo mass in Fig. 1, measured from ELUCID SAM catalog at z = 0. We can see that for central galaxies, it is more abundant with less massive host halo and less Hi mass. It is easy to understand because less massive halos are more abundant, and Hi mass positively correlates to the halo mass if the halo mass is not more massive than ∼ 10 12 M /h.
Together with the number of galaxies represented by different color, we have also shown the mean value and the standard deviation in Fig.1. We can see the mean Hi mass of the central galaxies clearly has a peak around host halo mass M halo = 10 11.6 M /h. Below that, Hi mass positively correlates to the halo mass and beyond that peak, Hi mass drops with larger halo mass until M halo ∼ 10 13 M /h. The mean Hi mass of satellite galaxies depends weakly on the halo mass. With larger halo mass, the number of satellite galaxies rises up quickly, and the standard deviation of their Hi mass also increase. Therefore, the total Hi mass in a dark matter halo consists of the contribution of the central galaxies and satellite galaxies differently for different halo masses.
The Hi mass function measured from the ELUCID SAM catalog at z = 0 is shown in Fig.2. We have shown the mass function of satellite galaxies and central galaxies separately, as well as the total Hi mass function. The Schechter function given in Martin et al. (2010) is shown in red curve for comparison. The SAM Hi mass function can well fit the Schechter function at M HI > 10 8 h −1 M but deviate from it at low mass end. Overall, the SAM catalog still provides a reasonably similar mass function to that fitted from ALFALFA observation.   The Hi mass function measured from the ELUCID SAM catalog at z = 0 is shown. The result of central (satellte) galaxy is shown in blue (orange), and the total result is shown in green. The red curve shows the Schechter function provided by Martin et al. (2010), which represents the fitting function of ALFALFA observation. Overall, the Hi mass function between the SAM and observation is consistent.

Hi Mass Halo Mass relation
For generating mock 21-cm map, the data product from N-body simulations or fast simulations is usually dark matter halo catalog. The most popular and simplest information in the halo catalog are the position, the velocity and the mass of the halo. Therefore, it is most useful if we can find an Hi mass-halo mass relation, providing Hi mass as a function of halo mass. We can see in Fig. 3 that, for low mass halos, the major contribution of Hi mass is from the central galaxies, and for high mass halos, the major contribution comes from satellite galaxies. Although in massive halos, there are more Hi mass in central galaxies than in every individual satellite galaxies on average, the number of satellite galaxies is so huge that they still contribute the most Hi mass in massive halos. It is very important to consider the satellite galaxy contribution of Hi mass in massive halos, which is ignored in HIR4 mock (Asorey et al. 2020). We have also com- In the upper(lower) panel, we show the Hi mass-halo mass relation at z = 0(0.66). The total Hi mass is shown in blue curve, the contribution from central galaxies is shown in orange and the contribution from satellite galaxies is shown in green. No matter z = 0 or 0.66, the central galaxy always contributes the most Hi mass in low mass halos and satellite galaxies contribute the most Hi mass in high mass halos. We have also shown the observed Hi mass-halo mass relation in red dots at z = 0 (Guo et al. 2020) for comparison.
pared the Hi mass-halo mass relation to the observational results (Guo et al. 2020) at z = 0, the results from ELUCID galaxy catalog fits relatively well to the observation. We first fitted the measured Hi mass-halo mass relation from ELUCID SAM catalog using a linear function plus a Gaussian function with free parameter a, b, c, d and f , then we get the parameters as functions of redshift z, by assuming that the parameters are linearly related with redshift in 0 < z < 0.66. Then we find the following analytical formula, which can be easily applied to halo catalogs with the information of halo mass, log 10 (m Hi ) = a log 10 (m) where z is the redshift, m Hi is the Hi mass in unit of M /h, m is the halo mass in unit of M /h. The halo mass definition is based on FoF (Friend-of-Friend) halo finding algorithm, the linking length is 0.2 times the mean particle separation distance. The fitting result is shown in Fig. 4 Notice that the Hi mass-halo mass relation is very similar between z = 0 and z = 0.66, therefore, the assumption that the parameters are linearly related with redshift in 0 < z < 0.66 is reasonable.

The Halo Mass Threshold
When building the mock Hi distribution, the mass resolution and box size of the halo catalog we use is usually limited. We need to know what is the effect of limited mass resolution on the Hi power spectrum. Such effect can be quantified by the bias of Hi gas in the mass limited sample. The Hi gas traces the matter distribution in a biased way, we usually use the bias parameter b Hi to describe the difference between Hi distribution and the matter distribution. Similarly, the dark matter halos also trace matter distribution in a biased way, and the halo bias is b. It has been addressed in HIR4 (Asorey et al. 2020) that the bias of Hi power spectrum can be calculated as the weighted average of halo bias, where b(m, z) represents the halo bias, m Hi is the Hi mass and m is the halo mass. The halo bias can be calculated by empirical formula calibrated by N-body simulations (Tinker et al. 2010). Notice that the halo bias used here is scale independent, which means the resulted Hi bias model is also scale independent. In other words, if the Hi bias is actually scale dependent, our bias model cannot provide a proper description. With the high resolution ELUCID simulation, we can shed light on whether a scale independent bias model is a proper description or not. We selected the galaxies according to their host halo mass with the mass threshold m > 5 × 10 9 M /h, 1 × 10 10 M /h, 1 × 10 11 M /h, 1 × 10 12 M /h. In Fig. 5, we show the Hi power spectrum of the selected samples and the full sample on the left panel, while on the right panel we show the ratio of Hi power spectrum of the selected samples and the full sample. The linear matter power spectrum calculated using Eisenstein & Hu (1998) is shown for reference in cyan color on the left panel.
Considering the RSD effect, the positions of the galaxies (x) are mapped to redshift space (s) using the plane-parallel approximation where v los is the peculiar velocity of the galaxies along the lineof-sight. The redshift space results are shown in dashed lines in Fig. 5. Notice that we have deducted the shot noise power for the illustration and calculation of the power spectrum ratio.
Only when the ratio of power spectrum is close to a constant, we can safely use the scale independent halo bias to calculate the Hi bias. Considering the scale of 1 degree in the sky, which corresponds to k ∼ 0.5hMpc −1 at z ∼ 0.12, the power spectrum ratio of m > 10 11 M /h from k ∼ 0 to k ∼ 0.5hMpc −1 is close to a constant, but in the m > 10 12 M /h case, the ratio cannot be well represented by a constant. Therefore, we conclude that the halo mass resolution for making a mock 21-cm intensity map requires at least ∼ 10 11 M /h. The finer the mass resolution, the better.

Repopulate Hi Mass
We use the Hi mass-halo mass relation introduced in Sec. 3.3 to repopulate the main halos in ELUCID simulation at z = 0 with Hi gas, in order to test the performance of our Hi HOD model. We have generated three kinds of Hi maps and calculated their power spectrum, 1 Use the galaxy (halo) positions in real space, treat the galaxies (halos) as point mass and calculate the Hi distribution in the SAM (HOD) catalog, 2 Use the galaxy (halo) positions in redshift space, treat the galaxies (halos) as point mass and calculate the Hi distribution in the SAM (HOD) catalog, 3 Use the galaxy (halo) positions in redshift space, treat the galaxies (halos) as extended mass due to Hi velocity dispersion inside galaxies (halos) and calculate the Hi distribution in the SAM (HOD) catalog.
We modelled the Hi velocity dispersion inside galaxies at z = 0 as,  Pk HOD /Pk SAM real red red2 Fig. 6. On the left, we show the Hi power spectrum, the solid lines represent the power spectrum measured from the SAM galaxy catalog, the dashed lines represent the power spectrum measured from the HOD model. "real" means real space point mass power spectrum, "red" means redshift space point mass power spectrum and "red2" means redshift space inside-galaxy Hi distribution modelled power spectrum (Zhang et al. 2020). On the right, we show the power spectrum ratio between the results of HOD model and SAM galaxy catalog.
we do not consider the Hi velocity dispersion inside galaxies (halos) when making the mock 21-cm intensity map using N-body simulations. Zhang et al. (2020) has pointed out that with the proper consideration of such Hi velocity dispersion, we can have a more realistic description of Hi mass distribution in redshift space. Therefore, we have made these three kinds of Hi maps for testing. The results are shown in Fig. 6

Horizon Run 4 Catalog
Horizon Run 4 (HR4 in short) is a cosmological N-body simulation with cubic box size of 3150Mpc/h and 6300 3 number of particles (Kim et al. 2015). HR4 simulation adopted WMAP5 ΛCDM cosmology, where Ω m = 0.26, Ω Λ = 0.74, h = 0.72, σ 8 = 0.79, n s = 0.96 (Dunkley et al. 2009). The halos are identified by Friend-of-Friend (FoF) algorithm with linking length equal to 0.2 times the particle mean distance in HR4 simulation, at least 30 particles are required for identifying a halo, which means the minimum halo mass is 2.7 × 10 11 M /h. We use the halo light cone catalog and the mock galaxy light cone catalog to generate our mock 21-cm intensity map. Usually, in order to generate a full sky light cone catalog, we will use the periodic boundary conditions of the simulations to pile up simulation boxes together to get a deep light cone. This technique will introduce some duplicate structures. However, the HR4 simulation box size is so large that there is no periodic structures in z < 0.6 in the full sky light cone catalog.
The redshift range of BINGO is 0.127 < z < 0.449, therefore, there is no duplicate structures in our mock map. In summary, there are four steps in making the mock map, 1 Populate Hi gas in the galaxy (halo) catalog and calculate the Hi density distribution ρ Hi 2 Calculate the brightness temperature where ρ c is the critical density of the Universe (Villaescusa-Navarro et al. 2018), 3 Re-scale the map according to required Ω Hi , 4 Smooth the brightness temperature map with a Gaussian kernel FWHM = 40 arcmin, which is due to the beam of BINGO.

Abundance Matching
Halo Abundance Matching (HAM) and Sub-Halo Abundance Matching (SHAM) have been widely used to link the dark matter halos to galaxies in observations. Under the assumption that more luminous galaxies (with more Hi mass) are populated inside more massive halos (sub-halo), we can populate the halos in the simulations with galaxies, which is adopted in HR4 galaxy catalog generation (Hong et al. 2016). In HR4 galaxy catalog, the position, velocity and an "artificial mass-like quantity", which is also used to estimate the galaxy luminosity with abundance matching, are provided. We use the Schechter form Hi mass function measured from the ALFALFA catalog ( where m Hi is in units of M . We also calculated the abundance of galaxy "mass" in the redshift range 0.127 < z < 0.449 in the HR4 galaxy light cone catalog. By matching the abundance of galaxy "mass" and the Hi mass function, we get a Hi massgalaxy "mass" relation. Therefore, we can populate Hi mass into the galaxy catalog. Notice that the measurement of Hi mass function from ALFALFA is at z = 0, which is not compatible with the light cone catalog. It is just an approximated way of populating Hi mass into galaxies. Then we sliced the catalog into redshift bins of 10 MHz each. We calculate the density of Hi in each redshift using HEALPix pixelization (Gorski et al. 2005) with N side = 256. Finally, using Eq. (5), we calculate the brightness temperature distribution. This set of mock map is labelled as HAM mock map in the following sections. In Fig. 7, we show the mollview of brightness temperature map at 990 MHz -1000 MHz bin, in units of mK, as an example. The map represents the Hi distribution in real space.

HOD Populate
Using abundance matching method has the following shortage, 1 Use the Hi mass function at z = 0 to match all galaxies at different redshift is not fair, the evolution with redshift is ignored and the Hi mass is mismatched, 2 The assumption that galaxy with more "mass" has more Hi mass may not hold, 3 It is not convenient to calculate the bias from theory to compare with the mock map, Therefore, we need to use the HOD approach to obtain the mock map.
Using the FoF halo light cone catalog from HR4 simulation, we can use the Hi mass-halo mass relation described by Eq. (1) to populate Hi mass into the halo catalog. Then we follow the same process described in the last subsection, separate the catalog into 30 bins and pixelize the full sky map using HEALPix with N side = 256. We get the brightness temperature map and show it in Fig. 8. This mollview map is from 990 MHz to 1000 MHz, representing the real space distribution of Hi gas, the same as Fig. 7. Comparing these two maps, we can see that different ways of populating Hi gas into galaxies (halos) leads to different results. Though the map looks very similar at large scale, the brightness temperature fluctuation in the HAM map is much higher than that in the HOD map. This illustrates different bias of Hi distribution. With abundance matching method, Hi gas is more clustered than that in the HOD method. More detailed discussion will be provided in Sec. 6. With intensity mapping, we can tell different ways of Hi clustering and occupation in halos.

Redshift Space Distortion
The intensity map we observe is binned according to the frequency, which is determined by the redshift of the 21-cm emission source. Thus, we are observing in redshift space rather than real space. We need to consider not only the redshift of emission sources due to the Hubble flow, but also consider the redshift contribution from the peculiar velocity of the sources. The positions of the 21-cm emission sources are different in redshift space from that in real space. This is known as RSD effect, which consists of Kaiser effect in large scale (Kaiser 1987) and Fingerof-God (FoG) effect in small scale (Jackson 1972). We considered two possible sources of the RSD effect in the mock map, 1 The peculiar velocity of galaxies (halos) as a point source, which mainly contributes to Kaiser effect, 2 The velocity dispersion inside the halos, which mainly contributes to FoG effect.
In Fig. 9, we show the mollview mock map of the brightness temperature adopting the point source assumption in the upper panel, while the result of additionally considering velocity dispersion is shown in the lower panel. Under the point source assumption it is easy to adopt the RSD effect, since the redshift of the source can be calculated by, where z is the redshift of the point source calculated according to its distance to the observer and hubble flow, which is provided in the light cone catalog, v r is the peculiar velocity along the lineof-sight direction and c is the speed of light. Then we can bin the sources according to their redshift z , calculate the Hi density with N side = 256 HEALPix pixelization and calculate the corresponding brightness temperature. Comparing to the real space mock map, after considering the RSD effect, some of the sources are moved away from its real space frequency range, and some of them are moved into the frequency bin. Thus, we can see slightly different large scale structures between the real space mock map and the redshift space mock map. The mock map considering velocity dispersion of Hi gas in halos is built following the steps, 1 Calculate the velocity dispersion of Hi gas in each halo according to the halo mass m and real redshift z, which is fitted from IllustrisTNG simulation (Villaescusa-Navarro et al. 2018) and adopted in Zhang et al. (2020) (notice that this is a fitting formula using multiple data points at different redshifts, the calculated result at z = 0 is not exactly the same as Eq. 4), 2 Calculate the corresponding comoving distance variance By doing so, we have assumed that the mass distribution of Hi in one halo in redshift space follows Gaussian distribution. This process can provide a rough estimation of FoG effect in the mock map.

Rescale the Map
The average brightness temperatureT b can be calculated from the mock maps at different redshifts, it also represents the Hi abundance Ω Hi since the mean temperature is proportional to it. The relation betweenT b and Ω Hi is In the mock map, since we have neglected the contribution from halos whose mass is smaller than the resolution, we may have missed some of the Hi gas in small halos. On the other hand, we may have also overestimated the Hi mass in halos comparing to the real observation by various reasons, such as overestimating the Hi mass-halo mass relation. Therefore we need to rescale the maps to meet the requirement of Ω Hi . We can rescale the maps according to the required Ω Hi value, by applying T b = T b Ω Hi Ω Hi .

Angular Power Spectra
We would like to address three questions in this subsection, 1 What is the effect of different ways of linking Hi gas with galaxies (halos)? 2 What is the effect of RSD in angular power spectrum? 3 Can we identify the above difference with BINGO telescope?
In Fig. 13, we have shown the zoom-in view of a 10 × 10 degree 2 brightness temperature map of the full sky mock map from 990 MHz to 1000 MHz. We have compared the map generated using HOD method (labeled as HOD), abundance matching method (labeled as HAM) and re-scale the HAM map to the same Ω Hi as the HOD map (labeled as Re-HAM). We have also shown the map smoothed by a Gaussian kernel with FWHM = 40 arcmin.The same color map has been given in each row for a fair comparison. The T b contrast in HAM map is clearly much higher than that in HOD map, which can be more easily identified in the Re-HAM map. In the HAM map, there are more Hi gas inside clusters than in the HOD map. Since the bias of more massive halos is higher than less massive halos, we will expect that with the same Ω Hi , the angular power spectrum of the mock map generated by abundance matching method will be higher than the map generated using HOD.
In Fig. 10, our expectation is verified. Re-HAM and HOD share the same Ω Hi , however the angular power spectrum of the Re-HAM mock map is much higher than the HOD mock map. The HAM mock shows very close angular power spectrum to the HOD mock by coincidence. We can also notice that at large l (l > 200), the real space angular power spectrum of the HAM mock is slighter larger than the HOD mock. This phenomena can also be noticed in Fig. 6. We know that in the HOD mock, the Hi gas are assumed to locate at the center of the halo, without any substructures. While in the mocks that contain subhalos, the subhalo distribution inside main halos leads to additional clustering at small scale, which is known as one-halo term. This can be seen when comparing the HR4 HAM mock and HR4 HOD Cl(mK 2 ) HOD real HAM real Re-HAM real Fig. 10. The comparison of angular power spectrum among HOD mock, HAM mock and Re-HAM mock in real space. With the same Ω Hi , Re-HAM mock has higher angular power spectrum than the HOD mock, which means it has a higher bias. Though very close at large scale, HAM mock has higher angular power spectrum than the HOD mock at small scale due to the existence of substructures in HAM mock. Cl(mK 2 ) HOD real HOD red HOD red2 Fig. 11. The comparison of angular power spectrum between the HOD mock in real space and redshift space. The power in redshift space is higher than that in real space at large scale due to Kaiser effect. No matter we treat the Hi gas in halos as point sources ("red") or extended sources due to internal velocity dispersion ("red2") in redshift space, the difference between them is not significant. mock in Fig. 10. It can also be seen when comparing ELUCID SAM mock and ELUCID HOD mock in Fig. 6. Therefore, it is expected that at small scale, the real space power spectrum of HOD mock has smaller amplitude than that of the HAM mock or SAM mock in real space. However, when considering RSD effect, since the subhalos in halos are also moving randomly in the frame of the center of mass of the main halo, there will be FoG effect that can suppress the power spectrum at small scale. Thus, in the "red" mock, the higher power from one-halo term and the lower power from FoG counter each other, the power spectrum between HOD mock and SAM from ELUCID simulation are getting closer. Furthermore, in the "red2" mock, the FoG effect is taken into account in the HOD mock, the difference between HOD and SAM are revealed again. In Fig. 11, we have shown the comparison between the HOD real space mock (labeled as "real"), the HOD redshift space point source mock (labeled as "red") and HOD redshift FoG effect in-A&A proofs: manuscript no. bingo06_mock Cl(mK 2 ) HOD red HOD red2 HAM red Re-HAM red Fig. 12. We show the angular power spectrum of the HOD mock, HAM mock and Re-HAM mock in redshift space. Re-HAM mock is clearly higher due to a higher bias. The difference between HAM mock and HOD mock due to one-halo term and the difference between "red" mock and "red2" mock due to FoG effect at small scale are not significant because of smoothing introduced by Gaussian kernel FWHM = 40 arcmin.
cluded mock (labeled as "red2"). The angular power spectrum difference between the "real" mock and "red" mock is clear at large scale due to Kaiser effect, but the difference between "red" mock and "red2" mock due to FoG effect at small scale is not significant. This is because of the bin width of the mock being 10 MHz, that corresponds to more than 20 Mpc/h, while the typical velocity dispersion is no larger than 1000 km/s, which roughly corresponds to 10 Mpc/h in redshift space. Therefore, the power spectrum suppression introduced by the velocity dispersion is averaged out by the large bin width. In Fig. 12, we have compared the angular power spectrum of the mocks smoothed by Gaussian kernel FWHM = 40 arcmin, which is the beam size of BINGO. After smoothing, the difference at small scale is easily covered. The difference of bias can still be clearly identified by comparing HOD mock and Re-HAM mock, which share the same Ω Hi value. Therefore, we can answer the questions raised in the beginning of this subsection, 1 Different ways of linking Hi gas with galaxies (halos) lead to different bias and one-halo term, the bias difference can be seen in the total power of angular power spectrum and the difference of one-halo term can be seen in the small scale, 2 Kaiser effect leads to higher angular power spectrum at large scale. FoG effect leads to suppression of power spectrum at small scale, but cannot be discriminated due to 10 MHz bin size, 3 Gaussian kernel smoothing eliminates the difference at small scale, but the bias difference shown in all scales and the Kaiser effect shown in large scales can still be easily identified.

Pixel Histogram
We have shown that the brightness temperature fluctuation in the HAM mock is larger than in the HOD mock in Fig. 13. More quantitatively, we calculated the histogram of the pixels in the unsmoothed mock map. In the upper panel of Fig. 14 On the other hand, the lognormal distribution was widely used to simulate dark matter tracer fields (Xavier et al. 2016), which is the 21-cm intensity map in our case. We calculated the histogram of the pixels from 0 mK to 0.4 mK and fitted with the lognormal distribution using the stats routine in scipy (Virtanen et al. 2019). The result is shown in the lower panel of Fig. 14. We found that the lognormal distribution can reasonably represent the measured results from our mock map. This means that using lognormal distribution to simulate the brightness temperature map is reasonable, but with clear limitations as well.
The difference between the pixel histogram of the HOD mock and the HAM mock shows that the way how the Hi gas is distributed inside halos can largely affect the intensity map, not only in the angular power spectrum, but also in other statistical measurements.

Lognormal Fitting
The lognormal distribution can provide a reasonable fitting to the pixel histogram of the mock map, both in redshift space and in real space. Furthermore, after we considered the 40 arcmin Gaussian kernel smoothing for the mock map, we found that the lognormal function can fit the pixel histogram very well. We illustrated one slice as an example in Fig. 15. We fitted the histogram by a lognormal distribution using the stats routine in scipy (Virtanen et al. 2019). For every redshift bin, we fitted the histogram by lognormal distribution and showed the parameters in Fig. 15. The lognormal probability density function is where the x variable is given by T b in our case. We have used the parabola function, f (z) = a + bz + cz 2 , to fit these three parameters of the lognormal distribution as a function of redshift z. The fitted results are given in Tab. 1.
Given these parameters, one can generate a similar mock map as ours using lognormal distribution (Xavier et al. 2016). It will be helpful for further study and comparison.

Effect of Bin Size
As discussed in Sec. 6.1, different ways of populating Hi gas in galaxies (halos) will lead to different bias and one-halo term, and  the bin width of 10 MHz will largely wipe out the effect of FoG. One step further, it is interesting to ask the following question, if we use different bin width, for example 2 MHz or 5 MHz, will the conclusion be different? Can we see more significant effect of FoG or different small scale power introduced by one-halo term with smaller bin width? To answer this question, we have generated the mock map with HOD and HAM method in 2 MHz bin width and 5 MHz bin width, following the same steps as we used for generating the 10 MHz bin width mock map.
In Fig. 16, we show the comparison of angular power spectrum among the mocks with 10 MHz, 5 MHz and 2 MHz. For a more fair comparison, we take the frequency range from 990 MHz to 1000 MHz, which consists of 1 bin in the 10 MHz mock, 2 bins in the 5 MHz mock and 5 bins in the 2 MHz mock. We take the average of the angular power spectrum measured from Article number, page 11 of 16 A&A proofs: manuscript no. bingo06_mock these 1 (2, 5) bin(s) of 10 (5, 2) MHz mocks. We have also amplified the results to keep the curves of different bin width separated for better illustration. We can clearly see the higher angular power spectrum of HAM mocks compared to the HOD mocks at small scale (l > 100). The angular power spectrum of HOD mocks comparing to HAM mocks are lower by more than 10% at l > 300. It is also obvious that the difference hardly depends on the bin width. It is well expected, since no matter if the difference of small scale power is due to scale dependent bias or one-halo term, they are all isotropic. The width of frequency bin is expected to play little role in the power spectrum. In contrast, the FoG effect is anisotropic. The FoG effect is pointing toward the line-of-sight direction. We expect to see stronger suppression due to FoG effect at small scale with smaller bin width, since there is less smoothing. The result is shown in Fig. 17, which is exactly what we have expected. We have taken point source assumption in the "red" mock, and take the velocity dispersion into account in the "red2" mock. Therefore, the FoG effect is modelled in the "red2" sample. We can see that the FoG effect is very clear in the 2 MHz mock, weaker in the 5 MHz mock and almost indistinguishable in the 10 MHz mock. We can conclude that smaller frequency bin width is very effective in identifying the FoG effect in 21-cm intensity map-ping observations. However, we should also notice that with smaller bin width, the noise level is also higher. Therefore, it is an open question whether we can identify the FoG effect in the real observations by BINGO.

Comparison with UCLCl
If we take the mock maps and the measured angular power spectrum as real observation, we need a fast model prediction to compare with the map and see whether we can fit it to the observation. It is essential for parameter fitting and cosmological parameter constraints. UCLCl (Unified Cosmological Library for C 's) is a package that calculates the angular power spectra (McLeod et al. 2017;Loureiro et al. 2019) from the power spectrum of density fluctuations through where W l (k) is the window function that absorbs all the processes involved in the evolution, convolution and the projection effects, and ∆ 2 (k) is the dimensionless matter power spectrum. By using the linear perturbation calculation of the matter power spectrum, we can get the linear angular power spectrum. By using the halofit nonlinear correction (Takahashi et al. 2012), we can take the one-halo term nonlinear effect into account. We have calculated the angular power spectrum using the UCLCl package with the cosmological parameters consistent with the HOD mock. In Fig. 18, we have compared the results of the UCLCl and our HOD mock measurements in redshift space. The point mass assumption was taken in this HOD mock for comparison. The HOD mock results are higher than the linearly calculated ones, but lower than the halofit corrected results at small scale. At large scale, the UCLCl results are consistent with the HOD mock except for the obvious fluctuations due to cosmic variance. At lower redshift (higher frequency), we can see larger small scale nonlinear correction than at higher redshift (lower frequency). There are two reasons for that, one is that at higher redshift, the one-halo term plays a smaller role than at lower redshift, the other one is that at lower redshift, the shell of Hi distribution that we can measure is closer to us, which means that at lower redshifts we can see more details from small scales with the same angular resolution when compared to higher redshift shells.
If we take the HOD mock as a benchmark, we can see that the linear case have underestimated the angular power spectrum at small scales and the halofit correction has overestimated the angular power spectrum. At higher redshift, the UCLCl calculation is closer to the measurements of the HOD mock. At l < 80, the UCLCl linearly calculated results, which adopted the linear bias assumption, can fit to the HOD mock from 1250 MHz to 1260 MHz. From 960 MHz to 970 MHz, the UCLCl linear calculation can fit the HOD mock result up to l < 200. This comparison clearly tells us the shortage of our modelling for the 21-cm intensity mapping angular power spectrum using UCLCl package. Without further improvement in the modelling at small scale and better understanding about the Hi distribution in halos, it is hard to use the full information of angular power spectrum at l > 80 at the lowest redshift and l > 200 at the highest redshift for BINGO.

Hi Bias and Abundance
The overall amplitude of the angular power spectra of 21-cm intensity mapping is closely related to Ω Hi and b Hi . Therefore, observations can give constraints on these parameters. The detailed   Fig. 15. In the upper left panel, we have shown the histogram of the pixels of the brightness temperature distribution of the redshift space mock map from 990 MHz to 1000 MHz smoothed by a 40 arcmin Gaussian kernel. After smoothing, the histogram can be very well fitted by the lognormal distribution. The fitted parameters, s, loc, and scale, are shown as a function of the redshift (frequency) in the other three panels. For each of them we have also fitted a parabola function, which turns out to be a good fit. The lognormal fitted parameters are shown in red color for the HOD mock and blue color for the HAM mock, crosses for the real space mock and dots for the redshift space mock. The parabola fitting curves are shown in red color for the HOD mock and blue color for the HAM mock, solid lines for the real space mock and dashed lines for the redshift space mock.
forecast is provided in our companion paper VII (Costa et al. 2021). Here we provide the "answer" of the mock map for the future parameter fitting and pipeline building. We calculate the Ω Hi in each mock map by adding up all the Hi mass in each redshift bin, and then divided by the volume of each bin. We show Ω Hi as a function of frequency (redshift) in Fig. 19. The measurement of HOD mock is shown in dotted lines and the HAM mock is shown in dashed lines. The population of Hi mass is different between the real space mock map and the redshift space mock map in the same frequency bin due to RSD effect. We have illustrated the result of real space mock and redshift space mock in different colors, but clearly there is little difference between the real space and redshift space. The RSD has little impact in the Ω Hi in each bin. We have provided a parabola fitting to the real space mocks, which gives (HOD)Ω Hi = 2.7 × 10 −4 + 1.0 × 10 −4 z − 8.0 × 10 −5 z 2 , for the HOD mock and (HAM)Ω Hi = 2.5 × 10 −4 − 4.0 × 10 −5 z − 7.0 × 10 −5 z 2 , for the HAM mock, shown in the black solid lines. For the HOD mock, we use the publicly available python package COLOSSUS (Diemer 2018) to calculate the bias of each halo using the empirical relation provided by Tinker et al. (2010), and by considering the following relation in real space we estimated the bias of the HAM mock. We show the bias b Hi as a function of frequency (redshift) in Fig. 20 the Hi gas distribution, related to the galaxy formation and evolution.
There are three advantages of using the fitting formula: it is useful for reproducing the map with other methods, such as using Gaussian realization or FLASK, which can easily generate large amount of realizations for constructing covariance matrix; it can largely reduce the number of free parameters about Ω Hi and b Hi with fitting formula, when using maps at many different redshifts; this fitting formula is independent of the number of bins, it is useful when we need to change the number of bins.

Summary and Discussion
In this paper, we have introduced the method of building mock 21-cm intensity maps using N-body simulation data. In summary, we have achieved the following points, 1 We have used the ELUCID N-body simulation and its semianalytical galaxy catalog to study the Hi halo occupation distribution. We have concluded that for halos less massive than Cl(mK 2 ) red 10Mhz red 5Mhz red 2Mhz red2 10Mhz red2 5Mhz red2 2Mhz Fig. 17. The comparison of "red" (point source assumption) HOD mock and "red2" (velocity dispersion induced FoG effect modelled) HOD mock is shown with 10 (5, 2) MHz in red (blue,green) color. The "red" result are shown in thick solid lines and the "red2" results are shown in faint solid lines. Notice that for better illustration, we have artificially amplified the curves by 0.5 (1, 2) times for 10 (5, 2) MHz results. The FoG suppression at small scale is more clear with smaller bin width.  10 12 M /h, the major contribution of Hi gas is from the central galaxies. While for halos more massive than 10 13 M /h, the major contribution of Hi gas is from the satellite galaxies. Neglecting the Hi contribution in satellite galaxies is not a good approximation. 2 We have provided a good fitting formula of Hi mass-halo mass relation, make it easier to populate Hi gas into dark matter halos. The fitting function is given by Eq.
(1). 3 The HOD method was verified using the ELUCID simulation by comparing to the SAM galaxy catalog, within 10% accuracy, the power spectrum of the HOD repopulated Hi distribution is correct up to k = 0.4hMpc −1 . 4 The mass resolution of the N-body simulation can lead to scale dependent bias. With the halo mass limit down to  Fig. 19. We show the Ω Hi as a function of frequency (redshift) in our mock map, the dotted lines show the result of the HOD mock and the dashed curves show the result of HAM mock. "real" means the real space mock, "red" means the redshift space mock under point source assumption. The black solid lines show the parabola fitting results to the real space mock. around 10 11 M /h, the linear bias assumption, which assume a constant bias number, is still valid. The linear bias assumption is not valid with 10 12 M /h halo mass limit. This sets a resolution requirement for N-body simulation to be used to build Hi mock map. 5 Using HR4 simulation, we have built the full-sky 21-cm intensity mock maps covering from 960 MHz to 1260 MHz, covering the range of BINGO (from 980 MHz to 1260 Mhz). Both the HOD method and abundance matching (HAM) method are applied to get two kinds of mocks. 6 The redshift space distortion (RSD) was taken into consideration, we introduced the method to consider both the point source redshift space distortion and the method to consider the velocity dispersion to account for finger-of-god (FoG) effect. 7 The lognormal distribution can provide a reasonable fitting to the pixel histogram of the mock map, if the map is smoothed by 40 arcmin Gaussian kernel, the fitting will be very good.
8 We have provided the lognormal fitting parameters for our mock maps. A further parabola fitting to the lognormal parameters as functions of redshift was done. The fitting results are given. It is easy to realize a similar brightness temperature distribution with lognormal distribution using our parameters. 9 We have discussed about the effect of Gaussian beam smoothing, RSD and the difference of HOD mock and HAM mock by comparing the angular power spectrum. Different bias and Kaiser effect can be discriminated. The nonlinear effect at small scale and FoG effect will not be important if the 40 arcmin Gaussian smoothing is taken into account. 10 The effect of different frequency bin sizes was discussed.
FoG effect can be more easily identified with smaller bin size. Different bin size has little impact in telling the difference between HOD mock and HAM mock. 11 The theoretical calculation was done using UCLCl package and compared to our mock map measurements. We found that the UCLCl results can fit to the mock at large scale, but failed at small scale. The lower the redshift, the worse. They fit well at l < 80 for the highest frequency bin and at l < 200 for the lowest frequency bin of BINGO. The linear theory calculation has underestimated the angular power spectra at small scale and the halofit nonlinear correction has overestimated it. We need a better modelling to make full use of the whole range of angular power spectra. 12 We have provided the measurement of Ω Hi and b Hi as a function of frequency (redshift) in the HOD mock and HAM mock. The parabola function can provide very good fitting to the measurements. The fitted functions in each case are listed in Sec. 6.6.
Beside constraining cosmological parameters, 21-cm intensity mapping is expected to constrain the astrophysical processes related to Hi gas. There is no conclusion about the best method, HAM and HOD, two different methods applied in this paper for populating Hi gas in haloes, can be regarded as valid tests. Since we can distinguish HAM from HOD using 21-cm intensity mapping by constraining the bias parameter, it proves the ability of discriminating the Hi population in haloes. The mock we have generated here will serve as the source signal to test the BINGO data analysis pipeline. This will include parameter fitting, BAO signal fitting and many other methods developed for intensity mapping.
The methodology we have introduced here will be useful in building mocks based on fast simulations such as COLA-HALO (Koda et al. 2016). After generating hundreds of mocks, we can measure the covariance matrix, which is essential for estimating the cosmological parameters. On the other hand, this methodology can also be applied for other 21-cm intensity mapping experiments such as SKA. Next step, we will combine the foreground, noise and mask of BINGO to generate a more realistic mock map for testing our data analysis pipeline. In the future, we will also include more effects that might change the map, such as point radio sources, lensing effects, gravitational redshift etc. By making more and more realistic mock maps, we can improve the understanding about intensity mapping and provide better data challenge for the pipeline.