Bayesian Model Selection with Future 21cm Observations of The Epoch of Reionisation

TL;DR

This paper uses Bayesian statistics and the Multinest algorithm to evaluate the ability of current and future 21cm observations to distinguish between different models of the Epoch of Reionisation, highlighting observational requirements and challenges.

Contribution

It applies Bayesian model selection to EoR scenarios with 21cm data, assessing observational capabilities and analyzing parameter redundancies within models.

Findings

LOFAR unlikely to enable model selection even with long observations

HERA needs 61 dipoles over 1080 hours for model selection

SKA can distinguish models with only 324 hours of data

Abstract

We apply Bayesian statistics to perform model selection on different reionisation scenarios via the Multinest algorithm. Initially, we recover the results shown by 21CMMC for the parameter estimation of 21cmFAST models. We proceed to test several toy models of the Epoch of Reionisation (EoR) defined in contrasting morphology and scale. We find that LOFAR observations are unlikely to allow model selection even with long integration times. HERA would require 61 dipoles to perform the same analysis in 1080 hours, and becomes comparable to the SKA with 217 dipoles. We find the SKA requires only 324 hours of observation to conclusively distinguish between our models. Once model selection is achievable, an analysis of observational priors is performed finding that neutral fraction checks at specific redshifts add little to no inference. We show the difficulties in model selection at the level of distinguishing fiducial parameters within a model or distinguishing galaxies with a constant versus power law mass-to-light ratio. Finally, we explore the use of the Savage-Dickey density ratio to show the redundancy of the parameter Rmfp within 21cmFAST.

Figures

Tables

Table 1. Table 1: A summary of the different telescope specifications used in this paper. Note that of the SKA’s 512 stations we only simulate using the ‘central’ 296 stations, where sensitivity to the EoR frequencies is dominant.

Parameter	LOFAR	HERA	SKA-Central
Number of dipole	48	331	296
Stations
Station Diameter [m]	31	14	35
Collecting Area [$m^2$]	35,762	50,953	492,602
$T_R$[K]	140	100	$40+\frac{T_{\mathrm{sky}}}{10}$
Observing Range	[110,240]	[50,250]	[50,350]
[MHz]
Observation Time	1080	1080	1080
[hrs]
Scan Type	1-hr Track	Drift	1-hr Track

Table 2. Table 2: This table summarises the 4 main EoR models considered in this work, see Watkinson & Pritchard ( 2014 ) for more details. It is worth noting here that filtering (F) the MHR models changes the scale of reionisation from local to global .

Model	Physical motivation	How is ionisation defined?	Resulting type of reionisation
FZH	Over-dense regions collapse to form stars, UV radiation is therefore the driving force behind the ionised bubbles. The size of these bubbles dictates the neutral fraction and therefore the brightness temperature.	Over-density increases above the critical barrier (${\delta }_{\mathrm{crit}}$).	Global Inside-Out
Inverted (Inv) FZH	Over-dense regions remain neutral due to a high recombination rate, hence reionisation begins in underdense regions. These then grow and dominate as hard radiation becomes the dominant radiation background and the IGM becomes ionised.	Over-density decreases below the critical barrier.	Global Outside-In
MHR	Stars are formed as above, ionising the circum-galactic gas. Ionisation then proceeds a long underdense regions because the atoms recombine faster than they are ionised in the dense regions. Eventually the background radiation is dense enough to dominate.	A pixel-by-pixel implementation of an under density threshold.	Local Outside-In
Inverted (Inv) MHR	Stars are formed as above, ionising the circum-galactic gas. Over-dense regions become ionised. It is easy to implement rather than being physically viable since even dense regions that do not have a line of sight to a radiation source will become ionised.	A pixel-by-pixel implementation of an over density threshold.	Local Inside-Out

Table 3. Table 3: The MAP parameter values to 1 ​ σ 1 𝜎 1\sigma (standard deviation) for the toy models fit to all three observational prior checks (without 21cm data). Note that the large errors (compared to each respective parameter) hint that this data will not be constraining in the context of distinguishing toy EoR scenario (discussed further in Section 6.4 ).

	$\zeta$	$R_{\mathrm{mfp}}$	${\log }_{10}\left[T_{\mathrm{vir}}\right]$
FZH	95.4 $\pm$ 61.9	12.3 $\pm$ 4.3	5.17 $\pm$ 0.29
Inv FZH	71.8 $\pm$ 68.5	14.9 $\pm$ 5.4	4.72 $\pm$ 0.34
F MHR	413 $\pm$ 289	9.3 $\pm$ 4.3	5.80 $\pm$ 0.31
F Inv MHR	268 $\pm$ 290	18.3 $\pm$ 4.4	5.61 $\pm$ 0.34
MHR	404 $\pm$ 293	-	5.72 $\pm$ 0.33
Inv MHR	799 $\pm$ 291	-	5.88 $\pm$ 0.33

Table 4. Table 4: A summary of the uniform prior ranges used for the 21cmFAST and toy models’ Evidence values. For a detailed description of these models and their parameters see Section 3.1 .

Model	$\zeta$	$R$	$\mathrm{Log}\left[T_{\mathrm{vir}}\right]$	$\alpha$
3pFZH	[5.,250.]	[5.,20.]	[4.0,5.3]	-
4pFZH	[5.,250.]	[5.,20.]	[4.0,5.3]	[-3.,3.]
InvFZH	[5.,5000.]	[0.1,10.]	[4.,7.]	-
MHR	[5.,4000.]	-	[4.,7.]	-
Inv MHR	[5.,1200.]	-	[4.,6.]	-
F MHR	[5.,1000.]	[0.1,10.]	[4.,6.]	-
F Inv MHR	[5.,1200.]	[0.1,10.]	[4.,6.]	-

Table 5. Table 5: The MAP parameter values to 1 ​ σ 1 𝜎 1\sigma for the 5 simulations performed in Section 6.1 (two of which are shown in Figure 5 , with fiducial results marked on). ζ eff subscript 𝜁 eff \zeta_{\rm eff} is the actual ionisation efficiency that the galaxies produce, via the implementation of the ζ 𝜁 \zeta parameter (via Equation 22 ). Note that in the 3pFZH cases ( α = 0 𝛼 0 \alpha=0 ) use the relationship of Equation 4 , while for 4pFZH we use Equation 6 ( α ≠ 0 𝛼 0 \alpha\neq 0 ). For the fiducial values, f1 and f2, σ ζ subscript 𝜎 𝜁 \sigma_{\zeta} has been approximated by using 25% of each value as suggested in the original 21cmFAST paper (Mesinger et al., 2011 ) ; this is motivated by the discrepancies between semi-numerical and full radiative transfer simulations. Deviations from recovering the exact fiducial values are due to the observational priors inputted in the likelihood.

	$\zeta$	${\log }_{10}\left[T_{\mathrm{vir}}\right]$	$\alpha$	${\zeta }_{\mathrm{eff}}$
3pFZHf1	18.7 $\pm$ 2.0	4.366 $\pm$ 0.09	(0.0)	18.7
4pFZHf1	5.8 $\pm$ 27.6	4.148 $\pm$ 0.16	0.36 $\pm$ 0.52	6.6
f1	20. $\pm$ 5.	4.447 $\pm$ 1.1	(0.0)	20.
3pFZHf2	34.9 $\pm$ 3.2	4.747 $\pm$ 0.05	(0.0)	34.9
3pFZHf2	12.2 $\pm$ 4.2	4.576 $\pm$ 0.13	(0.4)	66.1
4pFZHf2	5.3 $\pm$ 27.3	4.473 $\pm$ 0.15	0.624 $\pm$ 0.40	7.8
f2	15. $\pm$ 3.8	4.699 $\pm$ 1.2	0.40 $\pm$ 0.10	39.6

Table 6. Table 6: The MAP values ( ± 1 ​ σ plus-or-minus 1 𝜎 \pm 1\sigma ) for the toy model simulations shown in Figures 8 , 9 , and 10 when observing for 1080 hours with HERA-331.

Model	$\zeta$	R	$\mathrm{Log}\left[T_{\mathrm{vir}}\right]$
FZH	18.8 $\pm$ 2.0	14.4	4.38 $\pm$ 0.09
Inv FZH	1680 $\pm$ 1080	1.05	6.16 $\pm$ 0.1
MHR	1690 $\pm$ 470		6.02 $\pm$ 0.04
Inv MHR	71. $\pm$ 15		5.12 $\pm$ 0.07
F MHR	432 $\pm$ 80	4.65	5.56 $\pm$ 0.06
F Inv MHR	203 $\pm$ 170	4.35	5.3 $\pm$ 0.1

Equations

$$\delta T_{\mathrm{b}}$$

$$\times {\left(\frac{1+z}{10}\frac{0.15}{{\Omega }_{\mathrm{M}}{h}^2}\right)}^{\frac{1}{2}}\left(\frac{{\Omega }_{\mathrm{b}}{h}^2}{0.023}\right)mK.$$

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

$$f_{\mathrm{coll}}=\frac{1}{{\rho }_M}{\int }_{M_{\mathrm{vir}}}^{\infty }m\frac{dn}{dm}dm,$$

