Fundamental uncertainty levels of 21cm power spectra from a delay analysis

TL;DR

This paper quantifies the sample variance in 21cm power spectrum measurements from EoR experiments, highlighting its dependence on array design and implications for data analysis.

Contribution

It provides the first detailed analysis of sample variance levels in delay-transform based 21cm power spectrum estimators for redundant arrays.

Findings

Sample variance can reach 20-30% for individual baselines.

Combining multiple baselines reduces variance to 1-3%.

Results inform observing strategies and parameter constraints.

Abstract

Several experimental efforts are underway to measure the power spectrum of 21cm fluctuations from the Epoch of Reionization (EoR) using low-frequency radio interferometers. Experiments like the Hydrogen Epoch of Reionization Array (HERA) and Murchison Widefield Array Phase II (MWA) feature highly-redundant antenna layouts, building sensitivity through redundant measurements of the same angular Fourier modes, at the expense of diminished UV coverage. This strategy limits the numbers of independent samples of each power spectrum mode, thereby increasing the effect of sample variance on the final power spectrum uncertainty. To better quantify this effect, we measure the sample variance of a delay-transform based power spectrum estimator, using both analytic calculations and simulations of flat-spectrum EoR-like signals. We find that for the shortest baselines in HERA, the sample variance can reach as high as 20%, and up to 30% for the wider fields-of-view of the MWA. Combining estimates from all the baselines in a HERA- or MWA-like 37 element redundant hexagonal array can lower the variance to $1-3$% for some Fourier modes. These results have important implications for observing and analysis strategies, and suggest that sample variance can be non-negligible when constraining EoR model parameters from upcoming 21cm data.

Tables

Table 1. Table 1: Summary of simulation parameters.

${\sigma }_0$	31mK
Bandpass	100 – 130 MHz
Channel Width	78.125 kHz
Duration	24 Hours
Integration Time	11 seconds
Latitude	-30.72∘
Longitude	21.428∘
MWA beam FWHM	31.5∘
HERA beam FWHM	12.13∘
Field of View	110∘

Table 2. Table 2: Results at 8 hours and 24 hours for MWA and HERA beams, for single baselines, and for the spherically-binned power spectra of the 37 element hex. The binning used is discussed in section   7.3 . The maximum uncertainty across all k 𝑘 k bins is shown for the 37hex binned case, so no uncertainty is quoted.

	MWA		HERA
	8 hrs	24 hrs	8 hrs	24 hrs
14.6m	43.86% $\pm$ 2.94%	26.22% $\pm$ 2.13%	34.01% $\pm$ 2.15%	19.89% $\pm$ 1.50%
140m	16.55% $\pm$ 0.92%	9.57% $\pm$ 0.69%	21.23% $\pm$ 1.12%	12.37% $\pm$ 0.92%
37hex	3.41%	2.04%	3.22%	1.86%

Equations

$$\tilde{V}(\mathbf{b},\tau ,t)$$

$$=\int \int A(\hat{s})T(\hat{s},\nu ,t)e^{-2\pi i(\mathbf{b}\cdot \hat{s}\nu /c-\nu \tau )}d\Omega d\nu$$

$$\hat{P}(k_{b,\tau },t)$$

$$k_{b,\tau }$$

$${\hat{P}}_{\text{avg}}(k_{b,\tau };N_t)=\frac{1}{N_t}\sum _{i=0}^{N_t}\hat{P}({\mathbf{k}}_{\tau },t_i)$$

$$\hat{P}(k_{b,\tau },N_t;b)=\frac{\Phi }{N_t}\sum _{n=1}^{N_t}|\tilde{V}(b,\tau ,t_n){|}^2$$

$${\hat{P}}_N\equiv \hat{P}(k_{b,\tau },N_t;b)={\mathbf{v}}_{\tau }^{†}\mathrm{A}{\mathbf{v}}_{\tau }$$

$$\text{Var}[{\mathbf{v}}^{†}\mathrm{A}\mathbf{v}]=\text{Tr}[(\mathrm{AC}{)}^{†}(\mathrm{AC})]+2{\mathbf{\mu }}^{†}\mathrm{ACA}\mathbf{\mu }$$

$$\text{Var}[{\hat{P}}_{N_t}]={\left(\frac{\Phi }{N}\right)}^2\sum _{i,j}^{N_t}|{\mathrm{C}}_{ij}{|}^2$$

$${\mathrm{C}}_{ij}=\frac{P(k_{b,\tau })}{X^{2}Y}{\left(\frac{2k_B}{{\lambda }^2}\right)}^2$$

$$\times \exp \left[2\pi i\frac{\nu }{c}(\hat{s}-{\mathrm{R}}_{ij}\hat{s})\cdot \mathbf{b}\right]d\Omega d\nu$$

$$C_{ij}=\int d\nu {\Theta }_{\nu }$$

$${\Theta }_{\nu }$$

$$=\int d\Omega A_{\nu }^{\ast }(\hat{s})A({\mathrm{R}}_{ij}\hat{s})\exp \left[-2\pi i\frac{\nu }{c}\hat{s}\cdot {\mathbf{b}}^{\prime }\right]$$

$$\text{where }{\mathbf{b}}^{\prime }=(\mathcal{I}-{\mathrm{R}}_{ij}^T)\mathbf{b}$$

$${\Theta }_{\nu }=\int A_{\nu }^{\ast }(\vec{l})A_{\nu }({\mathrm{\Gamma }}_{ij}\vec{l})\exp \left[-2\pi i\frac{\nu }{c}\vec{l}\cdot {\vec{b}}^{\prime }\right]$$

$${\Theta }_{\nu }\simeq \int {\tilde{A}}^{\ast }\left(\frac{\nu }{c}{\vec{b}}^{\prime }-\vec{u}\right)\tilde{A}({\mathrm{\Gamma }}_{ij}^T\vec{u})d^{2}u$$

$${\gamma }_n={\Phi }^2|C_{ij}{|}^2\text{ such that }n=|i-j|$$

$$\text{Var}[P_k]=\frac{1}{N}\left(2\sum _{n=0}^{N-1}{\gamma }_n-{\gamma }_0\right)-\frac{2}{N^2}\sum _{n=0}^{N-1}n{\gamma }_n$$

$$\frac{1}{N}\sum _{n=0}^{N-1}{\gamma }_n$$

$$\frac{1}{N^2}\sum _{n=0}^{N-1}n{\gamma }_n$$

$$\gamma (t)={\gamma }_0\exp \left(\frac{-t^2}{2w_c^2}\right)$$

$$\text{Var}[{\hat{P}}_k(T)]={\gamma }_0^2$$

$$-2{\left(\frac{w_c}{T}\right)}^2\left[1-e^{-T^2/2w_c^2}\right])$$

$$q(l)=\frac{1}{l}\sum _{j=1}^l|x_j{|}^2\text{ for }l\le N_t$$

$$f(q;k,\alpha )=\frac{(\alpha q{)}^{k/2-1}\exp (-(\alpha q)/2)}{2^{k/2}\Gamma (k/2)}$$

$$\frac{l}{k}\to l_c\text{ as }l\to N_t$$

$${\hat{t}}_c=k\Delta t$$

$$V_b(\nu ,t_i)$$

$$=\sum _{n\in \text{FoV}}\omega T({\hat{s}}_n,t_i,\nu )A({\hat{s}}_n)e^{-2\pi i{\hat{s}}_n\cdot \mathbf{b}\nu /c}$$

