# Gravitational lensing of galaxy clustering

Brandon Buncher  Department of Physics, University of Illinois Urbana-
Champaign  Gilbert Holder  Department of Physics, University of Illinois
Urbana-Champaign  Selim C. Hotinli  Perimeter Institute for Theoretical
Physics

###### Abstract

We investigate lensing reconstruction using the clustered galaxy distribution
as a source field, using both the traditional cosmic microwave background
quadratic estimator and a shear-only estimator. We calculate the expected
signal-to-noise ratio of the cross power spectrum of such reconstructions with
cosmic shear measurements for an LSST-like galaxy survey. Modeling the galaxy
field as a Gaussian random field, we find that there is substantial clustering
signal in the source field at angular scales substantially smaller than those
typically used by CMB reconstructions. The expected signal-to-noise for cross-
correlations in LSST from cosmic shear is $\sim$60 in the presence of shape
noise, while cross-correlating with a sample-variance limited mass map would
have signal-to-noise in the hundreds. This type of cross-correlation could be
used as a way to identify systematic errors in lensing studies and is just one
example of many possible higher order correlations in galaxy surveys that may
contain substantial cosmological information.

]}

##  1 Introduction

There is large scale structure (LSS) in the universe, as evidenced by the
spatial clustering of galaxies (Peebles, 1980). This LSS deflects light as it
propagates through the inhomogeneous universe (Bartelmann & Schneider, 2001).
These deflections can occasionally lead to spectacular Einstein rings and
multiple images, but usually will lead to only small distortions to the shapes
and brightnesses of background objects. This has been measured in galaxy
surveys through the changes in the detailed shapes of galaxies (Luo et al.,
2018; Gatti et al., 2021), and has been measured in the distortion to the
characteristic clustering patterns in the cosmic microwave background (Aghanim
et al., 2020; Omori et al., 2023; Madhavacheril et al., 2023; Ade et al.,
2015; BICEP/Keck Collaboration et al., 2022).

In this paper we investigate the gravitational lensing of the clustering of
background galaxies by foreground matter, analogous to the case of
gravitational lensing of the CMB. A related method was recently discussed in
detail by Nistane et al. (2022) using an estimator that combined this
information with information about the magnification bias to increase the
signal-to-noise ratio (SNR). We determine the viability of a different (shear-
based) estimator and investigate the prospects for measuring this effect in
cross-correlation with cosmic shear measurements.

In the case of CMB lensing, the workhorse tool that is widely used is a
quadratic estimator Okamoto & Hu (2003), which amounts to multiplying together
two CMB maps that each have their own filtering. The filtering of the maps is
chosen to create an unbiased, minimum variance estimator, which optimizes the
SNR of the resulting lens reconstruction. If this lens reconstruction is then
correlated with a galaxy-based LSS tracer, the result could be expressed as a
3-pt function, as it involves products of two CMB maps and an LSS map. As we
will show below, for the similar case of lensing of the galaxy clustering
distribution by foreground matter, the reconstruction of the foreground
lensing potential will involve two versions of the background galaxy number
count field, and then cross-correlation with an LSS tracer such as cosmic
shear can be expressed as a galaxy-galaxy-shear 3-pt function. If we instead
used the foreground galaxy counts, then there is a 3-pt function that would be
quadratic in the background galaxy number counts and linear in foreground
galaxy number counts. In general, the three fields going into these estimates
can be from diverse LSS tracers: the background 2-pt function could be cross-
correlations between different LSS tracers, and the foreground could be taken
to be any LSS tracer of the foreground matter. In this work, we only look at
the simplest case of using the auto-correlation of galaxies in the background
being lensed by foreground matter in a way that leads to a measured shear of
the same background material, mainly to keep the accounting simple.

This work is closely related to other work that has looked at the lensing of
the clustering of LSS tracers. Prior work on lensing reconstruction has
primarily used the CMB as a light source (Okamoto & Hu, 2003). However,
modifications of the quadratic estimator have been applied to non-CMB sources,
particularly spectral line intensity maps (LIM). The closest works to ours use
the standard quadratic estimator to perform lensing reconstructions from
measurements of the Lyman-$\alpha$ forest (Metcalf et al., 2020), the cosmic
infrared background (CIB) (Schaan et al., 2018; Foreman et al., 2018), and
21cm lines (Pen, 2004; Zahn & Zaldarriaga, 2006; Lu & Pen, 2008; Agustin
Lozano Torres & Schaefer, 2022; Pourtsidou & Metcalf, 2014, 2015; Chakraborty
& Pullen, 2019); other 21cm reconstruction methods include (Lu et al., 2010;
Jalilvand et al., 2019; Zhang & Pen, 2006; Hilbert et al., 2007; Metcalf &
White, 2007, 2009).

An important point of focus in each of these studies has been the non-
Gaussianity of the background field being lensed. It has been shown (Schaan et
al., 2018; Foreman et al., 2018; Feng & Holder, 2019) that there can be
substantial biases in the recovered lensing signal if the non-Gaussianity is
ignored. Nevertheless, for this preliminary investigation we will disregard
the non-Gaussianity of the field and its uncertainties in what follows, so the
results presented below will no doubt require future work to study these
effects. These corrections arise because non-linear gravitational collapse
leads to an intrinsic background of spatial variation in the statistics of
small-scale clustering, the magnitude of which depends on factors including
the spatial scales probed, redshift, and survey parameters (Foreman et al.,
2018). However, some lensing estimators that are optimal for a Gaussian source
field have been shown to be resistant to non-Gaussian noise. Foreman et al.
(2018) constructed a bias-hardened estimator that included contributions from
the density field trispectrum that corrected for mild gravitational
nonlinearities, while Lu et al. (2010) calibrated their estimator with N-body
simulations to create a set of lensing weights that are optimal for non-
Gaussian fields.

In this work, we closely follow the procedures described by Schaan & Ferraro
(2019), who present a method for lensing reconstruction that is less affected
by foreground contamination and non-Gaussianities. In Schaan & Ferraro (2019),
the authors construct a convergence field hybrid quadratic estimator for a
Gaussian CMB source field. The hybrid estimator is comprised of an inverse-
variance weighted linear combination of a magnification estimator and a shear
estimator. It used the shear-only component for small scales and both the
convergence and shear components for large scales. The authors found that the
hybrid estimator outperformed the SNR of the standard quadratic estimator and
that the shear component was less affected by foreground contamination and
non-Gaussianities (Schaan & Ferraro, 2019).

Lu & Pen (2008) present a similar hybrid estimator for a 21 cm source field.
The authors combined the magnification and shear estimators at all spatial
scales. They show that, for a Gaussian source field, their estimator is
equivalent to the 3D optimal quadratic deflection estimator introduced by Zahn
& Zaldarriaga (2006), but was less affected by gravitational nonlinearities
when reconstructions were performed using a non-Gaussian source field. Like
Schaan & Ferraro (2019), the authors found that the hybrid estimators were
also less susceptible to foreground contamination.

