Gravitational lensing of cosmological 21 cm emission

Abstract

We investigate the feasibility of measuring weak gravitational lensing using 21-cm intensity mapping with special emphasis on the performance of the planned Square Kilometre Array (SKA). We find that the current design for SKA_Mid should be able to measure the evolution of the lensing power spectrum at z ∼ 2–3 using this technique. This will be a probe of the expansion history of the Universe and gravity at a unique range in redshift. The signal to noise is found to be highly dependent on evolution of the neutral hydrogen fraction in the Universe with a higher H i density resulting in stronger signal. With realistic models for this, SKA Phase 1 should be capable of measuring the lensing power spectrum and its evolution. The dependence of signal to noise on the area and diameter of the array is quantified. We further demonstrate the applications of this technique by applying it to two specific coupled dark energy models that would be difficult to observationally distinguish without information from this range of redshift. We also investigate measuring the lensing signal with 21 cm emission from the Epoch of Reionization (EoR) using SKA and find that it is unlikely to constrain cosmological parameters because of the small survey size, but could provide a map of the dark matter within a small region of the sky.

gravitational lensing: weak, cosmology: theory, dark energy

1 INTRODUCTION

21 cm cosmology is a new and exciting area of research with a great deal of potential. Future radio telescopes like the SKA¹ (Square Kilometre Array) will give us access to previously unexplored epochs of the Universe, such as the Epoch of Reionization (EoR), the end of the Dark Ages and the beginning of the Cosmic Dawn (Dewdney 2013).

Another important development in the field has been the H i intensity mapping technique in which the distribution of galaxies can be measured without detecting individual galaxies. Instead, the 21 cm emission is treated as a continuous (unresolved) background, much like the cosmic microwave background (CMB), but extended in the frequency (redshift) dimension. Several groups are planning to use this technique to measure the baryon acoustic oscillations scale at redshifts of order of unity (Chang et al. 2008, 2010; Seo et al. 2010; Masui, McDonald & Pen 2010; Ansari et al. 2012; Battye et al. 2012; Chen 2012; Pober et al. 2013; Smoot & Debono 2014; Bull et al. 2014), an important probe of dark energy (DE).

By combining results from different kinds of surveys (CMB, galaxies, 21 cm), which probe complimentary redshift ranges, we will be able to span a very large observational volume and more precisely investigate the acceleration of the Universe. Of particular relevance to this paper is the possibility of extending the cosmological probes to redshifts between those accessible with galaxy redshift surveys and the high redshift of the CMB. By studying the evolution of expansion and structure formation across a wide range of redshift we can investigate whether DE or modified gravity effects are present at a higher redshift than the standard cosmological model predicts (Copeland, Sami & Tsujikawa 2006; Clifton et al. 2012).

The 21 cm radiation provides an excellent source for gravitational lensing studies. Earlier work (Zahn & Zaldarriaga 2005; Metcalf & White 2009) had shown that if the EoR is at redshift z ∼ 8 or later, a large radio array such as SKA could measure the lensing convergence power spectrum and constrain the standard cosmological parameters – however, these studies assumed a much larger survey area than is currently planned for the largest such experiment, SKA_Low (see Section 3 for details). The authors extended the Fourier-space quadratic estimator technique, which was first developed by Hu (2001) for CMB lensing observations to three-dimensional observables, i.e. the 21-cm intensity field I(θ, z). These studies did not consider 21 cm observations from redshifts after reionization when the average H i density in the Universe is much smaller.

In Pourtsidou & Metcalf (2014), we extended the aforementioned studies to redshifts after reionization, but before those probed by galaxy surveys in the visible bands and showed that lensing can be measured using the H i intensity mapping technique. Here, we develop this concept further. In Section 2, we present a preliminary overview of our study, which spans two distinct regimes (EoR and z ∼ 2–5), and give the formal form of the lensing estimator and lensing reconstruction noise. In Section 3, we study how well the current SKA design will be able to map the lensing convergence at typical EoR redshifts. In Section 4, we investigate weak lensing intensity mapping with an SKA-like interferometer array at redshifts z ∼ 2–3. Different models for the H i mass function and different array configurations are considered. We concentrate on the ability of different SKA phases to probe the lensing signal and calculate the corresponding signal-to-noise (S/N) predictions. In Section 5, we show that these measurements can be used to differentiate between some specific and novel interacting DE models. We conclude in Section 6. Throughout the paper, we adopt a flat Λ cold dark matter (ΛCDM) Universe with the Planck cosmological parameter values (Planck Collaboration XVI 2014), unless otherwise stated.

2 PRELIMINARIES

We will investigate the lensing of 21 cm emission in two distinct regimes. The first is from the EoR and the second is from z ∼ 2 to 5. In the first case, the neutral hydrogen fraction is high and H i gas is not restricted to individual galaxies. In this case, the H i distribution will be approximated as a Gaussian random field that roughly follows the distribution of matter. More sophisticated models are under investigation, but will not be considered here. The biggest source of noise for 21 cm observations at these high redshifts comes from foreground contamination. SKA_Low is being planned to probe this regime. Current plans call for it to survey a 5° × 5° sky area with a frequency range of 40–250 MHz corresponding to redshifts z ∼ 5–35 (Dewdney 2013).

At lower redshift, the H i fraction is much lower and essentially all the H i gas is in discrete galaxies. In this regime, we treat the H i distribution as consisting of discrete sources that are clustered according to the standard CDM paradigm, but also exhibit random Poisson fluctuations. This is the standard way of modelling the distribution of galaxies in redshift surveys. We find that this Poisson or shot noise contribution to the galaxy clustering is important for measuring lensing from 21 cm after (and perhaps during) reionization. The evolution in the H i fraction as a function of redshift is a matter of some debate and speculation and has an important effect on the expected lensing S/N. We will address this uncertainty in Section 4 by adopting several different models.

The quadratic lensing estimator for the CMB or a single redshift slice of the 21 cm emission works by essentially measuring differences in the power spectrum in different regions of the sky which results in correlations in Fourier (or spherical harmonic) modes that would not exist otherwise. For the CMB, the local power spectrum becomes anisotropic because of shearing and the acoustic peaks are slightly shifted to larger (smaller) angular scales by gravitational magnification (demagnification).

The effect of lensing can be divided into a shearing which can make the local power spectrum anisotropic and magnification which isotropically scales the local power spectrum. It can be shown that an isotropic magnification cannot be measured (with a quadratic estimator) if the power spectrum is scale free (C_ℓ ∝ ℓ⁻²). It can also be shown that the shear cannot be measured if the power spectrum is constant (see Bucher et al. 2012 for a nice demonstration of this). As we will see, the projected matter power spectrum in the CDM model is approximately a constant at large scales (small ℓ) and C_ℓ ∝ ℓ⁻² at small scales so to the extent that the H i distribution follows the dark matter (DM) distribution, the quadratic estimator picks up the magnification on large scales and the shear on small scales and some combinations in between.

After reionization the remaining H i resides in discrete galaxies. The galaxies are clustered in the same way as matter on large scales (modulo a bias factor), but on small scales their discreteness enters into the power spectrum as shot noise or Poisson noise. The power spectrum of this component is flat (for equal luminosity sources, C_ℓ = 1/η where η is the number density on the sky) so only the magnification can be measured. Magnifying (demagnifying) a region of the sky reduces (increases) the number density of sources. The result is that the power spectrum will have features on the scale of the lensing.

Unlike the CMB, the 21 cm emission will be observable at many different redshifts. The 3D lensing estimator effectively stacks redshift slices. Since the lensing is coherent for different slices at the same angular position, but the 21 cm emission is statistically independent, the residual effect can be attributed to lensing. This will be made more explicit in the next section.

2.1 The quadratic lensing estimator

Here, we will give a more explicit description of the lensing estimator and noise. The advantage of 21 cm lensing is that one is able to combine information from multiple redshift slices, and that the 21 cm signal extends to far smaller angular scales than CMB fluctuations (Loeb 2004), which are suppressed by Silk damping at ℓ > 3000. The intensity fluctuations, I(θ, z), are expressed in discrete Fourier space in the radial direction owing to the finite band width (wave vector |$k_\parallel = \frac{2\pi }{{\cal L}}j$|⁠, where |${\cal L}$| is the depth of the observed volume) and in continuous Fourier space perpendicular to the line of sight (wave vector |$\boldsymbol {k}_{\perp }=\boldsymbol {l}/{\cal D}$| where |${\cal D}$| the angular diameter distance to the source redshift and |$\boldsymbol {l}$| is the dual of the angular coordinate on the sky). See Appendix A for the connection between estimators in discrete and continuous Fourier space. Considering modes with different j independent, an optimal estimator can be found by combining the individual estimators for different j modes without mixing them.

The quadratic estimator for the lensing potential Ψ is thus of the form:

\begin{equation} \hat{\Psi }(\boldsymbol {L})=\sum _j \sum _{{\boldsymbol \ell }}g({\boldsymbol \ell },\boldsymbol {L},j) \; \tilde{I}_{{\boldsymbol \ell },j}\tilde{I}^*_{{\boldsymbol \ell }-\boldsymbol {L},j}, \end{equation}

(1)

where |$\tilde{I}_{{\boldsymbol \ell },j}$| is the discrete Fourier transform (DFT) of the lensed intensity field |$\tilde{I}({ \boldsymbol \theta },x_\parallel )$|⁠. The form of the filter |$g({\boldsymbol \ell },\boldsymbol {L},j)$| and the expected noise in the estimator depends on whether the source is discrete or continuous. In Appendices A through C, we derive the cases that are used in the following sections.

The mean observed brightness temperature at redshift z due to the average H i density is |$\bar{T}(z)=180 \, \Omega _{\rm H\,\small {I}}(z) \, h\frac{(1+z)^2}{H(z)/H_0} \; {\rm mK}$| (Battye et al. 2012). The angular power spectrum of the displacement field (the displacement of points on the sky caused by gravitational lensing) for sources at a redshift z_s can be calculated with

\begin{equation} C_{\rm L}^{\delta \theta \delta \theta } = \displaystyle\frac{9 \Omega ^2_{\rm m} H^3_0}{L(L+1) c^3} \int _0^{z_{\rm s}}\mathrm{d}z \, P(k=L/{\cal D}(z),z) \displaystyle\frac{[W(z)]^2}{a^2E(z)}, \end{equation}

(2)

(Kaiser 1992) in the weak field limit where |$W(z)=({\cal D}(z_{\rm s})-{\cal D}(z))/{\cal D}(z_{\rm s})$|⁠, E(z) = H(z)/H₀ and P(k) is the power spectrum of matter. The convergence power spectrum is simply |$C^{\kappa \kappa }_{\rm L} = (L^2/4) \, C_{\rm L}^{\delta \theta \delta \theta }$|⁠. Note that the Limber approximation (Limber 1953) has been used in deriving equation (2). This approximation is quite accurate even for large scales ℓ ∼ 10, as the resulting formula involves an integral over the line of sight which effectively averages out the effect of the approximation, see Schmidt et al. (2008) for a detailed discussion.

The lensing reconstruction noise includes the thermal noise of the array which – assuming a uniform telescope distribution – is calculated using the formula

\begin{equation} C^{\rm N}_\ell = \displaystyle\frac{(2\pi )^3 T^2_{\rm sys}}{B t_{\rm obs} f^2_{\rm cover} \ell _{\rm max}(\nu )^2} \,, \end{equation}

(3)

where T_sys is the system temperature, B is the chosen frequency window, t_obs the total observation time, D_tel the diameter (maximum baseline) of the core array, ℓ_max(λ) = 2πD_tel/λ is the highest multipole that can be measured by the array at frequency ν (wavelength λ) and f_cover is the total collecting area of the core array A_coll divided by π(D_tel/2)² (Furlanetto, Oh & Briggs 2006).