$${\sigma }^2(r)=\frac{{\sigma }_0^{2}d{V}_0}{dV(r)}$$

$$s_k(t)=\frac{\sqrt{⟨|\hat{P}(k;N_t)-P_{\text{theory}}{|}^2⟩}}{P_{\text{theory}}}$$

$$t_c\approx t/M=t{\left[s_k(t)\right]}^2$$

$$\Delta t_c\approx 2t{s}_k(t)\Delta s_k(t)$$

Full text

††pubyear: 2019††pagerange: Fundamental uncertainty levels of 21cm power spectra from a delay analysis–Fundamental uncertainty levels of 21cm power spectra from a delay analysis

1 Introduction

The Epoch of Reionization (EoR) was the last global phase transition of the universe, when the first UV emitting sources formed and carved out bubbles of ionization in the previously-neutral intergalactic medium. The lack of bright sources and abundance of absorbing gas makes the EoR, and the dark ages that preceded it, inaccessible to traditional high redshift survey techniques. The hyperfine “spin-flip” transition of neutral hydrogen, which emits a photon of rest-wavelength 21cm, offers one of the most promising observational probes of the EoR (Furlanetto et al., 2006; Morales, 2010). The brightness of 21cm emission/absorption is tightly bound to the ionization state, density, and temperature of the gas, and as a line transition at cosmological distances, the redshift directly maps to distance. This technique therefore has the potential to directly map the three-dimensional structure of hydrogen during the EoR. The low opacity of hydrogen to 21cm means that this signal is both extremely weak and unobstructed by intervening gas, with brightness contrasts expected to be of order $\mathrm{mK}$ (Pritchard & Loeb, 2008)

The wealth of information promised by 21cm observations has spurred a renewed interest in low-frequency radio astronomy over the last decade. However, the expected weakness of the signal presents a considerable challenge. Directly imaging structures of the EoR will likely require an instrument of the scale of the upcoming Square Kilometre Array (SKA; Mellema et al. 2013; Furlanetto et al. 2006), so the current generation instruments are seeking statistical measures of the EoR signal, especially the power spectrum. Experiments like the GMRT (Paciga et al., 2013, 2011), MWA (Beardsley et al., 2016; Bowman et al., 2013), PAPER (Parsons et al., 2010; Ali et al., 2015; Cheng et al., 2018), and LOFAR (Patil et al., 2017), have placed upper limits on the power spectrum amplitude at several redshifts.

Choosing the optimal interferometer design to maximize power spectrum sensitivity requires balancing the effects of antenna positions, receiver element, and computational requirements of the design. The most significant factors that need to be reduced are thermal noise and foreground power (Franzen et al., 2016; Yatawatta et al., 2013; Bernardi et al., 2010; Jelić et al., 2008), so much work has focused on these effects. Short baselines are more sensitive to diffuse structures like the EoR signal, and a random antenna distribution typically improves imaging capability, which helps foreground modeling and subtraction (Lidz et al., 2008). This motivated the Phase I layout of the Murchison Widefield Array (MWA), in Murchison Radio Observatory in Western Australia, which consisted of a random distribution of 128 elements, featuring a dense core but with several baselines reaching as far as a few kilometers. On the other hand, a highly-redundant array can be made more sensitive to line of sight Fourier modes with relatively few antennas, because redundant baselines directly sample the same angular Fourier modes ($k_{\perp}$) (Dillon & Parsons, 2016; Parsons et al., 2012a). This has the additional benefit of improving calibration through redundant techniques (Liu et al., 2010).

With these considerations, many current EoR experiments feature highly redundant layouts, fine spectral resolution, and drift-scanning, zenith-pointing antennas. The Hydrogen Epoch of Reionization Array (HERA), under construction in the Karoo desert of South Africa, comprises a densely packed hexagonal grid of 14.6$\mathrm{m}$ parabolic dishes with dual-polarization feeds at the prime focus (DeBoer et al., 2017). It is sensitive to a bandpass of 100 – 200 $\mathrm{MHz}$ with 1024 frequency channels (97 $\mathrm{kHz}$ channel width). When finished, it will have 350 dishes, 30 of which are outriggers placed farther away from the compact core to improve imaging capabilities.

The Phase II deployment of the MWA is a hybrid design (Wayth et al., 2018) with 54 of the original core tiles of Phase I left in place, and the remaining 74 tiles re-purposed into two 37 element hexagonal grids with a 14 $\mathrm{m}$ spacing. Each “tile” is a phased array of 16 dual-polarization dipole antennas acting as a single element in the correlator. This hybrid design allows for a combination of calibration and foreground removal techniques to be applied, and provides data with well-tested hardware for comparison with HERA.

The main trade-off to these redundant designs is that for more antennas there are fewer independent samples of each Fourier mode, due to the redundancy of many measurements. Statistical quantities such as the power spectrum have an intrinsic uncertainty, called sample or cosmic variance, due to the limited size of the volume surveyed by an instrument. For a Gaussian signal, this is proportional to the power spectrum amplitude over the number of independent samples $N$ (see McQuinn et al. 2006). This holds true for non-Gaussian signals when $N$ is sufficiently large, due to the central limit theorem. The root-mean-square (RMS) error is the square root of the sample variance, and drops off as $1/\sqrt{N}$. Assuming no correlation between neighboring times, the number of independent samples of each $k_{\parallel}$ mode is just the number of integration times, but real observations do have correlations in time. Each radio interferometric measurement (visibility) is an integral over the primary beam across the sky, and visibilities at neighboring times will integrate over overlapping regions. The effective number of independent samples in a given survey volume is decreased by this correlation, so the sample variance decreases more slowly.

Previous sensitivity estimates have argued that sample variance is negligible except in a handful of k-modes, and so subsequent analysis has mostly ignored the issue in light of the more significant issues of thermal noise and foreground contamination. Since foreground sources tend to be smooth spectrum, it is possible to isolate foreground power to low-$k_{\parallel}$ modes, since $k_{\parallel}$ is proportional to the Fourier dual of frequency. This power tends to spread to higher $k_{\parallel}$ modes for longer baselines, forming the so-called foreground wedge (Datta et al., 2010; Morales et al., 2012; Hazelton et al., 2013; Thyagarajan et al., 2015). For this reason, many published power spectrum limits use data from short baselines only. For example, the power spectrum limits from PAPER (Ali et al., 2015; Cheng et al., 2018) were made using only the three shortest baselines in the array.

Ultimately, measurements of 21cm power spectra from the cosmic dawn will be used to constrain cosmological and astrophysical model parameters using such methods as Markov Chain Monte Carlo (MCMC) (Greig & Mesinger, 2015) and machine learning (Schmit & Pritchard, 2017; Gillet et al., 2019). Greig & Mesinger (2015) found that a 25% uncertainty in the measured power spectrum can increase uncertainty in model parameter estimates by a factor of a few, significantly degrading the constraints on the properties of galaxies during the EoR (although sampling the 21cm power spectrum over a wide range of redshifts can reduce the effect). Such uncertainty must be understood and accounted for in the era of precision cosmology.

In this paper, we argue that the sample variance of single-baseline power spectrum estimation does not fall below this 25% threshold for a typical data volume. In other words: even if all foreground and noise power can be properly mitigated, the fundamental uncertainty of the estimated 21cm power spectrum will be too high for it to be used for precision cosmology. Through instrument simulations of a mock EoR signal, we explore how this sample variance is affected by baseline length and beam width, and how it can be mitigated by combining data from multiple non-redundant baselines.

