Full-Sky Lensing Reconstruction of 21 cm Intensity Maps

TL;DR

This paper develops a full-sky formalism for reconstructing matter density perturbations from 21 cm intensity maps, highlighting the importance of accounting for sky curvature in future large-scale cosmological analyses.

Contribution

It extends the lensing reconstruction formalism to full-sky surveys using Spherical Fourier-Bessel expansion, providing predictions and comparisons with flat-sky approaches.

Findings

Agreement with flat-sky predictions at small scales

Significant deviations at large scales due to curvature

SNR for lensing signal increases by 107% with the new formalism

Abstract

Weak gravitational lensing of the 21 cm radiation is expected to be an important cosmological probe for post-reionization physics. We investigate the reconstruction of the matter density perturbations using a quadratic minimum variance estimator. The next generation of line intensity mapping (LIM) surveys such as HIRAX and CHIME will cover a larger sky fraction, which requires one to account for the curvature in the sky. Thus, we extend the plane-parallel flat-sky formalism for lensing reconstruction to account for a full-sky survey using the Spherical Fourier-Bessel (SFB) expansion. Using the HIRAX 21 cm survey as a basis, we make predictions for lensing-reconstruction noise in our formalism and compare our results with the predictions from the plane-parallel formalism. We find agreement with the plane-parallel noise power spectrum at small scales and a significant deviation at scales $L\lesssim \ell_{\rm res}-k_{\rm eq}R$ where $R$ is the radius of the shell volume, $k_{\rm eq}$ is the wavenumber for matter-radiation equality, and $\ell_{\rm res}$ is the angular resolution scale. Furthermore, we derive the SFB flat-sky reconstruction noise and compare it with the full-sky SFB case as well as the plane-parallel case, finding minor deviations from the full-sky noise due to sphericity. We also determine that, in the absence of non-Gaussian statistics of the intensity field but accounting for foregrounds, the signal-to-noise ratio (SNR) for $C_\ell^{\phi\phi}$ using our SFB estimator increases by 107%. This shows that accounting for the curved sky in LIM weak lensing will be crucial for large-scale cosmology.

Tables

Table 1. Table 1: HI clustering bias values used for our forecasts based on predictions from Castorina & Villaescusa-Navarro ( 2017 ) .

z	$b_{\mathrm{HI}}$
$1.38\le z\le 1.48$	1.6
$1.48\le z\le 1.59$	1.7
$1.59\le z\le 1.71$	1.8
$1.71\le z\le 1.84$	1.9
$1.84\le z\le 1.99$	2.0
$1.99\le z\le 2.16$	2.0
$2.16\le z\le 2.35$	2.2
$2.35\le z\le 2.57$	2.4

Table 2. Table 2: Signal-to-noise ratio (SNR) of C ℓ ϕ ​ ϕ superscript subscript 𝐶 ℓ italic-ϕ italic-ϕ C_{\ell}^{\phi\phi} for a HIRAX 21-cm intensity map at redshift 2 over a redshift bin size Δ ​ z = 0.16 Δ 𝑧 0.16 \Delta z=0.16 over redshift range 1.38 < z < 2.57 1.38 𝑧 2.57 1.38<z<2.57 and LSST galaxies in the redshift range 0 < z < 5 0 𝑧 5 0<z<5 . Note that these SNR values are overestimated in that they do not account for bias-hardening to remove the effect of nonlinear clustering in the source maps or the nonlinear clustering in the fiducial power spectrum, though we do limit the nonlinear scales used to k max = 0.5 ​ h subscript 𝑘 max 0.5 ℎ k_{\rm max}=0.5h /Mpc. We also show the increase in the SNR when replacing the plane-parallel estimator with the SFB full-sky estimator, finding that the SNR changes significantly.

z	Plane-Parallel	SFB	% Increase
$1.38\le z\le 1.48$	1.22	3.9	219
$1.48\le z\le 1.59$	1.27	3.2	151
$1.59\le z\le 1.71$	1.29	2.7	109
$1.71\le z\le 1.84$	1.25	2.1	68
$1.84\le z\le 1.99$	1.25	2.0	60
$1.99\le z\le 2.16$	1.02	1.5	47
$2.16\le z\le 2.35$	0.88	1.3	47
$2.35\le z\le 2.57$	0.96	1.1	14
$1.38\le z\le 2.57$	3.51	7.3	107

Table 3. Table 3: Signal-to-noise ratio (SNR) of C ℓ κ ​ g superscript subscript 𝐶 ℓ 𝜅 𝑔 C_{\ell}^{\kappa g} for a HIRAX 21-cm intensity map over redshift range 1.38 < z < 2.57 1.38 𝑧 2.57 1.38<z<2.57 and LSST galaxies in the redshift range 0 < z < 5 0 𝑧 5 0<z<5 . Note that the increase in the SNR for the SFB case is low overall, although the increase can be as high as 33% for the lowest redshifts.

z	Plane-Parallel	SFB	% Increase
$1.38\le z\le 1.48$	23	31	33
$1.48\le z\le 1.59$	24	27	12.5
$1.59\le z\le 1.71$	25	26	3.8
$1.71\le z\le 1.84$	25	26	3.8
$1.84\le z\le 1.99$	26	26	3.8
$1.99\le z\le 2.16$	24	25	4.16
$2.16\le z\le 2.35$	23	25	8.6
$2.35\le z\le 2.57$	25.1	25.3	0.7
$1.38\le z\le 2.57$	69	76	10.1

Equations

$$\tilde{T}(\hat{\text{r}},r)=T(\hat{\text{r}}+\nabla \varphi ,r).$$

$$T(\mathbf{k})$$

$$T(\mathbf{r})$$

$$⟨T(\mathbf{k})T^{\ast }({\mathbf{k}}^{\prime })⟩=(2\pi {)}^3{\delta }_D^3(\mathbf{k}-{\mathbf{k}}^{\prime })P_T(k)$$

$$\varphi (\hat{\mathbf{r}})=\frac{2}{c^2}{\int }_0^{{\chi }_s}d\chi \frac{{\chi }_s-\chi }{{\chi }_s\chi }\Phi ({\chi }_s\hat{\mathbf{r}},z),$$

$${\nabla }^2\Phi =\frac{3}{2}{\Omega }_{m,0}{H}_0^2(1+z)\delta ,$$

$$C_{\ell }^{\varphi \varphi }=\frac{9}{c^4}{\Omega }_{m,0}^2{H}_0^4\frac{2}{\pi }\int \frac{dk}{k^2}[W_{\ell }^{\varphi }(k){]}^{2}P(k),$$

$$W_{\ell }^{\varphi }(k)={\int }_0^{{\chi }_s}d\chi \frac{{\chi }_s-\chi }{{\chi }_s\chi }(1+z)D(z)j_{\ell }(k\chi ),$$

$$N_L^{\varphi \varphi }={\left[\sum _{j\ge j_{\min }}\int \frac{d^2\ell }{(2\pi {)}^2}\frac{f_{\varphi }^2(\mathbf{L},\mathbf{\ell },k_{\parallel ,j})}{C_{\ell }^{\mathrm{tot}}(k_{\parallel ,j})C_{|\mathbf{L}-\mathbf{\ell }|}^{\mathrm{tot}}(k_{\parallel ,j})}\right]}^{-1},$$

$$f_{\varphi }(\mathbf{L},\mathbf{\ell },k_{\parallel })$$

$$C_{\ell }(k_{\parallel })$$

$$C_{\ell }^{\mathrm{tot}}(k_{\parallel })$$

$$r(z)={\int }_0^z\frac{d{z}^{\prime }}{H(z^{\prime })},$$

$$r=r_{\mathrm{ref}}-\frac{c(1+z_{\mathrm{ref}}{)}^2}{{\nu }_{\mathrm{rest}}H(z_{\mathrm{ref}})}(\nu -{\nu }_{\mathrm{ref}}).$$

$$T_{\ell m}(k)=\sqrt{\frac{2}{\pi }}\int d^{3}r{j}_{\ell }(kr)Y_{\ell m}^{\ast }(\hat{\mathbf{r}})T(\mathbf{r}),$$

$$T(\mathbf{r})=\sqrt{\frac{2}{\pi }}\sum _{\ell m}\int dk{k}^2{j}_{\ell }(kr)Y_{\ell m}(\hat{\mathbf{r}})T_{\ell m}(k).$$

$$\int dr{r}^2{j}_{\ell }(kr)j_{\ell }(k^{\prime }r)=\frac{\pi }{2}\frac{{\delta }_D(k-k^{\prime })}{k^2}.$$

$$⟨T_{\ell m}(k)T_{{\ell }^{\prime }{m}^{\prime }}^{\ast }(k^{\prime })⟩=P_T(k)\frac{{\delta }_D(k-k^{\prime })}{k^2}{\delta }_{\ell {\ell }^{\prime }}{\delta }_{m{m}^{\prime }},$$

$$⟨T_{\ell m}(k)T_{{\ell }^{\prime }{m}^{\prime }}^{\ast }(k^{\prime })⟩=C_{\ell }(k,k^{\prime }){\delta }_{\ell {\ell }^{\prime }}{\delta }_{m{m}^{\prime }}.$$

$$C_{\ell }(k_1,k_2)$$

$$W_{\ell }^0(k,k^{\prime })$$

$$W_{\ell }^r(k,k^{\prime })$$

$$j_{\ell }^{\prime \prime }(x)$$

$$S_{\ell }(k)\equiv 2{\pi }^2{\left[\int d^{3}r{\varphi }^2(r)j_{\ell }^2(kr)\right]}^{-1}\frac{{\sum }_m|T_{\ell m}^{\mathrm{meas}}(k){|}^2}{2\ell +1},$$

$$T_{\ell mn}={\int }_{r_{\min }\le r\le r_{\max }}{d}^{3}r{j}_{\ell }(k_{\ell n}r)Y_{\ell m}^{\ast }(\hat{\mathbf{r}})T(\mathbf{r}),$$

$$s_{\ell n}\simeq \ell +n\pi \frac{r_{\max }}{\Delta R},$$

$${\int }_{r_{\min }}^{r_{\max }}dr{r}^2{j}_{\ell }(k_{\ell n}r)j_{\ell }(k_{\ell n^{\prime }}r)\simeq \frac{{\delta }_{n{n}^{\prime }}}{{\tau }_{\ell n}},$$

$$\sum _n{\tau }_{\ell n}{j}_{\ell }(k_{\ell n}r)j_{\ell }(k_{\ell n}{r}^{\prime })=\frac{\delta (r-r^{\prime })}{r^2},$$

$${\tau }_{\ell n}^{-1}$$

$$T(\mathbf{r})=\sum _{\ell mn}{\tau }_{\ell n}{j}_{\ell }(k_{\ell n}r)Y_{\ell m}(\hat{\mathbf{r}})T_{\ell mn}.$$

