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

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) and lensing convergence ($\kappa$). 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 $\sim$50-70\%; we present the angular and redshift dependence of the effect, which is a weak function of these variables. We then calculate S/N of $\kappa$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 $\kappa\Delta$$T_\mathrm{HI}$ correlation functions in a sufficient number of redshift bins, constraints on cosmology and HI bias will be possible


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 latetime 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).
★ E-mail: anut.sangka@port.ac.ukTo 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 (LRG), 5.7 (ELG), and 11.1 (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 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  g ) have focused on the impact of foreground removal (Chapman et al. 2012;Cunnington et al. 2019b;Padmanabhan et al. 2020;Cunnington et al. 2020;Spinelli et al. 2022).The study by Cunnington et al. (2019b) shows that the foreground removal affects 2-point function characteristics, especially when the redshift resolution is broad, as is the case in optical photometric surveys.
There are also numerous optical surveys measuring gravitational lensing shear () which distorts the shape of galaxy images; this is sensitive to density fluctuations of all the matter present along a line of sight, whether baryonic or dark matter.It is therefore of interest to consider the viability of the cross-correlation  −  HI , which will be able to be studied using a combination of lensing and IM surveys (Abbott et al. 2018;Baxter et al. 2019;Hu & Jain 2004;The CHIME Collaboration et al. 2022;CHIME Collaboration et al. 2022;Cunnington et al. 2022).The density projection along the unperturbed light ray trajectory, also known as 'lensing convergence'  can be considered instead of  as both share the same statistical properties.The 2-point functions between the pairs of  and HI could improve cosmological constraints and break degeneracies such as that between HI bias ( HI ) and clustering amplitude.
However, removing the HI foreground potentially affects these 2-point statistics, as the foreground removal effectively subtracts large-scale radial modes to which lensing is sensitive.In this paper we will calculate the cross-correlation function between convergence and 21cm intensity mapping, and will explore whether the foreground subtraction significantly hampers the cross-correlation measurement.We also explore whether the foreground removal impacts the viability of cosmological constraints from HI-HI and −HI correlations.
To achieve this, we will present theoretical and simulation approaches for calculating the -HI signal.We will then consider the effect of foreground removal on the signal, showing that the impact is significant (approximately a factor of 2 in signal reduction) but not lethal.We will then use the Fisher information matrix to make cosmological parameter forecasts for ideal and realistic surveys (including cut sky and the inclusion of telescope-specific noise), deploying the cross-correlation between convergence and intensity mapping, always including the effect of foreground removal.We discuss lensing convergence and HI simulation catalogues in Section 2, including modeling of the 2-point functions.We describe the HI foreground removal and its effect on -HI 2-point functions in Section 3. We present our Fisher forecasts for surveys in Section 4, effects of instrumental noise in Section 5 and present our conclusions in Section 6.

𝜅-HI 2-POINT STATISTICS: THEORY AND SIMULATIONS
In this section we discuss the relevant 2-point statistics.We shall start with theoretical calculations of 2-point functions of lensing convergence () and neutral hydrogen intensity maps (HI) in subsection 2.1.We will then discuss the generation of  catalogues and HI modelling from simulations of the matter overdensity .
The comparison between theoretical calculations and simulations is shown in subsection 2.2.The simulated HI maps will be used in the next section 3 where the foreground removal will be discussed.