2 Power Spectrum Sensitivity

Many previous measures of instrument sensitivity estimated the power in Fourier space of thermal noise and residual foreground contamination, and made simplified approximations of sample independence to calculate the variance in each Fourier mode (e.g. Beardsley et al., 2013; Pober et al., 2014). The approximate integration time per Fourier mode is done by dividing the UV plane into bins of size given by the effective antenna aperture in wavelengths and calculating the motion of baselines with time and frequency. The variance is estimated by modeling thermal noise and residual foreground power and dividing through by this integration time. Sample variance is estimated as a component proportional to the expected power spectrum amplitude.

Beardsley et al. (2013) used a method like this to predict that the MWA Phase I would be capable of detecting the EoR power spectrum, but relied on optimistic assumptions as to the degree which residual foreground power contaminated the signal (Dillon et al., 2015). Pober et al. (2014) estimated the sensitivity of HERA in this way by combining the thermal noise uncertainty, from a formula presented in Pober et al. (2013), and residual foreground power from three foreground models.

Thyagarajan et al. (2013) considered the contributions of thermal noise, residual foreground power, and sample variance to the uncertainty of MWA Phase I observations in comparing two observing modes: a long integration on a single field, or (for an equal total observing time) fewer observations of multiple fields. As a phased array, each MWA tile is capable of pointing in discrete directions. This allows for a “drift and shift” observing strategy, where the array can approximately track a field by phasing to a position just west of it, letting the field drift through, then “shifting” the pointing west again. Thyagarajan et al. (2013) also took into account a more realistic model of the antenna power pattern (beam) in UV space, providing a more realistic look at how primary beam effects (especially sidelobe sources) contribute to foreground power leakage. Sample variance was assumed to drop off as one over the number of fields, i.e., for perfectly independent fields.

Trott (2014) performed a similar analysis for MWA Phase I, but instead of looking at discrete fields also directly considered the effects of covariance between visibilities in a constantly-evolving drift scan. In this case, the independence of observing fields is harder to determine. Trott (2014) found that the choice of ideal observing strategy depended on the available observing time and survey volume. The drift-scanning mode reduced sample variance the most for a given observing time, as one might expect since drift-scanning covers a larger area, but did not reduce thermal noise as significantly. For a drift scanning approach, they estimate a best signal to noise ratio of 9.6 (or an RMS error of about 10%) for the power spectrum amplitude.

This paper takes an exact approach to study the effects of sample variance on realistic 21cm observations for redundant arrays. By simulating the observations of an idealized interferometer over 24 hours of observing time, we can directly measure the variance of the measured power spectra and the covariance of delay-transformed visibilities. We simulate observations of a full sky model with a known power spectrum, and measure the sample variance of power spectra estimated from the simulated datasets.

We ignore the effects of foregrounds and thermal noise in this work, assuming that foreground power has been sufficiently removed through sky modeling techniques. Thermal noise can in principle be reduced as much as desired by including data from multiple days. The sample variance ultimately is the most fundamental uncertainty, because we are limited to one sky.

We introduce our simplified delay spectrum estimator, similar to that used by PAPER and will be used by HERA, and discuss its expected statistics in section 3. Section 4 further explores the statistics of the estimator through controlled numerical tests, examining how visibility correlation relates to the distribution of estimator values. We describe our full simulation code and sky model in section 5. Section 6 describes our main questions and the simulations used to address them. In section 7 we present the results of several sets of simulations, exploring the correlations of visibilities with time and the sample variance of power spectrum estimates.

3 The Delay Spectrum

The delay transform was introduced in (Parsons & Backer, 2009; Parsons et al., 2012b; Parsons et al., 2014) as an inverse Fourier transform along the frequency axis into “delay” space.

[TABLE]

Equation 2 gives the delay-transform of the visibility measured by baseline $\boldsymbol{b}$, with primary antenna beam $A$, at time $t$. $T(\hat{s},\nu,t)$ is the fluctuating component of the sky brightness temperature at frequency $\nu$ which has mean 0. Equation 2 is explicitly written here in instrument (topocentric) coordinates, such that the sky moves while the primary beam and fringe terms are fixed in time. The unit vector $\hat{s}$ represents a position on the sky in this topocentric frame.

Delay modes $\tau$ approximately correspond with Fourier modes parallel to the line of sight $k_{\parallel}$, while the baseline length $|\boldsymbol{u}|$ corresponds to perpendicular modes $k_{\perp}$. We can therefore estimate the power spectrum by cross-multiplying delay-transformed visibilities for redundant baselines (Parsons et al., 2014).

[TABLE]

The Fourier mode $k_{b,\tau}$ probed by this estimator depends on the baseline length $u$ (in wavelengths) as well as the delay mode $\tau$. An average is taken over baselines $\alpha,\beta$ that are in the same redundant group $G$. $X$ and $Y$ are a cosmological scaling factors with units of length per angle and length per frequency, respectively. $\lambda$ is the observed wavelength and $k_{B}$ is Boltzmann’s constant. $\Omega_{\text{pp}}$ is the integral of the primary beam squared, and $B$ is the bandwidth of the observation (see Appendix B of Parsons et al. 2014).

Cosmological isotropy requires that the power spectrum be independent of direction, allowing us to get a better estimate by averaging samples in time:

[TABLE]

where $t_{i}$ are the separate integration times of the instrument, of which we’re averaging together delay spectra from $N_{t}$ time samples.

For a single baseline ($\alpha=\beta$), we can combine eqs. 3 and 4 into a single expression.

[TABLE]

where $\Phi=(\lambda^{2}/2k_{B})^{2}X^{2}Y/B\Omega_{pp}$ is a scalar combining the cosmological and volumetric scaling factors of eq. 3. Note that cross-multiplying data for a single baseline will lead to a bias in the presence of thermal noise. Cross-multiplying redundant baselines, which measure the same signal but have independent noise realizations, avoids this bias. We can also avoid this bias by cross-multiplying measurements of the same LSTs across multiple nights. The analysis presesnted in this work doesn’t include thermal noise, so we can use eq. 5 as an unbiased estimator.

3.1 Statistics of $\widehat{P}(k)$

The delay-transformed visibilities may be written as a measurement vector $\boldsymbol{v}_{\tau}=(v_{0},\ldots,v_{N})$, with $v_{n}=\widetilde{V}(\tau,t_{n})$. For a Gaussian sky signal, $\boldsymbol{v}_{\tau}$ follows a complex multivariate gaussian distribution, with mean vector $\boldsymbol{\mu}=0$ and covariance matrix $\mathrm{C}=\left\langle\boldsymbol{v}_{\tau}\boldsymbol{v}_{\tau}^{\dagger}\right\rangle$. $\mathrm{C}$ is an $N_{t}\times N_{t}$ symmetric, real, positive semi-definite matrix.

We may write our estimator (5) as

[TABLE]

