# Direct Optimal Mapping Image Power Spectrum and its Window Functions

Zhilei Xu MIT Kavli Institute, Massachusetts Institute of Technology,
Cambridge, MA  Honggeun Kim MIT Kavli Institute, Massachusetts Institute of
Technology, Cambridge, MA  Department of Physics, Massachusetts Institute of
Technology, Cambridge, MA  Jacqueline N. Hewitt MIT Kavli Institute,
Massachusetts Institute of Technology, Cambridge, MA  Department of Physics,
Massachusetts Institute of Technology, Cambridge, MA  Kai-Feng Chen MIT Kavli
Institute, Massachusetts Institute of Technology, Cambridge, MA  Department of
Physics, Massachusetts Institute of Technology, Cambridge, MA  Nicholas S.
Kern MIT Kavli Institute, Massachusetts Institute of Technology, Cambridge, MA
Department of Physics, Massachusetts Institute of Technology, Cambridge, MA
NHFP Hubble Fellow  Elizabeth Rath  MIT Kavli Institute, Massachusetts
Institute of Technology, Cambridge, MA  Department of Physics, Massachusetts
Institute of Technology, Cambridge, MA  Ruby Byrne  Cahill Center for
Astronomy and Astrophysics, California Institute of Technology, Pasadena, CA
Adélie Gorce  Department of Physics and McGill Space Institute, McGill
University, 3600 University Street, Montreal, QC H3A 2T8, Canada  Zachary E.
Martinot  Department of Physics and Astronomy, University of Pennsylvania,
Philadelphia, PA  Joshua S. Dillon Department of Astronomy, University of
California, Berkeley, CA  Bryna J. Hazelton Department of Physics, University
of Washington, Seattle, WA  eScience Institute, University of Washington,
Seattle, WA  Adrian Liu Department of Astronomy, University of California,
Berkeley, CA  Department of Physics and McGill Space Institute, McGill
University, 3600 University Street, Montreal, QC H3A 2T8, Canada  Miguel F.
Morales Department of Physics, University of Washington, Seattle, WA  Zara
Abdurashidova  Department of Astronomy, University of California, Berkeley, CA
Tyrone Adams  South African Radio Astronomy Observatory, Black River Park, 2
Fir Street, Observatory, Cape Town, 7925, South Africa  James E. Aguirre
Department of Physics and Astronomy, University of Pennsylvania, Philadelphia,
PA  Paul Alexander  Cavendish Astrophysics, University of Cambridge,
Cambridge, UK  Zaki S. Ali  Department of Astronomy, University of California,
Berkeley, CA  Rushelle Baartman  South African Radio Astronomy Observatory,
Black River Park, 2 Fir Street, Observatory, Cape Town, 7925, South Africa
Yanga Balfour  South African Radio Astronomy Observatory, Black River Park, 2
Fir Street, Observatory, Cape Town, 7925, South Africa  Adam P. Beardsley
School of Earth and Space Exploration, Arizona State University, Tempe, AZ
Department of Physics, Winona State University, Winona, MN  Gianni Bernardi
INAF-Istituto di Radioastronomia, via Gobetti 101, 40129 Bologna, Italy
Department of Physics and Electronics, Rhodes University, PO Box 94,
Grahamstown, 6140, South Africa  South African Radio Astronomy Observatory,
Black River Park, 2 Fir Street, Observatory, Cape Town, 7925, South Africa
Tashalee S. Billings  Department of Physics and Astronomy, University of
Pennsylvania, Philadelphia, PA  Judd D. Bowman School of Earth and Space
Exploration, Arizona State University, Tempe, AZ  Richard F. Bradley  National
Radio Astronomy Observatory, Charlottesville, VA  Philip Bull Jodrell Bank
Centre for Astrophysics, University of Manchester, Manchester, M13 9PL, United
Kingdom  Department of Physics and Astronomy, University of Western Cape, Cape
Town, 7535, South Africa  Jacob Burba Jodrell Bank Centre for Astrophysics,
University of Manchester, Manchester, M13 9PL, United Kingdom  Steven Carey
Cavendish Astrophysics, University of Cambridge, Cambridge, UK  Chris L.
Carilli National Radio Astronomy Observatory, Socorro, NM 87801, USA  Carina
Cheng  Department of Astronomy, University of California, Berkeley, CA  David
R. DeBoer Radio Astronomy Lab, University of California, Berkeley, CA  Eloy de
Lera Acedo  Cavendish Astrophysics, University of Cambridge, Cambridge, UK
Matt Dexter  Radio Astronomy Lab, University of California, Berkeley, CA  Nico
Eksteen  South African Radio Astronomy Observatory, Black River Park, 2 Fir
Street, Observatory, Cape Town, 7925, South Africa  John Ely  Cavendish
Astrophysics, University of Cambridge, Cambridge, UK  Aaron Ewall-Wice
Department of Astronomy, University of California, Berkeley, CA  Department of
Physics, University of California, Berkeley, CA  Nicolas Fagnoni  Cavendish
Astrophysics, University of Cambridge, Cambridge, UK  Randall Fritz  South
African Radio Astronomy Observatory, Black River Park, 2 Fir Street,
Observatory, Cape Town, 7925, South Africa  Steven R. Furlanetto Department of
Physics and Astronomy, University of California, Los Angeles, CA  Kingsley
Gale-Sides  Cavendish Astrophysics, University of Cambridge, Cambridge, UK
Brian Glendenning  National Radio Astronomy Observatory, Socorro, NM  Deepthi
Gorthi Department of Astronomy, University of California, Berkeley, CA
Bradley Greig School of Physics, University of Melbourne, Parkville, VIC 3010,
Australia  Jasper Grobbelaar  South African Radio Astronomy Observatory, Black
River Park, 2 Fir Street, Observatory, Cape Town, 7925, South Africa  Ziyaad
Halday  South African Radio Astronomy Observatory, Black River Park, 2 Fir
Street, Observatory, Cape Town, 7925, South Africa  Jack Hickish  Radio
Astronomy Lab, University of California, Berkeley, CA  Daniel C. Jacobs School
of Earth and Space Exploration, Arizona State University, Tempe, AZ  Austin
Julius  South African Radio Astronomy Observatory, Black River Park, 2 Fir
Street, Observatory, Cape Town, 7925, South Africa  MacCalvin Kariseb  South
African Radio Astronomy Observatory, Black River Park, 2 Fir Street,
Observatory, Cape Town, 7925, South Africa  Joshua Kerrigan Department of
Physics, Brown University, Providence, RI  Piyanat Kittiwisit Department of
Physics and Astronomy, University of Western Cape, Cape Town, 7535, South
Africa  Saul A. Kohn Department of Physics and Astronomy, University of
Pennsylvania, Philadelphia, PA  Matthew Kolopanis School of Earth and Space
Exploration, Arizona State University, Tempe, AZ  Adam Lanman  Department of
Physics, Brown University, Providence, RI  Paul La Plante Department of
Astronomy, University of California, Berkeley, CA  Department of Physics and
Astronomy, University of Pennsylvania, Philadelphia, PA  Anita Loots  South
African Radio Astronomy Observatory, Black River Park, 2 Fir Street,
Observatory, Cape Town, 7925, South Africa  David Harold Edward MacMahon
Radio Astronomy Lab, University of California, Berkeley, CA  Lourence Malan
South African Radio Astronomy Observatory, Black River Park, 2 Fir Street,
Observatory, Cape Town, 7925, South Africa  Cresshim Malgas  South African
Radio Astronomy Observatory, Black River Park, 2 Fir Street, Observatory, Cape
Town, 7925, South Africa  Keith Malgas  South African Radio Astronomy
Observatory, Black River Park, 2 Fir Street, Observatory, Cape Town, 7925,
South Africa  Bradley Marero  South African Radio Astronomy Observatory, Black
River Park, 2 Fir Street, Observatory, Cape Town, 7925, South Africa  Andrei
Mesinger Scuola Normale Superiore, 56126 Pisa, PI, Italy  Mathakane Molewa
South African Radio Astronomy Observatory, Black River Park, 2 Fir Street,
Observatory, Cape Town, 7925, South Africa  Tshegofalang Mosiane  South
African Radio Astronomy Observatory, Black River Park, 2 Fir Street,
Observatory, Cape Town, 7925, South Africa  Steven G. Murray School of Earth
and Space Exploration, Arizona State University, Tempe, AZ  Abraham R. Neben
MIT Kavli Institute, Massachusetts Institute of Technology, Cambridge, MA
Department of Physics, Massachusetts Institute of Technology, Cambridge, MA
Bojan Nikolic  Cavendish Astrophysics, University of Cambridge, Cambridge, UK
Hans Nuwegeld  South African Radio Astronomy Observatory, Black River Park, 2
Fir Street, Observatory, Cape Town, 7925, South Africa  Aaron R. Parsons
Department of Astronomy, University of California, Berkeley, CA  Nipanjana
Patra Department of Astronomy, University of California, Berkeley, CA
Samantha Pieterse  South African Radio Astronomy Observatory, Black River
Park, 2 Fir Street, Observatory, Cape Town, 7925, South Africa  Nima Razavi-
Ghods  Cavendish Astrophysics, University of Cambridge, Cambridge, UK  James
Robnett  National Radio Astronomy Observatory, Socorro, NM 87801, USA  Kathryn
Rosie  South African Radio Astronomy Observatory, Black River Park, 2 Fir
Street, Observatory, Cape Town, 7925, South Africa  Peter Sims Department of
Physics and McGill Space Institute, McGill University, 3600 University Street,
Montreal, QC H3A 2T8, Canada  Craig Smith  South African Radio Astronomy
Observatory, Black River Park, 2 Fir Street, Observatory, Cape Town, 7925,
South Africa  Hilton Swarts  South African Radio Astronomy Observatory, Black
River Park, 2 Fir Street, Observatory, Cape Town, 7925, South Africa
Nithyanandan Thyagarajan Commonwealth Scientific and Industrial Research
Organisation (CSIRO), Space & Astronomy, P. O. Box 1130, Bentley, WA 6102,
Australia  Pieter van Wyngaarden  South African Radio Astronomy Observatory,
Black River Park, 2 Fir Street, Observatory, Cape Town, 7925, South Africa
Peter K. G. Williams Center for Astrophysics, Harvard & Smithsonian,
Cambridge, MA  American Astronomical Society, Washington, DC  Haoxuan Zheng
Department of Physics, Massachusetts Institute of Technology, Cambridge, MA