$$T_{\ell mn}=\sqrt{\frac{2}{\pi }}\int dk{k}^2{W}_{\ell n}(k)T_{\ell m}(k),$$

$$W_{\ell n}(k)$$

Full text

††pubyear: 2018††pagerange: Full-Sky Lensing Reconstruction of 21 cm Intensity Maps–E

1 Introduction

Weak gravitational lensing by large-scale structure (LSS) (Blanchard & Schneider, 1987; Linder, 1988; Cole & Efstathiou, 1989) (see Bartelmann & Maturi (2017) for a recent review article) is a powerful and ubiquitous probe of cosmology. As an unbiased tracer of LSS it can be used to measure the distribution of matter (Chang et al., 2018) as well as the growth of structure (Huterer et al., 2015), an important probe of gravity. The two common tracers of the gravitational lensing potential are shear deformations in the shapes of galaxies (Prat et al., 2018; Hikage et al., 2018) and deflections in the cosmic microwave background (CMB) that sample both the convergence and shear fields (Planck Collaboration et al., 2018; Omori et al., 2017; Sherwin et al., 2017). A recent addition to the lensing tracers is the cosmic infrared background (CIB) (Schaan et al., 2018), in which the authors derive the optimal flat-sky lensing estimator for the CIB while accounting for the CIB trispectrum and self-lensing by low-redshift perturbations. The lensing potential also serves as a foreground for the $B$-mode CMB polarization through the lensing of the $E$-mode polarization field (Kamionkowski & Kovetz, 2016).

Line intensity mapping (IM) (Scott & Rees, 1990; Madau et al., 1997; Suginohara et al., 1999; Chang et al., 2008; Wyithe et al., 2008) has emerged as a new tracer of large-scale structure, in which the aggregate line emission from star-forming galaxies and the intergalactic medium (IGM) is tabulated into three-dimensional maps and has both astrophysical and cosmological applications. The most promising emission line being considered is the 21 cm line from neutral hydrogen, but several other lines are also being considered, such as CO, CII, and Ly$\alpha$ (see Kovetz et al. (2017) for the latest review on these lines). LIM has lately been considered as a promising new probe of the lensing potential. LIM lensing has mostly been considered in terms of the 3D Fourier-space quadratic estimator derived in Zahn & Zaldarriaga (2006) while numerous other estimators have also been considered (Mandel & Zaldarriaga, 2006; Metcalf & White, 2007; Lu & Pen, 2008; Metcalf & White, 2009; Lu et al., 2010; Pourtsidou & Metcalf, 2014, 2015; Romeo et al., 2018) LIM has significant advantages over traditional lensing tracers. For one, unlike CMB fluctuations, 21-cm fluctuations do not experience Silk damping on small scales (Hu & White, 1997); they are limited instead by the Jeans length of the gas which occurs at much smaller scales (Loeb & Zaldarriaga, 2004), providing more power to sample the lens field. Also, the CMB has only one source plane, while intensity maps will have multiple source planes which can sample the same lenses, greatly increasing the potential to constrain the gravitational perturbation field. Shear maps from galaxy lensing can also have multiple source planes, but these maps tend to use photometric redshifts which are more uncertain and limit the number of source planes that can be applied. Intensity maps are also optimized to reach redshifts all the way out through the reionization epoch and beyond, greatly extending the reach of lensing surveys. One disadvantage of intensity maps for lensing compared to CMB maps is that while the CMB is nearly Gaussian, intensity maps at lower redshifts will be non-Gaussian due to nonlinear clustering. However, recent work has determined various alternatives that can reduce this effect, including bias-hardened estimators (Foreman et al., 2018) and accounting for the trispectra (Schaan et al., 2018). Also, the foreground contamination for 21-cm intensity maps is orders of magnitude larger than for the CMB, which will degrade the ability to map the lensing potential.

The formalism for weak lensing of intensity maps using a quadratic estimator such as in Zahn & Zaldarriaga (2006) is typically written in terms of the plane-parallel formalism, where the 3D intensity maps are transformed into Fourier modes in terms of $\mathbf{k}=(k_{\parallel},\mathbf{k}_{\perp})$, where $k_{\parallel}$ is the component of the wavevector $\mathbf{k}$ along the line of sight, and $\mathbf{k}_{\perp}$ is the transverse component. This formalism is ideal for surveys over small areas such that the line of sight does not change significantly over the survey. However, some upcoming 21-cm surveys are currently being designed to survey very large areas. For example, the HIRAX 21-cm survey (Newburgh et al., 2016) will probe redshifts $0.8<z<2.5$ and map all the Southern Sky. In this limit, we cannot assume the line-of-sight is constant over the survey. In fact, we expect that Fourier analysis of intensity maps in this limit is not optimal, and that a new basis is required that accounts for the curvature of the sky.

A basis that serves this function is the spherical Fourier-Bessel basis (SFB). In this basis, the 3D map is decomposed into modes labeled by $(k,\ell)$, This basis has a rich history in cosmological theory and analysis, particularly for galaxy clustering and weak lensing (Binney & Quinn, 1991; Fisher et al., 1994; Fisher et al., 1995; Lahav et al., 1994; Zaroubi et al., 1995; Heavens & Taylor, 1995; Castro et al., 2005; LSST Dark Energy Science Collaboration, 2012; Rassat & Refregier, 2012; Shapiro et al., 2012; Leistedt et al., 2012; Pratten & Munshi, 2013; Yoo & Desjacques, 2013). Recently, Liu et al. (2016) proposed analyzing intensity maps in the SFB basis, and Castorina & White (2018) showed that SFB moments of 3D cosmological maps can be used directly to estimate the matter power spectrum in Fourier space. This paper will describe a powerful application for the SFB basis, namely weak gravitational lensing. We expect the SFB and plane-parallel lensing formalisms to deviate due to distorted correlations of the intensity field in mostly angular modes ($\ell\sim kR$), scales at which lensing can help answer cosmological questions. One such concern is constraining gravity, which has significant large-scale effects and for which lensing is an unbiased probe. These being large-scale effects imply that a SFB formalism is vital to modeling the effects in data correctly.

In this paper, we construct the formalism for estimating the gravitational potential from 3D intensity maps over the full sky, assuming Gaussian statistics in the intensity field. After reviewing the weak lensing estimator in the plane-parallel formalism as well as the spherical Fourier-Bessel (SFB) decomposition, we derive a minimum-variance quadratic estimator for the spherical harmonic moments of the weak lensing potential given the SFB series of a 3D intensity map. With this estimator, we then derive an expression for the noise bias of the angular power spectrum of the lensing potential $N_{\ell}^{\phi\phi}$ in this formalism. We also present the $N_{\ell}^{\phi\phi}$ prediction using the SFB flat-sky limit. We then compare both SFB predictions for $N_{\ell}^{\phi\phi}$ to the prediction from Zahn & Zaldarriaga (2006) using the plane-parallel formalism, assuming instrumental properties of the HIRAX 21 cm survey as a test case.

We find that both SFB predictions agree with the plane-parallel prediction at small scales but deviate dramatically at scales $L$ for the lensing potential where the estimator is sourced by some scales that are mostly angular in the intensity map with high signal-to-noise ratios. This corresponds to $L\lesssim\ell_{\rm res}-k_{\rm eq}R$ where $R$ is the radius of the shell volume, $k_{\rm eq}$ is the peak of the matter power spectrum set by the wavenumber that enters the cosmic horizon at matter-radiation equality, and $\ell_{\rm res}$ is the angular scale where the instrument beam increases the noise in angular harmonic space by 25%. Note that the full-sky and flat-sky SFB predictions deviate from the plane-parallel prediction at the same scale, while the full-sky SFB prediction deviates slightly from the flat-sky prediction for $\ell\lesssim 10$, similar to the full-sky vs. flat-sky deviation for CMB lensing. We then determine the sensitivity of $C_{\ell}^{\phi\phi}$ and $C_{\ell}^{\kappa g}$, the convergence-galaxy angular cross-power spectrum, for the HIRAX 21-cm intensity mapping survey (Newburgh et al., 2016) over a redshift range $1.38\leq z\leq 2.57$ using our new formalism, assuming the galaxy sample that will be produced by the Large Synoptic Survey Telescope (LSST) (LSST Dark Energy Science Collaboration, 2012). With a range of $100<L_{\rm dev}<200$ for the scale where $N_{\ell}^{\phi\phi}$ deviates from the plane-parallel prediction, we find that the SNR for $C_{\ell}^{\phi\phi}$ using our SFB estimator increases by 107%, while $C_{\ell}^{\kappa g}$ with galaxies from the upcoming Large Synoptic Survey Telescope only has a modest increase of 10% with much higher deviations if considering just large angular scales. Note that these results will change once non-Gaussian intensity fields due to nonlinear clustering are taken into account; we leave this for future work. These results suggest that at large angular scales intensity mapping could be a much more powerful probe of large-scale structure than previously thought, which could have many applications including improving measurements of gravity on the largest scales.

The plan of our paper is as follows: in Section 2 we review the form of the gravitational lensing potential and the noise bias from intensity map estimates in the plane-parallel formalism. In Section 3 we review the SFB formalism including our formalism for power spectrum predictions on a spherical shell. We present our full-sky $N_{\ell}^{\phi\phi}$ derivation in Section 4 and compare this result to the plane-parallel prediction in Section 5. We conclude in Section 6. We assume cosmological parameters consistent with the Planck 2015 cosmology (Planck Collaboration et al., 2016), namely $\Omega_{b}h^{2}=0.02225$, $\Omega_{m}h^{2}=0.1419$, $h=0.676$, $\Omega_{\nu}h^{2}=0.00064$, $N_{\rm eff}=3.046$, $n_{s}=0.964$, and $A_{s}=2.2\times 10^{-9}$.

2 Angular Power Spectrum of Lensing Potential

Let the 21 cm intensity field be $T(\hat{\textbf{r}},r)$ and the lensed field be $\widetilde{T}(\hat{\textbf{r}},r)$. For a lensing potential $\phi(\hat{\textbf{r}})$,

[TABLE]

We will use the conventions where the Fourier transform is given by

[TABLE]

the covariance between Fourier modes is written as

[TABLE]

where $P_{T}(k)$ is the power spectrum of the brightness temperature fluctuations and $k=|\mathbf{k}|$ is the mode wavenumber.

The lensing potential, which describes the weak lensing of extragalactic photons due to gravity, has been given extensively in the literature in the form

[TABLE]

where $\chi$ and $\chi_{s}$ are the comoving distances of the lenses and source plane. $\Phi$ is the gravitational potential, given in the Poisson equation as

[TABLE]

with $\Omega_{m,0}$ being the matter density today and $\delta$ being the matter overdensity. Inserting this into Eq. 4 and expanding $\delta$ in Fourier components, we can write the angular power spectrum for $\phi$ as

[TABLE]

where $P(k)$ is the matter power spectrum,

[TABLE]

and $D(z)$ is the growth function.