Modeling the 2-point functions
In this subsection, we describe the modeling of the 2-point functions.We begin by considering how to calculate the observable quantities, namely weak lensing convergence  and HI temperature fluctuations Δ HI .We will then turn to the angular cross-power spectra.We denote  as the power spectra between  fields, HI  HI  as the cross-power spectra between HI fields, and HI as the crosspower between  and HI.The dummy indices  and  refer to the th and th redshift bins.We will calculate the lensing convergence in an arbitrary direction on the sky n using the Born approximation, projecting the matter overdensity  along an unperturbed ray direction.This can be computed by (Bartelmann & Schneider 2001) where  is comoving distance, Ω m is the matter density parameter at the present epoch,  0 is the Hubble parameter today and the subscript  refers to the source plane.For lensing of distributed sources in redshift bins , the integrand is modified by including a source distribution, so that the integration now becomes where the lensing weight is given by where    () is the lensing source number density, and n  is its average in the th redshift bin.
HI will be a biased tracer of matter overdensity, so we write Δ HI ( n, ) = THI () HI ()( n, ), where  HI () is the HI bias at a given redshift  and THI () is the average temperature.The projected temperature fluctuation at the th redshift bin is then where where   HI () is the HI source number density, and n HI is its average in the th redshift bin.Battye et al. (2013) show that for a given redshift , THI () can be estimated by where  () =  ()/ 0 is the dimensionless Hubble function at redshift .The HI density parameter could be approximated to be Ω HI ℎ = 2.45 × 10 −4 (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) Constraining the HI bias  HI () will be discussed later in section 4.
Using the Limber approximation, the angular power spectra   (ℓ) are given by where   ℓ+1/2  , ( ) is the matter power spectrum (LoVerde & Afshordi 2008).We compute the nonlinear power spectrum using the Boltzmann code CAMB (Lewis & Bridle 2002) with the Halofit extension to nonlinear scales (Takahashi et al. 2012).
Since the ray tracing simulations by Takahashi et al. (2017) which we use below adopt a comoving bin size Δ = 150 ℎ −1 Mpc (see section 2.2), we choose a radial selection function for   HI (( ))/ n HI as a normal distribution around a central comoving position with 3  = 150 ℎ −1 Mpc.This   corresponds to the frequency bandwidth (Δ) selected.In practice, the frequency range and bandwidth will depend on the particular radio telescope being used; for example, BINGO (Baryon Acoustic Oscillations in Neutral Gas Observations) has operational frequency from 960 MHz to 1260 MHz (Battye et al. 2013;Wuensche & the BINGO Collaboration 2019) and MeerKAT's frequency bandwidth ranges from 900-1185 MHz and 580-1000 MHz for L-band and UHF-band, respectively (Wang et al. 2021;Cunnington et al. 2022) .
We show examples for the first time of calculations of the autoand cross-power for HI and convergence in Fig. 1 using Eq. 8.As expected, the auto-signal depends on the radial HI width   , while the cross-power is insensitive to this.

Lensing Convergence and HI Intensity Maps
The full-sky gravitational lensing mock catalogues by Takahashi et al. (2017) have been used throughout this work.They are based on a multiple-lens ray-tracing approach through N-body cosmological simulations.The datasets include weak lensing maps (convergence, shear and rotation data) up to redshift 5.3, and halo catalogues.The catalogues provide 108 realisations of N-body simulations, 35 of which are used in this research (due to storage limitations).The N-body simulations were produced with periodic boundary conditions following dark matter gravitational evolution without baryonic processes.14 simulation boxes of side length  = 450, 900, 1350, ..., 6300 ℎ −1 Mpc are nested to represent a region of the Universe in which lensing occurs; each box contains 2048 3 particles.The  fields are obtained by tracing the light ray path through planes with separation 150 ℎ −1 Mpc.By calculating the Jacobian matrix  along the light path, the lensing convergence , shear lensing  1,2 and rotation angle  can be obtained, via 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 1 and DESI2 , the cross-correlation between the lensing convergence and the HI intensity map is limited by the lower angular resolution of HI intensity maps expected with real radio telescopes.Therefore the resolution is reduced to NSIDE of 512; this is not only appropriate for our 2-point function measurements but also decreases the storage space requirement and computational time.
We will first consider a convergence map at a specific optical lensing catalogue source redshift, which we choose as  ≈ 0.78, for which the lensing will significantly occur at the redshift of an intensity map at redshift≃ 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  HI from halo catalogues, we assume that HI is a biased tracer of the total matter overdensity field (, )(see Eq. 4 and 5), where  HI is a HI bias.For instance the parametric form for  HI adopted by Cunnington et al. (2019b) is Since the neutral hydrogen signal is measured as the surface brightness temperature, we shall refer to the HI intensity map as the temperature fluctuation Δ HI : We apply this equation to the overdensity map obtained from Takahashi et al. ( 2017) catalogues to create HI intensity maps.Fig. 2 shows the uncleaned and cleaned intensity maps from one realisation for the zoom-in patch with area 5 × 5 square degrees.
As we are interested in the 2-dimensional projection of cosmological fields on the sky, together with their power spectra, it is convenient to describe these fields Θ( n, ) in spherical harmonics: where   ℓ ( n) and  ℓ () are spherical harmonics and their coefficients respectively (Pratten et al. 2016;Castro et al. 2005;Heavens 2003).Θ( n, ) represents an arbitrary cosmological field; in this work it can be either lensing convergence or HI temperature fluctuations.The angular power spectrum is then an average of  ℓ over  modes: where  and  stand for the cosmological fields at given redshifts  1 and  2 , respectively.The cross power-spectrum for HI and lensing  can be easily measured via HEALPIX's anafast routine especially if the data is for the full sky (however, if the data has missing regions or a cut sky, pseudo- ℓ methods are required (Brown et al. 2005;Upham et al. 2019)).
Using this routine, we obtain cross-power measurements for the HI and  fields.We measure the cross-power spectra of 35 realisations and evaluate their mean; we show the results in Fig. 3.Here the lensing convergence is measured at the central redshift 0.78 and HI is measured at the central redshift 0.3.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   HI (( ))/ n HI successfully match the simulations.Due to the lens shell approximation of the ray tracing code, the measured  ℓ is slightly affected at very high ℓ (see red line on Fig. 3).

