2D GENUS TOPOLOGY OF 21-CM DIFFERENTIAL BRIGHTNESS TEMPERATURE DURING COSMIC REIONIZATION

Abstract

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.

1. INTRODUCTION

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. SIMULATIONS OF REDSHIFTED 21-CM SIGNALS

2.1 N-body Simulations

We completed a 20483 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−1Mpc scale. The cubic box size is 66 h−1Mpc 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 20483 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 108M⊙ (corresponding to 30 particles). At each epoch, we calculate matter density fields on 20483 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 2563 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 < 109M⊙) and high-mass (M > 109M⊙) 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 5123 particles on a 5123 mesh on a 64 h−1Mpc box using GOTPM, resulting in somewhat poorer halomass resolution of 2 × 109M⊙. We created a binneddown density field on a 1283 mesh, similarly, on which radiative transfer is calculated.

We use the results from the high resolution (20483 N-body) run for two cosmic reionization simulations (Case 3 and 4), and those from the low resolution (5123 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.

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” 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 C2-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γ ≡ fesc f⋆ Ni over simulation, where fesc 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γMtot 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 (108 ≲ M/M⊙ ≳ 109) and high-mass (M/M⊙ ≳ 109). 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 H2-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, 105 ≲ M/M⊙ ≲ 108, 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 H2 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 (n1) and lower (n0) states are determined by the spin temperature TS, such that

Several pumping mechanisms determine TS as follows:

where TCMB, Tα, TK, yα, and yc 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 > TCMB, or in absorption if TS < TCMB. This is quantified by the 21-cm differential brightness temperature (e.g., Morales & Hewitt 2004)

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)

δ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,

We simply generate 3D maps of δTb on the radiative transfer grid over all redshifts using Eq. 5. Note that in the limit TS ≫ TCMB, δ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).

 

