# A Detection of Cosmological 21 cm Emission from CHIME in Cross-correlation
with eBOSS Measurements of the Lyman-$\alpha$ Forest

Mandana Amiri Department of Physics and Astronomy, University of British
Columbia, Vancouver, BC, Canada  Kevin Bandura Department of Computer Science
and Electrical Engineering, West Virginia University, Morgantown WV, USA
Center for Gravitational Waves and Cosmology, West Virginia University,
Morgantown, WV, USA  Arnab Chakraborty Department of Physics, McGill
University, Montreal, QC, Canada  Trottier Space Institute, McGill University,
Montreal, QC, Canada  Matt Dobbs Department of Physics, McGill University,
Montreal, QC, Canada  Trottier Space Institute, McGill University, Montreal,
QC, Canada  Mateus Fandino Department of Physical Sciences, Thompson Rivers
University, Kamloops, BC, Canada  Simon Foreman Department of Physics, Arizona
State University, Tempe, AZ, USA  Hyoyin Gan Department of Physics and
Astronomy, University of British Columbia, Vancouver, BC, Canada  Mark Halpern
Department of Physics and Astronomy, University of British Columbia,
Vancouver, BC, Canada  Alex S. Hill Department of Computer Science, Math,
Physics, and Statistics, University of British Columbia-Okanagan, Kelowna, BC,
Canada  Dominion Radio Astrophysical Observatory, Herzberg Astronomy &
Astrophysics Research Centre, National Research Council Canada, Penticton, BC,
Canada  Gary Hinshaw Department of Physics and Astronomy, University of
British Columbia, Vancouver, BC, Canada  Carolin Höfer Department of Physics
and Astronomy, University of British Columbia, Vancouver, BC, Canada  Kapteyn
Astronomical Institute, University of Groningen, Groningen, The Netherlands
T.L. Landecker Dominion Radio Astrophysical Observatory, Herzberg Astronomy &
Astrophysics Research Centre, National Research Council Canada, Penticton, BC,
Canada  Zack Li Canadian Institute for Theoretical Astrophysics, Toronto, ON,
Canada  Joshua MacEachern Department of Physics and Astronomy, University of
British Columbia, Vancouver, BC, Canada  Kiyoshi Masui MIT Kavli Institute for
Astrophysics and Space Research, Massachusetts Institute of Technology,
Cambridge, MA, USA  Department of Physics, Massachusetts Institute of
Technology, Cambridge, MA, USA  Juan Mena-Parra Dunlap Institute for Astronomy
and Astrophysics, University of Toronto, Toronto, ON, Canada  David A. Dunlap
Department of Astronomy & Astrophysics, University of Toronto, Toronto, ON,
Canada  Nikola Milutinovic Department of Physics and Astronomy, University of
British Columbia, Vancouver, BC, Canada  Arash Mirhosseini Department of
Physics and Astronomy, University of British Columbia, Vancouver, BC, Canada
Laura Newburgh Department of Physics, Yale University, New Haven, CT,
USAUniversity  Anna Ordog Department of Computer Science, Math, Physics, and
Statistics, University of British Columbia-Okanagan, Kelowna, BC, Canada
Dominion Radio Astrophysical Observatory, Herzberg Astronomy & Astrophysics
Research Centre, National Research Council Canada, Penticton, BC, Canada
Sourabh Paul Department of Physics, McGill University, Montreal, QC, Canada
Trottier Space Institute, McGill University, Montreal, QC, Canada  Ue-Li Pen
Canadian Institute for Theoretical Astrophysics, Toronto, ON, Canada  Tristan
Pinsonneault-Marotte Department of Physics and Astronomy, University of
British Columbia, Vancouver, BC, Canada  Alex Reda Department of Physics, Yale
University, New Haven, CT, USAUniversity  J. Richard Shaw Department of
Physics and Astronomy, University of British Columbia, Vancouver, BC, Canada
Seth R. Siegel Department of Physics, McGill University, Montreal, QC, Canada
Trottier Space Institute, McGill University, Montreal, QC, Canada  Perimeter
Institute for Theoretical Physics, Waterloo, ON, Canada  Keith Vanderlinde
David A. Dunlap Department of Astronomy & Astrophysics, University of Toronto,
Toronto, ON, Canada  Dunlap Institute for Astronomy and Astrophysics,
University of Toronto, Toronto, ON, Canada  Haochen Wang MIT Kavli Institute
for Astrophysics and Space Research, Massachusetts Institute of Technology,
Cambridge, MA, USA  Department of Physics, Massachusetts Institute of
Technology, Cambridge, MA, USA  D. V. Wiebe Department of Physics and
Astronomy, University of British Columbia, Vancouver, BC, Canada  Dallas Wulf
Department of Physics, McGill University, Montreal, QC, Canada  Trottier Space
Institute, McGill University, Montreal, QC, Canada T. Pinsonneault-Marotte
<EMAIL_ADDRESS>
###### Abstract

We report the detection of 21 cm emission at an average redshift $\bar{z}=2.3$
in the cross-correlation of data from the Canadian Hydrogen Intensity Mapping
Experiment (CHIME) with measurements of the Lyman-$\alpha$ forest from eBOSS.
Data collected by CHIME over 88 days in the $400-500$ MHz frequency band
($1.8<z<2.5$) are formed into maps of the sky and high-pass delay filtered to
suppress the foreground power, corresponding to removing cosmological scales
with $k_{\parallel}\lesssim 0.13\ \text{Mpc}^{-1}$ at the average redshift.
Line-of-sight spectra to the eBOSS background quasar locations are extracted
from the CHIME maps and combined with the Lyman-$\alpha$ forest flux
transmission spectra to estimate the 21 cm-Lyman-$\alpha$ cross-correlation
function. Fitting a simulation-derived template function to this measurement
results in a $9\sigma$ detection significance. The coherent accumulation of
the signal through cross-correlation is sufficient to enable a detection
despite excess variance from foreground residuals $\sim 6-10$ times brighter
than the expected thermal noise level in the correlation function. These
results are the highest-redshift measurement of 21 cm emission to date, and
set the stage for future 21 cm intensity mapping analyses at $z>1.8$.

††software:  bitshuffle (Masui et al., 2015), CAMB (Lewis et al., 2000), caput
(Shaw et al., 2020c), ch_pipeline (Shaw et al., 2020d), cora (Shaw et al.,
2020b), Cython (Behnel et al., 2011), draco (Shaw et al., 2020e), driftscan
(Shaw et al., 2020a), hankl (Karamanis & Beutler, 2021), h5py (Collette et
al., 2021), HDF5 (The HDF Group, 1997), HEALPix (Górski et al., 2005), healpy
(Zonca et al., 2019), Matplotlib (Hunter, 2007), mpi4py (Dalcin & Fang, 2021),
NumPy (Harris et al., 2020), OpenMPI (Gabriel et al., 2004), pandas (pandas
development team, 2020; McKinney, 2010), peewee (Leifer, C., 2021), SciPy
(Virtanen et al., 2020), Skyfield (Rhodes, 2019),

##  1 Introduction

