# Improved Modeling of the Kinematic Sunyaev-Zel’dovich Projected-Fields
signal and its Cosmological Dependence

Raagini Patki  Department of Astronomy, Cornell University, Ithaca, NY 14853,
USA.  Nicholas Battaglia  Department of Astronomy, Cornell University, Ithaca,
NY 14853, USA.  Simone Ferraro  Lawrence Berkeley National Laboratory, One
Cyclotron Road, Berkeley, CA 94720, USA.  Berkeley Center for Cosmological
Physics, Department of Physics, University of California, Berkeley, CA 94720,
USA.

###### Abstract

Over the past decade, the kinematic Sunyaev-Zel’dovich (kSZ) effect has
emerged as an observational probe of the distribution of baryons and velocity
fields in the late Universe. Of the many ways to detect the kSZ, the
‘projected-fields kSZ estimator’ has the promising feature of not being
limited to galaxy samples with accurate redshifts. The current theoretical
modeling of this estimator involves an approximate treatment only applicable
at small scales. As the measurement fidelity rapidly improves, we find it
necessary to move beyond the original treatment and hence derive an _improved_
theoretical model for this estimator without these previous approximations. We
show that the differences between the predicted signal from the two models are
scale-dependent and will be significant for future measurements from the
Simons Observatory and CMB-S4 in combination with galaxy data from WISE or the
Rubin Observatory, which have high forecasted signal-to-noise ratios ($>100$).
Thus, adopting our improved model in future analyses will be important to
avoid biases. Equipped with our model, we explore the cosmological dependence
of this kSZ signal for future measurements. With a Planck prior, residual
uncertainty on $\Lambda$CDM parameters leads to $\sim 7\%$ marginalized
uncertainties on the signal amplitude, compared to a sub-percent level
forecasted with a fixed cosmology. To illustrate the potential of this kSZ
estimator as a cosmological probe, we forecast initial constraints on
$\Lambda$CDM parameters and the sum of neutrino masses, paving the way for
jointly fitting both baryonic astrophysics and cosmology in future analyses.

##  I Introduction

Upcoming surveys such as the Simons Observatory (SO) [1] and CMB-S4 [2] will
map the cosmic microwave background (CMB) at arcminute-resolution, allowing
the precise measurement of secondary CMB anisotropies which arise due to the
large-scale structure (LSS) at low redshifts. The kinematic Sunyaev-Zel’dovich
(kSZ) effect [3, 4, 5] is one of the dominant sources of secondaries at
arcminute scales. It is the Doppler boosting of CMB photons as they Compton-
scatter off of free electrons moving with a non-zero line-of-sight (LOS)
velocity. The resulting shift in the observed CMB temperature is proportional
to the free electron momentum along the LOS, thus making the kSZ a key probe
of both astrophysics (e.g. [6, 7, 8, 9, 10]) and cosmology (e.g. [11, 12, 13,
14, 15]).

The kSZ effect was first detected in 2012 in Atacama Cosmology Telescope (ACT)
data [16] using the pairwise momenta of galaxies (e.g. [17]). This method was
subsequently used for detections in data from Planck [18], the South Pole
Telescope [19], and with updated ACT DR5 maps [20]. A velocity-weighted
stacking approach [21, 22] has also been used for measurements in Planck [23]
and ACT maps [24, 25]. While there are no measurements with the recently
proposed velocity reconstruction method (e.g. [26]) presently, it is
forecasted to become a promising cosmological probe (e.g. [27, 28]). A
limitation of these kSZ estimators is that they need accurate spectroscopic
redshifts, typically available only for smaller galaxy samples. If photometric
data is used instead, their signal-to-noise ratios (SNR) decrease sharply [29,
30], despite the increase in the number of galaxies.

In this paper, we focus on the ‘projected-fields’ estimator for the kSZ [31,
32], which has the advantage that it only needs a statistical redshift
distribution of the LSS tracers, rather than individual redshifts. This
distribution can be obtained from photometric galaxy surveys such as Wide-
field Infrared Survey Explorer (WISE) [33] and the future Vera C. Rubin
Observatory (VRO) [34], which are less expensive and map larger volumes. First
proposed in [31], this estimator extracts the kSZ signal that is associated
with a given tracer population of the LSS, from the CMB temperature map. Since
the LOS velocity of ionized gas (and the corresponding kSZ signal) is equally
likely to be positive or negative, a foreground-cleaned and filtered CMB map
is first squared in real space, and then cross-correlated with a projected
tracer field. Thus, this kSZ2-LSS cross-correlation is the simplest such non-
trivial estimator that does not require 3D information.

Following the initial formulation by [31], [32] further studied this estimator
and suggested an approximate theoretical model for it at small scales
(hereafter ‘D05’). The kSZ2-LSS cross-correlation was first detected in data
from Planck and WMAP with WISE galaxies [35, 6], and then with unWISE galaxies
[36]. The measured signal amplitude ($A_{\mathrm{kSZ}^{2}}$) was used to
address the ‘missing baryon’ problem by inferring the baryon fraction $f_{b}$,
which is proportional to the free electron density at the corresponding galaxy
redshifts. Recently, [37] forecasted future limits on astrophysical parameters
governing the baryon density profile with this signal.

While baryonic effects are largely dominant at small scales, a wider range of
scales is affected by the cosmological dependence of the signal. This
motivates us to do the first derivation of an ‘improved’ theoretical model for
the projected-fields estimator that is accurate across these larger scales,
while the D05 model used so far is only applicable at small scales. Moreover,
we show that there are $\sim$15$\%$ level differences between the signal
predicted by our model and the D05 model, which are significant for future CMB
surveys that have been forecasted to yield large SNR for the projected-fields
estimator [6, 37]. Since these differences are scale-dependent, as well as
large even at small scales, the resulting change in the shape of the signal
would likely impact _any_ parameter inference with this estimator.

This improved theoretical model enables us to explore the cosmological
dependence of the signal, paving the way for future analyses to jointly fit
both astrophysical and cosmological parameters. We study the cosmological
sensitivity of the improved signal for the first time, assuming a $\Lambda$CDM
cosmology with best-fit Planck-2018 parameters [38] as our baseline fiducial
model. Since the kSZ-induced anisotropy is proportional to the LOS velocity of
free electrons, it is an unbiased probe of the total underlying matter density
and the growth rate of LSS. Massive neutrinos directly impact the clustering
and growth rate of LSS (e.g. [39, 40]). Therefore, we additionally vary the
sum of neutrino masses ($\Sigma m_{\nu}$) about a fiducial value of 60 meV, as
an illustration of the potential sensitivity of kSZ measurements using this
estimator to physics beyond the Standard Model.

The rest of this paper is organized as follows: In Section II, we review the
theoretical formalism of the estimator, its existing approximate model, and
our ‘improved’ model, providing its rigorous derivation in the Appendix. In
Section III, we describe the survey specifications and our numerical
implementation. We compare the predicted signal from our model and the D05
model for future measurements in Section IV. In Section V, we present initial
Fisher forecasts jointly fitting for $\Lambda$CDM parameters, $\Sigma
m_{\nu}$, $A_{\mathrm{kSZ}^{2}}$, and the galaxy bias, for future SO and
CMB-S4 data combined with WISE and a VRO-like experiment. We show the effect
of marginalizing over cosmological parameters and the impact of residual
foregrounds present in cleaned CMB maps on these forecasts, and discuss future
directions. We conclude the paper in Section VI.

##  II Theoretical Formalism

###  II.1 The Projected-Fields kSZ estimator

The kSZ effect leads to a shift in the CMB temperature anisotropies in a
direction $\mathbf{\hat{n}}$ on the sky,
$\Theta^{\mathrm{kSZ}}\left(\mathbf{\hat{n}}\right)\equiv\Delta
T^{\mathrm{kSZ}}/T_{\mathrm{CMB}}\left(\mathbf{\hat{n}}\right)$, (taking c =
1):

$\Theta^{\mathrm{kSZ}}\left(\mathbf{\hat{n}}\right)=-\int d\eta\hskip
1.42271ptg(\eta)\hskip 2.84544pt\mathbf{p_{e}}\cdot\mathbf{\hat{n}},$ (1)

where $\eta(z)$ is the comoving distance to redshift $z$, and
$\mathbf{p_{e}}=(1+\delta_{e})\mathbf{v}_{e}$ is the electron momentum field,
with $\delta_{e}$ being the electron overdensity, and $\mathbf{v}_{e}$ the
peculiar velocity of electrons. $g(\eta)=e^{-\tau}d\tau/d\eta$ is the
visibility function, and $\tau$ is the optical depth to Thomson scattering.
Since $d\tau/d\eta=\sigma_{\mathrm{T}}\,n_{e}/(1+z)$, where
$\sigma_{\mathrm{T}}$ is the Thomson scattering cross-section and $n_{e}$ is
the free electron number density,

$\Theta^{\mathrm{kSZ}}\left(\mathbf{\hat{n}}\right)=-\sigma_{\mathrm{T}}\int\frac{d\eta}{1+z}e^{-\tau}n_{e}\left(\mathbf{\hat{n}},\eta\right)\mathbf{v_{e}}\left(\mathbf{\hat{n}},\eta\right)\cdot\mathbf{\hat{n}}.$
(2)

In this work, we consider galaxies as a tracer of the LSS. However, the same
formalism can be applied equivalently to other tracers of the underlying
matter density field, such as quasars, 21-cm fluctuations (e.g. [41, 10]), or
weak lensing convergence (e.g. [37]). The projected galaxy overdensity field
$\delta_{g}(\mathbf{\hat{n}})$ is:

$\delta_{g}(\mathbf{\hat{n}})=\int_{0}^{\eta_{\max}}d\eta\hskip
1.9919ptW^{g}(\eta)\hskip 1.9919pt\delta_{m}(\eta\mathbf{\hat{n}},\eta),$ (3)

where $\delta_{m}(\eta\mathbf{\hat{n}},\eta)$ is the 3D matter overdensity
field. $\eta_{\max}$ is the maximum comoving distance of the galaxy sample,
and its projection kernel $W^{g}(\eta)=b_{g}p_{s}(\eta)$, where
$p_{s}(\eta)\propto dn/d\eta$ is the redshift distribution of the number of
galaxies, normalized to have a unit integral. Here, we assume a linear galaxy
bias $b_{g}$ for simplicity.

We apply a Wiener filter $F(\ell)$ to foreground-cleaned CMB blackbody
temperature maps in harmonic space, in order to select angular scales
dominated by the kSZ contribution [6, 35] :

$F(\ell)=\frac{C_{\ell}^{\mathrm{kSZ}}}{C_{\ell}^{\mathrm{tot}}}.$ (4)

Recently, [37] have proposed a slightly more optimal filter, where
$C_{\ell}^{\mathrm{kSZ}}$ in the numerator of Eq.(4) is replaced by
$\sqrt{C_{\ell}^{\mathrm{kSZ}}}$. While this is slightly more optimal, here we
adopt a classical Wiener filter as in Eq.(4) for comparison with previous
works that adopt this choice.

Here, $C_{\ell}^{\mathrm{kSZ}}$ is the theoretical kSZ power spectrum, while
$C_{\ell}^{\mathrm{tot}}$ is the total power spectrum, which includes the
primary CMB, kSZ, and integrated Sachs-Wolfe (ISW) contributions, detector
noise, and residual foregrounds. In addition, telescopes observe the CMB
through a finite beam $b(\ell)$,

$b(\ell)=\mathrm{exp}\left(-\frac{1}{2}\ell(\ell+1)\frac{\theta_{\mathrm{FWHM}}^{2}}{8\,\mathrm{ln}\,2}\right),$
(5)

where $\theta_{\mathrm{FWHM}}$ is the full width at half maximum in radians,
for a Gaussian beam. Therefore, overall, the filtered CMB anisotropies
$\Theta_{f}(\ell)$ are related to the true CMB anisotropies $\Theta(\ell)$ as:

$\Theta_{f}(\mathbf{\ell})=F(\ell)b(\ell)\Theta(\mathbf{\ell})\equiv
f(\ell)\Theta(\mathbf{\ell}),$ (6)

where $f(\ell)\equiv F(\ell)b(\ell)$. Further details about our numerical
evaluation of the filters are given in Section III.1, and they are shown in
Figure 1.

Free electrons at low-redshift are equally likely to be moving towards or away
from us along the LOS. Due to this $\mathbf{v_{e}}\rightarrow-\mathbf{v_{e}}$
parity symmetry, the naive cross-correlation between $\Theta^{\mathrm{kSZ}}$
and $\delta_{g}$ is expected to vanish at small scales, (while the
contribution from linear ISW is detectable at $\ell<100$ [6]). Hence, we
follow [31, 32] and square the filtered CMB map in real space, before cross-
correlating it with $\delta_{g}$:

$\langle\Theta_{f}^{2}(\mathbf{\ell})\,\delta_{g}(\mathbf{\ell^{^{\prime}}})\rangle=(2\pi)^{2}\delta_{D}(\mathbf{\ell}+\mathbf{\ell^{^{\prime}}})C_{\ell}^{\mathrm{kSZ}^{2}\times\delta_{g}}.$
(7)

Using the Limber (Kaiser) approximation [42, 43], the angular power spectrum
of the $\mathrm{kSZ}^{2}$-galaxy cross-correlation above (i.e. the projected-
fields estimator) can be simplified in harmonic space and written as [31, 32]:

$C_{\ell}^{\mathrm{kSZ}^{2}\times\delta_{g}}=\int_{0}^{\eta_{\max}}\frac{d\eta}{\eta^{2}}W^{g}(\eta)g^{2}(\eta)\mathcal{T}\left(j=\frac{\ell}{\eta},\eta\right),$
(8)

where the ‘triangle power spectrum’ $\mathcal{T}$ is

$\mathcal{T}(j,\eta)=\int\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}f(q\eta)f(|\mathbf{j}+\mathbf{q}|\eta)\,B_{p_{\mathbf{\hat{n}}}p_{\mathbf{\hat{n}}}\delta}(\mathbf{q},-\mathbf{j}-\mathbf{q},\mathbf{j}).$
(9)

Here, the ‘hybrid bispectrum’
$B_{p_{\mathbf{\hat{n}}}p_{\mathbf{\hat{n}}}\delta}$ is the three-point
correlation in Fourier space between two electron momenta along the LOS
($p_{\mathbf{\hat{n}}}=\mathbf{p}\cdot\mathbf{\hat{n}}$), and one fractional
overdensity. $\mathcal{T}(j,\eta)$ is an integral over all (closed) triangles
lying in a constant-redshift plane (corresponding to $\eta(z)$) which have one
side of length $j$, where the wavevector $j\approx\ell/\eta$ follows from the
Limber approximation. $\mathcal{T}$ compresses the rich information contained
in $B_{p_{\mathbf{\hat{n}}}p_{\mathbf{\hat{n}}}\delta}$ into what is needed to
compute $C_{\ell}^{\mathrm{kSZ}^{2}\times\delta_{g}}$.

###  II.2 Existing Approximate Model for the Hybrid Bispectrum

$B_{p_{\mathbf{\hat{n}}}p_{\mathbf{\hat{n}}}\delta}$ is the key quantity that
must be modeled accurately enough, in order to interpret the measured signals
appropriately. The electron momentum field is:
$\mathbf{p}(\mathbf{x})=(1+\delta(\mathbf{x}))\mathbf{v}(\mathbf{x})$ (we drop
the subscript $e$ for clarity), which transforms into Fourier space as:
$\mathbf{p}(\mathbf{k})=\mathbf{v}(\mathbf{k})+\int
d^{3}k^{\prime}\mathbf{v}(\mathbf{k}^{\prime})\delta(\mathbf{k}-\mathbf{k}^{\prime})/(2\pi)^{3}$.
In linear theory, velocities $\mathbf{v}(\mathbf{x})$ are purely gradient, and