Assuming that the temperature (i.e. the 21 cm emissivity) is Gaussian distributed, we find that the lensing reconstruction noise N(L) is (see Zahn & Zaldarriaga 2005 and Appendix A)

\begin{eqnarray} &&{N(L,\nu ) =}\nonumber\\ &&{\quad\left[\sum _{j=1}^{j_{\rm max}} \displaystyle\frac{1}{L^4}\int \displaystyle\frac{\mathrm{d}^2\ell }{(2\pi )^2} \displaystyle\frac{[\boldsymbol {l} \cdot \boldsymbol {L} C_{\ell ,j}+\boldsymbol {L} \cdot (\boldsymbol {L}-\boldsymbol {l}) C_{|\ell -L|,j}]^2}{2 C^{\rm tot}_{\ell ,j}C^{\rm tot}_{|\boldsymbol {l}-\boldsymbol {L}|,j}}\right]^{-1},} \end{eqnarray}

(4)

where |$C^{\rm tot}_{\ell ,j}=C_{\ell ,j}+C^{\rm N}_\ell$| and |$C_{\ell ,j}=[\bar{T}(z)]^2P_{\ell ,j}$| with P_{ℓ, j} the underlying DM power spectrum. One can also consider stacking multiple bands ν and then the noise is reduced as N_L = 1/∑_ν[N(L, ν)]^{− 1}. The Gaussian approximation formalism will be used in the next section to investigate the possibility of measuring the lensing signal with 21 cm emission from the EoR.

3 LENSING STUDIES USING 21 CM RADIATION FROM THE EoR

Observing the EoR through 21 cm has been one of the primary science goals for several large radio telescope projects including the SKA. The collecting area and resolution of LOFAR² or MWA³ are unlikely to be high enough to observe gravitational lensing as we will see. For this reason, we concentrate on the SKA_Low instrument here which is the part of SKA designed for observing the EoR and the Cosmic Dawn.

The current SKA plans call for a 25 square degree survey⁴ with SKA_Low (Dewdney 2013). The convergence power spectrum measured in this size field will be dominated by sample variance. This will make it difficult to measure the cosmological parameters through linear growth in the matter power spectrum at a competitive level with these observations even if the S/N in the lensing map is very high. This is not true of the SKA_Mid at lower redshift where the survey area will be much larger – see Section 4.

However, it might be possible to map the lensing convergence within the 25 square degree EoR survey area with high fidelity. This would allow us to actually ‘see’ the distribution of DM in a typical region of the sky, something that is only possible with galaxy lensing around very atypical, large galaxy clusters.

For the EoR, the convergence (or, equivalently, the displacement field) estimator and the corresponding lensing reconstruction noise are calculated assuming that the temperature (brightness) distribution is Gaussian, and consider the reionization fraction |$f_{\rm H\,\small {I}}=1$| until the Universe is rapidly and uniformly reionized at redshift z_s. Before reionization, this is probably a very good approximation. Reionization is not expected to be homogeneous and may extend over a significant redshift range so for some period there will be ionized regions that grow and intersect until they fill almost all of space when reionization is complete. Because the contrast in the brightness temperature is larger during this period, detecting and mapping the 21 cm emission will be easier, but it will change the power spectrum and increase the fourth order moments of the brightness temperature (Metcalf & White 2009). It is not yet clear how this will affect the detectability of lensing. To assess this will require numerical simulations which we plan to present in a future publication. Note that a semi-analytic model used to describe the patchy regime in Zahn & Zaldarriaga (2005) showed an increase in noise due to a decrease of the fluctuation level of the 21 cm signal on the smallest resolved scales and the contribution of the connected four-point function which acts as a sample variance term correlating different k-modes. Of particular importance is the non-linear structure effect during the EoR. We will be able to account for these effects once more details about the reionization process are known, but under the assumption that these features appear in higher resolution than that of our considered interferometers there is no need for the connected four-point contribution to the variance of the quadratic lensing estimator to be taken into account (Kovetz & Kamionkowski 2013).

At redshift z_s ∼ 8, we assume the SKA1 Baseline Design (Dewdney 2013) parameters of A_coll ≃ 0.3 km² with maximum baseline D_tel = 4 km, while for SKA_Mid Phase 2 (SKA2), we can consider A_coll ≃ 1.2 km². The estimated lensing noise is shown in Fig. 1 along with the estimated signal. Here, |$C^{\delta \theta \delta \theta }_{\rm L}$| is the convergence field power spectrum at z_s = 8 and N_L the lensing reconstruction noise assuming a reionization fraction |$f_{\rm H\,\small {I}}=1$|⁠. Note that at such high redshifts, the system temperature T_sys is dominated by galactic synchrotron radiation and can be approximated by (Dewdney 2013)

\begin{equation} T_{\rm sys} = 60 \times \left(\displaystyle\frac{\nu }{300 \,{\rm MHz}}\right)^{-2.55} \, {\rm K}. \end{equation}

(5)

Figure 1.

The lensing displacement field power spectrum, |$C^{\delta \theta \delta \theta }_{\rm L}$|⁠, for sources at z = 8 is shown as a solid black line and lensing reconstruction noise N_L as dashed/dotted lines. The blue dotted curve is for the SKA1 Baseline Design with 108 MHz frequency bins around z = 8 spanning the redshift range z ≃ 6.5–11. The red dashed line is for SKA2 and the same frequency bins. The vertical line is approximately the lowest L accessible with a** 5-by-5 degree field. Where the noise curves are below C_L, typical fluctuations in the lensing deflection should be recoverable in a map.

For SKA_Low, we can consider a 1000-h observation time and multiple stacked bands ν with B = 8 MHz.⁵ Here, we note that the noise N(L, ν) converges quickly with j (as going to higher values of j the signal decays quickly below the thermal noise level), so that j_max ∼ 20 at z ∼ 8.

It can be seen in Fig. 1 that the noise for at least SKA2 should be well below the power spectrum over a large range in angular scales. This means that typical modes will be measurable with high S/N and a high fidelity map is possible. The noise crosses the power spectrum at l ∼ 800 for SKA2 which corresponds to an angular resolution of ∼10 arcmin. This measurement greatly benefits from the larger collecting area that will come with Phase 2 of the SKA. Such a high fidelity map could be cross-correlated with many other observables from lower redshift to determine how gas and stars follow DM as a function of redshift or to test the predictions of general relativity.

4 WEAK LENSING WITH H i INTENSITY MAPPING AFTER REIONIZATION

After reionization, H i gas exists only in discrete galaxies. This introduces some differences in the calculation of the lensing estimator. As will be seen, the noise per mode will increase from the EoR case, but because planned 21-cm intensity mapping surveys will cover a much larger fraction of the sky the sample variance is much lower and a very good measurement of the convergence power spectrum and its redshift evolution should be possible. The S/N is critically dependent on the array parameters, the area surveyed and the evolution of the H i mass fraction with redshift.

4.1 Measuring lensing with the SKA

In previous work (Pourtsidou & Metcalf 2014), we extended the 21-cm intensity mapping lensing method further, taking into account the discreteness and clustering of galaxies, and investigated the possibility of measuring lensing at intermediate redshifts without resolving (in angular resolution) or even identifying individual sources. Here, we perform an improved analysis of the S/N expected from an SKA-like interferometer, and demonstrate the performance of SKA_Mid Phases 1 and 2. A formal treatment of how a quadratic lensing estimator can be derived for discrete sources is included in Appendices B and C.

The expected measurement error in the power spectrum |$C_{\rm L}^{\delta \theta \delta \theta }$| in a band of width ΔL is given by

\begin{equation} \Delta C_{\rm L}^{\delta \theta \delta \theta } =\sqrt{\displaystyle\frac{2}{(2L+1)\Delta L f_{\rm {sky}}}}\left(C_{\rm L}^{\delta \theta \delta \theta }+N_{\rm L}\right). \end{equation}

(6)

Here, N_L is the derived δθ estimate reconstruction noise (see Appendix C for the detailed derivation and formulas for this case) which involves the underlying DM power spectrum, the H i density |$\Omega _{\rm H\,\small {I}}(z)$| as well as the H i mass moments up to fourth order and, of course, the thermal noise of the array |$C^{\rm N}_\ell$|⁠, equation (3). The system temperature in this case is the sum of the receiver and sky noise (Dewdney 2013)

\begin{equation} T_{\rm sys} = \left[40+1.1\times 60 \times \left(\displaystyle\frac{\nu }{300 \,{\rm MHz}}\right)^{-2.55}\right] \, {\rm K}. \end{equation}

(7)

The moments of the H i mass function (or 21-cm luminosity function) enter into the estimator's noise expressions (C15) through (C20). Assuming that the H i mass function follows the Schechter function (which is an excellent fit to the local Universe data – see for example, the HIPASS survey results in Zwaan et al. 2003), we can write

\begin{equation} \displaystyle\frac{{\rm d}n(M,z)}{{\rm d}M}{\rm d}M = \phi ^\star (z) \left(\displaystyle\frac{M}{M^\star (z)}\right)^\alpha {\rm exp}\left[-\displaystyle\frac{M}{M^\star (z)}\right]\displaystyle\frac{{\rm d}M}{M^\star (z)}. \end{equation}

(8)

Damped Lyman α systems observations are used to measure |$\Omega _{\rm H\,\small {I}}(z)$| up to z ∼ 5, while no measurements of ϕ^⋆(z) and M^⋆(z) are available other than in the local Universe. However, for this form of the mass function we have the relation

\begin{equation} \Omega _{\rm H\,\small {I}}(z) = 4.9 \times 10^{-4} \displaystyle\frac{\phi ^\star (z)}{\phi ^\star (z=0)}\displaystyle\frac{M^\star (z)}{M^\star (z=0)}. \end{equation}

(9)

At redshift z = 0, we use the values α = −1.3, M^⋆ = 3.47 h⁻²10⁹ M_⊙, ϕ^⋆ = 0.0204 h³ Mpc⁻³ reported from the HIPASS survey (Zwaan et al. 2003). We will initially consider the most conservative scenario, i.e. a no evolution model that we will call Model A. Model A was also used in Pourtsidou & Metcalf (2014). Here, we will also show results using a few H i evolution models.

Let us first concentrate on source redshift z_s ∼ 2. In the recently published SKA1 Baseline Design (Dewdney 2013), SKA_Mid operates in the frequency range ∼350–3050 MHz divided in three bands. Our chosen z_s corresponds to frequency ν_s = 473 MHz (Band 1). Following the proposed parameters for SKA_Mid, we consider a 2-yr observation time, f_sky ∼ 0.7, and we further choose B = 20 MHz and ΔL = 36. From equation (B6), we see that the measurement errors increase as |$1/\sqrt{f_{\rm sky}}$| and the fact that SKA_Mid will cover a large fraction of the sky contributes to the high S/N significance we predict. An increase in the observation time would reduce the thermal noise of the interferometer and would also lead to an increase of the overall S/N levels. A similar effect comes from increasing the bandwidth B of the observation – however, note that we avoid choosing very wide bands since in that case there would be non-negligible correlations within the band. It is not clear whether these correlations would just lead to an increase of the noise only or if they would also increase the lensing signal. We plan to investigate this issue in a future publication.