Emission from the hyperfine transition of neutral hydrogen (HI) at 21 cm can
be used to efficiently map the large-scale structure (LSS) of the Universe
over most of its history, a technique known as Hydrogen intensity mapping. In
this approach, HI contained in galaxies or the inter-galactic medium is
detected in aggregate – integrated over the relatively coarse angular
resolution afforded by telescopes at radio wavelengths, where redshifted 21 cm
emission is measured. Such telescopes can be designed to cost-effectively
observe large sky areas at high sensitivity, and digital receivers enable
broad bandwidths to be sampled with fine frequency resolution, which directly
maps to redshift due to the monochromaticity of the signal (Peterson et al.,
2009).

Over the last decade, instruments have been built with the aim of measuring
this signal in the late Universe (e.g. CHIME, MeerKAT) as well as at higher
redshift (HERA, EDGES, SARAS, PAPER, among others). Published results from the
field include upper limits (Abdurashidova et al., 2022), detections in cross-
correlation (Pen et al., 2009; Chang et al., 2010; Masui et al., 2013;
Anderson et al., 2018; Tramonte & Ma, 2020; Li et al., 2021; Wolz et al.,
2022; CHIME Collaboration et al., 2022a; Cunnington et al., 2023), and a first
detection in auto-correlation (Paul et al., 2023). The main challenge for all
of these efforts is separating the faint 21 cm signal from the extremely
bright foreground emission of our galaxy and radio point sources, which can in
principle be achieved by making use of their different spectral properties,
but is complicated by instrumental effects and radio-frequency interference
(RFI). Until foreground separation methods mature to a point where the signal
becomes dominant, cross-correlation with external surveys has proven to be an
effective way of mitigating residual foregrounds and obtaining interesting
scientific constraints, as demonstrated by the recently reported detection
from CHIME in cross-correlation with SDSS galaxy surveys (CHIME Collaboration
et al., 2022a).

The Canadian Hydrogen Intensity Mapping Experiment (CHIME) (CHIME
Collaboration et al., 2022b) is a compact interferometer composed of four
cylindrical reflectors instrumented with a total of 1024 dual-polarisation
antennas. It is located at the Dominion Radio Astrophysical Observatory (DRAO)
in Penticton, Canada, where it observes the entire sky visible at latitude
$\sim 49$∘, operating as a driftscan telescope over 1024 channels in the
$400-800$ MHz band (corresponding to 21 cm emission redshifted by
$2.5>z>0.8$). CHIME has been operating continuously since 2017.

The redshift range available to CHIME overlaps with measurements of the
Lyman-$\alpha$ forest from eBOSS (du Mas des Bourboux et al., 2020), allowing
for a cross-correlation analysis to be carried out in the high-redshift end of
the CHIME band. At these redshifts, the quasar catalogue that was stacked on
in CHIME Collaboration et al. (2022a) is sparser, and there is more
statistical weight in the Lyman-$\alpha$ data. Lyman-$\alpha$ forest
measurements are characterized by absorption of quasar light by HI clouds
along the line of sight whereas the 21 cm signal is in emission, so we expect
the cross-correlation with CHIME maps to be negative at small separations,
which provides a unique signal that is difficult to mimic any other way. This
is a feature that has been highlighted by Carucci et al. (2017), who also
emphasize the usefulness of such a cross-correlation to mitigate foreground
contamination and break degeneracies in a power spectrum analysis.

Several studies have found evidence that various classes of galaxies are
correlated with Lyman-$\alpha$ absorption by the IGM up to several tens of
megaparsecs (Mukae et al., 2017; Momose et al., 2021; Liang et al., 2021), but
these studies have been limited to small numbers of galaxies and/or absorption
systems. A joint analysis of CHIME and eBOSS Lyman-$\alpha$ forest data will
provide a huge statistical improvement on both sides, opening a new window on
the connection between HI-rich galaxies and low-density HI in the surrounding
IGM, and thus shedding light on the role of HI in galaxy evolution (e.g. the
strength of HI inflows from the IGM onto denser systems) when the cosmic star
formation rate is at its peak.

In this work, we present a detection of 21 cm emission in the cross-
correlation of CHIME data and eBOSS Lyman-$\alpha$ forest measurements, but
leave the interpretation of the signal (including modelling of its
cosmological and astrophysical implications) for future work. We describe the
data used in this analysis in Section 2 and the CHIME processing pipeline in
Section 3. Section 4 explains the cross-correlation method we use. A model for
the signal and its inclusion in simulations are described in Section 5. The
results of the analysis are presented in Section 6, and some validation tests
in Section 7. Finally, we conclude in Section 8.

##  2 Data

###  2.1 eBOSS

The Lyman-$\alpha$ forest measurements we use in our analysis are those from
the eBOSS sixteenth data release (SDSS DR16) (du Mas des Bourboux et al.,
2020). The data products consist of fractional flux transmission along the
line of sight to backlight quasars, defined as

$\delta_{q}(\lambda)=\frac{f_{q}(\lambda)}{\bar{F}(\lambda)C_{q}(\lambda)}-1\text{,}$
(1)

where $f_{q}$ is the measured flux, $C_{q}$ is the unabsorbed quasar continuum
and $\bar{F}$ is the mean transmission. The product
$\bar{F}(\lambda)C_{q}(\lambda)$ was modelled and fit to the data as explained
in du Mas des Bourboux et al. (2020). Apart from interpolating onto the CHIME
frequency band (Section 4.1), no additional processing is performed on the
catalogue provided by eBOSS.

Figure 1: Histogram of redshifts of eBOSS Lyman-$\alpha$ forest samples and
their overlap with the redshift-range probed by CHIME.

There are 210005 background quasars in the eBOSS Lyman-$\alpha$ catalogue,
with redshifts $2.1<z<4$. The CHIME data go down to 400 MHz, corresponding to
21 cm radiation redshifted at $z=2.5$, which overlaps with a total of $\sim
3\times 10^{7}$ spectral samples in the eBOSS measurements (see Figure 1). We
include the full eBOSS spectra in this total and the analysis, including the
regions blueward of the Lyman-$\beta$ wavelength ($\sim 102$ nm) at the quasar
redshift (referred to as the Ly$\beta$ region in the eBOSS papers). This part
of the spectrum is treated separately in their analysis because it can include
absorption from both transition levels, but since we are cross-correlating
with an external dataset, there is no concern of contamination from Ly$\beta$
correlations. Damped Lyman-$\alpha$ systems are identified and masked in the
eBOSS data products, so they do not contribute to the cross-correlation with
CHIME. The quasars are distributed throughout the SDSS North and South
Galactic Cap fields (NGC and SGC), which are fully contained within the sky
area observed by CHIME.

###  2.2 CHIME

The CHIME data that was used in this work is composed of nighttime
observations from 88 days spanning the calendar year 2019. This amounts to
$\sim 1000$ hours of total observation time, distributed approximately
uniformly over the sidereal day. The full frequency range observed with CHIME
is $400-800$ MHz but we use only the lowest quarter of the band, i.e.
$400-500$ MHz or $z>1.8$, that has overlap with the Lyman-$\alpha$ forest
measurements. This is a similar set of observations to what was used in CHIME
Collaboration et al. (2022a) to detect 21 cm emission in cross-correlation
with SDSS galaxy catalogues (refer to that paper and citations therein for a
general description of the CHIME data). However, that analysis considered data
from the upper part of the frequency band, $585-800$ MHz, whereas this work is
the first report on the high-redshift end.