HI FOREGROUND REMOVAL AND ITS EFFECT ON 𝜅HI 2-POINT FUNCTIONS
The HI signal is small compared to its foregrounds such as freefree thermal emission, extragalactic radio sources and Galactic synchrotron.For example, the synchrotron ( sync ) emission temperature which can be modelled by  sysnc ∝ (1 + ) 2.7 [K] (Platania et al. 1998;Smoot & Debono 2017), is approximately 3 to 4 orders of magnitude larger than  HI 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 Δ LoS HI ( n).Since the calculation of lensing involves integration along the light path (Eq.3), which will have a contribution from longwavelength radial modes, the HI foreground removal is a concern for the existence of the HI cross-correlation (i.e.we have just removed such modes from the HI signal).In this section we therefore  seek to ascertain the degree to which the HI cross-correlation survives foreground removal.
Here we follow the method for foreground removal emulation by Cunnington et al. (2019b).The cleaned intensity map Δ clean HI can be approximated as where Δ orig HI ( n, ) is the uncleaned signal in direction n at redshift .Δ LoS HI ( n) is defined by: so that Δ LoS HI ( n) 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 Δ LoS HI ( n).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/  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/  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): where ,  0 ,  1 and  2 are set to 1, 0.67, 0.18 and 0.05, respectively.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 ℎ −1 Mpc, which are our main interest.
We measure Δ LOS HI with two choices of maximum redshift,  max = 1 and 3.  max = 3 corresponds to futuristic HI-galaxy surveys (Square Kilometre Array Cosmology Science Working Group et al. 2020).On the other hand  max = 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 (Δ HI Δ HI ), and the cross-power spectra between HI and  (Δ HI ).We compare the signal of removed and unremoved Δ HI , 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 ℓ > 10, the foreground removal does not erase the Δ HI power spectrum; the signal is scaled down by a factor ( clean ) which is close to constant over a range of ℓ modes from 10 to 1000.Hence, in section 4, when cosmological constraints following foreground removal are considered, the estimation of cosmological parameters is based on the signal where ℓ > 10 .We describe the mean signal drop  clean by: From Fig. 4, we see that the higher the maximum redshift of the survey in which we remove the LOS signal, the less is the effect on the cleaned cross-correlation signal, as more radial modes are preserved in the removal process.Fig. 4 and 5 further indicate that the signal in the HI 2-point correlations drops by approximately the same factor  clean across a wide redshift range, if we remove the background noise up to a particular redshift  max , when crosscorrelating to the  field at a fixed redshift.Fig. 5 also implies that  clean ( HI ,   ,  max1 ) <  clean ( HI ,   ,  max2 ) if the maximum redshifts  max1 >  max2 .

FISHER FORECAST
In previous sections, we have presented the theoretical 2-point statistics for the HI-lensing cross-correlation, and have examined the impact of HI foreground removal on the cross-power spectrum.The results indicate that foreground removal reduces the 2-point statistics by a modest factor.
Here we begin the exploration of HI correlations as a tool for cosmological constraints.In particular we will make a Fisher information matrix forecast of this correlation in the case of low instrument noise (but including our foreground subtraction model); this will assess the best-case capacity of this probe to constrain cosmology, when one is dominated by large-scale structure fluctuations in the HI and lensing fields.We will then examine more realistic cases with cut sky and the inclusion of instrumental noise.