$$\zeta =$$

$${\zeta }_0=30\left(\frac{f_{\mathrm{esc}}}{0.3}\right)\left(\frac{f_{\star }}{0.05}\right)\left(\frac{N_{\gamma }}{4000}\right)\left(\frac{2}{1+n_{\mathrm{rec}}}\right)$$

$$\zeta =\{\begin{array}{ll}{\zeta }_0{\left(\frac{T_{\mathrm{vir}}}{T_{\mathrm{vir}}^{\mathrm{feed}}}\right)}^{\alpha },& \text{if }{T}_{\mathrm{vir}}\ge T_{\mathrm{vir}}^{\mathrm{feed}}\\ 0,& \text{otherwise}.\end{array}$$

$${\Delta }_{21}(k{)}^2$$

$${\delta }_{21}(\mathbf{x},z)\equiv \frac{\delta T_b(\mathbf{x},z)}{\overline{\delta T_b}(z)}-1.$$

$${\Delta }_{\mathrm{N}}^2(\mathbf{k})\approx X^{2}Y\frac{{\mathbf{k}}^3}{2{\pi }^2}\frac{{\Omega }_{\mathrm{Eff}}}{2t_{\mathrm{int}}}{T}_{\mathrm{sys}}^2.$$

$$m_{\mathrm{ion}}=\zeta m_{\mathrm{gal}}.$$

$$f_{\mathrm{coll}}=\text{erfc}\left\{\frac{{\delta }_{\mathrm{crit}}-\delta }{\sqrt[]{2[{\sigma }^2(m_{\min })-{\sigma }^2(m)]}}\right\}.$$

$$\delta \ge {\delta }_{\mathrm{crit}}-\sqrt[]{2[{\sigma }_{\min }^2(m_{\min })-{\sigma }^2(m)]}{\text{erfc}}^{-1}\left(1-\frac{1}{\zeta }\right).$$

$$f_{\mathrm{coll}}^{\prime }=$$

$$=$$

$$\frac{i}{N_{\mathrm{p}}}=1-{\bar{x}}_{\mathrm{HI}}.$$

$$P(\mathbf{\Theta }|D,M)=\frac{P(D|\mathbf{\Theta },M)P(\mathbf{\Theta }|M)}{P(D|M)},$$

$${\mathcal{B}}_{12}=\frac{P(D|M_1)}{P(D|M_2)}=\frac{{\mathcal{Z}}_1}{{\mathcal{Z}}_2}.$$

$${\mathcal{Z}}_i=\int P(D\mid {\Theta }_1,...,{\Theta }_N,M_i)P({\Theta }_1,...,{\Theta }_N|M_i)d^N\Theta$$

$$=\int \mathcal{L}(\mathbf{\Theta })\Pi (\mathbf{\Theta })d\mathbf{\Theta }$$

$$={\int }_0^1\mathcal{L}(X)dX.$$

$$\displaystyle\mathcal{B}(\Theta_{*})\equiv\frac{\mathcal{Z}(\Theta_{*})}{\mathcal{Z}}{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}..............}$$

$$\text{vspace}0.2cm=\frac{\int \mathcal{L}(\mathbf{\Theta })\Pi (\mathbf{\Theta })d{\Theta }^{N-1}}{\mathcal{Z}}=\int \mathcal{P}(\mathbf{\Theta })d{\Theta }^{N-1},$$

$${\zeta }_{\mathrm{eff}}=\frac{{\int }_{M_{\mathrm{vir}}}^{\infty }dm\frac{dn}{dm}m\zeta (m)}{{\int }_{M_{\mathrm{vir}}}^{\infty }dm\frac{dn}{dm}m}$$

Full text

1 Introduction

The Epoch of Reionisation (EoR) is the most recent phase change of the Universe. During this process, the hydrogen (HI) in the intergalactic medium (IGM) that has not formed stars and galaxies becomes heated by these very objects and eventually becomes ionised. We expect to uncover astrophysics involving radiative transfer, star formation, X-ray heating and the percolation of ionised hydrogen (HII) bubbles. Due to a lack of observation however, we remain in the dark.

An absorption signal in the Cosmic Microwave Background (CMB) radiation at 78 MHz was claimed to be observed for the first time by the EDGES team (Bowman et al., 2018) in February 2018. Around the same time an upper bound of $(79.6~{}\rm mk)^{2}$ on the EoR power spectrum has been achieved by the LOFAR team with a 13 hour observation of the NCP field within the redshift range $z=9.6-10.6$ at $k=0.053~{}h~{}\rm c\rm Mpc^{-1}$ (Patil et al., 2017).

Although these limits on the power spectrum are still several orders of magnitude larger than the signal, this paves the way for the first observations of the EoR with a new generation of radio interferometer. The Square Kilometre Array (SKA) will be capable of tomographically mapping primordial hydrogen from $z=6$ to $27$ (Mellema et al., 2013; Koopmans et al., 2015). The aim of this work is to explore the application of Bayesian model selection in the context of EoR models to explore the possibility for ruling out specific reionisation scenarios.

The redshifted 21cm line is set to be a powerful tool for analysing the high redshift universe (Furlanetto, 2006a). This atomic line originates from the spin-flip transition of the coupled magnetic moments between the proton and the electron in a hydrogen atom, releasing a 21cm (lab frame) wavelength photon which is then redshifted by the Hubble flow. The CMB acts as the initial backlight for the transition. By observing the signal at different frequencies, a tomographical image of how the primordial hydrogen evolves can be constructed. This methodology allows us to trace the formation of the first stars and galaxies. For a full discussion of how physically informative the 21cm signal could be, see e.g. Pritchard & Loeb (2012).

The current cosmological picture is bound by two sources of observational data: The CMB, at $z\sim 1100$ and high-redshift quasars z\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\sim}\hss}\raise 1.0pt\hbox{<}}7. The Big Bang model suggests that protons and electrons form separately. By $z\sim 1090$ the expanding universe has cooled enabling $99.9999\%$ of the universe’s baryonic matter to form hydrogen atoms (Peebles, 1968; Zeldovich & Sunyaev, 1969). Currently this is thought to remain the case until stars start to form $\sim 150\rm Myr$ after the Big Bang. The latter source of data includes measuring the Gunn-Peterson trough (Gunn & Peterson, 1965) within quasar spectra. This shows that by $z\sim 6$, $99.96\pm 0.03\%$ of the hydrogen in the IGM is ionised (Fan et al., 2006). In between these two milestones the universe re-ionises the hydrogen left over from recombination. A defining feature of EoR toy models is the morphology of reionisation (Furlanetto et al., 2004; Miralda-Escudé et al., 2000). In this work we consider two possibilities:

• •

Inside-Out - Where low density gas is ionised last, due to star formation being most dominant in dense regions of gas. The density field and radiation field correlate with position.

• •

Outside-In - Where high density gas is ionised last, due to the recombination rate of dense gas dominating over radiation. The density field and radiation field anti-correlate with position (in the IGM).

The other is the scale on which these behaviours are associated with:

• •

Local - Implementation of the ionisation criteria is performed pixel by pixel in the simulation. i.e. with no dependency on the surrounding area.

• •

Global - A scale dependency is included when implementing the ionisation criteria. Typically the surrounding region is influenced up to the mean free path of a photon.

These descriptions are approximate. In reality reionisation will progress with a variety of these behaviours at different times and scales, as is observed in detailed simulations (Finlator et al., 2018). For a detailed description of the astrophysics involved in the EoR, please see e.g. Loeb & Furlanetto (2013) or Dayal & Ferrara (2018).

Data is now becoming available from instruments such as LOFAR (Patil et al., 2017), MWA (Dillon et al., 2015), PAPER (Ali et al., 2015), and HERA (DeBoer et al., 2017); all of which hope to detect the 21cm power spectrum (Chapman et al., 2012). In the meantime mock data is created by simulations. The convention for finding the best fit parameters for these models is done using Markov Chain Monté Carlo (MCMC) algorithms (Hassan et al., 2017; Greig & Mesinger, 2015). Our sampler of choice is called Multinest (Feroz et al., 2009), a nested sampling algorithm. This subset of the MCMC family focuses on calculating the Bayesian Evidence, a quantity which measures how well a model fits a data set - see e.g. Sivia & Skilling (2006) for a comprehensive introduction to MCMC statistical analyses. Unfortunately full radiative transfer simulations are too computationally expensive for the number of samples required for an MCMC analysis111Typically $\sim 10^{4}$ samples guarantees convergence in this context. (Zahn et al., 2011). Instead semi-numerical simulations, such as 21cmFAST (Mesinger & Furlanetto, 2007; Mesinger et al., 2011) are used. These use approximate techniques, often combined with observations, to save on computation time - whilst producing consistent results to the full radiative transfer codes (Fialkov et al., 2012; Santos et al., 2010). Interestingly the application of machine learning algorithms, such as Convolutional Neural Networks, are beginning to show promise in reducing the computation times of these simulations even further (Gillet et al., 2018; Shimabukuro & Semelin, 2017; La Plante & Ntampaka, 2018; Kern et al., 2017; Hassan et al., 2019). With the first light on the SKA due around 2024, next generation radio interferometers are sure to reveal many secrets from the unseen universe. We are entering into an age of precision astrophysics and cosmology, particularly if synergies between telescopes such as SPICA222http://www.spica-mission.org/ and JWST333https://www.jwst.nasa.gov/ can be achieved. Sophisticated analysis techniques, such as Bayesian Inference, are becoming a necessity.