##  3 Processing

The processing that the CHIME data undergoes prior to cross-correlation with
the Lyman-$\alpha$ forest is largely the same as described in the first
detection paper (Section 3 of CHIME Collaboration et al. 2022a). In this
section we summarise the main steps and highlight improvements that were made
to the pipeline since the publication of that analysis as well as aspects of
the processing that had to be modified to accommodate the lower frequency
band.

###  3.1 Sidereal Stacking Pipeline

The final product of the CHIME processing pipeline is a set of visibilities
for unique baselines measuring the sky on a grid of right ascension spanning
the sidereal day. To produce this “sidereal stack”, observations from
individual days are collected and passed through daily processing that
includes RFI flagging, gain calibration and regridding from local time to a
fixed grid in sidereal time. Each day is inspected visually via a set of data
quality metrics to identify and flag days that appear corrupted. The remaining
days are averaged together to form the final sidereal stack that is used in
this analysis. Further details on each stage of this processing can be found
in CHIME Collaboration et al. (2022a).

####  3.1.1 Improved Flagging

Since the publication of the first CHIME detection, some improvements to the
processing pipeline have been made. These include:

##### Narrowband gain errors

The calibration algorithms employed in the real-time pipeline attempt to
identify and mask RFI-like features in the underlying data prior to estimating
the frequency- and feed-dependent gains. However, these algorithms are not
entirely robust to RFI and we have found that transient RFI occasionally
biases the resulting gain estimates. This is true for both the “digital
gains”, which are updated semi-annually and whose purpose is to minimize
quantization noise when truncating the data streams to 4 bits (real) + 4 bits
(imaginary), and the “calibration gains”, which are updated daily and whose
purpose is to correct instability in the analog receiver chain. The offending
RFI is narrowband, usually with a bandwidth less than a single 390.625 kHz
frequency channel, and results in spikes in the estimated gain as a function
of frequency. These spikes are imprinted on the foregrounds when the gains are
applied to the visibility data, leaking foreground power to high delays (the
conjugate axis to frequency, where smooth foregrounds would otherwise be
confined to low delays).

We have developed a new method to search the archived gains for these
narrowband artifacts and then mask the corresponding frequencies and times in
the visibility data. For each time sample in the visibility data, we calculate
the product of the digital gain and calibration gain that was applied to each
feed and frequency by the real-time pipeline, average the amplitude over
feeds, and then perform a search for narrowband features along the frequency
axis. The search algorithm applies an aggressive, high-pass filter along the
frequency axis, masks the frequency channel that is the largest outlier,
reapplies the filter accounting for the newly masked data, and iterates until
all unmasked frequencies are less than $6\sigma$, where $\sigma$ is an
estimate of the standard deviation of the noise. This procedure masked 6% of
the data that was considered for this analysis.

##### Decorrelation events

On rare occasions, communication errors during the corner-turn operation can
cause certain data streams in the correlator to become misaligned with all
other data streams. This results in one quarter of the feeds decorrelating
with all other feeds for 64 of the 1024 frequency channels, which persists
until the data streams are re-aligned by restarting the correlator. The 64
decorrelated frequency channels are uniformly spaced across the band and leak
significant foreground power to high delays. We now perform an automated
search for these decorrelation events and exclude any frequency channel and
time that is affected from further analysis. In total, 0.07% of the data
considered for this analysis was masked for this reason, corresponding to a
single decorrelation event that persisted for the majority of one sidereal
day.

##### Excessively small weights

The real-time RFI flagging excises corrupt samples within the time
integrations that are eventually saved to disk; for integrations where this
excision is nearly total, the data saved to disk can be extremely noisy. A
similar effect can occur due to packet loss or other digital errors in the
X-engine. These integrations are assigned very low noise weights in subsequent
analysis, but if the issues are confined to very narrow frequency bands or
small subsets of baselines, the low weights were not properly incorporated in
previous versions of the processing. This has been corrected using an extra
flagging step in the offline data pipeline, which flags 0.5% of all time-
frequency samples used in this analysis.

####  3.1.2 Thermal Gain Correction

Section 3.2.4 of CHIME Collaboration et al. (2022a) describes a stage of the
data processing that uses measurements of the outside temperature to correct
common-mode thermal variations in the amplitude of the gains of the analog
receiver chains. We recently discovered two errors in the resulting instrument
stability as reported in CHIME Collaboration et al. (2022a). First, the
pre(post)-correction stability of $0.8\%$ ($0.5\%$) was reported as the
standard deviation in fractional power, but it is actually the standard
deviation in fractional voltage. Second, the software that was deployed in the
pipeline mistakenly applied the inverse of the thermal correction, which
amplified the thermal variations instead of mitigating them. We have fixed the
software error and applied the correct factor to the data used for this
analysis. As a result, the estimated stability for the analysis presented in
CHIME Collaboration et al. (2022a) is $2.7\%$ (standard deviation in
fractional power), while the estimated stability for this analysis is $1.0\%$.

####  3.1.3 Impact on Previously-reported Results

The changes to the CHIME processing pipeline described in this section do not
result in an improved signal-to-noise on the detection of cosmological 21 cm
emission in cross-correlation with the eBOSS galaxy catalogues that was
reported in CHIME Collaboration et al. (2022a). In fact, the noise level in
the more recent revision of the processed data is found to be elevated
compared to what was measured in the previous one, but this has no qualitative
impact on the results reported in CHIME Collaboration et al. (2022a).
Investigating the impact of these changes and continued iteration on the
processing pipeline is currently ongoing.

###  3.2 Ringmaps

The way CHIME observes the sky – a primary beam that is a few degrees wide in
the East-West direction but spans the South and North horizons, illuminating a
regularly-spaced grid of feeds as the sky drifts overhead – is amenable to a
particular map-making method, which produces “ringmaps”. At every time sample,
corresponding to a specific right ascension (RA) on meridian, the North-South
(NS) baselines are phased and combined so as to form a grid of synthesized
beams within the instantaneous field of view. As the sky drifts overhead, each
formed beam records its intensity along a ring of constant declination. The
image generated in this way is called a ringmap. An example is displayed in
the top panel of Figure 2. A detailed description of this process is given in
CHIME Collaboration et al. (2022a).

For the cross-correlation analysis, we generate ringmaps from the sidereal
stacks at each polarisation, and find the angular pixel that contains the
position of each quasar in the eBOSS catalog. The spectrum from this pixel is
what will be correlated with the corresponding Lyman-$\alpha$ forest spectrum.

Figure 2:  _(First)_ Ringmap generated from inter-cylinder baselines ($\gtrsim
22$ m East-West distance) at 476.56 MHz for the XX polarisation. _(Second)_
Ringmap generated from visibilities filtered using DAYENU with a 200 ns delay
cut. _(Third)_ The same as the second panel but with an outlier mask derived
from a $6\sigma$ threshold based on the expected noise. _(Fourth)_ A noise
realisation at the level expected in the CHIME data, delay-filtered in the
same way as in the previous panels. In all panels, vertical dashed lines
indicate the region that was used to estimate the delay power spectrum in
Figure 4.