The estimator for the lensing potential in an intensity mapping context was first derived in Zahn & Zaldarriaga (2006) within a plane-parallel formalism. In order to account for a finite volume, $k_{\parallel}$ is discretized as $k_{\parallel,j}=2\pi j/\Delta R$ with $j=1,2,...$ and $\Delta R$ being the depth of the survey. In addition, $\mathbf{k}_{\perp}$ was rewritten as $\mathbf{k}_{\perp}=\bm{\ell}/R$ with $R$ being the distance to the source. In this formalism, the noise of the measured auto-power spectrum of the lensing potential $\phi$ at angular scale $\ell$ is given by

[TABLE]

where

[TABLE]

where $\beta=f/b$ is the growth rate divided by the linear bias and $j_{\rm min}$ denotes the minimum $k_{\parallel,j}$ that is not contaminated by foregrounds. Note that $C_{\ell}^{\rm N}$ is the noise power spectrum which will be defined in Sec. 5

3 Spherical Fourier Bessel Formalism

In this section we review the formalism for the SFB basis including a continuous $k$ and discrete $k$ values, after which we will derive the covariance of the discrete moments in terms of the power spectrum.

In the LIM context, a cosmology is assumed in order to relate the radial distance $r$ to redshift $z$ in the form

[TABLE]

where $H(z)$ is the Hubble rate in terms of the relative densities of matter and dark energy. The redshift $z$ is then related to the observed frequency of the survey line, such that $\nu=\nu_{\rm rest}/(1+z)$. In the linearized form, $z$ can be directly related to $r$ over a narrow ranged centered at $r_{\rm ref}(z_{\rm ref})$ according to

[TABLE]

3.1 Spherical Fourier Bessel Transform Review

In Liu et al. (2016), the SFB transform in three dimensions is given in the form

[TABLE]

where $k=|\mathbf{k}|$ is the mode wavenumber, $r=|\mathbf{r}|$ is the distance along the line-of-sight, $\mathbf{\hat{r}}=\mathbf{r}/r$ is the unit direction vector, $\ell$ and $m$ are the degree and order, respectively, of the spherical harmonic function $Y_{\ell m}(\mathbf{\hat{r}})$, and $j_{\ell}(x)$ is the spherical Bessel function of the first kind. The inverse of this transform is

[TABLE]

Any function defined on $\mathcal{R}^{3}$ can be written in this basis. A useful relation in this context is the identity

[TABLE]

It is shown in Liu et al. (2016) that for the case with an isotropic power spectrum in the temperature fluctuations $P_{T}(k)=b^{2}\bar{T}^{2}P(k)$ where $b$ is the clustering bias and $\bar{T}$ is the mean brightness temperature, the covariance between $T_{\ell m}(k)$ modes can be written as

[TABLE]

such that all the modes are independent. L16 also asserted that if rotational invariance of the fluctuations $T(\mathbf{r})$ is broken, e.g. redshift-space distortions, then the covariance between modes changes to

[TABLE]

In recent work (Castorina & White, 2018), the SFB power spectrum including redshift-space distortions (RSD) has been shown to be given in the form

[TABLE]

where the RSD parameter $\beta=f/b$ is the growth rate divided by the clustering bias and $W_{\ell}^{0}$ and $W_{\ell}^{r}$ are window functions given by

[TABLE]

where the second derivative of $j_{\ell}(x)$ can be written as

[TABLE]

Note that strictly speaking, $r$ is in redshift-space, not real-space, but this should have a small effect once the small survey window along the line of sight, which is also in redshift space, is taken into account; thus we will not mention this further.

For a partial volume all-sky survey defined by the window function $\phi(r)=1$ inside the survey and 0 outside, $W_{\ell}^{0}$ and $W_{\ell}^{r}$ in Eq. 3.1 are modulated by $\phi(r)$ in the integrand, in which case $W_{\ell}^{0}$ will no longer be a delta function. For an all-sky survey, L16 showed this leads to an estimator for $P(k)$ defined as

[TABLE]

where the $\ell$ modes can be combined to construct a minimum-variance estimator for $P(k)$.

3.2 Spherical Fourier Bessel Series

The spherical Fourier Bessel (SFB) transform can be discretized into a spherical Fourier Bessel series in a finite volume. This has typically been performed for a spherical volume (Fisher et al., 1995; Leistedt et al., 2012; Pratten & Munshi, 2013); however, more realistic surveys are constructed over spherical shells. Thus we will define our survey volume as a spherical shell with inner radius $r_{\rm min}$ and outer radius $r_{\rm max}$ when modeling a transformation using the spherical Fourier Bessel series.

For the sphere, the wavenumber $k$ in the transform must have discrete values to satisfy a boundary condition at $T(r=r_{\rm max})$. For a shell, a boundary condition at $T(r=r_{\rm min})$ must also be satisfied. In our paper, we use Dirichlet boundary conditions $T=0$ at both boundaries. Modeling the radial behavior as $j_{\ell}(kr)$ as in the SFB transform makes it impossible to satisfy both boundary conditions. However, we instead choose $k$ values that satisfy $j_{\ell}(kr_{\rm max})=0$ exactly, then picking values of $k$ that most closely satisfy $j_{\ell}(kr_{\rm min})=0$. This discrete SFB formalism is inspired by a full treatment for the shell which will be presented in Samushia (2019) and to which we had early access, in which $j_{\ell}(kr)$ is replaced by a linear combination of the spherical Bessel functions of the first and second kind, $j_{\ell}(kr)+Ay_{\ell}(kr)$, such that the boundary conditions at both $r_{\rm min}$ and $r_{\rm max}$ can be satisfied exactly.

We define moments for the map as $T_{\ell mn}$ given by

[TABLE]

where $k_{\ell n}=s_{\ell n}/r_{\rm max}$ and $s_{\ell n}$ is the $n^{\rm th}$ root of both $j_{\ell}(x)$ and $j_{\ell}(ax)$ where $a=r_{\rm min}/r_{\rm max}$, such that $j_{\ell}(s_{\ell n})=0\,\forall n$, while approximately satisfying $j_{\ell}(as_{\ell n})=0$. Note that for $a=0$ we attain the spherical case in which $j_{\ell}(as_{\ell n})=0$ is automatically satisfied for $\ell>0$, leaving only the boundary condition $j_{\ell}(s_{\ell n})=0\,\forall n$. In this case we define $q_{\ell n}$ as the $n^{\rm th}$ root of only $j_{\ell}(x)$.

Here we derive the lowest $k_{\ell n}$ allowed in our basis as well as the spacing between consecutive $k_{\ell n}$s. We use the approximation for the roots of the $j_{\ell}(x)$ given by $q_{\ell n}\simeq\ell+n\pi$ for $n=1,2,3,...$ (Leistedt et al., 2012). We will also only consider the case in which the volume is a thin shell, where $\Delta R=r_{\rm max}-r_{\rm min}\ll r_{\rm max}$ or equivalently $1-a\ll 1$. First we consider the case where $\ell\ll n\pi$ to find the lowest $s_{\ell 1}=q_{\ell n}$ by setting $n$. In this case we can approximate $q_{\ell n}\simeq n\pi$. Since the roots of $j_{\ell}(x)$ are approximately $\pi$ apart, and we need both $q_{\ell n}$ and $aq_{\ell n}$ to be roots, we can pick these roots to have the shortest possible difference of $\pi$. This leads to the condition that $n=1/(1-a)=r_{\rm max}/\Delta R$. Thus the first value for $s_{\ell n}$ is $s_{\ell 1}=\ell+\pi r_{\rm max}/\Delta R$. To get subsequent $s_{\ell n}$s, given an $s_{\ell n}$ we can always define $q_{\ell n^{\prime}}^{\prime}=s_{\ell n}+n^{\prime}\pi$ as the set of subsequent roots of $j_{\ell}(x)$. We can then repeat the previous exercise showing that the next root of both $j_{\ell}(x)$ and $j_{\ell}(ax)$ is a distance $\pi r_{\rm max}/\Delta R$ away. Thus in the case $\ell\ll n\pi$ we can set

[TABLE]

where $n=1,2,3,...$ and the spacing between $k_{\ell n}$s is $\Delta k=\pi/\Delta R$.

Next we consider the case $\ell\gtrsim n\pi$. Now the first roots for $j_{\ell}(x)$ and $j_{\ell}(ax)$ are approximately $\ell$ and $a\ell$, respectively. We then want to find a common root. Since $a<1$ we should find the subsequent roots for $j_{\ell}(ax)$; since $\ell\gg\pi$, one of these roots will probably be equal to the first root of $j_{\ell}(x)$, namely $x\simeq\ell$ where $\ell\gg 1$. Thus we set $\ell=a(\ell+n\pi)$, giving us $n=\ell(1-a)/(\pi a)$. Substituting this in $q_{\ell n}$ we find $s_{\ell 1}=\ell/a$. Getting subsequent $s_{\ell n}$s for this case is not trivial. Unlike the case where $\ell\ll n\pi$, we cannot set all subsequent roots of both $j_{\ell}(x)$ and $j_{\ell}(ax)$ to be multiples of $\pi$ away from $s_{\ell 1}$ since they are inherently offset. In this regime the $k$-spacing can be much smaller than $\pi/\Delta R$. Also in this regime the $x$ values that satisfy both $j_{\ell}(x)$ and $j_{\ell}(ax)$ tend to be well in between the $q_{\ell n}$s, causing the chosen $s_{\ell n}$s not to satisfy $j_{\ell}(as_{\ell n})$ very well. This is unavoidable as long as our basis functions only consist of $j_{\ell}$s. However, at high enough $n$ the offset mentioned earlier will become negligible, and from then on the $k_{\ell n}$s will be spaced as $\Delta k=\pi/\Delta R$.

Note that in both cases since $s_{\ell 1}=k_{\ell 1}r_{\rm max}>\ell$, the SFB basis does not probe modes $kr_{\rm max}<\ell$. This is similar to L16 in which modes $kr_{\rm max}<\ell$ have extremely low signal-to-noise for power spectrum measurements and correspond to modes $k_{\perp}>k$ in the plane-parallel limit, which is unphysical. In addition, the boundary condition alone does not account for only using an intensity map within the shell; we fully account for this in our theoretical derivations by only integrating Eq. 21 in the region $r_{\rm min}\leq r\leq r_{\rm max}$.

The orthogonal conditions for the shell’s SFB series are

[TABLE]

and

[TABLE]

where

[TABLE]

Note that since $j_{\ell}(k_{\ell n}r_{\rm min})$ sometimes differs from zero a fair amount for $\ell\gg\pi n$, the integral in Eq. 23 and the sum in Eq. 24 are only approximately diagonal. This allows the series for $T(\mathbf{r})$ to be expanded as

[TABLE]

