driftscan [Shaw et al., 2020],
emcee [Foreman-Mackey et al., 2013],
GetDist [Lewis, 2019],
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],
networkx [Hagberg et al., 2008],
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],

§ DELAY POWER SPECTRUM ESTIMATION VIA GIBBS SAMPLING

Delay power spectra
[For clarity, we will use delay spectrum to refer only to the direct Fourier transform of a frequency intensity or flux spectrum. The delay power spectrum will refer only to the variance of this quantity. Though the intensity and flux are both second-order statistics of the electric field and thus are power-like quantities in a physical sense, we do not think this is ambiguous anywhere in this text.]
measure the power at different time lags observed within a frequency spectrum, and are an extremely powerful tool for investigating instrumental effects as well as the frequency structure of radio emission from the sky (see <ref> for an example). Superficially estimating a delay power spectrum involves taking a Fourier transform of a frequency spectrum and estimating the resulting power in the time domain. However, in the presence of interference that causes certain frequencies to be masked out, and a large dynamic range between the power at different delays, significant care must be taken to avoid mixing of power between different delays. There are several existing strategies for dealing with this such as using a CLEAN-like algorithm in delay space [Parsons et al., 2014], and Least-squares Spectral Analysis <cit.>.

To understand the challenges involved, consider a noisy observation $\vf$ of a frequency spectrum with length $N_\mathrm{f}$. This is related to an underlying delay spectrum $\vd$ by

\begin{equation}
\label{eq:def_fi}
\vf = \mF \vd + \vn
\end{equation}

with noise $\vn$ and where the delay spectrum is assumed to be drawn from the input delay power spectrum $D[\tau_a]$:

\begin{equation}
D[\tau_a] = \left\langle \big| d_a \big|^2 \right\rangle\ .
\end{equation}

The noise is described by covariance $\mN$, and for the moment we treat the noise as being uniform except for entirely missing frequencies which we give infinite noise. We write this as $\mN^{-1} = \sigma^{-2} \mM$ where $\mM$ is a diagonal masking matrix with ones for included frequencies and zeros for missing frequencies. The matrix $\mF$, with  $F_{ab} = e^{-2 \pi j \tau_a \nu_b} / N_\mathrm{f}^{1/2}$, is unitary and performs a discrete Fourier transform from the time (delay) domain to the conjugate frequency domain.

A first attempt to estimate the delay power spectrum might start by simply applying the mask $\mM$ to the observed frequency spectrum and performing an inverse Fourier transform,

\begin{equation}
\hat{\vd}_\text{inv} = \mF^\hconj \mM \vf \ .
\end{equation}

With this estimate of the delay spectrum, we can then infer the delay power spectrum $D[\tau_a]$ by using a variance over $N_\mathrm{obs}$ observations, indexed by $i$:

\begin{equation}
\label{eq:D_samplevar}
\hat{D}[\tau_a] = \frac{1}{N_\mathrm{obs}} \sum_{i=1}^{N_\mathrm{obs}} \lv \hat{d}^i_a \rv^2 \ ,
\end{equation}

with $\hat{d}^i_a$ being any estimator for $d_a$ such as $\hat{\vd}_\text{inv}$ defined above (later we will introduce additional estimators). However, this procedure generates significant leakage between delay channels, with a delay spread function given by the matrix $\mF^\hconj \mM \mF$. In the case of random masking of $N_\text{masked}$ single frequencies, it can be shown that this gives leakage at the level of $\sim N_\text{masked} \sum_a D[\tau_a] / N_\mathrm{f}^2 $ uniformly across delays, and we are not able to see any structure in the delay power spectrum below this level.

To improve this, we could modify the delay spectrum estimate by deconvolving the delay spread function by its pseudo-inverse, or equivalently use a maximum-likelihood estimator[We note that this is similar to LSSA, though LSSA considers more general cases such as irregular sampling, and typically restricts the range of delays being solved for to minimise correlations and leakage.]

\begin{equation}
\hat{\vd}_\text{ml} = (\mF^\hconj \mN^{-1} \mF)^+ \mF^\hconj \mN^{-1} \, \vf \ .
\end{equation}

However, as $\mF$ is unitary, the pseudo-inverse $(\mF^\hconj \mN^{-1} \mF)^+$ is equal to $\mF^\hconj (\mN^{-1})^+ \mF$, and as the noise matrix is diagonal with zeros where samples are masked, $(\mN^{-1})^+ \mN^{-1} = \mM$; together, these imply that $\hat{\vd}_\text{ml} = \hat{\vd}_\text{inv}$. In words, the maximum-likelihood estimator is exactly equivalent to inverse Fourier transforming the masked frequency spectra.

Another option is using a Wiener filter instead of a maximum-likelihood type filter:

\begin{equation}
\label{eq:vd_wiener}
\hat{\vd}_\text{w} = (\mD^{-1} + \mF^\hconj \mN^{-1} \mF)^{-1} \mF^\hconj \mN^{-1} \vf \ ,
\end{equation}

where $D_{ab} = D[\tau_a] \delta_{ab}$ is the covariance matrix of the delay spectrum signal. By providing information about the distribution of power at various delays, the filter can distinguish delays related to true signal in the masked frequency spectra, resulting in delay spectra with significantly lower leakage and hence cleaner power spectra. However, constructing this requires that we already know the delay power spectrum $D[\tau_a]$, which is the quantity that we are trying to estimate. A close enough guess may minimise the leakage enough to produce accurate delay power spectrum estimates, but there is no knowing in advance if this is the case.

A resolution to this is to jointly solve for both the delay spectrum and the delay power spectrum, a problem which is tractable by Gibbs sampling [Geman & Geman, 1984], a Markov Chain Monte Carlo technique for drawing samples from a joint distribution where the conditional distributions are easily sampled. In particular, we draw inspiration from techniques used for power spectrum estimation of the Cosmic Microwave Background (e.g. Eriksen et al., 2004, Wandelt et al., 2004).

We want to infer both the delay spectrum $\vd$ and the delay power spectrum (equivalent to the diagonal matrix $\mD$) by drawing samples from the joint probability distribution $\calP(\vd, \mD \mid \vf)$. Gibbs sampling allows us to do that by alternately drawing from the conditional distributions $\calP(\vd \mid \mD, \vf)$ and $\calP(\mD \mid \vd, \vf)$; the ensuing set of samples will eventually converge to the joint distribution, and we can take the mean over $\mD$ samples as an estimate of the delay power spectrum. We now describe how to sample from each conditional distribution.

Starting with $\calP(\vd \mid \mD, \vf)$, we can use Bayes' theorem to write

\begin{equation}
\calP(\vd \mid \mD, \vf) \propto \calP(\vf \mid \vd, \mD) \, \calP(\vd \mid \mD) \ .
\end{equation}

The first term on the right hand side is the likelihood function for the frequency spectrum, which for Gaussian noise can be written as $\calP(\vf \mid \vd, \mD) = \calG_C(\vf - \mF \vd, \mN)$, where $\calG_C$ is a circularly symmetric complex Gaussian distribution:

\begin{equation}
\calG_C(\vz, \mC) = \frac{1}{\lv \pi \mC \rv} e^{-\vz^\hconj \mC^{-1} \vz} \ .
\end{equation}

We will also model the conditional prior distribution for the delay spectrum as Gaussian, with $\calP(\vd \mid \mD) = \calG_C(\vd, \mD)$. Combining these together and grouping the terms in $\vd$, we find that the conditional distribution is

\begin{equation}
\calP(\vd \mid \mD, \vf) = \calG_C(\vd - \hat{\vd}_\text{w}, \mC)\ ,
\end{equation}

