Modelling the 21 cm Signal From the Epoch of Reionization and Cosmic Dawn
Abstract
Studying the cosmic dawn and the epoch of reionization through the redshifted 21 cm line are among the major science goals of the SKA1. Their significance lies in the fact that they are closely related to the very first stars in the universe. Interpreting the upcoming data would require detailed modelling of the relevant physical processes. In this article, we focus on the theoretical models of reionization that have been worked out by various groups working in India with the upcoming SKA in mind. These models include purely analytical and seminumerical calculations as well as fully numerical radiative transfer simulations. The predictions of the 21 cm signal from these models would be useful in constraining the properties of the early galaxies using the SKA data.
Keywords: intergalactic medium – cosmology: theory – dark ages, reionization, first stars – diffuse radiation – largescale structure of Universe – methods: numerical –methods: statistical
1 Introduction
One of the major science goals of the SKA is to study the redshifted 21 cm signal of neutral hydrogen from the cosmic dawn and the epoch of reionization (Mellema et al., 2013; Koopmans et al., 2015; Carilli, 2015). The cosmic dawn refers to a period when the first stars formed in the universe, while the epoch of reionization is when the HI in the intergalactic medium was being ionized by the UV radiation from the first stars (see, e.g., Barkana and Loeb, 2001). It is believed that this is process which extended over redshift ranges (Mitra et al., 2015). The redshifted 21 cm radiation is possibly the most promising method of detecting the distribution of the HI at these redshifts (Furlanetto et al., 2006; Pritchard and Loeb, 2012). The signal will contain information about how this process occurred, and will indirectly help in studying the properties of the first stars and the surrounding medium at those early epochs.
Detailed models of reionization are a crucial ingredient for interpreting the data and constraining the EoR. The difficulty in constructing models at such high redshifts is that there are practically no observational constraints on the physics of galaxy formation. As a result it is often not possible to have a good idea about the number of ionizing photons available. The models thus assume simple prescriptions to assign ionizing luminosities to dark matter haloes, and try to constrain the free parameters using existing observations.
This article is meant to highlight the advances made in modelling the EoR and cosmic dawn by members of the Indian community keeping the upcoming SKA in mind. It turns out that groups working in related areas have contributed to different types of modelling, starting from purely analytical ones to complex radiative transfer simulations. The main features of these models and their main results are summarised in the following sections.
2 Present constraints on reionization history
The main components of building a reionization model are as follows:

The abundance (and locations) of dark matter haloes form the first step as the sources of ionizing photons (galaxies) form within these haloes. In analytical models, one can use the standard forms of the halo mass function given by, e.g., Press and Schechter (1974); Bond et al. (1991); Sheth and Tormen (1999), while in the body simulations, the masses and locations of these haloes are usually obtained by group finder algorithms, e.g., the Friendsoffriend (Davis et al., 1985).

Given the dark matter haloes, one needs to work out the physical processes related to galaxy formation, stellar radiation and escape of ionizing photons. All these are highly uncertain at , hence most reionization models tend to assume some prescription which relates the number of ionizing photons to the halo mass. The simplest of these assume the number of ionizing photons in the IGM produced by a halo of mass to be given by
(1) where is an unknown parameter and all other symbols have their usual meanings. Physically is a combination of starforming efficiency, the number of ionizing photons produced by stars and the fraction of ionizing photons that escape from the host galaxy into the IGM. The above relation can also be written in terms of the number of photons in the IGM per unit time per unit comoving volume
(2) where is the comoving number density of hydrogen and is the mass fraction in collapsed haloes (called the collapsed fraction) that are forming stars. The collapsed fraction is related to the mass function by
(3) where is the minimum mass of haloes that can form stars. The value of is decided by the cooling efficiency of the gas in collapsed haloes, e.g., if the gas contains only atomic hydrogen, it is unable to cool at virial temperatures lower than K while the presence of molecules can push the limit to K.