Since the galaxies are being treated as point sources, we wish to include only frequency modes that do not resolve the internal structure of the galaxies. Assuming a typical velocity dispersion for a galaxy at z = 2 to be around 200 km s⁻¹, we find j_max = 63. However, the noise has converged at j ∼ 40, which is the value we use in our numerical simulations (see Fig. C1 for a demonstration of the convergence with j). Keeping these values constant, we present a S/N contour plot (Fig. 2) at multipole L = 100 for z_s = 2 on the (A_coll, D_tell) parameter space, and show three SKA_Mid performance cases: SKA_Mid Phase 1 (SKA 1) with 50 per cent sensitivity (SKA0), SKA1 and SKA2. We also present the displacement field power spectrum and measurement errors for z_s = 2 and 3 using the SKA2 specifications in Fig. 3. These results are not sensitive to the exact value of ℓ_min. For example, the S/N values remain practically unchanged if we use a very small ℓ_min = 10 or the one corresponding to the current SKA1 field of view (∼1 deg), ℓ_min = 180.

Figure 2.

The S/N at L = 100 in a ΔL = 36 for various array configurations. Sources are at z_s = 2. The contour lines are labelled with the (S/N) values. The area under the thick blue solid line is excluded, since it corresponds to f_cover > 1. SKA0, SKA1 and SKA2 are shown. The no-evolution in |$\Omega _{\rm H\,\small {I}}$| (Model A) scenario is used.

Figure 3.

Displacement field power spectrum for z_s = 2 (solid black line) and z_s = 3 (dashed magenta line) and the corresponding measurement errors using the SKA2 specifications and Model A for the H i mass function.

This technique should enable us to measure the lensing power spectrum at source redshifts well beyond those accessible with more traditional weak lensing surveys based on the shearing of individual galaxy images. The noise is expected to be small enough that evolution in the lensing power spectrum between z_s = 2 and 3 will be clearly detectable as illustrated in Fig. 3.

4.2 Importance of the redshift evolution in the H i mass function

It is useful to illustrate how our results depend on a possible evolution of the H i mass function with redshift. We will follow Zhang & Pen (2006). The authors construct three different evolution models for the H i mass function. Observations of damped Lyman α systems and Lyman α limit systems measure |$\Omega _{\rm H\,\small {I}}$| from z = 0 to z ∼ 5 (see for example, Peroux et al. 2003).These observations found that |$\Omega _{\rm H\,\small {I}}$| increases by a factor of 5 towards z ∼ 3 and then decreases towards higher redshift. Combining equation (9) and these observations, one can put constraints on the product of ϕ^⋆ and M^⋆. Either an increase in ϕ^⋆ or M^⋆ increases the detectability of lensing, so the no evolution model is considered a conservative choice.

Following (Zhang & Pen 2006), the evolution of |$\Omega _{\rm H\,\small {I}}$| will be approximated as

\begin{equation} g(z) \equiv \displaystyle\frac{\Omega _{\rm H\,\small {I}}(z)}{\Omega _{\rm H\,\small {I}}(z=0)} = (1+z)^{2.9} \, {\rm exp}(-z/1.3). \end{equation}

(10)

This function is plotted in Fig. 4 and the constraint for the redshift evolution of the product ϕ^⋆M^⋆ is

\begin{equation} \phi ^\star (z)M^\star (z)=\phi ^\star (z=0)M^\star (z=0)g(z). \end{equation}

(11)

Here, we will investigate the following evolution scenarios (Zhang & Pen 2006):

• Model (B): No evolution in M^⋆(z). ϕ^⋆(z) = ϕ^⋆(z = 0)g(z).

• Model (C): No evolution in ϕ^⋆(z). M^⋆(z) = M^⋆(z = 0)g(z).

• Model (D): ϕ^⋆(z)/ϕ^⋆(z = 0) = M^⋆(z)/M^⋆(z = 0) = g(z)^1/2.

Figure 4.

The evolution of |$\Omega _{\rm H\,\small {I}}(z)$| with redshift in Zhang & Pen (2006).

We also note that H i mass function modelling has been performed by other authors (Abdalla, Blake & Rawlings 2010), where the evolution of the |$\Omega _{\rm H\,\small {I}}(z)$| function agrees with Zhang & Pen (2006) until z ∼ 3 and so would make very little difference to our predictions.

4.2.1 Model B

Model B assumes no evolution in M^⋆. The shot-noise terms in the estimator noise are changed in this case. The displacement field power spectrum and measurement errors for Model B are shown in Fig. 5, and the S/N contour at L = 100 is shown in Fig. 6 as a function of array parameters. Notice the improvement of the S/N across all of the parameter space from the no-evolution model – the lensing signal can be detected even with SKA0.

Figure 5.

Displacement field power spectrum for z_s = 2 (solid black line) and z_s = 3 (dashed magenta line) and the corresponding measurement errors using the SKA2 specifications and Model B for the H i mass function.

Figure 6.

The S/N at L = 100 in a ΔL = 36 for various array configurations. Model B for the evolution of the H i mass function is used. Sources are at z_s = 2. The contour lines are labelled with the (S/N) values. The area under the thick blue solid line is excluded, since it corresponds to f_cover > 1. SKA0, SKA1 and SKA2 are shown.

4.2.2 Model C

Model C assumes no evolution in ϕ^⋆. |$\Omega _{\rm H\,\small {I}}$| will increase, but the shot-noise terms will stay the same as the no-evolution ones (Model A), because the integrals do not get affected by M^⋆. Hence, the S/N is higher than Model A by a factor of ∼1.5 at L = 100 using the SKA2 specifications, but lower than Model B by a factor of ∼0.5 at L = 100 using the SKA2 specifications.

4.2.3 Model D

Model D assumes that both ϕ^⋆ and M^⋆ are evolving as ϕ^⋆(z)/ϕ^⋆(z = 0) = M^⋆(z)/M^⋆(z = 0) = g(z)^1/2. The S/N is higher than Model C by a factor of ∼1.5 at L = 100 using the SKA2 specifications, but lower than Model B by a factor of ∼0.7 at L = 100 using the SKA2 specifications.

We see that for all models where the H i density increases with redshift the lensing S/N is greater than in the no-evolution model. Furthermore, we learned that there is a stronger dependence on the evolution of ϕ^⋆ than M^⋆. This is also expected, since an increase in ϕ^⋆ increases the detectability of lensing – at this point we should stress that in our case the contributions of the Poisson moments contribute both to the signal as well as the noise of our estimator, so the dependence on ϕ^⋆ is crucial. Model B, in which ϕ^⋆ is solely responsible for the increase in |$\Omega _{\rm H\,\small {I}}$| is the most optimistic scenario, but the most conservative no-evolution Model A is probably less realistic.

Previously, we used the S/N calculations at a single multipole number (L = 100) to demonstrate the capabilities of various arrays (see Figs 2 and 6). A more comprehensive expression would involve contributions from many multipoles. We therefore define a new (S/N) for the detection of the lensing signal as

\begin{equation} {\rm S/N} = \left(\sum _{L_{\rm b}} \displaystyle\frac{C^2_{L_{\rm b}}}{\Delta C^2_{L_{\rm b}}}\right)^{1/2}, \end{equation}

(12)

where the sum is over deflection field power spectrum band powers (using ΔL = 36 as before) and |$\Delta C_{L_{\rm b}}$| is the forecasted error on each band power. We calculate (S/N) as a function of collecting area and maximum baseline, with |$L_{b_{\rm min}}=18$| and |$L_{b_{\rm max}}=990$|⁠. We show the S/N contour plots for Models A and B for the H i mass functions in Figs 7 and 8, respectively.

Figure 7.

The S/N from equation (12) for various array configurations. Sources are at z_s = 2 and we use Model A for the H i mass function. The contour lines are labelled with the (S/N) values. The area under the thick blue solid line is excluded, since it corresponds to f_cover > 1. SKA0, SKA1 and SKA2 are shown.

Figure 8.

The S/N from equation (12) for various array configurations. Sources are at z_s = 2 and we use Model B for the H i mass function. The contour lines are labelled with the (S/N) values. The area under the thick blue solid line is excluded, since it corresponds to f_cover > 1. SKA0, SKA1 and SKA2 are shown.

As we have already stressed, Model A is the most conservative, no-evolution scenario one can consider for modelling the H i mass function. From Figs 3 and 7, we can deduce that in order to get a good measurement of the lensing signal up to multipole number L ∼ 1000 we need a S/N of about 25. Then from Fig. 8, we can see that SKA1 assuming Model B does almost as well as SKA2 assuming no evolution (sources at z_s = 2). This shows how crucial the H i mass function evolution is for our method and the SKA's science goals. In Fig. 9, we show the results of our calculations using the SKA1 specifications and Model B – as expected from the above discussion, the measurement errors are similar to the SKA2 no-evolution case.

Figure 9.

Displacement field power spectrum for z_s = 2 and the corresponding measurement errors using the SKA1 specifications and Model B for the H i mass function.

5 CONSTRAINING INTERACTING DARK ENERGY MODELS

In this section, we demonstrate the constraining power of 21-cm lensing measurements using the intensity mapping method with an instrument like the SKA. The problems of the nature and evolution of DE and DM are arguably the most important ones in modern cosmology, and ambitious future missions like EUCLID⁶ are dedicated to the exploration and mapping of the dark Universe by investigating the evolution of cosmic structures out to redshifts ∼1.5–2, a period over which DE is thought to dominate the Universe's expansion.

An important question is whether DE is indeed a cosmological constant, like all the available data from CMB measurements (Planck Collaboration XVI 2014), the Hubble constant (Riess et al. 2009) luminosity and distance at high redshift with supernovae Ia (Kowalski et al. 2008) and BAO surveys (Lampeitl et al. 2009) seem to suggest (however, this success comes with the price of the cosmological constant and coincidence problems), an evolving scalar field like quintessence, or a modification of gravity. To differentiate between the plethora of available models, we might have to go even deeper in redshift space z > 1.5, and this is where 21 cm radiation becomes a unique probe of cosmology.

Interactions in the dark sector, such as a non-gravitational coupling between DM and DE (for example a coupling between a quintessence field ϕ playing the role of DE and the matter sector, see Amendola 2000), can cause modifications to the background evolution of the matter density and the Hubble parameter, as well as changes in the evolution of structure growth with respect to ΛCDM. However, as we will demonstrate below, the background effects are very difficult to probe as alternative DE and modified gravity models can be very successful in mimicking the ΛCDM background. Therefore, we need to study linear perturbations and quantify the effects of the various models in the CMB and matter power spectra. Modifications in the evolution of structure growth in interacting DE models leave distinctive signatures in the weak lensing signal, and measurements such as those presented in Fig. 3, could be used to constrain such models.

In a recent study, Pourtsidou, Skordis & Copeland (2013) constructed three general classes (types) of models of DE in the form of a scalar field ϕ coupled to CDM. The first class, Type 1, is a generalization of the coupled quintessence (CQ) model suggested by Amendola (2000). In such models, the CDM energy density is not conserved separately (like in the uncoupled case) but there is a non-zero coupling current which represents the decay of DM to DE or vice versa. In the specific Amendola (2000) model, the evolution of the background CDM density |$\bar{\rho }_{\rm c}$| is found to be

\begin{equation} \bar{\rho }_{\rm c} = \bar{\rho }_{{\rm c},0}a^{-3}{\rm e}^{\alpha \phi }. \end{equation}

(13)

We therefore see that the CDM density evolves differently than the uncoupled quintessence case and ΛCDM, which means that for positive coupling parameter α (decay of DM to DE) the CDM density has to be higher in the past in order to reach the same value today.

This model and its variants have been extensively studied in the literature, and constraints on the strength of the coupling parameter (roughly <0.1 assuming it is constant) have been derived using its effects on the CMB and matter power spectra, the growth of structure and the weak lensing signal (e.g. Caldera-Cabral, Maartens & Schaefer 2009; Amendola et al. 2011; De Bernardis et al. 2011; Tarrant et al. 2012; Xia 2013, and references therein).