The SFB transform basis and the SFB series basis both contain all the information on $T(\mathbf{r})$ within the radial range $r_{\rm min}<r<r_{\rm max}$. The discrete $k_{\ell n}$’s Nyquist sample the map in such a way that continuous $k$’s do not contain more (or less) information with a finite correlation length $\Delta k=\pi/\Delta R$.

We can write the moments for the SFB series $T_{\ell mn}$ in terms of the moments for the SFB transform $T_{\ell m}(k)$ by inserting Eq. 13 in Eq. 21 to find

[TABLE]

where $W_{\ell n}(k)$ is given by

[TABLE]

For $k_{\ell n}r_{\rm max}\gg\ell$, $W_{\ell n}(k)$ is a highly-peaked function at $k=k_{\ell n}$, and for all cases $W_{\ell n}(k_{\ell n})=1/\tau_{\ell n}$.

With these expressions, the covariance among modes $T_{\ell mn}$ can be written, assuming an intrinsic, isotropic power spectrum, as

[TABLE]

where

[TABLE]

Note that the even in the isotropic $P(k)$ case the covariance matrix $C_{\ell nn^{\prime}}$ is only approximately diagonal in $n$. This is just due to the survey volume being a finite region. This effect also occurs for the SFB transform in a finite volume. Inserting Eq. 17, we can write the above equation in the form

[TABLE]

where

[TABLE]

Adding instrumental noise and effects from the instrumental beam and spectral window functions, this expression augments to

[TABLE]

where we define $W_{\ell nn^{\prime}}^{AB}\equiv W_{\ell nn^{\prime}}^{A}W_{\ell}^{B}$ such that

[TABLE]

where

[TABLE]

and $\Delta r$ is the comoving distance corresponding to the spectral resolution and $W_{\ell}^{B}=e^{\ell^{2}\sigma_{b}^{2}}$, where $\sigma_{b}$ is the angular resolution. This instrumental noise is given by

[TABLE]

where $z(r)$ is the redshift given comoving distance $r$, $\sigma_{N}(r)$ is the instrumental noise given as a function of $r$ which is related to frequency by Eq. 11, and $\Delta(r)$ is the voxel size given by

[TABLE]

where $\Delta\Omega$ and $\nu$ are the solid angle of the angular pixels and $\Delta\nu$ is the size of the frequency channels. These instrumental effects are derived in Appendix B. Note that since $\sigma(r)$ and $\Delta(r)$ are smooth functions over the shell, we can take them out of the integral setting them with the values $\sigma(R)$ and $\Delta(R)$ where $R=(r_{\rm min}+r_{\rm max})/2$. With this approximation the integral can be performed analytically, allowing us to set

[TABLE]

3.3 SFB Flat-Sky Approximation

Here we consider a flat-sky approximation to the SFB basis in order to facilitate comparisons with the plane-parallel lensing formalism. The flat-sky limit is characterized by $\theta\rightarrow 0$, i.e., small angles around the pole in angular coordinates on the sphere $(\theta,\varphi)$, which translates into $\ell\gg 1$ in harmonic space. In the flat-sky coordinates, angular quantities are described by $\bm{\ell}=[\ell\cos(\varphi_{\ell}),\ell\sin(\varphi_{\ell})]$, where $\ell$ is the continuum approximation of the Spherical Harmonic index. In this approximation, sums transition into integrals, for which we use $\int d^{2}\ell=\int_{0}^{\infty}\ell d\ell\int_{0}^{2\pi}d\varphi_{\ell}\approx\sum_{\ell}\ell\int_{0}^{2\pi}d\varphi_{\ell}$.Inspired by the formalism from Castro et al. (2005) but converting to an SFB series formalism, we define the SFB flat-sky series transform as

[TABLE]

and

[TABLE]

We note the conversion from SFB full-sky moments to the flat-sky moments (and vice versa) are

[TABLE]

Further, one can derive the covariance

[TABLE]

It is evident that the SFB flat-sky basis is not strictly equivalent to the plane-parallel basis, given by

[TABLE]

with a power spectrum $\left\langle T^{\rm obs}(k_{||,j},\bm{\ell})T^{*,\rm obs}(k_{||,j^{\prime}},\bm{\ell}^{\prime})\right\rangle=(2\pi)^{2}\delta^{2}(\bm{\ell}-\bm{\ell}^{\prime})\delta_{jj^{\prime}}(R^{2}\Delta R)^{2}\tilde{C}_{\ell}(k_{||,j})$ given by

[TABLE]

where

[TABLE]

However, it’s useful to consider when the approximations

[TABLE]

where $C_{\ell}(k_{||})=P_{T}(k_{||},\ell)/(R^{2}\Delta R)$ are true. The approximation $\tilde{C}_{\ell}(k_{||,j})\simeq C_{\ell}(k_{||,j})$ is justified in that $W(k_{||}^{\prime},k_{||,j})$ is highly peaked at $k_{||}^{\prime}=k_{||,j}$ such that

[TABLE]

This has allowed previous treatments of the weak lensing of intensity maps to set the power spectrum as $C_{\ell}(k_{||,j})$ and not $\tilde{C}_{\ell}(k_{||,j})$.

In the SFB case the approximation $\tau_{\ell n}C_{\ell n}\simeq R^{2}\Delta RC_{\ell}(k_{||,j})$ is only satisfied if the integrand of $C_{\ell n}$ is highly peaked, such that

[TABLE]

The condition that the integrand be highly peaked is satisfied when $\ell\ll k_{\ell n}R$ and $\ell\ll kR$ such that the asymptotic expansion of $j_{\ell}(x)$ can be applied and $W_{\ell n}(k)$ approaches a sinc function peaked at $k=k_{\ell n}$ of the form

[TABLE]

The condition $\ell\ll k_{\ell n}R$ is generally true for most $k_{\ell n}$, and $W_{\ell n}(k)$ tends to not get any other peaks even for $kR\lesssim\ell$. Thus we can remove $P_{T}(k)$ from the integral which leaves us with $\tau_{\ell n}C_{\ell n}\simeq P_{T}(k_{\ell n})$ for $k_{\ell n}R\gg\ell$. However, for $k_{\ell n}R\sim\ell$, $W_{\ell n}(k)$ is not so highly peaked, causing mode mixing and giving us $\tau_{\ell n}C_{\ell n}<P_{T}(k_{\ell n})$.

Further, we wish to examine how $C_{\ell n}$ behaves at $\ell>k_{\rm eq}R$. Particularly, it is necessary to highlight the contribution of the integral from the region $kR<\ell$. In this regime, $j_{\ell}(kr)\approx\frac{(kr)^{\ell}}{(2\ell+1)!!}\sim\left(\frac{kr}{2\ell}\right)^{\ell}$, and thus,

[TABLE]

since the radial integrand here is bounded by $1$ and only contributes to suppression, the exact nature of which is irrelevant here. Thus we have shown that we may safely ignore the contribution from $kR<\ell$. Following the arguments for $k_{\ell n}R\gg\ell$, we have,

[TABLE]

Indeed, this may be corroborated from Fig. 1 for the specifications of the survey used in this paper (for which $P_{T}(k)$ is the matter-power spectrum). Particularly for the case $k_{\ell n}R>\ell>k_{\rm eq}R$, we have $P_{T}(k)\sim 1/k^{2}$, and given that $k_{\ell n}\sim\ell$, we may expect that $\tau_{\ell n}C_{\ell n}\sim\frac{1}{\ell^{2}}$, which is indeed suppressed.

4 Full Sky weak lensing reconstruction

Here we derive the lensing estimator for a 21 cm intensity map, applying a formalism similar to that derived in Okamoto & Hu (2003) for full-sky CMB lensing. We start with a 21 cm signal, perturbed weakly by a lensing potential. Expanding the lensing equation to first order produces

[TABLE]

Taking the SFB transform of this equation for the perturbation in the 21 cm field, $\delta\widetilde{T}$, produces its SFB moments

[TABLE]

where $W$ is defined as

[TABLE]

where $a\equiv r_{\rm min}/r_{\rm max}$ and

[TABLE]

In order to obtain this, we have used the well known identity:

[TABLE]

Note that $W_{\ell\ell^{\prime}}^{nn^{\prime}}$ serves as a coupling matrix that determines which unlensed modes source the lensed modes. In order to determine the behavior of this coupling, we use the asymptotic approximation for $j_{\ell}(kr)$ at large $\ell$ in Appendix D given by

[TABLE]

where $k_{||}=\sqrt{k_{\ell n}^{2}-(\ell/R)^{2}}$ and $\varphi_{k}={\rm arctan}(\ell/Rk_{||})$. Now, upon integrating $W_{\ell\ell^{\prime}}^{nn^{\prime}}$ using this approximation, and after some reorganization and letting the (un)primed $k_{||}$ correspond to the (un)primed $n$ and $\ell$, one gets

[TABLE]

Note that $k_{||}=2\pi j/\Delta R$ produces

[TABLE]

This shows that the coupling $W_{\ell\ell^{\prime}}^{nn^{\prime}}$ makes modes with the same $k_{||}$ source each other in weak lensing for the SFB formalism, just as in the plane-parallel formalism. Computationally, we have found that this conclusion is true for all $\ell$ and all but the first few of $k_{\ell n}$ which tend to couple with the first few $k_{\ell^{\prime}n^{\prime}}$’s, and we have found that neglecting these extra couplings do not affect our final results. Thus, in our analysis while we perform the full integral in Eq. 55, we choose to only use the entries in the coupling matrix that correspond to $k_{||}=k_{||}^{\prime}$, or $\sqrt{k_{\ell n}^{2}-(\ell/R)^{2}}=\sqrt{k_{\ell^{\prime}n^{\prime}}^{2}-(\ell^{\prime}/R)^{2}}$ in order to dramatically reduce the computational effort. Since this equality is never satisfied exactly, we instead find the $n^{\prime}$ that most satisfies it for every $(\ell,n,\ell^{\prime})$. Also, when $\varphi_{k}\ell\simeq\varphi_{k^{\prime}}\ell^{\prime}$ one can use the asymptotic expansion to show that $\sqrt{\tau_{\ell n}\tau_{\ell^{\prime}n^{\prime}}}\sim 2k_{\ell n}k_{\ell^{\prime}n^{\prime}}/\Delta R$, which motivates the approximation $W_{\ell\ell^{\prime}}^{nn^{\prime}}\simeq\delta_{\ell\ell^{\prime}}^{nn^{\prime}}/\sqrt{\tau_{\ell n}\tau_{\ell^{\prime}n^{\prime}}}$ where $\delta_{\ell\ell^{\prime}}^{nn^{\prime}}=1$ when the modes have the same $k_{||}$ and 0 otherwise. This is only accurate at the 20% level, so we do not use it for our main results; however, it is accurate for comparisons between different surveys where computing the full integrals would be very computationally expensive.

Using Eq. 54, it is possible to have the covariance matrix of the lensed signal for which we take the ensemble average over a fixed potential

[TABLE]