This paper is structured as follows: Section 2 introduces the base of this work, 21CMMC, including a summary of the recipes used for simulating the 21cm data; Section 3 highlights 21cmNest - the additions made to 21CMMC including the toy models; Section 4 is an introduction to the Bayesian Evidence and how it can be used for model selection. The results are split into seven subsections: Initially, we show agreement between the 21CMMC and Multinest in parameter estimation (Section 6.1), including cross checks of the methodology; then we perform model selection with LOFAR, HERA and SKA; in the context of HERA we successfully distinguish between EoR scenarios (Section 6.3) and show the redundancy of average neutral fraction checks as a way of ruling out toy models (Section 6.4); Finally, we attempt to distinguish two parameterisations of 21cmFAST (Section 6.7) - including the use of nested models in establishing redundant parameters with the Savage-Dickey density ratio (Section 6.7.1).

The main aims of this work look at the feasibility of ruling out toy models with the three selected observatories - These can be neatly summarised as the following three questions:

• •

How long would observations need to be to effectively distinguish between the toy models with LOFAR? (Section 6.2)

• •

Given a fixed observing timescale (1080 hours) how many dipoles are necessary to effectively distinguish between the toy models with HERA? (Section 6.5)

• •

What is the minimum observation time needed to effectively distinguish between the toy models with SKA? (Section 6.6).

We assume a $\Lambda$CDM cosmology throughout this work with the following parameters: ($\Omega_{M}$, $\Omega_{\Lambda}$, $\Omega_{b}$, $n_{S}$, $\sigma_{8}$, $H_{0}$) $=$ (0.27, 0.73, 0.046, 0.96, 0.82, 70 km s*-1* Mpc*-1*) (Planck Collaboration et al., 2016a).

2 The State of the Art: 21CMMC

21CMMC is a Bayesian parameter estimation tool for the EoR (Greig & Mesinger, 2015). It combines Cosmohammer (Akeret et al., 2013) - a framework for using the affine invariant MCMC algorithm Emcee (Foreman-Mackey et al., 2013); with 21cmFAST (Mesinger et al., 2011) - a numerical implementation of the analytic model described by Furlanetto, Zaldarraiga, and Hernquist (Furlanetto et al., 2004), referred to as FZH.

A streamlined version of 21cmFAST is used, where evolved matter density ($\delta$) and velocity ($dv_{\rm r}/dr$) fields are pre-calculated for the IGM. The 21cm signal intensity is directly proportional to the differential brightness temperature $\delta T_{\rm b}$ which can be approximated (Furlanetto et al., 2006) as

[TABLE]

where $x_{\rm HI}$ is the hydrogen neutral fraction; $T_{\rm CMB}$ is the CMB temperature; $T_{\rm s}$ is the spin temperature; and all other parameters follow the convention of a $\Lambda$CDM cosmology. In this original version of 21CMMC the simplifying assumption of $T_{s}\gg T_{\rm CMB}$ is used to ease computation (see Greig & Mesinger (2017a) for the inclusion of $T_{\rm s}$ and X-ray heating - which is also omitted here). This is valid during most of reionisation $z\sim[6,10]$ (Furlanetto, 2006b; Baek, S. et al., 2010; Chen & Miralda-Escudé, 2004) and is known as the post-heating regime.

This model uses the excursion set formalism in order to identify HII regions. A smoothing scale ($R$) is iterated down from the photon mean free path, $R_{\rm mfp}$ to the size of the cell. This parameter is called the photon mean free path for historical reasons. It is more accurately a photon’s mean horizon since only the photons whose path are stopped by ionising HI are important in the EoR - a simple way to account for recombinations in the ionised IGM444The newest 21CMMC version (Park et al., 2018) computes inhomogeneous recombinations explicitly.. What’s important is the number of ionising photons verses the number of HI per smoothing scale. For a more complete discussion please see Greig & Mesinger (2015, 2017a). The defined criteria for ionisation is

[TABLE]

with x being position and $\bar{M}_{\rm min}$ as the minimum virial halo mass defined by the atomic cooling threshold. The collapse fraction is dependent on the virial temperature through the virial mass. We assume the relationship derived in Barkana & Loeb (2001) where $M_{\rm vir}\sim T_{\rm vir}^{\frac{2}{3}}$. $f_{\rm coll}$ is the fraction of matter collapsed into dark matter halos (Press & Schechter, 1974; Bond et al., 1991; Sheth & Tormen, 1999) and is calculated by the following integral of the Press-Schechter halo mass function (a numerical correction factor to return on average the Sheth-Tormen result for ellipsoidal dark matter halo collapse is applied)

[TABLE]

where $\rho_{M}$ is the matter density of the chosen cosmology. $\zeta$ is the ionising efficiency of galaxies at high redshift, performed via a step function

[TABLE]

where $T_{\rm vir}$ is the virial temperature and

[TABLE]

to encompass general astrophysical properties of galaxies: A‘turn-on’ mass is defined by the aforementioned atomic cooling threshold, $T_{\rm vir}^{\rm feed}=10^{4}K$ (Barkana & Loeb, 2001); $N_{\gamma}$ is the number of ionising photon per baryon within a star, for $\rm popII$ stars only - $N_{\gamma}\approx 4000$ which is currently assumed to be the case during the EoR (Barkana & Loeb, 2005); $n_{\rm rec}$ is the typical number of recombinations per hydrogen atom, assumed to be $n_{\rm rec}\approx 1$ (Sobacchi & Mesinger, 2014) for a photon-starved end point for the EoR; $f_{\star}$ is the fraction of galactic gas in stars, which is taken to be $0.05$ (Ferrara & Loeb, 2013); and $f_{\rm esc}$ is the fraction of ionising photons escaping into the IGM. Since both $f_{\star}$ and $f_{\rm esc}$ are observationally uncertain (Gnedin et al., 2008; Wise & Cen, 2009; Ferrara & Loeb, 2013). We will return to these in Section 5, when we discuss the parameter priors.

21CMMC features the option to relax the single galaxy population assumed so far i.e. to depart from a constant $\zeta$. This is done by the addition of a 4th parameter, $\alpha$, enabling a power law relationship between $\zeta$ and $T_{\rm vir}$ as,

[TABLE]

During the MCMC analysis the likelihood is defined to be a $\chi^{2}$ statistic between a mock 21cm observed power spectra and a theoretically calculated one, which is combined with observational priors (2.1). The mock power spectrum uses a fiducial choice of parameters while the theoretical models’ best fit parameters are found via the MCMC analysis. This is done for $8$ k bins from a foreground corruption limit of $k=0.15~{}\rm Mpc^{-1}$ to a shot noise limit of $k=1.0~{}\rm Mpc^{-1}$ - this is maintained in this work unless specified otherwise.

The spherically averaged power spectrum, $\Delta_{21}^{2}$, is chosen to be the aforementioned likelihood statistic and is defined via the brightness temperature as

[TABLE]

The bar denotes a spatial average, while the angular brackets denote a k-space average. This likelihood computation chain is summarised in Figure 1.

2.1 Observational data priors

Depending on the choice of redshift range 21CMMC implements up to three observational checks:

• •

The Thomson optical depth of the IGM was observed by Planck to be: $\tau=0.058\pm 0.012$, assuming an instantaneous EoR (Planck Collaboration et al., 2016b). During each MCMC call, 21cmFAST can produce estimates of $\tau$ by interpolating the neutral hydrogen fraction across the desired redshifts. This allows testing between the observed and simulated values of $\tau$.

• •

The location of the Gunn-Peterson trough in QSO Ly-$\alpha$ forest data: McGreer et al. (2015) show that the EoR must be $>90\%$ complete by $z=5.9$. Therefore for all $z\leq 5.9$, $x_{\rm HI}=0$, and for $z>5.9$ the estimated neutral HI fraction is tested against an half-Gaussian with mean $\bar{x}_{\rm HI}=0.06$ and variance $\sigma^{2}=0.05$.

• •

The Red Ly-$\alpha$ damping wing allows the estimation of the $x_{\rm HI}$ surrounding QSO ULASJ1120+0641 (Greig et al., 2017). Provided $z=7.08$ is included in the redshift interpolation range, the neutral fraction at this redshift is checked against the observational measurement $\bar{x}_{\rm HI}=0.4^{+0.41}_{-0.32}$ at $2\sigma$.

All tests are performed with a $\chi^{2}$ between $x_{\rm HI}$ (and $\tau$) which are then combined into the likelihood linearly as: $ln\mathcal{L}=\chi^{2}_{21\rm cm}+\chi^{2}_{\rm Planck}+\chi^{2}_{\rm McGreer}+\chi^{2}_{\rm Greig}$.

2.2 Telescope sensitivity with 21cmSense