Here, we will investigate the weak lensing convergence signatures of models belonging to the new classes of coupled DE theories, namely Type 2 and Type 3, constructed in Pourtsidou et al. (2013). Type 2 models involve both energy and momentum transfer between the two components of the dark sector (i.e. like Type 1 models but with a distinctively different coupling mechanism), while Type 3 is a pure momentum transfer theory. The displacement field power spectrum for these models is calculated using (Kaiser 1992)

\begin{eqnarray} &&{C_{\rm L}^{\delta \theta \delta \theta } = \displaystyle\frac{9}{L(L+1) c^3} \int _0^{z_{\rm s}}\mathrm{d}z \, [\Omega _{\rm m}(z)]^2 P(k=L/{\cal D},z) [W(z)]^2} \nonumber\\ &&{\qquad\quad\ \,\times\, (1+z)^{-4}[H(z)]^3 ,} \end{eqnarray}

(14)

which simplifies to equation (2) for ΛCDM by setting |$\frac{\Omega _{\rm m}(z)}{\Omega _{\rm m}}=\frac{H^2_0}{a^3[H(z)]^2}$|⁠.

5.1 Type 2 models of CDM coupled to DE

Type 2 models are classified via the Lagrangian function (Pourtsidou et al. 2013)

\begin{equation} L(n,Y,Z,\phi ) = F(Y,\phi ) + f(n,Z), \end{equation}

(15)

where n is the particle number density, |$Y=\frac{1}{2}\nabla _\mu \phi \nabla ^\mu \phi$| is used to construct a kinetic term for the quintessence field ϕ and Z = u^μ∇_μϕ plays the role of a direct coupling of the fluid velocity to the gradient of the scalar field.

We can proceed by making further assumptions. First of all, we want to consider coupled quintessence models, hence we can write F = Y + V(ϕ). Secondly, since we choose the scalar field to be coupled to CDM, the function f is separable, i.e. f = n h(Z).

The details of the field and fluid equations for the Type 2 class of theories under consideration can be found in Pourtsidou et al. (2013). An important quantity that appears in the cosmological equations is the function K(Z) = h_z/(h − Zh_z), with h_z = dh/dZ. We therefore see that we cannot simply choose h to be proportional to Z, i.e. h(Z) = βZ with β the (constant) coupling parameter, as in that case K(Z) diverges.

Choosing h(Z) = Z^β, we find that the background CDM density |$\bar{\rho }_{\rm c}$| evolves as

\begin{equation} \bar{\rho }_{\rm c}=\bar{\rho }_{{\rm c},0}a^{-3}\bar{Z}^\beta , \end{equation}

(16)

where |$\bar{Z}=-\dot{\bar{\phi }}/a$| (note that the dot denotes derivatives with respect to conformal time). The fact that the CDM density depends on the time derivative |$\dot{\bar{\phi }}$| instead of |$\bar{\phi }$| itself is a notable difference from the Type 1 (Amendola) case. Note also that the specific form of the coupling function Z^β means that in order to have meaningful solutions |$\bar{Z}$| has to be positive throughout the cosmological evolution. In Pourtsidou et al. (2013), the authors derive the background and perturbed Klein–Gordon equations for the evolution of the quintessence field, as well as the perturbed CDM equations for the evolution of the density contrasts and velocities. With these equations at hand, and considering a single exponential potential for the quintessence field, we use a modified version of the camb code (Lewis, Challinor & Lasenby 2000) to study the observational signatures of the chosen coupled model. In Figs 10 and 11, we show the Hubble parameter evolution and the CMB TT power spectra for these models together with ΛCDM. We also construct the displacement field power spectrum (Fig. 12) and compare it with the ΛCDM prediction in Fig. 3. Note that the cosmology of each model evolves to the Planck cosmological parameter values at z = 0 (Ade et al. 2013) and the matter power spectra P(k) have been normalized to CMB fluctuations.

Figure 10.

The Hubble parameter versus redshift for ΛCDM (solid black line) and the Types 2 and 3 models under consideration (red dashed line for Type 2, blue dotted dashed line for Type 3), together with expansion history measurements from Moresco et al. (2012a), which combine the observational constraints on the Hubble parameter from Simon, Verde & Jimenez (2005), Stern et al. (2010) and Moresco et al. (2012b).

Figure 11.

The CMB TT spectra for ΛCDM (solid black line) and the Types 2 and 3 models under consideration (red dashed line for Type 2, blue dotted dashed line for Type 3), together with WMAP9 measurements.

Figure 12.

Displacement field power spectrum for ΛCDM (solid black line) compared with Type 2 (red dashed line) and Type 3 (blue dotted dashed line) interacting DE models. Sources are at z_s = 2. The circle points with black error bars are the same as in Fig. 3, i.e. using the most conservative scenario (no evolution in the H i mass function), Model A. Assuming a more optimistic scenario like Model B (magenta error bars) would further improve the constraints.

From Fig. 12, we see that the Type 2 model with a coupling parameter β = 0.02 would be excluded by the proposed lensing/intensity mapping observations. That is because in this case there is energy transfer from DM to DE making the DM density larger in the past compared to the non-interacting case for fixed Ω_m today, the growth is increased, the gravitational potential is higher and the convergence power spectrum is enhanced with respect to ΛCDM. Another important point demonstrated in Fig. 13 is that the Type 2 model linear growth factor is strongly scale dependent, and the difference between the growth of the two scales persists up to the source redshift z = 2. Note that there is also an effect on the Hubble parameter evolution (see Fig. 10) and the CMB power spectrum (see Fig. 11), but subtler than the effect on the lensing signal.

Figure 13.

The ratio of the linear matter power spectrum in the Type 2 coupled-DE model to the one in ΛCDM for two characteristic scales as a function of redshift.

5.2 Type 3 models of CDM coupled to DE

Type 3 models are classified via the Lagrangian function (Pourtsidou et al. 2013)

\begin{equation} L(n,Y,Z,\phi ) = F(Y,Z,\phi ) + f(n). \end{equation}

(17)

We can consider a CQ function of the form F = Y + V(ϕ) + h(Z), and choose h(Z) = γZ². Note that this model is physically acceptable for γ < 1/2, otherwise it suffers from ghost and strong coupling problems (Pourtsidou et al. 2013).

As we have already mentioned, Type 3 is a pure momentum transfer theory. The background CDM density |$\bar{\rho }_{\rm c}$| evolves as in the uncoupled case and the energy conservation equation remains uncoupled even at the linear level. Following the same procedure, using the equations derived in Pourtsidou et al. (2013) and the modified version of camb, we construct the displacement field power spectrum and compare it with the ΛCDM prediction in Fig. 3.

From Fig. 12, we see that the Type 3 model with a coupling parameter γ = 0.05 would be excluded. The signature of this model is growth suppression and hence a decrease in the lensing signal. Again, we note that there is also an effect on the Hubble parameter evolution (see Fig. 10) and the CMB power spectrum (see Fig. 11), but subtler than the effect on the lensing signal.

6 CONCLUSIONS

H i intensity mapping is an innovative technique which can be utilized to map the large-scale structure of the Universe in three dimensions. In this paper, we investigated the potential of an SKA-like interferometer to perform accurate measurements of the lensing signal over a wide range of redshifts using the intensity mapping method. Extending the work in Pourtsidou & Metcalf (2014), we presented our technique for measuring gravitational lensing in H i observations after reionization. Other methods have been proposed based on counting the number of detected galaxies in 21 cm emission (Zhang & Pen 2005, 2006; Zhang & Yang 2011). In such methods, the clustering of galaxies and the shot noise from Poisson fluctuations in the number counts contribute purely to the noise in the lensing estimator. The corresponding errors can be controlled to ∼10–20 per cent level for intermediate source redshifts z ∼ 1, but for z_s ≳ 2, the reconstruction is not successful (Zhang & Yang 2011). Our estimator takes into account both the clustering and discreetness of galaxies, which contribute to both the noise and to improving the signal resulting in a significant improvement in performance.

Our calculations showed that the lensing convergence can be measured with high fidelity at redshifts z ∼ 2–3, hence it can be used for tomographic studies along many redshift (or frequency) bins in order to map, for example, the evolution of the growth function at intermediate redshifts, i.e. higher than those of galaxy shear surveys. We also extended our calculations to include a possible evolution of the H i mass function, showing that the less likely no-evolution scenario is the most conservative one and our claim that a high-S/N measurement can be achieved using this method is robust. To quantify the power of this technique, we performed an optimization study using an SKA-like interferometer and showed that SKA will be able to deliver very good measurements of the lensing signal and its redshift evolution. We also demonstrated the constraining power of our technique by using it to distinguish specific interacting DE models. Our results confirm that the 21 cm technique can be used to compete (and complement) with future galaxy surveys, and SKA in particular can be used to investigate the evolution of DE and DM in order to test the standard ΛCDM paradigm against exotic models with possible observational signatures in higher redshifts than the ones probed by galaxy lensing surveys.

Another exciting prospect is the study of the reionization epoch – we showed that the low-frequency SKA instrument can map the lensing convergence and allow us to ‘see’ the distribution of DM in a typical region of the sky, something that is only possible with galaxy lensing around very atypical, large galaxy clusters. This would provide a great opportunity to correlate visible objects with mass and test the DM paradigm.

This research is part of the project GLENCO, funded under the Seventh Framework Programme, Ideas, Grant Agreement no. 259349. The authors would like to thank Leon Koopmans, Jonathan Pritchard and Mario Santos for useful discussions.

1

www.skatelescope.org

2

www.lofar.org

3

www.mwatelescope.org

4

Here, we note that suggestions for a three-tier SKA_Low survey have been recently put forward (Koopmans & Pritchard, private communication). These include an additional 200 square degree survey in 10 × 100 h as well as a 2000 square degree survey in 100 × 10 h. Note that these numbers assume one beam only, and a multibeam request (N_beam = 5) has also been put forward. This would further increase the sky area surveyed in a given amount of time and would greatly enhance the SKA_Low science output, adding the exciting possibility of measuring the cosmological parameters (see McQuinn et al. 2006 and Metcalf & White 2009).

5

As discussed in Metcalf & White (2009), the quadratic lensing estimator is optimized for the case where the statistical properties of the 21-cm radiation signal and noise are constant within a band, and an observation bandwidth of a few MHz is small enough so that this assumption is justified.

6

http://sci.esa.int/euclid/

APPENDIX A: LENSING ESTIMATOR FOR A CONTINUOUS SOURCE

Here, we derive the least squares, or minimum variance, quadratic estimator on a discrete grid for a source that can be approximated as a Gaussian random field. This will be the discrete version of the estimator found by Zahn & Zaldarriaga (2005). The discrete version is necessary for calculations on a finite grid and useful in the deviation of the discrete source case (see Appendix B). This appendix will also serve to introduce some notation.

Figure C1.

Displacement field power spectrum for z_s = 2 and the corresponding lensing reconstruction noise N(L) for different values of j_max using the SKA2 specifications and Model B for the H i mass function.

We write the DFT of the intensity field |$I(\boldsymbol x)$| as

\begin{equation} I_{\boldsymbol k} = \frac{\Omega _{\rm s}}{N_\perp N_\parallel } \sum _{\boldsymbol x}{\rm e}^{{\rm i}\boldsymbol k \cdot \boldsymbol x} I(\boldsymbol x), \end{equation}

(A1)

where |$\boldsymbol k=({\boldsymbol \ell },j)$|⁠, |$\boldsymbol x=(\theta ,z)$| and Ω_s = Θ_s × Θ_s for a square survey geometry. We also have