where $f_{\ell\ell_{1}\ell^{\prime}}^{nn^{\prime}}=M_{\ell\ell^{\prime}}^{nn^{\prime}}F_{\ell^{\prime}\ell_{1}\ell}+M_{\ell^{\prime}\ell}^{n^{\prime}n}F_{\ell\ell_{1}\ell^{\prime}}$, and $M$ is defined as

[TABLE]

which can be written as an integral by using the integral forms of $W_{\ell\ell^{\prime}}^{nn^{\prime}}$ and $C_{\ell nn^{\prime}}$ and the orthogonality condition from Eq. 24 to find

[TABLE]

where

[TABLE]

and

[TABLE]

We also have the property $M_{\ell\ell}^{nn^{\prime}}=C_{\ell nn^{\prime}}$.

In the following work, $C_{\ell nn^{\prime}}$ is the unlensed $T-T$ covariance, without any instrumental noise. Also, here we clarify that we assume that the temperature maps have already had their beam and spectral window profiles deconvolved. For the discussion that follows, note that the covariance with instrumental noise is denoted as

[TABLE]

where $C_{\ell nn^{\prime}}^{N}=N_{\ell nn^{\prime}}/W_{\ell nn^{\prime}}^{AB}$.

Now, we are prepared to consider a general quadratic estimator for the lensing potential. We use the following form, inspired by Okamoto & Hu (2003),

[TABLE]

where $A_{\ell}$ and $g_{\ell_{1}\ell_{2}\ell}^{n_{1}n_{2}}$ are to be determined. This is the most general quadratic estimator of the projected potential and we provide a proof of this fact in Appendix  C. We require an unbiased estimator such that it produces the potential upon taking an ensemble average over many realizations of the 21 cm signal. To achieve this, we set the normalization condition

[TABLE]

This condition straightforwardly produces the constraint

[TABLE]

Finally, we require that the yet undetermined $g_{\ell_{1}\ell_{2}\ell}^{n_{1}n_{2}}$ minimize the variance of this estimator. To achieve this, we first compute the covariance of the estimator

[TABLE]

In order to minimize the variance, we require the functional derivative of $A_{\ell}$

[TABLE]

Minimising the variance thus produces

[TABLE]

where $S$ is defined to be

[TABLE]

One will notice that $g_{\ell_{1}\ell_{2}\ell}^{n_{1}n_{2}}$ is taken to be real in this equation. It follows straightforwardly from the parity properties of both $\phi_{\ell m}$ and $T_{\ell mn}$, and the properties of the 3j-symbol that in fact $g_{\ell_{1}\ell_{2}\ell}^{n_{1}n_{2}}{}^{*}=(-1)^{\ell_{1}+\ell_{2}+\ell}g_{\ell_{1}\ell_{2}\ell}^{n_{1}n_{2}}$. However, the presence of $F_{\ell_{1}\ell\ell_{2}}$ in Eq. 72, due to the properties of the 3j-symbol, forces the reality of $g_{\ell_{1}\ell_{2}\ell}^{n_{1}n_{2}}$, so we make this change early in hindsight for convenience. With this in mind, we find that Eq. 72, through the definition of $A_{\ell}$ in Eq. 69, leads to the equality

[TABLE]

Using (column permutation) properties of the 3j-symbol once again, one can derive that $g_{ab\ell}^{n_{1}n_{2}}=(-1)^{a+b+\ell}g_{ba\ell}^{n_{2}n_{1}}$, where the phase factor is again omitted in hindsight for the purpose of convenience. Thus, we get

[TABLE]

For convenience, we temporarily define $\xi_{ab\ell}^{n_{1}d}$ as

[TABLE]

We solve for $\xi$ first, the equation for which is now

[TABLE]

Now, this is linear so we simply invert it. To be specific, for fixed $a$, $(C^{N}_{a}{}^{-1})_{nn^{\prime}}$ is the $(n,n^{\prime})$ element of the inverse of the (symmetric) matrix $C_{a}$, i.e., the $N\times N$ matrix with elements $(C_{a})_{nn^{\prime}}=C_{ann^{\prime}}$. Solving this now presents

[TABLE]

Proceeding identically as above to solve to solve for $g$ presents

[TABLE]

This completes our derivation of the minimum variance weighting for the lensing potential. Further, we have for the lensing potential

[TABLE]

We assume a fiducial survey, so we have for the noise power spectrum of the lensing potential the expression

[TABLE]

Additionally, we derive the flat-sky limit of our estimator in order to draw a rough comparison with the estimator of Zahn & Zaldarriaga (2006). The flat-sky noise power spectrum is:

[TABLE]

where $\bar{f}_{\ell_{1}\ell\ell_{2}}^{n_{1}n_{2}}=M_{\ell_{1}\ell_{2}}^{n_{1}n_{2}}\bm{\ell}\cdot\bm{\ell_{1}}+M_{\ell_{2}\ell_{1}}^{n_{2}n_{1}}\bm{\ell}\cdot\bm{\ell_{2}}$ and $\bm{\ell_{2}}=\bm{\ell}-\bm{\ell_{1}}$. In Appendix A, we present a full derivation of the flat-sky limit of our estimator and the above expression.

We perform our forecast by approximating the covariance matrix as diagonal, $C_{\ell nn^{\prime}}^{\rm tot}=C_{\ell n}^{\rm tot}\delta_{nn^{\prime}}$. This is justified since for most $k_{\ell n}$, $W_{\ell n}(k)$ and $W^{r}_{\ell n}(k)$) are highly peaked and the $k_{\ell n}$ values are separated enough to have little window function overlap. For $k_{\ell n}R\sim\ell$, we find that while $C_{\ell nn^{\prime}}$ is still diagonal in the non-RSD case ($\lvert C_{\ell nn^{\prime}}\rvert<10\%$ of $C_{\ell nn}$), we see that RSD can cause the off-diagonal terms to be significant. However, we expect that this should not affect the overall behavior of our results. We leave this investigation for future work. Considering the diagonal sets the quantity $M_{\ell\ell^{\prime}}^{nn^{\prime}}=\tau_{\ell n}W_{\ell\ell^{\prime}}^{nn^{\prime}}C_{\ell n}$ while trivially rewriting the sum

[TABLE]

Note that applying the approximation $\tau_{\ell n}C_{\ell n}\simeq R^{2}\Delta RC_{\ell}(k_{||,j})$ to the SFB flat-sky expression for $N_{\ell}^{\phi\phi}$ will give you the plane-parallel expression in Eq. 8. This suggests that the plane parallel and SFB (flat and full sky) formalisms will deviate when this approximation breaks down.

We already know the approximation breaks down for $\ell\sim k_{\ell n}R$. If the approximation were always valid, then in the low instrumental noise regime the sums in $N_{\ell}^{\phi\phi}$ would stays constant at low $\ell$ where the $k_{\ell n}$’s are approximately $\ell$-independent and would increase when the $k_{\ell n}$’s are approximately $\ell$-dependent around $\ell=2\pi R/\Delta R$. This is what occurs in the plane-parallel formalism and gives the same result, in particular a constant $N_{\ell}^{\phi\phi}$ at low $\ell$. However, because the approximation breaks down, at $k_{\ell n}\sim\ell$ the situation changes where the autocorrelation of modes with the same value of $k_{\ell n}$ but different $\ell$’s differ, causing the sums in $N_{\ell}^{\phi\phi}$ to increase for $\ell_{1}>k_{\rm eq}R$ and $\ell_{2}>\ell_{1}$ where $C_{\ell n}$ decreases with $\ell$. In addition, we expect that the sum in the lensing noise power spectrum $N_{\ell}^{\phi\phi}$ will terminate when the instrumental noise power spectrum $C_{\ell n}^{N}$ starts to increase due to finite angular resolution. Thus, we expect the SFB prediction for $N_{\ell}^{\phi\phi}$ to deviate from the plane-parallel prediction when there are a significant number of angular modes $k_{\rm eq}R<\ell<\ell_{\rm res}$, where modes $\ell>\ell_{\rm res}$ are beam smeared. In particular, the modes $T(\bm{\ell}_{2})$ that are correlated with $T(\bm{\ell}_{1})$ due to $\phi(\bm{L})$ are constrained by the range $|\bm{L}-\bm{\ell}_{1}|<\bm{\ell}_{2}<\bm{L}+\bm{\ell}_{1}$ such that for small enough $\bm{L}$, $\bm{\ell}_{1}$ and $\bm{\ell}_{2}$ will be close enough together and in the right range so that both $\tau_{\ell_{1}n_{1}}C_{\ell_{1}n_{1}}$ and $\tau_{\ell_{2}n_{2}}C_{\ell_{2}n_{2}}$ will deviate from $R^{2}\Delta RC_{\ell}(k_{||,j})$ while still being contributing to the sum. Specifically, we set $L_{\rm dev}$ as the scale at which for $L<L_{\rm dev}$ the SFB and plane-parallel predictions for $N_{\ell}^{\phi\phi}$ deviate, and we define it by setting the condition $\bm{\ell}_{\rm 2,max}=\bm{L}+\bm{\ell}_{1}<\ell_{\rm res}\forall\,L<L_{\rm dev}$ when $\bm{\ell}_{1}>k_{\rm eq}R$. Simplifying these conditions gives us the deviation scale

[TABLE]

and this equation seems to work well when we set $\ell_{\rm res}$ to be the angular mode where the instrumental noise is 25% higher than the value at $\ell=0$ (see Sec. 5). Note that this is correspondingly where the terms for the sums in $N_{\ell}^{\phi\phi}$ decrease by 50%. As a physical interpretation, we expect that for low angular resolution (low $\ell_{\rm res}$) or low curvature (high $R$), we cannot distinguish modes from different directions in harmonic space, such that the lensing noise reduces to the plane-plane parallel case.

5 Comparison between estimators

Here we compare our full-sky and flat-sky expressions for $N_{\ell}^{\phi\phi}$ to that from Zahn & Zaldarriaga (2006). Similar to the forecast in Foreman et al. (2018), hereafter F18, we assume the survey to be over the redshift range $1.38\leq z\leq 2.57$. For the 21-cm signal, we assume the same form used in F18 relevant at low-redshifts given as

[TABLE]

where $\Omega_{\rm HI}=5\times 10^{-4}$ (Masui et al., 2013). We write the power spectrum ${P_{T}(k,z)=b_{\rm HI}^{2}(z)T^{2}(z)P(k,z)}$, where ${P(k,z)}$ is the matter power spectrum computed from CAMB (Lewis et al., 2000) with the nonlinear clustering regime computed from HALOFIT (Smith et al., 2003), and $b_{\rm HI}$ is the clustering bias of HI emitters set based on predictions from Castorina & Villaescusa-Navarro (2017) to the values in Table 1. Note that we also include the anisotropy from redshift-space distortions (RSD), specifically the Kaiser effect. For the SFB power spectra, the Kaiser effect is given explicitly in Eq. 31, while for the plane-parallel power spectra this is accounted by adding an extra factor $(1+\beta\mu^{2})^{2}$ to the power spectrum $P(k)$ where $\beta=f_{g}/b_{\rm HI}$ and the growth rate $f_{g}(z)=[\Omega_{m}(z)]^{0.55}$. The SFB power spectrum for the signal is shown in Fig. 1.

