---
abstract: 'We introduce [<span style="font-variant:small-caps;">21CMMC</span>]{}: a parallelized, Monte Carlo Markov Chain analysis tool, incorporating the epoch of reionization ([EoR]{}) seminumerical simulation [<span style="font-variant:small-caps;">21CMFAST</span>]{}. [<span style="font-variant:small-caps;">21CMMC</span>]{} estimates astrophysical parameter constraints from 21 cm [EoR]{} experiments, accommodating a variety of [EoR]{} models, as well as priors on model parameters and the reionization history. To illustrate its utility, we consider two different [EoR]{} scenarios, one with a single population of galaxies (with a mass-independent ionizing efficiency) and a second, more general model with two different, feedback-regulated populations (each with mass-dependent ionizing efficiencies). As an example, combining three observations ($z=8$, 9 and 10) of the 21 cm power spectrum with a conservative noise estimate and uniform model priors, we find that interferometers with specifications like the Low Frequency Array/Hydrogen Epoch of Reionization Array (HERA)/Square Kilometre Array 1 (SKA1) can constrain common reionization parameters: the ionizing efficiency (or similarly the escape fraction), the mean free path of ionizing photons and the log of the minimum virial temperature of star-forming haloes to within 45.3/22.0/16.7, 33.5/18.4/17.8 and 6.3/3.3/2.4 per cent, $\sim1 \sigma$ fractional uncertainty, respectively. Instead, if we optimistically assume that we can perfectly characterize the [EoR]{} modelling uncertainties, we can improve on these constraints by up to a factor of $\sim$few. Similarly, the fractional uncertainty on the average neutral fraction can be constrained to within $\lesssim10$ per cent for HERA and SKA1. By studying the resulting impact on astrophysical constraints, [<span style="font-variant:small-caps;">21CMMC</span>]{} can be used to optimize (i) interferometer designs; (ii) foreground cleaning algorithms; (iii) observing strategies; (iv) alternative statistics characterizing the 21 cm signal; and (v) synergies with other observational programs.'
author:
- |
    Bradley Greig$^{1}$[^1] & Andrei Mesinger$^{1}$\
    $^1$Scuola Normale Superiore, Piazza dei Cavalieri 7, I-56126 Pisa, Italy
bibliography:
- '21cmMC.bib'
title: '21CMMC: an MCMC analysis tool enabling astrophysical parameter studies of the cosmic 21 cm signal'
---

galaxies: high-redshift – intergalactic medium – cosmology: theory – dark ages, reionization, first stars – diffuse radiation – early Universe

Introduction
============

The epoch of reionization ([EoR]{}) describes the period of enlightenment following the cosmic dark ages, when density perturbations grow into the first astrophysical sources (stars and galaxies). These sources produce ultraviolet (UV) ionizing photons, which escape into the intergalactic medium ([IGM]{}) eventually reionizing the pervasive neutral hydrogen fog.

The [EoR]{} is rich in astrophysical information, probing the formation and evolution of structure in the Universe, the nature of the first stars and galaxies and their impact on the [IGM]{} [see e.g. @Barkana:2007p2929; @Loeb:2013p2936; @Zaroubi:2013p2976]. Currently the astrophysics of the [EoR]{} are poorly understood. However, a wave of upcoming observations should produce a rapid advance in our understanding. At the forefront of these will be the detection of the redshifted 21 cm spin-flip transition of neutral hydrogen [see e.g. @Furlanetto:2006p209; @Morales:2010p1274; @Pritchard:2012p2958] from several dedicated radio interferometers.

The first-generation 21 cm experiments such as the Low Frequency Array (LOFAR; @vanHaarlem:2013p200 [@Yatawatta:2013p2980])[^2], the Murchison Wide Field Array (MWA; @Tingay:2013p2997)[^3] and the Precision Array for Probing the Epoch of Reionization (PAPER; @Parsons:2010p3000)[^4] are just now coming online. These have enabled progress in understanding and characterizing the relevant systematics and foregrounds, and may even yield a detection of the 21 cm power spectrum (PS) during the advanced stages of reionization [e.g. @Chapman:2012p344; @Mesinger:2014p244].

In the near future, second-generation experiments, the Square Kilometre Array (SKA; @Mellema:2013p2975)[^5] and the Hydrogen Epoch of Reionization Array (HERA; @Beardsley:2014p1529)[^6] will feature significant increases in collecting area and overall sensitivity. In addition to statistical measurements of the 21 cm PS, these second-generation experiments should provide the first tomographic maps of the 21 cm signal from the [EoR]{} [@Mellema:2013p2975], deepening our understanding of reionization-era physics.

However, we are faced with a fundamental question, *what exactly can we learn from these observations?* Qualitatively, we possess several key insights. Reionization driven by relatively massive galaxies should result in larger, more uniform ionized regions [e.g. @McQuinn:2007p1665; @Iliev:2012p112] and the large-scale PS should peak at the midpoint of the [EoR]{} [e.g. @Lidz:2008p1744], unless reionization is driven by X-rays [@Mesinger:2013p1835]. Abundant, absorption systems should result in an extended reionization [e.g. @Ciardi:2006p3030] characterized by small, disjoint ionized regions [e.g. @Alvarez:2012p1930]. Although these studies provide valuable intuition, they do not quantify the constraints and degeneracies among astrophysical parameters.

Therefore, it is imperative to develop an [EoR]{} analysis tool to quantitatively assess what astrophysical information can be gleaned from the cosmic 21 cm signal. This will serve to guide current and future 21 cm experiments in optimizing their observing strategies and interferometer designs to maximize their [EoR]{} scientific returns. In recent years, several authors have explored the benefits of this approach by performing [EoR]{} parameter searches. However, these approaches have been restricted to (i) sampling a finite, fixed grid of either simple analytic models (@Choudhury:2005p2859 [@Barkana:2009p116]) or 3D seminumerical simulations [@Mesinger:2012p1131; @Zahn:2012p1156; @Mesinger:2013p1835]; (ii) performing a Fisher Matrix analysis [@Pober:2014p35]; or (iii) Monte Carlo Markov Chain (MCMC) sampling only the reionization history in an analytic framework [@Harker:2012p2856; @Morandi:2012p2857; @Patil:2014p2858]

Here, we improve on these approaches by developing a public 21 cm [EoR]{} analysis tool, named [<span style="font-variant:small-caps;">21CMMC</span>]{}. [<span style="font-variant:small-caps;">21CMMC</span>]{} uses an optimized version of the well-tested, public simulation tool, [<span style="font-variant:small-caps;">21CMFAST</span>]{} [@Mesinger:2007p122; @Mesinger:2011p1123], to generate 3D tomographic maps of the 21 cm brightness temperature, statistically comparing these to a (mock) observation. It quantifies the constraints and degeneracies of astrophysical parameters governing the EoR, within a full Bayesian MCMC framework.

The remainder of this paper is organized as follows. In Section \[sec:21CMMC\], we outline the development and implementation of [<span style="font-variant:small-caps;">21CMMC</span>]{}. In Section \[sec:SGP\], we highlight its performance at providing astrophysical constraints from a popular theoretical model of the [EoR]{} containing a single population of ionizing galaxies with a mass-independent ionizing efficiency. To further emphasize the strength and flexibility of [<span style="font-variant:small-caps;">21CMMC</span>]{} in Section \[sec:DPL\], we then estimate the recovery of astrophysical parameters from a more generalized [EoR]{} model containing two ionizing galaxy populations characterized by their mass-dependent ionizing efficiency. Finally, in Section \[sec:conclusion\] we summarize the performance of [<span style="font-variant:small-caps;">21CMMC</span>]{} and finish with our closing remarks. Throughout this work, we adopt the standard set of $\Lambda$CDM cosmological parameters: ($\Omega_{\rm m}$, $\Omega_{\rm \Lambda}$, $\Omega_{\rm b}$, $n$, $\sigma_{8}$, $H_{0}$) =(0.27, 0.73, 0.046, 0.96, 0.82, 70 ${\rm km\, s^{-1} \, Mpc^{-1}}$), measured from the *Wilkinson Microwave Anisotropy Probe* [@Bennett:2013p3136] which are consistent with the latest results from Planck [@PlanckCollaboration:2014p3099].