where $\mathrm{A}=(\Phi/N_{t})\mathcal{I}_{N}$ for $(N_{t}\times N_{t})$ identity matrix $\mathcal{I}_{N}$. We let $\widehat{P}_{N}$ stand for eq. 5 to simplify notation, with the understanding that it refers to a single baseline and delay mode, likewise dropping the $\tau$ subscript on $\boldsymbol{v}$. In this notation, $\widehat{P}_{N}$ is recognizable as a real quadratic form. In the case that visibilities are uncorrelated and have similar noise, eq. (5) follows a $\chi^{2}$ distribution with $N_{t}$ degrees of freedom. More generally, this quadratic form may be shown to be $\chi^{2}$ distributed with degrees of freedom $d=\text{rank}(\mathrm{AC})$ under the condition that the matrix $\mathrm{A}\mathrm{C}$ is idempotent (Mathai & Provost, 1992; Searle, 1971). We find experimentally that eq. 5 is well-described by a $\chi^{2}$ distribution (see section 4), despite a lack of idempotency of $\mathrm{AC}$. The underlying distribution of eq. 5 is not important to our results, but is of interest for interpreting the measured uncertainty. We will fit $\chi^{2}$ distributions to the measured power spectra.

The variance of the estimator may be found directly from the covariance matrix using Theorem 1 of Chapter 5 of Searle (1971). The proof presented there readily extends to the complex case with the result:

[TABLE]

For the delay-transformed visibilities, $\boldsymbol{\mu}=0$, so this expression reduces to the trace of the squared covariance matrix with the scalar factors. Recalling that $C_{ij}$ is Hermitian, it follows that the variance becomes the sum of the squares of all entries:

[TABLE]

For the case of uncorrelated visibilities with variances $\sigma^{2}$, this reduces to $\text{Var}[\widehat{P}_{N}]=\Phi^{2}\sigma^{2}/N$, as expected. To understand the variance of our estimator, therefore, we need to look at the covariance between visibility measurements at different times.

3.2 Visibility covariance

A full expression for the covariance between two delay-transformed visibilities has been worked out in Zhang et al. (2018) in detail. For a single baseline and drift-scanning visibilities, their expression reduces to the form

[TABLE]

Equation 9 assumes the power spectrum is approximately constant across the bandpass (i.e. ignoring light-cone effects). $\mathrm{R}_{ij}$ is a three-dimensional rotation matrix describing the motion of source from time $t_{i}$ to $t_{j}$ in the observer’s topocentric frame (e.g., altitude/azimuth).111A notational note – The subscripts on $\mathrm{R}_{ij}$ are not indexes of the matrix itself. Note that in the case of $t_{i}=t_{j}$, $\mathrm{R}_{ij}=\mathcal{I}$ and we recover the estimator (3) from eq. 9.

We can manipulate eq. 9 further to gain some better intuition. Let $\Theta_{\nu}$ represent the sky integral, such that,

[TABLE]

Then

[TABLE]

In this form, we can see that $\Theta_{\nu}$ is given by the overlap of the beam between the pointings, weighted by a fringe term. The fringe width is given by length of the baseline projected onto the displacement vector of the sky, $|\boldsymbol{b}^{\prime}|$.

If we choose a fixed phase center $\hat{s}_{0}$, we can rewrite eq. 10 under a flat-sky approximation. We define $\vec{b}^{\prime}$ as the projection of $\boldsymbol{b}^{\prime}$ onto the plane orthogonal to $\hat{s}_{0}$, and $\vec{l}=\hat{s}-(\hat{s}\cdot\hat{s}_{0})\hat{s}_{0}$ is the projection of the sky position vector into the sky plane.222We denote 2D vectors with overhead arrows. Finally, we denote by $\mathrm{\Gamma}_{ij}$ the transformation matrix of $\vec{l}$ corresponding with the 3D transformation $R_{ij}$.

In these terms, eq. 10 may be written

[TABLE]

In this form, the covariance is a 2D Fourier transform of a product of the primary beams, which may be expressed as the convolution of the Fourier-transformed primary beams, which we call beam kernels:

[TABLE]

This quantity depends only on the baseline, beam kernel width, and time between visibilities (lag). Since beam kernels are compact in Fourier space, the correlation is likewise limited in range, and can be characterized by a correlation time $t_{c}$. For our purposes, we define $t_{c}$ to be the full-width at half maximum of the correlation vs. lag. Note that in the case of nearly-Gaussian primary beam functions, the convolution in eq. 12 will give a Gaussian function, which can be characterized entirely by its FWHM.

The width of beam kernel, typically equal to the antenna size in wavelengths, scales inversely with the primary beam width on the sky. From eq. 12 we therefore expect that measurements from shorter baselines and narrower beams should stay correlated longer than those from long baselines and wide beams. This may seem counter-intuitive, so it’s worth some further discussion. One may think of the primary (antenna) beam as selecting a “patch” of sky that is integrated into a single visibility. In this picture, visibilities which sample overlapping parts of the sky should be correlated, and so the correlation time should be the time it takes the sky to pass through the primary beam. However, the exponential term in eq. 10 complicates this interpretation, and makes it more sensible to think about the Fourier-transformed formula eq. 12, which is more in line with what previous sensitivity analyses have done. We will show in section 7.1, however, that the Fourier-space interpretation of the correlation is not always valid.

Fig. 1a shows an example covariance matrix obtained from simulated data. The simulations used to make this data are discussed in section 5. This is the correlation in time for a single baseline, averaged over delay, and shows the expected symmetry due to the isotropy of the sky signal. If we bin the covariance matrices for each delay mode in time separation, or lag, we get Fig. 1b, which we call the correlation function. This shows that the central stripe in the covariance matrix has a nearly Gaussian shape.

3.3 Theoretical sample variance

We can construct an analytic function for the variance, eq. 8, under some assumptions for the visibility correlation. We first assume that the time covariance matrix is symmetric and constant parallel to the diagonal (also known as a symmetric Toeplitz matrix). This matrix symmetry is equivalent to the assumption that the EoR signal is isotropic since the correlation matrix depends only on the separation in time index,

[TABLE]

The factor of $\Phi$, from eq. 5, is included to absorb the cosmological and unit scalars into the definition of $\gamma_{n}$. We can express the sum in eq. 8 by counting the number of matrix entries at each $n$.

[TABLE]

To continue, we will treat $\gamma_{n}$ as a discretization of a continuous correlation function $\gamma(t)$. Assuming the time is discretized into time steps of length $\delta t$, such that $t_{n}=n\Delta t$, with a total time $T=N\delta t$, the sums in eq. 13 become:

[TABLE]

Approximating the correlation function as a Gaussian with FWHM $t_{c}$, we can evaluate these integrals. Letting $w_{c}=t_{c}/2.355$, the correlation function is

[TABLE]

Carrying out the integrals in eq. 13 with this Gaussian, we’re left with the following for the variance of the power spectrum estimator at averaging time $T$:

[TABLE]

Erf$[\cdot]$ is the error function. This expression holds only under the assumption that $\Delta t\ll T$. In section 9 we show that this approximate expression fits with the sample variances measured from simulated data.

4 Toy Model

To better understand the relationship between the distribution of $\widehat{P}_{N}$ and the correlation of visibility data, we ran a series of simple tests using mock data with a known correlation. We generate an ensemble of $M$ arrays of complex gaussian variates $\boldsymbol{x}=(x_{0},\ldots,x_{N_{t}})$, each with a length $N_{t}$, and convolve each with a normalized gaussian function to introduce a correlation. We then take progressively longer averages of the absolute square of these arrays to mimic the form of the power spectrum estimator (5).

[TABLE]

For each averaging length $l$, this gives us an ensemble of $M$ averages $\{q(l)\}$. We can fit a scaled $\chi^{2}$ distribution to this set using maximum likelihood estimation. The distribution function has the form