We have found one other work that uses galaxy number counts as a source field
for weak lensing reconstruction in a similar manner. In (Nistane et al.,
2022), the authors demonstrate that measurements of the magnification bias can
be combined with the information from the background galaxy density field to
reconstruct the convergence power spectrum using a quadratic estimator at a
high SNR. The magnification bias has been well-studied both through
measurements and forecasts (Scranton et al., 2005; Ponce et al., 2013; Camera
et al., 2014; von Wietersheim Kramsta et al., 2021; Maartens et al., 2021);
however, for the purposes of this work, we will ignore its effects in order to
focus on the particular information that is carried in the distortion of the
clustering field. The effects of galaxy clustering and magnification bias on
weak lensing reconstruction were also examined in (Ma et al., 2023). Noise
from intrinsic galaxy clustering is much larger than the lensing signal, so
the authors developed a multi-band linear estimator for the convergence field
that minimizes the noise introduced by intrinsic clustering. While this work
does not use galaxy number counts as a source field, it supports the
conclusions of (Nistane et al., 2022) that magnification bias provides
valuable information for weak lensing reconstruction.

We perform projections of the SNR for a survey with parameters consistent with
the Legacy Survey of Space and Time at the Vera C. Rubin Observatory (LSST)
(Ivezić et al., 2019). To estimate the expected properties of LSST galaxies we
will use the CosmoDC2 simulations Korytov et al. (2019); for consistency, in
our calculations we will also use the cosmology assumed for that simulation.
This parameter set is similar to those of the best fit WMAP-7 parameter set
(Komatsu et al., 2011): $\omega_{\textnormal{cdm}}=0.1109$,
$\omega_{\textnormal{b}}=0.02258$, $n_{s}=0.963$, $h=0.71$, $\sigma_{8}=0.8$,
and $w=-1.0$. We investigate $r$-band magnitude cuts of $m_{r}<23$ and
$m_{r}<25$, chosen to correspond roughly with the Dark Energy Survey (DES)
limiting magnitude (Abbott et al., 2022) and the 5$\sigma$ single visit
limiting magnitude of the LSST survey (LSST Science Collaborations and LSST
Project, 2009), respectively. The LSST 5$\sigma$ 10 year coadd. dataset would
increase the limiting magnitude to around $m_{r}<27$. However, we find that
information from very faint galaxies was not needed to obtain a high SNR using
our method. We expect that the full LSST 5$\sigma$ 10 year coadd. dataset
would increase the SNR.

The next section lays out the theoretical framework for our calculations, and
the following section describes the computational setup. We then present the
expected SNR for an LSST-like survey and a DES-like survey, and then close
with a discussion.

##  2 Theory

In this section we describe the theoretical framework, including a brief
overview of gravitational lensing, a discussion of lensing estimation using
quadratic estimators, and a method for calculating the expected lensing power
spectrum for a screen of sources at a given redshift.

###  2.1 Lensing basics

We can observe the galaxy number density $n$ at redshift $z$ and locations
$\vec{\theta}$ to be

$n(\vec{\theta})=n^{u}(\vec{\theta}+\vec{\alpha}),$ (1)

where $\vec{\alpha}$ is the gravitational deflection and the superscript $u$
indicates the field that would be observed in the absence of gravitational
lensing (the unlensed field). We here ignore the effects of magnification
bias, where selection effects for the galaxy sample are affected by the
lensing magnification.

In the weak lensing regime (where the deflection $\alpha$ is small), the
lensed number density can be approximated via a Taylor series:

$n(\vec{\theta})\approx n^{u}(\vec{\theta})+\nabla
n^{u}(\vec{\theta})\cdot\vec{\alpha}.$ (2)

The deflection angle can be expressed as the gradient of a deflection
potential, $\vec{\alpha}=\nabla\psi$. The deflection potential $\psi$ is the
line-of-sight projection of the Newtonian gravitational potential $\Phi$
(Bartelmann & Schneider, 2001):

$\displaystyle\psi=\frac{2}{c^{2}}\frac{D_{\textnormal{ds}}}{D_{\textnormal{d}}D_{\textnormal{s}}}\int\Phi\left(D_{\textnormal{d}}\,\vec{\theta},z\right)dz,$

where $D_{\textnormal{d}}$ is the distance between the observer and the lens,
$D_{\textnormal{s}}$ is the distance between the observer and the source, and
$D_{\textnormal{ds}}$ is the distance between the source and the lens.

The effect of weak lensing is either to expand or contract the apparent solid
angle of a galaxy or cluster, known as “convergence,” or to systematically
stretch or squeeze them in orthogonal directions, known as “shear.” This
affects the perceived clustering of galaxies, leading to alterations to the
observed galaxy power spectrum. The convergence field $\kappa$ and two
components of the shear field $\gamma_{+}$ and $\gamma_{\times}$ can be
expressed as a transformation of the deflection potential $\psi$:

$\displaystyle\kappa=-\frac{1}{2}\left(\partial_{x}^{2}+\partial_{y}^{2}\right)\psi=-\frac{1}{2}\nabla^{2}\psi,$
$\displaystyle\gamma_{+}=-\frac{1}{2}\left(\partial_{x}^{2}-\partial_{y}^{2}\right)\psi,$
(3) $\displaystyle\gamma_{\times}=-\partial_{x}\partial_{y}\psi$

for some orthogonal angular directions $\hat{x}$ and $\hat{y}$ where here we
ignore sky curvature and approximate the sky as being a Cartesian grid.

A value of $\kappa>0$ represents image dilation (an increase in solid angle of
the object), while $\kappa<0$ represents contraction. A value of
$\gamma_{+}>0$ indicates stretching in the horizontal ($\hat{x}$) direction,
while $\gamma_{+}<0$ represents stretching in the vertical ($\hat{y}$)
direction. A value of $\gamma_{\times}>0$ represents stretching along the line
$y=-x$, while $\gamma_{\times}<0$ represents stretching along the line $y=x$.
Because each of the lensing observables are constructed from derivatives of
the deflection potential, this means that they may be reconstructed from one
another.

###  2.2 Lensing reconstruction with quadratic estimators

The foreground lensing deflection can be reconstructed using a “quadratic
estimator”; multiplying Eq. (2) by $\nabla n^{u}(\vec{\theta})$ leaves the
second term proportional to the square of the galaxy density field gradient
multiplied by the deflection (Okamoto & Hu, 2003). This is known as a
quadratic estimator because it uses the square of the observed field to
reconstruct the lensing deflection. An optimal quadratic estimator (i.e. a
quadratic estimator that minimizes the noise in the reconstructed field) can
be constructed that takes into account the noise in the background field; a
derivation of this estimator can be found in the appendix of many papers,
including in Schaan & Ferraro (2019), which we follow closely in this work.

On large scales, the main effect of the lensing distortion is to simply change
the observed small-scale clustering pattern of the background field. The
lensing effect therefore breaks statistical isotropy: the small-scale power
spectrum of the background field will be observed to vary from place to place
on the sky because of the variations in the large-scale lensing effects.

In the limit of large-scale lensing of small-scale structure, the optimal
quadratic estimator can be shown to be a mix of a convergence estimator and a
shear estimator (Bucher et al., 2012). On small scales, this separation
between the estimators breaks down, but it has been shown that the shear
estimator is more robust to non-Gaussian foregrounds and experimental non-
idealities (Lu & Pen, 2008; Schaan & Ferraro, 2019).