###### Abstract

The key to detecting neutral hydrogen during the epoch of reionization (EoR)
is to separate the cosmological signal from the dominating foreground
radiation. We developed direct optimal mapping (Xu et al., 2022) to map
interferometric visibilities; it contains only linear operations, with full
knowledge of point spread functions from visibilities to images. Here we
present an FFT-based image power spectrum and its window functions based on
direct optimal mapping. We use noiseless simulation, based on the Hydrogen
Epoch of Reionization Array (HERA) Phase I configuration, to study the image
power spectrum properties. The window functions show $<10^{-11}$ power leakage
from the foreground-dominated region into the EoR window; the 2D and 1D power
spectra also verify the separation between the foregrounds and the EoR.
Furthermore, we simulated visibilities from a $uv$-complete array and
calculated its image power spectrum. The result shows that the foreground–EoR
leakage is further suppressed below $10^{-12}$, dominated by the tapering
function sidelobes; the 2D power spectrum does not show signs of the horizon
wedge. The $uv$-complete result provides a reference case for future 21 cm
cosmology array designs.

21-cm lines (690), Early universe (435), Radio interferometry (1346),
Reionization (1383)

##  1 Introduction

Recent cosmological observations have established the standard cosmological
model — $\Lambda$CDM cosmology (Bennett et al., 1996; Riess et al., 1998;
Bennett et al., 2013; Hinshaw et al., 2013; Planck Collaboration et al., 2016,
2020). The cosmic microwave background (CMB) measures the early universe at
redshift $\sim$1100 (Fixsen et al., 1996; Hu & Dodelson, 2002; Staggs et al.,
2018); the late universe measurement gives statistical properties of the
universe below redshift 10 (Anderson et al., 2014; Alam et al., 2017; Abbott
et al., 2022; More et al., 2023). However, we have hardly observed the
universe from redshift 1100 to 10, including dark ages, cosmic dawn, and the
epoch of reionization (EoR) (Furlanetto et al., 2006; Pritchard & Loeb, 2012).

With measurements from the initial ($z\sim 1100$) and final ($z<10$)
conditions, the $\Lambda$CDM cosmology describes the initial perturbations’
growth from linear to non-linear until the ignition of first stars and
galaxies (Hogan & Rees, 1979; Madau et al., 1997). Our goal is to observe this
process by measuring radiation from the dominant baryonic content in the
universe — neutral hydrogen. The atomic hyperfine structure of neutral
hydrogen leads to emission or absorption of 21 cm radiation in the restframe.
The 21 cm radiation, redshifted, is observed at different frequency channels,
tracing the distribution of the baryonic matter during EoR, cosmic dawn, and
dark ages. Direct observations covering the above epochs are the ultimate goal
of 21 cm cosmology. The 21 cm observations have the potential to measure a
wide range of cosmological history, testing the current cosmological model
(Mao et al., 2008). The result will reveal more details of the cosmological
evolution and the universe’s contents, including dark matter and dark energy.

In practice, wavelengths of high-redshift 21 cm radiation are redshifted to
$>$ 2 m beyond $z\sim 8.5$; correspondingly, $>$ 230 m baselines are required
to achieve $\sim 0.5$∘ angular resolutions. Therefore, high-redshift 21 cm
measurements are conducted with interferometers for angular resolutions
(Parsons et al., 2010; Tingay et al., 2013; van Haarlem et al., 2013; DeBoer
et al., 2017). Interferometers are also cost-effective in achieving large
collecting areas, which eventually determines the measurement sensitivity.
Using interferometers, we aim to measure the distribution of neutral hydrogen
radiation across the sky with power spectra. Image-based power spectra are
measured in both CMB (Bennett et al., 2013; Planck Collaboration et al., 2020)
and other 21 cm experiments (Dillon et al., 2014, 2015; Beardsley et al.,
2016; Trott et al., 2016; Patil et al., 2017; Barry et al., 2019; Li et al.,
2019; Rahimi et al., 2021). Here we present one power spectrum estimator based
on direct optimal mapping (Xu et al., 2022, hereafter X22).

In this paper, we study the image power spectrum based on one interferometer
array — the Hydrogen Epoch of Reionization Array (HERA) (Pober et al., 2014;
DeBoer et al., 2017) — and calculate the image power spectrum and its window
functions. The frequency coverage provides a tomographic measurement of the
Epoch of Reionization (EoR), constraining the evolution of neutral hydrogen
spin temperature and ionization fraction during EoR (Furlanetto et al., 2006;
Morales & Wyithe, 2010; Pritchard & Loeb, 2012; Liu & Shaw, 2020). These
measurements also constrain sources for the reionization process, including UV
and X-ray properties of high-redshift galaxies (Ewall-Wice et al., 2016; Greig
et al., 2016; Kern et al., 2017).

Section 2 introduces the image power spectrum concept and presents the
mathematical formalism of the image power spectrum and its window functions.
In Section 3, we first introduce the image cube, constructed with direct
optimal mapping, and present the results for the window functions, 2D power
spectrum, and 1D power spectrum. In Section 4, we further explore the
potential of this image power spectrum by simulating visibilities from a
$uv$-complete array and investigating its image power spectrum. We discuss
future work in Section 5 before concluding in Section 6. We use the WMAP nine-
year cosmology (Bennett et al., 2013) throughout this paper.

##  2 Image Power Spectrum

The key challenge for 21 cm cosmology is to detect the EoR signals with the
presence of foregrounds, which have intensity at least four orders of
magnitude brighter than the cosmological signals (Liu et al., 2014).
Fortunately, the foreground and the EoR signals have different frequency
properties with pixels: foregrounds have smooth frequency spectra while the
EoR signals have both smooth and fast-changing components. One essential
question is how well different analyses can separate EoR from smooth
foregrounds, and detect EoR in the region of Fourier space expected to be free
of foreground emission, called the _EoR window_.

The raw data from an interferometer are the cross-correlations among antenna
pairs, called visibilities. Visibilities are measured at different frequency
channels, among which the instrument intrinsically has different spatial
responses, called instrumental chromaticity. This instrumental chromaticity
converts spatial structures into the frequency axis, causing the visibility
frequency spectrum to contain not only the pixel frequency spectrum but also
the features leaked in from spatial structures. If we calculate the frequency
Fourier Transform of the visibilities, we see not only the frequency spectrum
of sky pixels but also their spatial structures, mostly from the foregrounds.
The latter component contributes to the ‘wedge’ feature in the 2D power
spectrum (Datta et al., 2010; Parsons et al., 2012, 2012; Liu et al., 2014,
2014). The technique of measuring the power spectrum by Fourier Transforming
the visibility’s frequency axis is called the delay power spectrum. (Parsons
et al., 2012).

Alternatively, part of the instrumental chromaticity can be accounted for by
analyzing frequency spectra of _pixels_ , instead of visibilities. Therefore,
this approach is called image power spectrum (Morales & Hewitt, 2004; Liu &
Tegmark, 2011; Dillon et al., 2013). Analyzing sky pixels mitigates the space-
to-frequency leakage from the visibilities. Morales et al. (2019) investigated
the two approaches and concluded that, assuming uniform $uv$ coverage, the
image power spectrum achieves lower foreground leakage into the wedge than the
delay power spectrum.

However, working with pixels requires forming a 3D image cube from the
visibilities. Fundamentally, the imaging process presents the visibilities in
a linear combination that minimize the instrumental chromaticity. How much
instrumental chromaticity can be accounted for with imaging depends on how
much information we possess to reconstruct the sky. If we had perfect
knowledge of the primary beam, its effect could be eliminated in imaging,
removing the primary beam chromaticity; if we had complete $uv$-coverage, we
could reconstruct each sky pixel without degeneracy, removing chromaticity
from interferometric grating lobes.

Our image power spectrum is based on direct optimal mapping (X22). In contrast
to other imaging algorithms (Högbom, 1974; Clark, 1980; Cornwell, 2008; Rau &
Cornwell, 2011), our mapping method does not include iterative processes. In
addition, we map the visibility data at their original $uv$ locations without
gridding (Sullivan et al., 2012). The algorithm involves only linear
operations, robustly tracking the transfer function in the process; both the
point spread function and the pixel covariance matrix can be accurately
calculated during the mapping process. We then calculate the power spectra
from the image cube, and the point spread function is used to calculate the
power spectrum window functions.

In practice, interferometers rarely achieve complete and uniform $uv$
coverage, meaning that the observations do not contain sufficient information
to recover pixel intensities without ambiguity. This ambiguity shows up as
grating lobes in images, which keeps us from achieving the ideal spectral
leakage mitigation. Now the questions are how much mitigation one can achieve
by image power spectrum with a realistic interferometer, and how much more one
can achieve with a $uv$-complete interferometer. This paper addresses the two
questions in the following sections.

###  2.1 Power Spectrum Definition

We first present the power spectrum formula and discuss the symmetry
assumptions for constructing a power spectrum and binning a 3D power spectrum
into 2D and 1D.

With a 3D image cube, we calculate its 3D discrete Fourier Transform (DFT)
following the convention in Mesinger & Furlanetto (2007):

$\widetilde{m}(\mathbf{k})=\frac{V}{N}\sum m(\mathbf{r})e^{-i\mathbf{k\cdot
r}},$ (1)

where $m(\mathbf{r})$ contains the intensity of the image cube and
$\mathbf{r}$ is the 3D spatial vector, $\widetilde{m}(\mathbf{k})$ is the DFT
of the image cube and $\mathbf{k}$ is the 3D wavenumber vector, $V$ is the
physical size of the image cube (in the unit of Mpc${}^{3}\cdot$h-3), and $N$
is the number of voxels.

The power spectrum is calculated by squaring the DFT result and dividing out
the total volume of the image cube (Liu & Shaw, 2020)

$P(\mathbf{k})=\frac{\langle\widetilde{m}(\mathbf{k})^{*}\cdot\widetilde{m}(\mathbf{k})\rangle}{V}.$
(2)