where $\mC^{-1} = \mD^{-1} + \mF^\hconj \mN^{-1} \mF$. Thus, the mean of the conditional distribution is just the Wiener filter of <ref>, with the standard covariance. Although drawing from this can be done by solving for the mean, followed by inversion and factorization of $\mC^{-1}$ to add a random fluctuation, it is more efficiently done by constructing

\begin{equation}
\mC^{-1} \, \vd = \mF^\hconj \mN^{-1} \vf + \mD^{-1/2} \vw_1 + \mF^\hconj \mN^{-1/2} \vw_2
\label{eq:specsample}
\end{equation}

where $\vw_1$ and $\vw_2$ are standard Gaussian random samples, and then solving for $\vd$ [Jewell et al., 2004].

The conditional distribution for the delay power spectrum is more straightforward. We wish to calculate the conditional distribution $\calP(\mD \mid \vd, \vf)$, which is independent of $\vf$ as all the information about $\mD$ is contained within $\vd$. Using a flat prior on the elements of $\mD$, and the prior $\calP(\vd \mid \mD)$ we find that $\calP(\mD \mid \vd, \vf) \propto \calG_C(\vd, \mD)$.
Assuming that $\mD$ is diagonal we can rewrite this in terms of the sample variance estimates $\hat{D}[\tau_a]$ for each delay $\tau_a$ which are sufficient statistics for the diagonal elements of $\mD$ itself, $D[\tau_a]$. The sample variance $\hat{D}[\tau_a]$ has a chi-squared distribution,

\begin{equation}
N_\mathrm{obs} \frac{\hat{D}[\tau_a]}{D[\tau_a]} \sim \chi^2(N_\mathrm{obs}) \ ,
\label{eq:powspecchi2}
\end{equation}

and so we can draw samples from $\calP(\mD \mid \vd) \propto \calP(\hat{D}[\tau_a] \mid D[\tau_a])$ by drawing a standard chi-squared deviate for each delay $x_a \sim \chi^2(N_\mathrm{obs})$ and setting the new sample for $D[\tau_a]$ to be $N_\mathrm{obs} \hat{D}[\tau_a] / x_a$.

Our practical implementation of this algorithm is as follows:

* Pick a set of data whose delay spectra are expected to be similar enough that we can average over them. For computing delay spectra from visibilities, this might consist of all RA samples for individual baselines (after stacking over redundant copies). For the map delay spectrum described in sec:filter, we choose the set of RA samples at each polarization and declination. As above, we use $N_\mathrm{obs}$ to denote the size of this set.

* Apply an apodization window to each frequency spectrum, if desired. A Nuttall window is used sec:filter.

* Choose an initial guess $D_0[\tau_a]$ for the delay power spectrum. We use a white spectrum with amplitude 10 in this work.

* Loop over the following steps until convergence has been achieved:

* For each element $i$ of the set of $N_\mathrm{obs}$ spectra, draw the $n$th delay spectrum sample $\vd_n^i \leftarrow \calP(\vd \mid \vf, \mD_{n})$ using <ref>. Note that each $\vd_n^i$ is a delay spectrum with $N_\mathrm{f}$ elements.

* Draw the $(n+1)$th delay power spectrum sample $D_{n+1}[\tau_a] \leftarrow \calP(D[\tau_a] \mid \{\vd_n^i\}_{i=1}^{N_\mathrm{obs}})$, using the $N_\mathrm{obs}$ delay spectra drawn at step $n$ to compute $\hat{D}[\tau_a]$ in <ref>.

* Take the average of the converged samples, after removing burn-in and performing any necessary thinning. In sec:filter, we halt after 100 samples, and take the median over the final 50 samples as an estimate of the delay power spectrum.

In summary, the Gibbs sampling approach is a statistically well-motivated technique that iteratively deconvolves the delay spectra, uses them to update a delay power spectrum, and uses this to improve the next deconvolution round.

In <ref> we apply the various estimators discussed above  to a synthetic dataset with high dynamic range in delay space and a realistic frequency mask. We clearly see that the Gibbs sampling based estimator is able to accurately recover the input spectrum, while the naïve inverse-Fourier and Wiener estimators show various degrees of discrepancy.

To test the performance of the delay power spectrum estimation techniques discussed in Appendix <ref>, we generate a set of random delay spectra with a true delay power spectrum (black line) consisting of very high power at low delays, and a plateau of low power outside this region. We Fourier transform these into frequency spectra and apply the CHIME RFI mask used in this analysis (sec:rfi). In orange we show the direct inverse estimate, which has significant leakage at the $\sim 10^{-2}$ level. The Wiener filter estimate (green solid) produces a much closer estimate, correcting most of the leakage effects, at the expense of needing a good starting guess. If generated with a poor initial delay power spectrum (high power over twice the range of delays as the true power spectrum) the estimate is significantly worse (green dashed). The Gibbs sampler produces an estimate that is much closer to the true power spectrum (red). Like any MCMC scheme, attention must be paid to the convergence of the chain. We used 100 samples and derived our estimate from the median of the last half. The full chain is shown in gray and can be seen converging from a poor starting guess of a flat delay power spectrum to the true power spectrum over $\sim 30$ samples.

§ ESTIMATING THE PRIMARY BEAM BY DECONVOLVING A MODEL FOR THE POINT SOURCE SKY

In this appendix, we describe the algorithm that is used to directly reconstruct the average primary beam pattern of the CHIME antennas.  First, a model for the radio emission from extragalactic point sources is constructed from measurements made by other telescopes.  The specfind v2 table [Vollmer et al., 2010] in the Vizier database is queried for flux measurements of all known sources between declinations -4085.  For each source, all available measurements of the flux are fit to a power-law with frequency

\begin{align}
s(\nu) & = a \left(\frac{\nu}{\SI{600}{\mega\hertz}}\right)^{\gamma} \ ,
\end{align}

where the amplitude $a$ and exponent $\gamma$ are allowed to float.  The fit is done by performing a weighted linear regression of the logarithm of the flux to the logarithm of the frequency.  The uncertainties provided in the specfind v2 table are used to construct inverse variance weights.  These uncertainties are 20 of the measured flux [Vollmer et al., 2005], and the power-law model is in general a good fit given these large uncertainties.  Only sources with $s(\SI{600}{\mega\hertz}) > \SI{15}{\milli\jansky}$ that have at least one measurement on either side of the CHIME band are included in the sky model.  There are 97941 sources in total that meet these criteria.  All of the sources have at least 3 flux measurements, with 6 flux measurements on average.

Our model for the visibility measured by baseline $\vec{b}$ at frequency $\nu$ and local Earth rotation angle $\phi$ is then given by

\begin{align}
\label{eq:sky_model}
S(\vec{b}, \nu, \phi) & = \sum_{i} s_{i}(\nu) \ e^{j 2 \pi \nu \vec{b} \cdot \vec{\hat{n}}(\theta_{i}, \phi - \phi_{i}) / c} \ \delta\left(\phi, \phi_{i}\right)
\end{align}

where $s_{i}(\nu)$ is a power-law model for the flux of the $i$'th source, $\vec{\hat{n}'}(\theta_{i}, \phi - \phi_{i})$ is the unit vector pointing in the direction of the $i$'th source and is given by <Ref> with $\phi_{i}$ and $\theta_{i}$ denoting the source's right ascension and declination in CIRS coordinates, and

\begin{align}
\delta\left(\phi, \phi_{i}\right) & = \begin{cases}
1 & |\phi - \phi_{i}| < \frac{1}{2} \Delta \phi \\
0 & \mbox{otherwise} \\
\end{cases}
\end{align}

with $\Delta \phi = \SI{0.0879}{\degree}$ denoting the sample spacing of the data in local sidereal angle.  The sum in <Ref> runs over all sources.

