Simultaneously constraining the astrophysics of reionisation and the epoch of heating with 21CMMC

TL;DR

This paper extends the 21CMMC framework to simultaneously model the epoch of reionisation and epoch of heating, enabling precise constraints on early galaxy properties from 21 cm observations with upcoming radio telescopes.

Contribution

The authors develop an extended Bayesian inference framework that models both reionisation and heating, improving parameter estimation accuracy from 21 cm data.

Findings

HERA and SKA can constrain ionising and X-ray source properties to 1-10% precision.

The ionisation history can be constrained within a few percent.

Ignoring the spin temperature saturation can bias parameter recovery by 3-10 sigma.

Abstract

The cosmic 21 cm signal is set to revolutionise our understanding of the early Universe, allowing us to probe the 3D temperature and ionisation structure of the intergalactic medium (IGM). It will open a window onto the unseen first galaxies, showing us how their UV and X-ray photons drove the cosmic milestones of the epoch of reionisation (EoR) and epoch of heating (EoH). To facilitate parameter inference from the 21 cm signal, we previously developed 21CMMC: a Monte Carlo Markov Chain sampler of 3D EoR simulations. Here we extend 21CMMC to include simultaneous modelling of the EoH, resulting in a complete Bayesian inference framework for the astrophysics dominating the observable epochs of the cosmic 21 cm signal. We demonstrate that second generation interferometers, the Hydrogen Epoch of Reionisation Array (HERA) and Square Kilometre Array (SKA) will be able to constrain ionising and X-ray source properties of the first galaxies with a fractional precision of order $\sim1$-10 per cent (1$\sigma$). The ionisation history of the Universe can be constrained to within a few percent. Using our extended framework, we quantify the bias in EoR parameter recovery incurred by the common simplification of a saturated spin temperature in the IGM. Depending on the extent of overlap between the EoR and EoH, the recovered astrophysical parameters can be biased by $\sim3-10\sigma$.

Tables

Table 1. Table 1: Summary of telescope parameters we use to compute sensitivity profiles (see text for further details).

Parameter	HERA	SKA
Telescope antennae	331	512
Diameter (m)	14	35
Collecting Area (m2)	50 953	492 602
$T_{\mathrm{rec}}$ (K)	100	0.1Tsky + 40
Bandwidth (MHz)	8	8
Integration time (hrs)	1000 (drift)	1000 (tracked)

Table 2. Table 2: A summary of the astrophysical parameters used for our two mock observations, the corresponding electron scattering optical depth, τ e subscript 𝜏 e \tau_{\rm e} , as well as the adopted prior range for 21CMMC . See text for additional details.

	$\zeta$	$R_{\mathrm{mfp}}$	$T_{\mathrm{vir}}^{\min }$	log10$\left(T_{\mathrm{vir}}^{\min }\right)$	log10()	$E_0$	${\log }_{10}\left(N_{HI}\right)$	${\alpha }_X$	${\tau }_e$
Source Model		[Mpc]	[K]	[K]	[erg s-1 $M_{\odot }^{-1}$ yr]	[keV]	[cm-2]
faint galaxies	30.0	15.0	$5\times {10}^4$	4.70	40.0	0.5	21.45	1.0	0.059
bright galaxies	200.0	15.0	$3\times {10}^5$	5.48	40.0	0.5	21.45	1.0	0.056
Allowed range	[10, 250]	[5, 25]	[${10}^4$, ${10}^6$]	[4.0, 6.0]	[38.0, 42.0]	[0.1, 1.5]	[19.3,23.0]	[-0.5, 2.5]

Table 3. Table 3: Summary of the median recovered values (and associated 16th and 84th percentile errors) for the six parameter astrophysical model describing the EoR and EoX, ζ 𝜁 \zeta , R mfp subscript 𝑅 mfp R_{\rm mfp} , log 10 ( T vir min ) subscript superscript 𝑇 min vir (T^{\rm min}_{\rm vir}) , log 10 ( L X <  2 ​ keV / SFR subscript 𝐿 X 2 keV SFR L_{\rm X\,<\,2\,keV}/{\rm SFR} ), E 0 subscript 𝐸 0 E_{0} and α X subscript 𝛼 𝑋 \alpha_{X} . We assume a total 1000hr integration time with both the SKA and HERA. Our fiducial mock observation, corresponding to the faint galaxies model, assumes ( ζ 𝜁 \zeta , R mfp subscript 𝑅 mfp R_{\rm mfp} , log 10 ( T vir min ) subscript superscript 𝑇 min vir (T^{\rm min}_{\rm vir}) , log 10 ( L X <  2 ​ keV / SFR subscript 𝐿 X 2 keV SFR L_{\rm X\,<\,2\,keV}/{\rm SFR} ), E 0 subscript 𝐸 0 E_{0} , α X subscript 𝛼 𝑋 \alpha_{X} ) = (30, 15, 4.7, 40.0, 0.5, 1.0). We also provide the recovered biased constraints when we ignore the IGM spin temperature fluctuations (i.e. T S ≫ T CMB much-greater-than subscript 𝑇 S subscript 𝑇 CMB T_{\rm S}\gg T_{\mathrm{CMB}} ; see Section  5.1 ).

faint galaxies				Parameter
Model/instrument	$\zeta$	$R_{\mathrm{mfp}}$	log10$\left(T_{\mathrm{vir}}^{\min }\right)$	log10()	$E_0$	${\alpha }_X$
		[Mpc]	[K]	[erg s-1 $M_{\odot }^{-1}$ yr]	[keV]
Full X-ray heating
HERA 331	34.69$\begin{array}{c}+9.48\\ -6.50\end{array}$	14.82$\begin{array}{c}+2.47\\ -2.29\end{array}$	4.78$\begin{array}{c}+0.11\\ -0.11\end{array}$	40.08$\begin{array}{c}+0.14\\ -0.12\end{array}$	0.52$\begin{array}{c}+0.08\\ -0.05\end{array}$	1.01$\begin{array}{c}+0.82\\ -0.85\end{array}$
SKA	33.99$\begin{array}{c}+7.18\\ -5.13\end{array}$	14.56$\begin{array}{c}+2.58\\ -2.29\end{array}$	4.76$\begin{array}{c}+0.07\\ -0.07\end{array}$	40.08$\begin{array}{c}+0.10\\ -0.10\end{array}$	0.52$\begin{array}{c}+0.07\\ -0.08\end{array}$	0.87$\begin{array}{c}+0.86\\ -0.85\end{array}$
Ignoring $T_S$ ($T_S\gg T_{\mathrm{CMB}}$)
HERA 331	18.98$\begin{array}{c}+7.81\\ -2.80\end{array}$	16.80$\begin{array}{c}+3.36\\ -2.67\end{array}$	4.44$\begin{array}{c}+0.24\\ -0.13\end{array}$	-	-	-
SKA	15.51$\begin{array}{c}+4.97\\ -2.20\end{array}$	16.99$\begin{array}{c}+3.75\\ -4.17\end{array}$	4.29$\begin{array}{c}+0.17\\ -0.11\end{array}$	-	-	-

Table 4. Table 4: The same as Table  3 except now for the bright galaxies model. Our bright galaxies mock observation assumes ( ζ 𝜁 \zeta , R mfp subscript 𝑅 mfp R_{\rm mfp} , log 10 ( T vir min ) subscript superscript 𝑇 min vir (T^{\rm min}_{\rm vir}) , log 10 ( L X <  2 ​ keV / SFR subscript 𝐿 X 2 keV SFR L_{\rm X\,<\,2\,keV}/{\rm SFR} ), E 0 subscript 𝐸 0 E_{0} , α X subscript 𝛼 𝑋 \alpha_{X} ) = (200, 15, 5.48, 40.0, 0.5, 1.0).

bright galaxies	Parameter
Model/instrument	$\zeta$	$R_{\mathrm{mfp}}$	log10$\left(T_{\mathrm{vir}}^{\min }\right)$	log10()	$E_0$	${\alpha }_X$
		[Mpc]	[K]	[erg s-1 $M_{\odot }^{-1}$ yr]	[keV]
Full X-ray heating
HERA 331	188.59$\begin{array}{c}+32.66\\ -29.73\end{array}$	15.21$\begin{array}{c}+1.75\\ -1.65\end{array}$	5.47$\begin{array}{c}+0.03\\ -0.03\end{array}$	40.00$\begin{array}{c}+0.08\\ -0.08\end{array}$	0.50$\begin{array}{c}+0.08\\ -0.07\end{array}$	0.95$\begin{array}{c}+0.84\\ -0.86\end{array}$
SKA	177.32$\begin{array}{c}+33.62\\ -30.85\end{array}$	16.33$\begin{array}{c}+2.91\\ -2.08\end{array}$	5.47$\begin{array}{c}+0.02\\ -0.02\end{array}$	40.00$\begin{array}{c}+0.07\\ -0.07\end{array}$	0.50$\begin{array}{c}+0.07\\ -0.08\end{array}$	0.95$\begin{array}{c}+0.91\\ -0.82\end{array}$
Ignoring $T_S$ ($T_S\gg T_{\mathrm{CMB}}$)
HERA 331	25.50$\begin{array}{c}+12.42\\ -4.68\end{array}$	13.15$\begin{array}{c}+1.25\\ -1.19\end{array}$	4.55$\begin{array}{c}+0.22\\ -0.13\end{array}$	-	-	-
SKA	15.65$\begin{array}{c}+3.47\\ -2.22\end{array}$	11.17$\begin{array}{c}+2.06\\ -1.67\end{array}$	4.18$\begin{array}{c}+0.13\\ -0.11\end{array}$	-	-	-

Equations

$$\delta T_{\mathrm{b}}(\nu )$$

$$\zeta f_{\mathrm{coll}}(\mathbf{x},z,R,{\bar{M}}_{\min })\ge 1,$$

$$J(\text{x},E,z)=\frac{(1+z{)}^3}{4\pi }{\int }_z^{\infty }d{z}^{\prime }\frac{cdt}{d{z}^{\prime }}{ϵ}_{\mathrm{X}}{e}^{-\tau }.$$

