1 Introduction↩︎

The Epoch of Reionization (EoR) is an astrophysically complex era that has yet to be fully explored [1], [2]. Although the midpoint of reionization has been constrained by cosmic microwave background experiments [3]–[5], and the reionization history has also been constrained from various quasar sightline observations [6]–[12], recent observations of high-redshift galaxies—especially the discovery of numerous Lyman-\(\alpha\) emitters (LAEs) deep in reionization—have raised questions about the detailed processes of reionization [13]–[15]. A deeper understanding of the EoR will not only provide rich astrophysical insights into the properties of the first galaxies but also yield significant cosmological implications [16]–[20].

The redshifted 21-cm line from neutral hydrogen provides a direct and comprehensive probe of this otherwise opaque period in cosmic history [21], [22]. The current generation of large radio interferometers, such as the Hydrogen Epoch of Reionization Array (HERA, [23], [24]), has already been setting stringent limits on the 21-cm power spectrum—a statistical measurement of the spatial fluctuations in the 21-cm signal—and has placed important constraints on the properties of the intergalactic medium (IGM) during reionization [25], [26]. Continued observations with current experiments such as HERA and the upcoming Square Kilometre Array (SKA, [27]) will soon reach sufficient sensitivity to detect the 21-cm auto-spectrum [28].

Meanwhile, large ensembles of high-redshift galaxies discovered by ground- and space-based instruments [29]–[35] enable an alternative route to detecting the 21-cm signal via cross-correlation [36]–[38]. Under an “inside-out” reionization scenario, the 21-cm signal anticorrelates with the galaxy field: overdense regions around galaxies reionize first, while the 21-cm signal continues to trace the surrounding neutral regions [39], [40]. Thanks to the high signal-to-noise nature of optical observations, such an anticorrelation could be easier to detect. This is evident from the fact that most cosmological measurements of 21-cm fluctuations in the post-reionization universe have been made through correlating with galaxy surveys [41]–[44]. If detected, this cross-correlation will serve as a crucial sanity check for any future 21-cm auto-spectrum detection.

Here, we study the prospect of detecting a cross-correlation signal by stacking 21-cm image cubes around LAEs. Although many studies have investigated cross-correlating 21-cm data with galaxies [36]–[38], [45]–[48], most focused on correlation functions or cross power spectra. A stacking signal contains less information as full power spectra, detecting such a signal during the EoR still has significant astrophysical implications, because the signal strength is directly sensitive to the global \(\mathrm{HI}\) fraction [49].

In this work, we improve on previous studies by accounting for additional observational and theoretical complexities. We utilize a direct optimal mapping framework [50], which allows us to accurately quantify the statistical properties of these 21-cm image cubes. A realistic foreground-filtering algorithm [51] is applied to quantify signal loss from foreground mitigation. On the theory side, we use the radiation-magneto-hydrodynamics simulations thesan [52]–[55] to derive a signal template. The full radiative transfer calculations adopted within thesan provide a robust physical connection between galaxies and ionized bubbles [40], [56]–[59], including analyses specifically targeting LAE populations [55], [60]. This allows us to account for the complex correlation between the optical properties of LAEs and their surrounding IGM as observed in radio.

This paper is organized as follows. In Sec. 2 we discuss the procedure to generate foreground-filtered 21-cm image cubes and their statistical properties. A theory template for the cross-correlation signal inferred from simulations is given in Sec. 3. The prospects for a cross-correlation detection and its cosmological implications are discussed in Sec. 4. Conclusions are given in Sec. 5.

2 Foreground filtered 21-cm Maps↩︎

In this work, we adopt the direct optimal mapping¹ framework developed in [50] to generate 21-cm image cubes. This formalism is particularly beneficial for our application, as the statistical properties of the images are well understood. Here, we provide a brief overview of direct optimal mapping in Sec. 2.1. We describe how we filter foreground contamination from the data in Sec. 2.2. The statistical properties of these maps are presented in Sec. 2.3.

2.1 Direct Optimal Mapping↩︎

The most natural data product from a radio interferometer is the correlation of voltages measured by any two antennas, i.e., the visibility, as a function of frequency \(\nu\), \[\label{eq:measurement95eq} V(\mathbf{b}_{ij}, \nu) = \int \mathrm{d}\Omega\,I(\hat{\mathbf{s}}, \nu) B_{ij}(\hat{\mathbf{s}}, \nu) \exp\left(-\frac{i2\pi\nu}{c}\mathbf{b}_{ij}\cdot\hat{\mathbf{s}}\right)\, .\tag{1}\] Here, \(i\) and \(j\) are antenna indices; \(I\) is the brightness temperature of the sky; \(B_{ij}\) is the cross-power beam; \(\mathbf{b}_{ij}\) is the baseline vector; and \(\hat{\mathbf{s}}\) is the unit vector on the sky over which we integrate.

We can discretize Eq. 1 and describe the relation between the interferometric data \(\boldsymbol{\mathbf{d}}\) and the sky \(\boldsymbol{\mathbf{m}}\) using a linear system. For a given time and frequency, we write, \[\boldsymbol{\mathbf{d}} = \boldsymbol{\mathbf{A}}\boldsymbol{\mathbf{m}} + \boldsymbol{\mathbf{n}}\, .\] Here, \(\boldsymbol{\mathbf{d}}\) is a vector with a dimension equal to the number of baselines, and \(\boldsymbol{\mathbf{m}}\) is a vector with a dimension equal to the number of discretized sky pixels. \(\boldsymbol{\mathbf{n}}\) represents instrumental noise and can also absorb other uncertainties such as discretization error². The design matrix \(A\) is written as \[\label{eq:amat} A_{mn} \equiv B_{m}(\hat{\mathbf{s}}_n)\exp\left(-\frac{i2\pi\nu}{c}\mathbf{b}_{m}\cdot\hat{\mathbf{s}}_n\right)\, ,\tag{2}\] where the index \(m\) runs over the baseline axis and \(n\) runs over the sky pixel axis. We note that we absorb the area element \(\Delta\Omega\) from Eq. 1 into the sky vector \(\boldsymbol{\mathbf{m}}\). While \(I(\hat{\mathbf{s}}, \nu)\) has a unit of specific intensity (e.g., \(\left[\mathrm{Jy}/\mathrm{Sr}\right]\)), we choose to work with \(\boldsymbol{\mathbf{m}}\) representing the flux density from each pixel (e.g., in units of \(\left[\mathrm{Jy}\right]\)).