The following identical operations are then performed on the sidereal visibilities $V$ and sky model $S$.  First we arrange the baselines onto a 2D grid and then beamform in the $\vec{\hat{y}}$ direction using <Ref>.  The weights used in the beamformer are given by

\begin{align}
\label{eq:weights_ns_beamcal}
w_{xy}^{p}(\nu, \phi) & = W_{\alpha}\left(\frac{1}{2} \left[1 + \frac{\nu y}{\nu_{\rm min} y_{\rm max}}\right]\right)
\end{align}

where $W_{\alpha}$ denotes the Dolph-Chebyshev window, $\alpha = \SI{60}{\decibel}$ is the peak-to-sidelobe ratio, $\nu_{\rm min} = \SI{587.5}{\mega\hertz}$ is the minimum frequency examined, and $y_{\rm max} = 255$ corresponds to the maximum baseline distance in the $\vec{\hat{y}}$ direction.

The window function in <Ref> will result in a frequency-independent synthesized beam in the $\hat{\theta}$ direction that has a $\mbox{FWHM} = \SI{0.385}{\degree}$ and sidelobes that are $\lesssim 10^{-3}$ of the peak amplitude.  The Dolph-Chebyshev window minimizes the main lobe width for a given number of baselines and equiripple peak-to-sidelobe ratio.  It will degrade the point-source sensitivity relative to the inverse variance weighting scheme discussed in <ref>, however the loss of sensitivity is not problematic for beam calibration because it relies on a foreground signal that is $\gtrsim 500$ times brighter than the noise.  The low equiripple sidelobes help to ensure that each formed beam is sensitive to the primary beam at a narrow range of declinations.

The argument of the window function is scaled with frequency so that the synthesized beam in the $\hat{\theta}$ direction is frequency independent.  Essentially the resolution at every frequency is degraded to the resolution at the lowest frequency.  This ensures that all frequencies are sensitive to the same declinations, so that any errors in our sky model are not further modulated by a frequency-dependent synthesized beam pattern.

A schematic representation of the construction of the primary beam model used in our stacking analysis. Externally-measured spectra of 97941 radio point sources are propagated into mock visibilities, which, after several transformations, are cross-correlated with CHIME observations to construct a beam transfer function that assumes that the sky is solely composed of these sources. Further transformations are applied to minimize sensitivity to this assumption and remove artifacts.

Next we multiply the hybrid beamformed visibilities by a cosine-tapered window that is unity for $\SI{125}{\degree} < \phi < \SI{255}{\degree}$ and transitions to zero over a span of 15.  This restricts our attention to a relatively quiet portion of the radio sky, avoiding sharp features in the Galactic emission that are present in the data but not in our model, and also avoiding regions of the sky contaminated by Cygnus A and Casseopia A in the sidelobes, which this technique is unable to account for properly.  This range of $\phi$ also coincides with the range covered by the eBOSS NGC field.  The $m$-mode transform is then taken.

The model for the primary beam is obtained by cross-correlating the sky model and the visibilities in $m$-mode space,

\begin{align}
\tilde{B}^{p}_{xm}(\nu, \theta) & = \frac{\tilde{S}^{*}_{xm}(\nu, \theta) \ \tilde{V}^{p}_{xm}(\nu, \theta)}{|\tilde{S}_{xm}(\nu, \theta)|^{2} + \sigma^{p}_{x}(\nu)^{2}} \ ,
\end{align}

where $\tilde{S}$, $\tilde{V}$, and $\tilde{B}$ denote the $m$-mode transform of the hybrid beamformed visibilities for the sky model, data, and beam transfer function, respectively, and $\sigma$ is an estimate of the noise in $\tilde{V}$.

The $m$-mode transform of the beam transfer function, $\tilde{B}$, is multiplied by a cosine-tapered mask to remove any $m$-modes that cannot originate from the sky near meridian.  This mask is unity for $m_{{\rm center}, x} \pm 0.75 m_{\rm width}$ and then smoothly transitions to zero by $m_{{\rm center}, x} \pm m_{\rm width}$ (see <Ref> and <Ref>).  The inverse $m$-mode transform is then calculated to obtain our estimate of the beam transfer function $B$ for each east-west baseline separation $x$.  The beam transfer function is then “fringestopped”, or in other words, is multiplied by the complex conjugate of the exponential term in <Ref>, to recover $|A^{p}(\nu, \theta, \phi)|^{2}$, which we will refer to as the power beam.  Note that this method yields 4 distinct estimates of the power beam, one for each east-west baseline separation.

We find that the resulting estimate of the power beam exhibits small scale variations along the declination axis that are highly correlated as a function of frequency and hour angle.  We suspect that these variations are due to errors in the flux of the sources in the sky model, and remove them as follows.  At each declination, the logarithm of the power beam at $\phi = \SI{0}{\degree}$ is fit to a fourth order polynomial in frequency.  This logarithmic polynomial model is then high-pass filtered along the $\hat{\theta}$ direction so that only variations on scales $\lesssim \SI{3}{\degree}$ are preserved.  The power beam at each declination is then divided by the exponential of the high-pass filtered, logarithmic polynomial model.

The uncertainty in the power beam is estimated at each frequency and declination by examining the variance at large hour angle ($0.087 \leq |\cos{\theta} \sin{\phi}| \leq 0.42$).  This uncertainty varies significantly as a function of declination based on the brightness of the sources at that declination.  We apply a 2D Savitzky-Golay filter in ($\nu$, $\theta$) space to low-pass filter the beam model.  For each ($\nu$, $\theta$, $\phi$) a 4th-order Chebyshev polynomial in both $\nu$ and $\theta$ is fit to a small window centered on that location.  The best-fit polynomial model is evaluated at that location to obtain the low-pass filtered version of the beam model.  The variance in the beam model at large hour angle is used to estimate the weights in the fit and properly account for the declination dependent uncertainties.  The size of the window changes between three distinct values based on the declination and frequency in order to retain features in the beam at progressively smaller scales as one moves to lower declinations.  In addition to smoothing the beam, the low-pass filter interpolates the beam to the majority of the frequencies that have been masked because of missing data or RFI.

Even after applying the 2D smoothing operation, there are still sharp features in the beam along the frequency axis that we believe originate from unflagged RFI present in the sidereal visibilities.  These sharp features will leak foreground power to small-spectral scales when the beam model is deconvolved from the data.  To address this, at each $(\theta, \phi)$ we apply an 8th order low-pass Butterworth filter along the frequency axis.  The cutoff used  for the low-pass filter is declination dependent in order to retain what we suspect are actual features of the beam.  The cutoff ranges from 125200.

The final estimate of power beam is obtained from a weighted average of the estimate from baselines with a 44 and 66 east-west component.  The baselines with 0 east-west component are contaminated by diffuse Galactic emission, which is not present in our sky model, and also coupled noise that varies slowly as a function of Earth-rotation angle and thus appears at low-$m$'s that overlap with the range of $m$ at which the meridian sky fringes. The baselines with 22 east-west components are also contaminated by coupled noise, albeit to a lesser extent.

This technique is currently unable to measure the primary beam accurately at hour angles greater than $\approx \SI{2.0}{\degree}$, where the first-order approximation for the geometric phase given in <Ref> begins to break down and an additional term that depends on the north-south baseline distance, declination, and hour angle becomes relevant.  The phase due to this term will be equal to the first-order phase at a new “effective” declination given by

\begin{eqnarray}
\theta_{\rm eff}' & = \arcsin{\left( \cos{\Lambda} \sin{\theta'} - \sin{\Lambda} \cos{\theta'} \cos{\left(\phi - \phi'\right)}\right)} + \Lambda \ .
\end{eqnarray}