$${ϵ}_{\mathrm{X}}(\text{x},E_{\mathrm{e}},z^{\prime })=\frac{L_{\mathrm{X}}}{\mathrm{SFR}}\left[{\rho }_{\mathrm{crit},0}{\Omega }_b{f}_{\ast }(1+{\delta }_{\mathrm{nl}})\frac{d{f}_{\mathrm{coll}}(z^{\prime })}{dt}\right],$$

$$L_{\mathrm{X}<2\mathrm{keV}}/\mathrm{SFR}={\int }_{E_0}^{2\mathrm{keV}}d{E}_e{L}_{\mathrm{X}}/\mathrm{SFR}.$$

$$\zeta =30\left(\frac{f_{\mathrm{esc}}}{0.12}\right)\left(\frac{f_{\ast }}{0.05}\right)\left(\frac{N_{\gamma /b}}{4000}\right)\left(\frac{1.5}{1+n_{\mathrm{rec}}}\right)$$

$$M_{\mathrm{vir}}^{\min }$$

$${\Delta }_{\mathrm{N}}^2(k)\approx X^{2}Y\frac{k^3}{2{\pi }^2}\frac{{\Omega }^{\prime }}{2t}{T}_{\mathrm{sys}}^2,$$

$$\delta {\Delta }_{\mathrm{T}+\mathrm{S}}^2(k)={\left(\sum _i\frac{1}{({\Delta }_{\mathrm{N},i}^2(k)+{\Delta }_{21}^2(k){)}^2}\right)}^{-\frac{1}{2}},$$

Full text

1 Introduction

The 21 cm spin-flip transition of neutral hydrogen encodes a treasure trove of cosmological and astrophysical information (see e.g. Gnedin & Ostriker, 1997; Madau et al., 1997; Shaver et al., 1999; Tozzi et al., 2000; Gnedin & Shaver, 2004; Furlanetto et al., 2006; Morales & Wyithe, 2010; Pritchard & Loeb, 2012). The signal is expressed as the offset of the 21 cm brightness temperature, $\delta T_{\rm b}(\nu)$, relative to the cosmic microwave background (CMB) temperature, $T_{\rm CMB}$ (e.g. Furlanetto et al., 2006):

[TABLE]

where $x_{\mathrm{H\,{\scriptscriptstyle I}}{}}$ is the neutral fraction, $T_{\rm S}$ is the gas spin temperature, $\delta_{\rm nl}\equiv\rho/\bar{\rho}-1$ is the gas 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$, where $\nu_{0}$ is the 21 cm frequency.

Because equation (1) depends on the ionisation and temperature as a function of space and time, the 21 cm signal can provide insight into the formation, growth and evolution of structure in the Universe, the nature of the first stars and galaxies and their impact on the physics of the intergalactic medium (IGM; e.g. Barkana & Loeb 2007; Loeb & Furlanetto 2013; Zaroubi 2013). The most widely-studied of these properties is the ionisation state: in the first billion years, the Universe transitioned from being nearly fully neutral to being nearly fully ionised. This epoch of reionisation (EoR) was driven by the percolation of H II regions generated by the ionising photons escaping from the first galaxies. Sourcing the $x_{\mathrm{H\,{\scriptscriptstyle I}}{}}$ term in equation (1), the EoR should be evidenced by a rise and fall of large-scale fluctuations in the 21 cm brightness temperature (e.g. Lidz et al. 2008).

The second astrophysical term in equation (1) is the IGM spin temperature, $T_{\rm S}$. The spin temperature is thought to be regulated by the first sources of X-rays, which can heat the IGM from its post thermal decoupling adiabat to temperatures far above the CMB temperature. While the IGM is still adiabatically cooling, the $(1-T_{\rm CMB}/T_{\rm S})$ term in equation (1) can be of order $-200$. This large dynamic range means that the spatial fluctuations in temperature during this epoch of heating (EoH) can provide the strongest 21 cm signal, more than an order of magnitude larger than that during the EoR (e.g. Mesinger & Furlanetto, 2007; Pritchard & Furlanetto, 2007; Baek et al., 2010; Santos et al., 2010; McQuinn & O’Leary, 2012; Mesinger et al., 2013). The 21 cm signal can therefore be a powerful probe of high-energy processes in the first galaxies. The most likely of these X-ray sources are high mass X-ray binaries (HMXBs; Power et al. 2009; Mirabel et al. 2011; Fragos et al. 2013a; Power et al. 2013) and/or the hot interstellar medium (ISM) within the first galaxies (e.g. Oh, 2001; Pacucci et al., 2014). However, other alternative scenarios have been put forth including metal-free (Pop-III) stars (Xu et al., 2014), mini-QSOs (e.g Madau et al., 2004; Yue et al., 2013; Ghara et al., 2016), dark matter annihilation (e.g. Cirelli et al., 2009; Evoli et al., 2014; Lopez-Honorez et al., 2016) or cosmic rays (Leite et al., 2017).

Numerous 21 cm experiments are currently underway, attempting to detect the cosmic 21 cm signal. These fall into two broad categories. The first are large-scale interferometers, seeking to detect spatial 21 cm fluctuations. These include the Low Frequency Array (LOFAR; van Haarlem et al. 2013; Yatawatta et al. 2013), the Murchison Wide Field Array (MWA; Tingay et al. 2013), the Precision Array for Probing the Epoch of Reionisation (PAPER; Parsons et al. 2010), the Square Kilometre Array (SKA; Mellema et al. 2013) and the Hydrogen Epoch of Reionisation Array (HERA; DeBoer et al. 2017). The second class are single dipole or small compact array global-sky experiments seeking the volume averaged global 21 cm signal. These include the Experiment to Detect the Global EoR Signature (EDGES; Bowman & Rogers 2010), the Sonda Cosmológica de las Islas para la Detección de Hidrógeno Neutro (SCI-HI; Voytek et al. 2014), the Shaped Antenna Measurement of the Background Radio Spectrum (SARAS; Patra et al. 2015), Broadband Instrument for Global HydrOgen ReioNisation Signal (BIGHORNS; Sokolowski et al. 2015), the Large Aperture Experiment to detect the Dark Ages (LEDA; Greenhill & Bernardi 2012; Bernardi et al. 2016), and the Dark Ages Radio Explorer (DARE; Burns et al. 2012).

As a first step in preparation for the wealth of data expected from these 21 cm experiments, we developed a publicly available Monte Carlo Markov Chain (MCMC) EoR analysis tool 21CMMC111https://github.com/BradGreig/21CMMC (Greig & Mesinger, 2015). This is the first EoR analysis tool to sample 3D reionisation simulations (using 21CMFAST; Mesinger & Furlanetto 2007; Mesinger et al. 2011) within a fully Bayesian framework for astrophysical parameter space exploration and 21 cm experiment forecasting. In this introductory work, we adopted the common simplifying assumption of a saturated spin temperature: $T_{\rm S}\gg T_{\rm CMB}$. However, the applicability of this saturated limit is dependent on the poorly-known strength and spectral shape of the X-ray background in the early Universe. If the spin temperature is not fully saturated during the EoR, this could result in sizeable biases in the inferred EoR parameters (e.g. Watkinson & Pritchard, 2015).

Beyond its impact on EoR parameter recovery, the IGM spin temperature also encodes a wealth of information on the high-energy processes in the early Universe, as mentioned above. Second generation interferometers, HERA and SKA, have the bandwidth and sensitivity to easily probe temperature fluctuations during the EoH.

In this work, we extend 21CMMC to include a full treatment of the EoH, retaining our ability to perform on-the-fly sampling of 3D reionisation simulations. This updated 21CMMC is capable of jointly exploring the astrophysics of both the EoR and EoH, allowing us to maximise the scientific return of upcoming second-generation telescopes.

We note that recent studies suggest that machine learning can be a viable alternative to on-the-fly MCMC sampling of 21CMFAST. Shimabukuro & Semelin (2017) used an artificial neural network to predict astrophysical parameters, with an accuracy of $\sim$ tens of percent. This approach is fast, though producing Bayesian confidence limits becomes less straightforward. Alternately, Kern et al. (2017) bypassed the on-the-fly sampling of 3D simulations by using an emulator trained on the 21 cm power spectrum. An emulator can be used in an MCMC framework, and is orders of magnitude faster at parameter forecasting compared to a direct sampling of 3D simulations (such as 21CMMC). This comes at the cost of \mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{<}} ten per cent in power spectrum accuracy over most of the parameter space, when the learning is performed on $\sim 10^{4}$ 21CMFAST training samples (higher precision can be obtained by increasing the size of the training set). Future work will test emulator accuracy on high-order summary statistics.

The remainder of this paper is organised as follows. In Section 2 we summarise 21CMMC and the associated 21CMFAST simulations used to generate 3D realisations of the cosmic 21 cm signal, outlining the updated astrophysical parameterisation to model the EoH. In Section 3 we introduce our mock observations to be used in our 21 cm experiment forecasting, and present the forecasts and associated discussions in Section 4. We then explore the impact of assuming the saturated IGM spin temperature limit in Section 5, before summarising the improvements to 21CMMC and finishing with our closing remarks in Section 6. Unless stated otherwise, we quote all quantities in comoving units and adopt the cosmological parameters: ($\Omega_{\Lambda}$, $\Omega_{\rm M}$, $\Omega_{b}$, $n$, $\sigma_{8}$, $H_{0}$) = (0.69, 0.31, 0.048, 0.97, 0.81, 68 km s*-1* Mpc*-1*), consistent with recent results from the Planck mission (Planck Collaboration XIII, 2016).

2 21CMMC

In this section, we provide a short summary of the main aspects of 21CMMC, before delving into the modelling of the cosmic 21 cm signal during the EoR and EoH in Section 2.1 and developing an intuition about the full parameter set in Section 2.2.

21CMMC is a massively parallel MCMC sampler for exploring the astrophysical parameter space of reionisation. It was built using a modified version of the easy to use python module cosmohammer (Akeret et al., 2013), which employs the emcee python module developed by Foreman-Mackey et al. (2013) based on the affine invariant ensemble sampler of Goodman & Weare (2010). At each proposed step in the computation chain, 21CMMC performs a new, independent 3D reionisation of the 21 cm signal, using an optimised version of the publicly available 21CMFAST simulation code (Mesinger & Furlanetto, 2007; Mesinger et al., 2011) for computational efficiency. Using a likelihood statistic (fiducially the PS), it compares the model against a mock observation generated from a larger simulation with a different set of initial conditions. For further details, we refer the reader to (Greig & Mesinger, 2015).