The Python code 21cmSense555https://github.com/jpober/21cmSense, 21cmSense uses heavily the Astrophysical Interferometry Python code AIPY https://github.com/AaronParsons/aipy. (Pober, 2016) is used throughout this work to calculate the telescope noise profiles. In each U-V bin the thermal noise is modelled as:

[TABLE]

Where $X^{2}Y$ converts the solid angle of observing bandwidth to a comoving distance, $\Omega_{\rm Eff}$ represents the effective beam solid angle as derived in Parsons et al. (2014). $t_{\rm int}$ is the total integration time (per observed frequency). $T_{\rm sys}$ is the sum of the sky temperature (Thompson et al., 2017) and the receiver temperature $T_{\rm R}$. Cosmic variance is included assuming it is a Gaussian error - we can express the total uncertainty by spherically averaging an inversely weighted summation across all the k modes. Hence in Sections 6.2 and 6.6 (when $t_{\rm int}$ is altered) we expect the telescope noise to scale roughly as $t_{\rm int}^{-0.5}$. The specifications of LOFAR are taken from (van Haarlem et al., 2013). We use the 48 core high-band antenna stations only. HERA’s dipoles are structured in a filled hexagon (Dillon & Parsons, 2016; Beardsley et al., 2015). For the SKA we follow the specifications of LOW Phase 1 with 512 stations666https://www.skatelescope.org/key-documents/. The middle 224 stations are randomly distributed in a circular core ($\sim$ 400m in radius) with the remaining 288 stations split in three spiral arms extending outwards. For all telescopes, unless stated otherwise, we assume 6 hours of observing per day for 180 days for $t_{\rm int}$. The cosmological bandwidth is assumed to be $8~{}\rm MHz$ throughout. The telescope details are summarised in Table 1. Fundamental concerns with detecting the 21cm signal are due to foregrounds eg. Pober (2015). We adopt the ‘Moderate’ foreground setting within 21cmSense777The other options are: Pessimistic - Where baselines are added incoherently and no k modes are included in the horizon or the buffer zone; or Optimistic - All k modes in the primary field of view are used.. This means that all baselines are added coherently but no k modes from the horizon or a buffer zone (chosen to be the default $0.1h\rm Mpc^{-1}$) are used. See e.g. Pober et al. (2013, 2014a, 2014b) for more detail.

3 21CMNest

The primary objective of this work888https://binnietom@bitbucket.org/binnietom/21cmnest_1.0.git is to apply Bayesian model comparison to EoR scenarios. This has been achieved with two alterations to the 21CMMC framework. Firstly, we add various toy models (detailed in Section 3.1); and secondly we replace the Emcee sampler (Foreman-Mackey et al., 2013) (from Cosmohammer) with pyMultinest (Section 4).

3.1 Flexibility and toy EoR models

Since there is still uncertainty in the detailed physics of the EoR we explore several toy models. These toy models are explained below. Initially we use the same four models as Watkinson & Pritchard (2014), summarised in Table 2. See Figure 2 for a slice of brightness temperature from each models’ box at $z=8$ for $\bar{x}_{\rm HI}=0.5$. Figure 3 shows the corresponding 21cm power spectra and EoR ionisation histories for these slices. Although Table 2 contains physical motivations for all the models, we stress here that only the FZH and MHR models are physically liable - the others are included to test the proof-of-concept for Bayesian model selection. Reionisation may begin with the local outside-in scenario but global inside-out reionisation will eventually dominate assuming UV radiation drives reionisation (Garaldi et al., 2019). Given we understand the fluctuations of the IGM density on linear (Mpc) scales, this is easy to show by comparing the recombination and ionisation rates with e.g. 21cmFAST. Table 3 shows the maximum a posteriori (MAP) parameters for the toy models fit to the Greig, McGreer and Planck priors detailed in Section 2.1. Before 21cm data is included the allowed EoR parameter space is ample (Greig & Mesinger, 2017b).

3.1.1 FZH

The 21cmFAST simulation is a numerical implementation of the model developed by Furlanetto, Hernquist and Zaldarriagan (FZH), as summarised in Section 2. This assumes999The simplest scenario (Barkana & Loeb, 2001). that the mass of the collapsed object is relatable to the mass of the ionised region by the ionising efficiency ($\zeta$) as follows,

[TABLE]

Using the Press-Schechter collapse fraction, $f_{\rm coll}$ (Equation 3) can be written as

[TABLE]

This implies that the barrier for ionisation can be rewritten as a density constraint instead:

[TABLE]

Note that the factor of 2 in Equations 10 and 11 comes from a symmetry argument. For every excursion that reaches the critical barrier, only half of these will collapse to form virialised objects due to the nature of random walks.

Numerically, FZH is applied by smoothing around each pixel from large to small scale. If Equation 11 is achieved before the filter scale reaches the pixel size is met the central pixel is flagged as ionised. Since this is a global prescription, not all regions are purely ionised or neutral. After the filtering scale has reached the pixel size 21cmFAST sets the neutral fraction of $x_{\rm HI}=1-\zeta f_{\rm coll-per-pixel}$ to each pixel that is yet to be fully ionised101010This includes rounding to $0~{}if~{}x_{\rm HI}<0$ and $1~{}if~{}x_{\rm HI}>1$..

3.1.2 Inverted (Inv) FZH

Inverting the ionisation threshold (Equation 11) allows us to produce an outside-in model sensitive to global scales. Within this model underdense regions will be ionised first to form bubbles, the edge of which will remain in equilibrium111111Similar to a Strömgren sphere but on galactic scales. until the background ionising radiation increases. Eventually this will dominate over the recombination rate in overdense regions anticipating the completion of reionisation. Firstly we redefine $\zeta$ as $\zeta^{\prime}$ in order to keep the parameter framework within 21CMMC unchanged. $\zeta^{\prime}$ is now the background ionising radiation efficiency. Next we must turn to the critical barrier definition, which we choose to be $\delta_{\rm crit}^{\prime}=-\delta_{\rm crit}$. Assuming Gaussian perturbations for our regions of over-density, as is done for the FZH models (see the form of Equation 3 becomes 10), all we need to do is calculate the fraction of matter that satisfies $\delta<\delta_{\rm crit}^{\prime}$. In other words a new collapse fraction $f_{\rm coll}^{\prime}$, below which HI structures are sparse enough to be ionised. This can be done by changing the limits of the original $f_{\rm coll}$ integral as is shown below

[TABLE]

Now we have an analogous ionisation threshold, $f_{\rm coll}^{\prime}\zeta^{\prime}\geq 1.$, that can be implemented easily within the 21CMMC framework. As before, any partially ionised pixels left after using the excursion set formalism are assigned $x_{\rm HI}=1-\zeta^{\prime}f_{\rm coll}^{\prime}$.

3.1.3 MHR

For an alternate method of implementing ionisation we turn to the work of Miralda-Escudé, Haenelt and Rees referred to as MHR (Miralda-Escudé et al., 2000). This allows us to a produce local outside-in model of reionisation - underdense regions of the IGM are ionised first due to a recombination rate that is dependent on the density of the gas. The idea is that the denser a region of gas, the faster it can recombine and therefore the dense regions of the IGM are ionised last. Circumgalactic gas, within the halo, is ionised first allowing ionising photons to seep into the IGM. HII regions then expand in the directions of low gas density and the EoR is defined to end when these regions overlap. Typically this is achieved with a background of X-ray photons, due to their large mean free paths. Firstly, to define the aforementioned ionisation threshold we prescribe a neutral fraction $\bar{x}_{\rm HI}=1-f_{\rm coll}\zeta$. Secondly, we order the pixels by HI density. We then ensure consistency between the $x_{\rm HI}$ by finding the $i$’th pixel (of total $N_{\rm p}$) that satisfies

[TABLE]

Once the i’th pixel has been identified, we can set the ionisation definition as a function of density to be $\delta<\delta_{i}$. In summary the only deviation from 21cmFAST is the definition of ionisation via a new density threshold. Prescribing the $\bar{x}_{\rm HI}$ directly is referred to as the 1-parameter (1p) MHR model.

3.1.4 Inverted (Inv) MHR

Here we produce a local inside-out model of reionisation. The procedure follows as above, but with the greater than sign reversed to produce an over-density threshold: $\delta>\delta_{j}$. The choice of the j’th pixel using: $j=N_{\rm p}-i$, ensures (again) that the neutral fraction remains consistent between simulations.

3.1.5 Filtered (F) MHR

In order to test 21cmNest more thoroughly any skewing effect from extra parameters (see Section 4.2) in the Bayesian Evidence, a third parameter is added to the basic MHR model. This takes the form of a Gaussian k-space filter performed across the density field. The integrated volume of this Gaussian equates to that of a real space top-hat with radius $R_{\rm Filter}$. This smoothing has a tendency to reduce the range of values within the density field, which can have a large effect on the number and density of ionised structures. We expect this added flexibility to relax the local scenario, allowing the filtered model to fit the scale of the fiducial model. This means producing a global model when fitting for an FZH fiducial power spectrum, but being a redundant parameter when fitting for an MHR fiducial power spectrum (Section 6.1.2). The F Inv MHR proves to be the most challenging toy model to test, since it shares morphology and scale with FZH. In the F MHR model the background ionising efficiency, $\zeta$, is pushed below that of the MHR model. With the opposite being true for the F Inv MHR. As for the filter scale itself, this tends to the mean sizes of HII bubbles.