We assume the noise properties of the HIRAX survey (Newburgh et al., 2016) given in Foreman et al. (2018) to forecast the noise power spectra. In particular, we assume the survey covers half the sky with an angular resolution of 11 arcmin at $z=2$. For the plane-parallel formalism, the noise power spectrum for HIRAX is given by

[TABLE]

where $T_{\rm sys}(\nu)$ is instrumental thermal noise, $B$ is the bandwidth corresponding to the redshift window, $A_{\rm pix}$ and $t_{\rm pix}$ are the angular area and observing time per angular pixel, and $W(\ell)$ is the beam window function. For all these quantities we use the same values as in F18. For the SFB noise power spectrum, we set $C_{\ell n}^{N}$ in terms of $C_{\ell}^{N}$ based on Eq. 38 as

[TABLE]

Using this formalism, we plot $C_{\ell n}$ and $C_{\ell n}^{N}$ in Fig. 1. We see that for the SFB formalism, the signal is dominant over the noise up to $k_{\ell n}\simeq 0.7$ $h$/Mpc when the noise begins to exceed the signal. In particular, the beam window function eliminates scales greater than around $\ell=1000$ and the spectral window function starts to matter around $k_{\ell n}\simeq 2$ $h$/Mpc.

In our $N_{\ell}^{\phi\phi}$ forecasts we consider cuts to avoid the contamination from continuum foregrounds. In the plane-parallel (PP) basis for a single-dish survey the contamination is confined to low-$k_{\parallel}$ modes. It was shown in L17 that in the SFB basis low-$k$ modes are contaminated over the range $\sim 1/\Delta R$ except at large $\ell$ where the contamination reaches out to higher $\ell$. In Appendix E we explain this behavior in that the contours of constant contamination are equivalent to contours of constant $k_{||}=\sqrt{k^{2}-(\ell/R)^{2}}$. Thus in order to simulate equivalent foreground cuts for both the plane-parallel and SFB formalisms, we set $k_{\rm min}=\sqrt{k_{\rm||,min}^{2}+(\ell/R)^{2}}$, where $k_{\rm||,min}=2\pi j_{\rm min}/\Delta R$ and $j_{\rm min}$ is the minimum $j$ allowed by foreground cuts in the PP result.

We also simulate the removal of nonlinear modes by setting $k_{\rm max}=0.5h$/Mpc. Note that we do not use bias-hardened estimators in either formalism as introduced in F18, which will be needed to remove correlations due to nonlinear clustering. We leave this for future work.

We plot forecasts for $N_{\ell}^{\phi\phi}$ for the full-sky SFB formalism in Fig. 2 along with forecasts assuming the flat-sky SFB formalism we derived and the flat-sky, PP formalism from Zahn & Zaldarriaga (2006). We simulate reasonable foreground cuts by setting $j_{\rm min}=4$. We see that the PP and SFB forecasts agree at small scales, but deviate dramatically at $\ell_{\rm dev}\simeq 144$ and differing by as much as a factor of 26 at $\ell=2$. This would imply that $C_{\ell}^{\phi\phi}$ at the largest scales could be more accessible than previously thought. Note that with R(z=2)=3594 Mpc/$h$ and $k_{\rm eq}=0.02h$/Mpc for our fiducial cosmology and $\ell_{\rm res}=205$, we would guess a value $\ell_{\rm dev}=\ell_{\rm res}-k_{\rm eq}R=133$ which is close to where the plane-parallel and SFB curves deviate in Fig. 2. We also see that the deviation between the SFB flat and full sky cases is small and takes place for $\ell\lesssim 10$ with a magnitude similar to the CMB flat and full sky cases. The SFB full-sky vs flat-sky split location $\ell\sim 10$ is a result of the difference between $F_{\ell_{1}L\ell_{2}}$ and $\bm{L}\cdot\bm{\ell_{1}}$ which leads to a fractional difference for $N_{\ell}^{\phi\phi}$ of $1/(L+1)^{2}$, implying a >1% difference for $\ell\lesssim 10$.

5.1 Varying survey conditions

In this subsection we consider the effects of varying foreground cuts, the survey shell volume, and the angular resolution of the survey. Note that in all our examples we only compare $N_{\ell}^{\phi\phi}$ curves using SFB flat-sky instead of SFB full-sky since they are equal in most of the $\ell$-range and SFB flat-sky is much faster to compute. Also, based on our argument in the previous section we do not expect the SFB full-sky vs flat-sky split location $\ell\sim 10$ to vary with survey parameters. First, we consider how $N_{\ell}^{\phi\phi}$ is affected by the level of foreground cuts. We plot the SFB flat-sky forecasts with no foreground cut and with various values for $j_{\rm min}$ between 1 and 10 in Fig. 3.

First we find that although some modes are excised for $j_{\rm min}=1$, $N_{\ell}^{\phi\phi}$ is not affected by this cut ($\lesssim 2$% difference). Next, we find we would gain at most a 50% reduction in the noise by applying less stringent cuts. Finally, we see that in the limit of very stringent cuts, $j_{\rm min}\sim 10$, all the modes with $k_{\ell n}\sim\ell/R$ are removed such that $k_{\ell n}\gg\ell/R$ which is equivalent to $k_{\ell n}\sim k_{||}$ and that the result approaches a form similar to PP. It appears that $j_{\rm min}$ plays a role not in altering the behavior of $L_{\rm dev}$ but in fact the slope of $N_{\ell}^{\phi\phi}$ around $L_{\rm dev}$. To see why this is the case, note that the summand for $N_{\ell}^{\phi\phi}$ is roughly constant for low $\ell$s but is strongly suppressed at high $\ell$s ($\ell>\ell_{\rm res}$). While this suppression is immune to the choice of $j_{\rm min}$, the low $\ell$ behavior is suppressed by the increase of $j_{\rm min}$. Thus increasing $j_{\rm min}$ naturally produces a lower slope. Most importantly, we find that the reduction in the noise at low $\ell$ is present for reasonable levels of foreground cutting and will be prudent to consider.

Next we consider how $N_{\ell}^{\phi,\phi}$ changes if we change the survey volume redshift parameters $z$ and $\Delta z$. Note that we reduce the computational effort by using the approximation $W_{\ell\ell^{\prime}}^{nn^{\prime}}\simeq\delta_{\ell\ell^{\prime}}^{nn^{\prime}}/\sqrt{\tau_{\ell n}\tau_{\ell^{\prime}n^{\prime}}}$ for all $N_{\ell}^{\phi\phi}$ curves. In Fig. 4 we plot the PP and SFB $N_{\ell}^{\phi\phi}$ curves for redshifts $z=1.43$ and $z=2.46$. Based on the corresponding values of $\ell_{\rm res}$ and $R$ (see Fig. 4 for actual values), our formula would predict $\ell_{\rm dev}=182$ and 105, which are close to the apparent locations in the figure. Note that the curve for $z=2.46$ is lower because $T_{\rm sys}$ is higher leading to higher $C_{\ell}^{N}$. In Fig. 5 we plot the PP and SFB $N_{\ell}^{\phi\phi}$ curves for redshifts $z=2$ with both $\Delta z=0.16$ and 0.32. In the figure, $\ell_{\rm dev}$ is the same for both curves, which is consistent with our model since $\ell_{\rm res}$ is independent of $\Delta R$. We do see the noise decrease for $\Delta z=0.32$ since the bandwidth $B$ is longer.

Finally we consider what happens when we reduce the angular resolution in Fig. 6. We model this for the HIRAX survey by reducing the number of dishes in both orthogonal directions of the array ($n_{x}$ and $n_{y}$) by half from 32 to 16. Based on the corresponding values of $\ell_{\rm res}$ (see Fig. 6 for actual values), our formula would predict $\ell_{\rm dev}=134$ and 28 for the higher and lower resolution cases, respectively. The predicted $\ell_{\rm dev}$ value for the lower resolution case is a bit higher than the apparent location in the figure, which is more like $\ell\sim 15$. It is not obvious what is causing the difference; it could be that for a changing angular resolution our rule of a 25% increase in $C_{\ell n}^{N}$ is not adequate, but it is more likely that it just appears worse than the other examples because $\ell_{\rm dev}$ is so low in this case. Either way, the result is still of the right order, and it still confirms the directional behavior we predicted in our $\ell_{\rm dev}$ model. We also see that the overall noise for lower angular resolution is higher because the total number of beams $n_{\rm beam}\sim 4n_{x}n_{y}$ is inversely proportional to the noise $C_{\ell}^{N}$.

5.2 Effects to lensing power spectrum measurements

With our prediction for $N_{\ell}^{\phi\phi}$, we can now forecast the uncertainty on $C_{\ell}^{\phi\phi}$, given by

[TABLE]

where $f_{\rm sky}=0.5$ is the assumed sky fraction mapped by the HIRAX survey and we take the limit $C_{\ell}^{\phi\phi}\ll N_{\ell}^{\phi\phi}$. We find that the uncertainty for $C_{\ell}^{\phi\phi}$ in the SFB formalism decreases significantly in comparison to the plane-parallel case at scales where spherical geometry becomes relevant. In Table 2, we list the SNR values of $\phi\phi$ autocorrelation for the Plane-Parallel case and our SFB case over the redshift range $1.38\leq z\leq 2.57$, which corresponds to the lower half of the HIRAX frequency range. Particularly, we find an increase of 107% for the SNR. We note, however, that the SNR values are likely overestimated due to the fact that we have not used bias-hardened estimators for this calculation, as well as that nonlinear clustering was not included in the fiducial power spectrum. We expect to derive relevant results for a bias-hardened full-sky estimator for $\Phi$ in future work.

Following the example in F18, we also consider the uncertainty on the angular cross-power spectrum $C_{\ell}^{\kappa g}$, where $\kappa=-\nabla^{2}\phi/2$ is the lensing convergence and $g$ represents a galaxy distribution. For cross-correlations, $\kappa$ is generally favored since it directly traces the matter perturbations. The expression for $C_{\ell}^{\kappa g}$ is given by

[TABLE]

where

[TABLE]

$b_{g}$ is the galaxy clustering bias, and $f(z)$ is the redshift distribution. The uncertainty on $C_{\ell}^{\kappa g}$ is given by

[TABLE]

where $N_{\ell}^{\kappa\kappa}=[\ell(\ell+1)/2]^{2}N_{\ell}^{\phi\phi}$, $C_{\ell}^{gg}$ is the galaxy angular auto-power spectrum given by

[TABLE]