2.1 Modelling the cosmic 21 cm signal with 21CMFAST

In this work, we use an optimised version of the publicly-available version of 21CMFAST_v1.1222https://github.com/andreimesinger/21cmFAST. 21CMFAST employs approximate but efficient modelling of the underlying astrophysics of the reionisation and heating epochs. The resulting 21 cm PS during the EoR have been found to match those of computationally expensive radiative transfer simulations to with tens of per cent on the scales relevant to 21 cm interferometry, \mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{>}}1 Mpc (Zahn et al., 2011). We refer the reader to Mesinger & Furlanetto (2007); Mesinger et al. (2011) for explicit details on the semi-numerical approach, and only provide a summary below.

21CMFAST produces a full, 3D realisation of the 21 cm brightness temperature field, $\delta T_{\rm b}$ (c.f equation 1) which is dependent on the ionisation, density, velocity and IGM spin temperature fields. The evolved IGM density field at any redshift is obtained from an initial high resolution linear density field which is perturbed using the Zel’dovich approximation (Zel’dovich, 1970) before being smoothed onto a lower resolution grid.

The ionisation field is then estimated from this evolved IGM density field using the excursion-set approach (Furlanetto et al., 2004). The time-integrated number of ionising photons333 This includes ionisations from both UV and X-ray sources. While X-rays can produce some level of pre-reionisation (e.g. Ricotti & Ostriker, 2004; Dijkstra et al., 2012; McQuinn, 2012; Mesinger et al., 2013) their predominant contribution is in pre-heating the IGM before reionisation (e.g. McQuinn & O’Leary, 2012). is compared to the number of neutral atoms within regions of decreasing radius, $R$. These regions are computed from a maximum photon horizon, $R_{\rm mfp}$ down to the individual pixel resolution of a single voxel, $R_{\rm cell}$. A voxel at coordinates $(\boldsymbol{x},z)$ within the simulation volume is then tagged as fully ionised if,

[TABLE]

where $f_{\rm coll}(\boldsymbol{x},z,R,\bar{M}_{\rm min})$ is the fraction of collapsed matter residing within haloes more massive than $\bar{M}_{\rm.min}$ (Press & Schechter, 1974; Bond et al., 1991; Lacey & Cole, 1993; Sheth & Tormen, 1999) and $\zeta$ is an ionising efficiency describing the conversion of mass into ionising photons (see Section 2.2.1). Partial ionisations are included for voxels not fully ionised by setting their ionised fractions to $\zeta f_{\rm coll}(\boldsymbol{x},z,R_{\rm cell},\bar{M}_{\rm min})$.

Since 21 cm observations use the CMB as a background source, the IGM spin temperature has to be coupled to the kinetic gas temperature for the signal to be detected. This coupling is achieved through either collisional coupling or the Ly$\alpha$ background from the first generation of stars (so-called Wouthuysen-Field (WF) coupling (Wouthuysen, 1952; Field, 1958)). To compute the spin temperature, 21CMFAST solves for the evolution of the ionisation, temperature and impinging Ly$\alpha$ background in each voxel444 Note that within this work we do not vary the soft UV spectra of the first sources driving this epoch. This WF coupling epoch will be extremely challenging to detect in comparison to the EoR and the EoH. Nevertheless, we will return to this in future work. For more specifics regarding the computation of the Ly$\alpha$ background, we refer the reader to Section 3.2 of (Mesinger et al., 2011). These depend on the angle-averaged specific intensity, $J(\textbf{x},E,z)$, (in erg s*-1* keV*-1* cm*-2* sr*-1*), computed by integrating the comoving X-ray specific emissivity, $\epsilon_{\rm X}(\textbf{x},E_{e},z^{\prime})$ back along the light-cone:

[TABLE]

Here, $e^{-\tau}$ corresponds to the probability that a photon emitted at an earlier time, $z^{\prime}$, survives until $z$ owing to IGM attenuation (see Eq. 16 of Mesinger et al. 2011) and the comoving specific emissivity is evaluated in the emitted frame, $E_{\rm e}=E(1+z^{\prime})/(1+z)$, with

[TABLE]

where the quantity in square brackets is the star-formation rate (SFR) density along the light-cone, with $\rho_{\rm crit,0}$ being the current critical density and $f_{\ast}$ the fraction of galactic baryons converted into stars.

The quantity $L_{\rm X}/{\rm SFR}$ is the specific X-ray luminosity per unit star formation escaping the galaxies (in units of erg s*-1* keV*-1* $M^{-1}_{\odot}$ yr). The specific luminosity is taken to be a power law in photon energy, $L_{\rm X}\propto E^{-\alpha_{X}}$, with photons below some threshold energy, $E_{0}$, being absorbed inside the host galaxy.555We note that in the computation of the heating and ionisation rates in the default version of 21CMFAST, these integrals were performed out to infinity, which could result in divergent behaviour for $\alpha_{\rm X}\leq 0$. Here we adopt an upper limit of 10 keV for computing the rate integrals. This choice is arbitrary and purely for numerical convenience, as the EoH only depends on the SED below $\sim 2$ keV. Actual SEDs of course do not diverge but turn over at high energies: $\sim$ 10 – 100 keV (e.g. Lehmer et al. 2013, 2015).

Instead of the number of X-ray photons per stellar baryon used in the default version of 21CMFAST, here we normalise the X-ray efficiency in terms of an integrated soft-band ($<2$ keV) luminosity per SFR (in erg s*-1* $M^{-1}_{\odot}$ yr):

[TABLE]

This parametrisation is both physically-motivated (harder photons have mean free paths longer than the Hubble length and so do not contribute to the EoH; e.g. McQuinn 2012; Das et al. 2017), and easier to directly compare with X-ray observations of local star-forming galaxies.

2.2 The 21CMMC astrophysical parameter set

In the previous section, we outlined the semi-numerical approach, using 21CMFAST, to obtain the 3D realisations of the 21 cm brightness temperature field. This parameterisation yields six free parameters to be sampled within our Bayesian framework. In this section, we provide more detailed descriptions for each of these parameters, providing physical intuition for their impact on the IGM through X-ray heating or ionisation and defining their allowed ranges.