A simple but sufficient estimator for the true sky \(\boldsymbol{\mathbf{m}}\) can be formed as \[\label{eq:mapping95from95vis} \boldsymbol{\mathbf{\hat{m}}} \equiv \boldsymbol{\mathbf{D}}\boldsymbol{\mathbf{A}}^\dagger \boldsymbol{\mathbf{N}}^{-1}\boldsymbol{\mathbf{d}}\, ,\tag{3}\] where \(\boldsymbol{\mathbf{D}}\) is a normalization matrix and \(\boldsymbol{\mathbf{N}} \equiv \langle n n^\dagger\rangle\) is the noise covariance. \(\boldsymbol{\mathbf{\hat{m}}}\) satisfies the Fisher–Neyman criterion as long as \(\boldsymbol{\mathbf{D}}\) is invertible [61]. Throughout this work, we assume that the noise covariance in visibility space follows the form \[\label{eq:noise95matrix} N_{ij} = \frac{\sigma^2_\mathrm{rms}(\nu, t)}{n_{\boldsymbol{\mathbf{b}}_i}} \delta_{ij}\, ,\tag{4}\] where \(\sigma_\mathrm{rms}\) is given by the radiometer equation [62], \[\label{eq:radiometer} \sigma_\mathrm{rms} \equiv \frac{T_\mathrm{sys}}{\sqrt{\Delta\nu\Delta t}}\, .\tag{5}\] Here, \(n_{\boldsymbol{\mathbf{b}}_i}\) is the number of redundancy for each baseline group, \(T_\mathrm{sys}\) is the system temperature which is often estimated from the antenna’s auto-correlation, \(\Delta \nu\) is the correlator channel width, and \(\Delta t\) is the correlator integration time. For HERA, \(\Delta \nu = 122\,\mathrm{kHz}\) and \(\Delta t = 9.6\,\mathrm{s}\).

2.2 Imaging Delay-filtered Visibilities↩︎

One of the major barriers for 21-cm cosmology is the presence of bright foreground emission combined with complex instrumental responses. For a radio interferometer, foreground contamination is usually confined to a region of Fourier space known as the foreground wedge [63]–[70]. While there exists a rich literature in explicitly modeling and subtracting the foregrounds, foreground subtraction remains challenging in practice, especially in the presence of uncertainties in instrument response and systematic effects. Therefore, the most conservative method of mitigating foreground is to filter out all modes within the foreground wedge.

Ideally, for the visibility measured by each baseline \(\boldsymbol{\mathbf{b}}\), the foreground contamination should predominantly reside within delay \(|\tau| \leq |\boldsymbol{\mathbf{b}}|/c\), where \(c\) is the speed of light and \(\tau\) is the Fourier dual of the frequency for each baseline. Here, we adopt a foreground filtering method first developed in [51], utilizing a set of basis functions known as the discrete prolate spheroidal sequence (DPSS, [71]). The smooth foreground component in each visibility is removed by fitting these basis functions that are localized in Fourier space (within \(|\tau| \leq |\boldsymbol{\mathbf{b}}|/c\)). For details of this procedure, we refer the reader to 6. For the purpose of this work, we simply treat foreground filtering as a linear operation \(\boldsymbol{\mathbf{\mathcal{O}^\mathrm{fil}}}\), in which the filtered visibility is obtained via \[\label{eq:fg95filt} V^\mathrm{fil}(\mathbf{b}, \nu_i) = \sum_j\mathcal{O}^\mathrm{fil}_{ij}(\mathbf{b}) V(\mathbf{b}, \nu_j)\, .\tag{6}\] This foreground-filtering method has proven successful in real-world applications [26], [72]–[74] and is particularly beneficial in dealing with data with gaps, often introduced by radio frequency interference [75]. While residual foreground may persist due to systematic effects [76]–[78], we defer these to future studies. For the remainder of this work, we assume foreground contamination is completely removed by this delay-filtering procedure.

2.3 Noise Properties↩︎

Because the mapping and foreground filtering procedures outlined in Sec. 2.1 and 2.2 are both linear, a key advantage of this framework is that the noise properties of the image cubes can be easily modeled. Here, we discuss a couple of important statistical properties of our image cubes, including the choice of normalization, the optimal time averaging procedure, and both frequency-frequency and pixel-pixel correlations. We note that while the mathematical framework presented here is generic, the plots shown in this section require specifying an array layout. For all results this section, we use the portion of the HERA array commissioned as of 2022 as an example. We refer the reader to Sec. 4.1 for more details.

First, we discuss the choice of a sensible normalization convention. We note that the (pixel-pixel) noise covariance for each map is \[\label{eq:noise95img} \begin{align} \boldsymbol{\mathbf{N}}^\mathrm{img} &\equiv \langle \boldsymbol{\mathbf{\hat{m}}} \boldsymbol{\mathbf{\hat{m}}}^\dagger\rangle \\ &= \boldsymbol{\mathbf{D}}\boldsymbol{\mathbf{A}}^\dagger \boldsymbol{\mathbf{N}}^{-1} \langle \boldsymbol{\mathbf{n}} \boldsymbol{\mathbf{n}}^\dagger\rangle \boldsymbol{\mathbf{N}}^{-1} \boldsymbol{\mathbf{A}} \boldsymbol{\mathbf{D}}^\dagger \\ &= \frac{1}{2} \boldsymbol{\mathbf{D}}\boldsymbol{\mathbf{A}}^\dagger \boldsymbol{\mathbf{N}}^{-1} \boldsymbol{\mathbf{A}} \boldsymbol{\mathbf{D}}^\dagger\, , \end{align}\tag{7}\] where the factor of \(1/2\) arises because we only take the real part of the image. Without the normalization matrix \(\boldsymbol{\mathbf{D}}\), the noise variance in each sky pixel \(\boldsymbol{\mathbf{\hat{s}}}_i\) is proportional to