As a result, bright sources will exhibit a “U” shape track in the hybrid beamformed visibilities as they move out of the meridian beam centered on their true declination $\theta'$ and into meridian beams centered on more northern declinations at $\theta_{\rm eff}'$.  The recovered primary beam will be attenuated at large hour angles by a factor $\synth^{\hat{\theta}, p}(\nu, \theta, \theta_{\rm eff}', \phi) \ / \ \synth^{\hat{\theta}, p}(\nu, \theta, \theta', \phi)$.  In the main lobe of the primary beam the attenuation is less than 6 for the polarisations, frequencies, and declinations considered in this work, but in the sidelobes it quickly becomes significant.  We are actively exploring extensions to this algorithm that are capable of recovering the side lobes as well.

§ STACKING ON LOGNORMAL GALAXY DENSITY REALIZATIONS

In sec:skymap_generation, we made the following statement: if simulated galaxy catalogs are drawn from lognormal realizations of the galaxy density $\deltag$, and correlated Gaussian-distributed HI maps are stacked on the resulting galaxy positions, the measured stacking signal is the same as it would be if the galaxy catalogs were drawn from Gaussian realizations of the galaxy density. In this appendix, we justify this statement. We will make use of the following Gaussian integrals: if $\vec{\delta}$, $\vec{\alpha}$ are $n$-component vectors and $\mC$ is a symmetric, positive-definite $n\times n$ matrix, we can write

\begin{align}
\label{eq:gaussianint1}
\int d^n\vec{\delta}\, e^{-\frac{1}{2} \vec{\delta}^T \mC^{-1} \vec{\delta} + \vec{\alpha} \cdot \vec{\delta}}
&= \sqrt{ (2\pi)^{n} \det \mC }\, e^{ \frac{1}{2} \vec{\alpha}^T \mC \vec{\alpha} }\ , \\
\label{eq:gaussianint2}
\int d^n\vec{\delta}\, \delta_i\, e^{-\frac{1}{2} \vec{\delta}^T \mC^{-1} \vec{\delta} + \vec{\alpha} \cdot \vec{\delta}}
&= \sqrt{ (2\pi)^{n} \det \mC }\, e^{ \frac{1}{2} \vec{\alpha}^T \mC \vec{\alpha} }
\sum_j C_{ij} \alpha_j\ .
\end{align}

Consider an idealized version of the stacking analysis, in which we average the HI overdensity $\deltaHI$ at a 3d separation $\vr$ from the location of each of $N$ galaxies in a catalog:
S() = 1/N ∑_i=1^N (_i+) .

Suppose that $\deltaHI$ has Gaussian statistics, while the galaxy overdensity $\deltag$, from which the galaxy positions are drawn, is lognormal, related to a Gaussian field $\deltaG$ by
1+() = e^() + μ() .

Also, let $C_{ab}(\vx,\vx') \equiv \langle \delta_a(\vx) \delta_b(\vx') \rangle$, where $a, b \in \{ {\rm G}, {\rm HI} \}$. Following the standard procedure for lognormal fields, we fix $\mu(\vx)$ in <ref> such that $\langle 1+\deltag(\vx) \rangle = 1$. We can compute the relevant ensemble average by defining $\vec{\delta}$ to be $\deltag$ evaluated at a finite number of points, and writing

\begin{align*}
\left\langle 1+\deltag(\vx) \right\rangle
&= \int d^n\vec{\delta}\, \left[ 1+\deltag(\vx) \right] \calP(\vec{\delta}) \\
&= \int d^n\vec{\delta}\, e^{\deltaG(\vx) + \mu(\vx)}
\frac{1}{\sqrt{(2\pi)^n \det \mC_{\rm GG} }}
e^{ -\frac{1}{2} \vec{\delta}^T \mC_{\rm GG}^{-1} \vec{\delta} } \\
&= e^{\mu(\vx) + \frac{1}{2} C_{\rm GG}(\vx, \vx) }\ ,
\numberthis
\end{align*}

where in the second equality, we substituted <ref> and the pdf for a Gaussian random field, and in the final equality we used <ref>. Setting this to unity implies that $\mu(\vx) = -(1/2) C_{\rm GG}(\vx, \vx)$.

We wish to show that stacked HI overdensity in <ref> approaches the same result whether the galaxy positions are drawn from the lognormal field in <ref> or from the Gaussian field $\deltaG$ itself. To do so, we first consider an ensemble average over galaxy positions in the catalog, keeping the underlying fields ($\deltag$ and $\deltaHI$) fixed:
⟨S() ⟩_cat
= 1/N ∑_i=1^N  ⟨(_i+) ⟩_cat
= 1/N ∑_i=1^N ∫_V d^3_i  (_i + )   (_i | ) ,

with the integral evaluated over the survey volume $V$.
The probability distribution function of galaxy positions, $\calP(\vx | \deltag)$, is given by $\calP(\vx | \deltag) = V^{-1} [1+\deltag(\vx)]$. By combining Eqs. (<ref>)-(<ref>) and the expressions for $\calP(\vx | \deltag)$ and $\mu(\vx)$, we obtain
⟨S() ⟩_cat = 1/V ∫_V d^3
(+) e^() - 1/2 C_GG(, ) .

We now take the ensemble average of <ref> over the density fields as well as the catalog positions. Defining $\vec{\delta}$ to contain both $\deltaHI$ and $\deltaG$, and letting $\mC$ represent the joint covariance, we have

\begin{align*}
\left\langle \deltaHI(\vx+\vr) e^{\deltaG(\vx)} \right\rangle_{\rm fields}
&= \int d^n\vec{\delta}\, \deltaHI(\vx+\vr) e^{\deltaG(\vx)}
\frac{1}{\sqrt{(2\pi)^n \det \mC }}
e^{ -\frac{1}{2} \vec{\delta}^T \mC^{-1} \vec{\delta} } \\
&= e^{ \frac{1}{2} C_{\rm GG}(\vx, \vx) } C_{\rm HI, G}(\vx+\vr, \vx)\ ,
\numberthis
\end{align*}

where we used <ref> in the second equality. Combining this with <ref> and assuming that the fields have translation-invariant statistics, we arrive at
⟨S() ⟩_cat, fields = C_HI, G(, 0⃗) ,
which is also what we would obtain if the galaxy positions were drawn directly from $\deltaG$ itself.

§ TEMPLATE CALCULATION

In this appendix, we discuss the challenge of calculating the signal templates for arbitrary parameter combinations and the approach we take in this work. Other than the frequency bias parameter $\Dnu$, the parameters described in sec:freeparameters affect the properties of the underlying large-scale structure, or the or tracer density fields. That suggests that one way of calculating the template is to produce a realisation of the field and a correlated tracer catalog given a set of parameters, and then simulate a CHIME timestream from the field using a model for the instrumental transfer function, repeat the analysis procedure done to the actual data (flagging, filtering and map making), and finally stacking the output on the mock catalog. By repeating this procedure and averaging the results, we can estimate the expected signal.

Unfortunately, a full Monte-Carlo of this procedure is challenging, as even a single iteration requires around 900 core-hours of compute time (dominated by the timestream generation from input sky maps).
Instead we utilize the ergodic principle. We have an overlapping volume of $\gtrsim \SI{10}{\giga\parsec^3\per\h^3}$ (covered by the eBOSS quasar sample), but the stacking is probing scales $\sim \SI{10}{\mega\parsec\per\h}$. This gives many quasi-independent regions of that size within the volume, and so on those scales we expect the volume average to approach the ensemble average, or equivalently, that averaging over independent mock source catalogs drawn from a single large-scale structure realisation should give the same as averaging over completely independent large-scale structure realisations. Though this naive picture will break down on larger scales where the cosmic variance contribution is significant, we find that the cosmic variation in the stack signal is small:
around 0.8% for the LRG sample, 0.2% for the ELGs and 0.3% for QSOs (estimated by comparing the zero-lag amplitude for stacks drawn from distinct LSS realisations).
This is no more than the variation between single catalogs drawn from the same LSS realisation
($\sim 0.7\%$ for LRGs and QSOs, $\sim 1.2\%$ for ELGs)
although we average over a sufficiently large number of catalogs to reduce this contribution to well below the cosmic variance level.

