21 cm Intensity Mapping with the DSA-2000

Line-intensity mapping is an emerging approach for characterizing large volumes of the Universe (Battye et al. 2004; Chang et al. 2008; Wyithe et al. 2008; Kovetz et al. 2017). It involves measuring the statistical properties of emission across different angular scales, enabling surveys of large swaths of sky with relatively low-resolution instruments. While line-intensity mapping experiments aim to measure a variety of emission lines—including of CO, C ii, and Lyα—of particular interest is the narrow hyperfine transition line of hydrogen, called "21 cm emission" due to its rest-frame wavelength of 21 cm.

21 cm emission is the principal tracer of neutral hydrogen gas (H i) throughout the Universe, and measurement of the signal across redshift would constrain cosmological evolution and galactic dynamics. As the signal is quite faint and contaminated by bright foreground emission, it is most accessible at low redshifts. In an era of large galaxy surveys—such as the Sloan Digital Sky Survey (York et al. 2000; Blanton et al. 2017), the extended Baryon Oscillation Spectroscopic Survey (Dawson et al. 2016), the Dark Energy Spectroscopic Instrument (DESI; DESI Collaboration et al. 2016; Moon et al. 2023), Euclid (Euclid Collaboration et al. 2022), and the Legacy Survey of Space and Time from the Rubin Observatory (Ivezić et al. 2019)—21 cm intensity mapping has the advantage that it is not biased toward massive galaxies and instead probes the density field contribution from faint, unresolved galaxies. Additionally, mapping the narrow 21 cm emission line provides excellent redshift resolution beyond that achievable with galaxy surveys (Bull et al. 2015; Santos et al. 2015; Villaescusa-Navarro et al. 2018; Chen et al. 2021).

In recent years, the near-redshift 21 cm signal has been measured in cross-correlation with galaxy surveys with GBT (Chang et al. 2010; Masui et al. 2013; Wolz et al. 2021), Parkes (Anderson et al. 2018), CHIME (CHIME Collaboration et al. 2023), and MeerKAT (Cunnington et al. 2023); CHIME Collaboration et al. (2024) detected the 21 cm signal in correlation with the Lyα forest. Paul et al. (2023) reported the first-ever detection of the 21 cm autocorrelation power spectrum with the MeerKAT telescope. The detection corresponds to a redshift range of 0.32 ≤ z ≤ 0.44 and power spectrum modes of 0.48 h Mpc⁻¹ ≤ k ≤ 11.21 h Mpc⁻¹.

This paper explores measurement of the 21 cm power spectrum with the DSA-2000, a forthcoming 2000-element radio array to be built in Nevada (Hallinan et al. 2019). The DSA-2000 will operate at frequencies of 0.7–2 GHz and is optimized for fast survey speeds. It will leverage the novel Radio Camera imaging pipeline for real-time imaging, allowing for efficient processing of fully cross-correlated interferometric data. As a multipurpose observatory, the DSA-2000 has wide-reaching science applications, including fast radio burst detection and localization, weak lensing studies (Connor et al. 2022), multimessenger astronomy, pulsar timing with NANOGrav (Wahl et al. 2022), and more. We show that it will be a powerful instrument for near-redshift 21 cm intensity mapping.

Intensity-mapping experiments require both exquisite sensitivity and excellent systematics control. Sensitivity is crucial for detecting the faint signal, and in this paper we explore the expected sensitivity of the DSA-2000 to the 21 cm power spectrum signal. We show that, due to its many antennas and fast survey speeds, the DSA-2000 will have the necessary sensitivity to measure the 21 cm power spectrum across a large range of power spectrum modes with just a few minutes of data. In a forthcoming paper, we predict the resulting constraints on cosmological parameters (N. Mahesh et al. 2024, in preparation).

However, sensitivity alone is not sufficient for intensity mapping, as systematics can quickly overwhelm the signal and prevent separation of the astrophysical foregrounds. For this reason, instrument characterization through calibration is crucial. While several experiments in the field have the needed sensitivity to, in principle, measure the 21 cm power spectrum, the DSA-2000 is a uniquely calibratable instrument, owing once again to its many antennas and its pseudo-random array layout. Furthermore, its imaging capabilities will enable excellent foreground characterization and separation with minimal spectral leakage. This will make the DSA-2000 a powerful addition to the field of low-redshift 21 cm intensity mapping.

The DSA-2000 is a survey instrument that builds upon its pathfinders, the DSA-10 and DSA-110 (Kocz et al. 2019; Ravi et al. 2022). Over the course of 4 month observing seasons, the DSA-2000 will survey the entire sky above a decl. of −30°, delivering images at a spatial resolution of 35 and with a sensitivity of 2 μJy beam⁻¹ (Hallinan et al. 2019). It will operate as a fast survey complement to the Square Kilometre Array Mid and the Next Generation Very Large Array, with approximately 10 times the survey speed of each of those forthcoming instruments.

The DSA-2000 achieves this combination of resolution, sensitivity, and survey speed with an unprecedented 2000 antenna dishes, pictured in Figure 1. This is enabled by two key technological advances. First, innovative ambient-temperature low-noise amplifiers allow for inexpensive dish construction (Weinreb & Shi 2021). Next, advances in graphics processing unit technology in recent years will enable real-time imaging of the array's over 2 million visibilities through a data processing pipeline dubbed the Radio Camera.

Zoom In Zoom Out Reset image size

Each of the DSA-2000 antennas will be constructed from a 5 m hydroformed aluminum dish. The antennas will be fully steerable in both azimuth and elevation, differing from the pathfinder DSA-10 and DSA-110 dishes that are steerable in elevation only. A quad-ridge horn feed at the focus of each dish will allow for wideband dual-polarization measurements and will house the low-noise amplifier described in Weinreb & Shi (2021). Signals will be transported via RF-over-fiber to a central processing building, where they will be digitized, correlated, and imaged.

The Radio Camera data processing pipeline relies on the DSA-2000's exquisite imaging fidelity, owing to its many antennas and well-behaved point-spread function (PSF). The pseudo-random array configuration, pictured in Figure 2, is designed to produce a PSF with very low-amplitude sidelobes, plotted in Figure 3. This obviates the need for computationally costly visibility-based deconvolution (Connor et al. 2022). It also means that the DSA-2000 will have excellent uv coverage for power spectrum analyses such as 21 cm intensity mapping and that it will exhibit minimal spectral structure, reducing a dominant systematic for 21 cm analyses.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

While the DSA-2000 has a myriad of scientific applications, this paper concerns detection of the power spectrum of 21 cm emission. In this section, we define the predicted signal at the redshifts of interest.

The 21 cm power spectrum at low redshift is given by

$$\begin{eqnarray}&&{P}_{21}(k,z)={[{T}_{21}(z)]}^{2}{b}^{2}{P}_{\mathrm{mat}}(k,z),\end{eqnarray} \tag{ 1 }$$

where T₂₁(z) is the mean 21 cm brightness temperature at redshift z, P_mat(k, z) is the matter power spectrum, and b is the bias factor for H i halos. We represent T₂₁(z) as

$$\begin{eqnarray}&&{T}_{21}(z)\approx 0.084\,\mathrm{mK}\displaystyle \frac{{\left(1+z\right)}^{2}h}{\sqrt{{{\rm{\Omega }}}_{m}{\left(1+z\right)}^{3}+{{\rm{\Omega }}}_{{\rm{\Lambda }}}}}\displaystyle \frac{{{\rm{\Omega }}}_{B}}{0.044}\displaystyle \frac{{f}_{{\rm{H}}\,{\rm\small{I}}}(z)}{0.01}\end{eqnarray} \tag{ 2 }$$