\[\label{eq:map95pixel95variance} \begin{align} \left(\boldsymbol{\mathbf{A}}^\dagger \boldsymbol{\mathbf{N}}^{-1} \boldsymbol{\mathbf{A}}\right)_{ii} &= |B(\boldsymbol{\mathbf{\hat{s}}}_i)|^2 \sum_k N^{-1}_{kk} \\ &=\frac{n_\mathrm{bl}|B(\boldsymbol{\mathbf{\hat{s}}}_i)|^2}{\sigma_\mathrm{rms}^2}\, , \end{align}\tag{8}\] where we have assumed cross power beams are the same across all baselines. Here, \(n_\mathrm{bl}\) is the total number of baselines in the array and \(\sigma_\mathrm{rms}\) is given by Eq. 5 . Therefore, a natural choice for \(\boldsymbol{\mathbf{D}}\) so that the noise property is uniform across all pixels is to have \(D_{ij} \propto \delta_{ij}/B(\boldsymbol{\mathbf{\hat{s}}}_i)\).

On the signal side, our image estimator \(\boldsymbol{\mathbf{\hat{m}}}\) relates to the true sky \(\boldsymbol{\mathbf{m}}\) via \[\langle\boldsymbol{\mathbf{\hat{m}}}\rangle = \boldsymbol{\mathbf{D}}\boldsymbol{\mathbf{A}}^\dagger \boldsymbol{\mathbf{N}}^{-1}\boldsymbol{\mathbf{A}}\boldsymbol{\mathbf{m}}\, .\] In particular, in this work we focus on the response of the image estimator to a point source in the sky. While reionization bubbles usually span several megaparsecs, current radio instruments such as HERA do not have the sensitivity to spatially resolve them. Therefore, we treat these bubbles as point sources and focus on obtaining the cross-correlation signal along the frequency direction, where radio interferometers have exquisite resolution. For a point source centered on a sky pixel \(\boldsymbol{\mathbf{\hat{s}}}_i\), we have \[\label{eq:pt95source95signal} \langle{\hat{m}_i}\rangle = D_{ii}\frac{n_\mathrm{bl}|B(\boldsymbol{\mathbf{\hat{s}}}_i)|^2}{\sigma_\mathrm{rms}^2} f_i\, ,\tag{9}\] where \(f_i\) is the flux of the point source and we have assumed the normalization matrix \(\boldsymbol{\mathbf{D}}\) is diagonal. Hence, another choice of \(\boldsymbol{\mathbf{D}}\) so that the peak flux of any point source is preserved is to have \(D_{ij} = \delta_{ij} \sigma^2_\mathrm{rms} / n_\mathrm{bl} / |B(\boldsymbol{\mathbf{\hat{s}}}_i)|^2\). However, as shown in Eq. 8 , this normalization increases the noise variance for pixels farther from the pointing center.

In this work, we choose the following normalization matrix, \[\label{eq:normalization} D_{ij} = \frac{\sigma_\mathrm{rms}^2}{n_\mathrm{bl} |B(\boldsymbol{\mathbf{\hat{s}}}_i)|}\delta_{ij}\, .\tag{10}\] Plugging this into Eq. 9 , such a normalization gives rise to a beam-weighted sky. As discussed above, this choice of normalization also yields a uniform noise variance across all pixels. A clear advantage of this approach is that the mapping between visibilities to maps according to Eq. 3 does not depend on our knowledge of the primary beam, since the beam factor cancels between Eq. 2 and Eq. 10 . In reality, beam modeling can be complex and uncertain. This approach ensures that uncertainties in beam modeling do not propagate into calculations of noise statistics but are instead contained entirely within the signal modeling.

Because the current generation of radio interferometers lacks the sensitivity to detect individual ionized bubbles, it is crucial to correctly average the data along different axes to increase sensitivity. This can be done by averaging image cubes from different times or by stacking different pixels that contain ionized bubbles. However, extra care is required in the presence of pixel-pixel correlations and time-dependent noise properties.

First of all, we consider coherent time averaging by mapping the same pixel with visibilities obtained at different times. For a tracking telescope, the noise variance decreases as \(1/t_\mathrm{obs}\). However, for a transit telescope like HERA, as the source moves across the antenna’s primary beam, the accumulated signal-to-noise ratio is different depending on the location of the source. In other words, a visibility measured at a local sidereal time far from the source’s right ascension contributes little signal to the image. Combining Eq. 9 and 8 , the signal-to-noise ratio of a point source at \(\boldsymbol{\mathbf{\hat{s}}}_i\) with flux \(f_i\) observed at time \(t\) is \[\label{eq:pt95source95SN} \mathcal{S}/\mathcal{N}(t) = \frac{\sqrt{2n_\mathrm{bl}}}{\sigma_\mathrm{rms}}|B(\boldsymbol{\mathbf{\hat{s}}}_i; t)|f_i \, ,\tag{11}\] which is independent from the choice of the normalization matrix. Because of this, although we adopt a normalization that ensures uniform noise levels over time, simply averaging visibilities uniformly across different times does not maximize the signal-to-noise ratio. The signal is diluted as the source moves away from zenith. Here, we consider two averaging schemes \(w_i(t)\): a uniform weighting scheme where \(w_i(t) \propto 1\) and an optimal weighting scheme in which \(w_i(t) \propto |B(\boldsymbol{\mathbf{\hat{s}}}_i; t)|\). To assess how time-averaging helps increase the signal-to-noise ratio, we define an effective integration time \(t_\mathrm{eff}\), given a period of observation \(t_\mathrm{obs}\), to be \[\label{eq:t95eff} \begin{align} &t_\mathrm{eff}(t_\mathrm{obs}) =\\ &\left(\frac{1}{\mathcal{S}/\mathcal{N}(t_0)}\frac{\int_{t_0 - t_\mathrm{obs}/2}^{t_0 + t_\mathrm{obs}/2} w_i(t) \langle{\hat{m}_i}\rangle(t)~\mathrm{d}t}{\sqrt{\int_{t_0 - t_\mathrm{obs}/2}^{t_0 + t_\mathrm{obs}/2} w_i(t)N^\mathrm{img}_{ii}(t)~\mathrm{d}t}}\right)^2\, , \end{align}\tag{12}\] where \(t_0\) is the time where the source transits the zenith. For a tracking telescope with uniform weighting, \(t_\mathrm{eff}(t_\mathrm{obs}) = t_\mathrm{obs}\). The top panel of 1 shows \(t_\mathrm{eff}\) for the two weighting schemes considered above for a HERA-like transit experiment. We see that the signal-to-noise ratio under the optimal weighting scheme (red) saturates after a given time and reaches a higher signal-to-noise ratio compared to the naive uniform weighting scheme (blue). The signal-to-noise ratio for the uniform weighting scheme drops after a certain time as the averaging process contributes more noise than signal. Although an optimal weighting scheme achieves around \(10\%\) higher signal-to-noise ratio from averaging, it requires precise knowledge of the instrument’s primary beam. In this work, we adopt the more conservative and practical uniform-averaging scheme for our forecast. Therefore, for each night of observations, we assume that we can coherently average the image for around 45 minutes to reach an effective average time of roughly 25 minutes. The latter quantity is often referred to as the beam crossing time. We note that for simplicity, here we have assumed that the instrument is stable as a function of time, i.e., the number of available baselines, the system temperature, and the antenna’s primary beam are all not a function of time.