The Fisher matrix
The Fisher information matrix is a useful tool to estimate the expected uncertainty in cosmological parameters for forthcoming experiments (Heavens 2003;Tegmark et al. 1997).Assuming that the model parameters   are distributed by a multivariate Gaussian likelihood , the Fisher matrix can be calculated as where L = − ln .The Fisher matrix can be used to obtain the minimum uncertainty (  ) in parameter estimation due to the Cramér-Rao inequality (Mendez et al. 2014;Kamionkowski et al. 2011), which is equivalent to a 68% confidence level.For a dataset where the uncertainties are Gaussian, the Fisher matrix can be calculated by (Tegmark et al. 1997) where  denotes the covariance matrix of the data,   ≡  −1  , , the derivative data matrix    ≡  ,   ,  +  ,    , , and  is an expectation value of the data vector .The comma symbol means the partial derivative operator with respect to the parameter,  , ≡  /  .Note that all derivatives are performed at the maximum likelihood point.
As we expect only small changes in the covariance matrix COV(  ) 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 We calculate the cross-power spectra   (ℓ) by using Eq. 8 with Planck 2018 cosmological parameters (Planck Collaboration et al. 2018).We calculate the covariance matrices of , HI and HIHI from measured cross power-spectra of 35 realisations of the N-body simulation by Takahashi et al. (2017).All the HI temperature fluctuation maps which we use take into account foreground removal.We also calculate the correlation matrices CORR   = COV   / √︁ COV  COV   and show these in Fig. 6; these do not indicate significant correlations between ℓ bins.