\begin{equation} I(\boldsymbol x) = \frac{1}{\Omega _{\rm s}}\sum _{\boldsymbol k}{\rm e}^{-{\rm i}\boldsymbol k \cdot \boldsymbol x}I_{\boldsymbol k}. \end{equation}

(A2)

Here, N_∥ is the number of cells in the direction parallel to the line of sight and N_⊥ is the number of cells on a plane perpendicular to the line of sight. We use a flat sky or small angle approximation.

The correlation between discrete modes is

\begin{equation} \langle I_{\boldsymbol k}I^*_{\boldsymbol k^{\prime }}\rangle = \Omega _{\rm s}^2 \delta ^K_{kk^{\prime }} P_{\boldsymbol k}, \end{equation}

(A3)

where |$P_{\boldsymbol k}$| is the discrete power spectrum. The Kronecker delta implies statistical homogeneity. The discrete power spectrum is related to the continuous one, P(k), by

\begin{equation} P_{\boldsymbol k} = \frac{P(k)}{V_{\rm s}} = \frac{P(k)}{\Omega _{\rm s} {\cal D}^2 {\cal L}} = \frac{1}{\Omega _{\rm s}}C_{\ell ,j}. \end{equation}

(A4)

Here, |${\cal D}$| is the comoving angular size distance to the source volume from the observer. The angular Fourier coordinate is |$\ell = {\cal D} \boldsymbol k_\perp$| and j denotes the DFT coordinate in the radial direction. The comoving length of the source volume in the radial direction is |$\cal L$|⁠. This serves to define the angular power spectrum in angle and radial coordinates C_{ℓ, j}. So in terms of the angular power spectrum

\begin{equation} \langle I_{\boldsymbol k}I^*_{\boldsymbol k^{\prime }}\rangle = \Omega _{\rm s} \; C_{\ell ,j} \; \delta ^{K}_{\ell ,\ell ^{\prime }}\delta ^{K}_{jj^{\prime }}. \end{equation}

(A5)

Gravitational lensing causes the observed emission field to be inhomogeneous over a region of the sky with coherent deflection. The observed intensity after lensing is

\begin{equation} \tilde{I}({ \boldsymbol \theta },x_\parallel ) = I\left( \boldsymbol \theta - {\boldsymbol \alpha }(\boldsymbol \theta ), x_\parallel \right) \end{equation}

(A6)

\begin{eqnarray} \nonumber\\ \hphantom{\tilde{I}({ \boldsymbol \theta },x_\parallel ) }\simeq I\left( \boldsymbol \theta , x_\parallel \right) - {\boldsymbol \alpha }(\boldsymbol \theta ) \cdot \nabla _\theta I\left( \boldsymbol \theta , x_\parallel \right) \end{eqnarray}

(A7)

\begin{equation} \hphantom{\tilde{I}({ \boldsymbol \theta },x_\parallel ) }\simeq I\left( \boldsymbol \theta , x_\parallel \right) + \nabla _\theta \Psi (\boldsymbol \theta ) \cdot \nabla _\theta I\left( \boldsymbol \theta , x_\parallel \right), \end{equation}

(A8)

where |${\boldsymbol \alpha }(\boldsymbol \theta )$| is the deflection caused by lensing and |$\Psi (\boldsymbol \theta )$| is the lensing potential. The deflection field is a potential field to very good approximation. In Fourier space, this becomes

\begin{equation} \tilde{I}_{{\boldsymbol \ell },j} \simeq I_{{\boldsymbol \ell },j} + \frac{1}{\Omega _{\rm s}} \sum _{{\boldsymbol \ell }^{\prime }} {\boldsymbol \ell }^{\prime } \cdot ({\boldsymbol \ell } - {\boldsymbol \ell }^{\prime }) I_{{\boldsymbol \ell }^{\prime },j} \Psi _{{\boldsymbol \ell } - {\boldsymbol \ell }^{\prime }}. \end{equation}

(A9)

From this, we can find the correlation between modes to first order

\begin{equation} \langle \tilde{I}_{{\boldsymbol \ell },j}\tilde{I}^*_{{\boldsymbol \ell }-\boldsymbol {L},j^{\prime }}\rangle = \Omega _{\rm s} C_{\ell ,j} \delta ^{K}_{\boldsymbol {L},0}\delta ^K_{jj^{\prime }} + \delta ^K_{jj^{\prime }}[{\boldsymbol \ell }\cdot \boldsymbol {L}C_{\ell ,j}+\boldsymbol {L}\cdot (\boldsymbol {L}-{\boldsymbol \ell }) C_{|\ell -L|,j}]\Psi _{\boldsymbol L}. \end{equation}

(A10)

We seek to construct a quadratic estimator of the form

\begin{equation} \hat{\Psi }(\boldsymbol {L})=\sum _j \sum _{{\boldsymbol \ell }}g({\boldsymbol \ell },\boldsymbol {L},j) \; \tilde{I}_{{\boldsymbol \ell },j}\tilde{I}^*_{{\boldsymbol \ell }-\boldsymbol {L},j}. \end{equation}

(A11)

Putting equation (A10) into (A11) and requiring the estimator to be unbiased gives the constraint

\begin{equation} \sum _j \sum _{{\boldsymbol \ell }}g({\boldsymbol \ell },\boldsymbol {L},j)[{\boldsymbol \ell }\cdot \boldsymbol {L}C_{\ell ,j}+\boldsymbol {L}\cdot (\boldsymbol {L}-\boldsymbol L) C_{|\ell -L|,j}]=1. \end{equation}

(A12)

The variance of this estimator is

\begin{eqnarray} {\cal V} &=& \displaystyle\sum _j \sum _{j^{\prime }} \sum _{{\boldsymbol \ell }} \sum _{{\boldsymbol \ell }^{\prime }} g({\boldsymbol \ell },\boldsymbol {L},j) g^*({\boldsymbol \ell }^{\prime },\boldsymbol {L},j^{\prime })\langle I_{{\boldsymbol \ell },j}I^*_{|{\boldsymbol \ell }-\boldsymbol {L}|,j}I^*_{{\boldsymbol \ell }^{\prime },j^{\prime }}I_{|{\boldsymbol \ell }^{\prime }-\boldsymbol {L}|,j^{\prime }}\rangle \nonumber\\ &=&\displaystyle2\Omega ^2_{\rm s} \sum _j \sum _{{\boldsymbol \ell }}g^2({\boldsymbol \ell },\boldsymbol {L},j)C^{\rm tot}_{\ell ,j}C^{\rm tot}_{|{\boldsymbol \ell }-\boldsymbol {L}|,j}, \end{eqnarray}

(A13)

where in calculating the fourth-order correlations the result for a Gaussian random field has been used and

\begin{equation} C^{\rm tot}_{\ell ,j} = C_{\ell ,j} + C^{\rm N}_{\ell }, \end{equation}

(A14)

with |$C^{\rm N}_{\ell }$| the thermal noise of the array.

With the standard Lagrangian multiplier technique |${\cal V}$| can be minimized subject to the constraint (A12). The result is

\begin{equation} g(\tilde{{\boldsymbol \ell }},\boldsymbol {L},\tilde{j})= A_{\rm R} \frac{[\tilde{{\boldsymbol \ell }} \cdot \boldsymbol {L}C_{\tilde{\ell },\tilde{j}}+\boldsymbol {L}\cdot (\boldsymbol {L}-\tilde{{\boldsymbol \ell }}) C_{|\tilde{\ell }-L|,\tilde{j}}]}{\Omega ^2_{\rm s}C^{\rm tot}_{\tilde{\ell },\tilde{j}}C^{\rm tot}_{|\tilde{{\boldsymbol \ell }}-\boldsymbol {L}|,\tilde{j}}} ,\quad A_{\rm R} =\left[\sum _j \sum _{{\boldsymbol \ell }} \frac{[{\boldsymbol \ell }\cdot \boldsymbol {L}C_{\ell ,j}+\boldsymbol {L}\cdot (\boldsymbol {L}-{\boldsymbol \ell }) C_{|\ell -L|,j}]^2}{\Omega ^2_{\rm s}C^{\rm tot}_{\ell ,j}C^{\rm tot}_{|{\boldsymbol \ell }-\boldsymbol {L}|,j}} \right]^{-1}. \end{equation}

(A15)

The variance of this estimator is

\begin{equation} {\cal V} = 2 \left[\sum _j \sum _{{\boldsymbol \ell }} \frac{[{\boldsymbol \ell }\cdot \boldsymbol {L}C_{\ell ,j}+\boldsymbol {L}\cdot (\boldsymbol {L}-{\boldsymbol \ell }) C_{|\ell -L|,j}]^2}{\Omega ^2_{\rm s}C^{\rm tot}_{\ell ,j}C^{\rm tot}_{|{\boldsymbol \ell }-\boldsymbol {L}|,j}}\right]^{-1}. \end{equation}

(A16)

Now, we can verify the agreement of this estimator with previous work by going to continuous Fourier space using

\begin{equation} \sum _{{\boldsymbol \ell }} \rightarrow \Omega _{\rm s} \int \frac{\mathrm{d}^2\ell }{(2\pi )^2} \end{equation}

(A17)

that gives

\begin{equation} {\cal V} = \Omega _{\rm s} \left[\sum _j \int \frac{\mathrm{d}^2\ell }{(2\pi )^2} \frac{[{\boldsymbol \ell }\cdot \boldsymbol {L}C_{\ell ,j}+\boldsymbol {L}\cdot (\boldsymbol {L}-{\boldsymbol \ell }) C_{|\ell -L|,j}]^2}{2 C^{\rm tot}_{\ell ,j}C^{\rm tot}_{|{\boldsymbol \ell }-\boldsymbol {L}|,j}}\right]^{-1}, \end{equation}

(A18)

which with the identification Ω_s → (2π)²δ(0) is exactly the Zahn & Zaldarriaga (2005) result.

The behaviour of the lensing estimator and noise in the Gaussian approximation has been analysed in Zahn & Zaldarriaga (2005). Higher values of the discretized parallel wave vector k_∥ (i.e. higher values of j) mean higher values of the three-dimensional power spectrum P(k) for each C_ℓ value. P(k) is monotonically decreasing on all scales of interest and therefore going to higher values of j the effect of C_{ℓ, j} becomes negligible and the signal decays quickly below the thermal noise level. Hence, only a few modes contribute to the estimator (Zahn & Zaldarriaga 2005).

APPENDIX B: LENSING ESTIMATOR FOR UNCLUSTERED DISCRETE SOURCES

Here, we consider a collection of sources with random positions and derive a lensing estimator. The expansion used in the previous section to describe the action of lensing on the surface brightness (equation A7) is not formally valid here because the deflection angle may be large compared to the size of individual sources. For this reason, we derive the estimator by a different method in this case which also gives some insight into how the estimator works.

The fluctuation in the number of sources with luminosity L in a cell i in real space will be

\begin{equation} \delta n^L_i = n^L_i - \bar{n}_i, \end{equation}

(B1)

where the average cell occupation is

\begin{equation} \bar{n}_i = \bar{\eta }\delta V, \end{equation}

(B2)

|$\bar{\eta }$| being the average number density and δV the volume of a cell.

The surface brightness fluctuation δS(i) can be written as

\begin{equation} \delta S(i)=\sum _{\rm L} \delta n^L_i L. \end{equation}

(B3)

The continuous limit in luminosity will be taken later. The DFT of the above expression considering the full 3D case with cells i = (i_⊥, i_∥) (i_⊥ represents the two dimensions perpendicular to the line of sight) is