(Pober et al. 2013b). Here, f_{H I}(z) is the mass fraction of H i with respect to overall cosmological baryon content. Table 1 presents the cosmological parameter values used throughout this paper.

We calculate the matter power spectrum P_mat(k, z) at redshift z = 0.5 and k = 0–11.8 h Mpc⁻¹ with the Code for Anisotropies in the Microwave Background (CAMB; Lewis et al. 2000). ³ We define the normalized theoretical power spectrum P_theory(k, z), with units mK², as

$$\begin{eqnarray}&&{P}_{\mathrm{theory}}(k,z)=\displaystyle \frac{{k}^{3}}{2{\pi }^{2}}{P}_{21}(k,z).\end{eqnarray} \tag{ 3 }$$

Figures 5, 6, and 8 plot P_theory(k, z) in black. The dashed black line indicates extrapolated values at high k.

In this section, we quantify the sensitivity of the DSA-2000 by comparing the predicted statistical uncertainty of its measurement to the expected 21 cm power spectrum defined above in Section 3. This analysis determines the measurement capabilities of the DSA-2000 in the absence of systematics. Adequate sensitivity is a necessary but not sufficient condition for achieving the measurement, and Section 5 discusses anticipated sources of systematic error.

Our analysis is based on a delay spectrum approach discussed in, for example, Parsons et al. (2012) and Kolopanis et al. (2019), and described in detail below in Section 4.1. This approach offers a simple, analytic estimate of the impact of thermal noise on the 21 cm power spectrum.

We consider three noise sources: thermal noise, sample variance, and shot noise. The first benefits from a low system temperature, T_sys, and can be reduced through long time integrations, whereas the sample variance and shot noise are mitigated by measuring large swaths of the sky.

Our analysis is simplified by omitting the effects of synthesis rotation. Synthesis rotation exploits the Earth's rotation to improve the instrument's uv sampling. The degree of synthesis rotation depends on the survey strategy used, and we neglect it here in order to provide a more general estimate of the instrument's performance. Furthermore, our analysis omits consideration of techniques commonly used in the field to reduce noise bias, such as correlating interleaved time steps or removing baseline autocorrelations (see Morales et al. 2023 for an overview of those techniques). We use the cylindrically averaged 1D power spectrum as our figure of merit, and we therefore do not consider nonisotropic effects, such as the "fingers of God," that break symmetry between modes parallel and perpendicular to the line of sight.

We consider the effect of foreground contamination and use an aggressive foreground avoidance strategy that masks power spectrum modes likely to be affected by foreground emission. However, in Section 6 we note that the DSA-2000 may be able to use a less stringent foreground avoidance approach that recovers some of these masked modes.

The sensitivity of a power spectrum analysis depends on the data volume that contributes to each measurement bin. Our result therefore scales with the chosen redshift and power spectrum bin sizes. Larger bins produce more sensitive measurements but sacrifice measurement resolution. In this section, we use linearly spaced power spectrum bins of widths of either Δk = 0.1 h Mpc⁻¹ (Figures 5, 6, and 8) or Δk = 0.03 h Mpc⁻¹ (Figures 9, 10, and 12). We use the full frequency band of 0.7–1.43 GHz, corresponding to a redshift bin of 0 ≤ z ≤ 1.

All software underlying this analysis is open source and available on GitHub ⁴ (Byrne 2024).

4.1. The Delay Spectrum Analysis

Here, we detail the data analysis pipeline that underlies our thermal sensitivity analysis, presented below in Section 4.2.

Intensity-mapping analyses broadly fall into two categories: delay spectrum analyses and imaging power spectrum analyses (for a detailed comparison of these methods, see Morales et al. 2019 or Liu & Shaw 2020). An imaging power spectrum analysis requires a reconstructed uv plane, created by gridding visibilities and analogous to the angular Fourier transform of the reconstructed sky image. In contrast, the delay spectrum analysis bypasses gridding. Instead, the visibilities are Fourier transformed across frequency, squared, and then binned to produce a power spectrum estimate. This approach has the benefit that it is simple, computationally efficient, and enables analytic error propagation, and this paper therefore uses a delay spectrum analysis approach to derive the measurement sensitivity. In practice, however, we expect intensity mapping with the DSA-2000 to use an imaging power spectrum approach, as it better synergizes with the Radio Camera imaging pipeline and does not require storing the raw visibilities. While the noise properties of the imaging and delay spectrum approaches are not identical, we expect that our delay-spectrum-based sensitivity estimate is a reasonable approximation of the noise for any intensity-mapping analysis.

This delay spectrum analysis consists of three steps: Fourier transforming the visibilities (Section 4.1.1), implementing foreground avoidance to remove the contaminating foreground signal (Section 4.1.2), and finally combining the visibilities to form the power spectrum estimate (Section 4.1.3). In Section 4.2, we describe how we propagate the expected thermal noise through this pipeline.

4.1.1. Fourier Transform

We represent the visibility formed by correlating antennas i and j as v_ij (f), where f denotes frequency. The first step in the delay spectrum analysis involves performing a discrete Fourier transform across frequency for each visibility, which we define as follows:

$$\begin{eqnarray}&&{\tilde{v}}_{{ij}}(\eta )=\displaystyle \sum _{f}W(f){v}_{{ij}}(f){e}^{-2\pi i\eta f}.\end{eqnarray} \tag{ 4 }$$

Here, η represents delay, the Fourier dual of frequency with units of time, and the tilde (${}^{\mathop{}\limits^{\sim }}$) denotes the Fourier-transformed quantity. W(f) is the antialiasing window function (e.g., Blackman–Harris or Tukey). To preserve power, the window function is normalized such that

$$\begin{eqnarray}&&\displaystyle \sum _{f}{W}^{2}(f)={N}_{\mathrm{freq}},\end{eqnarray} \tag{ 5 }$$

where N_freq is the number of frequency channels.

In this analysis, we neglect the effects of flagging and assume evenly spaced frequency channels, such that a simple discrete Fourier transform suffices. The introduction of frequency-dependent data excision ("flagging"), for mitigation of radio-frequency interference (RFI) and other data contaminants, could result in deleterious spectral mode-mixing (Wilensky et al. 2022), if not handled correctly. At minimum, calculating the delay spectrum following data excision requires an approach that accounts for gaps in the spectral data, such as the Lomb–Scargle periodogram (see Lomb 1976; Scargle 1982, and VanderPlas 2018 for a description of the technique; Jacobs et al. 2016 describes its application to 21 cm cosmology with the Murchison Widefield Array). Beyond that, other algorithms may result in better spectral performance with reduced mode-mixing. These include iterative deconvolution techniques such as the CLEAN algorithm (Högbom 1974; Parsons et al. 2014); Markov Chain Monte Carlo (MCMC) algorithms such as Gibbs sampling (CHIME Collaboration et al. 2023; Kennedy et al. 2023); and frequency inpainting with Gaussian Process Regression (kriging; Mertens et al. 2018; Trott et al. 2020).

4.1.2. Foreground Avoidance

Foreground mitigation strategies for 21 cm analyses fall into two categories: foreground subtraction and foreground avoidance. Foreground subtraction entails modeling and subtracting bright foreground sources, while foreground avoidance removes foreground power by masking the contaminated power spectrum modes. The latter approach works because the foregrounds are inherently spectrally smooth, owing to their synchrotron and free–free emission mechanisms, and therefore inhabit a compact region in power spectrum space. While many analyses in the field employ some combination of foreground subtraction and avoidance, it is infeasible to model the foreground emission with the fidelity necessary to rely on foreground subtraction alone. We therefore adopt a foreground avoidance approach in this analysis.