The lensing maps are reconstructed using two combinations of the background
galaxy density field along with a Fourier space lensing kernel
$g_{{\bm{{L}}},{\bm{{\ell}}}}$ that captures the lensing-related correlations.
We use a notation where ${\bm{{L}}}$ refers to an angular mode in the
foreground lensing potential, while ${\bm{{\ell}}}$ refers to angular modes in
the background source that is being lensed.

The standard quadratic estimator is a minimum variance quadratic estimator
that is commonly used for reconstruction of the CMB lensing deflection
potential $\psi$. The QE is calculated using the Fourier transform (spherical
harmonics decomposition) of the background field. Throughout this paper, we
use the convention that $\bm{{\ell}}$ indexes the multipole of the Fourier
transform of the background field and $\bm{{L}}$ indexes the multipole of the
Fourier transform of the lensing field. Traditionally, the standard quadratic
estimator uses the background CMB temperature field $\Theta$ and its power
spectrum $C_{\ell}^{\Theta}$ to calculate the deflection field; in our case,
we replace the observed temperature field $\hat{\Theta}$ with the observed
galaxy number density $n$. The standard quadratic estimator performs the
reconstruction by calculating correlations between the observed number density
field at different locations $n_{\bm{{\ell}}}$ and
$n_{{\bm{{\ell}}}-{\bm{{L}}}}$; this convolution is weighted by the lensing
kernel $g_{{\bm{{\ell}}},{\bm{{L}}}}$. The reconstruction of the observed
lensing deflection potential in Fourier space $\psi({\bm{{L}}})$ using the
galaxy number density field is calculated (Lewis & Challinor, 2006)

$\displaystyle\psi({\bm{{L}}})=N_{\bm{{L}}}\int\frac{d^{2}{\bm{{\ell}}}}{2\pi}n_{\bm{{\ell}}}n^{*}_{{\bm{{\ell}}}-{\bm{{L}}}}g^{\psi}_{{\bm{{\ell}}},{\bm{{L}}}}.$
(4)

Here, the reconstruction noise $N_{\bm{{L}}}$, which is the variance of the
QE, is defined as (Lewis & Challinor, 2006)

$\displaystyle\frac{1}{N_{\bm{{L}}}^{\psi}}=\int\frac{d^{2}{\bm{{\ell}}}}{(2\pi)^{2}}\left[\left({\bm{{L}}}-{\bm{{\ell}}}\right)\cdot{\bm{{L}}}C_{\left|{\bm{{\ell}}}-{\bm{{L}}}\right|}^{0}+{\bm{{\ell}}}\cdot{\bm{{L}}}C_{\ell}^{0}\right]g_{{\bm{{\ell}}},{\bm{{L}}}}^{\psi},$
(5)

where $C_{\bm{{\ell}}}^{0}$ is the power spectrum of $n$ and
$C_{|{\bm{{\ell}}}-{\bm{{L}}}|}^{0}$ is the power spectrum of the shifted
number density field. The lensing kernel $g_{{\bm{{\ell}}},{\bm{{L}}}}$ for
the lensing potential is defined as (Lewis & Challinor, 2006)

$\displaystyle
g_{{\bm{{\ell}}},{\bm{{L}}}}^{\psi}=\frac{({\bm{{L}}}-{\bm{{\ell}}})\cdot{\bm{{L}}}C_{\left|{\bm{{\ell}}}-{\bm{{L}}}\right|}^{0}+{\bm{{\ell}}}\cdot{\bm{{L}}}C_{\ell}^{0}}{2C_{\ell}^{\textnormal{tot}}C_{\left|{\bm{{\ell}}}-{\bm{{L}}}\right|}^{\textnormal{tot}}}$
(6)

where $C_{\bm{{\ell}}}^{0}$ is the power spectrum of the background field
being lensed and $C_{\ell}^{\textnormal{tot}}$ is the same power spectrum plus
noise. For CMB lensing, this noise will be a combination of detector noise and
foregrounds, while for our purposes, this noise will be the shot noise of the
background field.

At large scales, the standard quadratic estimator can be broken down into a
combination of a magnification-only estimator and shear-only estimator. The
estimator is broken down by dividing the standard quadratic estimator weight
function into the contributions from the monopole ($m=0$)
$g^{0}_{{\bm{{L}}},{\bm{{\ell}}}}$ (the contribution from convergence) and
quadrupole ($m=2$) $g^{2}_{{\bm{{L}}},{\bm{{\ell}}}}$ (the contribution from
shear). Higher order terms are neglected in this work.

The lensing kernels for the convergence field are defined as (Bucher et al.,
2012; Schaan & Ferraro, 2019)

$\displaystyle
g_{{\bm{{L}}},{\bm{{\ell}}}}^{\kappa,0}=\frac{C_{\ell}^{0}}{2\left(C_{\ell}^{\textnormal{tot}}\right)^{2}}\frac{\textnormal{d}\ln\ell^{2}C_{\ell}^{0}}{\textnormal{d}\ln\ell}$
(7) $\displaystyle
g_{{\bm{{L}}},{\bm{{\ell}}}}^{\kappa,2}=\cos\left(2\theta_{{\bm{{L}}},{\bm{{\ell}}}}\right)\frac{C_{\ell}^{0}}{2\left(C_{\ell}^{\textnormal{tot}}\right)^{2}}\frac{\textnormal{d}\ln
C_{\ell}^{0}}{\textnormal{d}\ln\ell}.$ (8)

In our case, we will model the noise as being from the shot noise in the
galaxy number counts. In principle, this could be considered to be part of the
signal that would contribute to the convergence estimator (Eq. (7)).

The magnification ($m\,=\,0$) and shear ($m\,=\,2$) estimators of the
convergence field are defined by (Schaan & Ferraro, 2019)

$\displaystyle\hat{\kappa}^{m}_{\bm{{L}}}=\frac{\int\frac{d^{2}{\bm{{\ell}}}}{(2\pi)^{2}}g^{\kappa,m}_{{\bm{{L}}},{\bm{{\ell}}}}n_{\bm{{\ell}}}n_{{\bm{{L}}}-{\bm{{\ell}}}}}{\frac{2{\bm{{L}}}}{L^{2}}\int\frac{d^{2}{\bm{{\ell}}}}{(2\pi)^{2}}g^{\kappa,m}_{{\bm{{L}}},{\bm{{\ell}}}}\left[{\bm{{\ell}}}C_{\ell}^{0}+({\bm{{L}}}-{\bm{{\ell}}})C_{\left|{\bm{{L}}}-{\bm{{\ell}}}\right|}^{0}\right]},$
(9)

where $n_{\bm{{\ell}}}$ is the Fourier transform of the galaxy number density
field. The reconstruction noise $N_{\bm{{L}}}^{m}$ is the variance of each
estimator:

$\displaystyle
N_{\bm{{L}}}^{\kappa,m}=\frac{\int\frac{d^{2}{\bm{{\ell}}}}{(2\pi)^{2}}g^{m}_{{\bm{{L}}},{\bm{{\ell}}}}\left(g^{\kappa,m}_{{\bm{{L}}},{\bm{{\ell}}}}+g^{\kappa,m}_{{\bm{{L}}},{\bm{{L}}}-{\bm{{\ell}}}}\right)C_{\ell}^{\textnormal{tot}}C_{\left|{\bm{{L}}}-{\bm{{\ell}}}\right|}^{\textnormal{tot}}}{\left\\{\frac{2{\bm{{L}}}}{L^{2}}\int\frac{d^{2}{\bm{{\ell}}}}{(2\pi)^{2}}g^{\kappa,m}_{{\bm{{L}}},{\bm{{\ell}}}}\left[{\bm{{\ell}}}C_{\ell}^{0}+\left({\bm{{L}}}-{\bm{{\ell}}}\right)C_{\left|{\bm{{L}}}-{\bm{{\ell}}}\right|}^{0}\right]\right\\}^{2}},$
(10)

where again the $m\,=\,0$ mode is the magnification noise and the $m\,=\,2$
mode is the shape noise.

We will be comparing our lensing reconstructions using the lensing of the
galaxy correlation function to what can be obtained from a cosmic shear survey
that measures the shapes of individual galaxies as a proxy for the
gravitational shear. We thus consider the contribution of shape noise,
$N_{\gamma}$, to lensing reconstructions based on cosmic shear measurements;
$N_{\gamma}$ represents the intrinsic variance of the galaxy ellipticities
that act as noise in the shear measurements for each galaxy in a cosmic shear
survey. As is specified in LSST Science Collaborations and LSST Project
(2009), we take shape noise to be

$\displaystyle N_{\gamma}\equiv\frac{\gamma_{\textnormal{rms}}}{\sqrt{n}},$
(11)

where $\gamma_{\textnormal{rms}}=0.28$ is the shape noise per galaxy expected
for an LSST-like survey and $n$ is the average number of galaxies per
steradian.

###  2.3 The angular power spectra of LSST galaxy subsamples

The quadratic estimators we use require the angular power spectra of each
galaxy subsample, $C_{\ell}$. In a real survey these power spectra will be
measured in the data, but for this forecasting study we need to calculate the
expected power spectra. We use CosmoDC2 mock galaxy catalogs (Korytov et al.,
2019) to guide the inputs for the calculation. The calculation is performed
using the Limber equation (Limber, 1953) assuming a basic linear bias model:

$\displaystyle
C_{\ell}(z)=\frac{H_{0}^{2}}{c^{2}}\int\frac{1}{n^{2}}\left(\frac{d{n}}{d{z}}\right)^{2}b^{2}(z)P_{\delta}\left(\frac{\ell}{\chi},z\right)\,\frac{d\chi}{\chi^{2}},$
(12)

where $\chi$ is the comoving distance, $H_{0}$ is Hubble’s constant,
$\frac{d{n}}{d{z}}$ is the galaxy density distribution, $b(z)$ is the linear
galaxy bias, which represents the ratio of the mean galaxy overdensity to the
mean dark matter overdensity, and $P_{\delta}\left(k,z\right)$ is the total
matter power spectrum at wavenumber $k$ and redshift $z$.

(a) $\frac{d{n}}{d{z}}$

(b) $b(z)$

(c) $C_{\ell}$

Figure 1: (a) $\frac{d{n}}{d{z}}$ for the CosmoDC2 cutout. A fit was performed
using the form in Eq. (13) (solid line) to smooth the function for use in the
calculation of $C_{\ell}$, $C_{L}$, and the SNR. The points used for
performing the fit are shown as dots. (b) $b(z)$ for the CosmoDC2 cutout. (c)
The galaxy angular power spectrum $C_{\ell}$ in $z$-bins of width $\Delta
z=0.1$ for $m_{r}<25$ and $m_{r}<23$. The best-fit $C_{\ell}$ is represented
by solid lines, while the shot noise floor $N_{s}=\frac{1}{n(z)}$
corresponding to each $z$-bin and magnitude cut is indicated by the dotted
lines. The $z$-bins shown here are the ones with the maximum
$\textnormal{SNR}_{\Delta}(L,z)$.

We used a 25 $\textnormal{deg}^{2}$ cutout of CosmoDC2 (Korytov et al., 2019)
to calculate the tomographic galaxy number density distribution
$\frac{d{n}}{d{z}}$ and the linear galaxy bias $b(z)$ . We divided the
CosmoDC2 cutout into 12 redshift bins of equal width with $0<z<3$ and
calculated $n$ within each of these bins, which we used to estimate
$\frac{d{n}}{d{z}}$. To smooth the distribution, we parameterized
$\frac{d{n}}{d{z}}$ as

$\displaystyle\frac{d{n}}{d{z}}=A\,\left(\frac{z}{z_{0}}\right)^{m}\,\exp\left[-\left(z/z_{0}\right)^{\alpha}\right].$
(13)

The form of this fit was chosen to be similar to (Mandelbaum et al., 2018, p.
48); however, because our redshift range was larger than theirs, we introduced
the additional free paramter $m$. We found that for $m_{r}<23$,
$A=16.77\;\textnormal{arcmin}^{-2}$, $z_{0}=0.29$, $m=1.54$, and
$\alpha=1.08$, while for $m_{r}<25$, $A=46.12\;\textnormal{arcmin}^{-2}$,
$z_{0}=1.09$, $m=0.75$, and $\alpha=2.04$.

The shot noise $N_{s}=\frac{1}{n}$ within each $z$-bin was calculated by
integrating $\frac{d{n}}{d{z}}$ (converting $\textnormal{arcmin}^{-2}$ to sr)
over that $z$-bin; that is, for a $z$-bin with $z_{0}<z<z_{1}$,

$\displaystyle
N_{s}(z)=\left(\int_{z_{0}}^{z_{1}}\frac{d{n}}{d{z}}dz\right)^{-1}.$

The matter power spectrum $P_{\delta}\left(k,z\right)$ is modeled using
halofit (Smith et al., 2003). Using $\frac{d{n}}{d{z}}$ and
$P_{\delta}\left(k,z\right)$, we use Eq. (12) setting $b(z)=1$ to calculate
the angular matter power spectrum.

To calculate $b(z)$, we first estimate the sample galaxy angular power
spectrum $\hat{C}_{\ell}$ (which includes shot noise) directly from the
CosmoDC2 cutout using the 12 $z$-bins we used to calculate
$\frac{d{n}}{d{z}}$. We then estimate the bias within each redshift bin by
minimizing the residual of a weighted least-squares measure between the sample
galaxy power spectrum $\hat{C}_{\ell}$ and the halofit power spectrum (plus
noise) $C_{\ell}+N_{s}$:

$\displaystyle
R(b)=\sum_{\ell}\frac{\left(\hat{C}_{\ell}(z)-\left[b^{2}C_{\ell}(z)+N_{s}(z)\right]\right)^{2}}{\Delta
C_{\ell}^{2}(z)},$ (14) $\displaystyle\Delta
C_{\ell}^{2}(z)=\frac{2}{(2\ell+1)}\left(b^{2}C_{\ell}(z)+N_{s}(z)\right)^{2}.$

Any required $b(z)$ for redshifts not exactly coinciding with the $z$-bins
used for this calibration is obtained by linear interpolation.