###  3.3 Lower Band Processing

####  3.3.1 Beam Calibration

When ringmaps are generated, a beam model is used to deconvolve the effect of
the main lobe of the primary beam. The model that was used in CHIME
Collaboration et al. (2022a) was derived from data by fitting to the expected
flux of a large number of radio point sources. This model has since been
extended to the lower band considered in this analysis, using the same methods
that were described in detail in CHIME Collaboration et al. (2022a). Slices
along each axis of the beam model are presented in Figure 3, for both
polarisations.

Figure 3: Slices of the 3D beam model, from left to right, at constant hour
angle, declination, and frequency. The top row is the power beam for the X
polarisation, and the bottom row is Y. The beam pattern has been projected
onto sky coordinates HA and declination, which is why it broadens
significantly near the NCP. Note that the declination axis has been extended
beyond 90∘to reflect the sensitivity of the beam to regions of the sky on
either side of the NCP.

####  3.3.2 Frequency Flagging

Although RFI is flagged algorithmically at multiple stages in the real-time
and daily processing, low-level events which are fainter than the foregrounds
are very hard to detect. This contamination becomes apparent after a high-pass
delay filter has been applied to remove most of the foreground power, and is
strong enough to obscure the cosmological signal. We identify affected
frequencies by iteratively filtering and flagging outlier frequencies until
none remain. The flagging was performed by visually inspecting maps generated
from the filtered data at every frequency, where the spatial information helps
identify bad cases. In the end, 47% of the data in the $400-500$ MHz band was
flagged in this way, in addition to persistent RFI bands that are statically
masked, such that a total of 65% of the band was masked in this analysis.

We attributed this contamination to RFI, but instrumental miscalibration such
as the narrowband gain errors described above could have a similar effect. As
the calibration and flagging parts of the daily pipeline continue to improve,
we may be able to recover some of these flagged frequencies.

####  3.3.3 Delay Filter

In this analysis we suppress the foregrounds by high-pass filtering the
$400-500$ MHz band in delay. The filter is implemented using the DAYENU method
to account for the masked regions of the frequency band (Ewall-Wice et al.,
2021). We set the lower edge of the delay pass-band at 200 ns, which was
chosen to encompass the bulk of foreground power at all declinations in the
delay power spectrum, as illustrated in Figure 4. This power spectrum was
evaluated using the Gibbs sampling method described in Appendix A of CHIME
Collaboration et al. (2022a) and estimating the variance across right
ascensions 110-270∘, a region corresponding roughly to the eBOSS NGC field,
where the majority of the Lyman-$\alpha$ quasars are found. Shifting this
threshold by 50 or 100 ns in either direction did not appear to improve the
signal-to-noise, but no further optimisation was attemped. Filtering delays
below 200 ns corresponds to removing large scales with $k_{\parallel}\lesssim
0.13\ \text{Mpc}^{-1}$, at the average redshift $\bar{z}=2.3$ ($\sim 430$
MHz). Figure 2 shows maps at a single frequency and polarisation before and
after the delay filter. Pixels near bright foreground features tend to remain
corrupted after filtering, so we also derive a mask that excises any pixel
above a $6\sigma$ threshold determined by the expected thermal noise. This
threshold is what was used in CHIME Collaboration et al. (2022a) and is
sufficient but necessary to enable a detection of the 21 cm signal.

Figure 4: Delay power spectrum of a CHIME ringmap evaluated across right
ascensions between 110-270∘, at fixed declination. The colour scale represents
power in units of noise power. Vertical dashed lines show the lower bounds of
the DAYENU filter pass-band. Significant contamination exists at all delays.
Some of the horizontal features can be attributed to bright point sources in
the maps, but the vertical lines are of unknown origin.

##  4 Cross-correlation Method

The cross-correlation function of the 21 cm temperature $T_{21}$ and the
Lyman-$\alpha$ forest relative flux transmission $\delta_{\alpha}$ is defined
at a given redshift as

$\xi({\boldsymbol{\mathbf{r}}},{\boldsymbol{\mathbf{r}}}^{\prime})=\left\langle
T_{21}\left({\boldsymbol{\mathbf{r}}}\right)\delta_{\alpha}\left({\boldsymbol{\mathbf{r^{\prime}}}}\right)\right\rangle\text{,}$
(2)

where ${\boldsymbol{\mathbf{r}}}$ and ${\boldsymbol{\mathbf{r}}}^{\prime}$ are
three-dimensional comoving positions and the angle brackets indicate an
ensemble average. Assuming statistical isotropy, the correlation function only
depends on the magnitude of the separation in comoving space:
$\xi({\boldsymbol{\mathbf{r}}},{\boldsymbol{\mathbf{r}}}^{\prime})=\xi(\left\lvert{\boldsymbol{\mathbf{r}}}-{\boldsymbol{\mathbf{r}}}^{\prime}\right\rvert)$.
However, measurements of these fields are located on our past lightcone – in
redshift space and on the sphere – denoted by
$(z,\hat{{\boldsymbol{\mathbf{n}}}})$, so this is where we will evaluate the
correlation function.

In aggregate, the line-of-sight measurements from eBOSS and CHIME contain
three-dimensional information on the two fields, i.e. one could evaluate the
correlations along the redshift ($z$) and angular
($\hat{{\boldsymbol{\mathbf{n}}}}$) coordinates. However, CHIME’s resolving
power is much greater in frequency than it is on the sky, so we expect the
signal-to-noise to be stronger along this axis. Given the significant added
complexity of a 3-d analysis, in this work we only consider correlations along
the line-of-sight.

Here, we will describe the procedure used to compute an estimate of the radial
cross-correlation function in redshift space, corresponding to

$\xi\left(\Delta z\right)=\left\langle\int
dz~{}T_{21}^{i}\left(z\right)\delta_{\alpha}^{i}\left(z+\Delta
z\right)\right\rangle_{i}\text{,}$ (3)

where the angle brackets denote the average over individual line-of-sight
observations labelled by $i$. We describe how redshifts are computed for each
observable in Sec. 4.1, and how the cross-correlation is computed in Sec. 4.2.

###  4.1 Frequency Remapping

The eBOSS Lyman-$\alpha$ forest spectra measurements are arranged in
wavelength bins uniformly spaced in $\ln\lambda$ and spanning $\sim 360-700$
nm, whereas the CHIME data is uniformly spaced in frequency in the $400-800$
MHz band. In order to correlate the two, it will be necessary to map these
onto a common grid in redshift. Given that the resolution of the CHIME data is
coarser by a factor of $\sim 1.2-1.5$ in the redshift range where they
overlap, and that it will be convenient to work with the CHIME data directly,
we choose to map the Lyman-$\alpha$ data onto the CHIME band. The mapping is
defined by matching the frequency of Lyman-$\alpha$ absorption to the observed
frequency of 21 cm emission for the corresponding redshift, a relationship
given by the ratio of their rest wavelengths

$\displaystyle\nu^{\prime}=\frac{\lambda_{\alpha}}{\lambda_{21}}\nu\text{.}$
(4)

We use the reverse Lanczos interpolation method described in Section 3 of
CHIME Collaboration et al. (2022a) to regrid the Lyman-$\alpha$ forest spectra
onto the CHIME band. The regions of the spectra that are empty (redshifts
larger than that of the background quasar, or otherwise masked) are set to
zero in the regridded data.