To aid in this discussion, we provide Figure 1. For each parameter, denoted (i) - (vi), we provide the redshift evolution of the corresponding average neutral fraction (top panel), average 21 cm brightness temperature contrast (middle panel) and the amplitude of the 21 cm PS at $k=0.15$ Mpc*-1* (bottom panel)666Our choice of $k$-modes was motivated by the fact that the large-scale power in the 21 cm PS is the main discriminant between astrophysical models (e.g. McQuinn et al., 2007; Greig & Mesinger, 2015), Moreover, these $k$-modes are detectable by upcoming experiments (i.e. they are sufficiently small to avoid being contaminated by foregrounds yet large enough to obtain high signal to noise). ($\bar{\delta T^{2}_{\rm b}}\Delta^{2}_{21}(k,z)\equiv k^{3}/(2\pi^{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)/\bar{\delta T_{\rm b}}(z)-1$). In each set of panels, we vary the astrophysical parameter in question across its full allowed span, holding the other five parameters fixed.

The solid, black curve in each panel corresponds to our ‘fiducial’ faint galaxies model (see Section 3 for further details). Before discussing the impact of each parameter, it is instructive to note the general features exhibited by the fiducial model. The EoR history in the top panel has a monotonic evolution with a midpoint around $z\sim 8$. The global 21 cm signal in the middle panel shows a deep absorption trough, corresponding to when the X-ray heating rate surpasses the adiabatic cooling rate (start of the EoH). This is followed by a small emission peak, corresponding to the onset of the EoR (e.g. Furlanetto, 2006). The height of the peak is determined by the relative overlap of the EoR and EoH; if the overlap is strong, reionisation can proceed in a cold IGM, with the global signal never switching to emission.

The large-scale power in the bottom panel shows the characteristic three-peaked structure, driven by fluctuations in the Ly$\alpha$ background (WF coupling); IGM gas temperature (EoH), and the ionisation fraction (EoR), from right to left, respectively (e.g. Pritchard & Furlanetto, 2007; Baek et al., 2010). The redshift position of each peak therefore traces the timing of each epoch, while the amplitude of the 21 cm power traces the level of fluctuations determined by the typical source bias and/or X-ray SED (e.g. Mesinger et al., 2013). Troughs on the other hand correspond to the transitions between epochs, when the corresponding cross terms of the 21 cm PS cause the large-scale power to drop (e.g. Lidz et al., 2007; Pritchard & Furlanetto, 2007; Mesinger et al., 2016).

2.2.1 UV ionising efficiency, $\zeta$

The UV ionising efficiency of high-$z$ galaxies (Equation 2) can be expressed as

[TABLE]

where, $f_{\rm esc}$ is the fraction of ionising photons escaping into the IGM, $f_{\ast}$ is the fraction of galactic gas in stars, $N_{\gamma/b}$ is the number of ionising photons produced per baryon in stars and $n_{\rm rec}$ is the typical number of times a hydrogen atom recombines. While only the product of equation (6) is required for generating the ionisation field, we provide some plausible values for each of the terms. We adopt $N_{\gamma}\approx 4000$ as is typically expected from Population II stars (e.g. Barkana & Loeb, 2005). For high-z galaxies, both $f_{\star}$ and $f_{\rm esc}$ are extremely uncertain, with observations and simulations placing plausible values within the vicinity of $\sim 10$ per cent for $f_{\star}$ (e.g. Behroozi & Silk, 2015; Sun & Furlanetto, 2016) and $f_{\rm esc}$ (e.g Paardekooper et al., 2015; Xu et al., 2016; Kimm et al., 2017). Finally, we adopt $n_{\rm rec}\sim 0.5$, similar to those found in the ’photon-starved’ reionisation models of (Sobacchi & Mesinger, 2014), which are consistent with emissivity estimates from the Ly$\alpha$ forest (e.g. Bolton & Haehnelt, 2007; McQuinn et al., 2011). It is important to note that the $f_{\ast}$ appearing here, is equivalent to the $f_{\ast}$ in equation (4). That is, throughout all cosmic epochs, we adopt the same, constant value of 5 per cent for the fraction of galactic gas in stars. In future work, we will relax this assumption.

In this work, we adopt a flat prior over the fiducial range of $\zeta\in[10,250]$, which is a notable extension over the original upper limit of $\zeta=100$ chosen in Greig & Mesinger (2015). The extended range provides the flexibility to explore physically plausible models where reionisation is driven by bright, rare galaxies (see Section 3). In panel (i) of Figure 1 we highlight the impact of varying $\zeta$. As expected, $\zeta$ has a strong impact on the EoR and an almost negligible impact on the EoH. As $\zeta$ is increased, the EoR peak shifts to earlier redshifts and the width of the EoR peak reduces (i.e. shorter duration). For extremely large values of $\zeta$, the EoR and EoH peaks begin to merge, resulting in a larger amplitude EoR peak, sourced by the contrast between the cold IGM patches present in the early EoH stages and the H II regions.

2.2.2 Maximum ionising photon horizon within ionised regions, $R_{\rm mfp}$

The physical size of H II regions is regulated by the distance ionising photons can propagate into the IGM. This depends on the abundance of photon sinks (absorption systems such as Lyman limit systems) and the corresponding recombinations of these systems. When the H II regions start approaching the typical separation of the photon sinks, an increasing fraction of ionising photons are used to balance recombinations, and the EoR can slow down (e.g. Furlanetto & Oh, 2005; Furlanetto & Mesinger, 2009; Alvarez & Abel, 2012).

The details of this process can be complicated (e.g. Sobacchi & Mesinger 2014); however a common simplification in semi-numerical approaches is to adopt a maximum horizon for the ionising photons within the ionised IGM, $R_{\rm mfp}$777 For historical context, we adopt ‘$R_{\rm mfp}$’ to denote this effective horizon set by sub-grid recombinations. The sub-grid models of Sobacchi & Mesinger (2014) do not translate to a constant, uniform value of $R_{\rm mfp}$, and further the mean free path is typically larger than $R_{\rm mfp}$ (since the former is an instantaneous quantity while the latter depends on the cumulative contributions of sub-grid recombinations). However, these authors find that the approximation of a constant $R_{\rm mfp}$ can reproduce the 21 cm PS of their fiducial model to $\sim 10$ per cent. Future versions of 21CMMC will include the Sobacchi & Mesinger (2014) sub-grid recombination model, eliminating $R_{\rm mfp}$ as a free parameter. , which is implemented as the maximum filtering scale in the excursion-set EoR modelling (see Section 2.1). Motivated by recent sub-grid recombination models (Sobacchi & Mesinger, 2014), we adopt a flat prior over $R_{\rm mfp}\in[5,25]$ cMpc.

In panel (ii) of Figure 1, we show the impact of $R_{\rm mfp}$ over cosmic history. As mentioned above, $R_{\rm mfp}$ only becomes important in the advanced stages of reionization, when the typical H II region scale approaches $R_{\rm mfp}$. The result is a delay of the late stages of the EoR for small values of $R_{\rm mfp}$. Limiting the photon horizon with a decreasing $R_{\rm mfp}$ ($\lesssim 15$ Mpc) is also evidenced by a drop in the large-scale 21 cm PS (e.g. McQuinn et al., 2007; Alvarez & Abel, 2012; Mesinger et al., 2012; Greig & Mesinger, 2015). Values of $R_{\rm mfp}>15$ Mpc have little impact on the 21 cm PS (for this combination of the other astrophysical parameters), as the clustering of the ionising sources becomes the dominant source of power (note the same PS amplitude for $R_{\rm mfp}=15$ and 25 Mpc).

2.2.3 Minimum virial temperature of star-forming haloes, $T^{\rm min}_{\rm vir}$

We define the minimum threshold for a halo hosting a star-forming galaxy to be its virial temperature, $T^{\rm min}_{\rm vir}$, which is related to the halo mass via, (e.g. Barkana & Loeb, 2001)

[TABLE]

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_{c}=18\pi^{2}+82d-39d^{2}$ where $d=\Omega^{z}_{\rm m}-1$. The choice of $T^{\rm min}_{\rm vir}$ acts as a step-function cut-off to the UV luminosity function. Below $T^{\rm min}_{\rm vir}$, it is assumed that internal feedback mechanisms such as supernova or photo-heating suppress the formation of stars. Above, efficient star formation overcomes internal feedback, enabling these haloes to produce ionising photons capable of contributing to the EoR and EoH. We shall consider a flat prior across $T^{\rm min}_{\rm vir}\in[10^{4},10^{6}]$ K within this work.

Our lower limit, $T^{\rm min}_{\rm vir}\approx 10^{4}$ K is motivated by the minimum temperature for efficient atomic line cooling. In principle, $T^{\rm min}_{\rm vir}$ can be as low as $\approx 10^{2}$ K in the presence of radiative cooling (Haiman et al., 1996; Tegmark et al., 1997; Abel et al., 2002; Bromm et al., 2002), however, star formation within these haloes is likely inefficient (a few stars per halo; e.g. Kimm et al. 2017) and can quickly ($z>20$) be suppressed by Lyman-Werner or other feedback processes well before the EoR (Haiman et al., 2000; Ricotti et al., 2001; Haiman & Bryan, 2006; Mesinger et al., 2006; Holzbauer & Furlanetto, 2012; Fialkov et al., 2013). Our upper limit of $T^{\rm min}_{\rm vir}\approx 10^{6}$ K, is roughly consistent with the host halo masses of observed Lyman break galaxies at $z\sim 6$–8, as estimated with the abundance matching technique (e.g. Kuhlen & Faucher-Giguère 2012; Barone-Nugent et al. 2014).

Note that within the 21CMFAST framework, $T^{\rm min}_{\rm vir}$ is important both in the EoR and EoH, determining the ionisation field ($f_{\rm coll}$, equation 2) and the specific emissivity requisite for the X-ray heating background ($\frac{df_{\rm coll}}{dt}$, equation 4) respectively. This implies that the efficient star-forming galaxies responsible for reionisation are the same galaxies which host the sources responsible for X-ray heating (i.e. the physics of star formation drives both the X-ray heating and ionisation fields). As a result, $T^{\rm min}_{\rm vir}$ affects the timing of both the EoR and EoH.

This is evident in panel (iii) of Figure 1. The EoR and EoH milestones are pushed to lower redshifts for an increasing $T^{\rm min}_{\rm vir}$. For a suitably large choice of $T^{\rm min}_{\rm vir}$ (e.g. $5\times 10^{5}$ K), the timing of the EoH and WF-coupling peaks can be sufficiently delayed to result in their overlap. In addition to the timing, $T^{\rm min}_{\rm vir}$, also impacts the amplitude of the fluctuations. Rarer (more biased) galaxies, corresponding to larger values of $T^{\rm min}_{\rm vir}$, are evidenced by more 21 cm power at a given epoch (e.g. McQuinn et al., 2007).

2.2.4 Integrated soft-band luminosity, $L_{\rm X\,<\,2\,keV}/{\rm SFR}$

The efficiency of X-ray heating is driven by the total integrated soft-band ($<$ 2 keV) luminosity per SFR (equation 5) escaping into the IGM, which normalises the emergent specific emissivity produced by the X-ray sources within the first galaxies. Within this work, we adopt a flat prior over the range ${\rm log}_{10}(L_{\rm X\,<\,2\,keV}/{\rm SFR})\in[38,42]$. This range is conservatively selected888 It is also roughly consistent with the limits proposed by Fialkov et al. (2017) using observations of the unresolved cosmic X-ray background (e.g. Lehmer et al., 2012) and the upper limits on the measured 21 cm PS from PAPER-64 (e.g. Ali et al., 2015; Pober et al., 2015). to be one to two orders of magnitude broader than the distribution seen in local populations of star forming galaxies (Mineo et al., 2012a), and their stacked Chandra observations (Lehmer et al., 2016). It also encompasses values at high-redshifts predicted by population synthesis models (Fragos et al., 2013b).

In Panel (iv) of Figure 1 we highlight the impact of $L_{\rm X\,<\,2\,keV}/{\rm SFR}$. For very high values of $L_{\rm X\,<\,2\,keV}/{\rm SFR}$, the EoH commences prior to the completion of WF-coupling. As a result, no strong absorption feature in $\bar{\delta T_{\rm b}}$ is observed, and the Ly$\alpha$-EoH peaks in the 21 cm PS merge. In addition to heating, such high X-ray luminosities can also substantially ionise the EoR (at the $\sim 10$-20 per cent level). In this case, the EoR can complete earlier (top panel).

At the other end of the range, extremely in-efficient X-ray heating (dashed curve) results in a delayed EoH, with the EoR and EoH overlapping. In this ‘cold reionisation’ scenario (Mesinger et al., 2013), reionisation proceeds in an IGM which is significantly colder than the CMB. The resultant large temperature contrasts between the ionised H II regions and cold neutral IGM can yield an extremely large 21 cm PS amplitude ($\sim 10^{3}$ mK2).999Though the specific model shown here exceeds the PAPER-64 upper limits, for other choices of the EoR and EoH parameters, this same choice of $L_{\rm X\,<\,2\,keV}/{\rm SFR}=10^{38}$ erg s${}^{-1}\,M^{-1}_{\odot}$ yr can result in a 21 cm PS below these upper limits.

2.2.5 X-ray energy threshold for self-absorption by the host galaxies, $E_{0}$

The soft X-rays produced within galaxies can be absorbed by the intervening ISM, thus not being able to escape and contribute to heating the IGM. The impact that this ISM absorption has on the emergent X-ray SED depends on the ISM density and metallicity. Early galaxies responsible for the EoH are expected to be less polluted by metals than local analogues. Indeed, Das et al. (2017) find that the emergent X-ray SED from simulated high-redshift galaxies can be well approximated by a metal free ISM with a typical column density of ${\rm log_{10}}(N_{\rm\mathrm{H\,{\scriptscriptstyle I}}{}}/{\rm cm^{2}})=21.40\begin{subarray}{c}+0.40\\ -0.65\end{subarray}$. As the opacity of metal free gas is a steep function of energy, these authors find that the 21 cm PS from the emerging X-ray SED can be well approximated using the common assumption of a step function attenuation of the X-ray SED below an energy threshold, $E_{0}$. In this work, we adopt a flat prior over $E_{0}\in[0.1,1.5]$ keV which corresponds to ${\rm log_{10}}(N_{\rm\mathrm{H\,{\scriptscriptstyle I}}{}}/{\rm cm^{2}})\in[19.3,23.0]$101010The conversion to column densities is computed assuming a unity optical depth for a pristine, metal free, neutral ISM (i.e. the contribution from He II is negligible)..

Panel (v) of Figure 1 highlights the impact of $E_{0}$. As $L_{\rm X\,<\,2\,keV}/{\rm SFR}$ defines the total soft-band luminosity, as we increase (decrease) $E_{0}$ we effectively harden (soften) the spectrum of emergent X-ray photons. Since the absorption cross section scales as $\propto E^{-3}$, we would naively expect smaller values of $E_{0}$ to result in more efficient heating, shifting the minimum in the global signal to higher redshifts. However, this evolution is slightly reversed for very low values, E_{0}\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{<}}0.5 keV, since the energy of these very soft photons is continually deposited in a limited volume surrounding the first galaxies; thus the volume-averaged global signal during the EoH is slightly delayed for these very soft SEDs. This highlights the need to properly model the spatially-dependent IGM heating, even when predicting the average signal.