Plots of $\frac{d{n}}{d{z}}$, $b(z)$, and $C_{\ell}$ can be seen in Figure 1.
For many CMB lensing reconstructions, a typical $\ell_{\textnormal{max}}$ in
the CMB lensing field is around $3,\\!000$. We will here use the galaxy number
density field at full resolution, treating the shot noise as the limiting
source of noise.

We will assume that the angular power spectrum for the LSST-like survey is
measured over $20,\\!000$ $\textnormal{deg}^{2}$ field
($f_{\textnormal{sky}}=0.48$), which matches the projected sky coverage of
LSST (LSST Science Collaborations and LSST Project, 2009). For the DES-like
survey, we will use a sky coverage of $5,\\!000$ $\textnormal{deg}^{2}$
(Abbott et al., 2022)

###  2.4 The theoretical convergence power spectrum

The angular power spectrum for the convergence field at some comoving distance
$\chi_{s}$ is given by (Schaan et al., 2018)

$\displaystyle
C_{L}(z)=\frac{9}{4}\frac{H_{0}^{4}}{c^{4}}\Omega_{m,0}^{2}\int_{0}^{\chi_{s}}\frac{\chi^{2}\left(\chi_{s}-\chi\right)^{2}}{a^{2}\chi_{s}^{2}}P_{\delta}\left(\frac{L}{\chi},z\right)\frac{d\chi}{\chi^{2}}.$
(15)

Lensing reconstruction using galaxy density fields as a source has great
potential due to the massive amounts of galaxy position data available from
existing surveys such as DES. This data has not been used for lensing
reconstruction, so an efficient reconstruction method could be used to extract
information from this otherwise untapped data source. This algorithm could
also be used with new surveys such as LSST which provide unprecedented depth
and precision in galaxy position measurements (Ivezić et al., 2019; Connolly
et al., 2014; Korytov et al., 2019). This information could be used to
complement CMB lensing measurements. As an independent measurement, this would
help improve the statistical power of observations and identify systematic
errors.

When calculating $C_{L}$ and $N_{L}$ for $m_{r}<25$, we used all bins with
$0<z<3$. However, for the DES-like model ($m_{r}<23$), we performed these
calculations only for $0<z<1.6$, as it was difficult to get reliable
clustering measurements at such low number density. Shear measurements from
DES Y3 extended only to redshift $z<1$ (Abbott et al., 2022).

##  3 Results

In this section, we aim to determine the viability of using galaxy number
density fields as a source for reconstructing lensing fields using either the
standard quadratic estimator or a shear-only estimator. We will primarily
focus on the LSST-like survey, then briefly discuss the equivalent
calculations for a DES-like survey.

###  3.1 Convergence power spectra

The convergence spectrum $C_{L}^{\kappa\kappa}$ (solid lines) and
reconstruction noise $N_{L}$ (quadratic estimator: dashed lines, shear
estimator: dotted lines) are shown in Figure 2 for two choices of magnitude
cuts and $z$-bins with bin width $\Delta z=0.1$. This plot shows these
quantities when integration was performed to
$\ell_{\textnormal{max}}=70,\\!000$. The noise is substantially higher than
the signal at all $L$ for all choices of $z$, both for the traditional
quadratic estimator and the shear estimator.

The reconstruction noise is roughly constant for small $L$; however, the
integration kernels (Eq. (7) and (8)) used to reconstruct the power spectra
become less responsive to the signal at high $L$. This lack of response to
signal manifests as a large spike in the noise as described in (Schaan &
Ferraro, 2019). When integrating to much larger $\ell$ than is typically used
for CMB lensing reconstructions (Schaan & Ferraro, 2019), the divergence of
$N_{L}$ occurs at values of $L$ well beyond where the shape noise already
dominates the cross-spectrum between the lensing reconstruction and cosmic
shear measurements.

For verification, we calculate the total SNR of the autocorrelation function,
defined as

$\displaystyle\left(\frac{S}{N}\right)^{2}_{a}=\sum_{L,z}\frac{f_{\textnormal{sky}}(2L+1)}{2}\left(\frac{C_{L}^{\kappa\kappa}(z)}{C_{L}^{\kappa\kappa}(z)+N_{L}^{\kappa\kappa}(z)}\right)^{2}$
(16)

This is equivalent to Eq. (4.1) in (Nistane et al., 2022). For the survey
parameters we use, the signal is much smaller than the noise for all $L$ as
seen in Figure 2. Because of this, the SNR of the autocorrelation function
will be small. We find that, for the LSST-like survey we use, the total SNR of
the autocorrelation function for the standard quadratic estimator is
$\left(\frac{S}{N}\right)_{a}=5.25$ when performing integration to
$\ell=70,\\!000$. Including small scales in the source field has a large
effect. If we cut off our integration at $\ell=3000$, the total SNR of the
autocorrelation function for the standard quadratic estimator is only 0.52.
This is also substantially (two orders of magnitude) smaller than the total
signal to noise found in Nistane et al. (2022), as expected, where most of the
signal to noise was coming from the inclusion of magnification bias.

Figure 2: The convergence power spectra used for our reconstructions for
$m_{r}<23$ (orange line) and $m_{r}<25$ (purple line). These consist of the
convergence spectrum $C_{L}^{\kappa\kappa}$ (solid) and reconstruction noise
$N_{L}$ (quadratic estimator: dashed, shear estimator: dotted). The $z$-bins
used represent the bin with the maximum $\textnormal{SNR}_{\Delta}(L,z)$. The
noise for shear estimator diverges at very large $L$, leading to the spike
near $L=20,\\!000$.

###  3.2 Signal-to-noise ratio for Cross-Spectra

The lensing potential can be estimated in a wide variety of other ways. The
most direct comparison would be with a measurement of the cosmic shear of
galaxies at the same redshift as the source plane used for the reconstruction,
but one could also estimate the expected lensing from adding up the
contributions from the observed LSS tracers along the line of sight. We define
the differential SNR for the cross spectrum between a different measurement of
the lensing field, with total power spectrum $C_{L}+N_{\gamma}$ and the
reconstructed signal, with power spectrum $C_{L}+N_{L}$, within an $L-z$-bin
with widths $\Delta z$ and $\Delta L$ as

$\displaystyle\left(\frac{S}{N}\right)_{\Delta}(L,z)=\frac{C_{L}(z)}{\Delta
C_{L}^{ij}(z)},$ (17)

where for our purposes (Knox, 1995)

$\displaystyle\left(\Delta
C_{L}^{ij}\right)^{2}=\frac{1}{(2L+1)}\frac{\left(C_{L}^{2}+C_{L}^{i}C_{L}^{j}\right)}{f_{\textnormal{sky}}\,\Delta
L},$ $\displaystyle\hat{C}_{L}^{i}=C_{L}+N_{i},$ (18)
$\displaystyle\hat{C}_{L}^{j}=C_{L}+N_{j},$

where $C_{L}$ is the convergence power spectrum.

The cross spectra $\textnormal{SNR}_{\Delta}(L,z)$ (Eq. (17)) for $m_{r}<25$
are shown in Figure 3. The colored background shows the quadratic estimator
SNR, while the black contours show the shear estimator SNR.