Another way to increase the signal-to-noise ratio for cross-correlation is by stacking 21-cm image cubes around known LAEs. If each line of sight is completely independent, then stacking around \(N\) galaxies reduces the noise by \(\sqrt{N}\). However, as seen from Eq. 7 , there exists non-trivial pixel-pixel correlations due to instrumental response. The noise variance for each pixel is the same, \[\label{eq:pixel95variance} \sigma^2=\frac{\sigma^2_\mathrm{rms}}{2\,n_\mathrm{bl}}\, ,\tag{13}\] thanks to our choice of the normalization matrix. The correlation between two pixels \(i\) and \(j\) is given by \[\label{eq:pixel95corr} \rho_{ij} \equiv \frac{N^\mathrm{img}_{ij}}{\sqrt{N^\mathrm{img}_{ii}N^\mathrm{img}_{jj}}}\, .\tag{14}\] The resulting noise variance after averaging two lines of sight is \[\label{eq:avg95variance} \frac{1+\rho_{ij}}{2}\sigma^2\, .\tag{15}\] The bottom panel of 1 shows the agreement between an analytical calculation and numerical noise simulations. The solid black line is calculated using Eq. 15 , while the dashed red line is obtained from \(500\) realizations of noise simulations. For each noise realization, we generate images at different pointings (constant declination, separated by right ascension) from noise-only visibilities according to Eq. 5 . Each pointing is then time-averaged for \(45\) minutes to reach maximum sensitivity. We see that pixels within the array’s synthesis beam are particularly correlated with each other. This is well characterized by the diffraction limit scale \(1.22\,\lambda_\mathrm{obs} / b_\mathrm{max}\). In this plot, the observed wavelength \(\lambda_\mathrm{obs}\approx 1.67\,\mathrm{meter}\) traces 21-cm lines at \(z\sim7\), and the longest baseline length is at \(b_\mathrm{max}\approx265\,\mathrm{meter}\). Therefore, if we stack \(N\) galaxies that are pairwise separated by at least the diffraction limit, we can safely assume that the noise level decreases by \(\sqrt{N}\).

So far, we have focused on the image properties at one frequency. While images at different frequency bins should be independent, the foreground filtering procedure described in Sec. 2.2 introduces frequency-frequency correlations in each single-baseline visibility. Here, we investigate how optimal mapping propagates visibility-space correlations into image-space correlations.

Following Eq. 6 , after foreground filtering, the visibility measured by a given baseline is correlated between different frequency channels through \[\begin{align} \left(\boldsymbol{\mathbf{C}}^\mathrm{filt}_{\boldsymbol{\mathbf{b}}}\right)_{ij} &\equiv \langle V^\mathrm{fil}(\mathbf{b}, \nu_i) V^\mathrm{fil}(\mathbf{b}, \nu_j)^*\rangle \\ &= \sum_{mn}\mathcal{O}^\mathrm{fil}_{im}(\mathbf{b}) V(\mathbf{b}, \nu_m)\mathcal{O}^\mathrm{fil}_{jn}(\mathbf{b})^* V(\mathbf{b}, \nu_n)^* \\ &= \left(\boldsymbol{\mathbf{\mathcal{O}}}^\mathrm{fil}(\mathbf{b})\boldsymbol{\mathbf{C}}_{\boldsymbol{\mathbf{b}}}\boldsymbol{\mathbf{\mathcal{O}}}^\mathrm{fil}(\mathbf{b})^\dagger\right)_{ij} \, , \end{align}\] where we assume that the data covariance \(\boldsymbol{\mathbf{C}}_{\boldsymbol{\mathbf{b}}}\) is noise dominant, \[\left(\boldsymbol{\mathbf{C}}_{\boldsymbol{\mathbf{b}}}\right)_{ij}\equiv \frac{\sigma^2_\mathrm{rms}(\nu_i)}{n_{\boldsymbol{\mathbf{b}}}} \delta_{ij}\, .\] As our foreground-filtering method does not correlate between different baselines, the corresponding frequency-frequency covariance for a given line of sight \(\boldsymbol{\mathbf{\hat{\mathbf{s}}}}_n\) is \[\label{eq:img95covariance} \begin{align} &\left(\boldsymbol{\mathbf{C}}^\mathrm{img}_{\boldsymbol{\mathbf{\hat{\mathbf{s}}}}_n}\right)_{ij} \equiv \langle \hat{m}_n(\nu_i) \hat{m}_n(\nu_j)^*\rangle \\ =& \sum_k \left(\frac{n_{\boldsymbol{\mathbf{b}}_k}}{n_\mathrm{bl}}\right)^2 \exp\left(\frac{i2\pi(\nu_i-\nu_j)}{c}\boldsymbol{\mathbf{b}}_k\cdot\boldsymbol{\mathbf{\hat{\mathbf{s}}}}_n\right)\left(\boldsymbol{\mathbf{C}}^\mathrm{filt}_{\boldsymbol{\mathbf{b}}}\right)_{ij}\, . \end{align}\tag{16}\] The image-space frequency-frequency covariance is simply a weighted sum of the visibility-space covariance across all baselines.

The top panel of 2 shows the frequency-frequency correlations, \(\boldsymbol{\mathbf{C}}^\mathrm{img}_{\boldsymbol{\mathbf{\hat{\mathbf{s}}}}_n}/\sigma^2\), across the entire frequency band where we use to filter the foreground. Although ionized bubbles are fairly localized in frequency space, we still choose a wide frequency band because this minimizes signal attenuation during foreground filtering [51], [79]. For longer baselines, this introduces long-range correlations as more line-of-sight modes are filtered (e.g., the dashed green line in the bottom panel of 2). However, since the image-space correlation is a linear combination of visibility-space correlations across all baselines weighted by the redundancy of each baseline group, it is dominant by the behavior of shorter baselines. The bottom panel of 2 shows that there is in fact no significant correlation in the image space beyond its immediate neighboring frequency channels.