Spectrally smooth foregrounds inhabit a wedge-shaped region of 2D power spectrum space, dubbed the "foreground wedge" and discussed extensively in the literature (Morales et al. 2012, 2019; Trott et al. 2012; Vedantham et al. 2012; Hazelton et al. 2013; Pober et al. 2013a; Thyagarajan et al. 2013; Liu et al. 2014a, 2014b; Dillon et al. 2015). The extent of the foreground wedge depends on the largest angle at which the instrument detects off-axis foreground emission. This angle is set nominally by the instrument's field of view. However, as instruments may have significant nonzero sensitivity beyond the field of view extent, many analyses in the field take a conservative approach to foreground avoidance and omit wedge modes down to the horizon.

We represent the foreground mask with a binary function M_ij (η) where values of 0 correspond to masked modes. We define this function as

$$\begin{eqnarray}{M}_{{ij}}(\eta )=\left\{\begin{array}{ll}0, & | \eta | -{\rm{\Delta }}\eta /2\leqslant (\sin \omega )\,| {{\boldsymbol{b}}}_{{ij}}| /c\\ 1, & | \eta | -{\rm{\Delta }}\eta /2\gt (\sin \omega )\,| {{\boldsymbol{b}}}_{{ij}}| /c.\end{array}\right.\end{eqnarray} \tag{ 6 }$$

Here, ∣ b _ij ∣ is the length of the baselines formed by antennas i and j in units of meters and c represents the speed of light. $\sin \omega$ sets the extent of the wedge, where ω corresponds to the maximum angle from the pointing center where instrument is sensitive to foregrounds. Setting ω = 90° corresponds to a horizon cut. We assume the delay axis is discretized into evenly spaced bins of width Δη, and η corresponds to the value at the center of the bin. The term Δη/2 ensures that any delay bin containing foreground contamination will be masked.

4.1.3. Binning and Averaging

The next step in the analysis pipeline is combining the Fourier-transformed visibilities into power spectrum bins and averaging their squared magnitudes. Because the DSA-2000 is a nonregular array, visibilities cannot be coherently averaged before squaring. The resulting reconstructed power spectrum is given by

$$\begin{eqnarray}&&\hat{P}(k)=\displaystyle \frac{{\sum }_{{ij}\eta \in k}{M}_{{ij}}(\eta )| {\tilde{v}}_{{ij}}(\eta ){| }^{2}}{{\sum }_{{ij}\eta \in k}{M}_{{ij}}(\eta )},\end{eqnarray} \tag{ 7 }$$

where the ${}^{\hat{}}$ symbol indicates the empirically estimated value. Here, the sum is taken over all values that contribute to power spectrum bin k. M_ij (η) represents the binary masking function described above in Section 4.1.2, so the denominator of this expression simply counts the number of contributing values.

To determine which values contribute to each bin, we must convert the baseline length and delay η into cosmological units. See Appendix A for a description of that conversion and Table 1 for values of cosmological parameters used. In terms of the conversion factors C_∥ and C_⊥ defined in Appendix A, a visibility mode ${\tilde{v}}_{{ij}}(\eta )$ contributes to power spectrum bin k if

$$\begin{eqnarray}&&k-{\rm{\Delta }}k/2\lt \sqrt{{\left({C}_{\perp }| {{\boldsymbol{b}}}_{{ij}}| f/c\right)}^{2}+{C}_{\parallel }^{2}{\eta }^{2}}\lt k+{\rm{\Delta }}k/2.\end{eqnarray} \tag{ 8 }$$

Here, Δk is the bin width, c is the speed of light, ∣ b _ij ∣ is the baseline length, and f is the frequency, which we define as the central frequency of the band.

4.2. Thermal Noise Propagation

We can now propagate the expected thermal noise through the analysis steps outlined above in Section 4.1.

The thermal noise of a visibility is quantified by the root mean square (rms):

$$\begin{eqnarray}&&\mathrm{rms}[{v}_{{ij}}(f)]=\displaystyle \frac{{\lambda }^{2}{T}_{\mathrm{sys}}}{{A}_{e}\sqrt{{\rm{\Delta }}f\tau }},\end{eqnarray} \tag{ 9 }$$

where λ is the observed wavelength, T_sys is the system temperature, Δf is the channel width, τ is the integration time, and A_e is the effective collecting area (Morales & Wyithe 2010). We approximate the effective collecting area as

$$\begin{eqnarray}&&{A}_{e}=\displaystyle \frac{\pi }{4}{D}_{a}^{2}{e}_{a},\end{eqnarray} \tag{ 10 }$$

where D_a is the antenna diameter and e_a is the aperture efficiency. A table of the antenna parameter values used in this paper is provided in Table 2.

Table 2. Instrument Parameter Values Used throughout This Paper

Variable	Description	Value
Tsys	System temperature	25 K
Da	Antenna diameter	5 m
Nants	Number of antennas	2048
Nbls	Number of baselines	2,098,176
Δf	Frequency channel width	130.2 kHz
ea	Aperture efficiency	0.7
${f}_{\min }$	Minimum frequency	0.7 GHz
${f}_{\max }$	Maximum frequency	1.43 GHz
FoV	Field of view at 1.06 GHz	30.0 deg2
ΔΘ	Field-of-view diameter at 1.06 GHz	618

Notes. While the DSA-2000 measures frequencies up to 2 GHz, the maximum frequency used in this paper refers to the 21 cm rest-frame frequency. Here, we define the limits of the FoV at the 5% beam level, where we have modeled the antennas as circular apertures with Airy disk beams.

Download table as:  ASCIITypeset image

If we assume the thermal noise on each visibility is circularly Gaussian distributed, then we can define the thermal noise variance of each of the real and imaginary components of the visibilities (denoted ${\mathfrak{R}}[{v}_{{ij}}(f)]$ and ${\mathfrak{I}}[{v}_{{ij}}(f)]$, respectively) as

$$\begin{eqnarray}\begin{array}{rcl}{\sigma }_{\mathrm{vis}}^{2} & = & {\mathrm{Var}}_{\mathrm{therm}}\left({\mathfrak{R}}[{v}_{{ij}}(f)]\right)={\mathrm{Var}}_{\mathrm{therm}}\left({\mathfrak{I}}[{v}_{{ij}}(f)]\right)\\ & = & \displaystyle \frac{1}{2}{\left(\mathrm{rms}[{v}_{{ij}}(f)]\right)}^{2}.\end{array}\end{eqnarray} \tag{ 11 }$$

Here, Var_therm denotes the thermal noise variance. If we further assume that the noise is independent on each frequency channel and propagate the thermal noise variance through the Fourier transform operation in Equation (4), we get that

$$\begin{eqnarray}&&{\mathrm{Var}}_{\mathrm{therm}}\left({\mathfrak{R}}[{\tilde{v}}_{{ij}}(\eta )]\right)={\mathrm{Var}}_{\mathrm{therm}}\left({\mathfrak{I}}[{\tilde{v}}_{{ij}}(\eta )]\right)={\sigma }_{\mathrm{vis}}^{2}\displaystyle \sum _{f}{W}^{2}(f).\end{eqnarray} \tag{ 12 }$$

From the normalization of the window function in Equation (5), this reduces to

$$\begin{eqnarray}&&{\mathrm{Var}}_{\mathrm{therm}}\left({\mathfrak{R}}[{\tilde{v}}_{{ij}}(\eta )]\right)={\mathrm{Var}}_{\mathrm{therm}}\left({\mathfrak{I}}[{\tilde{v}}_{{ij}}(\eta )]\right)={N}_{\mathrm{freq}}{\sigma }_{\mathrm{vis}}^{2}.\end{eqnarray} \tag{ 13 }$$

We thereby find that the noise on the Fourier-transformed visibilities is independent of the window function used.

The Fourier transform is a linear operation, so it preserves the circular Gaussian noise distribution profile. It follows that the squared quantity has noise variance

$$\begin{eqnarray}&&{\mathrm{Var}}_{\mathrm{therm}}\left(| {\tilde{v}}_{{ij}}(\eta ){| }^{2}\right)=4{N}_{\mathrm{freq}}^{2}{\sigma }_{\mathrm{vis}}^{4}.\end{eqnarray} \tag{ 14 }$$

From Equation (7), and assuming constant noise across all baselines, the estimated power spectrum has noise variance

$$\begin{eqnarray}&&{\mathrm{Var}}_{\mathrm{therm}}\left[\hat{P}(k)\right]=\displaystyle \frac{4{N}_{\mathrm{freq}}^{2}{\sigma }_{\mathrm{vis}}^{4}{\sum }_{{ij}\eta \in k}{M}_{{ij}}^{2}(\eta )}{{\left[{\sum }_{{ij}\eta \in k}{M}_{{ij}}(\eta )\right]}^{2}}.\end{eqnarray} \tag{ 15 }$$

M_ij (η) can take only values of 0 and 1, so this is equivalent to

$$\begin{eqnarray}&&{\mathrm{Var}}_{\mathrm{therm}}\left[\hat{P}(k)\right]=\displaystyle \frac{4{N}_{\mathrm{freq}}^{2}{\sigma }_{\mathrm{vis}}^{4}}{{N}_{\mathrm{vis}}(k)},\end{eqnarray} \tag{ 16 }$$

where the denominator N_vis(k) = ∑_{ij η∈k} M_ij (η) is simply the number of measurements that contribute to bin k. Figure 4 plots the number of samples as a function of k_⊥ and k_∥.

Zoom In Zoom Out Reset image size

Here, we have omitted any treatment of the visibilities' polarization and assumed a per-polarization analysis. Because the DSA-2000 antennas are dual-polarized, and the 21 cm signal is unpolarized, the XX and YY visibilities can be averaged to produce an estimate of the Stokes I signal. ⁵ Assuming independent thermal noise realizations from the two polarizations, the result is a factor-of-two improvement in the variance of the combined visibilities, producing a factor-of-four improvement in the thermal variance of the power spectrum. We therefore define the thermal variance as

$$\begin{eqnarray}&&{\mathrm{Var}}_{\mathrm{therm}}\left[\hat{P}(k)\right]=\displaystyle \frac{{N}_{\mathrm{freq}}^{2}{\sigma }_{\mathrm{vis}}^{4}}{{N}_{\mathrm{vis}}(k)}.\end{eqnarray} \tag{ 17 }$$

Figure 5 plots the predicted thermal noise levels for different integration times. We use linearly spaced k bins of width Δk = 0.1 h Mpc⁻¹ and plot results from two foreground avoidance strategies: excluding all foreground power within the field of view (up to the dashed white line in Figure 4, denoted with circular markers in Figure 5) and excluding modes up to the horizon threshold (the solid white line in Figure 4, denoted with "X" markers in Figure 5).

Zoom In Zoom Out Reset image size

4.3. Sample Variance

The sample variance refers to the variance in the measurement due to the finite volume of space measured. ⁶ We measure finite regions on the sky across a limited redshift range, and in this section we estimate the resulting sample variance.

The first step in estimating the sample variance is to calculate the expected sample variance at each point in 3D power spectrum space, denoted k . We define a signal S( k ) drawn from a complex circular Gaussian distribution with mean zero, and we define the power spectrum as P( k ) = ∣S( k )∣². If the signal variance is

$$\begin{eqnarray}&&{\sigma }^{2}={\mathrm{Var}}_{\mathrm{samp}}\left({\mathfrak{R}}\left[S({\boldsymbol{k}})\right]\right)={\mathrm{Var}}_{\mathrm{samp}}\left({\mathfrak{I}}\left[S({\boldsymbol{k}})\right]\right),\end{eqnarray} \tag{ 18 }$$

then

$$\begin{eqnarray}&&\langle P({\boldsymbol{k}})\rangle =\langle | S({\boldsymbol{k}}){| }^{2}\rangle ={\sigma }^{2},\end{eqnarray} \tag{ 19 }$$

where the brackets 〈〉 denote the expectation value. The variance of the power spectrum is

$$\begin{eqnarray}&&{\mathrm{Var}}_{\mathrm{samp}}\left[P({\boldsymbol{k}})\right]=\langle {P}^{2}({\boldsymbol{k}})\rangle -\langle P({\boldsymbol{k}}){\rangle }^{2}={\sigma }^{4},\end{eqnarray} \tag{ 20 }$$

from which it follows that

$$\begin{eqnarray}&&{\mathrm{Var}}_{\mathrm{samp}}\left[P({\boldsymbol{k}})\right]=\langle P({\boldsymbol{k}}){\rangle }^{2}.\end{eqnarray} \tag{ 21 }$$

We now need to convert this estimate of the point-by-point sample variance of the power spectrum signal into the sample variance of the binned power spectrum measured by our instrument. This requires calculating the characteristic correlation length of the measured signal and the volume of power spectrum space measured.

The correlation length across delay is given by

$$\begin{eqnarray}&&{\rm{\Delta }}\eta =\displaystyle \frac{1}{{f}_{\max }-{f}_{\min }}.\end{eqnarray} \tag{ 22 }$$

For the ${f}_{\min }$ and ${f}_{\max }$ values given in Table 2, this corresponds to Δη = 1.37 × 10⁻⁹ s. Using the coordinate transformation given by Appendix A and the cosmological parameters from Table 1, this equates to a correlation length in k_∥ of 2.74 × 10⁻³ h Mpc⁻¹.

To estimate the uv plane correlation length for a single snapshot image, we use the instrument's expected field of view of 30.0 deg², corresponding to a field-of-view diameter ΔΘ = 618. For each individual observation, the correlation length in the uv plane is

$$\begin{eqnarray}&&{\rm{\Delta }}u=\displaystyle \frac{90^\circ }{{\rm{\Delta }}{\rm{\Theta }}},\end{eqnarray} \tag{ 23 }$$

equivalent to half-wavelength spacing for a horizon-to-horizon observation or approximately 14.6 wavelengths for the DSA-2000's field of view. Transforming into cosmological units (see Appendix A), this becomes a k_⊥ correlation length of 1.53 × 10⁻² h Mpc⁻¹.

We can define a correlation volume in 3D power spectrum space as

$$\begin{eqnarray}&&{\rm{\Delta }}V={\left({\rm{\Delta }}u\right)}^{2}{\rm{\Delta }}\eta .\end{eqnarray} \tag{ 24 }$$

In cosmological units, we estimate that the power spectrum signal is correlated across volumes of 6.43 × 10⁻⁷ h³ Mpc⁻³.

Next, we need to estimate the total volume of 3D power spectrum space measured by our experiment, V(k). Naively, this is equivalent to the volume of a spherical shell, but the calculation is complicated by the finite extent of the delay modes calculated and the foreground wedge, which mean that some power spectrum modes are excluded from the measurement. See Appendix B for the calculation of V(k). Note that here we assume full sampling of the uv plane. Our sample variance calculation does not consider the effect of "holes" in the uv plane, where incomplete sampling would mean that the instrument measures a smaller effective volume of the 3D power spectrum space. However, this is a good approximation for the DSA-2000, which has near-complete uv coverage.

For a given k bin, the number of independent samples that contribute is given by

$$\begin{eqnarray}&&{N}_{\mathrm{samp}}=\displaystyle \frac{V(k){N}_{\mathrm{fields}}}{{\rm{\Delta }}V},\end{eqnarray} \tag{ 25 }$$

where N_fields is the number of independent fields measured. The sample variance is then

$$\begin{eqnarray}&&{\mathrm{Var}}_{\mathrm{samp}}\left[\hat{P}(k)\right]=\displaystyle \frac{{\rm{\Delta }}V}{V(k){N}_{\mathrm{fields}}}{P}_{\mathrm{theory}}^{2}(k,z),\end{eqnarray} \tag{ 26 }$$

where P_theory(k, z) is the predicted power spectrum from Equations (1) and (3).

Figure 6 plots the sample variance for field-of-view foreground avoidance (circles) and the more conservative horizon foreground avoidance ("X"s). We find that the sample variance is reduced for large survey areas, and that the thermal noise, plotted in Figure 5, is the dominant source of measurement noise.

Zoom In Zoom Out Reset image size

4.4. Shot Noise

An additional source of noise in the measured signal is shot noise, which emerges from the discrete nature of the galaxy halos. The density of galaxies in each volume element follows Poisson statistics, and the resulting noise contribution to the 21 cm signal depends on the galaxies' H i masses.

A scale-invariant estimate of the shot noise is given by

$$\begin{eqnarray}&&{P}_{\mathrm{shot}}(z)={V}_{\mathrm{box}}\displaystyle \frac{\left\langle {\sum }_{i}{\left({M}_{{\rm{H}}\,{\rm\small{I}}}^{i}\right)}^{2}\right\rangle }{\langle {\sum }_{i}({M}_{{\rm{H}}\,{\rm\small{I}}}^{i}){\rangle }^{2}},\end{eqnarray} \tag{ 27 }$$

where V_box is the sampled volume, i indexes all H i sources within that volume, and ${M}_{{\rm{H}}\,{\rm\small{I}}}^{i}$ denotes the sources' H i mass (Spinelli et al. 2020; Chen et al. 2021). Evaluating this expression requires simulations of galaxy masses and halo occupation distributions, and the precise value depends on details of the simulation. Results from the Millennium II simulation, described in Boylan-Kolchin et al. (2009) and Spinelli et al. (2020), give P_shot(z = 0) = 46 Mpc³/h³ and P_shot(z = 1) = 61 Mpc³/h³; we infer P_shot(z = 0.5) ≈ 53.5 Mpc³/h³

In units of mK⁴, the shot noise contribution to the signal variance in 3D power spectrum space is

$$\begin{eqnarray}&&{\mathrm{Var}}_{\mathrm{shot}}\left[P({\boldsymbol{k}})\right]={\left(\displaystyle \frac{{k}^{3}}{2{\pi }^{2}}{[{T}_{21}(z)]}^{2}{P}_{\mathrm{shot}}(z)\right)}^{2},\end{eqnarray} \tag{ 28 }$$

where T₂₁(z), the 21 cm brightness temperature, is given by Equation (2). As with the sample variance, averaging independent samples reduces the shot noise. We then get that the variance on the spherically averaged 1D power spectrum is

$$\begin{eqnarray}&&{\mathrm{Var}}_{\mathrm{shot}}\left[\hat{P}(k)\right]=\displaystyle \frac{{\rm{\Delta }}V}{V(k){N}_{\mathrm{fields}}}{\left(\displaystyle \frac{{k}^{3}}{2{\pi }^{2}}{[{T}_{21}(z)]}^{2}{P}_{\mathrm{shot}}(z)\right)}^{2},\end{eqnarray} \tag{ 29 }$$

where, as in Section 4.3, ΔV is the correlation volume (Equation (24)) and V(k) is the total sampled volume in 3D power spectrum space (see Appendix B). N_fields is the number of independent fields measured.

Figure 7 plots the shot noise for different total fields of view and two foreground avoidance strategies. We find that shot noise becomes a dominant source of noise for large k.

Zoom In Zoom Out Reset image size

4.5. Combined Sensitivity Estimate

We define the combined signal variance as

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{Var}\left[\hat{P}({\rm{k}})\right] & = & {\mathrm{Var}}_{\mathrm{therm}}\left[\hat{P}(k)\right]\\ & & +\,{\mathrm{Var}}_{\mathrm{samp}}\left[\hat{P}(k)\right]+{\mathrm{Var}}_{\mathrm{shot}}\left[\hat{P}(k)\right].\end{array}\end{eqnarray} \tag{ 30 }$$

This quantifies the noise on the measurement but does not account for systematic error (see Section 5).

Figure 8 plots the combined variance as 1σ error bars, assuming that a field-of-view foreground wedge exclusion provides sufficient foreground mitigation. The left panel plots the error bars for a 15 minute snapshot observation, covering 30 deg² of the sky. The right panel includes 720 hr of data across 1700 deg², corresponding to a full season of observations at all local sidereal times. We find that, with just 15 minutes of data, the DSA-2000 has the sensitivity to produce a 5σ detection of the 21 cm power spectrum on modes k = 0.32–15.32 h Mpc⁻¹ with Δk = 0.1 h Mpc⁻¹ resolution. With 720 hr of data, a 5σ detection is achievable on modes of k = 0.03–35.12 h Mpc⁻¹.

Zoom In Zoom Out Reset image size

4.6. Baryon Acoustic Oscillations

Of particular interest for low-redshift 21 cm measurement are the Baryon Acoustic Oscillation (BAO) features, detectable as wiggles in the power spectrum signature. In Figures 9, 10, and 12, we have divided the predicted power spectrum signal by a smoothed signal to highlight the BAO features, plotted in black. The smoothed signal is simply the predicted power spectrum convolved with a Gaussian with FWHM of 0.06 h Mpc⁻¹, equal to the approximate peak-to-peak extent of each BAO wiggle. To resolve these features, we require better resolution than the Δk = 0.1 h Mpc⁻¹ resolution used above. We use a resolution of Δk = 0.03 h Mpc⁻¹, equivalent to the peak-to-trough extent of the features, and Figure 9 plots the resulting 1σ error bars for the 720 hr and 1700 deg² survey, assuming foreground wedge avoidance out to the field of view limit. We find that the nominal DSA-2000 array design, plotted in Figure 2, will not have the sensitivity to measure the BAO wiggles. This is because the DSA-2000 has a dearth of short baselines that sample the low-k modes of interest for BAO measurement.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 9 assumes zenith-pointed observations. Off-axis observations can sample lower k modes because of baseline foreshortening. In Figure 10, we plot error bars resulting from a simulation of observations taken at a zenith angle of 60° toward the East. Although the off-axis pointing reduces the noise on these low k modes, it is not sufficient to enable detection of the BAO wiggles.

Measurement of BAO wiggles would require supplementing the nominal DSA-2000 layout from Figure 2 with a more densely packed core subarray. We explore the impact of a 200-antenna core on the sensitivity of measurement of the BAO wiggles. We simulate randomly placed antennas within a radius of 50 m with a minimum antenna separation of 3 m, plotted in Figure 11. The resulting 1σ error bars on the BAO wiggles are plotted in Figure 12. The additional 200 close-packed core antennas enable the low-k sensitivity needed to measure the BAO features.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Section 4, we showed that the DSA-2000 has the necessary sensitivity to measure 21 cm power spectrum across a wide range of angular scales. However, sensitivity alone is not sufficient for enabling 21 cm intensity mapping, as systematics in the analysis can quickly overwhelm the faint signal. These systematics include calibration error, spectral mode-mixing, and contamination from RFI.

Calibration precision is a leading challenge for 21 cm intensity-mapping experiments, as frequency-dependent gain errors as small as one part in 10⁴–10⁵ can swamp the faint 21 cm signal with foreground contamination (Barry et al. 2016). The DSA-2000 will be uniquely calibratable, due to its many antennas and configuration optimized for imaging. An approximately 2000-antenna array produces over 2 million visibility measurements at each frequency and time, leading to an extremely constrained calibration problem. Calibratability is further enhanced by the DSA-2000's pseudo-random configuration. Byrne et al. (2019) shows that nonregular arrays are more resilient to calibration error than regular arrays such as CHIME, CHORD (Vanderlinde et al. 2019), or HERA (De Boer et al. 2014), where antennas are located on a repeating rectangular or hexagonal grid. This is because each visibility produced by a nonregular array captures an independent mode of the sky signal, and by extension, an independent mode of the sky model used in calibration. This holds true even when regular arrays are calibrated with redundant calibration (Byrne et al. 2019; for background on redundant calibration, see Wieringa 1992; Liu et al. 2010, and Dillon et al. 2020). Furthermore, the DSA-2000's excellent imaging capabilities make it well-suited to self-calibration, where iterative calibration and imaging steps improve calibration precision (Wieringa 1992). It could also benefit from unified calibration, a technique akin to self-calibration that is particularly effective for either regular arrays or nonregular arrays with dense uv coverage, such as the DSA-2000 (Byrne et al. 2021).

Precision calibration requires good knowledge of the direction-dependent beam response, as low-level beam errors introduce deleterious spectral calibration error (Orosz et al. 2019; Joseph et al. 2020). As a result, there has been significant investment in modeling the beam response of interferometers' antennas, and in extreme cases the gains must be calculated independently for different directions on the sky using sophisticated direction-dependent calibration algorithms (Bhatnagar et al. 2008; Yatawatta 2012, 2015; van Weeren et al. 2016; Patil et al. 2017; Kenyon et al. 2018). The DSA-2000's uniform, fully steerable dishes are ringed with a screen, pictured in Figure 1. The screen lowers the system noise by mitigating ground spillover and reduces the antennas' mutual coupling, which can lead to nonuniform beam responses and spectral systematics (Kern et al. 2019). Since the dishes are fully steerable, the beam response can be directly measured by rastering across a bright source or employing beam holography techniques (Berger et al. 2016).

While the DSA-2000 is in principle exquisitely calibratable, achieving the optimal calibration precision is complicated by the fact that it will be infeasible to recalibrate the raw data. The Radio Camera imaging approach reduces data in real time, enabling storage of the processed images rather than the raw visibilities. This is a departure from previous 21 cm intensity-mapping experiments, which typically achieve calibration precision through iterative processing of the same visibility data set. For the DSA-2000, the calibration pipeline must achieve the needed precision from the outset. To achieve the DSA-2000's full potential, the Radio Camera's calibration pipeline must incorporate state-of-the-art precision calibration techniques. Pipeline development should be informed by the strict calibration precision requirements for 21 cm intensity mapping.

Apart from its benefit to calibration, the DSA-2000's pseudo-random configuration, optimized for imaging, reduces systematics related to spectral mode-mixing. The DSA-2000 has an extremely well-behaved PSF with low sidelobe amplitudes (see Figure 3). This design was motivated by imaging performance requirements, as it enables accurate image-plane deconvolution (Connor et al. 2022). However, it also benefits spectral analyses such as 21 cm intensity mapping. PSF sidelobes are intrinsically frequency dependent, and this gives rise to spectral leakage of the foregrounds into the foreground wedge, as discussed above in Section 4.1.2. The amount of spectral leakage scales with the amplitude of the PSF sidelobes, meaning that the DSA-2000 would couple less foreground power into the foreground wedge than other 21 cm intensity-mapping experiments. This could be further enhanced with foreground subtraction, leveraging the DSA-2000's high-fidelity imaging capabilities. Further study is needed to explore the predicted amplitude of the DSA-2000's foreground wedge and determine if foreground wedge modes could be recovered, something that has thus far been infeasible in the field.

RFI contamination is a principal limitation for 21 cm analyses (Barry et al. 2019; Wilensky et al. 2020). Even when the affected channels are flagged, RFI flags themselves produce spectral contamination of the signal (Offringa et al. 2019; Chakraborty et al. 2022; Wilensky et al. 2022; Pagano et al. 2023). The state of Nevada has among the lowest population densities in the country, and the DSA-2000 site selection has been driven by the relative RFI environment of the candidate sites. While intensive RFI studies have informed the site selection, further study is required to quantify the expected signal contamination level and its impact on the DSA-2000 science cases. While RFI will certainly pose a challenge for 21 cm intensity mapping, the DSA-2000 has the distinct advantage of enabling deep statistical RFI detection algorithms from its many antennas. This will enhance the capability of state-of-the-art RFI excision algorithms such as AOFlagger (Offringa et al. 2015) and SSINS (Wilensky et al. 2019).

The forthcoming DSA-2000 has a myriad of impactful science applications, and its many antennas, low system temperature, and fast survey speeds provide the sensitivity needed for 21 cm intensity mapping on a wide range of angular scales.

In Section 4, we explore the effect of thermal noise, sample variance, and shot noise on observations made by the DSA-2000. We find that a season of observations, amounting to 720 hr of data, has the necessary sensitivity to produce a 5σ detection of the 21 cm power spectrum at z ≈ 0.5 and scales of k = 0.03–35.12 h Mpc⁻¹, with resolution of Δk = 0.1h Mpc⁻¹. Furthermore, most of these modes are measurable with far less data. A 15 minute snapshot observation has the sensitivity to produce this detection at scales of k = 0.32–15.32 h Mpc⁻¹.

The DSA-2000 is not the only instrument in the field with the necessary sensitivity for 21 cm intensity mapping. However, it features many more antennas than other near-redshift 21 cm intensity-mapping experiments, including CHIME, CHORD, MeerKAT, and HIRAX. The sheer number of antennas, along with its pseudo-random array configuration and fully steerable dishes, makes the DSA-2000 uniquely resilient to the systematic effects that limit other 21 cm intensity-mapping experiments. Namely, the DSA-2000 has the potential to be calibrated more precisely than other instruments, reducing the spectral calibration error that is a dominant systematic for many 21 cm analyses.

Measurement of the 21 cm power spectrum with the DSA-2000 will lead to better constraints on cosmological and astrophysical parameters, which we forecast in a detailed Fisher matrix analysis in N. Mahesh et al. (2024, in preparation). In general terms, the signal most directly traces the H i density fraction, Ω_{H I} b (Bull et al. 2015; Santos et al. 2016; Pourtsidou et al. 2017; Vanderlinde et al. 2019; Liu & Shaw 2020). It can also constrain the Hubble constant H₀, the dark energy density Ω_Λ, the scalar spectral index n_s , the matter fluctuation parameter σ₈, neutrino mass, and the angular diameter distance (Bull et al. 2015; Villaescusa-Navarro et al. 2015; Pourtsidou et al. 2017; Padmanabhan et al. 2019). When combined with other tracers, such as optical galaxy surveys, 21 cm intensity mapping overcomes the surveys' bias toward high-mass galaxies and captures faint, large-scale structure in galaxy filaments. This produces better constraints on astrophysical parameters, including the H i halo mass slope β, the H i mass function cutoff, the scale-dependent optical galaxy bias, and the velocity dispersion σ_v (Padmanabhan et al. 2020; Paul et al. 2023), with implications for understanding the baryon cycle and gas accretion (Sun et al. 2019; Walter et al. 2020). This makes 21 cm intensity mapping with the DSA-2000 highly complementary to upcoming galaxy surveys with DESI, Euclid, and the Rubin Observatory. Intensity-mapping measurements will also synergize with the DSA-2000's own surveys of resolved H i sources.

Measurements of the BAO wiggles is a particularly compelling goal for near-redshift 21 cm intensity mapping, as they would provide strong constraints on cosmic expansion history and H₀ (CHIME Collaboration et al. 2023). The DSA-2000's nominal design (Figure 2) is not optimized for measuring low k, where we expect to find strong BAO signatures. As we show in Section 4.6, supplementing this design with a dense core, such as the 200-antenna random core plotted in Figure 11, would supply additional short baselines and enable measurement of the BAO wiggles, as plotted in Figure 12. The feasibility of constructing such a core remains to be seen.

The emergence of many-element radio arrays, enabled by advances in radio interferometric hardware and computing systems, ushers in a new era for 21 cm intensity mapping. Arrays such as the DSA-2000, with pseudo-random array configurations and well-behaved PSFs, are particularly suited for these measurements. In the coming years, the DSA-2000 could become a valuable addition to the field of near-redshift 21 cm intensity mapping.

This research received support through the generosity of Eric and Wendy Schmidt by recommendation of the Schmidt Sciences program. It is based upon work supported by the National Science Foundation under Award No. 2303952. Part of this research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. Thank you to Jonathan Pober, Ari Cukierman, Yuping Huang, and the JPL Radio Cosmology Journal Club for their valuable input to this paper.

Throughout this paper, we represent the mapping between instrumental and cosmological variables with the parameters C_∥ and C_⊥, defined as

$$\begin{eqnarray}&&{k}_{\parallel }={C}_{\parallel }\eta ,\end{eqnarray} \tag{ A1 }$$

where k_∥ are the line-of-sight power spectrum modes and η is delay, and

$$\begin{eqnarray}&&{{\boldsymbol{k}}}_{\perp }={C}_{\perp }{\boldsymbol{u}},\end{eqnarray} \tag{ A2 }$$

where k _⊥ = (k_x , k_y ) are the power spectrum modes perpendicular to the line of sight and u are coordinates in the uv plane. In this appendix, we define C_∥ and C_⊥.

We transform between η and k_∥ via the relationship

$$\begin{eqnarray}&&{k}_{\parallel }\approx \displaystyle \frac{2\pi {f}_{21}E(z)}{{D}_{{\rm{H}}}{\left(1+z\right)}^{2}}\eta \end{eqnarray} \tag{ A3 }$$

(Hogg 1999). Here, f₂₁ is the 21 cm frequency and D_H is the Hubble distance, defined as D_H = c/H₀ where H₀ is the Hubble constant. The expression E(z) is given by

$$\begin{eqnarray}&&E(z)=\sqrt{{{\rm{\Omega }}}_{m}{\left(1+z\right)}^{3}+{{\rm{\Omega }}}_{k}{\left(1+z\right)}^{2}+{{\rm{\Omega }}}_{{\rm{\Lambda }}}}\end{eqnarray} \tag{ A4 }$$

(Morales & Hewitt 2004).

Perpendicular to the line of sight, we transform from the baseline vector to cosmological modes via the relationship

$$\begin{eqnarray}&&{{\boldsymbol{k}}}_{\perp }=\displaystyle \frac{2\pi }{{D}_{M}(z)}{\boldsymbol{u}},\end{eqnarray} \tag{ A5 }$$

where D_M (z) is the transverse comoving distance (Morales & Hewitt 2004). u is the uv coordinate of a given baseline, given by u = b /λ where b is the vector separation of the antenna pair and λ is the wavelength, which we define at the center of the bandpass.

From Hogg (1999), the transverse comoving distance is given by

$$\begin{eqnarray}{D}_{M}=\left\{\begin{array}{ll}{D}_{{\rm{H}}}\displaystyle \frac{1}{\sqrt{{{\rm{\Omega }}}_{k}}}\sinh (\sqrt{{{\rm{\Omega }}}_{k}}{D}_{C}/{D}_{{\rm{H}}}), & \mathrm{for}\ {{\rm{\Omega }}}_{k}\gt 0\\ {D}_{C}, & \mathrm{for}\ {{\rm{\Omega }}}_{k}=0\\ {D}_{{\rm{H}}}\displaystyle \frac{1}{\sqrt{{{\rm{\Omega }}}_{k}}}\sin (\sqrt{{{\rm{\Omega }}}_{k}}{D}_{C}/{D}_{{\rm{H}}}), & \mathrm{for}\ {{\rm{\Omega }}}_{k}\lt 0\end{array}\right..\end{eqnarray} \tag{ A6 }$$

D_C , the comoving distance, is

$$\begin{eqnarray}&&{D}_{C}={D}_{{\rm{H}}}{\int }_{0}^{z}\displaystyle \frac{{dz}^{\prime} }{E(z^{\prime} )}.\end{eqnarray} \tag{ A7 }$$

This integral is not analytically tractable, so we evaluate it numerically with the scipy.integrate function.

The sample variance calculation presented in Section 4.3 requires an estimate of the total volume of 3D power spectrum space measured by the experiment.

The theoretical maximum volume that contributes to a given power spectrum bin is given by the volume of a spherical shell:

$$\begin{eqnarray}\begin{array}{rcl}V(k) & = & {\int }_{k-{\rm{\Delta }}k/2}^{k+{\rm{\Delta }}k/2}{dr}\,{r}^{2}{\int }_{0}^{\pi }d\phi {\int }_{-1}^{1}d(\cos \theta )\\ & = & \,2\pi {k}^{2}\,{\rm{\Delta }}k+\displaystyle \frac{\pi }{6}{\left({\rm{\Delta }}k\right)}^{3}.\end{array}\end{eqnarray} \tag{ B1 }$$

Here, θ corresponds to the polar angle, ϕ corresponds to the azimuthal angle, and r corresponds to the radius. We integrate the azimuthal angle from 0 to π only to account for the fact that the uv plane is Hermitian, so only half the plane contributes independent values. While this is the maximum volume that could contribute to each measurement, in practice we need to exclude regions of power spectrum space outside the range of delay modes measured, as established by the instrument's frequency resolution, and within the foreground wedge, which is excluded to mitigate foreground contamination (see Section 4.1.2).

The delay range measured is given by $-{\eta }_{\max }\leqslant \eta \leqslant {\eta }_{\max }$, where ${\eta }_{\max }=\tfrac{1}{2{\rm{\Delta }}f}$. Transforming into cosmological units, this means that measured values of k_∥ fulfill

$$\begin{eqnarray}&&| {k}_{\parallel }| \leqslant \displaystyle \frac{{C}_{\parallel }}{2{\rm{\Delta }}f}.\end{eqnarray} \tag{ B2 }$$

Using the instrument's frequency resolution of Δf = 162.5 kHz, as specified in Table 2, and calculating C_∥ from Appendix A, we find that we measure values ∣k_∥∣ ≤ 6.01 h Mpc⁻¹. Converting into spherical coordinates, where ${k}_{\parallel }=r\cos \theta$, we get that

$$\begin{eqnarray}&&| \cos \theta | \leqslant \displaystyle \frac{{C}_{\parallel }}{2{\rm{\Delta }}f\,r}.\end{eqnarray} \tag{ B3 }$$

Next, we can account for the effects of foreground avoidance. In Section 4.1.2, we note that the only measurements that contribute to the power spectrum estimate fulfill

$$\begin{eqnarray}&&| \eta | \gtrsim \sin \omega \,| {{\boldsymbol{b}}}_{{ij}}| /c,\end{eqnarray} \tag{ B4 }$$

where ∣ b _ij ∣ is the baseline length and $\sin \omega$ determines the extent of the foreground wedge cut. This is equivalent to

$$\begin{eqnarray}&&| \eta | \gtrsim \sin \omega \,\sqrt{{u}^{2}+{v}^{2}}/f\end{eqnarray} \tag{ B5 }$$

or

$$\begin{eqnarray}&&| {k}_{\parallel }| \gtrsim \displaystyle \frac{\sin \omega \,{C}_{\parallel }}{{C}_{\perp }f}\,| {k}_{\perp }| ,\end{eqnarray} \tag{ B6 }$$

where C_∥ and C_⊥ are the unit conversion parameters defined in Appendix A. Letting f = 1.06 GHz, the central frequency, the coefficient $\tfrac{{C}_{\parallel }}{{C}_{\perp }f}=1.83$. We now convert this inequality into spherical coordinates, using the relationships ${k}_{\parallel }=r\cos \theta$ and ${k}_{\perp }=r\sin \theta$. We then require that

$$\begin{eqnarray}&&| \cos \theta | \gtrsim \displaystyle \frac{\sin \omega \,{C}_{\parallel }}{{C}_{\perp }f}\sin \theta .\end{eqnarray} \tag{ B7 }$$

Using the identity $\sin \theta =\sqrt{1-{\cos }^{2}\theta }$ and solving for $| \cos \theta |$, we get that

$$\begin{eqnarray}&&| \cos \theta | \gtrsim \displaystyle \frac{\sin \omega \,{C}_{\parallel }}{\sqrt{{\cos }^{2}\omega \,{C}_{\parallel }^{2}+{C}_{\perp }^{2}{f}^{2}}}.\end{eqnarray} \tag{ B8 }$$

We evaluate V(k), the volume of 3D power spectrum space measured, by expressing these inequalities as integration limits. Because the integral is symmetric about k_∥ = 0 (or equivalently, $\cos \theta =0$), we can evaluate the integral for positive $\cos \theta$ and account for the negative region with a factor of 2. When $k+\tfrac{{\rm{\Delta }}k}{2}\leqslant \tfrac{{C}_{\parallel }}{2{\rm{\Delta }}f}$, we need to only consider the foreground wedge cut, as the spherical shell lies fully within the range of delay values measured. We then get that

$$\begin{eqnarray}&&V(k)=2{\int }_{k-{\rm{\Delta }}k/2}^{k+{\rm{\Delta }}k/2}{dr}\,{r}^{2}{\int }_{0}^{\pi }d\phi {\int }_{\tfrac{\sin \omega \,{C}_{\parallel }}{\sqrt{{\cos }^{2}\omega \,{C}_{\parallel }^{2}+{C}_{\perp }^{2}{f}^{2}}}}^{1}d(\cos \theta ).\end{eqnarray} \tag{ B9 }$$

When $k+\tfrac{{\rm{\Delta }}k}{2}\gt \tfrac{{C}_{\parallel }}{2{\rm{\Delta }}f}$, we need to account for both the foreground wedge and the delay range limit. For $k-\tfrac{{\rm{\Delta }}k}{2}\geqslant \tfrac{{C}_{\parallel }}{2{\rm{\Delta }}f}$, the integral becomes

$$\begin{eqnarray}&&V(k)=2{\int }_{k-{\rm{\Delta }}k/2}^{k+{\rm{\Delta }}k/2}{dr}\,{r}^{2}{\int }_{0}^{\pi }d\phi {\int }_{\tfrac{\sin \omega \,{C}_{\parallel }}{\sqrt{{\cos }^{2}\omega \,{C}_{\parallel }^{2}+{C}_{\perp }^{2}{f}^{2}}}}^{\tfrac{{C}_{\parallel }}{2{\rm{\Delta }}f\,r}}d(\cos \theta ),\end{eqnarray} \tag{ B10 }$$

and when $k-\tfrac{{\rm{\Delta }}k}{2}\lt \tfrac{{C}_{\parallel }}{2{\rm{\Delta }}f}\lt k+\tfrac{{\rm{\Delta }}k}{2}$, we evaluate

$$\begin{eqnarray}\begin{array}{rcl}V(k) & = & 2{\int }_{k-{\rm{\Delta }}k/2}^{\tfrac{{C}_{\parallel }}{2{\rm{\Delta }}f}}{dr}\,{r}^{2}{\int }_{0}^{\pi }d\phi {\int }_{\tfrac{\sin \omega \,{C}_{\parallel }}{\sqrt{{\cos }^{2}\omega \,{C}_{\parallel }^{2}+{C}_{\perp }^{2}{f}^{2}}}}^{1}d(\cos \theta )\\ & & +\,2{\int }_{\tfrac{{C}_{\parallel }}{2{\rm{\Delta }}f}}^{k+{\rm{\Delta }}k/2}{dr}\,{r}^{2}{\int }_{0}^{\pi }d\phi {\int }_{\tfrac{\sin \omega \,{C}_{\parallel }}{\sqrt{{\cos }^{2}\omega \,{C}_{\parallel }^{2}+{C}_{\perp }^{2}{f}^{2}}}}^{\tfrac{{C}_{\parallel }}{2{\rm{\Delta }}f\,r}}d(\cos \theta ).\end{array}\end{eqnarray} \tag{ B11 }$$

Evaluating these integrals with Mathematica software, we get that

$$\begin{eqnarray}V(k)=\left\{\begin{array}{ll}\left[2\pi {k}^{2}\,{\rm{\Delta }}k+\displaystyle \frac{\pi }{6}{\left({\rm{\Delta }}k\right)}^{3}\right]\left(1-\displaystyle \frac{\sin \omega \,{C}_{\parallel }}{\sqrt{{\cos }^{2}\omega \,{C}_{\parallel }^{2}+{C}_{\perp }^{2}{f}^{2}}}\right), & k+\displaystyle \frac{{\rm{\Delta }}k}{2}\leqslant \displaystyle \frac{{C}_{\parallel }}{2{\rm{\Delta }}f}\\ \pi k\,{\rm{\Delta }}k\displaystyle \frac{{C}_{\parallel }}{{\rm{\Delta }}f}-\left[2\pi {k}^{2}\,{\rm{\Delta }}k+\displaystyle \frac{\pi }{6}{\left({\rm{\Delta }}k\right)}^{3}\right]\displaystyle \frac{\sin \omega \,{C}_{\parallel }}{\sqrt{{\cos }^{2}\omega \,{C}_{\parallel }^{2}+{C}_{\perp }^{2}{f}^{2}}}, & k-\displaystyle \frac{{\rm{\Delta }}k}{2}\geqslant \displaystyle \frac{{C}_{\parallel }}{2{\rm{\Delta }}f}\\ \displaystyle \frac{\pi }{8}\displaystyle \frac{{C}_{\parallel }}{{\rm{\Delta }}f}{\left(2k+{\rm{\Delta }}k\right)}^{2}-\displaystyle \frac{\sin \omega \,{C}_{\parallel }}{\sqrt{{\cos }^{2}\omega \,{C}_{\parallel }^{2}+{C}_{\perp }^{2}{f}^{2}}}\left[2\pi {k}^{2}{\rm{\Delta }}k+\displaystyle \frac{\pi }{6}{\left({\rm{\Delta }}k\right)}^{3}\right]+\displaystyle \frac{\pi }{12}{\left({\rm{\Delta }}k-2k\right)}^{3}-\displaystyle \frac{\pi }{24}\displaystyle \frac{{C}_{\parallel }^{3}}{{\left({\rm{\Delta }}f\right)}^{3}}, & k-\displaystyle \frac{{\rm{\Delta }}k}{2}\lt \displaystyle \frac{{C}_{\parallel }}{2{\rm{\Delta }}f}\lt k+\displaystyle \frac{{\rm{\Delta }}k}{2}.\end{array}\right.\end{eqnarray} \tag{ B12 }$$