Using Eqs. (17) and (18), we can obtain cumulative signal to noise for any
collection of $z_{j}$ and $M_{i}$ bins by adding them in quadrature:
$\textnormal{SNR}_{\textnormal{cumulative}}=\sum_{ij}\left(\frac{\textnormal{S}}{\textnormal{N}}\right)_{\Delta}(L_{i},z_{j})$.

Including all $L$ bins in the sum for a single $z$ bin leads to the total
signal to noise for that $z$ bin. This is shown as the curve on the right side
of Figure 3. Similarly, instead of summing over all $L$ bins for a given $z$
bin we can sum over all $z$ bins for a fixed $L$ bin to get the total signal
to noise in a given $L$ bin, shown as the curve at the top of Figure 3.

Both without (top panel) and with (bottom panel) shape noise,
$\textnormal{SNR}_{\Delta}(L,z)$ peaked near $z=1.15$, reflecting the fact
that the lensing signal is boosted by increased source distance while the
noise amplitude increases at higher redshift as $\frac{d{n}}{d{z}}$ decreases,
the latter of which falls off rapidly past $z\approx 1.3$. Without shape
noise, the shear estimator cross-spectrum $\textnormal{SNR}_{\Delta}(L,z)$ is
peaked at lower $L$ than the quadratic estimator, indicating that the shear
estimator is less sensitive to information at smaller scales. The increase in
the shear estimator SNR near $L=40,\\!000$ represents the information that can
be obtained beyond the noise spike near $L=10,\\!000$. With shape noise, the
$z$ and $L$ dependence of the quadratic estimator cross-spectrum is nearly
identical to the shear estimator for $\textnormal{SNR}_{\Delta}(L,z)$, just
with a larger magnitude.

(a) $m_{r}<25$, $\gamma_{\textnormal{rms}}=0$

(b) $m_{r}<25$, $\gamma_{\textnormal{rms}}=0.28$

Figure 3: The differential SNR per $L-z$-bin (Eq. (17)) for $\Delta\log
L=0.08$, $\Delta z=0.1$ for the cross spectrum. The colored image shows the
SNR for the standard quadratic estimator, while the black contours show the
SNR for the shear estimator. The inset plots show the differential SNR per
$L-z$-bin summed in quadrature over $z$ (top) and $L$ (right). We performed
these calculations using 20,000 $\textnormal{deg}^{2}$ of the sky. In the
upper plot, the sudden increase in the shear estimator SNR near $L=30,\\!000$
corresponds with the region in $L$ beyond which the response of the shear
estimator to the signal vanishes.

Figure 4: Cumulative SNR in a cross-spectrum between a reconstruction and
measurements of cosmic shear as a function of $L$, integrated over all
redshifts, for (left) the standard quadratic estimator and (right) the shear
estimator as a function of $L$ for various $\Delta z$ values. Colored lines
with circle markers represent the SNR for $m_{r}<25$, while the black line
with white diamond markers shows the SNR for $m_{r}<23$. The top row shows the
SNR without shape noise, while the bottom row shows it with LSST-like shape
noise; the order of the $\Delta z$ values from top to bottom remains the same.
The total SNR increases as $\Delta z$ decreases.

###  3.3 $\Delta z$ dependence

When performing tomographic measurements, it is important to ensure that each
$z$-bin contains a sufficient number of galaxies to prevent shot noise from
dominating the signal. However, making $\Delta z$ as small as possible will
generally improve reconstruction accuracy because $C_{\ell}$ will vary less
within the bin as well as having simply more measurements of the line of sight
lensing potential. While a wide $z$-bin will generally contain more galaxies
(hence lowering the shot noise and shape noise), $C_{\ell}$ will also be
smaller due to projection effects from the larger integration range in $z$.
Unfortunately, depending on the magnitude cut, survey depth, photo-$z$
uncertainties, and other factors, the minimum width of $z$-bins will be
limited. Thus, we will assess how our results change when varying $\Delta z$.

The cumulative SNR for all $z$-bins combined for an LSST-like survey for
various $\Delta z$ values is shown in Figure 4 (colored lines with circle
markers). This SNR was calculated by subdividing the full $0<z<3$ range into
bins with four different widths. For a given $L-z$-bin, neither $C_{L}$ nor
$N_{L}$ vary with $\Delta z$ because the effects of galaxy number counts on
shot noise are exactly canceled by the projection effects from performing the
integration in Eqs. (12) and (15). From Eq. (17), and using the fact that the
reconstructions are always strongly noise-dominated, it can be seen that the
signal to noise in a single $L,z$ bin will be unchanged by the size of the $z$
bin in the absence of shape noise. Because smaller bins allow for more bins in
redshift, the integrated SNR will increase as $1/\sqrt{\Delta z}$ in the
absence of shape noise.

However, the shape noise $N_{\gamma}$ does change with $\Delta z$ because it
depends on the number of galaxies within a given $z$-bin. When we include
shape noise, the SNR even for a single $L,z$ bin does vary with bin width. At
small $L$, the signal in the shear-based reconstruction still dominates over
shape noise, and at high $L$, the shape noise overtakes the signal at larger
$L$ and the cumulative SNR as a function of $L$ flattens out for all $\Delta
z$.

###  3.4 DES-like survey

The Dark Energy Survey already provides a massive amount of galaxy data that
could currently be used for lens reconstruction. To determine the
effectiveness of our method using this dataset, we replicated the above
calculations for a DES-like survey. We assumed a $5,\\!000$
$\textnormal{deg}^{2}$ DES-like survey with a maximum $r$-band magnitude of 23
and used a redshift range of $0<z<1.6$ (Abbott et al., 2022). While the SNR
can theoretically be calculated for higher $z$, our method for calculating
$b(z)$ was not valid for $z\gtrapprox 1.6$ because the shot noise in the mock
catalogs completely dominated the signal. Early estimates for the rms shape
noise (not including measurement error) for a DES-like survey were
$\gamma_{\textnormal{rms}}=0.25$ (Albrecht et al., 2006), while the standard
shear noise estimate (which includes both shape noise and measurement error
$\sigma_{m}$) is $\gamma_{\textnormal{rms}}^{2}+\sigma_{m}^{2}=0.255^{2}$
(Chang et al., 2013); as a conservative estimate and to match the parameters
of the LSST-like survey, we will estimate the rms shear noise as
$\gamma_{\textnormal{rms}}=0.28$ when calculating the SNR.

Figure 2 shows $C_{L}$ and $N_{L}$ for the full $5,\\!000$
$\textnormal{deg}^{2}$ DES-like survey for the $z$-bin that maximized
$\textnormal{SNR}_{\Delta}(L,z)$. The increased shot noise increased the
reconstruction noise $N_{L}$ by as much as an order of magnitude for some
$z$-bins relative to the LSST-like survey. There, we again see that the signal
is much smaller than the noise for all $L$.