$N^{gg}=1/\overline{n}_{g}$ is the galaxy shot noise related to the areal number density $\overline{n}_{g}$, $f_{\rm sky}^{\kappa g}$ is the sky fraction observed by both the LIM lensing and galaxy surveys, and we take the limit $C_{\ell}^{\kappa\kappa}\ll N_{\ell}^{\kappa\kappa}$. As in F18, we assume the specifications from the Large Synoptic Survey Telescope (LSST) (LSST Dark Energy Science Collaboration, 2012) which are listed in Sec. 4.4 of F18. As in the case for $C_{\ell}^{\phi\phi}$, we find that the uncertainty for $C_{\ell}^{\kappa g}$ in the SFB formalism decreases significantly at the same scales. In Table 3, we list the SNR values of the $\kappa g$ autocorrelation for the Plane-Parallel case and our SFB case. We find an increase of only 10%, mainly since the SNR is dominated by small scales where the PP noise prediction is already in parity with the signal $C_{\ell}^{\phi\phi}$ such that reducing the noise at large scales doesn’t help much. However, for lower redshift bins the increase is as high as 33% and all the redshift bins should have higher increases when you just consider angular scales $\ell<100$ which may be relevant for certain cosmological applications.

We expect the difference in these estimators to affect the science yield for upcoming 21-cm lensing surveys. However, this will be limited by wide-angle systematics in 21-cm intensity maps. CMB lensing maps from the Planck satellite are currently limited to $\ell>8$ due to systematics, and we expect this issue to be worse for 21-cm intensity maps. However, it may be possible for a futuristic survey with a significant increase in the number of detectors can achieve a high enough angular resolution to mitigate these systematic effects. Even if this is not the case, it would only reduce the utility of the SFB full sky estimator since the SFB flat-sky predictions differ from PP at even fairly small angular scales $\ell\sim 100$.

One application that can benefit from a full-sky estimator is constraints on modified gravity through measurements of $E_{G}$ (Pullen et al., 2015) using 21-cm lensing (Pourtsidou, 2016). In particular, modes $\ell\lesssim 100$ could be relevant in these constraints and thus it would be beneficial to consider $E_{G}$ constraints with these new estimators. Also, this estimator is not only applicable to 21-cm lensing, but any emission map lensing that can be mapped over large areas and multiple redshifts. This estimator will not readily be useful for upcoming intensity maps of galactic emission lines such as CO, CII, or Ly$\alpha$ since they tend to map over small areas. However, it may be beneficial for some science cases to sacrifice sensitivity to map a large area, in which case a large-scale lensing survey for these lines may be feasible, particularly if a futuristic survey could map larger areas while maintaining current sensitivities.

Finally, while we only construct an estimator for a 2D lensing field $\phi_{\ell m}$, this work could be readily extended to probe $\phi_{\ell m}(k)$, the SFB transform of the 3D lensing field as a function of the radial distance to the source. 3D lensing in SFB coordinates have already been considered and measured for galaxy lensing surveys (Castro et al., 2005; Kitching et al., 2014), and in principle this could also be applied to LIM lensing as well. A major advantage of this approach would be that while $C_{\ell}^{\kappa g}$, the angular lensing-convergence-galaxy cross-power spectrum, is a projection of the cross-power spectrum $P^{\nabla^{2}\Phi g}(k)$, the SFB cross-power spectrum $C_{\ell}^{\kappa g}(k,k^{\prime})$ is directly proportional to $P^{\nabla^{2}\Phi g}(k)$. This would especially be advantageous for measurements of $E_{G}$, in that the $E_{G}$ estimator would no longer require line-of-sight integrals over the numerator and denominator, eliminating systematics associated with the redshift distribution as well as naturally probing scale-dependent growth. We leave this for future work.

6 Conclusions

In this paper, we considered full sky effects on LIM lensing, particularly considering HIRAX, an upcoming survey for the 21 cm line. We have conducted our analysis using the spherical Fourier-Bessel series expansion of three-dimensional functions. Isotropy broken by weak lensing induces correlations which can be used to reconstruct the lensing potential by using a quadratic estimator. For large scales in particular, curved geometry becomes relevant through the extra correlations due to the three dimensional nature of the LIM signal, which is an additional benefit on top of the use of information at multiple redshifts.

In our forecasts of lensing reconstruction noise for the predicted 21 cm power spectra from HIRAX we include both our SFB Full-Sky and Flat-Sky estimators, in which we notice a significant drop in lensing noise for large scales ($\ell\lesssim 100$), scales past which we also notice a deviation between the SFB and PP lensing noise, which is suggestive of the key role of geometry in this effect. In addition, in our forecasts of $C_{\ell}^{\kappa g}$, where we assume galaxy measurements from the future LSST survey, we notice a similar trend at roughly the same angular scales. Conversely, at sufficiently small angular scales ($\ell\geqslant 100$), the difference between flat and full sky estimators is negligible, as expected. The corresponding SNR reflect this difference from the Plane-Parallel treatment of 3D weak lensing. Indeed, we find an increase in SNR by 107% for $C_{\ell}^{\phi\phi}$ and 10% for $C_{\ell}^{\kappa g}$, both in relation to the corresponding Plane-Parallel case.

Considering our formula for the angular scale at which the SFB and PP predictions deviate, $L_{\rm dev}=\ell_{\rm res}-k_{\rm eq}R$, we expect this effect to be relavant mainly for low-redshift LIM surveys. For high-redshift surveys, not only does $k_{\rm eq}R$ increase, but $\ell_{\rm res}$ decreases since the angular resolution is lower for larger observed wavelengths. Both of these effects would lower $L_{\rm dev}$. Another effect is that $L_{\rm dev}$ is dependent on cosmic parameters through both $k_{\rm eq}$ and $R(z)$. Whether $L_{\rm dev}$ could be measured with enough precision to serve as a standard ruler we will leave as an exercise for the reader.

Indeed, properly accounting for the spherical geometry shows promising benefits at the relevant scales. In particular, such a treatment is useful when one considers effects which dominate at large scales. Further, one may also consider full sky effects on reconstruction through polarization of the 21 cm line as well, which would be expected to demonstrate similar trends across large scales.

Further, we note that we have accounted for the fact that the survey is over a thin shell by applying an approximate boundary condition at the minimum radius of the survey volume. While this has worked reasonably well, in order to properly do a spherical-shell calculation of the lensing potential one would have to use Bessel functions of both first and second type as was done for the galaxy power spectrum in Samushia (2019), which naturally account for the fewer modes in a spherical shell. Finally, one should note that we have approximated $C_{\ell nn^{\prime}}\approx C_{\ell n}\delta_{nn^{\prime}}$ in our forecasts, which is reasonable for the non-RSD case. However, this is not strictly the case if one includes RSD (particularly for $k_{\ell n}R\sim\ell$). Properly factoring RSD effects would thus require using the full $C_{\ell nn^{\prime}}$. Regardless, we expect that the overall trends of our results should be immune to this detail. We leave a full investigation of these effects for future work.

Acknowledgements

We wish to thank S. Foreman, T. Kitching, B. Leidstedt, A. Liu, and L. Samushia for helpful comments and discussions. We also thank E. Schaan for reviewing an earlier manuscript, and we also thank the referee for their comments. AP was supported by NASA ROSES NNH17ZDA001N-ATP Grant No. 80NSSC18K10148. PC was supported by the Dean’s Undergraduate Research Fund at New York University.

Appendix A Quadratic Estimator for Lensing in the Flat-Sky Limit

Here, we derive the quadratic estimator of the lensing potential in the flat-sky limit, where we require to move from the SFB moments of the 21 cm signal to its Fourier moments in the limit of small angles. We use the following approximation from Okamoto & Hu (2003) for a function defined on a sphere in order to relate it with its harmonic moments on the sphere

[TABLE]

We apply this dictionary to Eq. 67 to approximate the flat-sky estimator for the lensing potential as

[TABLE]

We require an additional approximation in the limit $\ell,\ell_{1},\ell_{2}\gg 1$ (Okamoto & Hu, 2003)

[TABLE]

Thus, we extract the geometric factors in $f_{\ell_{1}\ell\ell_{2}}^{n_{1}n_{2}}$ and $g_{\ell_{1}\ell_{2}\ell}^{n_{1}n_{2}}$ and make the following redefinitions in terms of corresponding barred quantities

[TABLE]

where

[TABLE]

We will specify the approximations for $\bar{g}_{\ell_{1}\ell_{2}\ell}^{n_{1}n_{2}}$ later since they require additional machinery.

We use the formula for the integral of three spherical harmonics,

[TABLE]

and make use of the approximation (Okamoto & Hu, 2003)

[TABLE]

This allows us to simplify Eq (LABEL:eq:phiell1) to

[TABLE]

where $\bm{\ell_{2}}\equiv\bm{\ell}-\bm{\ell_{1}}$ and the same for the primed case. Following the same procedure as for the full-sky case, we take the covariance of the lensing estimator, which is given as

[TABLE]

Thus, the variance is

[TABLE]

We are required to approximate $A(\ell)$ as well, which we do below. We expand Eq. 69 to obtain

[TABLE]

We use the formula for the integral of three spherical harmonics in Eq. 112, which in the limit $\ell,\ell_{1},\ell_{2}\gg 1$ produces

[TABLE]

Further, note that

[TABLE]

where we have used Eq. 113. Using this, we can obtain the result for $A(\ell)$

[TABLE]

Similarly, we arrive at the flat-sky expression for $\bar{g}^{n_{1}n_{2}}_{\ell_{1}\ell_{2}\ell}$

[TABLE]

Finally, assuming a fiducial survey allows us to obtain the noise power spectrum

[TABLE]

Appendix B Intensity Mapping Instrumental Effects

B.1 Angular and Radial Window Functions

Just like for the CMB, a line intensity map signal will be convolved by the instrument. Yet while a CMB map is only convolved by the beam in two dimensions, an intensity map will also be convolved along the line-of-sight due to the spectral window in frequency space. In this section, we derive the angular (beam) and radial (spectral) window functions that can be applied to the LIM signal. Throughout this derivation, any integrals over $r$ will be understood to be over the range $r_{\rm min}<r<r_{\rm max}$.

We start by defining a window function in real space $D(\mathbf{r}^{\prime}|\mathbf{r})$ such that the observed map $T^{\rm obs}(\mathbf{r})$ is a convolution of the true map $T(\mathbf{r})$ with the window function, assuming no instrumental noise or masking, according to

[TABLE]

Next, we assume that the angular and radial window functions are separable such that $D(\mathbf{r}^{\prime}|\mathbf{r})=A(r^{\prime}|r)B(\mathbf{\hat{r}}^{\prime}|\mathbf{\hat{r}})$ where $A(r^{\prime}|r)$ and $B(\mathbf{\hat{r}}^{\prime}|\mathbf{\hat{r}})$ are the radial and angular window functions, respectively. This allows us to rewrite the convolution as

[TABLE]

It can be shown that the SFB moments of the observed map can be written in the form

[TABLE]

