Mon. Not. R. Astron. Soc. 000, 1–13 (2009) Printed 20 July 2009 (MN LATEX style file v2.

2)

Measuring the History of Cosmic Reionization using the

21-cm PDF from Simulations
arXiv:0907.2932v2 [astro-ph.CO] 20 Jul 2009

Kazuhide Ichikawa1,2,3, Rennan Barkana1,4,5 ⋆ , Ilian T. Iliev6,7, Garrelt Mellema8,

Paul
1
R. Shapiro9
Institute for Cosmic Ray Research, University of Tokyo, Kashiwa 277-8582, Japan
2 Department of Physics and Astronomy, University College London, Gower Street, London, WC1E 6BT, U.K.
3 Department of Micro Engineering, Kyoto University, Kyoto 606-8501, Japan
4 Division of Physics, Mathematics and Astronomy, California Institute of Technology, Mail Code 130-33, Pasadena, CA 91125, USA
5 Raymond and Beverly Sackler School of Physics and Astronomy, Tel Aviv University, Tel Aviv 69978, Israel
6 Universität Zürich, Institut für Theoretische Physik, Winterthurerstrasse 190, CH-8057 Zürich, Switzerland
7 Astronomy Centre, Department of Physics & Astronomy, University of Sussex, Brighton BN1 9QH, UK
8 Department of Astronomy, AlbaNova University Center, Stockholm University, SE 10691 Stockholm, Sweden
9 Department of Astronomy, University of Texas, Austin, TX 78712-1083, USA ;

Texas Cosmology Center, The University of Texas at Austin, TX 78712, USA

20 July 2009

ABSTRACT
The 21-cm PDF (i.e., distribution of pixel brightness temperatures) is expected to be
highly non-Gaussian during reionization and to provide important information on the
distribution of density and ionization. We measure the 21-cm PDF as a function of
redshift in a large simulation of cosmic reionization and propose a simple empirical
fit. Guided by the simulated PDF, we then carry out a maximum likelihood analysis
of the ability of upcoming experiments to measure the shape of the 21-cm PDF and
derive from it the cosmic reionization history. Under the strongest assumptions, we
find that upcoming experiments can measure the reionization history in the mid to
late stages of reionization to 1 − 10% accuracy. Under a more flexible approach that
allows for four free parameters at each redshift, a similar accuracy requires the lower
noise levels of second-generation 21-cm experiments.
Key words: galaxies:high-redshift – cosmology:theory – galaxies:formation

1 INTRODUCTION Low Frequency Array (www.lofar.org), the Giant Me-

trewave Radio Telescope (gmrt.ncra.tifr.res.in), and the
The earliest generations of stars are thought to have trans-
Precision Array to Probe the Epoch of Reionization
formed the universe from darkness to light and to have reion-
(astro.berkeley.edu/∼dbacker/eor/).
ized and heated the intergalactic medium (Barkana & Loeb
2001). Knowing how the reionization process happened is Studies of statistics of the 21-cm fluctuations have
a primary goal of cosmologists, because this would tell us focused on the two-point correlation function (or power
when the early stars formed and in what kinds of galaxies. spectrum) of the 21-cm brightness temperature. This is
The clustering of these galaxies is particularly interesting true both for analytical and numerical studies and anal-
since it is driven by large-scale density fluctuations in the yses of the expected sensitivity of the new experiments
dark matter (Barkana & Loeb 2004). While the distribution (Bowman, Morales, & Hewitt 2006; McQuinn et al. 2006).
of neutral hydrogen during reionization can in principle be The power spectrum is the natural statistic at very high red-
measured from maps of 21-cm emission by neutral hydro- shifts, as it contains all the available statistical information
gen, upcoming experiments are expected to be able to de- as long as Gaussian primordial density fluctuations drive the
tect ionization fluctuations only statistically (for reviews see, 21-cm fluctuations. More generally, the power spectrum is
e.g., Furlanetto et al. 2006; Barkana & Loeb 2007). Current also closely related to the directly observed radio visibili-
observational efforts include the Murchison Widefield Ar- ties. Now, during reionization the hydrogen distribution is a
ray (MWA, www.haystack.mit.edu/ast/arrays/mwa/), the highly non-linear function of the distribution of the underly-
ing ionizing sources. This follows most simply from the fact
that the H I fraction is constrained to vary between 0 and
⋆ E-mail: barkana@wise.tau.ac.il (RB) 1, and this range is fully covered in any scenario driven by
2 Ichikawa, Barkana, Iliev, Mellema, Shapiro
stars, in which the intergalactic medium is sharply divided 2 BASIC METHOD
between H I and H II regions. The resulting non-Gaussianity
In this section we develop our basic statistical method for
(Bharadwaj & Ali 2005) raises the possibility of using com-
fitting the PDF. While the statistical approach is general, for
plementary statistics to measuring additional information
concreteness we develop it within the context of a simple toy
that is not directly derivable from the power spectrum (e.g.,
model for the PDF. We also use this toy, double-Gaussian
Furlanetto et al. 2004; Saiyad-Ali et al. 2006).
model in order to get a crude quantitative intuition on how
hard it is to measure the 21-cm PDF. We note that we fol-
Numerical simulations have recently begun to reach the
low to some degree Oh et al. (2009), who considered such a
large scales (of order 100 Mpc) needed to capture the evolu-
double-Gaussian toy model and made a signal-to-noise study
tion of the intergalactic medium (IGM) during reionization
of this model with their analysis method.
(Iliev et al. 2006b; Mellema et al. 2006b; Zahn et al. 2007;
Santos et al. 2008). These simulations account accurately for
gravitational evolution and the radiative transfer of ionizing 2.1 A Toy Model for the PDF
photons, but still crudely for gas dynamics and star forma-
tion. Analytically, Furlanetto et al. (2004) used the statis- It is useful to have a simple PDF example on which to de-
tics of a random walk with a linear barrier to model the velop and test our methods. We present here a simplified
H II bubble size distribution during the reionization epoch. toy model that captures the main qualitative features of
Schematic approximations were developed for the two-point the PDF as seen in the simulations (and shown later in
correlation function (Furlanetto et al. 2004; McQuinn et al. the paper) during the central stage of reionization, when
2005), but recently Barkana (2007) developed an accurate, the cosmic ionization fraction x̄i ∼ 0.3 − 0.6. The PDF
self-consistent analytical expression for the full two-point at this stage has a sharp peak at a differential brightness
distribution within the Furlanetto et al. (2004) model, and temperature (defined as the difference between the actual
in particular used it to calculate the 21-cm correlation func- brightness temperature and the temperature of the cosmic
tion. microwave background at the same frequency) of Tb = 0 mK
corresponding to fully ionized pixels, and another peak at
Tb ∼ 20 mK corresponding to mostly neutral pixels, with
Noting the expected non-Gaussianity and the impor-
a rapidly declining probability at values above 20 mK, and
tance of additional statistics, Furlanetto et al. (2004) also
a smooth probability density in between the peaks that is
calculated the one-point probability distribution function
lower than the height of either peak. In the observations,
(PDF) of the 21-cm brightness temperature at a point. The
this physical PDF is convolved with a broad Gaussian cor-
PDF has begun to be explored in numerical simulations as
responding to the thermal noise, resulting in both positive
well (Ciardi & Madau 2003; Mellema et al. 2006b). Some of
and negative values of Tb . In the limit when we approxi-
the additional information available in the PDF can be cap-
mate both peaks as delta functions and neglect the physical
tured by its skewness (Wyithe & Morales 2007; Harker et al.
PDF at other points, the observed PDF becomes a sum of
2009). Barkana & Loeb (2008) have also considered the dif-
two Gaussians with equal standard deviations σ. While cer-
ference PDF, a two-dimensional function that generalizes
tainly highly simplified, this model does capture the main
both the one-point PDF and the correlation function and
question (relevant especially for low signal-to-noise data, i.e.,
yields additional information beyond those statistics.
when σ ≫ 20 mK) of whether it is at all possible to tell apart
the two peaks and not confuse them with a convolved single
Recently, Oh et al. (2009) have quantitatively consid-
peak (i.e., a single Gaussian).
ered the ability of upcoming experiments to determine the
Thus, we consider two Gaussian distributions with
cosmic reionization history from maximum likelihood fitting
equal standard deviation σ (where σ represents the mea-
of the 21-cm PDF. They specifically used mixture modeling
surement noise level). In the toy model we use a dimension-
of the PDF. In this paper we develop a method for statistical
less s as the dependent variable (which represents Tb in the
analysis of the PDF that is simpler and more efficient (allow-
real PDF). The Gaussian representing the reionized pixels
ing, in particular, binning of the PDF). We use our method
is centered at s = 0, while the neutral pixels are represented
to present a quantitative analysis of whether upcoming and
by a Gaussian centered at s = sG . The fraction of the total
future experiments can measure the detailed shape of the
probability contained in the first Gaussian is α. The total
21-cm PDF and derive from it the cosmic reionization his-
distribution is therefore
tory. In section 2 we develop our basic statistical method for
fitting the 21-cm PDF, and test it on a simple toy model for p(s) = αG(s, σ) + (1 − α)G(s − sG , σ) , (1)
the PDF. We then measure and follow the evolution of the
where
PDF in a large N-body and radiative transfer simulation of
1
cosmic reionization; since previous analytical models of the G(x, σ) ≡ √ exp(−x2 /2σ 2 ) . (2)
PDF differ qualitatively from the PDF in the simulation, 2πσ
here we simply fit the simulated PDF with an empirical, Since in the real case only differences in Tb can be mea-
four-parameter model (section 3). Finally, we present the sured, and not the absolute Tb (which is dominated by fore-
expected accuracy of reconstructing the 21-cm PDF and the grounds), in the toy model we assume that the absolute s
cosmic reionization history based on the simulated PDF, ei- cannot be measured. A simple practical way to do this is to
ther with strict assumptions that lead to one free parameter always measure s with respect to its average value according
at each redshift (section 4), or with a more flexible approach to the PDF of s; we do this separately in each model and in
that allows for four free parameters (section 5). We summa- each simulated data set, and thus only compare the relative
rize our conclusions in section 6. distributions between each model and each data set.
 Measuring Reionization using the 21-cm PDF 3