4 Bayesian Inference

Bayes’ theorem can be written as

[TABLE]

where $\mathbf{\Theta}$ are the parameters of a given model $M$, and $D$ is the data set to be fit against given an hypothesis. In this context the hypothesis is implied via the choice of model under consideration. $P(\mathbf{\Theta}\mid D,M)$ is the posterior parameter distribution, $\mathcal{P}$; $P(D\mid\mathbf{\Theta},M)$ is the likelihood $\mathcal{L}$; the prior knowledge $\Pi$ is $P(\mathbf{\Theta}\mid M)$; and the Bayesian Evidence121212Alternatively known as the model likelihood or marginal likelihood., $\mathcal{Z}$ is defined to be $P(D\mid M)$. Bayesian inference for parameter estimation is centred around the likelihood, as is the case for Cosmohammer within 21CMMC. In most MCMC algorithms (including 21CMMC) the Evidence is omitted because normalising the posterior is a redundant overhead in parameter estimation and the integral in Equation 18 is typically difficult to evaluate.

We now turn to the Bayesian Evidence, $\mathcal{Z}$, as this takes centre stage for model selection. By evaluating a $\mathcal{Z}_{i}$ for the i’th model, $M_{i}$, we can define what is called the Bayes factor as follows:

[TABLE]

Bayes factors can be interpreted as the odds on which model is thought to be correct, given the data. We can expand $\mathcal{Z}_{i}$ as

[TABLE]

Nested sampling shortcuts this difficult calculation of $\mathcal{Z}_{i}$ by transforming the integral to one-dimension (Equation 19), the fraction of prior volume ($X$) which reduces per sampling iteration. For more specific information on nested sampling see (Skilling, 2004). We discuss the interpretation of the Bayes factor in Section 4.2. We choose this machinery over simpler options, such as the Akaike Information Criterion (Akaike, 1974) or Bayesian Information Criterion (Schwarz, 1978), due to their lack of consistency - e.g. on WMAP data (Magueijo & Sorkin, 2007; Liddle, 2007).

The use of Multinest over other nested sampling algorithms is due to the low number of dimensions in our selected models ($\leqslant 4$). As dimensionality increases the acceptance fraction geometrically approaches [math] until the ellipsoidal rejection sampling within Multinest becomes inefficient. At higher dimension ($\sim 10$), the competing algorithm Polychord (Handley et al., 2015) might be a better choice.

In this work, 2000 live points131313Live points are used to calculate iso-likelihood contours within nested sampling algorithms. were found to be enough to produce consistent results within the statistical error bars produced by Multinest. The Multinest-obtained posterior distributions plotted throughout this text have been checked for convergence using NestCheck (Higson et al., 2018)141414https://github.com/ejhigson/nestcheck see the documentation for more details.. Dynamic nested sampling algorithms such as Dynesty (Higson et al., 2017) result in a maximum live point number of 1500 (to two s.f.) when run in this context. We use this as confirmation that nothing has been missed by our choice of 2000 live points within Multinest. The errors plotted in the Bayes factor figures in this work are those produced by Multinest, arisen purely from the numerical uncertainties within the algorithm (e.g. using the trapezium rule) - see Skilling (2004) for more details. It is worth emphasising that these errors do not reflect the sensitivity to Prior choice or the observational measurements.

4.1 The Savage-Dickey density ratio

One convenient approach to analysing the degeneracy of parameters is via the Savage-Dickey density ratio (Dickey, 1971; Verdinelli & Wasserman, 1995). Since this can only be applied to models that are nested within each other we hold one parameter, $\Theta_{*}$, constant and calculate the Bayesian Evidence. In this notation $d^{N-1}\Theta$ refers to the differential $d\Theta^{N}$ without the $d\Theta_{*}$ term (e.g. $d\Theta_{1}...d\Theta_{N}=d^{N-1}\Theta d\Theta_{*}$). Comparing the Bayesian Evidence factor of each of these $\mathcal{Z}(\Theta_{*})$ to $\mathcal{Z}$ (as in Equation 17) gives us a parameter dependent Bayes factor

[TABLE]

which, as shown above (and in Section 6.7.1), reflects the shape of the 1-dimensional posterior for the parameter $\Theta_{*}$ doubling as a cross check for calculating the Evidence integral.

4.2 The Jeffreys’ scale and potential pitfalls of the Bayes factor

The largest counter argument to using Bayesian statistics is the possibility of skewing results from a naive choice of prior distribution151515See Jaynes (2003) for a very well illustrated example in which results exactly reflect choices of prior distribution.. One of the many attractive properties Bayesian model selection has to offer is its inherent inclusion of Occam’s Razor - more parameters should not be used than are necessary. Because of this if any parameter prior region is chosen naively the added hyper-volume can skew the result via Equation 18. We therefore try to keep to observationally constrained information, particularly when considering the prior ranges. By doing so we can approach this as a way of rewarding good parameters’ rather than penalising redundant ones (Nesseris & García-Bellido, 2013).

The Jeffreys’ scale is conventionally used as a rough guide when treating this dilemma - since large Evidences must be obtained before conclusions are drawn. We use the Jeffreys’ scale throughout this work to ease interpretation, keeping in mind that the analytic odds are the real conclusions. We adopt the following version (Hobson et al., 2009) with three levels of significance:

• •

Strong - $\mathcal{B}_{12}>150$ then model 1 has outperformed model 2 and is objectively better at describing the data in hand.

• •

Moderate - $10<\mathcal{B}_{12}<150$ then the two models are likely to be distinguishable by this method, but care must be taken in assuring no skew has been introduced.

• •

Weak - $\mathcal{B}_{12}<10$ then the two models are likely to be indistinguishable by this method.

5 Considerations of sensitivity to the parameter priors

Throughout this work we apply only uniform (uninformative) prior distributions across all parameter ranges. Log priors have been tested (not shown) to produce the same values for $\mathcal{Z}$ (within error bars) for the 3pFZHf1, 3pFZHf2 and 3pFZH(with $\alpha=0.4$)f2 models with the same ranges of prior. We follow similar assumptions of prior boundaries taken from the original 21CMMC. The only differences from the initial paper are at the higher end of the $\zeta$ parameter, which we increase from $[5.,100.]$ to $[5.,250.]$ for the 21cmFast models. The upper end of this parameter is under debate due to the lack of observational data surrounding $f_{\star}$ and $f_{\rm esc}$ (see Section 2). The original choice of $f_{\rm esc}$ is in the range $[0.05,1]$ implying that $f_{\star}=0.05$. This is motivation by a theoretical modelling performed by Dayal et al. (2014) and is confirmed to be adequate by the observational constraints performed in Gorce et al. (2017). The latter of these two papers sums up the problem nicely: ‘the escape fraction is related to too many astrophysical parameters to allow us to use a complete and fully satisfactory model’. An updated parameterisation of $\zeta$ has been established (Park et al., 2018). We defer use of this to future work where we will analyse (similarly to Section 6.4) likelihood contributions via the galaxy luminosity function. The atomic cooling threshold mentioned in Section 2 is used as the lower boundary for $T_{\rm vir}$. The upper limit is set to be consistent with the observations of Lyman break galaxies and the cooling thresholds for ionised gas. These are limited due to feedback mechanisms as galaxies above this threshold are too small to maintain star formation (Fialkov et al., 2014). This results in the prior on $T_{\rm vir}[K]$ being $[10^{4},2\times 10^{5}]$. When relaxing this upper bound to find the posterior peaks of our toy models, the exponential drop off of the halo-mass function means there are insignificant changes to the Evidence. $R_{\rm mfp}$ is motivated from a simulation stand point. The mean free path of photons should depend only on the instantaneous recombination rate of hydrogen in the IGM; hence this scale must be captured. Typically the bubbles in a simulation such as 21cmFast will be of the order $10~{}h^{-1}~{}\rm Mpc$ and so we would expect this to be roughly where this value should lie (see Figure 16(b) in Section 6.7.1 for more discussion). To allow some flexibility in the scales in which the physics can proceed this is set to $[5.,20.]~{}h^{-1}\rm~{}Mpc$ (as in 21CMMC). We note that nothing is gained by extending the prior beyond this.

To test the impact of these choices we explore the skewing of the results with widened prior ranges and with $\delta$ function priors (blue and red points respectively, Figure 15 - Section 6.7 contains the full discussion). The widened prior values are changed as follows: $\zeta\sim[0.,1000.]$; $R_{\rm mfp}\sim[0.,30.]$; $\rm Log[T_{\rm vir}]\sim[3.5,6.0]$; and $\alpha\sim[-3.0,3.0]$.