3. 2D GENUS OF 21-CM DIFFERENTIAL BRIGHTNESS TEMPERATURE

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

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, g2D changes. Therefore, g2D 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

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 g2D,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, g2D,eff ≡ g2D − 2f, where f is the areal fraction of hot spots on the sky. g2D,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 δTb, th to construct 2D genus curve g2D(δTb, 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. RESULTS

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

and the transverse angle

where 𝜈0 ≡ 1420MHz is the rest frame frequency of the line, DA(z) is the angular diameter distance, Lbox is the comoving length of our simulation box, c is the speed of light and H0 ≡ 100 h kms−1Mpc−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

where N∆ν ≡ ∆L(∆ν) /lcell 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,

is the comoving length along LOS corresponding to ∆ν, and lcell 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 (Lbox/4). These maps will be almostmutually independent, because Lbox/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

with a full width at half maximum (FWHM) ∆θ = . We also use a compensated Gaussian filter given by

with a FWHM given by

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 WCG as an extreme case for “dirty” beams, clearly distinct from WG, because this filter is somewhat pathological as it suppresses both small and large scale features. When , WCG(θ) < 0 such that a data field convolved with WCG(θ) 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 ≫ TCMB.

Fig. 1 shows how the actual 2D δTb field (top rows) will be observed in different filtering schemes (middle: WG; bottom: WCG). 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. WC 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.

Fig. 1.— 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 WCG, 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, WCG changes the topology of the 21-cm signal. For most of our analysis, therefore, we just use WG. 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 × 109M⊙) halos, which collapse much later than low-mass (108 ≲ 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.

Fig. 2.— 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.

Fig. 3.— 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 ℓ ∼ 104, or the comoving length scale x ∼ 6 − 7 h−1Mpc, while f125_1000S shows a peak at ℓ ∼ 3 × 104, or x ∼ 2 h−1Mpc, and f125_125S shows a long uniform tail from ℓ ∼ 3×104 to ℓ ∼ 4 × 104. The behavior of 2D angular power spectra indicates stronger merger of (otherwise) individual bubbles into larger scale at x ∼ 6 − 7 h−1Mpc in f5, f250 and f125_125S than in f125_1000S. Even though f125_125S has a tail up to ℓ ∼ 4×104 (x ∼ 1 h−1Mpc), 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 × 104 (x ∼ 2 h−1Mpc), is shown at all redshifts in f125_1000S, due to the clustering of small halos of the highest efficiency.

Fig. 4.— 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 g2D(δ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 WG is used in all cases, except for f125_1000S to which WCG is also applied to understand the impact of beam shape on g2D.

4.3.1 2D Genus of Density Fluctuation

Is there any useful template genus curve to be compared with g2D(δTb)? There is such one indeed, which is g2D of an artificial field of δTb which reflects the density fluctuation only. Let us denote this quantity by

Distribution of δTb,n will be Gaussian in the linear regime (|δ| ≪ 1), which will provide well-defined, analytical g2D,Gauss(δTb,n). Even when g2D(δTb,n) starts to deviate from g2D,Gauss(δTb,n) due to nonlinear evolution of high density peaks, g2D(δ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.

g2D(δ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 g2D,Gauss, which is symmetrical around ν = 0, while high-density peaks (on high ν) that evolved non- linearly deviates from the Gaussian distribution to also make g2D deviate from g2D,Gauss. Because we will use filtered δTb maps, g2D of filtered δTb,n will be the template to be used throughout this paper. Filtering makes both the amplitude and width of g2D(δTb,n) shrink (Fig. 5). Nevertheless, even after filtering, regions with positive g2D correspond to overdense regions and negative g2D to underdense regions, because the meandensity point (where the g2D(δTb,n) curve crosses the δTb axis) is almost unshifted. This will work as a perfect indicator, if comparison is made with g2D(δTb), about how different the δTb field is from the underlying density field. g2D(δ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.

Fig. 5.— 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 g2D and Impact of Beam Shape

We found that the evolution of g2D(δ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 g2D(δTb) clearly shows how reionization proceeds in time. The left panel of Fig. 6 shows WG-filtered g2D(δTb,n) of the underlying matter density field (blue dashed line) and g2D(δTb) of the differential brightness temperature field of HI gas (red line). In the earliest phase of reionization (xv = 0.04), g2D(δTb) is much smaller than g2D(δ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 g2D(δ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 δTb,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 g2D 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 g2D(δ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 g2D 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 WCG on g2D is demonstrated on the right panel of Fig. 6. Note that the Fourier transform of WCG is proportional to k2σ2 exp(−k2σ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 WCG, δTbs 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 g2D(δ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 WCG than that of WG (Gaussian smoothing only), the interpretation becomes less transparent. We also note that WCG filters out almost all scales except for the filtering scale, while WG filters out only those scales smaller than the filtering scale.

Fig. 6.— 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 g2D(δTb) is somewhat lost with WCG. 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 g2D(δTb) and g2D(δTb,n) of all cases at xv = 0.4. The impact of small-mass (108 ≲ M/M⊙ ≲ 109) 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 g2D 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.

Fig. 7.— 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, g2D 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 g2D curves are mapped too deep in (newly filtered) underdense regions. Nevertheless, the amplitude of g2D is the largest and the rightmost wing of g2D is located at the smallest δTb in f125_1000S in all varying filtering scales, just as when ∆ν = 0.2′MHz is applied.

Fig. 8.— Same as Fig. 7, but with 2′ FWHM Gaussian beam applied.

Fig. 9.— Same as Fig. 7, but with 3′ FWHM Gaussian beam applied.

Fig. 10.— 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 g2D 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 f5, 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 g2D(δ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.

Fig. 11.— 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).

 

5. SUMMARY AND DISCUSSION

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 g2D of the underlying density fluctuation, g2D(δTb,n), clearly shows the evolution of the reionization process. g2D(δTb) is found to deviate quickly from the g2D 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.