2.2 Maximum Likelihood and the C-Statistic 2.3 Results for the Toy Model
In this subsection we develop our basic statistical method For the toy, double-Gaussian model, the parameters we wish
for fitting the PDF, referring to the above toy model as an to fit to mock data sets are sG and α. Note that we assume
example for the PDF. In general, we can create mock data that σ is known, as we expect that the level of thermal noise
sets by randomly generating Np values of s from a given per pixel will be known in the 21-cm experiments, given
p(s) distribution, and we can then try to estimate the best- the known array properties and the measured foreground
fitting parameters with a maximum likelihood method. For level. We perform 1000 Monte Carlo for each input model,
a given mock observed PDF, as given by the Np generated and thus obtain the full distribution of reconstructed model
values of s, we wish to find the best-fitting model PDF p(s) parameters. In order to develop intuition on how hard it is
by maximizing the likelihood L that the Np values si (i = 1, to measure the PDF, we define a parameter η that captures
2, . . . , Np ) came from p(s). Since the different values si are a simplistic notion of the total signal-to-noise ratio:
independent, this probability (apart from fixed ∆s factors)
is simply
Np
Y
L= p(si ) . (3) sG
 p
η≡ Np , (7)
i=1 σ
Now, it is standard to replace the problem of maximizing
the likelihood L with a minimization of −2 ln L, which in
this case is p
motivated by sG as a measure for the signal and σ/ Np
Np
X as a measure for the effective noise after Np measurements
− 2 ln L = −2 ln p(si ) . (4) with noise σ in each. Of course, the ability to detect the
i=1 two separate peaks also depends on α, in that values close
In comparing the data to a potential model, we bin the to 0 or 1 make one of the peaks insignificant. For a fixed α,
values of s in order to have a manageable number of bins though, we might naively expect that the accuracy of the
(NB = 1000) even when Np is very large. This is justified as reconstructed values of sG and α would not change with the
long as the bin width is much smaller than any s-scale that input value of sG , as long as we change Np so as to keep the
we hope to resolve in the PDF. We have explicitly checked combination η fixed.
that using NB = 1000 bins (with the C-statistic, see below) To test this, we fix the input α = 0.4 and sG = 1,
gives the same results as applying equation (4) directly, even and vary σ and Np together so as to keep η fixed. We test
for the largest values of Np that we use in this paper. Now, η = 400 and 4000, values comparable to those expected in
when the expected (according to a model p(s)) number of the real experiments discussed later in the paper. The Monte
points nexp,j in each bin j is large (i.e., nexp,j ≫ 1), the Carlo results are summarized in Figure 1. The results show
√
actual number nj has a standard error of nexp,j , and we that the parameters can be accurately reconstructed as long
can find the best-fitting model by minimizing a standard χ2 as the signal-to-noise per sample (or per pixel in real data)
statistic: sG /σ > 1. As long as this is the case, the relative error in
NB
sG and α is no worse than 4% (η = 400) or 0.4% (η =
X (nj − nexp,j )2 4000), and the average reconstructed values are essentially
χ2 = . (5)
nexp,j unbiased. However, once sG /σ drops below unity (i.e., σ > 1
j=1
in this particular case), the errors increase rapidly with σ, so
However, in modeling the PDF we often wish to include that for η = 400 reconstruction is impossible when σ = 10
a wide range of s, including some bins where the model prob- (i.e., both the bias and spread are of order unity) , and for
ability density p(s) is very low. When nexp,j is small, the χ2 η = 4000 the errors increase when σ = 4 to a 5% relative
distribution with its assumption of a Gaussian distribution spread in α.
for each nj severely underestimates the fluctuations in nj The reason for these increasing errors is parameter de-
compared to the correct Poisson distribution. Thus, equa- generacy, as illustrated in Figure 2 for η = 400. While for
tion (5) can lead to major errors if nexp,j ≪ 1 in any bin. In σ = 1 the reconstructed parameter distribution is fairly sym-
this situation, the correct statistic to use is the C-statistic metrical about the input values of sG and α, resembling a
(Cash 1979), derived from the Poisson distribution just as standard error ellipse, larger σ values produce a stretched
the χ2 statistic is derived from the Gaussian distribution. error contour that displays a clear (partial) degeneracy be-
The C-statistic is defined as tween the parameters sG and α. Intuitively, when sG /σ ≪ 1
NB the PDF consists of a narrow input signal (two peaks sep-
arated by sG , in the case of the toy model) convolved with
X
C=2 (nexp,j − nj ln nexp,j ) . (6)
j=1
a broad Gaussian of width σ. The result is a broad Gaus-
sian of width σ, with small bumps (distortions). Apparently
Note that the Poisson distribution also has a factor of nj ! these small bumps can be produced with very different pa-
in the denominator, which results in an additional ln nj ! rameter combinations, resulting in a degeneracy that leads
term within the sum in equation (6), but this term does not to a large uncertainty when fitting models. While we have
depend on the model parameters (which enter only through considered here a simple toy model, a similar degeneracy is
nexp,j ) and can thus be dropped from the minimization. encountered with the real 21-cm PDF, as discussed below.
4 Ichikawa, Barkana, Iliev, Mellema, Shapiro
3 THE 21-CM PDF IN SIMULATIONS
3.1 Numerical Simulation
In this paper we utilize a large-scale N-body and radia-
tive transfer simulation of cosmic reionization following the
methodology first presented in Iliev et al. (2006b). The cos-
mological structure formation and evolution is followed with
a particle-mesh N-body code called PMFAST (Merz et al.
2005). These N-body results then provide the evolving den-
sity field of the IGM, as well as the location and mass of
all the halo sources, as input to a separate radiative trans-
fer simulation of inhomogeneous reionization. The latter is
performed with the C2 −Ray (Conservative, Causal Ray-
Tracing) code, a regular-grid, ray-tracing, radiative transfer
and nonequilibrium chemistry code (Mellema et al. 2006a).
The ionizing radiation is ray-traced from every source cell
to every grid cell at a given timestep using a method of
short characteristics. C2 −Ray is designed to be explicitly
photon-conserving in both space and time, which ensures an
accurate tracking of ionization fronts, independently of the
spatial and time resolution. This is true even for grid cells
which are very optically thick to ionizing photons and time
Figure 1. For each model parameter x reconstructed in each steps long compared to the ionization time of the atoms,
Monte Carlo trial, we show the bias in the average (i.e., the en- which results in high efficiency. The code has been tested
semble average p hxi minus the input value xin ) and the standard against analytical solutions (Mellema et al. 2006a), and di-
deviation σx = hx2 i − hxi2 . We consider the model parameters rectly compared with other radiative transfer methods on a
sG (solid curves, input value 1) and α (dashed curves, input value standardized set of benchmark problems (Iliev et al. 2006a,
0.4), as a function of the noise level (i.e., width of each Gaussian) 2009).
σ, where η is held fixed at 400 (left panels) or 4000 (right panels). We simulated the ΛCDM universe with 16243 dark mat-
ter particles of mass 2.2 × 107 M⊙ , in a comoving simulation
volume of (100 h−1 Mpc)3 . This allowed us to resolve (with
100 particles or more per halo) all halos of mass 2.2×109 M⊙
and above. The radiative transfer grid has 2033 cells. The
H-ionizing photon luminosities per halo in our cosmic reion-
ization simulations are assigned as follows. A halo of mass
M is assumed to have converted a mass M ·(Ωb /Ωm )·f∗ into
stars, where f∗ is the star formation efficiency. Halo catalogs
are discrete in time, because N-body density fields are stored
every ∆t ∼ 20 Myrs and the corresponding halo catalogs are
produced at the same time. If each source forms stars over
a period of time ∆t and each stellar nucleus1 produces Ni
ionizing photons per stellar lifetime and is used only once
per ∆t, and if a fraction fesc of these photons escape into
the IGM, then the ionizing photon number luminosity of a
halo of mass M is given by
Ni · fesc · f∗ · M (Ωb /Ωm )
Qi = , (8)
∆t · µmH
where mH is the mass of a hydrogen atom and µ = 1.22
so that µmH is the mean mass per nucleus. In this model,
stars are produced in a burst, and they keep radiating with
a fixed Qi for ∆t ≃ 20 Myrs. We choose here a specific case,
first presented (and labeled f250) in Iliev et al. (2007) and
Figure 2. Distribution of reconstructed model parameters sG further discussed in Iliev et al. (2008). In this scenario, halos
and α in 1000 Monte Carlo simulations. The input parameter are assumed to host relatively low efficiency emitters, with
values are sG = 1 and α = 0.4. We vary σ keeping η = 400 fixed, fγ ≡ f∗ fesc Ni = 250 (corresponding, e.g., to Pop II stars
so that the number of samples is Np = 160, 000 σ2 . Different
with a Salpeter IMF).
panels cover different x ranges, but all x axes are shown on the
same scale for easy comparison. In the σ 6 2 panels, small tick The simulation we use in this work assumes a flat
marks are at 0.75 and 1.25.