3 Signal Modeling↩︎

To accurately forecast and interpret the cross-correlation signal, we use the radiation-magneto-hydrodynamic simulations thesan to derive a signal template which takes into account the selection effects of LAEs. In the following, we give a brief overview of the thesan simulation in Sec. 3.1. Observed properties of LAEs are discussed in Sec. 3.2. Combining these, we present the modeling template for the cross-correlation signal in Sec. 3.3.

3.1 Simulations↩︎

Thesan [52], [53], [55] is a suite of large (95.5 comoving Mpc per side) radiation-magneto-hydrodynamic cosmological simulations run down to \(z = 5.5\), which model reionization by self-consistently combining on-the-fly radiative transfer and realistic galaxy formation modeling from IllustrisTNG [80]–[86]. Here, we adopt the fiducial simulation thesan-1, which resolves dark matter to \(3.1 \times 10^6\, \text{M}_\odot\) and baryonic matter to \(5.8 \times 10^5\, \text{M}_\odot\). Atomic cooling halos are therefore marginally resolved down to masses of \(M_{\rm halo} \gtrsim 10^8\, h^{-1}\, \text{M}_\odot\). Thesan uses the efficient quasi-Lagrangian code arepo-rt [87], [88], an extension of the moving mesh code arepo [89], [90], with additional physics required to self-consistently model reionization. It solves the fluid dynamics equations on an adaptive unstructured Voronoi mesh produced by approximately following the flow of the gas. Gravity calculations utilize a hybrid Tree-PM approach, which splits the force into short- and long-range contributions [91]. The radiation transport equations are solved using a moment-based approach assuming the M1 closure relation [92], [93], with the spectrum discretized in three frequency ranges to accurately capture non-equilibrium photoionization and photoheating from stellar and AGN sources for primordial gas. A reduced speed of light approximation is used with an effective value of \(0.2\,c\), and a birth cloud escape fraction of \(0.37\) is employed to match constraints for the global reionization history. Data products from the thesan simulations are publicly available online for community use [94].

To obtain the properties of the 21-cm field from thesan, we model its brightness temperature via [21] \[\label{eq:21cm95brightness} \begin{align} \delta T_b \approx& 27 \,{\rm mK} (1+\delta_b)x_{\rm HI}\left(1-\frac{T_{\rm CMB}(\nu)}{T_{\rm spin}}\right)\left(\frac{\Omega_b h^2}{0.023}\right) \\ &\times\sqrt{\left(\frac{1+z}{10}\right)\left(\frac{0.15}{\Omega_m h^2}\right)}\left(\frac{H(z)/(1+z)}{\mathrm{d} v_\|/\mathrm{d} r_\|}\right)\,, \end{align}\tag{17}\] where \(\delta_b\) is the baryon overdensity field, \(x_{\rm HI}\) is the fraction of hydrogen that is neutral, \(T_{\rm CMB}(\nu)\) is the CMB temperature at frequency \(\nu\), \(\mathrm{d} v_\|/\mathrm{d} r_\|\) is the gradient of the proper velocity along the line-of-sight direction, and \(T_{\rm spin}\) is defined as the ratio of the occupancy of the spin-1 and spin-0 ground states of the neutral hydrogen: \[\begin{align} \frac{n_1}{n_0} = 3\exp\left(-T_* / T_{\rm spin}\right)\text{ with T_*=0.0681 K}\,. \end{align}\] For the redshift range of interest for this work, it is safe to assume that \(T_\mathrm{spin} \gg T_\mathrm{CMB}\) and ignore the term \((1-T_\mathrm{CMB}/T_\mathrm{spin})\). The remaining quantities in Eq. 17 are obtained using gas properties sampled on a \(512^3\) regular Cartesian grid [94]³. These quantities are binned in redshift space to take into account redshift-space distortions. We note that thesan was run under the cosmological parameters from [95]. In particular, we have \(h = 0.6774\), \(\Omega_m = 0.3089\), \(\Omega_\Lambda = 0.6911\), \(\Omega_b = 0.0486\), \(\sigma_8 = 0.8159\), and \(n_s =0.9667\). All cosmological calculations in this work assume the same cosmology.

3.2 LAE Properties↩︎

Because the neutral IGM is optically thick at the Ly\(\alpha\) wavelength, it is crucial to take into account the correlation between observed LAEs and their surrounding IGM. Here, we describe how we obtain the observed Ly\(\alpha\) properties through different lines of sight for each LAE. These observed properties are empirically calibrated to ensure that the resulting Ly\(\alpha\) luminosity functions match the observations. Details of this process will be described in the upcoming work [96].

First, intrinsic properties of the Ly\(\alpha\) emission from galaxies are calculated directly from thesan. The intrinsic Ly\(\alpha\) luminosity \(L_{\alpha,\rm int}\) incorporates contributions from recombinations, collisional excitations, and unresolved HII regions [55]. The frequency-dependent Ly\(\alpha\) transmission as the photons pass through the local IGM is also accurately captured through an effective absorption treatment with continuous Doppler shifting, i.e. \(\mathcal{T}_{\rm IGM}(\nu) = \exp[-\tau(\nu)]\), extracting sightlines with the colt code [97]–[99]. To account for unresolved galaxy-scale phenomena including dust, outflows, and other effects from the interstellar and circumgalactic medium, an idealized model for a Ly\(\alpha\) point source surrounded by an expanding or contracting gas cloud is applied and calibrated to observational constraints on the Ly\(\alpha\) luminosity functions at \(z=5.7\) and \(z=6.6\) [30]–[32]. The resulting fit is used to calculate the escape fraction \(f_{\rm esc}\) and calibrate the spectral profile for each galaxy. Together with the sightline-dependent IGM transmission \(\mathcal{T}_{\rm IGM}\), these quantities are combined to derive observed Ly\(\alpha\) luminosities \(L_{\alpha,\rm obs}\) and equivalent widths (EW) for each sightline from a galaxy. The observed luminosity is calculated as \(L_{\alpha,\rm obs} = f_{\rm esc} \times \mathcal{T}_{\rm IGM} \times L_{\alpha,\rm int}\) and the equivalent width is calculated as \({\rm EW} = L_{\alpha,\rm obs} / L_{\lambda,\mathrm{cont}}\), where \(L_{\lambda,\mathrm{cont}}\) is the specific luminosity of the stellar continuum surrounding the Ly\(\alpha\) emission line.

