Journal of The Korean Astronomical Society. 2014. Apr, 47(2): 49-67
Copyright © 2014, null
  • Received : January 16, 2014
  • Accepted : February 25, 2014
  • Published : April 30, 2014
Export by style
Cited by
About the Authors
Department of Physics, Korea Institute for Advanced Study, Seoul 130-722,,
Department of Earth Science Education, Chosun University, Gwangju 501-759,
Department of Physics, Korea Institute for Advanced Study, Seoul 130-722,,
Center for Advanced Computation, Korea Institute for Advanced Study, Seoul 130-722,
Astronomy Centre, Department of Physics & Astronomy, Pevensy II bldg, University of Sussex, Falmer, Brighton BN1 9QH,
Department of Astronomy and Oskar Klein Centre, AlbaNova, Stockholm University, SE-106 91 Stockholm,

A novel method to characterize the topology of the early-universe intergalactic medium during the epoch of cosmic reionization is presented. The 21-cm radiation background from high redshift is analyzed through calculation of the 2-dimensional (2D) genus. The radiative transfer of hydrogen- ionizing photons and ionization-rate equations are calculated in a suite of numerical simulations under various input parameters. The 2D genus is calculated from the mock 21-cm images of high-redshift universe. We construct the 2D genus curve by varying the threshold differential brightness temperature, and compare this to the 2D genus curve of the underlying density field. We find that (1) the 2D genus curve reflects the evolutionary track of cosmic reionization and (2) the 2D genus curve can discriminate between certain reionization scenarios and thus indirectly probe the properties of radiation-sources. Choosing the right beam shape of a radio antenna is found crucial for this analysis. Square Kilometre Array (SKA) is found to be a suitable apparatus for this analysis in terms of sensitivity, even though some deterioration of the data for this purpose is unavoidable under the planned size of the antenna core.
The epoch of cosmic reionization (EOR) commences with the birth of first astrophysical, nonlinear objects such as the first stars and miniquasars. These sources of radiation create individual H II bubbles by emit- ting hydrogen-ionizing photons. This is then followed by the end of cosmic reionization when all percolating individual H II bubbles fully merge with one another and almost all hydrogen atoms in the universe become ionized.
Direct observations of EOR have not yet been made, even though several indirect observations imply that (1) the end of reionization occurred at a redshift z ≃ 6.5 (e.g., Fan et al. 2006), (2) the intergalactic medium remained at high temperature during and after EOR (Hui & Haiman 2003), and (3) this epoch started no later than z ≃ 11 (Komatsu et al. 2009). Theories suggest that the growth of H II bubbles proceeded inhomogeneously (in patchy way) due to the clustered distribution of radiation sources, and the global ionized fraction 〈 x 〉 increased monotonically in time (e.g., Iliev et al. 2007; McQuinn et al. 2007; Shin et al. 2008, a few models predict, however, non-monotonic increase in the global ionized fraction 〈 x 〉 due to the possible recombination after regulated star formation, followed by emergence of stars of higher tolerance to photoheating, e.g., Cen 2003; Wyithe & Loeb 2003).
Observation of the redshifted 21-cm line from neutral hydrogen is one of the most promising methods for the direct detection of EOR, and the Cosmic Dark Ages that precedes this epoch as well. Both temporal and spatial fluctuations in the 21-cm signal are believed to be strong in general, thus easy to detect, during EOR due to the inhomogeneous growth of H II bubbles and relatively weak foregrounds at high frequencies. The strongest signal will come from the late dark ages, or very early EOR, when the spin temperature will be much lower than the cosmic microwave background radiation (CMB) due to the Ly-α coupling to yet unheated intergalactic medium (IGM), even though stronger foregrounds at lower frequencies will be an obstacle to observation (e.g., Pritchard & Loeb 2008). There are several large radio interferometer arrays which aim to detect 21-cm signal from EOR. These projects include Giant Meterwave Radio Telescope (GMRT), Murchison Widefield Array (MWA), Precision Array for Probing the Epoch of Reionization (PAPER: Parsons et al. 2010), 21 Centimeter Array § (21CMA; formerly known as PaST), the LOw Frequency ARray (LOFAR), and Square Kilometre Array (SKA; for EOR detection strategy by SKA, see Mellema et al. 2013).
The lack of direct observations of EOR results in poor constraints on the history of cosmic reionization. Theoretical predictions for the history of cosmic reionization are made by semi-analytical calculations or fully numerical simulations. These studies select a variety of input parameters, among which the most important one is the properties of sources of radiation. The mock 21-cm data produced by such studies will be compared to future observations to constrain, for example, emissivity of high-redshift sources of radiation.
The patchy, 3-dimensional 21-cm radiation background can be analyzed in various ways such as the 3D power spectrum, 2D power spectrum, distribution of H II bubble size, cross-correlation of ionized fraction and overdensity, etc (e.g., Iliev et al. 2007; Zahn et al. 2006; Iliev et al. 2013). Each analysis method adds to the capability to discriminate between different reionization scenarios. They are usually complementary to each other, allowing one to understand the underlying physics better the more tools are available. Therefore, it is always favorable to have as many different tools for data analysis as possible.
In this paper, we characterize the geometrical property of the distribution of neutral IGMand H II bubbles by using the topology of 21-cm differential brightness temperature field. For this purpose we measure the 2D genus of a series of snapshots of the high-redshift universe that are predicted in different models of reionization. By varying the threshold value for the differential brightness temperature, we construct 2D genus curves at different redshifts under different reionization models. These models are simulated by self-consistent calculation of radiative transfer and rate equations over our simulation box.
Recently, similar methods for studying topology of high-redshift IGM have been suggested by Gleser et al. (2006), Lee et al. (2008), and Friedrich et al. (2011). They calculated the 3D genus using either the neutral (Gleser et al. 2006; Lee et al. 2008) or ionized (Friedrich et al. 2011) fraction of IGM. These studies show that the 3D genus of the underlying neutral (or ionized) fraction reflects the evolutionary stages of cosmic reionization. The 3D genus is also found to be useful in discriminating between different reionization scenarios. Albeit the similarity of our work to these studies, there are several fundamental differences. First, we calculate the 2D genus instead of the 3D genus. 2D sky analysis of a given redshift universe, which is possible in the 21cm observations, is beneficial especially when reionization proceeds rapidly, because then 3D analysis will suffer from the light-cone effect (e.g., Datta et al. 2012; La Plante et al. 2013) by mixing different redshift information along the line of sight. Second, we use an observable quantity, the 21-cm differential brightness temperature, which is more directly applicable to real data observed by future radio antennae, when calculating the 2D genus. Third, we explicitly seek for the possibility to apply our method to the data to be observed by SKA, with appropriate choice on the filtering scales for the 21-cmdifferential brightness temperature.
We employ the 2D genus method developed by Melott et al. (1989). The method of genus topology was originally designed to test mainly for the Gaussian random-phase nature of the primordial density field in 3D (Gott et al. 1986; Hamilton et al. 1986; Gott et al. 1987) or in 2D (Melott et al. 1989). The 2D case has been studied for a variety of cosmological data sets: on redshift slices (Park et al. 1992; Colley 1997; Hoyle et al. 2002), on sky maps (Gott et al. 1992; Park et al. 2001), on the CMB (Coles 1989; Gott et al. 1990; Smoot et al. 1992; Kogut 1993; Kogut et al. 1996; Colley et al. 1996; Park et al. 1998, 2001), and on the HI gas distribution in galaxy disks (Kim & Park 2007; Chepurnov et al. 2008).
Because 2D genus is strongly affected by the number and distribution of H II bubbles as well as the distribution of underlying density field, our method may provide an insight about how reionization proceeded and was affected by properties of sources. We also seek for the possibility to apply our method to observed data from SKA, which so far has the highest sensitivity proposed.
The plan of our paper is as follows. In Section 2 we describe our numerical simulations and the basic procedure to calculate the 21-cm radiation background. In Section 3 we describe how the 2D genus is calculated. In Section 4, our simulations results are analyzed and possibility of using SKA observation results are investigated. We summarize our result and discuss future prospects in Section 5.
- 2.1N-body Simulations
We completed a 2048 3 particle simulation in a concordance ACDM model with WMAP 5-year parameters (Dunkley et al. 2009); Ω m = 0.258, Ω A = 0.742, Ω b = 0.044, ns = 0.96, 𝜎 8 = 0.79, and h = 0.719, where Ω m , Ω A , Ω b are the density parameters due to matter, cosmological constant, and baryon, respectively. Here, ns is the slope of the Harrison-Zeldovich power spectrum, and 𝜎 8 is the root-mean-square (rms) fluctuation of the density field smoothed at 8 h −1 Mpc scale. The cubic box size is 66 h −1 Mpc in a side length. The simulation uses GOTPM, a hybrid PM+Tree N -body code (Dubinski et al. 2004; Kim et al. 2009). The initial perturbation is generated with a random Gaussian distribution on a 2048 3 mesh at zi = 500. The force resolution scale is set to f ε = 3.2 h −1 kpc in the comoving scale. The total of 12,000 snapshots, uniformly spaced in the scale factor, are created. The time step is predetermined so that the maximum particle displacement in each time step is less than the force resolution scale, f ε .
We then extract friend-of-friend (FoF) halos at 84 epochs: 24 epochs with five million year interval from z = 40 to 20, and 60 epochs with ten million year separation between z = 20 and z = 7. The connection length is set to 0.2 times the mean particle separation to identify cosmological halos. Using this method, we find all halos with mass above 10 8 M (corresponding to 30 particles). At each epoch, we calculate matter density fields on 2048 3 mesh using the Triangular-Shaped-Cloud (TSC) scheme (Efstathiou et al. 1985). As the radiative transfer through the IGM does not need to be run at too high resolution, we further bin down the original density fields to 256 3 after subtracting the contribution of collapsed dark halos. We transform this dark matter density to baryonic density by assuming that in the IGM baryons follow the dark matter with the mean cosmic abundance. From the list of collapsed halos, we form a “source catalogue” by recording the total mass of low-mass ( M < 10 9 M ) and high-mass ( M > 10 9 M ) halos in those radiative transfer cells containing halos.
We also ran a similar N-body simulation with the same cosmological parameters but a different configuration of 512 3 particles on a 512 3 mesh on a 64 h −1 Mpc box using GOTPM, resulting in somewhat poorer halomass resolution of 2 × 10 9 M . We created a binneddown density field on a 128 3 mesh, similarly, on which radiative transfer is calculated.
We use the results from the high resolution (2048 3 N -body) run for two cosmic reionization simulations (Case 3 and 4), and those from the low resolution (512 3 N -body) run for the other two (Case 1 and 2). For each cosmic reionization simulation, we choose a fixed time step ∆ t . These parameters are listed in Table 1 .
Simulation parameters. Cases 1 and 2 have only high-mass sources. Cases 3 and 4 have both high-mass and low-mass sources. Numbers in the “Model” arefγ’s for high-mass halos, and if any follows, arefγ’s for low-mass halos (with “S” representing “self-regulated”). Correspondinggγ’s are also listed. Radiative transfer is calculated on the number of meshes listed under “RT mesh.” Note that even though case 1 has the samegγas cases 3 and 4, its minimum massMmindoes not reach the low-end (109M⊙) of the other two, thus making its reionization history end much later.
PPT Slide
Lager Image
Simulation parameters. Cases 1 and 2 have only high-mass sources. Cases 3 and 4 have both high-mass and low-mass sources. Numbers in the “Model” are fγ’s for high-mass halos, and if any follows, are fγ’s for low-mass halos (with “S” representing “self-regulated”). Corresponding gγ’s are also listed. Radiative transfer is calculated on the number of meshes listed under “RT mesh.” Note that even though case 1 has the same gγ as cases 3 and 4, its minimum mass Mmin does not reach the low-end (109M) of the other two, thus making its reionization history end much later.
- 2.2 Reionization Simulations
We performed a suite of cosmic reionization simulations based upon the density field of IGM and the source catalogue compiled from the N-body simulations, described in Section 2.1. We used C 2 -Ray (Mellema et al. 2006a) to calculate the radiative transfer of H-ionizing photons from each source of radiation to all points in the simulation box and the change of ionization rate at each point, simultaneously. The spatial resolution of the binned-down density field is the actual radiative-transfer resolution, which is described in Table 1 for each case. At each N -body simulation step, we produce corresponding 3-dimensional map of ionized fractions.
History of cosmic reionization may depend on various physical properties of radiative source, such as the initial mass function (IMF) and inclusion of Population III stars (e.g., Ahn et al. 2012) or quasi-stellar objects (QSOs) as sources (e.g., Iliev et al. 2005; Mesinger et al. 2013). We vary physical properties of sources by changing the parameter f γ ≡ f esc f Ni over simulation, where f esc is the escape fraction of ionizing photons, f is the star formation efficiency, and Ni is the number of ionizing photons emitted per stellar baryon during its lifetime ∆ t (e.g., Iliev et al. 2007). f γ determines the source property in such a way that f γ M tot photons are emitted from each grid cell during ∆ t , where Mtot is the mass of halos inside the cell. ∆ t is also used as the time-step for finite-differencing radiative transfer and rate equations. If ∆ t varies among simulations, fixed f γ will generate different numbers of photons accumulated. Therefore, it is sometimes preferable to use a quantity somewhat blind to ∆ t . g γ ≡ f γ/(∆ t /10Myr), as defined in Friedrich et al. (2011), is such a quantity, which is basically the emissivity of halos.
We also adopt two different types of halos: low-mass (10 8 M/M ≳ 10 9 ) and high-mass ( M/M ≳ 10 9 ). Both types are “atomically-cooling” halos, where collisional excitation of atomic hydrogen Lyα line predominantly initiates cooling of baryons, reaching the atomic-cooling temperature limit of ∼ 8000K, and is followed subsequently by H 2 -cooling into protostellar cloud. One main difference exists in their vulnerability to external radiation: low-mass halos, when exposed to ionizing radiation, the Jeans mass after photoheating overshoots their virial mass such that the gas collapse is prohibited, while high-mass halos have mass large enough to be unaffected by such photoheating. We thus adopt a simple “self-regulation” criterion: when a grid cell obtains x > 0.1, we fully suppress star formation inside low-mass halos in the cell (e.g., Iliev et al. 2007). We also use in general different f γ’s for different types, when we include low-mass halos. Cases 1 and 2 did not allow us to implement low-mass halos due to the limited mass resolution, while cases 3 and 4 did, thanks to the high mass resolution.
Finally, we note that halos in even smaller mass range, 10 5 M/M ≲ 10 8 , or minihalos, are not considered in this paper. This can be justified because their impact on cosmic reionization, especially at late stage, are believed to be negligible (e.g., Haiman et al. 2000; Haiman & Bryan 2006). Nonetheless, the earliest stage of cosmic reionization must have started with minihalo sources, most probably regulated by inhomogeneous H 2 Lyman-Werner (LW) radiation background (Ahn et al. 2009; see also Dijkstra et al. 2008 for the impact of LW background on growth of supermassive black holes). A self-consistent simulation of cosmic reionization including minihalos and LW feedback has been achieved recently, showing that about ∼ 20% of global ionization of universe can be completed by minihalo sources at high redshifts (Ahn et al. 2012).
2.3 Calculation of 21-cm Differential Brightness Temperature
The hyperfine splitting of the ground state of hydrogen leads to a transition with excitation temperature T = 0.068K. The relative population of upper ( n 1 ) and lower ( n 0 ) states are determined by the spin temperature TS , such that
PPT Slide
Lager Image
Several pumping mechanisms determine TS as follows:
PPT Slide
Lager Image
where T CMB , T α , T K , y α , and y c are the cosmic microwave background (CMB) temperature, Lyα color temperature, the kinetic temperature of gas, the Lyα coupling coefficient, and the collisional coupling coefficient, respectively (Purcell & Field 1956; Field 1959).
21-cm radiation is observed in emission against CMB if TS > T CMB , or in absorption if TS < T CMB . This is quantified by the 21-cm differential brightness temperature (e.g., Morales & Hewitt 2004)
PPT Slide
Lager Image
if 21-cm lines from neutral hydrogen at redshift z is transferred toward us. The optical depth τ is given by (see e.g., Iliev et al. 2002 for details)
PPT Slide
Lager Image
δTb becomes appreciable enough to be detected by sensitive radio antennae only when TS TCMB or TS TCMB . The former limit is believed to be reached quite early in the history of cosmic reionization, because X-ray sources can easily heat up the gas over cosmological scales with small optical depth, such that TS TK TCMB (e.g., Ciardi & Madau 2003). There exists another limiting case where TS is determined only by collisional pumping ( y α ≪ 1), which is usu- ally the case in the Cosmic Dark Ages as studied by, for instance, Shapiro et al. (2006) and Kim & Pen (2010). The latter limit is reached when strong Lyα background is built up while the IGM is still colder than CMB temperature, such that TS TK TCMB . In this paper, we simply assume early heating of IGM such that TS TCMB . Eq. 3 and 4 then give the limiting form of δTb ,
PPT Slide
Lager Image
We simply generate 3D maps of δTb on the radiative transfer grid over all redshifts using Eq. 5. Note that in the limit TS T CMB , δTb ∝ (1 + δ )(1 − x ) at a given redshift. This limit, because δTb then becomes proportional to the underlying matter density, allows one to perform cosmology without any bias, by generating 3D maps of matter density distribution as long as x ≪ 1. Of course, patchy evolution of H II bubbles will be strongly imprinted on this δTb map as well, when the reionization process becomes mature (e.g., Shapiro et al. 2013).
In 3D, the genus of a single connected surface is identical to the number of handles it has. In 2D, if a field of a variable T is given, the 2D genus for a threshold value, Tth , is given by
PPT Slide
Lager Image
where N +( Tth ) and N _( Tth ) are the number of connected regions with T > Tth (hot spots) and T < Tth (cold spots), respectively. In general, as one varies Tth , g 2D changes. Therefore, g 2D is a function of Tth , and its functional form is called the 2D genus curve. In our case, we calculate the 2D genus curve of the field of 21-cm differential brightness temperature δTb , by varying the threshold value δTb , th.
The 2D genus curve of a field can work as a characteristic of the field. A useful template for this is a Gaussian random field, because its 2D genus is given analytically by the relation
PPT Slide
Lager Image
PPT Slide
Lager Image
is the deviation from the average 〈 T 〉 in units of standard deviation σT . One can then compare the 2D genus curve of a given field to g 2D,Gauss to see how close the field is to a Gaussian random field, for example.
The 2D genus is also equal to the integral of the curvature along the contours of T = Tth , divided by 2π. This is because when a curvature is integrated along a closed contour around a hot spot (cold spot), its value is 2π (−2π). When a 2D field is pixelized, a contour of T = Tth is composed of a series of line segments with turns occurring at vertices. On our rectangular 2D grid, every vertex shares four pixels (except for vertices on the edge and the corner). Under these conditions, the CONTOUR2D program is able to calculate the 2D genus by counting the turning of a contour observed at each vertex of four pixels in an image (Melott et al. 1989): 1/4 is contributed to the total genus from each vertex with 1 hot pixel and 3 cold pixels, −1/4 from each vertex with 3 hot pixel and 1 cold pixels, and zero otherwise. We do not calculate the modified 2D genus, g 2D,eff g 2D − 2 f , where f is the areal fraction of hot spots on the sky. g 2D,eff is appropriate when applied to a full-sky field (Gott et al. 1990; Colley & Gott 2003; Gott et al. 2007), while g2D is more appropriate for a relatively small, restricted region on the sky. Our simulation box is large enough to capture the typical size of H II bubbles (∼ 10 − 20Mpc comoving) during most phases of the cosmic reionization, but still too small to cover the whole sky and capture fluctuations in much larger scales.
Our aim is to obtain 2D genus curve from the frequency- and angle-averaged differential brightness temperature at given frequencies. In order to compare different reionization models, we will choose those frequencies such that different reionization simulations yield the same global ionized fraction. We will then vary the threshold δT b, th to construct 2D genus curve g 2D ( δT b, th ) for each case.
The amplitude of 2D genus is roughly proportional to the field of view at a fixed redshift. Throughout this paper, the 2D genus we present is normalized to the field of view of size (1 degree) 2 .
- 4.1 Mock 21-cm Sky Maps
Radio signals, observed by radio antennae, should be integrated over the angle, frequency and time in order to achieve the sensitivity required for specific scientific goals. Observing signals from high redshifts are challenging, especially, because these signals are very weak while the required angle- and frequency-resolutions are relatively high. Therefore, very intensive integration in time (∼ 1000 hours) is required for EOR observations in general.
We therefore choose the beam size and the frequency bandwidth which we find adequate for generating distinctive 2D genus curves from different reionization scenarios, under the assumption of very long time integration (∼ 1000 hours or more). We generated mock 21-cm emission sky maps for our four reionization scenarios, f5, f250, f125_125S and f125_1000S, averaging δTb on our computational grid over the angle and the frequency. First, note that our box size corresponds to the frequency bandwidth
PPT Slide
Lager Image
and the transverse angle
PPT Slide
Lager Image
where 𝜈 0 ≡ 1420MHz is the rest frame frequency of the line, DA ( z ) is the angular diameter distance, L box is the comoving length of our simulation box, c is the speed of light and H 0 ≡ 100 h kms −1 Mpc −1 is the Hubble parameter at present. The numerator in the last parenthesis in Eq. 10 is a fitting formula for the line-of-sight comoving distance r ( z ) in units of Mpc, under the cosmological parameters used in this paper.
We first calculate bare differential brightness temperature on our computational grid by using Eq. 5. We then integrated them over frequency with chosen ∆ν’s. We also consider the Doppler shift due to peculiar velocity. Along the line-of-sight (LOS), the frequency-averaged δTb then becomes
PPT Slide
Lager Image
where N ∆ν ≡ ∆ L (∆ ν ) / l cell is the number of cells corresponding to ∆ ν at a given redshift, fi is the fraction of a cell i entering the band after the Doppler shift,
PPT Slide
Lager Image
is the comoving length along LOS corresponding to ∆ ν , and l cell is the comoving length of a unit cell. For a statistical purpose, we generate four 21-cm emission maps for each LOS direction from a single simulation box by periodically shifting the box by a quarter of box size ( L box /4). These maps will be almostmutually independent, because L box /4 ≫ ∆ L (∆ ν ). We also take three different directions for LOS ( x -, y -, and z -direction) by rotating the box, such that we have in total 12 (almost) independent 21-cm sky maps from a single simulation box at a given redshift.
Because the shape of the beam varies over different antennae, we choose two different types of beams to average the signal in angle: Gaussian and compensated Gaussian. We first convolve the 21-cm emission map from simulations with a Gaussian filter
PPT Slide
Lager Image
with a full width at half maximum (FWHM) ∆ θ =
PPT Slide
Lager Image
. We also use a compensated Gaussian filter given by
PPT Slide
Lager Image
with a FWHM given by
PPT Slide
Lager Image
where LambertW(e/2) ≈ 0.685.
A Gaussian filter is widely used in literature to approximate actual beams. In general, an interferometer will not retain the average (large scale) signal, but only record the fluctuations between minimum and maximum angular scales (depending on the largest and smallest baselines). As the largest scales to which the upcoming interferometers are sensitive are larger than our simulated images, a reasonable approach would be to subtract the mean signal from our images. However, for the determination of the genus this does not make a difference other than shifting the threshold value, so we retain the average signal here. A compensated Gaussian filter roughly mimics the beam of a compact interferometer (e.g., Mellema et al. 2006b), even though real beams have sometimes much more complicated shape depending on the actual configuration of the antenna. We just take W CG as an extreme case for “dirty” beams, clearly distinct from W G , because this filter is somewhat pathological as it suppresses both small and large scale features. When
PPT Slide
Lager Image
, W CG ( θ ) < 0 such that a data field convolved with W CG ( θ ) may change its sign from point to point. 21-cm signals filtered this way would show seemingly unnatural H II bubble feature, or in a worse case, even absorption signals ( δTb < 0), even though Eq. 5 implies that δTb > 0 when TS T CMB .
Fig. 1 shows how the actual 2D δTb field (top rows) will be observed in different filtering schemes (middle: W G ; bottom: W CG ). Three different epochs were chosen to represent early (volume weighted ionized fraction xv = 0.04), middle ( xv = 0.4), and late ( xv = 0.9) stages. W C generates reasonably filtered δTb maps at all three epochs. At very early ( xv = 0.04) and late ( xv = 0.9) stages, fluctuation in δTb is dominated by 1 + δ and 1 − x , respectively, while at the middle ( xv = 0.4), both 1 + δ and 1 − x contribute to the fluctuation.
PPT Slide
Lager Image
— 21-cm maps at redshifts z = 15.645, 12.184 and 11.325 for f125_1000S. Top row: Slices which are frequencyaveraged into one band of ∆𝜈 = 0.2MHz, showing ionized regions (yellow) and neutral regions (blue), superimposed on the density field (light and dark for high-density and low-density regions, respectively). Middle row: Corresponding 2D 21-cm differential brightness temperature, averaged by a Gaussian beam of 𝜈𝜃 = 1′. Bottom row: The same as middle, but with a compensated Gaussian beam of the same 𝜈𝜃.
The compensated Gaussian beam adds ripple structures to the emission signals. When a uniform field confined within a certain boundary is convolved with W CG , this ripple becomes visible only along the boundary. As the ionized regions grow, therefore, a compensated Gaussian beam produces negative through of width similar to the beam FWHM at the boundaries of bubbles. In general, W CG changes the topology of the 21-cm signal. For most of our analysis, therefore, we just use W G . Further comparison of these two filters will be made in Section 4.3.2.
- 4.2 Evolution, Sensitivity, Bubble Distribution, and Power Spectrum
The evolution of global volume-weighted ionized fraction xv varies significantly over different reionization scenarios, as seen in Fig. 2 . First, the end of reionization in f5 and f250 cases occurs relatively later than f125_125S and f125_1000S cases, because these are tied only to high-mass ( M > 2 × 10 9 M ) halos, which collapse much later than low-mass (10 8 M/M ≲ 109) halos. Second, f5 and f250 also show steeper evolution in xv , because there is no chance for slow evolution of total luminosity from self-regulation of small-mass halos.
PPT Slide
Lager Image
— Left: Evolution of the mean volume ionized fraction as a function of redshift. Middle: The mean differential brightness temperature with 1′ FWHM Gaussian beam and 0.2MHz bandwidth as a function of the mean volume ionized fraction. Right: rms fluctuation of the differential brightness temperatures filtered by the same beam and bandwidth as a function of the mean volume ionized fraction.
It is in general a better practice to compare different models at a fixed global ionized fraction than at a fixed redshift, partly because there are still too much freedom in the exact epoch of the end of reionization, the time-integrated optical depth to Thompson scattering of CMB photons, etc. More importantly, a fixed global ionized fraction among different models is reached by roughly the same number of H-ionizing photons accumulated from the beginning of EOR. We therefore use xv as the time indicator of the global evolution throughout this paper.
Using xv as the time indicator, we first observe similar trends in the evolution of δTb among different models. The mean differential brightness temperature 〈 δTb 〉 starts from 〈 δTb 〉 ≈ 30 − 40mK and gradually de- creases in time, as seen in Fig. 2 . The root-mean-square (rms) fluctuation of δTb , δTb,rms , reaches the maximum (∼ 8−11mK) at xv ≈ 0.4−0.6 ( Fig. 2 ), when the signal is filtered with ∆ ν = 0.2MHz and ∆ θ = 1′.
Albeit the similar trends, detailed analysis will be needed to discriminate between different models and ultimately probe the properties of sources. The probability distribution function (pdf) of bubble size is a useful tool to understand the impact of these source properties. The bubble size is associated with the mass spectrum of halos and the ionizing efficiency f γ, simply because it roughly reflects the total number of Hionizing photons emitted into the bubble. In order to find the pdf of bubble size, we use a hybrid method. First, the size of a bubble is determined by the method in Zahn et al. (2006): the bubble size is the maximum radius of a sphere from each simulation cell inside which the ionized fraction is over 90%. This way, every cell in the box is associated with a bubble with a certain size. We then use the void-finding method by Hoyle & Vogeley (2002), which has been extensively used for finding cosmological voids which are mutually exclusive in space. The bubbles are sorted in size from the largest to the smallest. The largest bubble is considered an isolated one. We then move to the next largest bubble and check if there is any overlap in volume with the fist one. We iterate this over the bubble list until the overlap becomes less than 10%, which then registers another isolated bubble. This process is further iterated and we obtain the full list of mutually exclusive bubbles (to the extent of 10% overlap).
In Fig. 3 , we show the pdf of bubble size obtained this way. We can observe clear distinction of f5 and f250 from f125_125S and f125_1000S. In f5 and f250, a given ionized fraction is reached only by high-mass halos. For example, to reach xv ∼ 0.6, the evolutionary stage of cosmological structure formation should enter a highly nonlinear phase because there is only one species (highmass halos) responsible for such high xv . In contrast, in f125_125S and f125_1000S, both of the two different species (high- and low-mass halos) contribute to, for example, xv ∼ 0.6. Therefore, the structure formation will be in less nonlinear stage than the former cases. Correspondingly, we can expect stronger clustering of sources for f5 and f250, and weaker for f125_125S and f125_1000S. Correspondingly, merger of bubbles would be stronger for the former and weaker for the latter, respectively. This is finally reflected in the relative contribution to a given xv from largest bubbles, as seen in Fig. 3 . Large bubbles dominate xv in f5 and f250, while there are relatively smaller contribution to xv in f125_125S and f125_1000S.
PPT Slide
Lager Image
— Left: The probability distribution functions of bubble size. We determined the size of a bubble using the method in Zahn et al. (2006), and select mutually exclusive bubbles using the method in Hoyle & Vogeley (2002). Right: The volume-weighted probability distribution functions of bubble size.
We also obtain 2D angular power spectra of the neutral fraction, 1 − x , for further analysis. This also provides an insight on the clustering scale of sources of radiation, or H II bubbles, through the location of peaks in the power spectrum. We follow the convention used in CMB analysis and constructed [ℓ(ℓ+1) C ℓ/2π] 1/2 in spherical harmonics, where ℓ roughly corresponds to 2π/ θ . See Fig. 4 . At xv = 0.04, it is clearly seen how cosmic reionization commences in each scenario. f5 and f250 show peaks in ℓ(ℓ+1)Cℓ at ℓ ∼ 10 4 , or the comoving length scale x ∼ 6 − 7 h −1 Mpc, while f125_1000S shows a peak at ℓ ∼ 3 × 10 4 , or x ∼ 2 h −1 Mpc, and f125_125S shows a long uniform tail from ℓ ∼ 3×10 4 to ℓ ∼ 4 × 10 4 . The behavior of 2D angular power spectra indicates stronger merger of (otherwise) individual bubbles into larger scale at x ∼ 6 − 7 h −1 Mpc in f5, f250 and f125_125S than in f125_1000S. Even though f125_125S has a tail up to ℓ ∼ 4×10 4 (x ∼ 1 h −1 Mpc), the smaller f γ(= 125) in small halos than that of f125_1000S ( f γ = 1000) makes these small halos less efficient. This is true at all redshifts: the power spectrum of f125_125S is hardly distinguishable from those of f5 and f250 any time except for higher power in the smallest length scale (highest ℓ). In contrast, a distinctive clustering scale, ℓ ∼ 3 × 10 4 (x ∼ 2 h −1 Mpc), is shown at all redshifts in f125_1000S, due to the clustering of small halos of the highest efficiency.
PPT Slide
Lager Image
— 2D angular power spectra of neutral fraction 1 − x.
In short, both bubble pdf and 2D angular power spectrum of neutral fraction may work as useful tools for understanding the nature of different reionization scenarios. In Section 4.3, when we analyze 2D genus curves of different cases, we will address these properties to see whether 2D genus curves created from possible observations can bear similar fundamental nature as well.
- 4.3 Genus Properties
We generate g 2D ( δTb ) from selected epochs at which xv = 0.04, 0.4, 0.6 and 0.9 as described in Section 3. The base 2D field of δTb is averaged in frequency with ∆ ν = 0.2, 1 and 2MHz, and in angle with ∆ θ = 1′, 2′ and 3′. We vary ∆ ν and ∆ θ to quantify the competing effects of changing resolution and sensitivity, and also to make comparison with characteristics of future radio antennae. In angle-averaging, only W G is used in all cases, except for f125_1000S to which W CG is also applied to understand the impact of beam shape on g 2D .
- 4.3.1 2D Genus of Density Fluctuation
Is there any useful template genus curve to be compared with g 2D ( δTb )? There is such one indeed, which is g 2D of an artificial field of δTb which reflects the density fluctuation only. Let us denote this quantity by
PPT Slide
Lager Image
Distribution of δTb,n will be Gaussian in the linear regime (| δ | ≪ 1), which will provide well-defined, analytical g 2D,Gauss ( δTb,n ). Even when g 2D ( δTb,n ) starts to deviate from g 2D ,Gauss( δTb,n ) due to nonlinear evolution of high density peaks, g 2D ( δTb,n ) will be used to indicate overdense ( δ > 0) and underdense ( δ < 0) regions. Fig. 5 implies that δTb,n > 0 corresponds to δ > 0, while δTb,n < 0 to δ < 0.
g 2D ( δTb,n ), plotted in Fig. 5 , show how matter density evolves. At z ≈ 10 − 15, there already exist extended tails into high δTb,n , because these correspond to high-density peaks that have evolved into nonlinear regime. Note that purely Gaussian random field generates g 2D,Gauss , which is symmetrical around ν = 0, while high-density peaks (on high ν ) that evolved non- linearly deviates from the Gaussian distribution to also make g 2D deviate from g 2D,Gauss . Because we will use filtered δTb maps , g 2D of filtered δTb,n will be the template to be used throughout this paper. Filtering makes both the amplitude and width of g 2D ( δTb,n ) shrink ( Fig. 5 ). Nevertheless, even after filtering, regions with positive g 2D correspond to overdense regions and negative g 2D to underdense regions, because the meandensity point (where the g 2D ( δTb,n ) curve crosses the δTb axis) is almost unshifted. This will work as a perfect indicator, if comparison is made with g 2D ( δTb ), about how different the δTb field is from the underlying density field. g 2D ( δTb,n ) can be obtained accurately for any redshift as long as cosmological parameters are known to high accuracy, which is the case in this era of precision cosmology.
PPT Slide
Lager Image
— 2D genus defined by the threshold differential brightness temperature within ∼ 30′ × 30′ square (corresponding to our box size) in the sky assuming all the gas was neutral, filtered with 1′ FWHM Gaussian beam and 0.2MHz bandwidth (red). In comparison, we also show the 2D genus of the differential brightness temperature with the same condition but from one-cell thick slice in our simulation (blue). g2D(𝛿Tb,n) of the bare N-body field is thus subject to smoothing over the size of the cell. This result is based upon the high-resolution N-body simulation (used for cases 3 and 4).
- 4.3.2 Evolution of g2Dand Impact of Beam Shape
We found that the evolution of g 2D ( δTb ) contained both generic and model-dependent features. We describe the generic features based on f125_1000S case. We will describe the model-dependent features in Section 4.3.3.
We find that the change of g 2D ( δTb ) clearly shows how reionization proceeds in time. The left panel of Fig. 6 shows W G -filtered g 2D ( δTb,n ) of the underlying matter density field (blue dashed line) and g 2D ( δTb ) of the differential brightness temperature field of HI gas (red line). In the earliest phase of reionization ( xv = 0.04), g 2D ( δTb ) is much smaller than g 2D ( δTb,n ) at high temperature thresholds ( δTb,n ≳ 50mK). It is because the highest density peaks have been ionized by the sources inside them, and are dropped out of the neutral gas distribution. It can be also noted that the amplitudes of both the maximum (at δTb ∼ 38mK) and minimum (at δTb ∼ 33mK) of g 2D ( δTb ) are higher than those of 2D ( δTb,n ). Birth of new islands and peninsulas made of H I regions — some peninsulas may appear as islands at certain δT b,th — is responsible for the former, while birth of new lakes made of H II regions is responsible for the latter. Birth of islands and peninsulas occurs at mildly overdense regions ( δTb,n ∼ 38mK) because bubbles are clustered such that some of them merge with one another to fully or partly surround neutral regions. Note that a neutral region identified as an island in 2D may well be a cross-section of a vast neutral region in 3D.
When the universe enters the middle phase of reionization ( xv ≈ 0.4 − 0.6), almost all overdense regions ( δTb,n ≳ 30mK) have been ionized to form larger bubbles. This is reflected in the relatively low amplitude trough of g2D(δTb), which appears in very low-density ( δTb,n ≳ 20mK) regions. Now many more islands appear in mildly underdense ( δTb,n ∼ 25mK) regions because larger and more clustered bubbles penetrate further into low-density IGM and are more efficient in forming new islands. As time passes from xv = 0.4 to xv = 0.6, the amplitude of g 2D decreases as bubbles merge with each other and neutral clumps disappear. Further penetration of bubbles into lower-density IGM from xv = 0.4 to xv = 0.6 is also reflected in the maximum value of δTb for nonzero g 2D ( δTb ).
Finally, in the late stage ( xv = 0.9), all ionized bubbles have been connected and the last surviving neutral clumps exist to give positive g 2D at δTb ∼ 20mK. Note that these clumps exist only in underdense regions, which is a clear indication of the fact that reionization proceeds in an inside-out fashion: high-density regions will be ionized first due to the proximity to sources of radiation, and low-density regions will be ionized later.
Now, we briefly investigate the impact of beam shapes. The impact of W CG on g 2D is demonstrated on the right panel of Fig. 6 . Note that the Fourier transform of W CG is proportional to k 2 σ 2 exp(− k 2 σ 2 /2). From the convolution theorem, it is obvious that the fil- tered field becomes the WG-smoothed, negative Lapla- cian (−∇ 2 ) of the field. Note that Laplacian is equiv- alent to the divergence of the gradient. Therefore, if the δTb field is filtered by W CG , δTb s in those regions with the steepest spatial gradient in x or δ will correspond to extrema. We also note that the sign of the Laplacian depends on the morphology of x : if H II bubbles form in the sea of neutral gas, gradients of δTb diverge (−∇ 2 δTb < 0), while if neutral clumps remain in the sea of ionized gas, gradients of δTb converge (−∇ 2 δTb > 0). An interesting feature in g 2D ( δTb ) is indeed observed to move from regions with δTb < 0 at the relatively early ( xv ∼ 0.4) stage to regions with δTb > 0 at the late ( xv ∼ 0.9) stage. Nevertheless, because the topology of the 21-cm signal is processed further (Gaussian smoothing and Laplacian) in the case of W CG than that of W G (Gaussian smoothing only), the interpretation becomes less transparent. We also note that W CG filters out almost all scales except for the filtering scale, while W G filters out only those scales smaller than the filtering scale.
PPT Slide
Lager Image
— Left: 2D genus of 𝛿Tb for f125_1000S, filtered with 1′ FWHM Gaussian beam and 0.2MHz bandwidth (solid), compared to the template 2D genus out of the underlying density field (dashed), also filtered the same way. Right: Same as left, but with a compensated Gaussian beam applied.
The nice interpretative power of g 2D ( δTb ) is somewhat lost with W CG . Based on this, we suggest that any artificial effects from “dirty” beams should be removed or minimized in real observations. Only when the effective beam shape becomes something close to WG, general analyses including our genus method will obtain their full potential.
- 4.3.3 Model-dependent Feature and Required Sensitivity
We finally investigate whether our 2D genus analysis can discriminate between different reionization scenarios. Ultimately, if this turns out true, we may be able to probe properties of radiation sources at least indirectly, because these properties determine how cosmic reionization proceeds.
Our 2D genus analysis seems very promising, indeed, in probing source properties. The top panel in Fig. 7 shows g 2D ( δTb ) and g 2D ( δTb,n ) of all cases at xv = 0.4. The impact of small-mass (10 8 M/M ≲ 10 9 ) halos is observed as follows. In f125_1000S, almost all (filtered) overdense regions have been ionized at this epoch due to efficient, small-mass halos. In contrast, when there are no small-mass halos (f5, f250) or if small-mass halos are not as efficient (f125 125S), some fraction of overdense regions still remains neutral at xv ∼ 0.4. These regions are mildly nonlinear, which are ionized due to the almost on-the-spot existence of efficient small-mass halos in f125_1000S, while they mostly remain neutral or only partially ionized in f5, f250 and f125_125S due to absence or low-efficiency of small-mass halos ( Fig. 10 ). The relatively larger amplitude of g 2D in f125_1000S is the reflection of the fact that there are many more small individual islands and bubbles, as explained in Section 4.3.2.
PPT Slide
Lager Image
— 2D genus curves of the differential brightness temperatures at xv = 0.4, filtered with 1′ FWHM Gaussian beam and with 0.2MHz (top), 1MHz (middle), 2MHz (bottom) bandwidth, compared to those of the matter densities (dashed).
The impact of beam size and bandwidth is shown in Figs. 7 9 . As the applied ∆ ν increases, g 2D curve shrinks in width and moves toward left. The interpretation we made on f125_1000S with ∆ ν = 0.2MHz and ∆ θ = 1′ somewhat weakens in the sense that g 2D curves are mapped too deep in (newly filtered) underdense regions. Nevertheless, the amplitude of g 2D is the largest and the rightmost wing of g 2D is located at the smallest δTb in f125_1000S in all varying filtering scales, just as when ∆ ν = 0.2′MHz is applied.
PPT Slide
Lager Image
— Same as Fig. 7, but with 2′ FWHM Gaussian beam applied.
PPT Slide
Lager Image
— Same as Fig. 7, but with 3′ FWHM Gaussian beam applied.
PPT Slide
Lager Image
— How 2D genus depends on reionization scenarios: from left to right, threshold 𝛿Tb’s (contour lines) are selected which give the minimum g2D(δTb), maximum g2D(δTb), and zero g2D(δTb,n), respectively, together with the snapshot of unfiltered ionized fraction (far right). xv = 0.4 in both (top: f125_1000S; bottom: f125_125S) cases. See arrows in Fig. 7 for contours a, b and c.
Since g 2D is found a useful tool for understanding the evolution (Section 4.3.2) and even discriminating between reionization scenarios, we need to ask whether this can be achieved in real observations. To estimate the required sensitivity, we show the degree of fluctuation in filtered δTb in Fig. 11 . δTb , rms becomes maximum around the middle stage of reionization in all cases. Roughly speaking, δTb , rms gives the required rms sensitivity limit of a radio antenna in each configuration, which is about a few mK.
Following Iliev et al. (2003), we estimate the sen- sitivity limit of SKA for the filtering scale treated in this paper, to figure out feasibility of our analysis on the real data of 21-cm background. We assume the core size to be ∼ 1 km for SKA. We also assume, just for the sake of sensitivity estimation, that xv ∼ 0.5 is reached at z ∼ 8. Even though our choice of input parameters (properties of sources) made reionization end much earlier than the usually believed epoch of end of reionization ( z ∼ 6.5), except for f 5, we may imagine that similar distinctive features among different models would still exist if we tuned these parameters to make these models achieve the half-ionized state at z ∼ 8. Fig. 11 shows the estimated sensitivity under these conditions. Due to the strong dependence of sensitivity on θ , ∼ ∆ θ −2 , increase in the beam size quickly makes our analysis feasible. This comes at a price that values of g 2D ( δTb ) shrinks further and the range of δTb shifts more to the left. At any rate, ∼ 1000 − 10000 hours of integration at δν ∼ 1−2 MHz and δθ ∼ 2′−3′ seems possible to allow our analysis.
PPT Slide
Lager Image
— Top row: Sensitivity limit of SKA in terms of bandwidth ∆𝜈, with varying integration time (100, 1000 and 10000 hours). Also plotted are the required sensitivity limits (cross) for useful results, at ∆𝜈 = 0.2, 1 and 2 MHz. Each column corresponds to varying beam size ∆𝜃’s (top labels). The rest: Evolution of rms differential brightness temperature for varying ∆𝜃 (each column) and ∆𝜈 (each row).
We calculated the 21-cm radiation background from high redshift using a suite of structure formation and radiative transfer simulations with varying properties of sources of radiation. These properties of sources are specified by the spectrum of halo masses capable of hosting sources of radiation and the emissivity of hydrogen-ionizing photons. Assuming the prereionization heating of IGM, we calculated the 21-cm radiation background in the saturated limit, TS TCMB , such that its differential brightness temperature is proportional to the underlying density and neutral fraction.
In order to understand the topology of the highredshift IGM, we developed a method of 2D genus topology applicable to the 21-cm radiation background we calculated. Basically, this method calculates the 2D genus, the difference between the number of hot spots and cold spots, under a given threshold differential brightness temperature. We construct 2D genus curves at different redshifts for different scenarios, by varying this threshold values.
We found that the 2D genus curve g2D( δTb ), if compared to g 2D of the underlying density fluctuation, g 2D( δTb,n ), clearly shows the evolution of the reionization process. g2D( δTb ) is found to deviate quickly from the g 2D of a Gaussian random field, qualitatively in accord with the findings by Iliev et al. (2006) and in contradiction to the findings by Shin et al. (2008). Shin et al. (2008) find that the non-Gaussianity of the H II region is small even when the universe is 50% ionized. It is also shown that the reionization proceeds in an inside-out fashion: high-density regions gets ionized first and low-density regions gets ionized later.
We also showed that our 2D genus method can be used to discriminate between various reionization scenarios, thus probing properties of sources indirectly. It seems most effective in discriminating the mass spectrums of halos which host sources of radiation. A hybrid mass spectrum with different emissivity in different species (f125_1000S in this paper) stands out from cases with a single species (f5, f250) or a case with equal-emissivity among different species (f125_125S).
Crucial ingredients needed for this analysis is the beam shape and sensitivity. We tested two different beams, Gaussian and compensated Gaussian. The Gaussian beam, even after the field is filtered and degraded, leaves discernible imprint on the 2D genus curve such that different reionization scenarios can be discriminated. In contrast, the compensated Gaussian beam somewhat loses the interpretative power that the Gaussian beam had. Assuming that the compensated Gaussian beam represents the so-called dirty beams, our analysis may be applied to full potential only when those artificial effects from dirty beams are removed.
We predict that SKA will be able to produce data suitable for this analysis, when 1000 − 10000 hours of integration is performed with δν ∼ 2 − 3MHz and δθ ∼ 2′ − 3′ at observing frequency of ∼ 150MHz. Therefore, our method seems promising. Note, however, that a direct link from the observed 2D genus curve to true properties of sources will be still farfetched, when we adopt such low-resolution filters to increase sensitivity. Moreover, there still exist too many uncertainties, such as the matter power spectrum in small scales, which is relevant to formation and evolution of sources of reionization. We hope to see many more useful constraints on cosmic reionization come from various other, direct or indirect, observations. In the future, we will explore more reionization scenarios to strengthen the potential of our 2D genus analysis.
One caveat of our work is that the simulation box is relatively small, while ≳ 100 Mpc box is favored to provide a more statistically reliable result. While 2D genus should be dominated by the most abundant H II bubbles, which are also the smallest bubbles, the cosmic variance guarantees that some of the bubble size, at late stages of reionization, will become larger than the box we used in this work. We have indeed simulated reionization in a very large box (425/ h Mpc) and observed variation of observational properties (e.g., 3D power spectrum) in certain range of scales and reionization stages (Iliev et al. 2013). We did not yet analyze the minihalo-included simulation (Ahn et al. 2012), which forms another important class of models with very small H II regions, with the 2D genus either. We thus plan to apply our 2D genus analysis to these new results (more model dependency and bigger box size) and provide a more statistically reliable result which can be used to analyze future 21-cm observations.
KA was supported in part by NRF grant funded by the Korean governmentMEST (No. 2012R1A1A1014646), and KA also acknowleges the generous support from Chosun University for KA’s sabbatical (research) year. We thank Korea Institute for Advanced Study for providing computing resources (KIAS Center for Advanced Computation Linux Cluster System) for this work. ITI was supported by the Science and Technology Facilities Council [grant numbers ST/F002858/1 and ST/I000976/1]; and The Southeast Physics Network (SEPNet). GM was supported by Swedish Research Council grant 2012-4144.
Adler R. J. 1981 The Geometry of Random Fields Wiley New York
Ahn K. , Shapiro P. R. , Iliev I. T. , Mellema G. , Pen U. L. 2009 The Inhomogeneous Background of H2 Dissociating Radiation During Cosmic Reionization ApJ 695 1430 -    DOI : 10.1088/0004-637X/695/2/1430
Ahn K. , Iliev I. T. , Shapiro P. R. , Mellema G. , Koda J. , Mao Y. 2012 Detecting the Rise and Fall of the First Stars by Their Impact on Cosmic Reionization ApJ 756 L16 -    DOI : 10.1088/2041-8205/756/1/L16
Cen R. 2003 Implications of WMAP Observations On the Population III Star Formation Processes ApJ 591 L5 -    DOI : 10.1086/377068
Chepurnov A. , Gordon J. , Lazarian A. , Stanimirovic S. 2008 Topology of Neutral Hydrogen within the Small Magellanic Cloud ApJ 688 1021 -    DOI : 10.1086/591655
Ciardi B. , Madau P. 2003 Probing Beyond the Epoch of Hydrogen Reionization with 21 Centimeter Radiation ApJ 596 1 -    DOI : 10.1086/377634
Coles P. 1989 The Clustering of LocalMaxima in Random Noise MNRAS 238 319 -    DOI : 10.1093/mnras/238.2.319
Colley W. N. 1997 Two Dimensional Topology of Large Scale Structure in the Las Campanas Redshift Survey ApJ 489 471 -    DOI : 10.1086/304799
Colley W. N. , Gott J. R. , Park C. 1996 Topology of COBE Microwave Background Fluctuations MNRAS 281 L82 -    DOI : 10.1093/mnras/281.4.L82
Colley W. N. , Gott J. R. I. 2003 Genus Topology of the Cosmic Microwave Background from WMAP MNRAS 344 686 -    DOI : 10.1046/j.1365-8711.2003.06907.x
Datta K. K. , Mellema G. , Mao Y. 2012 Light-Cone Effect on the Reionization 21-cm Power Spectrum MNRAS 424 1877 -    DOI : 10.1111/j.1365-2966.2012.21293.x
Dijkstra M. , Haiman Z. , Mesinger A. , Wyithe S. 2008 Fluctuations in the High-Redshift Lyman- Werner Background: Close Halo Pairs as the Origin of Supermassive Black Holes MNRAS 391 1961 -    DOI : 10.1111/j.1365-2966.2008.14031.x
Dubinski J. , Kim. J. , Park C. , Humble R. 2004 GOTPM: A Parallel Hybrid Particle-Mesh Treecode New Astron 9 111 -    DOI : 10.1016/j.newast.2003.08.002
Dunkley J. , Spergel D. N. , Komatsu E. , Hinshaw G. , Larson D. , Nolta M. R. , Odegard N. , Page L. , Bennett C. L. , Gold B. , Hill R. S. , Jarosik N. , Weiland J. L. , Halpern M. , Kogut A. , Limon M. , Meyer S. S. , Tucker G. S. , Wollack E. , Wright E. L. 2009 Five-Year Wilkinson Microwave Anisotropy Probe (WMAP) Observations: Likelihoods and Parameters from theWMAP data ApJS 180 306 -    DOI : 10.1088/0067-0049/180/2/306
Efstathiou G. , Davis M. , Frenk C. S. , White S. D. M. 1985 Numerical Techniques For Large Cos- mological N-Body Simulations ApJS 57 241 -    DOI : 10.1086/191003
Fan X. , Strauss M. A. , Becker R. H. , White R. L. , Gunn J. E. , Knapp G. R. , Richards G. T. , Schneider D. P. , Brinkmann J. , Fukugita M. 2006 Constraining the Evolution of the Ionizing Back- ground and the Epoch of Reionization with z - 6 Quasars II: A Sample of 19 Quasars AJ 132 117 -    DOI : 10.1086/504836
Field G. B. 1959 The Spin Temperature of Intergalactic Neutral Hydrogen ApJ 129 536 -    DOI : 10.1086/146653
Friedrich M.M. , Mellema G. , Alvarez M. A. , Shapiro P. R. , Iliev I. T. 2011 Topology and Sizes of H II Regions during Cosmic Reionization MNRAS 413 1353 -    DOI : 10.1111/j.1365-2966.2011.18219.x
Gleser L. , Nusser A. , Benedetta C. , Desjacques V. 2006 The Morphology of Cosmological Reionization by means of Minkowski Functionals MNRAS 370 1329 -    DOI : 10.1111/j.1365-2966.2006.10556.x
Gott J. R. , Park C. , Juszkiewicz R. , Bies W. E. , Bennett D. P. , Bouchet F. R. , Stebbins A. 1990 Topology of Microwave Background Fluctuations: Theory ApJ 352 1 -    DOI : 10.1086/168511
Gott. J. R. , Mao. S. , Park C. , Lahav. O. 1992 The Topology of Large-Scale Structure. V. - Two Dimensional Topology of Sky Map ApJ 385 26 -    DOI : 10.1086/170912
Gott J. R. , Colley W. N. , Park C. G. , Park C. , Mugnolo C. 2007 Genus Topology of the Cosmic Microwave Background from the WMAP 3-Year Data MNRAS 377 1668 -    DOI : 10.1111/j.1365-2966.2007.11730.x
Gott J. R. , Choi Y.-Y. , Park C. , Kim J. 2009 3D Genus Topology of Luminous Red Galaxies ApJ 695 L45 -    DOI : 10.1088/0004-637X/695/1/L45
Gott. J. R. , Melott A. L. , Dickinson M. 1986 The Sponge-Like Topology of Large-Scale Structure in the Universe ApJ 306 341 -    DOI : 10.1086/164347
Gott J. R. , Weinberg D. H. , Melott A. L. 1987 A Quantitative Approach to the Topology of Large- Scale Structure ApJ 319 1 -    DOI : 10.1086/165427
Haiman Z. , Abel T. , Rees M. J. 2000 The Radiative Feedback of the First Cosmological Objects ApJ 534 11 -    DOI : 10.1086/308723
Haiman Z. , Bryan G. L. 2006 Was Star-Formation Suppressed in High-Redshift Minihalos? ApJ 650 7 -    DOI : 10.1086/506580
Hamilton A. J. S. , Gott J. R. , Weinberg D. H. 1986 The Topology of the Large-Scale Structure of the Universe ApJ 309 1 -    DOI : 10.1086/164571
Hoyle F. , Vogeley M. S. 2002 Voids in the PSCz Survey and the Updated Zwicky Catalog ApJ 566 641 -    DOI : 10.1086/338340
Hoyle F. , Vogeley. M. S. , Gott J. R. 2002 Two- Dimensional Topology of the 2dF Galaxy Redshift Survey ApJ 570 44 -    DOI : 10.1086/339582
Hui L. , Haiman Z. 2003 The Thermal Memory of Reionization History ApJ 596 9 -    DOI : 10.1086/377229
Ikeuchi S. 1986 Baryon Clump within an Extended Dark Matter Ap&SS 118 509 -    DOI : 10.1007/BF00651178
Iliev I. T. , Shapiro P. R. , Ferrara A. , Martel H. 2002 On the Direct Detectability of the Cosmic Dark Ages: 21-cm Emission from Minihalos ApJ 572 123 -    DOI : 10.1086/341869
Iliev I. T. , Scannapieco E. , Martel H. , Shapiro P. R. 2003 Non-Linear Clustering during the Cosmic Dark Ages and its Effect on the 21-cm Back- ground from Minihaloes MNRAS 341 81 -    DOI : 10.1046/j.1365-8711.2003.06410.x
Iliev I. T. , Shapiro P. R. , Raga A. C. 2005 Minihalo Photoevaporation during Cosmic Reionization: Evaporation Times and Photon Consumption Rates MNRAS 361 405 -    DOI : 10.1111/j.1365-2966.2005.09155.x
Iliev I. T. , Mellema G. , Pen U.-L. , Merz H. , Shapiro P. R. , Alvarez M. A 2006 Simulating Cosmic Reionization at Large Scales – I The Geometry of Reionization 369 1625 -
Iliev I. T. , Mellema G. , Shapiro P. R. , Pen U. L. 2007 Self-Regulated Reionization MNRAS 376 534 -    DOI : 10.1111/j.1365-2966.2007.11482.x
Iliev I. T. , Mellema G. , Ahn K. , Shapiro P. R. , Mao Y. , Pen U. L. 2013 Simulating Cosmic Reionization: How Large a Volume is Large Enough?
Kim S. , Park C. 2007 Topology of H I Gas Dis- tribution in the Large Magellanic Cloud ApJ 663 244 -    DOI : 10.1086/518470
Kim J. , Park C. , Gott J. R. , Dubinski J. 2009 The Horizon Run N-body Simulation: Baryon Acoustic Oscillations and Topology of Large Scale Structure of the Universe ApJ 701 1547 -    DOI : 10.1088/0004-637X/701/2/1547
Kim J. , Pen U.-L. 2010 Redshifted 21-cm Signals in the Dark Ages
Kogut A. 1993 Topology of the COBE-DMR First Year Sky Map BAAS 183 121.03 -
Kogut A. , Banday A. J. , Bennett C. L. , Gorski K. , Hinshaw G. , Smoot G. F. , Wright E. L. 1996 Tests for Non-Gaussian Statistics in the DMR Four- Year Sky Maps ApJ 464 L29 -    DOI : 10.1086/310078
Komatsu E. , Dunkley J. , Nolta M. R. , Bennett C. L. , Gold B. , Hinshaw G. , Jarosik N. , Larson D. , Limon M. , Page L. , Spergel D. N. , Halpern M. , Hill R. S. , Kogut A. , Meyer S. S. , Tucker G. S. , Weiland J. L. , Wollack E. , Wright E. L. 2009 Five-Year Wilkinson Microwave Anisotropy Probe (WMAP) Observations: Cosmological Interpretation ApJS 180 330 -    DOI : 10.1088/0067-0049/180/2/330
La Plante P. , Battaglia N. , Natarajan A. 2013 Reionization on Large Scales IV: Predictions for the 21 cm Signal Incorporating the Light Cone Effect
Lee K.-G. , Cen R. , Gott J. R. , Trac H. 2008 The Topology of Cosmological Reionization ApJ 675 8 -    DOI : 10.1086/525520
Madau P. , Meiksin A. , Rees M. J. 1997 21-cm Tomography of the Intergalactic Medium at High Redshift ApJ 475 429 -    DOI : 10.1086/303549
McQuinn M. , Lidz A. , Zahn O. , Dutta S. , Hernquist L. , Zaldarriaga M. 2007 The Morphology of H II Regions during Reionization MNRAS 377 1043 -    DOI : 10.1111/j.1365-2966.2007.11489.x
Mellema G. , Iliev I. T. , Alvarez, M A. , Shapiro P. R. 2006 C2-Ray: A New Method for Photon-Conserving Transport of Ionizing Radiation New Astron. 11 374 -    DOI : 10.1016/j.newast.2005.09.004
Mellema G. , Iliev I. T. , Pen U. L. , Shapiro P. R. 2006 Simulating Cosmic Reionization at Large Scales II: the 21-cm Emission Features and Statistical Signals MNRAS 372 679 -    DOI : 10.1111/j.1365-2966.2006.10919.x
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. , Jelic 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. , Zaroubi S. 2013 Reionization and the Cosmic Dawn with the Square Kilometre Array Experimental Astronomy 36 235 -    DOI : 10.1007/s10686-013-9334-5
Melott A. L. , Cohen A. P. , Hamilton A. J. S. , Gott J. R. , Weinberg D. H. 1989 Topology of Large Scale Structure. 4. Topology in Two-Dimensions ApJ 345 618 -    DOI : 10.1086/167935
Mesinger A. , Ferrara A. , Spiegel D. S. 2013 Signatures of X-Rays in the Early Universe MNRAS 431 621 -    DOI : 10.1093/mnras/stt198
Morales M. F. , Hewitt J. 2004 Toward Epoch of Re-Ionization Measurements with Wide-Field LO- FAR Observations ApJ 615 7 -    DOI : 10.1086/424437
Park C. , Gott J. R. , Melott A. L. , Karachentsev I. D. 1992 The Topology of Large Scale Structure. 6. Slices of the Universe ApJ 387 1 -    DOI : 10.1086/171055
Park C. , Colley W. N. , Gott J. R. , Ratra B. , Spergel D. N. , Sugiyama N 1998 CMB Anisotropy Correlation Function and Topology from Simulated Maps for MAP ApJ 506 473 -    DOI : 10.1086/306259
Park C. , Gott J. R. 2001 Topology of the Galaxy Distribution in the Hubble Deep Fields ApJ 553 33 -    DOI : 10.1086/320640
Park C. G. , Park C. , Ratra B. , Tegmark M. 2001 Gaussianity of Degree-Scale Cosmic Microwave Background Anisotropy Observations ApJ 556 582 -    DOI : 10.1086/321591
Park C. , Kim Y. R. 2010 Large-Scale Structure of the Universe as a Cosmic Standard Ruler ApJ 715 L185 -    DOI : 10.1088/2041-8205/715/2/L185
Parsons A. R. , Backer D. C. , Foster G. S. , Wright M. C. H. , Bradley R. F. , Gugliucci N. E. , Parashare C. R. , Benoit E. E. , Aguirre J. E. , Jacobs D. C. , Carilli C. L. , Herne D. , Lynch M. J. , Manley J. R. , Werthimer D. J. 2010 The Precision Array for Probing the Epoch of Re-Ionization: Eight Station Results AJ 139 1468 -    DOI : 10.1088/0004-6256/139/4/1468
Pritchard J. R. , Loeb A. 2008 Evolution of the 21cm Signal throughout Cosmic History Phys. Rev. D 78 103511 -    DOI : 10.1103/PhysRevD.78.103511
Purcell E. M. , Field G. B. 1956 Influence of Col- lisions upon Population of Hyperfine States in Hydrogen ApJ 124 542 -    DOI : 10.1086/146259
Rees M. J. 1986 Lyman Absorption Lines in Quasar Spectra-Evidence for Gravitationally-Confined Gas in Dark Minihaloes MNRAS 218 25 -    DOI : 10.1093/mnras/218.1.25P
Shapiro P. R. , Ahn K. , Alvarez M. , Iliev I. T. , Martel H. , Ryu D. 2006 The 21 cm Background from the Cosmic Dark Ages: Minihalos and the In- tergalactic Medium before Reionization ApJ 646 681 -    DOI : 10.1086/504972
Shapiro P. R. , Mao Y. , Iliev I. T. , Mellema G. , Datta K. K. , Ahn K. , Koda J. 2013 Will Nonlinear Peculiar Velocity and Inhomogeneous Reion- ization Spoil 21 cm Cosmology from the Epoch of Reionization? Phys. Rev. Lett. 110 151301 -    DOI : 10.1103/PhysRevLett.110.151301
Shin M.-S. , Trac H. , Cen R. 2008 Cosmological H II Bubble Growth during Reionization ApJ 681 756 -    DOI : 10.1086/588247
Lubin P. , Mather J. , Meyer S. S. , Moseley S. H. , Murdock T. , Rokke L. , Silverberg R. F. , Tenorio L. , Weiss R. , Wilkinson D. T. 1992 Structure in the COBE Differential Microwave Radiometer First Year Maps ApJ 396 L1 -    DOI : 10.1086/186504
Wyithe J. S. B. , Loeb A. 2003 Was the Universe Reionized by Massive Population-III Stars? ApJ 588 L69 -    DOI : 10.1086/375682
Zahn O. , Lidz A. , McQuinn M. , Dutta S. , Hernquist L. , Zaldarriaga M. , Furlanetto S. R. 2006 Simulations and Analytic Calculations of Bubble Growth During Hydrogen Reionization ApJ 654 12 -