For spatially homogeneous signals, like cosmological signals from the early
universe, all cross terms from
$\langle\widetilde{m}(\mathbf{k^{\prime}})^{*}\cdot\widetilde{m}(\mathbf{k})\rangle$,
when $\mathbf{k^{\prime}}\neq\mathbf{k}$, average to zero. Our observable
universe is one realization from its underlying power spectrum. The angle
brackets indicate that the true power spectrum is defined as the average over
sufficient realizations of the universe. The angle brackets only apply to the
variables with randomness. For a deterministic variable, the angled value
equals itself, say the image cube volume $\langle V\rangle=V$.

The power spectrum is originally expressed in 3D $\mathbf{k}$ space. We bin
the 3D power spectrum to 2D (or 1D) under different symmetry assumptions. For
the 2D binning, the two dimensions perpendicular to the line-of-sight (LOS)
are circularly binned into $k_{\bot}$, assuming the two perpendicular
dimensions share the same statistical properties. Together with
$k_{\parallel}$, parallel to LOS, we form the 2D $k_{\bot}-k_{\parallel}$
power spectrum, which is the natural space for foreground avoidance (Datta et
al., 2010; Morales et al., 2012; Parsons et al., 2012). For the 1D binning, we
assume isotropy of the signal, and all three dimensions are spherically binned
into one dimension; the 3D vector $\mathbf{k}$ is condensed to only its length
$k$. Since the cosmological signals are believed to be homogeneous and
isotropic at large scales, the 1D binning is often used to describe the power
spectrum of the cosmological signals.

Homogeneity and isotropy in cosmological signals ensure the validity of the
above steps; however, when the foregrounds are involved, which can be neither
homogeneous nor isotropic, one should be aware of the underlining assumptions
as we reduce power spectrum dimensions.

###  2.2 Power Spectrum Estimator

We form an estimated image cube with the direct optimal mapping algorithm.
Here we quickly review the method and refer to more details presented in X22.
We first create the data model of the visibilities as

$\mathbf{d}_{f_{i}}=\mathbf{A}_{f_{i}}\mathbf{m}_{f_{i}}+\mathbf{n}_{f_{i}},$
(3)

where $\mathbf{d}_{f_{i}}$ is the measured data at the frequency channel
$f_{i}$, $\mathbf{m}_{f_{i}}$ is the true sky emission at the frequency
channel, $\mathbf{A}_{f_{i}}$ is the measurement matrix which describes how
the interferometer integrates sky emission to produce visibilities, and
$\mathbf{n}_{f_{i}}$ represents noise in the visibilities.

With this data model, the optimized sky recovery is described in a linear
operation

$\mathbf{\hat{m}}_{f_{i}}=\mathbf{D}_{f_{i}}\mathbf{A}_{f_{i}}^{\dagger}\mathbf{N}_{f_{i}}^{-1}\mathbf{d}_{f_{i}},$
(4)

where $\mathbf{\hat{m}}_{f_{i}}$ is the estimated sky map,
$\mathbf{D}_{f_{i}}$ is a normalization matrix, and $\mathbf{N}_{f_{i}}$ is
the noise covariance matrix calculated as
$\mathbf{N}_{f_{i}}=\langle\mathbf{n}_{f_{i}}^{\dagger}\mathbf{n}_{f_{i}}\rangle_{n}$.
The $\langle...\rangle_{n}$ operator means that we are averaging over
different visibility noise realizations. Combining the above two equations, we
define the point spread function (PSF) matrix $\mathbf{P}_{f_{i}}$

$\displaystyle\langle\mathbf{\hat{m}}_{f_{i}}\rangle_{n}$ $\displaystyle=$
$\displaystyle\mathbf{D}_{f_{i}}\mathbf{A}_{f_{i}}^{\dagger}\mathbf{N}_{f_{i}}^{-1}(\mathbf{A}_{f_{i}}\mathbf{m}_{f_{i}}+\langle\mathbf{n}_{f_{i}}\rangle_{n})$
(5) $\displaystyle=$
$\displaystyle(\mathbf{D}_{f_{i}}\mathbf{A}_{f_{i}}^{\dagger}\mathbf{N}_{f_{i}}^{-1}\mathbf{A}_{f_{i}})\,\mathbf{m}_{f_{i}}$
$\displaystyle\equiv$ $\displaystyle\mathbf{P}_{f_{i}}\,\mathbf{m}_{f_{i}},$

where we use the property that noise averages to zero
($\langle\mathbf{n}_{f_{i}}\rangle_{n}=0$) in the above derivation, and the
last line defines the PSF matrix as
$\mathbf{P}_{f_{i}}=\mathbf{D}_{f_{i}}\mathbf{A}_{f_{i}}^{\dagger}\mathbf{N}_{f_{i}}^{-1}\mathbf{A}_{f_{i}}$.
The map pixels are chosen to be smaller than the synthesized beam to ensure
that the measurement resolution is limited by the instrument. The subscript
$f_{i}$ in all the variables indicates that they are for one frequency
channel. Sky images are generated for each frequency channel, and frequency
channels are converted into LOS distances. Combining the 2D sky images at
different LOS distances yields a 3D image cube — $\mathbf{\hat{m}}$. It is
related to the true 3D image cube $\mathbf{m}$ by the 3D PSF matrix
$\mathbf{P}$. We aggregate the 2D image PSF matrix $\mathbf{P}_{f_{i}}$ to
form the 3D PSF matrix with $\mathbf{P}_{f_{i}}$ as block-diagonal matrices

$\mathbf{P}=\mathrm{diag}(\mathbf{P}_{f_{1}},...,\mathbf{P}_{f_{i}},...,\mathbf{P}_{f_{N}}).$
(6)

Here we do not consider frequency-frequency correlations. The aggregated
$\mathbf{P}$ maps the true 3D image cube to the estimation

$\langle\mathbf{\hat{m}}\rangle_{n}=\mathbf{Pm},$ (7)

where the 3D image cube $\mathbf{m}$ and its estimation $\mathbf{\hat{m}}$ are
flattened to column vectors.

The image cube estimation is then tapered along the three dimensions. The
tapering apodizes hard boundaries, suppressing ringing in Fourier space.
Meanwhile, we treat all the power within each voxel as a single point source
(Liu & Tegmark, 2011; Dillon et al., 2013), and we will correct for the voxel
window function later. With the Fourier Transform convention in Equation 1 and
2, we define the quadratic estimator of the image cube to be

$\hat{q}_{\alpha}=\frac{V}{N^{2}}\mathbf{\hat{m}}^{\dagger}\mathbf{R}^{\dagger}\mathbf{E}_{\alpha}\mathbf{R}\mathbf{\hat{m}},$
(8)

where $\mathbf{R}$ is the 3D tapering function of the image cube. Then, we
define

$\mathbf{E}_{\alpha}=\mathbf{c}^{\dagger}_{\alpha}\mathbf{c}_{\alpha},$ (9)

where $\mathbf{c}_{\alpha}$ is the 3D discrete Fourier Transform operator at
$\mathbf{k}_{\alpha}$

$\mathbf{c}_{\alpha}=(e^{-i\mathbf{k}_{\alpha}\cdot\mathbf{r}_{1}},...,e^{-i\mathbf{k}_{\alpha}\cdot\mathbf{r}_{j}},...,e^{-i\mathbf{k}_{\alpha}\cdot\mathbf{r}_{N}}).$
(10)

Elements within $\mathbf{c}_{\alpha}$ cover all voxels $\mathbf{r}_{i}$ within
the image cube; multiplying $\mathbf{c}_{\alpha}$ on an image cube vector
gives the Fourier Transform of the image cube at $\mathbf{k}_{\alpha}$. On the
contrary, the conjugate transpose operator $\mathbf{c}_{\alpha}^{\dagger}$ is
the inverse 3D discrete Fourier Transform

$\mathbf{c}_{\alpha}^{\dagger}=(e^{i\mathbf{k}_{\alpha}\cdot\mathbf{r}_{1}},...,e^{i\mathbf{k}_{\alpha}\cdot\mathbf{r}_{j}},...,e^{i\mathbf{k}_{\alpha}\cdot\mathbf{r}_{N}})^{T}.$
(11)

Neither $\mathbf{c}_{\alpha}$ nor $\mathbf{c}_{\alpha}^{\dagger}$ contains DFT
normalization factors; instead, the normalization factors are explicitly
represented in the leading factors of Equation 8.

Coming back to the quadratic estimator $\hat{q}_{\alpha}$, plugging in the
$\mathbf{m}$ and $\mathbf{\hat{m}}$ relation of Equation 7, the quadratic
estimator is written as

$\langle\hat{q}_{\alpha}\rangle_{n}=\frac{V}{N^{2}}\mathbf{m}^{\dagger}\mathbf{P}^{\dagger}\mathbf{R}^{\dagger}\mathbf{E}_{\alpha}\mathbf{R}\mathbf{P}\mathbf{m}.$
(12)

With the true sky brightness, we define the voxel covariance
$\mathbf{C\equiv\langle m\,\mathbf{m}^{\dagger}\rangle}$. The
$\langle...\rangle$ operator indicates we are averaging over different sky
realization, as initially defined in Equation 2. This covariance matrix is
associated with the power spectrum, and can be expressed as a linear
combination of the band powers (Dillon et al., 2013; Liu & Shaw, 2020)

$\mathbf{C\equiv\langle
m\,\mathbf{m}^{\dagger}\rangle}=\frac{1}{V}\mathbf{\sum_{\beta}}p_{\beta}\mathbf{E_{\beta}}\cdot|\Phi(\mathbf{k_{\beta}})|^{2},$
(13)

The above equation adds band powers across the $\mathbf{k}$-space, and the
normalization factor $1/V$ is the $\mathbf{k}$-space resolution for the
summation. In addition, the $\Phi(\mathbf{k_{\beta}})$ function is the voxel
window function at the $\mathbf{k}_{\beta}$ position. We use the 3D boxcar
profile as the real-space voxel profile, therefore, $\Phi(\mathbf{k_{\beta}})$
is expressed as the product of three sinc functions:

$\Phi(\mathbf{k})\equiv j_{0}\left(\frac{k_{x}\Delta
x}{2}\right)j_{0}\left(\frac{k_{y}\Delta
y}{2}\right)j_{0}\left(\frac{k_{z}\Delta z}{2}\right),$ (14)

where $j_{0}(x)=\mathrm{sin}(x)/x$, $k_{x},k_{y},k_{z}$ are the three
components of the $\mathbf{k}$ vector, and $\Delta x,\Delta y,\Delta z$ are
the real-space resolution of the image cube for each dimension (Dillon et al.,
2013).

