Homogeneous vs Biased IGM: Impact on Reionization

This paper considers the impact of large scale biasing of the IGM on reionization. The two simplest but extreme scenarios for IGM biasing are: an unbiased IGM which has a constant density and an IGM with density equal to the collapsed matter density. In this work, the relationship between the IGM density and the collapsed matter density is defined through an IGM bias parameter. The two extreme scenarios of homogeneous and perfectly biased IGM are produced for two extreme values of this bias parameter. It is found that, for the same level of reionization (i.e., for same global neutral hydrogen fraction). one could get very different 21 cm brightness temperature distributions for different values of this bias parameter. These distributions could give an order of magnitude more or less power as compared to the uniform case. It is also found that there exists a critical value for the IGM bias parameter for which there could be a near washout of the structure in the 21 cm brightness temperature distribution (i.e., zero power or a nearly uniform 21 cm brightness temperature distribution). To address the problem, a new method of generating 21 cm brightness temperature maps is used. The method uses the results of n-body simulations and then employs ray tracing to obtain the 21 cm brightness temperature maps. Towards the end, a prescription for the IGM bias parameter is given. This is derived within the framework of the Press-Schechter theory.


INTRODUCTION
Epoch of Reionization (EoR) is considered as one of the most difficult challenges in modern cosmology (see e.g., Shapiro and Giroux, 1987;Gnedin and Ostriker, 1997; Barkana and Loeb, 2001, Barkana and Loeb, 2005;Furlanetto et al., 2009;Zaroubi, 2013;Cooray et al., 2019). Analysis of the EoR is demanding on observational (see e.g., Tozzi et al., 2000;Bernardi et al., 2009;Harker et al., 2010;Parsons et al., 2010;Hazelton et al., 2013;Thyagarajan et al., 2015;Dillon and Parsons, 2016;Thyagarajan et al., 2016;Mertens et al., 2018;Barry et al., 2019;Greig et al., 2020) as well as theoretical (see e.g., Madau et al., 1997;Ciardi et al., 2000;Benson et al., 2006;Furlanetto et al., 2006;Mesinger et al., 2011;Giri et al., 2018;Qin et al., 2018;La Plante and Ntampaka, 2019;Park et al., 2019;Liu and Shaw, 2020) fronts. The possibility of observing EoR in the radio window using the 21 cm radiation was suggested a long time ago (Wouthuysen, 1952;Field, 1959b;Sunyaev and Zeldovich, 1975). The 21 cm line emitted by the neutral hydrogen in the IGM gets redshifted to much larger wavelengths and a continuous distribution of neutral hydrogen gives rise to a continuum of emission in the radio window (for reviews, see Fan et al., 2006;Furlanetto et al., 2006;Choudhury, 2009;Pritchard and Loeb, 2012). This paper explores the impact of large scale biasing of the IGM on the 21 cm power spectrum from EoR. Many earlier works have considered a clumpy IGM during the EoR (see e.g., Kohler et al., 2007;Pawlik et al., 2009;McQuinn et al., 2011;Raičević and Theuns, 2011;Finlator et al., 2012;Emberson et al., 2013;Jeeson-Daniel et al., 2014;So et al., 2014;Sobacchi and Mesinger, 2014;Trombetti and Burigana, 2014;Mao et al., 2020). Such an IGM is granular or clustered on small scales and is considered on subgrid level. Its eventual impact on the EoR power spectrum is considered by some authors (see e.g., Mao et al., 2020). This paper addresses large scale clumping or large scale biasing of the IGM. By large scale biasing one means how the distribution of the IGM is in the whole simulation box. Many numerical approaches of simulating EoR assume IGM to be having uniform density far from halo boundaries (see e.g., Mellema et al., 2006;Thomas and Zaroubi, 2008;Thomas et al., 2009;Ghara et al., 2015). While the early versions of semi-analytic methods determine the neutral hydrogen fraction, x HI , based on excursion sets arguments and then assume that the hydrogen density follows the dark matter density (see e.g., Choudhury et al., 2009;Mesinger et al., 2011). Now, the hydrogen in the IGM could have a more complicated dependence on the local density of the massive haloes. So the regions of the universe that have more or less massive halos also have more or less density of the background IGM. The simplest way to consider this is through following expression: In this work, α is termed as the IGM bias parameter and it notifies the biasing of the IGM with respect to the collapsed massive haloes 1 . For α 0, the IGM is uniform which means that the uncollapsed gas is distributed uniformly in the Inter-Galactic regions. For α 1, the IGM is biased in the way that its density is proportional to the collapsed matter density: ρ IGM ∝ (1 + δ collapsed massive haloes ). These are the two extremes the models suggested by Eq. 1 take. Note that as 〈1 + α · δ〉 1, introduction of alpha does not change the value of mean or global amount of neutral hydrogen. It only changes its distribution. It is seen in this work that the two extreme distributions, given by α 0 and α 1, could have very different power spectra. As 21 cm power spectrum is one the prime observables of the EoR, this work suggests that the distribution of the IGM with respect to the collapsed matter should be studied in more depth in order to gain a proper understanding of the epoch.
The paper is organized as follows: The next section discusses the method of generating the 21 cm maps. The method is ideally suited to address the problem of impact of a biased IGM on EoR. The later section discusses the main results of the paper. A prescription for the IGM bias parameter within the Press-Schechter framework is then derived. This is followed by the conclusion of the work in the last section.

