The feasibility of weak lensing and 21cm intensity mapping cross-correlation measurements

Anut Sangka1,2,3 and David Bacon1
1Institute of Cosmology and Gravitation, University of Portsmouth, Portsmouth PO1 3FX, UK
2National Astronomical Research Institute of Thailand, Chiangmai 50180, Thailand,
3 Department of Physics, Faculty of Science, Ubon Ratchatani University, Ubon Ratchatani 34190, Thailand
E-mail: anut.sangka@port.ac.uk
(Accepted XXX. Received YYY; in original form ZZZ)
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 z<1𝑧1z<1italic_z < 1) and lensing convergence (κ𝜅\kappaitalic_κ). 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 similar-to\sim50-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 κ𝜅\kappaitalic_κ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 κΔTHI𝜅Δsubscript𝑇HI\kappa\Delta{T}_{\mathrm{HI}}italic_κ roman_Δ italic_T start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT 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 Universe
pubyear: 2023pagerange: The feasibility of weak lensing and 21cm intensity mapping cross-correlation measurementsThe feasibility of weak lensing and 21cm intensity mapping cross-correlation measurements

1 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, 7.1σ7.1𝜎7.1\sigma7.1 italic_σ (LRG), 5.7σ5.7𝜎5.7\sigma5.7 italic_σ (ELG), and 11.1σ11.1𝜎11.1\sigma11.1 italic_σ (QSO).  Cunnington et al. (2022) have detected the correlated clustering between MeerKAT measurements of HI and galaxies from the WiggleZ Dark Energy Survey at 7.7σ7.7𝜎7.7\sigma7.7 italic_σ 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 (ΔHIδgsubscriptΔHIsubscript𝛿g\Delta_{\mathrm{HI}}\delta_{\mathrm{g}}roman_Δ start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT italic_δ start_POSTSUBSCRIPT roman_g end_POSTSUBSCRIPT) 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 (γ𝛾\gammaitalic_γ) 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 γδHI𝛾subscript𝛿HI\gamma-\delta_{\mathrm{HI}}italic_γ - italic_δ start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT, 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’ κ𝜅\kappaitalic_κ can be considered instead of γ𝛾\gammaitalic_γ as both share the same statistical properties. The 2-point functions between the pairs of κ𝜅\kappaitalic_κ and HI could improve cosmological constraints and break degeneracies such as that between HI bias (bHIsubscript𝑏HIb_{\mathrm{HI}}italic_b start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT) 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 κlimit-from𝜅\kappa-italic_κ -HI correlations.

To achieve this, we will present theoretical and simulation approaches for calculating the κ𝜅\kappaitalic_κ-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 κ𝜅\kappaitalic_κ-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 κ𝜅\kappaitalic_κ-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 (κ𝜅\kappaitalic_κ) and neutral hydrogen intensity maps (HI) in subsection 2.1. We will then discuss the generation of κ𝜅\kappaitalic_κ catalogues and HI modelling from simulations of the matter overdensity δ𝛿\deltaitalic_δ. 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 κ𝜅\kappaitalic_κ and HI temperature fluctuations ΔTHIΔsubscript𝑇HI\mathrm{\Delta}T_{\mathrm{HI}}roman_Δ italic_T start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT. We will then turn to the angular cross-power spectra. We denote κκ𝜅𝜅\kappa\kappaitalic_κ italic_κ as the power spectra between κ𝜅\kappaitalic_κ fields, HIiHIj as the cross-power spectra between HI fields, and κ𝜅\kappaitalic_κHI as the cross-power between κ𝜅\kappaitalic_κ and HI. The dummy indices i𝑖iitalic_i and j𝑗jitalic_j refer to the i𝑖iitalic_ith and j𝑗jitalic_jth redshift bins. We will calculate the lensing convergence in an arbitrary direction on the sky n^^𝑛\hat{n}over^ start_ARG italic_n end_ARG using the Born approximation, projecting the matter overdensity δ𝛿\deltaitalic_δ along an unperturbed ray direction. This can be computed by (Bartelmann & Schneider, 2001)

κ(χs,n^)=3ΩmH022c20χs𝑑χχ(χχ)χδ(n^,χ)a(χ),𝜅subscript𝜒𝑠^𝑛3subscriptΩmsubscriptsuperscript𝐻202superscript𝑐2superscriptsubscript0subscript𝜒𝑠differential-dsuperscript𝜒superscript𝜒𝜒𝜒𝜒𝛿^𝑛superscript𝜒𝑎superscript𝜒\kappa(\chi_{s},\hat{n})=\frac{3\Omega_{\mathrm{m}}H^{2}_{0}}{2c^{2}}\int_{0}^% {\chi_{s}}d\chi^{\prime}\frac{\chi^{\prime}(\chi-\chi)}{\chi}\frac{\delta(\hat% {n},\chi^{\prime})}{a(\chi^{\prime})},italic_κ ( italic_χ start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT , over^ start_ARG italic_n end_ARG ) = divide start_ARG 3 roman_Ω start_POSTSUBSCRIPT roman_m end_POSTSUBSCRIPT italic_H start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG start_ARG 2 italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG ∫ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_χ start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT end_POSTSUPERSCRIPT italic_d italic_χ start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT divide start_ARG italic_χ start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ( italic_χ - italic_χ ) end_ARG start_ARG italic_χ end_ARG divide start_ARG italic_δ ( over^ start_ARG italic_n end_ARG , italic_χ start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ) end_ARG start_ARG italic_a ( italic_χ start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ) end_ARG , (1)

where χ𝜒\chiitalic_χ is comoving distance, ΩmsubscriptΩm\Omega_{\rm m}roman_Ω start_POSTSUBSCRIPT roman_m end_POSTSUBSCRIPT is the matter density parameter at the present epoch, H0subscript𝐻0H_{0}italic_H start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT is the Hubble parameter today and the subscript s𝑠sitalic_s refers to the source plane. For lensing of distributed sources in redshift bins i𝑖iitalic_i, the integrand is modified by including a source distribution, so that the integration now becomes

κi(n^)=0𝑑χqκi(χ)δ(n^,χ),superscript𝜅𝑖^𝑛superscriptsubscript0differential-dsuperscript𝜒subscriptsuperscript𝑞𝑖𝜅superscript𝜒𝛿^𝑛superscript𝜒\kappa^{i}(\hat{n})=\int_{0}^{\infty}d\chi^{\prime}q^{i}_{\kappa}(\chi^{\prime% })\delta(\hat{n},\chi^{\prime}),italic_κ start_POSTSUPERSCRIPT italic_i end_POSTSUPERSCRIPT ( over^ start_ARG italic_n end_ARG ) = ∫ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ∞ end_POSTSUPERSCRIPT italic_d italic_χ start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT italic_q start_POSTSUPERSCRIPT italic_i end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_κ end_POSTSUBSCRIPT ( italic_χ start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ) italic_δ ( over^ start_ARG italic_n end_ARG , italic_χ start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ) , (2)

where the lensing weight is given by

qκi(χ)=3ΩmH022c20χsδ(n^,χ)a(χ)χ𝑑χχχχnsi(z(χ))dzdχn¯si,subscriptsuperscript𝑞𝑖𝜅𝜒3subscriptΩmsubscriptsuperscript𝐻202superscript𝑐2superscriptsubscript0subscript𝜒𝑠𝛿^𝑛superscript𝜒𝑎superscript𝜒superscriptsubscript𝜒differential-dsuperscript𝜒𝜒𝜒𝜒subscriptsuperscript𝑛𝑖𝑠𝑧superscript𝜒𝑑𝑧𝑑superscript𝜒subscriptsuperscript¯𝑛𝑖𝑠q^{i}_{\kappa}(\chi)=\frac{3\Omega_{\mathrm{m}}H^{2}_{0}}{2c^{2}}\int_{0}^{% \chi_{s}}\frac{\delta(\hat{n},\chi^{\prime})}{a(\chi^{\prime})}\int_{\chi}^{% \infty}d\chi^{\prime}\frac{\chi-\chi}{\chi}\frac{n^{i}_{s}(z(\chi^{\prime}))% \frac{dz}{d\chi^{\prime}}}{\bar{n}^{i}_{s}},italic_q start_POSTSUPERSCRIPT italic_i end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_κ end_POSTSUBSCRIPT ( italic_χ ) = divide start_ARG 3 roman_Ω start_POSTSUBSCRIPT roman_m end_POSTSUBSCRIPT italic_H start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG start_ARG 2 italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG ∫ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_χ start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT end_POSTSUPERSCRIPT divide start_ARG italic_δ ( over^ start_ARG italic_n end_ARG , italic_χ start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ) end_ARG start_ARG italic_a ( italic_χ start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ) end_ARG ∫ start_POSTSUBSCRIPT italic_χ end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ∞ end_POSTSUPERSCRIPT italic_d italic_χ start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT divide start_ARG italic_χ - italic_χ end_ARG start_ARG italic_χ end_ARG divide start_ARG italic_n start_POSTSUPERSCRIPT italic_i end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT ( italic_z ( italic_χ start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ) ) divide start_ARG italic_d italic_z end_ARG start_ARG italic_d italic_χ start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT end_ARG end_ARG start_ARG over¯ start_ARG italic_n end_ARG start_POSTSUPERSCRIPT italic_i end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT end_ARG , (3)

where nsi(z)subscriptsuperscript𝑛𝑖𝑠𝑧n^{i}_{s}(z)italic_n start_POSTSUPERSCRIPT italic_i end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT ( italic_z ) is the lensing source number density, and n¯sisubscriptsuperscript¯𝑛𝑖𝑠\bar{n}^{i}_{s}over¯ start_ARG italic_n end_ARG start_POSTSUPERSCRIPT italic_i end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT is its average in the i𝑖iitalic_ith redshift bin.

HI will be a biased tracer of matter overdensity, so we write ΔTHIΔsubscript𝑇HI\mathrm{\Delta}T_{\mathrm{HI}}roman_Δ italic_T start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT (n^,z)^𝑛𝑧(\hat{n},z)( over^ start_ARG italic_n end_ARG , italic_z ) = T¯HI(z)bHI(z)δ(n^,z)subscript¯𝑇HI𝑧subscript𝑏HI𝑧𝛿^𝑛𝑧\bar{T}_{\mathrm{HI}}(z)b_{\mathrm{HI}}(z)\delta(\hat{n},z)over¯ start_ARG italic_T end_ARG start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT ( italic_z ) italic_b start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT ( italic_z ) italic_δ ( over^ start_ARG italic_n end_ARG , italic_z ), where bHI(z)subscript𝑏HI𝑧b_{\mathrm{HI}}(z)italic_b start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT ( italic_z ) is the HI bias at a given redshift z𝑧zitalic_z and T¯HI(z)subscript¯𝑇HI𝑧\bar{T}_{\mathrm{HI}}(z)over¯ start_ARG italic_T end_ARG start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT ( italic_z ) is the average temperature. The projected temperature fluctuation at the i𝑖iitalic_ith redshift bin is then

ΔTHIi(n^)=0χi𝑑χqHIi(χ)δ(n^,χ),Δsuperscriptsubscript𝑇HI𝑖^𝑛superscriptsubscript0subscript𝜒𝑖differential-dsuperscript𝜒subscriptsuperscript𝑞𝑖HIsuperscript𝜒𝛿^𝑛superscript𝜒\mathrm{\Delta}T_{\mathrm{HI}}^{i}(\hat{n})=\int_{0}^{\chi_{i}}d\chi^{\prime}q% ^{i}_{\mathrm{HI}}(\chi^{\prime})\delta(\hat{n},\chi^{\prime}),roman_Δ italic_T start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_i end_POSTSUPERSCRIPT ( over^ start_ARG italic_n end_ARG ) = ∫ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_χ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT end_POSTSUPERSCRIPT italic_d italic_χ start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT italic_q start_POSTSUPERSCRIPT italic_i end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT ( italic_χ start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ) italic_δ ( over^ start_ARG italic_n end_ARG , italic_χ start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ) , (4)

where

qHIi(χ)=T¯HI(z(χ))bHIi(χ)nHIi(z(χ))dzdχn¯HIi,subscriptsuperscript𝑞𝑖HI𝜒subscript¯𝑇HI𝑧𝜒subscriptsuperscript𝑏𝑖HI𝜒subscriptsuperscript𝑛𝑖HI𝑧𝜒𝑑𝑧𝑑𝜒subscriptsuperscript¯𝑛𝑖HIq^{i}_{\mathrm{HI}}(\chi)=\bar{T}_{\mathrm{HI}}(z(\chi))b^{i}_{\mathrm{HI}}(% \chi)\frac{n^{i}_{\mathrm{HI}}(z(\chi))\frac{dz}{d\chi}}{\bar{n}^{i}_{\mathrm{% HI}}},italic_q start_POSTSUPERSCRIPT italic_i end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT ( italic_χ ) = over¯ start_ARG italic_T end_ARG start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT ( italic_z ( italic_χ ) ) italic_b start_POSTSUPERSCRIPT italic_i end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT ( italic_χ ) divide start_ARG italic_n start_POSTSUPERSCRIPT italic_i end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT ( italic_z ( italic_χ ) ) divide start_ARG italic_d italic_z end_ARG start_ARG italic_d italic_χ end_ARG end_ARG start_ARG over¯ start_ARG italic_n end_ARG start_POSTSUPERSCRIPT italic_i end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT end_ARG , (5)

where nHIi(z)subscriptsuperscript𝑛𝑖HI𝑧n^{i}_{\mathrm{HI}}(z)italic_n start_POSTSUPERSCRIPT italic_i end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT ( italic_z ) is the HI source number density, and n¯HIisubscriptsuperscript¯𝑛𝑖HI\bar{n}^{i}_{\mathrm{HI}}over¯ start_ARG italic_n end_ARG start_POSTSUPERSCRIPT italic_i end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT is its average in the i𝑖iitalic_ith redshift bin.

Battye et al. (2013) show that for a given redshift z𝑧zitalic_z, T¯HI(z)subscript¯𝑇HI𝑧\bar{T}_{\mathrm{HI}}(z)over¯ start_ARG italic_T end_ARG start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT ( italic_z ) can be estimated by