###  4.2 Correlation Function Estimation

We will write the discrete samples in frequency as
$T_{21}[n]=T_{21}(n\Delta\nu+400\mathrm{MHz})$, with integers $n$ labelling
samples at the CHIME frequency resolution, $\Delta\nu=390.625$ kHz. Both the
Lyman-$\alpha$ and CHIME spectra have noise estimates associated with every
sample so we use these to weight the samples and improve signal-to-noise. With
inverse variance weights, $w=\sigma^{-2}$, we estimate the cross-correlation
function as

$\xi[n]=W^{-1}\sum_{i}^{M}\sum_{j}^{N}w_{21}^{i}[j]\,w_{\alpha}^{i}[j+n]\,T_{21}^{i}[j]\,\delta_{\alpha}^{i}[j+n]\>\text{,}$
(5)

where $i$ labels lines of sight, $n$ indexes separation in frequency, and
$W^{-1}$ is a normalisation term given by

$W=\sum_{i}^{M}\sum_{j}^{N}w_{21}^{i}[j]\,w_{\alpha}^{i}[j+n]\text{.}$ (6)

$M=210005$ is the total number of lines of sight and $N\leq 255$ is the number
of samples along the frequency axis being averaged together at that
separation.

##  5 Signal Model and Simulations

To allow us to interpret the results of this analysis and assess detection
significance, we generate simulations of correlated 21 cm and Lyman-$\alpha$
forest measurements. These are run through the same processing pipeline and
analysis as the data in order to forward model the cross-correlation result.
Although the models described below are parametric, we do not attempt to fit
them to the data. Instead, the purpose of these models is to generate a
reasonably realistic template to compare to our measurement and assess the
S/N.

###  5.1 Large-scale Structure Simulation

A detailed description of the model used to generate maps of large-scale
structure in the CHIME redshift range is given in Section 5 of CHIME
Collaboration et al. (2022a). Very briefly, a realisation of Gaussian
fluctuations on the lightcone is generated at the desired redshifts given a
non-linear matter power spectrum evaluated using the halo model prediction
from Mead et al. (2021).

###  5.2 CHIME Signal Model

The HI linear bias, Fingers of God effect, and 21 cm brightness temperature
are modelled to generate a map of 21 cm temperature. See Section 5 of CHIME
Collaboration et al. (2022a) for the model definitions. In this work we fix
the parameters to the fiducial values defined there. Synthetic CHIME
observations are derived from these maps, and directly substituted for the
real data in the cross-correlation analysis.

###  5.3 Lyman-$\alpha$ Signal Model

In order to model the Lyman-$\alpha$ forest measurements, we need a
prescription for converting the matter density fluctuation field
$\delta_{m}(z,\hat{{\boldsymbol{\mathbf{n}}}})$ generated by our LSS
simulation into an optical depth to Lyman-$\alpha$ photons. A commonly used
and straightforward approach is to model the mildly non-linear baryon field as
log-normally distributed (Bi & Davidsen, 1997) and the Lyman-$\alpha$
absorption by the fluctuating Gunn-Peterson approximation (FGPA) (Farr et al.,
2020).

The log-normal transform is performed so as to preserve the variance of the
original field $\sigma_{m}^{2}$:

$1+\delta_{LN}(z,\hat{{\boldsymbol{\mathbf{n}}}})=\exp\left[\delta_{m}(z,\hat{{\boldsymbol{\mathbf{n}}}})-\frac{\sigma_{m}^{2}(z)}{2}\right]\text{.}$
(7)

This step also ensures that the density is strictly non-negative.

The optical depth $\tau$ is proportional to the number density of HI atoms,
which the FGPA assumes is related to the density by a power law

$\tau(z,\hat{{\boldsymbol{\mathbf{n}}}})=\tau_{0}(z)\left(1+\delta_{LN}(z,\hat{{\boldsymbol{\mathbf{n}}}})\right)^{\alpha}\text{.}$
(8)

We adopt the parameter values

$\tau_{0}(z)=0.3\left(\frac{1+z}{3.4}\right)^{4.5},\quad\alpha=1.6\text{,}$
(9)

quoted in Cieplak & Slosar (2016) and Seljak (2012).

These steps are illustrated in the top three panels of Figure 5 for a
realisation of LSS along a given line of sight.

Rather than use the FGPA, it might be more self-consistent to derive the
optical depth from the HI density that was modelled for the 21 cm temperature
directly. This is expected to be part of a future effort to more carefully
model the Lyman-$\alpha$ signal and perhaps derive constraints on the model
parameters. We emphasize that for the purpose of this work, the simplest
approach was taken. However, we will show that a simple FGPA-based model is
sufficient to describe our measurements at the current signal-to-noise level,
so we leave more detailed modelling to future analyses.

###  5.4 Mock Lyman-$\alpha$ Forest Catalogue

From the Lyman-$\alpha$ fractional flux transmission simulation we extract
individual spectra along the lines of sight at a number of sky locations,
chosen to be those of the quasars in the eBOSS catalogue. Each spectrum is
also masked so that regions of missing data – redshifts larger than the
backlight quasar and instrumental masks – match those in the corresponding
eBOSS data. Gaussian random noise is optionally added to the simulated spectra
at a level consistent with the noise variance recorded alongside the data.
This results in a mock catalogue with exactly the same quasar sample as the
data, and thus the same sensitivity to cross-correlation with CHIME, but with
a synthetic Lyman-$\alpha$ forest signal. The bottom panel of Figure 5 shows a
mock spectrum generated in this way.

Figure 5: Fluctuations along a given line of sight in redshift, at subsequent
steps of the procedure for generating simulated Lyman-$\alpha$ forest spectra.
From top to bottom: the matter fluctuations from the LSS simulation; the
result of the log-normal transform; the Lyman-alpha fractional flux
transmission; the final spectrum masked and with noise added to emulate the
real data.

It must be noted that these simulations neglect a number of effects that would
be necessary to include in order to produce realistic simulations of the eBOSS
spectra (and which _are_ included in more sophisticated simulations, e.g. Farr
et al. (2020)). Among the most significant are

  * •

Redshift-space distorsions.

  * •

Metal and Lyman-$\beta$ line absorption.

  * •

High column density absorbers (which include damped Lyman-$\alpha$ absorber
systems)

  * •

Biases due to the modelling of the quasar spectra continuum.

A careful accounting of these will be necessary for future work that will aim
to constrain the physics of the Lyman-$\alpha$ forest.

##  6 Results

The cross-correlation functions of CHIME and the eBOSS Lyman-$\alpha$ forest
are presented in Figure 6, where the two polarisations have been combined to
approximate total intensity (they are shown separately in Figure 7). Note that
the cross-polar response of the telescope beam has not been calibrated, so
this approximation may include leakage from polarised emission as well. An
estimate of the contamination from noise and residual foregrounds was derived
by cross-correlating the CHIME data against permutations of the Lyman-$\alpha$
forest spectra, i.e. exchanging the spectrum measurements between lines of
sight so that they are uncorrelated with the CHIME observations but maintain
exactly the same distribution on the sky (this is described in more detail
below). The measurements show a clear excess over the background around the
zero lag bin, with the negative sign characteristic of this correlation. Also
shown in this figure is a signal template derived from the result of cross-
correlating the simulated datasets described in the previous section. Only the
amplitude of the template was fit to the data, as explained below, and it
appears to be broadly consistent with the measurement.

