# Epoch of Reionization modelling and simulations for SKA

###### Abstract

In this chapter we provide an overview of the current status of the simulations and modelling of the Cosmic Dawn and Epoch of Reionization. We discuss the modelling requirements as dictated by the characteristic scales of the problem and the SKA instrumental properties and the planned survey parameters. Current simulations include most of the relevant physical processes. They can follow the full nonlinear dynamics and are now reaching the required scale and dynamic range, although small-scale physics still needs to be included at sub-grid level. However, despite a significant progress in developing novel numerical methods for efficient utilization of current hardware they remain quite computationally expensive. In response, a number of alternative approaches, particularly semi-analytical/semi-numerical methods, have been developed. While necessarily more approximate, if appropriately constructed and calibrated on simulations they could be used to quickly explore the vast parameter space available. Further work is still required on including some physical processes in both simulations and semi-analytical modelling. This hybrid approach of fast, approximate modelling calibrated on numerical simulations can then be used to construct large libraries of reionization models for reliable interpretation of the observational data.

Epoch of Reionization modelling and simulations for SKA

Suman Majumdar, Garrelt Mellema

Department of Astronomy & Oskar Klein Centre, AlbaNova, Stockholm University, SE-106 91 Stockholm, Sweden

\abstract@cs

## 1 Introduction

The Cosmic Dark Ages, Cosmic Dawn and Epoch of Reionization, span altogether about a billion years between the CMB last scattering surface at redshift and the complete ionisation of hydrogen in the inter-galactic medium (IGM) at . This phase is quite distant from us and difficult to study, remaining one of the last poorly understood epochs of the Universe. During this time, the very first stars and galaxies formed, then gradually ionized the intergalactic medium and enriched it with metals, thereby laying the foundations for the cosmic structures we see today. Better knowledge of these epochs is therefore of key importance for understanding galaxy formation and evolution.

The main obstacle to further progress is the scarcity of observational data, which currently mostly probes the tail-end of reionization (both Ly- source surveys and IGM absorption lines probe low neutral fractions, due to the high optical depth to resonant scattering of such radiation by even small amounts of neutral hydrogen) or are integral measures of its history (Cosmic Microwave Background, CMB, optical depth and polarization; kinetic Sunyaev-Zel’dovich effect, kSZ; Near Infrared Background, NIRB). The redshifted 21-cm signal promises to provide full 3D tomographic observations of the intergalactic medium throughout, and possibly even before reionization. The first generation of experiments, currently ongoing, will most likely provide just a statistical detection of the signal, e.g. power spectra, variance and/or probability distribution functions (PDFs) and their higher moments. In contrast, SKA should be able to also perform imaging and 3D mapping, as well as go much deeper in redshift, due to its far superior sensitivity, thereby completely transforming reionization research.

The shortage of observational constraints has meant that simulations have played and continue to play a larger role than in other areas, since we need to rely on modelling to understand even the basic properties of the expected signals. This in turn steers the design of the observational experiments targeted at this science. Furthermore, better understanding of the expected 21-cm signatures and their cross-correlations with other observations is important for reliably confirming any detection, since the signal is weak compared to the foregrounds which are several orders of magnitude stronger. Separating the signal and the foregrounds will not be trivial (see e.g. Chapman et al. 2015, reference: PoS(AASKA14)005) and subtraction may leave contaminating residuals mimicking the signal. The simulations will have a further important role to play in the interpretation of any detections. The reason for this is the disconnect in scales between the objects we want to study, i.e. the first galaxies, which are very small and faint, thus largely impossible to observe directly, and the actual 21-cm and other signals we will detect, which typically relate to large-scale patterns such as the reionization patchiness. The latter are caused by large number of clustered, individually dim sources. Therefore, the interpretation of any signals in terms of the nature and properties of the first galaxies is non-trivial and requires detailed modelling.