T¯HI(z)=44μK(ΩHI(z)h2.45×104)(1+z)2E(z),subscript¯𝑇HI𝑧44𝜇KsubscriptΩHI𝑧2.45superscript104superscript1𝑧2𝐸𝑧\bar{T}_{\mathrm{HI}}(z)=44\mu\mathrm{K}\bigg{(}\frac{\mathrm{\Omega}_{\mathrm% {HI}}(z)h}{2.45\times 10^{-4}}\bigg{)}\frac{(1+z)^{2}}{E(z)},over¯ start_ARG italic_T end_ARG start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT ( italic_z ) = 44 italic_μ roman_K ( divide start_ARG roman_Ω start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT ( italic_z ) italic_h end_ARG start_ARG 2.45 × 10 start_POSTSUPERSCRIPT - 4 end_POSTSUPERSCRIPT end_ARG ) divide start_ARG ( 1 + italic_z ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG italic_E ( italic_z ) end_ARG , (6)

where E(z)=H(z)/H0𝐸𝑧𝐻𝑧subscript𝐻0E(z)=H(z)/H_{0}italic_E ( italic_z ) = italic_H ( italic_z ) / italic_H start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT is the dimensionless Hubble function at redshift z𝑧zitalic_z. The HI density parameter could be approximated to be ΩHIh=2.45×104subscriptΩHI2.45superscript104\mathrm{\Omega}_{\mathrm{HI}}h=2.45\times 10^{-4}roman_Ω start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT italic_h = 2.45 × 10 start_POSTSUPERSCRIPT - 4 end_POSTSUPERSCRIPT  (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)

ΩHI(z)=0.00048+0.00039z0.000065z2.subscriptΩHI𝑧0.000480.00039𝑧0.000065superscript𝑧2\mathrm{\Omega}_{\mathrm{HI}}(z)=0.00048+0.00039z-0.000065z^{2}.roman_Ω start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT ( italic_z ) = 0.00048 + 0.00039 italic_z - 0.000065 italic_z start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT . (7)

Constraining the HI bias bHI(z)subscript𝑏HI𝑧b_{\mathrm{HI}}(z)italic_b start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT ( italic_z ) will be discussed later in section 4.

Using the Limber approximation, the angular power spectra CXY()superscript𝐶𝑋𝑌C^{XY}(\ell)italic_C start_POSTSUPERSCRIPT italic_X italic_Y end_POSTSUPERSCRIPT ( roman_ℓ ) are given by

CXY()=𝑑χqX(χ)qY(χ)χ2Pδ(+1/2χ,z(χ)),superscript𝐶𝑋𝑌differential-d𝜒subscript𝑞𝑋𝜒subscript𝑞𝑌𝜒superscript𝜒2subscript𝑃𝛿12𝜒𝑧𝜒C^{XY}(\ell)=\int d\chi\frac{q_{X}(\chi)q_{Y}(\chi)}{\chi^{2}}P_{\delta}\bigg{% (}\frac{\ell+1/2}{\chi},z(\chi)\bigg{)},italic_C start_POSTSUPERSCRIPT italic_X italic_Y end_POSTSUPERSCRIPT ( roman_ℓ ) = ∫ italic_d italic_χ divide start_ARG italic_q start_POSTSUBSCRIPT italic_X end_POSTSUBSCRIPT ( italic_χ ) italic_q start_POSTSUBSCRIPT italic_Y end_POSTSUBSCRIPT ( italic_χ ) end_ARG start_ARG italic_χ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG italic_P start_POSTSUBSCRIPT italic_δ end_POSTSUBSCRIPT ( divide start_ARG roman_ℓ + 1 / 2 end_ARG start_ARG italic_χ end_ARG , italic_z ( italic_χ ) ) , (8)

where Pδ(+1/2χ,z(χ))subscript𝑃𝛿12𝜒𝑧𝜒P_{\delta}\bigg{(}\frac{\ell+1/2}{\chi},z(\chi)\bigg{)}italic_P start_POSTSUBSCRIPT italic_δ end_POSTSUBSCRIPT ( divide start_ARG roman_ℓ + 1 / 2 end_ARG start_ARG italic_χ end_ARG , italic_z ( italic_χ ) ) 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 ΔχΔ𝜒\rm\Delta{\chi}roman_Δ italic_χ = 150 h1superscript1h^{-1}italic_h start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPTMpc (see section 2.2), we choose a radial selection function for nHIi(z(χ))/n¯HIisubscriptsuperscript𝑛𝑖HI𝑧𝜒subscriptsuperscript¯𝑛𝑖HIn^{i}_{\mathrm{HI}}(z(\chi))/\bar{n}^{i}_{\mathrm{HI}}italic_n start_POSTSUPERSCRIPT italic_i end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT ( italic_z ( italic_χ ) ) / over¯ start_ARG italic_n end_ARG start_POSTSUPERSCRIPT italic_i end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT as a normal distribution around a central comoving position with 3σχ3subscript𝜎𝜒3\sigma_{\chi}3 italic_σ start_POSTSUBSCRIPT italic_χ end_POSTSUBSCRIPT = 150 h1superscript1h^{-1}italic_h start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPTMpc. This σχsubscript𝜎𝜒\sigma_{\chi}italic_σ start_POSTSUBSCRIPT italic_χ end_POSTSUBSCRIPT corresponds to the frequency bandwidth (Δν)Δ𝜈(\Delta\nu)( roman_Δ italic_ν ) 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) .

We show examples for the first time of calculations of the auto- and cross-power for HI and convergence in Fig. 1 using Eq. 8. As expected, the auto- signal depends on the radial HI width σχsubscript𝜎𝜒\sigma_{\chi}italic_σ start_POSTSUBSCRIPT italic_χ end_POSTSUBSCRIPT, while the cross-power is insensitive to this.

Refer to caption
Figure 1: Power spectra Csubscript𝐶C_{\ell}italic_C start_POSTSUBSCRIPT roman_ℓ end_POSTSUBSCRIPT for HI and cross-power between HI and convergence, for radial HI width σχsubscript𝜎𝜒\sigma_{\chi}italic_σ start_POSTSUBSCRIPT italic_χ end_POSTSUBSCRIPT = 50 and 150 h1superscript1h^{-1}italic_h start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPTMpc. The effect of the width is less important for the cross-correlation. Different σχsubscript𝜎𝜒\sigma_{\chi}italic_σ start_POSTSUBSCRIPT italic_χ end_POSTSUBSCRIPT corresponds to different frequency bandwidths, ΔνΔ𝜈\Delta{\nu}roman_Δ italic_ν of the radio data.

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 L𝐿Litalic_L = 450, 900, 1350, …, 6300 h1superscript1h^{-1}italic_h start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPTMpc are nested to represent a region of the Universe in which lensing occurs; each box contains 20483superscript204832048^{3}2048 start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT particles. The κ𝜅\kappaitalic_κ fields are obtained by tracing the light ray path through planes with separation 150 h1superscript1h^{-1}italic_h start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPTMpc. By calculating the Jacobian matrix A𝐴Aitalic_A along the light path, the lensing convergence κ𝜅\kappaitalic_κ, shear lensing γ1,2subscript𝛾12\gamma_{1,2}italic_γ start_POSTSUBSCRIPT 1 , 2 end_POSTSUBSCRIPT and rotation angle ω𝜔\omegaitalic_ω can be obtained, via