Figure 4 shows the full-survey SNR for $m_{r}<23$, $\Delta z=0.1$ (black line
with white diamond markers). While the total SNR with shape noise is
substantially lower than for the LSST-like survey (about 5.6 for the quadratic
estimator, 4.6 for the shear estimator for $\Delta z=0.1$) due to the
increased $N_{L}$, it remains above unity. This indicates that it may be
possible to use DES galaxy number density data as a source field for lensing
reconstruction; however, a more thorough analysis of additional error
contributions is required. In addition, we found that there was substantial
signal for $z>1.6$, indicating that more information could be obtained from a
DES-like survey if the bias calculations could be performed at higher $z$.

##  4 Conclusion

We demonstrated that a galaxy number density field can be used as a source for
convergence field reconstruction. We performed our reconstructions with a
modification of the standard quadratic estimator and a shear-only estimator
using a background galaxy density field as a source distribution. We
demonstrated that this method can be used to reconstruct the convergence power
spectrum $C_{L}$ and we estimated the SNR in cross-correlation that can be
obtained from an LSST-like survey within various tomographic bins at two
magnitude cuts, ranging from barely detectable with $m_{r}<23$ to over 50 for
$m_{r}<25$.

The total SNR using both our standard quadratic estimator and shear estimator
increases with galaxy number density, indicating that increasing the depth and
sensitivity of surveys will improve the SNR within each tomographic bin. We
expect that using the full LSST 5$\sigma$ 10 year coadd. dataset would
increase the SNR. However, we still expect to get a high SNR when using only
single snapshots from an LSST-like survey. The SNR in cross-correlation
continues to increase as the redshift resolution increases, indicating that
accurate photometric redshifts will be of great value.

Whether this reconstruction method will be useful for DES data is unclear. The
total SNR is greater than unity for both estimators, but the inclusion of
additional error sources will decrease the SNR. The information that we
obtained in our forecast from the DES-like survey is also somewhat limited by
the $z=1.6$ cutoff imposed in our analysis. If the galaxy clustering can be
measured at higher redshift, we expect that more information might be
obtained. Future work should recalculate the SNR for a DES-like survey beyond
this $z=1.6$ cutoff.

The analysis we performed used only reconstruction noise and shape noise as
sources of error for $C_{L}$; it did not include some error sources associated
with survey measurements. These include foreground contamination and
photometric redshift and other measurement uncertainties. Future work should
recalculate the SNR with these error sources included. This would provide
information about the effectiveness of our method for performing convergence
field reconstruction from observed galaxy survey data.

A substantial concern that we have neglected for this analysis is the non-
Gaussianity of the galaxy number density counts on small scales, exactly the
scales where we are obtaining most of our lensing information. Correcting for
non-Gaussianities would help further improve our estimator. The non-
Gaussianities and strong nonlinearities at small scales in galaxy density
fields may make higher-order moments of the galaxy density field (the
bispectrum, trispectrum, etc.) non-negligible and introduce additional noise.
Line intensity mapping (LIM) source fields are highly nonlinear because they,
like the galaxy number counts, trace the underlying matter distribution.
Reconstruction methods designed for 21 cm LIM (Pen, 2004; Zahn & Zaldarriaga,
2006; Foreman et al., 2018; Lu et al., 2010) may be applicable for
reconstruction using a galaxy density field as a source. Future work should
evaluate the impact of non-Gaussianities and foreground contamination on the
SNR of reconstructions using our estimator. We would evaluate if, like the
estimators described by Lu & Pen (2008) and Schaan & Ferraro (2019), the shear
estimator is less susceptible to foreground contamination and non-
Gaussianities than the standard quadratic estimator when using a galaxy
density source field. The impact of non-Gaussianities should be compared with
reconstructions using the methods of Foreman et al. (2018) and Lu et al.
(2010) for a galaxy density field background because these estimators are
designed to reduce non-Gaussian noise (or, for the case of Lu et al. (2010),
are optimal for non-Gaussian fields).

##  5 Software Used

Our calculations of the angular power spectra, convergence spectra, and
reconstruction noise were performed using code based off of LensQuEst (Schaan
& Ferraro, 2019). CosmoDC2 data (Korytov et al., 2019) was used for various
calculations and was obtained using GCRCatalogs (Mao et al., 2018) and CatSim
(Connolly et al., 2014).

##  6 Acknowledgements

This work was supported by the Brand & Monica Fortner Chair in Theoretical
Astrophysics and by the National Science Foundation Graduate Research
Fellowship Program Grant No. DGE — 1746047.

This work utilizes resources supported by the National Science Foundation’s
Major Research Instrumentation program, grant #1725729, as well as the
University of Illinois at Urbana-Champaign (Kindratenko et al., 2020).

This research was supported in part by Perimeter Institute for Theoretical
Physics. Research at Perimeter Institute is supported by the Government of
Canada through the Department of Innovation, Science and Economic Development
Canada and by the Province of Ontario through the Ministry of Research,
Innovation and Science.