[<span style="font-variant:small-caps;">21CMMC</span>]{} {#sec:21CMMC}
========================================================

Prior to showcasing the performance of [<span style="font-variant:small-caps;">21CMMC</span>]{}[^7], we first outline the various ingredients that contribute to the overall analysis tool. First, in Section \[sec:21CMFAST\], we summarize the included [EoR]{} simulation code [<span style="font-variant:small-caps;">21CMFAST</span>]{}. We then discuss the likelihood statistic in Section \[sec:likelihood\]. Next, in Section \[sec:MCMC\] we briefly describe the specifics of the MCMC driver and in Section \[sec:TNP\] we discuss the telescope sensitivities. Finally, in Section \[sec:pipeline\] we summarize the full [<span style="font-variant:small-caps;">21CMMC</span>]{} pipeline.

[<span style="font-variant:small-caps;">21CMFAST</span>]{} {#sec:21CMFAST}
----------------------------------------------------------

[<span style="font-variant:small-caps;">21CMMC</span>]{} employs a streamlined version of [<span style="font-variant:small-caps;">21CMFAST</span>]{}v1.1, a publicly available semi-numerical simulation code specifically designed for enabling astrophysical parameter searches. It employs approximate but efficient methods for [EoR]{} physics, and is accurate compared to computationally expensive radiative transfer (RT) simulations on scales relevant to 21 cm interferometry, $\geq 1$ Mpc [@Zahn:2011p1171].

Specifically, [<span style="font-variant:small-caps;">21CMFAST</span>]{} generates the [IGM]{} density, velocity, source and ionization fields by first creating a 3D realization of the linear density field within a cubic volume and then evolving the density field using the Zel’dovich approximation [@Zeldovich:1970p2023] before smoothing on to a lower resolution grid. Following this, [<span style="font-variant:small-caps;">21CMFAST</span>]{} estimates the ionization field by comparing the time-integrated number of ionizing photons to the number of baryons within spherical regions of decreasing radius, $R$, following the excursion-set approach described in @Furlanetto:2004p123. A cell at coordinates $(\boldsymbol{x},z)$ within the simulation volume is then tagged as fully ionized if, $$\begin{aligned}
\label{eq:ioncrit}
\zeta f_{\rm coll}(\boldsymbol{x},z,R,\bar{M}_{\rm min}) \geq 1,\end{aligned}$$ where $f_{\rm coll}(\boldsymbol{x},z,R,\bar{M}_{\rm min})$ is the fraction of collapsed matter within a spherical radius $R$ residing within haloes larger than $\bar{M}_{\rm min}$ [@Press:1974p2031; @Bond:1991p111; @Lacey:1993p115; @Sheth:1999p2053] and $\zeta$ is an ionizing efficiency describing the conversion of mass into ionizing photons (see Section \[sec:Zeta\]). Partial ionization is additionally included for cells not fully ionized by setting the ionized fraction of a cell smoothed at the minimum smoothing scale, $R_{\rm cell}$ to $\zeta f_{\rm coll}(\boldsymbol{x},z,R_{\rm cell},\bar{M}_{\rm min})$.

It is common to characterize the [EoR]{} with just three parameters: (i) the ionizing efficiency of high-$z$ galaxies; (ii) the mean free path of ionizing photons; (iii) the minimum virial temperature hosting star-forming galaxies. These parameters are somewhat empirical and in most cases are averaged over redshift and/or halo mass dependences. Nevertheless, they offer the flexibility to describe a wide variety of EoR signals, and can have a straightforward physical interpretation. Below, we summarize each of these three parameters and their physical origins.

### Ionizing efficiency, $\zeta$ {#sec:Zeta}

The ionizing efficiency of high-$z$ galaxies (equation \[eq:ioncrit\]) can be expressed as $$\begin{aligned}
 \label{eq:Zeta}
 \zeta = 30\left(\frac{f_{\rm esc}}{0.3}\right)\left(\frac{f_{\star}}{0.05}\right)\left(\frac{N_{\gamma}}{4000}\right)
 \left(\frac{2}{1+n_{\rm rec}}\right)
 \end{aligned}$$ where, $f_{\rm esc}$ is the fraction of ionizing photons escaping into the IGM, $f_{\star}$ is the fraction of galactic gas in stars, $N_{\gamma}$ is the number of ionizing photons produced per baryon in stars and $n_{\rm rec}$ is the typical number of times a hydrogen atom recombines. While our reionization model only depends on the product of equation (\[eq:Zeta\]), we provide reasonable estimates for the individual factors. The choice of $N_{\gamma}\approx4000$ is expected from Population II stars [e.g. @Barkana:2005p1934], while $f_{\star}$ and $f_{\rm esc}$ are extremely uncertain in high-$z$ galaxies [e.g. @Gnedin:2008p1935; @Wise:2009p1968; @Ferrara:2013p1248], though our choices are consistent with observations of galaxy luminosity functions [e.g. @Robertson:2013p2123; @Cooke:2014p2861; @Dayal:2014p2876]. Finally, $n_{\rm rec} \sim 1$ is consistent with the models of [@Sobacchi:2014p1157] which result in a ‘photon-starved’ end to reionization, consistent with emissivity estimates from the [Ly$\alpha$]{} forest [e.g. @Bolton:2007p3273; @McQuinn:2011p3293]. We allow $\zeta$ to vary within the range $\zeta\in[5,100]$. Since, $\zeta$ is empirically defined unlike $f_{\rm esc}$ which can be observationally constrained at lower redshift, we will often provide both when recovering our parameter constraints. For our range of $\zeta$, this corresponds to $f_{\rm esc}\in[0.05,1]$ using our fiducial values adopted in equation (\[eq:Zeta\]). In this work, we will first consider a constant $\zeta$ (Section \[sec:SGP\]) before considering a more generalized halo mass-dependent model in Section \[sec:DPL\].

### Mean free path of ionizing photons within ionized regions, $R_{\rm mfp}$

The propagation of ionizing photons through the [IGM]{} strongly depends on the abundances and properties of absorption systems (Lyman limit as well as more diffuse systems), which are below the resolution limits of EoR simulations. These systems act as photon sinks, roughly dictating the maximum scales to which [H[II]{}]{} bubbles can grow around ionizing galaxies. In EoR modelling this effect is usually parametrized with an implicit maximum horizon for ionizing photons (corresponding to the maximum filtering scale in excursion-set EoR models). Physically, this maximum horizon is taken to correspond to the mean free path of ionizing photons within ionized regions, $R_{\rm mfp}$ [e.g. @OMeara:2007p3360; @Prochaska:2009p3339; @Songaila:2010p3348; @McQuinn:2011p3293]. Motivated by recent sub-grid recombination models [@Sobacchi:2014p1157], here we explore the range $R_{\rm mfp}\in[5, 20]$ cMpc[^8]. While this range is relatively narrow, we point out that this parameter only becomes important once the size of the ionized regions become larger than $R_{\rm mfp}$ [e.g. @McQuinn:2007p1665; @Alvarez:2012p1930; @Mesinger:2012p1131]. Expanding the allowed range of $R_{\rm mfp}$ we consider in this work would therefore not greatly modify our conclusions.

### Minimum virial temperature of star-forming haloes, $T^{\rm Feed}_{\rm vir}$ {#sec:Mthresh}

Throughout, we choose to define the minimum threshold for a halo hosting a star-forming galaxy to be in terms of its virial temperature, which regulates processes important for star formation: gas accretion, cooling and retainment of supernovae outflows. The virial temperature is related to the halo mass via, [e.g. @Barkana:2001p1634] $$\begin{aligned}
M_{\rm min} &=& 10^{8} h^{-1} \left(\frac{\mu}{0.6}\right)^{-3/2}\left(\frac{\Omega_{\rm m}}{\Omega^{z}_{\rm m}}
\frac{\Delta_{\rm c}}{18\upi^{2}}\right)^{-1/2} \nonumber \\
& & \times \left(\frac{T_{\rm vir}}{1.98\times10^{4}~{\rm K}}\right)^{3/2}\left(\frac{1+z}{10}\right)^{-3/2}M_{\sun},\end{aligned}$$ where $\mu$ is the mean molecular weight, $\Omega^{z}_{\rm m} = \Omega_{\rm m}(1+z)^{3}/[\Omega_{\rm m}(1+z)^{3} + 
\Omega_{\Lambda}]$, and $\Delta_{\rm c} = 18\upi^{2} + 82d - 39d^{2}$ where $d = \Omega^{z}_{\rm m}-1$. Typically, $T_{\rm vir}\approx10^{4}$ K has been adopted in the literature as the minimum temperature when efficient atomic cooling occurs, corresponding to a DM halo mass of $\approx10^{8}M_{\sun}$ at $z\sim10$. In principle, this could be lower as the first stars are likely hosted within haloes with $M_{\rm halo}\approx 10^{6}$–$10^{7}$ $M_{\sun}$ [@Haiman:1996p2144; @Abel:2002p2149; @Bromm:2002p2153]. However, star formation within these haloes is likely inefficient (a few stars per halo) and can be quickly ($z>20$) suppressed by Lyman–Werner or other feedback processes well before the [EoR]{} [@Haiman:2000p2155; @Ricotti:2001p2160; @Haiman:2006p2169; @Mesinger:2006p2171; @Holzbauer:2012p2890; @Fialkov:2013p2903].

Here, we allow for a critical temperature threshold, $T^{\rm Feed}_{\rm vir}$, which can even be larger than $10^{4}$ K. Below $T^{\rm Feed}_{\rm vir}$ (but above $10^{4}$ K), haloes are still sufficiently large to produce low-mass galaxies; however, their potential wells are not deep enough to retain supernovae-driven outflows that quench star formation and their corresponding contribution to reionization [e.g. @Springel:2003p2176]. In this work, we allow $T^{\rm Feed}_{\rm vir}$ to vary within the range $T^{\rm Feed}_{\rm vir}\in[10^{4},2\times10^{5}]$ K, with the lower limit set by the atomic cooling threshold and the upper limit consistent with observed high-$z$ Lyman break galaxies (see for example Fig. \[fig:UVLum\]).

Likelihood statistic {#sec:likelihood}
--------------------

The offset of the 21 cm brightness temperature relative to the [CMB]{} temperature, $T_{\gamma}$ [e.g. @Furlanetto:2006p209], can be written as, $$\begin{aligned}
\delta T_{\rm b}(\nu) &=& \frac{T_{\rm S} - T_{\gamma}}{1+z}\left(1-{\rm e}^{-\tau_{\nu_0}}\right) \nonumber \\
&\approx& 27x_{{\ifmmode\mathrm{H\,{\scriptscriptstyle I}}\else{}H\,{\scriptsize I}\fi}{}}(1+\delta_{\rm nl})\left(\frac{H}{{\rm d}v_{\rm r}/{\rm d}r+H}\right)
\left(1 - \frac{T_{\gamma}}{T_{\rm S}}\right) \nonumber \\
& & \times \left(\frac{1+z}{10}\frac{0.15}{\Omega_{\rm m}h^{2}}\right)^{1/2}
\left(\frac{\Omega_{\rm b}h^{2}}{0.023}\right)~{\rm mK},\end{aligned}$$ where $T_{\rm S}$ is the gas spin temperature, $\tau_{\nu_0}$ is the optical depth at the 21 cm frequency, $\nu_0$, $\delta_{\rm nl}(\boldsymbol{x},z)$ is the evolved (Eularian) overdensity, $H(z)$ is the Hubble parameter, ${\rm d}v_{\rm r}/{\rm d}r$ is the gradient of the line of sight component of the velocity and all quantities are evaluated at redshift $z = \nu_{0}/\nu - 1$. Throughout this work, we include the effects of peculiar velocities and assume $T_{\rm S} \gg T_{\gamma}$ (likely appropriate for the advanced stages of reionization, corresponding to an average neutral fraction of $\bar{x}_{{\ifmmode\mathrm{H\,{\scriptscriptstyle I}}\else{}H\,{\scriptsize I}\fi}} \lesssim 0.8$, e.g. @Mesinger:2013p1835) as the explicit computation of the spin temperature within [<span style="font-variant:small-caps;">21CMFAST</span>]{} is computationally expensive; however, in future we will remove this assumption.

Although [<span style="font-variant:small-caps;">21CMFAST</span>]{} produces 3D maps of the 21 cm brightness temperature, we require a goodness of fit statistic in order to quantitatively compare an [EoR]{} realization with the mock observation. We choose the common spherically averaged 21 cm PS, $\Delta^{2}_{21}(k,z) \equiv k^{3}/(2\upi^{2}V)\,\delta \bar{T}^{2}_{\rm b}(z)\,\langle |\delta_{21}(\boldsymbol{k},z)|^{2}\rangle_{k}$ where $\delta_{21}(\boldsymbol{x},z) \equiv \delta T_{\rm b}(\boldsymbol{x},z)/\delta \bar{T}_{\rm b}(z) -1$, to be this statistic. It is important to note that while we restrict our default analysis to $\Delta^{2}_{21}(k,z)$, in principle any statistical measure of the cosmic 21 cm signal could be used in our framework.

The MCMC driver {#sec:MCMC}
---------------

An MCMC based [EoR]{} tool is a computationally challenging prospect since each step in the computation chain requires a [<span style="font-variant:small-caps;">21CMFAST</span>]{} realization. For reference, computing a 128$^{3}$ 21 cm realization at a single redshift takes $\sim3$s using a single core. This is efficient, however, it is a case of diminishing returns if we attempt to improve the computation time using the inbuilt multithreading of [<span style="font-variant:small-caps;">21CMFAST</span>]{}. Therefore, constructing an MCMC analysis tool which samples [<span style="font-variant:small-caps;">21CMFAST</span>]{} utilizing existing serial MCMC algorithms such as Metropolis–Hastings (@Metropolis:1953p2224 [@Hastings:1970p97]), which require multiple long computation chains, would be too computationally inefficient.

In recent years, several new, more efficient and parallel MCMC algorithms have been developed. Once such algorithm is the affine invariant ensemble sampler developed by @Goodman:2010p843. This algorithm has since been modified and improved before being released as the publicly available <span style="font-variant:small-caps;">python</span> module [<span style="font-variant:small-caps;">emcee</span>]{} by @ForemanMackey:2013p823. Recently, @Akeret:2012p842 discussed the virtues of adopting these computational algorithms for applications of cosmological parameter sampling. These authors released an easy to implement, massively parallel, publicly available <span style="font-variant:small-caps;">python</span> module [<span style="font-variant:small-caps;">cosmohammer</span>]{} which has already been extensibly adopted for a wide variety of cosmological application. We therefore choose to build [<span style="font-variant:small-caps;">21CMMC</span>]{} using a modified version of [<span style="font-variant:small-caps;">cosmohammer</span>]{} which performs our MCMC reionization parameter sampling by calling [<span style="font-variant:small-caps;">21CMFAST</span>]{} at each step in the computation chain to compute the log-likelihood.

Telescope noise profiles {#sec:TNP}
------------------------

To illustrate the performance of [<span style="font-variant:small-caps;">21CMMC</span>]{} for current and future 21 cm experiments, we compute the telescope sensitivities using the <span style="font-variant:small-caps;">python</span> module [[21]{}<span style="font-variant:small-caps;">cmsense</span>]{}[^9][@Pober:2013p41; @Pober:2014p35]. In this section, we outline only the specifics and assumptions we make in producing our realistic telescope noise profiles, summarizing the key telescope parameters in Table \[tab:TelescopeParams\], and defer the reader to the more detailed discussions within @Parsons:2012p95 and @Pober:2013p41 [@Pober:2014p35].

The thermal noise PS is computed at each $uv$-cell according to the following [e.g. @Morales:2005p1474; @McQuinn:2006p109; @Pober:2014p35], $$\begin{aligned}
 \label{eq:NoisePS}
\Delta^{2}_{\rm N}(k) \approx X^{2}Y\frac{k^{3}}{2\upi^{2}}\frac{\Omega^{\prime}}{2t}T^{2}_{\rm sys},\end{aligned}$$ where $X^{2}Y$ is a cosmological conversion factor between observing bandwidth, frequency and comoving distance units, $\Omega^{\prime}$ is a beam-dependent factor derived in @Parsons:2014p781, $t$ is the total time spent by all baselines within a particular $k$ mode and $T_{\rm sys}$ is the system temperature, the sum of the receiver temperature, $T_{\rm rec}$, and the sky temperature $T_{\rm sky}$. For all telescope configurations considered in this work, we assume drift scanning with a total synthesis time of 6 h per night. While HERA can only operate in drift scan mode, LOFAR and SKA can perform tracked scans. In the future, we will investigate various observing strategies; however, drift scanning suffices here for our like-to-like comparisons.

Of great concern for 21 cm experiments is the estimation and removal of the bright foreground emission. Most notable of these are the chromatic effects due to how the interferometric array’s $uv$ coverage depends on frequency. It has been shown that these effects do not lead to mode-mixing across all frequencies, but instead localize the majority of the foreground noise into a foreground ‘wedge’ within cylindrical 2D $k$-space [@Datta:2010p2792; @Morales:2012p2828; @Parsons:2012p2833; @Trott:2012p2834; @Vedantham:2012p2801; @Thyagarajan:2013p2851; @Liu:2014p3465; @Liu:2014p3466]. This results in a relatively pristine observing window where the cosmological 21 cm signal is dominated only by thermal noise. However, it is still uncertain where to define the transition from the foreground ‘wedge’ to the [EoR]{} observing window.

Moreover, the summation over redundant baselines within an array configuration can reduce uncertainties in the thermal noise estimates [@Parsons:2012p95]. However, the gain depends on whether they are perfectly or partially coherent [@Hazelton:2013p1481].

In @Pober:2014p35, these authors consider numerous foreground removal strategies, including a variety of treatments for the summation over redundant baselines and how to define the transition from the foreground ‘wedge’. We defer the reader to their work for more detailed discussions, and highlight that throughout this work, we adopt what they refer to as their ‘moderate’ foreground removal. For this, we assume a fixed location of the foreground ‘wedge’ to extend $\Delta k_{\parallel} = 0.1 \,h$ Mpc$^{-1}$ beyond the horizon limit [@Pober:2014p35] and compute our estimates of the thermal noise PS sampling only $k$-modes that fall within the [EoR]{} observing window. This strategy further assumes the coherent addition of all redundant baselines.

  Parameter                  LOFAR   HERA             SKA           
  ------------------------- ------- ------- ----------------------- --
  Telescope antennas          48      331             866           
  Diameter (m)               30.8     14              35            
  Collecting area (m$^2$)    35762   50953          833190          
  $T_{\rm rec}$ (K)           140     100    0.1T$_{\rm sky}$ + 40  
  Bandwidth (MHz)              8       8               8            
  Integration time (hrs)     1000    1000            1000           

  : Summary of telescope parameters we use to compute sensitivity profiles (see the text for further details).[]{data-label="tab:TelescopeParams"}

We include another level of conservatism to our errors, applying a strict $k$-mode cut to our 21 cm PS. We sample only modes above $|k|$ = 0.15 Mpc$^{-1}$. Following the moderate foreground removal strategy of @Pober:2014p35, it effectively turns out that there is limited information available below this $k$-mode; however, we still apply this strict cut on top of the estimated total noise PS. Modes with $|k| \geq$ 0.15 Mpc$^{-1}$ can still be affected by the foreground ‘wedge’ due to its spherically asymmetric structure, hence we apply both cuts.

Since equation (\[eq:NoisePS\]) is computed for each individual $k$-mode, the sample variance of the cosmological PS can be easily included to produce an estimate of the total noise by performing an inverse-weighted summation over all the individual modes [@Pober:2013p41], $$\begin{aligned}
 \label{eq:T+S}
\delta\Delta^{2}_{\rm T+S}(k) = \left(\sum_{i}\frac{1}{(\Delta^{2}_{{\rm N},i}(k) + \Delta^{2}_{21}(k))^{2}}\right)^{-1/2},\end{aligned}$$ where $\delta\Delta^{2}_{\rm T+S}(k)$ is the total uncertainty from thermal noise and sample variance in a given $k$-mode and $\Delta^{2}_{21}(k)$ is the cosmological 21 cm PS. Here, we assume Gaussian errors for the cosmic-variance term, which is a good approximation on large scales.

Within this work, we restrict our analysis to three telescope designs (LOFAR, HERA and the SKA). In @Pober:2014p35, these authors qualitatively compare the merits of a variety of telescope designs through a measurement of the overall sensitivity of the 21 cm PS. To perform this direct comparison, these authors set a variety of design specific parameters to be equivalent. Instead, we perform our comparison using the exact design specifics of each experiment as outlined below. Note, the difference is marginal, however using system specific parameters we better mimic the expected overall sensitivities of the modelled telescope. For all, we model $T_{\rm sky}$ using the frequency dependent scaling $T_{\rm sky} = 60\left(\frac{\nu}{300~{\rm MHz}}\right)^{-2.55}~{\rm K}$ [@Thompson2007].

-   *LOFAR:* we model LOFAR using the antennas positions listed in @vanHaarlem:2013p200. Following the approach of @Pober:2014p35, we focus solely on the Netherlands core for the purposes of detecting the 21 cm PS assuming that all 48 high-band antennas (HBA) stations can be correlated separately maximizing the total number of redundant baselines. We assume $T_{\rm rec} =140$ K consistent with @Jensen:2013p1389.\

-   *HERA:* we follow the design specifics outlined in @Beardsley:2014p1529, where they propose a final design including 331 antennas. Each antennas is 14m in diameter closely packed into a hexagonal configuration to maximize the total number of redundant baselines [@Parsons:2012p95]. While the total observing area for HERA is an order of magnitude lower than the SKA, its design is optimized specifically for PS measurements allowing for comparable sensitivity [@Pober:2014p35][^10]. However, our 331 antennas telescope is smaller than the larger 547 antennas instrument proposed by these authors, therefore, our design performs more as an intermediate instrument. For HERA, we model the total system temperature as $T_{\rm sys} = 100 + T_{\rm sky}~{\rm K}$.\

-   *SKA:* we mimic the current design brief for the SKA-low Phase 1 instrument outlined in the SKA System Baseline Design document[^11]. Specifically, SKA-low Phase 1 includes a total of 911 35 m antennas stations, with 866 stations normally distributed in radius within a core extending out to a radius of $\sim$2 km and a further 45 stations arranged into three spiral arms each with 15 stations extending to $\sim50$ km. For our SKA design we include only the core stations, as the spiral arms add little to the PS sensitivity. As in @Pober:2014p35, we model the central core of SKA as a hexagon to maximize baseline redundancy, however our core extends out to a radius of $\sim300$m which includes a total of 218 telescope antennas. We then randomly place the telescope antennas normalising to ensure we contain 433 and 650 telescopes within a radii of 600 and 1000 m and randomly place the remaining 216 antennas at a radii larger than 1000 m. The total SKA system temperature is modelled as outlined in the SKA System Baseline Design, $T_{\rm sys} = 1.1T_{\rm sky} + 40~{\rm K}$.

The [<span style="font-variant:small-caps;">21CMFAST</span>]{} pipeline {#sec:pipeline}
-----------------------------------------------------------------------

Eventually [<span style="font-variant:small-caps;">21CMMC</span>]{} will be performed on actual observations. However for this proof-of-concept, we first create mock observations of the cosmic [EoR]{} signal and generate sensitivities for various instruments[^12]. Using [<span style="font-variant:small-caps;">21CMFAST</span>]{}, we perform a high-resolution, large-box realization of our fiducial [EoR]{} model, simulated within a 500$^{3}$ cMpc$^{3}$ box with an initial grid of 1536$^{3}$ voxels smoothed on to a 256$^{3}$ grid. We then calculate the 21 cm PS from this simulation at several different redshifts and compute the total noise contribution from thermal and sample variance using [[21]{}<span style="font-variant:small-caps;">cmsense</span>]{}. In addition to the thermal noise and sample variance, we add an additional uncertainty due to EoR modelling, accounting for, e.g. differences in the implementation of RT [e.g. @Iliev:2006p3393] and errors in semi-numerical approximations [e.g. @Zahn:2011p1171]. In other words, even if we have a well behaved Universe characterized exactly by our 3 or 5 parameter models, the simulated PS will be different from the true one due to the numerical implementation. We assume for this a constant multiplicative error of 25 per cent, guided by the observed 10–30 per cent difference at the relevant scales between [<span style="font-variant:small-caps;">21CMFAST</span>]{} and the cosmological RT simulations in @Zahn:2011p1171. For each mock observation of the 21 cm PS, we add this 25 per cent error in quadrature with the noise contribution from equation (\[eq:T+S\]). In principle this would be a correlated error, however, at this stage a constant multiplicative factor is sufficient. Below, we include some results excluding this error, showing the maximum information which can be extracted from the signal assuming a perfect characterization of modelling uncertainties.

To sample the reionization parameter space in [<span style="font-variant:small-caps;">21CMMC</span>]{}, we perform a smaller, lower resolution [EoR]{} simulation. Using a different set of initial conditions to the mock observations, we simulate a 250$^{3}$ cMpc$^{3}$ volume with an initial grid of 768$^{3}$ voxels smoothed down on to a 128$^{3}$ grid. Our two simulations are both sampled at the same physical resolution ($\sim2$ Mpc), and resampled by the same factor of 6 from the initial density grid to ensure convergence between the recovered 21 cm PS from the mock observation and the PS from the smaller 128$^{3}$ simulations. To account for the impact of shot noise on the smallest spatial modes (largest $k$) within these low-resolution simulations, we apply an additional cut at $k\geq1.0~{\rm Mpc}^{-1}$.

At each step in the MCMC algorithm we compute the 21 cm PS corresponding to that point in our astrophysical parameter space using the 3D realization from [<span style="font-variant:small-caps;">21CMFAST</span>]{}. We compute its likelihood using the $\chi^{2}$ statistic between $k=0.15$ and 1.0 Mpc$^{-1}$ to assess whether to accept or reject the new proposed set of [EoR]{} parameters.

Finally, we briefly discuss the improvements in our method to sample the astrophysical [EoR]{} parameter space, relative to previous works. These previous approaches have been restricted to either sampling a fixed, course grid of reionization models [e.g. @Barkana:2009p116; @Mesinger:2012p1131; @Zahn:2012p1156] or the fundamental assumption of Gaussian errors in Fisher Matrix applications (@Pober:2014p35; although new methods exist to overcome this assumption, e.g. @Joachimi:2011p3193 [@Sellentin:2014p3227]). Instead, with a Bayesian MCMC framework, we are able to more efficiently sample the [EoR]{} parameter space and by construction obtain direct estimates of the individual probability distribution functions (PDFs) of the [EoR]{} parameters, removing the *a priori* assumptions on the intrinsic probability distribution (i.e. Gaussian). Moreover, by sampling over 10$^{5}$ realizations within [<span style="font-variant:small-caps;">21CMMC</span>]{} this allows us to sample many features within the [EoR]{} parameter space not visible in course-grid models. Additionally, this allows for the better characterization of parameter degeneracies, which is a problem for Fisher Matrix analyses.

Single Ionizing Galaxy Population {#sec:SGP}
=================================

To illustrate the use of [<span style="font-variant:small-caps;">21CMMC</span>]{}, we first adopt the popular three-parameter EoR model, discussed above. This model is characterized by a single population of efficient star-forming galaxies hosted by haloes above our defined critical temperature $T^{\rm Feed}_{\rm vir}$. The ionizing efficiency of this model is a step function, $$\begin{aligned}
 \label{eq:zeta_0}
\zeta(T_{\rm vir}) = \left\{ \begin{array}{cl} 
      \zeta_{0}  & {T_{\rm vir} \geq T^{\rm Feed}_{\rm vir}}\\
      0 & {T_{\rm vir} < T^{\rm Feed}_{\rm vir}}.
\end{array} \right.\end{aligned}$$ In this model, a constant fraction of the host halo mass goes into the production of ionizing photons. To illustrate the performance of [<span style="font-variant:small-caps;">21CMMC</span>]{}, we choose $\zeta_{0} = 30$, $R_{\rm mfp} = 15$ Mpc and $T^{\rm Feed}_{\rm vir} = 3\times10^{4}$ K for our mock observation. Our chosen values are roughly consistent with observational data, though they merely serve for illustrative purposes.

Building physical intuition {#sec:intuition}
---------------------------

### UV luminosity function {#sec:UVLum}

It is first constructive to illustrate how this [EoR]{} model would appear with respect to current observations. In Fig. \[fig:UVLum\], we provide a schematic picture of our single ionizing galaxy population compared to the observed non-ionizing ($\sim$ 1500 Å) UV luminosity function at $z=8$ of @Bouwens:2014p1528. To obtain the UV luminosity of our galaxy population, we use the Schechter function parameter values described in @Kuhlen:2012p1506 and abundance matching to estimate the host halo mass as a function of the UV magnitude, $M(M_{\rm UV})$. Fitting this expression as a power-law, we convert our $T^{\rm Feed}_{\rm vir}$ into a UV magnitude which produces a sharp drop in the UV luminosity function at $M_{\rm UV} = -12.6$ ($M\sim10^{8.78}\,M_{\sun}$). By construction, the non-ionizing UV luminosity function matches the observational data and we extrapolate down into the faint-end until $T^{\rm Feed}_{\rm vir}$, below which no ionizing photons are produced/escape the galaxy. In our illustrative model, this drop, tracing fundamental galaxy formation physics, is well beyond the sensitivity limits even for JWST [e.g. @Salvaterra:2011p2908], showcasing the power of our approach.

The UV luminosity function with an appropriate dust model can constrain the star formation rate, which (approximately) collapses the uncertainty in the ionizing emissivity to just that in the escape fraction, $f_{\rm esc}$. In other words, the total ionizing luminosity, $L_{\rm ion}$ is roughly proportional to the product of the non-ionzing UV luminosity ($L_{\rm UV}$) and $f_{\rm esc}$: $L_{\rm ion}(M) \propto f_{\rm esc}(M)L_{\rm UV}(M)$. Our fiducial choice of a constant $\zeta$ implies $L_{\rm ion}\propto M$, from which we can obtain an expression for $L_{\rm UV}(M)$ using our power-law expression for $M(M_{\rm UV})$ to recover $L_{\rm UV}(M)\propto M^{0.68}$ for the faint end of the luminosity function. Therefore, our illustrative, single population model at $z=8$ roughly translates to $f_{\rm esc}\propto M^{0.32}$.

### How does the 21 cm PS depend on the [EoR]{} parameters? {#sec:variation}

In Fig. \[fig:3param\_variation\], we highlight how each individual [EoR]{} parameter affects the 21 cm PS, and also provide a visual representation of the total noise profiles for each telescope. In all panels, we show the mock $z=9$ observation for our fiducial [EoR]{} parameters (solid black curve) which results in $\bar{x}_{{\ifmmode\mathrm{H\,{\scriptscriptstyle I}}\else{}H\,{\scriptsize I}\fi}{}} = 0.71$. In order to facilitate the discussion, we decompose the 21 cm PS to first order (neglecting the contribution from redshift-space distortions and spin temperature fluctuations), $$\begin{aligned}
 \label{eq:decomp}
\Delta^{2}_{21}(k) &\propto& \bar{x}^{2}_{{\ifmmode\mathrm{H\,{\scriptscriptstyle I}}\else{}H\,{\scriptsize I}\fi}{}}\Delta^{2}_{\delta\delta}(k) -2\bar{x}_{{\ifmmode\mathrm{H\,{\scriptscriptstyle I}}\else{}H\,{\scriptsize I}\fi}}\left(1-\bar{x}_{{\ifmmode\mathrm{H\,{\scriptscriptstyle I}}\else{}H\,{\scriptsize I}\fi}}\right)\Delta^{2}_{\delta {\rm x}}(k) \nonumber \\
& & + (1-\bar{x}_{{\ifmmode\mathrm{H\,{\scriptscriptstyle I}}\else{}H\,{\scriptsize I}\fi}})^{2}\Delta^{2}_{\rm xx}(k),\end{aligned}$$ where $\Delta^{2}_{\delta\delta}$, $\Delta^{2}_{\delta {\rm x}}$ and $\Delta^{2}_{\rm xx}$ are the matter density, density-ionization and ionization PS respectively and fluctuations in the ionized fraction are defined as $\delta_{\rm x} = (\bar{x}_{{\ifmmode\mathrm{H\,{\scriptscriptstyle II}}\else{}H\,{\scriptsize II}\fi}{}} - \langle \bar{x}_{{\ifmmode\mathrm{H\,{\scriptscriptstyle II}}\else{}H\,{\scriptsize II}\fi}{}} \rangle)/\langle \bar{x}_{{\ifmmode\mathrm{H\,{\scriptscriptstyle II}}\else{}H\,{\scriptsize II}\fi}{}} \rangle)$. Although @Lidz:2007p1929 showed that this first-order approximation provides a rather poor description for the 21 cm PS, for our qualitative discussion below equation (\[eq:decomp\]) sufficiently highlights the major components to the total power.

Before reionization gets well underway ($\bar{x}_{{\ifmmode\mathrm{H\,{\scriptscriptstyle I}}\else{}H\,{\scriptsize I}\fi}{}} > 0.9$) $\Delta^{2}_{21}$ closely follows the matter density PS, $\Delta^{2}_{\delta\delta}$. Initially, cosmic [H[II]{}]{} bubbles grow around biased density peaks hosting galaxies. This causes the large-scale power to drop due to the increased (negative) contribution from the cross-correlation term, $\Delta^{2}_{\delta {\rm x}}$. As time evolves ($\bar{x}_{{\ifmmode\mathrm{H\,{\scriptscriptstyle I}}\else{}H\,{\scriptsize I}\fi}{}}$ decreases), the [H[II]{}]{} bubbles become more numerous (with the increasing population of ionizing sources) and begin to overlap, increasing the contribution from the ionization PS, $\Delta^{2}_{\rm xx}$. Eventually, this term begins to dominate, resulting in an overall increase in the total power in $\Delta^{2}_{21}$ across all scales. At $\bar{x}_{{\ifmmode\mathrm{H\,{\scriptscriptstyle I}}\else{}H\,{\scriptsize I}\fi}{}} = 0.5$, the large-scale power peaks. As the [H[II]{}]{} bubbles grow and overlap, the peak power in $\Delta^{2}_{21}$ shifts to larger physical scales (smaller $k$) as $\Delta^{2}_{\rm xx}$ is sensitive to the physical size of the [H[II]{}]{} regions. At the same time, the power on small scales begins to drop as the number of small, isolated [H[II]{}]{} regions decreases.

This generic PS evolution during the EoR is evident in the first three panels of Fig. \[fig:3param\_variation\] (see also, e.g. @McQuinn:2007p1665 [@Lidz:2008p1744; @Friedrich:2011p1816]). A larger value of $\zeta_{0}$ corresponds to a more advanced EoR at $z=9$, as galaxies are more efficient at producing ionizing photons. Similarly, a smaller value of $T^{\rm Feed}_{\rm vir}$ results in a more advanced EoR, as smaller, more abundant galaxies can efficiently form stars. The progress of the EoR is less sensitive to variation in $R_{\rm mfp}$ (second panel); instead this maximum horizon for the ionizing photons produces a prominent ‘knee’ in the PS, suppressing coherent ionization structure on large scales (e.g. @Alvarez:2012p1930 [@Mesinger:2012p1131]).

The PS is not however entirely determined by $\bar{x}_{{\ifmmode\mathrm{H\,{\scriptscriptstyle I}}\else{}H\,{\scriptsize I}\fi}{}}$. In the final panel of Fig. \[fig:3param\_variation\] we illustrate the impact of our [EoR]{} parameters at a fixed [IGM]{} neutral fraction of $\bar{x}_{{\ifmmode\mathrm{H\,{\scriptscriptstyle I}}\else{}H\,{\scriptsize I}\fi}{}} = 0.71$ (i.e. all at the equivalent stage of reionization). For example, if the EoR is driven by more biased, larger mass haloes (larger $T^{\rm Feed}_{\rm vir}$), their relative dearth must be compensated for by a higher ionizing efficiency. This results in an EoR morphology with larger, more isolated cosmic [H[II]{}]{} patches [e.g. @McQuinn:2007p1665]. $\Delta^{2}_{\rm xx}$ is more peaked on large scales, with a decrease in $\Delta^{2}_{\delta {\rm x}}$. As a result, the 21 cm PS is higher, with a more prominent large-scale ‘bump’ tracing the typical [H[II]{}]{} region size. Also explicitly evident in the last panel is the prominent ‘knee’ feature, imprinted by small values of $R_{\rm mfp}$. Here, for $R_{\rm mfp} = 3$ Mpc (blue curve), the ionization PS peaks on much smaller scales, with a sharp drop in large-scale power compared to the other curves computed with $R_{\rm mfp} = 15$ Mpc. Therefore, the 21 cm PS contains information on [*both*]{} the EoR epoch (i.e. $\bar{x}_{{\ifmmode\mathrm{H\,{\scriptscriptstyle I}}\else{}H\,{\scriptsize I}\fi}{}}$), as well as [*independent*]{} (i.e. at fixed $\bar{x}_{{\ifmmode\mathrm{H\,{\scriptscriptstyle I}}\else{}H\,{\scriptsize I}\fi}{}}$) information on the properties of galaxies and the IGM.

In all panels, we provide the estimated 1$\sigma$ errors on the 21 cm PS for each of the considered 21 cm experiments, including the additional 25 per cent modelling uncertainty. It is immediately obvious that LOFAR will have difficulty discriminating among various models, having to rely on the largest scales close to the foreground-dominated regime. HERA and the SKA on the other hand should be able to recover meaningful constraints on [EoR]{} parameters, with the large-scales limited by the additional 25 per cent modelling uncertainty we include. On smaller spatial scales, SKA performs considerably better than HERA.

Single epoch observation of the 21 cm PS {#sec:HERA_SKA}
----------------------------------------

21 cm experiments have coverage over a large bandwidth (e.g. 50–350 MHz for the SKA allowing for observations to $z\lesssim28$). However, foregrounds and high data rates can limit the coverage to narrower instantaneous bandwidths. Hence, we begin by analysing an observation at a single redshift (assuming our fiducial bandwidth of 8 MHz), before moving on to a broader bandpass observation. We restrict our analysis to only second-generation 21 cm experiments, HERA and the SKA, since with a single bandwidth the first-generation instruments will only achieve a marginal detection at best (e.g. @Mesinger:2014p244 [@Pober:2014p35]; see also Fig. \[fig:3param\_variation\]).

In Fig. \[fig:2D\_HERA\_SKA\_singlez\], we compare the outputs of [<span style="font-variant:small-caps;">21CMMC</span>]{} for both HERA (red curve) and SKA (blue curve) for an assumed single $z=9$ observation of our fiducial 21 cm PS with a total integration time of 1000 h. In this figure, and for the remainder of this work we assume uniform priors on all recovered parameters within their allowed ranges outlined in Section \[sec:21CMFAST\]. The dashed vertical lines denote the mock observation values. Across the diagonal panels of this figure we provide the 1D marginalized PDFs for each of our [EoR]{} parameters. In the top-right panel, we provide the marginalized 1D PDFs of the [IGM]{} neutral fraction. Additionally, we choose to renormalize all 1D PDFs to have the peak probability equal to unity to better visually emphasize the shape and width of the recovered distribution. In the lower-left corner of Fig. \[fig:2D\_HERA\_SKA\_singlez\], we show the 2D joint marginalized likelihood distributions. In each panel, the cross denotes the location of our fiducial [EoR]{} parameters. For each, we construct both the 1$\sigma$ and 2$\sigma$ contours (thick and thin lines, respectively), by computing a smoothed 2D histogram of the entire MCMC sample and estimating the likelihood contours which enclose 68 and 95 per cent of the data sample, respectively.

For all [EoR]{} parameters the recovered PDFs are centred around their input fiducial parameters, highlighting the strong performance of [<span style="font-variant:small-caps;">21CMMC</span>]{} at recovering our input model. This is despite the asymmetric distributions of the recovered PDFs, emphasizing the strength of the MCMC approach within [<span style="font-variant:small-caps;">21CMMC</span>]{}. Comparing these two 21 cm experiments, we note the significant improvement achievable with the increased sensitivity of the SKA. This can easily be seen in the case of both $\zeta_{0}$ and $T^{\rm Feed}_{\rm vir}$, where the widths of the 1D PDFs, as well as the 2D joint likelihood contours, are noticeably tighter for SKA than for HERA.

The improved constraints from the SKA relative to HERA can be best explained by referring to Fig. \[fig:3param\_variation\]. We observe that while the sensitivity at large scales for both SKA and HERA are comparable, only the SKA is capable of probing the available small-scale information. Given the degenerate nature of the [EoR]{} parameters, by accessing this additional shape information tighter constraints are achievable. While it is not sufficient to break the degeneracies it does reduce their amplitude, substantially limiting the allowed [EoR]{} parameter space. It is important to note that with only a single redshift observation, the uncertainty is large and the likelihood is non-Gaussian distributed. Nevertheless, our Bayesian MCMC framework captures parameter constraints and degeneracies, recovering the actual shape of the likelihood distribution (without the Gaussian assumptions of the Fisher matrix formalism; @Pober:2014p35).

  ----------------------- ------------------------------------ ---------------------------------- -------------------------------------- -- ---------------------------------------------------------------------------------------- -- --
  Instrument                                                               Parameter                                                         $\bar{x}_{{\ifmmode\mathrm{H\,{\scriptscriptstyle I}}\else{}H\,{\scriptsize I}\fi}{}}$     
  (single-$z$)                        $\zeta_{0}$                        $R_{\rm mfp}$             log$_{10}$$(T^{\rm Feed}_{\rm vir})$                                             $z = 9$                                             
  HERA 331                 $38.38\substack{+22.43 \\ -11.33}$   $14.34\substack{+3.72 \\ -4.95}$     $4.57\substack{+0.36 \\ -0.32}$                                    $0.73\substack{+0.10 \\ -0.21}$                                 
  HERA 331 (with prior)    $38.37\substack{+20.86 \\ -11.53}$   $13.88\substack{+3.76 \\ -4.39}$     $4.65\substack{+0.32 \\ -0.35}$                                    $0.75\substack{+0.09 \\ -0.13}$                                 
  SKA                       $32.11\substack{+8.32 \\ -6.47}$    $13.78\substack{+3.87 \\ -4.09}$     $4.50\substack{+0.20 \\ -0.21}$                                    $0.71\substack{+0.07 \\ -0.09}$                                 
  ----------------------- ------------------------------------ ---------------------------------- -------------------------------------- -- ---------------------------------------------------------------------------------------- -- --

In the $\zeta_{0}$-${\rm log}_{10}T^{\rm Feed}_{\rm vir}$ plane for HERA in lower-left panel of Fig. \[fig:2D\_HERA\_SKA\_singlez\], we observe at the 95 per cent confidence level a multimodal distribution for $\zeta_{0}$ and $T^{\rm Feed}_{\rm vir}$ as highlighted by the lower streaks for high $\zeta_{0}$ and low $T^{\rm Feed}_{\rm vir}$. This region of parameter space is capable of reionizing the [IGM]{} earlier than our fiducial [EoR]{} model, due to a brighter population of lower mass galaxies (see Section \[sec:variation\]). This is highlighted by overlaying the $\bar{x}_{{\ifmmode\mathrm{H\,{\scriptscriptstyle I}}\else{}H\,{\scriptsize I}\fi}{}}$ isocontours for the [EoR]{} parameter space, where it is clear that these two streaks reproduce different [IGM]{} neutral fractions. Correspondingly, these models result in a smaller, secondary feature around $\bar{x}_{{\ifmmode\mathrm{H\,{\scriptscriptstyle I}}\else{}H\,{\scriptsize I}\fi}{}} = 0.4$ in the [IGM]{} neutral fraction PDF. Models in this region of the [EoR]{} parameter space have similar large-scale 21 cm power, but different small-scale structure which HERA cannot differentiate. The SKA, due to its higher sensitivity on small scales does not exhibit the same behaviour. Alternatively, these degeneracies can be broken with additional redshift observations which constrain the redshift evolution of the large-scale power (see below).

In order to quantitatively assess the performance of [<span style="font-variant:small-caps;">21CMMC</span>]{} we choose to report the median value of each [EoR]{} parameter and the associated 16th and 84th percentiles. This choice follows on from the fact that the recovered 1D PDFs do not need to be normally distributed within the Bayesian MCMC framework, as observed in the case of HERA in Fig. \[fig:2D\_HERA\_SKA\_singlez\]. In Table \[tab:3param\_summary\_singlez\] we summarize the recovered [EoR]{} parameter constraints for each of the 21 cm experiments from a single epoch observation of the 21 cm PS.

For HERA, we note that while both $R_{\rm mfp}$ and $T^{\rm Feed}_{\rm vir}$ (also $\bar{x}_{{\ifmmode\mathrm{H\,{\scriptscriptstyle I}}\else{}H\,{\scriptsize I}\fi}{}}$), are recovered to within 5 per cent of our fiducial parameters, in the case of $\zeta_{0}$ we over predict the expected value by close to 30 per cent. This is not surprising given the mild positive asymmetry observed in the recovered 1D PDF for $\zeta_{0}$. Since we choose to report the median recovered value for all [EoR]{} parameters, whenever we recover an asymmetric distribution our reported value can be shifted away from the input fiducial value. Despite the skewness, the 16th and 84th percentiles should for the most part enclose our fiducial values.

In the case of the SKA, we are able to recover all [EoR]{} parameters values at less than 10 per cent deviation from their fiducial values. While the errors on $R_{\rm mfp}$ remain relatively consistent with those for HERA, we observe a reduction of close to a factor of 2 (3) for the lower (upper) 1$\sigma$ bounds on $\zeta_{0}$ and a 30 per cent reduction on the 1$\sigma$ errors on $T^{\rm Feed}_{\rm vir}$. We note a reduction also of approximately 25 per cent on the recovered 1$\sigma$ errors for the [IGM]{} neutral fraction.

Single redshift observation with a neutral fraction prior {#sec:NFPrior}
---------------------------------------------------------

In the previous section, we noted that the same large-scale power could occur at different [IGM]{} neutral fractions. This degeneracy can be broken with increased small-scale sensitivity (as in the case of SKA). Here, we investigate whether we can improve upon the available constraining power from a single redshift observation by instead including a prior on the [IGM]{} neutral fraction from observations at the tail end of the reionization epoch. For example, such constraints can be obtained from the spectra of high-$z$ quasars such as ULAS J1120+0641 [@Bolton:2011p1063; @Mortlock:2011p1049] or from the observed drop in Lyman-$\alpha$ emission from high-$z$ galaxies [e.g. @Caruana:2014p1585; @Pentericci:2014p1531; @Dijkstra:2014p3264; @Mesinger:2015p1584]. Here, we consider a purely illustrative conservative lower limit on the [IGM]{} neutral fraction at $z=7$ of 5 per cent. We caution however that priors on different epochs can make our results more model dependent. For example, our fiducial [EoR]{} model does not allow for a redshift evolution in the reionization parameters. This is straightforward to account for, however, we postpone it to future work.

In Fig. \[fig:2D\_HERA\_prior\], we compare the recovered [EoR]{} parameter constraints for HERA from the same 1000 h observation of the 21 cm PS at $z=9$, following the inclusion of the [IGM]{} neutral fraction prior at $z=7$. Additionally, in Table \[tab:3param\_summary\_singlez\], we again provide the recovered median and 1$\sigma$ errors for our [EoR]{} parameters. After the inclusion of the $z=7$ neutral fraction prior we observe only a very marginal improvement in either the median recovered values or their 1$\sigma$ errors and this is equivalently shown by the slightly narrower widths of the recovered 1D PDFs.

However, in the case of the 2D joint likelihood distributions we are able to observe an improvement in the allowed [EoR]{} parameter space. Specifically, focusing again on the $\zeta_{0}$-${\rm log}_{10}T^{\rm Feed}_{\rm vir}$ plane, while the 1$\sigma$ errors are only marginally improved, the 2$\sigma$ contours are notably reduced. Most importantly, the [IGM]{} neutral fraction prior removes the second degenerate streak which marginally allowed a significantly different reionization history in the absence of the prior when observed with HERA. This can clearly be seen from the removal of the previously noted secondary bump at $\bar{x}_{{\ifmmode\mathrm{H\,{\scriptscriptstyle I}}\else{}H\,{\scriptsize I}\fi}{}} = 0.4$ in the [IGM]{} neutral fraction PDF. By imposing that the [IGM]{} must at least be 5 per cent neutral at $z=7$, this rules out reionization scenarios which would otherwise predict an early end to the reionization epoch from a single redshift observation of the 21 cm PS. This highlights the importance of the additional leverage gained from probing the reionization history at improving the [EoR]{} reionization constraints, especially in the absence of strong telescope sensitivity at small scales.

The level of improvement achievable from an [IGM]{} neutral fraction prior on a single epoch measurement of the 21 cm PS depends on numerous factors. First, imposing a harder prior than our conservative choice would allow increased constraining power on the allowed [EoR]{} parameter space. Secondly, with decreasing telescope sensitivity, a neutral fraction prior becomes more important, as knowledge of the EoR history compensates for poor single epoch PS coverage. It is feasible that with improved foreground removal strategies, a marginal detection of the 21 cm PS could be achieved with first generation instruments [@Pober:2014p35]. Therefore, in this scenario, a neutral fraction prior would be of vital importance for improving the available [EoR]{} constraining power.

Finally, due to the flexibility of [<span style="font-variant:small-caps;">21CMMC</span>]{} we emphasize that it is straightforward to additionally include direct priors on our model [EoR]{} parameters motivated by current and upcoming observations and theoretical advances. Given the observed degeneracies between the [EoR]{} parameters within this model, it is clear that including priors on our model parameters would notably improve the overall constraining power.

  -------------------------------------------------- ------------------------------------ ----------------------------------- -------------------------------------- --------------------------------- ---------------------------------------------------------------------------------------- ---------------------------------- --
  Instrument                                                                                           Parameter                                                                                        $\bar{x}_{{\ifmmode\mathrm{H\,{\scriptscriptstyle I}}\else{}H\,{\scriptsize I}\fi}{}}$                                     
  (multiz-$z$)                                                   $\zeta_{0}$                      $R_{\rm mfp}$ (Mpc)          log$_{10}$$(T^{\rm Feed}_{\rm vir})$               $z = 8$                                                      $z = 9$                                                       $z = 10$              
  LOFAR                                               $45.21\substack{+26.27 \\ -18.06}$   $12.98\substack{+ 4.44\\ -5.50}$      $4.73\substack{+0.34 \\ -0.32}$      $0.57\substack{+ 0.20\\ -0.21}$                              $0.76\substack{+0.13 \\ -0.15}$                               $0.88\substack{+ 0.07\\ -0.10}$   
  HERA 331                                             $30.17\substack{+7.60 \\ -5.35}$    $15.13\substack{+ 2.92\\ -3.01}$      $4.49\substack{+ 0.15\\ -0.16}$      $0.49\substack{+ 0.05\\ -0.05}$                              $0.71\substack{+ 0.04\\ -0.04}$                               $0.84\substack{+ 0.03\\ - 0.03}$  
  SKA                                                  $30.04\substack{+ 5.65\\ -4.48}$    $15.22\substack{+ 2.82 \\ -2.91}$     $4.48\substack{+ 0.11\\ - 0.12}$     $0.49\substack{+0.04 \\ -0.04}$                              $0.71\substack{+ 0.03\\ -0.03}$                               $0.84\substack{+ 0.02\\ - 0.02}$  
  HERA 331 (w/o 25 per cent modelling uncertainty)     $28.78\substack{+ 2.54\\ -1.98}$    $14.33\substack{+ 1.12 \\ -0.93}$     $4.47\substack{+ 0.06\\ - 0.06}$     $0.47\substack{+0.02 \\ -0.02}$                              $0.69\substack{+ 0.02\\ -0.02}$                               $0.82\substack{+ 0.01\\ - 0.01}$  
  SKA (w/o 25 per cent modelling uncertainty)          $30.31\substack{+ 1.70\\ -1.62}$    $15.47\substack{+ 1.41 \\ -1.20}$     $4.48\substack{+ 0.04\\ - 0.04}$     $0.48\substack{+0.01 \\ -0.01}$                              $0.71\substack{+ 0.01\\ -0.01}$                               $0.83\substack{+ 0.01\\ - 0.01}$  
  -------------------------------------------------- ------------------------------------ ----------------------------------- -------------------------------------- --------------------------------- ---------------------------------------------------------------------------------------- ---------------------------------- --

Combining multiple redshift observations {#sec:multiz}
----------------------------------------

Here, we consider combining multiple epoch observations of the 21 cm PS, assuming a broader instantaneous bandwidth is available for the [EoR]{} observations. However, this comes at the cost of increased model dependence. For example, our three empirical [EoR]{} parameters are redshift-independent, whereas physically, we would expect each parameter to vary with redshift across the reionization epoch. In this sense, our astrophysical parameters can be thought of as ‘effective’ EoR properties, averaging over any redshift evolution. In future, we will explicitly include a generalized redshift evolution for [EoR]{} parameters (e.g. @Barkana:2009p116 [@Mesinger:2012p1131; @Zahn:2012p1156]).

To combine all the available information across the multiple epochs and differentiate between the sampled [EoR]{} parameters, we compute the total maximum likelihood fit of the [EoR]{} model, obtained by linearly combining the individual $\chi^{2}$ statistics from each redshift. In principle, we could instead perform a weighted summation to estimate the total maximum likelihood; however, in neglecting to do so our errors are more conservative.

In Fig. \[fig:2D\_LOFAR\_HERA\_SKA\_multiz\], we compare the recovered 1D marginalized PDFs and 2D marginalized likelihood contours for our three fiducial [EoR]{} parameters ($\zeta_{0}$, $R_{\rm mfp}$ and log$_{10}$$(T^{\rm Feed}_{\rm vir})$). For this, we assume three independent 1000 h observations of the 21 cm PS at $z=8$, 9 and 10 for LOFAR (turquoise), HERA (red) and SKA (blue) each with a 1000 h integration time. In Fig. \[fig:LOFAR\_HERA\_SKA\_multiz\_NF\] we provide the recovered 1D PDFs of the [IGM]{} neutral fraction from each individual epoch. Note, for our fiducial [EoR]{} model, the corresponding [IGM]{} neutral fractions are $\bar{x}_{{\ifmmode\mathrm{H\,{\scriptscriptstyle I}}\else{}H\,{\scriptsize I}\fi}{}}$ = 0.48, 0.71 and 0.83, respectively. Finally, in Table \[tab:3param\_summary\] we provide the recovered median, 16th and 84th percentiles for each [EoR]{} parameter and the corresponding [IGM]{} neutral fraction.

From Fig. \[fig:2D\_LOFAR\_HERA\_SKA\_multiz\], we observe the 1D PDFs and joint 2D likelihood contours to be almost indistinguishable between HERA and the SKA, with the SKA performing only marginally better. Furthermore, comparing with Fig. \[fig:2D\_HERA\_SKA\_singlez\], the widths of both the 1D PDFs and 2D likelihood contours have significantly reduced, drastically in the case of HERA. This emphasizes the strength of combining multiple epoch observations of the 21 cm PS (under the simple assumption they can be linearly combined). Quantitatively, we find little difference between the recovered constraints, with both instruments recovering median values of our [EoR]{} parameters to better than 1.5 per cent. Due to the improved constraining power from sampling the reionization history the recovered PDFs are now normally distributed. Therefore, we can estimate the total 1$\sigma$ fractional error for each parameter. For SKA (HERA), we find errors of 16.7 (22.0) on $\zeta_{0}$, 17.8 (18.4) on $R_{\rm mfp}$ and 2.4 (3.3) per cent on log$_{10}$$(T^{\rm Feed}_{\rm vir})$. Additionally, across all three observed epochs, both HERA and SKA tightly recover the expected [IGM]{} neutral fraction of our fiducial [EoR]{} model (see Fig. \[fig:LOFAR\_HERA\_SKA\_multiz\_NF\]).

Importantly, while for a single epoch observation, the SKA notably outperforms HERA, the improvements available when probing the reionization history can make these two instruments comparable. By combining multiple epoch observations, these two experiments are able to pick up the same evolution in the 21 cm PS as the power rises and falls during reionization (see, e.g. @Lidz:2008p1744 and our Section \[sec:variation\]), compensating for their different sensitivity at small scales (for the case of these EoR models). Therefore, while the additional small-scale power can aid significantly in improving the constraining power from an individual observation, it only marginally improves the constraining power across multiple epochs.

For LOFAR, we find that by jointly combining multiple epoch observations we are now capable of recovering reionization constraints. As expected though, the quality of these constraints is drastically reduced compared to the second-generation experiments. Additionally, we recover broad marginalized 1D PDFs, which encompass the majority of the [EoR]{} parameter space and their asymmetric nature complicates a quantitative comparison with the second-generation experiments. We find LOFAR can only marginally rule out (at $\sim 2\sigma$) either a rapid reionization ($\bar{x}_{{\ifmmode\mathrm{H\,{\scriptscriptstyle I}}\else{}H\,{\scriptsize I}\fi}{}} < 0.2$ at $z=8$) or a delayed/extended reionization ($\bar{x}_{{\ifmmode\mathrm{H\,{\scriptscriptstyle I}}\else{}H\,{\scriptsize I}\fi}{}} > 0.9$ at $z=8$). However, this is still encouraging for the first generation of 21 cm experiments given our conservative foreground and sensitivity assumptions.

Impact of the [EoR]{} modelling uncertainty {#sec:ModUncert}
-------------------------------------------

Thus far we have discussed the available constraints on our [EoR]{} parameters including a 25 per cent [EoR]{} modelling uncertainty, which dominates on large scales (see Fig. \[fig:3param\_variation\]). Characterization of these uncertainties (e.g. by calibrating to more detailed RT simulations) will be essential to maximizing the achievable scientific gains from upcoming second-generation interferometers. To emphasize this, in Fig. \[fig:2D\_HERA\_SKA\_multiz\_ModUncert\] we perform the same analysis as in the previous section, instead now removing the 25 per cent error for HERA and the SKA, corresponding to a perfect characterization of the [EoR]{} modelling uncertainties. We also list the corresponding constraints in Table \[tab:3param\_summary\]. As expected, the resulting fractional errors are dramatically improved by a factor of 2–3 to 8.1/5.5, 7.2/8.1 and 1.4/0.9 per cent for HERA/SKA (for $\zeta_{0}$, $R_{\rm mfp}$ and log $T^{\rm Feed}_{\rm vir}$, respectively). The fact that the parameter constraints improve dramatically and in a comparable manner for both HERA and the SKA shows that the redshift evolution of large-scale power dominates the constraining power for this EoR model. Unlike for narrow-band, single-$z$ observations in Fig. \[fig:2D\_HERA\_SKA\_singlez\], the additional leverage provided by small scales does not dramatically improve the parameter constraints if the redshift evolution of the large-scale power can be constrained. We caution however that this conclusion is dependent on the EoR model.

We also note that these values without the modelling error are in reasonable agreement with the roughly comparable[^13] constraints recovered by @Pober:2014p35, suggesting the Fisher matrix approach provides a reasonably accurate description of available constraints for this simple EoR parametrization.

Two galaxy populations with mass-dependent ionizing efficiencies {#sec:DPL}
================================================================

In the previous section, we highlighted the strength of [<span style="font-variant:small-caps;">21CMMC</span>]{} at providing constraints on our model parameters from a simple single ionizing source population. However, physically we expect galaxy formation to be more complicated. Therefore, in this section, we consider a more flexible [EoR]{} model which allows for both a halo mass-dependent ionizing efficiency and multiple ionizing source populations.

that for our fiducial choice of $\beta=-1/3$ we roughly obtain the scaling, $f_{\rm esc}\propto M^{-0.18}$ (see the text for further discussions). \[fig:UVLum\_DPL\]

Generalized five parameter [EoR]{} model
----------------------------------------

Rather than assuming a step-function for the ionizing efficiency as in the previous section, here we allow it to vary with the host halo potential both for the bright, star-forming galaxies and the faint, feedback limited galaxies, $$\begin{aligned}
 \label{eq:zeta_M}
\zeta(T_{\rm vir}) = \zeta_{0}\left\{ \begin{array}{cl} 
      \left(\frac{T_{\rm vir}}{T^{\rm Feed}_{\rm vir}}\right)^{\beta}  & {T_{\rm vir} \geq T^{\rm Feed}_{\rm vir}}\\
      \left(\frac{T_{\rm vir}}{T^{\rm Feed}_{\rm vir}}\right)^{\alpha} & {10^{4}\,{\rm K} \leq T_{\rm vir} < T^{\rm Feed}_{\rm vir}}\\ 
      0 & {T_{\rm vir} < 10^{4}\,{\rm K}}.
\end{array} \right.\end{aligned}$$ Here, $\zeta_{0}$, is simply the normalization of the ionization efficiency and retains the same definition as in Section \[sec:Zeta\]. When $\beta = 0$ and $\alpha\rightarrow \infty$, this expression is equivalent to the step function ionizing efficiency we considered earlier. A mass-dependent ionizing efficiency has previously been investigated to describe the ionizing galaxy population; however, these only considered a single population of galaxies described by a single power-law index [e.g. @Furlanetto:2006p1647; @McQuinn:2007p1665; @Paranjape:2014p121]. A sharp mass-scale separating feedback limited and non-limited galaxies is expected in some theoretical models [e.g. @Finlator:2011p3408; @Raicevic:2011p3440].

Our generalized [EoR]{} model contains a total of five parameters, the same three [EoR]{} parameters as described in Section \[sec:21CMFAST\] ($\zeta_{0}$, $R_{\rm mfp}$ and $T^{\rm Feed}_{\rm vir}$) and the two new power-law indices describing the mass-dependent ionizing efficiency, $\alpha$ and $\beta$. For our fiducial mock observation, we choose to adopt $\alpha = 4/3$ and $\beta = -1/3$, both within the allowed range, $\alpha$, $\beta\in[-4/3,10/3]$. If instead we define equation (\[eq:zeta\_M\]) with respect to mass, these power-law indices become $\alpha = 2$ and $\beta = -1/2$ (since $M\propto T^{3/2}_{\rm vir}$). Our specific choice of $\alpha$ and $\beta$ for our fiducial model is purely illustrative, and allows a strong mass-dependent scaling for the feedback limited galaxies, with a flatter slope contribution from the high-mass galaxies. In the local Universe, observations show that for low-mass galaxies the star formation efficiency may increase with host halo mass, $f_{\star} \propto M^{2/3}$ [@Kauffmann:2003p1701]. On the other hand, our framework also allows small-mass galaxies to have even higher ionizing efficiencies than the high-mass galaxies. Such models [e.g. @Alvarez:2012p3266] could be motivated by a much higher escape fraction in small galaxies [e.g. @Wise:2014p3267]. For our remaining [EoR]{} parameters, we adopt $\zeta_{0}$ = 70, ${\rm log}_{10}T^{\rm Feed}_{\rm vir} = 4.7$ and $R_{\rm mfp}=15$ Mpc.

Again, we provide an illustrative example of our generalized [EoR]{} model by constructing a schematic picture of the expected $z=8$ UV luminosity function in Fig. \[fig:UVLum\_DPL\] following the same approach as in Section \[sec:21CMFAST\]. We convert our new critical feedback temperature, $T^{\rm Feed}_{\rm vir}$ into a UV magnitude (total halo mass) of $M_{\rm UV} = -13.7$ ($M\sim10^{9.11}\,M_{\sun}$), and denote this and the atomic cooling threshold by the vertical dashed lines. For the massive (bright) galaxies, denoted by $\beta = -1/3$, we use the same fits to the UV luminosity function in @Kuhlen:2012p1506. As before, this fit for these bright galaxies can be translated to $f_{\rm esc} \propto M^{-0.18}$ (see the discussion in Section \[sec:UVLum\]).

  ------------- ------------------------------------ ---------------------------------- --------------------------------- ---------------------------------- -------------------------------------- --------------------------------- ---------------------------------------------------------------------------------------- ---------------------------------- --
  Instrument                                                                                        Parameter                                                                                                                          $\bar{x}_{{\ifmmode\mathrm{H\,{\scriptscriptstyle I}}\else{}H\,{\scriptsize I}\fi}{}}$                                     
  (multi-$z$)               $\zeta_{0}$                     $R_{\rm mfp}$ (Mpc)                     $\alpha$                           $\beta$                log$_{10}$$(T^{\rm Feed}_{\rm vir})$               $z = 8$                                                      $z = 9$                                                       $z = 10$              
  HERA 331       68.28$\substack{+20.27 \\ -25.96}$   14.99$\substack{+2.89 \\ -2.36}$   1.10$\substack{+1.28 \\ -0.72}$   -0.31$\substack{+0.39 \\ -0.45}$     4.74$\substack{+0.42 \\ -0.22}$      0.28$\substack{+0.04 \\ -0.03}$                              0.55$\substack{+0.04 \\ -0.04}$                               0.74$\substack{+0.04 \\ -0.04}$   
  SKA            70.20$\substack{+18.62 \\ -25.68}$   14.87$\substack{+2.79 \\ -2.13}$   1.11$\substack{+1.33 \\ -0.70}$   -0.33$\substack{+0.37 \\ -0.43}$     4.73$\substack{+0.44 \\ -0.21}$      0.28$\substack{+0.03 \\ -0.03}$                              0.55$\substack{+0.03 \\ -0.03}$                               0.74$\substack{+0.03 \\ -0.03 }$  
  ------------- ------------------------------------ ---------------------------------- --------------------------------- ---------------------------------- -------------------------------------- --------------------------------- ---------------------------------------------------------------------------------------- ---------------------------------- --

Constraints from combining multiple redshift observations
---------------------------------------------------------

Given the increased parameter degeneracies in this generalized five parameter [EoR]{} model, we will focus only on multiple epoch observations of the 21 cm PS. As in Section \[sec:multiz\], we assume three independent 1000 h observations of the 21 cm PS at $z=8$, 9 and 10; however, only for HERA and the SKA. In Fig. \[fig:1DPDFs\_HERA\_SKA\_DPL\] we show the recovered 1D marginalized PDFs for each of the five [EoR]{} model parameters from HERA (red curves) and SKA (blue curves) and in Fig. \[fig:2D\_HERA\_SKA\_DPL\] we show the corresponding 2D joint marginalized likelihood contours and the resulting 1D PDFs of the [IGM]{} neutral fraction at each observed epoch. Finally, in Table \[tab:5param\_summary\], we report the recovered median values and the corresponding 16th and 84th percentiles of our [EoR]{} model parameters.

For the majority of these [EoR]{} parameters, we recover rather broad PDFs, noting mild to strong asymmetries emphasizing the strength of the degeneracies. This is to be expected with more model parameters. As an example, focusing on the $\alpha$-${\rm log}_{10}T^{\rm Feed}_{\rm vir}$ panel for a decreasing $T^{\rm Feed}_{\rm vir}$ (increasingly narrower mass range for the feedback limited galaxies), any strong suppression (high $\alpha$) can mimic the same behaviour as our fiducial model. This is due to the fact that for our assumed value of $\alpha = 4/3$ the strong mass suppression in the ionizing efficiency only allows low-mass galaxies closest to $T^{\rm Feed}_{\rm vir}$ to provide any meaningful contribution to reionization.

Despite these degeneracies, we are still able to recover meaningful constraints on these EoR parameters. For example, the contribution of low-mass galaxies to reionization as well as a well-defined feedback scale are picked up in the marginalized 1D $T^{\rm Feed}_{\rm vir}$ PDFs. Such fundamental issues in galaxy formation are unlikely to be addressed with direct observations, even with JWST.

Quantitatively, both HERA and the SKA are able to recover median values for the five [EoR]{} parameters to within 5 per cent of their fiducial values with $\alpha$ being the exception ($\sim20$ per cent). The corresponding percentile errors for $\alpha$ are additionally rather broad, however, this was to be expected from the qualitative analysis above (i.e. a strong mass dependence of $\alpha$ is hard to constrain due to the narrow range in mass of these galaxies which can contribute to reionization). Finally, from Fig. \[fig:2D\_HERA\_SKA\_DPL\], we see that the PDFs of the [IGM]{} neutral fraction are tightly centred around the expected values from our fiducial model, and their quantitative errors are equivalent to our simple three parameter [EoR]{} model. This indicates that we have extracted all the available information from sampling the reionization history for both our 21 cm experiments. In principle, we could improve these constraints by (i) adding measurements of the 21 cm PS in additional bands (e.g. at $z=7$); (ii) applying priors on to our [EoR]{} model parameters; (iii) considering alternative statistics of the 21 cm signal in combination with the 21 cm PS; or (iv) reducing the intrinsic modelling uncertainty.

It is important to note, that the discussion of the degeneracies and constraints of these [EoR]{} parameters is specific to our chosen fiducial model for the mock observation. Due to the sharp break in the power-law indices describing the two populations of ionizing galaxies and our choice for $T^{\rm Feed}_{\rm vir}$, the contribution to reionization from the high-mass (efficient star-forming) galaxy population can be relatively well constrained (as evident by the recovered PDF for $\beta$). Such a model is motivated by theoretical predictions of a sharp feedback scale [e.g. @Finlator:2011p3408; @Raicevic:2011p3440]. If however, we select a shallower break in the power-law indices, we effectively spread out the allowed range of masses of the ionizing galaxies (i.e. a shallower, positive $\alpha$ allows an increased contribution from the lower mass galaxies as they are not as efficiently suppressed). To illustrate this pessimistic case, in Fig. \[fig:1DPDFs\_HERA\_SKA\_DPL\_alternative\] we consider ($\zeta_{0}$, $R_{\rm mfp}$, $\alpha$, $\beta$, $T^{\rm Feed}_{\rm vir}$) = (30, 15 Mpc, 2/5, 1/15, 4.8). As expected, we now recover much broader PDFs. Most notably, with the shallower break in the power-law indices, it becomes increasingly difficult to constrain $T^{\rm Feed}_{\rm vir}$ due to the smoother transition between the ionizing efficiency of the two galaxy populations. This is understandable, as in this model the two galaxy populations resemble a single population.

Conclusion {#sec:conclusion}
==========

The [EoR]{} is rich in astrophysical information, probing the nature of the first stars and galaxies, as well as their impact on the [IGM]{}. However, our current understanding of EoR astrophysics is poorly understood. Our understanding will rapidly improve in the near future with the detection of the cosmic 21 cm signal from a broad range of dedicated 21 cm experiments. However, we are faced with the fundamental question, *what exactly can we learn from these observations?*

To this end, we developed [<span style="font-variant:small-caps;">21CMMC</span>]{}[^14], a new, massively parallel MCMC analysis tool specifically designed to recover constraints on any astrophysical parameters during the [EoR]{} from observations of the cosmic 21 cm signal. [<span style="font-variant:small-caps;">21CMMC</span>]{} utilizes existing astrophysical simulation codes and parameter estimation techniques to ensure its efficiency and flexibility. It combines the astrophysical parameter sampler [<span style="font-variant:small-caps;">cosmohammer</span>]{} [@Akeret:2012p842], calling the seminumerical 21 cm simulation code [<span style="font-variant:small-caps;">21CMFAST</span>]{} [@Mesinger:2007p122; @Mesinger:2011p1123] at each point in the sampled [EoR]{} astrophysical parameter space. This allows sampling over $10^{5}$ individual realizations of the cosmic 21 cm signal to recover robust [EoR]{} parameter constraints, without common assumptions of Gaussian errors.

Throughout this work, we have emphasized several key strengths of [<span style="font-variant:small-caps;">21CMMC</span>]{}, outlining its applicability to a broad range of studies of the 21 cm signal during the [EoR]{}. In summary, we

-   showed by utilizing a Bayesian MCMC framework, we can efficiently and accurately recover robust astrophysical constraints by recovering the marginalized 1D PDFs, irrespective of any model degeneracies. This removes the fundamental limitations and assumptions inherent in previous [EoR]{} astrophysical parameter studies such as sampling fixed grids [e.g. @Choudhury:2005p2859; @Barkana:2009p116; @Mesinger:2012p1131; @Zahn:2012p1156], Fisher Matrix approaches [@Pober:2014p35] or simplistic analytic modelling of the [EoR]{} [@Harker:2012p2856; @Morandi:2012p2857; @Patil:2014p2858];

-   highlighted that it is generalizable to any theoretical model of the [EoR]{}. We showcased two such models: (i) a popular, three parameter [EoR]{} model driven by a single population of star-forming galaxies characterized by a constant (mass-independent) ionizing efficiency; and (ii) a generalized five parameter [EoR]{} model containing two ionizing galaxy populations each characterized by a unique mass-dependent power-law ionizing efficiency;

-   highlighted the available constraining power from both single redshift and combining multiple redshift observations of the ionization 21 cm PS for three current and proposed 21 cm experiments, LOFAR, HERA and the SKA;

-   recovered for LOFAR/HERA/SKA fractional errors of 45.3/22.0/16.7, 33.5/18.4/17.8 and 6.3/3.3/2.4 per cent on the ionizing efficiency (escape fraction), mean free path of ionizing photons and the minimum halo mass hosting star-formation from combining three independent 1000hr observations of the 21 cm PS at $z=8$, 9 and 10 assuming no priors and a conservative noise estimate;

-   showed that by removing the EoR modelling uncertainty from our fiducial analysis, these constraints can be improved by up to a factor of a $\sim$few. This result motivates further work on characterizing EoR modelling errors in order to maximize science returns from second-generation instruments;

-   quantified the improvements in constraining power on the [EoR]{} parameters from accessing the additional information from the cosmic reionization history through the application of physically motivated priors on the [IGM]{} neutral fraction at the tail-end of reionization;

-   emphasized that [<span style="font-variant:small-caps;">21CMMC</span>]{} additionally allows for the inclusion of individual priors on the [EoR]{} model parameters motivated by current and upcoming observations and theoretical advances. Furthermore, while we have focused solely on the 21 cm PS, it is easily adaptable to consider any statistical measure of the cosmic 21 cm signal.

Finally, due its applicability to a broad range of [EoR]{} studies, [<span style="font-variant:small-caps;">21CMMC</span>]{} could be an important analysis tool for the 21 cm [EoR]{} community. For example, it can be used to quantify how foregrounds and other contaminants can hinder the recovery of EoR astrophysics. Moreover, it can serve to guide designs and observing strategies for current and future 21 cm experiments to maximize their scientific returns.

Acknowledgements {#acknowledgements .unnumbered}
================

We thank the anonymous referee for their helpful suggestions. We are grateful to Jonathan Pober, Adrian Liu and Adam Lidz for valuable comments on a draft version of this manuscript. Additionally, we thank Jonathan Pober for making his telescope sensitivity code [[21]{}<span style="font-variant:small-caps;">cmsense</span>]{} publicly available and also Joel Akeret and Sebastian Seehars for the development and public release of [<span style="font-variant:small-caps;">cosmohammer</span>]{}.

[^1]: E-mail: bradley.greig@sns.it

[^2]: http://www.lofar.org/

[^3]: http://www.mwatelescope.org/

[^4]: http://eor.berkeley.edu/

[^5]: https://www.skatelescope.org

[^6]: http://reionization.org

[^7]: [<span style="font-variant:small-caps;">21CMMC</span>]{} will be made available at http://homepage.sns.it/mesinger/21CMMC.html. Interested users can obtain a pre-release beta version by contacting the authors.

[^8]: The ionization structure of these subgrid models of recombinations inside unresolved systems does not in fact directly translate to a constant, uniform value of $R_{\rm mfp}$. These authors find that it takes a long time for cosmic [H[II]{}]{} regions to approach the Str[ö]{}mgren limit, with their expansion only slowing more and more as reionization progresses. In other words, the mean free path (which is only defined in terms of the [*instantaneous*]{} recombination rate) is almost always [*larger*]{} than the actual size of the [H[II]{}]{} regions (which depend on the [*cumulative*]{} number of recombinations), even at the late stages of reionization. The resulting ionization fields have a dearth of ionization structure on scales smaller than the instantaneous mean free path. Moreover, the spatial correlations between the sources and sinks of ionizing photons results in sizeable spatial and temporal variations in the mean free path. Nevertheless, an effective, uniform maximum horizon for ionizing photons of $\sim 10$ Mpc can reproduce the 21 cm power spectra of the @Sobacchi:2014p1157 simulations to within tens of percent on relevant scales, at the epochs studied in that work $0.25 \lesssim \bar{x}_{{\ifmmode\mathrm{H\,{\scriptscriptstyle I}}\else{}H\,{\scriptsize I}\fi}{}} \lesssim 0.75$ (Sobacchi, private communication; Mesinger & Greig, in preparation). For historical context, we use the variable ‘$R_{\rm mfp}$’ to denote this effective horizon set by subgrid recombinations, but caution that the actual mean free path is larger than this parameter, as discussed above.

[^9]: https://github.com/jpober/21cmSense

[^10]: These authors found a factor of 2.5 difference between their SKA and HERA designs for their estimated total sensitivity on a measurement of the 21 cm PS at $z=9.5$. However, this factor was a result of a numerical error in their calculation of their SKA sensitivities. Accounting for this, their HERA design has in fact only a marginally reduced sensitivity relative to the SKA for their pessimistic and moderate foreground scenarios (Pober, private communication)

[^11]: http://www.skatelescope.org/wp-content/uploads/2012/07/SKA-TEL-SKO-DD-001-1BaselineDesign1.pdf

[^12]: At present, [<span style="font-variant:small-caps;">21CMMC</span>]{} computes only co-evolving data cubes of the [EoR]{}. In future, we intend to improve the simulation pipeline to mimic a more realistic observational pipeline, taking into account several effects such as finite light-cones, foreground cleaning etc.

[^13]: A like-to-like comparison is difficult for several reasons. First, these authors jointly combine four different epoch observations of the 21 cm PS compared to our three. Since instrument sensitivity increases with decreasing redshift, improved constraining power is available from their additional $z=7$ observation (although it is not a significant improvement as three redshift measurements will have already strongly constrained the reionization history). Secondly, we consider only a 331 telescope antenna for HERA, unlike their 547 telescope design, while for the SKA there is the previously noted numerical error in their total sensitivities.

[^14]: http://homepage.sns.it/mesinger/21CMMC.html