Figure 6: _(Top)_ Cross-correlation function of CHIME and eBOSS Lyman-$\alpha$
forest data, for the combined X and Y polarisations. An estimate of the
standard deviation of the background, as described in Section 6.1, is plotted
as a shaded region. The black dashed line is a template derived from
simulations with an amplitude fit to the data. _(Bottom)_ Residuals normalised
by the estimated background level.

Figure 7: Measured cross-correlation functions for each polarisation, along
with the corresponding template fit to that polarisation.

###  6.1 Lyman-$\alpha$ Forest Permutations

In order to characterise the background fluctuations in the cross-correlation
estimates, we would like to generate a catalogue of Lyman-$\alpha$ forest
measurements with the same properties as the eBOSS one but which is
uncorrelated with the CHIME data. Noise, foreground residuals, and any other
systematic effects that are not correlated with the signal will remain present
in the resulting correlation functions, and we can measure their power. We do
this by drawing random permutations of the Lyman-$\alpha$ forest spectra,
while keeping the sample of line-of-sight directions and quasar redshifts
fixed. 1000 such permutations were generated and correlation functions
evaluated for all of them. The standard deviation over permutations at every
spectral separation is shown as the shaded region in Figure 6. Note that
correlations between frequency separations were not considered (i.e. the
covariance was assumed to be diagonal).

Also note that this procedure keeps the selection of lines-of-sight in the
CHIME data fixed, such that if by chance it included data which deviates
significantly from the rest it would be present in every permutation and
potentially bias the estimated variance. To check for this, we generate
another set of permutations where the catalogue positions have also been
rotated by $\sim 3$∘(approximately one CHIME primary beam width at these
wavelengths) in a random direction. We find that the measured variance in the
background fluctuations from these rotated permutations is indistinguishable
from the unrotated set within the sample noise. Thus, we do not expect that
the variance computed from the unrotated permutations is subject to
significant bias due to the fixed set of lines of sight.

###  6.2 Detection Significance

To assess the significance of the detection, we derive a signal template from
the simulations. Synthetic CHIME and Lyman-$\alpha$ forest data were generated
without adding noise and run through the same filtering and cross-correlation
analysis. The resulting correlation functions, the dashed lines in Figure 7,
have a shape that is compact in frequency separation, a result of the delay
filter removing larger scales. These were fit to the cross-correlation
measurements by varying a single amplitude parameter and minimising the
$\chi^{2}$. Finding the $\Delta\chi^{2}=1$ bounds in parameter space provides
an estimate for the standard deviation of the amplitude constraint. The fitted
template and residuals are presented in Figure 6 for the combined
polarisations. We quote the significance of the detection as the amplitude of
the template in units of its standard deviation from the fit in Table 1.

The fitted amplitude of the template to match the measured signal is a factor
of $\sim 4$, i.e. the measurement is about 4 times brighter than the
simulation. Given that the CHIME instrument model used in the simulations is a
simple approximation – in particular, the model for the beam response which is
essential to accurately normalise the signal amplitude – we are not confident
in interpreting this as a discrepancy between the physics of the model and the
observations. However, future work that refines both the instrumental
calibration and modelling should allow a useful amplitude to be constrained.

Table 1: Detection significance of correlation function template fits. The number of degrees of freedom used to compute $\chi_{\nu}$ is 254. measurement | $\mathbf{\chi^{2}}$ | $\mathbf{\chi^{2}_{\nu}}$ | significance
---|---|---|---
X pol | 217.1 | 0.85 | $6.4\sigma$
Y pol | 295.8 | 1.16 | $6.5\sigma$
combined | 238.5 | 0.94 | $9.1\sigma$

##  7 Validation

In this section we describe checks that were performed to validate the
detection and background estimate.

###  7.1 Impact of Selection Function on Catalogue Permutations

To estimate the level of background fluctuations in the correlation function
measurements, random permutations of the Lyman-$\alpha$ forest catalogue were
drawn and correlated against the CHIME data, as described in the previous
section. By construction, this method preserves the distribution of background
quasars on the sky and in redshift separately, but if there exist correlations
between these two axes in the selection function, they will be erased by our
procedure. We check for the presence of such a correlation by binning the
background quasars by their angular position onto a HEALPix (Górski et al.,
2005) grid of ‘nside’ 32 and computing the average redshift in each bin. The
standard deviation across non-empty bins is $\Delta z<0.08$, or about 3% of
the average redshift. For comparison, we shuffled the positions of the
quasars, and repeated the binning calculation, resulting in a typical standard
deviation of just under 3%. There doesn’t appear to be more variation in the
redshifts with location on the sky than would be expected for an uncorrelated
distribution and we conclude that selection effects are unlikely to be
important.

###  7.2 Even/Odd Days

For this test, the full set of days of CHIME data that went into the sidereal
stack were divided into two sets according to an even/odd split in
chronological order, and each was stacked in the same way as the full set (see
Section 2). The 21 cm and foreground signals are the same on every sidereal
day, so they are expected to be common to the even and odd stacks, but thermal
noise and RFI should be uncorrelated between the two. By taking the difference
(and dividing by a factor of two), the common foreground residuals should
cancel and what is left can be compared to our expectation for thermal noise.
Any remaining excess over the noise can be attributed to day-to-day
variations, e.g. RFI or changes to the instrumental response that were not
entirely captured in the daily calibration.

Figure 8: Correlation function evaluated on the difference of even and odd
days, with the standard deviation derived from noise realisations as a shaded
region

The correlation functions evaluated on the even/odd difference are shown in
Figure 8, along with the expected noise level. Not surprisingly, there is no
visible excess around the zero lag bin, which was confirmed by fitting the
signal template to this data, resulting in amplitude parameters that are
consistent with zero within $1\sigma$. Note that the background fluctuations
appear to have a magnitude in excess of the noise. To make this comparison
more clearly, the correlation function was evaluated for 1000 permutations of
the Lyman-$\alpha$ forest spectra for the even-odd difference, as well as for
a CHIME noise realisation. In all cases, the same set of permutations was
used. Figure 9 shows the variance across permutations for these two cases
along with the regular data. The ratio between the variance in the regular
analysis and the even-odd difference indicates the contribution of residual
foregrounds to the background fluctuations in the correlation function. These
account for about twice as much variance as the even-odd difference, which in
turn is about 3-5 times as large as the noise. Note however that the variance
is skewed by significant non-Gaussian tails in the distributions of the pixel
values.

Figure 9: Variance computed across correlation functions estimated on 1000
permutations of the Lyman-$\alpha$ forest spectra. In blue is the full CHIME
data, orange is the difference of even and odd days and green is a noise
realisation. Note that the minimum does not occur at zero separation because
the maximal overlap of the highest weighted regions in the CHIME and eBOSS
spectra is not achieved there. Also, the full range of frequency separations
is shown here, whereas previous plots were cropped to better show the signal.