Cosmological Constraints for Single-Slice Cross-correlations
In this section, the cosmological constraint viability of  and HI cross-correlations is explored.We first start with the simplest observational configuration, considering only one redshift slice of HI and .We further assume that the  HI () behaves as in Eq. 17.
We use the Planck 2018 cosmological parameters as the fiducial cosmology (Planck Collaboration et al. 2018).The fiducial cosmological parameters are ℎ = 0.67, Ω m = 0.3,  8 = 0.82, Ω k = 0, Ω Λ = 0.7,  = 0.06, and   = 0.96.To make a covariance matrix of crosspower spectra for the Fisher matrix (Eq.22), we combine ℓ modes into 15 bins; each bin contains 101 ℓ modes with 11 ≤ ℓ ≥ 1527 and averages over 35 realisations.We first consider the 3×2 functions for a joint analysis of (0.78)(0.78),Δ HI (0.3)Δ HI (0.3) and (0.78)ΔHI (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×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  − (green) and HI-HI(grey) provide relatively weak constraints, while 2x2pt and particularly 3x2pt are more promising, with fewto ten-per cent constraints available on parameters in this low noise case.The zoom-in version of 3 × 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  or Δ HI Δ HI constraints alone.Therefore in the next section, we will examine a joint likelihood between more redshift bins, and where the HI bias ( HI ()) is taken into account.

HI bias and multi-redshift bin joint likelihood analysis
It is well known that there is a degeneracy between galaxy bias, Ω m and  8 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 Δ HI Δ HI measures  2 HI (), Δ HI additionally measures  HI ().Combining the cross-bin intensity mapping power spectra with the Δ HI cross-spectra we can therefore tighten our constraints on bias and cosmological parameters.
In this section we consider two different bias models, with distinct parameter sets.The first, more restricted model explores bias amplitude variation via the  parameter in Eq. 17, setting the rest of the parameters to the best fit values (Cunnington et al. 2019b).In the second model,  0 ,  1 and  2 are the bias parameters with  set to equal 1.We include these parameters when evaluating the Fisher matrices (Eq.22).Note that both  HI models are scale-invariant and depend only on .We will consider both full-sky and 300 deg 2 surveys to explore the viability of HIHI and HI in cosmological constraints.
For the full-sky case, as we include more parameters for  HI , we also examine more redshift bins for both HI and  to obtain the best possible results.We consider the redshift range for Δ HI 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).Ta-ble 1 shows the central redshifts we consider for Δ HI and  bins; the width of each bin is 150 ℎ −1 Mpc, Δ ≈ 0.05.Table 1 lists both HI and  central redshifts.As we have 16  HI , we shall refer to '16-HIHI' which corresponds to 16 pairs of HI auto-correlation functions; we cross-correlate HI intensity maps to  fields at   = 0.44, 0.78 and 1.77 respectively.We refer to 16-HIHI+1- as the joint analysis for 16-HIHI and HI 2-point statistics at which   = 0.44.We add further   bins and label joint data as 16-HIHI+2- and 16-HIHI+3-.We calculate both the futuristic case where HI can be measured with high ℓ max ≥ 1000 and the current state of art where 100 < ℓ max < 400.
We further calculate the figure of merit (FoM) for the Ω m −  8 constraint.The FoM is the inverse of the area of the Ω m −  8 contours; in this case we calculate the FoM at 95% confidence level.Fig. 8 shows the FoM of the Ω m −  8 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  < 1 HI slices are available).The green, red and black dots in Fig. 8 show the FoM for ℓ max = 1530 when  we further add the cross-correlations between consecutive HI slices and the  slices at  = 0.44, 0.78 and 1.78, respectively.We see that these cross-correlations significantly improve the FoM, and appear to be converging to a maximal constraint when including all slices.
We now present Fisher forecast results for the first bias model, using the redshift bins in Table 1.Fig. 10 shows the utility of multiredshift bin power spectra measurements with HI and .With the appropriate redshift bin size of HI (Δ HI ) for ℓ max = 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  redshift slices, the constraints are improved significantly especially for Ω m and .However, this In contrast the second model considers quadratic parameters;  2 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 also makes the contours more elliptical, as there are remaining degeneracies among parameters.The uncertainties on parameters are measured at 95% confidence level and reported in Table 2.We next consider the current state of the art case, where since the typical angular resolution ∼ 1 • , we set ℓ max = 375 (we choose this particular value as it is convenient to consider the Δℓ bins as 15 bins with Δℓ = 25).Fig. 9 shows the cumulative Figure of Merit in this case, while Fig. 11 illustrates the likelihood contours for cosmological parameters with this maximum multipole.Comparing with Fig. 10, where ℓ max = 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  bins to 16-HIHI (green shade) in ℓ max = 375.We are therefore seeing that by joining HI 2-point statistics to HIHI auto-correlations, we can improve cosmological constraints significantly.
For bias model 2, we find that the second order coefficient of HI bias  2 is very poorly constrained.Marginalising over this parameter does not significantly affect the other cosmological parameter constraints (ℎ 0 , Ω m ,  8 and   ).We set this parameter to 0.05 following Cunnington et al. (2019b).Fig. 12 shows the likelihood constraints for this model.We can see that by adding more  slices, the ℎ 0 constraints do not improve much but the improvement in the Ω m −  8 constraint can be easily noticed.The uncertainties on our HI bias models and cosmological parameters are reported in Figure 10.The likelihood contours for our multi-bin analysis with HI bias model 1 with ℓ max = 1530; the contours show 68% and 95% confidence levels.The cosmological parameter uncertainties at 95% are measured and reported in Table 2.  Constraints on cosmological parameters and  HI () for our second bias model, for 2-point functions at 68% and 95% levels of confidence (for our 16HIHI redshift slice case).The HI correlation functions do not significantly improve the ℎ 0 ,   ,  0 and  1 constraints.For this figure, we marginalised over  2 as it is poorly constrained.However, we see significant improvement in Ω m and  8 constraints.The parameter uncertainties at 95% level are reported in Table 2.
Table 2. Comparing the parameter constraints from both  HI models (see Table 2), we notice that model 1 gives slightly better (but very comparable) cosmological constraints.

The effect of sky coverage
Now let us consider the effect of sky coverage area on 2-point statistics of both HIHI auto and HI cross angular power spectra.We now consider a combined survey area 300 deg 2 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ℓ ) as an estimator of  ℓ .
Suppose the survey footprint of the observations can be expressed using the weight function  ( n).Normalised by the sky factor  sky (the fraction of the sky covered by the data), the weight moments are given by: where   represents the -th moment of weighting.The power spectrum of the window function is: For a spin-0 field  ( n) weighted by  ( n), a spherical harmonic coefficient ãℓ can be expressed as (Kim & Naselsky 2010;Kim 2011): where we approximate the integration over sky factor by the summation over pixel area with the surface density Ω  .The pseudo power spectrum estimator, Cℓ , is then: Similarly for spin-2 fields ( 1 ,  2 ), we can obtain the coefficients, ã±2,ℓ , by: where Similarly to the full-sky case, ãE ℓ and ãB ℓ are then: The pseudo power spectra for E and B modes are: The pseudo power spectra Cℓ and true  ℓ are related by the modemode coupling resulting from masking ( ℓℓ ′ ): This kernel depends solely on the geometry of a cut-sky W ℓ and plays a crucial role in the pseudo- ℓ method.For details concerning the mode-mode coupling, see Hivon et al. (2002); Alonso et al. (2019).We utilise NaMaster (Alonso et al. 2019), which is a software package to calculate pseudo- ℓ for any spin fields, to evaluate HIHI and HI pseudo- ℓ for 300 deg 2 of our simulations above, within a masked region RA = [0, 30] deg and dec = +[0,10] deg, for the same redshift bins (see Table 1).We choose the same ℓ bins as before up to ℓ < 375 with  beam = 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 deg 2 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).

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.