In recent years two basic approaches emerged in modelling the reionization patchiness - either by direct numerical simulations (e.g. Gnedin & Ostriker, 1997; Ciardi et al., 2000; Gnedin, 2000a; Ricotti, Gnedin & Shull, 2002; Iliev et al., 2006, 2014) or through semi-analytical/semi-numerical modelling (e.g. Furlanetto, Zaldarriaga & Hernquist, 2004; Zahn et al., 2007; Alvarez et al., 2009; Mesinger & Furlanetto, 2007a; Choudhury, Haehnelt & Regan, 2009; Santos et al., 2010; Battaglia et al., 2013), along with some intermediate methods involving simplified simulations (Kohler, Gnedin & Hamilton, 2007; Thomas & et al., 2009). These different lines of attack are in practice highly complimentary and each has its advantages as well as disadvantages, as discussed in more detail in the next section.

## 2 Overview of current modelling

### 2.1 Radiative transfer simulations

The simulations of cosmic reionization fall in two broad types, as a consequence of the huge range of scales involved, which in practice cannot be covered in a single simulation on current computer hardware. The first type consists of small volume, high-resolution simulations which can be used to study in detail the formation of early galaxies and their radiative and supernova feedback on the gas with fully-coupled radiative hydrodynamic simulations. Such simulations cannot capture a large enough volume to represent the global reionization process or its observational signatures (e.g. Gnedin & Ostriker, 1997; Ciardi et al., 2000; Gnedin, 2000a; Ricotti, Gnedin & Shull, 2002). The second type consists of large-scale simulations, which instead follow volumes sufficiently large to study the global evolution, but lack resolution to directly resolve the small-scale physics (e.g. Iliev et al., 2006; Trac & Cen, 2007; Iliev et al., 2014) (see Figs. 1 and 2). This type of simulations decouple the radiative transfer, star formation and feedback processes from the underlying structure formation. The latter is done first, typically using large-scale N-body simulations, which is then followed up by detailed radiative transfer and non-equilibrium chemistry simulations. This approach allows for much larger dynamic range, and thus the larger volumes required for understanding the global reionization process and its observational signatures. In this case the detailed physical processes at unresolved scales have to be included using subgrid recipes, which are themselves based on high-resolution simulations or other modelling (e.g. Iliev et al., 2007; McQuinn et al., 2007; Ahn et al., 2012, 2014). Both types of simulations have played a major role in the significant recent advances in this area. We now have a good understanding of the characteristic scales of the process (Friedrich et al., 2011; Iliev et al., 2014), the physical processes affecting the patchiness and the various observational features which might help to discriminate between different EoR models (e.g. Wise & Abel, 2008; Baek et al., 2009; Ahn et al., 2012; Iliev et al., 2012; Alvarez & Abel, 2012; O’Leary & McQuinn, 2012; Park et al., 2013; Jensen et al., 2013; Sobacchi & Mesinger, 2014).

### 2.2 Semi-numerical simulations

\pdfdestnameref-seminanalyt:sect fith

Some of the early approaches to characterise the 21cm signal from the epoch of reionization have relied on more simple analytical models, mostly based on the excursion-set approach. These can be very useful to quickly generate the signal, in particular the power spectrum of 21-cm brightness temperature fluctuations (Furlanetto, Zaldarriaga & Hernquist, 2004; Furlanetto, McQuinn & Hernquist, 2006; Sethi, 2005). The speed and ease of use, allows for straightforward tests of the ability of a given experiment to constraint cosmological and astrophysical parameters (Santos & Cooray, 2006; McQuinn et al., 2006; Mao et al., 2008) and more generally to rapidly explore the vast parameter space available. The analytical models have also been useful to understand the possible contributions to the 21 cm signal at high redshifts (Barkana & Loeb, 2005; Pritchard, 2007), although it is at low redshifts () that they seem to provide a better description of the 21-cm 2-point correlation function (Santos et al., 2008). However, these models have several issues in properly dealing with the spatial distribution of the reionization process, such as bubble overlap and ignore complicated astrophysics during reionization. Such complexities are better handled in full (but expensive) numerical simulations, as was discussed above.