[TABLE]

where $\alpha$ is a scaling parameter and $k$ is the degree of freedom parameter, and $\Gamma(k)$ is the gamma function. When defined from a sum of squares of unit variance variables, the $\chi^{2}$ distribution has a mean equal to its degrees of freedom. In the form (19), the scaling parameter reflects that $x$ is not unit variance and that $q$ is an average, not a sum.

We expect that the degrees of freedom in the fits should be related to the correlation time, and we do indeed find that

[TABLE]

where $l_{c}$ is the full-width at half maximum of the convolution kernel. This is shown in Fig. 2. The scaling parameter $\alpha$ converges to $\alpha\rightarrow\sigma^{2}/(k-1)$, where $\sigma^{2}$ is the variance of $x$.

To make for a better comparison with our simulations, we define a time step $\Delta t=11$ $\sec$ and use $N_{t}=7854$ samples as described in Table 1. The convolution kernel width is set in terms of time as $l_{c}=t_{c}/\Delta_{t}$ for $t_{c}\in\{$20$\min$, 30$\min$, 50$\min$, 100\min$$\}. For these cases, we try to recover the input $t_{c}$ through the estimator

[TABLE]

where $k$ is the fitted degrees of freedom. Fig. 2 shows the result.

For the simple, short-range correlations of the delay-transformed visibilities we can still expect the $\widehat{P}_{N}$ to be $\chi^{2}$ distributed, and that we will be able to relate the correlation time. However, the convergence of $\hat{t}_{c}$ to the true value is slow, and depends on the number of samples going into each fit and how the correlation time compares with the total time available (24 hours). For longer correlation times and fewer samples, the convergence is worse.

Overall, we interpret the results of this toy model to mean that for the baseline types, beam widths, and averaging lengths under consideration, the distribution of $\widehat{P}_{N}$ will be best described by a $\chi^{2}$ distribution with relatively few degrees of freedom. The degrees of freedom are connected to the correlation time, but it may not be possible to reliably estimate $t_{c}$ from fitting the degrees of freedom. For the rest of this work, we use more robust methods to estimate the correlation time and compare with the results of a $\chi^{2}$ fit for consistency.

5 Simulation

Simulations were carried by directly evaluating the radio interferometry measurement equation in sky coordinates (Thompson et al., 2017), using an analytic and frequency-independent primary beam function. The simulator code is publicly available.333https://github.com/RadioAstronomySoftwareGroup/healvis We consider the ideal case of no foregrounds or thermal noise, and that our baselines perfectly measure the Stokes-I visibilities without any polarization leakage. By making these ideal assumptions, we can examine the effects of sample variance in isolation.

We integrate the measurement equation numerically across a spherical sky with uniform pixel areas $\omega$.

[TABLE]

The field of view (FoV) in eq. 21 is a circular selection of the sky centered on the zenith extended out to a chosen angular distance. For foreground simulations, this should be set to extend to the horizon, (see Thyagarajan et al. 2015), but since the choice of FoV increases the number of points to be integrated (especially near the horizon) increased FoV comes at the expense of increased runtime. We choose the FoV diameter to be $110^{\circ}$, at which distance the widest beam used in this analysis drops to 0.01% of its maximum.

For each time step, the pointing center is calculated using astropy coordinate frame transformations (Astropy Collaboration et al., 2013; Price-Whelan et al., 2018). The zenith angle and azimuth coordinates of each pixel in the FoV are calculated relative to this pointing center, ensuring that the sky rotates correctly throughout the simulation. Each pixel can therefore be thought of as small, uniform-brightness patch whose positions are in the international celestial reference frame (ICRS). This integration scheme ignores variation in the beam and fringe across the pixel, treating the beam and fringe as having a uniform value across each pixel equal to their values at the center of the pixel. Comparisons between simulations of different resolutions showed minimal effect from these assumptions.

We use a simple Gaussian function, $A(\theta,\phi)=\exp\left(-\theta^{2}/2\sigma^{2}\right)$ as our primary beam, where $\sigma$ sets the beam width and $\theta$ is the zenith angle. From here on, “beam width” refers to the FWHM of this Gaussian function ($\approx 2.355\sigma$). This is chosen for computational simplicity and to allow us to explore the effects of varied beam width directly. More realistic beam models, especially those obtained from calibrator sources or EM modeling, are important when modeling foregrounds, because bright sources moving through the beam and its sidelobes can introduce a spectral structure (Pober et al., 2016). Since the signal we analyze here is isotropic, there won’t be any drastic changes in integrated flux due to, say, a bright source near the edge of beam dropping into a null at lower frequency. We therefore don’t expect these subtle beam effects to be important here.

We note that, though well-defined in sky coordinates, this beam function does not directly correspond with a physical antenna on the ground, such that we are ignoring the potential impacts of beam chromaticity. Within the chosen bandpass, the effects of baseline chromaticity and the signal structure are much more important. Future work will refine these results with more realistic, chromatic beams.

5.1 Mock signal model

Instrument simulation tools are generally designed to read sky models in instrument coordinates (angle and frequency). Since frequency maps to distance, the survey volume forms a spherical shell with thickness corresponding with the bandwidth. This shell is divided into voxels with perpendicular area $r_{\nu}^{2}\omega$ and thickness $\delta r_{\nu}=|r(\nu+\delta\nu)-r(\nu)|$, where $r_{\nu}$ is the comoving distance corresponding with frequency $\nu$, $\omega$ is the pixel solid angle, and $\delta\nu$ is the frequency channel width. Within the HEALPix pixellization scheme, the Nside parameter sets the pixel solid angle to $\omega=4\pi/(12\times$Nside${}^{2})$ $\mathrm{sr}$, such that all pixels have the same area and pixel centers are evenly spaced (Gorski et al., 2005).

For our sky model, we select a brightness for each voxel sampled from a Gaussian distribution with frequency-dependent variance $\sigma^{2}(\nu)$. A space with equal voxel volumes has a flat power spectrum with amplitude $\sigma^{2}dV$, since there are no spatial correlations among the values. The comoving voxel volume changes with distance for two reasons. The radial thickness of the voxel changes because the frequency channels are evenly spaced and are nonlinearly related to comoving distance. In addition, the transverse comoving size of the voxel grows quadratically with distance because the solid angle remains the same. To account for this, we scale the variance at each frequency to ensure that the power spectrum amplitude is the same across the band.

[TABLE]

Here, $dV(\nu)=r_{\nu}^{2}\omega\delta r_{\nu}$ is the comoving voxel volume as a function of frequency. Quantities subscripted [math] refer to a reference frequency, chosen to be the center of the bandpass.

The sky model used here has a flat power spectrum, $P(k)=\sigma_{0}^{2}dV_{0}$, which has some advantages for further analysis, and is simple to generate in large quantities. The problem of generating a more realistic EoR model involves establishing spatial correlations that are isotropic and homogeneous within the spherical shell geometry, and is not trivial. It will be of interest to simulated data from EoR models that take into account both the spherical nature of the sky and the expected non-gaussian features of the late-stage EoR signal, especially since these nonlinearities can increase uncertainty on power spectrum measurements (Mondal et al., 2015). We leave full-sky modeling of EoR signals to future work.

6 Key Questions

We ran three sets of simulations to explore the following questions:

What is the relationship between correlation time and baseline length / beam width?