A=[1κγ1γ2ωγ2+ω1κ+γ1].𝐴matrix1𝜅subscript𝛾1subscript𝛾2𝜔subscript𝛾2𝜔1𝜅subscript𝛾1\displaystyle A=\begin{bmatrix}1-\kappa-\gamma_{1}&-\gamma_{2}-\omega\\ -\gamma_{2}+\omega&1-\kappa+\gamma_{1}\end{bmatrix}.italic_A = [ start_ARG start_ROW start_CELL 1 - italic_κ - italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT end_CELL start_CELL - italic_γ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT - italic_ω end_CELL end_ROW start_ROW start_CELL - italic_γ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT + italic_ω end_CELL start_CELL 1 - italic_κ + italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT end_CELL end_ROW end_ARG ] . (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 z𝑧absentz\approxitalic_z ≈ 0.78, for which the lensing will significantly occur at the redshift of an intensity map at redshift0.3similar-to-or-equalsabsent0.3\simeq 0.3≃ 0.3. 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 MHIsubscript𝑀HIM_{\mathrm{HI}}italic_M start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT from halo catalogues, we assume that HI is a biased tracer of the total matter overdensity field δ(θ,z)𝛿𝜃𝑧\delta(\theta,z)italic_δ ( italic_θ , italic_z )(see Eq. 4 and 5),

δHI(n^,z)THI(n^,z)T¯HI(z)T¯HI(z)=bHI(z)δ(n^,z),subscript𝛿HI^𝑛𝑧subscript𝑇HI^𝑛𝑧subscript¯𝑇HI𝑧subscript¯𝑇HI𝑧subscript𝑏HI𝑧𝛿^𝑛𝑧\delta_{\mathrm{HI}}(\hat{n},z)\equiv\frac{T_{\mathrm{HI}}(\hat{n},z)-\bar{T}_% {\mathrm{HI}}(z)}{\bar{T}_{\mathrm{HI}}(z)}=b_{\mathrm{HI}}(z)\delta(\hat{n},z),italic_δ start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT ( over^ start_ARG italic_n end_ARG , italic_z ) ≡ divide start_ARG italic_T start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT ( over^ start_ARG italic_n end_ARG , italic_z ) - over¯ start_ARG italic_T end_ARG start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT ( italic_z ) end_ARG start_ARG over¯ start_ARG italic_T end_ARG start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT ( italic_z ) end_ARG = italic_b start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT ( italic_z ) italic_δ ( over^ start_ARG italic_n end_ARG , italic_z ) , (10)

where bHIsubscript𝑏HIb_{\mathrm{HI}}italic_b start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT is a HI bias. For instance the parametric form for bHIsubscript𝑏HIb_{\mathrm{HI}}italic_b start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT adopted by  Cunnington et al. (2019b) is

bHI(z)=0.67+0.18z+0.05z2.subscript𝑏HI𝑧0.670.18𝑧0.05superscript𝑧2b_{\mathrm{HI}}(z)=0.67+0.18z+0.05z^{2}.italic_b start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT ( italic_z ) = 0.67 + 0.18 italic_z + 0.05 italic_z start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT . (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 ΔTHIΔsubscript𝑇HI\mathrm{\Delta}T_{\mathrm{HI}}roman_Δ italic_T start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT :

ΔTHI(n^,z)=THI(n^,z)T¯HI(z)=T¯HI(z)bHI(z)δ(n^,z).Δsubscript𝑇HI^𝑛𝑧subscript𝑇HI^𝑛𝑧subscript¯𝑇HI𝑧subscript¯𝑇HI𝑧subscript𝑏HI𝑧𝛿^𝑛𝑧\mathrm{\Delta}T_{\mathrm{HI}}(\hat{n},z)=T_{\mathrm{HI}}(\hat{n},z)-\bar{T}_{% \mathrm{HI}}(z)=\bar{T}_{\mathrm{HI}}(z)b_{\mathrm{HI}}(z)\delta(\hat{n},z).roman_Δ italic_T start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT ( over^ start_ARG italic_n end_ARG , italic_z ) = italic_T start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT ( over^ start_ARG italic_n end_ARG , italic_z ) - over¯ start_ARG italic_T end_ARG start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT ( italic_z ) = over¯ start_ARG italic_T end_ARG start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT ( italic_z ) italic_b start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT ( italic_z ) italic_δ ( over^ start_ARG italic_n end_ARG , italic_z ) . (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 ×\times× 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 Θ(n^,z)Θ^𝑛𝑧\Theta(\hat{n},z)roman_Θ ( over^ start_ARG italic_n end_ARG , italic_z ) in spherical harmonics:

Θ(n^,z)==0m=m=am(z)Ym(n^),Θ^𝑛𝑧superscriptsubscript0superscriptsubscript𝑚𝑚subscript𝑎𝑚𝑧subscriptsuperscript𝑌𝑚^𝑛\Theta(\hat{n},z)=\sum_{\ell=0}^{\infty}\sum_{m=-\ell}^{m=\ell}a_{\ell m}(z)Y^% {m}_{\ell}(\hat{n}),roman_Θ ( over^ start_ARG italic_n end_ARG , italic_z ) = ∑ start_POSTSUBSCRIPT roman_ℓ = 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ∞ end_POSTSUPERSCRIPT ∑ start_POSTSUBSCRIPT italic_m = - roman_ℓ end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_m = roman_ℓ end_POSTSUPERSCRIPT italic_a start_POSTSUBSCRIPT roman_ℓ italic_m end_POSTSUBSCRIPT ( italic_z ) italic_Y start_POSTSUPERSCRIPT italic_m end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_ℓ end_POSTSUBSCRIPT ( over^ start_ARG italic_n end_ARG ) , (13)

where Ym(n^)subscriptsuperscript𝑌𝑚^𝑛Y^{m}_{\ell}(\hat{n})italic_Y start_POSTSUPERSCRIPT italic_m end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_ℓ end_POSTSUBSCRIPT ( over^ start_ARG italic_n end_ARG ) and am(z)subscript𝑎𝑚𝑧a_{\ell m}(z)italic_a start_POSTSUBSCRIPT roman_ℓ italic_m end_POSTSUBSCRIPT ( italic_z ) are spherical harmonics and their coefficients respectively (Pratten et al., 2016; Castro et al., 2005; Heavens, 2003). Θ(n^,z)Θ^𝑛𝑧\Theta(\hat{n},z)roman_Θ ( over^ start_ARG italic_n end_ARG , italic_z ) 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 amsubscript𝑎𝑚a_{\ell m}italic_a start_POSTSUBSCRIPT roman_ℓ italic_m end_POSTSUBSCRIPT over m𝑚mitalic_m modes:

CXY()=amX(z1)amY(z2),superscript𝐶𝑋𝑌delimited-⟨⟩subscriptsuperscript𝑎𝑋𝑚subscript𝑧1subscriptsuperscript𝑎𝑌𝑚subscript𝑧2C^{XY}(\ell)=\langle a^{X}_{\ell m}(z_{1})a^{Y*}_{\ell m}(z_{2})\rangle,italic_C start_POSTSUPERSCRIPT italic_X italic_Y end_POSTSUPERSCRIPT ( roman_ℓ ) = ⟨ italic_a start_POSTSUPERSCRIPT italic_X end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_ℓ italic_m end_POSTSUBSCRIPT ( italic_z start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ) italic_a start_POSTSUPERSCRIPT italic_Y ∗ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_ℓ italic_m end_POSTSUBSCRIPT ( italic_z start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ) ⟩ , (14)

where X𝑋Xitalic_X and Y𝑌Yitalic_Y stand for the cosmological fields at given redshifts z1subscript𝑧1z_{1}italic_z start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT and z2subscript𝑧2z_{2}italic_z start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT, respectively.

The cross power-spectrum for HI and lensing κ𝜅\kappaitalic_κ 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-Csubscript𝐶C_{\ell}italic_C start_POSTSUBSCRIPT roman_ℓ end_POSTSUBSCRIPT methods are required (Brown et al., 2005; Upham et al., 2019)).

Using this routine, we obtain cross-power measurements for the HI and κ𝜅\kappaitalic_κ 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 𝐂𝐎𝐕𝐂𝐎𝐕\mathbf{COV}bold_COV(CXYsuperscript𝐶𝑋𝑌C^{XY}italic_C start_POSTSUPERSCRIPT italic_X italic_Y end_POSTSUPERSCRIPT) of 2-point statistics from 35 realisations. The error bars are the square root of the diagonal elements of 𝐂𝐎𝐕𝐂𝐎𝐕\mathbf{COV}bold_COV(CXYsuperscript𝐶𝑋𝑌C^{XY}italic_C start_POSTSUPERSCRIPT italic_X italic_Y end_POSTSUPERSCRIPT) of the estimators. The correlation matrix for 𝐂𝐎𝐕𝐂𝐎𝐕\mathbf{COV}bold_COV(CXYsuperscript𝐶𝑋𝑌C^{XY}italic_C start_POSTSUPERSCRIPT italic_X italic_Y end_POSTSUPERSCRIPT) 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 nHIi(z(χ))/n¯HIisubscriptsuperscript𝑛𝑖HI𝑧𝜒subscriptsuperscript¯𝑛𝑖HIn^{i}_{\mathrm{HI}}(z(\chi))/\bar{n}^{i}_{\mathrm{HI}}italic_n start_POSTSUPERSCRIPT italic_i end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT ( italic_z ( italic_χ ) ) / over¯ start_ARG italic_n end_ARG start_POSTSUPERSCRIPT italic_i end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT successfully match the simulations. Due to the lens shell approximation of the ray tracing code, the measured Csubscript𝐶C_{\ell}italic_C start_POSTSUBSCRIPT roman_ℓ end_POSTSUBSCRIPT is slightly affected at very high \ellroman_ℓ (see red line on Fig. 3).

Refer to caption
Refer to caption
Refer to caption
Figure 2: Top: the uncleaned intensity map, ΔTHIorigΔsuperscriptsubscript𝑇HIorig\Delta{T}_{\mathrm{HI}}^{\mathrm{orig}}roman_Δ italic_T start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_orig end_POSTSUPERSCRIPT, at z=0.3𝑧0.3z=0.3italic_z = 0.3 from an example realisation. The fluctuations were measured by assuming bHI(z)subscript𝑏HI𝑧b_{\mathrm{HI}}(z)italic_b start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT ( italic_z ) (see Eq. 12 and 17). Middle: the foreground-removed intensity map, ΔTHIcleanΔsuperscriptsubscript𝑇HIclean\Delta{T}_{\mathrm{HI}}^{\mathrm{clean}}roman_Δ italic_T start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_clean end_POSTSUPERSCRIPT, at the same redshift. The foregrounds were removed by eliminating radial long wavelength modes up to redshift zmax=1subscript𝑧max1z_{\mathrm{max}}=1italic_z start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT = 1. The NSIDEs of the fluctuation maps is reduced from 4096 to 512 to match the resolution of forthcoming radio surveys. Bottom: residual map of cleaned and uncleaned maps. Each of these detail maps has area 5 ×\times× 5 square degrees (a small patch of the entire sky maps).
Refer to caption
Figure 3: Comparison between theoretical Csubscript𝐶C_{\ell}italic_C start_POSTSUBSCRIPT roman_ℓ end_POSTSUBSCRIPT (see Eq. 8) and measured Csubscript𝐶C_{\ell}italic_C start_POSTSUBSCRIPT roman_ℓ end_POSTSUBSCRIPT from our simulations. Here the lensing convergence is measured at central redshift 0.78 and HI is measured at central redshift 0.3.

3 HI foreground removal and its effect on κ𝜅\kappaitalic_κ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 (Tsyncsubscript𝑇syncT_{\mathrm{sync}}italic_T start_POSTSUBSCRIPT roman_sync end_POSTSUBSCRIPT) emission temperature which can be modelled by Tsysnc(1+z)2.7[K]proportional-tosubscript𝑇sysncsuperscript1𝑧2.7delimited-[]KT_{\mathrm{sysnc}}\propto(1+z)^{2.7}[\mathrm{K}]italic_T start_POSTSUBSCRIPT roman_sysnc end_POSTSUBSCRIPT ∝ ( 1 + italic_z ) start_POSTSUPERSCRIPT 2.7 end_POSTSUPERSCRIPT [ roman_K ]  (Platania et al., 1998; Smoot & Debono, 2017), is approximately 3 to 4 orders of magnitude larger than THIsubscript𝑇HIT_{\mathrm{HI}}italic_T start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT 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 ΔTHILoS(n^)Δsuperscriptsubscript𝑇HILoS^𝑛\mathrm{\Delta}T_{\mathrm{HI}}^{\mathrm{LoS}}(\hat{n})roman_Δ italic_T start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_LoS end_POSTSUPERSCRIPT ( over^ start_ARG italic_n end_ARG ).

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 κ𝜅\kappaitalic_κ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 κ𝜅\kappaitalic_κHI cross-correlation survives foreground removal.

Here we follow the method for foreground removal emulation by Cunnington et al. (2019b). The cleaned intensity map ΔTHIcleanΔsuperscriptsubscript𝑇HIclean\mathrm{\Delta}T_{\mathrm{HI}}^{\mathrm{clean}}roman_Δ italic_T start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_clean end_POSTSUPERSCRIPT can be approximated as

ΔTHIclean(n^,z)=ΔTHIorig(n^,z)ΔTHILoS(n^),Δsuperscriptsubscript𝑇HIclean^𝑛𝑧Δsuperscriptsubscript𝑇HIorig^𝑛𝑧Δsuperscriptsubscript𝑇HILoS^𝑛\mathrm{\Delta}T_{\mathrm{HI}}^{\mathrm{clean}}(\hat{n},z)=\mathrm{\Delta}T_{% \mathrm{HI}}^{\mathrm{orig}}(\hat{n},z)-\mathrm{\Delta}T_{\mathrm{HI}}^{% \mathrm{LoS}}(\hat{n}),roman_Δ italic_T start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_clean end_POSTSUPERSCRIPT ( over^ start_ARG italic_n end_ARG , italic_z ) = roman_Δ italic_T start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_orig end_POSTSUPERSCRIPT ( over^ start_ARG italic_n end_ARG , italic_z ) - roman_Δ italic_T start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_LoS end_POSTSUPERSCRIPT ( over^ start_ARG italic_n end_ARG ) , (15)

where ΔTHIorig(n^,z)Δsuperscriptsubscript𝑇HIorig^𝑛𝑧\mathrm{\Delta}T_{\mathrm{HI}}^{\mathrm{orig}}(\hat{n},z)roman_Δ italic_T start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_orig end_POSTSUPERSCRIPT ( over^ start_ARG italic_n end_ARG , italic_z ) is the uncleaned signal in direction n^^𝑛\hat{n}over^ start_ARG italic_n end_ARG at redshift z𝑧zitalic_z. ΔTHILoS(n^)Δsuperscriptsubscript𝑇HILoS^𝑛\mathrm{\Delta}T_{\mathrm{HI}}^{\mathrm{LoS}}(\hat{n})roman_Δ italic_T start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_LoS end_POSTSUPERSCRIPT ( over^ start_ARG italic_n end_ARG ) is defined by:

ΔTHILoS(n^)=1NziT¯HI(zi)bHI(zi)δ(n^,zi),Δsuperscriptsubscript𝑇HILoS^𝑛1subscript𝑁𝑧subscript𝑖subscript¯𝑇HIsubscript𝑧𝑖subscript𝑏HIsubscript𝑧𝑖𝛿^𝑛subscript𝑧𝑖\mathrm{\Delta}T_{\mathrm{HI}}^{\mathrm{LoS}}(\hat{n})=\frac{1}{N_{z}}\sum_{i}% \bar{T}_{\mathrm{HI}}(z_{i})b_{\mathrm{HI}}(z_{i})\delta(\hat{n},z_{i}),roman_Δ italic_T start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_LoS end_POSTSUPERSCRIPT ( over^ start_ARG italic_n end_ARG ) = divide start_ARG 1 end_ARG start_ARG italic_N start_POSTSUBSCRIPT italic_z end_POSTSUBSCRIPT end_ARG ∑ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT over¯ start_ARG italic_T end_ARG start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT ( italic_z start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) italic_b start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT ( italic_z start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) italic_δ ( over^ start_ARG italic_n end_ARG , italic_z start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) , (16)

so that ΔTHILoS(n^)Δsuperscriptsubscript𝑇HILoS^𝑛\mathrm{\Delta}T_{\mathrm{HI}}^{\mathrm{LoS}}(\hat{n})roman_Δ italic_T start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_LoS end_POSTSUPERSCRIPT ( over^ start_ARG italic_n end_ARG ) 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 ΔTHILoS(n^)Δsuperscriptsubscript𝑇HILoS^𝑛\mathrm{\Delta}T_{\mathrm{HI}}^{\mathrm{LoS}}(\hat{n})roman_Δ italic_T start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_LoS end_POSTSUPERSCRIPT ( over^ start_ARG italic_n end_ARG ). 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 1/f1𝑓1/f1 / italic_f 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 1/f1𝑓1/f1 / italic_f 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):

bHI(z)=α(b0+b1z+b2z2),subscript𝑏HI𝑧𝛼subscript𝑏0subscript𝑏1𝑧subscript𝑏2superscript𝑧2b_{\mathrm{HI}}(z)=\alpha(b_{0}+b_{1}z+b_{2}z^{2}),italic_b start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT ( italic_z ) = italic_α ( italic_b start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT + italic_b start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT italic_z + italic_b start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT italic_z start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) , (17)

where α𝛼\alphaitalic_α, b0subscript𝑏0b_{0}italic_b start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT, b1subscript𝑏1b_{1}italic_b start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT and b2subscript𝑏2b_{2}italic_b start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT 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 h1superscript1h^{-1}italic_h start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPTMpc, which are our main interest.

We measure ΔTHILOSΔsuperscriptsubscript𝑇HILOS\mathrm{\Delta}T_{\mathrm{HI}}^{\mathrm{LOS}}roman_Δ italic_T start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_LOS end_POSTSUPERSCRIPT with two choices of maximum redshift, zmaxsubscript𝑧maxz_{\mathrm{max}}italic_z start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT = 1 and 3. zmaxsubscript𝑧maxz_{\mathrm{max}}italic_z start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT = 3 corresponds to futuristic HI-galaxy surveys (Square Kilometre Array Cosmology Science Working Group et al., 2020). On the other hand zmaxsubscript𝑧maxz_{\mathrm{max}}italic_z start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT = 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 (ΔTHIΔTHIΔsubscript𝑇HIΔsubscript𝑇HI\Delta{T}_{\mathrm{HI}}\Delta{T}_{\mathrm{HI}}roman_Δ italic_T start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT roman_Δ italic_T start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT), and the cross-power spectra between HI and κ𝜅\kappaitalic_κ (κΔTHI𝜅Δsubscript𝑇HI\kappa\Delta{T}_{\mathrm{HI}}italic_κ roman_Δ italic_T start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT). We compare the signal of removed and unremoved κΔTHI𝜅Δsubscript𝑇HI\kappa\Delta{T}_{\mathrm{HI}}italic_κ roman_Δ italic_T start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT, 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 >1010\ell>10roman_ℓ > 10, the foreground removal does not erase the κΔTHI𝜅Δsubscript𝑇HI\kappa\Delta{T}_{\mathrm{HI}}italic_κ roman_Δ italic_T start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT power spectrum; the signal is scaled down by a factor (Acleansubscript𝐴cleanA_{\mathrm{clean}}italic_A start_POSTSUBSCRIPT roman_clean end_POSTSUBSCRIPT) which is close to constant over a range of \ellroman_ℓ 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 >1010\ell>10roman_ℓ > 10 . We describe the mean signal drop Acleansubscript𝐴cleanA_{\mathrm{clean}}italic_A start_POSTSUBSCRIPT roman_clean end_POSTSUBSCRIPT by:

Aclean(zHI,zκ,zmax)κΔTHIuncleanedκΔTHIcleaned10<<1500.subscript𝐴cleansubscript𝑧HIsubscript𝑧𝜅subscript𝑧maxsubscriptdelimited-⟨⟩𝜅Δsuperscriptsubscript𝑇HIuncleaned𝜅Δsuperscriptsubscript𝑇HIcleaned101500A_{\mathrm{clean}}(z_{\mathrm{HI}},z_{\kappa},z_{\mathrm{max}})\equiv\bigg{% \langle}\frac{\kappa\Delta{T}_{\mathrm{HI}}^{\mathrm{uncleaned}}}{\kappa\Delta% {T}_{\mathrm{HI}}^{\mathrm{cleaned}}}\bigg{\rangle}_{10<\ell<1500}.italic_A start_POSTSUBSCRIPT roman_clean end_POSTSUBSCRIPT ( italic_z start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT , italic_z start_POSTSUBSCRIPT italic_κ end_POSTSUBSCRIPT , italic_z start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT ) ≡ ⟨ divide start_ARG italic_κ roman_Δ italic_T start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_uncleaned end_POSTSUPERSCRIPT end_ARG start_ARG italic_κ roman_Δ italic_T start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_cleaned end_POSTSUPERSCRIPT end_ARG ⟩ start_POSTSUBSCRIPT 10 < roman_ℓ < 1500 end_POSTSUBSCRIPT . (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 κ𝜅\kappaitalic_κHI 2-point correlations drops by approximately the same factor Acleansubscript𝐴cleanA_{\mathrm{clean}}italic_A start_POSTSUBSCRIPT roman_clean end_POSTSUBSCRIPT across a wide redshift range, if we remove the background noise up to a particular redshift zmaxsubscript𝑧maxz_{\mathrm{max}}italic_z start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT, when cross-correlating to the κ𝜅\kappaitalic_κ field at a fixed redshift. Fig. 5 also implies that Aclean(zHI,zκ,zmax1)<Aclean(zHI,zκ,zmax2)subscript𝐴cleansubscript𝑧HIsubscript𝑧𝜅subscript𝑧max1subscript𝐴cleansubscript𝑧HIsubscript𝑧𝜅subscript𝑧max2A_{\mathrm{clean}}(z_{\mathrm{HI}},z_{\kappa},z_{\mathrm{max1}})<A_{\mathrm{% clean}}(z_{\mathrm{HI}},z_{\kappa},z_{\mathrm{max2}})italic_A start_POSTSUBSCRIPT roman_clean end_POSTSUBSCRIPT ( italic_z start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT , italic_z start_POSTSUBSCRIPT italic_κ end_POSTSUBSCRIPT , italic_z start_POSTSUBSCRIPT max1 end_POSTSUBSCRIPT ) < italic_A start_POSTSUBSCRIPT roman_clean end_POSTSUBSCRIPT ( italic_z start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT , italic_z start_POSTSUBSCRIPT italic_κ end_POSTSUBSCRIPT , italic_z start_POSTSUBSCRIPT max2 end_POSTSUBSCRIPT ) if the maximum redshifts zmax1subscript𝑧max1z_{\mathrm{max1}}italic_z start_POSTSUBSCRIPT max1 end_POSTSUBSCRIPT > zmax2subscript𝑧max2z_{\mathrm{max2}}italic_z start_POSTSUBSCRIPT max2 end_POSTSUBSCRIPT.

Refer to caption
Figure 4: The ratio between cleaned and uncleaned κΔTHI𝜅Δsubscript𝑇HI\kappa\Delta{T}_{\mathrm{HI}}italic_κ roman_Δ italic_T start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT power spectra. Two maximum redshifts (zmax)z_{\mathrm{max}})italic_z start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT ) for foreground removal are considered; zmax=1subscript𝑧max1z_{\mathrm{max}}=1italic_z start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT = 1 corresponds to current and imminent radio dishes meanwhile zmax=3subscript𝑧max3z_{\mathrm{max}}=3italic_z start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT = 3 represents a future SKA survey.
Refer to caption
Figure 5: The average ratio Acleansubscript𝐴cleanA_{\mathrm{clean}}italic_A start_POSTSUBSCRIPT roman_clean end_POSTSUBSCRIPT over >1010\ell>10roman_ℓ > 10 modes, as a function of redshift of HI slice used in the cross-correlation.

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 κ𝜅\kappaitalic_κ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 θisubscript𝜃𝑖\theta_{i}italic_θ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT are distributed by a multivariate Gaussian likelihood L𝐿Litalic_L, the Fisher matrix can be calculated as

𝑭ij<2θiθj>,subscript𝑭𝑖𝑗expectationsuperscript2subscript𝜃𝑖subscript𝜃𝑗\boldsymbol{F}_{ij}\equiv\bigg{<}\frac{\partial^{2}\mathcal{L}}{\partial\theta% _{i}\partial\theta_{j}}\bigg{>},bold_italic_F start_POSTSUBSCRIPT italic_i italic_j end_POSTSUBSCRIPT ≡ < divide start_ARG ∂ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT caligraphic_L end_ARG start_ARG ∂ italic_θ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ∂ italic_θ start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT end_ARG > , (19)

where =lnL𝐿\mathcal{L}=-\ln{L}caligraphic_L = - roman_ln italic_L. The Fisher matrix can be used to obtain the minimum uncertainty (σisubscript𝜎𝑖\sigma_{i}italic_σ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT) in parameter estimation due to the Cramér-Rao inequality (Mendez et al., 2014; Kamionkowski et al., 2011),

σi𝑭ii1,subscript𝜎𝑖subscriptsuperscript𝑭1𝑖𝑖\sigma_{i}\geqslant\sqrt{\boldsymbol{F}^{-1}_{ii}},italic_σ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ⩾ square-root start_ARG bold_italic_F start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_i italic_i end_POSTSUBSCRIPT end_ARG , (20)

which is equivalent to a 68%percent\%% confidence level. For a dataset where the uncertainties are Gaussian, the Fisher matrix can be calculated by (Tegmark et al., 1997)

𝑭ij=12Tr[𝑨i𝑨j+𝑪1𝑴ij]subscript𝑭𝑖𝑗12Trdelimited-[]subscript𝑨𝑖subscript𝑨𝑗superscript𝑪1subscript𝑴𝑖𝑗\boldsymbol{F}_{ij}=\frac{1}{2}\mathrm{Tr}[\boldsymbol{A}_{i}\boldsymbol{A}_{j% }+\boldsymbol{C}^{-1}\boldsymbol{M}_{ij}]bold_italic_F start_POSTSUBSCRIPT italic_i italic_j end_POSTSUBSCRIPT = divide start_ARG 1 end_ARG start_ARG 2 end_ARG roman_Tr [ bold_italic_A start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT bold_italic_A start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT + bold_italic_C start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT bold_italic_M start_POSTSUBSCRIPT italic_i italic_j end_POSTSUBSCRIPT ] (21)

