The feasibility of weak lensing and 21cm intensity mapping cross-correlation measurements
Abstract
One of the most promising probes to complement current standard cosmological surveys is the HI intensity map, i.e. the distribution of temperature fluctuations in neutral hydrogen. In this paper we present calculations of the 2-point function between HI (at redshift ) and lensing convergence (). We also construct HI intensity maps from N-body simulations, and measure 2-point functions between HI and lensing convergence. HI intensity mapping requires stringent removal of bright foregrounds, including emission from our galaxy. The removal of large-scale radial modes during this HI foreground removal will reduce the HI-lensing cross-power spectrum signal, as radial modes are integrated to find the convergence; here we wish to characterise this reduction in signal. We find that after a simple model of foreground removal, the cross-correlation signal is reduced by 50-70%; we present the angular and redshift dependence of the effect, which is a weak function of these variables. We then calculate S/N of HI detection, including cases with cut sky observations, and noise from radio and lensing measurements. We present Fisher forecasts based on the resulting 2-point functions; these forecasts show that by measuring correlation functions in a sufficient number of redshift bins, constraints on cosmology and HI bias will be possible.
keywords:
Radio lines: General, Gravitational Lensing: Weak, Large Scale Structure of Universe1 Introduction
The clustering of matter in the Universe provides an important insight into the origins and evolution of the cosmic structure. Inflation predicts that early structure formation generates a near-Gaussian random field in overdensity; evolution due to gravity causes late-time large-scale structures to exhibit non-Gaussian features. Two point statistics of the density field at different redshifts capture information about the evolution of structures, and correlation functions between different pairs of cosmological probes can precisely constrain cosmological parameters (Abbott et al., 2018; Tröster et al., 2022; Pandey et al., 2021; Fang et al., 2022; Upham et al., 2019). Two dimensional surveys of the cosmic microwave background (CMB) have been effectively carried out through the last few decades (Planck Collaboration et al., 2018; Hinshaw et al., 2013). The complement to this is deep sky observations of the 3-dimensional galaxy and dark matter fields. While conventional optical and infrared surveys have high angular resolution, long integration times are needed for these to obtain precise redshifts via spectroscopy. In contrast, photometric surveys provide faster redshift capture but less radial resolution (Fernandez-Soto et al., 2001).
To complement the low radial resolution of optical photometric surveys, alternative techniques with higher radial resolution are desirable; radio 21cm intensity mapping is a rapidly developing candidate for this purpose. Unlike most optical surveys, this technique does not measure the brightness of individual objects, but focuses on the larger-scale fluctuations in intensity of the 21cm radio signal from neutral hydrogen (HI). The temperature fluctuations can be used as a tracer for the underlying cosmic density field. This intensity mapping is a complementary technique to a photometric survey, with excellent redshift resolution but lower angular resolution (Bull et al., 2015). Hence, combining HI and optical surveys is potentially valuable, as the two techniques compensate for each other’s limitations (Cunnington et al., 2019b; Square Kilometre Array Cosmology Science Working Group et al., 2020).
Recently, HI intensity mapping techniques have been actively developed (Santos et al., 2010; Harker et al., 2010; Mao et al., 2008; The CHIME Collaboration et al., 2022; Wolz et al., 2016; Cunnington et al., 2022). The Canadian Hydrogen Intensity Mapping Experiment (CHIME) (CHIME Collaboration et al., 2022) has provided a detection of HI via cross-correlations with three probes of Large-Scale Structure (LSS), namely luminous red galaxies (LRG), emission line galaxies (ELG), and quasars (QSO) from the eBOSS clustering catalogs at high significant levels, (LRG), (ELG), and (QSO). Cunnington et al. (2022) have detected the correlated clustering between MeerKAT measurements of HI and galaxies from the WiggleZ Dark Energy Survey at significance. Intensity mapping is therefore on its way to becoming an independent observational probe, providing useful information from low to high redshifts, via future surveys with radio telescopes such as MeerKAT (Pourtsidou, 2017; Pourtsidou et al., 2017; Spinelli et al., 2022) and the Square Kilometre Array, SKA (Santos et al., 2015).
The major challenge for the intensity mapping technique is that the foreground signals are much stronger than the cosmic HI brightness temperature, especially due to the galactic plane synchrotron radiation (Spinelli et al., 2018; Su et al., 2018; Switzer et al., 2013). Hence several studies of 2-point functions between HI and optical () have focused on the impact of foreground removal (Chapman et al., 2012; Cunnington et al., 2019b; Padmanabhan et al., 2020; Cunnington et al., 2020; Spinelli et al., 2022). The study by Cunnington et al. (2019b) shows that the foreground removal affects 2-point function characteristics, especially when the redshift resolution is broad, as is the case in optical photometric surveys.
There are also numerous optical surveys measuring gravitational lensing shear () which distorts the shape of galaxy images; this is sensitive to density fluctuations of all the matter present along a line of sight, whether baryonic or dark matter. It is therefore of interest to consider the viability of the cross-correlation , which will be able to be studied using a combination of lensing and IM surveys (Abbott et al., 2018; Baxter et al., 2019; Hu & Jain, 2004; The CHIME Collaboration et al., 2022; CHIME Collaboration et al., 2022; Cunnington et al., 2022). The density projection along the unperturbed light ray trajectory, also known as ’lensing convergence’ can be considered instead of as both share the same statistical properties. The 2-point functions between the pairs of and HI could improve cosmological constraints and break degeneracies such as that between HI bias () and clustering amplitude.
However, removing the HI foreground potentially affects these 2-point statistics, as the foreground removal effectively subtracts large-scale radial modes to which lensing is sensitive. In this paper we will calculate the cross-correlation function between convergence and 21cm intensity mapping, and will explore whether the foreground subtraction significantly hampers the cross-correlation measurement. We also explore whether the foreground removal impacts the viability of cosmological constraints from HI-HI and HI correlations.
To achieve this, we will present theoretical and simulation approaches for calculating the -HI signal. We will then consider the effect of foreground removal on the signal, showing that the impact is significant (approximately a factor of 2 in signal reduction) but not lethal. We will then use the Fisher information matrix to make cosmological parameter forecasts for ideal and realistic surveys (including cut sky and the inclusion of telescope-specific noise), deploying the cross-correlation between convergence and intensity mapping, always including the effect of foreground removal. We discuss lensing convergence and HI simulation catalogues in Section 2, including modeling of the 2-point functions. We describe the HI foreground removal and its effect on -HI 2-point functions in Section 3. We present our Fisher forecasts for surveys in Section 4, effects of instrumental noise in Section 5 and present our conclusions in Section 6.
2 -HI 2-point statistics: theory and simulations
In this section we discuss the relevant 2-point statistics. We shall start with theoretical calculations of 2-point functions of lensing convergence () and neutral hydrogen intensity maps (HI) in subsection 2.1. We will then discuss the generation of catalogues and HI modelling from simulations of the matter overdensity . The comparison between theoretical calculations and simulations is shown in subsection 2.2. The simulated HI maps will be used in the next section 3 where the foreground removal will be discussed.
2.1 Modeling the 2-point functions
In this subsection, we describe the modeling of the 2-point functions. We begin by considering how to calculate the observable quantities, namely weak lensing convergence and HI temperature fluctuations . We will then turn to the angular cross-power spectra. We denote as the power spectra between fields, HIiHIj as the cross-power spectra between HI fields, and HI as the cross-power between and HI. The dummy indices and refer to the th and th redshift bins. We will calculate the lensing convergence in an arbitrary direction on the sky using the Born approximation, projecting the matter overdensity along an unperturbed ray direction. This can be computed by (Bartelmann & Schneider, 2001)
(1) |
where is comoving distance, is the matter density parameter at the present epoch, is the Hubble parameter today and the subscript refers to the source plane. For lensing of distributed sources in redshift bins , the integrand is modified by including a source distribution, so that the integration now becomes
(2) |
where the lensing weight is given by
(3) |
where is the lensing source number density, and is its average in the th redshift bin.
HI will be a biased tracer of matter overdensity, so we write = , where is the HI bias at a given redshift and is the average temperature. The projected temperature fluctuation at the th redshift bin is then
(4) |
where
(5) |
where is the HI source number density, and is its average in the th redshift bin.
Battye et al. (2013) show that for a given redshift , can be estimated by
(6) |
where is the dimensionless Hubble function at redshift . The HI density parameter could be approximated to be (Battye et al., 2013). However, throughout this research we shall follow the fitting formula for the SKA-MID I by Square Kilometre Array Cosmology Science Working Group et al. (2020)
(7) |
Constraining the HI bias will be discussed later in section 4.
Using the Limber approximation, the angular power spectra are given by
(8) |
where is the matter power spectrum (LoVerde & Afshordi, 2008). We compute the nonlinear power spectrum using the Boltzmann code CAMB (Lewis & Bridle, 2002) with the Halofit extension to nonlinear scales (Takahashi et al., 2012).
Since the ray tracing simulations by Takahashi et al. (2017) which we use below adopt a comoving bin size = 150 Mpc (see section 2.2), we choose a radial selection function for as a normal distribution around a central comoving position with = 150 Mpc. This corresponds to the frequency bandwidth selected. In practice, the frequency range and bandwidth will depend on the particular radio telescope being used; for example, BINGO (Baryon Acoustic Oscillations in Neutral Gas Observations) has operational frequency from 960 MHz to 1260 MHz (Battye et al., 2013; Wuensche & the BINGO Collaboration, 2019) and MeerKAT’s frequency bandwidth ranges from 900-1185 MHz and 580-1000 MHz for L-band and UHF-band, respectively (Wang et al., 2021; Cunnington et al., 2022) .
2.2 Lensing Convergence and HI Intensity Maps
The full-sky gravitational lensing mock catalogues by Takahashi et al. (2017) have been used throughout this work. They are based on a multiple-lens ray-tracing approach through N-body cosmological simulations. The datasets include weak lensing maps (convergence, shear and rotation data) up to redshift 5.3, and halo catalogues. The catalogues provide 108 realisations of N-body simulations, 35 of which are used in this research (due to storage limitations). The N-body simulations were produced with periodic boundary conditions following dark matter gravitational evolution without baryonic processes. 14 simulation boxes of side length = 450, 900, 1350, …, 6300 Mpc are nested to represent a region of the Universe in which lensing occurs; each box contains particles. The fields are obtained by tracing the light ray path through planes with separation 150 Mpc. By calculating the Jacobian matrix along the light path, the lensing convergence , shear lensing and rotation angle can be obtained, via
(9) |
The convergence maps were created in the HEALPIX scheme with NSIDE of 4096 (Górski et al., 2005), which contain 200 megapixels. While this resolution is appropriate to study nonlinear structure and matches forthcoming galaxy surveys such as EUCLID 111https://www.euclid-ec.org/ and DESI 222https://www.desi.lbl.gov/, the cross-correlation between the lensing convergence and the HI intensity map is limited by the lower angular resolution of HI intensity maps expected with real radio telescopes. Therefore the resolution is reduced to NSIDE of 512; this is not only appropriate for our 2-point function measurements but also decreases the storage space requirement and computational time.
We will first consider a convergence map at a specific optical lensing catalogue source redshift, which we choose as 0.78, for which the lensing will significantly occur at the redshift of an intensity map at redshift. This particular choice of redshift allows us to compare our results to current and forthcoming optical and radial surveys (Baxter et al., 2019; Square Kilometre Array Cosmology Science Working Group et al., 2020; Euclid Collaboration, 2019; Pourtsidou et al., 2017; Santos et al., 2015). We will then extend to multiple lensing planes (see Table 1).
We turn now to generating our IM maps. Crucially, we will emulate removal of the IM foreground by removing the radial temperature fluctuations on large scales. The foreground removal will be discussed in detail in Section 3.
First we need to make the pre-foreground-removal IM maps. Instead of calculating the individual HI masses from halo catalogues, we assume that HI is a biased tracer of the total matter overdensity field (see Eq. 4 and 5),
(10) |
where is a HI bias. For instance the parametric form for adopted by Cunnington et al. (2019b) is
(11) |
Since the neutral hydrogen signal is measured as the surface brightness temperature, we shall refer to the HI intensity map as the temperature fluctuation :
(12) |
We apply this equation to the overdensity map obtained from Takahashi et al. (2017) catalogues to create HI intensity maps. Fig. 2 shows the uncleaned and cleaned intensity maps from one realisation for the zoom-in patch with area 5 5 square degrees.
As we are interested in the 2-dimensional projection of cosmological fields on the sky, together with their power spectra, it is convenient to describe these fields in spherical harmonics:
(13) |
where and are spherical harmonics and their coefficients respectively (Pratten et al., 2016; Castro et al., 2005; Heavens, 2003). represents an arbitrary cosmological field; in this work it can be either lensing convergence or HI temperature fluctuations. The angular power spectrum is then an average of over modes:
(14) |
where and stand for the cosmological fields at given redshifts and , respectively.
The cross power-spectrum for HI and lensing can be easily measured via HEALPIX’s anafast routine especially if the data is for the full sky (however, if the data has missing regions or a cut sky, pseudo- methods are required (Brown et al., 2005; Upham et al., 2019)).
Using this routine, we obtain cross-power measurements for the HI and fields. We measure the cross-power spectra of 35 realisations and evaluate their mean; we show the results in Fig. 3. Here the lensing convergence is measured at the central redshift 0.78 and HI is measured at the central redshift 0.3. Fig. 3 also displays a comparison between theoretical 2-point statistics and the measurements from the mock catalogues. We then measure the covariance matrices () of 2-point statistics from 35 realisations. The error bars are the square root of the diagonal elements of () of the estimators. The correlation matrix for () are shown in Fig. 6.
We see that the measurements from simulations agree very well with our theory curves on this plot, which indicates that our theoretical calculation and selection function successfully match the simulations. Due to the lens shell approximation of the ray tracing code, the measured is slightly affected at very high (see red line on Fig. 3).
3 HI foreground removal and its effect on HI 2-point functions
The HI signal is small compared to its foregrounds such as free-free thermal emission, extragalactic radio sources and Galactic synchrotron. For example, the synchrotron () emission temperature which can be modelled by (Platania et al., 1998; Smoot & Debono, 2017), is approximately 3 to 4 orders of magnitude larger than at low redshift. Thus, 21cm foreground removal is a major challenge for HI cosmology. Several studies suggest that the foreground spectrum appears to be smooth in the radial direction (Cunnington et al., 2019a, b; Shaw et al., 2014). This is equivalent to being present in the long radial wavelengths in Fourier space. We therefore remove such modes in the line of sight background temperature fluctuations .
Since the calculation of lensing involves integration along the light path (Eq. 3), which will have a contribution from long-wavelength radial modes, the HI foreground removal is a concern for the existence of the HI cross-correlation (i.e. we have just removed such modes from the HI signal). In this section we therefore seek to ascertain the degree to which the HI cross-correlation survives foreground removal.
Here we follow the method for foreground removal emulation by Cunnington et al. (2019b). The cleaned intensity map can be approximated as
(15) |
where is the uncleaned signal in direction at redshift . is defined by:
(16) |
so that is the mean surface brightness temperature fluctuation along the entire line of sight. This is an initial very approximate model of Principal Component Analysis (PCA) foreground removal, as most dominant components are included in the line of sight expectation temperature fluctuations . It is worth mentioning that this blind foreground removal technique assumes the smoothness of the foreground. However, this smoothness can be hampered by non-smooth features of the beam, e.g. beamwidth of the radio dish, and some oscillating features in all bands of MeerKAT. A simple dependence of the beam could generate artificial HI signals. This leads to the conclusion in Spinelli et al. (2022) that it is fundamental to develop accurate beam deconvolution algorithms and test data post-processing steps carefully before cleaning. This topic of beam deconvolution is beyond the scope of our research; here we shall assume that the behaviour is sufficiently small. For more sophisticated foreground cleaning methods we encourage the reader to explore e.g. Cunnington et al. (2023).
In this work we adopt the same bias model as Cunnington et al. (2019a):
(17) |
where , , and are set to 1, 0.67, 0.18 and 0.05, respectively. Cunnington et al. (2019a) obtained this parameter set by investigating HI as a biased tracer of the large-scale structure via HI intensity map and optical galaxy number density cross-correlations (see Eq. 39 from Cunnington et al. (2019a)). We use this as a fiducial model since the HI redshift range in our work is similar to Cunnington et al. (2019b, a). Note that in this model, we solely account for the redshift evolution of HI bias and assume any transverse scale dependence of the bias is negligible. Martin et al. (2012) shows that this is a good approximation for scales > 10 Mpc, which are our main interest.
We measure with two choices of maximum redshift, = 1 and 3. = 3 corresponds to futuristic HI-galaxy surveys (Square Kilometre Array Cosmology Science Working Group et al., 2020). On the other hand = 1 is an approximate limit for HI maps with SKA1-MID and MeerKAT (Square Kilometre Array Cosmology Science Working Group et al., 2020; Cunnington et al., 2022).
We use these intensity maps with removed foreground to calculate the auto-power spectra of the intensity map (), and the cross-power spectra between HI and (). We compare the signal of removed and unremoved , resulting in Fig. 4 and Fig. 5. From Fig. 4, we note that foreground removal strongly affects the signal on large scales. However we find that on smaller scales, at , the foreground removal does not erase the power spectrum; the signal is scaled down by a factor () which is close to constant over a range of modes from 10 to 1000. Hence, in section 4, when cosmological constraints following foreground removal are considered, the estimation of cosmological parameters is based on the signal where . We describe the mean signal drop by:
(18) |
From Fig. 4, we see that the higher the maximum redshift of the survey in which we remove the LOS signal, the less is the effect on the cleaned cross-correlation signal, as more radial modes are preserved in the removal process. Fig. 4 and 5 further indicate that the signal in the HI 2-point correlations drops by approximately the same factor across a wide redshift range, if we remove the background noise up to a particular redshift , when cross-correlating to the field at a fixed redshift. Fig. 5 also implies that if the maximum redshifts > .
4 Fisher forecast
In previous sections, we have presented the theoretical 2-point statistics for the HI-lensing cross-correlation, and have examined the impact of HI foreground removal on the cross-power spectrum. The results indicate that foreground removal reduces the 2-point statistics by a modest factor.
Here we begin the exploration of HI correlations as a tool for cosmological constraints. In particular we will make a Fisher information matrix forecast of this correlation in the case of low instrument noise (but including our foreground subtraction model); this will assess the best-case capacity of this probe to constrain cosmology, when one is dominated by large-scale structure fluctuations in the HI and lensing fields. We will then examine more realistic cases with cut sky and the inclusion of instrumental noise.
4.1 The Fisher matrix
The Fisher information matrix is a useful tool to estimate the expected uncertainty in cosmological parameters for forthcoming experiments (Heavens, 2003; Tegmark et al., 1997). Assuming that the model parameters are distributed by a multivariate Gaussian likelihood , the Fisher matrix can be calculated as
(19) |
where . The Fisher matrix can be used to obtain the minimum uncertainty () in parameter estimation due to the Cramér-Rao inequality (Mendez et al., 2014; Kamionkowski et al., 2011),
(20) |
which is equivalent to a 68 confidence level. For a dataset where the uncertainties are Gaussian, the Fisher matrix can be calculated by (Tegmark et al., 1997)
(21) |
where denotes the covariance matrix of the data, , the derivative data matrix + , and is an expectation value of the data vector . The comma symbol means the partial derivative operator with respect to the parameter, . Note that all derivatives are performed at the maximum likelihood point.
As we expect only small changes in the covariance matrix under a modest change in cosmological parameters (see Sec. 2 and 3), the first term on the right hand side in Eq. 21 will be negligible. Then the Fisher matrix can be written
(22) |
We calculate the cross-power spectra by using Eq. 8 with Planck 2018 cosmological parameters (Planck Collaboration et al., 2018). We calculate the covariance matrices of , HI and HIHI from measured cross power-spectra of 35 realisations of the N-body simulation by Takahashi et al. (2017). All the HI temperature fluctuation maps which we use take into account foreground removal. We also calculate the correlation matrices and show these in Fig. 6; these do not indicate significant correlations between bins.
4.2 Cosmological Constraints for Single-Slice Cross-correlations
In this section, the cosmological constraint viability of and HI cross-correlations is explored. We first start with the simplest observational configuration, considering only one redshift slice of HI and . We further assume that the behaves as in Eq. 17. We use the Planck 2018 cosmological parameters as the fiducial cosmology (Planck Collaboration et al., 2018). The fiducial cosmological parameters are = 0.67, = 0.3, = 0.82, = 0, = 0.7, = 0.06, and = 0.96. To make a covariance matrix of cross-power spectra for the Fisher matrix (Eq. 22), we combine modes into 15 bins; each bin contains 101 modes with and averages over 35 realisations. We first consider the 32 functions for a joint analysis of , and , where the numbers in brackets are the central redshifts. We choose these central redshifts as examples of current HI and lensing surveys’ central redshifts. The ‘22’ functions refer to the same combination but exclude the weak lensing-HI cross power spectrum.
Fig. 7 shows the joint likelihood obtained via Fisher matrices (see Eq. 22). We see that single redshift slice correlations of (green) and HI-HI(grey) provide relatively weak constraints, while 2x2pt and particularly 3x2pt are more promising, with few- to ten-per cent constraints available on parameters in this low noise case. The zoom-in version of pt functions is shown on the right hand side of Fig. 7; these likelihood contours, which include the cross-correlation, show a significant improvement in cosmological constraints compared to or constraints alone. Therefore in the next section, we will examine a joint likelihood between more redshift bins, and where the HI bias () is taken into account.
4.3 HI bias and multi-redshift bin joint likelihood analysis
It is well known that there is a degeneracy between galaxy bias, and in parameter constraints, since these parameters all affect the amplitude of the power spectrum (see Eq. 8). However, they contribute differently to the evolution of the power spectrum with time; hence by measuring the power spectra in various redshifts we can break the degeneracies between them. From Eq. 8, we can see that while measures , additionally measures . Combining the cross-bin intensity mapping power spectra with the cross-spectra we can therefore tighten our constraints on bias and cosmological parameters.
In this section we consider two different bias models, with distinct parameter sets. The first, more restricted model explores bias amplitude variation via the parameter in Eq. 17, setting the rest of the parameters to the best fit values (Cunnington et al., 2019b). In the second model, , and are the bias parameters with set to equal 1. We include these parameters when evaluating the Fisher matrices (Eq. 22). Note that both models are scale-invariant and depend only on . We will consider both full-sky and 300 deg2 surveys to explore the viability of HIHI and HI in cosmological constraints.
For the full-sky case, as we include more parameters for , we also examine more redshift bins for both HI and to obtain the best possible results. We consider the redshift range for which would be measured by pre-SKA and SKA-MID experiments (Square Kilometre Array Cosmology Science Working Group et al., 2020; Pourtsidou et al., 2017; Santos et al., 2015). Table 1 shows the central redshifts we consider for and bins; the width of each bin is 150 Mpc, 0.05. Table 1 lists both HI and central redshifts. As we have 16 , we shall refer to ’16-HIHI’ which corresponds to 16 pairs of HI auto-correlation functions; we cross-correlate HI intensity maps to fields at = 0.44, 0.78 and 1.77 respectively. We refer to 16-HIHI+1- as the joint analysis for 16-HIHI and HI 2-point statistics at which = 0.44. We add further bins and label joint data as 16-HIHI+2- and 16-HIHI+3-. We calculate both the futuristic case where HI can be measured with high and the current state of art where .
We further calculate the figure of merit (FoM) for the constraint. The FoM is the inverse of the area of the contours; in this case we calculate the FoM at 95 confidence level. Fig. 9 shows the FoM of the constraints. The blue dots show the FoM from 16HIHI, where we consecutively add HI auto-correlations for the redshift bins in the order listed in Table 1. We see that all redshift bins contribute to an improved signal, with a nearly linearly increasing contribution (for this experiment, we are assuming that only HI slices are available). The green, red and black dots in Fig. 9 show the FoM for when we further add the cross-correlations between consecutive HI slices and the slices at = 0.44, 0.78 and 1.78, respectively. We see that these cross-correlations significantly improve the FoM, and appear to be converging to a maximal constraint when including all slices.
0.02 | 0.44 |
0.08 | 0.78 |
0.13 | 1.77 |
0.18 | |
0.24 | |
0.29 | |
0.35 | |
0.41 | |
0.47 | |
0.54 | |
0.60 | |
0.68 | |
0.75 | |
0.83 | |
0.91 | |
0.99 |
We now present Fisher forecast results for the first bias model, using the redshift bins in Table 1. Fig. 11 shows the utility of multi-redshift bin power spectra measurements with HI and . With the appropriate redshift bin size of HI () for , we can achieve tight cosmological constraints (in this low noise case) which are comparable to the optical and CMB probes such as (Abbott et al., 2019; Alam et al., 2017; Planck Collaboration et al., 2018; Abbott et al., 2018). By including more redshift slices, the constraints are improved significantly especially for and . However, this also makes the contours more elliptical, as there are remaining degeneracies among parameters. The uncertainties on parameters are measured at 95 confidence level and reported in Table 2.
We next consider the current state of the art case, where since the typical angular resolution , we set = 375 (we choose this particular value as it is convenient to consider the bins as 15 bins with = 25). Fig. 9 shows the cumulative Figure of Merit in this case, while Fig. 11 illustrates the likelihood contours for cosmological parameters with this maximum multipole. Comparing with Fig. 11, where we can see that there is a substantial difference in the 16-HIHI contours. However, we notice the significant improvement in parameter constraints when we add 3 bins to 16-HIHI (green shade) in = 375. We are therefore seeing that by joining HI 2-point statistics to HIHI auto-correlations, we can improve cosmological constraints significantly.
Parameters | HI bias model 1 | HI bias model 2 |
---|---|---|
0.02 | 0.02 | |
0.01 | 0.02 | |
0.03 | 0.04 | |
0.04 | 0.05 | |
0.04 | - | |
- | 0.04 | |
- | 0.03 |
For bias model 2, we find that the second order coefficient of HI bias is very poorly constrained. Marginalising over this parameter does not significantly affect the other cosmological parameter constraints (, , and ). We set this parameter to 0.05 following Cunnington et al. (2019b). Fig. 12 shows the likelihood constraints for this model. We can see that by adding more slices, the constraints do not improve much but the improvement in the constraint can be easily noticed. The uncertainties on our HI bias models and cosmological parameters are reported in Table 2. Comparing the parameter constraints from both models (see Table 2), we notice that model 1 gives slightly better (but very comparable) cosmological constraints.
4.4 The effect of sky coverage
Now let us consider the effect of sky coverage area on 2-point statistics of both HIHI auto and HI cross angular power spectra. We now consider a combined survey area 300 deg2 of lensing and HI observations, comparable to current pathfinder intensity mapping surveys. As this area is much smaller than the full-sky case, we therefore now need to use the pseudo angular power () as an estimator of .
Suppose the survey footprint of the observations can be expressed using the weight function . Normalised by the sky factor (the fraction of the sky covered by the data), the weight moments are given by:
(23) |
where represents the -th moment of weighting. The power spectrum of the window function is:
(24) |
For a spin-0 field weighted by , a spherical harmonic coefficient can be expressed as (Kim & Naselsky, 2010; Kim, 2011):
(25) |
where we approximate the integration over sky factor by the summation over pixel area with the surface density . The pseudo power spectrum estimator, , is then:
(26) |
Similarly for spin-2 fields (), we can obtain the coefficients, , by:
(27) |
where
(28) |
Similarly to the full-sky case, and are then:
(29) |
(30) |
The pseudo power spectra for E and B modes are:
(31) |
The pseudo power spectra and true are related by the mode-mode coupling resulting from masking ():
(32) |
This kernel depends solely on the geometry of a cut-sky and plays a crucial role in the pseudo- method. For details concerning the mode-mode coupling, see Hivon et al. (2002); Alonso et al. (2019).
We utilise NaMaster (Alonso et al., 2019), which is a software package to calculate pseudo- for any spin fields, to evaluate HIHI and HI pseudo- for 300 deg2 of our simulations above, within a masked region RA = [0, 30] deg and dec = +[0,10] deg, for the same redshift bins (see Table 1). We choose the same bins as before up to with deg convolution, and calculate the cut-sky covariance matrix using NaMaster. Fig. 13 shows the cosmological constraints feasible of this scenario. Comparing the results to the full-sky case (Fig. 11), we can see that the statistical incompleteness due to having a cut-sky survey reduces the feasibility; with 300 deg2 sky, we cannot detect the HIHI cosmic signal. However, incorporating the cross-correlation with weak lensing improves the significance of the joint statistics (Fig. 13).
5 Instrument Noise
In this section we consider the instrument noise for both lensing and HI surveys for the current state of art case. We will consider the expected thermal noise for a single-dish survey for HI measurement.
5.1 Single Dish Thermal Noise
We begin by examining the noise on 3D measured power spectra for the HI auto-correlation (). This discussion is based on the works of (Battye et al., 2013; Santos et al., 2015; Bigot-Sazy et al., 2015). Subsequently, we derive the root mean square (rms) thermal noise expected on IM maps from the power spectra, which we will use to assess the effect of realistic noise on the ability to detect the lensing-HI cross-correlation.
The expected uncertainty () on the power spectrum can be estimated by evaluating its expected second moment. By calculating the ratio between and averaging over the radial wavenumber bin size , we can estimate the uncertainty in the IM map measurements. Following this procedure the error on can be estimated by the following expression (Seo et al., 2010; Feldman et al., 1994; Battye et al., 2013):
(33) |
where is the volume of survey, is windows function, is the average temperature at given redshift and is the pixel volume. In this work we shall ignore the contribution of shot noise and only consider the contribution from pixel (thermal) noise . These parameters are discussed in detail in the following paragraphs.
It is essential to pick the proper size of to optimise the viability of a single dish. As we focus on the cosmic signal the acoustic scale should be an aim. That means we require , where is the wavenumber of acoustic scale. The volume of the survey can be computed by
(34) |
where
(35) |
where we assume a flat universe. The window function is set by the instrument’s specification. As we map the from multiple redshift bins, we can ignore the contribution of radial directions in . However for the angular direction, this is not the case. Importantly, the angular resolution of the radio beam will define the noise level. We can model by:
(36) |
where is the full width half max of angular resolution. Realistically, depends on frequency (). However, in this analysis we assume the variation of is small and negligible.
The pixel volume, , can be determined by
(37) |
where is the central redshift and corresponds to frequency width and is the pixel solid angle.
In radio astronomy we normally measure the signal in terms of power. The antenna temperature then generates the thermal noise; the pixel noise can be approximated by (Santos et al., 2015; Seo et al., 2010),
(38) |
where represents the observation time per pointing, (approximately 1) signifies the efficiency of the telescope, meaning that almost no signal is lost when radio radiation is transmitted to the antenna. represents the system temperature, which includes
(39) |
where we ignore the contribution from the Earth’s atmosphere. is the spill over from ground radiation (approximately 3K), 2.73K and galactic temperature 25K (408 MHz/)2.75K (Square Kilometre Array Cosmology Science Working Group et al., 2020). The observing time per pointing relates to the total observation by = , where is an angular pixel size.
As is the rms thermal noise, by definition its square is the power per pixel volume (. Therefore, the 3D noise power spectrum () is then (Battye et al., 2013; Santos et al., 2015)
(40) |
where
(41) |
If we have dishes where each dish has beams, we can take less time for each pointing area. The noise power spectrum is then reduced to
(42) |
To determine the optimisation of a survey strategy, we can estimate the suitable and that minimises for acoustic scale . This acoustic scale can be estimated by following the work of Blake & Glazebrook (2003); Battye et al. (2013),
(43) |
where is the overall amplitude which can be marginalised. The subscript ’ref’ refers to reference cosmological parameters.
We now calculate the thermal noise of a MeerKAT-like instrument = . In this case we consider a single dish telescope consisting of 64 dishes of 13.5m diameter, operating in UHF, L and S-band. The MeerKAT pilot survey by Wang et al. (2021) focuses on L-band from 856 to 1712 MHz with 4096 frequency channels. This pilot survey has 10.5 hours observation time with approximately 200 deg2 observation field (Wang et al., 2021; Cunnington et al., 2022). The summary statistics of this MeerKAT pilot survey are listed in Table 3. Wang et al. (2021) show that for integrated frequency channels, they can achieve thermal noise mK.
0.2 MHz | |
[200,250] | |
[mK] | |
ttot | 10.5 hours |
[0.3885, 0.4623] | |
64 | |
16 mK | |
87500 | |
1.48 deg | |
200 deg2 |
If we use Eq. (38) and (39) together with Table 3, the expected for a single frequency channel of MeerKAT pilot survey is
(44) |
where we assume each dish has equal efficiency =1 and consider only a single frequency channel = 0.2 MHz. If we consider the whole frequency range like (Wang et al., 2021) the 2mK.
We now consider the case where the total observation time = 1000 hours and 250 frequency channels with = 0.2 MHz.
5.2 S/N of
In section 5.1 we have estimated the rms thermal noise for MeerKAT-like surveys similar to the current state-of-the-art. In this section we explore the future case, where the observation time can take longer than MeerKAT’s pilot survey, and we assume a full-sky survey to estimate the best possible S/N for weak lensing-intensity mapping (HI) 2-point statistics. The estimate from Eq. 44 is 15mK for one frequency channel, which is based on the specification of the current MeerKAT survey in Table 3 and Eq. 38; to detect the cross-correlation we should find ways to reduce the as far as possible.
We first model the S/N of HI. We can consider the zero lag noise level for , i.e. where and are measured in the same pixel. There is no reason why the statistical noise of should be correlated with . The noise for the cross-correlation will therefore be proportional to the product of and , where is HI thermal noise which can be estimated by Eq. 44. The rms noise for weak lensing can be estimated by
(45) |
where is the variance of intrinsic galaxy ellipticities and is galaxy number per pixel. For KiDS and DES-like surveys, 0.3. The KiDS DR4 effective galaxy number density is = 0.325 arcmin-2 for the whole redshift range (Heymans et al., 2021; Giblin et al., 2021). For a pixel size 0.252 deg2, we find 0.03.
The S/N of then can be estimated by
(46) |
and we find that the rms lensing convergence signal is similar to the rms noise on 0.25 deg scales (Gatti et al., 2021; Amon et al., 2022).
To estimate the signal of HI intensity mapping, we first recall that the HI brightness temperature fluctuations can be expressed by:
(47) |
where is an average temperature over angular position () for a given and is HI bias. As the power spectrum is the power of the temperature fluctuation , the square root of per volume is effectively the root mean square of the HI true signal (),
(48) |
We consider a survey similar to the MeerKAT pilot survey (0.39 < z < 0.46) over a moderately thick = 0.075 redshift bin, with 0.252 deg2 pixel size, 4000 Mpc3. Assuming the foregrounds and redshift space distortions have been appropriately dealt with, we can use
(49) |
We assume an effective redshift of the survey = 0.42, =1 and = 0.07 mK (Cunnington et al., 2022; Wang et al., 2021; Santos et al., 2015). In this case, the estimation of HI rms is then
(50) |
where we estimate at the = 0.1 Mpc scale. This estimation generally agrees with Table I in Santos et al. (2015) with slightly better signal rms, because Santos et al. (2015) uses smaller channel bins than the thick redshift bin we have here.
is approximately 15 mK for a single frequency channel given Table 3. If we stack over 200 channels and assume 10 hours observing time, then
(51) |
Then using equation (46), the estimation of S/N for HI 2-point statistics for KiDS-like lensing surveys and MeerKAT for = 0.1 Mpc with a pixel size 0.252 deg2 covering 200 deg2 sky, is then
(52) |
This means that by the current state of art, we expected to observe more instrument noise than cosmic signal from HI cross correlations at zero lag. By including cross-correlations at different angular separations, we expect a higher total signal-to-noise.
Note that this estimation is based on the pilot survey by MeerKAT which only contains 64 dishes of 13m diameter and only observes for 10hr; for the full operation of MeerKAT or SKA-Mid, we will have more dishes, longer observation time and more frequency channels. If we assume a longer observation time such as 1000 hours (as is recommended by (Zhang et al., 2023)) and increase the number of frequency channels to 250, the estimation of is then
(53) |
Now we have better S/N by one order of magnitude for a 3000 pixel sky. (If we use the pilot survey footprint ( 200 deg2), then S/N 2.4 for zero lag.)
We confirm this calculation by generating Gaussian random fields for the instrument noise for both and HI fields and add these noise maps using (Eq. 44 and 45) to the simulations described in section 4. We utilise the full sky maps with NSIDE=128. The pixel number for this resolution is = 196608. Hence the estimation of S/N for this configuration is then,
(54) |
However, when we take the LoS foreground subtraction into account, the signal would be reduced by a factor of 3. This means
(55) |
This estimation indicates that we would expect to observe the true signal at zero lag for 1000 hour exposure and large sky surveys such as SKA.
5.3 Fisher Analysis
We now consider the feasibility of HI and HIHI 2-point statistics for cosmological constraints in the presence of these current noise levels. In the S/N analysis of HI (sec. 5.2), we only considered the case where MHz. However to compare results in sec. 4.3, we will adjust to match in Table. 1.
We generate Gaussian white noise fields for both and HI using the previous subsection’s calculated amplitudes. We note again that in this analysis, we consider only the full-sky case. The marginalisations of cosmological parameters are illustrated by Fig. 14. Comparing this result to the no-noise case (see Fig. 11) we can see that the cosmological feasibility of HIHI and HI are reduced significantly. We also show the resulting FoM in Fig. 15 and 2-constraints in Table 4. We note that the constraint shows a precise estimation, although it is degenerate with . Considering Fig. 15, we can see that the maximum of FoM when considering instrument noise is one order of magnitude less than the no-noise case (see Fig. 9). Nevertheless, all FoM plots (see Fig. 9, 9, and 15) indicate that by combining HI likelihoods with HIHI, we can significantly enhance the feasibility of cosmological constraints from the HI intensity maps.
Parameters | HI bias model 1 |
---|---|
0.28 | |
0.09 | |
0.32 | |
0.17 | |
0.34 |
6 Conclusions
In this paper we have studied the 2-point statistics of lensing convergence and HI intensity mapping.
We first presented the theoretical framework for calculating convergence-intensity mapping cross-correlations. Next, by using realisations from an N-body simulation we have emulated HI intensity maps, and have shown that their cross-correlation with convergence maps from these simulations agree with our theoretical cross-correlation calculations.
We proceeded to study the effect of HI foreground removal on the 2-point functions. We model the effect of foreground removal by removing the mean along each line-of-sight, which effectively represents the largest radial mode. from our HI maps, following the method of Cunnington et al. (2019b). We then measure the post-removal cross-power; we find that the foreground removal modestly reduces the HI power spectrum signal, by a factor . In the case of forthcoming HI experiments that will measure HI at < 1, is approximately 2.5 for our catalogues.
In the following section, we utilised the Fisher matrix formalism to forecast best-case cosmological constraints for the convergence-HI probe, for the maximal case of full sky and subdominant telescope noise, but while including foreground removal. We calculated the Fisher matrix for , HIHI and HI 2-point functions by using the measured covariance matrices from Sec 3.
We find that a single redshift slice of the HI intensity map and can constrain cosmological parameters for known bias (see Fig. 7), but when is a further parameter (or several), the few-slice point functions do not sufficiently constrain the cosmological parameters compared to current cosmological surveys such as Planck and DES (Planck Collaboration et al., 2018; Abbott et al., 2019).
Hence, several cross-bin correlations are required in order for this probe to be of interest. In Sec 4.3, we have explored the use of several redshift bins for HI and convergence, together with the effect of on cosmological constraints. We consider both the current state of art where , and the futuristic case where > 1000. Both cases show that a set of 2-point functions constrain the uncertainty in cosmological parameters to a comparable level with current experiments. All FoMs show that by including the cross-correlation of a lensing survey with the 21cm signal, we can improve the HI auto constraints.
We then examined the impact of a cut-sky survey. In this analysis, we evaluated the 2-point statistics based on a 300 deg2 observed patch of sky. Due to the statistical incompleteness, detecting cosmic signals becomes marginal in this context.
In sec. 5.1, we explore the instrument noise affecting lensing and HI galaxy surveys. The thermal noise of a single-dish survey was calculated. In this study, we focus on instruments similar to MeerKAT for radio observations and KIDs-like surveys for optical counterparts. Our analysis demonstrates the feasibility of detecting the HI cross-correlation, provided we have sufficient sky coverage and long exposure times for the radio measurements.
Even though we have shown positive results for 2-point statistics between the field and HI intensity map, there are important caveats that remain to be explored further:
-
•
In this study we have created HI intensity maps based on the assumption that they are linearly biased in relation to overdensity . A more realistic study should construct HI intensity maps by assigning HI mass to simulated halo catalogues. Also in future work the generated HI maps will be compared in detail to real data.
-
•
We have approximated the foreground cleaning by removing the mean fluctuation along the line of sight, which effectively represents the largest radial mode. More detailed emulation of the foreground cleaning will be studied in future.
-
•
The cosmological and constraint predictions have been obtained by using the idealised Fisher matrix analysis; for real data, Markov Chain Monte Carlo methods are required to deal with non-Gaussian likelihoods and realistic degeneracies between parameters.
-
•
Only the CDM model is considered in these parameter constraints. Further extensions (e.g. wCDM) should be considered in further work.
-
•
Only 2-point statistics have been explored in this research. However the study by Schmit et al. (2019) shows that combining the bispectrum and power spectrum can reduce the error of cosmological parameters by an order of magnitude compared to Planck.
-
•
This study has only explored the low-redshift 2-point functions. The high redshift probes at the time of the Epoch of Reionization and the Cosmic Microwave Background are not taken into account. It is an interesting matter for future work to consider whether 2-point statistics between the EOR HI and CMB weak lensing can also be measured.
In conclusion, -HI cross-correlations are an intriguing additional probe for cosmology, which are not destroyed by foreground removal. This probe will be available for measurement with forthcoming HI and lensing surveys this decade.
Acknowledgements
We acknowledge studentship support from the Thai Royal Government for AS. We thank Tzu-Ching Chang (JPL), Yu-Wei Liao (ASIAA), Steve Cunnington (University of Manchester) and Utane Sawangwit (National Astronomical Research Institute of Thailand) for extremely helpful discussions. DB was supported by STFC Consolidator Grant ST/S000550/1. Calculations were made using Chanlawan High Performance Computer at the National Astronomical Research Institute of Thailand.
Data Availability
In this work we use the N-Body simulations and weak lensing catalogues from Takahashi et al. (2017). The catalogues can be obtained at http://cosmo.phys.hirosaki-u.ac.jp/takahasi/allsky_raytracing. The measured 2-point functions and Python codes for calculations in sections 3 and 4 are available in Github, at https://github.com/AnutUoP/2020s-KappaHI. Other data will be shared on reasonable requests to the corresponding author.
References
- Abbott et al. (2018) Abbott T. M. C., et al., 2018, Phys. Rev. D, 98, 043526
- Abbott et al. (2019) Abbott T. M. C., et al., 2019, Astrophys. J., 872, L30
- Alam et al. (2017) Alam S., et al., 2017, MNRAS, 470, 2617
- Alonso et al. (2019) Alonso D., Sanchez J., Slosar A., LSST Dark Energy Science Collaboration 2019, MNRAS, 484, 4127
- Amon et al. (2022) Amon A., et al., 2022, Phys. Rev. D, 105, 023514
- Bartelmann & Schneider (2001) Bartelmann M., Schneider P., 2001, Phys. Rep., 340, 291
- Battye et al. (2013) Battye R. A., et al., 2013, MNRAS, 434, 1239
- Baxter et al. (2019) Baxter E. J., et al., 2019, Phys. Rev., D99, 023508
- Bigot-Sazy et al. (2015) Bigot-Sazy M.-A., et al., 2015, Monthly Notices of the Royal Astronomical Society, 454, 3240
- Blake & Glazebrook (2003) Blake C., Glazebrook K., 2003, ApJ, 594, 665
- Brown et al. (2005) Brown M. L., Castro P. G., Taylor A. N., 2005, MNRAS, 360, 1262
- Bull et al. (2015) Bull P., Ferreira P. G., Patel P., Santos M. G., 2015, ApJ, 803, 21
- CHIME Collaboration et al. (2022) CHIME Collaboration et al., 2022, arXiv e-prints, p. arXiv:2202.01242
- Castro et al. (2005) Castro P. G., Heavens A. F., Kitching T. D., 2005, Phys. Rev. D, 72, 023516
- Chapman et al. (2012) Chapman E., et al., 2012, MNRAS, 423, 2518
- Cunnington et al. (2019a) Cunnington S., Harrison I., Pourtsidou A., Bacon D., 2019a, MNRAS, 482, 3341
- Cunnington et al. (2019b) Cunnington S., et al., 2019b, MNRAS, 488, 5452
- Cunnington et al. (2020) Cunnington S., Pourtsidou A., Soares P. S., Blake C., Bacon D., 2020, MNRAS, 496, 415–433
- Cunnington et al. (2022) Cunnington S., et al., 2022, arXiv e-prints, p. arXiv:2206.01579
- Cunnington et al. (2023) Cunnington S., et al., 2023, MNRAS, 523, 2453
- Euclid Collaboration (2019) Euclid Collaboration 2019, arXiv e-prints, p. arXiv:1910.09273
- Fang et al. (2022) Fang X., Eifler T., Schaan E., Huang H.-J., Krause E., Ferraro S., 2022, MNRAS, 509, 5721
- Feldman et al. (1994) Feldman H. A., Kaiser N., Peacock J. A., 1994, ApJ, 426, 23
- Fernandez-Soto et al. (2001) Fernandez-Soto A., Lanzetta K. M., Chen H.-W., Pascarelle S. M., Yahata N., 2001, The Astrophysical Journal Supplement Series, 135, 41
- Gatti et al. (2021) Gatti M., et al., 2021, MNRAS, 504, 4312
- Giblin et al. (2021) Giblin B., et al., 2021, A&A, 645, A105
- Górski et al. (2005) Górski K. M., et al., 2005, ApJ, 622, 759
- Harker et al. (2010) Harker G., et al., 2010, MNRAS, 405, 2492
- Heavens (2003) Heavens A., 2003, MNRAS, 343, 1327
- Heymans et al. (2021) Heymans C., et al., 2021, A&A, 646, A140
- Hinshaw et al. (2013) Hinshaw G., et al., 2013, ApJS, 208, 19
- Hivon et al. (2002) Hivon E., Górski K. M., Netterfield C. B., Crill B. P., Prunet S., Hansen F., 2002, ApJ, 567, 2
- Hu & Jain (2004) Hu W., Jain B., 2004, Phys. Rev., D70, 043009
- Kamionkowski et al. (2011) Kamionkowski M., Smith T. L., Heavens A., 2011, Phys. Rev. D, 83, 023007
- Kim (2011) Kim J., 2011, A&A, 531, A32
- Kim & Naselsky (2010) Kim J., Naselsky P., 2010, A&A, 519, A104
- Lewis & Bridle (2002) Lewis A., Bridle S., 2002, Phys. Rev. D, 66, 103511
- LoVerde & Afshordi (2008) LoVerde M., Afshordi N., 2008, Phys. Rev. D, 78, 123506
- Mao et al. (2008) Mao Y., et al., 2008, Phys. Rev. D, 78, 023529
- Martin et al. (2012) Martin A. M., Giovanelli R., Haynes M. P., Guzzo L., 2012, The Astrophysical Journal, 750, 38
- Mendez et al. (2014) Mendez R. A., Silva J. F., Orostica R., Lobos R., 2014, PASP, 126, 798
- Padmanabhan et al. (2020) Padmanabhan H., Refregier A., Amara A., 2020, MNRAS, 495, 3935–3942
- Pandey et al. (2021) Pandey S., et al., 2021, arXiv e-prints, p. arXiv:2105.13545
- Planck Collaboration et al. (2018) Planck Collaboration et al., 2018, arXiv e-prints, p. arXiv:1807.06205
- Platania et al. (1998) Platania P., Bensadoun M., Bersanelli M., De Amici G., Kogut A., Levin S., Maino D., Smoot G. F., 1998, ApJ, 505, 473
- Pourtsidou (2017) Pourtsidou A., 2017, arXiv e-prints, p. arXiv:1709.07316
- Pourtsidou et al. (2017) Pourtsidou A., Bacon D., Crittenden R., 2017, MNRAS, 470, 4251
- Pratten et al. (2016) Pratten G., Munshi D., Valageas P., Brax P., 2016, Phys. Rev. D, 93, 103524
- Santos et al. (2010) Santos M. G., Ferramacho L., Silva M. B., Amblard A., Cooray A., 2010, MNRAS, 406, 2421
- Santos et al. (2015) Santos M., et al., 2015, in Advancing Astrophysics with the Square Kilometre Array (AASKA14). p. 19 (arXiv:1501.03989)
- Schmit et al. (2019) Schmit C. J., Heavens A. F., Pritchard J. R., 2019, MNRAS, 483, 4259
- Seo et al. (2010) Seo H.-J., Dodelson S., Marriner J., Mcginnis D., Stebbins A., Stoughton C., Vallinotto A., 2010, ApJ, 721, 164
- Shaw et al. (2014) Shaw J. R., Sigurdson K., Pen U.-L., Stebbins A., Sitwell M., 2014, The Astrophysical Journal, 781, 57
- Smoot & Debono (2017) Smoot G. F., Debono I., 2017, A&A, 597, A136
- Spinelli et al. (2018) Spinelli M., Bernardi G., Santos M. G., 2018, MNRAS, 479, 275–283
- Spinelli et al. (2022) Spinelli M., Carucci I. P., Cunnington S., Harper S. E., Irfan M. O., Fonseca J., Pourtsidou A., Wolz L., 2022, MNRAS, 509, 2048
- Square Kilometre Array Cosmology Science Working Group et al. (2020) Square Kilometre Array Cosmology Science Working Group et al., 2020, Publ. Astron. Soc. Australia, 37, e007
- Su et al. (2018) Su H., et al., 2018, MNRAS, 479, 4041–4055
- Switzer et al. (2013) Switzer E. R., et al., 2013, MNRAS: Letters, 434, L46–L50
- Takahashi et al. (2012) Takahashi R., et al., 2012, ApJ, 761, 152
- Takahashi et al. (2017) Takahashi R., et al., 2017, ApJ, 850, 24
- Tegmark et al. (1997) Tegmark M., Taylor A. N., Heavens A. F., 1997, ApJ, 480, 22
- The CHIME Collaboration et al. (2022) The CHIME Collaboration et al., 2022, arXiv e-prints, p. arXiv:2201.07869
- Tröster et al. (2022) Tröster T., et al., 2022, A&A, 660, A27
- Upham et al. (2019) Upham R. E., Whittaker L., Brown M. L., 2019, MNRAS, 491, 3165–3181
- Wang et al. (2021) Wang J., et al., 2021, MNRAS, 505, 3698
- Wolz et al. (2016) Wolz L., Tonini C., Blake C., Wyithe J. S. B., 2016, MNRAS, 458, 3399
- Wuensche & the BINGO Collaboration (2019) Wuensche C. A., the BINGO Collaboration 2019, in Journal of Physics Conference Series. p. 012002 (arXiv:1803.01644), doi:10.1088/1742-6596/1269/1/012002
- Zhang et al. (2023) Zhang M., Li Y., Zhang J.-F., Zhang X., 2023, MNRAS, 524, 2420