We explore the validity of eq. 9 by running a suite of short simulations with a variety of baseline lengths and beam widths, and calculate the correlation times from them. 2. 2.

How does the RMS error of single baseline delay spectra drop off with increasing averaging time for HERA-like and MWA-like beams?

Here we wish to know how well a single baseline can measure the EoR power spectrum through a delay transform method. 3. 3.

What is the lowest RMS error that can be achieved using all the elements of a 37 element hexagon, like MWA Phase II?

Combining data from multiple independent baselines can improve the sensitivity. Within the basic delay spectrum formalism, we are limited, however, to combining baselines of similar lengths, which correspond to the same $k_{\perp}$ mode. With all antennas in the 37-element array, how much does this binning improve the RMS error?

Table 1 summarizes the common parameters used by the simulations. Each simulation used the same bandpass, channel width, and array location. The latitude of the array does affect how baselines move through Fourier space over time, which may impact results. We use the HERA array position for all simulations, which is within 5∘ of latitude from the MWA site, so we expect our results apply to both cases. We use the MWA channel width, which is the finer of the two instruments, and use a 30 $\mathrm{MHz}$ bandpass. Most delay analyses focus on a fiducial 10 $\mathrm{MHz}$ band, to avoid redshift evolution effects, but for our purposes choosing a wider bandwidth increases the number of $k_{\parallel}$ probed.

For our sky models, we choose an Nside of 128 and a band-center pixel brightness variance \sigma_{0}^{2}=(31$$\mathrm{mK}$$)^{2}, taken from the expected EoR signal at $z\sim 13$, the bottom of our bandpass (Pritchard & Loeb, 2008). Since we’re looking at the relative error, the actual signal amplitude is inconsequential. The choice of Nside is a compromise with available computational resources. The simulator had 120GB of memory available, which was sufficient for a shell with an Nside of 256 (resolution of about 13$$\mathrm{\SIUnitSymbolArcminute}) and a 30MHz bandpass, or an Nside of 512 (resolution 6.8$$\mathrm{\SIUnitSymbolArcminute}) with a 10MHz bandpass. Both higher resolutions come at the expense of significantly longer run times, and since many tasks involved running hundreds of simulations we choose a modest resolution of Nside=128, which corresponds with an angular resolution of about 27.8$\mathrm{\SIUnitSymbolArcminute}$. This is still small relative to the smallest interferometric fringes in the simulations, and comparisons with higher resolution simulations showed no noticeable impact on results.

7 Results

7.1 Correlation time vs. Baseline and Beam Width

The first set of simulations was done to explore the effect of baseline length and beam width on visibility covariance. For 75 beam widths ranging from 15∘ to 50∘ we simulated an array of 50 coplanar baselines pointing East, with lengths evenly spaced from 20 to 100 $\mathrm{m}$. These simulations were repeated for 50 separate sky realizations, and run for 2 hours worth of simulated time.

For each simulated dataset, we delay-transform the visibility array and compute the covariance by cross-multiplying and averaging over delay. By treating the delays as an ensemble axis this way, we implicitly assume that all visibilities follow the same distribution at all delays, which we will show later to be a fair assumption due to the flatness of the power spectrum. This cross-multiplication gives us an $N_{t}\times N_{t}$ covariance matrix which is symmetric across both diagonals. We then bin the covariance matrix in lag, producing a correlation function that has a nearly Gaussian shape.

Fig. 3 shows the results of this set of simulations. As expected from eq. 9, correlation time decreases rapidly with increasing baseline length. For most baseline lengths, correlation time generally decreases with beam width. However, for beam widths less than 25∘ for baselines shorter than 50 $\mathrm{m}$, correlation time increases with beam width. This is unexpected, and counter to our interpretation in section 3.2. This is also the regime in which we find the shortest HERA and MWA baselines.

There are several factors at play that may be causing this change in behavior at short baselines and narrow beams. As discussed in section 3.2, if we ignore the fringe term in eq. 10 then we would expect the correlation time to increase with beam width. The fringe width is determined by the baseline length, and for short baselines will grow more slowly with time than for a long baseline. It may be that, for short baselines and narrow beams, the decoherence due to narrow beams has a greater effect than the fringe term.

Another possible reason for this turnover may come from our use of a Gaussian primary beams. For a physical antenna, the beam kernel width is approximately equal to the antenna diameter in wavelengths. A Gaussian beam has a Gaussian beam kernel, which does not have a finite width, and therefore would correspond with an antenna of infinite diameter and a diminishing response toward its edge. For short baselines, it may be that the antennas are effectively overlapping under this assumption. However, simulations with more realistic MWA and HERA beam models have hinted at this behavior on short baselines, and so we do believe this is a real effect.

The main takeaway here is that our intuition for correlation, as the overlap of beam kernels in Fourier space, typically leads to the behavior we expect: longer baselines and wider beams give us a shorter correlation time, and so visibilities will decorrelate faster with time and sample variance will drop off faster. However, the shortest baselines in HERA and the MWA do fall into the counter-intuitive regime where increasing beam width increases correlation time. We will recall this behavior on short baselines when interpreting the results of the 24 hour simulations exploring sample variance dropoff.

7.2 Sample Variance

The second set of simulations was done to measure the root-mean-square (RMS) error of the power spectrum estimator (5) for a range of observation times for individual HERA- and MWA-like baselines. By “observation time,” we mean the total contiguous simulated time that is averaged together in the estimator eq. 5. The RMS error $s_{k}(t)$, with different observation times $t$, as a fraction of power spectrum amplitude, for each $k=k_{b\tau}$ mode, is defined as

[TABLE]

where $P_{\text{theory}}=\sigma_{0}^{2}dV_{0}$ is the input power spectrum amplitude.

The MWA/HERA beam widths used here are set based on the true HERA/MWA beam widths at the lowest frequency of the selected bandpass, where they are widest. This is a conservative approximation, since the covariance is smaller for wider beams, so any error from the lack of beam chromaticity is likely to cause us to underestimate sample variance. The HERA beam width is set to $12.13^{\circ}$, which is the FWHM of an Airy disk pattern for a circular aperture of diameter 14.6$\mathrm{m}$ (the HERA dish diameter) at 100$\mathrm{MHz}$. Estimating a good effective beam width for the MWA is less straightforward, since each MWA tile is a phased array of 16 dipole antennas and so has a more complicated response than a parabolic dish. We choose a FWHM of $31.4^{\circ}$, extrapolating from the $23^{\circ}$ at 137$\mathrm{MHz}$ estimated in Neben et al. (2016).

We simulate 100 sky realizations with each beam width and with two baseline configurations. The first is an equilateral triangle of 14.6$\mathrm{m}$ baselines, one oriented East, the others pointing NW and NE. The second configuration is three long baselines: (136,143) points East and is 102.2$\mathrm{m}$, (2,136) points NW and is 139.2$\mathrm{m}$, and (2,143) points NE and is 115.9$\mathrm{m}$. These baselines were selected from the current 65 element layout of HERA, and are plotted in Fig. 4.

For each baseline, we obtain one delay spectrum per time step by doing an inverse FFT along the frequency axis, taking the absolute square, and applying the cosmological scaling factors as in (3). We want to measure the variance of progressively longer averages of these power spectra (i.e., vary $N_{t}$ in the sum of (5)). We consider a set of observation times $t\in T_{\text{avg}}$, which range from one time sample (11s) to the full day (24 hours).