An intermediate approach that has become more popular are the so called semi-numerical 21cm simulations which try to merge the speed and ease of use of analytical models with the ability to follow the spatial evolution of reionization in more detail as provided by numerical simulations. The basic feature of these semi-numerical simulations is that they start from a random Gaussian field realization of the cosmological density field, or a full N-body simulation data, but then replace the time-consuming radiative transfer (RT) with an excursion-set approach (Furlanetto, Zaldarriaga & Hernquist, 2004) that determines if a region is ionized by comparing the number of ionizing photons to the number of atoms (plus recombinations) inside it. As in conventional cosmological RT simulations, this algorithm can be applied to discrete halo source fields, obtained either from traditional N-body simulations (e.g. Zahn et al., 2007; Choudhury, Haehnelt & Regan, 2009; Zahn et al., 2011) or through faster methods involving excursion-set formalism and perturbation theory (Mesinger & Furlanetto, 2007a; Santos et al., 2010; Mesinger, Furlanetto & Cen, 2011; Alvarez et al., 2009; Geil & Wyithe, 2008). The latter approach results in some reduction in accuracy (see Fig. 3; Zahn et al., 2011; Majumdar et al., 2014), but substantially increases the achievable dynamic range. Some of these approaches have also been extended to simulating the earlier epochs of the 21cm signal: the Cosmic Dawn/Dark Ages (Mesinger, Furlanetto & Cen, 2011; Santos et al., 2010). During the Dark Ages, the 21cm signal is governed by soft UV and X-ray photons, which can have mean free paths of hundreds of Mpc, requiring large volumes, which makes it a challenge to implement this in a complete radiative transfer/gas dynamics algorithm.

Semi-numerical simulations of reionization have been extensively tested against numerical simulations, with good agreement on moderate to large scales (Mpc), relevant for most SKA science (Mesinger, Furlanetto & Cen, 2011; Zahn et al., 2011; Majumdar et al., 2014). The approximations start to break down however on non-linear scales, and care should be taken by comparing against more accurate methods when extending their application to any non-standard problems. Nevertheless, their speed and efficiency allows for rapid exploration of the large-parameter space of uncertainties, ushering in an exciting era of 21cm astrophysical parameter studies (e.g. Mesinger, Ewall-Wice & Hewitt, 2014; Pober et al., 2013).

The analysis of Zahn et al. (2011) has been largely focused on the comparison of the spherically average power spectrum of the 21-cm signal from different simulations. Their study suggested that these semi-numerical techniques can mimic the 21-cm power spectrum from that of a radiative transfer simulation with accuracy at length scales relevant for the present and future surveys. This study was extended in more recent works (Mesinger, Furlanetto & Cen, 2011; Majumdar et al., 2014) to include the comparison of EoR history, the redshift space anisotropies in the power spectrum and the morphology of the ionization fields (Fig. 3). Those analyses suggest that the halo based semi-numerical simulations can predict the history of EoR with an accuracy of when compared with radiative transfer simulations for length scales , whereas the same accuracy for the Press-Schechter based technique is at same length scale range. They have also shown that the anisotropies in the power spectrum of the 21-cm signal due to redshift space distortions can also be mimicked by both type of the semi-numerical simulations with an accuracy of 85% or more, which is well within the noise uncertainties of hr LOFAR or hr SKA1-LOW observation at 150 MHz, provided that the actual peculiar velocities are used to introduce the redshift space distortions rather than using a perturbative technique. It is also found that the cross-correlation between the ionization fields from the halo based semi-numerical technique and the radiative transfer simulations is more than during almost the entire period of EoR, for length scales , whereas the same correlation between the conditional Press-Schechter formalism and the radiative transfer simulations, while good in some regimes, can reach values as low as , specifically at late stages of EoR for the same length scale range.

These comparisons however, are also limited as the reionization scenario/model that has been compared in these cases is relatively simplistic in nature. The reionization history and the nature of the 21-cm signal from this epoch may differ depending on various other factors. One such important factor is the effect of radiative feedback on the star formation in low mass halos ( in mass). This will most likely affect their star formation efficiency, although the details remain unclear (Couchman & Rees, 1986; Rees, 1986; Efstathiou, 1992; Thoul & Weinberg, 1995, 1996; Gnedin, 2000b; Kitayama et al., 2000; Dijkstra et al., 2004; Hoeft et al., 2006; Okamoto, Gao & Theuns, 2008). There has been some effort to include this in numerical (Iliev et al., 2007, 2012; Ahn et al., 2014) as well as in semi-numerical simulations (Sobacchi & Mesinger, 2013). Similarly, another such important issue is the effect of enhanced recombination in the Lyman Limit Systems and other small-scale structures, which are not resolvable in any of these simulations and thus require detailed sub-grid modelling (Ciardi et al. 2006; Choudhury, Haehnelt & Regan 2009; Sobacchi & Mesinger 2014, Koda et al. in prep., Shukla et al., in prep.). Also, the effect of X-ray heating on the spin temperature evolution at the very early stages of reionization (Mesinger, Ferrara & Spiegel, 2013; Ghara, Choudhury & Datta, 2014) is another important factor which has not been included in these comparison studies.