1 Note that we defined this number per atomic nucleus rather

than per baryon in stars.
 Measuring Reionization using the 21-cm PDF 5
(Ωk = 0) ΛCDM cosmology. The simulation is based on
the WMAP 3-year results, which derived the parameters
(Ωm , ΩΛ , Ωb , h, σ8 , n) = (0.24, 0.76, 0.042, 0.73, 0.74, 0.95)
(Spergel et al. 2007). Here Ωm , ΩΛ , and Ωb are the total
matter, vacuum, and baryonic densities in units of the crit-
ical density, σ8 is the root-mean-square density fluctuation
on the scale of 8h−1 Mpc linearly extrapolated to the present,
and n is the power-law index of the primordial power spec-
trum of density fluctuations.

3.2 The Simulated 21-cm PDF

During cosmic reionization, we assume that there are suffi-
cient radiation backgrounds of X-rays and of Lyα photons so
that the cosmic gas has been heated to well above the cosmic
microwave background temperature and the 21-cm level oc-
cupations have come into equilibrium with the gas tempera-
ture. In this case, the observed 21-cm differential brightness
temperature (i.e., relative to the cosmic microwave back-
ground) is independent of the spin temperature and, for our
assumed cosmological parameters, is given by (Madau et al.
1997)
Figure 3. The global progress of cosmic reionization in the sim-
 r ulation, as a function of the redshift z. Bottom panel: we show
Ωb h 0.3 1+z
  