METHOD FOR GENERATING 21 CM MAPS
The method used in this work (Raut, 2019) borrows ideas of semianalytic methods (Furlanetto et al., 2004;Choudhury et al., 2009;Mesinger et al., 2011) and extensively numerical methods (Razoumov and Scott, 1999;Abel et al., 1999;Nakamoto et al., 2001). The parameters that are borrowed from semi-analytic approach are: M min (minimum mass of collapsed haloes), n ion (number of ionizing photons deposited in the IGM per collapsed baryon in stars) and n rec (number representing the ratio of ionization rate to recombination rate). This method has similarities with the hybrid approaches of Trac and Cen (2007); Lee et al. (2008); Zheng et al. (2011), but it is computationally less demanding as it emphasizes upon the scales of interests (0.1Mpc −1 (k(1.0Mpc −1 ). Upcoming radio telescopes like HERA (DeBoer et al., 2017) 2 and SKA1-Low (Koopmans et al., 2015) 3 will be more sensitive to these modes (Kulkarni et al., 2016). The fiducial values of reionization parameters used in this work are: M min 10 8 M ⊙ , n rec 1 and x HI {0.4, 0.6, 0.75} for z {7, 8, 9} (denotes the mean or global value of neutral hydrogen fraction.) The details of the method used for generating 21 cm maps is x HI as follows: • The first step is to do N-body simulations and compute overdensity for each grid cell. I have used Illustris-1 simulations (Genel et al., 2014) for this step. The cosmological parameters used are wmap-9 (Hinshaw et al., 2013): Ω M 0.2726, Ω b 0.0456, σ 8 0.821 and h 0.704. The box-size is 106.5Mpc and gridsize is 64 3 . This paper uses a method that is a small modification to the one used in the earlier work. Instead of using excursions sets (Press and Schechter, 1974;Bardeen et al., 1986;Bond et al., 1991;Lacey and Cole, 1993)to find out the collapsed fraction f coll , one uses the total halo distribution itself to find the extent of collapsed structure. The M min provided by above n-body simulations is about 10 8 M ⊙ and so this approach is justified. Second step assigns ionizing photons to each grid cell based on following prescription (Furlanetto et al., 2004).
• Here, f esc is the escape fraction of the ionizing photons, f p is the star formation efficiency and N c/b is the number of ionizing photons produced per baryon in stars. n c is the number of ionizing photons produced in a grid cell which has collapse fraction f coll .
• The next step is to split all the photons, N c cells n c , into a large number ( ∼ 10 7 ) of photon-rays. Each photon ray carries a weight indicating the amount of photons it carries. Initially most of the photon rays have equal weights. Each photon ray is assigned a random position and direction inside the grid cells. Then each photon ray is propagated in small steps, i.e. 0.05 of grid cell size. Note that the number of photon-rays originating in a grid cell is proportional to ∼ 10 7 count of that cell.
• The next step is to implement reionization. At each step of propagation of photon-rays, one computes the total weights of the photon-rays arriving in the grid cells. This total is compared with the amount of neutral hydrogen present in the grid cells. The photon-ray weights and neutral hydrogen counts are decremented at every step accounting for reionization. It is assumed at this step that each photon-ray can ionize amount of neutral hydrogen equal to the photon-ray weight (i.e., 100% ionization efficiency or a very large cross section for ionization reaction).
• The next step implements recombination. After the photonray weights and neutral hydrogen content of each grid cell is revised, the number of neutral hydrogen atoms that were ionized in that particular step, n 1 , is computed. This number is computed for each grid cell separately. Before implementation of the next step of photon-ray propagation, r · n 1 number of neutral hydrogen atoms are added back to each grid cell to account for recombinations. Here r n rec /(1 + n rec ) and is the ratio of recombination rate to ionization rate.
• The final step is to terminate the photon-ray propagation and assign brightness temperature to each grid cell. The first part is implemented when the total of photon-ray weights falls below 0.1% of its initial value. For the second part, we use the wellknown expression (Wouthuysen, 1952;Field, 1959a;Madau et al., 1997;Furlanetto et al., 2006),

THE BIASED IGM ANALYSIS
This section discusses the effect of biased IGM on 21 cm brightness temperature distribution. The biasing of the IGM is modelled through Eq. 1 in the next two subsections.

Modification of the 21 cm Brightness Temperature Distribution With the IGM Bias Parameter
The 21 cm maps are generated for three values of α: 0, 0.5, and 1. As suggested earlier, α 0 corresponds to an IGM that is distributed uniformly in the whole simulation box, while, α 1 corresponds to an IGM that is distributed proportional to the extent of the collapsed structures. This way of assigning IGM to each grid cell could be over simplified and actual distribution of the IGM could be more complicatedly biased with respect to the collapsed matter (see the next section).
Here, value of one is used for the recombination parameter, n rec (Sobacchi and Mesinger, 2014) and value of 10 8 M ⊙ is used for the minimum halo mass M min (Barkana and Loeb, 2001). The n ion values were chosen in such a way so as to produce global neutral hydrogen fraction of 0.4, 0.6, and 0.75 for redshifts 7, 8, and 9 respectively. This set of values of neutral hydrogen fraction is consistent with the CMBR and LyAlpha data (Kulkarni et al., 2016;Raut and Choudhury, 2018). To get the n ion values for each redshift, we solve the following equation (Y He is the Helium fraction): In this work, the global neutral hydrogen fraction, denoted by x HI , is considered as a parameter instead of n ion . Once n rec and 〈f coll 〉 values are fixed, above equation uniquely determines the n ion value for a given x HI .
The reionization data cubes are generated as discussed in the previous section. The 21 cm brightness temperature slices (midway through the simulation box) are shown in Figure 1. This is done for three values of the IGM bias parameter (α 0, 0.5, 1) and for three values of redshift (z 7, 8, 9).
As seen from Figure 1, the topology of the 21 cm brightness temperature distribution is changed as the IGM bias parameter is changed. For α 0, it is outside-in with the 21 cm bright regions are surrounding the 21 cm dark regions. While for α 1, it is inside-out with the 21 cm bright regions are embedded inside the 21 cm dark regions. This is expected as gas in the voids will remain neutral and hence 21 cm bright for the uniform IGM or α 0 case (The hydrogen near collapsed structure will be ionized while the hydrogen away from collapsed structures will remain neutral.) On the other hand, the voids will be devoid of gas and hence 21 cm dark for α 1 case. As expected, for the intermediate value of α 0.5, the distribution is also intermediate and is close to being uniform.
These conclusions are fortified if one analyzes the 21 cm brightness temperature probability distributions. As seen from Figure 2, the probability distributions are opposite for α 0 and α 1 case. For α 0, there are more bright regions as compared to the dark ones and for α 1, there are more dark regions are compared to the bright ones. This is because, voids, which occupy a larger volume, are brighter for α 0 case and darker for α 1 case. For the intermediate case of α 0.5, the plots indicate that the distributions are more or less uniform as they have a single peak and most of the areas in the probability distribution curves are covered by a small range of δT b values.
One can also study the impact of the variation of the IGM bias parameter on the power spectrum. To obtain the power spectrum, the standard procedure is followed. First, the 21 cm brightness temperature distribution (δT b ≡ δ 21 ) is Fourier transformed and dimensionless power spectrum is obtained using following equations: The spherically averaged dimensionless power spectrum is the obtained by averaging Δ 2 (k) over all possible angles Here μ k /k is the cosine of the angle that wavevector k makes with the LOS or parallel direction and ϕ is the azimuthal angle.
As seen from Figure 3, for z 7, the power spectrum is nearly same for the two extreme scenarios of α 0 and α 1. One should take note that although the power spectrum is almost the same, the topologies are opposite. For the other two redshifts, though, there is significant enhancement of power for the case of α 1. This could be because of combination of larger neutral hydrogen fraction and more isolation of collapsed structures as compared to the z 7 case. These two, combined, are producing more contrast in

Raut
Biased IGM and EoR δT b distribution for the α 1 case as compared to the α 0 case. There is also a notable fall in the power spectrum for the α 0.5 case for redshifts z 7 and z 8. As mentioned before, this is because the intermediate value of the IGM bias parameter drives the δT b distribution toward uniformity. This means that distributions of the remnant neutral hydrogen and the collapsed structures are almost perfectly inversely related to each other. Thus the product x HI · (1 + δ b ) that appears in the δT b expression is almost constant. This is giving a nearly constant brightness temperature distribution and very small power spectrum. This behaviour is generic and seems to appear for every redshift for some α (see the discussion in the next subsection).

The Critical Value of the IGM Bias Parameter and the Uniform Distribution
The severe degradation of the 21 cm brightness temperature power spectrum is one of the main results of this paper and hence some further analysis is justified. Figure 4 shows the dependence of the power for three typical modes for three redshifts as a function of α. As mentioned before, these three modes are representative of the scales that would be measured by upcoming and future experiments. As seen from the figure, behavior of all the modes is similar: the plots are U-shaped. The power goes to almost zero for each mode for some value of α. Notably, this value of α seems to depend on the redshift alone and not on the scales. This means that all the scales hit minima at the same value of α. This observation indicates that there is indeed a near washout of δT b fluctuations in the entire simulation box.
A critical value of the IGM bias parameter that causes a near uniform brightness temperature distribution is possibly a characteristic of the reionization process. To check this all three reionization parameters are varied from their fiducial values. This is done for the z 7 case. As mentioned earlier, the fiducial parameter values are M min 10 8 M ⊙ , n rec 1 and x HI 0.4. Figure 5 shows the variation of the three k-modes as a function of the IGM bias parameter for six different combinations of reionization parameters. The central region of the plots indicate the changed parameter while the other two  parameters are kept unchanged from their fiducial values. For all the cases, one can see that there is either a critical value of the IGM bias parameter or a range of values of the IGM bias parameter for which there is a low signal. As the signal is nearly zero, this work warrants a careful study of distribution of IGM with respect to the collapsed structures. It should also be noted that the topology of the 21 cm brightness temperature distribution would be opposite on two sides of the critical value of the IGM bias parameter.

PRESCRIPTION FOR IGM DISTRIBUTION BASED ON PRESS-SCHECHTER THEORY
The IGM bias parameter being a constant throughout the box is ideally not possible but it is the next step to the other two extreme scenarios of it being equal to one or zero. These two extreme scenarios form basis of the EoR simulations done by many earlier works (see e.g., Mellema et al., 2006;Thomas and Zaroubi, 2008;Choudhury et al., 2009;Thomas et al., 2009;Mesinger et al., 2011;Ghara et al., 2015). A possible next step to a constant value for the IGM bias parameter is an IGM bias parameter that is a function of the overdensity parameter. This could be predicted within the framework of the Press-Schecheter theory. The following assumptions are employed to obtain the IGM distribution within this framework: • The results of the Press-Schechter theory are exact, • Baryons follow dark matter, • All the baryons that are not part of galaxies (i.e., massive collapsed haloes) are considered as being part of the IGM. 4  Let, δ IGM ≡ δ lmh ≡ δ low mass haloes (7) δ cmh ≡ δ collapsed massive haloes It is shown in the Supplementary Appendix that α δ IGM /δ cmh for a grid cell with overdensity δ and collapse fraction f coll is given by: Here f coll is average or global value of the collapse fraction. Figure 6 shows the dependence of the IGM bias parameter obtained using above equation.
This prescription should be able to give a good zeroth level baryon distribution. Environmental effects, baryonic physics and other astrophysical effects would give rise to modulations to this distribution. Using this prescription for the IGM density, one can now obtain the 21 cm brightness temperature distribution and it is shown in Figure 7.
The under dense grid cells are more in number than the over dense grid cells and they dominate the overall behaviour. For z 9, the average value of the IGM bias parameter is about 0.22 and is close to the critical value. As a result, the distribution seen is also nearly uniform. One can find the power spectrum of the 21 cm brightness temperature distribution for the Press-Schechter prescription and compare it to the two extreme scenarios of the uniform and the perfectly biased IGM. The results are shown in Figure 8. As seen from the figure, there is degradation of power for all three cases compared to the other two extreme scenarios. It can also be seen that for z 9 there is a severe loss of power for the large scales as the average value of the IGM bias parameter is close to the critical value.

DISCUSSION AND CONCLUSION
In this work we discussed how an IGM that is having varied bias with respect to the massive halo population could impact the 21 cm brightness temperature distribution. It is found that the distribution as well as the power spectrum are significantly affected with the difference in the IGM bias parameter α (Eq. 1), while, the global neutral hydrogen fraction remains constant. It is also found that for some critical value of the IGM bias parameter, the distribution becomes almost uniform causing a near washout of the 21 cm power spectrum. This critical value of the IGM bias parameter is different for different redshifts. Thus, it is essential to estimate the large scale distribution of the IGM. In this work, we also saw that this distribution could be estimated within the framework of the Press-Schechter theory. It was also observed that for z 9, the framework yielded a nearly uniform distribution. The degradation in the structure of 21 cm brightness temperature distribution occurs because there are more ionizing photons in the regions where there is more neutral hydrogen. This situation can also occur for standard scenarios (i.e. α 0 or α 1) if one has an n ion that is density dependent. In this work n ion has been assumed to be globally constant. But, consider case of α 1 for which neutral hydrogen follows the collapsed matter. If one has larger n ion for over dense regions and smaller n ion for under dense regions then one can arrive at a situation in which the extent of neutral hydrogen in each grid cell is about the same. And so the 21 cm brightness temperature distribution is nearly uniform.
Many early large scale structure studies assume that the baryons follow dark matter or ρ b ∝ ρ dm . It is found in many works that, based on theory (see e.g. Cen and Ostriker, 1999;Springel et al., 2001;He et al., 2005;Zentner et al., 2005;Crain et al., 2007;Diemand et al., 2007;Giocoli et al., 2008;Moster et al., 2018) and observations (see e.g. Fukugita et al., 1998;Ettori, 2003;Dai et al., 2010;McGaugh et al., 2010;Eckert et al., 2016;Chan, 2019), the above assumption need not be valid and more work is needed to get a conclusion on the exact dependence. As suggested by these works, the problem of finding exact dependence of the hydrogen content on the massive collapsed structures is a difficult one. In the realistic scenarios, where baryons do not strictly follow the dark matter, there could be more structure in the 21 cm brightness temperature distribution. This would possibly give more power on small scales. Baryonic physics could alter the bias parameter but the average value over the entire box is likely to remain close to the one obtained using the assumption of baryons following dark matter.
There are possibly two future directions to the problem posed in this work: One could assume that the baryon fraction parameter, f b Ω b /Ω M , to be dependent on halo mass and find the its correlation to the halo mass based on observations (see e.g., Padmanabhan et al., 2016) or simulations (see e.g., Genel et al., 2014;Crain et al., 2015;Pillepich et al., 2018;Diemer et al., 2019). One could also take a more direct approach and use the high resolution simulations themselves to use the hydrogen content of various cosmic structures (Martizzi et al., 2019). A very high resolution simulation could directly yield the hydrogen density for the whole simulation box.
The author would like to thank the Illustris collaboration for their N-body simulations. The author would also like to thank anonymous referees whose comments helped in improving the quality of the manuscript.
Both, δ cmh and δ lmh , can be obtained as function of true overdensity δ and hence the IGM bias parameter α can be obtained as a function of overdensity of collapsed massive haloes δ cmh .

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
The author confirms being the sole contributor of this work and has approved it for publication.