## 3 SKA simulations

### 3.1 Basic Simulation Requirements

The volume and resolution required for proper numerical modelling of the reionization process are dictated by both its intrinsic characteristic scales (especially the typical sizes of ionized and neutral patches) and the instrumental parameters (FOV, beam size and bandwidth). The simulation volume should be large enough to faithfully sample all relevant scales, while the numerical resolution should be such that every beam/bandwidth is sampled with a sufficient number of elements so as to avoid discretization effects and other numerical artifacts.

The nominal SKA1-LOW Epoch of Reionization survey at the frequencies relevant for reionization is proposed to have FOV (FWHM) of (at MHz) to 10 degrees (at MHz), which corresponds to Mpc-1.5 Gpc, with maximum resolution of arcmin, or roughly Mpc. Therefore simulations of the full SKA1-LOW FOV with sufficient resolution should have volume of at least 500 Mpc and grid sizes of at least several thousand cells per side. These simulation parameters are now just becoming achievable on current hardware. These are of course only very rough estimates and they only relate to the 21-cm sky maps and statistical studies (e.g. power spectra, or PDFs). For other purposes, for example 21-cm absorption studies against bright radio galaxies (see e.g. Ciardi et al. (2015) reference: PoS(AASKA14)006), a much higher simulation resolution will be required, corresponding to frequency channels of kHz, or cells of kpc comoving. This is unrealistic when combined with Mpc volumes unless adaptive-mesh refinement (AMR) techniques are employed, or else smaller simulation volumes have to be used for such work.

Apart from the purely instrumental requirements discussed above, the simulation parameters are also dictated by the characteristic scales of the reionization itself and the various processes that occur at different stages and will be of interest to study. Before the first UV sources form, the 21-cm fluctuations are dictated by the underlying density fluctuations, the temperature of the IGM, as well as additional effects like baryon-dark matter displacement (Tseliakhovich & Hirata, 2010, see also Maio et al. 2015, reference: PoS(AASKA14)009) and star formation suppression in minihaloes due to Lyman-Werner bands photons. Once the first ionizing sources appear, they propagate ionization fronts into the IGM initially forming small, local HII regions, which then quickly percolate and merge locally into larger ones. The same first sources also produce copious amounts of soft-UV radiation and likely some X-rays, which heat the IGM gas and decouple its spin temperature from the CMB, making it detectable. Outside dense regions/very high redshifts, the 21-cm line decoupling from the CMB also requires sufficiently strong Lyman- background. All these different processes impose some characteristic scales on the 21-cm fluctuations. These are typically long-range modulations ranging from tens of Mpc (Lyman-), to Mpc (baryon-dark matter displacement and Lyman-Werner) and up to hundreds of Mpc (X-rays). Later-on, likely below redshift when the reionization becomes more widespread and is driven by larger, atomically-cooling halos, the 21-cm fluctuations due to these early backgrounds become less significant, baryon-dark matter displacement decreases and the dominant process determining the 21-cm fluctuations become the ionization patchiness. Its typical scales and geometry depend on the abundance and clustering of the main sources driving the process. Low-mass galaxies are exponentially more abundant, but are also liable to get their star formation suppressed by the Jeans mass filtering due to photoheathing, as well as other radiative and mechanical feedback mechanisms. The details of these processes are still not firmly established and are currently subject of active research.