3.3 Signal Template↩︎

To derive a signal template for the stacked 21-cm spectrum around LAEs, we utilize 768 lines of sight (in Healpix directions) from each galaxy in thesan. Based on the observed LAE properties along each line of sight, we can implement the selection of any galaxy survey and derive the signal template by stacking the 21-cm spectra only when we can observe an LAE.

As an example, we consider we have a sample of LAEs at \(z_{\mathrm{LAE}}\sim7\). The leftmost panel of 3 shows the 21-cm brightness temperature in thesan at this redshift with LAEs marked in white stars. What does the stacked 21-cm spectrum around these LAEs look like? The right two panels of 3 show how the result differs if we select LAEs based on their intrinsic versus observed properties. Here, a positive value in the \(x\)-axis indicates the direction toward the observer.

We note that the prescriptions in thesan give rise to a particular model of reionization history. To generalize our signal template to account for different reionization scenarios, we calculate the stacked 21-cm spectra with LAEs at various snapshots with different \(\bar{x}_\mathrm{HI}\). Following Eq. 17 , we scale the resulting 21-cm spectra by \(\sqrt{(1+z_{\mathrm{LAE}})/(1+z_{\mathrm{snap}})}\) where \(z_{\mathrm{snap}}\) is the redshift of each snapshot. Hence, each curve in the right two panels of 3 corresponds to a stacked 21-cm spectrum at \(z\sim7\) assuming a different global neutral fraction. Here, we also show the difference between stacking 21-cm spectra around intrinsically bright LAEs (left) versus observed LAEs (right). In both cases, the brightness temperature dips around the center as the IGM are mostly ionized there, except for a small emission peak from the neutral hydrogen within the galaxies. However, if we stack around intrinsically bright LAEs, the absorption troughs do not go all the way to zero, especially when the IGM is more opaque (higher \(\bar{x}_\mathrm{HI}\)). This is because not every LAE resides in an ionized bubble, whereas the observed LAEs are guaranteed to be surrounded by a more transparent IGM. A major feature of this result is that the amplitude of the stacked 21-cm spectrum becomes a direct tracer of the global neutral fraction \(\bar{x}_\mathrm{HI}\). This coincides with the finding in [49] as the amplitude of the stacked spectrum is approximately equivalent to the two-point correlation function between 21-cm and galaxies at very small scales. Moreover, we see that the absorption profiles are largely symmetric in the left panel. The asymmetry in the rightmost panel arises because observed LAEs preferentially reside in the back side of ionized regions. Observationally speaking, while stacking around a sample of LAEs yields spectra like those on the right-hand side, we can obtain signal that bear more resemblances to the template on the left-hand side if we stack around galaxies detected through other emission lines such as [OIII].

Here, we choose a Ly\(\alpha\) luminosity threshold of \(10^{42}\,\mathrm{erg}/\mathrm{s}\) as an example. At \(z\sim7\), this roughly corresponds to a survey with a flux limit of \(10^{-18}\,\mathrm{erg}/\mathrm{s}/\mathrm{cm}^2\). This can be achieved by large ground-based spectroscopy [100]–[102] and is approximately an order of magnitude deeper than what the Roman grism survey will achieve [103]. Changing the selection criteria in either Ly\(\alpha\) luminosity or the equivalent width modifies the morphology of our signal, but the variation is not significant given the sensitivity of current 21-cm experiments.

4 Forecast↩︎

4.1 Setup↩︎

In this work, we focus on forecasting the detectability of a cross-correlation signal with HERA. HERA is a 350-element radio interferometer located in the Karoo desert in South Africa. In particular, we consider only the 320 elements that form a compact core array with dishes that almost touch each other. The maximum inter-antenna distance is 292 meter and the shortest baseline is 14.6 meter. The antenna configuration can be seen in 4.

Currently, science data are being taken by a subset of commissioned antennas, while new antennas are continuously being added. In particular, 172 antennas marked in green in 4 have been taking data since 2022. We use this subset of antennas to form a conservative forecast to investigate whether a cross-correlation detection is possible with data that have already been taken by HERA to date.

To simulate the signal observed in real HERA images, we take the stacked signal template derived in Sec. 3.3 and generate simulated visibility following Eq. 1 . In this work, we assume that each antenna has an Airy beam profile, \[B\left(\hat{\mathbf{s}}(\theta, \phi); \nu\right) = \left[\frac{2J_1(2\pi\nu a \sin\theta/c)}{2\pi\nu a \sin\theta/c}\right]^2\, ,\] where \(J_1\) is the Bessel function of the first kind, \(\theta\) is the zenith angle, and \(a\) is the aperture radius which we set to be six meters to mimic an underillumination HERA dish [23], [104], [105]. The visibility from each baseline is then filtered according to the procedure outlined in Sec. 2.2 and 6. The foreground-filtered visibilities are then combined and map to the image space following Eq. 3 in which the normalization is chosen to be as Eq. 10 . One important feature is that this entire procedure—from signal to mock image—is linear. Hence, mapping a stacked signal is equivalent to mapping individual galaxies and stacking them afterward. We denote the observed stacked signal as \(s_\mathrm{obs}(\nu)\).

Another source of uncertainty in the signal comes from the redshift of the LAEs. In order to perfectly align and stack the 21-cm spectra, we need to know precisely the redshift of these LAEs, at least to the precision that matches the frequency resolution of the radio instrument. At redshift 7, the frequency resolution of HERA corresponds to a redshift uncertainty \(\sigma_z\approx0.005\). This matches well with the uncertainties provided by a space-based grism or a large ground-based spectroscopy. Here, we consider three different redshift uncertainties, \(\sigma_z = 0.001, 0.01, 0.1\). These redshift uncertainties are incorporated by perturbing the location of the galaxies when deriving the stacked signal template in Sec. 3.3.