The amplitude of the PS is very highly dependent on $E_{0}$. Softer SEDs result in very inhomogeneous heating, with PS amplitudes larger by up to an order of magnitude (e.g. Pacucci et al. 2014). The amplitude of the EoH peak consistently decreases for an increasing $E_{0}$, as the harder SEDs make the EoH more homogeneous. Eventually for $E_{0}>0.7$ keV, no EoH peak occurs, as the large mean free path of the X-ray photons result in an inefficient, relatively uniform heating of the IGM (e.g. Fialkov et al., 2014).

2.2.6 X-ray spectral index, $\alpha_{\rm X}$

The spectral index, $\alpha_{\rm X}$, describing the emergent spectrum from the first galaxies hosting X-ray sources depends on what is assumed to be the dominant process producing X-ray photons. In this work, we adopt a fiducial flat prior of $\alpha_{\rm X}\in[-0.5,2.5]$ which should encompass the most relevant high energy X-ray SEDs that describe the first galaxies (e.g. HMXBs, host ISM, mini-quasars, SNe remnants etc.; see for example McQuinn 2012; Pacucci et al. 2014).

Finally, in panel (vi) of Figure 1, we illustrate the impact of $\alpha_{\rm X}$. For an increasing $\alpha_{\rm X}$, more soft X-ray photons are produced (as the soft-band luminosity is kept fixed) resulting in more efficient X-ray heating. This results in an increase in the temperature fluctuations, driving a larger amplitude of the EoH peak in the 21 cm PS. Conversely, for an inverted spectral index, $\alpha_{\rm X}=-0.5$, the emergent X-ray photons are spectrally harder, producing inefficient uniform X-ray heating, and limited temperature fluctuations. Thus, the EoH peak in the 21 cm PS for $\alpha_{\rm X}=-0.5$ (dashed curve) is suppressed. However, over our adopted range of priors, the absorption by the host galaxy, $E_{0}$, is much more potent in hardening/softening the SED than the spectral index, $\alpha_{\rm X}$.

2.3 Telescope noise profiles

In order to provide astrophysical parameter forecasting, we must model the expected noise of the 21 cm experiments. Within this work, we focus solely on the 21 cm PS. To generate the sensitivity curves for 21 cm interferometer experiments, we use the python module 21cmsense111111https://github.com/jpober/21cmSense(Pober et al., 2013, 2014). Below, we summarise the main aspects and assumptions required to produce telescope noise profiles.

Firstly, the thermal noise PS is calculated at each gridded $uv$-cell according to the following (e.g. Morales, 2005; McQuinn et al., 2006; Pober et al., 2014),

[TABLE]

where $X^{2}Y$ is a cosmological conversion factor between the observing bandwidth, frequency and comoving distance, $\Omega^{\prime}$ is a beam-dependent factor derived in Parsons et al. (2014), $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}$. 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}$ (Thompson et al., 2007).

The sample variance of the cosmological 21 cm PS can easily be combined with the thermal noise to produce the total noise PS using an inverse-weighted summation over all the individual modes (Pober et al., 2013),

[TABLE]

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 (mock observation). Here we assume Gaussian errors for the cosmic-variance term, which is a good approximation on large-scales.

The largest primary uncertainty for 21 cm experiments is dealing with the bright foreground contamination. However, for the most part these bright foregrounds are spectrally smooth and have been shown to reside within a confined region of cylindrical 2D $k$-space known as the ‘wedge’ (Datta et al., 2010; Vedantham et al., 2012; Morales et al., 2012; Parsons et al., 2012b; Trott et al., 2012; Hazelton et al., 2013; Thyagarajan et al., 2013; Liu et al., 2014a, b; Thyagarajan et al., 2015b, a; Pober et al., 2016). Outside of this ‘wedge’ we are left with a relatively pristine observing window where the cosmic 21 cm signal is only affected by the instrumental thermal noise. At present, the location of the boundary separating this observing window and the ‘wedge’-like feature is uncertain.

Within 21cmsense three foreground removal strategies are provided (Pober et al., 2014), “optimistic”, “moderate” and “pessimistic”. We defer the reader to this work for further details on these scenarios, highlighting here that we choose to adopt the “moderate” scenario. This entails 21 cm observations only within the pristine 21 cm window (i.e. avoiding the “wedge”), with the wedge location defined to extend $\Delta k_{\parallel}=0.1\,h$ Mpc*-1* beyond the horizon limit (Pober et al., 2014). Furthermore, this scenario includes the coherent summation over all redundant baselines within the array configuration allowing the reduction of thermal noise (Parsons et al., 2012a).

Within this work, we focus on the two second generation 21 cm interferometer experiments capable of simultaneously measuring the EoR and the EoH, namely the SKA and HERA. Below (and in Table 1) we summarise the specific design features and assumptions required to model the theoretical noise of both instruments:

HERA: We follow the design specifics outlined in Beardsley et al. (2014) with a core design consisting of 331 dishes121212Note, the final fully funded design for HERA consists of 350 dishes, 320 in the core and 30 outriggers (DeBoer et al., 2017). For all intents and purpose for this work, the difference between a 320 and 331 core layout is negligible.. Each dish is 14m in diameter closely packed into a hexagonal configuration to maximise the total number of redundant baselines (Parsons et al., 2012a). We model the total system temperature as $T_{\rm sys}=100+T_{\rm sky}~{}{\rm K}$. HERA will operate in a drift-scanning mode, for which we assume a total 1080 hr observation, spread across 180 nights at 6 hours per night.

• 2.

SKA: We use the latest131313 In Greig & Mesinger (2015) we used the original design prior to the re-baselining (the 50 per cent reduction in the number of antennae dipoles), while in Greig et al. (2015) we investigated several design layouts to maximise the sensitivity specifically for the 21 cm PS following the re-baselining. The latest design for the SKA results in reduced sensitivity when compared to both our previous works, therefore we caution comparisons between the SKA within this work and our previous studies. available design for SKA-low Phase 1, using the telescope positions provided in the most recent SKA System Baseline Design document141414http://astronomers.skatelescope.org/wp-content/uploads/2016/09/SKA-TEL-SKO-0000422_02_SKA1_LowConfigurationCoordinates-1.pdf. Specifically, SKA-low Phase 1 includes a total of 512 35m antennae stations randomly distributed within a 500m core radius. 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}$. For the SKA we adopt a single, deep 1000 hr tracked scan of a cold patch on the sky.

It is non trivial to perform a like-for-like comparison between the two experiments, as HERA intends to perform a rotational synthesis drift scan, whereas the SKA intends to conduct track scanned observations151515Additionally, HERA is a dedicated 21 cm experiment specifically designed for a 21 cm PS measurement, while the SKA is a multidisciplinary experiment, with detection of the 21 cm signal only one of many key science goals. Moreover, the instrument layout and design is tailored towards 3D tomographic imaging (e.g. Mellema et al., 2013) rather than the 21 cm PS.. These two strategies result in considerably different noise PS. A single tracked field with the SKA will have considerably lower thermal noise than HERA owing to the deeper integration time, therefore SKA will be superior on small scales (large $k$) important for imaging. On the other hand, by observing numerous patches of the sky rotating through the zenith pointing field of view of HERA per night, sample variance can be better mitigated compared to the single tracked field of the SKA (i.e. HERA will have reduced noise on large scales, small $k$). For the most part, the strongest constraints on the astrophysical parameters come from the large-scale (small $k$) modes of the 21 cm PS, therefore using the PS as the likelihood statistic will favour the approach of HERA over that of SKA. We quantify these claims further in the Appendix, showing the noise power spectra at different redshifts.

Further complicating a direct comparison is the fact that the SKA is planning a tiered survey. For simplicity, in this introductory work we only consider a single, deep 1000 hr observation (though we will return to this in the future). However, we could have considered either the intermediate $10\times 100$ hr or wide and shallow $100\times 10$ hr strategies (e.g. Greig et al., 2015). These latter two surveys concede thermal noise sensitivity (from the reduced per field integration time) for an increased sample variance sensitivity by surveying multiple fields. In terms of the simplified three parameter EoR model considered in Greig & Mesinger (2015), the single, deep 1000 hr observation recovered the largest uncertainties on the astrophysical parameters relative to the intermediate or wide and shallow surveys (i.e. a single, tracked field was the worst performed strategy).

3 Mock cosmic signal