$\mathbf{v}(\mathbf{k})=i\frac{faH\delta(\mathbf{k})}{k}\mathbf{\hat{k}},$
(10)

where $a,f$, and $H$ are the scale factor, linear growth rate, and Hubble
parameter at the corresponding redshift, respectively. In the Limber [42] or
‘weak coupling’ [44] approximations, only those $\mathbf{p}(\mathbf{k})$ modes
with $\mathbf{k}$ perpendicular to the LOS contribute to the kSZ effect (Eq.
1). Therefore, the linear effect from the first term
$\mathbf{p}\sim\mathbf{v}$ vanishes [31]. However, the second term
$\mathbf{p}\sim\delta\mathbf{v}$ can have a non-zero “curl” or transverse
component ($\equiv\mathbf{p}_{\perp}$) that has velocity perpendicular to its
wavevector [5]. To capture the full kSZ effect including nonlinear scales, we
follow [32] and only consider the $\delta\mathbf{v}$ term while substituting
the non-linear density field in Eq.(10) for the velocity. [45] first proposed
this phenomenologically, and showed that it agrees excellently with their
simulations.

Under the Limber approximation, the power spectrum of the LOS component
$p_{\mathbf{\hat{n}}}$ is $P_{p_{\mathbf{\hat{n}}}}=(1/2)P_{p_{\perp}}$ [46,
47]. Similarly, our hybrid bispectrum
$B_{p_{\mathbf{\hat{n}}}p_{\mathbf{\hat{n}}}\delta}=(1/2)B_{p_{\perp}p_{\perp}\delta}$.
The fully-nonlinear improved expression for $P_{p_{\mathbf{\hat{n}}}}$ was
derived in [46] (see Eq.(6) therein); they show that it simplifies to
$P_{p_{\mathbf{\hat{n}}}}(k)\approx(1/3)v_{\mathrm{rms}}^{2}P_{\delta\delta}(k)$
in the high-$k$ (i.e. small scales) limit, where $v_{\mathrm{rms}}^{2}$ is the
3D velocity dispersion and $P_{\delta\delta}$ is the nonlinear matter power
spectrum. The original paper [31] suggested an analogous ansatz in the
high-$k_{1},k_{2},k_{3}$ limit:

$B_{p_{\mathbf{\hat{n}}}p_{\mathbf{\hat{n}}}\delta}^{\mathrm{approx}}(\mathbf{k_{1}},\mathbf{k_{2}},\mathbf{k_{3}})\approx\frac{1}{3}v_{\mathrm{rms}}^{2}B_{\delta\delta\delta}(\mathbf{k_{1}},\mathbf{k_{2}},\mathbf{k_{3}}),$
(11)

where $B_{\delta\delta\delta}$ is the nonlinear matter bispectrum. [32] (D05)
provided a rough justification for the application of this ansatz for
estimating $C_{\ell}^{\mathrm{kSZ}^{2}\times\delta_{g}}$ at small scales (i.e.
high-$\ell$); we outline it below. We refer to Eq.(11) as the ‘D05’ or
‘approximate’ model for calculating
$C_{\ell}^{\mathrm{kSZ}^{2}\times\delta_{g}}$. While all previous works using
this estimator assume this model, it makes some approximations that limit its
validity for future surveys; we address them in our _improved_ model described
in the next subsection.

By definition, we have (similar to Eq. 7):

$\displaystyle\langle
p_{\mathbf{\hat{n}}}(\mathbf{k_{1}})p_{\mathbf{\hat{n}}}(\mathbf{k_{2}})\delta(\mathbf{k_{3}})\rangle=(2\pi)^{3}\delta_{D}(\mathbf{k_{1}}+\mathbf{k_{2}}+\mathbf{k_{3}})$
$\displaystyle
B_{p_{\mathbf{\hat{n}}}p_{\mathbf{\hat{n}}}\delta}(\mathbf{k_{1}},\mathbf{k_{2}},\mathbf{k_{3}}),$
(12)

where the Dirac delta function enforces that the three wavevectors must form a
closed triangle, and $\langle
p_{\mathbf{\hat{n}}}p_{\mathbf{\hat{n}}}\delta\rangle=(1/2)\langle
p_{\perp}p_{\perp}\delta\rangle$. Since $\mathbf{p}\sim\delta\mathbf{v}$,
$\langle p_{\perp}p_{\perp}\delta\rangle$ is a five-point function
$\langle\delta\mathbf{v}\delta\mathbf{v}\delta\rangle$, whose Wick contraction
consists of $\binom{5}{2}=10$ different terms that are products of two- and
three- point functionals (Table 1), as first explained by D05. The connected
five-point term is assumed to be negligible due to weak non-Gaussianity (as
demonstrated for $P_{p_{\mathbf{\hat{n}}}}$ in simulations [46]), and the
‘tree-level’ expansion is sufficient. Six out of the ten terms turn out to be
zero due to the parity of the functionals.

Terms | Geometric scaling
---|---
$\langle v^{i}(\mathbf{k})v^{j}(\mathbf{k}^{\prime})\rangle\langle\delta(\mathbf{k}_{1}-\mathbf{k})\delta(\mathbf{k}_{2}-\mathbf{k}^{\prime})\delta(\mathbf{k}_{3})\rangle$ | 1
$\langle v^{i}(\mathbf{k})\delta(\mathbf{k}_{1}-\mathbf{k})\rangle\langle v^{j}(\mathbf{k}^{\prime})\delta(\mathbf{k}_{2}-\mathbf{k}^{\prime})\delta(\mathbf{k}_{3})\rangle$ | 0
$\langle v^{i}(\mathbf{k})\delta(\mathbf{k}_{2}-\mathbf{k}^{\prime})\rangle\langle\delta(\mathbf{k}_{1}-\mathbf{k})v^{j}(\mathbf{k}^{\prime})\delta(\mathbf{k}_{3})\rangle$ | $k/k_{2}$
$\langle\delta(\mathbf{k}_{2}-\mathbf{k}^{\prime})v^{j}(\mathbf{k}^{\prime})\rangle\langle\delta(\mathbf{k}_{1}-\mathbf{k})v^{i}(\mathbf{k})\delta(\mathbf{k}_{3})\rangle$ | 0
$\langle\delta(\mathbf{k}_{1}-\mathbf{k})v^{j}(\mathbf{k}^{\prime})\rangle\langle v^{i}(\mathbf{k})\delta(\mathbf{k}_{2}-\mathbf{k}^{\prime})\delta(\mathbf{k}_{3})\rangle$ | $k/k_{1}$
$\langle v^{i}(\mathbf{k})\delta(\mathbf{k}_{3})\rangle\langle\delta(\mathbf{k}_{1}-\mathbf{k})\delta(\mathbf{k}_{2}-\mathbf{k}^{\prime})v^{j}(\mathbf{k}^{\prime})\rangle$ | 0
$\langle\delta(\mathbf{k}_{3})v^{j}(\mathbf{k}^{\prime})\rangle\langle\delta(\mathbf{k}_{1}-\mathbf{k})\delta(\mathbf{k}_{2}-\mathbf{k}^{\prime})v^{i}(\mathbf{k})\rangle$ | 0
$\langle\delta(\mathbf{k}_{1}-\mathbf{k})\delta(\mathbf{k}_{2}-\mathbf{k}^{\prime})\rangle\langle v^{i}(\mathbf{k})v^{j}(\mathbf{k}^{\prime})\delta(\mathbf{k}_{3})\rangle$ |  $[-k+(k_{1}$ or $k_{2})]/k_{3}$
$\langle\delta(\mathbf{k}_{2}-\mathbf{k}^{\prime})\delta(\mathbf{k}_{3})\rangle\langle v^{i}(\mathbf{k})v^{j}(\mathbf{k}^{\prime})\delta(\mathbf{k}_{1}-\mathbf{k})\rangle$ | 0
$\langle\delta(\mathbf{k}_{1}-\mathbf{k})\delta(\mathbf{k}_{3})\rangle\langle v^{i}(\mathbf{k})v^{j}(\mathbf{k}^{\prime})\delta(\mathbf{k}_{2}-\mathbf{k}^{\prime})\rangle$ | 0
Table 1: The ten terms in the Wick contraction of $\langle
p_{\perp}p_{\perp}\delta\rangle\sim\langle\delta\mathbf{v}\delta\mathbf{v}\delta\rangle$
which contribute to
$B_{p_{\mathbf{\hat{n}}}p_{\mathbf{\hat{n}}}\delta}(\mathbf{k_{1}},\mathbf{k_{2}},\mathbf{k_{3}})$,
as first explained in D05 [32].

D05 stated the following ansatz for the contribution to
$B_{p_{\perp}p_{\perp}\delta}^{\mathrm{approx}}(\mathbf{k_{1}},\mathbf{k_{2}},\mathbf{k_{3}})$
from the leading order (or ‘usual’) term (first row in Table 1) among the four
non-zero terms:

$\int\frac{d^{3}k}{\left(2\pi\right)^{3}}[1-\mu_{1}\mu_{2}]\,P_{vv}(k)B_{\delta\delta\delta}(\mathbf{k_{1}}-\mathbf{k},\mathbf{k_{2}}+\mathbf{k},\mathbf{k_{3}}),$
(13)

where $\mu_{1}\equiv\mathbf{\hat{k_{1}}}\cdot\mathbf{\hat{k}}$ and
$\mu_{2}\equiv\mathbf{\hat{k_{2}}}\cdot\mathbf{\hat{k}}$. The kSZ effect is
dominant at small scales in the CMB, and thus the filters (e.g. Fig. 1) peak
at high-$\ell$. D05 then argued that for this small-scale regime of the
$C_{\ell}^{\mathrm{kSZ}^{2}\times\delta_{g}}$ estimator, only those triangles
with high-$k_{1},k_{2},k_{3}$ need to be considered for
$B_{p_{\mathbf{\hat{n}}}p_{\mathbf{\hat{n}}}\delta}(\mathbf{k_{1}},\mathbf{k_{2}},\mathbf{k_{3}})$
in Eq.(9). In this limit, their ansatz above simplifies to:

$\displaystyle
B_{p_{\mathbf{\hat{n}}}p_{\mathbf{\hat{n}}}\delta}(\mathbf{k_{1}},\mathbf{k_{2}},\mathbf{k_{3}})\approx\frac{1}{2}[1-\mu_{12}]v_{\mathrm{rms}}^{2}B_{\delta\delta\delta}(\mathbf{k_{1}},\mathbf{k_{2}},\mathbf{k_{3}}),$

which is close to the [31] ansatz (Eq. 11 above), where
$\mu_{12}\equiv\mathbf{\hat{k_{1}}}\cdot\mathbf{\hat{k_{2}}}$. The three other
non-zero terms listed in Table 1 turn out to have geometric scalings $\propto
1/k_{i}$ (where $i=1,2,$ or 3). Since D05 chose to confine to the
high-$k_{1},k_{2},k_{3}$ limit, they argued that these additional ‘extra’
terms can be neglected.

###  II.3 Our Improved Model

While the framework and initial basic assumptions for the D05 model are well-
founded, they make a series of further approximations that we find to be
inaccurate. We list these inaccuracies here and describe how our rigorously
derived _improved_ model addresses and improves upon them:

  * •

By definition, the calculation of triangle power spectra $\mathcal{T}(j,\eta)$
given by Eq.(9) is a convolution, in which the shape of the hybrid bispectrum
$B_{p_{\mathbf{\hat{n}}}p_{\mathbf{\hat{n}}}\delta}(\mathbf{q},-\mathbf{j}-\mathbf{q},\mathbf{j})$
gets modified by the filtering product
$f(q\eta)f(|\mathbf{j}+\mathbf{q}|\eta)$. Therefore, even if one restricts to
the high-$\ell$ regime for $C_{\ell}^{\mathrm{kSZ}^{2}\times\delta_{g}}$, the
convolution leads to a mixing of angular scales, and it is not sufficient to
only consider the high-$k_{1},k_{2},k_{3}$ limit. Note that here,
$(\mathbf{k_{1}},\mathbf{k_{2}},\mathbf{k_{3}})\equiv(\mathbf{q},-\mathbf{j}-\mathbf{q},\mathbf{j})$.

  * •

In fact, when restricted to the high-$\ell$ regime of
$C_{\ell}^{\mathrm{kSZ}^{2}\times\delta_{g}}$, Eq.(8) only dictates that
$k_{3}$ must be high (where $k_{3}\equiv j\approx\ell/\eta$). Consequently, we
expect that two types of triangle configurations of
$B_{p_{\mathbf{\hat{n}}}p_{\mathbf{\hat{n}}}\delta}(\mathbf{k_{1}},\mathbf{k_{2}},\mathbf{k_{3}})$
contribute significantly to the overall estimator at small scales: (1)
Triangles whose all 3 sides are large i.e. high-$k_{1},k_{2},k_{3}$, and (2)
Squeezed triangles whose one side is small and two sides (including
$\mathbf{k_{3}}$) are large.

Both of the above considerations imply that even for astrophysical
applications of the projected-fields estimator (e.g. constraining the baryon
profile, which largely affects high-$\ell$), it is important to account for
all triangle configurations of
$B_{p_{\mathbf{\hat{n}}}p_{\mathbf{\hat{n}}}\delta}$ (especially the squeezed
triangles), and not just those with high-$k_{1},k_{2},k_{3}$. In the absence
of an assumed limiting scale regime,

  * •

A generalized expression must be used for the leading order term of
$B_{p_{\perp}p_{\perp}\delta}$. Assuming only the Limber approximation, we re-
derived the improved expression for $P_{p_{\perp}}$ as given in [46] (Eq. 6
therein). After this verification, we analogously derive a rigorous expression
for the ‘usual’ term’s contribution to $B_{p_{\perp}p_{\perp}\delta}$ (denoted
by $B_{p_{\perp}p_{\perp}\delta}^{\mathrm{usual}}$):

$\int\frac{d^{3}k}{\left(2\pi\right)^{3}}\left[\sqrt{1-\mu_{1}^{2}}\sqrt{1-\mu_{2}^{2}}\right]P_{vv}(k)B_{\delta\delta\delta}(\mathbf{k_{1}}-\mathbf{k},\mathbf{k_{2}}+\mathbf{k},\mathbf{k_{3}}).$
(14)

Comparing this to the corresponding expression for D05 in Eq.(13), we see that
only the geometrical factors (in square brackets) are different, while the
rest is the same. Mathematically, our rigorously derived factor
$\left[\sqrt{1-\mu_{1}^{2}}\sqrt{1-\mu_{2}^{2}}\right]$ is always less than or
equal to D05’s $[1-\mu_{1}\mu_{2}]$, with equality occurring only when
$\mu_{1}=\mu_{2}$.

  * •

We also include the contributions to $B_{p_{\perp}p_{\perp}\delta}$ from the
three non-zero ‘extra’ terms.

The $B_{p_{\perp}p_{\perp}\delta}^{\mathrm{extra-1}}$ and
$B_{p_{\perp}p_{\perp}\delta}^{\mathrm{extra-2}}$ terms are symmetric and are
effectively the same as $B_{p_{\perp}p_{\perp}\delta}^{\mathrm{usual}}$ with
an additional scaling of $(k/k_{1})$ and $(k/k_{2})$, respectively. The
$B_{p_{\perp}p_{\perp}\delta}^{\mathrm{extra-3}}$ term has a different
functional form, and a rough geometric scaling: $([-k+(k_{1}$ or
$k_{2})]/k_{3})$, where either $k_{1}$ or $k_{2}$ appear in the numerator
depending on how the expression is simplified (both ways are equivalent).
Unlike D05’s claim, the extra$-3$ term does not have a pure $(k/k_{3})$
scaling, since it is associated with the $\delta$ mode (while extra$-1$ and
extra$-2$ terms are associated with $p_{\mathbf{\hat{n}}}$).