In the context of the toy models, we have expanded the priors in the plots to find the peaks of the likelihood distributions - this is to maximise their Evidence values for the sake of argument. Values obtained by integrating the full likelihood distributions are overwhelming compared to the penalties of using these large priors (large compared to FZH). The conclusions obtained from any Bayes factors are not challenged by any change in prior hyper-volume that is performed in this work. Any false Jeffreys’ scale conclusions can be disputed by physical motivations. The Evidences shown in this paper use the priors given in Table 4 - in general as long as the width of the prior captures the peak of the likelihood, Bayes factor conclusions are repeatable.

Finally we have included two reparameterisations of the MHR models (1p MHR and F MHR, both detailed in Section 3.1 with the results discussed in Section 6.1.3) as well as nested versions of the 3pFZH model (Section 6.7.1) to quantitatively asses the impact of Occam’s razor implicit in calculating the Evidence.

6 Results

Unless specified otherwise all results are tested against an FZH power spectrum with one of two fiducial sets of parameters as the simulated data, namely: f1 $[\zeta=20.,~{}R_{\rm mfp}=15.~{}\rm Mpc,~{}T_{\rm vir}=30000~{}\rm K]$; and f2 $[\zeta=15.,~{}R_{\rm mfp}=15.,~{}T_{\rm vir}=50000,~{}\alpha=0.4]$, as is done in the original 21CMMC.

In this work, the 21cmFast models is the collective name we use for FZH simulations, with 3p referring to those with $\alpha=0$ (therefore implying Equation 4) and 4p to any with $\alpha\neq 0$ (Equation 6). The toy models are FZH, Inv FZH, MHR, Inv MHR, F MHR and F Inv MHR collectively. For example, the notation 3pFZHf1 refers to a 3 parameter FZH model fitting for the f1 fiducial power spectrum.

The noise profile of HERA-331 with 1080 hours of observation is used throughout this work, unless specified otherwise. Sections 6.2, 6.5, and 6.6 respectively summarise the results that follow when simulating with LOFAR-48, multiple HERA configurations, and the SKA-512. The error on the f1 power spectrum for LOFAR, HERA-61, HERA-331 and the SKA at $z=9$ are shown in Figure 4. The SKA has longer baselines than HERA and is more sensitive at small scales. Due to HERA-331 having many more small baselines it has greater large scale sensitivity making it comparable to the SKA at model selection via the 21cm power spectrum.

6.1 Direct comparisons of 21cmNest vs 21cmMC

We first demonstrate our ability to reproduce the results of 21CMMC. These are done in comparison to two fiducial power spectra across three redshifts ($z~{}=~{}8,9,10$). As is done in Greig & Mesinger (2015), we attempt to recover the f1 and f2 fiducial parameter sets (specified in the previous section) with the 21cmFast models. Figure 5 shows the overlaid posterior distributions from 21CMMC and 21cmNest for the 3pFZHf1 and 3pFZHf2 results. A full summary of all runs shown in Table 5. The 4pFZHf1, 4pFZHf2 and 3pFZH(with $\alpha=0.4$)f2 exhibit similarly good agreement to Figure 5(a) however these are not shown.

Multinest is a more efficient algorithm for model selection purposes than Emcee is. Less likelihood evaluations are necessary due to the weighting of points calculated by the change in prior volume at each step of the nested sampling algorithm. In comparison, correlation in the MCMC chain must be removed for algorithms, such as MCEvidence (Heavens et al., 2017), that approximate the Evidence using the points in the chain - this requires an order of magnitude more points than nested sampling algorithms do. At low dimensionality (up to 4 is considered in this work) Multinest performs as well as Emcee for parameter estimation purposes also, as can be seen by the matching posteriors in Figure 5. None of the EoR models we consider consist of likelihood distributions with sharp edges that would typically cause trouble for the ellipsoidal rejection sampling used in Multinest (GAMBIT Scanner Workgroup et al., 2017). The MCMC chains produced by Emcee are run until the passing of a Gelmin-Rubin (Gelman & Rubin, 1992) convergence test. We used these settings for 21CMMC and Multinest: walkersRatio=16, burninIterations=250, sampleIterations=3000, threadCount=20, reuseBurnin=False; n_live_points=2000, max_iter=0, multimodal=True, evidence_tolerance=0.5, sampling_efficiency=0.3 respectively. On a 20 core machine, Emcee obtained 39000 likelihood evaluations in 32 hours while Multinest produced 50031 samples in 15 hours. For the problem at hand Multinest is faster, but the scaling to higher dimensions needs exploring within the EoR context.

6.1.1 Cross check: $\zeta_{\rm eff}$ the effective ionising efficiency

Figure 5 shows agreement between the two statistical algorithms well in Figure 5(a). In plot 5(b) however both exhibit a posterior peaked at a value that is different from the fiducial (shown by the pale blue lines). This is because the FZH model prescribes an effective $\zeta_{\rm eff}$ that is calculated as

[TABLE]

(note the denominator can be rewritten $f_{\rm coll}~{}\rho_{\rm M}$ from Equation 10). Since the f2 set was created with $\alpha\neq 0$ it is impossible (by definition) for the 3p model to retrieve the exact fiducial values that satisfy this mock data. This is evident in 5(b). The consistency between the $\zeta$ of 3pf1 and the $\zeta_{\rm eff}$ of 4pf1 in Table 5 shows how the added parameter allows a compensation between $\zeta$ and $\rm log_{10}[T_{\rm vir}]$. Therefore the standard deviations of $\zeta$ ($\sigma_{\zeta}$) are much larger when $\alpha$ is allowed to vary (the 4p models in Table 5). In particular note that the $\zeta_{\rm eff}$ agrees for all of the f2 simulations and all of the f1 simulations within $\sigma_{\zeta}$ (rather than for the inputted $\zeta$).

6.1.2 Cross check: recovering an Inv MHR fiducial power spectrum with the toy models

To further test our methodology, a fiducial power spectrum (mock data) has been generated by the Inv MHR model with $[\zeta=30,\rm Log[T_{\rm vir}]=4.5]$. This showed that different input models can be recovered - proof of concept for Bayesian model selection.

6.1.3 Cross Check: Do alternate parameterisations of MHR skew the Bayes factor?

Testing for the Inv MHR fiducial as in the previous section, we assess the choice of prior width with the reparameterisations of MHR. We can recover the Inv MHR mock data power spectrum using the F Inv MHR and 1p Inv MHR models detailed in Section 2.

The 1p Inv MHR is a simplified version of Inv MHR, and in particular has a vastly reduced prior space (the parameters’ volume has been reduced from $\sim 2000$ to $1$). But this should significantly reduce the predictive power of the Inv MHR model for two reasons: without the use of $f_{\rm coll}$ there is no intrinsic redshift dependence via the cosmological parameters required to calculate Equation 3; and mainly, we expect the neutral fraction to decrease with redshift as reionisation progresses.

On the other hand with F MHR by smoothing the density field via the addition of a third parameter $R_{\rm Filter}$, the filtered models can adapt to the scaling of the ionised bubbles. When fitting for f1 this complexity makes the F MHR and F Inv MHR harder to distinguish than their unfiltered counterparts. In this context however (when fitting mock data from Inv MHR), increasing the flexibility of this model increases the prior volume without any improvement to the fit (because the fiducial power spectrum was made without filtering). The 1p Inv MHR and F Inv MHR models achieve Bayes factors (odds) of $10^{-25}$ and $10^{-29}$ respectively. This illustrates that the added parameter is penalised because it is redundant in this context (for comparison 3pFZH obtains $10^{-28}$; Inv MHR has $\mathcal{B}=1$, by definition).

6.2 Can LOFAR-48 perform model selection?

Next we ask the question: what length of LOFAR observations is needed to perform Bayesian model selection with the toy EoR models?

Up until now we have used 1080 hours observing time with all the telescopes considered, as shown in Table 1. Each observation is performed with a 1 hour tracked scan of the f1 power spectrum with 6 hours of observing time per day, corresponds to 6 different fields. The Bayes factors and best fit power spectra ($z=8$) are shown for each toy model in Figure 6. Due to the large error bars on each power spectrum weak results are obtained for all except the InvFZH model. For the InvFZH model moderate results are obtained because of penalties from the observational priors (see Section 11(a) for further discussion). The breadth of power spectrum that are capable of fitting within the grey error region of the LOFAR-48 mock signal is reflected in the width of the posteriors in plot 7(a) - the posterior distributions for $R_{\rm mfp}$ and $Log_{10}[T_{\rm vir}]$ have not peaked at the fiducial values.

We now vary the observation times. LOFAR achieves161616These are a very ambitious trajectory considering LOFAR has observed a total of $\sim 1300$ hours of the NCP field in its 7 years of EoR activity. signal to noise of 1, 2, 3, 7, and 11 with 1080, 2160, 4320, 10800, and 21600 hours of integration time respectively at $z=8$ (to the nearest integer). A total of 21600 hours of observing time is necessary to obtain strong conclusions for the toy models considered in this work - Figure 7(c). The posteriors shown in Figure 7(b) now show peaks that agree with HERA-331 in Figure 5(a). In summary model selection is unlikely to be achievable with LOFAR because of the limited k-range in which it is sufficiently sensitive.