Tb = T̃b Ψ; T̃b = 23.7 mK , (9) the mass-weighted ionized fraction x̄i (solid curve) and the corre-
0.032 Ωm 8
sponding neutral fraction x̄n = 1 − x̄i (dashed curve). Top panel:
with Ψ = xn [1 + δ], where xn is the neutral hydrogen frac- we show the cosmic mean 21-cm differential brightness tempera-
tion and δ is the relative density fluctuation. Under these ture T̄b in the simulation (solid curve), and the mean Tb expected
conditions, the 21-cm fluctuations are thus determined by for a neutral universe of uniform density (dotted curve). Also
indicated in each panel are the 26 output redshifts used in the
fluctuations
R in Ψ. We denote the PDF by p(Tb ), normalized analysis below (points).
so that p(Tb )dTb = 1.
To calculate the 21-cm PDF, we smooth the 21-cm emis-
sion intensity over our full simulation volume with a cubical
top-hat filter (sometimes referred to as “boxcar” averaging) drops more slowly with Tb than the Gaussian fit (more on
of a pre-determined size Rpix . We then assemble the PDF of the fitting function below); this results from the non-linear
the resulting values over a fine grid, much finer than Rpix . growth of density fluctuations.
This effectively smooths out the fluctuations in the PDF and As reionization gets under way, the high-density tail
yields a smooth function, but we note that there is still a real drops off and (coincidentally) approaches the Gaussian
sample (or “cosmic”) variance limit on the accuracy of our shape, as high-density pixels are more likely to be partially
simulated PDF, resulting from the limited number of inde- or fully ionized and thus have their Tb reduced. When x̄i
pendent volumes of size Rpix within our simulation box. We reaches a fraction of a percent, the still fairly Gaussian PDF
use Rpix = 5 h−1 Mpc, 10 h−1 Mpc, and 20 h−1 Mpc, yield- develops a significant low-Tb tail which is roughly exponen-
ing a number of independent volumes equal to 8000, 1000, tial (i.e., linear in the plot of log of the PDF). This tail cor-
and 125, respectively. The analogous results for the first- responds to pixels that are substantially ionized, i.e., where
year WMAP cosmology were previously presented, for a few a large fraction of the pixel volume partially overlaps one
redshifts only, in Mellema et al. (2006b) (with a similarly- or more ionized bubbles. Soon afterward, a significant peak
defined ionized fraction PDF shown in Iliev et al. (2006b)). can be seen near Tb = 0 mK, corresponding to fully ion-
Figure 3 shows the overall progress of reionization as a ized pixels (i.e., pixels in which the hydrogen in the IGM
function of redshift in the simulation. We calculate the PDF has been fully ionized, but there may remain a small bit
at 26 redshifts spanning a global mass-weighted ionization of high-density neutral gas). Near the mid-point of reion-
fraction x̄i from 6 × 10−6 to 99.0%, with the cosmic mean ization (x̄i = 50%), there is still a half-Gaussian peak (at
21-cm differential brightness temperature T̄b ranging from Tb ∼ 20 mK), i.e., with a Gaussian drop-off towards higher
36.5 mK to 0.27 mK. Of course, we assume that T̄b itself is Tb , now with a nearly flat exponential tail towards lower
not directly observable, due to the bright foregrounds. The Tb , and a prominent peak at Tb = 0 mK. The peak at zero
main goal of the PDF analysis is to reconstruct x̄i vs. z using increasingly dominates towards the end of reionization, as
the Tb fluctuations as captured in the PDF at each redshift. most pixels become fully ionized, but there remains an ex-
We show the measured simulation PDFs for various red- ponential tail out to higher Tb , with a cutoff (at Tb ∼ 20
shifts and Rpix = 5 h−1 Mpc in Figure 4. The PDF starts out mK).
close to Gaussian at high redshift, when the ionized volume The PDFs are shown for Rpix = 10 h−1 Mpc and
−1
is negligible and the density fluctuations on the scale Rpix 20 h Mpc in Figure 5. The qualitative evolution of the PDF
are fairly linear and thus give a Gaussian PDF. There is also throughout reionization is similar to the Rpix = 5 h−1 Mpc
a clear skewness, seen particularly in a high-density tail that case, but the PDF is narrower for larger Rpix since the 21-
6 Ichikawa, Barkana, Iliev, Mellema, Shapiro

Figure 5. Same as Figure 4 but for cubic pixels of size 10 h−1 Mpc
(left panel) or 20 h−1 Mpc (right panel). In the right panel we show
only the simulated PDFs, and the × marks the peak of the PDF
right after the midpoint of reionization (z = 8.79).

the same model used for the other PDFs. Thus, in this pa-
per we focus on the two smaller values of Rpix and only
present fits to the corresponding PDFs. Observations of the
PDF are most promising during the central stage of reion-
ization, when the PDF has two significant, well-separated
peaks rather than a single narrow peak (as is the case ei-
ther very early or very late in reionization). This two-peak
regime covers x̄i ∼ 30 − 90% for Rpix = 5 h−1 Mpc, but only
x̄i ∼ 75−95% for Rpix = 10 h−1 Mpc, because of the rarity of
fully ionized pixels in the latter case. However, even without
Figure 4. The 21-cm PDF in 5 h−1 Mpc cubic pixels, shown ver-
a strong peak at zero, the extended nearly flat (exponential)
sus the differential brightness temperature Tb . We show log10 of
the PDF, which itself is expressed in units of 1/mK. We show the part of the PDF during reionization helps in measuring the
PDF obtained from the simulation (alternating solid and dotted PDF, as we find below.
curves) and our best fits to them (alternating long-dashed and
short-dashed curves). The 26 redshifts (see Figure 3) range from
z = 15.729 (top) to 7.460 (bottom). The highest-redshift PDF is
shown at its actual value, corresponding to the labels at the top 3.3 The GED Model Fit to the Simulated PDF
of the y-axis; each subsequent PDF is shifted vertically down by
a factor of 10 in the PDF. The × mark points (where Tb equals Previous analytical models of the PDF do not describe our
the best-fit TL ) on three simulated PDFs: early in reionization simulated PDFs well. While the Gaussian at high redshift
(z = 10.08, x̄i = 0.156), right after the midpoint (z = 8.79, and the Tb = 0 mK Delta function at the end of reioniza-
x̄i = 0.530), and late in reionization (z = 7.75, x̄i = 0.948); these tion are obvious, the precise shape at intermediate redshifts
points mark the 12-redshift range that is used in the fitting of seems to depend on the precise topology of the ionized bub-
mock data in the following sections.
bles and the geometry of their overlap with the cubic pixels.
Here we take an empirical approach based on our numer-
ical simulation. Thus we use a Gaussian + Exponential +
cm fluctuations are smaller when smoothed on larger scales. Delta function (GED) model for the PDF p(Tb ). The Dirac
Also, for larger Rpix there are fewer pixels in the peak near Delta function is centered at zero, and is connected with an
Tb = 0 mK since it is more difficult to fully ionize large pix- exponential to the Gaussian. The model depends on four in-
els. The PDFs for Rpix = 20 h−1 Mpc are not so reliable, as dependent parameters: TG (center of Gaussian), σG (width
they are measured based on only 125 independent volumes. of Gaussian), cG (height of Gaussian peak) and TL (transi-
Also, their shapes differ significantly from the PDFs in the tion point between the exponential and the Gaussian). Our
smaller pixels, and so they cannot be successfully fitted with GED model is thus:
 Measuring Reionization using the 21-cm PDF 7

 p1 (Tb ) = PD δD (Tb ) + a exp(λTb ) , 0 6 Tb 6 TL observe. The fits are also quite good at the highest redshifts,