where 𝑪𝑪\boldsymbol{C}bold_italic_C denotes the covariance matrix of the data, 𝑨isubscript𝑨𝑖\boldsymbol{A}_{i}bold_italic_A start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT \equiv 𝑪1𝑪,i\boldsymbol{C}^{-1}\boldsymbol{C}_{,i}bold_italic_C start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT bold_italic_C start_POSTSUBSCRIPT , italic_i end_POSTSUBSCRIPT, the derivative data matrix 𝑴ijsubscript𝑴𝑖𝑗\boldsymbol{M}_{ij}bold_italic_M start_POSTSUBSCRIPT italic_i italic_j end_POSTSUBSCRIPT \equiv 𝝁,i𝝁,jT\boldsymbol{\mu}_{,i}\boldsymbol{\mu}^{T}_{,j}bold_italic_μ start_POSTSUBSCRIPT , italic_i end_POSTSUBSCRIPT bold_italic_μ start_POSTSUPERSCRIPT italic_T end_POSTSUPERSCRIPT start_POSTSUBSCRIPT , italic_j end_POSTSUBSCRIPT + 𝝁,j𝝁,iT\boldsymbol{\mu}_{,j}\boldsymbol{\mu}^{T}_{,i}bold_italic_μ start_POSTSUBSCRIPT , italic_j end_POSTSUBSCRIPT bold_italic_μ start_POSTSUPERSCRIPT italic_T end_POSTSUPERSCRIPT start_POSTSUBSCRIPT , italic_i end_POSTSUBSCRIPT, and 𝝁𝝁\boldsymbol{\mu}bold_italic_μ is an expectation value of the data vector 𝒙𝒙\boldsymbol{x}bold_italic_x. The comma symbol means the partial derivative operator with respect to the parameter, 𝝁,i\boldsymbol{\mu}_{,i}bold_italic_μ start_POSTSUBSCRIPT , italic_i end_POSTSUBSCRIPT \equiv 𝝁/θi𝝁subscript𝜃𝑖\partial\boldsymbol{\mu}/\partial\theta_{i}∂ bold_italic_μ / ∂ italic_θ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT. Note that all derivatives are performed at the maximum likelihood point.

As we expect only small changes in the covariance matrix 𝐂𝐎𝐕(CXY)𝐂𝐎𝐕superscript𝐶𝑋𝑌\mathbf{COV}(C^{XY})bold_COV ( italic_C start_POSTSUPERSCRIPT italic_X italic_Y end_POSTSUPERSCRIPT ) 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

𝑭ij=XYCXYθiT𝐂𝐎𝐕(CXY)1CXYθj.subscript𝑭𝑖𝑗subscript𝑋𝑌superscriptsuperscript𝐶𝑋𝑌subscript𝜃𝑖𝑇𝐂𝐎𝐕superscriptsuperscript𝐶𝑋𝑌1superscript𝐶𝑋𝑌subscript𝜃𝑗\boldsymbol{F}_{ij}=\sum_{XY}\frac{\partial C^{XY}}{\partial\theta_{i}}^{T}% \mathbf{COV}(C^{XY})^{-1}\frac{\partial C^{XY}}{\partial\theta_{j}}.bold_italic_F start_POSTSUBSCRIPT italic_i italic_j end_POSTSUBSCRIPT = ∑ start_POSTSUBSCRIPT italic_X italic_Y end_POSTSUBSCRIPT divide start_ARG ∂ italic_C start_POSTSUPERSCRIPT italic_X italic_Y end_POSTSUPERSCRIPT end_ARG start_ARG ∂ italic_θ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT end_ARG start_POSTSUPERSCRIPT italic_T end_POSTSUPERSCRIPT bold_COV ( italic_C start_POSTSUPERSCRIPT italic_X italic_Y end_POSTSUPERSCRIPT ) start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT divide start_ARG ∂ italic_C start_POSTSUPERSCRIPT italic_X italic_Y end_POSTSUPERSCRIPT end_ARG start_ARG ∂ italic_θ start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT end_ARG . (22)

We calculate the cross-power spectra CXY()superscript𝐶𝑋𝑌C^{XY}(\ell)italic_C start_POSTSUPERSCRIPT italic_X italic_Y end_POSTSUPERSCRIPT ( roman_ℓ ) by using Eq. 8 with Planck 2018 cosmological parameters (Planck Collaboration et al., 2018). We calculate the covariance matrices of κκ𝜅𝜅\kappa\kappaitalic_κ italic_κ, κ𝜅\kappaitalic_κ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 𝐂𝐎𝐑𝐑ij=𝐂𝐎𝐕ij/𝐂𝐎𝐕ii𝐂𝐎𝐕jjsubscript𝐂𝐎𝐑𝐑𝑖𝑗subscript𝐂𝐎𝐕𝑖𝑗subscript𝐂𝐎𝐕𝑖𝑖subscript𝐂𝐎𝐕𝑗𝑗\mathbf{CORR}_{ij}=\mathbf{COV}_{ij}/\sqrt{\mathbf{COV}_{ii}\mathbf{COV}_{jj}}bold_CORR start_POSTSUBSCRIPT italic_i italic_j end_POSTSUBSCRIPT = bold_COV start_POSTSUBSCRIPT italic_i italic_j end_POSTSUBSCRIPT / square-root start_ARG bold_COV start_POSTSUBSCRIPT italic_i italic_i end_POSTSUBSCRIPT bold_COV start_POSTSUBSCRIPT italic_j italic_j end_POSTSUBSCRIPT end_ARG and show these in Fig. 6; these do not indicate significant correlations between \ellroman_ℓ bins.

Refer to caption
Refer to caption
Refer to caption
Figure 6: The correlation matrices of CXYsuperscriptsubscript𝐶𝑋𝑌C_{\ell}^{XY}italic_C start_POSTSUBSCRIPT roman_ℓ end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_X italic_Y end_POSTSUPERSCRIPT measured from 35 realisations for the cleaned HI signal at central redshift z=0.3𝑧0.3z=0.3italic_z = 0.3 and κ𝜅\kappaitalic_κ signal at z=0.78𝑧0.78z=0.78italic_z = 0.78. Left: κκ𝜅𝜅\kappa\kappaitalic_κ italic_κ, middle: ΔTHIΔTHIΔsubscript𝑇HIΔsubscript𝑇HI\mathrm{\Delta}T_{\mathrm{HI}}\mathrm{\Delta}T_{\mathrm{HI}}roman_Δ italic_T start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT roman_Δ italic_T start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT and right: κΔTHI𝜅Δsubscript𝑇HI\kappa\mathrm{\Delta}T_{\mathrm{HI}}italic_κ roman_Δ italic_T start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT.

4.2 Cosmological Constraints for Single-Slice Cross-correlations

In this section, the cosmological constraint viability of κ𝜅\kappaitalic_κ and HI cross-correlations is explored. We first start with the simplest observational configuration, considering only one redshift slice of HI and κ𝜅\kappaitalic_κ. We further assume that the bHI(z)subscript𝑏HI𝑧b_{\mathrm{HI}}(z)italic_b start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT ( italic_z ) 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 hhitalic_h = 0.67, ΩmsubscriptΩm\Omega_{\mathrm{m}}roman_Ω start_POSTSUBSCRIPT roman_m end_POSTSUBSCRIPT = 0.3, σ8subscript𝜎8\sigma_{8}italic_σ start_POSTSUBSCRIPT 8 end_POSTSUBSCRIPT = 0.82, ΩksubscriptΩk\Omega_{\mathrm{k}}roman_Ω start_POSTSUBSCRIPT roman_k end_POSTSUBSCRIPT = 0, ΩΛsubscriptΩΛ\Omega_{\Lambda}roman_Ω start_POSTSUBSCRIPT roman_Λ end_POSTSUBSCRIPT = 0.7, τ𝜏\tauitalic_τ = 0.06, and nssubscript𝑛𝑠n_{s}italic_n start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT = 0.96. To make a covariance matrix of cross-power spectra for the Fisher matrix (Eq. 22), we combine \ellroman_ℓ modes into 15 bins; each bin contains 101 \ellroman_ℓ modes with 11152711152711\leq\ell\geq 152711 ≤ roman_ℓ ≥ 1527 and averages over 35 realisations. We first consider the 3×\times×2 functions for a joint analysis of κ(0.78)κ(0.78)𝜅0.78𝜅0.78\kappa(0.78)\kappa(0.78)italic_κ ( 0.78 ) italic_κ ( 0.78 ), ΔTHI(0.3)ΔTHI(0.3)Δsubscript𝑇HI0.3Δsubscript𝑇HI0.3\mathrm{\Delta}T_{\mathrm{HI}}(0.3)\mathrm{\Delta}T_{\mathrm{HI}}(0.3)roman_Δ italic_T start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT ( 0.3 ) roman_Δ italic_T start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT ( 0.3 ) and κ(0.78)ΔTHI(0.3)𝜅0.78Δsubscript𝑇HI0.3\kappa(0.78)\mathrm{\Delta}T_{\mathrm{HI}}(0.3)italic_κ ( 0.78 ) roman_Δ italic_T start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT ( 0.3 ), 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 ‘2×\times×2’ 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 κκ𝜅𝜅\kappa-\kappaitalic_κ - italic_κ(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 3×2323\times 23 × 2 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 κκ𝜅𝜅\kappa\kappaitalic_κ italic_κ or ΔTHIΔsubscript𝑇HI\mathrm{\Delta}T_{\mathrm{HI}}roman_Δ italic_T start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPTΔTHIΔsubscript𝑇HI\mathrm{\Delta}T_{\mathrm{HI}}roman_Δ italic_T start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT constraints alone. Therefore in the next section, we will examine a joint likelihood between more redshift bins, and where the HI bias (bHI(z)subscript𝑏HI𝑧b_{\mathrm{HI}}(z)italic_b start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT ( italic_z )) is taken into account.

Refer to caption
Refer to caption
Figure 7: Left: Likelihood contours for the dataset described in section 4.2; contours show 68%percent\%% and 95%percent\%% confidence levels. Right: zoom-in of likelihood contours of 3×\times×2-point functions and their marginalisations from the left panel.

4.3 HI bias and multi-redshift bin joint likelihood analysis

It is well known that there is a degeneracy between galaxy bias, ΩmsubscriptΩm\Omega_{\mathrm{m}}roman_Ω start_POSTSUBSCRIPT roman_m end_POSTSUBSCRIPT and σ8subscript𝜎8\sigma_{8}italic_σ start_POSTSUBSCRIPT 8 end_POSTSUBSCRIPT 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 ΔTHIΔTHIΔsubscript𝑇HIΔsubscript𝑇HI\mathrm{\Delta}T_{\mathrm{HI}}\mathrm{\Delta}T_{\mathrm{HI}}roman_Δ italic_T start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT roman_Δ italic_T start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT measures bHI2(z)superscriptsubscript𝑏HI2𝑧b_{\mathrm{HI}}^{2}(z)italic_b start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ( italic_z ), κΔTHI𝜅Δsubscript𝑇HI\kappa\mathrm{\Delta}T_{\mathrm{HI}}italic_κ roman_Δ italic_T start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT additionally measures bHI(z)subscript𝑏HI𝑧b_{\mathrm{HI}}(z)italic_b start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT ( italic_z ). Combining the cross-bin intensity mapping power spectra with the κΔTHI𝜅Δsubscript𝑇HI\kappa\mathrm{\Delta}T_{\mathrm{HI}}italic_κ roman_Δ italic_T start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT 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 α𝛼\alphaitalic_α parameter in Eq. 17, setting the rest of the parameters to the best fit values (Cunnington et al., 2019b). In the second model, b0subscript𝑏0b_{0}italic_b start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT, b1subscript𝑏1b_{1}italic_b start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT and b2subscript𝑏2b_{2}italic_b start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT are the bias parameters with α𝛼\alphaitalic_α set to equal 1. We include these parameters when evaluating the Fisher matrices (Eq. 22). Note that both bHIsubscript𝑏HIb_{\mathrm{HI}}italic_b start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT models are scale-invariant and depend only on z𝑧zitalic_z. We will consider both full-sky and 300 deg2 surveys to explore the viability of HIHI and κ𝜅\kappaitalic_κHI in cosmological constraints.

For the full-sky case, as we include more parameters for bHIsubscript𝑏HIb_{\mathrm{HI}}italic_b start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT, we also examine more redshift bins for both HI and κ𝜅\kappaitalic_κ to obtain the best possible results. We consider the redshift range for ΔTHIΔsubscript𝑇HI\mathrm{\Delta}T_{\mathrm{HI}}roman_Δ italic_T start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT 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 ΔTHIΔsubscript𝑇HI\mathrm{\Delta}T_{\mathrm{HI}}roman_Δ italic_T start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT and κ𝜅\kappaitalic_κ bins; the width of each bin is 150 h1superscript1h^{-1}italic_h start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPTMpc, ΔzΔ𝑧\Delta{z}roman_Δ italic_z \approx 0.05. Table 1 lists both HI and κ𝜅\kappaitalic_κ central redshifts. As we have 16 zHIsubscript𝑧HIz_{\mathrm{HI}}italic_z start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT , we shall refer to ’16-HIHI’ which corresponds to 16 pairs of HI auto-correlation functions; we cross-correlate HI intensity maps to κ𝜅\kappaitalic_κ fields at zκsubscript𝑧𝜅z_{\kappa}italic_z start_POSTSUBSCRIPT italic_κ end_POSTSUBSCRIPT = 0.44, 0.78 and 1.77 respectively. We refer to 16-HIHI+1-κ𝜅\kappaitalic_κ as the joint analysis for 16-HIHI and κ𝜅\kappaitalic_κHI 2-point statistics at which zκsubscript𝑧𝜅z_{\kappa}italic_z start_POSTSUBSCRIPT italic_κ end_POSTSUBSCRIPT = 0.44. We add further zκsubscript𝑧𝜅z_{\kappa}italic_z start_POSTSUBSCRIPT italic_κ end_POSTSUBSCRIPT bins and label joint data as 16-HIHI+2-κ𝜅\kappaitalic_κ and 16-HIHI+3-κ𝜅\kappaitalic_κ. We calculate both the futuristic case where HI can be measured with high max1000subscriptmax1000\ell_{\mathrm{max}}\geq 1000roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT ≥ 1000 and the current state of art where 100<max<400100subscriptmax400100<\ell_{\mathrm{max}}<400100 < roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT < 400.

We further calculate the figure of merit (FoM) for the Ωmσ8subscriptΩmsubscript𝜎8\Omega_{\mathrm{m}}-\sigma_{8}roman_Ω start_POSTSUBSCRIPT roman_m end_POSTSUBSCRIPT - italic_σ start_POSTSUBSCRIPT 8 end_POSTSUBSCRIPT constraint. The FoM is the inverse of the area of the Ωmσ8subscriptΩmsubscript𝜎8\Omega_{\mathrm{m}}-\sigma_{8}roman_Ω start_POSTSUBSCRIPT roman_m end_POSTSUBSCRIPT - italic_σ start_POSTSUBSCRIPT 8 end_POSTSUBSCRIPT contours; in this case we calculate the FoM at 95%percent\%% confidence level. Fig. 9 shows the FoM of the Ωmσ8subscriptΩmsubscript𝜎8\Omega_{\mathrm{m}}-\sigma_{8}roman_Ω start_POSTSUBSCRIPT roman_m end_POSTSUBSCRIPT - italic_σ start_POSTSUBSCRIPT 8 end_POSTSUBSCRIPT 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 z<1𝑧1z<1italic_z < 1 HI slices are available). The green, red and black dots in Fig. 9 show the FoM for max=1530subscriptmax1530\ell_{\mathrm{max}}=1530roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT = 1530 when we further add the cross-correlations between consecutive HI slices and the κ𝜅\kappaitalic_κ slices at z𝑧zitalic_z = 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.