These calculations use 21cmSense’s moderate foreground settings (Section 2.2) with LOFAR for ease of comparison with the other telescopes. In principle LOFAR aims to observe within the foreground wedge (i.e. in a k range that extends lower than $k=0.15~{}\rm Mpc^{-1}$), meaning these results are conservative estimates. The results published in Patil et al. (2017) are within $k=[0.05,~{}0.13]~{}\rm Mpc^{-1}$. We repeated our 1080 hour LOFAR analyses within this k range to obtain parameter posteriors that are comparable to the 21600 hour observations in the previous k range ($[0.15,~{}1.]~{}\rm Mpc^{-1}$). In the context of model selection, the results become moderate but not strong in ruling out F Inv MHR ($\mathcal{B}=148$) and so model selection would still be challenging for LOFAR. Ideally LOFAR will be operational without any foreground contamination (the optimum 21cmSense setting) in a total possible k range of $k=[0.05,~{}1.0]~{}\rm Mpc^{-1}$. If these criteria are achievable, LOFAR can obtain $\mathcal{B}\approx 120000$ and is likely able to perform model selection. We stick with our conservative estimates to have a consistent analysis between the telescopes.

6.3 Model selection of reionisation scenarios with HERA

Here we attempt to distinguish the toy models against the f1 power spectrum (mock data generated with 3pFZH). The obtained parameter posteriors and power spectra are shown in Figures 8 and 9 respectively. As mentioned in Section 4.2 the prior ranges on these plots have been expanded to show the peaks of the likelihood distributions. Note that adding the filter scale to the MHR models smooths the density field, producing over-densities that are fewer in number but larger on all scales. The background ionising efficiency $\zeta$ must compensate by being bigger in the F Inv MHR case (plot 8(d) has a larger $\zeta$ than 8(f)) so that these over-densities can be ionised before recombining again. The opposite is true for the MHR and its filtered counterpart (plot 8(e) has a smaller $\zeta$ than 8(c)); the background ionising radiation is less efficient when ionising a smoother density field. A qualitative reflection can be interpreted from the sizes of the bubbles (signal-less regions), the smaller and more numerous bubbles require lots of faint small galaxies. By construction, the MHR models contain a degeneracy between the two main parameters ($\zeta$ and $\rm Log_{10}[T_{\rm vir}]$). This can be observed by the length of the arc in the corresponding 2d posterior plots.

The toy models predict vastly different reionisation morphologies and can all be distinguished strongly on the Jeffrey’s scale, as shown by the Bayes factors in Figure 10. The F Inv MHR proves to be the hardest to distinguish. FZH and F Inv MHR share both morphology and scale of reionisation (refer to Figure 2 for slices of the brightness temperature). In contrast, Inv MHR has the inside-out morphology, but is local rather than global in scale - this explains why it is rejected more heavily than its filtered counterpart in Figure 10. The filter scale added to the MHR models can clearly increase the flexibility of our local based models, resolving this issue. Our local inside out model proves a challenging test only when the local constraint is relaxed by smoothing the density field. The relative similarity of Inv FZH, MHR and F MHR (red points) suggest that if the morphology is different, there is little compensation achievable by the scale of reionisation.

6.4 QSO data verses the 21cm power spectrum with HERA

In this section we discuss the observational prior contribution towards the total Bayesian Evidence values calculated for the HERA-331. Since the observational priors are linearly combined into the likelihood (Figure 1) they are easily separable. By separating out the 21cm power spectrum contribution from the QSO data it becomes very apparent how little quantitative impact this observational data has on differentiating between models (the lack of deviation of the blue points from the black dashed line, Figure 11(a)). Note that these blue points have been run separately from the simulation with only the 21cm power spectrum and can be fit easily by all of our toy models (Table 3). All of the toy models posses a region in their parameter space that agrees with the constraints of the observational priors. To break this model degeneracy the 21cm power spectrum must be used. An example of this can be seen from the MAP parameters of the Inv FZH model - the values obtained differ significantly (by a few $\sigma$, see Figure 8) depending on what contributions to the likelihood are included. The MAP parameters found by the 21cm power spectrum can be heavily penalised by those found by the observational prior checks, and vice versa. The discrepancies of the black (full likelihood, Figure 11(a)) Evidence points compared with the red (only the 21cm Likelihood) show that these models are un-physical as they are penalised by these observational prior checks. Despite this the F Inv MHR model still proves the hardest to distinguish from FZH (with odds $\sim 10^{7}:1$), hence a better way of harnessing the non-21cm data is required if a cross check is to be performed. How the observational prior checks are so poorly constraining can be visualised for the McGreer and Greig neutral fraction checks (Figure 11(b)). The plotted points at the simulated redshifts are interpolated to see if they fit within the error bars on the observed neutral fractions (shown in gold). Since the error bars on these neutral fraction values are large it is difficult for them to add any constraining power to our toy models of the EoR. Similarly small likelihood contributions come from the Planck prior which utilised the optical depth. This quantitative approach agrees with the discussion in Greig & Mesinger (2017b). When recovering f1, we indirectly observe this via $\sigma$ values that are of a comparable size to that obtained parameter value (see Table 6). This should be interpreted such that other data171717This includes more numerous and precise neutral fraction measurements from quasars., as well as the 21cm power spectrum, is necessary to break model degeneracies for EoR inference.

A recent publication (Park et al., 2018) attempts to harness more information via galaxy luminosity functions, in which the parameterisation of $\zeta_{0}$ in Equation 5 is expanded via power laws in $f_{\rm esc}$ and $f_{*}$. Equation 6 is no longer used due to its inadequacy at reproducing observed UV luminosity functions. We defer a model selection analysis including this improved framework to future work.

6.5 Varying the number of HERA dipoles

In this section we illustrate what number of dipoles are required for HERA to attain the model selection capabilities of LOFAR or the SKA. As mentioned previously HERA is organised in a hexagonal structure with redundant baselines aligned for the purpose of noise reduction, for more details on this see Dillon & Parsons (2016). HERA achieves signal to noise of 3, 9, 18, 28, 39, and 50 with 19, 61, 127, 217, 331, and 469 dipoles respectively at $z=8$ on the f1 power spectrum.

As shown in figure 12 HERA-19 (the lowest dipole configuration) performs similarly to LOFAR. Model selection between the models we consider becomes strong with HERA-61. HERA-217 is comparable to SKA (see Section 6.6).

6.6 How well will the SKA-512 perform at model selection?

The analyses of the SKA-512 are shown in this section for a 1 hour tracked scan of the f1 power spectrum with 6 hours of observing time per day (split across 6 fields). We repeat the analyses that have been done for LOFAR-48 in Section 6.2. At $z=8$ SKA achieves signal to noise of 4, 7, 10, 14, and 22 with 108, 216, 324, 540, and 1080 hours of integration time.

Similarly to Section 6.2, we calculate the Bayes factors between 3pFZHf1 and the F Inv MHRf1, with decreasing the observation times. we assume the moderate foreground model for the wedge (and a buffer, Section 2.2) within 21cmSense. The toy EoR models we consider can be separated out with 324 hours of observing with the SKA. Figure 13 shows the Bayes factors obtained for the toy models considered. The F Inv MHR is used because it is the toy model that has been hardest to distinguish from 3pFZH in Figure 10.

In comparison to HERA, SKA has less repeated baselines and therefore has a marginally lower signal to noise. This is compensated for by having a wider range of baselines, meaning a fuller interferometric U-V coverage. Despite this SKA does not achieve tighter posterior contours on parameters, when considering the power spectrum - this might not remain the case when considering higher order statistics such as the bispectrum (Watkinson et al., 2019). For a direct comparison of HERA-331 and SKA’s observational error on the power spectrum see Figure 4.

In conclusion, SKA performs comparably at EoR model selection purposes as HERA-217.

6.7 Model selection on the number of parameters in 21cmFAST

Having shown that HERA & SKA can easily determine the most suitable reionisation scenario we now explore refining the number of parameters. Namely testing the ability to determine the use of a 3 or 4 parameter FZH model (detailed in Section 2). We then proceed to inspect the effect of skewing these results with a naive choice of prior.

Figure 15 shows the Bayes factors for 3p and 4p FZH models tested with $\delta$ function, used (Table 4) and widened priors in red, black and blue respectively. This figure therefore serves both aforementioned purposes.

Primarily, the black triangles show that the different parameterisations of 21cmFast considered score weakly on the Jeffrey’s scale with HERA-331. In other words the power spectrum from this 1080 hour observation is not sensitive to changes of the fiducial parameters - nor the change of a constant $\zeta$ to a power law in halo mass (Equations 4 and 6). Even when increasing the SKA’s observation times (not shown) we score weak results recovering the f1 and f2 power spectra with the 3p and 4p FZH models. Applying larger values of $\alpha$ makes these models more distinguishable however this is because the observational priors become dominant in the Evidence. Increasing $\alpha$ increases the ionising efficiency of the galaxies (Equation 5) and therefore the EoR will finish too early if $\alpha$ is largely positive.