Having outlined our astrophysical model to describe the EoR and EoH, we now introduce our mock observations of the cosmic 21 cm signal. It is impractical to vary all available astrophysical parameters when creating mock observations, therefore, following Mesinger et al. (2016) we take two extreme choices for $T^{\rm min}_{\rm vir}$. This parameter characterises both the timing of the epochs and the typical bias of the dominant galaxies, thus encoding the largest variation in the 21 cm signal (e.g panel (iii) of Figure 1). Specifically, we adopt $T^{\rm min}_{\rm vir}=5\times 10^{4}$ K (faint galaxies) and $T^{\rm min}_{\rm vir}=3\times 10^{5}$ K (bright galaxies). These choices approximately match the evolution of the cosmic SFR density inferred from extrapolating the observed luminosity functions161616 In the future, we will provide a more generalised parameterisation of the ionising source model by allowing each individual constituent of $\zeta$ (equation 6) to vary with halo mass (Park et al., in prep.). This will enable additional flexibility in the source modelling, while at the same time allowing high-$z$ galaxy luminosity functions to be applied as priors to the source model prescription. (Bouwens et al., 2015) down to a UV magnitude of $M_{\rm UV}=-10\,(-17)$ for the faint galaxies (bright galaxies) galaxies model (see figure 5 of Das et al. 2017).

In order to match the latest constraints on the electron scattering optical depth, $\tau_{\rm e}$, from Planck ($\tau_{\rm e}=0.058\pm 0.012$; Planck Collaboration XLVII 2016), we adopt $\zeta=30$ for the faint galaxies model. This corresponds to the fiducial parameter set outlined in equation (6). For the biased, rarer galaxies in the bright galaxies model, we adopt $\zeta=200$. We select a fiducial $R_{\rm mfp}=15$ Mpc for both models. In Table 2 we summarise the adopted astrophysical parameter set for each of the two models, while also providing the corresponding optical depth, $\tau_{\rm e}$.

We adopt the same X-ray source model for both the bright galaxies and faint galaxies models. We assume HMXBs to be the dominant X-ray heating source within the first galaxies, and use the results of Das et al. (2017) to describe the emergent X-ray SED171717 Note that the typical star-forming haloes used by Das et al. (2017), e.g. $\lesssim 3\times 10^{8}\,M_{\odot}$, are smaller than the corresponding masses set by our minimum on $T^{\rm min}_{\rm vir}$ (e.g. $10^{4}$ K). Therefore, we have assumed that the intrinsic attenuation is independent of the host galaxy mass. This, however, need not be the case. . This corresponds to an energy threshold, $E_{0}=0.5$ keV (${\rm log_{10}}(N_{\rm\mathrm{H\,{\scriptscriptstyle I}}{}}/{\rm cm^{2}})=21.5$) and an X-ray spectral index of $\alpha_{X}=1.0$. Finally, we assume a HMXB soft band luminosity of $L_{\rm X\,<\,2\,keV}/{\rm SFR}=10^{40}$ erg s${}^{-1}\,M^{-1}_{\odot}$ yr, which is consistent with estimates from the HMXB population synthesis models of Fragos et al. (2013b).

In Figure 2 we provide the cosmic evolution of our two EoR source models: faint galaxies (solid curve) and bright galaxies (dashed curve). The top, middle and bottom panels are the global reionisation history, global 21 cm signal, and the 21 cm PS at $k=0.15$ (black) and $k=0.5$ Mpc*-1* (red) respectively. These show qualitatively the same behaviour as we have discussed in the previous section. These mock observations are generated from simulations with a volume of 6003 Mpc3, on a 4003 grid smoothed down from the high-resolution initial conditions generated on a 24003 grid.

By construction, these two EoR models have a similar reionisation history, consistent with the latest observations. The faint galaxies model has a somewhat more extended EoR, as the less biased DM halos which host the dominant galaxies form slower. Unlike for the observationally-constrained EoR history, we allow the EoH history to be different in the two models. In other words, we take the same X-ray luminosity per SFR. As a result, the faint galaxies galaxy model has an earlier EoH (see the middle panel), governed by the formation of the first $T^{\rm min}_{\rm vir}=5\times 10^{4}$ K structures.

Due to this, these two models will be interesting for our astrophysical parameter forecasting for both HERA and the SKA. With a decreasing sensitivity with increasing redshift, the location of the EoR and EoH peaks will impact the resultant significance of a detection of the 21 cm signal. For example, in the bottom panel of Figure 2, it is evident that the bright galaxies model exhibits additional 21 cm power relative to the faint galaxies model. Furthermore, this occurs at a lower redshift, where the instrumental noise is expected to be lower. Naively, one would therefore expect astrophysical parameters to be more tightly constrained with a mock bright galaxies signal.

4 21 cm forecasts

We now quantify astrophysical parameter constraints for each of these two models, for both second generation interferometers, HERA and the SKA. We first summarise our 21CMMC configuration, and then provide parameter constraints while discussing their implications.

4.1 21CMMC setup

As in Greig & Mesinger (2015), we use the 21 cm PS as the likelihood statistic to sample the astrophysical parameter space in 21CMMC. We adopt a modelling uncertainty of 20 per cent on the sampled 21 cm PS (not the mock observation of the 21 cm PS), which is added in quadrature with the total noise power spectrum (equation 8).181818A modelling uncertainty accounts for the inaccuracy of semi-numerical approaches such as 21CMFAST. Here we simply take an uncorrelated fixed percentage error, with a value consistent with the comparisons in Zahn et al. (2011). This constant error purely accounts for the observed fractional differences in the recovered 21cm PS between the simulations. For example, it does not account for any potential errors arising from the different algorithms used to compute the IGM spin temperature. Future work will seek to characterise this uncertainty, potentially even mediating it by comparing 21CMFAST to suites of RT simulations. As in Greig & Mesinger (2015) we compare the 21 cm PS at each redshift only over the reduced $k$-space range, $k=0.15\,$–$\,1.0$ Mpc*-1*, which we deem to be free of cosmic variance and shot-noise effects arising within the simulations.

We perform our forecasting using eight co-evolution191919The observed 3D 21 cm signal is spatially dependent on the 2D sky, with the third (line-of-sight) direction being frequency dependent. The line-of-sight axis encodes evolution along the light-cone, which can impact the observed 21 cm PS, when compared to predictions from co-evolution cubes (e.g. Datta et al., 2012, 2014; La Plante et al., 2014; Ghara et al., 2015). However, this effect is generally only pronounced on much larger scales than considered here (note that our lowest $k$-modes, $|k|=0.15$ Mpc*-1*, correspond to a modest $\Delta z\sim 0.1$ at EoR redshifts). We shall nevertheless return to this in future work, extending 21CMMC to operate directly on the light-cone.

redshifts ($z=6,\,7,\,8,\,9,\,11,\,13,\,15$ and 17). Note that our choice is relatively arbitrary, taken to span the EoH and EoR in both of our mock observations and is not driven by any computational or numerical reasons.

For the MCMC sampling, we use smaller boxes than used for the mock observations: 3003 Mpc3 volume on a 2003 grid smoothed down from a high-resolution 12003 grid. Both the mock observations and the sampled 21 cm PS have the same voxel resolution ($\sim 1.5$ Mpc per voxel). These box and cell sizes were selected to ensure per cent level convergence in the PS for a few randomly sampled astrophysical models.

Finally, it is instructive to provide some timing estimates for 21CMMC. We note that the public version of 21CMFAST is not optimised for our MCMC framework. Therefore, we heavily streamlined the computation for a single core implementation (since the advantage of 21CMMC is that a realisation of 21CMFAST is performed for each available CPU core). For the $300^{3}$ Mpc3 volume, 2003 voxel setup above, our streamlined version of 21CMFAST is $\sim 20\times$ faster on a single core202020 This boost in computational efficiency arises owing to the inclusion of additional interpolation tables to remove redundant calculations. Additionally, 21CMFAST is reduced to a single executable removing file I/O and limiting memory overheads. . For our astrophysical parameter forecasting, we perform $\sim 8\times 10^{4}$ 21CMFAST runs. As we are using the parallelised affine invariant ensemble sampler emcee, we can achieve $\sim 8\times 10^{4}$ samplings with the following setup: 400 walkers, each performed for 200 iterations212121We explored several combinations of the walker/iterations configuration to confirm our experimental setup was providing converged results for our MCMC.. Performed using 200 cores (i.e. 100 physical cores + 100 virtual cores) on a shared memory cluster, such a setup takes $\sim 8$ days.

4.2 Parameter recovery for the faint galaxies mock observation

In Figure 3 we present the astrophysical parameter constraints for our faint galaxies model for our assumed 1000hr observation with both HERA (blue) and the SKA (red). Across the diagonals we provide the 1D marginalised probability distribution functions (PDFs) for each of the six model astrophysical parameters. In the lower left half of the figure, we provide the 2D marginalised joint likelihood contours, with crosses denoting the input fiducial values, while the thick and thin contours denote the 68 (1$\sigma$) and 95 (2$\sigma$) percentiles, respectively. In the top right half of the figure, we provide the 2$\sigma$ marginalised constraints on the IGM neutral fraction, $\bar{x}_{\mathrm{H\,{\scriptscriptstyle I}}{}}$, with respect to the reionisation history of the mock observation (solid black curve). Note that to generate this figure, we interpolate the reionisation history between the marginalised distributions for $\bar{x}_{\mathrm{H\,{\scriptscriptstyle I}}{}}$ at the eight co-evolution redshifts, sampling at a rate of $\Delta z=0.4$. In Table 3 we provide the median and associated 16th and 84th percentiles for each of our six astrophysical parameters for both HERA and the SKA.

From the relatively narrow 1D PDFs, it is clear that both HERA and the SKA can simultaneously constrain the EoR and EoH to high accuracy. The only parameter for which we do not achieve strong constraints is the X-ray spectral index, $\alpha_{X}$. However, this is not overly surprising given the relatively small effect it has on the amplitude of the 21 cm PS over the entire allowed parameter range (see Figure 1). Folding in the 20 per cent modelling uncertainty and the instrumental noise, the relative difference in the 21 cm PS amplitude across the full allowed parameter range is roughly consistent at about 2$\sigma$. In terms of the reionisation history both HERA and the SKA recover comparably tight constraints on the IGM neutral fraction. At $1\sigma$, we recover constraints on $\bar{x}_{\mathrm{H\,{\scriptscriptstyle I}}{}}$ of the order of $\sim 5$ per cent.