Refer to caption
Figure 8: Figure of Merit for σ8Ωmsubscript𝜎8subscriptΩ𝑚\sigma_{8}-\Omega_{m}italic_σ start_POSTSUBSCRIPT 8 end_POSTSUBSCRIPT - roman_Ω start_POSTSUBSCRIPT italic_m end_POSTSUBSCRIPT constraints; the horizontal axis is the number of redshift bin pairs for cosmological constraints. We show cumulative FoM when including increasing numbers of HI auto-correlation redshift bins (green); then increasing numbers of cross-correlations with convergence bins (grey, red, blue). Here max=1530subscriptmax1530\ell_{\mathrm{max}}=1530roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT = 1530.
Refer to caption
Figure 9: Figure of Merit for σ8Ωmsubscript𝜎8subscriptΩ𝑚\sigma_{8}-\Omega_{m}italic_σ start_POSTSUBSCRIPT 8 end_POSTSUBSCRIPT - roman_Ω start_POSTSUBSCRIPT italic_m end_POSTSUBSCRIPT constraints; the horizontal axis is the number of redshift bin pairs for cosmological constraints. We show cumulative FoM when including increasing numbers of HI auto-correlation redshift bins (green); then increasing numbers of cross-correlations with convergence bins (grey, red, blue). Here max=375subscriptmax375\ell_{\mathrm{max}}=375roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT = 375.
Table 1: The central redshifts for intensity and lensing convergence maps in our multi-redshift bin analysis. For κΔTHI𝜅Δsubscript𝑇HI\kappa\mathrm{\Delta}T_{\mathrm{HI}}italic_κ roman_Δ italic_T start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT, we require that zHI<0.7zκsubscript𝑧HI0.7subscript𝑧𝜅z_{\mathrm{HI}}<0.7z_{\kappa}italic_z start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT < 0.7 italic_z start_POSTSUBSCRIPT italic_κ end_POSTSUBSCRIPT.
zHIsubscript𝑧HIz_{\mathrm{HI}}italic_z start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT zκsubscript𝑧𝜅z_{\kappa}italic_z start_POSTSUBSCRIPT italic_κ end_POSTSUBSCRIPT
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 κ𝜅\kappaitalic_κ. With the appropriate redshift bin size of HI (ΔzHIΔsubscript𝑧HI\Delta{z}_{\mathrm{HI}}roman_Δ italic_z start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT) for max=1530subscriptmax1530\ell_{\mathrm{max}}=1530roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT = 1530, 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 κ𝜅\kappaitalic_κ redshift slices, the constraints are improved significantly especially for ΩmsubscriptΩm\Omega_{\mathrm{m}}roman_Ω start_POSTSUBSCRIPT roman_m end_POSTSUBSCRIPT and α𝛼\alphaitalic_α. However, this also makes the contours more elliptical, as there are remaining degeneracies among parameters. The uncertainties on parameters are measured at 95%percent\%% confidence level and reported in Table 2.

We next consider the current state of the art case, where since the typical angular resolution 1similar-toabsentsuperscript1\sim 1^{\circ}∼ 1 start_POSTSUPERSCRIPT ∘ end_POSTSUPERSCRIPT, we set maxsubscriptmax\ell_{\mathrm{max}}roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT = 375 (we choose this particular value as it is convenient to consider the ΔΔ\Delta{\ell}roman_Δ roman_ℓ bins as 15 bins with ΔΔ\Delta{\ell}roman_Δ roman_ℓ = 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 max=1530subscriptmax1530\ell_{\mathrm{max}}=1530roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT = 1530 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 κ𝜅\kappaitalic_κ bins to 16-HIHI (green shade) in maxsubscriptmax\ell_{\mathrm{max}}roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT = 375. We are therefore seeing that by joining κ𝜅\kappaitalic_κHI 2-point statistics to HIHI auto-correlations, we can improve cosmological constraints significantly.

Refer to caption
Figure 10: The likelihood contours for our multi-bin analysis with HI bias model 1 with max=1530subscriptmax1530\ell_{\mathrm{max}}=1530roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT = 1530; the contours show 68%percent\%% and 95%percent\%% confidence levels. The cosmological parameter uncertainties at 95%percent\%% are measured and reported in Table 2.
Refer to caption
Figure 11: The likelihood contours for our multi-bin analysis with HI bias model 1 with max=375subscriptmax375\ell_{\mathrm{max}}=375roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT = 375; the contours show 68%percent\%% and 95%percent\%% confidence levels. Compared to Fig. 11, we can see the difference in 16-HIHI constraint. However, when combining the κHI𝜅HI\kappa\mathrm{HI}italic_κ roman_HI-correlation the constraints are similar to those in Fig. 11.
Table 2: Cosmological forecast from Fisher analysis for 16-HIHI + 3-κ𝜅\kappaitalic_κ correlations; the uncertainties on cosmological parameters and HI bias are quoted at 95%percent\%% confidence level. HI bias model 1 considers only the re-scaling parameter α𝛼\alphaitalic_α. In contrast the second model considers quadratic parameters; b2subscript𝑏2b_{2}italic_b start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT is poorly constrained, while all other parameters are able to be measured well (see Fig. 12).
Parameters HI bias model 1 HI bias model 2
Δh0Δsubscript0\Delta{h_{0}}roman_Δ italic_h start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ±plus-or-minus\pm± 0.02 ±plus-or-minus\pm± 0.02
ΔΩmΔsubscriptΩm\Delta{\Omega_{\mathrm{m}}}roman_Δ roman_Ω start_POSTSUBSCRIPT roman_m end_POSTSUBSCRIPT ±plus-or-minus\pm± 0.01 ±plus-or-minus\pm± 0.02
Δσ8Δsubscript𝜎8\Delta{\sigma_{8}}roman_Δ italic_σ start_POSTSUBSCRIPT 8 end_POSTSUBSCRIPT ±plus-or-minus\pm± 0.03 ±plus-or-minus\pm± 0.04
ΔnsΔsubscript𝑛𝑠\Delta{n_{s}}roman_Δ italic_n start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT ±plus-or-minus\pm± 0.04 ±plus-or-minus\pm± 0.05
ΔαΔ𝛼\Delta\alpharoman_Δ italic_α ±plus-or-minus\pm± 0.04 -
Δb0Δsubscript𝑏0\Delta{b_{0}}roman_Δ italic_b start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT - ±plus-or-minus\pm± 0.04
Δb1Δsubscript𝑏1\Delta{b_{1}}roman_Δ italic_b start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT - ±plus-or-minus\pm± 0.03
Refer to caption
Figure 12: Constraints on cosmological parameters and bHI(z)subscript𝑏HI𝑧b_{\mathrm{HI}}(z)italic_b start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT ( italic_z ) for our second bias model, for 2-point functions at 68% and 95% levels of confidence (for our 16HIHI redshift slice case). The κ𝜅\kappaitalic_κHI correlation functions do not significantly improve the h0subscript0h_{0}italic_h start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT, nssubscript𝑛𝑠n_{s}italic_n start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT, b0subscript𝑏0b_{0}italic_b start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT and b1subscript𝑏1b_{1}italic_b start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT constraints. For this figure, we marginalised over b2subscript𝑏2b_{2}italic_b start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT as it is poorly constrained. However, we see significant improvement in ΩmsubscriptΩm\Omega_{\mathrm{m}}roman_Ω start_POSTSUBSCRIPT roman_m end_POSTSUBSCRIPT and σ8subscript𝜎8\sigma_{8}italic_σ start_POSTSUBSCRIPT 8 end_POSTSUBSCRIPT constraints. The parameter uncertainties at 95%percent\%% level are reported in Table 2.

For bias model 2, we find that the second order coefficient of HI bias b2subscript𝑏2b_{2}italic_b start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT is very poorly constrained. Marginalising over this parameter does not significantly affect the other cosmological parameter constraints (h0subscript0h_{0}italic_h start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT, ΩmsubscriptΩm\Omega_{\mathrm{m}}roman_Ω start_POSTSUBSCRIPT roman_m end_POSTSUBSCRIPT, σ8subscript𝜎8\sigma_{8}italic_σ start_POSTSUBSCRIPT 8 end_POSTSUBSCRIPT and nssubscript𝑛𝑠n_{s}italic_n start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT). 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 κ𝜅\kappaitalic_κ slices, the h0subscript0h_{0}italic_h start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT constraints do not improve much but the improvement in the Ωmσ8subscriptΩmsubscript𝜎8\Omega_{\mathrm{m}}-\sigma_{8}roman_Ω start_POSTSUBSCRIPT roman_m end_POSTSUBSCRIPT - italic_σ start_POSTSUBSCRIPT 8 end_POSTSUBSCRIPT 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 bHIsubscript𝑏HIb_{\mathrm{HI}}italic_b start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT 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 κ𝜅\kappaitalic_κ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 (C~subscript~𝐶\tilde{C}_{\ell}over~ start_ARG italic_C end_ARG start_POSTSUBSCRIPT roman_ℓ end_POSTSUBSCRIPT) as an estimator of Csubscript𝐶C_{\ell}italic_C start_POSTSUBSCRIPT roman_ℓ end_POSTSUBSCRIPT.

Suppose the survey footprint of the observations can be expressed using the weight function W(n^)𝑊^𝑛W(\hat{n})italic_W ( over^ start_ARG italic_n end_ARG ). Normalised by the sky factor fskysubscript𝑓skyf_{\mathrm{sky}}italic_f start_POSTSUBSCRIPT roman_sky end_POSTSUBSCRIPT (the fraction of the sky covered by the data), the weight moments are given by:

fskywi=14π𝑑n^Wi(n^),subscript𝑓skysubscript𝑤𝑖14𝜋differential-d^𝑛superscript𝑊𝑖^𝑛f_{\mathrm{sky}}w_{i}=\frac{1}{4\pi}\int d\hat{n}W^{i}(\hat{n}),italic_f start_POSTSUBSCRIPT roman_sky end_POSTSUBSCRIPT italic_w start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT = divide start_ARG 1 end_ARG start_ARG 4 italic_π end_ARG ∫ italic_d over^ start_ARG italic_n end_ARG italic_W start_POSTSUPERSCRIPT italic_i end_POSTSUPERSCRIPT ( over^ start_ARG italic_n end_ARG ) , (23)

where wisubscript𝑤𝑖w_{i}italic_w start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT represents the i𝑖iitalic_i-th moment of weighting. The power spectrum of the window function is:

𝒲=12+1m|wm|2.subscript𝒲121subscript𝑚superscriptsubscript𝑤𝑚2\mathcal{W}_{\ell}=\frac{1}{2\ell+1}\sum_{m}|w_{\ell m}|^{2}.caligraphic_W start_POSTSUBSCRIPT roman_ℓ end_POSTSUBSCRIPT = divide start_ARG 1 end_ARG start_ARG 2 roman_ℓ + 1 end_ARG ∑ start_POSTSUBSCRIPT italic_m end_POSTSUBSCRIPT | italic_w start_POSTSUBSCRIPT roman_ℓ italic_m end_POSTSUBSCRIPT | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT . (24)

For a spin-0 field ζ(n^)𝜁^𝑛\zeta(\hat{n})italic_ζ ( over^ start_ARG italic_n end_ARG ) weighted by W(n^)𝑊^𝑛W(\hat{n})italic_W ( over^ start_ARG italic_n end_ARG ), a spherical harmonic coefficient a~msubscript~𝑎𝑚\tilde{a}_{\ell m}over~ start_ARG italic_a end_ARG start_POSTSUBSCRIPT roman_ℓ italic_m end_POSTSUBSCRIPT can be expressed as (Kim & Naselsky, 2010; Kim, 2011):

a~m=𝑑n^ζ(n^)W(n^)Ym(n^)Ωppζ(p)W(p)Ym(p),subscript~𝑎𝑚differential-d^𝑛𝜁^𝑛𝑊^𝑛subscriptsuperscript𝑌𝑚^𝑛subscriptΩ𝑝subscript𝑝𝜁𝑝𝑊𝑝subscriptsuperscript𝑌𝑚𝑝\tilde{a}_{\ell m}=\int d\hat{n}\zeta(\hat{n})W(\hat{n})Y^{*}_{\ell m}(\hat{n}% )\approx\Omega_{p}\sum_{p}\zeta(p)W(p)Y^{*}_{\ell m}(p),over~ start_ARG italic_a end_ARG start_POSTSUBSCRIPT roman_ℓ italic_m end_POSTSUBSCRIPT = ∫ italic_d over^ start_ARG italic_n end_ARG italic_ζ ( over^ start_ARG italic_n end_ARG ) italic_W ( over^ start_ARG italic_n end_ARG ) italic_Y start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_ℓ italic_m end_POSTSUBSCRIPT ( over^ start_ARG italic_n end_ARG ) ≈ roman_Ω start_POSTSUBSCRIPT italic_p end_POSTSUBSCRIPT ∑ start_POSTSUBSCRIPT italic_p end_POSTSUBSCRIPT italic_ζ ( italic_p ) italic_W ( italic_p ) italic_Y start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_ℓ italic_m end_POSTSUBSCRIPT ( italic_p ) , (25)

where we approximate the integration over sky factor by the summation over pixel area with the surface density ΩpsubscriptΩ𝑝\Omega_{p}roman_Ω start_POSTSUBSCRIPT italic_p end_POSTSUBSCRIPT. The pseudo power spectrum estimator, C~subscript~𝐶\tilde{C}_{\ell}over~ start_ARG italic_C end_ARG start_POSTSUBSCRIPT roman_ℓ end_POSTSUBSCRIPT, is then:

C~=12+1m|a~m|2.subscript~𝐶121subscript𝑚superscriptsubscript~𝑎𝑚2\tilde{C}_{\ell}=\frac{1}{2\ell+1}\sum_{m}|\tilde{a}_{\ell m}|^{2}.over~ start_ARG italic_C end_ARG start_POSTSUBSCRIPT roman_ℓ end_POSTSUBSCRIPT = divide start_ARG 1 end_ARG start_ARG 2 roman_ℓ + 1 end_ARG ∑ start_POSTSUBSCRIPT italic_m end_POSTSUBSCRIPT | over~ start_ARG italic_a end_ARG start_POSTSUBSCRIPT roman_ℓ italic_m end_POSTSUBSCRIPT | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT . (26)

Similarly for spin-2 fields (γ1,γ2subscript𝛾1subscript𝛾2\gamma_{1},\gamma_{2}italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , italic_γ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT), we can obtain the coefficients, a~±2,msubscript~𝑎plus-or-minus2𝑚\tilde{a}_{\pm 2,\ell m}over~ start_ARG italic_a end_ARG start_POSTSUBSCRIPT ± 2 , roman_ℓ italic_m end_POSTSUBSCRIPT, by:

a~±2,m=[γ~1(n^)±iγ~2(n^)]Ym±2(n^)𝑑n^,subscript~𝑎plus-or-minus2𝑚delimited-[]plus-or-minussubscript~𝛾1^𝑛𝑖subscript~𝛾2^𝑛subscriptsubscriptsuperscript𝑌𝑚plus-or-minus2^𝑛differential-d^𝑛\tilde{a}_{\pm 2,\ell m}=\int[\tilde{\gamma}_{1}(\hat{n})\pm i\tilde{\gamma}_{% 2}(\hat{n})]\,\,{{}_{\pm 2}}Y^{*}_{\ell m}(\hat{n})d\hat{n},over~ start_ARG italic_a end_ARG start_POSTSUBSCRIPT ± 2 , roman_ℓ italic_m end_POSTSUBSCRIPT = ∫ [ over~ start_ARG italic_γ end_ARG start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ( over^ start_ARG italic_n end_ARG ) ± italic_i over~ start_ARG italic_γ end_ARG start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ( over^ start_ARG italic_n end_ARG ) ] start_FLOATSUBSCRIPT ± 2 end_FLOATSUBSCRIPT italic_Y start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_ℓ italic_m end_POSTSUBSCRIPT ( over^ start_ARG italic_n end_ARG ) italic_d over^ start_ARG italic_n end_ARG , (27)

where

γ~1(n^)±iγ~2=W(n^)[γ1(n^)±iγ2(n^)].plus-or-minussubscript~𝛾1^𝑛𝑖subscript~𝛾2𝑊^𝑛delimited-[]plus-or-minussubscript𝛾1^𝑛𝑖subscript𝛾2^𝑛\tilde{\gamma}_{1}(\hat{n})\pm i\tilde{\gamma}_{2}=W(\hat{n})[\gamma_{1}(\hat{% n})\pm i\gamma_{2}(\hat{n})].over~ start_ARG italic_γ end_ARG start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ( over^ start_ARG italic_n end_ARG ) ± italic_i over~ start_ARG italic_γ end_ARG start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT = italic_W ( over^ start_ARG italic_n end_ARG ) [ italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ( over^ start_ARG italic_n end_ARG ) ± italic_i italic_γ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ( over^ start_ARG italic_n end_ARG ) ] . (28)

Similarly to the full-sky case, a~mEsubscriptsuperscript~𝑎E𝑚\tilde{a}^{\mathrm{E}}_{\ell m}over~ start_ARG italic_a end_ARG start_POSTSUPERSCRIPT roman_E end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_ℓ italic_m end_POSTSUBSCRIPT and a~mBsubscriptsuperscript~𝑎B𝑚\tilde{a}^{\mathrm{B}}_{\ell m}over~ start_ARG italic_a end_ARG start_POSTSUPERSCRIPT roman_B end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_ℓ italic_m end_POSTSUBSCRIPT are then:

a~E=(a~2,m+a~2,m)/2,superscript~𝑎Esubscript~𝑎2𝑚subscript~𝑎2𝑚2\tilde{a}^{\mathrm{E}}=-(\tilde{a}_{2,\ell m}+\tilde{a}_{-2,\ell m})/2,over~ start_ARG italic_a end_ARG start_POSTSUPERSCRIPT roman_E end_POSTSUPERSCRIPT = - ( over~ start_ARG italic_a end_ARG start_POSTSUBSCRIPT 2 , roman_ℓ italic_m end_POSTSUBSCRIPT + over~ start_ARG italic_a end_ARG start_POSTSUBSCRIPT - 2 , roman_ℓ italic_m end_POSTSUBSCRIPT ) / 2 , (29)
a~B=i(a~2,ma~2,m)/2.superscript~𝑎B𝑖subscript~𝑎2𝑚subscript~𝑎2𝑚2\tilde{a}^{\mathrm{B}}=i(\tilde{a}_{2,\ell m}-\tilde{a}_{-2,\ell m})/2.over~ start_ARG italic_a end_ARG start_POSTSUPERSCRIPT roman_B end_POSTSUPERSCRIPT = italic_i ( over~ start_ARG italic_a end_ARG start_POSTSUBSCRIPT 2 , roman_ℓ italic_m end_POSTSUBSCRIPT - over~ start_ARG italic_a end_ARG start_POSTSUBSCRIPT - 2 , roman_ℓ italic_m end_POSTSUBSCRIPT ) / 2 . (30)

The pseudo power spectra for E and B modes are:

C~E,B=12+1m|a~mE,B|2.subscriptsuperscript~𝐶EB121subscript𝑚superscriptsubscriptsuperscript~𝑎EB𝑚2\tilde{C}^{\mathrm{E,B}}_{\ell}=\frac{1}{2\ell+1}\sum_{m}|\tilde{a}^{\mathrm{E% ,B}}_{\ell m}|^{2}.over~ start_ARG italic_C end_ARG start_POSTSUPERSCRIPT roman_E , roman_B end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_ℓ end_POSTSUBSCRIPT = divide start_ARG 1 end_ARG start_ARG 2 roman_ℓ + 1 end_ARG ∑ start_POSTSUBSCRIPT italic_m end_POSTSUBSCRIPT | over~ start_ARG italic_a end_ARG start_POSTSUPERSCRIPT roman_E , roman_B end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_ℓ italic_m end_POSTSUBSCRIPT | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT . (31)

The pseudo power spectra C~subscript~𝐶\tilde{C}_{\ell}over~ start_ARG italic_C end_ARG start_POSTSUBSCRIPT roman_ℓ end_POSTSUBSCRIPT and true Csubscript𝐶C_{\ell}italic_C start_POSTSUBSCRIPT roman_ℓ end_POSTSUBSCRIPT are related by the mode-mode coupling resulting from masking (Msubscript𝑀superscriptM_{\ell\ell^{\prime}}italic_M start_POSTSUBSCRIPT roman_ℓ roman_ℓ start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT end_POSTSUBSCRIPT):

C~=MC.delimited-⟨⟩subscript~𝐶subscriptsubscript𝑀superscriptdelimited-⟨⟩subscript𝐶\langle\tilde{C}_{\ell}\rangle=\sum_{\ell}M_{\ell\ell^{\prime}}\langle C_{\ell% }\rangle.⟨ over~ start_ARG italic_C end_ARG start_POSTSUBSCRIPT roman_ℓ end_POSTSUBSCRIPT ⟩ = ∑ start_POSTSUBSCRIPT roman_ℓ end_POSTSUBSCRIPT italic_M start_POSTSUBSCRIPT roman_ℓ roman_ℓ start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT end_POSTSUBSCRIPT ⟨ italic_C start_POSTSUBSCRIPT roman_ℓ end_POSTSUBSCRIPT ⟩ . (32)

This kernel depends solely on the geometry of a cut-sky 𝒲subscript𝒲\mathcal{W}_{\ell}caligraphic_W start_POSTSUBSCRIPT roman_ℓ end_POSTSUBSCRIPT and plays a crucial role in the pseudo-Csubscript𝐶C_{\ell}italic_C start_POSTSUBSCRIPT roman_ℓ end_POSTSUBSCRIPT 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-Csubscript𝐶C_{\ell}italic_C start_POSTSUBSCRIPT roman_ℓ end_POSTSUBSCRIPT for any spin fields, to evaluate HIHI and κ𝜅\kappaitalic_κHI pseudo-Csubscript𝐶C_{\ell}italic_C start_POSTSUBSCRIPT roman_ℓ end_POSTSUBSCRIPT 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 \ellroman_ℓ bins as before up to <375375\ell<375roman_ℓ < 375 with θbeam=1subscript𝜃beam1\theta_{\mathrm{beam}}=1italic_θ start_POSTSUBSCRIPT roman_beam end_POSTSUBSCRIPT = 1 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).