The ionized regions growth continues throughout reionization (see Fig. 2 for an example). The characteristic patch sizes are dictated by the abundance, clustering and typical luminosities of the dominant ionizing sources, as well as various effects which modify the growth like recombinations in the IGM gas, Lyman-limit systems and other photon sinks. The intrinsic source clustering results from the statistics of the initial Gaussian random noise density fluctuations and depends strongly on redshift and the typical mass of the halos hosting the dominant sources. In CDM the cosmological structures form hierarchically, with the smallest ones forming first and then growing and merging over time to form larger ones. Consequently, low-mass galaxies likely dominated the ionizing photon output throughout reionization. Figure 4 (left) shows the relative abundance of halos of different mass vs. redshift. At very high redshifts all halos are rare () and only below redshift the low-mass galaxies () become somewhat more common (). Consequently, these halos are strongly clustered, with bias with respect to the underlying mass distribution well above 1 (Figure 4, right). This strong clustering yields quick percolation of the individually small HII regions around each relatively weak source and thus the rapid growth to much larger scales. Furthermore, there is a notable power in the density fluctuations at fairly large scales, which required large simulation volumes to account for properly (Mesinger & Furlanetto, 2007b; Iliev et al., 2014). This additional power means that local photon output is modulated significantly, which is not modelled correctly in volumes smaller than about per side, resulting in artificial suppression of the redshifted 21-cm fluctuations, among others. Different measures of the characteristic scales of the patchiness all yield typical peak sizes of tens of comoving Mpc (see Figure 5), which corresponds to tens of arcmin (see also Sobacchi & Mesinger, 2014). In summary, accurate modelling of all these physical processes again points to requirements of large, at least several hundred Mpc per side, volumes and sufficient resolution to reliably identify, either directly or sub-grid, the main sources and sinks driving the reionization process.

### 3.2 Computational Requirements

As mentioned above, reliable guidance and interpretation of the observations in terms of deriving the properties of the first galaxies requires the creation of large libraries of models sampling the available parameters space. At present a few full radiative transfer codes can be scaled up to the volumes and resolutions required for full SKA1 EoR simulation discussed in the previous section. This can only be achieved using massively-parallel approach (in some cases using accelerators like Graphical Processing Units, GPUs, and Intel Phi multicore processors) on the largest available computers and a single simulation of this size requires up to tens of millions of core-hours. Even with the expected advances in computing technology it is unlikely that a large number of such simulations could be practically performed. It is therefore most likely that more approximate methods like semi-numerical modelling would have to be employed, guided by fewer detailed full simulations.

Both approaches also require large amounts of data storage. The large N-body structure formation simulations (needed for precision and accuracy even in semi-numerical models, as explained in § 2.2), alone require up to several PB of storage per simulation, with further significant storage required for the reionization data itself. Accessing and using this data within the community will require also investment in databases and other efficient methods for data sharing.

## 4 Summary

We presented an overview of the current status and future prospects for simulations and modelling of the Epoch of Reionization, with focus on the specific requirements for SKA1-LOW. As discussed in detail, extensive modeling is particularly important for Cosmic Dawn and EoR science compared to other areas of study because as yet there is very limited direct observational data to guide us. Two basic approaches are available - full numerical simulations and semi-analytical/semi-numerical modeling. The simulations can handle the nonlinear dynamics, feedback effects and complex geometry, but are fairly computationally expensive. This drawback is being alleviated by innovative, efficient radiative transfer methods, some of which are able to efficiently utilise the latest Petascale computing facilities and specialist hardware (e.g. GPU accelerators). This now allows simulations which were impossible just a few years ago, with volumes and resolution matching both the intrinsic scales of the reionization process and the FOV and resolution expected for the SKA1 EoR experiment. Nonetheless, full simulations remain quite expensive and consequently a significant effort has been invested in the development of semi-numerical modelling and other faster modelling methods. These approaches are by nature more approximate, generally use simplified physics and are missing the non-linear effects. However, when they are carefully constructed they still can retain many of the key features, while at the same time introducing the option of constructing large libraries of models to be used for interpreting the data, which is not feasible for full simulations. Much additional work is still required on these models, as well as on the numerical simulations, particularly in adding important physical processes which are still missing and in calibration and verification of semi-analytical models against numerical simulations. We have also outlined the requirements for SKA1-LOW-specific simulations and provided estimates on the resources which will be required for this effort.