Figure 10 shows histograms of the pixels in the filtered and masked maps
across all frequencies that were used in the upper and lower bands, within a
region of the sky restricted roughly to the eBOSS NGC field. In the even-odd
difference of the lower band, the distribution is closer to Gaussian for
values $\lesssim 2\sigma$, but large tails remain. The histograms of the upper
band are much closer to Gaussian for $\lesssim 3\sigma$, especially in the
even-odd difference.

Figure 10: Histogram of pixel values in filtered and masked ringmaps across
all frequencies used in _(left)_ the lower band, 400-500 MHz, and _(right)_
the upper band, 585-800 MHz. The pixels are normalised to the expected noise
standard deviation, and a Gaussian curve with a matching amplitude is plotted
in orange for comparison. Rows show the two polarisations, and columns show
the difference of even and odd days and the average of all days. The effect of
the outlier mask can be seen as the sharp cutoff at $6\sigma$.

The analysis of the difference of even and odd days in the lower band suggests
that

  1. 1.

There is a large excess of variance in the maps above the thermal noise (by a
factor $\sim 6-10$).

  2. 2.

This excess is not fully accounted for by features that are repeated day-to-
day like celestial foregrounds (the variance is only $\sim$halved in the even-
odd difference).

  3. 3.

The excess above the noise is non-Gaussian, with large tails in the
distribution of map pixels.

  4. 4.

The distribution is significantly less Gaussian in the lower band than in the
upper band, especially in the even-odd difference.

A possible explanation for this is that the lower band is more contaminated
than the upper band by unflagged RFI – which does not repeat day-to-day and
contributes in a non-Gaussian-distributed way. Another contributing effect
could also be errors in the daily calibration, which would have a similar
signature. Both of these factors have been observed in other internal CHIME
analyses of the lower band. Identifying the source of this contamination and
removing it will likely require further investigation of the individual days,
before they are averaged into a sidereal stack. While these residuals are
concerning and finding their origin in order to remove them will enable more
sensitive analyses, their presence has not impeded the detection in cross-
correlation with the Lyman-$\alpha$ forest reported in this work.

##  8 Conclusion

This work reports on the first detection of 21 cm emission at redshifts
$z>1.5$, obtained by cross-correlating CHIME data from 88 days in the
$400-500$ MHz band with the DR16 eBOSS measurements of the Lyman-$\alpha$
forest. To do so, we extended the analysis methods described in CHIME
Collaboration et al. (2022a) to the lower quarter of the CHIME frequency band,
and developed a cross-correlation method to combine this data with the
Lyman-$\alpha$ forest as a tracer of large-scale structure.

This demonstrates that CHIME is able to detect emission from cosmological HI
in the entire redshift range available to it, but it also highlights the
difficulty of separating this signal from the extremely bright foregrounds.
Even in cross-correlation and after aggressive delay filtering, excess
variance from foregrounds and other sources is found to be about $6-10$ times
above thermal noise. This level of contamination nevertheless allows for a
$\sim 9\sigma$ detection. The shape of the measured correlation function is
dominated by the effect of the foreground filter, and although we have not
attempted to interpret it, we have found that a simple physical model agrees
with the measurement well enough to use this model to quantify the detection
significance. The amplitude of the cross-correlation should in principle tell
us something about the physics of the HI that sources the Lyman-$\alpha$
forest and its relationship to the high-HI-density systems primarily probed
with 21 cm intensity mapping, but interpreting its measurement will require
more sophisticated modelling of the Lyman-$\alpha$ and 21 cm signals. This
will be the subject of future work.

Improvements to the CHIME calibration and analysis pipelines, as well as the
inclusion of several more years of data that have yet to be processed, will
lead to future improvements in the constraining power of the analysis.
Significant residuals in day-to-day jacknives also suggest that the lower band
considered in this work is more contaminated than the upper band that was
reported on previously, and may especially benefit from improved RFI flagging.
In the short term, this might reduce the number of frequencies that need to be
flagged and allow for a less restrictive delay cut, enabling access to
information from scales larger than the $k_{\parallel}\gtrsim 0.13\
\text{Mpc}^{-1}$ that remain in this analysis. These improvements will
position CHIME’s dataset as a powerful source of cross- and auto-correlation
measurements.

## Acknowledgements

We thank Meiling Deng for her substantial contributions to the instrument
design and analysis of CHIME.

We thank the Dominion Radio Astrophysical Observatory, operated by the
National Research Council Canada, for gracious hospitality and expertise. The
DRAO is situated on the traditional, ancestral, and unceded territory of the
Syilx Okanagan people. We are fortunate to live and work on these lands.

CHIME is funded by grants from the Canada Foundation for Innovation (CFI) 2012
Leading Edge Fund (Project 31170), the CFI 2015 Innovation Fund (Project
33213), and by contributions from the provinces of British Columbia, Québec,
and Ontario. Long-term data storage and computational support for analysis is
provided by WestGrid, SciNet and the Digital Research Alliance of Canada, and
we thank their staff for flexibility and technical expertise that has been
essential to this work.

Additional support was provided by the University of British Columbia, McGill
University, and the University of Toronto. CHIME also benefits from NSERC
Discovery Grants to several researchers, funding from the Canadian Institute
for Advanced Research (CIFAR), from Canada Research Chairs, from the FRQNT
Centre de Recherche en Astrophysique du Québec (CRAQ) and from the Dunlap
Institute for Astronomy and Astrophysics at the University of Toronto, which
is funded through an endowment established by the David Dunlap family. This
material is partly based on work supported by the NSF through grants (2008031)
(2006911) and (2006548) and by the Perimeter Institute for Theoretical
Physics, which in turn is supported by the Government of Canada through
Industry Canada and by the Province of Ontario through the Ministry of
Research and Innovation.

We thank the Sloan Digital Sky Survey and eBOSS collaborations for publicly
releasing the Lyman-$\alpha$ forest catalogs used in this work. Funding for
the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan
Foundation, the U.S. Department of Energy Office of Science, and the
Participating Institutions. SDSSIV acknowledges support and resources from the
Center for High Performance Computing at the University of Utah. The SDSS
website is www.sdss.org.