5 shows the resulting stacked 21-cm signal as observed by HERA. These signals have been foreground-filtered; hence, the smooth component of the signal is removed. The three panels indicate the different signal strength if the global neutral fraction \(\bar{x}_\mathrm{HI}\) is \(0.86\), \(0.79\), or \(0.44\) at redshift \(\sim7\). The default reionization history in thesan predicts \(\bar{x}_\mathrm{HI} = 0.44\) at \(z\sim7\). A higher neutral fraction means bigger contrast between the average brightness temperature of the IGM and the ionized bubble (which has a brightness temperature around \(0\)). In each panel, the three different curves show the effect of redshift uncertainties on the stacked signal. The cross-correlation signal is maximized with minimal redshift uncertainties. We see that the signal almost vanishes for \(\sigma_z = 0.1\) (dotted blue), which is typical for a photometric redshift estimate. Spectroscopy confirmation of these LAEs is therefore necessary for a successful cross-correlation detection through stacking along the line-of-sight direction.

Throughout this work, we focus on forecasting the cross-correlation signal around redshift \(7\). This is the redshift where most LAEs are currently being identified on the ground due to sky lines. Searching for such a signal at higher redshift might be easier due to the increased signal strength. At the same time, a more opaque IGM makes it harder to identify a large sample of LAEs. With the launch of the James Web Space Telescope, we are starting to more evenly sample LAEs at even higher redshifts [35], [106], [107]. With upcoming wide-area grism surveys on Euclid [108] and Roman [109], it will soon be possible to probe cross-correlation signals across the entire reionization history.

We note that while this forecast focuses on HERA, the result presented here can be generalized to other experiments. As we do not attempt to map structures in the spatial direction, the exact layout of antennas is less relevant to our sensitivity forecast. This is seen in Eq. 13 as the noise level in a given pixel in the 21-cm image depends only on the total number of baselines in the experiment. For experiments with steerable antennas, the dilution of signal from antenna’s primary beam in Eq. 11 can be avoided by tracking a source. The forecast for a transit telescope therefore forms a conservative lower bound for an equivalent tracking experiment. For context, the full HERA core array has 51360 baselines, while the 172 antennas that are currently taking data form 14706 baselines.

4.2 Results↩︎

The main question we focus on in this work is: what are the observing resources required to make a significant cross-correlation detection through stacking 21-cm spectra around LAEs? Using the observed signal template \(s_\mathrm{obs}(\nu)\) derived in Sec. 4.1, we calculate the signal-to-noise ratio for the cross-correlation to be \[\mathrm{SNR}\mathrel{\vcenter{:}}=\sqrt{\sum_{ij} s_\mathrm{obs}(\nu_i) C^{-1}_{ij} s_\mathrm{obs}(\nu_j)}\, ,\] where the covariance matrix \(\boldsymbol{\mathbf{C}}\) is given in Eq. 16 . The covariance matrix takes into account the correlation between different frequency channels introduced by foreground filtering, and its variance is determined by the amount of observing resources used to reduce the noise level.

The noise level in the stacked 21-cm images is determined by three factors: the integration time of each nightly observation, the number of nights we can observe each object, and the number of objects (LAEs) available for stacking. Based on the results described in Sec. 2.3, the maximum sensitivity we can achieve around each galaxy occurs when it is observed for roughly 45 minutes as it transits the zenith. This gives an effective coherent averaging time of 25 minutes. Ideally, each source can be observed repeatedly every night for half the year. We present our forecast in terms of the number of nights and LAEs required to obtain a significant cross-correlation signal.

The top panel of 6 shows the minimum amount of resources required for a \(3\sigma\) cross-correlation detection with HERA. This assumes the signal to be the strongest as seen in 5. The green lines are for the current HERA layout, while the black lines are for the full HERA-320 (see 4). Lower redshift uncertainties also make the signal stronger and easier to be detected. Since 2022, three seasons of HERA Phase II data have been taken. With roughly \(540\) nights of available data, a cross-correlation detection starts to become possible with around \(43\) (\(90\)) LAEs assuming \(\sigma_z = 0.001~(0.01)\).

The fiducial reionization model adopted in thesan predicts \(\bar{x}_\mathrm{HI}\sim0.44\) at \(z\sim7\). In this case, around \(1000\) LAEs are required to detect the cross-correlation signal with the current HERA dataset. However, we note that the neutral fraction changes rapidly as we move to higher redshift while the properties of the 21-cm maps remain roughly the same. Moving to \(z\sim8\), we need about \(300\) LAEs to make a detection with the same dataset. In reality, as we have access to 21-cm information across a wide range of redshift, we can make inference at any redshift bin where we have a significant sample of LAEs.

Even in the event where there are not enough LAEs to make a detection, we can turn a non-detection with a given number of LAEs into an upper limit on the global neutral fraction. This is because the signal strength is proportional to the global neutral fraction \(\bar{x}_\mathrm{HI}\). A non-detection at a given sensitivity suggests that the neutral fraction must be lower than a certain level. The bottom panel of 6 shows an example of upper limits that can be derived from existing HERA data, if we do not measure a cross-correlation signal after stacking a given amount of spectroscopically confirmed LAEs. Such a result will be a unique constraint on the reionization history measured directly from the IGM.

5 Conclusion↩︎

We have presented a detailed analysis of a direct way to detect a cross-correlation signal between galaxies and the 21-cm field—through stacking 21-cm image cubes around Lyman-alpha emitters (LAEs). Detections of the cross-correlation signal are of increasing importance as 21-cm experiments approach the sensitivity to detect an auto power spectrum. Cross-correlations can be crucial to verify any 21-cm auto power spectra detection is cosmological in nature. Moreover, we have shown that even an upper limit on the cross-correlation signal provides additional cosmological information. Current limits on the 21-cm auto power spectra are only sensitive to the heating of the intergalactic medium (IGM). By combining with information provided by high-redshift galaxies, a simple cross-correlation through stacking can already provide additional constraints on the reionization history.

One of the foremost requirements to make any inference in 21-cm cosmology is to understand the statistical properties of the data. In the case of stacking, the 21-cm image cubes. In this work, we have derived the statistical properties of foreground-filtered 21-cm maps generated with a series of linear operations. We choose to work with a linear map-making and foreground-filtering algorithm to ensure that we can correctly estimate and propagate the statistical properties of our maps.

On the theory side, we have derived a signal template using state-of-the-art radiation-magneto-hydrodynamic cosmological simulations thesan. This provides a physically driven model that connects galaxies to their surrounding ionized bubbles. Our results also take into account the selection effect of LAEs. The correlation between the observed Lyman-alpha properties of a galaxy and its surrounding IGM is accounted for with radiative transfer modeling. Unresolved galaxy-scale phenomena are further calibrated to observational constraints on the Lyman-alpha luminosity functions at high redshift. An important feature we have confirmed is that the stacking signal is proportional to the averaged neutral fraction of the universe.