By taking the covariance of these moments, we find

[TABLE]

where

[TABLE]

and

[TABLE]

A simpler expression for the mode covariance arises if it is assumed that $A(r^{\prime}|r)$ and $B(\mathbf{\hat{r}}^{\prime}|\mathbf{\hat{r}})$ are both highly peaked, in which case $W_{\ell nn^{\prime}}^{A}$ and $W_{\ell}^{B}$ should be independent of pixel locations $\mathbf{r}_{1}$ and $\mathbf{r}_{2}$. In this case, the integrals in Eq. 132 can be performed, giving us

[TABLE]

For the analytical formulae, this representation for the radial and angular harmonic window functions will be retained. In Appendix B.2 we derive that the radial harmonic window function can be written as a sinc function in the form

[TABLE]

where

[TABLE]

where $\Delta r$ is the full width of the radial window.

B.2 Spectral Window Function

Here we derive the radial harmonic window function $W_{\ell nn^{\prime}}^{A}=A_{\ell}(k_{\ell n})A_{\ell}(k_{\ell n^{\prime}})$, where

[TABLE]

and $A(r^{\prime}|r)$ is given by

[TABLE]

where $\Delta r$ is the full width of the radial window. This is similar to the window assumed in L17, except we add an extra pre-factor to make it normalize conveniently. Technically the integration over $r$ should be over the range $r_{\rm min}<r<r_{\rm max}$; however, since the window function A(r’|r) is expected to be relatively narrow, we can extend the integration range to $r>0$ without losing accuracy. We are mainly going to consider the case where the modes are more radial than angular, $\ell\ll kr$. However, for the case where $\ell\gtrsim kr$, $j_{\ell}(kr)$ is a smooth, monotonic function that can be taken out of the integral compared to the highly-peaked $A(r^{\prime}|r)$, which itself is normalized. Thus, for $\ell\gtrsim kr$, $A_{\ell}(k)\simeq 1$.

For $\ell\ll kr$, we can use the approximation

[TABLE]

allowing us to write the numerator of Eq. 138, which we label $Q$, as

[TABLE]

Since $A(r^{\prime}|r)$ is cut off for $r^{\prime}>r+\Delta r/2$, it is fine to extend the integral out to infinity. We evaluate this integral using a similar trick as was performed in L17. First, we convert the sine function to an exponential with a derivative.

[TABLE]

Then, we define a new function $C(r^{\prime})=A(r^{\prime}+r|r)$ which can be written as

[TABLE]

Given this, we find

[TABLE]

where we define $x\equiv r^{\prime}-r$ and we can extend the integral out to $-\infty$ since $C(x<-\Delta r/2)=0$, allowing to then define the Fourier transform $\tilde{C}(k)$. Next, we make the argument that $C(x)$ is approximately an even function in the limit where $r^{\prime}\ll r$. Under this approximation, we can apply the Re operator to just the exponential, giving us

[TABLE]

Inserting this into Eq. 138, we can write

[TABLE]

Now we evaluate $\tilde{C}(k)$. Using the $r^{\prime}<<r$ approximation, which allows us to just evaluate the real part of the Fourier transform, we find

[TABLE]

From this, we know $\partial\tilde{C}/\partial k\propto\Delta r/r\ll 1$ such that we can approximate $A_{\ell}(k)\simeq r^{2}\tilde{C}(k)$, which gives us Eq. 137.

B.3 Instrumental Noise

Here we derive the instrumental noise contribution to the covariance of the SFB moments. In this analysis we will set each three-dimensional pixel, or voxel, to be labeled by an index such that the noise fluctuations are labeled $T_{i}^{N}$. We assume that the noise fluctuations are uncorrelated between voxels and that their distribution varies only with frequency, and thus the radial coordinate for a given line. This gives us a noise covariance matrix $\left\langle T_{i}^{N}T_{j}^{N}\right\rangle=\sigma_{N,i}^{2}\delta_{ij}$.

We begin by writing the SFB moment for the noise fluctuations as a sum over voxels

[TABLE]

where $\Delta_{i}$ is the volume of voxel $i$, given by $\Delta_{i}=\Delta\Omega\Delta r_{i}r_{i}^{2}$. Taking the covariance, we find

[TABLE]

Converting back from a sum to an integral, we can write the covariance, after integrating over angle, as

[TABLE]

where

[TABLE]

and $\Delta(r)$ is the size of the voxel as a function or $r$. The voxel depth is determined by the frequency bin size, such that we can write

[TABLE]

where

[TABLE]

Using the results from the previous section, we can now write

[TABLE]

where we define $W_{\ell n_{1}n_{2}}^{AB}\equiv W_{\ell n_{1}n_{2}}^{A}W_{\ell}^{B}$.

Appendix C Proof of Form of Quadratic Estimator

Here we prove that our lensing estimator is the most general quadratic estimator that can be constructed for this problem. Consider a general quadratic estimator for $\phi(\hat{n})$

[TABLE]

We may expand $\hat{\phi}$ and $\widetilde{T}$ in spherical harmonics and SFB series respectively to obtain

[TABLE]

Now, $G(\vec{x},\vec{y},\hat{n})$ depends on three directions, thus it may be expanded in tripolar spherical harmonics as

[TABLE]

One may integrate the angular part of Eq. 160 by using the orthogonality relations of spherical harmonics to obtain

[TABLE]

Now, one may make the re-definition,

[TABLE]

which leaves the estimator in the form

[TABLE]

But, $g_{\ell\ell^{\prime}L}^{mm^{\prime}M,nn^{\prime}}$ is to be determined by minimizing $\langle\hat{\phi}_{LM}\hat{\phi}^{*}_{L^{\prime}M^{\prime}}\rangle$ which we assume to be diagonal, i.e., $\propto\delta_{LL^{\prime}}\delta_{MM^{\prime}}$. This can only be achieved due to the orthogonality of the 3j-symbol, so $g$ cannot depend on $\{m,m^{\prime},M\}$. So we now have

[TABLE]

Finally, since we need an unbiased estimator, we need a suitable normalization which we denote as $A_{L}$. Thus we have the final form of the estimator as given in the main text.

Appendix D Approximation for Spherical Bessel Function at large $\ell$

Recall an integral representation of $j_{\ell}(x)$

[TABLE]

It is known that performing a saddle-point approximation (along with a suitable substitution) produces the following approximation for $j_{\ell}(x)$ for the region $x\leq\ell$ ($\ell\gg 1$)

[TABLE]

Similarly, one may use an integral representation for the spherical Hankel function $h_{\ell}^{(1)}(x)$ to saddle-point approximate it in the region $x\geq\ell$ ($\ell\gg 1$)

[TABLE]

where the real part of this expression gives the approximation for $j_{\ell}(x)$. We may consider the regime $kr\geq\ell$ in particular with $x\equiv kr$. Firstly, note,

[TABLE]

Now we may define $k_{||}^{2}=k^{2}-\ell^{2}/r_{1}^{2}=k^{2}-k_{\perp}^{2}$. So we have

[TABLE]

since $a^{2}=1+\frac{1}{N^{2}}$ where $N=\frac{R}{\Delta R}>>1$. In light of this, we have the following,

[TABLE]

and

[TABLE]

The correction term above is small since $k_{\perp}/k<1$ and $N>>1$, so we have

[TABLE]

In the $k_{\perp}-k_{||}$ plane we may now make the following redefinition $k_{\perp}+ik_{||}=ike^{-i\varphi_{k}}$. This makes the above expression simply $e^{i\varphi_{k}\ell-\frac{\pi\ell}{2}}$. Thus, the final expression for the Bessel function becomes

[TABLE]

Note that in the limit $kr>>\ell\,(k_{||}>>k_{\perp})$, $\varphi_{k}\approx 0$ and $k\approx k_{||}$, which recovers the well known approximation for the Bessel function in this regime.

Appendix E Foregrounds

Foreground contamination in the flat-sky case is rather straightforward, owing to the fact that the foreground power spectrum has a simple expression (L17) which can be derived straightforwardly for a $T(\vec{r})=T(\mathbf{\hat{r}})$ given so,

[TABLE]

for a survey at mean comoving radius $R$ and survey width $\Delta R$, particularly for a tophat window function (which is the case we are concerning ourselves with in the SFB formalism). In particular, this expression is independent of $k_{\perp}$ since $C^{\text{fg}}_{\ell}$ is actually constant, i.e., one only needs to impose a cutoff $k_{||,\text{min}}$. Below, we will derive an expression to determine the appropriate methods for foreground cutoffs for SFB, which has a nontrivial expression for foreground power. For instance, it was shown in L17 that the foreground cuts for SFB will be scale-dependent (which is not the case for Plane-Parallel). Here, we will show that the correct cuts in SFB coordinates are given in terms of $k_{||,\rm min}$ by $k^{2}_{\rm min}=k_{||,\rm min}^{2}+\ell^{2}/r_{\rm max}^{2}$. Note that any radial integrals will be considered to be from $r_{\rm min}$ to $r_{\rm max}$, but we will not explicitly write the limits.

To start, consider the foreground signal in SFB (L17),

[TABLE]

In particular, $\tau_{\ell n}C_{\ell n}^{\text{fg}}$ is some nontrivial function which is expected to produce a nontrivial $k_{\ell n,\text{min}}$ corresponding to the plane-parallel $k_{||,\text{min}}$

In particular we wish to examine $\tau_{\ell n}C_{\ell n}^{\text{fg}}$ in the regime $\ell>>1$ and $k_{\ell n}r\gtrsim\ell$ in order to obtain the correct $k_{\text{min}}$. Using the approximation for $j_{\ell}(kr)$ obtained in Appendix D we have

[TABLE]

where $k_{||,\ell n}^{2}\equiv k_{\ell n}^{2}-\ell^{2}/r_{\rm max}^{2}$ and $\zeta=\varphi_{k}\ell-\pi\ell/2$. Now since we have $k_{\ell n}r\approx\ell$, we have $k_{||,\ell n}\approx 0$ (to order 1/N). Thus, $\varphi_{k}\approx\pi/2$ so $\zeta\approx 0$. So we have

[TABLE]

However, note that in the limit $k_{||}\rightarrow 0$, particularly for the case $R\gg\Delta R$, the second term is subdominant. So to leading order, we may say

[TABLE]

Next we need the denominator which is just $W_{\ell\ell}(k,k)$, i.e.,

[TABLE]

since $k_{||}\approx 0$. To complete, we have

[TABLE]

However, we may average out the oscillatory part of this expression as in L17, which gives exactly that $\tau_{\ell n}C_{\ell n}\approx R^{2}\Delta RC_{\ell}(k_{||,\ell n})$ when $k_{\ell n}R\gtrsim\ell$. This analysis indicates that the correct contours should be $k_{\text{min}}^{2}=k_{||,\text{min}}^{2}+\frac{\ell^{2}}{r_{\rm max}^{2}}$, as shown in Fig. 9.