## References

  * Abdurashidova et al. (2022) Abdurashidova, Z., Aguirre, J. E., Alexander, P., et al. 2022, ApJ, 925, 221, doi: 10.3847/1538-4357/ac1c78
  * Anderson et al. (2018) Anderson, C. J., Luciw, N. J., Li, Y. C., et al. 2018, MNRAS, 476, 3382, doi: 10.1093/mnras/sty346
  * Behnel et al. (2011) Behnel, S., Bradshaw, R., Citro, C., et al. 2011, Computing in Science Engineering, 13, 31 , doi: 10.1109/MCSE.2010.118
  * Bi & Davidsen (1997) Bi, H., & Davidsen, A. F. 1997, 479, 523, doi: 10.1086/303908
  * Carucci et al. (2017) Carucci, I. P., Villaescusa-Navarro, F., & Viel, M. 2017, J. Cosmology Astropart. Phys, 2017, 001, doi: 10.1088/1475-7516/2017/04/001
  * Chang et al. (2010) Chang, T.-C., Pen, U.-L., Bandura, K., & Peterson, J. B. 2010, Nature, 466, 463, doi: 10.1038/nature09187
  * CHIME Collaboration et al. (2022a) CHIME Collaboration, Amiri, M., Bandura, K., et al. 2022a, Detection of Cosmological 21 cm Emission with the Canadian Hydrogen Intensity Mapping Experiment, arXiv, doi: 10.48550/arXiv.2202.01242
  * CHIME Collaboration et al. (2022b) —. 2022b, The Astrophysical Journal Supplement Series, 261, 29, doi: 10.3847/1538-4365/ac6fd9
  * Cieplak & Slosar (2016) Cieplak, A. M., & Slosar, A. 2016, 2016, 016, doi: 10.1088/1475-7516/2016/03/016
  * Collette et al. (2021) Collette, A., Kluyver, T., Caswell, T. A., et al. 2021, h5py/h5py: 3.5.0, 3.5.0, Zenodo, doi: 10.5281/zenodo.5585380
  * Cunnington et al. (2023) Cunnington, S., Li, Y., Santos, M. G., et al. 2023, MNRAS, 518, 6262, doi: 10.1093/mnras/stac3060
  * Dalcin & Fang (2021) Dalcin, L., & Fang, Y.-L. L. 2021, Computing in Science Engineering, 23, 47, doi: 10.1109/MCSE.2021.3083216
  * du Mas des Bourboux et al. (2020) du Mas des Bourboux, H., Rich, J., Font-Ribera, A., et al. 2020, ApJ, 901, 153, doi: 10.3847/1538-4357/abb085
  * Ewall-Wice et al. (2021) Ewall-Wice, A., Kern, N., Dillon, J. S., et al. 2021, MNRAS, 500, 5195, doi: 10.1093/mnras/staa3293
  * Farr et al. (2020) Farr, J., Font-Ribera, A., Bourboux, H. d. M. d., et al. 2020, 2020, 068, doi: 10.1088/1475-7516/2020/03/068
  * Gabriel et al. (2004) Gabriel, E., Fagg, G. E., Bosilca, G., et al. 2004, in Proceedings, 11th European PVM/MPI Users’ Group Meeting, Budapest, Hungary, 97–104
  * Górski et al. (2005) Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759, doi: 10.1086/427976
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Karamanis & Beutler (2021) Karamanis, M., & Beutler, F. 2021, arXiv e-prints, arXiv:2106.06331.  https://arxiv.org/abs/2106.06331
  * Leifer, C. (2021) Leifer, C. 2021, peewee/peewee, latest, Github.  https://github.com/coleifer/peewee
  * Lewis et al. (2000) Lewis, A., Challinor, A., & Lasenby, A. 2000, ApJ, 538, 473, doi: 10.1086/309179
  * Li et al. (2021) Li, L.-C., Staveley-Smith, L., & Rhee, J. 2021, Research in Astronomy and Astrophysics, 21, 030, doi: 10.1088/1674-4527/21/2/30
  * Liang et al. (2021) Liang, Y., Kashikawa, N., Cai, Z., et al. 2021, ApJ, 907, 3, doi: 10.3847/1538-4357/abcd93
  * Masui et al. (2015) Masui, K., Amiri, M., Connor, L., et al. 2015, Astronomy and Computing, 12, 181, doi: 10.1016/j.ascom.2015.07.002
  * Masui et al. (2013) Masui, K. W., Switzer, E. R., Banavar, N., et al. 2013, ApJ, 763, L20, doi: 10.1088/2041-8205/763/1/L20
  * McKinney (2010) McKinney, W. 2010, in Proceedings of the 9th Python in Science Conference, ed. Stéfan van der Walt & Jarrod Millman, 56 – 61, doi: 10.25080/Majora-92bf1922-00a
  * Mead et al. (2021) Mead, A. J., Brieden, S., Tröster, T., & Heymans, C. 2021, MNRAS, 502, 1401, doi: 10.1093/mnras/stab082
  * Momose et al. (2021) Momose, R., Shimasaku, K., Kashikawa, N., et al. 2021, ApJ, 909, 117, doi: 10.3847/1538-4357/abd2af
  * Mukae et al. (2017) Mukae, S., Ouchi, M., Kakiichi, K., et al. 2017, ApJ, 835, 281, doi: 10.3847/1538-4357/835/2/281
  * pandas development team (2020) pandas development team, T. 2020, pandas-dev/pandas: Pandas, latest, Zenodo, doi: 10.5281/zenodo.3509134
  * Paul et al. (2023) Paul, S., Santos, M. G., Chen, Z., & Wolz, L. 2023, arXiv e-prints, arXiv:2301.11943, doi: 10.48550/arXiv.2301.11943
  * Pen et al. (2009) Pen, U.-L., Staveley-Smith, L., Peterson, J. B., & Chang, T.-C. 2009, MNRAS, 394, L6, doi: 10.1111/j.1745-3933.2008.00581.x
  * Peterson et al. (2009) Peterson, J. B., Aleksan, R., Ansari, R., et al. 2009, in astro2010: The Astronomy and Astrophysics Decadal Survey, Vol. 2010, 234, doi: 10.48550/arXiv.0902.3091
  * Rhodes (2019) Rhodes, B. 2019, Skyfield: High precision research-grade positions for planets and Earth satellites generator.  http://ascl.net/1907.024
  * Seljak (2012) Seljak, U. 2012, 2012, 004, doi: 10.1088/1475-7516/2012/03/004
  * Shaw et al. (2020a) Shaw, J. R., Foreman, S., Nitsche, R., et al. 2020a, radiocosmology/driftscan: 20.5.0, v20.5.0, Zenodo, doi: 10.5281/zenodo.5949193
  * Shaw et al. (2020b) Shaw, J. R., Nitsche, R., Foreman, S., & Kefala, A. 2020b, radiocosmology/cora: 20.5.0, v20.5.0, Zenodo, doi: 10.5281/zenodo.5949195
  * Shaw et al. (2020c) Shaw, J. R., Masui, K., Nitsche, R., et al. 2020c, radiocosmology/caput: 20.10.0, v20.10.0, Zenodo, doi: 10.5281/zenodo.5846375
  * Shaw et al. (2020d) Shaw, J. R., Masui, K., Hincks, A. D., et al. 2020d, chime-experiment/ch_pipeline: 20.10.0, v20.10.0, Zenodo, doi: 10.5281/zenodo.5846379
  * Shaw et al. (2020e) Shaw, R., Nitsche, R., Siegel, S. R., et al. 2020e, radiocosmology/draco: 20.10.0, v20.10.0, Zenodo, doi: 10.5281/zenodo.5828373
  * The HDF Group (1997) The HDF Group. 1997, Hierarchical Data Format, version 5
  * Tramonte & Ma (2020) Tramonte, D., & Ma, Y.-Z. 2020, MNRAS, 498, 5916, doi: 10.1093/mnras/staa2727
  * 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
  * Wolz et al. (2022) Wolz, L., Pourtsidou, A., Masui, K. W., et al. 2022, MNRAS, 510, 3495, doi: 10.1093/mnras/stab3621
  * Zonca et al. (2019) Zonca, A., Singer, L., Lenz, D., et al. 2019, Journal of Open Source Software, 4, 1298, doi: 10.21105/joss.01298