2
p(Tb ) = (Tb − TG ) (10) where the simulated PDF is essentially Gaussian except for
 p2 (Tb ) = cG exp − 2 σ 2 , Tb > TL
the skewness. This skewness, though, affects mainly the tails
G
of the distribution; e.g., at the highest redshift (z = 15.729)
where δD (x) is the Dirac delta function. The quantities a and for Rpix = 5 h−1 Mpc, ∼ 60% of the total probability is con-
λ can be expressed in terms of the above four parameters by tained at Tb values above the peak of the PDF, i.e., the
requiring the exponential and Gaussian functions to connect high-density tail adds about 10% to the 50% of a symmetri-
smoothly at Tb = TL . The conditions p1 (TL ) = p2 (TL ) and cal Gaussian. As noted above, this high-density tail declines
p′1 (TL ) = p′2 (TL ) lead to with time due to ionization offsetting the high density of
TG − TL overdense pixels. Thus, the high-Tb tail becomes well fitted
λ = 2
, (11) by the Gaussian model once x̄i rises above a few percent. At
σG
  later times the cutoff becomes somewhat steeper than the
(TL − TG )2 Gaussian fit, especially for Rpix = 10 h−1 Mpc, but this only
a = cG exp − 2
− λTL . (12)
2 σG affects the insignificant tail end of the PDF at the highest
Also, PD is determined by the requirement of normalization; Tb . For instance, for Rpix = 10 h−1 Mpc at x̄i = 0.530, the
the total integrated probability is unity if PD = 1−PE −PG , tail beyond Tb = 23 mK (where the cutoff starts to differ
where significantly from the fit) contains only 0.2% of the total
Z TL
probability.
a Another small mismatch occurs when reionization gets
PE = p1 (Tb )dTb = [exp(λTL ) − 1] , (13)
+ǫ
λ significantly under way but is still fairly early. The transition
Z ∞ q
π

TL − TG
 region from a near-Gaussian to a near-exponential shape is
PG = p2 (Tb )dTb = cG σG erfc √ .(14) not well-fit at these times by our model, and as a result the
TL
2 2 σG
fit is significantly below the low-Tb , roughly linear (expo-
Note that the parameters PD , PE and PG represent the rel- nential) tail. This mismatch is significant in the range of x̄i
ative contribution to the total probability from the delta from a few percent up to ∼ 30%, and at these redshifts this
function, the exponential function, and the Gaussian func- exponential tail typically contains only a few percent of the
tion, respectively. total probability (up to 10%).
Using the GED model, we determine the values of TG , Figure 6 shows how our model parameters vary as cos-
σG , cG and TL as functions of redshift by fitting to the sim- mic reionization progresses. The Gaussian peak position TG
ulation PDFs for pixels of 5 h−1 Mpc and 10 h−1 Mpc. In ap- and height cG both decline with time due to the increas-
proaching this fitting, we note that we focus on the main ing ionization of even low-density pixels. At least a half-
features of the PDF, and not on the fine details. In partic- Gaussian is present until x̄i ∼ 60%, but after that TL > TG
ular, we do not worry about features that contain a small and only the Gaussian cutoff remains. The parameter σG
fraction of the total probability, or on the detailed PDF remains at a value of a few mK throughout reionization; it
shape on scales finer than several mK. This is justified since gives a measure of density fluctuations, initially purely and
the observations are difficult, and most likely will not be later together with some correlation with ionization. At the
sensitive to these fine details, at least in the upcoming 21- very end of reionization, TG → 0 and then σG and cG lose
cm experiments. In addition, our simulated PDF may not their usual meaning (e.g., cG becomes an indirect parame-
be reliable in its fine details, since we are using a single, lim- terization of the normalization of the exponential portion).
ited simulated volume, and more generally, numerical simu- Figure 7 shows the evolution of the probabilities PD
lations of reionization still lack a detailed demonstration of (representing the delta function), PE (exponential), and PG
convergence. (Gaussian), which together add up to unity. The Figure
Thus, we do not try to fit the detailed peak shape at shows how the 21-cm PDF is gradually transformed from
Tb = 0, but instead represent the total probability of that a Gaussian to a delta function, with the exponential dom-
region with the Delta function. In practice we only fit to the inating at intermediate times (mid to late reionization).
data beyond the lowest values of Tb , and then set the Delta Note that in the limit of infinite resolution, we would have
function contribution PD to get the correct overall normal- PD = x̄i . With a finite resolution, PD can be thought of as
ization. Specifically, for each PDF we first find Tb,h which is the cosmic ionized fraction smoothed at the observed res-
the highest value of Tb where p(Tb,h ) > 10−4 . We then only olution. In practice, converting observed values of PD , PE ,
fit to the data with Tb > 5 mK, if Tb,h > 20 mK, or to the and PG to the true x̄i requires some modeling.
data with Tb > Tb,h /4, if Tb,h < 20 mK. At redshifts where We also calculate the variance hTb2 i − hTb i2 from the
the simulation data do not have a Delta function feature, i.e., PDF both directly from the original simulation data and
there are no pixels near Tb = 0, we make a fit constrained from our GED model fits. We plot this in Figure 8 for two
by setting PD = 1 − PE − PG = 0; this is the case at the reasons. First, the plot shows that the GED model repro-
highest redshifts, namely z > 10.924 for Rpix = 5 h−1 Mpc duces the variance of the real PDF rather well, especially
and z > 9.034 for Rpix = 10 h−1 Mpc. where the upcoming measurements are more promising (i.e.,
Our GED model fits are shown along with the PDFs later in reionization). Second, the Figure illustrates a sym-
in Figs. 4 and 5. The fits are very good during the central metry in that the variance is maximum near the midpoint of
and late stages of reionization, except for the detailed shape reionization, and has lower values both before and after the
(which we do not try to fit) of the Tb = 0 peak which ex- midpoint; this symmetry helps explain a near-degeneracy
tends out to Tb ∼ 2 − 4 mK. These are the redshifts that that we sometimes find below, when we consider low signal-
we focus on in this paper, and which are most promising to to-noise data for which it is difficult to measure the detailed
8 Ichikawa, Barkana, Iliev, Mellema, Shapiro

p
Figure 6. Our best-fitting GED model parameters TG (solid Figure 8. Standard deviation hTb2 i − hTb i2 as a function of the
curve), TL (long-dashed curve), σG (short-dashed curve), and cG cosmic mass-weighted ionization fraction. We show this quantity
(dotted curve, different y-axis range) as functions of the cosmic for the original simulation data (solid curves) and from our GED
mass-weighted ionization fraction. They are obtained by fitting to model fits (dashed curves). We consider the PDF in boxes of
the simulated PDFs for pixels of size 5 h−1 Mpc (bottom panel) size 5 h−1 Mpc, 10 h−1 Mpc and 20 h−1 Mpc (top to bottom, only
or 10 h−1 Mpc (top panel). simulation data for the 20 h−1 Mpc case).

shape of the PDF, and the variance is a major part of what

can be measured.

4 MONTE-CARLO RESULTS WITH ONE