\begin{equation} \delta \tilde{S}(j) = \frac{\Omega _{\rm s}}{N_\perp N_\parallel }\sum _{i_\parallel } \sum _{i_\perp } \delta S(i){\rm e}^{{\rm i} \frac{2\pi i_\perp j_\perp }{N_\perp }} {\rm e}^{{\rm i} \frac{2\pi i_\parallel j_\parallel }{N_\parallel }} \end{equation}

(B4)

\begin{equation} \hphantom{\delta \tilde{S}(j)}= \frac{\Omega _{\rm s}}{N_\perp N_\parallel } \sum _{i_\parallel } \sum _{i_\perp } \sum _{\rm L} L\delta n^L_i {\rm e}^{{\rm i} \frac{2\pi i_\perp j_\perp }{N_\perp }} {\rm e}^{{\rm i} \frac{2\pi i_\parallel j_\parallel }{N_\parallel }}. \end{equation}

(B5)

For the average flux in a cell |$\bar{S,}$| we will have

\begin{equation} \bar{S} = \bar{n}_i \langle L\rangle = \bar{\eta }\delta V \langle L\rangle . \end{equation}

(B6)

Now let us look at correlations between modes:

\[ \langle \delta \tilde{S}(j)\delta \tilde{S}^{*}(j-m)\rangle = \frac{\Omega _{\rm s}^2}{(N_\perp N_\parallel )^2} \Bigg\langle \sum _{i_\parallel } \sum _{i_\perp } \sum _{\rm L} L\delta n^L_i {\rm e}^{{\rm i} \frac{2\pi i_\perp j_\perp }{N_\perp }} {\rm e}^{{\rm i} \frac{2\pi i_\parallel j_\parallel }{N_\parallel }} \sum _{i_\parallel ^{\prime }} \sum _{i_\perp ^{\prime }} \sum _{L^{\prime }} L^{\prime }\delta n^{L^{\prime }}_{i^{\prime }} {\rm e}^{-{\rm i} \frac{2\pi i_\perp ^{\prime } (j_\perp -m_\perp )}{N_\perp }} {\rm e}^{-{\rm i} \frac{2\pi i_\parallel ^{\prime } (j_\parallel -m_\parallel )}{N_\parallel }}\Bigg\rangle . \]

Since |$\langle \delta n^L_i\rangle =0$|⁠, only the i = i^′, L = L^′ terms contribute:

\begin{equation} \langle \delta \tilde{S}(j)\delta \tilde{S}^{*}(j-m)\rangle = \frac{\Omega _{\rm s}^2}{(N_\perp N_\parallel )^2} \Bigg\langle \sum _{i_\parallel } \sum _{i_\perp } \sum _{\rm L} L^2 (\delta n^L_i)^2 {\rm e}^{{\rm i} \frac{2\pi i_\perp m_\perp }{N_\perp }} {\rm e}^{{\rm i} \frac{2\pi i_\parallel m_\parallel }{N_\parallel }}\Bigg\rangle . \end{equation}

(B7)

The second moment of the number counts in a cell is given by a Poisson distribution |$\langle (\delta n_i)^2 \rangle = \overline{n}_i$|⁠. If a cell of fixed angular size is magnified by a factor μ_i, the galaxies within it will be a factor μ_i brighter. At the same time, the true volume of that cell will be a factor of 1/μ_i smaller and so the average number of galaxies will go down by the same factor. The result is

\begin{equation} \langle \delta \tilde{S}(j)\delta \tilde{S}^{*}(j-m)\rangle = \frac{\Omega _{\rm s}^2}{(N_\perp N_\parallel )^2} \Bigg\langle \sum _{i_\parallel } \sum _{i_\perp } \sum _{\rm L} \mu _{i_\perp } L^2 \bar{n}_i {\rm e}^{{\rm i} \frac{2\pi i_\perp m_\perp }{N_\perp }} {\rm e}^{{\rm i} \frac{2\pi i_\parallel m_\parallel }{N_\parallel }} \Bigg\rangle \end{equation}

(B8)

\begin{equation} \hphantom{\langle \delta \tilde{S}(j)\delta \tilde{S}^{*}(j-m)\rangle}=\frac{\Omega _{\rm s}^2}{(N_\perp N_\parallel )^2} \bar{\eta }\delta V \langle L^2\rangle \Bigg\langle \sum _{i_\perp } \mu _{i_\perp }{\rm e}^{{\rm i} \frac{2\pi i_\perp m_\perp }{N_\perp }} \Bigg\rangle \Bigg\langle \sum _{i_\parallel } {\rm e}^{{\rm i} \frac{2\pi i_\parallel m_\parallel }{N_\parallel }}\Bigg\rangle \end{equation}

(B9)

\begin{equation} \hphantom{\langle \delta \tilde{S}(j)\delta \tilde{S}^{*}(j-m)\rangle}= \frac{\Omega _{\rm s}^2}{N_\perp N_\parallel } \bar{\eta }\delta V \langle L^2\rangle \tilde{\mu }(m_\perp )\delta ^K_{m_\parallel }, \end{equation}

(B10)

where |$\tilde{\mu }(m_\perp )$| is the DFT of the magnification.

Dividing by the average flux |$\bar{S}$| from equation (B6) in order to obtain the dimensionless fluctuations, we find

\begin{eqnarray} \langle \Delta \tilde{S}(j) \Delta \tilde{S}^{*}(j-m_\perp )\rangle &=&\displaystyle\frac{\langle \delta \tilde{S}(j) \delta \tilde{S}^{*}(j-m_\perp )\rangle }{\bar{S}^2} \nonumber\\ &=&\displaystyle\frac{\Omega _{\rm s}^2}{N_\perp N_\parallel } \frac{1}{\bar{\eta }\delta V} \frac{\langle L^2\rangle }{\langle L\rangle ^2} \; \tilde{\mu }(m_\perp ) \end{eqnarray}

(B11)

\begin{equation} \hphantom{\langle \Delta \tilde{S}(j) \Delta \tilde{S}^{*}(j-m_\perp )\rangle} =\frac{\Omega _{\rm s}^2}{\bar{N}_{\rm g}}\frac{\langle L^2\rangle }{\langle L\rangle ^2} \tilde{\mu }(m_\perp ) \end{equation}

(B12)

\begin{equation} \hphantom{\langle \Delta \tilde{S}(j) \Delta \tilde{S}^{*}(j-m_\perp )\rangle}=\Omega _{\rm s} C^{\rm shot} \tilde{\mu }(m_\perp ), \end{equation}

(B13)

where we have used |$\delta V = \frac{V_{{\rm tot}}}{N_\perp N_\parallel }$| and |$\bar{N}_{\rm g}=\bar{\eta } V_{{\rm tot}}$|⁠. Equation (B13) defines C^shot. The moments of the luminosity function can now be calculated for a continuous distribution of luminosities.

We construct the quadratic estimator |$\hat{\mu }(m_\perp )$| as

\begin{equation} \hat{\mu }(m_\perp )= \frac{1}{\Omega _{\rm s} C^{\rm shot} }\frac{1}{N_\perp N_\parallel } \sum _{j_\perp } \sum _{j_\parallel } \Delta \tilde{S}(j) \Delta \tilde{S}^{*}(j-m_\perp ). \end{equation}

(B14)

In this case (with no thermal noise or clustering), the optimal filter (⁠|$g({\boldsymbol \ell },{\boldsymbol L},j)$| in the previous section) is a function of only |${\boldsymbol m}_\perp$| and its value is given by the requirement that |$\langle \hat{\mu }(m_\perp ) \rangle = \tilde{\mu }(m_\perp )$|⁠. The variance of this estimator is

\begin{eqnarray} &&{{\cal V} = \langle |\hat{\mu (m_\perp )}|^2\rangle = \langle \hat{\mu }(m_\perp )\hat{\mu }^{*}(m_\perp )\rangle} \nonumber\\ &&{\hphantom{{\cal V}}= \displaystyle\frac{1}{(N_\perp N_\parallel )^2} \frac{1}{(\Omega _{\rm s} C^{\rm shot})^2 }\sum _{j_\perp } \sum _{j_\parallel } \sum _{j_\perp ^{\prime }} \sum _{j_\parallel ^{\prime }} \langle \Delta \tilde{S}(j)\Delta \tilde{S}^{*}(j-m_\perp ) \Delta \tilde{S}^{*}(j^{\prime })\Delta \tilde{S}(j^{\prime }-m_\perp ) \rangle} \end{eqnarray}

(B15)

\begin{equation} \hphantom{{\cal V}}=\frac{\Omega _{\rm s}^2}{\bar{N}_{\rm g}}\frac{\langle L^4\rangle }{\langle L^2\rangle ^2}\left(1+3\frac{\bar{N}_{\rm g}}{N_\perp N_\parallel }\right) +2\Omega _{\rm s}^2 \frac{N_\perp N_\parallel -1}{(N_\perp N_\parallel )^2}. \end{equation}

(B16)

In deriving equation (B16), we have used the higher moments of the Poisson number counts in cells. Note that in the limit N_⊥, N_∥ → ∞ only the first term survives. In the notation of the previous section |$\bar{N}_{\rm g} = \overline{\eta } \Omega _{\rm s} {\cal D}^2 {\mathcal {L}}$|⁠, where |$\overline{\eta }$| is the density of galaxies. Note that for 21 cm emission the luminosity is proportional to the H i density so the luminosity moments are also the moments of the H i mass function.

In the weak lensing limit μ ≃ 1 + 2κ = 1 − ∇²Ψ. So the estimator for the lensing potential is related to the estimator for magnification by

\begin{equation} \hat{\Psi }({\boldsymbol L}) = |{\boldsymbol L}|^{-2} \hat{\mu }({\boldsymbol L}). \end{equation}

(B17)

As it turns out, this filter is the same one we would have found by plugging a constant power spectrum into the optimal Gaussian filter (A15), but the noise is different.

APPENDIX C: LENSING ESTIMATOR FOR CLUSTERED DISCRETE SOURCES

In this appendix, the results of the two previous appendices are combined to find a lensing estimator and its expected noise in the case of clustered point sources. For the two-point correlation function, we have

\begin{equation} \langle I(\boldsymbol x) I(\boldsymbol x^{\prime }) \rangle = \frac{1}{\bar{\eta }\delta V}\frac{\langle L^2\rangle }{\langle L\rangle ^2}\delta ^K_{\boldsymbol x\boldsymbol x^{\prime }} + \xi _{\boldsymbol x \boldsymbol x^{\prime }}, \end{equation}

(C1)

where the first part comes from the Poisson fluctuations in the number counts and the second from the clustering. In Fourier space

\begin{equation} \langle I_{\boldsymbol k} I^*_{\boldsymbol k^{\prime }}\rangle = \Omega _{\rm s} \; (C_{\ell ,j}+C^{\rm shot}) \; \delta ^{K}_{\ell ,\ell ^{\prime }}\delta ^{K}_{jj^{\prime }}. \end{equation}

(C2)

After lensing the observed correlations will be

\[ \langle \tilde{I}_{{\boldsymbol \ell },j}\tilde{I}^*_{{\boldsymbol \ell }-\boldsymbol {L},j^{\prime }}\rangle = \delta ^K_{jj^{\prime }} [{\boldsymbol \ell }\cdot \boldsymbol {L}C_{\ell ,j}+\boldsymbol {L}\cdot (\boldsymbol {L}-{\boldsymbol \ell }) C_{|\ell -L|,j}+L^2 \; C^{\rm shot} ]\; \Psi (\boldsymbol {L}) \]

to first order. We will not attempt to find an optimal estimator here. Instead, we use an estimator of the form