While this gives us a tractable method of computing the template for a given set of parameters, it still requires a costly timestream simulation for each set of values. To avoid this, we note that as both the process of observation and analysis are linear (other than data-derived RFI and bright pixel masking), if we can isolate individual terms in the cross-power spectrum description, they map to distinct contributions to the stacked signal.

For the moment we will fix the Fingers of God parameters $\alphaFoG{\sHI}$ and $\alphaFoG{\sg}$ as well as the non-linear power spectrum parameter $\alphaNL$, and to make the notation more compact we will define scaling parameters about the fiducial model: $\alphaOmega = \OmegaHI / \OmegaHI^\fid(\zeff)$, $\alphah = \bHI / \bHI^\fid(\zeff)$ and $\alphag = \bg / \bg^\fid(\zeff)$. With this we can rewrite <ref> as

\begin{equation}
\label{eq:crossps_sep}
P_{\sHI,\sg}(\alphaOmega, \alphah, \alphag, \Mten) =
\alphaOmega (\alphah \alphag \, P_{hg} + \alphah \, P_{hv} + \alphag \, P_{vg} + P_{vv}) + \Mten \, P_\mathrm{sn}
\end{equation}

where we have left the $k$, $\mu$ and $z$ dependence implicit. The power spectrum terms on the right hand side are

\begin{align}
P_{hg} & = \Tb^\fid(z_1) \bHI^\fid(z_1) \bg^\fid(z_2) P \\
P_{hv} & = \Tb^\fid(z_1) \bHI^\fid(z_1) f(z_2) \mu^2 P \\
P_{vg} & = \Tb^\fid(z_1) \bg^\fid(z_2) f(z_1) \mu^2 P \\
P_{vv} & = \Tb^\fid(z_1) f(z_1) f(z_2) \mu^4 P \\
P_\mathrm{sn} & = C_\sHI(z_1) 10^{10} \: \si{\Msolar}
\end{align}

\begin{equation}
\label{eq:P_sep}
P = \ls \PL(k) + \alphaNL \lp \PNL(k) - \PL(k) \rp \rs \DFoG{\sHI}(\alphaFoG{\sHI} \, k \mu; z_1) \: \DFoG{\sg}(\alphaFoG{\sg} \, k \mu; z_2) \; .
\end{equation}

The linearity of the simulation and analysis procedure means that the stack signal should be separable into distinct terms like <ref>. If we write the template generated by the given parameters as $s(\alphaOmega, \alphah, \alphag, \Mten)$, where we make the dependence on $\Delta\nu$ implicit, we find

\begin{equation}
\label{eq:template_sep}
s(\alphaOmega, \alphah, \alphag, \Mten) = \alphaOmega (\alphah \alphag \, s_{hg} + \alphah \, s_{hv} + \alphag \, s_{vg} + s_{vv}) + \Mten \, s_\mathrm{sn}
\end{equation}

where each $s_{xy}$ term is the stack signal corresponding to cross-power spectrum term $P_{xy}$ as defined above. However, as we cannot directly propagate a cross-power spectrum into a stack signal, we must determine these terms indirectly. This can be done by running simulations through with specific $\alpha$ parameters that generate known linear combinations of the $s_{xy}$. By choosing these simulated parameters judiciously we can easily invert these combinations to generate the individual $s_{xy}$ terms. One such choice is

\begin{align}
s_{hg} & = s(1, 1, 1, 0) - s(1, 1, 0, 0) - s(1, 0, 1, 0) + s(1, 0, 0, 0) \\
s_{hv} & = s(1, 1, 0, 0) - s(1, 0, 0, 0) \\
s_{vg} & = s(1, 0, 1, 0) - s(1, 0, 0, 0) \\
s_{vv} & = s(1, 0, 0, 0) \\
s_\mathrm{sn} & = s(0, 0, 0, 1) \; .
\end{align}

Each of the five unique combinations of parameters passed to $s(\alphaOmega, \alphah \alphag, \Mten)$ in the equation above requires a separate simulation to determine, but after that, we can use these modes and <ref> to determine the stacked template for any combination of parameters.

This scheme allows us to exactly treat the effects of the three linear parameters and the shot-noise contribution on the template. Incorporating the effect of the non-linear power spectrum shape is straightforward as the parameterisation used for $\alphaNL$ means that the output stack signal is a simple linear mixing of two terms, as it is in the cross power spectrum <ref>. This means that four more simulations can be used to generate templates at any $\alphaNL$. However, the template is a non-linear function of the Fingers of God parameters and so cannot be exactly generated in a finite number of modes.

To account for this we start by noting that the effect of the Fingers of God treatment we use (see sec:fingerofgod), at constant time and in comoving distance, is a convolution of the underlying fields along the line of sight. In a narrow enough interval in redshift, such that we could ignore evolutionary effects and the constant frequency spacing of our measurements maps to a constant separation in comoving distance, this effect commutes with the stacking and we could apply it via convolving a post-simulation $s_{xy}$ template to the desired $\alphaFoG{x}$ value, rather than needing to incorporate it into the simulations directly. However, as our sources are located over wide redshift intervals, evolution of the cosmological fields as well as the pairwise velocity dispersion $\sigma_P$ cannot be neglected, and in addition the RFI masking, redshift-dependent source number density, and sensitivity further break the stationarity of the radial axis.

However, even if there is no exact mapping from the Fingers of God effects into a convolution on the template modes, we can still attempt to find an effective one that is accurate in the vicinity of the fiducial Fingers of God parameters $\alphaFoG{x} = 1$. To do this we use use a transfer function in delay, $\tau$, the Fourier conjugate of the frequency separation $\Delta\nu$, of the form

\begin{equation}
D^\mathrm{eff}_{x}(\tau, \alphaFoG{x}) = \frac{2 + \tau^2 \sigma_\mathrm{eff}^2}{2 + \alphaFoG{x} \tau^2 \sigma_\mathrm{eff}^2} \; ,
\end{equation}

that will be applied to the templates with the fiducial Fingers of God strength.
To motivate this choice, we note that within a short redshift interval,

\begin{equation}
\tau \approx - \frac{1}{2 \pi \nuHI} \frac{c}{H(z)} (1 + z)^2 k_\parallel
\end{equation}

and so if we set

\begin{equation}
\sigma_\mathrm{eff} = 2 \pi \nu_0 \frac{H(z)}{c (1 + z)^2} \sigma_{\rm P}\
\end{equation}

and then apply the transfer function above to the template, the numerator would effectively undo the Lorentzian Fingers of God model with the fiducial $\alphaFoG{x} = 1$ and the denominator would reapply it with the desired $\alphaFoG{x}$. Thus we would have transformed the template mode from the fiducial to the the desired $\alphaFoG{x}$ parameter. To take into account the wide redshift range and non-stationarity, we estimate an effective smoothing width $\sigma_\mathrm{eff}$ by finding the value which minimises the template error at a higher $\alphaFoG{x} = 1.2$ compared to an exact simulation at the same value. This effective convolution approach is applicable over a wide range of values of $\alphaFoG{x}$, with errors at $\alphaFoG{x} = 0$ or $3$ of $\lesssim 0.5 \%$ and much smaller around the pivot $\alphaFoG{x} = 1$. Computationally this requires an additional eight simulations, one for each of the four $(\alphah$, $\alphag)$ combinations with a perturbed value of $\alphaFoG{\sHI} = 1.2$, and an additional four with $\alphaFoG{\sg} = 1.2$