We refer the reader to Appendix A for details of our rigorous derivation and a
comparison of the four different terms’ contributions in an illustrative case.
Overall, our improved model is given by
$B_{p_{\mathbf{\hat{n}}}p_{\mathbf{\hat{n}}}\delta}=(1/2)\left[B_{p_{\perp}p_{\perp}\delta}^{\mathrm{usual}}+B_{p_{\perp}p_{\perp}\delta}^{\mathrm{extra-1}}+B_{p_{\perp}p_{\perp}\delta}^{\mathrm{extra-2}}+B_{p_{\perp}p_{\perp}\delta}^{\mathrm{extra-3}}\right]$,
where the expressions for the four terms are given by Eqs. (A.2, 33, 34, and
A.3). Since our improved model for
$B_{p_{\mathbf{\hat{n}}}p_{\mathbf{\hat{n}}}\delta}$ does not make any further
approximations, it accurately predicts the theoretical signal
$C_{\ell}^{\mathrm{kSZ}^{2}\times\delta_{g}}$ at all scales of interest, which
allows us to use the projected-fields estimator for extracting cosmological
information.

##  III Methodology

###  III.1 Survey Specifications

We now review the specifications of the CMB and LSS experiments that we have
assumed in our analysis. The projected-fields estimator requires a CMB map
that has been ‘cleaned’ or separated from the other contaminating foregrounds
in the microwave sky using some multi-frequency component-separation
technique. Following the analysis in [35, 6] for the joint data from Planck
[38] and WMAP [48], we take the measured CMB temperature map that has been
cleaned [49] using the “local-generalized morphological component analysis”
(LGMCA) technique [50]. This map has a Gaussian beam $b(\ell)$ with
$\theta_{\mathrm{FWHM}}$ = 5 arcmin, and a white noise power spectrum
$N_{\ell}$ given by

$N_{\ell}=\Delta_{T}^{2}\,b(\ell)^{-2}=\Delta_{T}^{2}\,\mathrm{exp}\left(\ell(\ell+1)\frac{\theta_{\mathrm{FWHM}}^{2}}{8\,\mathrm{ln}\,2}\right),$
(15)

where the “effective” pixel noise level $\Delta_{T}=47\mu$K-arcmin includes
noise from residual foregrounds. In order to obtain the Wiener filter
$F(\ell)=C_{\ell}^{\mathrm{kSZ}}/C_{\ell}^{\mathrm{tot}}$ (Eq. 4), we compute
the theoretical kSZ power spectrum using a
template111https://github.com/nbatta/SILC/blob/master/data/
ksz_template_battaglia.csv derived from cosmological hydrodynamical
simulations [51]. Here, $C_{\ell}^{\mathrm{tot}}$ =
$\left(C_{\ell}^{TT}+C_{\ell}^{\mathrm{kSZ}}+N_{\ell}\right)$, where the
lensed primary CMB temperature power spectrum $C_{\ell}^{TT}$ is calculated
for our fiducial model using CAMB [52, 53]. The filter for Planck is computed
assuming its effective noise (Eq. 15), with an $\ell_{\mathrm{\max}}=3000$.

On the other hand, the upcoming Simons Observatory (SO) and the subsequent
CMB-S4 experiments will map the CMB at a high resolution of
$\theta_{\mathrm{FWHM}}\approx 1.4$ arcmin (or $\ell_{\mathrm{\max}}\approx
8000$). We forecast for their cleaned CMB maps that will be obtained using the
standard “Internal Linear Combination” (ILC) component-separation technique
[54, 55], which minimizes the total variance. While some earlier forecasts for
the kSZ (e.g.[13]) and the projected-fields estimator [32, 6] assumed
instrumental noise only, we use the simulated post-ILC noise curves for SO [1]
and CMB-S4 [2] respectively, which are publicly available
online222https://github.com/simonsobs/so_noise_models/tree/master/
LAT_comp_sep noise/v3.1.0; standard
ILC:deproj-0,333https://sns.ias.edu/jch/S4_190604d_2LAT_Tpol_default
_noisecurves.tgz. These ‘ILC’ noise power spectra include the detectors’ white
noise as well as noise due to any residual foregrounds present after component
separation, and are thus considerably larger (e.g. see Fig. 4 in [56]).

Additionally, we also present forecasts assuming detector-only noise for
CMB-S4 (denoted as ‘CMB-S4 (Det)’ henceforth) in Section V.3 to demonstrate
the major impact of residual foregrounds on these forecasts, and to connect
them with previous works. Throughout this paper, unless mentioned otherwise,
CMB-S4 refers to the CMB-S4 (ILC) case. We construct the standard Wiener
filters (Eq. 4) for SO and CMB-S4 using these noise curves, and the same
theoretical $C_{\ell}^{TT}$ and $C_{\ell}^{\mathrm{kSZ}}$ as for Planck above.
Our filters for these CMB experiments are shown in Figure 1, and their
specifications are summarized in Table 2.

CMB | $\theta_{\mathrm{FWHM}}$ | Noise $\Delta_{T}$ | Noise | $\ell_{\mathrm{\max}}$ | $f_{\mathrm{sky}}$
---|---|---|---|---|---
experiment | [arcmin] | [$\mu$K-arcmin] | description |  |
Planck | 5 | 47 | effectivea | 3000 | 0.7
SO | 1.4 | ILC | realisticb | 8000 | 0.4
CMB-S4 (Det) | 3 | 1 | detector-only | 8000 | 0.4
CMB-S4 (ILC) | 1.4 | ILC | realisticb | 8000 | 0.4
Table 2: Specifications for the forecasted CMB experiments. We assume the
aeffective noise (following [35, 6]) and bpost-ILC noise [1] for Planck and SO
respectively, which include residual foreground contributions and are thus
realistic. We consider two distinct noise models for CMB-S4: post-ILC noise
[2], and detector-only noise (as assumed for case 3 in [6]). Unless mentioned
otherwise, CMB-S4 refers to the CMB-S4 (ILC) case. Figure 1: Wiener filters
(Eq. 4) $F(\ell)$ for the CMB experiments studied in this work, before
multiplying with their corresponding beams as described in Table 2. The
normalization here is arbitrarily chosen such that each $F(\ell)$ peaks at 1,
and it does not affect any of the results.

We consider galaxies as a tracer of the LSS (Eq. 3) and include forecasts for
two different photometric galaxy surveys. The WISE mission [33] imaged the
entire sky in the infrared, collecting a catalog of more than 500 million
objects. We use the same criteria as [57] to select a sample of galaxies,
which were originally studied in [58]. Following the data analysis in [6, 35],
we similarly assume a mask with $f_{\mathrm{sky}}=0.447$ and 46.2 million
selected galaxies. The redshift distribution of these galaxies peaks at
$z\approx 0.3$ and extends up to $z=1$ [59] (as plotted in Fig. 10 of [6]).

The future VRO (formerly named LSST) [60] will take deep optical images of
about half of the sky, creating catalogs of an unprecedented number of
objects. The high-SNR “gold” sample of LSST (i.e. with $i<25.3$) will contain
around 4 billion galaxies over 20,000 deg2, based on empirical estimates from
previous surveys. We assume the predicted redshift distribution of these
galaxies that extends from $z=0$ to $z=3$, with a peak at $z\approx 0.7$ [34]
(Eq.(3.8) therein with $z_{0}\approx 0.3$). We refer the reader to Section 3.7
of [34] for further details.

We assume a constant fiducial linear galaxy bias $b_{g}=1.13$ for WISE
galaxies, the best-fit value determined from Planck (2015) CMB lensing maps
[6] and the Planck$\times$WISE projected-fields measurement [6, 35]. We allow
a redshift-dependent fiducial $b_{g}$ for VRO galaxies [34] given by the model
$b_{g}(z)=1+0.84z$, which has been estimated from simulations in [61].

###  III.2 Total Covariance Matrix

Our parameter forecasts depend crucially on the covariance matrix of
measurement of the projected-fields estimator. We follow [31, 6] and assume
Gaussianity, so that the total covariance matrix $M_{\ell\ell^{^{\prime}}}$ is
diagonal and written as:

$M_{\ell\ell^{^{\prime}}}=\frac{\delta_{\ell\ell^{^{\prime}}}}{(2\ell+1)f_{\mathrm{sky}}}\left[C_{\ell}^{\bar{T}^{2}\bar{T}^{2},f}C_{\ell}^{\delta_{g}\delta_{g}}+\left(C_{\ell}^{\mathrm{kSZ}^{2}\times\delta_{g}}\right)^{2}\right].$
(16)

Here, $f_{\mathrm{sky}}$ is the sky fraction over which the CMB and LSS
experiments overlap (with values given in Table 2), and for
$C_{\ell}^{\bar{T}\bar{T},f}\equiv f^{2}(\ell)C_{\ell}^{\mathrm{tot}}$,

$C_{\ell}^{\bar{T}^{2}\bar{T}^{2},f}\approx
2\int\frac{d^{2}\mathbf{L}}{(2\pi)^{2}}C_{\ell}^{\bar{T}\bar{T},f}C_{|\mathbf{\ell}-\mathbf{L}|}^{\bar{T}\bar{T},f}.$
(17)

The projected galaxy density power spectrum is: [43, 62]

$C_{\ell}^{\delta_{g}\delta_{g}}=\int_{0}^{\eta_{\max}}\frac{d\eta}{\eta^{2}}\,[W^{g}(\eta)]^{2}P_{\delta\delta}\left(k=\frac{\ell}{\eta},\eta\right)+\frac{1}{\bar{n}},$
(18)

where $\bar{n}$ is the projected number of galaxies per steradian, appearing
in the shot noise term $(1/\bar{n})$. The first term in brackets in Eq.(16) is
the contribution due to the actual measurement which includes the instrumental
(or ILC) noise of the CMB and LSS experiments. Meanwhile, the second term is
the corresponding cosmic variance \- an inherent statistical uncertainty in
the measurement that arises because only one realization among all possible
universes described by a certain model can be observed.

###  III.3 Numerical Implementation and the Effect of Massive Neutrinos

In order to calculate the projected-fields estimator
$C_{\ell}^{\mathrm{kSZ}^{2}\times\delta_{g}}$ given by Eq.(8), we equivalently
compute the LOS integral in the redshift $(z)$ variable, where
$\eta(z)\equiv\int_{0}^{z}dz[c/H(z)]$. We discretize the integral and use a
coarse bin size $\Delta z=0.1$ in the window functions $W^{g}(z)$ for the WISE
and VRO photometric surveys. The triangle power spectrum $\mathcal{T}$ (Eq. 9)
is calculated independently in each redshift bin, using the corresponding CMB
filter $f(\ell)$ and the theoretical hybrid bispectrum
$B_{p_{\mathbf{\hat{n}}}p_{\mathbf{\hat{n}}}\delta}$. We now summarize our
numerical computation of $B_{p_{\mathbf{\hat{n}}}p_{\mathbf{\hat{n}}}\delta}$
assuming our _improved_ model Eq.(36), which involves the computation of the
four terms given in Eqs. (A.6, A.8-10).

As noted in Section II.2, we follow [32] and subsequent works and assume the
simulation-based phenomenological substitution suggested by [45], in which we
substitute the non-linear matter density field into Eq.(10) for every velocity
factor appearing in the improved model equations. The resulting non-linear
matter power spectrum $P_{\delta\delta}$ is calculated [63] using CAMB [52,
53] for our assumed fiducial model, defined by:
$\\{H_{0},\Omega_{b}h^{2},\Omega_{c}h^{2},10^{9}A_{s},n_{s},\tau_{re}\\}+\\{\Sigma
m_{\nu}\\}$. We use the best-fit Planck-2018 values (Table 1 of [38]) for the
first six, which are the base parameters in $\Lambda$CDM cosmology: Hubble
constant, baryon density, cold dark matter density, amplitude and slope of the
primordial spectrum of metric fluctuations, and optical depth to reionization,
respectively.

We also consider an extension to this base model by varying the sum of
neutrino masses ($\Sigma m_{\nu}$) in Section V.2, which has been confirmed to
be non-zero by detections of flavor oscillations in solar and atmospheric
neutrinos (e.g. [64, 65]). These detections allow for a normal hierarchy of
neutrinos with a minimum $\Sigma m_{\nu}$ of $\sim 58$ meV, as well as an
inverted hierarchy ($m_{3}\ll m_{1}$ $\simeq$ $m_{2}$) with a minimum $\Sigma
m_{\nu}$ of $>100$ meV [66, 40]. We assume a normal hierarchy for the
neutrinos ($m_{1}$ $\simeq$ $m_{2}$ $\simeq$ 0 $\ll$ $m_{3}$) and a fiducial
value of 60 meV. Aside from affecting the expansion rate, importantly, massive
neutrinos suppress $P_{\delta\delta}$ at scales smaller than their free-
streaming scale, defined as
$k_{\mathrm{fs}}(z)=0.018\,\Omega_{m}^{1/2}(z)[m_{\nu}$/1 eV] h Mpc-1.(e.g.
[39, 40]).

The $C_{\ell}^{\mathrm{kSZ}^{2}\times\delta_{g}}$ estimator contains two
factors of the peculiar velocity field of the LSS, which is largely dominant
at linear and quasi-linear scales. After the phenomenological substitution
[45] using Eq.(10), the overall signal is proportional to the square of the
linear growth rate $f$ of matter perturbations. The presence of massive
neutrinos induces a small scale dependence in $f$ that is very well-modeled by
the fitting function [67]: $f(k,z)\approx\mu(k)\Omega_{m}^{\alpha}(z)$ with
$\alpha\approx 0.55$ for $\Lambda$CDM, and

$\mu(k)=1-A(k)\Omega_{\Lambda}f_{\nu}+B(k)f_{\nu}^{2}-C(k)f_{\nu}^{3},$ (19)

where $f_{\nu}=\Omega_{\nu}/\Omega_{m}$ is the fractional contribution of
neutrinos to the total matter density ($\approx$ 0.0045 in fiducial model),
and $\Omega_{\Lambda}$ is the dimensionless energy density of dark energy
($\approx 0.6847$ in fiducial model). The coefficients $A(k),B(k),C(k)$ are
obtained by interpolating the values given in Table II of [67].

To compute the non-linear matter bispectrum $B_{\delta\delta\delta}$, we
follow previous works [6, 35, 36] and use the fitting function from [68],
which has the form:

$B_{\delta\delta\delta}(k_{1},k_{2},k_{3})=\sum_{cyc}2F_{2}(k_{1},k_{2},k_{3})P_{\delta\delta}(k_{1})P_{\delta\delta}(k_{2}),$
(20)

where $\sum_{cyc}$ denotes a cyclic sum over $(k_{1},k_{2},k_{3})$ and
$P_{\delta\delta}$ is again the non-linear matter power spectrum. The kernel
$F_{2}$’s scale and redshift dependence is governed by parameters that are
fitted using $\Lambda$CDM-only N-body simulations. However, this formula is
likely general enough to be applicable at the first-order for certain
$\Lambda$CDM-extensions since the kernel is weakly dependent on cosmology
[68]. To compute the matter bispectrum, we use the matter power spectrum in
the presence of massive neutrinos in the fitting function above (Eq. 20). At
first-order, this is consistent with results from second-order perturbation
theory [69] and simulations [70, 71] which found that the suppression in the
bispectrum’s amplitude is almost twice that of the power spectrum.