Single Dish Thermal Noise
We begin by examining the noise on 3D measured power spectra for the HI auto-correlation ().This discussion is based on the works of (Battye et al. 2013;Santos et al. 2015;Bigot-Sazy et al. 2015).Subsequently, we derive the root mean square (rms) thermal noise expected on IM maps from the power spectra, which we will use to assess the effect of realistic noise on the ability to detect the lensing-HI cross-correlation.
The expected uncertainty (  ) on the power spectrum  can be estimated by evaluating its expected second moment.By calculating the ratio between   and  averaging over the radial wavenumber bin size Δ, we can estimate the uncertainty in the IM map measurements.Following this procedure the error on  can be estimated by the following expression (Seo et al. 2010;Feldman et al. 1994;Battye et al. 2013): where  sur is the volume of survey,  () is windows function, T () is the average temperature at given redshift  and  pix is the pixel volume.In this work we shall ignore the contribution of shot noise and only consider the contribution from pixel (thermal) noise  pix .These parameters are discussed in detail in the following paragraphs.
It is essential to pick the proper size of Δ to optimise the viability of a single dish.As we focus on the cosmic signal the acoustic scale should be an aim.That means we require Δ/  < 1, where   is the wavenumber of acoustic scale.The volume of the survey  sur can be computed by where where we assume a flat universe.The window function  () is set by the instrument's specification.As we map the  HI from multiple redshift bins, we can ignore the contribution of radial directions in  ().However for the angular direction, this is not the case.Importantly, the angular resolution of the radio beam will define the noise level.We can model  () by: where  FWHM is the full width half max of angular resolution.
Realistically,  FWHM depends on frequency ().However, in this analysis we assume the variation of  FWHM is small and negligible.The pixel volume,  pix , can be determined by where  c is the central redshift and Δ corresponds to frequency width (Δ) and Ω pix 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  pix can be approximated by (Santos et al. 2015;Seo et al. 2010), where   represents the observation time per pointing,  (approximately 1) signifies the efficiency of the telescope, meaning that almost no signal is lost when radio radiation is transmitted to the antenna. sys represents the system temperature, which includes where we ignore the contribution from the Earth's atmosphere. spl is the spill over from ground radiation (approximately 3K),  CMB ≈ 2.73K and galactic temperature  gal ≈ 25K (408 MHz/) 2.75K (Square Kilometre Array Cosmology Science Working Group et al. 2020).The observing time per pointing   relates to the total observation by   =  tot ( B ) 2 /Ω sur , where  B is an angular pixel size.
As  pix is the rms thermal noise, by definition its square is the power per pixel volume (  / pix ).Therefore, the 3D noise power spectrum (  ) is then (Battye et al. 2013;Santos et al. 2015) where If we have  d dishes where each dish has  b beams, we can take less time for each pointing area.The noise power spectrum is then reduced to To determine the optimisation of a survey strategy, we can estimate the suitable  FWHM and Ω sur that minimises  A / A for acoustic scale  A .This acoustic scale  A can be estimated by following the work of Blake & Glazebrook (2003); Battye et al. (2013), where  is the overall amplitude which can be marginalised.The subscript 'ref' refers to reference cosmological parameters.We now calculate the thermal noise of a MeerKAT-like instrument  MK pix =  T .In this case we consider a single dish telescope consisting of 64 dishes of 13.5m diameter, operating in UHF, L and S-band.The MeerKAT pilot survey by Wang et al. (2021) focuses on L-band from 856 to 1712 MHz with 4096 frequency channels.This pilot survey has 10.5 hours observation time with approximately ∼ 200 deg 2 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  T ≈ 2 mK.
If we use Eq. ( 38) and (39) together with Table 3, the expected  pix for a single frequency channel of MeerKAT pilot survey is  pix (Δ = 0.2; 10hr) ≈ 15 mK, (44) where we assume each dish has equal efficiency  =1 and consider only a single frequency channel Δ = 0.2 MHz.If we consider the whole frequency range like (Wang et al. 2021) the  T ≈ 2mK.We now consider the case where the total observation time  tot = 1000 hours and 250 frequency channels with Δ = 0.2 MHz.