FREE PARAMETER
In the rest of this paper, we present results for the expected
accuracy of reconstructing the 21-cm PDF itself and the cos-
mic reionization history from the PDF. To obtain these re-
sults, we assume that our simulation accurately reproduces
the real reionization process in the universe, and furthermore
we assume that our GED model introduced in the previous
section can be used as a substitute for the PDF from the sim-
ulation. In the future, more realistic simulations and more
elaborate PDF fits can be used instead, but the general idea
will be the same: as long as the overall signal-to-noise ratio
is low, it is essential to rely on simulations in order to both
reconstruct and interpret the observed PDF.
Of course, even if simulations perfectly predicted the
21-cm PDF for given inputs, various astrophysical scenar-
ios would give somewhat different ionizing source and sink
Figure 7. The derived probabilities PD (solid curve), PE (short- properties, and might yield a variety of possible PDFs. We
dashed curve) and PG (long-dashed curve) as functions of the leave the detailed exploration of this issue for future work,
cosmic mass-weighted ionization fraction. They are obtained by
and here assume that the simulated scenario matches reality,
fitting the GED model to the simulation PDFs for pixels of size
5 h−1 Mpc (bottom panel) or 10 h−1 Mpc (top panel). except that a small number of free parameters are allowed
to vary and must be reconstructed by trying to match the
observed PDF. In this section, we reconstruct reionization
from the PDF under the most optimistic assumption, where
we assume that the real PDF matches the simulated one as
a function of just a single parameter, the ionization frac-
tion x̄i . Thus, at each redshift, we find the value of x̄i that
 Measuring Reionization using the 21-cm PDF 9
best matches the observed PDF, assuming that the PDF pixels of comoving size rcom , we find
varies with x̄i as in the simulation. In practice we expect  −3  0.9
that x̄i is indeed the main parameter that determines the rcom 1+z
Np = 6.0 × 106 , (18)
PDF, but there should be some small additional dependence 5h−1 Mpc 10
on redshift and astrophysical inputs. In the next section we  −2.5  5.25
rcom 1+z
explore a more flexible approach which makes much weaker σN = 200 mK . (19)
5h−1 Mpc 10
assumptions.
Thus, here we wish to know how well a certain experi- In order to explore the dependence on the noise level, we also
ment can determine x̄i assuming this one-parameter model. consider various specifications with lower noise in the same
An experiment is specified by a total number of pixels Np field of view, e.g., 1/2 the noise we denote as MWA/2 (which
and a noise per pixel σN . We can simulate an observed PDF corresponds, e.g., to four-year data with the MWA), while
from such an experiment at a given input x̄i by generat- 1/10 the noise we denote as MWA/10 (which corresponds
ing Np data points from the PDF of equation (10) and to the regime of larger, second-generation 21-cm arrays).
adding to each noise generated from a Gaussian distribution Note that we include only Gaussian thermal noise, whose
with standard deviation σN . The resulting Monte-Carlo- magnitude is determined by the receiver’s system tempera-
generated “observed” PDF is then compared, via the C- ture, which in turn is set by the sky’s brightness tempera-
statistic of equation (6), to the model, which is equation (10) ture which is dominated by Galactic synchrotron emission
convolved with the Gaussian noise. This convolved function (Furlanetto et al. 2006). In particular this assumes perfect
q(Tb ) equals: foreground removal from the 21-cm maps; we leave an anal-
ysis of the effect of foreground residuals for future work.
q(Tb ) = PD G(Tb , σN ) + q1 (Tb ) + q2 (Tb ) , (15) We note the following conversions between comoving
distance and observational units of angular and frequency
where G is a Gaussian (eq. 2), and resolution:
−0.2 −0.5
1+z 1+z
 
5h−1 Mpc ≈ 2.6′
 
1 λ 2 σN
2
≈ 0.37 MHz .(20)
q1 (Tb ) = a exp + λTb × (16) 10 10
2 2
  2
  2
 The diffraction limit of the MWA is several arcminutes, but
λσN + Tb λσN − TL + Tb its frequency resolution will be around 10 kHz. In principle,
erf √ − erf √ ,
2σN 2σN this allows a measurement of the PDF in skinny boxes (thin-
ner in the redshift direction) rather than cubes. This would
 
1 σG (Tb − TG )2
q2 (Tb ) = cG exp − × (17) give us more points but with less signal in each, keeping the
2 σc 2σc2
 2 2
 overall signal-to-noise ratio about the same. By accessing
σN (TL − TG ) + σG (TL − Tb ) fluctuations on smaller scales, this skinny-box PDF would
erfc √ ,
2σc σG σN be somewhat broader than the cubic one but on the other
hand, our quantitative results for the toy model above sug-
where σc2 = σG 2
+ σN2
. As noted above, in this section we gest that decreasing the signal-to-noise ratio per pixel in
regard q(Tb ) as a one-parameter function of x̄i , taking TG , this way would have a strong tendency to introduce partial
σG , cG and TL to be functions of x̄i as shown in Figure 6. degeneracies. Thus, we do not expect this option to be pro-
For clarity we denote the input, real cosmic ionized fraction ductive (except in the cases when the errors in the cubic
simply x̄i , while the free parameter which is the output of PDF are very small), and focus here on the simplest case of
the fitting we denote x̄out
i . Note that we assume that the the 21-cm PDF measured in cubes.
experimental setup is sufficiently well characterized that σN At each redshift, we generate 1000 Monte Carlo in-
is known and need not be varied in the fitting. Also note stances of observed PDFs and minimize the C-statistic to
that while the various temperatures we have defined (Tb , TL , find the best-fitting model in each case. Results for MWA
and TG ) refer to the differential brightness temperature (i.e., and MWA/2 errors are plotted in Figure 9, which shows
0 mK refers to the absence of a cosmological signal), in prac- that for first-generation experiments the larger (10 h−1 Mpc)
tice, when the foregrounds as well as the cosmic microwave boxes are much more promising, since the lower noise σN (by
background are subtracted, the mean cosmological signal on a factor of ∼ 6) dominates despite the narrower PDF (com-
the sky will be removed as well, since there is no easy way pare Figs. 4 and 5) and smaller number of pixels Np (by
to separate out different contributions except through their a factor of 8). We note that lower noise is particularly im-
fluctuations. Thus, as in section 2.1, we assume that the portant in view of the partial degeneracy (demonstrated in
absolute Tb cannot be measured, and in our fitting always section 2.3 for the toy model) that arises when σN is greater
measure Tb with respect to its average value according to than the characteristic width of the intrinsic PDF. The par-
the PDF, both in each model and in each simulated data tial degeneracy is also apparent in comparing the MWA and
set. MWA/2 cases, where at some x̄i values, halving the errors
For the experimental specification, we adopt the crosses a degeneracy threshold and cuts the output uncer-
(rough) expected parameters for one-year observations of tainty in a non-linear fashion. We caution that cases that are
a single field of view with the MWA. We use the relations very near such a threshold may be susceptible to additional
for 21-cm arrays from the review by Furlanetto et al. (2006), numerical errors.
adopting a net integration time tint = 1000 hours, a collect- The same results are shown in Figure 10 in terms of
ing area Atot = 7 × 103 m2 , a field of view of π162 deg2 , and relative errors, making it easier to see and compare both
a total bandwidth ∆νtot = 6 MHz. Then assuming cubic small and large errors. Specifically, in terms of the various
10 Ichikawa, Barkana, Iliev, Mellema, Shapiro