To measure the statistics of the estimator, we need an ensemble of averaged power spectra. We build this ensemble by taking separate advantage of all parts of the dataset that independent of each other. For each observation time $t$, we partition the time array into non-overlapping segments of length $t$, such that $N_{\text{parts}}=\text{floor}(24\text{hrs}/t)$ where the “floor” function rounds its argument down to the nearest integer. The single time-step delay spectra are averaged within each segment, giving us $N_{\text{parts}}\times N_{\text{skies}}$ independent samples of power spectra averaged at a length $t$, where $N_{\text{skies}}=100$ is the number of sky realizations. We then measure the RMS error of each $P(k)$ across this ensemble, giving us the $s_{k}(t)$ of (23). This partitioning method naturally gives more samples at shorter averaging time, since there are more short non-overlapping segments within 24 hours than long. Using multiple sky realizations ensures that there are always at least $N_{\text{skies}}$ samples in the ensemble.

Fig. 5 shows the RMS of power spectrum estimates at different time observation times. The lines show the mean of $s_{k}(t)$ across all $k$, while the error shading show the standard error calculated across $k$.444Note that this is effectively the variance of the variance of the power spectrum, which is itself a variance. The discrete jumps in the shaded regions are due to the sharp changes in $N_{\text{part}}$ at certain observation times.

For the short baselines, the RMS drops off more quickly for the HERA beam (dashed) than for the MWA beam (solid). This is consistent with the conclusions of Fig. 3, since the 14.6$\mathrm{m}$ baselines are in the regime where increasing beam width increases correlation time, which means a slower drop-off. This behavior is reversed in the long baseline case. The black lines and grey shading are the corresponding results after binning power spectra in $k_{\perp}$. Averaging together the three baselines in each case drops the error by a factor of $\sim\sqrt{3}$, as expected for independent baselines.

The baseline length appears to have the most significant effect on RMS error, yet even with a 139$\mathrm{m}$ baseline the RMS error doesn’t get below $10\%$. The situation is much worse for the short baselines, which individually don’t reach below 19% for HERA beams or below 25% for MWA beams. Single-baseline estimators, particularly for short baselines, cannot be used reliably for precision cosmology.

7.3 Full array simulations

Results for individual baselines of specific lengths and orientations motivated us to determine the sample variance limits of a 37 element hexagonal array, i.e., the layout of the MWA Phase II East hexagon. Antenna elements are spaced by 14$\mathrm{m}$ in this configuration. Baselines are grouped by redundancy and one baseline is selected from each redundant group. We ran a series of simulations of 100 sky realizations, calculating the visibilities only for these selected baselines. Since there is no thermal noise in this simulation, visibilities from redundant baselines are identical. For the 37 element hexagonal array, there are 66 unique baseline types.

Each baseline samples the same range of $k_{\parallel}$ for a given $k_{\perp}$. Fig. 6 shows the RMS errors of single-baseline power spectrum estimates for the different baseline lengths in the array vs. averaging time. The discontinuities at 12 hours, 8 hours, and elsewhere are due to the sudden changes in the number of partitions available.

The spread of curves for the MWA beam, compared with the HERA beam curves, is due to the change in behavior observed in Fig. 3: for short baselines the narrower HERA beam has a shorter correlation time than the MWA, and vice versa for long baselines. The general trend of longer baselines having faster variance dropoff is clear for both instruments.

We can see from Fig. 6 that, for individual baselines in the 37-element hexagonal layout, $s_{k}$ stays above about 13% for MWA beams and above 15% for HERA beams. This is in the optimistic case that all 24 hours of LST can be used. With only 8 hours of data, which is how much time was used by Ali et al. (2015) for PAPER-64, the error for a single baseline of HERA is around 35%, and is upwards of 45% for the MWA.

By combining estimates from multiple baseline lengths, we can reduce this error dramatically. We can extend eq. 5 by binning separate $k_{b,\tau}$ modes. Each $k_{b,\tau}$ is a single delay mode for a single baseline, and coresponds with a point in the cylindrical power spectrum space $(k_{\perp},k_{\parallel})$. We can define bin in the spherical-averaged $k=\sqrt{k_{\parallel}^{2}+k_{\perp}^{2}}$, where $k_{\parallel}=\tau/Y$ and $k_{\perp}=b/(\lambda X)$, and average together $\widehat{P}(k_{b,\tau},N_{t};b)$ in each bin.

Choosing a set of $k$ bins is a subtle problem that depends on the expected shape of the power spectrum. For a highly-structured power spectrum, one needs to choose fine $k$ bins to resolve the structure. For a smoother spectrum, one can combine measurements into wider $k$ bins with little loss of information. Choosing wider bins means more measurements are combined, which can reduce sample variance, so there is a trade-off between the precision with which we can constrain the power spectrum shape and the per-bin sample variance. For a flat power spectrum model, like the ones we’re using, one could in principle average all $k$ modes with impunity, since there is no shape to constrain.

We choose to bin in the logarithmically spaced $k$ bins of the default configuration in 21cmFAST, which provide a representative range of scales and $k$ space resolution for describing the information in power spectra typical of the EoR (Mesinger et al., 2011). The bin edges are shown as contours in Fig. 7, over the power spectrum $P(k_{\perp},k_{\parallel})$ averaged 24 hours for a single sky realization.

Fig. 8 shows the RMS error in each $k$ bin, at different observation times and for the different instruments, using all baselines in the 37 element array for both HERA and MWA-like beams and both 8 and 24 long observations. The variance drops off with increasing $k$ due to the logarithmic bins, which become wider at high $k$ and therefore include more samples. We note that for all bins considered, the error is under 3.5%, and can be brought under a percent for the highest bins. Therefore, we expect that even with a relatively small number of antennas, a maximally redundant array can place high precision constraints on the power spectrum, as long as all baseline types are used in the measurement.

8 Comparing Correlation Time Measures

The correlation time $t_{c}$ specifies the approximate separation in time of independent samples of the power spectrum, and can in principle be calculated from the baseline length and primary beam model. This allows one to estimate the sample variance error for a given baseline and beam model directly from a formula like eq. 17, without running expensive simulations. In this section, we show that the correlation time as defined in section 7.1 gives the RMS error for longer averaging times, and is consistent with the number of degrees of freedom in a $\chi^{2}$ distribution, using the six baselines simulated in section 7.2. We also show that using the correlation time calculated here fits the RMS error vs. time plots, using eq. 17. The results are summarized in Fig. 9.

8.1 Binned covariance

The first $t_{c}$ measurement method is to calculate the covariance matrix, bin it in lag, and find the width of the central peak. Since we have $N_{\text{skies}}$ simulations in each configuration, we can calculate the covariance matrix separately for each delay mode. This is in contrast to the method in section 7.1, where we used the delay axis as an ensemble when computing covariance matrices. We switch to this method now because it is more consistent with the per-k RMS error measurement of section 7.2. The results of this method are plotted as triangles in Fig. 9. The points and error bars are the means and standard deviation, across $k$, of the measured correlation times.

We note that there is no technical reason why the FWHM of the covariance should set the correlation time. The correlation beyond the half-maximum distance is not negligible. For this reason, the FWHM is likely underestimating $t_{c}$. This explains why the triangles in Fig. 9 are a bit lower, though still within the uncertainty, of the results from the other methods.

8.2 RMS in the central limit