Similar to previous works [35, 6, 36], our numerical computation here
implicitly assumes that baryons trace the dark matter distribution (where we
have dropped the subscripts on overdensities $\delta$), which is true at large
scales ($\sim$Mpc), but not at smaller $\sim$kpc scales (e.g. [25, 72]).
Alternatively, the matter bispectrum could be computed within the Halo Model
framework [73], as done recently in [37] assuming the D05 model and using
their code
class_sz444https://github.com/borisbolliet/class_sz/releases/tag/v1.0.0. We
discuss this and other possible future directions in Section V.

The direct numerical evaluation of
$C_{\ell}^{\mathrm{kSZ}^{2}\times\delta_{g}}$ with our improved model is even
more time-consuming than the D05 model. While [37] used Fourier transforms to
accelerate the approximate computation, the same cannot be done for our model
since the integrands for $B_{p_{\mathbf{\hat{n}}}p_{\mathbf{\hat{n}}}\delta}$
are not ‘separable’. In contrast, we speed up by creating emulators for
$B_{p_{\mathbf{\hat{n}}}p_{\mathbf{\hat{n}}}\delta}$ in each redshift bin
using the ostrich library555https://github.com/dylancromer/ostrich [74], which
uses Principal Component Analysis and Gaussian Process Interpolation. The
resulting computation of a $C_{\ell}^{\mathrm{kSZ}^{2}\times\delta_{g}}$ is
accurate and relatively efficient, taking $\mathcal{O}$(10 s) on a laptop.

###  III.4 CMB Lensing Contribution

The projected-fields estimator involves squaring a (cleaned and) filtered CMB
temperature map. As a result, the actual measurement gets a contribution due
to the weak lensing of the CMB by intervening LSS, in addition to the kSZ
signal $C_{\ell}^{\mathrm{kSZ}^{2}\times\delta_{g}}$ modeled in Section II.
The leading-order lensing contribution is given by [6]:

$\displaystyle C_{\ell}^{\mathrm{lens}}\approx\frac{-2\ell
C_{\ell}^{\psi\delta_{g}}}{\left(2\pi\right)^{2}}\int_{0}^{\infty}dL^{{}^{\prime}}L^{{}^{\prime}2}f(L^{{}^{\prime}})C_{L^{{}^{\prime}}}^{TT}\int_{0}^{2\pi}d\phi
f(|\mathbf{L}^{{}^{\prime}}+\mathbf{\ell}|)\,\mathrm{cos}\phi,$ (21)

where $C_{\ell}^{\psi\delta_{g}}$ is the cross-power spectrum between the
lensing potential $\psi$ (see [75] for a review of CMB lensing) and the tracer
field, and $C_{\ell}^{TT}$ is the unlensed primary CMB power spectrum. The
lensing convergence, $\kappa(\theta)=-$$\nabla^{2}\psi(\theta)$, so that in
Fourier space,
$C_{\ell}^{\psi\delta_{g}}=(2/\ell(\ell+1))C_{\ell}^{\kappa\delta_{g}}$. We
thus compute $C_{\ell}^{\psi\delta_{g}}$ using the Limber approximation [43,
62]:

$C_{\ell}^{\kappa\delta_{g}}=\int_{0}^{\infty}\frac{d\eta}{\eta^{2}}\,W^{g}(\eta)W^{\kappa}(\eta)P_{\delta\delta}\left(k=\frac{\ell}{\eta},\eta\right),$
(22)

where $W^{\kappa}(\eta)$ is the CMB lensing kernel (e.g. [76]).

As seen from Eq.(21), $C_{\ell}^{\mathrm{lens}}$ vanishes if the beam is
infinitesimally small and no filter is applied [6]. For the projected-fields
measurements with Planck data [35, 36], the lensing contribution is much
higher than the kSZ signal [6]. In contrast, SO and CMB-S4 will have a much
smaller beam, and will cover $\ell\sim 3000-8000$ where the kSZ effect
dominates over the anisotropies due to lensing. Thus, unlike Planck,
$C_{\ell}^{\mathrm{kSZ}^{2}\times\delta_{g}}$ will dominate over
$C_{\ell}^{\mathrm{lens}}$ contribution in future measurements with SO and
CMB-S4, as shown in the left panel of Figure 2.

##  IV Results

Figure 2: (left panel): Projected-fields estimator signal as predicted by the
approximate ‘D05’ model [32] ($C_{\ell}^{\mathrm{approx}}$; dashed black
lines) and our improved model ($C_{\ell}^{\mathrm{improved}}$; solid colored
lines), for (top-to-bottom) Planck, SO, and CMB-S4(ILC) data cross-correlated
with WISE galaxies, and CMB-S4(ILC) data with VRO galaxies, respectively. The
$C_{\ell}^{\mathrm{lens}}$ contributions (dot-dashed lines) are shown for all
experiments except Planck$\times$WISE, for which $C_{\ell}^{\mathrm{lens}}$ is
much higher than the $C_{\ell}^{\mathrm{kSZ}^{2}\times\delta_{g}}$ [6].
(right panel): Ratios of the two models’ predicted signals,
($C_{\ell}^{\mathrm{improved}}/C_{\ell}^{\mathrm{approx}}$; solid colored
lines). The shaded regions show ratios between the statistical uncertainties
(Eq. 16) and the corresponding $C_{\ell}^{\mathrm{approx}}$, with a linear
binning of $\Delta\ell=289$. The relatively large statistical errors for
Planck$\times$WISE are beyond the axis range shown here [35], and are thus
omitted.

In this section, we compare the predictions of our rigorous improved model and
the approximate D05 model, for WISE galaxies ($\delta_{g}$) cross-correlated
with Planck, SO and CMB-S4 data, and for VRO galaxies cross-correlated with
CMB-S4 data. $C_{\ell}^{\mathrm{kSZ}^{2}\times\delta_{g}}$ assuming our
improved model is numerically computed as described in Section III.3. The
approximate signal is evaluated similarly, using Eq.(11) instead for the
hybrid bispectrum, where the matter bispectrum and non-linear
$v_{\mathrm{rms}}^{2}$ are consistently computed with the same assumptions as
above.

The first plot in the left panel of Figure 2 compares the models’ predictions
for Planck$\times$WISE, while the plot on the right shows the ratio between
the two. The differences between the predicted signals are up to a $\sim 5\%$
level and are thus negligible as compared to the large statistical errors for
the Planck$\times$WISE measurement [35]. We note that the numerical
simulations and the approximate model up to an $\ell_{\mathrm{max}}$ of
$\approx 3000$ (i.e. Planck’s resolution) agree to within $\sim$5-10$\%$ in
[6]. Thus, the differences we find from the improved model would not be
resolved in previous theory comparisons or impact previous measurements with
Planck [35]. Interestingly, the improved model predicts a higher signal for
most $\ell<3000$ and could explain some (but not all) of the excess
$A_{\mathrm{kSZ}^{2}}$ found for the green unWISE sample in [36].

On the other hand, the three lower rows of plots in Figure 2 imply that our
improvements in the theoretical model lead to $\sim 15\%$ level differences in
the predicted signal for the upcoming SO and CMB-S4 experiments, which have
enormous forecasted SNRs (see [6, 37], and our Section 23). The right panel of
Figure 2 shows that these differences are significantly larger than the
associated statistical uncertainties for these measurements (given by Eq. 16),
and hence must be taken into account.

From Figure 2, we also see that our improvements in the theoretical model have
the largest impact on the predicted
$C_{\ell}^{\mathrm{kSZ}^{2}\times\delta_{g}}$ at both the low-$\ell$ and
high-$\ell$ ends, despite the arguments in [31, 32] justifying the use of the
approximate model at small scales. From Section II.3, the ‘usual’ term in our
model for $B_{p_{\perp}p_{\perp}\delta}$ is always $\leq
B_{p_{\perp}p_{\perp}\delta}^{\mathrm{approx}}$, while the other 3 additional
terms may or may not contribute comparably depending on the triangle
configuration.

Therefore, the ratios of the overall
$C_{\ell}^{\mathrm{kSZ}^{2}\times\delta_{g}}$ shown here can be less than or
greater than 1 at different scales, and differ significantly at high-$\ell$s
due to the convolution in Eq.(9) and because squeezed triangles beyond the
high-$k_{1},k_{2},k_{3}$ limit contribute significantly to the signal.
Moreover, since these differences in the signal are scale-dependent, the
corresponding change in the shape of
$C_{\ell}^{\mathrm{kSZ}^{2}\times\delta_{g}}$ would also likely impact
parameter inferences. In summary, it is crucial to adopt our improved model in
_any_ analysis that applies the projected-fields estimator to upcoming high-
resolution CMB experiments, to ensure unbiased results.

##  V Parameter Forecasts

###  V.1 Cosmological Dependence

Not only is our improved model necessary to avoid biased inferences from
future measurements, but it also enables us to accurately study the
cosmological dependence of the $C_{\ell}^{\mathrm{kSZ}^{2}\times\delta_{g}}$
signal for the first time since it is not restricted to small scales. As noted
in Section III.3, this signal depends on the baryon profile, as well as other
astrophysics such as Halo Occupation Distribution (HOD) [37] and non-linear
bias for galaxies. Here, we focus on the cosmological dependence of the signal
and assume a linear galaxy bias $b_{g}$ for simplicity. The other
astrophysical parameters are assumed to be known externally.

So far, the two Planck measurements with this estimator [35, 6, 36] jointly
fit for $b_{g}$ and the amplitude of this kSZ signal, $A_{\mathrm{kSZ}^{2}}$.
It is defined as a single free parameter such that
$C_{\ell}^{\mathrm{kSZ}^{2}\times\delta_{g}}=A_{\mathrm{kSZ}^{2}}C_{\ell}^{\mathrm{kSZ}^{2}\times\delta_{g}}|_{\mathrm{fid}}$;
its fiducial value is 1. From Eq.(8),
$A_{\mathrm{kSZ}^{2}}\propto(f_{b}f_{\mathrm{free}})^{2}$, where $f_{b}$ and
$f_{\mathrm{free}}$ are the baryon and free electron fractions at the redshift
of the LSS tracers. While we use a redshift-dependent galaxy bias for VRO
(Section III.1), here, we vary the bias amplitude $b_{g}$ which can be
redefined so that its fiducial value matches that of WISE galaxies ($\approx
1.13$). The total measured signal is
$C_{\ell}^{\mathrm{tot}}=(C_{\ell}^{\mathrm{kSZ}^{2}\times\delta_{g}}+C_{\ell}^{\mathrm{lens}})$.
Since $C_{\ell}^{\mathrm{lens}}\propto b_{g}$ and
$C_{\ell}^{\mathrm{kSZ}^{2}\times\delta_{g}}\propto
b_{g}A_{\mathrm{kSZ}^{2}}$, there is a degeneracy between
$A_{\mathrm{kSZ}^{2}}$ and $b_{g}$, and some of it is broken by the
$C_{\ell}^{\mathrm{lens}}$ term [6].

We forecast the 1$\sigma$ statistical uncertainty on $A_{\mathrm{kSZ}^{2}}$,
denoted by $\Delta A_{\mathrm{kSZ}^{2}}$ $(=\Delta
A_{\mathrm{kSZ}^{2}}/A_{\mathrm{kSZ}^{2}}|_{\mathrm{fid}})$, for future
measurements with SO$\times$WISE, CMB-S4$\times$WISE, and CMB-S4$\times$VRO,
assuming post-ILC CMB noise. We use the Fisher matrix formalism (e.g.[77]),
where the Fisher matrix is given by:

$F_{ij}=\sum_{\ell,\ell^{{}^{\prime}}}\frac{\partial{C_{\ell}^{\mathrm{tot}}}}{\partial{\theta_{i}}}[M^{-1}]_{\ell\ell^{^{\prime}}}\frac{\partial{C_{\ell^{{}^{\prime}}}^{\mathrm{tot}}}}{\partial{\theta_{j}}},$
(23)

and the covariance matrix of the measurement $M_{\ell\ell^{^{\prime}}}$ is
given by Eq.(16). For this Fisher matrix, the marginalized $1\sigma$
uncertainty on each parameter $\theta_{i}$ is given by $\sqrt{[F^{-1}]_{ii}}$.

As in previous measurements, we first assume that the cosmology is fixed, so
the parameter set is
$\boldsymbol{\theta}_{\mathrm{base}}=\\{A_{\mathrm{kSZ}^{2}},b_{g}\\}$. Here,
partial derivatives in Eq.(23) are calculated directly due to the linear
dependence. We also set a 1$\%$ Gaussian prior on $b_{g}$, since it can be
determined by galaxy auto-correlation ($C_{\ell}^{\delta_{g}\delta_{g}}$) and
galaxy-lensing cross-correlation ($C_{\ell}^{\kappa\delta_{g}}$) external
measurements (e.g. [6, 78]). Figure 3 shows the resulting 68$\%$ C.L.
contours. Marginalized 1$\sigma$ uncertainties on $A_{\mathrm{kSZ}^{2}}$ are
at the sub-percent level for all 3 survey combinations (Table 3: top row),
corresponding to an SNR $(\equiv A_{\mathrm{kSZ}^{2}}/\Delta
A_{\mathrm{kSZ}^{2}})$ of about 113, 127, and 182 for SO$\times$WISE,
CMB-S4$\times$WISE, and CMB-S4$\times$VRO, respectively.

Figure 3: Contours show 68$\%$ confidence levels ($1\sigma$) for $\\{A_{\mathrm{kSZ}^{2}},b_{g}\\}$ for forecasted projected-fields kSZ measurements with SO$\times$WISE (green), CMB-S4$\times$WISE (blue), and CMB-S4$\times$VRO (red). Here, we assume realistic post-ILC noise for SO and CMB-S4, a 1$\%$ Gaussian prior on $b_{g}$, and fixed cosmological parameters. |  $\%\Delta A_{\mathrm{kSZ}^{2}}$: 1$\sigma$ $\%$ uncertainty on $A_{\mathrm{kSZ}^{2}}$
---|---
| SO $\times$ WISE | CMB-S4(ILC)$\times$ WISE | CMB-S4(ILC)$\times$ VRO
Fixed Cosmology | 0.9 | 0.8 | 0.6
Varying $\Lambda$CDM + Planck prior | 7.1 | 6.7 | 6.9
Varying $\Lambda$CDM + (low-$\ell$ Planck+high-$\ell$ CMB-S4) prior | 5.9 | 5.6 | 5.9
Varying $\Lambda$CDM + (low-$\ell$ LiteBIRD+high-$\ell$ CMB-S4) prior | 3.7 | 3.6 | 3.6
Table 3: Marginalized 1$\sigma$ % uncertainty on $A_{\mathrm{kSZ}^{2}}$ for
the three survey combinations considered, when cosmological $\Lambda$CDM
parameters are: fixed at their best-fit Planck-2018 values [38] (top row),
allowed to vary around their fiducial values with a known Planck prior (second
row), with a forecasted (low-$\ell$ Planck \+ high-$\ell$ CMB-S4) prior (third
row) and with a forecasted (low-$\ell$ LiteBIRD + high-$\ell$ CMB-S4) prior
(bottom row) [79]. In all cases here, we assume an external 1$\%$ Gaussian
prior on $b_{g}$, realistic post-ILC noise, and $\Sigma m_{\nu}$ fixed at 60
meV.

Next, we consider the cosmological dependence of this kSZ signal, and allow
the $\Lambda$CDM parameters to vary about their fiducial Planck-2018 values.
Thus, we compute kSZ Fisher matrices for the parameter set:

$\boldsymbol{\theta}_{\mathrm{min}}=\\{H_{0},\Omega_{b}h^{2},\Omega_{c}h^{2},10^{9}A_{s},n_{s}\\}+\\{A_{\mathrm{kSZ}^{2}},b_{g}\\},$
(24)

where we numerically compute partial derivatives with respect to $\Lambda$CDM
parameters using the same step-sizes as [80] and [81] (Table I therein), and
check for numerical stability. We apply an external prior based on current
Planck measurements of the CMB (TT,TE,EE+lowE; no CMB lensing) [38], obtained
using the code
FishLSS666https://github.com/NoahSailer/FishLSS/tree/master/FishLSS/ [79],
which we marginalize over $\tau_{re}$.

From Table 3, with this Planck prior, the marginalized uncertainties on
$A_{\mathrm{kSZ}^{2}}$ are around $\sim 7\%$ for all 3 experiments. Compared
to the fixed cosmology case, this pronounced increase in $\Delta
A_{\mathrm{kSZ}^{2}}$ arises due to the high sensitivity of
$C_{\ell}^{\mathrm{kSZ}^{2}\times\delta_{g}}$ to cosmology, particularly
$\Omega_{c}h^{2}$ and $A_{s}$, as we explain in detail in Appendix B. To
understand this intuitively, note that Eqs.(8), (10), and (20) imply a rough
scaling for the signal,

$C_{\ell}^{\mathrm{kSZ}^{2}\times\delta_{g}}\propto(faH/\eta)^{2}\,P_{\delta\delta}^{3}.$
(25)

Even if the bracketed growth factor above is neglected,
$P_{\delta\delta}\propto\sigma_{8}^{2-3}$, where $\sigma_{8}$ is the power
spectrum normalization. Hence,
$C_{\ell}^{\mathrm{kSZ}^{2}\times\delta_{g}}\propto\sigma_{8}^{6-7}$, as
suggested roughly in [31]. Moreover, $\sigma_{8}$ is a derived parameter
mainly dependent on $\Omega_{m}$ (and thus $\Omega_{c}h^{2}$) and $A_{s}$.
Current Planck constraints on $\sigma_{8}$ are $\sim 0.9\%$, so this scaling
contributes around $5.7\%$ to $\Delta A_{\mathrm{kSZ}^{2}}$.

We note that the strong cosmological dependence of the signal highlighted
above does not change the detection significance of the kSZ signal, which will
be very high, but it will limit the interpretation of the measured amplitude
in terms of astrophysical quantities (such as baryon fraction).

With a basic CMB prior combining information from low-$\ell$ ($\ell<30$)
Planck data and future high-$\ell$ ($30<\ell<5000$) CMB-S4 data [79], the
forecasted uncertainties on $\Delta A_{\mathrm{kSZ}^{2}}$ improve to $\sim
6\%$ (Table 3). The upcoming LiteBIRD satellite will map the CMB at $\ell<200$
with high sensitivity, allowing a cosmic variance-limited measurement of
$\tau_{re}$ [82]. Applying a forecasted (low-$\ell$ LiteBIRD + high-$\ell$
CMB-S4) prior on $\Lambda$CDM, $\Delta A_{\mathrm{kSZ}^{2}}$ improve further
to $\sim 3.5\%$, since $A_{s}$ and $\tau_{re}$ are highly degenerate in
primary CMB data, and $C_{\ell}^{\mathrm{kSZ}^{2}\times\delta_{g}}$ is
extremely sensitive to $A_{s}$ (Appendix B).

###  V.2 Fisher forecasts including Neutrino Masses

Figure 4: Fisher forecasts for marginalized errors on $\Lambda$CDM parameters,
$\Sigma m_{\nu}$ (in eV here), and $\\{A_{\mathrm{kSZ}^{2}},b_{g}\\}$ from
kSZ-only projected-fields constraints, without assuming _any_ priors, and with
realistic post-ILC noise for SO and CMB-S4. Contours show the 2D marginalized
68$\%$ confidence ellipses for SO$\times$WISE (green), CMB-S4$\times$WISE
(blue), and CMB-S4$\times$VRO (red).

Equipped with our improved model, we now illustrate the potential sensitivity
of this kSZ estimator to extensions of the Standard Model, by studying the
dependence of $C_{\ell}^{\mathrm{kSZ}^{2}\times\delta_{g}}$ on the sum of
neutrino masses. Here, we compute kSZ Fisher matrices for the parameters
$\boldsymbol{\theta}_{\mathrm{ext}}=\boldsymbol{\theta}_{\mathrm{min}}+\\{\Sigma
m_{\nu}\\}$. We note that there is no marginalization over higher-order bias
parameters or astrophysical parameters related to the halo profile performed
here. Therefore, when interpreting these forecasts they should only be taken
as representative of the parameter sensitivity of the kSZ2 observable, and
would only be attainable if the bias and astrophysical parameters were known
externally (for example from a joint analysis with other probes).

At first, we do not apply any external prior. Contours in Figure 4 show the
resulting 2D marginalized 1$\sigma$ uncertainties for each pair of parameters,
thus informing about their degeneracies for
$C_{\ell}^{\mathrm{kSZ}^{2}\times\delta_{g}}$. CMB-S4$\times$WISE gives
tighter constraints than SO$\times$WISE on all parameters except for $b_{g}$,
due to CMB-S4’s lower noise. The relative scale-dependence of
$C_{\ell}^{\mathrm{lens}}$ and $C_{\ell}^{\mathrm{kSZ}^{2}\times\delta_{g}}$
together determine the error on $b_{g}$; the relative amplitude of
$C_{\ell}^{\mathrm{lens}}$ is higher for SO$\times$WISE (Figure 2).

In all the forecasts presented here, we do not account for the cosmological
dependence of $C_{\ell}^{\mathrm{lens}}$, but only its linear dependence on
the galaxy bias. This is a reasonable assumption since a template for
$C_{\ell}^{\mathrm{lens}}$ can be obtained separately by cross-correlating the
galaxy sample with a CMB lensing map. However, as a check, we also performed
an alternative Fisher analysis including the dependence of
$C_{\ell}^{\mathrm{lens}}$ on all $\Lambda$CDM+$\\{\Sigma m_{\nu}\\}$
parameters. In that scenario, forecasted errors on all parameters shift within
a margin of $\sim 10\%$ for all 3 experiments; we omit the exact values here
for brevity.

| kSZ + Planck \+ 1% $b_{g}$ prior
---|---
| SO $\times$ WISE | CMB-S4(ILC)$\times$ WISE | CMB-S4(ILC)$\times$ VRO
$\sigma(\Sigma m_{\nu})$ [meV] | 168 | 100 | 254
$\%\Delta A_{\mathrm{kSZ}^{2}}$ | 10.9 | 8.1 | 20.7
Table 4: Forecasts for marginalized 1$\sigma$ errors on $\Sigma m_{\nu}$ (top
row) and on $A_{\mathrm{kSZ}^{2}}$ (written as a %; bottom row) from the
$C_{\ell}^{\mathrm{kSZ}^{2}\times\delta_{g}}$ kSZ signal jointly constraining
the $\boldsymbol{\theta}_{\mathrm{ext}}$ parameters, after applying our Planck
prior [79] on $\Lambda$CDM parameters only (given in Table 5). We also set an
external 1$\%$ Gaussian prior on $b_{g}$, and assume realistic post-ILC CMB
noise. 1$\sigma$ errors on $\Lambda$CDM parameters improve only marginally
over the Planck prior, and are thus omitted here.

From Figure 4, we show that the forecasted errors on
$\Omega_{c}h^{2},10^{9}A_{s},n_{s}$, and $b_{g}$ are smaller for
CMB-S4$\times$VRO than CMB-S4$\times$WISE, since
$C_{\ell}^{\mathrm{kSZ}^{2}\times\delta_{g}}$ with the former combination is
more sensitive to those parameters than the latter. However,
CMB-S4$\times$WISE constrains the other four parameters:
$\Omega_{b}h^{2},H_{0},\Sigma m_{\nu}$, and $A_{\mathrm{kSZ}^{2}}$ more
tightly. Note that when cosmology is assumed to be fixed (in Section V.1
above), the SNR for CMB-S4$\times$VRO is higher than CMB-S4$\times$WISE, and
the corresponding $\Delta A_{\mathrm{kSZ}^{2}}$ is smaller with VRO (Figure
3).

Once cosmological parameters are freed up, i.e. $\Lambda$CDM in Table 3 and
$\Lambda$CDM+$\\{\Sigma m_{\nu}\\}$ in Table 4, the uncertainty on
$A_{\mathrm{kSZ}^{2}}$ is higher with VRO than with WISE, even though VRO is a
deeper survey. [31] shows that differential contribution to
$C_{\ell}^{\mathrm{kSZ}^{2}\times\delta_{g}}$ comes largely from lower-
redshift bins peaking around $z\sim 0.5$, which reasonably matches the
$z$-range of WISE. For deeper surveys such as VRO that have a long tail at the
high-$z$ end of their distribution, using an optimal redshift-weighting in
$W_{g}$ for $C_{\ell}^{\mathrm{kSZ}^{2}\times\delta_{g}}$ may possibly improve
their SNR [6] and their cosmological constraints. Alternatively, a low-$z$
subsample of the galaxy survey may be chosen to improve such constraints; we
leave this direction for future work.

We now apply the same Planck prior [79] as the previous subsection on only
$\Lambda$CDM parameters, marginalizing over $\tau_{re}$, and include a 1%
Gaussian prior on $b_{g}$. Table 4 shows the resulting 1$\sigma$ %
uncertainties on $\Delta A_{\mathrm{kSZ}^{2}}$, which are higher than those in
Table 3 due to the additional marginalization over $\Sigma m_{\nu}$.
Constraints on $\Lambda$CDM parameters improve only slightly over the Planck
prior. Our results are robust to the choice of the assumed neutrino hierarchy.

We do not include _any_ external prior on $\Sigma m_{\nu}$ here, and forecast
a 1$\sigma$ constraint of $\sim 168$ meV and 100 meV on $\Sigma m_{\nu}$ with
the projected-fields kSZ estimator for SO$\times$WISE and CMB-S4$\times$WISE
respectively. In comparison, [13] forecasted a $\sigma(\Sigma m_{\nu})$ of
$\sim 220$ meV and 96 meV for a SO-like and CMB-S4-like experiment
respectively, by combining with spectroscopic galaxy surveys using the
pairwise-kSZ estimator, where they included a forecasted Planck prior on
$\Lambda$CDM as well as $\Sigma m_{\nu}$. We further discuss our
$\sigma(\Sigma m_{\nu})$ constraints in the context of other probes in the
next subsection.

In our initial analysis here, we demonstrate the potential of this estimator
as a probe of cosmology and of the sum of neutrino masses, through their
effect on the growth of structure. A more accurate analysis would include a
non-linear galaxy bias that also models the induced scale-dependence due to
$\Sigma m_{\nu}$ (e.g. [83, 84]), instead of the simple linear galaxy bias
$b_{g}$ assumed here. As discussed in Section III.3, the matter bispectrum too
could be numerically computed with a more accurate model, that also possibly
accounts for the effect of neutrinos more directly (e.g. [84, 70]).

###  V.3 Impact of Residual Foregrounds

| 1$\sigma$ errors on parameters
---|---
Parameter | Planck | CMB-S4(Det)$\times$WISE
| prior | kSZ-only | kSZ + Planck
$H_{0}$ [km/s/Mpc] | 0.57 | 3.20 | 0.34
$\Omega_{b}h^{2}$ | 0.00013 | 0.00430 | 0.00010
$\Omega_{c}h^{2}$ | 0.0013 | 0.0131 | 0.0008
$10^{9}A_{s}$ | 0.031 | 0.366 | 0.023
$n_{s}$ | 0.003 | 0.030 | 0.003
$\Sigma m_{\nu}$ [meV] | - | 136 | 38
$A_{\mathrm{kSZ}^{2}}$ | - | 0.411 | 0.025
$b_{g}$ | - | 0.011 | 0.007
Table 5: Forecasts for marginalized 1$\sigma$ errors on the
$\boldsymbol{\theta}_{\mathrm{ext}}$ parameters assuming detector-only noise
and the absence of residual foregrounds in CMB-S4 maps. The third column shows
kSZ-only 1$\sigma$ constraints with the
$C_{\ell}^{\mathrm{kSZ}^{2}\times\delta_{g}}$ signal from
CMB-S4(Det)$\times$WISE. We include a Planck prior [79] as described in
Section V.1 on $\Lambda$CDM only (second column), and the resulting kSZ+Planck
constraints are given (last column).

In the absence of accurate post-ILC noise models for SO and CMB-S4, earlier
Fisher forecasts for the kSZ (e.g. [31, 32, 6, 13]) often assumed detector-
only white noise. In order to connect with these previous works, in this
subsection, we forecast for a CMB-S4(Det)$\times$WISE measurement assuming
this noise model for CMB-S4 (Table 2).

Similar to Section V.1, we first keep the cosmology fixed and forecast an SNR
$(\equiv A_{\mathrm{kSZ}^{2}}/\Delta A_{\mathrm{kSZ}^{2}})$ of $\sim 345$,
consistent with the [6] forecast (Table II therein; case 3). Thus, comparing
with the CMB-S4(ILC)$\times$WISE case, we see that the presence of residual
foregrounds in post-ILC CMB maps diminishes the significance of this detection
by a factor of $\sim 3$.

We then vary the $\Lambda$CDM+$\\{\Sigma m_{\nu}\\}$ parameters as well, and
include the same Planck prior [79] as described in Section V.1 on
$\Lambda$CDM. The resulting kSZ+Planck marginalized errors improve upon the
Planck errors significantly, unlike the CMB-S4(ILC)$\times$WISE case. This
again shows the enormous impact of residual contaminants in post-ILC CMB maps.
In the idealistic case of detector-only noise, we forecast a $\sigma(\Sigma
m_{\nu})$ of 38 meV for CMB-S4(Det)$\times$WISE, compared to 100 meV for
CMB-S4(ILC)$\times$WISE. However, the projected-fields
$C_{\ell}^{\mathrm{kSZ}^{2}\times\delta_{g}}$ estimator does need cleaned CMB
maps. Multi-frequency techniques such as ILC remove the bulk of foregrounds,
but leave some residual contamination. Thus, the CMB-S4(ILC)$\times$WISE
forecasts are realistic for this future measurement.

Although we do not include _any_ prior on $\Sigma m_{\nu}$, our realistic
$\sigma(\Sigma m_{\nu})$ forecast is tighter than the 95% C.L. bound of 240
meV from Planck’s extended (TT,TE,EE+lowE)+lensing analysis. Assuming
detector-only noise, a TT,TE,EE+lowE)+lensing analysis for CMB-S4
optimistically forecasts a $\sigma(\Sigma m_{\nu})=73$ meV (e.g. [85]), which
is again insufficient alone for a significant $\Sigma m_{\nu}$ detection.
Combining it with upcoming BAO measurements from DESI [86] or with clustering
and shear data from VRO, would allow a 2-3$\sigma$ significant detection of
the normal hierarchy (e.g. [85, 66, 40]). The
$C_{\ell}^{\mathrm{kSZ}^{2}\times\delta_{g}}$ kSZ signal is a complementary
probe of $\Sigma m_{\nu}$ and can possibly be combined with these LSS probes
to improve upon their constraints.