Figure 9. Expected 1σ errors on reconstructing the cosmic mean Figure 10. Same as Figure 9 but showing relative errors (see
ionized fraction from the PDF, assuming just one free parameter. text), for better visibility of cases with small errors. We show
Specifically, for each input value of x̄i we show the output median f0 (absolute value shown, where negative values are open circles
(i.e., 50 percentile) x̄out
i as well as the 16 to 84 percentile range. and positive values are solid circles) , f+ (+ symbols), and f− (−
We consider MWA 1-yr errors (left panels) or MWA/2 (right pan- symbols). We consider MWA 1-yr errors (left panels) or MWA/2
els), for the PDF in 5 h−1 Mpc boxes (top panels) or 10 h−1 Mpc (right panels), for the PDF in 5 h−1 Mpc boxes (top panels) or
boxes (bottom panels). 10 h−1 Mpc boxes (bottom panels).

success in reconstructing the reionization history under the

percentile output ionization fractions (e.g., we denote the strict assumption of a single free parameter motivates us to
median by x̄out,50
i ), we show f0 = (x̄out,50
i /x̄i ) − 1 (the rel- consider in the following section a more flexible reconstruc-
ative difference between the median output value and the tion method.
true input value, representing the fractional bias of the re-
construction), f+ = (x̄out,84
i /x̄out,50
i ) − 1 (the relative dif-
ference between the 84% and median values, representive
5 MONTE-CARLO RESULTS WITH A
the fractional +1σ spread), and f− = 1 − (x̄out,16 i /x̄out,50
i )
FLEXIBLE FOUR-PARAMETER MODEL
(the relative difference between the 16% and median values,
representive the fractional −1σ spread). The Figure shows In the previous section we showed the expected accuracy of
that the reconstruction is typically unbiased within the er- reconstructing the cosmic reionization history from the 21-
rors (i.e., the 1σ range is significantly larger than the bias cm PDF, assuming the PDF shape is known as a function of
in the median), except for some points in the early stages of the cosmic mean ionized fraction. In this section we drop the
reionization. Only a little information is available with the latter assumption and present results for the expected accu-
PDF in the smaller boxes (except for a few redshifts with racy of reconstructing the detailed shape of the 21-cm PDF
MWA/2 errors); typically the error ranges are smaller near directly from the data. We focus on the regime of second-
the mid-point of reionization, partly due to the fact (see Fig- generation experiments, since the expected MWA errors do
ure 8) that the variance of the PDF suffices to distinguish not allow such a reconstruction. Even with the lower errors,
the mid-point of reionization from its two ends, but the early the PDF cannot be reconstructed parameter free, so we as-
and late stages are degenerate with each other in terms of sume that our four-parameter GED model from section 3.3
the variance. A rather good measurement of the reionization correctly describes the real intrinsic PDF (an assumption
history is expected with 10 h−1 Mpc boxes, in the mid to late which is explicitly true in our Monte-Carlo setup). Other-
stages of reionization, down to 1% errors in measuring the wise we do not assume any restrictions, and allow the four
cosmic mean ionized fraction (or even better with MWA/2 parameters of the model to vary freely when fitting (again
errors). When the errors are small, the measurement is un- by minimizing the C-statistic) to the noisy mock PDF data.
biased and has symmetric error bars. Specifically, we fit the four parameters TG , TL , σG , and
As shown in Figure 11, lower errors (approaching cG . We consider fitting the PDF in 5 h−1 Mpc boxes with
second-generation experiments) would avoid the degeneracy MWA/10 or MWA/20 errors. Figure 12 shows that signifi-
and allow a meaningful measurement of the cosmic reion- cant information can be reconstructed with MWA/10 errors,
ization history even with the PDF in the smaller boxes, but although the errors in the reconstructed parameters are usu-
10 h−1 Mpc boxes always give a more precisely measured out- ally fairly large (with particular failures at the early stage
put value by about an order of magnitude. The expected of reionization). The derived total probabilities of the GED
 Measuring Reionization using the 21-cm PDF 11

Figure 11. Same as Figure 10, but we consider MWA/4 errors Figure 12. Expected 1σ errors on reconstructing the PDF pa-
(left panels) or MWA/10 (right panels), for the PDF in 5 h−1 Mpc rameters assuming the four-parameter GED model, assuming
boxes (top panels) or 10 h−1 Mpc boxes (bottom panels). MWA/10 errors on the PDF in 5 h−1 Mpc boxes. We show the
16, 50, and 84 percentiles, as before, and also the assumed input
values (circles).
model components are shown in Figure 13; in particular,
the statistically significant measurement of the evolution of
PD (which is the cosmic ionized fraction smoothed over the
5 h−1 Mpc resolution) shows that significant information can
be extracted about the cosmic reionization history, even in
this more flexible fitting approach.
Since the errors on the reconstructed parameters with
MWA/10 noise are still mostly of order unity, we explored
further and found that MWA/20 is necessary to break most
of the degeneracies. Figure 14 shows that in this case the pa-
rameters can usually be reconstructed to 1 − 10% accuracy
(with symmetric error bars and insignificant bias). Specifi-
cally we show the four quantities PD , PG , TG and σG , which
together comprise a complete set that specifies the GED
model. Note that the measurement of PD is particularly pre-
cise, during the latter stages of reionization.
As in the previous section, the PDF in larger,
10 h−1 Mpc boxes, is easier to measure, due to the lower noise
per pixel. Thus, here we consider somewhat larger noise lev-
els, MWA/5 and MWA/10, with results shown in Figs. 15
and 16. Note that the last (highest x̄i ) point in TG is not
shown, since the input TG there is zero (see Figure 6), and
also, we show PD only during late reionization, where it is
non-zero (see Figure 7), and PE at earlier times. While the
errors are fairly large with MWA/5 errors, they reach the
Figure 13. Expected 1σ errors on reconstructing the derived
1 − 10% level with MWA/10, corresponding to a second- probabilities of the GED model, from a four-parameter fit to the
generation 21-cm experiment. PDF, assuming MWA/10 errors and 5 h−1 Mpc boxes. We show
the 16, 50, and 84 percentiles, as before, and also the assumed
input values (circles).
6 CONCLUSIONS
We have carried out a detailed quantitative analysis of
whether upcoming and future experiments can measure the
shape of the 21-cm PDF and derive from it the cosmic reion-
ization history. This is an important question since the PDF
12 Ichikawa, Barkana, Iliev, Mellema, Shapiro

Figure 14. Expected 1σ errors on reconstructing various quan- Figure 16. Expected 1σ errors on reconstructing various quan-
tities of the GED model, from a four-parameter fit to the PDF, tities of the GED model, from a four-parameter fit to the PDF,
assuming MWA/20 errors and 5 h−1 Mpc boxes. As in the previ- assuming MWA/10 errors and 10 h−1 Mpc boxes. We show the
ous section, we show the relative errors f0 (absolute value shown, relative errors f0 (absolute value shown, where negative values
where negative values are open circles and positive values are solid are open circles and positive values are solid circles) , f+ (+ sym-
circles) , f+ (+ symbols), and f− (− symbols). bols), and f− (− symbols).