Refer to caption
Figure 13: Constraints on cosmological parameters and bHI(z)subscript𝑏HI𝑧b_{\mathrm{HI}}(z)italic_b start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT ( italic_z ) for 2-point functions at 68% and 95% levels of confidence, for the small area 300 deg2 case. This plot indicates that the feasibility of κ𝜅\kappaitalic_κHI (pseudo) 2-point statistics in cosmological constraints is heavily affected by the statistical incompleteness due to having a small-sky survey.

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 (P𝑃Pitalic_P). 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 (σpsubscript𝜎𝑝\sigma_{p}italic_σ start_POSTSUBSCRIPT italic_p end_POSTSUBSCRIPT) on the power spectrum P𝑃Pitalic_P can be estimated by evaluating its expected second moment. By calculating the ratio between σpsubscript𝜎𝑝\sigma_{p}italic_σ start_POSTSUBSCRIPT italic_p end_POSTSUBSCRIPT and P𝑃Pitalic_P averaging over the radial wavenumber bin size ΔkΔ𝑘\Delta{k}roman_Δ italic_k, we can estimate the uncertainty in the IM map measurements. Following this procedure the error on P𝑃Pitalic_P can be estimated by the following expression (Seo et al., 2010; Feldman et al., 1994; Battye et al., 2013):

σpP=2(2π3)Vsur14πk2Δk¯(1+σpix2Vpix[T¯(z)]2W(k)2P),subscript𝜎𝑝𝑃22𝜋3subscript𝑉sur14𝜋superscript𝑘2¯Δ𝑘1subscriptsuperscript𝜎2pixsubscript𝑉pixsuperscriptdelimited-[]¯𝑇𝑧2𝑊superscript𝑘2𝑃\frac{\sigma_{p}}{P}=2\sqrt{\left(\frac{2\pi}{3}\right)V_{\text{sur}}\frac{1}{% 4\pi k^{2}\overline{\Delta k}}}\bigg{(}1+\frac{\sigma^{2}_{\text{pix}}V_{\text% {pix}}}{\left[\overline{T}(z)\right]^{2}W(k)^{2}P}\bigg{)},divide start_ARG italic_σ start_POSTSUBSCRIPT italic_p end_POSTSUBSCRIPT end_ARG start_ARG italic_P end_ARG = 2 square-root start_ARG ( divide start_ARG 2 italic_π end_ARG start_ARG 3 end_ARG ) italic_V start_POSTSUBSCRIPT sur end_POSTSUBSCRIPT divide start_ARG 1 end_ARG start_ARG 4 italic_π italic_k start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT over¯ start_ARG roman_Δ italic_k end_ARG end_ARG end_ARG ( 1 + divide start_ARG italic_σ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT start_POSTSUBSCRIPT pix end_POSTSUBSCRIPT italic_V start_POSTSUBSCRIPT pix end_POSTSUBSCRIPT end_ARG start_ARG [ over¯ start_ARG italic_T end_ARG ( italic_z ) ] start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_W ( italic_k ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_P end_ARG ) , (33)

where Vsursubscript𝑉surV_{\mathrm{sur}}italic_V start_POSTSUBSCRIPT roman_sur end_POSTSUBSCRIPT is the volume of survey, W(k)𝑊𝑘W(k)italic_W ( italic_k ) is windows function, T¯(z)¯𝑇𝑧\bar{T}(z)over¯ start_ARG italic_T end_ARG ( italic_z ) is the average temperature at given redshift z𝑧zitalic_z and Vpixsubscript𝑉pixV_{\mathrm{pix}}italic_V start_POSTSUBSCRIPT roman_pix end_POSTSUBSCRIPT is the pixel volume. In this work we shall ignore the contribution of shot noise and only consider the contribution from pixel (thermal) noise σpixsubscript𝜎pix\sigma_{\mathrm{pix}}italic_σ start_POSTSUBSCRIPT roman_pix end_POSTSUBSCRIPT. These parameters are discussed in detail in the following paragraphs.

It is essential to pick the proper size of ΔkΔ𝑘\Delta kroman_Δ italic_k 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 Δk/kA<1Δ𝑘subscript𝑘𝐴1\Delta k/k_{A}<1roman_Δ italic_k / italic_k start_POSTSUBSCRIPT italic_A end_POSTSUBSCRIPT < 1, where kAsubscript𝑘𝐴k_{A}italic_k start_POSTSUBSCRIPT italic_A end_POSTSUBSCRIPT is the wavenumber of acoustic scale. The volume of the survey Vsursubscript𝑉surV_{\mathrm{sur}}italic_V start_POSTSUBSCRIPT roman_sur end_POSTSUBSCRIPT can be computed by

Vsur=Ωsurzminzmax𝑑zdVdzdΩ,subscript𝑉sursubscriptΩsursubscriptsuperscriptsubscript𝑧maxsubscript𝑧mindifferential-d𝑧𝑑𝑉𝑑𝑧𝑑ΩV_{\mathrm{sur}}=\Omega_{\mathrm{sur}}\int^{z_{\mathrm{max}}}_{z_{\mathrm{min}% }}dz\frac{dV}{dzd\Omega},italic_V start_POSTSUBSCRIPT roman_sur end_POSTSUBSCRIPT = roman_Ω start_POSTSUBSCRIPT roman_sur end_POSTSUBSCRIPT ∫ start_POSTSUPERSCRIPT italic_z start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_z start_POSTSUBSCRIPT roman_min end_POSTSUBSCRIPT end_POSTSUBSCRIPT italic_d italic_z divide start_ARG italic_d italic_V end_ARG start_ARG italic_d italic_z italic_d roman_Ω end_ARG , (34)

where

dVdzdΩ=cr2(z)H0E(z),𝑑𝑉𝑑𝑧𝑑Ω𝑐superscript𝑟2𝑧subscript𝐻0𝐸𝑧\frac{dV}{dzd\Omega}=\frac{cr^{2}(z)}{H_{0}E(z)},divide start_ARG italic_d italic_V end_ARG start_ARG italic_d italic_z italic_d roman_Ω end_ARG = divide start_ARG italic_c italic_r start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ( italic_z ) end_ARG start_ARG italic_H start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT italic_E ( italic_z ) end_ARG , (35)

where we assume a flat universe. The window function W(k)𝑊𝑘W(k)italic_W ( italic_k ) is set by the instrument’s specification. As we map the δHIsubscript𝛿HI\delta_{\mathrm{HI}}italic_δ start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT from multiple redshift bins, we can ignore the contribution of radial directions in W(k)𝑊𝑘W(k)italic_W ( italic_k ). 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 W(k)𝑊𝑘W(k)italic_W ( italic_k ) by:

W(k)=exp[12k2r2(z)θFWHM8ln(2)],𝑊𝑘12superscript𝑘2superscript𝑟2𝑧subscript𝜃FWHM82W(k)=\exp\bigg{[}-\frac{1}{2}k^{2}r^{2}(z)\frac{\theta_{\mathrm{FWHM}}}{8\ln{(% 2)}}\bigg{]},italic_W ( italic_k ) = roman_exp [ - divide start_ARG 1 end_ARG start_ARG 2 end_ARG italic_k start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_r start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ( italic_z ) divide start_ARG italic_θ start_POSTSUBSCRIPT roman_FWHM end_POSTSUBSCRIPT end_ARG start_ARG 8 roman_ln ( 2 ) end_ARG ] , (36)

where θFWHMsubscript𝜃FWHM\theta_{\mathrm{FWHM}}italic_θ start_POSTSUBSCRIPT roman_FWHM end_POSTSUBSCRIPT is the full width half max of angular resolution. Realistically, θFWHMsubscript𝜃FWHM\theta_{\mathrm{FWHM}}italic_θ start_POSTSUBSCRIPT roman_FWHM end_POSTSUBSCRIPT depends on frequency (ν𝜈\nuitalic_ν). However, in this analysis we assume the variation of θFWHMsubscript𝜃FWHM\theta_{\mathrm{FWHM}}italic_θ start_POSTSUBSCRIPT roman_FWHM end_POSTSUBSCRIPT is small and negligible.

The pixel volume, Vpixsubscript𝑉pixV_{\mathrm{pix}}italic_V start_POSTSUBSCRIPT roman_pix end_POSTSUBSCRIPT, can be determined by

Vpix=ΩpixzcΔzzc+Δz𝑑zdVdzdΩ,subscript𝑉pixsubscriptΩpixsubscriptsuperscriptsubscript𝑧cΔ𝑧subscript𝑧cΔ𝑧differential-d𝑧𝑑𝑉𝑑𝑧𝑑ΩV_{\mathrm{pix}}=\Omega_{\mathrm{pix}}\int^{z_{\mathrm{c}}+\Delta z}_{z_{% \mathrm{c}}-\Delta{z}}dz\frac{dV}{dzd\Omega},italic_V start_POSTSUBSCRIPT roman_pix end_POSTSUBSCRIPT = roman_Ω start_POSTSUBSCRIPT roman_pix end_POSTSUBSCRIPT ∫ start_POSTSUPERSCRIPT italic_z start_POSTSUBSCRIPT roman_c end_POSTSUBSCRIPT + roman_Δ italic_z end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_z start_POSTSUBSCRIPT roman_c end_POSTSUBSCRIPT - roman_Δ italic_z end_POSTSUBSCRIPT italic_d italic_z divide start_ARG italic_d italic_V end_ARG start_ARG italic_d italic_z italic_d roman_Ω end_ARG , (37)

where zcsubscript𝑧cz_{\mathrm{c}}italic_z start_POSTSUBSCRIPT roman_c end_POSTSUBSCRIPT is the central redshift and ΔzΔ𝑧\Delta{z}roman_Δ italic_z corresponds to frequency width (Δν)Δ𝜈(\Delta{\nu})( roman_Δ italic_ν ) and ΩpixsubscriptΩpix\Omega_{\mathrm{pix}}roman_Ω start_POSTSUBSCRIPT roman_pix end_POSTSUBSCRIPT 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 σpixsubscript𝜎pix\sigma_{\mathrm{pix}}italic_σ start_POSTSUBSCRIPT roman_pix end_POSTSUBSCRIPT can be approximated by (Santos et al., 2015; Seo et al., 2010),

σpixTsysεtp2Δν,subscript𝜎pixsubscript𝑇sys𝜀subscript𝑡𝑝2Δ𝜈\sigma_{\mathrm{pix}}\approx\frac{T_{\mathrm{sys}}}{\varepsilon\sqrt{t_{p}2% \Delta{\nu}}},italic_σ start_POSTSUBSCRIPT roman_pix end_POSTSUBSCRIPT ≈ divide start_ARG italic_T start_POSTSUBSCRIPT roman_sys end_POSTSUBSCRIPT end_ARG start_ARG italic_ε square-root start_ARG italic_t start_POSTSUBSCRIPT italic_p end_POSTSUBSCRIPT 2 roman_Δ italic_ν end_ARG end_ARG , (38)

where tpsubscript𝑡𝑝t_{p}italic_t start_POSTSUBSCRIPT italic_p end_POSTSUBSCRIPT represents the observation time per pointing, ε𝜀\varepsilonitalic_ε (approximately 1) signifies the efficiency of the telescope, meaning that almost no signal is lost when radio radiation is transmitted to the antenna. Tsyssubscript𝑇sysT_{\mathrm{sys}}italic_T start_POSTSUBSCRIPT roman_sys end_POSTSUBSCRIPT represents the system temperature, which includes

Tsys=Trx+Tspl+TCMB+Tgal,subscript𝑇syssubscript𝑇rxsubscript𝑇splsubscript𝑇CMBsubscript𝑇galT_{\mathrm{sys}}=T_{\mathrm{rx}}+T_{\mathrm{spl}}+T_{\mathrm{CMB}}+T_{\mathrm{% gal}},italic_T start_POSTSUBSCRIPT roman_sys end_POSTSUBSCRIPT = italic_T start_POSTSUBSCRIPT roman_rx end_POSTSUBSCRIPT + italic_T start_POSTSUBSCRIPT roman_spl end_POSTSUBSCRIPT + italic_T start_POSTSUBSCRIPT roman_CMB end_POSTSUBSCRIPT + italic_T start_POSTSUBSCRIPT roman_gal end_POSTSUBSCRIPT , (39)

where we ignore the contribution from the Earth’s atmosphere. Tsplsubscript𝑇splT_{\mathrm{spl}}italic_T start_POSTSUBSCRIPT roman_spl end_POSTSUBSCRIPT is the spill over from ground radiation (approximately 3K), TCMBsubscript𝑇CMBT_{\mathrm{CMB}}italic_T start_POSTSUBSCRIPT roman_CMB end_POSTSUBSCRIPT \approx 2.73K and galactic temperature Tgalsubscript𝑇galT_{\mathrm{gal}}italic_T start_POSTSUBSCRIPT roman_gal end_POSTSUBSCRIPT \approx 25K (408 MHz/ν𝜈\nuitalic_ν)2.75K  (Square Kilometre Array Cosmology Science Working Group et al., 2020). The observing time per pointing tpsubscript𝑡𝑝t_{p}italic_t start_POSTSUBSCRIPT italic_p end_POSTSUBSCRIPT relates to the total observation by tpsubscript𝑡𝑝t_{p}italic_t start_POSTSUBSCRIPT italic_p end_POSTSUBSCRIPT = ttot(θB)2/Ωsursubscript𝑡totsuperscriptsubscript𝜃B2subscriptΩsurt_{\mathrm{tot}}(\theta_{\mathrm{B}})^{2}/\Omega_{\mathrm{sur}}italic_t start_POSTSUBSCRIPT roman_tot end_POSTSUBSCRIPT ( italic_θ start_POSTSUBSCRIPT roman_B end_POSTSUBSCRIPT ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT / roman_Ω start_POSTSUBSCRIPT roman_sur end_POSTSUBSCRIPT, where θBsubscript𝜃B\theta_{\mathrm{B}}italic_θ start_POSTSUBSCRIPT roman_B end_POSTSUBSCRIPT is an angular pixel size.

As σpixsubscript𝜎pix\sigma_{\mathrm{pix}}italic_σ start_POSTSUBSCRIPT roman_pix end_POSTSUBSCRIPT is the rms thermal noise, by definition its square is the power per pixel volume (PN/Vpix)P_{N}/V_{\mathrm{pix}})italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT / italic_V start_POSTSUBSCRIPT roman_pix end_POSTSUBSCRIPT ). Therefore, the 3D noise power spectrum (PNsubscript𝑃𝑁P_{N}italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT) is then  (Battye et al., 2013; Santos et al., 2015)