\begin{equation} \hat{\Psi }(\boldsymbol {L})= f(\boldsymbol {L}) \sum _j \sum _{{\boldsymbol \ell }} \tilde{I}_{{\boldsymbol \ell },j}\tilde{I}^*_{{\boldsymbol \ell }-\boldsymbol {L},j}, \end{equation}

(C3)

where |$f(\boldsymbol {L})$| is a normalization. This would be optimal were the sources unclustered. For an unbiased estimator, we require that |$\langle \hat{\Psi }(\boldsymbol {L})\rangle = \Psi (\boldsymbol {L})$| so we find

\begin{eqnarray} &&{f(\boldsymbol {L})=\bigg \lbrace \sum _j \sum _{{\boldsymbol \ell }}[{\boldsymbol \ell }\cdot \boldsymbol {L}C_{\ell ,j}+\boldsymbol {L}\cdot (\boldsymbol {L}-{\boldsymbol \ell }) C_{|\ell -L|,j}+L^2 \; C^{\rm shot}]\bigg \rbrace ^{-1}} \nonumber\\ &&{\hphantom{f(\boldsymbol {L})}=\,\bigg \lbrace \sum _j \sum _{{\boldsymbol \ell }}[{\boldsymbol \ell }\cdot \boldsymbol {L}C_{\ell ,j}+\boldsymbol {L}\cdot (\boldsymbol {L}-{\boldsymbol \ell }) C_{|\ell -L|,j}]+(N_\perp N_\parallel )L^2 \; C^{\rm shot}\bigg \rbrace ^{-1}.} \end{eqnarray}

(C4)

The variance of the estimator |${\cal V}=\langle \hat{\Psi }(\boldsymbol {L})\hat{\Psi }^*(\boldsymbol {L})\rangle$| is given by

\begin{equation} {\cal V}=[f(\boldsymbol {L})]^2 \sum _j \sum _{j^{\prime }} \sum _{{\boldsymbol \ell }} \sum _{{\boldsymbol \ell }^{\prime }} \langle I_{{\boldsymbol \ell },j}I^*_{{\boldsymbol \ell }-\boldsymbol {L},j}I^*_{{\boldsymbol \ell }^{\prime },j^{\prime }}I_{{\boldsymbol \ell }^{\prime }-\boldsymbol {L},j^{\prime }}\rangle , \end{equation}

(C5)

with

\begin{eqnarray} &&{\langle I_{{\boldsymbol \ell },j}I^*_{{\boldsymbol \ell }-\boldsymbol {L},j}I^*_{{\boldsymbol \ell }^{\prime },j^{\prime }}I_{{\boldsymbol \ell }^{\prime }-\boldsymbol {L},j^{\prime }}\rangle = \Omega _{\rm s} \displaystyle\frac{1}{\bar{\eta }^3}\frac{1}{(D^2{\cal L})^3}\frac{\langle L^4\rangle }{\langle L\rangle ^4}\left(1+3\frac{\bar{N}_{\rm g}}{N_\perp N_\parallel }\right)} \nonumber\\ &&{\quad+\, \Omega _{\rm s}\displaystyle\frac{1}{\bar{\eta }^2}\frac{1}{(D^2{\cal L})^2}\frac{\langle L^3\rangle }{\langle L\rangle ^3}\left(\frac{N_\perp N_\parallel -1}{N_\perp N_\parallel }\right) [C_{|\ell ^{\prime }-L|,j^{\prime }}+C_{\ell ^{\prime },j^{\prime }}+C_{|\ell -L|,j}+C_{\ell ,j}]} \nonumber\\ &&{\quad+\, \Omega _{\rm s}\displaystyle\frac{1}{\bar{\eta }^2}\frac{1}{(D^2{\cal L})^2}\frac{\langle L^2\rangle ^2}{\langle L\rangle ^4}\left(\frac{N_\perp N_\parallel -1}{N_\perp N_\parallel }\right) [C_{{\rm L},0}+C_{|\ell -\ell ^{\prime }|,|j-j^{\prime }|}+C_{|\ell +\ell ^{\prime }-L|,j+j^{\prime }}]} \nonumber\\ &&{\quad+\, \Omega ^2_{\rm s} \displaystyle\frac{1}{\bar{\eta }}\frac{1}{(D^2{\cal L})}\frac{\langle L^2\rangle }{\langle L\rangle ^2}\frac{(N_\perp N_\parallel -1)(N_\perp N_\parallel -2)}{(N_\perp N_\parallel )^2}} \nonumber\\ &&{\quad \times\, [C_{\ell ^{\prime },j^{\prime }}\delta ^K_{\ell ,\ell ^{\prime }}\delta ^K_{j,j^{\prime }}\delta ^K_{\ell ,\ell ^{\prime }}\delta ^K_{j,j^{\prime }}+C_{{\rm L}-\ell ^{\prime },-j^{\prime }} \delta ^K_{\ell ,L-\ell ^{\prime }}\delta ^K_{j,-j^{\prime }}\delta ^K_{\ell ,L-\ell ^{\prime }}\delta ^K_{j,-j^{\prime }}} \nonumber\\ &&{\quad+\, C_{\ell ^{\prime },j^{\prime }}\delta ^K_{L-\ell ,\ell ^{\prime }}\delta ^K_{j^{\prime },-j}\delta ^K_{L-\ell ,\ell ^{\prime }}\delta ^K_{j^{\prime },-j} +C_{\ell ^{\prime }-L,j^{\prime }}\delta ^K_{\ell ,\ell ^{\prime }}\delta ^K_{j,j^{\prime }}\delta ^K_{\ell ,\ell ^{\prime }}\delta ^K_{j,j^{\prime }}]} \nonumber\\ &&{\quad+\, \Omega ^2_{\rm s} \displaystyle\frac{(N_\perp N_\parallel -1)(N_\perp N_\parallel -2)(N_\perp N_\parallel -3)}{(N_\perp N_\parallel )^3}} \nonumber\\ &&{\quad\times\, [C_{\ell ,j}C_{|\ell ^{\prime }-L|,j^{\prime }}\delta ^K_{\ell ,\ell ^{\prime }}\delta ^K_{j,j^{\prime }}\delta ^K_{\ell ,\ell ^{\prime }}\delta ^K_{j,j^{\prime }} +C_{\ell ,j}C_{\ell ^{\prime },j^{\prime }}\delta ^K_{\ell ,L-\ell ^{\prime }}\delta ^K_{j,-j^{\prime }}\delta ^K_{\ell ,L-\ell ^{\prime }}\delta ^K_{j,-j^{\prime }}]. } \end{eqnarray}

(C6)

Here, we have assumed that the clustering and the noise are Gaussian.

We will now calculate the sum over |${\boldsymbol \ell },{\boldsymbol \ell }^{\prime },j,j^{\prime }$| in equation (C5) separately for each of the terms of the fourth moment (C6). We will also consider the limit N_∥ → ∞ to simplify the constant factors. Note that, to include the thermal noise of the array in our calculation, we just send |$C_{\ell ,j} \rightarrow C^{\rm tot}_{\ell ,j}$|⁠, with

\begin{equation} C^{\rm tot}_{\ell ,j} = C_{\ell ,j} + C^N \end{equation}

(C7)

with C^N the power spectrum of the thermal noise. We have

\begin{equation} {\cal V} = [f(\boldsymbol {L})]^2 \left({\cal I}_0+{\cal I}_1+{\cal I}_2+{\cal I}_3+{\cal I}_4\right), \end{equation}

(C8)

with

\begin{eqnarray} {\cal I}_0 &=& \sum _j \sum _{j^{\prime }} \sum _{{\boldsymbol \ell }} \sum _{{\boldsymbol \ell }^{\prime }} \Omega _{\rm s} \frac{1}{\bar{\eta }^3}\displaystyle\frac{1}{(D^2{\cal L})^3}\frac{\langle L^4\rangle }{\langle L\rangle ^4} \nonumber\\ &=& (N_\perp N_\parallel )^2 \; \Omega _{\rm s} \displaystyle\frac{1}{\bar{\eta }^3}\frac{1}{(D^2{\cal L})^3}\frac{\langle L^4\rangle }{\langle L\rangle ^4}. \end{eqnarray}

(C9)

\begin{eqnarray} {\cal I}_1 &=& \sum _j \sum _{j^{\prime }} \sum _{{\boldsymbol \ell }} \sum _{{\boldsymbol \ell }^{\prime }}\Omega _{\rm s}\frac{1}{\bar{\eta }^2}\displaystyle\frac{1}{(D^2{\cal L})^2}\frac{\langle L^3\rangle }{\langle L\rangle ^3} [C^{\rm tot}_{|\ell ^{\prime }-L|,j^{\prime }}+C^{\rm tot}_{\ell ^{\prime },j^{\prime }} + C^{\rm tot}_{|\ell -L|,j}+C^{\rm tot}_{\ell ,j}]\nonumber\\ &=&\Omega _{\rm s}\displaystyle\frac{1}{\bar{\eta }^2}\frac{1}{(D^2{\cal L})^2}\frac{\langle L^3\rangle }{\langle L\rangle ^3} \times (N_\perp N_\parallel )\sum _j\sum _{{\boldsymbol \ell }}[2C^{\rm tot}_{\ell ,j}+2C^{\rm tot}_{|\ell -L|,j}]. \end{eqnarray}

(C10)

\begin{eqnarray} {\cal I}_2 &=& \sum _j \sum _{j^{\prime }} \sum _{{\boldsymbol \ell }} \sum _{{\boldsymbol \ell }^{\prime }} \Omega _{\rm s}\frac{1}{\bar{\eta }^2}\frac{1}{(D^2{\cal L})^2}\frac{\langle L^2\rangle ^2}{\langle L\rangle ^4} [C^{\rm tot}_{L,0}+C^{\rm tot}_{|\ell -\ell ^{\prime }|,|j-j^{\prime }|} +\, C^{\rm tot}_{|\ell +\ell ^{\prime }-L|,j+j^{\prime }}] \nonumber\\ &=& \Omega _{\rm s}\frac{1}{\bar{\eta }^2}\frac{1}{(D^2{\cal L})^2}\frac{\langle L^2\rangle ^2}{\langle L\rangle ^4} (N_\perp N_\parallel )^2 C^{\rm tot}_{L,0} \nonumber\\ &&+\, \Omega _{\rm s}\frac{1}{\bar{\eta }^2}\frac{1}{(D^2{\cal L})^2}\frac{\langle L^2\rangle ^2}{\langle L\rangle ^4} \sum _j \sum _{j^{\prime }} \sum _{{\boldsymbol \ell }} \sum _{{\boldsymbol \ell }^{\prime }} [C^{\rm tot}_{|\ell -\ell ^{\prime }|,|j-j^{\prime }|}+C^{\rm tot}_{|\ell +\ell ^{\prime }-L|,j+j^{\prime }}]. \end{eqnarray}

(C11)