Equation 12 has scalars on both sides, so we add the trace operator on both
sides. With the cyclic property of the trace, we can circularly shift the
matrices on the right-hand side. Also using the above equation, we rewrite
Equation 12

$\displaystyle\langle\hat{q}_{\alpha}\rangle$ $\displaystyle=$
$\displaystyle\frac{V}{N^{2}}\mathrm{tr}[\langle\mathbf{m}^{\dagger}\mathbf{P}^{\dagger}\mathbf{R}^{\dagger}\mathbf{E}_{\alpha}\mathbf{R}\mathbf{P}\mathbf{m}\rangle]$
(15) $\displaystyle=$
$\displaystyle\frac{V}{N^{2}}\mathrm{tr}[\mathbf{P}^{\dagger}\mathbf{R}^{\dagger}\mathbf{E}_{\alpha}\mathbf{R}\mathbf{P}\langle\mathbf{m}\mathbf{m}^{\dagger}\rangle]$
$\displaystyle=$
$\displaystyle\frac{V}{N^{2}}\mathrm{tr}[\mathbf{P}^{\dagger}\mathbf{R}^{\dagger}\mathbf{E}_{\alpha}\mathbf{R}\mathbf{P}\frac{1}{V}\mathbf{\sum_{\beta}}p_{\beta}\mathbf{E_{\beta}}|\Phi(\mathbf{k_{\beta}})|^{2}]$
$\displaystyle=$
$\displaystyle\frac{1}{N^{2}}\mathbf{\sum_{\beta}}\mathrm{tr}[\mathbf{P}^{\dagger}\mathbf{R}^{\dagger}\mathbf{E}_{\alpha}\mathbf{R}\mathbf{P}\mathbf{E_{\beta}}|\Phi(\mathbf{k_{\beta}})|^{2}]p_{\beta}$
$\displaystyle\equiv$ $\displaystyle\sum_{\beta}H_{\alpha\beta}p_{\beta}.$

Please note that we are averaging over both noise realizations and sky
realizations above. Here we define the power spectrum response function
$\mathbf{H}$, which transfer the true sky power to our quadratic estimator.
Element-wise, $\mathbf{H}$ is written as

$\displaystyle H_{\alpha\beta}$ $\displaystyle=$
$\displaystyle\frac{|\Phi(\mathbf{k_{\beta}})|^{2}}{N^{2}}\mathrm{tr}(\mathbf{P^{\dagger}R^{\dagger}E_{\alpha}RPE_{\beta}})$
(16) $\displaystyle=$
$\displaystyle\frac{|\Phi(\mathbf{k_{\beta}})|^{2}}{N^{2}}\mathrm{tr}(\mathbf{P^{\dagger}R^{\dagger}c_{\alpha}^{\dagger}c_{\alpha}RPc_{\beta}^{\dagger}c_{\beta}}).$

Using the matrix trace cyclic property again, we can rewrite the above
equation

$\displaystyle H_{\alpha\beta}$ $\displaystyle=$
$\displaystyle\frac{|\Phi(\mathbf{k_{\beta}})|^{2}}{N^{2}}\mathrm{tr}(\mathbf{c_{\beta}P^{\dagger}R^{\dagger}c_{\alpha}^{\dagger}c_{\alpha}RPc_{\beta}^{\dagger}})$
(17) $\displaystyle=$
$\displaystyle\frac{|\Phi(\mathbf{k_{\beta}})|^{2}}{N^{2}}\mathrm{tr}[(\mathbf{c_{\alpha}RPc_{\beta}^{\dagger})^{\dagger}(c_{\alpha}RPc_{\beta}^{\dagger}})]$
$\displaystyle=$
$\displaystyle|\Phi(\mathbf{k_{\beta}})|^{2}\,\mathrm{tr}[(\frac{\mathbf{c_{\alpha}}}{\sqrt{N}}\mathbf{RP}\frac{\mathbf{c_{\beta}^{\dagger}}}{\sqrt{N}})^{\dagger}(\frac{\mathbf{c_{\alpha}}}{\sqrt{N}}\mathbf{RP}\frac{\mathbf{c_{\beta}^{\dagger}}}{\sqrt{N}})]$
$\displaystyle=$
$\displaystyle|\Phi(\mathbf{k_{\beta}})|^{2}\,|\frac{\mathbf{c_{\alpha}}}{\sqrt{N}}\mathbf{RP}\frac{\mathbf{c_{\beta}^{\dagger}}}{\sqrt{N}}|^{2},$

The last line indicates the normalization for DFT: we divide forward and
background DFT by $\sqrt{N}$. Equation 17 is the formula to calculate the
response matrix $\mathbf{H}$.

The $\mathbf{H}$ matrix links the quadratic estimator $\mathbf{\hat{q}}$ and
the true power spectrum $\mathbf{p}$ as

$\mathbf{\langle\hat{q}\rangle=Hp}.$ (18)

The quadratic estimator $\mathbf{\hat{q}}$ is an intermediate product of the
power spectrum estimation, and we the response matrix $\mathbf{H}$ to convert
it to the power spectrum estimator $\mathbf{\hat{p}}$. Equation 17 also shows
that $\mathbf{H}$ is determined by $\mathbf{RP}$ — the transfer function in
image space. Essentially, $\mathbf{H}$ and $\mathbf{RP}$ represent the same
relation in two domains.

###  2.3 Window Functions

Window functions describe, in power spectrum space, how the true sky power is
mapped into the measured power spectrum. Ideally, window functions are delta
functions, which is not always achievable. In practice, one measured band
power takes contribution from all parts of the true sky power. The
contribution is characterized by a set of weights, which is the window
function.

Continuing with Equation 18, if we could construct a matrix $\mathbf{M}$ to
invert $\mathbf{H}$, unambiguously recovering the underlining band power
$\mathbf{p}$ from $\mathbf{\hat{q}}$, we would equivalently obtain delta-
function-like window functions. In practice, $\mathbf{H}$ is generally not
invertible for interferometers because of their incomplete $uv$ coverage.
Therefore, different forms of $\mathbf{M}$ are constructed for various
purposes (Seljak, 1998; Tegmark et al., 2002; Liu & Tegmark, 2011; Ali et al.,
2015; HERA Collaboration et al., 2022; Kern & Liu, 2021). Here we sum the
columns for each row in the $\mathbf{H}$ matrix and divide out the sum for
each element in the corresponding row. Mathematically, we form a diagonal
matrix $\mathbf{M}$ with each element as the inverse of the row
sum111Analogous to normalizing the weights for a weighted sum.

$\mathbf{M}=\mathrm{diag}(1/\sum_{\beta}H_{1,\beta},...,1/\sum_{\beta}H_{i,\beta},...,1/\sum_{\beta}H_{N,\beta}).$
(19)

The full-rank $\mathbf{M}$ matrix takes our quadratic estimator
$\mathbf{\hat{q}}$ to the estimated band powers $\mathbf{\hat{p}}$ (Tegmark,
1997; Liu et al., 2014)

$\mathbf{\hat{p}=M\hat{q}}.$ (20)

With this definition of $\mathbf{M}$, the window function matrix is defined as
$\mathbf{W=MH}$, the normalized power is expressed as

$\mathbf{\langle\hat{p}\rangle=M\langle\hat{q}\rangle=MHp}=\mathbf{Wp}.$ (21)

The window function matrix $\mathbf{W}$ is fully defined by $\mathbf{H}$

$\displaystyle W_{\alpha\beta}$ $\displaystyle=$
$\displaystyle\sum_{\gamma}M_{\alpha\gamma}H_{\gamma\beta}$ (22)
$\displaystyle=$ $\displaystyle M_{\alpha\alpha}H_{\alpha\beta}$
$\displaystyle=$
$\displaystyle\frac{H_{\alpha\beta}}{\sum_{\delta}H_{\alpha\delta}}.$

The above equation also shows that weights in each row of $\mathbf{W}$ sums to
one. Despite being normalized, the weights’ distribution varies across
different power spectrum estimators. Accordingly, the measured power spectrum
depends on the window function; we do not expect different power spectrum
estimators to recover exactly the same power spectrum.

Window functions in the delay power spectrum were studied in Parsons et al.
(2012); Liu et al. (2014); Gorce et al. (2023); here, we investigate window
functions from this image power spectrum.

##  3 Image Power Spectrum Application

Equipped with the above formalism, we investigate the properties of the image
power spectrum. We map simulated noiseless visibilities to one image cube and
calculate its image power spectrum along with the window functions. We compare
the window functions with the delay spectrum ones to understand the difference
between the two estimators. We show the 2D/1D power spectrum to evaluate the
foreground separation.

###  3.1 The Image Cube

We use the HERA validation simulations (Aguirre et al., 2022), which is based
on HERA Phase I array configuration with 33 unflagged antennas (HERA
Collaboration et al., 2022, 2023). We choose three datasets with different sky
models, including EoR-only, foregrounds-only, and EoR+foregrounds. Please note
that the EoR signals in these simulations are boosted compared to the fiducial
model (Aguirre et al., 2022). We employ these three variations to analyze
their image power spectra. The simulation includes neither noise nor
systematics. We select the data from one time-integration when the zenith is
at 24.2∘ in R.A.; and we select a 17.5 MHz frequency range from 150.34 MHz to
167.84 MHz (with 0.098 MHz resolution). After Blackman-Harris tapering, the
effective bandpass is around 8 MHz, within which cosmological evolution is
insignificant. The frequency range and resolution are consistent with Band 2
defined in HERA Collaboration et al. (2022). The data contain both East-West
and North-South polarization, we only analyze the East-West polarization in
this paper.

We use the Fast Fourier Transform (FFT) algorithm to calculate the DFT, which
requires a regular grid. Therefore, we create a regular R.A./Dec. grid in the
sky-plane. However, the pixels do not have equal areas on a curved surface;
their solid angle changes with the pixel declination,

$\Delta\Omega_{\mathrm{px}}=\Delta\alpha\Delta\delta\,\mathrm{cos}(\delta),$
(23)