The final effect we need to apply is the frequency shift $\Dnu$. This is performed in Fourier space by phase rotating the delay transform of the template.

[Ahumada et al., 2020]
Ahumada, R., Prieto, C. A., Almeida, A., et al. 2020, , 249, 3,

[Alam et al., 2020]
Alam, S., Peacock, J. A., Kraljic, K., Ross, A. J., & Comparat, J.
2020, , 497, 581, 10.1093/mnras/staa1956

[Alam et al., 2021]
Alam, S., Aubert, M., Avila, S., et al. 2021, , 103, 083533,

[Albareti et al., 2017]
Albareti, F. D., Allende Prieto, C., Almeida, A., et al. 2017, ,
233, 25, 10.3847/1538-4365/aa8992

[Anderson et al., 2018]
Anderson, C. J., Luciw, N. J., Li, Y. C., et al. 2018, , 476,
3382, 10.1093/mnras/sty346

[Ando et al., 2019]
Ando, R., Nishizawa, A. J., Hasegawa, K., Shimizu, I., & Nagamine, K. 2019,
Monthly Notices of the Royal Astronomical Society, 484, 5389–5399,

[Baars et al., 1977]
Baars, J. W. M., Genzel, R., Pauliny-Toth, I. I. K., & Witzel, A.
1977, , 500, 135

[Battye et al., 2004]
Battye, R. A., Davies, R. D., & Weller, J. 2004, , 355, 1339,

[Behnel et al., 2011]
Behnel, S., Bradshaw, R., Citro, C., et al. 2011, Computing in Science
Engineering, 13, 31 , 10.1109/MCSE.2010.118

[Bera et al., 2019]
Bera, A., Kanekar, N., Chengalur, J. N., & Bagla, J. S. 2019, ,
882, L7, 10.3847/2041-8213/ab3656

[Blanton et al., 2017]
Blanton, M. R., Bershady, M. A., Abolfathi, B., et al. 2017, , 154,
28, 10.3847/1538-3881/aa7567

[Bull et al., 2015]
Bull, P., Ferreira, P. G., Patel, P., & Santos, M. G. 2015, , 803,
21, 10.1088/0004-637X/803/1/21

[Calderhead & Girolami, 2009]
Calderhead, B., & Girolami, M. 2009, Computational Statistics & Data
Analysis, 53, 4028–4045, 10.1016/j.csda.2009.07.025

[Castorina & Villaescusa-Navarro, 2017]
Castorina, E., & Villaescusa-Navarro, F. 2017, , 471, 1788,

[Chakraborty et al., 2021]
Chakraborty, A., Datta, A., Roy, N., et al. 2021, , 907, L7,

[Chang et al., 2010]
Chang, T.-C., Pen, U.-L., Bandura, K., & Peterson, J. B. 2010, ,
466, 463, 10.1038/nature09187

[Chang et al., 2008]
Chang, T.-C., Pen, U.-L., Peterson, J. B., & McDonald, P. 2008, ,
100, 091303, 10.1103/PhysRevLett.100.091303

[Chen et al., 2021]
Chen, Q., Meyer, M., Popping, A., et al. 2021a, ,
508, 2758, 10.1093/mnras/stab2810

[Chen et al., 2021]
Chen, S.-F., Vlah, Z., & White, M. 2021b, arXiv e-prints,

[20]
CHIME Collaboration. In Prep.

[CHIME Collaboration et al., 2022]
CHIME Collaboration, Amiri, M., Bandura, K., et al. 2022a,
arXiv e-prints, arXiv:2201.07869.

[CHIME Collaboration et al., 2022]
—. 2022b, arXiv e-prints, arXiv:2201.11822.

[CHIME/FRB Collaboration et al., 2018]
CHIME/FRB Collaboration, Amiri, M., Bandura, K., et al. 2018, ,
863, 48, 10.3847/1538-4357/aad188

[CHIME/Pulsar Collaboration et al., 2021]
CHIME/Pulsar Collaboration, Amiri, M., Bandura, K. M., et al. 2021,
, 255, 5, 10.3847/1538-4365/abfdcb

[Chowdhury et al., 2020]
Chowdhury, A., Kanekar, N., Chengalur, J. N., Sethi, S., &
Dwarakanath, K. S. 2020, , 586, 369, 10.1038/s41586-020-2794-7

[Cohen et al., 2007]
Cohen, A. S., Lane, W. M., Cotton, W. D., et al. 2007, The Astronomical
Journal, 134, 1245, 10.1086/520719

[Collette et al., 2021]
Collette, A., Kluyver, T., Caswell, T. A., et al. 2021, h5py/h5py: 3.5.0,
3.5.0,  Zenodo, 10.5281/zenodo.5585380

[Condon et al., 1998]
Condon, J. J., Cotton, W. D., Greisen, E. W., et al. 1998, , 115,
1693, 10.1086/300337

[Cosmic Visions 21 cm Collaboration et al., 2018]
Cosmic Visions 21 cm Collaboration, Ansari, R., Arena, E. J., et al.
2018, arXiv e-prints, arXiv:1810.09572.

[Crichton et al., 2021]
Crichton, D., Aich, M., Amara, A., et al. 2021, arXiv e-prints,

[Crighton et al., 2015]
Crighton, N. H. M., Murphy, M. T., Prochaska, J. X., et al. 2015,
, 452, 217, 10.1093/mnras/stv1182

[Dalcin & Fang, 2021]
Dalcin, L., & Fang, Y.-L. L. 2021, Computing in Science Engineering, 23, 47,

