Epoch of reionization parameter estimation with the 21-cm bispectrum
MMon. Not. R. Astron. Soc. , 000–000 (2021) Printed 5 February 2021 (MN L A TEX style file v2.2)
Epoch of reionization parameter estimation with the 21-cmbispectrum
Catherine A. Watkinson (cid:63) , Bradley Greig , , Andrei Mesinger School of Physics and Astronomy, Queen Mary University of London, Mile End Road, London E1 4NS, UK School of Physics, University of Melbourne, Parkville, VIC 3010, Australia ARC Centre of Excellence for All-Sky Astrophysics in 3 Dimensions (ASTRO 3D), University of Melbourne, VIC 3010, Australia Scuola Normale Superiore, Piazza dei Cavalieri 7, I-56126 Pisa, Italy
ABSTRACT
We present the first application of the isosceles bispectrum to MCMC parameter inferencefrom the cosmic 21-cm signal. We extend the MCMC sampler CMMC to use the fast bis-pectrum code, B I FFT , when computing the likelihood. We create mock 1000h observationswith SKA1-low, using PY O BS to account for uv-sampling and thermal noise. Assuming thespin temperature is much higher than that of the CMB, we consider two different reioniza-tion histories for our mock observations: fiducial and late-reionization. For both models wefind that bias on the inferred parameter means and 1- σ credible intervals can be substantiallyreduced by using the isosceles bispectrum (calculated for a wide range of scales and triangleshapes) together with the power spectrum (as opposed to just using one of the statistics). Wefind that making the simplifying assumption of a Gaussian likelihood with a diagonal covari-ance matrix does not notably bias parameter constraints for the three-parameter reionizationmodel and basic instrumental effects considered here. This is true even if we use extreme (un-likely) initial conditions which would be expected to amplify biases. We also find that usingthe cosmic variance error calculated with Monte-Carlo simulations using the fiducial modelparameters whilst assuming the late-reionization model for the simulated data also does notstrongly bias the inference. This implies we may be able to sparsely sample and interpolate thecosmic variance error over the parameter space, substantially reducing computational costs.All codes used in this work are publicly-available. Key words: methods: statistical – dark ages, reionization, first stars – intergalactic medium– cosmology: theory.
The Square Kilometre Array aims to detect the high-redshift 21-cm line of neutral hydrogen. It is projected to produce high pre-cision maps at a wide range of redshifts. These maps can be usedto infer the properties of early generations of stars and galaxies asthey influence the intergalactic medium (IGM) via coupling, heat-ing and ionizations (Dewdney 2016). The phase change in the Uni-verse’s ionization state induced by the latter process is called theEpoch of Reionization (EoR).Numerous studies have predicted great benefits from usinghigher-order statistics such as the bispectrum in our analysis ofsuch datasets. For example, Shimabukuro et al. 2016b; Majumdar (cid:63) Email: [email protected] The Square Kilometre Array and https://astronomers.skatelescope.org/wp-content/uploads/2016/05/SKA-TEL-SKO-0000002_03_SKA1SystemBaselineDesignV2.pdf et al. 2017; Watkinson et al. 2019b; Hutter et al. 2019 and Gorce& Pritchard 2019 show that, due to the non-Gaussian nature of thesignal, additional information is contained in higher order statis-tics, which unlike the power spectrum are sensitive to non-Gaussianstructure in a dataset. In particular, Shimabukuro et al. 2016b per-form a Fisher forecast and find that using the equilateral bispectrumin addition to the power spectrum substantially shrinks the crediblelimits of the parameters of a three-parameter EoR model comparedto those resulting from using the power spectrum alone.Furthermore, it appears that the error due to instrumentalnoise is not as large as one might naively expect; see for example,Yoshiura et al. 2015; Watkinson et al. 2019b, and Trott et al. 2019.This is because Gaussian distributed noise has zero bispectrum sothat it is only the statistical fluctuations of the noise bispectrum thatcontributes to the measured error on the bispectrum (Yoshiura et al.2015).The Fisher analysis of Shimabukuro et al. 2016b, whilst an im-portant first step towards understanding the improvements gained © 2021 RAS a r X i v : . [ a s t r o - ph . C O ] F e b C. A. Watkinson, B. Greig, A. Mesinger in performing parameter estimation with the bispectrum, likelyunderestimates the credible limits associated with each parame-ter. This is because a parameter’s covariance matrix is only accu-rately described by the inverse of the Fisher matrix if the errorson the measured quantities are perfectly Gaussian (i.e. the likeli-hood surface is Gaussian around the maximum likelihood point),which is not a given for even the 21-cm power spectrum. It hasalso been shown that the covariance predicted by a Fisher fore-cast, by the Cramer-Rao bound, provides the smallest possible at-tainable error, i.e. it provides a lower limit (Fisher 1935; Cramér1946; Rao 1945; Tegmark et al. 1997). In this paper we take thework of Shimabukuro et al. 2016b a step further by adding theisosceles bispectrum (in which we include a wide range of triangleconfigurations in addition to the equilateral) within a Monte-CarloMarkov Chain (MCMC) parameter estimation framework, buildingon the established 21-cm MCMC code CMMC (Greig & Mesinger2015a, 2017b; Park et al. 2018).Section 2 describes our bispectrum likelihood and the meth-ods used to simulate instrumental effects and measure the bispec-trum. In Section 3 we look at an idealised case with no instrumentaleffects or sample variance to see the maximal achievable improve-ment to the parameter constraints when combining the bispectrumand power spectrum. In Section 4.1 we compare analytic approxi-mations to the sample-variance error with the true sample-varianceerror calculated using Monte-Carlo (MC) methods. We will showin this section that assuming a sample-variance error that is a fixedpercentage of the statistics in any given bin is a very poor approx-imation, as is propagating the power-spectrum sample-variance er-ror onto the bispectrum assuming Gaussianity. In Section 4.2 wepresent our main analysis that include instrumental effects ( uv sam-pling & noise) and sample-variance. We will show in this sectionthat using the bispectrum in combination with the power spectrumreduces the bias (and in some cases the credible intervals) on all pa-rameters relative to that of the power-spectrum only analysis. Thisis true regardless of how likely is the realization of the "true" Uni-verse (i.e. if the initial conditions are outliers) or its reionizationhistory. CMMC
For the purposes of this analysis, we modify the latest versionof CMMC : an MC sampler of V (a python-wrapped,semi-numerical simulation of the 21-cm signal at high redshifts)(Murray et al. 2020). CMMC can be downloaded from https://github.com/21cmFAST/21CMMC , and is detailed in: Greig& Mesinger 2015a (which describes the first implementation thatused a three-parameter model for reionization), Greig & Mesinger2017b (which extends sampling to parameters responsible for heat-ing and Lyman- α coupling effects), and Park et al. 2018 (whichintroduces mass dependence to the star formation rates and es-cape fraction of ionizing radiation, as well as luminosity func-tions). The latest version of CMMC has the option of using ei-ther the EMCEE or Multinest samplers; here we use EMCEEwhich is an Affine-invariant, openMP-parallelized MCMC sam-pler (for more details see https://emcee.readthedocs.io/en/stable/ ) (Goodman & Weare 2010; Foreman-Mackeyet al. 2013). V is a standalone code for computing 3D real-izations of the 21-cm signal and its component fields. SamplingGaussian initial conditions, it uses Lagrangian perturbation the- ory to generate density and velocity fields (e.g. Bernardeau et al.2001); then using a combination of excursion set (Furlanetto et al.2004) and lightcone integration it generates ionization and temper-ature fields. We refer the interested reader to Mesinger & Furlan-etto 2007 and Mesinger et al. 2010 for details, as well as to theextensive documentation associated with the code itself availableat https://github.com/21cmfast/21cmFAST .For this demonstrative work, we use the simplest, three param-eter reionization model (as described in Greig & Mesinger 2015a),and assume the spin temperature exceeds the CMB temperature.We also compute our summary statistics from coeval cubes, insteadof lightcones. These choices keep the analysis time to a minimumfacilitating the ability to experiment with different aspects of theanalysis whilst still being informative. In future work, we will re-lax these assumptions.The parameters that we vary in our analysis are: • ζ = f esc f ∗ N γ/b (1+ n rec ) − which is the ionizing efficiencyof galaxies. Here f esc is the escape fraction of ionizing photons, N γ/b is the number of ionizing photons produced per baryon instars, and n rec is the cumulative number of IGM recombinationsper baryon. This is assumed to be a constant, and a region is deemedto be ionized if the collapsed fraction within that region is greaterthan or equal to ζ − . Increasing ζ therefore speeds up the EoR. • T vir is the minimum virial temperature needed for halos tohost star-forming galaxies (determined by cooling and feedbackmechanisms that allow star formation). Smaller T vir means star for-mation is possible in lower-mass halos that are less biased. Thusreducing T vir results in an earlier EoR, characterized by smaller,more uniformly-distributed cosmic HII regions. • R max defines the maximum distance a photon can travel in anionized IGM before it encounters a recombined atom. This effec-tive parameter can loosely be related to a characteristic mean freepath (c.f. Furlanetto & Oh 2005 and Sobacchi & Mesinger 2014).As R max is only relevant when it is smaller than the typical HIIregion size, reducing it extends the late stages of the EoR withoutimpacting the early stages.We make the assumption that the power spectrum and bis-pectrum measurements are independent (from each other and be-tween each k bin for the power spectrum or triangle configurationfor the bispectrum). We also assume independence of these statis-tics at each redshift. This allows us to approximate the total like-lihood using a simple sum over χ values. Specifically, we take ln L ( θ | d ) = − (cid:80) ij ( d ij − m ij ) / (2 σ ij ) where the indices denoteredshift and statistical bins, i.e. each ij corresponds to a the mea-surement of a single power spectrum or bispectrum bin (from thedata d ij or model m ij ) at one of the redshift bins under consid-eration. For the main results of this paper we pre-compute σ ij byforward simulating the fiducial model, each time varying the ini-tial seed of the simulation to account for sample variance error, andincluding a random realisation of instrumental noise. The standarddeviation we use in this study is calculated using 2000 such Monte-Carlo (MC) samples of the power spectrum and bispectrum in eachbin (although it is worth noting that the error estimate has mostlyconverged by 1000 iterations).We ignore the contribution to the power spectrum and bispec- A coeval cube is a datacube that has been simulated using a fixed cosmo-logical time throughout. A lightcone dataset is one in which the simulatedepoch evolves with frequency (or redshift), i.e. each slice along the z -axisrepresents a different cosmological time.© 2021 RAS, MNRAS000
For the purposes of this analysis, we modify the latest versionof CMMC : an MC sampler of V (a python-wrapped,semi-numerical simulation of the 21-cm signal at high redshifts)(Murray et al. 2020). CMMC can be downloaded from https://github.com/21cmFAST/21CMMC , and is detailed in: Greig& Mesinger 2015a (which describes the first implementation thatused a three-parameter model for reionization), Greig & Mesinger2017b (which extends sampling to parameters responsible for heat-ing and Lyman- α coupling effects), and Park et al. 2018 (whichintroduces mass dependence to the star formation rates and es-cape fraction of ionizing radiation, as well as luminosity func-tions). The latest version of CMMC has the option of using ei-ther the EMCEE or Multinest samplers; here we use EMCEEwhich is an Affine-invariant, openMP-parallelized MCMC sam-pler (for more details see https://emcee.readthedocs.io/en/stable/ ) (Goodman & Weare 2010; Foreman-Mackeyet al. 2013). V is a standalone code for computing 3D real-izations of the 21-cm signal and its component fields. SamplingGaussian initial conditions, it uses Lagrangian perturbation the- ory to generate density and velocity fields (e.g. Bernardeau et al.2001); then using a combination of excursion set (Furlanetto et al.2004) and lightcone integration it generates ionization and temper-ature fields. We refer the interested reader to Mesinger & Furlan-etto 2007 and Mesinger et al. 2010 for details, as well as to theextensive documentation associated with the code itself availableat https://github.com/21cmfast/21cmFAST .For this demonstrative work, we use the simplest, three param-eter reionization model (as described in Greig & Mesinger 2015a),and assume the spin temperature exceeds the CMB temperature.We also compute our summary statistics from coeval cubes, insteadof lightcones. These choices keep the analysis time to a minimumfacilitating the ability to experiment with different aspects of theanalysis whilst still being informative. In future work, we will re-lax these assumptions.The parameters that we vary in our analysis are: • ζ = f esc f ∗ N γ/b (1+ n rec ) − which is the ionizing efficiencyof galaxies. Here f esc is the escape fraction of ionizing photons, N γ/b is the number of ionizing photons produced per baryon instars, and n rec is the cumulative number of IGM recombinationsper baryon. This is assumed to be a constant, and a region is deemedto be ionized if the collapsed fraction within that region is greaterthan or equal to ζ − . Increasing ζ therefore speeds up the EoR. • T vir is the minimum virial temperature needed for halos tohost star-forming galaxies (determined by cooling and feedbackmechanisms that allow star formation). Smaller T vir means star for-mation is possible in lower-mass halos that are less biased. Thusreducing T vir results in an earlier EoR, characterized by smaller,more uniformly-distributed cosmic HII regions. • R max defines the maximum distance a photon can travel in anionized IGM before it encounters a recombined atom. This effec-tive parameter can loosely be related to a characteristic mean freepath (c.f. Furlanetto & Oh 2005 and Sobacchi & Mesinger 2014).As R max is only relevant when it is smaller than the typical HIIregion size, reducing it extends the late stages of the EoR withoutimpacting the early stages.We make the assumption that the power spectrum and bis-pectrum measurements are independent (from each other and be-tween each k bin for the power spectrum or triangle configurationfor the bispectrum). We also assume independence of these statis-tics at each redshift. This allows us to approximate the total like-lihood using a simple sum over χ values. Specifically, we take ln L ( θ | d ) = − (cid:80) ij ( d ij − m ij ) / (2 σ ij ) where the indices denoteredshift and statistical bins, i.e. each ij corresponds to a the mea-surement of a single power spectrum or bispectrum bin (from thedata d ij or model m ij ) at one of the redshift bins under consid-eration. For the main results of this paper we pre-compute σ ij byforward simulating the fiducial model, each time varying the ini-tial seed of the simulation to account for sample variance error, andincluding a random realisation of instrumental noise. The standarddeviation we use in this study is calculated using 2000 such Monte-Carlo (MC) samples of the power spectrum and bispectrum in eachbin (although it is worth noting that the error estimate has mostlyconverged by 1000 iterations).We ignore the contribution to the power spectrum and bispec- A coeval cube is a datacube that has been simulated using a fixed cosmo-logical time throughout. A lightcone dataset is one in which the simulatedepoch evolves with frequency (or redshift), i.e. each slice along the z -axisrepresents a different cosmological time.© 2021 RAS, MNRAS000 , 000–000 oR parameter estimation with the 21-cm bispectrum. trum for k modes that fall outside of the range . (cid:54) k (cid:54) . cMpc − . The lower k cut is motivated by avoiding modes that arelikely to suffer from corruption due to foreground leakage, and theupper cut excludes modes that will suffer from the effects of shotnoise (Greig & Mesinger 2015a). For the bispectrum this meansthat if any one of the three k -vectors that form a given triangleconfiguration fall outside of this range, then the configuration isexcluded from our likelihood calculation.We set our fiducial model parameter values as ζ = 30 . , log T vir = 4 . and R max = 15 . We also consider a late reion-ization model with ζ = 17 . , log T vir = 5 and R max = 10 .We initialise the core of CMMC to simulate coeval cubes at z =[6 . , , , , chosen to sample a range of ionized fractions, withour redshifts corresponding to x H I = [0 . , . , . , . forour fiducial model and x H I = [0 . , . , . , . for our latereionization model. Note that our late reionization model is notpicked as a realistic model, it is selected somewhat arbitrarily toprovide a test case that is quite different to the fiducial model. We use the same prior ranges as Greig & Mesinger 2015a, i.e. (cid:54) ζ (cid:54) , (cid:54) T vir (cid:54) and (cid:54) R max(bubble) (cid:54) .Our coeval cubes are and (256 Mpc) in dimension, chosento keep both sample variance and analysis time to acceptable levels(Iliev et al. 2014; Kaur et al. 2020). B I FFT - a fast code for measuring the bispectrum
The bispectrum is defined as the Fourier transform of the three-point correlation function (which measures excess probability as afunction of three points in real space). It can be written as (2 π ) B ( k , k , k ) δ D ( k + k + k ) = (cid:104) ∆( k )∆( k )∆( k ) (cid:105) , (1)where δ D ( k + k + k ) is the Dirac-delta function. According-ingly, the bispectrum is a function of three k vectors that form aclosed triangle, often referred to (as we will from here on) as a trian-gle configuration. It is necessary to perform some kind of averagingwhen measuring the bispectrum to beat down statistical noise. Asis common in bispectrum and power spectrum analysis, we chooseto perform spherical averaging, i.e. our bispectrum measurementsare functions of triangle shape and size only, not orientation.The bispectrum is the lowest order polyspectra that is sensitiveto non-Gaussian information, or structure, in a dataset. For a nicedescription of what real-space structures different k -space triangleconfigurations are sensitive to see Lewis 2011; Watkinson et al.2019b and Hutter et al. 2019 (see in particular Figure 1).Due to computational limitations, the bispectrum is often over-looked in forward-modeling frameworks. Naively, it requires mul-tiple nested loops to find the k -space pixels that form closed tri-angles of the desired shape and size. However, there are meth-ods that make the calculation tractable for many applications. Oneof these is to use Fast-Fourier Transforms to enforce the Dirac-delta function in equation 1 (Scoccimarro 2015; Sefusatti et al.2016). B I FFT is a python package that wraps a C implementation ofthe Fourier-transform bispectrum method, described in Watkinsonet al. 2017a and publicly available from https://bitbucket.org/caw11/bifft . It’s very fast, taking only a few secondsper triangle configuration on a MacBookPro (2.3GHz i9 intel core, The ionized fractions we quote are for our "standard" seed, which wediscuss in section 4.1. . This method is exten-sively described in Watkinson et al. 2017a and Watkinson et al.2021c.Throughout we will normalise out the amplitude of the bis-pectrum to isolate the non-Gaussian information: b ( k , k , k ) = B ( k , k , k ) (cid:112) ( k k k ) − P ( k ) P ( k ) P ( k ) . (2)Equation 2 is commonly applied in signal processing, see for exam-ple Hinich & Clay 1968; Kim & Powers 1978; Hinich & Messer1995 and Hinich & Wolinsky 2005. It has also been argued byBrillinger & Rosenblatt 1967 that Equation 2 is the optimal normal-isation for the bispectrum. The findings of Watkinson et al. 2019bsupport this claim for 21-cm datasets. Note also that the normalisedbispectrum is not a direct function of the power spectrum and so lin-early combining its likelihood contribution with that of the powerspectrum is not an unreasonable choice. uv sampling and noise generation with PY O BS In order to carry out our investigation we wrote PY O BS (whichcan be used as a bolt-on module for CMMC or V ) toapply uv sampling and add Gaussian random noise (with standarddeviation based on SENSE calculations) to a 21-cm brightness-temperature coeval simulation. The established code
SENSE outputs the noise andsample-variance error of the spherically-averaged power spec-trum as a function of k . PY O BS relies on an adapted version ofcalc_sense.py from SENSE which instead outputs a file con-taining the k x , k y , k z (in cMpc − ) corresponding to the instru-ment’s uv sampling and bandwidth associated with the simulationdimensions, along with the noise power spectrum associated witheach uv sample. SENSE is described extensively in Pober et al.2013a and Pober et al. 2014b. We assume optimistically that fore-grounds are fully removed and assume a track scan mode of oper-ation. On the first call to PY O BS , a maskfile of the same dimen-sion as the V simulation is created containing the noisepower in each pixel (the noise in pixels that are repeat samples arecombined coherently using inverse-covariance weighting) and ze-roed where there are no uv -samples. Once the uv -noise maskfileis written to file, PY O BS accesses it each time it is called, zeroesany unsampled pixels in the cosmological simulation and adds arandom sample of Gaussian noise to each pixel (based on the noisepower in the corresponding uv -noise maskfile pixel).By working in simulation co-ordinates (i.e. cMpc) and creat-ing the uv-noise maskfile on the first call, PY O BS is very quick,making it suitable for MC calculations, including calculating in-strumental error on any statistic (that is in itself also relativelyquick to compute). This approach is approximate in that it ignoresthe evolution of the uv sampling along the line of sight. It alsoignores the effect of the primary beam, effectively assuming thefield size is small enough to not be affected by this (which for thebox sizes simulated here is not unreasonable) or that the primarybeam as been corrected for. The SKA noise level produced by this PY O BS (using the central region from the current design for the PY O BS
21 can be used for lightcone data if it is chunked into cubes, butsince PY O BS
21 assumes a fixed redshift in translating the uv sampling ofthe instrument to simulation co-ordinates it is not the ideal tool for use withlightcones.© 2021 RAS, MNRAS , 000–000 C. A. Watkinson, B. Greig, A. Mesinger
SKA-Low phase 1 telescope model and assuming 1000 hours ob-servation time) is consistent with that predicted by Mellema et al.2013 and Koopmans et al. 2015. The SKA1-Low details and an-tenna locations used for our noise calculations are based on the lat-est SKA configuration coordinates (central region) and Dewdney2016. In this work we only consider the isosceles configuration as a func-tion of angle between k1 and k2, and for a range of scales. Ourrange of isosceles triangles span shapes from squeezed to stretched,and should therefore be able to pick up a large range of non Gaus-sian structures in the 21-cm maps. We refer the reader to Section3 of Watkinson et al. 2019b and Lewis 2011 for discussions ofthe types of structures that various configurations are sensitive to,as well as to the results of Majumdar et al. 2017 for verificationthat the isosceles configuration captures key features of reioniza-tion maps.In this section we compare the parameter constraints achievedwhen using the isosceles bipectrum (for k = k =[0 . − , . − , . − , . − ] ) and for θ/π = [0 . , . , . , . , . , . , . , . , . , . , . (where θ is the internal angle to k + k ), the power spectrum,and a combination of the two statistics. To do so we assume abest case scenario of negligible instrumental effects, perfect fore-ground removal, and negligible sample variance. In practice thisinvolves running analysis on the raw coeval cubes produced by V and assuming the same random seed for the data andmodel. We also include a modelling uncertainty of 15% of the mod-elled statistics, as is default in CMMC ).The corresponding corner plot for the three parameter modelis shown in Figure 1. Darker/lighter shading encloses 68%/95% ofthe credible limits. Different colors indicate different statistics usedfor computing the likelihood: (i) bispectrum is shown with grey; (ii)power spectrum is shown with red; and (iii) bispectrum + powerspectrum is shown with blue.Under these idealised conditions the power spectrum only(pspec-only) statistic results in tight, unbiased constraints, whichcan be seen in the bottom of Figure 1 where we plot the marginalstatistics, i.e. the marginalised posterior’s mean +/- the 68% up-per and lower credible limits. As in Greig & Mesinger 2015a, wesee a moderate degeneracy between the ionizing efficiency and theVirial temperature. This is because both parameters effect the tim-ing of reionization; for example, both a high virial temperature anda low ionizing efficiency will delay and slow the progress of reion-ization. The epoch of heating, ignored in this exploratory work,should break this degeneracy (e.g. Greig & Mesinger 2017b).We find the pspec-only statistic generally results in tighter The SKA antenna positions we use are given by the central region an-tenna positions of https://astronomers.skatelescope.org/wp-content/uploads/2016/09/SKA-TEL-SKO-0000422_02_SKA1_LowConfigurationCoordinates-1.pdf For all the statistics we consider, we disregard contribution from any k modes that fall outside of the range k f < k < k nyq where k f = 2 π/L is the fundamental k scale and L is the length of a side of the simulation,and k nyq = 1 . / . ∗ N ∗ k f where N is the resolution on a side. Aconsequence of these restrictions is that not all θ bins will be included forlarger values of k . constraints than the bispectrum only (bispec-only) statistic. Thistells us that even in the idealized scenario, the amplitude of the sig-nal is more informative than the non-Gaussian information alone(at least for our fiducial model). However, the credible intervals of R max are reduced by a factor of 0.47 relative to the pspec-only case(see also Shaw et al. 2020a). This is because R max (by applying ahard limit beyond which photons from a source will cease to beeffective at ionizing the IGM) induces structural features to whichthe bispectrum is particularly sensitive.When we combine the bispectrum with the power spectrum,the additional information from the non-Gaussianities in the mapsgreatly reduces the degeneracies of the credible limits for all the pa-rameters. This corresponds to a shrinkage of the credible intervalsby a factor of 0.70, 0.50, 0.60 for ζ , log( T vir ) , and R max respec-tively (with respect to those of pspec-only). Although we note thatthe marginalised posterior mean is closer to the truth for both ζ and log( T vir ) for the pspec-only case. This degree of improvement isroughly in agreement with Shimabukuro et al. 2016a who performa Fisher analysis using V in which they consider thesensitivity levels of LOFAR and MWA. Although, even in the bestcase scenario of a perfect observation, the degree of improvementis not as extreme as the Fisher analysis suggests. This is under-standable, since the inverse of the Fisher matrix only provides anestimate of smallest achievable credible limit. A major challenge to performing parameter estimation with 21-cmdata and simulations is correctly accounting for sample variance.Even at the level of the power spectrum this is difficult as the errordue to sample variance is dependent on the 21-cm signal itself, andtherefore the model parameters. This makes it a great challengeto model the sample-variance error using MC simulations as wehave here. One would need to effectively sample the full model pa-rameter space (which for the current most complex CMMC modelconsists of 17 astrophysical parameters, see Qin et al. 2020) at eachpoint performing at least several hundred, ideally thousands of sim-ulations with different initial conditions. This would realisticallyrequire the use of a machine-learning interpolation procedure tomake this tractable. You would also need to decide a-priori howyou are going to chop up your lightcone to measure your statisticsas a function of redshift (necessary to effectively capture the evo-lution of the signal with redshift using such summary statistics), orstore all the simulations to avoid being locked into any such choice(not a terribly practical option). It is therefore interesting to con-sider whether we might be able to approximate the sample-varianceerror using an analytic approach.We first consider whether using a constant sample varianceerror could suffice, which would allow us to use a simple fac-torial error term similar to the modelling error that is built into CMMC (intended to account for numerical inaccuracies in thesimulation code). Under this ansatz, the sample-variance error onthe power spectrum of a particular k bin, would be describedas ∆ sv P ( k ) = A P ( k ) where A is the error factor; by de-fault A = 0 . is used for the similar but distinct modelling er-ror in CMMC . We can estimate the true sample variance usinga Monte-carlo approach in which we vary the initial-condition’srandom seed and random-noise realisation assuming the fiducial © 2021 RAS, MNRAS , 000–000 oR parameter estimation with the 21-cm bispectrum.
15 20 25 30 351314151617 R m a x l o g ( T v i r ) log( T vir )
13 14 15 16 17 R max bispec-only (same seed, no inst, 15% mod err)pspec-only (same seed, no inst, 15% mod err)bispec+PS (same seed, no inst, 15% mod err) fiducial model with same seed l o g ( T v i r ) pspec-only bispec-only bispec+pspec1415 R m f p Figure 1.
Corner plot (top) for a likelihood based on the spherically-averaged isosceles bispectrum (bispec-only; grey), power spectrum (pspec-only; red), and power spectrum + bispectrum (bispec+pspec; blue). Thebottom plot shows the mean +/- 68% credible intervals for each pa-rameter. All assume a best-case scenario of no instrumental effects orforegrounds and use the same random seed for our models and data.In this and all the figures that follow, our simulations have dimen-sions of pixels and (256cMpc) and redshifts simulated are z =[6 . , , , . For our bispectrum likelihood, we use the isosceles trian-gle configurations for 11 linearly spaced θ bins and for k = k =[0 . − , . − , . − , . − ] (where θ is theinternal angle to k + k ). We see the power spectrum in such a case doesa good job of constraining the data but constraints are improved by the in-clusion of the bispectrum. model parameters. If we look at the ratio of this ‘true’ sample-variance error for the power spectrum with the power spectrummagnitude, i.e. ∆ sv P ( k ) /P ( k ) , as shown in Figure 2, we see that ∆ sv P ( k ) /P ( k ) (cid:47) . for all scales and redshifts considered.Therefore this is a reasonable choice for our fiducial model and wewill not get biased results as a result of this assumption. It is im-portant to emphasise that this may not be universally true for allmodels, as such further studies would be necessary to be able toinvestigate this point more deeply.If we perform the same exercise for the isosceles bispec-trum and plot ∆ sv B ( k , θ ) /B ( k , θ ) as we have in Figure Figure 2.
Plot of the ratio of power-spectrum sample variance (calculatedfrom repeated simulations) with the power spectrum for all the redshiftsconsidered. We see that simply assuming a fixed sample-variance error of15% (marked with the dotted horizontal line) is a reasonable choice and willnot result in any biasing of results. . B ( k , θ ) would severely underestimate the sample varianceand would undoubtedly impact on our results. It is also clear thatassuming any value for the constant sample-variance error cannotprovide a decent approximation to the true sample-variance errorcalculated using an MC approach.Assuming the signal is Gaussian, an estimate for the powerspectrum sample-variance error is given by ∆ SV ( k ) = ∆ ( k ) = k / (2 π ) P ( k ) / (cid:112) N ( k ) , where P ( k ) is the 21-cm brightness-temperature power spectrum and (cid:112) N ( k ) is the number of times aparticular mode has been sampled. Similarly, we can calculate thetheoretical bispectrum sample variance error assuming it is Gaus-sian distributed (as is often done in the case of Gaussian noise) as, [∆ sv B ( k , k , k )] = k f n V ∆ sv P ( k ) ∆ sv P ( k ) ∆ sv P ( k ) , (3)in this expression k f = 2 π/L is the fundamental k scale, V ≈ . π k k k ( s k f ) gives the number of fundamental trianglesin units of k , s k f is the binwidth, and n = 1 , , for gen-eral, isosceles and equilateral triangle configurations respectively(Scoccimarro et al. 1998b, 2004a; Liguori et al. 2010). We assume s = 1 to obtain the maximum possible estimate for the theoreticalsample-variance contribution to the bispectrum using this approxi-mation.The ratio of the MC sample-variance error to that calculatedusing Equation 3 is plotted in Figure 4 (where solid line correspondto z = 6 . , dot-dashed to z = 7 , dotted with triangles to z = 8 anddashed with circles to z = 9 ). It is clear that this approximationis orders of magnitude lower than the true sample variance for thisbox size and resolution. It is also clear there is no clean connectionbetween this theoretical sample variance and the true sample vari-ance. It is possible to improve on this theoretical approximation;for example, one can add the trispectrum contribution to the sam-ple variance of the power spectrum as per Shaw et al. 2020a, whichincludes a contribution from the non-Gaussianity of the data. Wedefer such extentions to this aspect of this study to future work.However, as is evident from the results of Section 3, Figure 4and the many works studying non-Gaussianity of the 21-cm signal,our signal is far from Gaussian. Various works including Mondal © 2021 RAS, MNRAS , 000–000 C. A. Watkinson, B. Greig, A. Mesinger
Figure 3.
Plot of the ratio of the isosceles bispectrum sample variance (cal-culated from repeated simulations) and the bispectrum for all the redshiftsand triangle configurations considered. We see that simply assuming a fixedsample-variance error of 15% would massively underestimate the samplevariance for most scales and redshifts considered. et al. 2016, Shaw et al. 2019b, and Shaw et al. 2020a have shownthat correctly accounting for this non-Gaussianity in the power-spectrum covariance will have a non-negligible effect on the re-sulting parameter constraints (provided large-scale measurementsare limited by thermal noise). Therefore, for the rest of the paperwe will use the MC estimated error estimates. s v t r u e B ( k , ) / s v t h e o r y B ( k , ) z = 6.32, k1 = k2 = 0.12z = 6.32, k1 = k2 = 0.3z = 6.32, k1 = k2 = 0.7z = 6.32, k1 = k2 = 0.98z = 7.02, k1 = k2 = 0.12z = 7.02, k1 = k2 = 0.3z = 7.02, k1 = k2 = 0.7z = 7.02, k1 = k2 = 0.98z = 8.05, k1 = k2 = 0.12z = 8.05, k1 = k2 = 0.3z = 8.05, k1 = k2 = 0.7z = 8.05, k1 = k2 = 0.98z = 9.00, k1 = k2 = 0.12z = 9.00, k1 = k2 = 0.3z = 9.00, k1 = k2 = 0.7z = 9.00, k1 = k2 = 0.98 Figure 4.
Ratio of the sample variance on the bispectrum as measured frombrute force repeat simulation to that measured from theory assuming thesignal is Gaussian. Solid line correspond to z = 6 . , dot-dashed to z = 7 ,dotted with triangles to z = 8 and dashed with circles to z = 9 . TheGaussian assumption for the sample variance is unable to even qualitativelycapture the features we see in the simulated sample variance. The initial conditions of our Universe can impact the outcome ofour parameter estimation. To quantify this, we choose a "standard"and an "extreme" model for our mock observations used for param-eter inference. Specifically, we use two different random seeds thatexhibit minimal and maximal χ from the mean of the signal, se-lected from among ∼
50 different realizations. In the analysis of thissection we use the MCMC estimated noise+sample variance error,but since we are using V for generating our mock ob-servations, we set the modelling error factor to A = 0 . . We showthe bispectrum of these two random seeds in Figure 5, we also plotin thin lines the full range of bispectrum produced in the repeatsampling we used to estimate the 1 σ sample-variance errors (whichare the error bars on each of our random seed bispectra). The plotsfrom top to bottom correspond to k = [0 . , . , . , . and z = [9 . , . , . , . . As can be seen from this plot, seed 6937 isour "extreme" seed and seed 54321 is our "standard" seed. For theinterested reader we have include the equivalent plots for the powerspectrum in Appendix A.The top plot of Figure 6 shows the resulting credible intervalswhen we use the standard seed and assume the parameters of ourfiducial model for our mock observed data. The forward model andmock observed data used for the analysis behind this plot both in-clude instrumental effects (i.e. uv sampling and noise). As before,the largest grey contour shows the bispectrum-only case, the redcontours the power-spectrum only case, and the blue contours thebispectrum + power spectrum case. For both models the true pa-rameters values (marked with the black dashed lines) lie within the95% credible intervals for all three combinations of statistic, how-ever for the fiducial model the power spectrum posterior is bimodal.Furthermore, there is more probability density in the mode that iscentred around different parameter values to the truth, leading tobiased marginal statistics (this can be seen from the marginalisedstatistics for this case which we show in the bottom plot of Fig-ure 6). The posterior for the bispectrum-only case has its probabil-ity density focused around the true parameter values for T vir and ζ ,but as with the power spectrum exhibits bias towards larger R max .All is saved by combining the power spectrum with the bispectrum, © 2021 RAS, MNRAS000
50 different realizations. In the analysis of thissection we use the MCMC estimated noise+sample variance error,but since we are using V for generating our mock ob-servations, we set the modelling error factor to A = 0 . . We showthe bispectrum of these two random seeds in Figure 5, we also plotin thin lines the full range of bispectrum produced in the repeatsampling we used to estimate the 1 σ sample-variance errors (whichare the error bars on each of our random seed bispectra). The plotsfrom top to bottom correspond to k = [0 . , . , . , . and z = [9 . , . , . , . . As can be seen from this plot, seed 6937 isour "extreme" seed and seed 54321 is our "standard" seed. For theinterested reader we have include the equivalent plots for the powerspectrum in Appendix A.The top plot of Figure 6 shows the resulting credible intervalswhen we use the standard seed and assume the parameters of ourfiducial model for our mock observed data. The forward model andmock observed data used for the analysis behind this plot both in-clude instrumental effects (i.e. uv sampling and noise). As before,the largest grey contour shows the bispectrum-only case, the redcontours the power-spectrum only case, and the blue contours thebispectrum + power spectrum case. For both models the true pa-rameters values (marked with the black dashed lines) lie within the95% credible intervals for all three combinations of statistic, how-ever for the fiducial model the power spectrum posterior is bimodal.Furthermore, there is more probability density in the mode that iscentred around different parameter values to the truth, leading tobiased marginal statistics (this can be seen from the marginalisedstatistics for this case which we show in the bottom plot of Fig-ure 6). The posterior for the bispectrum-only case has its probabil-ity density focused around the true parameter values for T vir and ζ ,but as with the power spectrum exhibits bias towards larger R max .All is saved by combining the power spectrum with the bispectrum, © 2021 RAS, MNRAS000 , 000–000 oR parameter estimation with the 21-cm bispectrum. k = 0 . , z = 9 . k = 0 . , z = 8 . k = 0 . , z = 7 . k = 0 . , z = 6 . Figure 5.
Here we plot with thin lines all 2000 bispectra used in estimat-ing the error due to sample variance for our simulation dimensions. Theplots from top to bottom correspond to k = [0 . , . , . , . and z = [9 . , . , . , . . We overplot the two random seeds used in ourparameter estimation analysis chosen from about 50 trial runs to minimise(54321) and maximise (6937) the reduced χ between them and the meanof the distribution of the thin lines shown by the thin lines in the plot. the marginal statistics of which do not suffer from bias on the in-ferred parameter values.If we now consider the results when we use the "extreme" seedfor generating our mock observed datasets, then we see that the95% credible intervals for all combinations of summary statisticstill contain the true model parameters for all parameters. However,they are in a lower probability region of the posterior than they were
20 30 401415161718 R m a x l o g ( T v i r ) log( T vir )
14 15 16 17 18 R max bispec-only (standard seed)pspec-only (standard seed)bispec + pspec (standard seed) fiducial model with instrumental noise & sampling error l o g ( T v i r ) pspec-only bispec-only bispec+pspec1618 R m f p Figure 6.
Corner plot (top) of credible intervals when using mock ob-served data generated using the fiducial model and the standard seedand the bispec-only (grey contours), the pspec-only (red contours), andbipsec+pspec (blue contours) as summary statistics in the likelihood. Theblacked dashed lines indicate the parameter values used to generate themock observed datasets for each model. The bottom plot shows the mean+/- 68% credible intervals for each parameter. All include the effects ofSKA-LOW (phase 1) uv sampling and noise, as well as sample variance,which we model the associated standard deviation using MC methods andusing the parameters of the fiducial model. Whilst all cases contain the truthwithin their 95% credible intervals, the posterior probability mass for thepspec-only case is concentrated in a different region of model parameterspace, resulting in biased marginal statistics. for the case of the more standard seed. This can be seen in Figure 7where the top plot shows the corner plot for the fiducial model withextreme seed. We see that for the case of the "extreme" seed, theweight is more evenly spread across the two posterior modes re-sulting in marginal statistics (which are summarised in the bottomplot of Figure 7) that are less biased than one might imagine fromexamining the credible intervals by eye. The bias of the marginalposterior’s mean is even reduced for the pspec-only case relativeto the results using the more standard seed for the mock observeddataset. Combining the bispectrum still improves the robustness ofthe results; however, the bias on the marginal statistics of R max is © 2021 RAS, MNRAS , 000–000 C. A. Watkinson, B. Greig, A. Mesinger
30 40 501415161718 R m a x l o g ( T v i r ) log( T vir )
14 15 16 17 18 R max bispec-only (extreme seed)pspec (extreme seed)bispec + pspec (extreme seed) fiducial model ("extreme" seed) with instrumental noise & sampling error l o g ( T v i r ) pspec-only bispec-only bispec+pspec1618 R m f p Figure 7.
Corner plot (top) of credible intervals when using mock observeddata generated using the fiducial model and the "extreme" seed for thebispec-only (grey contours), pspec-only (red contours), and bipsec+pspec(blue contours) as summary statistics in the likelihood. The blacked dashedlines indicate the parameter values used to generate the mock observeddatasets for each model. The bottom plot shows the mean +/- 68% credibleintervals for each parameter. All include SKA-LOW (phase 1) instrumentaleffects (assuming negligible primary beam effects) as well as sample vari-ance, which we model the associated standard deviation using MC methodsand using the parameters of the fiducial model. All cases contain the truthwithin their 95% credible intervals, albeit in a lower probability region ofthe posterior. not as reduced when the bispectrum and power spectrum are com-bined as it is for the more standard seed.As can be seen in the corner plot of Figure 8 (top), there ismuch less of an issue with bi-modality in the posterior for mockobserved data generated with the parameters of our late reioniza-tion model; clearly this region of parameter space is less generic(i.e. the model statistics are very distinct from those of other mod-els). The marginal statistics for this model are summarised in thebottom plot of 8. We see that for our late reionization model usingthe bispectrum in combination with the power spectrum still over-all reduces bias on the marginal statistics (although at the cost ofintroducing a small bias on the marginal statistics of ζ ) and shrinksthe credible intervals relative to those of either statistic alone. We
14 17 20 2381012141618 R m a x l o g ( T v i r ) log( T vir )
10 14 18 R max bispec-only (standard seed)pspec-only (standard seed)bispec + pspec (standard seed) Model 2 with instrumental noise & sampling error l o g ( T v i r ) pspec-only bispec-only bispec+pspec1012 R m f p Figure 8.
Corner plot (top) of credible intervals when using mock ob-served data generated using our late reionization model and the standardseed for the bispec-only (grey contours), the pspec-only (red contours),and bipsec+pspec (blue contours) as summary statistics in the likelihood.The blacked dashed lines indicate the parameter values used to generate themock observed datasets for each model. The bottom plot shows the mean +/-68% credible intervals for each parameter. All include SKA-LOW (phase1) instrumental effects (assuming negligible primary beam effects) as wellas sample variance, which we model the associated standard deviation us-ing MC methods and using the parameters of the fiducial model. Whilstall cases contain the truth within their 95% credible intervals, the poste-rior probability mass for the pspec-only case is concentrated in a differentregion of model parameter space, resulting in biased marginal statistics. found that even in test runs where we fixed the modelled initial con-ditions using the standard seed and the extreme seeds for the datathat the results for our late reionization model were still robust withno serious issue with biased results.Figure 9 shows the results for our late reionization modelwhen the extreme seed is used to generate the mock observeddata. As with the standard seed, the results for our late reioniza-tion model are less biased than they are for the fiducial model withthe 95% credible intervals of all combinations of statistic contain-ing the true model and with the combining of the bispectrum andpower spectrum improving the quality of the constraints. As wewill discuss further in the following paragraph, this is because this © 2021 RAS, MNRAS000
Corner plot (top) of credible intervals when using mock ob-served data generated using our late reionization model and the standardseed for the bispec-only (grey contours), the pspec-only (red contours),and bipsec+pspec (blue contours) as summary statistics in the likelihood.The blacked dashed lines indicate the parameter values used to generate themock observed datasets for each model. The bottom plot shows the mean +/-68% credible intervals for each parameter. All include SKA-LOW (phase1) instrumental effects (assuming negligible primary beam effects) as wellas sample variance, which we model the associated standard deviation us-ing MC methods and using the parameters of the fiducial model. Whilstall cases contain the truth within their 95% credible intervals, the poste-rior probability mass for the pspec-only case is concentrated in a differentregion of model parameter space, resulting in biased marginal statistics. found that even in test runs where we fixed the modelled initial con-ditions using the standard seed and the extreme seeds for the datathat the results for our late reionization model were still robust withno serious issue with biased results.Figure 9 shows the results for our late reionization modelwhen the extreme seed is used to generate the mock observeddata. As with the standard seed, the results for our late reioniza-tion model are less biased than they are for the fiducial model withthe 95% credible intervals of all combinations of statistic contain-ing the true model and with the combining of the bispectrum andpower spectrum improving the quality of the constraints. As wewill discuss further in the following paragraph, this is because this © 2021 RAS, MNRAS000 , 000–000 oR parameter estimation with the 21-cm bispectrum.
12 14 16 18 20 22101520 R m a x l o g ( T v i r ) log( T vir )
10 15 20 R max bispec-only (extreme seed)pspec-only (extreme seed)bispec + pspec (extreme seed) Model 2 ("extreme" seed) with instrumental noise & sampling error l o g ( T v i r ) pspec-only bispec-only bispec+pspec1015 R m f p Figure 9.
Corner plot (top) of credible intervals when using mock ob-served data generated using our late reionization model and the "extreme"seed for the bispec-only (grey contours), pspec-only (red contours), andbipsec+pspec (blue contours) as summary statistics in the likelihood. Theblacked dashed lines indicate the parameter values used to generate themock observed datasets for each model. The bottom plot shows the mean +/-68% credible intervals for each parameter. All include SKA-LOW (phase 1)instrumental effects (assuming negligible primary beam effects) as well assample variance, which we model the associated standard deviation usingMC methods and using the parameters of the fiducial model. Here the com-bining the bispectrum with the power spectrum still helps less in relievingbias in the credible intervals, although the marginal statistics are seen toreturn reasonable predictions of the true parameters. model is at a much early stage of the reionization process for whichdifferences between seeds are suppressed relative to that of the fidu-cial model.What is potentially important about the results of the our latereionization model analysis, is that we have used the standard-deviation due to sample variance as calculated for the fiducialmodel, rather than calculating it for the our late reionization modelparameters, i.e. we have seen no serious negative impact from as-suming sample variance is the same in both regions of model pa-rameter space, despite them being very different models. This islikely because the sample-variance error for the later-reionizationmodel is smaller or similar to that of the fiducial model because as
Figure 10.
Covariance matrix for the power spectrum for all bins and red-shifts considered here.We see that there is correlations between statisticalbins in different redshifts, most notably between z = 7 and z = 8 . the process of reionization is less advanced (the late-reionizaton’sionized fraction is only 0.7 at the lowest redshift we consider as op-posed to 0.3 in the fiducial model). In the later stages of reionization(in the regime of sparse neutral islands) the amplitude varies morebetween realisations, as can be seen by the trend further away fromthe theoretical sample-variance with decreasing redshift in Figure3. This result implies that one could use the sample-variance froma single well-chosen model for all regions in parameter space. Abetter, and still tractable, option would be to sparsely sample thesample-variance error in parameter space and use some form ofinterpolation to approximate the sample-variance error in other re-gions of parameter space. Whether or not this would be a sufficientapproximation, and whether this finding extends to the full covari-ance matrix should be addressed in future work.It is clear that using a diagonal covariance matrix and assum-ing independence between statistical bins are not disastrous as-sumptions in that the true parameters are constrained by the re-sulting parameter estimation analysis, even if we consider outlierdata. However, it will give stronger and more robust results tonot make such assumptions and to use a fully multivariate Gaus-sian likelihood that includes all correlations between the statisticalbins, statistics and redshifts. We have discussed the difficulty of ac-curately accounting for sample variance errors, it will equally bechallenging to capture correlations between redshifts, which canbe seen in Figure 10 where we plot the covariance matrix for thepower spectrum. These correlations would likely be less severe ifwe were working with chunks of lightcones, which is the more cor-rect thing to consider; however, it is unlikely that there would beno correlations whatsoever. It is also likely, that as the complexityof our forward model increases (necessary if we are to fully char-acterise the instrumental effects, foreground residuals, ionosphericeffects, unresolved RFI and polarisation leakage) that assuming amultivariate Gaussian form for the likelihood will be insufficient.A method to bypass all these issues would be to uselikelihood-free inference, which bypasses the need to ever pre- We do not plot the correlations between the power spectrum and bis-pectrum, because the amplitude contribution has been normalised out ofour bispectrum analysis; we therefore expect correlations between the twostatistics to be negligible.© 2021 RAS, MNRAS , 000–000 C. A. Watkinson, B. Greig, A. Mesinger calculate a covariance matrix since the likelihood (or posterior de-pending on the type of likelihood-free inference) is estimated us-ing forward simulations during the inference process. It also meansone need never explicitly write down a likelihood. This approachwill also be able to deal with cross-correlations of the cosmologicalsignal with the noise and foregrounds biasing parameter-inferenceresults, as seen in Nasirudin et al. 2020 who perform far more accu-rate and detailed forward-modelling than that attempted here (theyalso use a fully multivariate Gaussian likelihood). We will discussthe application of likelihood-free methods as applied to 21-cm ob-servations in Watkinson, Alsing, Greig & Mesinger (in prep).
In this work, we have added an isosceles bispectrum likelihoodmodule to the established CMMC code that assumes independencebetween all statistical bins and redshifts. We are able to make thisassumption by using a normalised version of the bispectrum inwhich the contribution of the power spectrum to the bispectrum hasbeen removed. In order to perform our analysis we use two newpublicly available codes B I FFT (to measure the bispectrum withsufficient speed) and PY O BS (to simulate uv sampling and ran-dom instrumental noise for coeval cubes).We generate various mock observations by varying astrophys-ical parameters as well as the random seed for initial conditions.We consider two sets of astrophysical parameters, that result in dif-ferent reionization histories: a fiducial model and a late reionizationmodel. We also consider two different random seeds: one chosen toproduce relatively standard bispectra (in terms of its χ comparedto the mean) and another to produce more extreme outlier bispectradata.Various approaches for handling the bispectrum sample-variance error term have also been considered. We find that thebispectrum sample-variance error cannot be effectively describedas a percentage of the bispectrum in a given bin, nor by propagat-ing the power spectrum sample-variance error onto the bispectrumunder the assumption of Gaussianity. We find that using the 1 σ error generated using Monte-Carlo methods for a simple 1D Gaus-sian likelihood is sufficient to constrain the parameters of the threeparameter model of reionization considered here. We also find thatusing the sample-variance error generated under our fiducial modelwhilst assuming simulated data from a late-reionization model hasno serious negative impact on our results. This is important as itimplies that we may be able to get away with a sparse samplingof the bispectrum sample-variance error as a function of parame-ter space combined with some form of interpolation to estimate theerror term at the unsampled points of parameter space.We find that combining the power spectra and the bispectrumin the likelihood can significantly reduce the bias away from theinput reionization parameters, for all of the mock observations andmodels considered here (see also Gazagnes et al. 2020). For thelate-reionization model we also see a reduction in the credible lim-its. These findings hold true even if we consider outlier mock ob-servations.Further work is needed to establish the improvements from us-ing the bispectrum in more complex models for reionization, suchas the mass-dependent parametrisation including spin temperaturefluctuations V model. It will also be important for fu-ture works that consider the issue of modelling the bispectrum sam-ple variance, to better understand its dependence on simulation res-olution and dimensions. It will also be necessary to get a better un- derstanding of how these results will be impacted by the inclusionof more levels of observational realism, as there has already beenindication that foreground residuals and observational effects willbe more problematic for the bispectrum (Watkinson et al. 2021c). ACKNOWLEDGEMENTS
CAW’s research is currently supported by a UK Research and In-novation Future Leaders Fellowship, grant number MR/S016066/1(PI Alkistis Pourtsidou). CAW also thanks Jonathan Pritchard forfinancial support during the early stages of this project via theEuropean Research Council under ERC grant number 638743-FIRSTDAWN, as well the ARC Centre for Excellence in All-SkyAstrophysics in 3D Visitor Fund. This research utilised QueenMary’s Apocrita HPC facility, supported by QMUL ITS Research.http://doi.org/10.5281/zenodo.438045. Parts of this research weresupported by the Australian Research Council Centre of Excel-lence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D),through project number CE170100013. AM acknowledges supportfrom the European Research Council (ERC) under the EuropeanUnion’s Horizon 2020 research and innovation programmes (AIDA–
The codes used to produce this work are all publicly available withlinks to access provided in the main text.
REFERENCES
Bernardeau F., Colombi S., Gaztanaga E., Scoccimarro R.,2001, Phys. Rep., 367, 1, arXiv:0112551, doi:10.1016/S0370-1573(02)00135-7Brillinger D., Rosenblatt M., 1967, in Bernard Harris ed., , Spectr.Anal. time Ser.. Wiley, New York, pp 189–232Cramér H., 1946, Mathematical Methods of Statistics (PMS-9).Princeton University Press, doi:10.1515/9781400883868Dewdney P., 2016, Technical report, SKA1 SYSTEM BASELINEDESIGN V2. SKA OfficeFisher R. A., 1935, J. R. Stat. Soc., 98, 39, doi:10.2307/2342435Foreman-Mackey D., Hogg D. W., Lang D., Goodman J.,2013, Publ. Astron. Soc. Pacific, 125, 306, arXiv:1202.3665,doi:10.1086/670067Furlanetto S. R., Oh S. P., 2005, MNRAS, 363, 1031,arXiv:0505065, doi:10.1111/j.1365-2966.2005.09505.xFurlanetto S. R., Zaldarriaga M., Hernquist L., 2004, Astrophys.J., 613, 1, arXiv:0403697, doi:10.1086/423025Gazagnes S., Koopmans L. V., Wilkinson M. H., 2020, arXiv,arXiv:2011.08260Goodman J., Weare J., 2010, Commun. Appl. Math. Comput. Sci.,5, 65, doi:10.2140/camcos.2010.5.65Gorce A., Pritchard J. R., 2019, arXiv:1903.11402Greig B., Mesinger A., 2015a, MNRAS, 449, 4246,doi:10.1093/mnras/stv571Greig B., Mesinger A., 2017b, MNRAS, 472, 2651,doi:10.1093/mnras/stx2118Hinich M. J., Clay C. S., 1968, Rev. Geophys., 6, 347,doi:10.1029/RG006i003p00347 © 2021 RAS, MNRAS , 000–000 oR parameter estimation with the 21-cm bispectrum. Hinich M., Messer H., 1995, IEEE Trans. Signal Process., 43,2130, doi:10.1109/78.414775Hinich M. J., Wolinsky M., 2005, J. Stat. Plan. Inference, 130,405, doi:10.1016/J.JSPI.2003.12.022Hutter A., Watkinson C. A., Seiler J., Dayal P., Sinha M., CrotonD. J., 2019, MNRAS, doi:10.1093/mnras/stz3139Iliev I. T., Mellema G., Ahn K., Shapiro P. R., Mao Y., Pen U.-L.,2014, MNRAS, p. 725, arXiv:1310.7463Kaur H. D., Gillet N., Mesinger A., , 2020, Minimum size of cos-mological 21-cm simulationsKim Y. C., Powers E. J., 1978, Phys. Fluids, 21, 1452,doi:10.1063/1.862365Koopmans L. et al., 2015, Proc. Adv. Astrophys. with Sq. Km.Array, AASKA14Lewis A., 2011, J. Cosmol. Astropart. Phys., 10, 1475,arXiv:1107.5431, doi:10.1088/1475-7516/2011/10/026Liguori M., Sefusatti E., Fergusson J. R., Shellard E.P. S., 2010, Adv. Astron., 2010, 64, arXiv:1001.4707,doi:10.1155/2010/980523Majumdar S., Pritchard J. R., Mondal R., Watkinson C. A.,Bharadwaj S., Mellema G., 2017, MNRAS, 476, 4007,arXiv:1708.08458, doi:10.1093/mnras/sty535Mellema G. et al., 2013, Exp. Astron., 36, 235, arXiv:1210.0197Mesinger A., Furlanetto S. R., 2007, Astrophys. J., 669, 663,doi:10.1086/521806Mesinger A., Furlanetto S., Cen R., 2010, MNRAS, 411, 955,arXiv:1003.3878, doi:10.1111/j.1365-2966.2010.17731.xMondal R., Bharadwaj S., Majumdar S., 2016, Mon. Not.R. Astron. Soc. Vol. 464, Issue 3, p.2992-3004, 464, 2992,arXiv:1606.03874, doi:10.1093/mnras/stw2599Murray S., Greig B., Mesinger A., Muñoz J., Qin Y., ParkJ., Watkinson C., 2020, J. Open Source Softw., 5, 2582,arXiv:2010.15121, doi:10.21105/joss.02582Nasirudin A., Murray S. G., Trott C. M., Greig B., Joseph R. C.,Power C., 2020, Astrophys. J., 893, 118, arXiv:2003.08552,doi:10.3847/1538-4357/ab8003Park J., Mesinger A., Greig B., Gillet N., 2018, Mon. Not.R. Astron. Soc. Vol. 484, Issue 1, p.933-949, 484, 933,arXiv:1809.08995, doi:10.1093/mnras/stz032Pober J. C. et al., 2013a, Astron. J., 145, 65, doi:10.1088/0004-6256/145/3/65Pober J. C. et al., 2014b, Astron. J., 782, 66, arXiv:1310.7031Qin Y., Mesinger A., Park J., Greig B., Muñoz J. B., 2020, MN-RAS, 495, 123, arXiv:2003.04442, doi:10.1093/mnras/staa1131Rao C., 1945, Bull. Calcutta Math. Soc., 37, 81Scoccimarro R., 2015, Phys. Rev. D, Vol. 92, Issue 8, id.083532,92, arXiv:1506.02729, doi:10.1103/PhysRevD.92.083532Scoccimarro R., Sefusatti E., Zaldarriaga M., 2004a, Phys. Rev.D, 69, 1550, arXiv:0312286, doi:10.1103/PhysRevD.69.103513Scoccimarro R., Colombi S., Fry J. N., Frieman J. A., Hivon E.,Melott A., 1998b, Astrophys. J., 496, 586, doi:10.1086/305399Sefusatti E., Crocce M., Scoccimarro R., CouchmanH., 2016, MNRAS, 460, 3624, arXiv:1512.07295,doi:10.1093/mnras/stw1229Shaw A. K., Bharadwaj S., Mondal R., 2019b, arXiv:1902.08706,doi:10.1093/mnras/stz1561Shaw A. K., Bharadwaj S., Mondal R., 2020a, MNRAS,arXiv:2005.06535, doi:10.1093/mnras/staa2090Shimabukuro H., Yoshiura S., Takahashi K., Yokoyama S., IchikiK., 2016a, p. 11, arXiv:1608.00372Shimabukuro H., Yoshiura S., Takahashi K., Yokoyama S.,Ichiki K., 2016b, MNRAS, 468, 1542, arXiv:1608.00372, doi:10.1093/mnras/stx530Sobacchi E., Mesinger A., 2014, MNRAS, pp 1662–1673,arXiv:1402.2298Tegmark M., Taylor A., Heavens A., 1997, Astrophys. Jour-nal, Vol. 480, Issue 1, pp. 22-35., 480, 22, arXiv:9603021,doi:10.1086/303939Trott C. M. et al., 2019, Publ. Astron. Soc. Aust., 36, e023,arXiv:1905.07161, doi:10.1017/pasa.2019.15Watkinson C. A., Majumdar S., Pritchard J. R., 2017a, MNRAS,472, 2436, arXiv:1705.06284, doi:10.1093/mnras/stx2130Watkinson C. A. et al., 2019b, MNRAS, 482, 2653,arXiv:1808.02372, doi:10.1093/mnras/sty2740Watkinson C. A., Trott C. M., Hothi I., 2021c, MNRAS, 501, 367,doi:10.1093/mnras/staa3677Yoshiura S., Shimabukuro H., Takahashi K., MomoseR., Nakanishi H., Imai H., 2015, MNRAS, 451, 266,doi:10.1093/mnras/stv855
APPENDIX A: THE SAMPLE VARIANCE OF THEPOWER SPECTRUM
In Figure A1 we plot the 2000 power-spectra realisations that weuse to calculate our sample-variance errors. We also overplot thetwo random seeds used for mock observed data in this study. Aswith the bispectrum the extreme seed (purple solid like) is morethan σ away from the more standard seed (blue dot-dashed line)for many bins, especially at the later stages of reionization, i.e. z (cid:54) .This paper has been typeset from a TEX/ L A TEX file prepared by theauthor. © 2021 RAS, MNRAS , 000–000 C. A. Watkinson, B. Greig, A. Mesinger k /Mpc468101214 P ( k ) z = 6.32167913284 z = 6 . k /Mpc101520253035 P ( k ) z = 7.02489125294 z = 7 . k /Mpc10152025303540 P ( k ) z = 8.04717039344 z = 8 . k /Mpc51015202530 P ( k ) z = 9.00285740683 z = 9 . Figure A1.
Here we plot with thin lines all 2000 power spectra used inestimating the error due to sample variance for our simulation dimensions.The plots from top to bottom correspond to z = [6 . , . , . , . . Weoverplot the two random seeds used in our parameter estimation analysischosen from about 50 trial runs to minimise (54321) and maximise (6937)the reduced χ between them and the mean of the distribution of the thinlines shown by the thin lines in the plot. © 2021 RAS, MNRAS000