Moreover, note that this estimator only requires the removal of foregrounds
that are correlated with the considered LSS tracers [6], so for WISE galaxies,
higher-$z$ foregrounds such as the CIB (whose bulk emission comes from
$z\gtrsim 1$ [87]) may be neglected. Also, [36] used a different
‘$\alpha$-cleaning’ method which maximizes their SNR. Therefore, such
alternative cleaning techniques that are more suited for the
$C_{\ell}^{\mathrm{kSZ}^{2}\times\delta_{g}}$ estimator may give parameter
constraints that are intermediate between the two CMB-S4(Det) and CMB-S4(ILC)
cases that we have considered here.

###  V.4 Future Directions

Apart from the future directions described earlier in this section, there are
a few key ways to build upon our analysis here. As noted in Section III.3,
baryon and dark matter density profiles are found to differ at small scales
($\lesssim$ kpc) due to feedback and energy injection, and neglecting this can
introduce systematics in the future high-resolution regime. Recently, [37]
forecasted constraints on the baryon profile using
$C_{\ell}^{\mathrm{kSZ}^{2}\times\delta_{g}}$, although they are limited in
the sense that they assume the approximate D05 model and a fixed cosmology.
Incorporating our improved model within the Halo Model approach of [37] that
uses a more realistic baryon profile (e.g. [88]) would thus enable improved
forecasts. The baryon profile can also be determined externally, for example,
with kSZ measurements using velocity-weighted stacking [25].

Another systematic is the galaxy bias that becomes non-linear at small scales,
while $\Lambda$CDM extensions such as massive neutrinos induce an additional
scale-dependence in it [83, 84]. Thus, a non-linear galaxy bias should be
included in future analyses for further accuracy. Alternatively, this effect
could be consistently accounted for via HOD parameters within the Halo Model
framework.

##  VI Conclusions

In this work, we have considered the projected-fields kSZ estimator which has
the distinct advantage of not requiring accurate spectroscopic redshifts. The
original formulation of this estimator [31, 32] has an approximate model (D05)
only applicable at small scales, which has been used to estimate the baryon
fraction from the signal amplitude $A_{\mathrm{kSZ}^{2}}$ with Planck and WISE
data [35, 36]. Here, we have presented a first rigorous derivation of an
_improved_ theoretical model for this estimator, that is thus accurate even at
larger scales where cosmological effects are significant. While the
differences between the predicted signal from our model and the D05 model are
negligible for previous Planck detections, they are significant at all scales
($\sim$10-15$\%$) for upcoming measurements with SO and CMB-S4 that are
forecasted to have high SNR ($>100$). Given the scale-dependence of these
differences, it is crucial to adopt our improved model in _any_ future
analysis with this estimator to ensure unbiased results.

Our improved model enables us to accurately study the cosmological dependence
of the $C_{\ell}^{\mathrm{kSZ}^{2}\times\delta_{g}}$ signal. While the
marginalized uncertainty on $A_{\mathrm{kSZ}^{2}}$ assuming fixed cosmology is
at the sub-percent level for future measurements, it increases to $\sim 7\%$
when $\Lambda$CDM parameters are freed up with a Planck prior, limiting the
signal’s interpretation. The detection significance remains high, and $\Delta
A_{\mathrm{kSZ}^{2}}$ will improve with reduced uncertainty on $\Lambda$CDM
parameters (especially $\tau_{re}$) from CMB-S4 and LiteBIRD. As an
illustration of the potential sensitivity of this signal to $\Lambda$CDM
extensions, our initial analysis forecasts a $\sigma(\Sigma m_{\nu})$ of $\sim
168$ meV and 100 meV for SO and CMB-S4 (with WISE, and assuming realistic CMB
noise), similar to corresponding constraints from pairwise-kSZ [13]. This
reaffirms the kSZ as a complementary probe of neutrino masses among other CMB
and LSS probes.

While it is important to marginalize over cosmology in any analysis with this
estimator, the numerical computation using the improved model will need to be
sped up for a full MCMC analysis in the future. Building upon our current
method of emulators, a faster implementation may be possible within the
cosmopower emulation framework [89, 90] in the future. Our work is essential
for future extended analyses jointly fitting both baryonic astrophysics and
cosmology with this kSZ estimator. This could be achieved by incorporating our
improved model and its cosmological dependence within the Halo Model framework
of [37]. We note that the formalism of this estimator, further developed here,
is directly applicable for analyses with any tracers of the underlying matter
density, such as quasars, weak lensing convergence (e.g. [37]), or 21-cm
fluctuations probing patchy reionization (e.g. [41, 10]).

###### Acknowledgements.

We thank Boris Bolliet, J. Colin Hill, Aleksandra Kusiak, Mathew
Madhavacheril, and Anirban Roy for their helpful feedback and suggestions on
this work. We thank Emily Moser and Dylan Cromer for their help in speeding up
our computation using the ostrich emulator library. R.P and N.B acknowledge
the support from NSF grant AST-1910021 and NASA grant 17-ATP17-0060 for this
work and N.B acknowledges additional support from NASA grants 21-ADAP21-0114
and 21-ATP21-0129. S.F. is supported by the Physics Division of Lawrence
Berkeley National Laboratory.

##  Appendix A Derivation of our improved Model

In Section II, we describe the framework for calculating the hybrid bispectrum
$B_{p_{\mathbf{\hat{n}}}p_{\mathbf{\hat{n}}}\delta}$, and summarize our
improved model and its features. Here, we outline its derivation.

###  A.1 A useful lemma

Following [46], each of the four non-zero terms in Table 1 contributing to
$B_{p_{\mathbf{\hat{n}}}p_{\mathbf{\hat{n}}}\delta}$ (=
$B_{p_{\perp}p_{\perp}\delta}$/2) can be written in index notation as:
$B_{p_{\perp}p_{\perp}\delta}(\mathbf{k_{1}},\mathbf{k_{2}},\mathbf{k_{3}})=\mathbf{\hat{p}}_{\perp;1}^{i}\mathbf{\hat{p}}_{\perp;2}^{j}B_{pp\delta}^{ij}(\mathbf{k_{1}},\mathbf{k_{2}},\mathbf{k_{3}})$,
where $\mathbf{p}_{\perp;\alpha}$ is the curl or transverse component of
$\mathbf{p}(\mathbf{k}_{\alpha})$ for $\alpha=1,2$. So, by definition,
$\mathbf{p}_{\perp;\alpha}$ must be perpendicular to the wavevector
$\mathbf{k}_{\alpha}$. As noted earlier (and in [46, 31, 32]), in the Limber
approximation, only those modes of momenta with wavevectors
$\mathbf{k}_{\alpha}$ perpendicular to the LOS direction $\mathbf{\hat{n}}$
contribute to the integrated kSZ effect. Moreover, as noted previously in [32]
under the Limber approximation, the orthogonality of
$\mathbf{p}_{\perp;\alpha}$ to the LOS gives us the effective relation:
$\mathbf{\hat{p}}_{\perp;\alpha}\approx\mathbf{\hat{k}}_{\alpha}\times\mathbf{\hat{n}}$,
upto a sign difference.

Suppose we define the LOS as the z-direction in 3D:
$\mathbf{\hat{n}}\equiv\mathbf{\hat{z}}$. For any unit wavevector
$\mathbf{\hat{m}}$, consider the dot product
$\mathbf{\hat{p}}_{\perp;\alpha}^{i}\mathbf{\hat{m}}^{i}=\mathbf{\hat{p}}_{\perp;\alpha}\cdot\mathbf{\hat{m}}$:

$\mathbf{\hat{p}}_{\perp;\alpha}\cdot\mathbf{\hat{m}}=(\mathbf{\hat{k}}_{\alpha}\times\mathbf{\hat{z}})\cdot\mathbf{\hat{m}}\hskip
42.67912pt\text{[i.e. scalar triple
product]}=\mathbf{\hat{m}}^{x}\mathbf{\hat{k}}_{\alpha}^{y}-\mathbf{\hat{m}}^{y}\mathbf{\hat{k}}_{\alpha}^{x}\approx|\mathbf{\hat{m}}\times\mathbf{\hat{k}}_{\alpha}|\hskip
56.9055pt\text{[Limber Approximation]}=\pm\sin(\theta_{\alpha})\hskip
51.21504pt\text{{\hbox{\theta_{\alpha} \equiv}} angle between
{\hbox{\mathbf{\hat{k}}_{\alpha}}} and
{\hbox{\mathbf{\hat{m}}}}}=\sqrt{1-\mu_{\alpha}^{2}}\hskip
51.21504pt{\mu_{\alpha}\equiv\mathbf{\hat{k}}_{\alpha}\cdot\mathbf{\hat{m}}},$
(26)

upto a sign difference. As a verification, we first apply this lemma to derive
the general expression for the power spectrum $P_{p_{\perp}}$, which matches
exactly with the one given in [46] (see Eq. (6) therein; note that the scaling
of the (negligible) connected four-point term also has an analogous form). We
then use this relation to calculate the geometric scalings of
$B_{p_{\perp}p_{\perp}\delta}$ below.

###  A.2 The usual term

For the leading order ‘usual’ term of $B_{p_{\perp}p_{\perp}\delta}$ given in
the first row of Table 1, we have

$\langle vv\rangle\langle\delta\delta\delta\rangle=\mathbf{\hat{p}}_{\perp
1}^{i}\mathbf{\hat{p}}_{\perp
2}^{j}\int\frac{d^{3}k}{\left(2\pi\right)^{3}}\int\frac{d^{3}k^{{}^{\prime}}}{\left(2\pi\right)^{3}}\\\
\langle
v^{i}(\mathbf{k})v^{j}(\mathbf{k^{\prime}})\rangle\langle\delta(\mathbf{k_{1}}-\mathbf{k})\delta(\mathbf{k_{2}}-\mathbf{k{{}^{\prime}}})\delta(\mathbf{k_{3}})\rangle$
(27)

We substitute the corresponding power spectrum and bispectrum in place of the
cross-correlations, similar to Eq.(12). Using the sampling property of the
Dirac delta function, we eliminate $\mathbf{k^{\prime}}$ and substitute it
with (-$\mathbf{k}$):

$\displaystyle\langle vv\rangle\langle\delta\delta\delta\rangle$
$\displaystyle=\left(2\pi\right)^{3}\delta_{D}(\mathbf{k_{1}}+\mathbf{k_{2}}+\mathbf{k_{3}})\times$
$\displaystyle\int\frac{d^{3}k}{\left(2\pi\right)^{3}}[\mathbf{\hat{p}}_{\perp
1}^{i}\mathbf{\hat{p}}_{\perp
2}^{j}\mathbf{\hat{k}}^{i}\widehat{(-\mathbf{k})}^{j}]P_{vv}(k)B_{\delta\delta\delta}(\mathbf{k_{1}}-\mathbf{k},\mathbf{k_{2}}+\mathbf{k},\mathbf{k_{3}})$
(28)

By canceling the Dirac delta enforcing closed triangles on both sides, we get

$\displaystyle
B_{p_{\perp}p_{\perp}\delta}^{\mathrm{usual}}(\mathbf{k_{1}},\mathbf{k_{2}},\mathbf{k_{3}})$
$\displaystyle=\int\frac{d^{3}k}{\left(2\pi\right)^{3}}[\mathbf{\hat{p}}_{\perp
1}^{i}\mathbf{\hat{p}}_{\perp
2}^{j}\mathbf{\hat{k}}^{i}\widehat{(-\mathbf{k})}^{j}]P_{vv}(k)B_{\delta\delta\delta}(\mathbf{k_{1}}-\mathbf{k},\mathbf{k_{2}}+\mathbf{k},\mathbf{k_{3}})$
(29)

Consider the orientation factor within square brackets in the equation above.
Let the angle between $\mathbf{\hat{k}}$ and $\mathbf{\hat{k}}_{\alpha}$ be
$\theta_{\alpha}$ for $\alpha$ = 1, 2. Then, by the lemma in Eq.(26), the
geometric scaling of the usual term is

$(\mathbf{\hat{p}}_{\perp;1}^{i}\mathbf{\hat{k}}^{i})(\mathbf{\hat{p}}_{\perp;2}^{j}\widehat{(-k)}^{j})=(\mathbf{\hat{p}}_{\perp;1}\cdot\mathbf{\hat{k}})(\mathbf{\hat{p}}_{\perp;2}\cdot\widehat{(-\mathbf{k})})=\sin(\theta_{1})\sin(\pi-\theta_{2})=\sin(\theta_{1})\sin(\theta_{2})=\sqrt{1-\mu_{1}^{2}}\sqrt{1-\mu_{2}^{2}}$
(30)

where $\mu_{\alpha}\equiv\mathbf{\hat{k}}_{\alpha}\cdot\mathbf{\hat{k}}$ for
$\alpha$= 1, 2. Thus, the usual term’s expression is

$\displaystyle
B_{p_{\perp}p_{\perp}\delta}^{\mathrm{usual}}(\mathbf{k_{1}},\mathbf{k_{2}},\mathbf{k_{3}})$
$\displaystyle=\int\frac{d^{3}k}{\left(2\pi\right)^{3}}\left[\sqrt{1-\mu_{1}^{2}}\sqrt{1-\mu_{2}^{2}}\right]P_{vv}(k)B_{\delta\delta\delta}(\mathbf{k_{1}}-\mathbf{k},\mathbf{k_{2}}+\mathbf{k},\mathbf{k_{3}})$
(31)

###  A.3 The three extra terms

Figure 5: The hybrid bispectrum $B_{p_{\perp}p_{\perp}\delta}$ at $z$ = 0.5
for equilateral triangles of side lengths $k$. (top): The 4 terms as predicted
by our theoretical model and their sum - the ‘improved’
$B_{p_{\perp}p_{\perp}\delta}$ (red), compared with the approximate D05
model’s prediction (black). (bottom): Ratios of these 4 terms and the improved
$B_{p_{\perp}p_{\perp}\delta}$ with respect to the approximate
$B_{p_{\perp}p_{\perp}\delta}$. In this equilateral case, extra-1 and extra-2
terms match.

We now calculate the three additional non-zero terms given in Table 1,
referring to each of them as the ‘extra$-i$’ terms, because we find their
geometric scalings to be $\propto 1/k_{i}$, for $i=1,2,3$. For the extra$-1$
term, we have:

$\displaystyle\langle\delta v\rangle\langle v\delta\delta\rangle$
$\displaystyle=\mathbf{\hat{p}}_{\perp 1}^{i}\mathbf{\hat{p}}_{\perp
2}^{j}\int\frac{d^{3}k}{\left(2\pi\right)^{3}}\int\frac{d^{3}k^{{}^{\prime}}}{\left(2\pi\right)^{3}}$
$\displaystyle~{}~{}~{}~{}~{}\langle\delta(\mathbf{k_{1}}-\mathbf{k})v^{j}(\mathbf{k^{\prime}})\rangle\langle
v^{i}(\mathbf{k})\delta(\mathbf{k_{2}}-\mathbf{k{{}^{\prime}}})\delta(\mathbf{k_{3}})\rangle$
(32)

We simplify the Dirac delta functions analogous to the usual term’s
calculation and eliminate the integration variable $\mathbf{k}$. On
substituting
$\widehat{(\mathbf{k}+\mathbf{k}_{1})}=(\mathbf{k}+\mathbf{k}_{1})/|\mathbf{k}+\mathbf{k}_{1}|$,
we get