where $\alpha,\delta$ are R.A. and Dec., $\Delta\Omega$ is the pixel solid
angle. From the above equation, we can see the pixel solid angle changes the
slowest around $\delta=0$ since
$d\Delta\Omega_{\mathrm{px}}/d\delta=\Delta\alpha\Delta\delta[-\mathrm{sin}(\delta)]=0$
at $\delta=0$. Therefore, we pixelize the sky with $\delta=0$ across the field
center to minimize the pixel size change. This is equivalent to rotating the
sky patch to the horizon for pixelization.

For each frequency channel, the visibilities are mapped with the direct
optimal mapping algorithm. The grid is chosen to center at 24.2∘ in R.A. and
-30.7∘ in Dec., with 0.5∘ resolution along both directions. The 0.5∘
resolution is chosen to be half of the size of the $\sim$1∘ synthesized beam.
The final map contains $32\times 16=512$ pixels, covering $16^{\circ}\times
8^{\circ}=128$ square degrees. The top panel of Figure 1 shows the map at the
central frequency channel — 159.04 MHz. The mapping repeats for the 180
frequency channels.

Figure 1: Original and tapered maps at 159.04 MHz. We use the direct optimal
mapping algorithm (X22) to map the noiseless simulation data to regular
R.A./Dec. grids. The original map consists of $32\times 16=512$ pixels,
covering $16^{\circ}\times 8^{\circ}=128$ square degrees. After tapering with
the Blackman-Harris function along three dimensions, the off-center signals
are highly attenuated to suppress ringing structures from Fourier Transform.
We select the center frequency for plotting where the frequency tapering
effect is minimal. In total, we map 180 maps in the frequency range from
150.34 MHz to 167.84 MHz.

The frequency channels are converted to redshifts and to comoving distances in
the unit of $\mathrm{Mpc\cdot h^{-1}}$. Uniform sampling in frequency leads to
nonuniform sampling in comoving distances, because of the nonlinear relation
between frequency and comoving distance. The distance intervals differ up to
$\sim$5.5% across the whole range, we chose the average as the resolution
along the LOS direction. In addition, comoving distances determine the
physical spacing in the sky-plane. Within our frequency range, the comoving
distance changes up to $\sim$3.2%; we use the average LOS comoving distance to
calculate the sky-plane resolution. Finally, the voxel size is calculated by
multiplying the physical resolutions along three dimensions, which is $\sim
3.5\times 10^{3}\,\mathrm{Mpc^{3}\cdot h^{-3}}$. Table 1 collects the related
parameters of the image cube. Please note the difference in comoving
resolutions between the LOS and the sky-plane: the frequency channels provide
much finer comoving resolution compared to the angular solution.

Table 1: Parameters for the Image Cube

Parameters | LOS | R.A. | Dec.
---|---|---|---
$N_{\mathrm{pixel}}$ | 180 | 32 | 16
Range | 17.5 MHz | 16∘ | 8∘
Resolution | 0.098 MHz | 0.5∘ | 0.5∘
Range ($\mathrm{Mpc\cdot h^{-1}}$) | 205.4 | 1772 | 886
Resolution ($\mathrm{Mpc\cdot h^{-1}}$) | 1.21 | 55.4 | 55.4

Note. — Image cube parameters shown in measurement units, including
angular/frequency resolution, and the related comoving distances. We use the
average comoving distance to calculate sky-plane resolution; we use the
average LOS interval as the LOS resolution.

###  3.2 Window Function Results

Figure 2: Window functions for image power spectrum and delay power spectrum.
The top panels show image window functions at three $(k_{\bot},k_{\parallel})$
points; the bottom panels show delay window functions at similar
$(k_{\bot},k_{\parallel})$ locations (Gorce et al., 2023). Locations of the
native $(k_{\bot},k_{\parallel})$ bins are marked with white plus signs (the
symbols are out of the plot range in the first column). All panels share the
same logarithmic scales in two axes; they also share the same color range from
1 to $10^{-12}$. The dashed gray lines shows the horizon wedge and the solid
white lines show 200 ns beyond the wedge. Both sets of window functions are
based on the same array configuration and observation parameters. The
$k_{\bot}-k_{\parallel}$ binning is different from the image and delay power
spectrum; white space is displayed where delay window functions do not have
data. All the window functions are sum-normalized. The top and right axes show
the corresponding baseline and delay values for $k_{\bot}$ and $k_{\parallel}$
respectively. The image window functions show lower power leakage along the
$k_{\parallel}$ direction ($<10^{-11}$) compared to the delay window
functions, which is critical for foreground avoidance; the spread along the
$k_{\bot}$ direction is from the limited sky patch of the image cube.

Now we look at the image power spectrum window functions for the HERA Phase I
configuration. Section 2.3 gives the theoretical formalism of the window
function; here we present the actual calculation.

The aggregated $\mathbf{P}$ matrix is a square matrix with
$(N_{\mathrm{pixel}}\times N_{\mathrm{freq.}})^{2}=(512\times 180)^{2}\approx
8\times 10^{9}$ elements. Storing and operating this matrix is the most
computationally expensive step in our image power spectrum analysis. The
tapering matrix $\mathbf{R}$ applies tapering functions along the three
dimensions of the original image cube. Here we use the Blackman-Harris
(4-term) tapering function along the three dimensions. A tapered map,
perpendicular to LOS, is shown in the bottom panel of Figure 1.

To calculate $\mathbf{H}$, we pick individual columns in the $\mathbf{RP}$
matrix, which are 3D vectors flattened to 1D. We reshape the 1D columns back
into 3D and calculate their 3D FFT. The result is then flattened to replace
the same column. The above process is repeated for all columns of
$\mathbf{RP}$. Correspondingly, inverse FFT is performed across all rows from
the previous result. Then, we square the individual elements of the matrix.
Finally, we calculate the pixel window function $|\Phi(\mathbf{k})|^{2}$ and
apply the result to all rows of the previous matrix to obtain $\mathbf{H}$.

Summing up columns in $\mathbf{H}$, we calculate the $\mathbf{M}$ matrix
according to Equation 19. The $\mathbf{M}$ matrix normalizes the quadratic
estimator to get the power spectrum estimator as in Equation 20. Finally, the
window function matrix $\mathbf{W}$ is calculated as $\mathbf{W=MH}$.

Figure 3: 2D image power spectra from EoR-only, foregrounds-only, and
EoR+foregrounds simulations. The three power spectra were calculated from
simulated visibilities. The power spectra share the same $k_{\bot}$,
$k_{\parallel}$ range as well as the color range; the dashed line shows the
horizon wedge and the solid line shows 200 ns beyond the wedge. The EoR power
spectrum shows power distributed across the whole space, with more power at
low-$k_{\parallel}$. The foregrounds’ power is concentrated within
low-$k_{\parallel}$ region, with an amplitude four orders of magnitude higher
than that from the EoR power spectrum. However, there is a window not
contaminated by the foregrounds in high-$k_{\bot}$ regions. Finally, the
EoR+foregrounds power spectrum shows the result with both sky components:
foreground features dominate low-$k_{\parallel}$ and EoR features prevail in
the high-$k_{\parallel}$ EoR window.

The window function matrix $\mathbf{W}$ is defined in 3D band powers; for each
$\mathbf{k}$ bin, the corresponding window function shows the contribution
from all the 3D $\mathbf{k}$ bins. The dimension of the window functions can
be reduced with symmetry, similar to the power spectrum. The two sky-plane
dimensions are circularly binned into $k_{\bot}$; the LOS direction is kept
unchanged as $k_{\parallel}$. Values from related 3D $\mathbf{k}$ bins are
averaged into ($k_{\bot}$, $k_{\parallel}$) space with equal weights. The
$k_{\parallel}$ direction is divided into 90 linear bins from 0 to 2.71
$\mathrm{h}\cdot\mathrm{Mpc}^{-1}$; the $k_{\bot}$ direction is divided into
16 linear bins from $7.4\times 10^{-3}$ to $7.8\times
10^{-2}\,\mathrm{h}\cdot\mathrm{Mpc}^{-1}$.

The top panels of Figure 2 show image power spectrum window functions at three
($k_{\bot},k_{\parallel}$) points. Along the $k_{\bot}$ direction, the window
functions show wide spreading. This is because we currently only consider a
small sky patch as shown in Figure 1. The small sky patch and the sky-plane
tapering functions generate the extended $k_{\bot}$ kernels.

Along the $k_{\parallel}$ direction, power leakage along the $k_{\parallel}$
direction is suppressed below $10^{-11}$. This is critical for separating EoR
from foregrounds: the smooth foregrounds only occupy low-$k_{\parallel}$
regions while the EoR signal takes up the entire $k_{\bot}-k_{\parallel}$
space; therefore, the detectability of the EoR signal lies in the assumption
that the foreground power can be contained only in low-$k_{\parallel}$
regions. The window functions in Figure 2 directly show how much power leaks
from low-$k_{\parallel}$ to high-$k_{\parallel}$.

The bottom panels of Figure 2 show window functions of the same datasets at
similar ($k_{\bot},k_{\parallel}$) locations from the delay power spectrum
(Gorce et al., 2023). In Gorce et al. (2023), the same Blackman-Harris
tapering is applied along the visibility frequency axis. The $k_{\parallel}$
leakage is approximately three orders of magnitude higher than the image power
spectrum, due to the spatial structures leaking into the frequency axis, as
discussed in Section 2. However, along the $k_{\bot}$ direction, the window
functions display a more compact kernel compared to the image power spectrum.
This is because visibilities integrate the whole sky and the delay spectrum
analyzes visibilities directly, equivalently measuring the full sky weighted
by the primary beam; compared to the small sky patch in the image power
spectrum, the delay spectrum window functions show more compact kernels along
$k_{\bot}$.

Because we cannot fully deconvolve the power spectrum response function, the
measured power spectrum receives contributions across all $\mathbf{k}$ bins
from the true band powers, whose distribution is presented in window
functions. The window functions show that the image and delay power spectrum
are different estimators of the true sky band powers.

###  3.3 2D Power Spectrum

The previous section examines the window functions from three representative
$(k_{\bot},k_{\parallel})$ bins, now we investigate the full 2D image power
spectra.