Secondly, Figure 15 contains the information obtained when attempting to skew our own results with prior choice (see Section 5 for the full discussion) - this is the only way one can be sure the results are prior independent. Namely the blue triangles have uniform priors expanded to reduce the posterior density (and therefore the Evidence). The red shapes have the prior volume reduced to a $\delta$ function for the opposite effect. These blue results use unphysical prior ranges and still fall short of the strong Bayes factor score; we consider this a success of the machinery. In other words if one chose wide uniform prior distributions out of naivety (a sensible approach for a first attempt), the result could be moderate. The user should look closer at the physics at hand and try to reduce the wide priors (chosen for the blue points) to those resembling Table 4 - which were chosen via physical reasoning and familiarisation with the problem (used for the black points). It is worth noting that the interpolated shape of the black points is translated down for the blue points - in the uniform case the parameter prior volume is factored out of the Evidence integral as a constant (Equation 17).

For the red shapes (using $\delta$ function priors), the maximum possible Evidence is obtained by choosing the MAP parameters, $\mathcal{Z}=\mathcal{L}(\theta_{\rm MAP})$, and is observed by the alignment of the red triangles in Figure 15. The starred red point (only above 3pf1) is calculated with the fiducial parameters. It is worse than this maximum because no consideration of the observational priors was made when choosing these values. Therefore there is a disagreement between the parameters that produce the fiducial power spectrum and those satisfying the observational priors. This was also observed when cross checking the methodology against the Inv MHR fiducial data set in Section 6.1.2 (and is discussed in Section 6.4). It is worth noting here that if a prior is chosen to be too narrow, the posterior plots will show close to uniform distributions and suspicion should be applied - either the parameters dynamic range is not expressed by the chosen prior (the case here), or the parameter is redundant (Section 6.7.1).

6.7.1 The Savage-Dickey density ratio: how useful are a model’s parameters?

Here we explore redundancy of parameters within models via the Bayesian Evidence (see Section 4.1 for the methodology). We consider each parameter of the 3pf1 simulation in turn ($R_{\rm mfp}$ and $\zeta$ are shown in Figure 16). Since the $R_{\rm mfp}$ parameter is actually the maximum smoothing scale (a top-hat in k-space via the excursion-set formalism to be precise, see Section 2) it is an easy target for this evaluation. Above a certain (model dependent) threshold the maximum smoothing scale will capture all of the structural detail within the model. Below this scale we expect the model to perform poorly in model selection (against a larger $R_{\rm mfp}$) due to a loss of large-scale structure. As is shown (Figure 16(a), the power spectrum at different $R_{\rm mfp}$) by the descending left hand side of the blue power spectra while the right hand side remains fixed. This corresponds to the Bayes factor results $\mathcal{B}(R_{\rm mfp})$ from Equation 20 is plotted in Figure 16(b). All values of $R_{\rm mfp}>10$ agree with the fiducial results in green - showing negligible inference is provided by varying $R_{\rm mfp}$ in this range. For a comparison $\zeta$ produces no degenerate power spectra and is therefore a useful parameter (Figure 16(c), the power spectrum at different $\zeta$). In the Bayes factor results, $\mathcal{B}(\zeta)$ from Equation 20 is plotted in Figure 16(d), this leads to the sharp peaked value that agrees with the fiducial result shown in green. When $\zeta$ is small, the power spectra show that the simulation attempts to compensate by increasing the number of small structures (the blue power spectrum is the largest at the small scale end of Figure 16(c)). On the other hand if galaxies have a high ionising efficiency the simulation responds by decreasing the size of structures on all scales. This can be identified from the shape of each black power spectrum being maintained as the amplitude decreases. Within all of $\zeta$’s prior range, there are dynamic changes to the simulation that are observable with the power spectrum.

The lines in Plots 16(b) and 16(d) rise above the green fiducial result due to the reduction in prior volume of the parameter being nested. In 16(d) the evidence rises above the moderate Jeffreys’ scale threshold - an even wider choice of uniform prior distribution for this parameter could have provided strong evidence however this should not be taken seriously. This draws light on the implicit nature of Occam’s razor within the Bayesian Evidence, as mentioned in Section 4.2. Similarly with $R_{\rm mfp}$, we are able to see the increase in Bayes factor as the number of parameters (and therefore prior volume) is reduced. Here we obtain an weak scoring on the Jeffrey’s scale due to a smaller prior distribution ($[5,20]\rm Mpc$ compared with $\zeta$’s $[5,100]$). If in doubt, comfort is found in the variability of the power spectra per parameter as this quantitatively reveals the dynamic range of that parameter - 16(a) and 16(c). Figures 16(b) and 16(d) reflect the posterior distributions of each parameter, as derived in Section 4.1.

7 Conclusions

Initially we show that Multinest can produce matching posterior distributions to those produced by Cosmohammer in 21CMMC (Section 6.1).

The literature warns of two potential heffalump traps within Bayesian model selection: these are the prior dependency of results and the overcompensation of redundant parameters via Occam’s Razor. We address carefully these issues in Sections 4.2 and 5 to show the validity of the Bayesian model selection we apply. For a rough interpretation of our results we adopt the Jeffreys scale to distinguish weak, moderate and strong conclusions.

This work considers toy EoR models as well as those in the original 21CMMC model. Namely: an inverted version of the original 3 parameter FZH (global inside-out) model used in 21cmFAST (within which the excursion set formalism solves for reionised bubbles); a simpler model MHR (local outside-in), in which ionisation is defined via a density criteria; both of these models have had their ionisation criterion mathematical inverted (Inv FZH, global outside-in; and Inv MHR, local inside-out); and finally a 3rd parameter is added to both the MHR and Inv MHR models in the form of a Gaussian density field filter - F MHR, global outside-in; and F Inv MHR (global inside-out). By using the specifications of HERA-331 (see Table 1) for 1080 hour observations we can distinguish these toy models strongly on the Jeffreys’ scale (Section 3.1).

By calculating Evidences for the separate likelihood contributions of 21CMMC - we quantitatively show that more QSO data is required for effective non-21cm observational checks than only that of neutral fractions (surrounging high $z$ quasars) and the optical depth (Section 6.4). In future work, we will quantitatively evaluate the inference capabilities of including luminosity function data - following the methodology of Park et al. (2018).

On a subtler level the considered version of 21CMMC contains the addition of a 4th parameter, enabling a mass dependent ionising efficiency, in other words relaxing the assumption of a constant mass to light ratio for the simulated galaxies. We show in Section 6.7 that these cannot be distinguished with 1080 hour observations as they score weak results on the Jeffreys’ scale, as does changing the fiducial parameters within these models.

Using the Savage-Dickey density ratio we proceeded to show the redundancy of the $R_{\rm mfp}>10h^{-1}~{}\rm Mpc$ in 21CMMC (Section 6.7.1) - hence justifying its omission from the newest parameterisations of 21CMMC and 21cmFAST in favour of inhomogeneous recombinations (Sobacchi & Mesinger, 2014).

The most difficult to distinguish of the toy models is F Inv MHR because it shares scale and morphology with FZH. We therefore use it as a bench mark for our primary questions on feasibility:181818Assuming the ‘moderate’ foreground wedge model within 21cmSense in all cases (Section 2.2).

• •

LOFAR-48 would struggle to perform model selection. We require 21600 hours of observation to strongly disfavour the toy models (with odds $\sim 400:1$).

• •

Using 1080 observing hours, HERA would require configurations of at least 61 dipoles in order to rule out any of the toy models with strong evidence on the Jeffreys’ scale (scoring odds of $\sim 600:1$).

• •

Finally, with the SKA-512: The strong rejection of the considered EoR toy models can be achieved quickly with 324 observing hours at $\sim 500:1$.

At observing the 21cm EoR power spectrum HERA-217 becomes comparable to the SKA. As can be seen in Figure 4, SKA is better at observing small scales while HERA dominates the sensitivity at large scale. Assuming the moderate noise settings within 21cmSense can be acheived, the redundant baseline instrumental method is very promising at performing model selection with the 21cm EoR power spectrum (Byrne et al., 2018; Dillon et al., 2019).

We have now set the scene for ruling out toy models of the EoR with Bayesian Model Selection. This work shows that modest 21cm experiments, such as HERA-61, are likely to observationally pin down the correct morphology and scale of reionisation. In future work we intend to look more closely at the level of precision that can be obtained from more involved models of the EoR using the newer versions of 21CMMC. An updated version of 21cmNest191919https://binnietom@bitbucket.org/binnietom/21cmnest_1.0.git will be released soon.

Once the desired reionisation scenario can be quantitatively chosen, parameter inference will be performed rigorously to infer the involved astrophysics.

Acknowledgments

The authors would like to especially thank Catherine Watkinson for her contributions to an earlier draft of this work. We would also like to thank Daniel Mortlock, Conor O’Riordan, Brad Greig, Andrei Mesinger, Emma Chapman, Roberto Trotta, Alan Heavens, Suman Majumdar and the anonymous referee for their helpful conversations and suggestions. We acknowledge use of Linux OS, C, C++, Fortran90, gnu, gsl, Overleaf, GitHub, Python, and the Python packages numpy, scipy, cornerplot, matplotlib. T Binnie thanks STFC for their studentship funding. JRP is pleased to acknowledge support from the European Research Council under ERC grant number 638743-FIRSTDAWN.