$B_{p_{\perp}p_{\perp}\delta}^{\mathrm{extra-1}}(\mathbf{k_{1}},\mathbf{k_{2}},\mathbf{k_{3}})\\\
=\int\frac{d^{3}k^{\prime}}{\left(2\pi\right)^{3}}\left[\sqrt{1-\mu_{1}^{2}}\sqrt{1-\mu_{2}^{2}}\frac{k^{\prime}}{|\mathbf{k^{\prime}}+\mathbf{k}_{1}|}\right]\\\
P_{\delta
v}(k^{\prime})B_{v\delta\delta}(\mathbf{k_{1}}+\mathbf{k^{\prime}},\mathbf{k_{2}}-\mathbf{k^{\prime}},\mathbf{k_{3}}),$
(33)

where
$\mu_{\alpha}\equiv\mathbf{\hat{k}}_{\alpha}\cdot\mathbf{\hat{k^{\prime}}}$
for $\alpha$ = 1, 2.

Similarly, we eliminate the integration variable $\mathbf{k^{\prime}}$ to
obtain the extra-2 term:

$B_{p_{\perp}p_{\perp}\delta}^{\mathrm{extra-2}}(\mathbf{k_{1}},\mathbf{k_{2}},\mathbf{k_{3}})\\\
=\int\frac{d^{3}k}{\left(2\pi\right)^{3}}\left[\sqrt{1-\mu_{1}^{2}}\sqrt{1-\mu_{2}^{2}}\frac{k}{|\mathbf{k}+\mathbf{k}_{2}|}\right]\\\
P_{v\delta}(k)B_{\delta
v\delta}(\mathbf{k_{1}}-\mathbf{k},\mathbf{k_{2}}+\mathbf{k},\mathbf{k_{3}}),$
(34)

where $\mu_{\alpha}\equiv\mathbf{\hat{k}}_{\alpha}\cdot\mathbf{\hat{k}}$ for
$\alpha$ = 1, 2. As expected, the expressions for the extra$-1$ and extra$-2$
terms above are symmetric in $\mathbf{k_{1}}$ and $\mathbf{k_{2}}$.

Figure 6: Power-law indices $\alpha_{\theta}$ for each of the $\Lambda$CDM
parameters
$\boldsymbol{\theta}=\\{H_{0},\Omega_{b}h^{2},\Omega_{c}h^{2},10^{9}A_{s},n_{s}\\}$,
showing the approximate dependence of
$C_{\ell}^{\mathrm{kSZ}^{2}\times\delta_{g}}$ on them as a function of $\ell$
for CMB-S4$\times$WISE (solid lines) and CMB-S4$\times$VRO (dotted lines).

The extra$-3$ term looks different from them, since $\mathbf{k_{3}}$
corresponds to $\delta_{g}$, and not a momentum. We do an analogous
calculation, eliminating the integration variable $\mathbf{k^{\prime}}$, and
substituting
$\widehat{(\mathbf{k}+\mathbf{k}_{3})}=(\mathbf{k}+\mathbf{k}_{3})/|\mathbf{k}+\mathbf{k}_{3}|$.
Since $\mathbf{k_{3}}=-\mathbf{k_{1}}-\mathbf{k_{2}}$,

$\displaystyle
B_{p_{\perp}p_{\perp}\delta}^{\mathrm{extra-3}}(\mathbf{k_{1}},\mathbf{k_{2}},\mathbf{k_{3}})$
$\displaystyle=\int\frac{d^{3}k}{\left(2\pi\right)^{3}}\bigg{[}-\sqrt{1-\mu_{1}^{2}}\sqrt{1-\mu_{2}^{2}}\frac{k}{|\mathbf{k}+\mathbf{k}_{3}|}+$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\sqrt{1-\mu_{1}^{2}}\sqrt{1-\mu_{12}^{2}}\frac{k_{1}}{|\mathbf{k}+\mathbf{k}_{3}|}\bigg{]}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}P_{\delta\delta}(|k-k_{1}|)B_{vv\delta}(\mathbf{k},\mathbf{k_{1}}+\mathbf{k_{2}}-\mathbf{k},\mathbf{k_{3}}),$
(35)

where $\mu_{\alpha}\equiv\mathbf{\hat{k}}_{\alpha}\cdot\mathbf{\hat{k}}$ for
$\alpha$ = 1, 2, and
$\mu_{12}\equiv\mathbf{\hat{k}}_{1}\cdot\mathbf{\hat{k}}_{2}$. Note that this
extra$-3$ term is generally negative at large scales and positive at small
scales since it has two parts. A similar expression is alternatively obtained
if we choose to eliminate the integration variable $\mathbf{k}$, where the
roles of $k_{1}$ and $k_{2}$ get exchanged.

Thus, the ‘improved’ hybrid bispectrum
$\left(B_{p_{\mathbf{\hat{n}}}p_{\mathbf{\hat{n}}}\delta}=\frac{1}{2}B_{p_{\perp}p_{\perp}\delta}\right)$
predicted by our derived improved model is:

$B_{p_{\mathbf{\hat{n}}}p_{\mathbf{\hat{n}}}\delta}=(1/2)\left[B_{p_{\perp}p_{\perp}\delta}^{\mathrm{usual}}+B_{p_{\perp}p_{\perp}\delta}^{\mathrm{extra-1}}+B_{p_{\perp}p_{\perp}\delta}^{\mathrm{extra-2}}+B_{p_{\perp}p_{\perp}\delta}^{\mathrm{extra-3}}\right].$
(36)

In Figure 5 (top panel), we show the individual contributions of the four
terms of $B_{p_{\perp}p_{\perp}\delta}$ (plotting the absolute value of
extra$-3$) and their sum - the improved model’s prediction for
$B_{p_{\perp}p_{\perp}\delta}$ (red), for the special case of equilateral
triangles at $z=0.5$. We also show the approximate D05 model’s prediction
(black). The bottom panel shows the ratios of these four terms and of their
sum with respect to the approximate prediction. For this particular case, the
models match reasonably at the smallest scales but differ at larger scales.
When the total estimator $C_{\ell}^{\mathrm{kSZ}^{2}\times\delta_{g}}$ is
computed, the CMB filters applied suppress $B_{p_{\perp}p_{\perp}\delta}$ at
large scales (small $k$), and $B_{p_{\perp}p_{\perp}\delta}$ contributions
from all triangle shapes are summed together.

##  Appendix B Understanding the Cosmological Dependence

In Section V.1, we show that the marginalized uncertainty on
$A_{\mathrm{kSZ}^{2}}$ goes from a sub-percent level to $\sim 7\%$ when
cosmology is allowed to vary. To understand this result roughly, we
approximately model the dependence of
$C_{\ell}^{\mathrm{kSZ}^{2}\times\delta_{g}}$ (denoted here as ‘$C_{\ell}$’)
on a parameter $\theta$ as a power-law scaling:
$C_{\ell}\propto\theta^{\alpha_{\theta}}$. We estimate the index
$\alpha_{\theta}$ for each of the $\Lambda$CDM parameters as:
$\alpha_{\theta}\approx\frac{(\Delta
C_{\ell}/C_{\ell})}{(\Delta\theta/\theta)}$, where $\Delta\theta$ are the same
step-sizes used in Section V.1 for computing partial derivatives, and $\Delta
C_{\ell}$ are the corresponding changes in the signal.

From Eq.(25), since $P_{\delta\delta}\propto A_{s}$, we expect that $C_{\ell}$
roughly scales as $A_{s}^{3}$, as seen in Figure 6 across scales.
$P_{\delta\delta}$ depends non-linearly on the total matter density today,
$\Omega_{m}$, while the linear growth rate $f$ roughly scales as
$\Omega_{m}^{0.55}$. Since $\Omega_{m}=\Omega_{c}+\Omega_{b}$, and the
fiducial value of $\Omega_{c}$ is about 5 times $\Omega_{b}$, $C_{\ell}$ has
the strongest dependence on $\Omega_{c}h^{2}$ with an $\alpha_{\theta}$ of
$\sim 4.5$-5, while $\alpha_{\Omega_{b}h^{2}}\approx 1$. Because $n_{s}$
itself is the slope of the primordial spectrum, its $\alpha_{\theta}$ varies
strongly with $\ell$ and $C_{\ell}$’s dependence on it is not well-modeled as
a power-law.

Given rough power-law scalings, each $\theta$ with a residual error of
$\sigma(\theta)$ contributes $\sim(|\alpha_{\theta}|\sigma(\theta)/\theta)$ to
the relative propagated error on $C_{\ell}$. Thus, if they were fully
uncorrelated, $A_{s}$ and $\Omega_{c}h^{2}$ together would imply a propagated
error of $\sim 9\%$ on $A_{\mathrm{kSZ}^{2}}$. While the actual uncertainty on
$A_{\mathrm{kSZ}^{2}}$ is lower due to degeneracies between $\Lambda$CDM
parameters, a majority of it arises from residual uncertainties on $A_{s}$ and
$\Omega_{c}h^{2}$ in the Planck prior. The dependence on $\Lambda$CDM
parameters of the $C_{\ell}^{\mathrm{kSZ}^{2}\times\delta_{g}}$ signal from
SO$\times$WISE is similar to that from CMB-S4$\times$WISE, although it differs
for the signal from CMB-S4$\times$VRO due to the $z$-dependent effects of
cosmology. We discuss this further in Section V.2.