[d'Amico et al., 2020]
d'Amico, G., Gleyzes, J., Kokron, N., et al. 2020, , 2020, 005,

[Dawson et al., 2016]
Dawson, K. S., Kneib, J.-P., Percival, W. J., et al. 2016, , 151,
44, 10.3847/0004-6256/151/2/44

[de Mattia et al., 2021]
de Mattia, A., Ruhlmann-Kleider, V., Raichoor, A., et al. 2021, ,
501, 5616, 10.1093/mnras/staa3891

[DeBoer et al., 2017]
DeBoer, D. R., Parsons, A. R., Aguirre, J. E., et al. 2017, , 129,
045001, 10.1088/1538-3873/129/974/045001

[Delhaize et al., 2013]
Delhaize, J., Meyer, M. J., Staveley-Smith, L., & Boyle, B. J. 2013,
, 433, 1398, 10.1093/mnras/stt810

[Dey et al., 2019]
Dey, A., Schlegel, D. J., Lang, D., et al. 2019, , 157, 168,

[Eriksen et al., 2004]
Eriksen, H. K., O'Dwyer, I. J., Jewell, J. B., et al. 2004, , 155,
227, 10.1086/425219

[Ewall-Wice et al., 2021]
Ewall-Wice, A., Kern, N., Dillon, J. S., et al. 2021, , 500,
5195, 10.1093/mnras/staa3293

[Foreman-Mackey et al., 2013]
Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, ,
125, 306, 10.1086/670067

[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

[Gelman & Meng, 1998]
Gelman, A., & Meng, X.-L. 1998, Statistical Science, 13, 163 ,

[Geman & Geman, 1984]
Geman, S., & Geman, D. 1984, IEEE Transactions on Pattern Analysis and Machine
Intelligence, PAMI-6, 721, 10.1109/TPAMI.1984.4767596

[Gil-Marín et al., 2020]
Gil-Marín, H., Bautista, J. E., Paviot, R., et al. 2020, ,
498, 2492, 10.1093/mnras/staa2455

[Gordon & Trotta, 2007]
Gordon, C., & Trotta, R. 2007, , 382, 1859,

[Górski et al., 2005]
Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, , 622,
759, 10.1086/427976

[Gould, 1994]
Gould, R. J. 1994, , 423, 522, 10.1086/173830

[Gregory et al., 1996]
Gregory, P. C., Scott, W. K., Douglas, K., & Condon, J. J. 1996,
, 103, 427, 10.1086/192282

[Gunn et al., 2006]
Gunn, J. E., Siegmund, W. A., Mannery, E. J., et al. 2006, , 131,
2332, 10.1086/500975

[Guo et al., 2020]
Guo, H., Jones, M. G., Haynes, M. P., & Fu, J. 2020, , 894, 92,

[Hagberg et al., 2008]
Hagberg, A. A., Schult, D. A., & Swart, P. J. 2008, in Proceedings of the 7th
Python in Science Conference, ed. G. Varoquaux, T. Vaught, & J. Millman,
Pasadena, CA USA, 11 – 15

[Handley & Lemos, 2019]
Handley, W., & Lemos, P. 2019, , 100, 023512,

[Harris et al., 2020]
Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585,
357, 10.1038/s41586-020-2649-2

[Hu et al., 2020]
Hu, W., Catinella, B., Cortese, L., et al. 2020, , 493, 1587,

[Hu et al., 2019]
Hu, W., Hoppmann, L., Staveley-Smith, L., et al. 2019, , 489,
1619, 10.1093/mnras/stz2038

[Hunter, 2007]
Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90,

[Ivanov et al., 2020]
Ivanov, M. M., Simonović, M., & Zaldarriaga, M. 2020, , 2020,
042, 10.1088/1475-7516/2020/05/042

[Jackson, 1972]
Jackson, J. C. 1972, Mon. Not. Roy. Astron. Soc., 156, 1P,

[Jeffreys, 1961]
Jeffreys, H. 1961, Theory of Probability, 3rd edn. (Oxford, England: Oxford)

[Jewell et al., 2004]
Jewell, J., Levin, S., & Anderson, C. H. 2004, , 609, 1,

[Jones et al., 2018]
Jones, M. G., Haynes, M. P., Giovanelli, R., & Moorman, C. 2018, Monthly
Notices of the Royal Astronomical Society, 477, 2,

[Kaiser, 1987]
Kaiser, N. 1987, Mon. Not. Roy. Astron. Soc., 227, 1

[Kanekar et al., 2016]
Kanekar, N., Sethi, S., & Dwarakanath, K. S. 2016, , 818, L28,

[Kaplan, 2005]
Kaplan, G. H. 2005, U.S. Naval Observatory Circulars, 179.

[Karamanis & Beutler, 2021]
Karamanis, M., & Beutler, F. 2021, arXiv e-prints, arXiv:2106.06331.

[Kass & Raftery, 1995]
Kass, R. E., & Raftery, A. E. 1995, Journal of the American Statistical
Association, 90, 773, 10.1080/01621459.1995.10476572

[Kokron et al., 2021]
Kokron, N., DeRose, J., Chen, S.-F., White, M., & Wechsler, R. H.
2021, , 505, 1422, 10.1093/mnras/stab1358

[Kovetz et al., 2019]
Kovetz, E., Breysse, P. C., Lidz, A., et al. 2019, , 51, 101.

[Lah et al., 2007]
Lah, P., Chengalur, J. N., Briggs, F. H., et al. 2007, , 376,
1357, 10.1111/j.1365-2966.2007.11540.x

[Landecker et al., 2000]
Landecker, T. L., Dewdney, P. E., Burgess, T. A., et al. 2000, ,
145, 509, 10.1051/aas:2000257

[Lang et al., 2016]
Lang, D., Hogg, D. W., & Schlegel, D. J. 2016, , 151, 36,

[Laurent et al., 2017]
Laurent, P., Eftekharzadeh, S., Le Goff, J.-M., et al. 2017, ,
2017, 017, 10.1088/1475-7516/2017/07/017

[Leifer, C., 2021]
Leifer, C. 2021, peewee/peewee, latest,  Github.

[Lewis, 2019]
Lewis, A. 2019.

[Lewis et al., 2000]
Lewis, A., Challinor, A., & Lasenby, A. 2000, , 538, 473,

[Li et al., 2020]
Li, J., Zuo, S., Wu, F., et al. 2020, Science China Physics, Mechanics,
and Astronomy, 63, 129862, 10.1007/s11433-020-1594-8

[Li et al., 2021]
Li, L.-C., Staveley-Smith, L., & Rhee, J. 2021, Research in Astronomy
and Astrophysics, 21, 030, 10.1088/1674-4527/21/2/30

[Liu & Tegmark, 2011]
Liu, A., & Tegmark, M. 2011, , 83, 103006,

[Lyke et al., 2020]
Lyke, B. W., Higley, A. N., McLane, J. N., et al. 2020, , 250, 8,

[Masui et al., 2015]
Masui, K., Amiri, M., Connor, L., et al. 2015, Astronomy and Computing,
12, 181, 10.1016/j.ascom.2015.07.002

[Masui et al., 2017]
Masui, K. W., Shaw, J. R., Ng, C., et al. 2017, arXiv e-prints,

[Masui et al., 2013]
Masui, K. W., Switzer, E. R., Banavar, N., et al. 2013, , 763,
L20, 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,

[Mead et al., 2021]
Mead, A. J., Brieden, S., Tröster, T., & Heymans, C. 2021, ,
502, 1401, 10.1093/mnras/stab082

[Merson et al., 2019]
Merson, A., Smith, A., Benson, A., Wang, Y., & Baugh, C. 2019,
, 486, 5737, 10.1093/mnras/stz1204

[Oh & Mack, 2003]
Oh, S. P., & Mack, K. J. 2003, , 346, 871,

[Paciga et al., 2013]
Paciga, G., Albert, J. G., Bandura, K., et al. 2013, , 433, 639,

[Padmanabhan, 2021]
Padmanabhan, H. 2021, arXiv e-prints, arXiv:2109.00003.

[pandas development team, 2020]
pandas development team, T. 2020, pandas-dev/pandas: Pandas, latest,  Zenodo,

[Parsons et al., 2014]
Parsons, A. R., Liu, A., Aguirre, J. E., et al. 2014, , 788, 106,

[Pen et al., 2009]
Pen, U.-L., Staveley-Smith, L., Peterson, J. B., & Chang, T.-C. 2009,
, 394, L6, 10.1111/j.1745-3933.2008.00581.x

[Perley & Butler, 2017]
Perley, R. A., & Butler, B. J. 2017, , 230, 7,

[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.

[Petit & Luzum, 2010]
Petit, G., & Luzum, B. 2010, IERS Technical Note, 36, 1.

[Rafiei-Ravandi et al., 2021]
Rafiei-Ravandi, M., Smith, K. M., Li, D., et al. 2021, , 922, 42,

[Raichoor et al., 2021]
Raichoor, A., de Mattia, A., Ross, A. J., et al. 2021, , 500,
3254, 10.1093/mnras/staa3336

[Rao et al., 2017]
Rao, S. M., Turnshek, D. A., Sardane, G. M., & Monier, E. M. 2017, Monthly
Notices of the Royal Astronomical Society, 471, 3428–3442,

[Rengelink et al., 1997]
Rengelink, R. B., Tang, Y., de Bruyn, A. G., et al. 1997, , 124,
259, 10.1051/aas:1997358

[Rhee et al., 2018]
Rhee, J., Lah, P., Briggs, F. H., et al. 2018, , 473, 1879,

[Rhee et al., 2016]
Rhee, J., Lah, P., Chengalur, J. N., Briggs, F. H., & Colless, M.
2016, , 460, 2675, 10.1093/mnras/stw1097

[Rhee et al., 2013]
Rhee, J., Zwaan, M. A., Briggs, F. H., et al. 2013, , 435, 2693,

[Rhodes, 2019]
Rhodes, B. 2019, Skyfield: High precision research-grade positions for
planets and Earth satellites generator.

[Richards et al., 2011]
Richards, G. T., Kruczek, N. E., Gallagher, S. C., et al. 2011, ,
141, 167, 10.1088/0004-6256/141/5/167

[Ross et al., 2020]
Ross, A. J., Bautista, J., Tojeiro, R., et al. 2020, , 498, 2354,

[Santos et al., 2005]
Santos, M. G., Cooray, A., & Knox, L. 2005, The Astrophysical Journal, 625,
575, 10.1086/429857

[Sarkar & Bharadwaj, 2019]
Sarkar, D., & Bharadwaj, S. 2019, , 487, 5666,

[Scoccimarro, 2004]
Scoccimarro, R. 2004, , 70, 083007, 10.1103/PhysRevD.70.083007

[Shaver et al., 1999]
Shaver, P. A., Windhorst, R. A., Madau, P., & de Bruyn, A. G. 1999,
, 345, 380.

[Shaw et al., 2020]
Shaw, J. R., Foreman, S., Nitsche, R., et al. 2020a,
radiocosmology/driftscan: 20.5.0, v20.5.0,  Zenodo,

[Shaw et al., 2014]
Shaw, J. R., Sigurdson, K., Pen, U.-L., Stebbins, A., & Sitwell, M.
2014, , 781, 57, 10.1088/0004-637X/781/2/57

[Shaw et al., 2015]
Shaw, J. R., Sigurdson, K., Sitwell, M., Stebbins, A., & Pen, U.-L.
2015, , 91, 083514, 10.1103/PhysRevD.91.083514

[Shaw et al., 2020]
Shaw, J. R., Masui, K., Nitsche, R., et al. 2020b,
radiocosmology/caput: 20.10.0, v20.10.0,  Zenodo,

[Shaw et al., 2020]
Shaw, J. R., Masui, K., Hincks, A. D., et al. 2020c,
chime-experiment/ch_pipeline: 20.10.0, v20.10.0,  Zenodo,

[Shaw et al., 2020]
Shaw, R., Nitsche, R., Foreman, S., & Kefala, A. 2020d,
radiocosmology/cora: 20.5.0, v20.5.0,  Zenodo, 10.5281/zenodo.5949195

[Shaw et al., 2020]
Shaw, R., Nitsche, R., Siegel, S. R., et al. 2020e,
radiocosmology/draco: 20.10.0, v20.10.0,  Zenodo,

[Shen et al., 2016]
Shen, Y., Brandt, W. N., Richards, G. T., et al. 2016, , 831, 7,

[Skilling, 2006]
Skilling, J. 2006, Bayesian Analysis, 1, 833 , 10.1214/06-BA127

[Smee et al., 2013]
Smee, S. A., Gunn, J. E., Uomoto, A., et al. 2013, , 146, 32,

[Smith et al., 2003]
Smith, R. E., Peacock, J. A., Jenkins, A., et al. 2003, , 341,
1311, 10.1046/j.1365-8711.2003.06503.x

[Spinelli et al., 2020]
Spinelli, M., Zoldan, A., De Lucia, G., Xie, L., & Viel, M. 2020, Monthly
Notices of the Royal Astronomical Society, 493, 5434–5455,

[Stevens et al., 2016]
Stevens, A. R. H., Croton, D. J., & Mutch, S. J. 2016, , 461, 859,

[Subrahmanya et al., 2017]
Subrahmanya, C. R., Manoharan, P. K., & Chengalur, J. N. 2017, Journal
of Astrophysics and Astronomy, 38, 10, 10.1007/s12036-017-9430-4

[Switzer et al., 2013]
Switzer, E. R., Masui, K. W., Bandura, K., et al. 2013, , 434,
L46, 10.1093/mnrasl/slt074

[Takahashi et al., 2012]
Takahashi, R., Sato, M., Nishimichi, T., Taruya, A., & Oguri, M.
2012, , 761, 152, 10.1088/0004-637X/761/2/152

[Tamone et al., 2020]
Tamone, A., Raichoor, A., Zhao, C., et al. 2020, , 499, 5527,

[The HDF Group, 1997]
The HDF Group. 1997, Hierarchical Data Format, version 5

[Tinker et al., 2008]
Tinker, J., Kravtsov, A. V., Klypin, A., et al. 2008, , 688, 709,

[Tramonte & Ma, 2020]
Tramonte, D., & Ma, Y.-Z. 2020, , 498, 5916,

[Tramonte et al., 2019]
Tramonte, D., Ma, Y.-Z., Li, Y.-C., & Staveley-Smith, L. 2019, ,
489, 385, 10.1093/mnras/stz2146

[Trott et al., 2016]
Trott, C. M., Pindor, B., Procopio, P., et al. 2016, , 818, 139,

[Vanderlinde et al., 2019]
Vanderlinde, K., Liu, A., Gaensler, B., et al. 2019, in Canadian Long
Range Plan for Astronomy and Astrophysics White Papers, Vol. 2020, 28,

[Vaníček, 1969]
Vaníček, P. 1969, , 4, 387, 10.1007/BF00651344

[Villaescusa-Navarro et al., 2018]
Villaescusa-Navarro, F., Genel, S., Castorina, E., et al. 2018, ,
866, 135, 10.3847/1538-4357/aadba0

[Virtanen et al., 2020]
Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17,
261, 10.1038/s41592-019-0686-2

[Vollmer et al., 2005]
Vollmer, B., Davoust, E., Dubois, P., et al. 2005, , 436, 757,

[Vollmer et al., 2010]
Vollmer, B., Gassmann, B., Derrière, S., et al. 2010, , 511,
A53, 10.1051/0004-6361/200913460

[Wandelt et al., 2004]
Wandelt, B. D., Larson, D. L., & Lakshminarayanan, A. 2004, , 70,
083511, 10.1103/PhysRevD.70.083511

[Wang et al., 2021]
Wang, Z., Chen, Y., Mao, Y., et al. 2021, The Astrophysical Journal, 907, 4,

[Wiese & Fuhr, 2009]
Wiese, W. L., & Fuhr, J. R. 2009, Journal of Physical and Chemical
Reference Data, 38, 565, 10.1063/1.3077727

[Wolz et al., 2017]
Wolz, L., Blake, C., & Wyithe, J. S. B. 2017, , 470, 3220,

[Wolz et al., 2022]
Wolz, L., Pourtsidou, A., Masui, K. W., et al. 2022, , 510, 3495,

[Wu et al., 2021]
Wu, F., Li, J., Zuo, S., et al. 2021, , 506, 3455,

[Wyithe & Loeb, 2008]
Wyithe, J. S. B., & Loeb, A. 2008, , 383, 606,

[Xi et al., 2020]
Xi, H., Staveley-Smith, L., For, B.-Q., et al. 2020, Monthly Notices of the
Royal Astronomical Society, 501, 4550, 10.1093/mnras/staa3931

[Zafar et al., 2013]
Zafar, T., Péroux, C., Popping, A., et al. 2013, Astronomy & Astrophysics,
556, A141, 10.1051/0004-6361/201321154

[Zarrouk et al., 2018]
Zarrouk, P., Burtin, E., Gil-Marín, H., et al. 2018, , 477,
1639, 10.1093/mnras/sty506

[Zhai et al., 2017]
Zhai, Z., Tinker, J. L., Hahn, C., et al. 2017, , 848, 76,

[Zonca et al., 2019]
Zonca, A., Singer, L., Lenz, D., et al. 2019, Journal of Open Source
Software, 4, 1298, 10.21105/joss.01298