## References

- Ahn et al. (2012) Ahn K., Iliev I. T., Shapiro P. R., Mellema G., Koda J., Mao Y., 2012, ApJL, 756, L16
- Ahn et al. (2014) Ahn K., Iliev I. T., Shapiro P. R., Srisawat C. B., 2014, ArXiv e-prints, 1407.2637
- Alvarez & Abel (2012) Alvarez M. A., Abel T., 2012, ApJ, 747, 126
- Alvarez et al. (2009) Alvarez M. A., Busha M., Abel T., Wechsler R. H., 2009, ApJL, 703, L167
- Baek et al. (2009) Baek S., di Matteo P., Semelin B., Combes F., Revaz Y., 2009, A.&A, 495, 389
- Barkana & Loeb (2005) Barkana R., Loeb A., 2005, ApJL, 624, L65
- Battaglia et al. (2013) Battaglia N., Trac H., Cen R., Loeb A., 2013, ApJ, 776, 81
- Choudhury, Haehnelt & Regan (2009) Choudhury T. R., Haehnelt M. G., Regan J., 2009, MNRAS, 394, 960
- Ciardi et al. (2000) Ciardi B., Ferrara A., Governato F., Jenkins A., 2000, MNRAS, 314, 611
- Ciardi et al. (2006) Ciardi B., Scannapieco E., Stoehr F., Ferrara A., Iliev I. T., Shapiro P. R., 2006, MNRAS, 366, 689
- Couchman & Rees (1986) Couchman H. M. P., Rees M. J., 1986, MNRAS, 221, 53
- Dijkstra et al. (2004) Dijkstra M., Haiman Z., Rees M. J., Weinberg D. H., 2004, ApJ, 601, 666
- Efstathiou (1992) Efstathiou G., 1992, MNRAS, 256, 43P
- Friedrich et al. (2011) Friedrich M. M., Mellema G., Alvarez M. A., Shapiro P. R., Iliev I. T., 2011, MNRAS, 413, 1353
- Furlanetto, McQuinn & Hernquist (2006) Furlanetto S. R., McQuinn M., Hernquist L., 2006, MNRAS, 365, 115
- Furlanetto, Zaldarriaga & Hernquist (2004) Furlanetto S. R., Zaldarriaga M., Hernquist L., 2004, ApJ, 613, 1
- Geil & Wyithe (2008) Geil P. M., Wyithe J. S. B., 2008, MNRAS, 386, 1683
- Ghara, Choudhury & Datta (2014) Ghara R., Choudhury T. R., Datta K. K., 2014, ArXiv e-prints
- Gnedin (2000a) Gnedin N. Y., 2000a, ApJ, 535, 530
- Gnedin (2000b) Gnedin N. Y., 2000b, ApJ, 542, 535
- Gnedin & Ostriker (1997) Gnedin N. Y., Ostriker J. P., 1997, ApJ, 486, 581
- Hoeft et al. (2006) Hoeft M., Yepes G., Gottlöber S., Springel V., 2006, MNRAS, 371, 401
- Iliev et al. (2014) Iliev I. T., Mellema G., Ahn K., Shapiro P. R., Mao Y., Pen U.-L., 2014, MNRAS, 439, 725
- Iliev et al. (2006) Iliev I. T., Mellema G., Pen U.-L., Merz H., Shapiro P. R., Alvarez M. A., 2006, MNRAS, 369, 1625
- Iliev et al. (2007) Iliev I. T., Mellema G., Shapiro P. R., Pen U.-L., 2007, MNRAS, 376, 534
- Iliev et al. (2012) Iliev I. T., Mellema G., Shapiro P. R., Pen U.-L., Mao Y., Koda J., Ahn K., 2012, MNRAS, 423, 2222
- Jensen et al. (2013) Jensen H. et al., 2013, MNRAS, 435, 460
- Kitayama et al. (2000) Kitayama T., Tajiri Y., Umemura M., Susa H., Ikeuchi S., 2000, MNRAS, 315, L1
- Kohler, Gnedin & Hamilton (2007) Kohler K., Gnedin N. Y., Hamilton A. J. S., 2007, ApJ, 657, 15
- Majumdar et al. (2014) Majumdar S., Mellema G., Datta K. K., Jensen H., Choudhury T. R., Bharadwaj S., Friedrich M. M., 2014, ArXiv e-prints, arXiv:1403.0941
- Mao et al. (2008) Mao Y., Tegmark M., McQuinn M., Zaldarriaga M., Zahn O., 2008, Phys Rev D, 78, 023529
- McQuinn et al. (2007) McQuinn M., Lidz A., Zahn O., Dutta S., Hernquist L., Zaldarriaga M., 2007, MNRAS, 377, 1043
- McQuinn et al. (2006) McQuinn M., Zahn O., Zaldarriaga M., Hernquist L., Furlanetto S. R., 2006, ApJ, 653, 815
- Mesinger, Ewall-Wice & Hewitt (2014) Mesinger A., Ewall-Wice A., Hewitt J., 2014, MNRAS; arXiv:1310.0465
- Mesinger, Ferrara & Spiegel (2013) Mesinger A., Ferrara A., Spiegel D. S., 2013, MNRAS, 431, 621
- Mesinger & Furlanetto (2007a) Mesinger A., Furlanetto S., 2007a, ApJ, 669, 663
- Mesinger & Furlanetto (2007b) Mesinger A., Furlanetto S., 2007b, ApJ, 669, 663
- Mesinger, Furlanetto & Cen (2011) Mesinger A., Furlanetto S., Cen R., 2011, MNRAS, 411, 955
- Okamoto, Gao & Theuns (2008) Okamoto T., Gao L., Theuns T., 2008, MNRAS, 390, 920
- O’Leary & McQuinn (2012) O’Leary R. M., McQuinn M., 2012, ApJ, 760, 4
- Park et al. (2013) Park H., Shapiro P. R., Komatsu E., Iliev I. T., Ahn K., Mellema G., 2013, ApJ, 769, 93
- Pober et al. (2013) Pober J. C., et al., 2013, The Astrophysical Journal Letters, 768, L36
- Pritchard (2007) Pritchard J. R., 2007, PhD thesis, California Institute of Technology
- Rees (1986) Rees M. J., 1986, MNRAS, 218, 25P
- Ricotti, Gnedin & Shull (2002) Ricotti M., Gnedin N. Y., Shull J. M., 2002, ApJ, 575, 33
- Santos et al. (2008) Santos M. G., Amblard A., Pritchard J., Trac H., Cen R., Cooray A., 2008, ApJ, 689, 1
- Santos & Cooray (2006) Santos M. G., Cooray A., 2006, Phys Rev D, 74, 083517
- Santos et al. (2010) Santos M. G., Ferramacho L., Silva M. B., Amblard A., Cooray A., 2010, MNRAS, 406, 2421
- Sethi (2005) Sethi S. K., 2005, MNRAS, 363, 818
- Sobacchi & Mesinger (2013) Sobacchi E., Mesinger A., 2013, MNRAS, 432, 3340
- Sobacchi & Mesinger (2014) Sobacchi E., Mesinger A., 2014, MNRAS, 440, 1662
- Thomas & et al. (2009) Thomas R. M., et al., 2009, MNRAS, 393, 32
- Thoul & Weinberg (1995) Thoul A. A., Weinberg D. H., 1995, ApJ, 442, 480
- Thoul & Weinberg (1996) Thoul A. A., Weinberg D. H., 1996, ApJ, 465, 608
- Trac & Cen (2007) Trac H., Cen R., 2007, ApJ, 671, 1
- Tseliakhovich & Hirata (2010) Tseliakhovich D., Hirata C., 2010, Phys Rev D, 82, 083520
- Wise & Abel (2008) Wise J. H., Abel T., 2008, ApJ, 685, 40
- Zahn et al. (2007) Zahn O., Lidz A., McQuinn M., Dutta S., Hernquist L., Zaldarriaga M., Furlanetto S. R., 2007, ApJ, 654, 12
- Zahn et al. (2011) Zahn O., Mesinger A., McQuinn M., Trac H., Cen R., Hernquist L. E., 2011, MNRAS, 414, 727