\begin{eqnarray} \nonumber\\ &&{{\cal I}_3 = \sum _j \sum _{j^{\prime }} \sum _{{\boldsymbol \ell }} \sum _{{\boldsymbol \ell }^{\prime }} \Omega ^2_s \frac{1}{\bar{\eta }}\frac{1}{(D^2{\cal L})}\frac{\langle L^2\rangle }{\langle L\rangle ^2}} \nonumber\\ &&{\qquad \times\, [C^{\rm tot}_{\ell ^{\prime },j^{\prime }}\delta ^K_{\ell ,\ell ^{\prime }}\delta ^K_{j,j^{\prime }}\delta ^K_{\ell ,\ell ^{\prime }}\delta ^K_{j,j^{\prime }}+C^{\rm tot}_{L-\ell ^{\prime },-j^{\prime }} \delta ^K_{\ell ,L-\ell ^{\prime }}\delta ^K_{j,-j^{\prime }}\delta ^K_{\ell ,L-\ell ^{\prime }}\delta ^K_{j,-j^{\prime }}} \nonumber\\ &&{\qquad+\,C^{\rm tot}_{\ell ^{\prime },j^{\prime }}\delta ^K_{L-\ell ,\ell ^{\prime }}\delta ^K_{j^{\prime },-j}\delta ^K_{L-\ell ,\ell ^{\prime }}\delta ^K_{j^{\prime },-j} +C^{\rm tot}_{\ell ^{\prime }-L,j^{\prime }}\delta ^K_{\ell ,\ell ^{\prime }}\delta ^K_{j,j^{\prime }}\delta ^K_{\ell ,\ell ^{\prime }}\delta ^K_{j,j^{\prime }}]} \nonumber\\ &&{\quad=\, \Omega ^2_{\rm s} \frac{1}{\bar{\eta }}\frac{1}{(D^2{\cal L})}\frac{\langle L^2\rangle }{\langle L\rangle ^2} \sum _j \sum _{{\boldsymbol \ell }} [2C^{\rm tot}_{\ell ,j}+2C^{\rm tot}_{|\ell -L|,j}].} \end{eqnarray}

(C12)

\begin{eqnarray} &&{{\cal I}_4 =\sum _j \sum _{j^{\prime }} \sum _{{\boldsymbol \ell }} \sum _{{\boldsymbol \ell }^{\prime }} \Omega ^2_{\rm s} [C^{\rm tot}_{\ell ,j}C^{\rm tot}_{|\ell ^{\prime }-L|,j^{\prime }}\delta ^K_{\ell ,\ell ^{\prime }}\delta ^K_{j,j^{\prime }}\delta ^K_{\ell ,\ell ^{\prime }}\delta ^K_{j,j^{\prime }}} \nonumber\\ &&{\qquad +\,C^{\rm tot}_{\ell ,j}C^{\rm tot}_{\ell ^{\prime },j^{\prime }}\delta ^K_{\ell ,L-\ell ^{\prime }}\delta ^K_{j,-j^{\prime }}\delta ^K_{\ell ,L-\ell ^{\prime }}\delta ^K_{j,-j^{\prime }}]} \nonumber\\ &&{\quad=\,\Omega ^2_{\rm s} \sum _j \sum _{{\boldsymbol \ell }} 2C^{\rm tot}_{\ell ,j}C^{\rm tot}_{|\ell -L|,j}. } \end{eqnarray}

(C13)

From appendices A and B, we can recognize equation (C9) as a pure Poisson clustering term and equation (C13) as a pure Gaussian clustering term. The other terms are from the interplay of these contributions.

To calculate these terms, it is convenient to go into continuous ℓ-space where the integrals can be done numerically. We have

\begin{eqnarray} &&{f(\boldsymbol {L})=\bigg \lbrace \sum _j \sum _{{\boldsymbol \ell }}[{\boldsymbol \ell }\cdot \boldsymbol {L}C_{\ell ,j}+\boldsymbol {L}\cdot (\boldsymbol {L}-{\boldsymbol \ell }) C_{|\ell -L|,j}+L^2 \; C^{\rm shot}]\bigg \rbrace ^{-1}} \nonumber\\ &&{\hphantom{f(\boldsymbol {L})}=\, \frac{1}{\Omega _{\rm s}}\bigg \lbrace \sum _j \int \frac{\mathrm{d}^2\ell }{(2\pi )^2}[{\boldsymbol \ell }\cdot \boldsymbol {L}C_{\ell ,j}+\boldsymbol {L}\cdot (\boldsymbol {L}-{\boldsymbol \ell }) C_{|\ell -L|,j}+L^2 \; C^{\rm shot}]\bigg \rbrace ^{-1}.} \end{eqnarray}

(C14)

\begin{eqnarray} &&{{\cal I}_0 = \sum _j \sum _{j^{\prime }} \sum _{{\boldsymbol \ell }} \sum _{{\boldsymbol \ell }^{\prime }} \Omega _{\rm s} \displaystyle\frac{1}{\bar{\eta }^3}\displaystyle\frac{1}{(D^2{\cal L})^3}\displaystyle\frac{\langle L^4\rangle }{\langle L\rangle ^4}} \nonumber\\ &&{\quad=\, (N_\parallel )^2\Omega ^3_{\rm s} \displaystyle\frac{1}{\bar{\eta }^3}\displaystyle\frac{1}{(D^2{\cal L})^3}\displaystyle\frac{\langle L^4\rangle }{\langle L\rangle ^4}\left(\int \displaystyle\frac{\mathrm{d}^2\ell }{(2\pi )^2}\right)^2,} \end{eqnarray}

(C15)

where

\begin{equation} \int \mathrm{d}^2\ell = \pi (\ell ^2_{\rm max}-\ell ^2_{\rm min}). \end{equation}

(C16)

\begin{eqnarray} &&{{\cal I}_1 =\Omega _{\rm s}\displaystyle\frac{1}{\bar{\eta }^2}\displaystyle\frac{1}{(D^2{\cal L})^2}\displaystyle\frac{\langle L^3\rangle }{\langle L\rangle ^3} \times (N_\perp N_\parallel )\sum _j\sum _{{\boldsymbol \ell }}[2C^{\rm tot}_{\ell ,j}+2C^{\rm tot}_{|\ell -L|,j}]} \nonumber\\ &&{\quad=\,\Omega ^3_{\rm s}\displaystyle\frac{1}{\bar{\eta }^2}\displaystyle\frac{1}{(D^2{\cal L})^2}\displaystyle\frac{\langle L^3\rangle }{\langle L\rangle ^3} N_\parallel \left(\int \displaystyle\frac{d^2\ell ^{\prime }}{(2\pi )^2}\right) \sum _j \int \displaystyle\frac{\mathrm{d}^2\ell }{(2\pi )^2} [2C^{\rm tot}_{\ell ,j}+2C^{\rm tot}_{|\ell -L|,j}].} \end{eqnarray}

(C17)

\begin{eqnarray} &&{{\cal I}_2 = \sum _j \sum _{j^{\prime }} \sum _{{\boldsymbol \ell }} \sum _{{\boldsymbol \ell }^{\prime }} \Omega _{\rm s}\displaystyle\frac{1}{\bar{\eta }^2}\displaystyle\frac{1}{(D^2{\cal L})^2}\displaystyle\frac{\langle L^2\rangle ^2}{\langle L\rangle ^4} [C^{\rm tot}_{L,0}+C^{\rm tot}_{|\ell -\ell ^{\prime }|,|j-j^{\prime }|}+ C^{\rm tot}_{|\ell +\ell ^{\prime }-L|,j+j^{\prime }}]} \nonumber\\ &&{\hphantom{{\cal I}_2}=\,\Omega ^3_{\rm s}\displaystyle\frac{1}{\bar{\eta }^2}\displaystyle\frac{1}{(D^2{\cal L})^2}\displaystyle\frac{\langle L^2\rangle ^2}{\langle L\rangle ^4} (N_\parallel )^2 \left(\int \displaystyle\frac{\mathrm{d}^2\ell }{(2\pi )^2}\right)^2 C^{\rm tot}_{L,0} +\Omega ^3_{\rm s}\displaystyle\frac{1}{\bar{\eta }^2}\displaystyle\frac{1}{(D^2{\cal L})^2}\displaystyle\frac{\langle L^2\rangle ^2}{\langle L\rangle ^4} \sum _j \sum _{j^{\prime }} \int \displaystyle\frac{\mathrm{d}^2\ell }{(2\pi )^2} \int \displaystyle\frac{\mathrm{d}^2\ell ^{\prime }}{(2\pi )^2} \times [C^{\rm tot}_{|\ell -\ell ^{\prime }|,|j-j^{\prime }|}+C^{\rm tot}_{|\ell +\ell ^{\prime }-L|,j+j^{\prime }}].} \nonumber \\ \end{eqnarray}

(C18)

We are not able to probe the k_∥ = 0 mode (its signal will be removed with foreground cleaning), so we ignore the C_{L, 0} contribution to the |${\cal I}_2$| part of the noise and all our plots are with j_min = 1.

\begin{eqnarray} &&{{\cal I}_3 = \Omega ^2_{\rm s} \displaystyle\frac{1}{\bar{\eta }}\displaystyle\frac{1}{(D^2{\cal L})}\displaystyle\frac{\langle L^2\rangle }{\langle L\rangle ^2} \sum _j \sum _{{\boldsymbol \ell }} [2C^{\rm tot}_{\ell ,j}+2C^{\rm tot}_{|\ell -L|,j}]} \nonumber\\ &&{\quad =\,\Omega ^3_{\rm s} \displaystyle\frac{1}{\bar{\eta }}\displaystyle\frac{1}{(D^2{\cal L})}\displaystyle\frac{\langle L^2\rangle }{\langle L\rangle ^2} \sum _j \int \displaystyle\frac{\mathrm{d}^2\ell }{(2\pi )^2} [2C^{\rm tot}_{\ell ,j}+2C^{\rm tot}_{|\ell -L|,j}].} \end{eqnarray}

(C19)

\begin{eqnarray} \nonumber \\ &&{{\cal I}_4 =\Omega ^2_{\rm s} \sum _j \sum _{{\boldsymbol \ell }} 2C^{\rm tot}_{\ell ,j}C^{\rm tot}_{|\ell -L|,j}} \nonumber\\ &&{\quad =\, \Omega ^3_{\rm s} \sum _j \int \displaystyle\frac{\mathrm{d}^2\ell }{(2\pi )^2} 2C^{\rm tot}_{\ell ,j}C^{\rm tot}_{|\ell -L|,j}.} \end{eqnarray}

(C20)

Here, we should note that foreground subtraction techniques will remove the first few k_∥ modes (Zahn & Zaldarriaga 2005; McQuinn et al. 2006), meaning that we will need to use some j_min > 1 for the lensing reconstruction. Removing only the first few modes would correspond to a relatively simple foreground contamination, while a large number would denote a much more complex problem. Here, we should note that in our case the S/N is not greatly affected by removing the few first modes, and we plan to investigate foreground contamination and subtraction techniques in future work.

The behaviour of the lensing estimator and noise for clustered discrete sources is more complicated than in the Gaussian approximation due to the complexity of the various contributing terms. As we have stressed in the main text, in our method the Poisson fluctuations contribute both to the signal and the noise of the estimator. More specifically, the contribution of the second moment C^shot (see equation C14) is crucial for obtaining a low lensing reconstruction noise level. Of the noise terms the dominant contribution comes from the fourth moment Poisson term |${\cal I}_0$|⁠, while the |${\cal I}_2$| term has the smallest contribution. It is also useful to note that, in a similar manner with the Gaussian estimator, the lensing reconstruction noise converges with increasing j, so that a relatively small number of modes contributes to the final estimator. For example, our calculations for the source redshift z = 2 allow us to use a maximum number of parallel modes j_max = 63, but the noise has already converged at j ∼ 40. To illustrate that, we show a plot demonstrating the convergence of the lensing reconstruction noise N(L) with j using the SKA2 specifications and Model B for the H i mass function. We see that for this case the noise converges fast (this naturally depends on the interferometer specifications, which determine the contribution of the instrument's thermal noise, as well as the evolution of the H i mass function which affects the signal and the Poisson moments). Also note that the shape of L²N(L) is nearly flat up to a scale where noise becomes important – this behaviour is qualitatively similar to the one of the Gaussian case (Zahn & Zaldarriaga 2005).