For a large number of independent samples, the estimator $\widehat{P}_{t}$ is approximately Gaussian-distributed under the central limit theorem. We therefore expect RMS of the power spectrum estimator should go roughly as $s_{k}(t)=1/\sqrt{M(t)}$, where $M$ is the number of independent time samples within the time $t$, for $t>>t_{c}$. From $M$ we can estimate $t_{c}$ as

[TABLE]

Equation 24b is the error on the estimate, propagating from the measured error $\Delta s_{k}(t)$. The squares in Fig. 9 show this estimate from the 24 hour averaging time, using the results in Fig. 5 with eqs. 24a and 24b.

8.3 $\chi^{2}$ distribution fit

The third method, discussed in section 4, is to fit a scaled $\chi^{2}$ distribution to the measured power spectra and treat the fitted degrees of freedom parameter $\kappa$ as the number of independent samples in the measurement, such that the correlation time is $t_{c}=t_{\text{avg}}/\kappa$555To address a possible confusion in terminology: In most cases, the phrase “$\chi^{2}$ fitting” refers to fitting model parameters to data using a chi-squared goodness-of-fit test. Here, we mean we are fitting the parameters of a $\chi^{2}$ distribution that describes the statistics of our data. We take the ensemble of $\widehat{P}_{N}$ measurements for a chosen observation time, and estimate the parameters of the scaled $\chi^{2}$ distribution separately for each $k$-mode using a maximum likelihood estimator. The square points on Fig. 9 show the mean of these per-k estimates, with error bars given by the standard error across $k$.

Two factors influence how well this method may be used. For shorter averaging times, the ensemble of $\widehat{P}_{N}$ measurements is much larger, making fewer samples to fit with the $\chi^{2}$ distribution. However, as Fig. 2 shows, for longer correlation times (e.g. the red curve), there is a risk of underestimating in this case. With much longer averaging times, there is the opposite problem. Since the other methods suggest the correlation times to be between 20 and 100 minutes, we look at estimates from 12 hours of averaging. Fig. 10 shows histograms of the averaged power spectra for three observation times, along with the fitted $\chi^{2}$ for that time. With longer averaging times, the degrees of freedom increases and the distributions become less skewed, as expected. We note that there is still a noticeable skew around 8 hours.

8.4 Theoretical RMS error

Fig. 11 shows the RMS vs. averaging time curves of Fig. 5, but with the theoretical sample variances for those baselines indicated with dashed black lines. This theoretical curve is obtained from eq. 17 by dividing out by $\gamma_{0}=P_{\text{theory}}$ and taking the square root. The $t_{c}=2.355w_{c}$ used here are the correlation times measured by $\chi^{2}$ distribution fitting (i.e. the squares in Fig. 9).

The good agreement shown in Fig. 11 indicates that the form of eq. 17 can be used to estimate the single-baseline sample variance at different averaging times, knowing only $t_{c}$, in the case that the correlation function is approximately Gaussian in time. The correlation time may be calculated directly from eq. 10 or eq. 9. It is thus possible to estimate the sample variance of an experimental setup without carrying out simulations. This may prove invaluable for future instrument design.

9 Discussion and Conclusions

We have used wide-field instrument simulations with Gaussian beams to estimate the statistical independence of delay-transformed visibilities over time, and calculate the uncertainty per-baseline for the simple delay-spectrum estimator of eq. 5. We have further looked further at the case of a 37 element hexagonal array, like the MWA Phase II, and measured the uncertainty of spherically-binned power spectrum estimates. This was done under the assumption that thermal noise and foreground power have been fully mitigated, and so represents a fundamental limit on the uncertainty that this estimator can reach.

Table 2 summarizes sample variance estimates in several of the simulated scenarios. The shortest baselines in HERA and the MWA cannot bring the sample variance below 20%, even in the ideal case that all 24 hours of LST are available. Spherically averaging power spectrum estimates across all baselines in the 37 element hexagonal array brings the RMS error under 3.5% for all cases.

As mentioned in section 7.3, choosing spherical-k bins for the power spectrum estimator does come with its own difficulties. In the case of a flat power spectrum, the only parameter we have to constrain is the amplitude, which does not depend on $k$. For a non-flat spectrum, however, binning in $k$ can potentially smooth over small scale features, which will propagate through to constraints on the model parameters. We chose the logarithmic bins of 21cmFAST as a typical binning used in 21cm analysis, but the actual choice of $k$ bins will affect the available improvement in sample variance.

We note also that these results largely relied on the Gaussianity of the tested signal. The EoR signal, especially at late times, is expected to be highly non-Gaussian due to the evolution of structure and appearance of ionized regions in the field. Mondal et al. (2015), and more recently Shaw et al. (2019), have shown that higher-order statistics such as the trispectrum enhance the errors of power spectrum estimates by 50% – 100% at redshifts up to $z\lesssim 10$. Future simulations with a more realistic EoR model will be able to examine the effect of these non-Gaussianities on realistic power spectrum estimates.

Nonetheless, our results bode well for HERA and the MWA’s usefulness for precision cosmology. As suggested by Greig & Mesinger (2015), getting the sample variance under 25% can be enough to constrain EoR parameters within a factor of a few, though the tolerable level of power spectrum uncertainty remains to be determined. As HERA construction continues, improvements in foreground modeling and removal will clean up more of Fourier space, allowing longer baselines to be used. Further, the split-hexagon final layout of HERA, discussed in Dillon & Parsons (2016), will add more non-redundant short baselines to the analysis, providing more independent samples of low $k_{\perp}$ modes. The MWA phase II also consists of another hexagonal array and a number of randomly-placed tiles, which greatly increases the numbers of available baselines. Most importantly, none of these sample variance uncertainties prevent the possibility of an EoR power spectrum detection, which is the most immediate goal of these instruments.

To summarize our main conclusions:

• •

Sample variance is an important factor on single baselines, and severely limits the precision of cosmological models that can be drawn from power spectrum measurements.

• •

Including measurements from longer baselines can help to mitigate the impact of sample variance, which drops dramatically with increasing baseline length. Combining measurements from multiple baseline lengths by binning across $k_{\perp}$ will further reduce the sample variance down to tolerable levels, while still avoiding the use of very long baselines.

• •

For a gaussian signal, the uncertainty on a delay spectrum estimator follows a scaled chi-squared distribution. The effective degrees of freedom of this distribution can cause a non-negligble skew, which may further bias parameter constraints by underestimating the error toward larger $P_{k}$ values. Analyses attempting to make cosmological parameter constraints should account for this potential skew.

• •

The sample variance of a Gaussian or near-Gaussian signal may be estimated from the time-correlation of visibilities. This correlation does not scale simply with beam width and baseline length, but can be evaluated analytically or numerically from knowledge of the primary beam shape and baseline vector.

• •

The time-correlation of delay-transformed visibilities depends on the primary beam width and the baseline length, but can behave in unexpected ways. The correlation time of short baselines will increase with increasing beam width for narrow beams, and decrease for wider beams.

10 Acknowledgements

This work was supported by National Science Foundation Grant Nos. 1636646 and 1613040, and NASA Grant 80NSSC18K0389. Thanks are owed to Miguel Morales, Cathryn Trott, Nithyanandan Thyagarajan, Michael Wilensky, Daniel Jacobs, Bryna Hazelton, Adrian Liu, and Bradley Greig for their useful advice and discussions. We thank the developers of scipy, matplotlib, astropy, healpy, and pyuvdata.