We first calculate the quadratic estimator $\mathbf{\hat{q}}$ with Equation 8.
We FFT and inverse FFT the tapered image cube and multiply the results
together element-wise, then we multiply the $V/N$ factor222The other $N$ in
the denominator is included in the FFT and iFFT normalization. to obtain the
3D quadratic estimator $\mathbf{\hat{q}}$. We use the matrix $\mathbf{M}$ to
normalize the quadratic estimator for the power spectrum estimator (Equation
20). After obtaining the normalized 3D power spectrum estimator, we bin it to
the $(k_{\bot},k_{\parallel})$ space by averaging the related 3D band power.
In the averaging step, we use equal weights. The power spectrum in the 2D
$(k_{\bot},k_{\parallel})$ space are shown in Figure 3.

The EoR-only power spectrum shows power across the entire
$k_{\bot}-k_{\parallel}$ space. Not much change is observed along the
$k_{\bot}$ direction while there is a clear decreasing trend from
low-$k_{\parallel}$ to high-$k_{\parallel}$. The peak of the 2D power spectrum
is around $10^{10}\mathrm{mK^{2}\cdot h^{-3}Mpc^{3}}$; the amplitude is highly
amplified in the simulation. The foregrounds-only power spectrum shows a
different situation: the power is constrained within the low-$k_{\parallel}$
regions, mostly within 200 ns beyond the horizon wedge. The peak of the
foregrounds-only power spectrum rises as high as $10^{14}\,\mathrm{mK^{2}\cdot
h^{-3}Mpc^{3}}$, four orders of magnitude higher than the simulated EoR-only
peak. However, its power does not leak into the high-$k_{\parallel}$ region
with $>10^{8}$ dynamic range. This is consistent with the window function
result above: the window functions show $<10^{-11}$ leakage from
low-$k_{\parallel}$ to high-$k_{\parallel}$, which ensures the overwhelming
power from the smooth foregrounds does not affect the EoR window.

Finally, we come to the power spectrum of EoR+foregrounds. Please note that
this power spectrum is not obtained by simply adding the previous two power
spectra; instead, we start from visibilities, form the image cubes, and
calculate the 2D power spectrum. In this power spectrum, the
low-$k_{\parallel}$ region is dominated by the foregrounds, resembling the
features in the foregrounds-only power spectrum; the high-$k_{\parallel}$
region is filled by EoR signals. The clear separation between the foregrounds
and the EoR further illustrates the detectability of the EoR with the image
power spectrum.

Figure 4: 1D power spectrum measured from three simulation datasets. The
foreground signal dominates in the low-$k$ regions; the boosted EoR signal
dominates in the high-$k$ region. Also shown is the EoR-only spectrum measured
by the delay power spectrum. Both the image power spectrum and the delay power
spectrum agree with the input power, represented in a dashed line. The
deviation of the measurement at high-$k$ is due to aliasing from DFT (Aguirre
et al., 2022). The deviation from a straight line in the mid-$k$ range is due
to the small sky patch selected for the image power spectrum.

###  3.4 1D Power Spectrum

The EoR power spectrum is often collapsed into 1D assuming it is homogeneous
and isotropic.333Please refer to Section 2.1 for discussions on symmetry and
power spectrum dimension reduction. We present the 1D power spectrum in this
section by binning the original 3D power spectrum. The binning is based on
equal-weight averaging, similar to the 2D power spectrum.

Figure 4 shows the 1D power spectrum from EoR-only, foregrounds-only, and
EoR+foregrounds simulations. The EoR+foregrounds 1D power spectrum closely
follows the foregrounds-only spectrum in low-$k$ and the EoR-only spectrum in
high-$k$. In high-$k$, the EoR signal is more than five orders of magnitude
stronger than that of the foregrounds. The foreground avoidance capability of
the image power spectrum is illustrated in this figure. Also plotted is the 1D
EoR-only power spectrum measured with the delay power spectrum, and the input
EoR power. The image power spectrum and delay power spectrum both recover the
input EoR power.

Figure 5: 1D Power Spectrum from the evolution of 21 cm structure (EOS)
simulation. We downsample and convolve the EOS image cube and measure the 1D
image power spectrum. The normalized 1D power spectrum agrees with the input
after being adjusted by six orders of magnitude from the normalization step.

In addition, we use the Evolution Of 21 cm Structure (EOS)444URL:
http://homepage.sns.it/mesinger/EOS.html simulation (Sobacchi & Mesinger,
2014) to further validate the image power spectrum 1D power spectrum. The EOS
project provides image cubes and also their input 1D power spectrum.
Therefore, we can calculate the 1D image power spectrum and compare it with
the input.

The EOS simulation provides image cubes with 1.6 Gpc (1.11 $\mathrm{Gpc\cdot
h^{-1}}$) and 1024 bins on each side. We use one image cube at redshift 5.76.
Along the x-axis, we downsample the EOS image cube into 32 bins; along the
y-axis, we also downsample the EOS image cube into 32 bins but select the
middle 16 bins; along the z-axis, we select the first 180 bins without
downsampling. We average the intensity among associated voxels during
downsampling. This gives us the $180\times 32\times 16$ image cube similar to
the DOM image cube. Then we convolve the downsampled EOS image cube with the
$\mathbf{P}$ matrix to get the instrument-convolved image cube.

We measure the image power spectrum and bin the 3D power spectrum into 1D.
Figure 5 shows the input 1D power spectrum and the normalized 1D power
spectrum from DOM. The measurement recovers the input 1D power spectrum across
the majority of the $k$, with some mismatch at the low-$k$ end. The
inconsistency at the low-$k$ end comes from the edge effect of the finite EOS
image cube. The finite EOS image cube contains no power beyond its size,
beyond the lower limit of $k_{\bot}$; the window functions average down the
power close to the $k_{\bot}$ lower limit. Also shown in Figure 5 is the 1D
power spectrum before normalization. The comparison shows the normalization
correctly rescales the power spectrum over six orders of magnitude.

##  4 Complete UV-coverage

Figure 6: Window functions for a $uv$-complete array. The three native $k$
bins, $k_{\bot}$/$k_{\parallel}$ range, and the color range are the same as in
Figure 2. With the complete $uv$ coverage, the power leakage from low
$k_{\parallel}$ regions is further suppressed by the improved $uv$ coverage
but not completely eliminated in the given color range. The spread along the
$k_{\bot}$ direction is also reduced from the complete $uv$ coverage.  Figure
7: 2D power spectrum from a $uv$-complete array. The $k_{\parallel}$ range and
the color range are identical to the ones in Figure 3. However, the $k_{\bot}$
range is extended given the enlarged $uv$ coverage. The extended part in
$k_{\bot}$ is shaded. The horizon wedge is not present with the complete $uv$
coverage, even after extending the $k_{\bot}$ range. The $uv$-complete 2D
power spectrum shows less power compared to the previous foregrounds+EoR
simulation. This is expected because the $uv$-complete simulation assumes
different sky models, points at a different sky patch, and covers a different
frequency band.

The HERA array configuration does not have complete $uv$ coverage (Dillon &
Parsons, 2016). The incomplete $uv$ coverage comes with ambiguity in imaging;
this is clearly demonstrated as grating lobes in the PSF (X22). The size and
location of the grating lobes change with frequency, leaking spatial power
into the frequency axis. Therefore, with a $uv$ incomplete array, we do not
achieve the full advantage of the image power spectrum (Morales et al., 2019).

To evaluate how much improvement we can expect if we had complete $uv$
coverage, we simulate visibilities from a $uv$-complete interferometer array
with the pyuvsim package (Lanman et al., 2019). Using the average wavelength
at 1.64 m as a reference, we first evenly divide half of the $uv$ space from
10 to 100 wavelengths into one-wavelength cells. We then fill each cell with a
baseline at a randomized location, giving us 15,604 visibilities in total.
This setup represents an array with complete $uv$ coverage from ten to one
hundred wavelengths at the resolution of one wavelength — 1.64 m. The $uv$
pattern is created at the reference wavelength and fixed in its physical size,
so it evolves with frequency: both the $uv$ coverage and resolution change
across the frequency band. This $uv$ pattern does not require any certain
pairs of antennas located at one reference wavelength apart, because the $uv$
coverage starts from 10 reference wavelengths. For example, antennas
representing two close baselines can come from four antennas — two pairs —
with each pair separated by at least 10 reference wavelengths; an array
designed to cover these $uv$ locations will have additional baselines beyond
the ones in the simulation, we do not consider them for this study. Adding the
additional baseline, equivalently adding additional information, will only
improve the results from this study.

After adopting a 14-m Airy profile for the primary beam, we use a diffuse sky
map (Byrne et al., 2022) and the GLEAM sources (Wayth et al., 2015; Hurley-
Walker et al., 2017, 2019) to simulate the visibilities for one time-
integration. During the integration, the array points at 0∘ R.A. and -26.7∘
Dec. No noise or mitigated systematics is added to the visibilities.

We use direct optimal mapping to map the visibilities to one image cube; each
sky image has the same range and resolution as the image cube in Section 3.1
except centering at the 0∘ R.A. and -26.7∘ Dec. We select 180 channels closest
to the frequency range in the previous analysis — 167.07 MHz to 181.40 MHz
with a resolution of 0.08 MHz. Please note that the $uv$-complete simulation
uses a different sky model, assumes a different primary beam model, centers on
a different sky patch, and measures a different frequency range compared to
the previous section; so the exact foreground map and power are different.
However, we can still investigate how the foreground power is distributed in
the 2D $k_{\bot}-k_{\parallel}$ space.

After the image cube is formed, we use the image power spectrum pipeline to
calculate its window functions and 2D power spectrum. Figure 6 shows window
functions at the three representative points as in Figure 2. Window functions
with complete $uv$-coverage show less $k_{\parallel}$-leakage from the ones in
Figure 2. However, we still see leakage along the $k_{\parallel}$ at the
$10^{-12}$ level. We find the $k_{\parallel}$ leakage disappears in the
current color scale if we use a more aggressive tapering function — 7-term
Blackman-Harris — along the LOS direction; in exchange, we get a larger kernel
along $k_{\parallel}$.

Figure 7 shows the 2D power spectrum similar to the ones in Figure 3. The
$uv$-complete 2D power spectrum extends farther out along the $k_{\bot}$
direction, thanks to the array’s larger $uv$ range; for comparison, the
extended part is shaded in Figure 7. The horizon wedge is generally more
prominent in high-$k_{\bot}$, so extending $k_{\bot}$ should reveal more of
the horizon wedge. We start to see the horizon wedge in Figure 3 with the HERA
array; however, even after extending to higher $k_{\bot}$ values, we do not
see the wedge in Figure 7 with the $uv$-complete array. This delivers an
interesting message: image power spectrum from a $uv$-complete array does not
see the horizon wedge.