S/N of 𝜅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  tot can take longer than MeerKAT's pilot survey, and we assume a full-sky survey to estimate the best possible S/N for weak lensingintensity mapping (HI) 2-point statistics.The estimate  T pix 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  T pix as far as possible.
We first model the S/N of HI.We can consider the zero lag noise level for ⟨Δ⟩, i.e.where  and Δ are measured in the same pixel.There is no reason why the statistical noise of  should be correlated with Δ.The noise for the cross-correlation will therefore be proportional to the product of  T and    , where  T is HI thermal noise which can be estimated by Eq. 44.The rms noise for weak

Figure 1 .
Figure 1.Power spectra  ℓ for HI and cross-power between HI and convergence, for radial HI width   = 50 and 150 ℎ −1 Mpc.The effect of the width is less important for the cross-correlation.Different   corresponds to different frequency bandwidths, Δ of the radio data.
Fig. 3 also displays a comparison between theoretical 2-point statistics and the measurements from the mock catalogues.We then measure the covariance matrices COV(  ) of 2-point statistics from 35 realisations.The error bars are the square root of the diagonal elements of COV(  ) of the estimators.The correlation matrix for COV(  ) are shown in Fig. 6.

Figure 2 .
Figure 2. Top: the uncleaned intensity map, Δ orig HI , at  = 0.3 from an example realisation.The fluctuations were measured by assuming  HI () (see Eq. 12 and 17).Middle: the foreground-removed intensity map, Δ clean HI , at the same redshift.The foregrounds were removed by eliminating radial long wavelength modes up to redshift  max = 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 × 5 square degrees (a small patch of the entire sky maps).

Figure 3 .
Figure 3.Comparison between theoretical  ℓ (see Eq. 8) and measured  ℓ from our simulations.Here the lensing convergence is measured at central redshift 0.78 and HI is measured at central redshift 0.3.

Figure 4 .
Figure 4.The ratio between cleaned and uncleaned Δ HI power spectra.Two maximum redshifts ( max ) for foreground removal are considered;  max = 1 corresponds to current and imminent radio dishes meanwhile  max = 3 represents a future SKA survey.

Figure 5 .
Figure 5.The average ratio  clean over ℓ > 10 modes, as a function of redshift of HI slice used in the cross-correlation.

Figure 6 .
Figure 6.The correlation matrices of   ℓ measured from 35 realisations for the cleaned HI signal at central redshift  = 0.3 and  signal at  = 0.78.Left:  , middle: Δ HI Δ HI and right: Δ HI .

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

Figure 8 .
Figure 8. Figure of Merit for  8 − Ω  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 = 1530.

Figure 9 .
Figure 9. Figure of Merit for  8 − Ω  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 = 375.

Figure 11 .
Figure 11.The likelihood contours for our multi-bin analysis with HI bias model 1 with ℓ max = 375; the contours show 68% and 95% confidence levels.Compared to Fig. 10, we can see the difference in 16-HIHI constraint.However, when combining the HI-correlation the constraints are similar to those in Fig. 10.

Figure
Figure12.Constraints on cosmological parameters and  HI () for our second bias model, for 2-point functions at 68% and 95% levels of confidence (for our 16HIHI redshift slice case).The HI correlation functions do not significantly improve the ℎ 0 ,   ,  0 and  1 constraints.For this figure, we marginalised over  2 as it is poorly constrained.However, we see significant improvement in Ω m and  8 constraints.The parameter uncertainties at 95% level are reported in Table2.

Figure 13 .
Figure13.Constraints on cosmological parameters and  HI () for 2-point functions at 68% and 95% levels of confidence, for the small area 300 deg 2 case.This plot indicates that the feasibility of HI (pseudo) 2-point statistics in cosmological constraints is heavily affected by the statistical incompleteness due to having a small-sky survey.

Table 1 .
The central redshifts for intensity and lensing convergence maps in our multi-redshift bin analysis.For Δ HI , we require that  HI < 0.7  .

Table 2 .
Cosmological forecast from Fisher analysis for 16-HIHI + 3- correlations; the uncertainties on cosmological parameters and HI bias are quoted at 95% confidence level.HI bias model 1 considers only the re-scaling parameter .