We recover no strong degeneracies with our astrophysical parameters, with mild degeneracies for $\zeta$-$T^{\rm min}_{\rm vir}$, $L_{\rm X\,<\,2\,keV}/{\rm SFR}$-$T^{\rm min}_{\rm vir}$ and $E_{0}$-$\alpha_{X}$. This is generally consistent with both Ewall-Wice et al. (2016) and Kern et al. (2017), whose mock observation most closely resembles our faint galaxies model. However, these authors find somewhat stronger $f_{\rm X}$-$E_{0}$ and $f_{\rm X}$-$\alpha_{X}$ degeneracies (where $f_{\rm X}$ corresponds to the number of X-ray photons per stellar baryon, and can thus be related to our $L_{\rm X\,<\,2\,keV}/{\rm SFR}$ for a given $\alpha_{X}$ and $E_{0}$). This discrepancy could arise due to (i) the inclusion of the 20 per cent modelling uncertainty which broadens our contours; (ii) the approximations made in those works [assumptions of Gaussian errors in Fisher matrices (Ewall-Wice et al., 2016) or modelling errors in an emulator method Kern et al. (2017)]; and/or (iii) our choice for the soft-band energy as a normalisation parameter (instead of the X-ray photon number), which can provide a more independent basis vector for the EoH evolution222222 In selecting the soft-band X-ray luminosity instead of a harder X-ray band (e.g. 0.5 - 8 keV), we have preferentially minimised the degeneracy between the X-ray luminosity and $\alpha_{X}$. Additionally, adopting a soft-band X-ray luminosity enables straightforward comparison with numerous observations of nearby galaxies (e.g. Tzanavaris & Georgantopoulos, 2008; Mineo et al., 2012b; Fragos et al., 2013a; Lehmer et al., 2015, 2016). In the near future, we expect observations of the intrinsic soft-band X-ray luminosity escaping the host galaxy to improve with the upcoming Athena telescope (Barcons et al., 2012), which will provide a soft-band effective area more than an order of magnitude larger that existing experiments (T. Dauser, private communication). This will provide stronger priors on the X-ray SED, even if 21cm observations themselves are less discriminatory. .

If we assume the 1D marginalised PDFs can be modelled by a normal distribution, which for the most part is reasonable at the 1$\sigma$ level (i.e. some tails begin to appear at $\sim 2\sigma$), we can provide some approximate fractional uncertainties for the astrophysical parameters. For the SKA (HERA), the 1$\sigma$ percentage errors are: $\zeta$ = 18 (24), $R_{\rm mfp}$ = 16 (16), log10$(T^{\rm min}_{\rm vir})$ = 1.4 (2.3), log10($L_{\rm X\,<\,2\,keV}/{\rm SFR}$) = 0.2 (0.3), $E_{0}$ = 17 (14) and $\alpha_{X}$ = 88 (73). The uncertainty on the EoR parameters is comparable to what we obtained in Greig & Mesinger (2015): $\zeta$ = 17 (22), $R_{\rm mfp}$ = 18 (18), log10$(T^{\rm min}_{\rm vir})$ = 2.4 (3.3). Therefore, despite increasing the model complexity by including the EoH, the relative constraints are comparable.

Ewall-Wice et al. (2016) quote fractional precisions on their six parameter model, with 1-2 per cent accuracy on the EoR parameters and 6 per cent on their EoH parameters. Their constraints are smaller than ours by about an order of magnitude for the EoR and factor of a few for the EoH. Approximately half of this difference can be attributed to the inclusion of the modelling uncertainty (see e.g. Greig & Mesinger, 2015). The remaining discrepancy can arise from either the fundamental assumptions in their Fisher matrix approach, or their larger number of redshift samples (more than a factor of two). In future, we will modify 21CMMC to directly work on the observed light cone, removing the necessity of an ad-hoc sampling of co-evolution cubes.

For our faint galaxies model, we find both HERA and the SKA will recover comparable parameter constraints. This is despite the significantly increased sensitivity achievable with the SKA, resulting in a larger total integrated signal to noise (S/N). As pointed out in Greig & Mesinger (2015), the S/N is not a reliable metric for predicting an instrument’s ability at parameter constraints, since model constraining power is biased towards large scales. This highlights the importance of using parameter forecasting as a figure of merit, instead of just the total S/N. We caution however that the SKA performance can be improved if one can better mediate modelling uncertainties. Indeed, the increased thermal noise sensitivity on small-scales is washed out by our assumed 20 per cent modelling uncertainty. Moreover, the SKA will be superior at tomography; using higher-order likelihood statistics should therefore favour the SKA over HERA.

4.3 Parameter recovery for the bright galaxies mock observation

In Figure 4 we present our 1 and 2D joint marginalised posterior distributions for each of the six astrophysical parameters for our bright galaxies model assuming a 1000hr observation with HERA (blue) and the SKA (red). Table 4 provides the median and associated 16th and 84th percentiles for each of our astrophysical model parameters. As in the previous section, we provide approximate fractional precisions on the model parameters assuming normally distributed marginal likelihoods. For the SKA (HERA), the 1$\sigma$ percent errors are: $\zeta$ = 17 (15), $R_{\rm mfp}$ = 16 (12), log10$(T^{\rm min}_{\rm vir})$ = 0.4 (0.6), log10($L_{\rm X\,<\,2\,keV}/{\rm SFR}$) = 0.2 (0.2), $E_{0}$ = 16 (17) and $\alpha_{X}$ = 80 (79).

The constraints are comparable for most of the astrophysical parameters held fixed across the two models (i.e. $R_{\rm mfp}$, log10($L_{\rm X\,<\,2\,keV}/{\rm SFR}$), $E_{0}$, $\alpha_{X}$). The largest difference is in $R_{\rm mfp}$, which shows tighter constraints relative to the faint galaxies model. At a given stage in the EoR, the H II regions in the bright galaxies model are larger and more isolated, sourced by the brighter, rarer, more-biased sources. Recombinations play a larger role in limiting the growth of such H II regions, whose characteristic sizes approach $R_{\rm mfp}$ earlier in reionisation. Since the signal is more sensitive to $R_{\rm mfp}$, a degeneracy between $\zeta$ and $R_{\rm mfp}$ emerges, with both now able to control the timing of the EoR: one can compensate for the slowing-down of the EoR due to recombinations (a smaller $R_{\rm mfp}$) by increasing the ionising efficiency, $\zeta$. This degeneracy leads to poorer overall constraints on $\zeta$ when compared to the faint galaxies model.

These improvements in the constraints of $R_{\rm mfp}$ are however only available with HERA. The marginalised PDFs for the SKA exhibit a noticeable tail towards increasing $R_{\rm mfp}$. The source of this tail can simply be attributed to the reduced large-scale sensitivity of the SKA relative to HERA owing to our adopted observing strategies within this work (see Section 2.3 and Figure 10).

5 Can we ignore the spin temperature in EoR parameter recovery?

The majority of studies that constrain the EoR with the cosmic 21 cm signal assume that the IGM spin temperature is saturated (i.e. $T_{\rm S}\gg T_{\rm CMB}$). Since the corresponding term in the brightness temperature saturates to unity in this regime (c.f. equation 1), assuming saturation greatly simplifies the computational load. However, the validity of such an approximation is highly dependent on the poorly-constrained relative efficiencies of ionising and X-ray sources in the first galaxies. Not properly taking into account the IGM spin temperature can bias the recovered EoR parameter constraints232323 In Ewall-Wice et al. (2016), these authors found that by incorrectly accounting for the EoH, their recovered uncertainties in the EoR parameters could be biased. They explored this by fitting their EoR parameters by either (i) assuming a fixed EoH model; or (ii) properly fitting and marginalising over the EoH model parameters. The fractional uncertainties from (i) were smaller than (ii) by almost a factor of two. This emphasises the importance of properly accounting for the EoH when performing astrophysical parameter recovery from the 21 cm signal. . Here, we quantify the impact of this using our two mock 21 cm observations. The two adopted extrema in $T^{\rm min}_{\rm vir}$ allow the exploration of vastly different levels of overlap between the EoR and EoH, enabling us to estimate the available span in the corresponding bias. Here we use 21CMMC to recover the fractional precision on the EoR parameters under the simplification of $T_{\rm S}\gg T_{\rm CMB}$, and compare these to the constraints obtained by properly modelling $T_{\rm S}$.

5.1 Faint galaxies

The 21 cm signal for the faint galaxies model transitions from absorption to emission at $z\lesssim 12$ (see e.g. the middle panel of Figure 2). At lower redshifts, the IGM spin temperature continues to increase and the $(1-T_{\rm CMB}/T_{S})$ factor approaches unity. Assuming the saturated limit ($T_{\rm S}\gg T_{\rm CMB}$) results in a fractional error in the power spectrum less than 10% when $(1-T_{\rm CMB}/T_{S})^{2}\gtrsim 0.9$, corresponding to $z\sim 7.3$ and an IGM neutral fraction of $\bar{x}_{\mathrm{H\,{\scriptscriptstyle I}}{}}\sim 0.4$ (cf. the similar faint galaxies model considered in Mesinger et al. 2016). Therefore, the saturated spin temperature approximation is only reasonable during the second half of reionisation in our faint galaxies model.

In Figure 5, we present the 1 and 2D joint marginalised posterior distributions for the EoR parameters, namely $\zeta$, $R_{\rm mfp}$ and log${}_{10}(T^{\rm min}_{\rm vir})$. For this comparison, we only consider the SKA (the results for HERA are nearly identical as shown in Table 3), and run 21CMMC on the redshifts spanning the EoR: $z=6,\,7,\,8$ and 9. Red curves correspond to the constraints on the EoR model parameters from fitting the full six parameter astrophysical model (i.e. including the EoH, Figure 3) whereas the purple curves correspond to the constraints obtained with the saturated spin temperature approximation. In the lower half of Table 3 we provide the median, 16th and 84th percentiles for the recoverer uncertainties on the EoR parameters.