This study decouples the effect from a $uv$-incomplete array and the intrinsic
properties from the image power spectrum. For an array with complete $uv$
coverage, we achieve the full potential of the image power spectrum, with the
2D power spectrum showing a mitigated horizon wedge, which is our initial
motivation to develop the image power spectrum; here we quantitatively
verified this conclusion with end-to-end simulation and analysis. Now, two
caveats are worth mentioning: we assumed a 14-m Airy beam in this study,
without considering the complexity of realistic primary beams (Kim et al.,
2022, 2023); we did not simulate noise, instrument systematics, or calibration
errors, which can introduce the wedge even in image power spectra (Morales et
al., 2019).

Murray & Trott (2018) conducted an investigation of different array layouts,
and concluded that radially randomized arrays show the wedge feature in the 2D
power spectrum. The power spectrum estimator in Murray & Trott (2018) is based
on $uv$-mapmaking, different from the image power spectrum in this paper.

##  5 Future Work

We have developed the image power spectrum pipeline and understand its window
function and power spectra with noiseless simulations. We will investigate the
effects of realistic noise, radio-frequency interference, realistic errors in
beam models (Kim et al., 2022, 2023), calibration accuracy, realistic antenna
position, and beam variations (Orosz et al., 2019; Aguirre et al., 2022). We
plan to apply this method on the HERA Phase I data (HERA Collaboration et al.,
2022, 2023). The image power spectrum provides an independent statistic of the
sky in addition to the delay power spectrum (HERA Collaboration et al., 2022,
2023).

Currently, we map only 512 pixels on the sky plane, limited by the size of the
$\mathbf{P}$ matrix and the calculation of the window function. After
understanding the window function, we do not need to calculate its full form;
instead, we only need the sum of each $\mathbf{k}$ bin for normalization,
which will dramatically reduce the RAM requirement of calculation. The
reduction will enable us to increase the size and resolution of the image
cube, enlarging the $k_{\bot}$ range. With an enlarged sky patch, we can
include more data at different local sidereal times (LSTs); we leave these
topics for future work.

##  6 Conclusion

The image power spectrum measures $k_{\parallel}$ directly through sky pixels,
mitigating spatial power leaking into the frequency axis through measuring
visibilities. This paper presents an FFT-based image power spectrum based on
the direct optimal mapping (X22).

We use noiseless simulation data from Aguirre et al. (2022) to explore the
image power spectrum. After obtaining the image cube with direct optimal
mapping, the image cube is tapered along the three axes before performing a 3D
FFT to calculate the power spectrum. We calculate the normalization factors
for the power spectrum and the power spectrum window functions. The window
functions show $<10^{-11}$ leakage along the $k_{\parallel}$ direction,
separating the foreground from EoR signals. The 2D and 1D power spectra also
show that the EoR and foreground signals are clearly separated in $\mathbf{k}$
space. However, we still see the horizon wedge in the 2D power spectrum.

To understand whether the wedge feature comes from the image power spectrum or
from HERA’s incomplete $uv$ coverage, we simulated visibilities from a
$uv$-complete array and analyzed the visibilities with our image power
spectrum pipeline. For the window functions, leakage along the $k_{\parallel}$
direction is less than $10^{-12}$, currently dominated by the tapering
function’s dynamic range; the 2D power spectrum does not show wedge features
either, even after extending the $k_{\bot}$ range. This study highlights the
potential utility of complete $uv$-coverage when evaluating future array
designs: a $uv$-complete array, which contains complete information for
imaging, does not show the horizon wedge in 2D image power spectra.

We plan to measure the image power spectrum using HERA data (DeBoer et al.,
2017; HERA Collaboration et al., 2022, 2023). The result will provide a
measurement that complements the delay spectrum, with different foreground
features in the 2D power spectrum.

## Acknowledgements