## References

  * The Simons Observatory Collaboration [2019] The Simons Observatory Collaboration, J. Cosmology Astropart. Phys 2019, 10.1088/1475-7516/2019/02/056 (2019), arXiv:1808.07445 [astro-ph.CO] .
  * The CMB-S4 Collaboration  _et al._ [2019] The CMB-S4 Collaboration, K. Abazajian, G. Addison, P. Adshead, Z. Ahmed, S. W. Allen, D. Alonso, M. Alvarez, A. Anderson, K. S. Arnold, C. Baccigalupi, K. Bailey, D. Barkats, D. Barron, P. S. Barry, _et al._ , arXiv e-prints , arXiv:1907.04473 (2019), arXiv:1907.04473 [astro-ph.IM] .
  * Sunyaev and Zeldovich [1972] R. A. Sunyaev and Y. B. Zeldovich, Comments on Astrophysics and Space Physics 4 (1972).
  * Sunyaev and Zeldovich [1980] R. A. Sunyaev and I. B. Zeldovich, ARA&A 18, 537 (1980).
  * Ostriker and Vishniac [1986] J. P. Ostriker and E. T. Vishniac, ApJ 306, L51 (1986).
  * Ferraro  _et al._ [2016] S. Ferraro, J. C. Hill, N. Battaglia, J. Liu, and D. N. Spergel, Phys. Rev. D 94, 123526 (2016), arXiv:1605.02722 [astro-ph.CO] .
  * Battaglia  _et al._ [2017] N. Battaglia, S. Ferraro, E. Schaan, and D. N. Spergel, J. Cosmology Astropart. Phys 2017, 10.1088/1475-7516/2017/11/040 (2017), arXiv:1705.05881 [astro-ph.CO] .
  * Battaglia  _et al._ [2019] N. Battaglia, J. C. Hill, S. Amodeo, J. G. Bartlett, K. Basu, J. Erler, S. Ferraro, L. Hernquist, M. Madhavacheril, M. McQuinn, T. Mroczkowski, D. Nagai, E. Schaan, R. Somerville, R. Sunyaev, M. Vogelsberger, and J. Werk, BAAS 51, 10.48550/arXiv.1903.04647 (2019), arXiv:1903.04647 [astro-ph.CO] .
  * Amodeo  _et al._ [2021] S. Amodeo, N. Battaglia, E. Schaan, S. Ferraro, E. Moser, and ACT Collaboration, Phys. Rev. D 103, 10.1103/PhysRevD.103.063514 (2021), arXiv:2009.05558 [astro-ph.CO] .
  * La Plante  _et al._ [2022] P. La Plante, J. Sipple, and A. Lidz, ApJ 928, 10.3847/1538-4357/ac5752 (2022), arXiv:2111.13717 [astro-ph.CO] .
  * Bhattacharya and Kosowsky [2008] S. Bhattacharya and A. Kosowsky, Phys. Rev. D 77, 10.1103/PhysRevD.77.083004 (2008), arXiv:0712.0034 [astro-ph] .
  * Mueller  _et al._ [2015a] E.-M. Mueller, F. de Bernardis, R. Bean, and M. D. Niemack, ApJ 808, 10.1088/0004-637X/808/1/47 (2015a), arXiv:1408.6248 [astro-ph.CO] .
  * Mueller  _et al._ [2015b] E.-M. Mueller, F. de Bernardis, R. Bean, and M. D. Niemack, Phys. Rev. D 92, 10.1103/PhysRevD.92.063501 (2015b), arXiv:1412.0592 [astro-ph.CO] .
  * Roncarelli  _et al._ [2018] M. Roncarelli, M. Baldi, and F. Villaescusa-Navarro, MNRAS 481, 2497 (2018), arXiv:1805.11607 [astro-ph.CO] .
  * Roncarelli  _et al._ [2017] M. Roncarelli, F. Villaescusa-Navarro, and M. Baldi, MNRAS 467, 985 (2017), arXiv:1702.00676 [astro-ph.CO] .
  * Hand  _et al._ [2012] N. Hand, G. E. Addison, E. Aubourg, N. Battaglia, E. S. Battistelli, D. Bizyaev, J. R. Bond, H. Brewington, J. Brinkmann, B. R. Brown, S. Das, K. S. Dawson, M. J. Devlin, J. Dunkley, R. Dunner, _et al._ , Phys. Rev. Lett. 109, 041101 (2012), arXiv:1203.4219 [astro-ph.CO] .
  * Ferreira  _et al._ [1999] P. G. Ferreira, R. Juszkiewicz, H. A. Feldman, M. Davis, and A. H. Jaffe, ApJ 515, L1 (1999), arXiv:astro-ph/9812456 [astro-ph] .
  * Planck Collaboration [2016] Planck Collaboration, A&A 586, A140 (2016), arXiv:1504.03339 [astro-ph.CO] .
  * Soergel  _et al._ [2016] B. Soergel, S. Flender, K. T. Story, and The DES and SPT Collaborations, Monthly Notices of the Royal Astronomical Society 461, 3172 (2016).
  * Calafut  _et al._ [2021] V. Calafut, P. A. Gallardo, E. M. Vavagiakis, S. Amodeo, S. Aiola, J. E. Austermann, N. Battaglia, E. S. Battistelli, J. A. Beall, R. Bean, J. R. Bond, E. Calabrese, S. K. Choi, N. F. Cothard, M. J. Devlin, _et al._ , Phys. Rev. D 104, 043502 (2021).
  * Ho  _et al._ [2009] S. Ho, S. Dedeo, and D. Spergel, arXiv e-prints , arXiv:0903.2845 (2009).
  * Shao  _et al._ [2011] J. Shao, P. Zhang, W. Lin, Y. Jing, and J. Pan, MNRAS 413, 628 (2011).
  * Hernández-Monteagudo  _et al._ [2015] C. Hernández-Monteagudo, Y.-Z. Ma, F. S. Kitaura, W. Wang, R. Génova-Santos, J. Macías-Pérez, and D. Herranz, Phys. Rev. Lett. 115, 191301 (2015).
  * Schaan  _et al._ [2016] E. Schaan, S. Ferraro, M. Vargas-Magaña, K. M. Smith, S. Ho, S. Aiola, N. Battaglia, J. R. Bond, F. De Bernardis, E. Calabrese, H.-M. Cho, M. J. Devlin, J. Dunkley, P. A. Gallardo, M. Hasselfield, and ACTPol Collaboration, Phys. Rev. D 93, 10.1103/PhysRevD.93.082002 (2016), arXiv:1510.06442 [astro-ph.CO] .
  * Schaan  _et al._ [2021] E. Schaan, S. Ferraro, S. Amodeo, N. Battaglia, S. Aiola, _et al._ (Atacama Cosmology Telescope Collaboration), Phys. Rev. D 103, 10.1103/PhysRevD.103.063513 (2021).
  * Smith  _et al._ [2018] K. M. Smith, M. S. Madhavacheril, M. Münchmeyer, S. Ferraro, U. Giri, and M. C. Johnson, arXiv e-prints 10.48550/arXiv.1810.13423 (2018).
  * Münchmeyer  _et al._ [2019] M. Münchmeyer, M. S. Madhavacheril, S. Ferraro, M. C. Johnson, and K. M. Smith, Phys. Rev. D 100, 10.1103/PhysRevD.100.083508 (2019).
  * Giri and Smith [2022] U. Giri and K. M. Smith, J. Cosmology Astropart. Phys 2022, 028 (2022), arXiv:2010.07193 [astro-ph.CO] .
  * Keisler and Schmidt [2013] R. Keisler and F. Schmidt, ApJ 765, 10.1088/2041-8205/765/2/L32 (2013).
  * Flender  _et al._ [2016] S. Flender, L. Bleem, H. Finkel, S. Habib, K. Heitmann, and G. Holder, ApJ 823, 10.3847/0004-637X/823/2/98 (2016), arXiv:1511.02843 [astro-ph.CO] .
  * Doré  _et al._ [2004] O. Doré, J. F. Hennawi, and D. N. Spergel, ApJ 606, 46 (2004), arXiv:astro-ph/0309337 [astro-ph] .
  * DeDeo  _et al._ [2005] S. DeDeo, D. N. Spergel, and H. Trac, arXiv e-prints , astro-ph/0511060 (2005), arXiv:astro-ph/0511060 [astro-ph] .
  * Wright  _et al._ [2010] E. L. Wright, P. R. M. Eisenhardt, A. K. Mainzer, M. E. Ressler, R. M. Cutri, T. Jarrett, J. D. Kirkpatrick, D. Padgett, R. S. McMillan, M. Skrutskie, S. A. Stanford, M. Cohen, R. G. Walker, J. C. Mather, D. Leisawitz, _et al._ , AJ 140, 1868 (2010).
  * The LSST collaboration [2009] The LSST collaboration, arXiv e-prints , arXiv:0912.0201 (2009), arXiv:0912.0201 [astro-ph.IM] .
  * Hill  _et al._ [2016] J. C. Hill, S. Ferraro, N. Battaglia, J. Liu, and D. N. Spergel, Phys. Rev. Lett. 117, 10.1103/PhysRevLett.117.051301 (2016), arXiv:1603.01608 [astro-ph.CO] .
  * Kusiak  _et al._ [2021] A. Kusiak, B. Bolliet, S. Ferraro, J. C. Hill, and A. Krolewski, Phys. Rev. D 104, 10.1103/PhysRevD.104.043518 (2021), arXiv:2102.01068 [astro-ph.CO] .
  * Bolliet  _et al._ [2023a] B. Bolliet, J. Colin Hill, S. Ferraro, A. Kusiak, and A. Krolewski, J. Cosmology Astropart. Phys 2023, 10.1088/1475-7516/2023/03/039 (2023a), arXiv:2208.07847 [astro-ph.CO] .
  * Planck Collaboration [2020] Planck Collaboration, A&A 641, 10.1051/0004-6361/201833910 (2020), arXiv:1807.06209 [astro-ph.CO] .
  * Lesgourgues and Pastor [2006] J. Lesgourgues and S. Pastor, Phys. Rep. 429, 307 (2006), arXiv:astro-ph/0603494 [astro-ph] .
  * Dvorkin  _et al._ [2019] C. Dvorkin, M. Gerbino, D. Alonso, N. Battaglia, S. Bird, A. Diaz Rivero, A. Font-Ribera, G. Fuller, M. Lattanzi, M. Loverde, J. B. Muñoz, B. Sherwin, A. Slosar, and F. Villaescusa-Navarro, BAAS 51, 10.48550/arXiv.1903.03689 (2019), arXiv:1903.03689 [astro-ph.CO] .
  * Ma  _et al._ [2018] Q. Ma, K. Helgason, E. Komatsu, B. Ciardi, and A. Ferrara, MNRAS 476, 4025 (2018), arXiv:1712.05305 [astro-ph.CO] .
  * Kaiser [1992] N. Kaiser, ApJ 388, 272 (1992).
  * Limber [1953] D. N. Limber, ApJ 117, 134 (1953).
  * Hu and White [1996] W. Hu and M. White, A&A 315, 33 (1996), arXiv:astro-ph/9507060 [astro-ph] .
  * Zhang  _et al._ [2004] P. Zhang, U.-L. Pen, and H. Trac, MNRAS 347, 1224 (2004), arXiv:astro-ph/0304534 [astro-ph] .
  * Ma and Fry [2002] C. P. Ma and J. N. Fry, Physical Review Letters 88, 2113011 (2002).
  * Hu [2000] W. Hu, The Astrophysical Journal 529, 12 (2000).
  * Bennett  _et al._ [2013] C. L. Bennett, D. Larson, J. L. Weiland, N. Jarosik, G. Hinshaw, N. Odegard, K. M. Smith, R. S. Hill, B. Gold, M. Halpern, E. Komatsu, M. R. Nolta, L. Page, D. N. Spergel, E. Wollack, _et al._ , ApJS 208, 10.1088/0067-0049/208/2/20 (2013), arXiv:1212.5225 [astro-ph.CO] .
  * Bobin, J.  _et al._ [2016] Bobin, J., Sureau, F., and Starck, J.-L., A&A 591, A50 (2016).
  * Bobin, J.  _et al._ [2013] Bobin, J., Starck, J.-L., Sureau, F., and Basak, S., A&A 550, A73 (2013).
  * Battaglia  _et al._ [2010] N. Battaglia, J. R. Bond, C. Pfrommer, J. L. Sievers, and D. Sijacki, ApJ 725, 91 (2010), arXiv:1003.4256 [astro-ph.CO] .
  * Lewis  _et al._ [2000] A. Lewis, A. Challinor, and A. Lasenby, ApJ 538, 473 (2000), arXiv:astro-ph/9911177 [astro-ph] .
  * Howlett  _et al._ [2012] C. Howlett, A. Lewis, A. Hall, and A. Challinor, J. Cosmology Astropart. Phys 2012, 10.1088/1475-7516/2012/04/027 (2012), arXiv:1201.3654 [astro-ph.CO] .
  * Tegmark  _et al._ [2003] M. Tegmark, A. de Oliveira-Costa, and A. J. S. Hamilton, Phys. Rev. D 68, 10.1103/PhysRevD.68.123523 (2003).
  * Eriksen  _et al._ [2004] H. K. Eriksen, A. J. Banday, K. M. Górski, and P. B. Lilje, ApJ 612, 633 (2004), arXiv:astro-ph/0403098 [astro-ph] .
  * Ferraro  _et al._ [2022] S. Ferraro, E. Schaan, and E. Pierpaoli, arXiv e-prints , arXiv:2205.10332 (2022).
  * Ferraro  _et al._ [2015] S. Ferraro, B. D. Sherwin, and D. N. Spergel, Phys. Rev. D 91, 10.1103/PhysRevD.91.083533 (2015), arXiv:1401.1193 [astro-ph.CO] .
  * Jarrett  _et al._ [2011] T. H. Jarrett, M. Cohen, F. Masci, E. Wright, D. Stern, D. Benford, A. Blain, S. Carey, R. M. Cutri, P. Eisenhardt, C. Lonsdale, A. Mainzer, K. Marsh, D. Padgett, S. Petty, _et al._ , ApJ 735, 112 (2011).
  * Yan  _et al._ [2013] L. Yan, E. Donoso, C.-W. Tsai, D. Stern, R. J. Assef, P. Eisenhardt, A. W. Blain, R. Cutri, T. Jarrett, S. A. Stanford, E. Wright, C. Bridge, and D. A. Riechers, AJ 145, 55 (2013), arXiv:1209.2065 [astro-ph.CO] .
  * Ivezić  _et al._ [2019] Ž. Ivezić, S. M. Kahn, J. A. Tyson, B. Abel, E. Acosta, R. Allsman, D. Alonso, Y. AlSayyad, S. F. Anderson, J. Andrew, and otheers, ApJ 873, 10.3847/1538-4357/ab042c (2019), arXiv:0805.2366 .
  * Weinberg  _et al._ [2004] D. H. Weinberg, R. Davé, N. Katz, and L. Hernquist, ApJ 601, 1 (2004), arXiv:astro-ph/0212356 [astro-ph] .
  * LoVerde and Afshordi [2008] M. LoVerde and N. Afshordi, Physical Review D 78, 10.1103/physrevd.78.123506 (2008).
  * Takahashi  _et al._ [2012] R. Takahashi, M. Sato, T. Nishimichi, A. Taruya, and M. Oguri, ApJ 761, 10.1088/0004-637X/761/2/152 (2012), arXiv:1208.2701 [astro-ph.CO] .
  * Fukuda  _et al._ [1998] Y. Fukuda, T. Hayakawa, E. Ichihara, K. Inoue, K. Ishihara, H. Ishino, Y. Itow, _et al._ (Super-Kamiokande Collaboration), Phys. Rev. Lett. 81, 1562 (1998).
  * SNO Collaboration [2002] SNO Collaboration, Phys. Rev. Lett. 89, 10.1103/PhysRevLett.89.011301 (2002).
  * Abazajian  _et al._ [2015] K. N. Abazajian, K. Arnold, J. Austermann, B. A. Benson, C. Bischoff, J. Bock, J. R. Bond, J. Borrill, E. Calabrese, J. E. Carlstrom, C. S. Carvalho, C. L. Chang, H. C. Chiang, S. Church, A. Cooray, _et al._ , Astroparticle Physics 63, 66 (2015), arXiv:1309.5383 [astro-ph.CO] .
  * Kiakotou  _et al._ [2008] A. Kiakotou, Ø. Elgarøy, and O. Lahav, Phys. Rev. D 77, 10.1103/PhysRevD.77.063005 (2008).
  * Gil-Marín  _et al._ [2012] H. Gil-Marín, C. Wagner, F. Fragkoudi, R. Jimenez, and L. Verde, J. Cosmology Astropart. Phys 2012, 10.1088/1475-7516/2012/02/047 (2012), arXiv:1111.4477 [astro-ph.CO] .
  * Levi and Vlah [2016] M. Levi and Z. Vlah, arXiv e-prints 10.48550/arXiv.1605.09417 (2016).
  * Ruggeri  _et al._ [2018] R. Ruggeri, E. Castorina, C. Carbone, and E. Sefusatti, J. Cosmology Astropart. Phys 2018, 10.1088/1475-7516/2018/03/003 (2018), arXiv:1712.02334 [astro-ph.CO] .
  * Coulton  _et al._ [2019] W. R. Coulton, J. Liu, M. S. Madhavacheril, V. Böhm, and D. N. Spergel, J. Cosmology Astropart. Phys 2019, 043 (2019), arXiv:1810.02374 [astro-ph.CO] .
  * Mallaby-Kay  _et al._ [2023] M. Mallaby-Kay, S. Amodeo, J. C. Hill, M. Aguena, S. Allam, O. Alves, J. Annis, N. Battaglia, E. S. Battistelli, E. J. Baxter, K. Bechtol, M. R. Becker, E. Bertin, J. R. Bond, D. Brooks, _et al._ , arXiv e-prints 10.48550/arXiv.2305.06792 (2023).
  * Seljak [2000] U. Seljak, MNRAS 318, 203 (2000), arXiv:astro-ph/0001493 [astro-ph] .
  * Cromer  _et al._ [2022] D. Cromer, N. Battaglia, H. Miyatake, and M. Simet, J. Cosmology Astropart. Phys 2022, 10.1088/1475-7516/2022/10/034 (2022), arXiv:2104.06925 [astro-ph.CO] .
  * Lewis and Challinor [2006] A. Lewis and A. Challinor, Phys. Rep. 429, 1 (2006), arXiv:astro-ph/0601594 [astro-ph] .
  * Sherwin  _et al._ [2012] B. D. Sherwin, S. Das, A. Hajian, G. Addison, J. R. Bond, D. Crichton, M. J. Devlin, J. Dunkley, M. B. Gralla, M. Halpern, J. C. Hill, A. D. Hincks, J. P. Hughes, K. Huffenberger, R. Hlozek, _et al._ , Phys. Rev. D 86, 10.1103/PhysRevD.86.083006 (2012).
  * Tegmark  _et al._ [1997] M. Tegmark, A. N. Taylor, and A. F. Heavens, ApJ 480, 22 (1997), arXiv:astro-ph/9603021 [astro-ph] .
  * Krolewski  _et al._ [2020] A. Krolewski, S. Ferraro, E. F. Schlafly, and M. White, J. Cosmology Astropart. Phys 2020, 10.1088/1475-7516/2020/05/047 (2020), arXiv:1909.07412 [astro-ph.CO] .
  * Sailer  _et al._ [2021] N. Sailer, E. Castorina, S. Ferraro, and M. White, J. Cosmology Astropart. Phys 2021, 10.1088/1475-7516/2021/12/049 (2021), arXiv:2106.09713 [astro-ph.CO] .
  * Madhavacheril  _et al._ [2017] M. S. Madhavacheril, N. Battaglia, and H. Miyatake, Phys. Rev. D 96, 10.1103/PhysRevD.96.103525 (2017), arXiv:1708.07502 [astro-ph.CO] .
  * Cromer  _et al._ [2019] D. Cromer, N. Battaglia, and M. S. Madhavacheril, Phys. Rev. D 100, 063529 (2019), arXiv:1903.00976 [astro-ph.CO] .
  * LiteBIRD Collaboration [2020] LiteBIRD Collaboration, in  _Space Telescopes and Instrumentation 2020: Optical, Infrared, and Millimeter Wave_, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 11443, edited by M. Lystrup and M. D. Perrin (2020) arXiv:2101.12449 [astro-ph.IM] .
  * Muñoz and Dvorkin [2018] J. B. Muñoz and C. Dvorkin, Physical Review D 98, 10.1103/physrevd.98.043503 (2018).
  * Chiang  _et al._ [2018] C.-T. Chiang, W. Hu, Y. Li, and M. LoVerde, Phys. Rev. D 97, 10.1103/PhysRevD.97.123526 (2018), arXiv:1710.01310 [astro-ph.CO] .
  * Mishra-Sharma  _et al._ [2018] S. Mishra-Sharma, D. Alonso, and J. Dunkley, Phys. Rev. D 97, 10.1103/PhysRevD.97.123544 (2018), arXiv:1803.07561 [astro-ph.CO] .
  * DESI Collaboration [2016] DESI Collaboration, arXiv e-prints , arXiv:1611.00036 (2016), 1611.00036 [astro-ph.IM] .
  * Addison  _et al._ [2013] G. E. Addison, J. Dunkley, and J. R. Bond, Monthly Notices of the Royal Astronomical Society 436, 1896 (2013).
  * Battaglia [2016] N. Battaglia, J. Cosmology Astropart. Phys 2016, 058 (2016), arXiv:1607.02442 [astro-ph.CO] .
  * Spurio Mancini  _et al._ [2022] A. Spurio Mancini, D. Piras, J. Alsing, B. Joachimi, and M. P. Hobson, MNRAS 511, 1771 (2022), arXiv:2106.03846 [astro-ph.CO] .
  * Bolliet  _et al._ [2023b] B. Bolliet, A. Spurio Mancini, J. C. Hill, M. Madhavacheril, H. T. Jense, E. Calabrese, and J. Dunkley, arXiv e-prints 10.48550/arXiv.2303.01591 (2023b).