## References

  * Abbott et al. (2022) Abbott, T. M. C., Aguena, M., Alarcon, A., et al. 2022, Phys. Rev. D, 105, 023520, doi: 10.1103/PhysRevD.105.023520
  * Ade et al. (2015) Ade, P. A. R., Arnold, K., Atlas, M., et al. 2015, Phys. Rev. D, 92, 123509, doi: 10.1103/PhysRevD.92.123509
  * Aghanim et al. (2020) Aghanim, N., et al. 2020, Astron. Astrophys., 641, A8, doi: 10.1051/0004-6361/201833886
  * Agustin Lozano Torres & Schaefer (2022) Agustin Lozano Torres, J., & Schaefer, B. M. 2022, arXiv e-prints, arXiv:2205.03713.  https://arxiv.org/abs/2205.03713
  * Albrecht et al. (2006) Albrecht, A., et al. 2006.  https://arxiv.org/abs/astro-ph/0609591
  * Bartelmann & Schneider (2001) Bartelmann, M., & Schneider, P. 2001, Phys. Rep., 340, 291, doi: 10.1016/S0370-1573(00)00082-X
  * BICEP/Keck Collaboration et al. (2022) BICEP/Keck Collaboration, :, Ade, P. A. R., et al. 2022, arXiv e-prints, arXiv:2203.16556, doi: 10.48550/arXiv.2203.16556
  * Bucher et al. (2012) Bucher, M., Carvalho, C. S., Moodley, K., & Remazeilles, M. 2012, Phys. Rev. D, 85, 043016, doi: 10.1103/PhysRevD.85.043016
  * Camera et al. (2014) Camera, S., Fedeli, C., & Moscardini, L. 2014, JCAP, 2014, 027, doi: 10.1088/1475-7516/2014/03/027
  * Chakraborty & Pullen (2019) Chakraborty, P., & Pullen, A. R. 2019, MNRAS, 488, 1828, doi: 10.1093/mnras/stz1781
  * Chang et al. (2013) Chang, C., Jarvis, M., Jain, B., et al. 2013, Monthly Notices of the Royal Astronomical Society, 434, 2121, doi: 10.1093/mnras/stt1156
  * Connolly et al. (2014) Connolly, A. J., Angeli, G. Z., Chandrasekharan, S., et al. 2014, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9150, Modeling, Systems Engineering, and Project Management for Astronomy VI, ed. G. Z. Angeli & P. Dierickx, 915014
  * Feng & Holder (2019) Feng, C., & Holder, G. 2019, arXiv e-prints, arXiv:1905.02084.  https://arxiv.org/abs/1905.02084
  * Foreman et al. (2018) Foreman, S., Meerburg, P. D., van Engelen, A., & Meyers, J. 2018, JCAP, 2018, 046, doi: 10.1088/1475-7516/2018/07/046
  * Gatti et al. (2021) Gatti, M., et al. 2021, Mon. Not. Roy. Astron. Soc., 504, 4312, doi: 10.1093/mnras/stab918
  * Hilbert et al. (2007) Hilbert, S., Metcalf, R. B., & White, S. D. M. 2007, MNRAS, 382, 1494, doi: 10.1111/j.1365-2966.2007.12615.x
  * Ivezić et al. (2019) Ivezić, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, ApJ, 873, 111, doi: 10.3847/1538-4357/ab042c
  * Jalilvand et al. (2019) Jalilvand, M., Majerotto, E., Durrer, R., & Kunz, M. 2019, JCAP, 2019, 020, doi: 10.1088/1475-7516/2019/01/020
  * Kindratenko et al. (2020) Kindratenko, V., Mu, D., Zhan, Y., et al. 2020, in Practice and Experience in Advanced Research Computing, PEARC ’20 (New York, NY, USA: Association for Computing Machinery), 41–48.  https://doi.org/10.1145/3311790.3396649
  * Knox (1995) Knox, L. 1995, Phys. Rev. D, 52, 4307, doi: 10.1103/PhysRevD.52.4307
  * Komatsu et al. (2011) Komatsu, E., Smith, K. M., Dunkley, J., et al. 2011, The Astrophysical Journal Supplement Series, 192, 18, doi: 10.1088/0067-0049/192/2/18
  * Korytov et al. (2019) Korytov, D., Hearin, A., Kovacs, E., et al. 2019, The Astrophysical Journal Supplement Series, 245, 26, doi: 10.3847/1538-4365/ab510c
  * Lewis & Challinor (2006) Lewis, A., & Challinor, A. 2006, Phys. Rep., 429, 1, doi: 10.1016/j.physrep.2006.03.002
  * Limber (1953) Limber, D. N. 1953, ApJ, 117, 134, doi: 10.1086/145672
  * LSST Science Collaborations and LSST Project (2009) LSST Science Collaborations and LSST Project. 2009, arXiv e-prints, arXiv:0912.0201.  https://arxiv.org/abs/0912.0201
  * Lu & Pen (2008) Lu, T., & Pen, U.-L. 2008, MNRAS, 388, 1819, doi: 10.1111/j.1365-2966.2008.13524.x
  * Lu et al. (2010) Lu, T., Pen, U.-L., & Doré, O. 2010, Phys. Rev. D, 81, 123015, doi: 10.1103/PhysRevD.81.123015
  * Luo et al. (2018) Luo, W., Yang, X., Lu, T., et al. 2018, The Astrophysical Journal, 862, 4, doi: 10.3847/1538-4357/aacaf1
  * Ma et al. (2023) Ma, R., Zhang, P., Yu, Y., & Qin, J. 2023, doi: 10.1093/mnras/stad3664
  * Maartens et al. (2021) Maartens, R., Fonseca, J., Camera, S., et al. 2021, JCAP, 2021, 009, doi: 10.1088/1475-7516/2021/12/009
  * Madhavacheril et al. (2023) Madhavacheril, M. S., Qu, F. J., Sherwin, B. D., et al. 2023, arXiv e-prints, arXiv:2304.05203, doi: 10.48550/arXiv.2304.05203
  * Mandelbaum et al. (2018) Mandelbaum, R., et al. 2018.  https://arxiv.org/abs/1809.01669
  * Mao et al. (2018) Mao, Y.-Y., Kovacs, E., Heitmann, K., et al. 2018, ApJS, 234, 36, doi: 10.3847/1538-4365/aaa6c3
  * Metcalf et al. (2020) Metcalf, R. B., Tessore, N., & Croft, R. A. C. 2020, A&A, 642, A122, doi: 10.1051/0004-6361/202038056
  * Metcalf & White (2007) Metcalf, R. B., & White, S. D. M. 2007, MNRAS, 381, 447, doi: 10.1111/j.1365-2966.2007.12212.x
  * Metcalf & White (2009) —. 2009, MNRAS, 394, 704, doi: 10.1111/j.1365-2966.2009.14401.x
  * Nistane et al. (2022) Nistane, V., Jalilvand, M., Carron, J., Durrer, R., & Kunz, M. 2022, Journal of Cosmology and Astroparticle Physics, 2022, 024, doi: 10.1088/1475-7516/2022/06/024
  * Okamoto & Hu (2003) Okamoto, T., & Hu, W. 2003, Phys. Rev. D, 67, 083002, doi: 10.1103/PhysRevD.67.083002
  * Omori et al. (2023) Omori, Y., Baxter, E. J., Chang, C., et al. 2023, Phys. Rev. D, 107, 023529, doi: 10.1103/PhysRevD.107.023529
  * Peebles (1980) Peebles, P. J. E. 1980, The large-scale structure of the universe
  * Pen (2004) Pen, U.-L. 2004, New Astron., 9, 417, doi: 10.1016/j.newast.2004.01.006
  * Ponce et al. (2013) Ponce, R., Sánchez, E., & Sevilla, I. 2013, in Highlights of Spanish Astrophysics VII, ed. J. C. Guirado, L. M. Lara, V. Quilis, & J. Gorgas, 459–459
  * Pourtsidou & Metcalf (2014) Pourtsidou, A., & Metcalf, R. B. 2014, MNRAS, 439, L36, doi: 10.1093/mnrasl/slt175
  * Pourtsidou & Metcalf (2015) —. 2015, MNRAS, 448, 2368, doi: 10.1093/mnras/stv102
  * Schaan & Ferraro (2019) Schaan, E., & Ferraro, S. 2019, Phys. Rev. Lett., 122, 181301, doi: 10.1103/PhysRevLett.122.181301
  * Schaan et al. (2018) Schaan, E., Ferraro, S., & Spergel, D. N. 2018, Phys. Rev. D, 97, 123539, doi: 10.1103/PhysRevD.97.123539
  * Scranton et al. (2005) Scranton, R., Ménard, B., Richards, G. T., et al. 2005, ApJ, 633, 589, doi: 10.1086/431358
  * Smith et al. (2003) Smith, R. E., Peacock, J. A., Jenkins, A., et al. 2003, MNRAS, 341, 1311, doi: 10.1046/j.1365-8711.2003.06503.x
  * von Wietersheim Kramsta et al. (2021) von Wietersheim Kramsta, M., Joachimi, B., van den Busch, J. L., et al. 2021, Monthly Notices of the Royal Astronomical Society, 504, 1452, doi: 10.1093/mnras/stab1000
  * Zahn & Zaldarriaga (2006) Zahn, O., & Zaldarriaga, M. 2006, ApJ, 653, 922, doi: 10.1086/508916
  * Zhang & Pen (2006) Zhang, P., & Pen, U.-L. 2006, Mon. Not. Roy. Astron. Soc., 367, 169, doi: 10.1111/j.1365-2966.2006.09971.x