This analysis utilized custom-built, publicly-accessible software by the HERA
Collaboration (https://github.com/Hera-Team) in addition to software built by
both HERA members and collaborators
(https://github.com/RadioAstronomySoftwareGroup), especially pyuvdata
(Hazelton et al., 2017). This analysis also relied on number of public, open-
source software packages, including numpy (Harris et al., 2020), scipy
(Virtanen et al., 2020), matplotlib (Hunter, 2007), and astropy (Astropy
Collaboration et al., 2018).

This material is based upon work supported by the National Science Foundation
under grants #1636646 and #1836019 and institutional support from the HERA
collaboration partners. This research is funded in part by the Gordon and
Betty Moore Foundation through Grant GBMF5212 to the Massachusetts Institute
of Technology. HERA is hosted by the South African Radio Astronomy
Observatory, which is a facility of the National Research Foundation, an
agency of the Department of Science and Innovation.

We acknowledge the use of the Ilifu cloud computing facility (www.ilifu.ac.za)
and the support from the Inter-University Institute for Data Intensive
Astronomy (IDIA; https://www.idia.ac.za).

J.S. Dillon gratefully acknowledges the support of the NSF AAPF award
#1701536. NK acknowledges support from NASA through the NASA Hubble Fellowship
grant #HST-HF2-51533.001-A awarded by the Space Telescope Science Institute,
which is operated by the Association of Universities for Research in
Astronomy, Incorporated, under NASA contract NAS5-26555. R. Byrne is supported
by National Science Foundation Award No. 2303952. P. Kittiwisit and M.G.
Santos acknowledge support from the South African Radio Astronomy Observatory
(SARAO; www.sarao.ac.za) and the National Research Foundation (Grant No.
84156). This result is part of a project that has received funding from the
European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation programme (Grant agreement No. 948764; P. Bull and
M.J. Wilensky). P. Bull and H. Garsden acknowledge support from STFC Grant
ST/T000341/1. Parts of this research were supported by the Australian Research
Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO
3D), through project number CE170100013. G. Bernardi acknowledges funding from
the INAF PRIN-SKA 2017 project 1.05.01.88.04 (FORECaST), support from the
Ministero degli Affari Esteri della Cooperazione Internazionale – Direzione
Generale per la Promozione del Sistema Paese Progetto di Grande Rilevanza
ZA18GR02 and the National Research Foundation of South Africa (Grant Number
113121) as part of the ISARP RADIOSKY2020 Joint Research Scheme, from the
Royal Society and the Newton Fund under grant NA150184 and from the National
Research Foundation of South Africa (grant No. 103424). E. de Lera Acedo
acknowledges the funding support of the UKRI Science and Technology Facilities
Council SKA grant. A. Liu acknowledges support from the New Frontiers in
Research Fund Exploration grant program, the Canadian Institute for Advanced
Research (CIFAR) Azrieli Global Scholars program, a Natural Sciences and
Engineering Research Council of Canada (NSERC) Discovery Grant and a Discovery
Launch Supplement, the Sloan Research Fellowship, and the William Dawson
Scholarship at McGill.

## References

  * Abbott et al. (2022) Abbott, T. M. C., Aguena, M., Alarcon, A., et al. 2022, Phys. Rev. D, 105, 023520, doi: 10.1103/PhysRevD.105.023520
  * Aguirre et al. (2022) Aguirre, J. E., Murray, S. G., Pascua, R., et al. 2022, ApJ, 924, 85, doi: 10.3847/1538-4357/ac32cd
  * Alam et al. (2017) Alam, S., Ata, M., Bailey, S., et al. 2017, MNRAS, 470, 2617, doi: 10.1093/mnras/stx721
  * Ali et al. (2015) Ali, Z. S., Parsons, A. R., Zheng, H., et al. 2015, ApJ, 809, 61, doi: 10.1088/0004-637X/809/1/61
  * Anderson et al. (2014) Anderson, L., Aubourg, É., Bailey, S., et al. 2014, MNRAS, 441, 24, doi: 10.1093/mnras/stu523
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., SipHocz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Barry et al. (2019) Barry, N., Beardsley, A. P., Byrne, R., et al. 2019, PASA, 36, e026, doi: 10.1017/pasa.2019.21
  * Beardsley et al. (2016) Beardsley, A. P., Hazelton, B. J., Sullivan, I. S., et al. 2016, ApJ, 833, 102, doi: 10.3847/1538-4357/833/1/102
  * Bennett et al. (1996) Bennett, C. L., Banday, A. J., Gorski, K. M., et al. 1996, ApJ, 464, L1, doi: 10.1086/310075
  * Bennett et al. (2013) Bennett, C. L., Larson, D., Weiland, J. L., et al. 2013, ApJS, 208, 20, doi: 10.1088/0067-0049/208/2/20
  * Byrne et al. (2022) Byrne, R., Morales, M. F., Hazelton, B., et al. 2022, MNRAS, 510, 2011, doi: 10.1093/mnras/stab3276
  * Clark (1980) Clark, B. G. 1980, A&A, 89, 377
  * Cornwell (2008) Cornwell, T. J. 2008, IEEE Journal of Selected Topics in Signal Processing, 2, 793, doi: 10.1109/JSTSP.2008.2006388
  * Datta et al. (2010) Datta, A., Bowman, J. D., & Carilli, C. L. 2010, ApJ, 724, 526, doi: 10.1088/0004-637X/724/1/526
  * DeBoer et al. (2017) DeBoer, D. R., Parsons, A. R., Aguirre, J. E., et al. 2017, PASP, 129, 045001, doi: 10.1088/1538-3873/129/974/045001
  * Dillon et al. (2013) Dillon, J. S., Liu, A., & Tegmark, M. 2013, Phys. Rev. D, 87, 043005, doi: 10.1103/PhysRevD.87.043005
  * Dillon & Parsons (2016) Dillon, J. S., & Parsons, A. R. 2016, ApJ, 826, 181, doi: 10.3847/0004-637X/826/2/181
  * Dillon et al. (2014) Dillon, J. S., Liu, A., Williams, C. L., et al. 2014, Phys. Rev. D, 89, 023002, doi: 10.1103/PhysRevD.89.023002
  * Dillon et al. (2015) Dillon, J. S., Neben, A. R., Hewitt, J. N., et al. 2015, Phys. Rev. D, 91, 123011, doi: 10.1103/PhysRevD.91.123011
  * Ewall-Wice et al. (2016) Ewall-Wice, A., Dillon, J. S., Hewitt, J. N., et al. 2016, MNRAS, 460, 4320, doi: 10.1093/mnras/stw1022
  * Fixsen et al. (1996) Fixsen, D. J., Cheng, E. S., Gales, J. M., et al. 1996, ApJ, 473, 576, doi: 10.1086/178173
  * Furlanetto et al. (2006) Furlanetto, S. R., Oh, S. P., & Briggs, F. H. 2006, Phys. Rep., 433, 181, doi: 10.1016/j.physrep.2006.08.002
  * Gorce et al. (2023) Gorce, A., Ganjam, S., Liu, A., et al. 2023, MNRAS, 520, 375, doi: 10.1093/mnras/stad090
  * Greig et al. (2016) Greig, B., Mesinger, A., & Pober, J. C. 2016, MNRAS, 455, 4295, doi: 10.1093/mnras/stv2618
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357–362, doi: 10.1038/s41586-020-2649-2
  * Hazelton et al. (2017) Hazelton, B. J., Jacobs, D. C., Pober, J. C., & Beardsley, A. P. 2017, The Journal of Open Source Software, 2, 140, doi: 10.21105/joss.00140
  * HERA Collaboration et al. (2022) HERA Collaboration, Abdurashidova, Z., Aguirre, J. E., et al. 2022, ApJ, 925, 221, doi: 10.3847/1538-4357/ac1c78
  * HERA Collaboration et al. (2023) HERA Collaboration, Abdurashidova, Z., Adams, T., et al. 2023, The Astrophysical Journal, 945, 124, doi: 10.3847/1538-4357/acaf50
  * Hinshaw et al. (2013) Hinshaw, G., Larson, D., Komatsu, E., et al. 2013, ApJS, 208, 19, doi: 10.1088/0067-0049/208/2/19
  * Hogan & Rees (1979) Hogan, C. J., & Rees, M. J. 1979, MNRAS, 188, 791, doi: 10.1093/mnras/188.4.791
  * Högbom (1974) Högbom, J. A. 1974, A&AS, 15, 417
  * Hu & Dodelson (2002) Hu, W., & Dodelson, S. 2002, ARA&A, 40, 171, doi: 10.1146/annurev.astro.40.060401.093926
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Hurley-Walker et al. (2017) Hurley-Walker, N., Callingham, J. R., Hancock, P. J., et al. 2017, MNRAS, 464, 1146, doi: 10.1093/mnras/stw2337
  * Hurley-Walker et al. (2019) Hurley-Walker, N., Hancock, P. J., Franzen, T. M. O., et al. 2019, PASA, 36, e047, doi: 10.1017/pasa.2019.37
  * Kern & Liu (2021) Kern, N. S., & Liu, A. 2021, MNRAS, 501, 1463, doi: 10.1093/mnras/staa3736
  * Kern et al. (2017) Kern, N. S., Liu, A., Parsons, A. R., Mesinger, A., & Greig, B. 2017, ApJ, 848, 23, doi: 10.3847/1538-4357/aa8bb4
  * Kim et al. (2022) Kim, H., Nhan, B. D., Hewitt, J. N., et al. 2022, ApJ, 941, 207, doi: 10.3847/1538-4357/ac9eaf
  * Kim et al. (2023) Kim, H., Kern, N. S., Hewitt, J. N., et al. 2023, ApJ, 953, 136, doi: 10.3847/1538-4357/ace35e
  * Lanman et al. (2019) Lanman, A., Hazelton, B., Jacobs, D., et al. 2019, The Journal of Open Source Software, 4, 1234, doi: 10.21105/joss.01234
  * Li et al. (2019) Li, W., Pober, J. C., Barry, N., et al. 2019, ApJ, 887, 141, doi: 10.3847/1538-4357/ab55e4
  * Liu et al. (2014) Liu, A., Parsons, A. R., & Trott, C. M. 2014, Phys. Rev. D, 90, 023018, doi: 10.1103/PhysRevD.90.023018
  * Liu et al. (2014) Liu, A., Parsons, A. R., & Trott, C. M. 2014, Phys. Rev. D, 90, 023019, doi: 10.1103/PhysRevD.90.023019
  * Liu & Shaw (2020) Liu, A., & Shaw, J. R. 2020, PASP, 132, 062001, doi: 10.1088/1538-3873/ab5bfd
  * Liu & Tegmark (2011) Liu, A., & Tegmark, M. 2011, Phys. Rev. D, 83, 103006, doi: 10.1103/PhysRevD.83.103006
  * Madau et al. (1997) Madau, P., Meiksin, A., & Rees, M. J. 1997, ApJ, 475, 429, doi: 10.1086/303549
  * Mao et al. (2008) Mao, Y., Tegmark, M., McQuinn, M., Zaldarriaga, M., & Zahn, O. 2008, Phys. Rev. D, 78, 023529, doi: 10.1103/PhysRevD.78.023529
  * Mesinger & Furlanetto (2007) Mesinger, A., & Furlanetto, S. 2007, ApJ, 669, 663, doi: 10.1086/521806
  * Morales et al. (2019) Morales, M. F., Beardsley, A., Pober, J., et al. 2019, MNRAS, 483, 2207, doi: 10.1093/mnras/sty2844
  * Morales et al. (2012) Morales, M. F., Hazelton, B., Sullivan, I., & Beardsley, A. 2012, ApJ, 752, 137, doi: 10.1088/0004-637X/752/2/137
  * Morales & Hewitt (2004) Morales, M. F., & Hewitt, J. 2004, ApJ, 615, 7, doi: 10.1086/424437
  * Morales & Wyithe (2010) Morales, M. F., & Wyithe, J. S. B. 2010, ARA&A, 48, 127, doi: 10.1146/annurev-astro-081309-130936
  * More et al. (2023) More, S., Sugiyama, S., Miyatake, H., et al. 2023, arXiv e-prints, arXiv:2304.00703, doi: 10.48550/arXiv.2304.00703
  * Murray & Trott (2018) Murray, S. G., & Trott, C. M. 2018, ApJ, 869, 25, doi: 10.3847/1538-4357/aaebfa
  * Orosz et al. (2019) Orosz, N., Dillon, J. S., Ewall-Wice, A., Parsons, A. R., & Thyagarajan, N. 2019, MNRAS, 487, 537, doi: 10.1093/mnras/stz1287
  * Parsons et al. (2012) Parsons, A., Pober, J., McQuinn, M., Jacobs, D., & Aguirre, J. 2012, The Astrophysical Journal, 753, 81, doi: 10.1088/0004-637X/753/1/81
  * Parsons et al. (2012) Parsons, A. R., Pober, J. C., Aguirre, J. E., et al. 2012, ApJ, 756, 165, doi: 10.1088/0004-637X/756/2/165
  * Parsons et al. (2010) Parsons, A. R., Backer, D. C., Foster, G. S., et al. 2010, AJ, 139, 1468, doi: 10.1088/0004-6256/139/4/1468
  * Patil et al. (2017) Patil, A. H., Yatawatta, S., Koopmans, L. V. E., et al. 2017, ApJ, 838, 65, doi: 10.3847/1538-4357/aa63e7
  * Planck Collaboration et al. (2016) Planck Collaboration, Adam, R., Ade, P. A. R., et al. 2016, A&A, 594, A1, doi: 10.1051/0004-6361/201527101
  * Planck Collaboration et al. (2020) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A1, doi: 10.1051/0004-6361/201833880
  * Pober et al. (2014) Pober, J. C., Liu, A., Dillon, J. S., et al. 2014, ApJ, 782, 66, doi: 10.1088/0004-637X/782/2/66
  * Pritchard & Loeb (2012) Pritchard, J. R., & Loeb, A. 2012, Reports on Progress in Physics, 75, 086901, doi: 10.1088/0034-4885/75/8/086901
  * Rahimi et al. (2021) Rahimi, M., Pindor, B., Line, J. L. B., et al. 2021, MNRAS, 508, 5954, doi: 10.1093/mnras/stab2918
  * Rau & Cornwell (2011) Rau, U., & Cornwell, T. J. 2011, A&A, 532, A71, doi: 10.1051/0004-6361/201117104
  * Riess et al. (1998) Riess, A. G., Filippenko, A. V., Challis, P., et al. 1998, AJ, 116, 1009, doi: 10.1086/300499
  * Seljak (1998) Seljak, U. 1998, ApJ, 503, 492, doi: 10.1086/306019
  * Sobacchi & Mesinger (2014) Sobacchi, E., & Mesinger, A. 2014, MNRAS, 440, 1662, doi: 10.1093/mnras/stu377
  * Staggs et al. (2018) Staggs, S., Dunkley, J., & Page, L. 2018, Reports on Progress in Physics, 81, 044901, doi: 10.1088/1361-6633/aa94d5
  * Sullivan et al. (2012) Sullivan, I. S., Morales, M. F., Hazelton, B. J., et al. 2012, ApJ, 759, 17, doi: 10.1088/0004-637X/759/1/17
  * Tegmark (1997) Tegmark, M. 1997, Phys. Rev. D, 55, 5895, doi: 10.1103/PhysRevD.55.5895
  * Tegmark et al. (2002) Tegmark, M., Hamilton, A. J. S., & Xu, Y. 2002, MNRAS, 335, 887, doi: 10.1046/j.1365-8711.2002.05622.x
  * Tingay et al. (2013) Tingay, S. J., Goeke, R., Bowman, J. D., et al. 2013, PASA, 30, e007, doi: 10.1017/pasa.2012.007
  * Trott et al. (2016) Trott, C. M., Pindor, B., Procopio, P., et al. 2016, ApJ, 818, 139, doi: 10.3847/0004-637X/818/2/139
  * van Haarlem et al. (2013) van Haarlem, M. P., Wise, M. W., Gunst, A. W., et al. 2013, A&A, 556, A2, doi: 10.1051/0004-6361/201220873
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Wayth et al. (2015) Wayth, R. B., Lenc, E., Bell, M. E., et al. 2015, PASA, 32, e025, doi: 10.1017/pasa.2015.26
  * Xu et al. (2022) Xu, Z., Hewitt, J. N., Chen, K.-F., et al. 2022, ApJ, 938, 128, doi: 10.3847/1538-4357/ac9053