In conclusion, our forecast suggests that we are in a position to place significant constraints on the reionization history with existing HERA data through stacking. A sample of around 50 (100) LAEs with redshift uncertainties of 0.001 (0.01) is sufficient to begin with. Around \(300\) (\(1000\)) LAEs are required to make a detection with existing HERA dataset if the neutral fraction is around \(0.6\) (\(0.4\)) at \(z\sim7\). Such a sample of LAEs could soon be available with upcoming space-based grism surveys. The prospect of a detection will improve in the meantime with more commissioned antennas and continuous observations. Significant improvements in detectability can also be achieved with advanced analysis techniques in inferring cosmological modes within the foreground wedge [110], [111], or through designing experiments that reduce the foreground wedge [112]. These methods can alleviate the amount of signal loss during foreground mitigation.

Lastly, we discuss some of the limitations of this work. On the theory side, our signal template is slightly limited by the box size of the simulation. At higher neutral fractions, the number of LAEs that pass the selection criteria is low, leading to a nosier signal template, as can be seen in 3. With larger boxes of radiation-hydrodynamic simulations soon to be available, we expect a significant improvement in comic variance uncertainties in our signal template. On the observation side, an important assumption we have made is that the 21-cm image cubes are free of systematic effects after foreground filtering. This might not be true in the real world due to imperfect or insufficient characterizations of the instrument, especially after significantly averaging down the data. Another source of systematic uncertainty that we omit is the contamination from surrounding emissions. The brightness temperature at a given pixel in our map is a convolution of all pixels on the sky with the point spread function (synthesized beam) of the array. If the point spread function varies smoothly as a function of frequency, this should only add smooth contamination to our signal and would be removed by foreground filtering. However, careful consideration of these systematic effects is necessary when handling real data.

6 Foreground Filtering with Discrete Prolate Spheroidal Sequence↩︎

In this appendix, we describe the detailed procedure for generating the linear operator \(\boldsymbol{\mathbf{\mathcal{O}^\mathrm{fil}}}\) that performs foreground filtering using the discrete prolate spheroidal sequence (DPSS). DPSS is the set of eigenvectors to the prolate matrix \(\boldsymbol{\mathbf{B}}\) where \[\label{eq:prolate95matrix} B_{ij} = \frac{\sin{2\pi T(\nu_i - \nu_j)}}{\pi (\nu_i - \nu_j)}\, .\tag{18}\] Here, \(T\) denotes the baseline-dependent delay range below which we want to filter out the smooth foreground component. In this work, for a baseline vector \(\boldsymbol{\mathbf{b}}\), we choose \(T=|\boldsymbol{\mathbf{b}}|/c\) to remove all modes within the foreground wedge. We note that in practice, one might choose a slightly larger \(T\) to also filter out foreground contamination into the EoR window due to various systematic effects.

To fit and filter the smooth foreground, we use all eigenvectors \(f_i(\nu)\) of Eq. 18 with eigenvalues \(\lambda_i \geq 10^{-12}\). The eigenvalues of the prolate matrix are always between \(0\) and \(1\) and denote how localize each eigenvector is within the given delay range \([-T, +T]\) (\(1\) being completely localized). A lower eigenvalue cut here ensures that we have a more complete basis, which reduces the amount of residual foreground. Although this also gives us some eigenvectors that could filter out signal in the high delay EoR window, the number of these eigenvectors is fairly limited [113] and the signal loss is carefully quantified in our work.

Once a set of basis \(\{f_i\}_{i=1}^N\) is chosen, we fit the smooth foreground by solving the linear system \[\boldsymbol{\mathbf{v}}_\mathrm{obs} = \boldsymbol{\mathbf{A}}\boldsymbol{\mathbf{\alpha}} + \boldsymbol{\mathbf{n}}\,,\] where \(\mathbf{v}_\mathrm{obs}\) is the observed visibility of a given baseline, \(A_{ij} = f_j(\nu_i)\) is the design matrix, \(\boldsymbol{\mathbf{\alpha}}\) is the DPSS coefficient we wish to solve, and \(\mathbf{n}\) is the instrumental thermal noise that can be modeled with Eq. 4 . The maximum likelihood estimator of \(\boldsymbol{\mathbf{\alpha}}\) is then \[\label{eq:b95ML} {\boldsymbol{\mathbf{\hat{\alpha}}}} = (\mathbf{A}^\dagger\mathbf{N}^{-1}\mathbf{A})^{+} \mathbf{A}^\dagger\mathbf{N}^{-1}\mathbf{v}_\mathrm{obs}\,,\tag{19}\] where \(\boldsymbol{\mathbf{M}}^+\) denotes the Moore–Penrose pseudo-inverse of a matrix \(\boldsymbol{\mathbf{M}}\). By subtracting the best-fit smooth foreground, we obtain the filtered visibility \(\boldsymbol{\mathbf{v}}^\mathrm{fil}\) \[\boldsymbol{\mathbf{v}}^\mathrm{fil} \mathrel{\vcenter{:}}= \boldsymbol{\mathbf{v}}_\mathrm{obs} - \boldsymbol{\mathbf{A}}\boldsymbol{\mathbf{\hat{\alpha}}} = \left[\boldsymbol{\mathbf{I}} - \boldsymbol{\mathbf{A}}(\mathbf{A}^\dagger\mathbf{N}^{-1}\mathbf{A})^{+} \mathbf{A}^\dagger\mathbf{N}^{-1}\right]\boldsymbol{\mathbf{v}}_\mathrm{obs}\,.\] Hence, \(\left[\boldsymbol{\mathbf{I}} - \boldsymbol{\mathbf{A}}(\mathbf{A}^\dagger\mathbf{N}^{-1}\mathbf{A})^{+} \mathbf{A}^\dagger\mathbf{N}^{-1}\right]\) is the linear foreground filter operator in Eq. 6 . In this work, we filter the visibility across a wide frequency range to reduce signal attenuation [51], [79]. For the results presented in this work, foreground is fit with data between \(108.08\) to \(234.30\) MHz which correspond to the full range of HERA data above the FM radio band.

References↩︎