during reionization is highly non-Gaussian, it directly pro-

vides important information such as the cosmic ionization
fraction at each redshift (though smoothed on the experi-
mental resolution scale), and is potentially a way to derive
the cosmic reionization history independently of the stan-
dard power spectrum analysis.
We developed a maximum-likelihood approach that
achieves maximum efficiency by minimizing the C-statistic
(eq. 6) applied to binned PDF data. We used a toy PDF
model of two Gaussians (eq. 1) to show that the simplistic
notion of signal-to-noise ratio (eq. 7) does not fully describe
the ability to extract the PDF out of noisy data. Instead,
once the noise per pixel rises above a few times the signal
(i.e., the width of the intrinsic PDF), the errors blow up due
to a strong degeneracy, even if the total signal-to-noise ratio
is kept fixed by increasing the number of pixels (Figs. 1 and
2).
We measured the 21-cm PDF as a function of redshift
in a large-scale N-body and radiative transfer simulation of
cosmic reionization (Figs. 4 and 5). The PDF starts out
close to Gaussian at high redshift, due to still-linear density
fluctuations, later develops an exponential tail at low Tb ,
Figure 15. Expected 1σ errors on reconstructing various quan- and finally becomes strongly peaked at zero towards the
tities of the GED model, from a four-parameter fit to the PDF, end of reionization. We empirically fit the PDF from the
assuming MWA/5 errors and 10 h−1 Mpc boxes. We show the rel- simulation with a four-parameter Gaussian + Exponential
ative errors f0 (absolute value shown, where negative values are + Delta function (GED) model (eq. 10, Figs. 6 and 7).
open circles and positive values are solid circles) , f+ (+ symbols), Assuming the simulations as a reliable guide for the
and f− (− symbols). evolution of the PDF, we quantitatively explored how well
parameters can be measured with two different approaches.
In the most optimistic approach, we assumed that the real
PDF matches the simulated one as a function of just a sin-
gle free parameter, the ionization fraction x̄i , and tried to
 Measuring Reionization using the 21-cm PDF 13
reconstruct this parameter from noisy mock data. We found Furlanetto S. R., Oh S. P., Briggs, F., 2006, Phys. Rep.,
that first-generation experiments (such as the MWA) are 433, 181
promising, at least if relatively large (10 h−1 Mpc) pixels are Harker, G. J. A., et al. 2009, MNRAS, 393, 1449
used along with their relatively low noise level per pixel. Iliev I. T., et al., 2006a, MNRAS, 371, 1057
Specifically, a rather good measurement of the reionization Iliev I. T., Mellema G., Pen U.-L., Merz H., Shapiro P. R.,
history is expected in the mid to late stages of reionization, Alvarez M. A., 2006b, MNRAS, 369, 1625
down to 1% errors in measuring the cosmic mean ionized Iliev I. T., Mellema G., Shapiro P. R., Pen U.-L. 2007,
fraction. MNRAS, 376, 534
We also considered reconstructing the cosmic reioniza- Iliev I. T., Mellema G., Pen U.-L., Bond J. R., Shapiro
tion history together with the PDF shape, all while assuming P. R. 2008, MNRAS, 384, 863
that the four-parameter GED model correctly describes the Iliev I. T., et al., 2009, arXiv0905.2920
real intrinsic PDF, but allowing the four parameters to vary Madau, P., Meiksin, A., & Rees, M. J. 1997, ApJ, 475, 429
freely when fitting mock data at each redshift. We found that McQuinn M., Furlanetto S. R., Hernquist L., Zahn O., Zal-
this flexible approach requires much lower noise levels, char- darriaga M., 2005, ApJ, 630, 643
acteristic of second-generation 21-cm experiments, to reach McQuinn M., Zahn O., Zaldarriaga M., Hernquist L.,
the level of 1 − 10% accuracy in measuring the parameters Furlanetto S. R., 2006, ApJ, 653, 815
of the 21-cm PDF. Mellema G., Iliev I. T., Alvarez M. A., Shapiro P. R., 2006a,
We note that cosmic reionization ends in our simulation New Astronomy, 11, 374
at redshift 7.5 (Fig. 3). If reionization in the real universe Mellema G., Iliev I. T., Pen U.-L., Shapiro P. R., 2006b,
ends later (e.g., closer to z = 6.5), then observations will MNRAS, 372, 679
be somewhat easier than we have assumed, due to the re- Merz H., Pen U.-L., Trac H., 2005, New Astronomy, 10,
duced foregrounds at lower redshift. On the simulation side, 393
further work is necessary to establish the numerical con- Oh S. P., Hansen M., Furlanetto S. R., Mesinger A., 2009,
vergence of the simulated 21-cm PDF during reionization, submitted
and to explore the dependence of the PDF on various astro- Saiyad-Ali, S., Bharadwaj, S., & Pandey, S. K., 2006, MN-
physical scenarios for the ionizing sources and sinks during RAS, 366, 213
reionization. This further effort is warranted since we have Santos M. G., Amblard A., Pritchard J., Trac H., Cen R.,
shown that the 21-cm PDF is a promising alternative to the Cooray A. 2008, ApJ, 689, 1
power spectrum which can independently probe the cosmic Spergel, D. N., et al., 2007, ApJS, 170, 377
reionization history. Wyithe S., Morales M., 2007, MNRAS, 379, 1647
Zahn O., Lidz A., McQuinn M., Dutta S., Hernquist L.,
Zaldarriaga M., Furlanetto S. R., 2007, ApJ, 654, 12

ACKNOWLEDGMENTS
We thank Oh, Hansen, Furlanetto, & Mesinger for provided
us with a draft of their paper far in advance of publica-
tion. RB is grateful for support from the ICRR in Tokyo,
Japan, the Moore Distinguished Scholar program at Cal-
tech, the John Simon Guggenheim Memorial Foundation,
and Israel Science Foundation grant 629/05. This study
was supported in part by Swiss National Science Foun-
dation grant 200021-116696/1, NSF grant AST 0708176,
NASA grants NNX07AH09G and NNG04G177G, Chandra
grant SAO TM8-9009X and Swedish Research Council grant
60336701.

REFERENCES
Barkana R., 2007, MNRAS, 376, 1784
Barkana, R., & Loeb, A. 2001, Phys. Rep., 349, 125
Barkana R., Loeb A., 2004, ApJ, 609, 474
Barkana R., Loeb A., 2007, Rep. Prog. Phys., 70, 627
Barkana, R., & Loeb, A. 2008, MNRAS, 384, 1069
Bharadwaj, S., & Pandey, S. K., 2005, MNRAS, 358, 968
Bowman J. D., Morales M. F., Hewitt J. N., 2006, ApJ,
638, 20
Cash W., 1979, ApJ, 228, 939
Ciardi B., Madau, P., 2003, ApJ, 596, 1
Furlanetto S. R., Zaldarriaga M., Hernquist L., 2004, ApJ,
613, 1