The third component is a description of the inhomogeneous baryonic density field. In analytical models, this could be described by the PDF of the baryonic overdensity (smoothed over the Jeans scale). Some standard forms that are used are the lognormal distribution (Bi and Davidsen, 1997; Choudhury et al., 2001a, b; Choudhury and Ferrara, 2005), or a fitting form motivated by the hydrodynamical simulations (MiraldaEscudé et al., 2000; Bolton and Becker, 2009). While simulating the baryonic densities, it is important to be able to resolve reasonably small scales so as to identify the dense optically thick systems. These regions, because of their high recombination rate, act as “sinks” of ionizing photons and can alter the distribution of ionized regions.

Finally, given the ionization sources and the clumpy baryonic field, one has to solve the transfer of radiation through the IGM accounting for all relevant physical processes. The cosmological radiative transfer equation
(4) where is the monochromatic specific intensity of the radiation field, is a unit vector along the direction of propagation of the radaiation, is the absorption coefficient and is the emissivity, has to be solved at every point in the sevendimensional space. However, the high dimensionality of the problem makes the solution of the complete radiative transfer equation well beyond our capabilities, particularly since we do not have any obvious symmetries in the problem and often need high spatial and angular resolution in cosmological simulations. Hence, the approach to the problem has been to use different numerical schemes and approximations (Iliev et al., 2006a).
In analytical models, one can simplify the problem by taking the global average of equation (4) under the assumption that the mean free path of photons is significantly smaller than the horizon size (Choudhury, 2009). The equation can be written in terms of the volume filling factor of ionized regions as
(5) where is the clumping factor and is the recombination coefficient at a temperature .
The uncertainties in the reionization models can, in principle, be constrained by comparing with existing data, in particular, the quasar absorption spectra blueward of the Ly emission line and the measurements of the Thomson scattering optical depth from the CMB observations. One such physically motivated semianalytical model (Mitra et al., 2011, 2012, 2015) which treats the as an unknown function of and constraints it using a Principal Component Analysis, predicts the evolution of and the average neutral fraction as shown in Figure 1. One can see that the present data sets imply that reionization begins at and is completed close to . There is still considerable uncertainty around the final stages of reionization as is shown by the width of the shaded regions around , and one expects the 21 cm experiments to play an important role in reducing the uncertainties.
3 Modelling the 21 cm signal from EoR
The basic principle which is central to the 21 cm experiments is the neutral hydrogen hyperfine transition line at a rest wavelength of 21 cm. This line, when redshifted, is observable at low frequencies ( MHz for ) against the CMB. The strength of the signal is quantified by the differential brightness temperature given by
(6) 
where is the spin temperature of the gas and
(7) 
The quantities and are the neutral hydrogen fraction and the baryonic overdensity respectively. The radiointerferometric observations would measure only the fluctuations in , e.g., the power spectrum defined as
(8) 
where is Fourier transform of . As can be seen from equation (6), the fluctuations in the signal are sourced by either fluctuations in the neutral hydrogen density field or in the spin temperature . Unless one is interested in the very early stages of reionization (cosmic dawn), the IGM can be taken to be sufficiently heated by Xrays and the coupled to the gas temperature through the Ly coupling (Wouthuysen, 1952; Field, 1959). In that case and hence the signal is simply , i.e., the signal will simply follow the neutral hydrogen distribution.
The above expressions do not account for the line of sight peculiar velocities of the gas which can make the power spectrum anisotropic. Other sources of anisotropies are the light cone effect and the AlcockâPaczynski effect. In that case, the power spectrum can be denoted as , where is the direction cosine between the line of sight and . It is expected that the first generation of 21 cm experiments would measure the spherically averaged power spectrum .
A related quantity of interest which can be measured from the observation is the multifrequency angular power spectrum (MAPS) defined as (Datta et al., 2007)
(9) 
where is the comoving distance to the redshift and has a magnitude . The quantity is essentially the twodimensional power spectrum on a plane at the distance from the observer.
3.1 Analytical models
Analytical models of reionization are based on modelling the size distribution of ionized bubbles around galaxies, which then can be extended to obtain the power spectrum. The simplest model would be to approximate the ionization field as a collection of nonoverlapping sphere of fixed radius (Bharadwaj and Ali, 2005; Datta et al., 2007). The values of and the number density of such bubbles can be chosen to obtain a particular ionized fraction at a given epoch. The quantity for different values of for such a model is shown in Fig 2. One can see that the signal is enhanced compared to the underlying dark matter fluctuations at angular scales larger than the bubble size. It peaks around a value of which is inversely proportional to the angular size of the bubbles. The frequency decorrelation function as function of is plotted in Fig 3 which shows that the fluctuations decorrelate quite rapidly particularly for large . This decorrelation can be useful in separating the cosmological signal from other astrophysical foregrounds which have relatively smooth frequency spectra.
These simple models of ionized bubbles do not account for the overlap and thus are valid only when the ionized fraction is very small. It is possible to account for such overlaps using the excursion set approach as was proposed by Furlanetto et al. (2004), which we refer to as FZH04. In this approach the condition for a spherical region of radius having a density contrast (linearly extrapolated to present epoch) to be fully ionized is given by a condition on the collapsed fraction in the region
(10) 
The condition for a region to be “selfionized” can be converted into a condition in terms of the density contrast , and then problem reduces to the one for a barrier crossing. An improvement to the above model was proposed by Paranjape and Choudhury (2014) by accounting for the fact that haloes would preferentially form near the density peaks (Musso and Sheth, 2012; Paranjape and Sheth, 2012), known as the excursion sets peak (ESP) model. Fig 4 shows the bubble size distribution for the two models FZH04 and ESP for fixed values of the ionized fraction . The main effect of the ESP model is to predict bubbles of relatively larger sizes, a fact which was confirmed while calculating the bubble distribution from seminumerical calculations. This implies that the power spectrum predicted from the ESP model would be very different from the earlier calculations, a fact which can play an important role in interpreting the 21 cm data.
3.2 Seminumerical models
Although the excursionset based analytical calculations provide a reasonable description of the growth of HII regions during reionization, it is difficult to incorporate various complexities into the analytical framework e.g., the nonsphericity of bubble during overlap, the effects of line of sight peculiar velocities, quantifying the cosmic variance. It is often useful to have realisations of the ionization field so that one can construct realistic radio maps and other quantities relevant to the observations. One can use the radiative transfer simulations to generate such maps, however, they are often not suited for exploring the parameter space.
A compromise has been proposed where it is possible to generate HI maps without running the full simulation. These seminumerical methods are based on the excursion set formalism discussed in the previous section and can be used for generating reasonably large volumes in a small amount of time (Mesinger and Furlanetto, 2007; Geil and Wyithe, 2008). One of these methods is based on generating the DM density field using the perturbation theory, calculating the collapsed fraction within each grid cell (of size and density contrast ) using the analytical expression for the conditional mass function and then generating the ionization field using the excursion set formalism (Mesinger et al., 2011). A slightly different approach is to run a full dark matter only body simulation and identify the haloes using a suitable groupfinder algorithm (Zahn et al., 2007; Choudhury et al., 2009). The method for generating the ionization field remains the same, i.e., using the excursion set formalism.
To understand how these seminumerical models work, let us assume that we have a realisation of the density field along with the collapsed fraction at each grid point of the simulation box. The basic idea is to compute the spherically averaged collapsed fraction for each grid point in the box for a wide range of values. If the selfionization condition (10) is satisfied for any , then the grid point is flagged as ionized. Points which fail to satisfy the condition are assigned a neutral fraction , where is the size of the grid cell typically set by the resolution of the map. The procedure is repeated for all points in the simulation volume. The method thus provides a realisation of the ionization field for a given value of .
One important physical process which is not taken into account in the above discussion is the recombination. A significant fraction of the ionizing photons may be lost in ionizing the recombined hydrogen atoms, which may lead to very different HI topology. In the simplest case where recombinations are taken to be spatially homogeneous, the effect can be absorbed in the definition of . In other words, one can rewrite the ionization condition as
(11) 
where is the average number of recombinations per hydrogen atom. Since is assumed to be independent of the spatial point under consideration, the above equation is simply the earlier ionization condition (10) with the redefinition .
In reality, however, the recombinations are not homogeneous. In fact, the high density regions will tend to recombine faster and will be able to “shield” themselves from the ionizing radiation. These regions would act as “sinks” of ionizing photons. One way of implementing the effect of such sinks into the seminumerical models is by including an additional condition for ionization (Choudhury et al., 2009)
(12) 
where is the number density of ionizing photons averaged over the region of radius , is the recombination timescale, is the timescale over which the recombination term has significant contribution with being the Hubble time and is the comoving size of the absorbing region, usually taken to be the local Jeans scale. The effect of these recombinations is that the reionization becomes outsidein once a substantial fraction of the IGM is ionized (in contrast to models without sinks where the process is always insideout). The presence of sinks predict smaller amplitude of fluctuations at scales Mpc accessible to first generation of 21 cm experiments.
This model of Choudhury et al. (2009) was used by Majumdar et al. (2013) to study the effect of peculiar velocities on the signal. It was proposed that the magnitude and nature of the ratio between the quadrupole and monopole moments of the power spectrum () can be a possible probe for the epoch of reionization. The same seminumerical models have also been used for showing that the 21 cm anisotropy is best measured by the quadrupole moment of the power spectrum, i.e., it evolves predictably as a function of (Majumdar et al., 2016).
Since the seminumerical simulations are quite fast in terms of their computing time, they are suitable for studying various statistical properties of the signal which may require large number of realisations. One such application is to study the effects of nonGaussianity of the ionization field on error predictions for the power spectrum (Mondal et al., 2015, 2016). Unlike the Gaussian case where SNR , being the number Fourier modes in a given bin, the SNR has upper limit which cannot be exceeded by increasing (Mondal et al., 2015) which can be seen in Fig 6. This could have severe implications for estimating the sensitivities for the future 21 cm experiments. In fact, one can see from Fig 7 that the error covariance is not diagonal any more and the coupling between different modes cannot be neglected for error estimations (Mondal et al., 2016).
3.3 Numerical simulations
The full complexities of the radiative transfer through the clumpy IGM can only be taken into account through the numerical solution of the radiative equation (4). However, radiative transfer simulations are still computationally extremely challenging and hence the equation is usually solved under reasonable approximations. One such algorithm has been implemented in the radiative transfer code “Conservative Causal Raytracing method” (CRAY) which works by tracing rays from all sources and iteratively solving for the time evolution of the ionized hydrogen fraction (Mellema et al., 2006; Iliev et al., 2006b). It turns out that the ionization fields generated by the CRAY have properties which are quite similar to the seminumerical calculations described in the previous section (Majumdar et al., 2014). Fig 5 shows the 21 cm power spectrum obtained from the CRAY compared with various seminumerical schemes. The match turns out to be quite reasonable thus showing that one can possibly use the fast seminumerical models to generate the 21 cm maps.
An alternate, relatively faster, method is to use the density field and haloes from a dark matter only body simulation, and postprocess with a spherically symmetric onedimensional radiative transfer algorithm (Thomas and Zaroubi, 2008; Thomas et al., 2009). The main idea is to generate spherically symmetric 21 cm patterns around individual sources (galaxies) and then account for the overlap of such regions suitably. This method has been implemented by Ghara et al. (2015a) whose code not only accounts for ionization of hydrogen and helium, but also the effects of Xray heating and Ly radiation. Inclusion of these effects allows one to estimate the spin temperature at different points selfconsistently and hence study the effect of fluctuations in . Some sample outputs from the code are shown in Fig 8 for three different models of treating the spin temperature fluctuations. The corresponding evolution of the 21 cm signal fluctuations at a scale Mpc is shown in Fig 9. If we concentrate on model C where the heating and the Ly coupling are calculated selfconsistently, we find that the evolution shows three distinct “peaks”. They arise from different physical processes, i.e., the one at the lowest redshift arises from fluctuations in , the second peak at corresponds to fluctuations in the Xray heating and the one at is from the Ly coupling fluctuations. Measuring these peaklike features could be quite important in understanding the physical processes at early times. It should however be kept in mind that various line of sight effects may affect the amplitude of these peaks and hence care should be maintained while interpreting the results (Ghara et al., 2015a, c).
4 Future outlook
There has been tremendous progress in modelling the physics of cosmic dawn and reionization in recent times. It is now possible to construct models that are consistent with all available observations related to reionization (Mitra et al., 2015). The prediction of the 21 cm signal is possibly a greater challenge and there are various approaches in modelling this.
The construction of the SKAlow phase 1 will allow us to probe the 21 cm cosmological signal with unprecedented sensitivity. Fig 10 compares the expected sensitivity of SKA1low with some of the existing facilities. As one can see, the noise variance for the SKA1low will be almost a order of magnitude better than any of the existing ones, thus allowing us to probe the high universe in much more detail.
Given this advancement, we require equally advanced techniques in modelling the signal so as to interpret the data accurately. Some of the directions in which the models discussed in the article can be improved are as follows:

Since the physical processes at high redshifts are uncertain, they need to be parametrized by a number of free parameters (and functions). The analytical calculations could be quite helpful in probing the space of these unknown parameters as they are fast. The challenge then would be to improve the analytical models and make them as accurate as possible, most likely by comparing with seminumerical and radiative transfer simulations. There also remain some conceptual issues with these excursion set based approaches, e.g., violation of photon conservation (Paranjape et al., 2015).

One important physical process which is still not accounted for satisfactorily is the effect of high density regions or sinks of ionizing photons. Appropriate selfconsistent and accurate treatments of these regions are required for analytical, seminumerical and radiative transfer simulations. This could turn out to be a major challenge for reionization models in the near future (Mesinger et al., 2015; Choudhury et al., 2015; Shukla et al., 2016).

In addition to the 21 cm experiments, there would be instruments at other wave bands which are likely to come up in the near future probing the high universe, e.g., JWST, TMT and so on. The models of galaxy formation as well as HI distribution should be sufficiently broad so as to explain all these observations simultaneously. This would be very important for understanding the EoR.

The sensitivities of the SKA1low would not only probe the power spectrum of HI fluctuations, but would also allow us to image the HI distribution. It is thus important to explore the new information one can obtain through these images.

With the improvement in the noise properties of the instrument, it becomes important to devise advanced statistical estimators which can be used for obtaining the relevant quantities of interest. Thus a better modelling of the system noise coupled with improved reionization models is required in the near future.
References
 Barkana and Loeb (2001) Barkana, R. and Loeb, A.: 2001, Phys. Rep. 349, 125
 Bharadwaj and Ali (2005) Bharadwaj, S. and Ali, S. S.: 2005, MNRAS 356, 1519
 Bi and Davidsen (1997) Bi, H. and Davidsen, A. F.: 1997, ApJ 479, 523
 Bolton and Becker (2009) Bolton, J. S. and Becker, G. D.: 2009, MNRAS 398, L26
 Bond et al. (1991) Bond, J. R., Cole, S., Efstathiou, G., and Kaiser, N.: 1991, ApJ 379, 440
 Carilli (2015) Carilli, C.: 2015, Advancing Astrophysics with the Square Kilometre Array (AASKA14) p. 171
 Choudhury (2009) Choudhury, T. R.: 2009, Current Science 97, 841
 Choudhury and Ferrara (2005) Choudhury, T. R. and Ferrara, A.: 2005, MNRAS 361, 577
 Choudhury et al. (2009) Choudhury, T. R., Haehnelt, M. G., and Regan, J.: 2009, MNRAS 394, 960
 Choudhury et al. (2001a) Choudhury, T. R., Padmanabhan, T., and Srianand, R.: 2001a, MNRAS 322, 561
 Choudhury et al. (2015) Choudhury, T. R., Puchwein, E., Haehnelt, M. G., and Bolton, J. S.: 2015, MNRAS 452, 261
 Choudhury et al. (2001b) Choudhury, T. R., Srianand, R., and Padmanabhan, T.: 2001b, ApJ 559, 29
 Datta et al. (2007) Datta, K. K., Choudhury, T. R., and Bharadwaj, S.: 2007, MNRAS 378, 119
 Davis et al. (1985) Davis, M., Efstathiou, G., Frenk, C. S., and White, S. D. M.: 1985, ApJ 292, 371
 Field (1959) Field, G. B.: 1959, ApJ 129, 536
 Furlanetto et al. (2006) Furlanetto, S. R., Oh, S. P., and Briggs, F. H.: 2006, Phys. Rep. 433, 181
 Furlanetto et al. (2004) Furlanetto, S. R., Zaldarriaga, M., and Hernquist, L.: 2004, ApJ 613, 1
 Geil and Wyithe (2008) Geil, P. M. and Wyithe, J. S. B.: 2008, MNRAS 386, 1683
 Ghara et al. (2015a) Ghara, R., Choudhury, T. R., and Datta, K. K.: 2015a, MNRAS 447, 1806
 Ghara et al. (2015b) Ghara, R., Choudhury, T. R., and Datta, K. K.: 2015b, ArXiv eprints
 Ghara et al. (2015c) Ghara, R., Datta, K. K., and Choudhury, T. R.: 2015c, MNRAS 453, 3143
 Iliev et al. (2006a) Iliev, I. T., Ciardi, B., Alvarez, M. A., Maselli, A., Ferrara, A., Gnedin, N. Y., Mellema, G., Nakamoto, T., Norman, M. L., Razoumov, A. O., Rijkhorst, E.J., Ritzerveld, J., Shapiro, P. R., Susa, H., Umemura, M., and Whalen, D. J.: 2006a, MNRAS 371, 1057
 Iliev et al. (2006b) Iliev, I. T., Mellema, G., Pen, U.L., Merz, H., Shapiro, P. R., and Alvarez, M. A.: 2006b, MNRAS 369, 1625
 Koopmans et al. (2015) Koopmans, L., Pritchard, J., Mellema, G., Aguirre, J., Ahn, K., Barkana, R., van Bemmel, I., Bernardi, G., Bonaldi, A., Briggs, F., de Bruyn, A. G., Chang, T. C., Chapman, E., Chen, X., Ciardi, B., Dayal, P., Ferrara, A., Fialkov, A., Fiore, F., Ichiki, K., Illiev, I. T., Inoue, S., Jelic, V., Jones, M., Lazio, J., Maio, U., Majumdar, S., Mack, K. J., Mesinger, A., Morales, M. F., Parsons, A., Pen, U. L., Santos, M., Schneider, R., Semelin, B., de Souza, R. S., Subrahmanyan, R., Takeuchi, T., Vedantham, H., Wagg, J., Webster, R., Wyithe, S., Datta, K. K., and Trott, C.: 2015, Advancing Astrophysics with the Square Kilometre Array (AASKA14) p. 1
 Majumdar et al. (2013) Majumdar, S., Bharadwaj, S., and Choudhury, T. R.: 2013, MNRAS 434, 1978
 Majumdar et al. (2016) Majumdar, S., Jensen, H., Mellema, G., Chapman, E., Abdalla, F. B., Lee, K.Y., Iliev, I. T., Dixon, K. L., Datta, K. K., Ciardi, B., Fernandez, E. R., Jelić, V., Koopmans, L. V. E., and Zaroubi, S.: 2016, MNRAS 456, 2080
 Majumdar et al. (2014) Majumdar, S., Mellema, G., Datta, K. K., Jensen, H., Choudhury, T. R., Bharadwaj, S., and Friedrich, M. M.: 2014, MNRAS 443, 2843
 Mellema et al. (2006) Mellema, G., Iliev, I. T., Alvarez, M. A., and Shapiro, P. R.: 2006, New A 11, 374
 Mellema et al. (2013) Mellema, G., Koopmans, L. V. E., Abdalla, F. A., Bernardi, G., Ciardi, B., Daiboo, S., de Bruyn, A. G., Datta, K. K., Falcke, H., Ferrara, A., Iliev, I. T., Iocco, F., Jelić, V., Jensen, H., Joseph, R., Labroupoulos, P., Meiksin, A., Mesinger, A., Offringa, A. R., Pandey, V. N., Pritchard, J. R., Santos, M. G., Schwarz, D. J., Semelin, B., Vedantham, H., Yatawatta, S., and Zaroubi, S.: 2013, Experimental Astronomy 36, 235
 Mesinger et al. (2015) Mesinger, A., Aykutalp, A., Vanzella, E., Pentericci, L., Ferrara, A., and Dijkstra, M.: 2015, MNRAS 446, 566
 Mesinger and Furlanetto (2007) Mesinger, A. and Furlanetto, S.: 2007, ApJ 669, 663
 Mesinger et al. (2011) Mesinger, A., Furlanetto, S., and Cen, R.: 2011, MNRAS 411, 955
 MiraldaEscudé et al. (2000) MiraldaEscudé, J., Haehnelt, M., and Rees, M. J.: 2000, ApJ 530, 1
 Mitra et al. (2011) Mitra, S., Choudhury, T. R., and Ferrara, A.: 2011, MNRAS 413, 1569
 Mitra et al. (2012) Mitra, S., Choudhury, T. R., and Ferrara, A.: 2012, MNRAS 419, 1480
 Mitra et al. (2015) Mitra, S., Choudhury, T. R., and Ferrara, A.: 2015, MNRAS 454, L76
 Mondal et al. (2016) Mondal, R., Bharadwaj, S., and Majumdar, S.: 2016, MNRAS 456, 1936
 Mondal et al. (2015) Mondal, R., Bharadwaj, S., Majumdar, S., Bera, A., and Acharyya, A.: 2015, MNRAS 449, L41
 Musso and Sheth (2012) Musso, M. and Sheth, R. K.: 2012, MNRAS 423, L102
 Paranjape and Choudhury (2014) Paranjape, A. and Choudhury, T. R.: 2014, MNRAS 442, 1470
 Paranjape et al. (2015) Paranjape, A., Choudhury, T. R., and Padmanabhan, H.: 2015, ArXiv eprints
 Paranjape and Sheth (2012) Paranjape, A. and Sheth, R. K.: 2012, MNRAS 426, 2789
 Press and Schechter (1974) Press, W. H. and Schechter, P.: 1974, ApJ 187, 425
 Pritchard and Loeb (2012) Pritchard, J. R. and Loeb, A.: 2012, Reports on Progress in Physics 75(8), 086901
 Sheth and Tormen (1999) Sheth, R. K. and Tormen, G.: 1999, MNRAS 308, 119
 Shukla et al. (2016) Shukla, H., Mellema, G., Iliev, I. T., and Shapiro, P. R.: 2016, MNRAS 458, 135
 Thomas and Zaroubi (2008) Thomas, R. M. and Zaroubi, S.: 2008, MNRAS 384, 1080
 Thomas et al. (2009) Thomas, R. M., Zaroubi, S., Ciardi, B., Pawlik, A. H., Labropoulos, P., Jelić, V., Bernardi, G., Brentjens, M. A., de Bruyn, A. G., Harker, G. J. A., Koopmans, L. V. E., Mellema, G., Pandey, V. N., Schaye, J., and Yatawatta, S.: 2009, MNRAS 393, 32
 Wouthuysen (1952) Wouthuysen, S. A.: 1952, AJ 57, 31
 Zahn et al. (2007) Zahn, O., Lidz, A., McQuinn, M., Dutta, S., Hernquist, L., Zaldarriaga, M., and Furlanetto, S. R.: 2007, ApJ 654, 12