PN=σpix2Vpix=r2yTsys2Ωsur2ε2ttot,subscript𝑃𝑁superscriptsubscript𝜎pix2subscript𝑉pixsuperscript𝑟2𝑦subscriptsuperscript𝑇2syssubscriptΩsur2superscript𝜀2subscript𝑡totP_{N}=\sigma_{\mathrm{pix}}^{2}V_{\mathrm{pix}}=r^{2}y\frac{T^{2}_{\mathrm{sys% }}\Omega_{\mathrm{sur}}}{2\varepsilon^{2}t_{\mathrm{tot}}},italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT = italic_σ start_POSTSUBSCRIPT roman_pix end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_V start_POSTSUBSCRIPT roman_pix end_POSTSUBSCRIPT = italic_r start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_y divide start_ARG italic_T start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_sys end_POSTSUBSCRIPT roman_Ω start_POSTSUBSCRIPT roman_sur end_POSTSUBSCRIPT end_ARG start_ARG 2 italic_ε start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_t start_POSTSUBSCRIPT roman_tot end_POSTSUBSCRIPT end_ARG , (40)

where

y=cH(z)1(1+z)2ν21.𝑦𝑐𝐻superscript𝑧1superscript1𝑧2subscript𝜈21y=cH(z)^{-1}\frac{(1+z)^{2}}{\nu_{21}}.italic_y = italic_c italic_H ( italic_z ) start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT divide start_ARG ( 1 + italic_z ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG italic_ν start_POSTSUBSCRIPT 21 end_POSTSUBSCRIPT end_ARG . (41)

If we have Ndsubscript𝑁dN_{\mathrm{d}}italic_N start_POSTSUBSCRIPT roman_d end_POSTSUBSCRIPT dishes where each dish has Nbsubscript𝑁bN_{\mathrm{b}}italic_N start_POSTSUBSCRIPT roman_b end_POSTSUBSCRIPT beams, we can take less time for each pointing area. The noise power spectrum is then reduced to

PN(Nd,Nb)=r2yTsys2Ωsur2ε2ttotNbNd.subscript𝑃𝑁subscript𝑁dsubscript𝑁bsuperscript𝑟2𝑦subscriptsuperscript𝑇2syssubscriptΩsur2superscript𝜀2subscript𝑡totsubscript𝑁bsubscript𝑁dP_{N}(N_{\mathrm{d}},N_{\mathrm{b}})=\frac{r^{2}yT^{2}_{\mathrm{sys}}\Omega_{% \mathrm{sur}}}{2\varepsilon^{2}t_{\mathrm{tot}}N_{\mathrm{b}}N_{\mathrm{d}}}.italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT ( italic_N start_POSTSUBSCRIPT roman_d end_POSTSUBSCRIPT , italic_N start_POSTSUBSCRIPT roman_b end_POSTSUBSCRIPT ) = divide start_ARG italic_r start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_y italic_T start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_sys end_POSTSUBSCRIPT roman_Ω start_POSTSUBSCRIPT roman_sur end_POSTSUBSCRIPT end_ARG start_ARG 2 italic_ε start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_t start_POSTSUBSCRIPT roman_tot end_POSTSUBSCRIPT italic_N start_POSTSUBSCRIPT roman_b end_POSTSUBSCRIPT italic_N start_POSTSUBSCRIPT roman_d end_POSTSUBSCRIPT end_ARG . (42)

To determine the optimisation of a survey strategy, we can estimate the suitable θFWHMsubscript𝜃FWHM\theta_{\mathrm{FWHM}}italic_θ start_POSTSUBSCRIPT roman_FWHM end_POSTSUBSCRIPT and ΩsursubscriptΩsur\Omega_{\mathrm{sur}}roman_Ω start_POSTSUBSCRIPT roman_sur end_POSTSUBSCRIPT that minimises δkA/kA𝛿subscript𝑘Asubscript𝑘A\delta{k}_{\mathrm{A}}/k_{\mathrm{A}}italic_δ italic_k start_POSTSUBSCRIPT roman_A end_POSTSUBSCRIPT / italic_k start_POSTSUBSCRIPT roman_A end_POSTSUBSCRIPT for acoustic scale kAsubscript𝑘Ak_{\mathrm{A}}italic_k start_POSTSUBSCRIPT roman_A end_POSTSUBSCRIPT. This acoustic scale kAsubscript𝑘Ak_{\mathrm{A}}italic_k start_POSTSUBSCRIPT roman_A end_POSTSUBSCRIPT can be estimated by following the work of  Blake & Glazebrook (2003); Battye et al. (2013),

P(k)Pref=1+Akexp[(k0.1hMpc1)1.4]sin(2πk/kA),𝑃𝑘subscript𝑃ref1𝐴𝑘superscript𝑘0.1superscriptMpc11.42𝜋𝑘subscript𝑘A\frac{P(k)}{P_{\mathrm{ref}}}=1+Ak\exp{\bigg{[}-\bigg{(}\frac{k}{0.1h\mathrm{% Mpc^{-1}}}\bigg{)}^{1.4}\bigg{]}}\sin{(2\pi k/k_{\mathrm{A}})},divide start_ARG italic_P ( italic_k ) end_ARG start_ARG italic_P start_POSTSUBSCRIPT roman_ref end_POSTSUBSCRIPT end_ARG = 1 + italic_A italic_k roman_exp [ - ( divide start_ARG italic_k end_ARG start_ARG 0.1 italic_h roman_Mpc start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT end_ARG ) start_POSTSUPERSCRIPT 1.4 end_POSTSUPERSCRIPT ] roman_sin ( 2 italic_π italic_k / italic_k start_POSTSUBSCRIPT roman_A end_POSTSUBSCRIPT ) , (43)

where A𝐴Aitalic_A 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 σpixMKsubscriptsuperscript𝜎MKpix\sigma^{\mathrm{MK}}_{\mathrm{pix}}italic_σ start_POSTSUPERSCRIPT roman_MK end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_pix end_POSTSUBSCRIPT = σTsubscript𝜎T\sigma_{\mathrm{T}}italic_σ start_POSTSUBSCRIPT roman_T end_POSTSUBSCRIPT. 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 similar-to\sim 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 σT2subscript𝜎T2\sigma_{\mathrm{T}}\approx 2italic_σ start_POSTSUBSCRIPT roman_T end_POSTSUBSCRIPT ≈ 2 mK.

Table 3: MeerKAT pilot survey specifications (Wang et al., 2021)
ΔνΔ𝜈\Delta\nuroman_Δ italic_ν 0.2 MHz
NΔνsubscript𝑁Δ𝜈N_{\Delta\nu}italic_N start_POSTSUBSCRIPT roman_Δ italic_ν end_POSTSUBSCRIPT [200,250]
Trxsubscript𝑇rxT_{\mathrm{rx}}italic_T start_POSTSUBSCRIPT roman_rx end_POSTSUBSCRIPT 7.5×103+103(ν[MHz]/10000.75)27.5superscript103superscript103superscript𝜈delimited-[]MHz10000.7527.5\times 10^{3}+10^{3}(\nu[\mathrm{MHz}]/1000-0.75)^{2}7.5 × 10 start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT + 10 start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT ( italic_ν [ roman_MHz ] / 1000 - 0.75 ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT [mK]
ttot 10.5 hours
z𝑧zitalic_z [0.3885, 0.4623]
Ndishsubscript𝑁dishN_{\mathrm{dish}}italic_N start_POSTSUBSCRIPT roman_dish end_POSTSUBSCRIPT 64
Tsyssubscript𝑇sysT_{\mathrm{sys}}italic_T start_POSTSUBSCRIPT roman_sys end_POSTSUBSCRIPT 16 ×103absentsuperscript103\times 10^{3}× 10 start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT mK
Npixsubscript𝑁pixN_{\mathrm{pix}}italic_N start_POSTSUBSCRIPT roman_pix end_POSTSUBSCRIPT 87500
θFHWMsubscript𝜃FHWM\theta_{\mathrm{FHWM}}italic_θ start_POSTSUBSCRIPT roman_FHWM end_POSTSUBSCRIPT 1.48 deg
ΩsursubscriptΩsur\Omega_{\mathrm{sur}}roman_Ω start_POSTSUBSCRIPT roman_sur end_POSTSUBSCRIPT 200 deg2

If we use Eq. (38) and (39) together with Table 3, the expected σpixsubscript𝜎pix\sigma_{\mathrm{pix}}italic_σ start_POSTSUBSCRIPT roman_pix end_POSTSUBSCRIPT for a single frequency channel of MeerKAT pilot survey is

σpix(Δν=0.2;10hr)15mK,subscript𝜎pixΔ𝜈0.210hr15mK\sigma_{\mathrm{pix}}(\Delta\nu=0.2;10\rm{hr})\approx 15\,\mathrm{mK},italic_σ start_POSTSUBSCRIPT roman_pix end_POSTSUBSCRIPT ( roman_Δ italic_ν = 0.2 ; 10 roman_h roman_r ) ≈ 15 roman_mK , (44)

where we assume each dish has equal efficiency ε𝜀\varepsilonitalic_ε =1 and consider only a single frequency channel ΔνΔ𝜈\Delta{\nu}roman_Δ italic_ν = 0.2 MHz. If we consider the whole frequency range like (Wang et al., 2021) the σTsubscript𝜎T\sigma_{\mathrm{T}}italic_σ start_POSTSUBSCRIPT roman_T end_POSTSUBSCRIPT \approx 2mK.

We now consider the case where the total observation time ttotsubscript𝑡tott_{\mathrm{tot}}italic_t start_POSTSUBSCRIPT roman_tot end_POSTSUBSCRIPT = 1000 hours and 250 frequency channels with ΔνΔ𝜈\Delta\nuroman_Δ italic_ν = 0.2 MHz.

5.2 S/N of κHI𝜅HI\kappa\mathrm{HI}italic_κ roman_HI

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 ttotsubscript𝑡tott_{\mathrm{tot}}italic_t start_POSTSUBSCRIPT roman_tot end_POSTSUBSCRIPT 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 (κ𝜅\kappaitalic_κHI) 2-point statistics. The estimate σpixTsubscriptsuperscript𝜎Tpix\sigma^{\mathrm{T}}_{\mathrm{pix}}italic_σ start_POSTSUPERSCRIPT roman_T end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_pix end_POSTSUBSCRIPT 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 σpixTsubscriptsuperscript𝜎Tpix\sigma^{\mathrm{T}}_{\mathrm{pix}}italic_σ start_POSTSUPERSCRIPT roman_T end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_pix end_POSTSUBSCRIPT as far as possible.

We first model the S/N of κ𝜅\kappaitalic_κHI. We can consider the zero lag noise level for κΔTdelimited-⟨⟩𝜅Δ𝑇\langle\kappa\Delta{T}\rangle⟨ italic_κ roman_Δ italic_T ⟩, i.e. where κ𝜅\kappaitalic_κ and ΔTΔ𝑇\Delta Troman_Δ italic_T are measured in the same pixel. There is no reason why the statistical noise of κ𝜅\kappaitalic_κ should be correlated with ΔTΔ𝑇\Delta{T}roman_Δ italic_T. The noise for the cross-correlation will therefore be proportional to the product of σTsubscript𝜎T\sigma_{\mathrm{T}}italic_σ start_POSTSUBSCRIPT roman_T end_POSTSUBSCRIPT and σensubscriptsuperscript𝜎𝑛𝑒\sigma^{n}_{e}italic_σ start_POSTSUPERSCRIPT italic_n end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT, where σTsubscript𝜎T\sigma_{\mathrm{T}}italic_σ start_POSTSUBSCRIPT roman_T end_POSTSUBSCRIPT is HI thermal noise which can be estimated by Eq. 44. The rms noise for weak lensing σensubscriptsuperscript𝜎𝑛𝑒\sigma^{n}_{e}italic_σ start_POSTSUPERSCRIPT italic_n end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT can be estimated by

σen=σengal,subscriptsuperscript𝜎𝑛𝑒subscript𝜎𝑒subscript𝑛gal\sigma^{n}_{e}=\frac{\sigma_{e}}{\sqrt{n_{\mathrm{gal}}}},italic_σ start_POSTSUPERSCRIPT italic_n end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT = divide start_ARG italic_σ start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT end_ARG start_ARG square-root start_ARG italic_n start_POSTSUBSCRIPT roman_gal end_POSTSUBSCRIPT end_ARG end_ARG , (45)

where σesubscript𝜎𝑒\sigma_{e}italic_σ start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT is the variance of intrinsic galaxy ellipticities and ngalsubscript𝑛galn_{\mathrm{gal}}italic_n start_POSTSUBSCRIPT roman_gal end_POSTSUBSCRIPT is galaxy number per pixel. For KiDS and DES-like surveys, σesubscript𝜎𝑒\sigma_{e}italic_σ start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT \approx 0.3. The KiDS DR4 effective galaxy number density is neffsubscript𝑛effn_{\mathrm{eff}}italic_n start_POSTSUBSCRIPT roman_eff end_POSTSUBSCRIPT = 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 σensubscriptsuperscript𝜎𝑛𝑒\sigma^{n}_{e}italic_σ start_POSTSUPERSCRIPT italic_n end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT \approx 0.03.

The S/N of κHIdelimited-⟨⟩𝜅HI\langle\kappa\mathrm{HI}\rangle⟨ italic_κ roman_HI ⟩ then can be estimated by

S/N=σHIσκσTσenNpix.𝑆𝑁subscript𝜎HIsubscript𝜎𝜅subscript𝜎Tsubscriptsuperscript𝜎𝑛𝑒subscript𝑁pixS/N=\frac{\sigma_{\mathrm{HI}}\sigma_{\kappa}}{\sigma_{\mathrm{T}}\sigma^{n}_{% e}}\sqrt{N_{\mathrm{pix}}}.italic_S / italic_N = divide start_ARG italic_σ start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT italic_σ start_POSTSUBSCRIPT italic_κ end_POSTSUBSCRIPT end_ARG start_ARG italic_σ start_POSTSUBSCRIPT roman_T end_POSTSUBSCRIPT italic_σ start_POSTSUPERSCRIPT italic_n end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT end_ARG square-root start_ARG italic_N start_POSTSUBSCRIPT roman_pix end_POSTSUBSCRIPT end_ARG . (46)

and we find that the rms lensing convergence signal σκsubscript𝜎𝜅\sigma_{\kappa}italic_σ start_POSTSUBSCRIPT italic_κ end_POSTSUBSCRIPT is similar to the rms noise σensubscriptsuperscript𝜎𝑛𝑒\sigma^{n}_{e}italic_σ start_POSTSUPERSCRIPT italic_n end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT 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 δTHI𝛿subscript𝑇HI\delta{T}_{\mathrm{HI}}italic_δ italic_T start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT can be expressed by:

δTHI=T¯HI(z)bHI(z)δm(θ,z),𝛿subscript𝑇HIsubscript¯𝑇HI𝑧subscript𝑏HI𝑧subscript𝛿m𝜃𝑧\delta{T}_{\mathrm{HI}}=\bar{T}_{\mathrm{HI}}(z)b_{\mathrm{HI}}(z)\delta_{% \mathrm{m}}(\theta,z),italic_δ italic_T start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT = over¯ start_ARG italic_T end_ARG start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT ( italic_z ) italic_b start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT ( italic_z ) italic_δ start_POSTSUBSCRIPT roman_m end_POSTSUBSCRIPT ( italic_θ , italic_z ) , (47)

where T¯HI(z)subscript¯𝑇HI𝑧\bar{T}_{\mathrm{HI}}(z)over¯ start_ARG italic_T end_ARG start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT ( italic_z ) is an average temperature over angular position (θ𝜃\thetaitalic_θ) for a given z𝑧zitalic_z and bHI(z)subscript𝑏HI𝑧b_{\mathrm{HI}}(z)italic_b start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT ( italic_z ) is HI bias. As the power spectrum PHIsubscript𝑃HIP_{\mathrm{HI}}italic_P start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT is the power of the temperature fluctuation δTHI𝛿subscript𝑇HI\delta T_{\mathrm{HI}}italic_δ italic_T start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT, the square root of PHIsubscript𝑃HIP_{\mathrm{HI}}italic_P start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT per volume Vsursubscript𝑉surV_{\mathrm{sur}}italic_V start_POSTSUBSCRIPT roman_sur end_POSTSUBSCRIPT is effectively the root mean square of the HI true signal (σHIsubscript𝜎HI\sigma_{\mathrm{HI}}italic_σ start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT),

σHI=PHI/Vsur.subscript𝜎HIsubscript𝑃HIsubscript𝑉sur\sigma_{\mathrm{HI}}=\sqrt{P_{\mathrm{HI}}/V_{\mathrm{sur}}}.italic_σ start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT = square-root start_ARG italic_P start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT / italic_V start_POSTSUBSCRIPT roman_sur end_POSTSUBSCRIPT end_ARG . (48)

We consider a survey similar to the MeerKAT pilot survey (0.39 < z < 0.46) over a moderately thick ΔzΔ𝑧\Delta{z}roman_Δ italic_z = 0.075 redshift bin, with 0.252 deg2 pixel size, Vsursubscript𝑉surV_{\mathrm{sur}}italic_V start_POSTSUBSCRIPT roman_sur end_POSTSUBSCRIPT similar-to\sim 4000 Mpc3h3superscript3h^{-3}italic_h start_POSTSUPERSCRIPT - 3 end_POSTSUPERSCRIPT. Assuming the foregrounds and redshift space distortions have been appropriately dealt with, we can use

PHI=T¯HI2bHI2Pmsubscript𝑃HIsubscriptsuperscript¯𝑇2HIsubscriptsuperscript𝑏2HIsubscript𝑃mP_{\mathrm{HI}}=\bar{T}^{2}_{\mathrm{HI}}b^{2}_{\mathrm{HI}}P_{\mathrm{m}}italic_P start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT = over¯ start_ARG italic_T end_ARG start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT italic_b start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT roman_m end_POSTSUBSCRIPT (49)

We assume an effective redshift of the survey zeffsubscript𝑧effz_{\mathrm{eff}}italic_z start_POSTSUBSCRIPT roman_eff end_POSTSUBSCRIPT = 0.42, bHIsubscript𝑏HIb_{\mathrm{HI}}italic_b start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT =1 and T¯HIsubscript¯𝑇HI\bar{T}_{\mathrm{HI}}over¯ start_ARG italic_T end_ARG start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT = 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

σHI(k=0.1,z=0.42)=5μK,subscript𝜎HIformulae-sequence𝑘0.1𝑧0.425𝜇K\sigma_{\mathrm{HI}}(k=0.1,z=0.42)=5\mu\mathrm{K},italic_σ start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT ( italic_k = 0.1 , italic_z = 0.42 ) = 5 italic_μ roman_K , (50)

where we estimate at the k𝑘kitalic_k = 0.1 h1superscript1h^{-1}italic_h start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPTMpc 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.

σTsubscript𝜎T\sigma_{\mathrm{T}}italic_σ start_POSTSUBSCRIPT roman_T end_POSTSUBSCRIPT is approximately 15 mK for a single frequency channel given Table 3. If we stack over 200 ΔνΔ𝜈\Delta\nuroman_Δ italic_ν channels and assume 10 hours observing time, then

σT1.1mK.subscript𝜎T1.1mK\sigma_{\mathrm{T}}\approx 1.1\,\mathrm{mK}.italic_σ start_POSTSUBSCRIPT roman_T end_POSTSUBSCRIPT ≈ 1.1 roman_mK . (51)

Then using equation (46), the estimation of S/N for κ𝜅\kappaitalic_κHI 2-point statistics for KiDS-like lensing surveys and MeerKAT for k𝑘kitalic_k = 0.1 h1superscript1h^{-1}italic_h start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT Mpc with a pixel size 0.252 deg2 covering similar-to\sim 200 deg2 sky, is then

S/N0.24.𝑆𝑁0.24S/N\approx 0.24.italic_S / italic_N ≈ 0.24 . (52)

This means that by the current state of art, we expected to observe more instrument noise than cosmic signal from κ𝜅\kappaitalic_κ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 σMKsubscript𝜎MK\sigma_{\mathrm{MK}}italic_σ start_POSTSUBSCRIPT roman_MK end_POSTSUBSCRIPT is then

σMK10000.01mK.subscriptsuperscript𝜎1000MK0.01mK\sigma^{1000}_{\mathrm{MK}}\approx 0.01\,\mathrm{mK}.italic_σ start_POSTSUPERSCRIPT 1000 end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_MK end_POSTSUBSCRIPT ≈ 0.01 roman_mK . (53)

Now we have better S/N by one order of magnitude for a 3000 pixel sky. (If we use the pilot survey footprint (similar-to\sim 200 deg2), then S/N \approx 2.4 for zero lag.)

We confirm this calculation by generating Gaussian random fields for the instrument noise for both κ𝜅\kappaitalic_κ 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 Npixsubscript𝑁pixN_{\mathrm{pix}}italic_N start_POSTSUBSCRIPT roman_pix end_POSTSUBSCRIPT = 196608. Hence the estimation of S/N for this configuration is then,

S/Nfull100022.𝑆subscriptsuperscript𝑁1000full22S/N^{1000}_{\mathrm{full}}\approx 22.italic_S / italic_N start_POSTSUPERSCRIPT 1000 end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_full end_POSTSUBSCRIPT ≈ 22 . (54)

However, when we take the LoS foreground subtraction into account, the signal would be reduced by a factor of similar-to\sim 3. This means

S/Nfull10007.S/N^{1000}_{\mathrm{full}}\longrightarrow\simeq 7.italic_S / italic_N start_POSTSUPERSCRIPT 1000 end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_full end_POSTSUBSCRIPT ⟶ ≃ 7 . (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 κ𝜅\kappaitalic_κHI and HIHI 2-point statistics for cosmological constraints in the presence of these current noise levels. In the S/N analysis of κ𝜅\kappaitalic_κHI (sec. 5.2), we only considered the case where Δνtotal=250×0.2Δsubscript𝜈total2500.2\Delta\nu_{\mathrm{total}}=250\times 0.2roman_Δ italic_ν start_POSTSUBSCRIPT roman_total end_POSTSUBSCRIPT = 250 × 0.2 MHz. However to compare results in sec. 4.3, we will adjust ΔνΔ𝜈\Delta{\nu}roman_Δ italic_ν to match zHIsubscript𝑧HIz_{\mathrm{HI}}italic_z start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT in Table. 1.

We generate Gaussian white noise fields for both κ𝜅\kappaitalic_κ 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 κ𝜅\kappaitalic_κHI are reduced significantly. We also show the resulting FoM in Fig. 15 and 2σ𝜎\sigmaitalic_σ-constraints in Table 4. We note that the ΩmsubscriptΩm\Omega_{\mathrm{m}}roman_Ω start_POSTSUBSCRIPT roman_m end_POSTSUBSCRIPT constraint shows a precise estimation, although it is degenerate with α𝛼\alphaitalic_α. 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. 99, and 15) indicate that by combining κ𝜅\kappaitalic_κHI likelihoods with HIHI, we can significantly enhance the feasibility of cosmological constraints from the HI intensity maps.

Refer to caption
Figure 14: The likelihood contours for our multi-bin analysis with HI bias model 1 with max=375subscriptmax375\ell_{\mathrm{max}}=375roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT = 375 for full-sky case including the instrument noise contributions; the contours show 68%percent\%% and 95%percent\%% confidence levels. Comparing to Fig. 11 where we ignore instrument noise contributions, we can see that h0subscript0h_{0}italic_h start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT is now poorly constrained.
Refer to caption
Figure 15: Figure of Merit for σ8Ωmsubscript𝜎8subscriptΩ𝑚\sigma_{8}-\Omega_{m}italic_σ start_POSTSUBSCRIPT 8 end_POSTSUBSCRIPT - roman_Ω start_POSTSUBSCRIPT italic_m end_POSTSUBSCRIPT constraints including instrument noise; the horizontal axis is the number of redshift bin pairs for cosmological constraints. We show cumulative FoM when including increasing numbers of HI auto-correlation redshift bins (green); then increasing numbers of cross-correlations with convergence bins (grey, red, blue). Here max=375subscriptmax375\ell_{\mathrm{max}}=375roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT = 375. Comparing this figure to Fig. 9, we find that the FoM is lower by 1 order of magnitude.
Table 4: Cosmological forecast from Fisher analysis including instrument noise contributions for 16-HIHI + 3-κ𝜅\kappaitalic_κ correlations; the uncertainties on cosmological parameters and HI bias are quoted at 95%percent\%% confident level. Here we consider only HI bias model 1.
Parameters HI bias model 1
Δh0Δsubscript0\Delta{h_{0}}roman_Δ italic_h start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ±plus-or-minus\pm± 0.28
ΔΩmΔsubscriptΩm\Delta{\Omega_{\mathrm{m}}}roman_Δ roman_Ω start_POSTSUBSCRIPT roman_m end_POSTSUBSCRIPT ±plus-or-minus\pm± 0.09
Δσ8Δsubscript𝜎8\Delta{\sigma_{8}}roman_Δ italic_σ start_POSTSUBSCRIPT 8 end_POSTSUBSCRIPT ±plus-or-minus\pm± 0.32
ΔnsΔsubscript𝑛𝑠\Delta{n_{s}}roman_Δ italic_n start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT ±plus-or-minus\pm± 0.17
ΔαΔ𝛼\Delta\alpharoman_Δ italic_α ±plus-or-minus\pm± 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 κ𝜅\kappaitalic_κHI power spectrum signal, by a factor Aclean(zHI,zκ,zmax)subscript𝐴cleansubscript𝑧HIsubscript𝑧𝜅subscript𝑧maxA_{\mathrm{clean}}(z_{\mathrm{HI}},z_{\kappa},z_{\mathrm{max}})italic_A start_POSTSUBSCRIPT roman_clean end_POSTSUBSCRIPT ( italic_z start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT , italic_z start_POSTSUBSCRIPT italic_κ end_POSTSUBSCRIPT , italic_z start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT ). In the case of forthcoming HI experiments that will measure HI at zmaxsubscript𝑧maxz_{\mathrm{max}}italic_z start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT < 1, Aclean(zHI<0.5,zκ=0.78,zmax=1)subscript𝐴cleanformulae-sequencesubscript𝑧HI0.5formulae-sequencesubscript𝑧𝜅0.78subscript𝑧max1A_{\mathrm{clean}}(z_{\mathrm{HI}}<0.5,z_{\kappa}=0.78,z_{\mathrm{max}}=1)italic_A start_POSTSUBSCRIPT roman_clean end_POSTSUBSCRIPT ( italic_z start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT < 0.5 , italic_z start_POSTSUBSCRIPT italic_κ end_POSTSUBSCRIPT = 0.78 , italic_z start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT = 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 κκ𝜅𝜅\kappa\kappaitalic_κ italic_κ, HIHI and κ𝜅\kappaitalic_κ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 κ𝜅\kappaitalic_κ can constrain cosmological parameters for known bias (see Fig. 7), but when bHIsubscript𝑏HIb_{\mathrm{HI}}italic_b start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT is a further parameter (or several), the few-slice 3×2323\times 23 × 2 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 bHIsubscript𝑏HIb_{\mathrm{HI}}italic_b start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT on cosmological constraints. We consider both the current state of art where max<400subscriptmax400\ell_{\mathrm{max}}<400roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT < 400, and the futuristic case where maxsubscriptmax\ell_{\mathrm{max}}roman_ℓ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT > 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 κ𝜅\kappaitalic_κ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 κ𝜅\kappaitalic_κ 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 δmsubscript𝛿m\delta_{\mathrm{m}}italic_δ start_POSTSUBSCRIPT roman_m end_POSTSUBSCRIPT. A more realistic study should construct HI intensity maps by assigning HI mass MHIsubscript𝑀HIM_{\mathrm{HI}}italic_M start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT 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 bHIsubscript𝑏HIb_{\mathrm{HI}}italic_b start_POSTSUBSCRIPT roman_HI end_POSTSUBSCRIPT 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 ΛΛ\Lambdaroman_Λ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, κ𝜅\kappaitalic_κ-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