It is clear from Figure 5 that the saturated spin temperature approximation significantly biases parameter constraints. The saturated temperature approximation prefers models with a smaller virial temperature and ionising efficiency compared to the ‘true’ ones. The marginalised 1D PDFs for $T^{\rm min}_{\rm vir}$ are offset by at least $\sim 3\sigma$. For $\zeta$, we recover more modest offsets of $\sim 1.5\sigma$, and for $R_{\rm mfp}$ we recover comparable constraints. In all cases, the recovered astrophysical constraints return larger uncertainties, owing to the lower marginalised likelihoods. For reference, the maximum likelihood (ML: $\mathcal{L}={\rm exp}(-\frac{1}{2}\chi^{2})$) is a factor of four lower than obtained with the full spin temperature modelling in Section 4.2.

To understand this bias, in Figure 6 we show the 21 cm PS for the mock observation (black curve), the ML 21 cm PS assuming $T_{\rm S}\gg T_{\rm CMB}$ (red curve) and the 21 cm PS adopting the ‘true’ EoR model parameters assuming $T_{\rm S}\gg T_{\rm CMB}$ (blue curve). Error bars denote the 20 per cent modelling uncertainty on the model 21 cm PS. The main impact of ignoring the spin temperature can be seen by comparing the blue and black curves. Since both curves are generated from the same EoR model parameters, any discrepancies arise from the spin temperature. The largest discrepancy, as expected, arises at the highest redshift, where the IGM is still undergoing the final stages of X-ray heating. As discussed in Mesinger et al. (2016), during the EoR the cosmic H I patches effectively have the same, uniform IGM spin temperature, resulting in fairly negligible temperature fluctuations for the faint galaxies model (see also Pober et al., 2015). However, the amplitude of the 21 cm PS is decreased by a factor of $(1-T_{\rm CMB}/T_{\rm S})^{2}$ (compare the blue and black curves at $z\sim$ 8–9). In order to mimic this $(1-T_{\rm CMB}/T_{\rm S})^{2}$ decrease in amplitude, a model which assumes $T_{\rm S}\gg T_{\rm CMB}$ will tend to prefer EoR models with intrinsically less power, i.e. those in which the sources are less biased, having a lower $T^{\rm min}_{\rm vir}$. A lower $T^{\rm min}_{\rm vir}$ implies more abundant ionising sources, which must be compensated for by decreasing $\zeta$ to attempt to recover the correct reionisation history. Therefore, the likelihood peaks at smaller virial temperatures and ionising efficiencies (c.f. the red curve, corresponding to $T^{\rm min}_{\rm vir}=10^{4.29}$ K and $\zeta=15.5$). The full evolution of the PS is unable to be completely reproduced, producing the factor of four lower ML compared to the complete model including the EoH242424Note, that the ML from our six parameter model was computed from eight different redshifts, whereas in the saturated limit we only considered four redshifts. As additional redshift data decreases the unnormalised ML value of a model, the relative differences here are even larger..

5.2 Bright galaxies

Already for the faint galaxies model we find significant biases in the recovered EoR parameters under the assumption that the spin temperature is saturated. In the bright galaxies model, the EoH and EoR overlap even more strongly, in effect maximising this bias in parameter recovery. We quantify this in Figure 7, in which we provide the marginalised distributions, and in the lower half of Table 4 where we provide the median, 16th and 84th percentiles for our recovered EoR model parameters.

For both $\zeta$ and $T^{\rm min}_{\rm vir}$ we recover marginalised constraints discrepant at $>10\sigma$. As in the faint galaxies model, $R_{\rm mfp}$ remains consistent to within $1\sigma$, but it instead prefers marginally lower values, $R_{\rm mfp}\sim 10$ Mpc. In the inset of Figure 7, the recovered reionisation history is discrepant at the $>5\sigma$ level, beyond $z\sim 7.5$. In Figure 8 we present the 21 cm PS from the mock bright galaxies model, the ML estimate and the ‘true’ EoR model parameters assuming the saturated spin temperature limit. It is immediately evident that all the constraining power for the ML arises from the 21 cm PS at $z=6$ and 7. At $z>7$, the saturated limit cannot reproduce the reduced amplitude of the mock 21 cm PS252525For reference, the ML is $\sim 10^{28}$ times lower when we assume $T_{\rm S}\gg T_{\rm CMB}$..

This is not surprising. As discussed in Section 4.3, the global averaged 21 cm brightness temperature contrast is still in absorption at $z>8$ (see Figure 2), indicating $T_{\rm S}\lesssim T_{\rm CMB}$ for a significant fraction of the simulation volume. At the same time, the IGM is 35 per cent ionised ($\bar{x}_{\mathrm{H\,{\scriptscriptstyle I}}{}}=0.65$) by $z=8$, indicating reionisation is well underway. This significant overlap of the EoR and EoH breaks the fundamental assumption of the saturated spin temperature limit. Furthermore, at $z=8$ the bright galaxies model is transitioning closely to the $T_{\rm S}\equiv T_{\rm CMB}$ limit, producing a precipitous drop in the 21 cm PS amplitude. Relative to the bright galaxies model in the saturated limit (blue curve) this is a factor of $\sim 40$ difference in the 21 cm PS amplitude. Under the saturated limit assumption, there is no avenue to mimic such behaviour, resulting in hugely discrepant astrophysical parameter constraints.

This highlights the importance of properly including the EoH. Incorrectly ignoring the EoH and associated IGM spin temperature fluctuations when interpreting a realistic observation could significantly bias the inferred EoR source model.

6 Conclusion

Detecting the cosmic 21 cm signal during the EoR and the EoH stands to reveal insights into the formation, growth and evolution of structure in the Universe. However, how do we interpret the underlying astrophysics once we have a detection? To aid this, we developed 21CMMC, a massively parallel Bayesian MCMC analysis tool, which performs full 3D reionisation simulations (using 21CMFAST) on the fly, for recovering EoR astrophysical parameter constraints. Second generation experiments such as HERA and the SKA, along with global 21 cm experiments, will be able to measure both the EoR and the pre-heating of the IGM by X-rays. Therefore, in order to facilitate simultaneous astrophysical forecasting of the EoR and EoH and consequently aid in the development and construction of the data reduction and signal extraction pipelines, in this work we extend 21CMMC into the epoch of X-ray heating.

We demonstrate that both HERA and the SKA will be able to simultaneously constrain the astrophysics of reionisation and the X-ray heating. We consider two models describing the mock observation: faint galaxies and bright galaxies. These are intended to encompass the physically plausible region of parameter space provided by extrapolating the faint end of the observed UV luminosity function.

Assuming a 1000 hr observation of the 21 cm PS at eight co-evolution redshifts we recover the fractional precision on our six parameter model describing the sources responsible for the EoR and the EoH. These parameters include an ionising source efficiency ($\zeta$), effective photon horizon ($R_{\rm mfp}$), minimum halo mass of ionising sources ($T^{\rm min}_{\rm vir}$), soft-band luminosity ($<2$ keV) of X-ray sources ($L_{\rm X\,<\,2\,keV}/{\rm SFR}$), minimum energy threshold for the attenuation of the X-rays by the host galaxy ISM ($E_{0}$) and the spectral index of the X-ray source SED ($\alpha_{X}$). Additionally assuming a 20 per cent modelling uncertainty in the power spectrum, we recover the parameters of the mock signals with the following percentage error (1$\sigma$):

• •

faint galaxies: $\zeta$ = 18 (24), $R_{\rm mfp}$ = 16 (16), log10$(T^{\rm min}_{\rm vir})$ = 1.4 (2.3), log10($L_{\rm X\,<\,2\,keV}/{\rm SFR}$) = 0.2 (0.3), $E_{0}$ = 17 (14) and $\alpha_{X}$ = 88 (73), respectively for the SKA (HERA).

• •

bright galaxies: $\zeta$ = 17 (15), $R_{\rm mfp}$ = 16 (12), log10$(T^{\rm min}_{\rm vir})$ = 0.4 (0.6), log10($L_{\rm X\,<\,2\,keV}/{\rm SFR}$) = 0.2 (0.2), $E_{0}$ = 16 (17) and $\alpha_{X}$ = 80 (79), respectively for the SKA (HERA).

Both the SKA and HERA perform equally as well at simultaneously constraining the astrophysics of reionisation and the epoch of X-ray heating due to their comparable sensitivities on the large scales (under our assumptions regarding the survey strategy) which most strongly discriminate between the astrophysical models.

With our expanded framework, we also quantify the impact of the common assumption of a saturated spin temperature, $T_{\rm S}\gg T_{\rm CMB}$, during the EoR. Our faint galaxies model has a relatively distinct EoH and EoR, for which we would typically expect $T_{\rm S}\gg T_{\rm CMB}$ to be a reasonable approximation during the bulk of reionisation. Nevertheless, even this modest overlap of epochs leads to biases in $\zeta$ and $T^{\rm min}_{\rm vir}$ of up to $\sim 3\sigma$, and a reduction in the ML by a factor of $>4$. The bright galaxies model on the other hand represents the extreme case of EoR and EoH overlap. For this model, the recovered constraints are discrepant at $>10\sigma$, with the ML under in the saturated limit being a factor of $\sim 10^{28}$ lower than in the full model. Therefore, adopting the saturated spin temperature approximation can significantly bias inferences on EoR parameters.

Acknowledgements

We thank Adrian Liu and Nicholas Kern for comments on a draft version of this work. This work was supported by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 638809 – AIDA – PI: Mesinger).

Appendix A Sensitivity Curves

In Figures 9 and 10 we present the mock 21 cm PS (black curves) for both of our EoR source models faint galaxies and bright galaxies, respectively. We show the 21 cm PS for each of the eight co-evolution redshifts we use in this work to achieve our astrophysical constraints, covering both the EoR and the EoH. For each, we provide the corresponding total noise PS for HERA (blue) and SKA (red) for our assumed 1000hr observational strategy (see Section 2.3 for more details). Finally, the black dashed curve is the 20 per cent modelling uncertainty, which we have applied to the mock 21 cm PS (it is applied to the sampled 21 cm PS not the mock observation in 21CMMC) to provide reference to the dominate source of error at any $k$-mode.

It is interesting to note that while most of the EoH parameters are tightly constrained for both mock observations, these constraints appear to be arise solely from a relatively narrow range in redshifts, which pick up the rise and fall of the large scale power. This emphasises the importance of observing the cosmic 21 cm signal in an extended frequency range, allowing us to pick up the major milestones in the signal which drive strong astrophysical parameter constraints.