# A new “temperature inversion” estimator to detect CMB patchy screening by
large-scale structure

Theo Schutt  0000-0002-7187-9628<EMAIL_ADDRESS>Kavli Institute for
Particle Astrophysics and Cosmology, 382 Via Pueblo Mall, Stanford, CA 94305,
USA  SLAC National Accelerator Laboratory, 2575 Sand Hill Road, Menlo Park, CA
94025, USA  Department of Physics, Stanford University, Stanford, CA 94305,
USA  Abhishek S. Maniyar  0000-0002-4617-9320 Kavli Institute for Particle
Astrophysics and Cosmology, 382 Via Pueblo Mall, Stanford, CA 94305, USA  SLAC
National Accelerator Laboratory, 2575 Sand Hill Road, Menlo Park, CA 94025,
USA  Department of Physics, Stanford University, Stanford, CA 94305, USA
Emmanuel Schaan  0000-0002-4619-8927 Kavli Institute for Particle Astrophysics
and Cosmology, 382 Via Pueblo Mall, Stanford, CA 94305, USA  SLAC National
Accelerator Laboratory, 2575 Sand Hill Road, Menlo Park, CA 94025, USA
William R. Coulton  0000-0002-1297-3673 Center for Computational Astrophysics,
Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA  Kavli Institute
for Cosmology Cambridge, Madingley Road, Cambridge CB3 0HA, UK  Nishant Mishra
0000-0002-9141-9792 Department of Astronomy, University of Michigan, Ann
Arbor, MI 48109, USA  Department of Physics, University of California,
Berkeley, CA 94720, USA

###### Abstract

Thomson scattering of cosmic microwave background (CMB) photons imprints
various properties of the baryons around galaxies on the CMB. One such
imprint, called patchy screening, is a direct probe of the gas density profile
around galaxies. It usefully complements the information from the kinematic
and thermal Sunyaev-Zel’dovich effects and does not require individual
redshifts. In this paper, we derive new estimators of patchy screening called
the “temperature inversion” (TI) and “signed” estimators, analogous to the
gradient inversion estimator of CMB lensing. Pedagogically, we clarify the
relation between these estimators and the standard patchy screening quadratic
estimator (QE). The new estimators trade optimality for robustness to biases
caused by the dominant CMB lensing and foreground contaminants, allowing the
use of smaller angular scales. We perform a simulated analysis to
realistically forecast the expected precision of patchy screening measurements
from four CMB experiments, ACT, SPT, Simons Observatory (SO) and CMB-S4,
cross-correlated with three galaxy samples from BOSS, unWISE and the simulated
Rubin LSST Data Challenge 2 catalog. Our results give further confidence in
the first detection of this effect from the ACT$\times$unWISE data in the
companion paper and show patchy screening will be a powerful observable for
future surveys like SO, CMB-S4 and LSST. Implementations of the patchy
screening QE and the TI and signed estimators are publicly available in our
LensQuEst111https://github.com/EmmanuelSchaan/LensQuEst and
ThumbStack222https://github.com/EmmanuelSchaan/ThumbStack software packages,
respectively.

##  I Introduction

While the cosmic microwave background (CMB) is an extremely powerful window
onto the initial conditions of our Universe, it also contains a wealth of
information on its late-time large-scale structure (LSS). Indeed, large-scale
structure properties including gas density, thermal and bulk velocities,
gravitational potential and more leave shadows on the CMB called secondary
anisotropies; see Aghanim2008Rev  for a review. One important class of
secondary anisotropies is due to Thomson scattering of the CMB photons by free
electrons in the ionized plasma (gas) around galaxies and clusters. This class
includes the thermal and kinematic Sunyaev-Zel’dovich (tSZ & kSZ) effects ZS69
; SZ72 ; SZ80 , which are measured with high precision and have already
yielded crucial insights into cosmology Salvati_2018 ; Bolliet_2018 ;
Soergel:2017ahb ; Amon_2022 ; Arico_2023.678A.109A  and the thermodynamics of
gaseous halos around galaxies and clusters (e.g., Schaan:2015uaa, ;
Schaan2021a, ; Amodeo2021, ; Schiappucci_2023, ).

Another effect in this class is the patchy screening effect (Ostriker_1986, ;
Vishniac:1987wm, ; DvorkinSmith2009, ; Feng_2018, ). This effect is perhaps
the simplest, most intuitive one: a free electron along a given line of sight
(LOS) scatters CMB photons out of that LOS, and replaces them with CMB photons
from the average blackbody distribution, imprinting a shadow. However, patchy
screening is smaller in amplitude than the more complex tSZ and kSZ effects
since they scale as

$\left\\{\begin{aligned} &\delta T^{\rm
tSZ}\propto\tau(v_{\text{thermal}}/c)^{2},&\left(v_{\text{thermal}}/c\right)^{2}&\sim
10^{-2},\\\ &\delta T^{\rm
kSZ}\propto\tau(v_{\text{bulk}}/c),&v_{\text{bulk}}/c&\sim 10^{-3},\\\ &\delta
T^{\rm screening}\propto\tau(\delta T_{0}/T_{0}),&\delta T_{0}/T_{0}&\sim
10^{-4},\\\ \end{aligned}\right.$ (1)

where $\tau\sim 10^{-3}$ for galaxy group-sized halos.

These scalings also highlight the complementarity of these different effects.
Like kSZ, patchy screening is independent of the gas temperature. It thus
scales linearly with halo mass, unlike tSZ, and is thus ideally suited for
constraining baryonic effects in galaxy lensing by measuring baryon density
profiles around the group-sized halos that dominate the matter power spectrum
Amodeo2021 .

Like tSZ, patchy screening does not require any individual redshift
information to be detected. While the kSZ effect can also be detected without
redshift information Dore:2003ex ; Hill:2016dta ; Kusiak_2021 , this is done
at the cost of squaring the kSZ, which destroys its linear scaling with gas
density and halo mass. Combining screening, kSZ and tSZ allows us to form a
complete picture of the gas thermodynamics around halos (gas density,
temperature). Screening measurements also enables breaking of the kSZ
degeneracy between galaxy gas content and galaxy velocity. Thus, by combining
them, one may extract the amplitude of large-scale peculiar velocities in the
Universe to test the neutrino masses Mueller:2014dba , the growth rate of
structure Sugiyama_2017.01.057S , general relativity Mueller:2014nsa , and
measure primordial non-Gaussianity 2019PhRvD.100h3508M . Finally, while
screening and kSZ also probe the epoch of reionization Hu:1999vq ;
DvorkinSmith2009 ; Natarajan2013.776.82N ; Roy2018.05.014R , we will not
consider this in this paper, focusing instead on the late-time signal from
galaxies and clusters.

Observationally, the patchy screening effect is unique in that it changes sign
depending on the sign of the local large-scale CMB temperature. This is unlike
extragalactic foregrounds like the cosmic infrared background (CIB) or the tSZ
effect, whose sign are fixed at any given frequency. This crucial feature is
analogous to kSZ, whose signal changes sign based on the sign of the local LOS
peculiar velocity. Like for kSZ, it allows us to cleanly distinguish screening
from these foregrounds. This property also differentiates screening from CMB
lensing, whose sign changes based on the local large-scale temperature
gradient, rather than temperature value.

Leveraging this crucial feature of patchy screening, we derive new estimators
dubbed “temperature inversion” (TI) and “signed” estimators, analogous to the
gradient inversion (GI) estimator of CMB lensing Horowitz2019.485.3919H ;
Hadzhiyska2019.100b3547H . By making use of the sign change of the screening
effect, these estimators are naturally robust to contamination from CIB, tSZ
and kSZ, as well as CMB lensing.

We carefully present the relation between these new estimators and the
standard quadratic estimator (QE) for patchy screening Hu:2001kj ;
DvorkinSmith2009 ; Su2011.1106.4313S . We show how the new estimators do not
outperform the standard QE in terms of raw statistical power when using same
multipole range. However, in its simplest form, the QE is biased by
foregrounds and CMB lensing. While CMB lensing contamination to the QE can be
removed via “bias hardening,” the situation may be more complex for
extragalactic foregrounds.

While we focus on detecting screening from the CMB temperature only, these
estimators generalize trivially to polarization by simply replacing
temperature with Stokes Q and U.

The paper is organized as follows. In Section II, we review the physical
origin of patchy screening. Section III covers the standard QE for patchy
screening and its application to reconstructing both large- and small-scale
optical depth $\tau$. Section IV introduces the new TI and signed estimators
and clarifies their relationship to the small-scale QE. In Section V, we study
the biases to the QE due to lensing and extragalactic foregrounds and explain
the TI estimator’s robustness to these biases. Section VI provides an overview
of our realistic forecast simulations and analysis and presents the resulting
patchy screening detection signal-to-noise ratios (SNRs) for TI and QE for a
set of current and future CMB and LSS surveys. We conclude with Section VII.

##  II Patchy screening of CMB temperature and polarization

Consider a line of sight where the CMB temperature $T^{0}$ deviates from the
sky averaged CMB temperature $\bar{T}^{0}$, by an amount $\delta
T^{0}=T^{0}-\bar{T}^{0}$. If an electron cloud with optical depth $\tau$ to
Thomson scattering is present along that line of sight, two things occur.
First, a fraction $1-e^{-\tau}$ of the CMB photons with temperature $T^{0}$ is
deflected out of the line of sight, leaving only a fraction $e^{-\tau}$ of the
incident photons at temperature $T^{0}$. Second, an equal fraction
$1-e^{-\tau}$ of CMB photons from all other directions is scattered into the
line of sight. Because this effect averages over all directions, as seen by
the electron cloud, these photons have temperature $\bar{T}^{0}$. Thus the
observed temperature $T$ along the line of sight, after Thomson screening, is:

$T=T^{0}e^{-\tau}+\bar{T}^{0}\left(1-e^{-\tau}\right).$ (2)

Considering temperature deviations with respect to the sky average CMB
temperature $\bar{T}^{0}$, we thus obtain:

$\delta T(\bm{x})=\delta T^{0}(\bm{x})e^{-\tau(\bm{x})},$ (3)

where $\bm{x}$ is the angular position on the sky. The screening effect
operates in the same way on temperature T, and Q and U polarizations. For ease
of notation, we will thus focus on temperature here, but note that all our
results apply to polarization as well.

In polarization, an additional “scattering” effect adds to this screening
effect, where the local CMB temperature quadrupole observed by the free
electrons is converted to linear polarization. This “scattering” occurs most
noticeably at the last scattering surface and around the epoch of
reionization. It also occurs at later times due to the ionized gas around
galaxies and clusters, an efect which has been called polarized Sunyaev-
Ze’dovich (pSZ) effect SZ80Rev ; Sazonov1999 ; Emritte2016 .

Going further, we decompose the optical depth field $\tau(\bm{x})$ into a
spatial average $\bar{\tau}$, well measured by Planck through the large-scale
scattering effect ($\bar{\tau}=0.054\pm 0.007$ Planck2020A&A…641A…6P ), and a
fluctuating part $\delta\tau$, which traces the inhomogenous projected
distribution of free electrons on the sky. Reabsorbing the mean $\bar{\tau}$
into the definition of $\delta T^{0}$ and assuming small fluctuations
$\delta\tau(\bm{x})\ll 1$, the screening effect reduces to an additive term:

$\delta T^{\text{screening}}(\bm{x})=-\delta T^{0}(\bm{x})\delta\tau(\bm{x}).$
(4)

This effect is illustrated in Fig. 1. Our goal in this paper is to propose new
“temperature inversion” and “signed” estimators in order to reconstruct the
optical depth fluctuation field from the observed temperature map. Please note
that in the discussion that follows, we use notations $\tau$ and $\delta\tau$
interchangeably, and they both represent the fluctuations on top of the mean
$\bar{\tau}$ field.

Figure 1: Top: 1 deg2 simulated map of a CMB realization unscreened (left) and
screened by 5’ FWHM 2D Gaussian $\tau$ profiles placed at CMASS galaxy
positions (right). The $\tau$ amplitude is greatly exaggerated for visual
effect.
Bottom: A 1D schematic illustrating the small angular scale temperature
deviations $\delta T_{\rm small}$ induced by the patchy screening effect
(green). Unlike any other foreground, these deviations are (anti-)correlated
with the large-scale primary CMB anisotropies $\delta T_{\rm large}$ (black).
The tSZ- and CIB-induced anisotropies always have the same sign at a given
frequency; lensing is correlated with $\nabla\delta T_{\rm large}$, not
$\delta T_{\rm large}$ itself; and kSZ-induced anisotropies change sign with
the gas’s LOS velocity. Thus, stacking the observed CMB temperature patches
around galaxies weighted by the sign of the large-scale temperature provides a
means to detect screening. Moreover, the patchy screening effect has greater
magnitude where $\delta T_{\rm large}$ is large, motivating weighting by
$\delta T_{\rm large}$ (and not just its sign) when stacking.

##  III Standard Quadratic Estimator

###  III.1 Review: derivation of the standard QE

The standard QE for $\tau$ Hu:2001kj ; DvorkinSmith2009 ; Su2011.1106.4313S
is generally derived starting from

$\delta T(\bm{x})=\delta
T^{0}(\bm{x})\left[1-\tau(\bm{x})+\mathcal{O}(\tau^{2}(\bm{x}))\right].$ (5)

In Fourier space with a flat-sky approximation, this becomes

$\delta T_{\bm{\ell}}\simeq\delta
T^{0}_{\bm{\ell}}-\int_{{\bm{\ell}}^{\prime}}\delta
T^{0}_{{\bm{\ell}}^{\prime}}\,\tau_{{\bm{\ell}}-{\bm{\ell}}^{\prime}},$ (6)

where
$\int_{{\bm{\ell}}^{\prime}}\equiv\int\frac{d^{2}{\bm{\ell}}^{\prime}}{(2\pi)^{2}}$.
Thus screening produces non-zero off-diagonal covariances for the Fourier
modes:

$\langle\delta T_{\bm{\ell}}\delta
T_{{\bm{L}}-{\bm{\ell}}}\rangle_{\text{fixed
}\tau}=f^{\tau}({\bm{\ell}},{\bm{L}}-{\bm{\ell}})\tau_{\bm{L}}\,,$ (7)

where

$f^{\tau}({\bm{\ell}},{\bm{L}}-{\bm{\ell}})=-\left(C^{0}_{\ell}+C^{0}_{|{\bm{L}}-{\bm{\ell}}|}\right).$
(8)

An elementary QE for $\tau_{L}$ is thus simply $-\delta T_{\bm{\ell}}\delta
T_{{\bm{L}}-{\bm{\ell}}}/\left(C^{0}_{\ell}+C^{0}_{|{\bm{L}}-{\bm{\ell}}|}\right)$,
with variance $C^{\text{total}}_{\ell}\
C^{\text{total}}_{|{\bm{L}}-{\bm{\ell}}|}\
/\left(C^{0}_{\ell}+C^{0}_{|{\bm{L}}-{\bm{\ell}}|}\right)$. The standard QE
for $\tau$ is therefore obtained from the inverse-variance linear combination
of the elementary QEs, giving

$\displaystyle\hat{\tau}^{\text{QE}}_{\bm{L}}$ $\displaystyle\equiv
N_{{\bm{L}}}\int_{\bm{\ell}}F^{\tau}({\bm{\ell}},{\bm{L}}-{\bm{\ell}})\delta
T_{\bm{\ell}}\delta T_{{\bm{L}}-{\bm{\ell}}}\,$ (9)
$\displaystyle=-\frac{\int_{\bm{\ell}}\frac{\left(C^{0}_{\ell}+C^{0}_{|{\bm{L}}-{\bm{\ell}}|}\right)}{C^{\text{total}}_{\ell}\
C^{\text{total}}_{|{\bm{L}}-{\bm{\ell}}|}}\ \delta T_{\bm{\ell}}\delta
T_{{\bm{L}}-{\bm{\ell}}}}{\int_{\bm{\ell}}\frac{\left(C^{0}_{\ell}+C^{0}_{|{\bm{L}}-{\bm{\ell}}|}\right)^{2}}{C^{\text{total}}_{\ell}\
C^{\text{total}}_{|{\bm{L}}-{\bm{\ell}}|}}}.$

where

$F^{\tau}({\bm{\ell}},{\bm{L}}-{\bm{\ell}})=\frac{f^{\tau}({\bm{\ell}},{\bm{L}}-{\bm{\ell}})}{C^{\text{total}}_{\ell}\
C^{\text{total}}_{|{\bm{L}}-{\bm{\ell}}|}}\,,$ (10)

and

$(N_{{\bm{L}}})^{-1}=\int_{\bm{\ell}}F^{\tau}({\bm{\ell}},{\bm{L}}-{\bm{\ell}})f^{\tau}({\bm{\ell}},{\bm{L}}-{\bm{\ell}})\,,$
(11)

where the weights $F^{\tau}({\bm{\ell}},{\bm{L}}-{\bm{\ell}})$ are determined
such that the variance of Eq. 9, $\mathcal{N^{\tau}_{\bm{L}}}$, given below,
is minimized.

$\mathcal{N^{\tau}_{\bm{L}}}=2(N_{\bm{L}})^{2}\int_{\bm{\ell}}\left(F^{\tau}({\bm{\ell}},{\bm{L}}-{\bm{\ell}})\right)^{2}C^{\text{total}}_{\ell}\
C^{\text{total}}_{|{\bm{L}}-{\bm{\ell}}|}\,.$ (12)

###  III.2 Intuition: large-scale VS small-scale $\tau$ reconstruction

The standard QE derived above (Eq. 9) can be understood more simply in two
limiting cases, where large-scale and small-scale Fourier modes are being
reconstucted.

In the former limit, where large-scale $\tau_{\bm{L}}$ Fourier modes are
reconstructed from small-scale temperature modes $\delta T_{\bm{\ell}}$, with
$L\ll\ell$, the standard QE is effectively looking at the local power spectrum
and searching for a change of the form

$C_{\ell}\leftarrow C_{\ell}e^{-2\delta\tau},\quad\text{i.e.},\quad\frac{d\ln
C_{\ell}}{d\tau}=-2.$ (13)

This is illustrated in Fig. 2.

In the latter limit, where we reconstruct the optical depth field
$\tau_{\bm{L}}$ on small scales, i.e., $L\gg\ell$, we can approximate the
power spectra as smoothly-varying functions:
$C^{0}_{|{\bm{L}}-{\bm{\ell}}|}\simeq C^{0}_{L}$ and
$C^{\text{total}}_{|{\bm{L}}-{\bm{\ell}}|}\simeq C^{\text{total}}_{L}$ to
lowest order. Furthermore, the CMB power spectrum is very red, such that
$C^{0}_{\ell}\gg C^{0}_{L}$. In this limit, the standard QE takes the
suggestive form

$\hat{\tau}^{\text{QE}}_{\bm{L}}\simeq-\frac{\int_{\bm{\ell}}\left(\frac{C^{0}_{\ell}}{C^{\text{total}}_{\ell}}\delta
T_{\bm{\ell}}\right)\delta
T_{{\bm{L}}-{\bm{\ell}}}}{\int_{\bm{\ell}}\frac{\left(C^{0}_{\ell}\right)^{2}}{C^{\text{total}}_{\ell}}},$
(14)

or in configuration space:

$1-\hat{\tau}^{\text{QE}}(\bm{x})\simeq\delta T(\bm{x})\ \frac{\delta
T^{\text{large-WF}}(\bm{x})}{\langle\left(\delta T^{\text{large-
WF}}\right)^{2}\rangle}.$ (15)

In this limit, the standard QE is thus simply multiplying the local small-
scale temperature $\delta T(\bm{x})$ by the Wiener-filtered large-scale
temperature $\delta T^{\text{large-WF}}$, and normalizing it by the expected
variance of the Wiener filter. The significance of this expression will become
clear in Sec. IV, where we connect it to the new TI and signed estimators.

Figure 2: Top: A large-scale optical depth fluctuation $\delta\tau$ spanning
10∘ in RA.
Middle: A $10^{\circ}\times 5^{\circ}$ simulated temperature map of a CMB
realization screened by the large-scale $\delta\tau$ fluctuation shown in the
top panel. The dotted (dashed) square corresponds to the area with suppressed
(enhanced) anisotropies with respect to the full-sky average optical depth
$\bar{\tau}$.
Bottom: The power spectra corresponding to the two regions where
$\delta\tau>0$ (dotted) and $\delta\tau<0$ (dashed). The power spectrum for
the CMB screened only by $\bar{\tau}$ (i.e. $\delta\tau=0$, solid) is shown
for comparison. The screening of the CMB by a large-scale mode of the optical
depth modulates the locally measured power spectrum as $C_{\ell}\leftarrow
C_{\ell}e^{-2\delta\tau}$. The standard QE for $\delta\tau_{\bm{L}}$ is
looking precisely for this effect when reconstructing from temperature modes
$\delta T_{\bm{\ell}}$ with $L\ll\ell$.

##  IV “Temperature-inversion” estimator & approximate small-scale QE

###  IV.1 New “temperature-inversion” estimator

The optical depth fluctuation field possesses power on very small scales since
it traces galaxy halos, whereas the primary CMB temperature has almost none
due to Silk damping. Denoting small and large scales with subscripts $S$ and
$L$ respectively, the screening effect on small scales is thus almost
exclusively sourced by the following configurations:

$\delta T^{\text{screening}}_{S}(\bm{x})=-\delta
T^{0}_{L}(\bm{x})\delta\tau_{S}(\bm{x}).$ (16)

This separation of scales is crucial: thanks to it, both the unscreened and
screened temperature fluctuations can be measured from the same temperature
map Natarajan2013.776.82N .

Indeed, we can low- and high-pass filter the observed temperature map to
estimate the large and small scale temperature fluctuations
$\delta\hat{T}_{L}(\bm{x})$ and $\hat{\delta T_{S}}(\bm{x})$. A natural, if
somewhat singular, estimator for $\tau$ is then obtained by simply taking the
ratio:

$\hat{\tau}^{\text{TI}}(\bm{x})\equiv-\frac{\delta\hat{T}_{S}(\bm{x})}{\delta\hat{T}_{L}(\bm{x})}.$
(17)

We call this the “temperature inversion” (TI) estimator in analogy with the
gradient inversion (GI) estimator of CMB lensing Horowitz2019.485.3919H ;
Hadzhiyska2019.100b3547H . A natural choice of low-pass filter in this case
would be a Wiener filter, the linear estimator for $\delta T^{0}_{L}(\bm{x})$
that has the lowest mean-squared error. A natural choice for the high-pass
filter would be an inverse-variance weighted average of the small scales, the
unbiased linear estimator of $\delta\hat{T}_{S}(\bm{x})$ with the lowest
variance. In general, we can define linear filters $W_{L}$ and $W_{S}$ such
that

$\left\\{\begin{aligned} &\hat{\delta T_{L}}(\bm{x})=\int
d\bm{x}^{\prime}W_{L}(\bm{x}-\bm{x}^{\prime}){\delta T}(\bm{x}),\\\
&\hat{\delta T_{S}}(\bm{x})=\int
d\bm{x}^{\prime}W_{S}(\bm{x}-\bm{x}^{\prime}){\delta T}(\bm{x}).\\\
\end{aligned}\right.$ (18)

In practice, we apply low- and high-pass filters (see Sec. VI.1) to the
temperature map to isolate $\delta T_{L}$ and $\delta T_{S}$, respectively. If
the shape of the $\tau$ profile for the objects of interest is known (e.g.
from kSZ, which has higher SNR), we can use a matched filter with this profile
to estimate the amplitude of the signal from the data. However, our filters,
with variable size, allow us to actually measure the scale dependence of the
electron density profile. When jointly analysed, including their covariance
matrix, they are as statistically optimal as the matched filter.

In App. A, we further derive the filter $W_{L}$ which minimizes the mean-
squared error in the presence of foreground contamination with known profiles
(a combination of matched filter and Wiener filter), but do not implement it
in what follows.

###  IV.2 Comparison to the approximate small-scale QE

In the limit where we are reconstructing small-scale optical depth
fluctuations, we have shown that the QE takes the approximate form

$1-\hat{\tau}^{\text{QE}}(\bm{x})\simeq\delta T(\bm{x})\ \frac{\delta
T^{\text{large-WF}}(\bm{x})}{\left<\left(\delta T^{\text{large-
WF}}\right)^{2}\right>},$ (19)

where $\delta T^{\text{large-WF}}$ is the Wiener filter estimator for the
unscreened CMB, given the map noise. This can be contrasted with our TI
estimator of the form

$1-\hat{\tau}^{\text{TI}}(\bm{x})\simeq\frac{\delta T(\bm{x})}{\delta
T^{\text{large-WF}}(\bm{x})}.$ (20)

In other words, in the small-scale limit, the two estimators are related via

$1-\hat{\tau}^{\text{QE}}(\bm{x})=\left[1-\hat{\tau}^{\text{TI}}(\bm{x})\right]\frac{\left(\delta
T^{\text{large-WF}}(\bm{x})\right)^{2}}{\left<\left(\delta T^{\text{large-
WF}}\right)^{2}\right>}.$ (21)

Thus, even if the large-scale temperature fluctuation $\delta T^{\text{large-
WF}}$ is measured with infinite precision, the QE suffers from an irreducible
variance due to the ratio between $(\delta T^{\text{large-WF}}(\bm{x}))^{2}$
and its mean. This is completely analogous to the GI estimator of CMB lensing.
In App. B, we expand on this comparison, showing that the SNR is higher for TI
than the QE when the large-scale temperature leg is measured at high SNR
(already true for Planck data) and the small-scale temperature leg is
dominated by the screening effect. This latter condition is likely not
satisfied; even in the absence of detector noise, and assuming that all non-
blackbody foregrounds have been subtracted, kSZ and CMB lensing will likely
dominate the noise unless they can be largely subtracted.

If TI does not outperform the QE in terms of SNR when considering the same set
of temperature Fourier modes, why consider it further? As we show in Secs. V.2
and V.4, TI is more robust to extragalactic foreground biases and naturally
robust to CMB lensing bias. This means that TI can use more temperature
Fourier modes (higher $\ell_{\text{max}}^{T}$) and thus potentially recover a
higher total SNR while still being as or more robust to biases.

###  IV.3 Stacked TI estimator for cross-correlation & comparison with small-
scale QE

For a single galaxy, the naïve TI estimator (Eq. (17)) is clearly singular for
lines of sight where the unscreened temperature fluctuation is zero. This
issue is automatically addressed when stacking around a set of galaxies, as we
now show. The large-scale temperature $\hat{\delta T_{L}}(\bm{x})$ can be
measured to high precision from the temperature map, such that we ignore the
noise on $\hat{\delta T_{L}}(\bm{x})$ for now. On the other hand, the small-
scale temperature $\hat{\delta T_{S}}(\bm{x})$ is typically affected by a non-
negligible noise $n_{S}$ (with variance $\sigma_{n_{S}}^{2}$), such that

$\hat{\tau}^{\text{TI}}(\bm{x})=\tau(\bm{x})+\frac{n_{S}(\bm{x})}{\hat{\delta
T_{L}}(\bm{x})}.$ (22)

In this case, the estimator is unbiased, and its noise is smallest along lines
of sight with the highest large-scale temperature, as expected intuitively.

To extract the mean optical depth profile of a sample of tracers (galaxies,
galaxy groups or clusters), we look for the unbiased linear combination of the
$\hat{\tau}^{\text{TI}}$ around each cluster that has the smallest noise
variance. The answer is the inverse-weighted average:

$\hat{\tau}^{\text{TI
stack}}\equiv\frac{\sum_{i}\hat{\tau}^{\text{TI}}_{i}\times\hat{\delta
T_{L}}_{i}^{2}/\sigma_{n_{S}}^{2}}{\sum_{i}\hat{\delta
T_{L}}_{i}^{2}/\sigma_{n_{S}}^{2}},$ (23)

which can be rewritten as

$\hat{\tau}^{\text{TI stack}}=-\frac{\sum_{i}\hat{\delta T}_{Si}\hat{\delta
T_{Li}}/\sigma_{n}^{2}}{\sum_{j}\hat{\delta T_{L}}_{j}^{2}/\sigma_{n}^{2}}.$
(24)

Conveniently, the inverse-variance weighting regularizes the singularity of
the naïve estimator, such that lines of sight where the unscreened CMB
temperature $\hat{\delta T^{0}}$ is estimated to vanish are given zero weight
in the sum.

Although the numerator of Eq. (24) is now quadratic in the temperature map,
the estimator as a whole is not since the denominator also involves a
quadratic combination of the data. As derived in App. C, the stacked QE takes
a very similar form in the limit of small-scale screening:

$\hat{\tau}^{\text{QE stack}}=-\frac{\sum_{i}\hat{\delta T}_{Si}\hat{\delta
T_{Li}}/\sigma_{n}^{2}}{\sum_{j}\sigma_{T_{L}}^{2}/\sigma_{n}^{2}}.$ (25)

The numerator is identical to that of the stacked TI, and the denominator only
differs by replacing $\hat{\delta T_{L}}^{2}_{j}$ with $\sigma_{T_{L}}^{2}$.
Thus, the ratio of the stacked QE to the stacked TI is simply

$\frac{\hat{\tau}^{\text{QE stack}}}{\hat{\tau}^{\text{TI
stack}}}=\frac{\sum_{j}\hat{\delta
T_{L}}_{j}^{2}/\sigma_{n}^{2}}{\sum_{j}\sigma_{T_{L}}^{2}/\sigma_{n}^{2}}.$
(26)

This ratio has mean 1, but its variance is non-zero, equal to
$2/N_{\text{galaxies}}$. Thus, in the limit of zero temperature noise, the
stacked QE retains an irreducible variance
$2\sigma_{\tau}^{2}/N_{\text{galaxies}}$, whereas the stacked TI yields a
noiseless reconstruction of the true $\tau$. While this variance vanishes as
$\propto 1/N_{\text{galaxies}}$ for large galaxy samples, the other sources of
variance (e.g., from detector noise) also share the same scaling, such that
this additional noise remains relevant for all sample sizes.

###  IV.4 “Signed estimator”: suboptimal to TI, but potentially even more
robust to foregrounds

The key intuition behind the derivation of the stacked TI estimator (Eq. (24))
is shown in Fig. 1: unlike all other imprints of galaxies on the CMB, the
screening signal changes sign and amplitude depending on those of the local
large-scale temperature $\delta T_{L}$. Discarding the amplitude information
and weighting only by the sign of $\delta T_{L}$ leads to a different
screening estimator. We name it the “signed” estimator:

$\hat{\tau}^{\rm sgn}\equiv\frac{\sum_{i}\operatorname{sgn}(\delta
T_{L,i})\delta T_{S,i}}{\sum_{i}|\delta T_{L,i}|},$ (27)

where $\delta T_{L,i}$ (or $\delta T_{S,i}$) is the large-scale (or small-
scale) temperature fluctuation at the position of galaxy $i$.

While the TI and signed estimators are both unbiased estimators of $\tau$,
their variances differ, such that

$\frac{\text{var}\ \hat{\tau}^{\rm sgn}}{\text{var}\ \hat{\tau}^{\rm
TI}}=\frac{\langle\left(\delta T_{L,i}\right)^{2}\rangle}{\langle|\delta
T_{L,i}|\rangle^{2}}=\frac{\pi}{2},$ (28)

for the Gaussian distributed large-scale temperatures $\delta T_{L}$. Thus,
the signed estimator has lower SNR than TI by a factor $\sqrt{\pi/2}=1.25$.

On the other hand, this signed estimator enables two interesting avenues for
controlling potential foregrounds, which may be worth the cost in SNR:

  * •

Sign balancing. We can discard CMB patches on galaxies with a certain sign to
ensure that $\sum_{i}\operatorname{sgn}(\delta T_{L,i})=0$. Doing so exactly
nulls any part of the foreground emission (e.g., CIB, tSZ) that is constant
across galaxies. This does not lead to any bias in the estimated screening
since this selection depends on the large-scale temperature only,
statistically independent from the screening signal. Furthermore, this
typically only requires discarding a fraction $\sim 1/N_{\text{gal}}$ of the
galaxy sample, completely negligible for the large galaxy catalogs of interest
here. This is described in more detail in App. D.

  * •

Sign thresholding. After masking bright point sources and clusters from the
temperature map, any residual foreground emission from our galaxies will be
small ($\sim\mu$K) compared to the typical large-scale temperature
fluctuations ($\sim 100\mu$K for the CMB). Thus, we may select for the stack
only the objects for which $|\delta T_{L,i}|$ is larger than a few times the
typical foreground signal (e.g., ($\sim 10\mu$K)). This ensures that residual
foregrounds are not large enough to affect the sign weighting used in the
signed estimator. This is in contrast with TI, where residual foregrounds in
$\delta T_{L,i}$, however small, remain present. This is described in more
detail in App. E.

Including tests on simulations for biases from extragalactic sources, this
estimator is presented in detail in Coulton_2023 , where we test for
foreground biases and validate these arguments by applying the signed
estimator to realistic, non-Gaussian sky simulations.

##  V Biases to QE & TI from CMB lensing and extragalactic foregrounds

###  V.1 CMB lensing bias to QE: lens hardening needed

From Eq. (9), we can see that the QE for $\tau$ is of the form

$\hat{\tau}^{\text{QE}}_{\bm{L}}\equiv-
N_{L}\int_{\bm{\ell}}\frac{\left(C^{0}_{\ell}+C^{0}_{|{\bm{L}}-{\bm{\ell}}|}\right)}{C^{\text{total}}_{\ell}\
C^{\text{total}}_{|{\bm{L}}-{\bm{\ell}}|}}\ \delta T_{\bm{\ell}}\delta
T_{{\bm{L}}-{\bm{\ell}}}\,.$ (29)

This has a similar form as that of the QE for the lensing potential:
$\hat{\phi}^{\text{QE}}_{\bm{L}}\propto\int_{\bm{\ell}}F({\bm{\ell}},{\bm{L}}-{\bm{\ell}})\delta
T_{\bm{\ell}}\delta T_{{\bm{L}}-{\bm{\ell}}}$, where
$F({\bm{\ell}},{\bm{L}}-{\bm{\ell}})$ are optimized weights for the pair of
multipoles $\delta T_{\bm{\ell}}$ and $\delta T_{{\bm{L}}-{\bm{\ell}}}$, such
that the variance of the estimator is minimized. Thus, Eq. (9) can be
interpreted as an estimator for the CMB lensing potential reconstruction with
non-optimized weights, and as such, has a non-zero response to lensing. We can
rewrite Eq. (29) as

$\langle\hat{\tau}^{\text{QE}}_{\bm{L}}\rangle={\tau}^{\text{QE}}_{\bm{L}}+N_{L}\underbrace{\left(\int_{\bm{\ell}}F^{\tau}({\bm{\ell}},{\bm{L}}-{\bm{\ell}})f^{\phi}_{TT}({\bm{\ell}},{\bm{L}}-{\bm{\ell}})\right)}_{\text{Response
of QE to lensing}}\phi({\bm{L}})\,,$ (30)

where we have used $\langle\delta T_{\bm{\ell}}\delta
T_{{\bm{L}}-{\bm{\ell}}}\rangle=f^{\phi}_{TT}({\bm{\ell}},{\bm{L}}-{\bm{\ell}})\phi({\bm{L}})$
with
$f^{\phi}_{TT}({\bm{\ell}},{\bm{L}}-{\bm{\ell}})=C^{0}_{\ell}({\bm{L}}\cdot{\bm{\ell}})+C^{0}_{|{\bm{L}}-{\bm{\ell}}|}({\bm{L}}\cdot({\bm{L}}-{\bm{\ell}}))$,
and $F^{\tau}({\bm{\ell}},{\bm{L}}-{\bm{\ell}})$ are given by Eq. (10). In
practice, the lensing bias in Eq. (30) is non-zero and can be dominant even in
cross-correlations with other tracers, as shown in Tab. 1 and Fig. 3 (see also
Fig. 1 of Namikawa_2021 ).

This large lensing bias to the QE is typically avoided via “lens hardening”
Namikawa_2013 . The idea behind the lens-hardened QE is simple. Just as we can
construct a $\hat{\tau}$ QE, we can also construct a lensing convergence
$\hat{\kappa}$ QE, which can then be used to subtract off the lensing bias in
Eq. (30). Therefore, the lens-hardened $\tau$ QE minimizes contamination from
lensing. However, this comes at an increased reconstruction noise cost when
compared with the standard QE variance $\mathcal{N^{\tau}_{\bm{L}}}$ given by
Eq. 12. For a Planck-like experiment, this noise cost is $\sim 0-40\%$
depending on scales Namikawa_2021 , whereas the reconstruction noise from
temperature-only (no polarization) QE increases by a factor of $\sim 2-3$
across all $\ell$ values for a CMB-S4-like experiment Roy_2023 .

Intuitively, this lensing bias to QE can be understood from Fig. 2 when
reconstructing large-scale optical depth fluctuations $\tau_{\bm{L}}$ from
small-scale temperature modes $T_{\bm{\ell}}$ with $L\ll\ell$. Both lensing
and screening cause variations in the local power spectrum from patch to
patch. While these variations are different, they have non-zero overlap.

The opposite regime, where we reconstruct small-scale optical depth
fluctuations $\tau_{\bm{L}}$ from temperature modes $T_{\bm{\ell}}$ and
$T_{{\bm{L}}-{\bm{\ell}}}$ with $L\gg\ell$, is different. In this limit, the
lensing term in Eq. (30) vanishes. Thus, the small-scale part of the QE is
unbiased with respect to lensing. This part of the QE is very similar to the
TI estimator, only with a slightly different normalization (see Eqs.
(19)-(21)), such that both are unbiased by lensing. We explore the robustness
of TI to CMB lensing in more detail in the next subsection.

###  V.2 TI is robust to CMB lensing

The lensed temperature map can be written, to first order, as

$\delta\tilde{T}(\bm{x})\approx\delta T(\bm{x})+\nabla\phi(\bm{x})\nabla\delta
T(\bm{x}),$ (31)

where $\tilde{T}$ is the lensed temperature field and $\phi$ is the lensing
potential. Thus, lensing by a galaxy of interest produces a local small-scale
dipole in the CMB map, whose direction and amplitude matches those of the
local large-scale unlensed CMB gradient. The spatial symmetry of this dipole
effect means that lensing is nulled by our azimuthally symmetric small-scale
temperature filter, typically a disk or ring aperture around the galaxy, which
only extracts the local small-scale monopole. Thus, to first order, lensing
from the galaxies in our sample does not bias the TI or signed estimators and
does not even contribute noise.

The situation is different for the lensing by other objects than our galaxies.
Since these small-scale dipoles are no longer aligned with our small-scale
apertures (rings or disks), lensing does enter the TI and signed estimators in
the form of noise. In fact, on small scales $\ell\gtrsim 5000$, the CMB power
spectrum is dominated by lensing, meaning that lensing may be the dominant
source of noise (other than kSZ, other foregrounds and detector noise
potentially present in the map).

In summary, TI and signed estimators null the lowest order bias due to lensing
$\sim\mathcal{O}\left(\phi g\right)$. At higher order in lensing
($\mathcal{O}\left(\phi^{2}\right)$), a bias to the TI and signed estimators
may be present if

  * •

lensing is non-negligible in the large-scale temperature map as well;

  * •

this lensing is not due to the galaxy we are centered on, such that the
corresponding dipole is not nulled by our small-scale aperture filter;

  * •

the lensing is still caused by objects correlated with the galaxy we are
centered on, e.g., two-halo term or other correlated galaxies in the sample.

We leave the exploration of these higher order bias terms in simulations to
future work.

###  V.3 Extragalactic foregrounds bias QE

Similar to CMB lensing, foregrounds like tSZ, kSZ and CIB cause mode coupling
of temperature anisotropies. As a result, they bias the QE for screening.
Again, similar to CMB lensing, this bias also affects cross-correlation of the
reconstructed $\tau$ field with other tracers. Fig. 3 shows the expected
cross-correlation of the $\tau$ field reconstructed from a QE for CMASS
galaxies using an Advanced ACT-like experiment with CMASS galaxies. Here, “tSZ
high” and “tSZ low” are two tSZ maps we simulate which give an upper and lower
bound respectively on this bias (see App. F for details). We can see that the
bias from CMB lensing is higher than the signal itself. The bias from tSZ
emission is not negligible either. While we only discuss the foreground bias
due to CMB lensing and tSZ here, foregrounds like CIB could also have non-
negligible bias to $\tau$ QE. Such foreground biases have been extensively
studied and shown to be important in the context of CMB lensing QE
vanEngelen_2014 ; Osborne_2014 . Similar conclusions hold for the $\tau$ QE as
well due to the similarity between the two estimators. It should also be noted
that these biases are likely bigger for lower noise CMB experiments as they
will have increased weights from small scales where the foregrounds are
dominant.

Figure 3:  The expected cross-correlation signal between BOSS CMASS galaxies
and the screening QE from ACT temperature (black) is compared to the CMB
lensing bias (red) and the tSZ foreground bias (cyan area) for a range of tSZ
models. The importance of these biases highlight the need for lens hardening
and foreground hardening for the $\tau$ QE. The corresponding biases to the
overall amplitude of the signal are shown in Tab. 1. The simulations used to
create these correlated lensing and tSZ maps are described in detail in App.
F.

In order to quantify the bias from foregrounds like CMB lensing and tSZ to the
$\tau$ QE in cross-correlation with galaxies, we define an amplitude parameter
$A$ for a given power spectrum measurement $\hat{C}_{L}$ such that

$\hat{C}_{L}=A\ C^{\rm theory}_{L}+n_{L}\,,$ (32)

where $C^{\rm theory}_{L}$ is the theory power spectrum and $n_{L}$ is noise.
The minimum-variance estimator for $A$, linear in the measured power spectrum,
is then

$\hat{A}=\frac{\sum_{L}\hat{C}_{L}C^{\rm
theory}_{L}/{\sigma_{L}}^{2}}{\sum_{L}{C^{\rm
theory}_{L}}^{2}/{\sigma_{L}}^{2}}\,.$ (33)

The nominal value of $\hat{A}$ is 1; however, due to biases from foregrounds,
the value differs from 1. Thus, bias from a foreground $X$ to cross-
correlation of the QE reconstructed $\tau$ field with galaxies $\delta_{g}$
($C_{L}^{\hat{\tau}\times\delta_{g}}$) can be calculated in the continuum
limit as

${\rm
bias}(C^{\hat{\tau}\times\delta_{g}})=\frac{\int\frac{d^{2}{\mathbb{L}}}{(2\pi)^{2}}\frac{{\rm
bias}(C_{L}^{X\times\delta_{g}})}{C_{L}^{\hat{\tau}\times\delta_{g}}}\frac{(C_{L}^{\hat{\tau}\times\delta_{g}})^{2}}{\sigma_{L}^{2}}}{{\int\frac{d^{2}{\mathbb{L}}}{(2\pi)^{2}}\frac{(C_{L}^{\hat{\tau}\times\delta_{g}})^{2}}{\sigma_{L}^{2}}}}\,,$
(34)

where ${\rm bias}(C_{L}^{X\times\delta_{g}})$ is the cross-correlation with
the galaxies of the foreground $X$ passed through a $\tau$ QE. Tab. 1 shows
that these foreground biases can be larger than the screening signal and
highly significant. Some form of foreground mitigation is therefore likely
required for the QE. This may include multi-frequency component separation,
masking, template deprojection, using shear only estimators, and bias
hardening with respect to point sources or specific profiles. Such approaches
and their challenges have been explored in the context of CMB lensing and
$\tau$ QE Namikawa_2013 ; Namikawa_2021 ; MacCrann_2023 ; Schaan_2019 ;
Sailer_2020  and similar results hold here too. We do not explore these
further in this paper.

| Galaxy
---
sample
| Fore-
---
ground
| SPT
---
| ACT
---
| SO
---
| CMB-S4
---
| CMASS
---
| $\kappa$
---
tSZ low
tSZ high
| 410.9%
---
-2.7%
-39.7%
| 346.9%
---
-2.5%
-28.4%
| 394%
---
-2.7%
-36.5%
| 433.5%
---
-2.8%
-44.2%
Table 1:  The amplitude of the cross-correlation of galaxies with the
screening QE suffers large fractional biases from CMB lensing and tSZ, using
$30<\ell<3000$. App. F.2 describes the “low” and “high” tSZ modeling in
detail.

###  V.4 TI is robust to extragalactic foregrounds

The TI and signed estimators are automatically robust to extragalactic
foreground biases, and this can be further enhanced by sign balancing and sign
thresholding for the signed estimator. We summarize the reasons here and
present more detailed derivations in App. D and E.

Unlike the patchy screening effect, foregrounds like the CIB, tSZ and kSZ do
not change sign depending on the sign of the local large-scale temperature
$\delta T_{L}$. In the absence of foreground contamination to the large-scale
leg, the weights in the TI and signed estimators, which change sign and have
zero mean value, automatically average out these foregrounds. For a finite
galaxy sample size, this foreground reduction may leave a residual bias $\sim
f_{S}/\sqrt{N_{\text{patch}}}$, where $f_{S}$ is the typical small scale
foreground emission, and $N_{\text{patch}}$ is the number of patches with
independent large-scale CMB temperature $\delta T_{L}$. Since the CMB is
coherent on degree scales, $N_{\text{patch}}$ is roughly equal to the area of
the CMB map in square degrees.

For the signed estimator, this residual can be further reduced by sign
balancing the weights, as described in Sec. IV.4. If the foreground signal is
identical for all galaxies in the sample, sign balancing exactly nulls the
foreground bias. While foreground emissions actually vary across galaxies,
sign balancing is still effective in realistic scenarios, as shown in App. D

In practice, foreground contamination may be present both in the large and
small-scale temperature legs. In this case, the TI and signed estimators may
suffer a bispectrum bias $\sim f_{L}f_{S}g_{S}$, where $f$ is the foreground
and $g$ the galaxy map. However, for the signed estimator, this bias can be
nulled exactly with the sign thresholding procedure presented in Sec. IV.4.
This procedure ensures that absolutely no foreground signal leaks into
$\text{sgn}\left(\delta T_{L}\right)$, thus removing any bispectrum foreground
bias.

##  VI SNR forecasts

We produce SNR forecasts for the TI and $\tau$ QE estimators using realistic
simulations for four CMB experiments, South Pole Telescope (SPT) Benson_2014 ;
2018SPIE10708E..03B ; Sobrin_2022 , Atacama Cosmology Telescope (ACT)
2016JLTP..184..772H , Simons Observatory (SO) SOScience  and CMB-S4
2019arXiv190704473A , cross-correlated with three galaxy survey samples, the
Baryonic Oscillation Spectroscopic Survey (BOSS) CMASS sample Dawson_2013 ;
Reid_2016 , the Wide-field Infrared Survey Explorer (WISE) unWISE “blue” and
“green” samples Wright_2010 ; Meisner_2018 , and the simulated Legacy Survey
of Space and Time (LSST) Data Challenge 2 (DC2) extragalactic sample DC2_2021
. Details of the sky fraction, galaxy number density and number of galaxies
overlapping the CMB survey footprints are listed in Tab. 2.

###  VI.1 Simulated analysis for TI and QE

Here we give an overview of the methods used to generate the simulations for
the TI and QE detection forecast. See App. G for further details. We create
realistic simulated $\tau$ signal maps by placing identical, isotropic
Gaussian $\tau$ profiles at CMASS galaxy positions as described in App. G.1.
The parameters of the $\tau$ profiles are matched to the best-fit values from
the ACT$\times$unWISE first measurement Coulton_2023 .

To calculate the covariance matrix, we simulate a set of 128 maps for each CMB
experiment, using the power spectra shown in Fig. 4. All of the spectra are
beamed and include noise contributions from the CMB, detector noise and–with
the exception of CMB-S4–extragalactic foreground contamination (e.g. lensing,
CIB, tSZ, kSZ, radio sources). More details may be found in App. G.2.

Figure 4: Power spectra used in the SNR forecasts for ACT (solid orange), SPT
(dashed green), SO (dash-dotted blue), and CMB-S4 (dotted purple). Details for
the spectra are given in App. G. The shaded bands show the $\ell$ ranges used
by the TI estimator: $30<\ell\lesssim 2075$ (grey) and $2425\lesssim\ell<5950$
(green) for large- and small-scale temperature measurements, respectively. See
Eqs. 83 and 84 for the exact $\ell$ filtering used. The QE uses $30<\ell<3000$
(black dotted line) as it is sensitive to foreground contamination at higher
$\ell$. Since the TI estimator is robust to foregrounds at high $\ell$, it can
take advantage of future experiments’ lower detector noise and smaller beam
widths.

To separate the large- and small-scale temperature fluctuations $\delta T_{L}$
and $\delta T_{S}$, we apply low- and high-pass filters (see Eqs. (83) and
(84)) to generate two separate maps. We have implemented the stacked TI and
signed estimators in the aperture photometry filtering and stacking software,
ThumbStack (Schaan2021a, ). We use this code to generate small postage stamps
centered on each galaxy position for both maps, apply annular (i.e. “ring”)
aperture photometry filters to measure the $\delta T_{S}$ radial profile for
each stamp, and stack the stamps according to Eq. (24), using each stamp’s
average $\delta T_{L}$ as weights.

We construct the covariance matrix taking into account the clustering of
galaxies and scaling by the proper Hartlap correction factor Hartlap2007 .
Details are provided in App. G.4.

We then measure the SNR for each CMB-galaxy sample combination, scaling to
each experiment’s full galaxy sample and sky fraction, $f_{\rm sky}$, shown in
Tab. 2. The resulting SNR values for the surveys are shown in Tab. 4.

For the QE of $\tau$, details of the SNR estimation for the cross-power
spectrum of the $\tau$ map reconstructed using QE and galaxy map $g$ are
presented in App. H. The SNR is calculated analytically as given in Eq. (88)
where the noise terms contain the total power spectrum of the reconstructed
$\tau$ map along with its reconstruction noise and the total galaxy power
spectrum, which includes both the clustering and shot noise contributions. The
shot noise levels for each galaxy sample are calculated using the number
density values provided in Tab. 2.

| Galaxy
---
sample
| $\bar{n}$
---
$\rm{arcmin}^{-2}$
| SPT
---
($f_{\rm sky}=0.036$)
| ACT/SO/CMB-S4
---
($f_{\rm sky}=0.315$)
| CMASS
---
| 0.012
---
| 65,600
---
| 568,776
---
| unWISE
---
| 1.24
---
| 6.8M
---
| 59M
---
| LSST
---
| 34.1
---
| 187M
---
| 1.62B
---
Table 2: Number of galaxies assumed to overlap with each CMB experiment used
in the SNR forecasts. The shot noise, $\bar{n}$, for each galaxy sample and
the sky fraction, $f_{\rm sky}$, for each CMB experiment are also listed. ACT,
SO, and CMB-S4 are assumed to have equal $f_{\rm sky}$.

###  VI.2 Galaxy clustering effects

As we go towards higher density galaxy samples, the shot noise of the sample
decreases by $\sqrt{N_{\rm gal}}$, but the clustering power spectrum still
adds noise to the covariance. Here we assess the impact of clustering effects
for the TI and QE forecasts.
For the QE, we complete the simulated analysis with and without the galaxy
clustering ($gg$) contribution at the power spectrum level, where the total
galaxy power spectrum is calculated with and without the clustering power
spectrum. The results of this study are summarized in Tab. 3.

|  |  $gg=0$ / $gg\neq 0$
---|---|---
| CMB
---
Expt.
| Estimator
---
| CMASS
---
|  unWISE
---
|  LSST
---
SPT |  | TI5950
---
QE6000
| 0.34 / 0.37
---
0.54 / 0.53
| 3.4 / 3.0
---
5.5 / 4.7
| 18 / 11
---
29 / 12
ACT |  | TI5950
---
QE6000
| 0.60 / 0.65
---
1.1 / 1.1
| 6.1 / 5.2
---
11 / 9.3
| 31 / 18
---
58 / 22
SO |  | TI5950
---
QE6000
| 0.84 / 0.92
---
1.4 / 1.4
| 8.5 / 7.5
---
14 / 12
| 45 / 25
---
73 / 29
CMB-S4 |  | TI5950
---
QE6000
| 1.5 / 1.7
---
2.7 / 2.7
| 16 / 13
---
28 / 25
| 82 / 49
---
145 / 65
Table 3:  Summary of the detection SNRs for TI and QE varying the noise
contribution due to galaxy clustering, $gg$. The first value in each pair
corresponds to the case with no clustering ($gg=0$) power spectrum
contribution to the total galaxy power spectrum; the second value includes it
($gg\neq 0$). For QE, we use $60<\ell<6000$. For TI, we use $30<\ell\lesssim
2075$ for $\delta T_{L}$ weights and $2425\lesssim\ell<5950$ for $\tau$
reconstruction. Ratios of the $gg=0$ and $gg\neq 0$ SNR values are visualized
in Fig. 5.

For the TI estimator, we compare the SNR values obtained with a covariance
matrix estimated from bootstrapping from a single map realization (thus
neglecting clustering effects) versus that estimated from the sample variance
over 128 map realizations.

As shown in Fig. 5, neglecting clustering effects results in overestimation of
the SNR for both the QE and TI. While the CMASS sample is sparse enough to not
show significant differences, using unWISE results in a moderate $\sim$20%
overestimated SNR. Ignoring clustering in the LSST sample is a more dramatic
oversight: $\sim$60% for TI and $\sim$150% for QE. The marked differences in
sensitivity to clustering effects in TI and QE, especially for an LSST-like
experiment, is unclear. We leave an investigation of this difference to future
work. Despite these uncertainties, it is clear clustering has a significant
effect on noise estimation. Thus, we include these effects in the present
forecasts and emphasize they must be properly accounted for in future
analyses.

We also study the effect of neglecting the cosmic variance of the signal
itself, i.e., $C_{\ell}^{\tau g}$, and find it is negligible for all the
experiments considered.

Figure 6 shows the comparison of the stacked high-pass filtered radial $\tau$
profile on the unWISE (orange) and CMASS (red) locations in our simulations
where the $\tau$ profile has been reconstructed using the TI esimator with an
ACT-like experiment. The input profile used when painting each galaxy in both
maps is shown for comparison in black. The average $\tau$ profile stacked on
CMASS galaxy locations matches the profile painted on an individual galaxy,
i.e., the CMASS sample is sparse enough for the two-halo term to have
negligible contribution. However, after stacking on the unWISE sample, there
is a slight excess in the amplitude of the profile, which reflects the two-
halo term contribution. This shows that the clustering term already starts
becoming important for the unWISE-like galaxy sample and will be even more
important for an LSST-like sample. A tentative detection of this two-halo term
might be possible using ACT$\times$unWISE data.

Figure 5: Ratio of the SNRs when accounting for galaxy shot noise only versus
galaxy shot noise and clustering. For all three galaxy samples, the difference
in SNR ratio between CMB experiments (same colors as Fig. 4) is relatively
small. However, the impact of clustering on SNR increases significantly for
higher density galaxy samples, especially LSST. When neglecting clustering in
the LSST sample, the SNR is overestimated by $\sim 60\%$ for TI (solid lines)
and $\sim 150\%$ for QE (dashed lines).  Figure 6: Comparison of the high-pass
filtered radial $\langle\tau\rangle$ profile when stacking on the unWISE (red)
and CMASS (orange) $\tau$ maps with TI estimator with ACT-like experiment. The
input profile used when painting each galaxy in both maps is shown for
comparison (black). Note the shape of the profiles dips below zero due to the
high-pass filtering applied. The $\tau$ profile stacked on CMASS galaxy
locations matches the profile painted on an individual galaxy, i.e., the CMASS
sample is sparse enough for the two-halo term to have negligible contribution.
However, after stacking on the unWISE sample, the profile amplitude shows a
slight excess, which reflects the two-halo term contribution. The error bars
on the CMASS and unWISE profiles correspond to the expected noise from ACT. A
tentative detection of this two-halo term may be possible with
ACT$\times$unWISE.

###  VI.3 Results

The forecast SNRs for TI, QE, and lens-hardened QE (LHQE), using realistic
simulations and accounting for clustering effects as described above, are
summarized in Tab. 4 and visualized in Fig. 7. Unlike for the clustering
study, we limit $\ell_{\rm max}$ to 3000 for the QE and LHQE to mitigate
foreground contamination. As TI is robust to foregrounds, we use scales up to
$\ell_{\rm max}=5950$.333This rather specific value corresponds to the most
conservative $\ell_{\rm max}$ of the input ILC power spectra used in the
simulations. As expected, we find that the LHQE degrades the SNR to some
extent ($\sim 10-15\%$) when compared to the standard QE as the noise in the
reconstruction of $\tau$ for a LHQE goes up. As this is the forecast SNR for
cross-correlation of $\tau$ and galaxies $C_{\ell}^{\tau g}$ rather than auto-
power spectrum of $\tau$ $C_{\ell}^{\tau\tau}$, effects of the noise penalty
for hardening against lensing are slightly less dramatic. Factors contributing
to the noise (denominator of Eq. (88)) are total galaxy and $\tau$ power
spectra. As expected, we find that for a given galaxy sample, as we go towards
more sensitive CMB experiments, effects of lens-hardening are felt slightly
more as the noise penalty is more for more sensitive experiments Roy_2023 . As
a result, SNR is affected more by lens-hardening for a CMB-S4-like experiment
than a ACT-like experiment. In the following discussion, we limit our
comparisons to TI and QE.

For current experiments like ACT and SPT in combination with CMASS and unWISE
galaxies, TI and QE give similar results. Due to a fairly high shot noise, the
cross-correlation with CMASS does not yield a significant SNR for any CMB
experiment, whereas the signal is detectable using the unWISE sample. We note
this result is consistent with Coulton_2023 , which uses ACT$\times$unWISE
data for the first detection of this effect.

For future CMB experiments, the TI estimator can attain higher SNR than the
QE, most dramatically when using an LSST-like galaxy sample (see Tab. 4 and
Fig. 7). This is due to the use of smaller scales of the CMB. Future CMB
experiments, particularly CMB-S4, have exquisite projected sensitivity and
beam size. With lower noise and a smaller beam, the TI estimator can extract
information from the small scales that QE cannot reliably use due to
foreground biases. Thus, the SNR values across all cross-correlations with
LSST are improved for TI compared to the QE.

| Galaxy
---
sample
| Estimator
---
| SPT
---
| ACT
---
| SO
---
| CMB-S4
---
| CMASS
---
| TI5950
---
QE3000
LHQE3000
| 0.37$\sigma$
---
0.34$\sigma$
0.30$\sigma$
| 0.65$\sigma$
---
0.87$\sigma$
0.78$\sigma$
| 0.92$\sigma$
---
0.96$\sigma$
0.84$\sigma$
| 1.7$\sigma$
---
1.1$\sigma$
0.94$\sigma$
| unWISE
---
| TI5950
---
QE3000
LHQE3000
| 3.0$\sigma$
---
2.9$\sigma$
2.6$\sigma$
| 5.2$\sigma$
---
7.3$\sigma$
6.6$\sigma$
| 7.5$\sigma$
---
8.1$\sigma$
7.2$\sigma$
| 13$\sigma$
---
9.2$\sigma$
8.0$\sigma$
| LSST
---
| TI5950
---
QE3000
LHQE3000
| 11$\sigma$
---
6.6$\sigma$
5.9$\sigma$
| 18$\sigma$
---
16$\sigma$
15$\sigma$
| 25$\sigma$
---
18$\sigma$
16$\sigma$
| 49$\sigma$
---
21$\sigma$
18$\sigma$
Table 4: Summary of the detection SNRs for TI and QE accounting for galaxy
clustering and using realistic $\ell$ ranges for each estimator. We do not
include cosmic variance on the $\tau g$ signal, which is negligible for all
experiments considered. For the standard QE and lens-hardened QE (LHQE), we
use $60<\ell<3000$ where these estimators have limited response to foreground
contamination. For TI, we use the same $\ell$ range as in Tab. 3. TI can use
higher $\ell$ modes as it is robust to foregrounds; the SNR for TI can thus
exceed QE because of this added information. All SNR values are visualized in
Fig. 7.  Figure 7: Visualization of the SNR values given in Tab. 4. The LHQE
(dotted lines) results in a SNR reduction compared to the standard QE (dashed
lines). For current experiments, the QE and TI (solid lines) give similar
SNRs. Cross-correlation with CMASS (orange circles) does not yield a
significant measurement for any CMB experiment, whereas the signal is
detectable using the unWISE sample (red diamonds). For future experiments, TI
can attain higher SNR than the QE, most dramatically when using the LSST
galaxy sample (blue pluses), because it can use smaller CMB scales.

##  VII Conclusions

In this paper, we have derived and studied new estimators for the patchy
screening of the CMB, called TI and signed estimators. These estimators are
modifications of the standard screening QE, analogous to the gradient
inversion estimators for CMB lensing.

The TI and signed estimators make use of a unique feature of patchy screening:
it is the only effect which changes sign based on the local large-scale
temperature. This allows TI and signed estimators to be automatically unbiased
to extragalactic foregrounds like the CIB, tSZ and kSZ, unlike the standard
QE. This foreground robustness can be further enhanced for the signed
estimator via “sign balancing” and “sign thresholding,” which further immunise
it against foregrounds at a negligible noise cost. For the QE, the
corresponding biases can be large and likely require some form of hardening at
some cost in SNR.

The TI and signed estimators are also robust to CMB lensing bias, unlike the
standard QE. Lens hardening remedies this issue for QE, again, with higher
noise penalty.

We also clarify the relationship between the TI estimator, signed estimator
and QE. We show that the QE reconstructs both small- and large-scale screening
information and present simple intuition about these two regimes. We show how
TI is a subtle modification to the small-scale-only QE, which in principle
avoids the noise floor above. Finally, the signed estimator is a modification
to TI which allows even more robustness to foregrounds (via sign balancing and
sign thresholding), at the cost of a factor of $\sqrt{\pi/2}\approx 1.25$
reduction in SNR. Even in the absence of detector noise, foregrounds, and
lensing, the $\tau$ QE has an irreducible source of variance and hits a noise
floor, whereas the TI and signed estimators could perform arbitrarily better
under such noiseless conditions. However, this noise floor is not reached in
practice. For a fixed set of Fourier scales included in the analysis, the
statistical SNR in the TI or signed estimators are similar to the QE but do
not outperform it. However, by being robust to lensing and foregrounds, TI and
signed estimators can be applied to more Fourier modes (pushing to smaller
scales), thus potentially increasing the overall SNR.

We perform realistic forecasts for current and upcoming CMB and LSS
experiments (see Tab. 4 and Fig. 7). Our results highlight the large noise
enhancement due to the clustering of galaxies once the galaxy number density
reaches LSST-like values. Further, our SNR forecast for ACT$\times$unWISE is
consistent with the measurement found in the companion paper Coulton_2023 ,
which is the first detection of this signal.

We have focused on measuring the patchy screening from the CMB temperature
only. The TI and signed estimators should work identically when substituting
temperature T with Stokes Q or U polarization maps. We leave this exploration
to future work.

In the future, patchy screening measurements around galaxies will be a
powerful probe of the gas density profile around them, highly complementary
with kSZ and tSZ effects. Joint analyses of these signals around the same
galaxies will yield valuable insights into the thermodynamical properties of
baryons within these halos and shed more light on the lesser understood
properties of feedback processes. This will have a direct impact on weak
galaxy lensing studies where baryonic effects are a major systematic. They
will also help in breaking the velocity-optical depth degeneracy in kSZ, thus
expanding its power as a probe of the large-scale velocities to constrain the
growth rate of structure, modifications to general relativity and local
primordial non-Gaussianity.

## Data Availability

The LensQuEst and ThumbStack software packages are publicly available under an
open source license.

## Acknowledgments

We thank Nicholas Battaglia, Simone Ferraro, Darby Kramer, Srinivasan
Raghunathan, Rashid Sunyaev, Alexander van Engelen and Ziang Yan for useful
discussions. We thank Joanna Dunkley, Lyman Page and Edward Wollack for
helpful comments on the manuscript.

This material is based upon work supported by the National Science Foundation
Graduate Research Fellowship under Grant No. DGE-2146755. TS thanks the LSSTC
Data Science Fellowship Program, which is funded by LSSTC, NSF Cybertraining
Grant #1829740, the Brinson Foundation, and the Moore Foundation; their
participation in the program has benefited this work. This work received
support from the U.S. Department of Energy under contract number DE-
AC02-76SF00515 to SLAC National Accelerator Laboratory. Some of the computing
for this project was performed on the Sherlock cluster. We thank Stanford
University and the Stanford Research Computing Center for providing
computational resources and support.

## References

  * (1) N. Aghanim, S. Majumdar, and J. Silk, Reports on Progress in Physics 71, 066902 (2008).
  * (2) Y. B. Zeldovich and R. A. Sunyaev, Ap&SS 4, 301 (1969).
  * (3) R. A. Sunyaev and Y. B. Zeldovich, Comments on Astrophysics and Space Physics 4, 173 (1972).
  * (4) R. A. Sunyaev and I. B. Zeldovich, MNRAS 190, 413 (1980).
  * (5) L. Salvati, M. Douspis, and N. Aghanim, A&A 614, A13 (2018).
  * (6) B. Bolliet, B. Comis, E. Komatsu, and J. F. Macías-Pérez, MNRAS 477, 4957 (2018).
  * (7) B. Soergel et al., Mon. Not. Roy. Astron. Soc. 478, 5320 (2018).
  * (8) A. Amon and G. Efstathiou, MNRAS 516, 5355 (2022).
  * (9) G. Aricò et al., A&A 678, A109 (2023).
  * (10) E. Schaan et al., Phys. Rev. D 93, 082002 (2016).
  * (11) E. Schaan et al., Phys. Rev. D 103, 063513 (2021).
  * (12) S. Amodeo et al., Phys. Rev. D 103, 063514 (2021).
  * (13) E. Schiappucci et al., Phys. Rev. D 107, 042004 (2023).
  * (14) J. P. Ostriker and E. T. Vishniac, ApJ 306, L51 (1986).
  * (15) E. T. Vishniac, Astrophys. J. 322, 597 (1987).
  * (16) C. Dvorkin and K. M. Smith, Phys. Rev. D 79, 043003 (2009).
  * (17) C. Feng and G. Holder, Phys. Rev. D 97, 123523 (2018).
  * (18) O. Doré, J. F. Hennawi, and D. N. Spergel, Astrophys. J. 606, 46 (2004).
  * (19) J. C. Hill et al., Phys. Rev. Lett. 117, 051301 (2016).
  * (20) A. Kusiak et al., Phys. Rev. D 104, 043518 (2021).
  * (21) E.-M. Mueller, F. de Bernardis, R. Bean, and M. D. Niemack, Phys. Rev. D 92, 063501 (2015).
  * (22) N. S. Sugiyama, T. Okumura, and D. N. Spergel, J. Cosmology Astropart. Phys 2017, 057 (2017).
  * (23) E.-M. Mueller, F. de Bernardis, R. Bean, and M. D. Niemack, Astrophys. J. 808, 47 (2015).
  * (24) M. Münchmeyer et al., Phys. Rev. D 100, 083508 (2019).
  * (25) W. Hu, Astrophys. J. 529, 12 (2000).
  * (26) A. Natarajan et al., ApJ 776, 82 (2013).
  * (27) A. Roy, A. Lapi, D. Spergel, and C. Baccigalupi, J. Cosmology Astropart. Phys 2018, 014 (2018).
  * (28) B. Horowitz, S. Ferraro, and B. D. Sherwin, MNRAS 485, 3919 (2019).
  * (29) B. Hadzhiyska, B. D. Sherwin, M. Madhavacheril, and S. Ferraro, Phys. Rev. D 100, 023547 (2019).
  * (30) W. Hu and T. Okamoto, Astrophys. J. 574, 566 (2002).
  * (31) M. Su et al., arXiv e-prints arXiv:1106.4313 (2011).
  * (32) R. A. Sunyaev and I. B. Zeldovich, ARA&A 18, 537 (1980).
  * (33) S. Y. Sazonov and R. A. Sunyaev, MNRAS 310, 765 (1999).
  * (34) M. Shehzad Emritte, S. Colafrancesco, and P. Marchegiani, J. Cosmology Astropart. Phys 2016, 031 (2016).
  * (35) Planck Collaboration et al., A&A 641, A6 (2020).
  * (36) W. R. Coulton et al., (in prep.).
  * (37) T. Namikawa et al., Phys. Rev. D 104, 063514 (2021).
  * (38) T. Namikawa, D. Hanson, and R. Takahashi, MNRAS 431, 609 (2013).
  * (39) A. Roy, A. van Engelen, V. Gluscevic, and N. Battaglia, ApJ 951, 50 (2023).
  * (40) A. van Engelen et al., ApJ 786, 13 (2014).
  * (41) S. J. Osborne, D. Hanson, and O. Doré, J. Cosmology Astropart. Phys 2014, 024 (2014).
  * (42) N. MacCrann et al., arXiv e-prints arXiv:2304.05196 (2023).
  * (43) E. Schaan and S. Ferraro, Phys. Rev. Lett. 122, 181301 (2019).
  * (44) N. Sailer, E. Schaan, and S. Ferraro, Phys. Rev. D 102, 063517 (2020).
  * (45) B. A. Benson et al., in Millimeter, Submillimeter, and Far-Infrared Detectors and Instrumentation for Astronomy VII, edited by W. S. Holland and J. Zmuidzinas (SPIE, Bellingham, Washington, 2014).
  * (46) A. N. Bender et al., in Millimeter, Submillimeter, and Far-Infrared Detectors and Instrumentation for Astronomy IX, Vol. 10708 of Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, edited by J. Zmuidzinas and J.-R. Gao (SPIE, Bellingham, Washington, 2018), p. 1070803.
  * (47) J. A. Sobrin et al., ApJS 258, 42 (2022).
  * (48) S. W. Henderson et al., Journal of Low Temperature Physics 184, 772 (2016).
  * (49) Simons Observatory Collaboration et al., J. Cosmology Astropart. Phys 2019, 056 (2019).
  * (50) K. Abazajian et al., arXiv e-prints arXiv:1907.04473 (2019).
  * (51) K. S. Dawson et al., AJ 145, 10 (2013).
  * (52) B. Reid et al., MNRAS 455, 1553 (2016).
  * (53) E. L. Wright et al., AJ 140, 1868 (2010).
  * (54) A. M. Meisner, D. Lang, and D. J. Schlegel, Research Notes of the American Astronomical Society 2, 1 (2018).
  * (55) LSST Dark Energy Science Collaboration (LSST DESC) et al., ApJS 253, 31 (2021).
  * (56) J. Hartlap, P. Simon, and P. Schneider, A&A 464, 399 (2007).
  * (57) A. Maniyar, M. Béthermin, and G. Lagache, A&A 645, A40 (2021).
  * (58) W. R. Coulton et al., arXiv e-prints arXiv:2307.01258 (2023).
  * (59) S. Raghunathan and Y. Omori, ApJ 954, 83 (2023).
  * (60) A. Krolewski, S. Ferraro, E. F. Schlafly, and M. White, J. Cosmology Astropart. Phys 2020, 047 (2020).
  * (61) A. Krolewski, S. Ferraro, and M. White, J. Cosmology Astropart. Phys 2021, 028 (2021).

##  Appendix A Optimal $\delta T_{L}$ estimator in the presence of foreground
bias

In this section, we derive estimators of the large-scale CMB temperature
around a given galaxy. To declutter the derivation, we simplify our notation
as:

$d_{\bm{\ell}}=T_{\bm{\ell}}+n_{\bm{\ell}}+Au_{\bm{\ell}},$ (35)

where $d_{\bm{\ell}}$ is the temperature data as a function of Fourier
multipole ${\bm{\ell}}$, $T_{\bm{\ell}}$ is the true CMB (with known power
spectrum $C_{\ell}$) to be estimated, $n_{\bm{\ell}}$ is the noise (with known
power spectrum $N_{\ell}$), and $Au_{\bm{\ell}}$ is a foreground emission,
with a known profile $u_{\bm{\ell}}$ (e.g., the beam for a point source) and
unknown amplitude $A$. We would like to build an estimator
$\hat{T}_{\bm{\ell}}$ of $T_{\bm{\ell}}$, which minimizes the mean squared
error $\langle\left(\hat{T}_{\bm{\ell}}-T_{\bm{\ell}}\right)^{2}\rangle$. We
restrict ourselves to estimators linear in the data, i.e. we wish to determine
the weights $W_{\bm{L}}$ for each ${\bm{\ell}}$ such that:

$\hat{T}_{\bm{\ell}}=\sum_{\bm{L}}W_{\bm{L}}d_{{\bm{\ell}}-{\bm{L}}}.$ (36)

In the absence of foregrounds, i.e., if $A$ is known to be zero, minimizing
the mean-squared error leads to the well-known Wiener filter:

$\hat{T}_{\bm{\ell}}=\frac{C_{\ell}}{C_{\ell}+N_{\ell}}d_{\bm{\ell}}.$ (37)

The Wiener filter can also be derived as the maximum a posteriori estimator
for $T_{\bm{\ell}}$, in the absence of foregrounds ($A=0$), as well as in the
presence of foregrounds for which the prior on the amplitude $A$ is flat. Such
foregrounds would have a random sign, which is typically not the case in our
problem.

We now consider the case of non-zero foreground and impose the additional
constraint that the estimator $\hat{T}_{\bm{\ell}}$ should have zero response
to the foreground amplitude:

$\frac{d\hat{T}_{\bm{\ell}}}{dA}=0.$ (38)

In what follows, we show that the linear estimator which minimizes the mean
squared error, while having zero response to the foreground emission, is
simply obtained by Wiener-filtering the data, after having subtracted the
matched-filtered foreground template.

We introduce a Lagrange multiplier $\lambda$, to minimize the mean-squared
error under the constraint

$\mathcal{L}=\langle\left(\hat{T}_{\bm{\ell}}-T_{\bm{\ell}}\right)^{2}\rangle-\lambda\frac{d\hat{T}_{\bm{\ell}}}{dA}.$
(39)

In terms of the weights $W_{\bm{L}}$, this loss function becomes

$\displaystyle\mathcal{L}$
$\displaystyle=\left(W_{0}-1\right)^{2}C_{\ell}+W_{0}^{2}N_{\ell}$ (40)
$\displaystyle+\sum_{{\bm{L}}\neq
0}W_{\bm{L}}^{2}\left(C_{{\bm{\ell}}-{\bm{L}}}+N_{{\bm{\ell}}-{\bm{L}}}\right)$
$\displaystyle+A^{2}\sum_{{\bm{L}},{\bm{L}}^{\prime}}W_{\bm{L}}W_{{\bm{L}}^{\prime}}u_{{\bm{\ell}}-{\bm{L}}}u_{{\bm{\ell}}-{\bm{L}}^{\prime}}$
$\displaystyle-\lambda\sum_{{\bm{L}}}W_{\bm{L}}u_{{\bm{\ell}}-{\bm{L}}}.$

Setting $\partial\mathcal{L}/\partial W_{\bm{L}}$ to zero gives

$W_{\bm{L}}=\frac{C_{\ell}}{C_{\ell}+N_{\ell}}\delta^{K}_{{\bm{L}},0}-\left(\alpha+\lambda\right)\frac{u_{{\bm{\ell}}-{\bm{L}}}}{C_{{\bm{\ell}}-{\bm{L}}}+N_{{\bm{\ell}}-{\bm{L}}}},$
(41)

where we have defined

$\alpha\equiv
A\sum_{{\bm{L}}^{\prime}}W_{{\bm{L}}^{\prime}}u_{{\bm{\ell}}-{\bm{L}}^{\prime}}.$
(42)

Setting $\partial\mathcal{L}/\partial\lambda$ to zero gives
$\sum_{{\bm{L}}}W_{{\bm{L}}}u_{{\bm{\ell}}-{\bm{L}}}=0$, which means that
$\alpha=0$ and lets us solve for $\lambda$ in Eq. (40):

$\lambda=\frac{C_{\ell}}{C_{\ell}+N_{\ell}}\
\frac{u_{\bm{\ell}}}{\left(\sum_{\bm{L}}\frac{u_{\bm{L}}^{2}}{C_{L}+N_{L}}\right)}.$
(43)

Plugging this into Eq. (40) yields the expressions for the weights:

$W_{\bm{L}}=\frac{C_{\ell}}{C_{\ell}+N_{\ell}}\left[\delta^{K}_{{\bm{L}},0}-\frac{\frac{u_{{\bm{\ell}}-{\bm{L}}}}{C_{{\bm{\ell}}-{\bm{L}}}+N_{{\bm{\ell}}-{\bm{L}}}}}{\left(\sum_{{\bm{L}}^{\prime}}\frac{u_{{\bm{L}}^{\prime
2}}}{C_{L^{\prime}}+N_{L^{\prime}}}\right)}u_{\bm{\ell}}\right].$ (44)

In conclusion, the estimator for $T_{\bm{\ell}}$, linear in the data
$d_{\bm{\ell}}$, which minimizes the mean squared error while having zero
response to the foreground emission is

$\hat{T}_{\bm{\ell}}=\frac{C_{\ell}}{C_{\ell}+N_{\ell}}\left[d_{\bm{\ell}}-\hat{A}u_{\bm{\ell}}\right],$
(45)

where $\hat{A}$ is the matched filter estimator for the foreground amplitude
$A$:

$\hat{A}\equiv\frac{\left(\sum_{{\bm{L}}}\frac{d_{{\bm{L}}}u_{{\bm{L}}}}{C_{L}+N_{L}}\right)}{\left(\sum_{{\bm{L}}^{\prime}}\frac{u_{{\bm{L}}^{\prime
2}}}{C_{L^{\prime}}+N_{L^{\prime}}}\right)}.$ (46)

As stated, this reduces to subtracting the foreground template
$\hat{A}u_{\bm{\ell}}$ from the data, where $\hat{A}$ is the matched filter
for the profile $u_{\bm{\ell}}$, and then Wiener-filtering the result.

##  Appendix B Variance comparison: small-scale QE vs. TI

We wish to get intuition about the noise properties of the QE and TI, in order
to understand in which case one may be a better estimator than the other.

We decompose the observed small and large-scale temperatures around a given
galaxy into signal and noise:

$\left\\{\begin{aligned} &\delta\hat{T}_{S}=-\delta T_{L}\tau+n_{S}\\\
&\delta\hat{T}_{L}=\delta T_{L}+n_{L}\\\ \end{aligned}\right..$ (47)

Here, the noise $n_{S}$ includes everything in the small-scale map other than
our screening effect, i.e. detector noise, foregrounds, lensing and the small,
Silk-damped primary CMB. The large-scale noise $n_{L}$ should be thought of as
very small, since the primary CMB is very well-measured on degree scales.

To get a sense of the relevant noise regimes, we introduce the variances
$\sigma^{2}_{T_{L}},\sigma^{2}_{\tau},\sigma^{2}_{n_{L}}$ and
$\sigma^{2}_{n_{S}}$ of $\delta\hat{T}_{L}$, the true $\tau$, $n_{L}$ and
$n_{S}$. In what follows, we assume $\sigma^{2}_{n_{L}}\ll\sigma^{2}_{T_{L}}$,
realistic for current and upcoming experiments, and that $\delta
T_{L},\tau,n_{L}$ and $n_{S}$ are statistically independent.

We can now express the TI and QE estimators in real space as:

$\left\\{\begin{aligned}
&\hat{\tau}^{\text{TI}}=-\delta\hat{T}_{S}/\delta\hat{T}_{L}\\\
&\hat{\tau}^{\text{QE}}=-\delta\hat{T}_{S}\delta\hat{T}_{L}/\sigma^{2}_{T_{L}}\\\
\end{aligned}\right.$ (48)

We then plug in the expressions for the observed quantities Eq. (47), in terms
of the true ones, into the QE and TI estimators. As expected, we obtain a
linear response in the true $\tau$ and a noise term:

$\left\\{\begin{aligned}
&\hat{\tau}^{\text{TI}}=\alpha^{\text{TI}}\tau+n^{\text{TI}}\\\
&\hat{\tau}^{\text{QE}}=\alpha^{\text{QE}}\tau+n^{\text{QE}}\\\
\end{aligned}\right..$ (49)

###  B.1 Irreducible noise in the QE response

For the TI estimator,

$\alpha^{\text{TI}}=\frac{\delta T_{L}}{\delta
T_{L}+n_{L}}=1-\frac{n_{L}}{\delta T_{L}}+\left(\frac{n_{L}}{\delta
T_{L}}\right)^{2}+...$ (50)

such that the mean response $\langle\alpha^{\text{TI}}\rangle$ differs from
unity by a bias
$\mathcal{O}\left(\sigma^{2}_{n_{L}}/\sigma^{2}_{T_{L}}\right)$. The variance
of the response $\alpha^{\text{TI}}$ is similarly small, with

$\text{var}\left(\alpha^{\text{TI}}\right)=\mathcal{O}\left(\sigma^{2}_{n_{L}}/\sigma^{2}_{T_{L}}\right).$
(51)

Thus, in practice, the bias and the variance of the TI response are negligibly
small, i.e. $\alpha^{\text{TI}}\simeq 1$.

In contrast, the linear response of the QE estimator is

$\alpha^{\text{QE}}=\frac{\delta T_{L}\left(\delta
T_{L}+n_{L}\right)}{\sigma^{2}_{T_{L}}}.$ (52)

This response is exactly unbiased, i.e. $\langle\alpha^{\text{QE}}\rangle=1$.
However, its variance does not go to zero in the limit of low noise:

$\text{var}\left(\alpha^{\text{TI}}\right)=2+\sigma^{2}_{n_{L}}/\sigma^{2}_{T_{L}}.$
(53)

Thus, if the noises $n^{\text{TI}}$ and $n^{\text{QE}}$ are negligible, the TI
estimator converges to a perfect reconstruction of the true $\tau$, whereas
the QE estimator converges to $\alpha^{\text{QE}}\tau$, where
$\alpha^{\text{QE}}$ is a random number with variance 2.

###  B.2 Comparable estimator noises

The noise $n^{\text{TI}}$ of the TI and QE estimators are

$\left\\{\begin{aligned} &n^{\text{TI}}=\frac{n_{S}}{\delta T_{L}+n_{L}}\\\
&n^{\text{QE}}=\frac{n_{S}\left(\delta
T_{L}+n+L\right)}{\sigma^{2}_{T_{L}}}\\\ \end{aligned}\right..$ (54)

Since $\langle 1/\delta T_{L}^{2}\rangle\gtrsim 1/\sigma^{2}_{T_{L}}$ by
convexity of $x\rightarrow 1/x$, we obtain

$\text{var}\left(n^{\text{TI}}\right)\gtrsim\frac{\sigma^{2}_{n_{S}}}{\sigma^{2}_{T_{L}}}\left(1+\frac{\sigma^{2}_{n_{L}}}{\sigma^{2}_{T_{L}}}\right)=\text{var}\left(n^{\text{QE}}\right).$
(55)

Thus, while the TI noise is larger than the QE noise, they are comparable.

###  B.3 Total variance: TI outperforms QE in the low-noise regime

To decide in which regime one estimator might be preferred over the other, we
take into account all the contributions to the estimator variance, including
the true signal $\tau$, the noise $n^{\text{TI/QE}}$ and the response
$\alpha^{\text{TI/QE}}$. This leads to

$\text{var}\left(\hat{\tau}^{\text{TI/QE}}\right)=\text{var}\left(\alpha^{\text{TI/QE}}\right)\sigma^{2}_{\tau}+\text{var}\left(n^{\text{TI/QE}}\right),$
(56)

such that we find:

$\left\\{\begin{aligned}
&\text{var}\left(\hat{\tau}^{\text{TI}}\right)=\frac{\sigma^{2}_{n_{L}}}{\sigma^{2}_{T_{L}}}\
\sigma^{2}_{\tau}+\mathcal{O}\left(\frac{\sigma^{2}_{n_{S}}}{\sigma^{2}_{T_{L}}}\right)\\\
&\text{var}\left(\hat{\tau}^{\text{QE}}\right)=\left(2+\frac{\sigma^{2}_{n_{L}}}{\sigma^{2}_{T_{L}}}\right)\
\sigma^{2}_{\tau}+\frac{\sigma^{2}_{n_{S}}}{\sigma^{2}_{T_{L}}}\\\
\end{aligned}\right.$ (57)

It thus appears that the TI is preferred over the QE when the “2” term is
relevant, i.e. when:

$\sigma^{2}_{n_{L}}\lesssim\sigma^{2}_{T_{L}}\text{ and
}\sigma^{2}_{n_{S}}\lesssim\sigma^{2}_{T_{L}}\sigma^{2}_{\tau}.$ (58)

The first condition is easily satisfied, since the large scale temperature is
measured at high SNR in current experiments. The second condition effectively
requires that the $n_{S}$ be smaller than the screening effect $\delta
T_{L}\tau$ in the measured small-scale temperature. However, the noise $n_{S}$
includes not only detector noise, which could in principle be made arbitrarily
small, but also foregrounds and lensing. Thus, going beyond QE with TI
requires significant foreground cleaning. Lensing is much larger than
screening in the small-scale temperature power spectrum, which violates
$\sigma^{2}_{n_{S}}\lesssim\sigma^{2}_{T_{L}}\sigma^{2}_{\tau}$. Thus, the TI
SNR should only outperform the QE’s if a large amount of delensing can be
achieved. We leave quantifying this required amount of delensing to future
work.

##  Appendix C Two versions of the approximate small-scale QE; same comparison
with TI

Here we give the expressions for the stacked TI and QE estimators. These show
that the key difference between TI and QE remains true for the stacked
estimators: the QE’s response to the true $\tau$ has an irreducible variance,
even in the absence of any noise.

As we show in the main text, the individual TI estimator for galaxy $i$ is
$\tau^{\text{TI}}_{i}\equiv\delta\hat{T}_{Si}/\delta\hat{T}_{Li}$, with
variance $\sigma^{2}_{n_{s}}/\delta\hat{T}_{Li}^{2}$, such that the inverse-
variance weighted average over all galaxies is given by Eq. (24):

$\hat{\tau}^{\text{TI stack}}=\frac{\sum_{i}\hat{\delta T}_{Si}\hat{\delta
T_{Li}}/\sigma_{n_{S}}^{2}}{\sum_{j}\hat{\delta
T_{L}}_{j}^{2}/\sigma_{n_{S}}^{2}}.$ (59)

Ignoring the small large-scale noise $n_{L}$, but taking into account the
small-scale noise $n_{S}$, this can be rewritten as

$\hat{\tau}^{\text{TI stack}}=\underbrace{1}_{\text{TI
response}}\times\tau+\underbrace{\frac{\sum_{i}n_{Si}\hat{\delta
T_{Li}}/\sigma_{n_{S}}^{2}}{\sum_{j}\hat{\delta
T_{L}}_{j}^{2}/\sigma_{n_{S}}^{2}}}_{\text{TI noise}},$ (60)

i.e. the response to the true $\tau$ is unity (again, for small $n_{L}$).

On the other hand, the small-scale QE estimator for a given galaxy $i$ takes
the form
$\tau^{\text{QE}}_{i}\equiv\delta\hat{T}_{Si}\delta\hat{T}_{Li}/\sigma^{2}_{T_{L}}$.
When quantifying the variance of $\tau^{\text{QE}}_{i}$, two options are
possible. A “local” variance estimates takes into account the local value of
$\delta\hat{T}_{Li}$, whereas a “global” variance estimates marginalizes over
$\delta\hat{T}_{Li}$, given its variance $\sigma^{2}_{T_{L}}$.

The global variance estimate $\sigma^{2}_{n_{S}}/\sigma^{2}_{T_{L}}$ is in the
spirit of the standard QE, as it preserves the quadratic nature of the
estimator, by not including the data itself in the weighting. In this case,
the stacked estimator is then:

$\displaystyle\hat{\tau}^{\text{QE ``global'' stack}}$
$\displaystyle=\frac{\sum_{i}\delta\hat{T}_{Si}\delta\hat{T}_{Li}/\sigma^{2}_{n_{S}}}{\sum_{j}\sigma^{2}_{T_{L}}/\sigma_{n_{S}}^{2}}$
(61)
$\displaystyle=\underbrace{\frac{\sum_{i}\delta\hat{T}_{Li}^{2}/\sigma^{2}_{n_{S}}}{\sum_{j}\sigma^{2}_{T_{L}}/\sigma_{n_{S}}^{2}}}_{\text{response,
var}\neq
0}\times\tau+\underbrace{\frac{\sum_{i}n_{Si}\delta\hat{T}_{Li}/\sigma^{2}_{n_{S}}}{\sum_{j}\sigma^{2}_{T_{L}}/\sigma_{n_{S}}^{2}}}_{\text{noise}}.$

If one were to instead weigh by the local variance
$\sigma^{2}_{n_{S}}/\sigma^{4}_{T_{L}}\delta\hat{T}_{Li}$, one obtains the
“local” stacked QE, no longer quadratic in the data:

$\displaystyle\hat{\tau}^{\text{QE ``local'' stack}}$
$\displaystyle=\frac{\sum_{i}\frac{\delta\hat{T}_{Si}}{\delta\hat{T}_{Li}}\frac{\sigma^{2}_{T_{L}}}{\sigma^{2}_{n_{S}}}}{\sum_{j}\frac{1}{\delta\hat{T}^{2}_{Lj}}\frac{\sigma^{4}_{T_{L}}}{\sigma^{2}_{n_{S}}}}$
(62)
$\displaystyle=\underbrace{\frac{\sum_{i}\sigma^{2}_{T_{L}}/\sigma^{2}_{n_{S}}}{\sum_{j}\frac{\sigma^{2}_{T_{L}}}{\delta\hat{T}_{Lj}^{2}}\frac{\sigma^{2}_{T_{L}}}{\sigma^{2}_{n_{S}}}}}_{\text{response,
var}\neq
0}\times\tau+\underbrace{\frac{\sum_{i}\frac{n_{Si}}{\delta\hat{T}_{Li}}\frac{\sigma^{2}_{T_{L}}}{\sigma^{2}_{n_{S}}}}{\frac{\sigma^{2}_{T_{L}}}{\delta\hat{T}_{Lj}^{2}}\frac{\sigma^{2}_{T_{L}}}{\sigma^{2}_{n_{S}}}}}_{\text{noise}}.$

Thus, both for the “local” and “global” weightings, the stacked QE estimator
differs from the stacked TI estimator. Furthermore, the response to the true
$\tau$ of the stacked QE estimators still has an irreducible variance, even as
the small-scale noise $\sigma^{2}_{n_{S}}$ goes to zero.

As discussed in the main text, the stacked global QE and stacked TI appear
very similar, differing only in their denominator. The fractional difference
between the denominators scales as $1/\sqrt{N_{\text{galaxies}}}$, such that
it becomes smaller and smaller with larger galaxy samples. However, this
difference is not irrelevant, because the variance of the estimators (stacked
QE and TI) also decreases as $1/\sqrt{N_{\text{galaxies}}}$. Thus, this
difference is actually significant.

##  Appendix D $\tau^{\rm sgn}$ estimator: Sign balancing helps with imperfect
foreground cancellation

“Sign balancing” refers to requiring that the number of positively and
negatively weighted galaxies in the stack is the same. Sign balancing is
effective as long as the contribution to the foreground emission comes from
$>0.1\%$ of the total number of galaxies. This appendix presents the evidence
for this statement.
The signed estimator $\hat{\tau}^{\rm sgn}$ for the foreground emission
becomes

$\hat{\tau}^{\rm
sgn}\equiv\frac{\sum_{i}\operatorname{sgn}(T_{L,i})f_{i}}{\sum_{i}|T_{L,i}|}$
(63)

where $T_{L,i}$ is the large-scale temperature at the position of galaxy $i$,
and $f_{i}$ is the small-scale temperature due to foreground emission from
galaxy $i$. Simplifying notation, we write

$\displaystyle\hat{\tau}^{\rm sgn}$ $\displaystyle=\frac{\frac{1}{n_{\rm
gal}}\sum_{i}\operatorname{sgn}(T_{L,i})f_{i}}{\frac{1}{n_{\rm
gal}}\sum_{i}|T_{L,i}|}$ (64) $\displaystyle\equiv\frac{T^{\rm N}}{T^{\rm
D}},$

where $n_{\rm gal}$ is the number of galaxies in the stack, and $T^{\rm N}$
and $T^{\rm D}$ are defined as the numerator and denominator of the r.h.s.,
respectively. We have the following relations:

$\left\\{\begin{aligned}
\operatorname{sgn}(T_{L,i})&\mathrel{\perp\\!\\!\\!\perp}|T_{L,j}|\\\
\operatorname{sgn}(T_{L,i})&\mathrel{\perp\\!\\!\\!\perp}f_{j}\\\
|T_{L,i}|&\mathrel{\perp\\!\\!\\!\perp}f_{j}\\\
|T_{L,i}|&\mathrel{\,\not\\!\perp\\!\\!\\!\perp}|T_{L,j}|\end{aligned}\right.\forall\,i,j,$
(65)

where $\mathrel{\perp\\!\\!\\!\perp}$ represents statistical independence of
the l.h.s. and r.h.s.. The last relation comes from the fact that galaxies may
be located in a region where the CMB temperature is highly correlated (i.e. in
$\sim$1 deg2 patches). Since we assume
$\operatorname{sgn}(T_{L,i})\mathrel{\perp\\!\\!\\!\perp}f_{j}$, there is no
mean bias from foregrounds, i.e. $\left<\hat{\tau}\right>=0$. However, for a
given realization we may find $\left<\hat{\tau}\right>\neq 0$ due to imperfect
cancellation among the $n_{\rm gal}$ sources. Here we derive the typical value
of $\hat{\tau}_{\rm bias}$, i.e. the RMS error of $\hat{\tau}_{\rm bias}$, to
estimate the bias due to imperfect cancellation. First, we simplify our
expression for $\operatorname{var}(\hat{\tau}_{\rm bias})$:

$\displaystyle\operatorname{var}(\hat{\tau}_{\rm bias})$
$\displaystyle=\operatorname{var}\left(\frac{1}{T^{\rm D}}T^{\rm N}\right)$
(66) $\displaystyle=\operatorname{var}\left(\frac{1}{T^{\rm
D}}\right)\operatorname{var}(T^{\rm N})+\left<\frac{1}{T^{\rm
D}}\right>^{2}\operatorname{var}(T^{\rm N})$
$\displaystyle\qquad+\underbrace{\left<T^{\rm
N}\right>^{2}}_{=0}\operatorname{var}\left(\frac{1}{T^{\rm D}}\right)$
$\displaystyle=\left(\operatorname{var}\left(\frac{1}{T^{\rm
D}}\right)+\left<\frac{1}{T^{\rm
D}}\right>^{2}\right)\operatorname{var}(T^{\rm N})$
$\displaystyle=\left<\frac{1}{(T^{\rm D})^{2}}\right>\operatorname{var}(T^{\rm
N}).$

Since $\operatorname{var}(\hat{\tau}_{\rm
bias})\propto\operatorname{var}(T^{\rm N})$, we can simply use the variance of
the numerator $T^{\rm N}\equiv\frac{1}{n_{\rm
gal}}\sum\operatorname{sgn}(T_{L,i})f_{i}$ for the remainder of the analysis.
We first determine the typical value of $T^{\rm N}$ in the absence of
fluctuations in foreground contamination, i.e.,
$f_{i}=\left<f_{i}\right>=f~{}\forall\,i$. Within a correlated $\sim 1$ deg2
patch of the CMB, sgn($T_{L,i}$) are identical. For $n_{\rm patch}$ patches,
we expect $\sqrt{n_{\rm patch}}$ excess patches of one sign (positive or
negative). Hence, the total excess number of galaxies with positive or
negative weight is

$n_{\rm gal,excess}=\sqrt{n_{\rm patch}}\left(\frac{n_{\rm gal}}{n_{\rm
patch}}\right)=\frac{n_{\rm gal}}{\sqrt{n_{\rm patch}}}.$ (67)

Therefore, the typical value of $T^{\rm N}$ is

$\displaystyle T^{\rm N}_{\rm typical}$ $\displaystyle=\frac{1}{n_{\rm
gal}}\left(n_{\rm gal,excess}f\right)$ (68) $\displaystyle=\frac{1}{n_{\rm
gal}}\left(\frac{n_{\rm gal}f}{\sqrt{n_{\rm patch}}}\right)$
$\displaystyle=\frac{f}{\sqrt{n_{\rm patch}}}.$

By “sign balancing,” i.e., ensuring $n_{\rm gal,excess}=0$ before stacking the
sources, the bias due to imperfect cancellation of sources is completely
negated in the regime where all sources have equal foreground emission.
We can relax this assumption and allow $f_{i}$ to be distinct for all $i$. We
define $\sigma_{f}^{2}\equiv\operatorname{var}(f_{i})$ and assume the source
weights are already sign balanced such that
$\sum_{i}\operatorname{sgn}(T_{L,i})=0$. In this case,

$\displaystyle{\rm var}(T^{\rm N})$ $\displaystyle=\frac{1}{n_{\rm
gal}^{2}}\sum_{i}\operatorname{sgn}(T_{L,i})^{2}\sigma_{f}^{2}$ (69)
$\displaystyle=\frac{\sigma_{f}^{2}}{n_{\rm gal}}.$

Thus,

$T^{\rm N}_{\rm typical}=\frac{\sigma_{f}}{\sqrt{n_{\rm gal}}}.$ (70)

We see sign balancing does not fully negate the imperfect cancellation caused
by foreground fluctuations. It is therefore only effective if the bias due to
mismatched positive/negative weights is much greater than the bias due to
scatter in foreground levels, i.e.,

$\displaystyle\frac{f}{\sqrt{n_{\rm patch}}}$
$\displaystyle\gg\frac{\sigma_{f}}{\sqrt{n_{\rm gal}}}$ (71)
$\displaystyle\Rightarrow\frac{\sigma_{f}}{f}$
$\displaystyle\ll\sqrt{\frac{n_{\rm gal}}{n_{\rm patch}}}$
$\displaystyle\ll\sqrt{\rho_{\rm gal}A_{\rm patch}},$

where $\rho_{\rm gal}$ is the galaxy number density and $A_{\rm patch}$ is the
area of the correlated temperature patch. $A_{\rm patch}\sim 1\,$deg2, and
$\rho_{\rm gal}\sim 1000\,$deg-2 for current galaxy surveys. Eq. 71 is not
satisfied if a small fraction of the galaxy sample dominates the foreground
signal. We use the following toy model to determine the minimum fraction of
foreground emitting galaxies required for sign balancing to be effective.
Consider the case where a fraction $x$ of the galaxy sample emits foregrounds,
and $1-x$ emits no foregrounds. The population of emitting galaxies must
contribute a total of $f/x$ foreground emission such that
$\left<f_{i}\right>=(f/x)*x+(0)*(1-x)=f$. The variance of the foreground
emission is

$\displaystyle{\rm var}(f_{i})$
$\displaystyle=\left<f_{i}^{2}\right>-\left<f_{i}\right>^{2}$ (72)
$\displaystyle=\left[\left(\frac{f}{x}\right)^{2}x+0^{2}(1-x)\right]-f^{2}$
$\displaystyle=f^{2}\left(\frac{1}{x}-1\right).$

Thus,

$\frac{\sigma_{f}}{f}=\sqrt{\frac{1}{x}-1}.$ (73)

Substituting this result into Eq. 71 and squaring both sides, we determine

$\displaystyle\frac{1}{x}-1$ $\displaystyle\ll\rho_{\rm gal}A_{\rm patch}$
(74) $\displaystyle\ll 1000$ $\displaystyle\Rightarrow x$ $\displaystyle\gg
10^{-3}.$

Thus, sign balancing is effective as long as the foreground emission comes
from $>0.1\%$ of the total galaxies. To ensure this condition is satisfied,
one can exclude clusters or high-mass groups from the galaxy sample, which
would dominate the foreground contamination.

##  Appendix E Mean foreground contamination comparison for stacked TI and
signed estimators

There are three regimes of foreground contamination for the TI and signed
estimators: (1) present small-scale foreground $f_{S}$ but negligible large-
scale foreground $f_{L}$, (2) present large-scale foreground $f_{L}$ but
negligible small-scale foreground $f_{S}$, and (3) non-negligible correlated
foregrounds $f_{S}$ and $f_{L}$. In the first two cases, $f_{S}$ (or $f_{L}$)
is uncorrelated with the temperature fluctuation $\delta T_{L}$ (or $\delta
T_{S}$) and thus introduces no bias to either the TI or signed estimators. We
study the third regime for both estimators below.
The stacked TI estimator including both large- and small-scale foreground
contamination is

$\hat{\tau}^{\rm TI}=-\frac{\sum(\delta T_{L}+f_{L})(\delta
T_{S}+f_{S})}{\sum(\delta T_{L}+f_{L})^{2}}.$ (75)

We assume $f_{L}\ll\sigma_{T_{L}}$ (for the typical values of
$\sigma_{T_{L}}\sim 100\,\mu$K and $f_{L}\sim 1\,\mu$K this is quite
accurate). For simplicity, we also assume all objects have equal $f_{L}$ and
$f_{S}$. The bias to the stacked TI estimator as a function of radial distance
$r$ in the 1D stacked photometry profile to first order in
$f_{L}/\sigma_{T_{L}}$ is

$\displaystyle\hat{\tau}_{\rm bias}^{\rm TI}(r)$
$\displaystyle=-\left<\frac{\sum f_{S}(r)f_{L}}{\sum\delta
T_{L}^{2}}\right>+\mathcal{O}\left(\left(\frac{f_{L}}{\sigma_{T_{L}}}\right)^{2}\right)$
(76) $\displaystyle=-f_{S}(r)f_{L}\left<\frac{N_{\text{object}}}{\sum\delta
T_{L}^{2}}\right>+\mathcal{O}\left(\left(\frac{f_{L}}{\sigma_{T_{L}}}\right)^{2}\right)$
$\displaystyle\simeq-f_{S}(r)f_{L}\frac{1}{\left<\delta
T_{L}^{2}\right>}+\mathcal{O}\left(\left(\frac{f_{L}}{\sigma_{T_{L}}}\right)^{2}\right)$
$\displaystyle\simeq-\frac{f_{S}(r)f_{L}}{\sigma_{T_{L}}^{2}}+\mathcal{O}\left(\left(\frac{f_{L}}{\sigma_{T_{L}}}\right)^{2}\right).$

Next, we evaluate the foreground bias for the signed estimator. The stacked
signed estimator including large- and small-scale foregrounds is

$\hat{\tau}^{\rm Sgn}=-\frac{\sum\mathrm{sgn}(\delta T_{L}+f_{L})(\delta
T_{S}+f_{S})}{\sum|\delta T_{L}+f_{L}|}.$ (77)

The bias to the stacked signed estimator to first order in
$f_{L}/\sigma_{T_{L}}$ is

$\displaystyle\hat{\tau}_{\rm bias}^{\rm sgn}(r)$
$\displaystyle=-\left<\frac{\sum\mathrm{sgn}(\delta
T_{L}+f_{L})f_{S}(r)}{\sum|\delta
T_{L}|}\right>+\mathcal{O}\left(\left(\frac{f_{L}}{\sigma_{T_{L}}}\right)^{2}\right)$
(78) $\displaystyle\simeq-f_{S}(r)\frac{\left<\sum\mathrm{sgn}(\delta
T_{L}+f_{L})\right>}{\left<\sum|\delta
T_{L}|\right>}+\mathcal{O}\left(\left(\frac{f_{L}}{\sigma_{T_{L}}}\right)^{2}\right).$

To determine the value of $\left<\sum\mathrm{sgn}(\delta T_{L}+f_{L})\right>$,
consider the case where $\delta T_{L}\sim\mathcal{N}(0,\sigma_{T_{L}})$ and
$f_{L}$ is a small, positive, fixed value. In this case, the distribution of
$\delta T_{L}+f_{L}$ is just a translated Gaussian, i.e., $\delta
T_{L}+f_{L}\sim\mathcal{N}(f_{L},\sigma_{T_{L}})$. Thus, the average value of
the signed estimator is

$\displaystyle\left<\sum\mathrm{sgn}(\delta T_{L}+f_{L})\right>$
$\displaystyle=A\int_{0}^{\infty}\mathrm{exp}\left[-\frac{(x-f_{L})^{2}}{2\sigma_{T_{L}}^{2}}\right]dx$
(79)
$\displaystyle\qquad-A\int_{-\infty}^{0}\mathrm{exp}\left[-\frac{(x-f_{L})^{2}}{2\sigma_{T_{L}}^{2}}\right]dx$
$\displaystyle=A\int_{-f_{L}}^{f_{L}}\mathrm{exp}\left[-\frac{x^{2}}{2\sigma_{T_{L}}^{2}}\right]dx$
$\displaystyle=\sqrt{\frac{2}{\pi}}\frac{f_{L}}{\sigma_{T_{L}}}+\mathcal{O}\left(\left(\frac{f_{L}}{\sigma_{T_{L}}}\right)^{2}\right),$

where $A\equiv\frac{1}{\sqrt{2\pi}\sigma_{T_{L}}}$.

Substituting the result from Eq. 79 into Eq. 78 and using the relation
$\left<|\mathcal{N}(\mu,\sigma)|\right>^{2}=\frac{2}{\pi}\left<\mathcal{N}(\mu,\sigma)^{2}\right>$,
we find

$\displaystyle\hat{\tau}_{\rm bias}^{\rm Sgn}(r)$
$\displaystyle\simeq-\sqrt{\frac{2}{\pi}}\frac{f_{L}}{\sigma_{T_{L}}}\frac{f_{S}(r)}{\left<|\delta
T_{L}|\right>}+\mathcal{O}\left(\left(\frac{f_{L}}{\sigma_{T_{L}}}\right)^{2}\right)$
(80)
$\displaystyle\simeq-\frac{f_{S}(r)f_{L}}{\sigma_{T_{L}}^{2}}+\mathcal{O}\left(\left(\frac{f_{L}}{\sigma_{T_{L}}}\right)^{2}\right).$

Thus, to first order in $f_{L}/\sigma_{T_{L}}$, the TI and signed estimators
have the same sensitivity to foreground contamination. For the signed
estimator, this bias is only caused by regions where $|\delta T_{L}|<f_{L}$.
As sources are weighted by $\delta T_{L}$ in the TI estimator, these regions
contribute very little to the stacked TI profile. For the signed estimator,
the mean bias can be completely mitigated by enforcing a minimum $\delta
T_{L}$ cut such that $|\delta T_{L}|>|f_{L}|$ for all sources.

The analysis in this appendix assumes no overlapping of sources, which would
contribute a non-trivial bias. A detailed analysis of the effect of
overlapping sources on the stacked estimators is left to future work.

##  Appendix F Simulations of correlated extragalactic foregrounds

###  F.1 Simulating a lensing map correlated with $\tau$ map

The deflection angle $\overrightarrow{\alpha}$ due to the gravitational
lensing of photons is given by

$\overrightarrow{\alpha}=\int\frac{d^{2}\overrightarrow{\theta}}{\pi}\kappa(\overrightarrow{\theta})\frac{\hat{\theta}}{\theta}\,$
(81)

where

$\kappa(\overrightarrow{\theta})=\frac{4\pi
G\Sigma(\overrightarrow{\theta})}{c^{2}a}\frac{d_{L}d_{LS}}{d_{S}}$ (82)

In the above equations, $\overrightarrow{\theta}$ are the 2D angles in the
plane of the sky, $\kappa$ is the convergence field, $G$ is Newton’s
gravitational constant, $c$ is the speed of light, $a$ is the scale factor,
$d_{L}$, $d_{LS}$, and $d_{S}$ are the comoving distance from observer to the
lens, between the lens and the source, and to the source plane respectively.
$\Sigma(\overrightarrow{\theta})$ is the comoving angular surface density. We
therefore need $\Sigma(\overrightarrow{\theta})$ to get a lensing map
$\kappa(\overrightarrow{\theta})$ correlated with the $\tau$ map. For this, we
first assume a Gaussian projected profile for the angular surface density
which we apply to the point-like galaxy maps. For CMASS galaxies with typical
mass $M=2\times 10^{13}M_{\odot}$ at $z=0.55$, the virial radius is $\approx
1.6^{\prime}$ Schaan2021a . If we assume that the $3\sigma$ span of the
Gaussian surface density profile should cover most of the profile and that all
the mass falls within this then we have $\sigma=1.6^{\prime}/3$ i.e. angular
surface density profile is a Gaussian with $\sigma\approx 0.53^{\prime}$. In
order to determine the normalization factor of this profile, we use the fact
that $\int d^{2}\overrightarrow{r}\Sigma(\overrightarrow{r})=N_{\rm gal}M$
where $N_{\rm gal}$ is the total number of galaxies, $M$ is their mass, and
$d^{2}\overrightarrow{r}={\rm pixel\>area(in\>steradians)}\times
d_{\chi}^{2}$. Here the pixel area corresponds to the area of a single pixel
of the map and $d_{\chi}$ is the average comoving distance to the galaxy
sample. After getting this normalization, we have
$\Sigma(\overrightarrow{\theta})$ which gives us
$\kappa(\overrightarrow{\theta})$ using Eq. 82. Since the $\tau$ map is
created using the same galaxy map as described in App. G, the $\kappa$ map and
$\tau$ map are correlated.

###  F.2 Simulating a tSZ map correlated with the $\tau$ map

We make some simple assumptions to simulate a tSZ map correlated with the
$\tau$ map. For example, we see from Fig. 8 of Schaan2021a  that the amplitude
of the tSZ profile of CMASS galaxies is around $-10\mu K\>{\rm arcmin}^{2}$ at
150 GHz. We then assume that the tSZ profile is Gaussian in nature centered at
each galaxy location. The profile from Fig. 8 of Schaan2021a  extends out to
$6^{\prime}$ and is a beam convolved profile where the beam full width half
maximum (FWHM) is $1.4^{\prime}$. Therefore, debeamed FWHM for our simulated
tSZ profiles is $\sqrt{6^{2}-1.4^{2}}\approx 5.8^{\prime}$. Similar to
lensing, in order to determine the correct normalization for the Gaussian
profile, we use the fact that $\int d^{2}\overrightarrow{\theta}T=N_{\rm
gal}\times-10\mu K\>{\rm arcmin}^{2}$ where $T$ is the tSZ profile in $\mu K$
and $\int d^{2}\overrightarrow{\theta}T=({\rm
pixel\>area\>in\>{arcmin}^{2}})\sum_{i}T_{i}$ where $\sum_{i}T_{i}$ is the
temperature profile summed over all the pixels. The normalization is therefore
determined such that this relation holds and we get our simulated tSZ map.

We then calculate the power spectrum of our tSZ map and compare it with the
expected power from galaxies with halo mass $5\times
10^{12}M_{\odot}<M<10^{14}M_{\odot}$ and redshift $0.4<z<0.7$ i.e. mass and
redshift range similar to that of CMASS galaxies. We use the halo model
predictions from Maniyar_2021  for this purpose. This prediction provides an
upper limit on the expected power from the CMASS galaxies. As expected, we
find that the tSZ power from our simulated map is below the halo model
predictions and it also falls down much quicker on small scales. We call this
map as ‘tSZ low’. In order to determine the upper limit for the tSZ bias to
the $\tau$ QE, we thus simulate another tSZ map whose power spectrum roughly
matches with the halo model predictions. For this, we simply convolve our
galaxy map by a Gaussian beam with $\sigma=0.88^{\prime}$ instead of
$\sigma=2.46^{\prime}$ for our first tSZ map, such that the shape of its power
spectrum roughly matches with that of halo model predictions. Then we multiply
it by a factor of 1.3 on map level to match the amplitude of the power
spectrum with predictions. Please note that these two numbers, namely the beam
with $\sigma=0.88^{\prime}$ and extra normalization factor of 1.3 have been
determined through trial and error to match the halo model prediction. This
gives us a second tSZ map which we use to determine the highest expected bias
from the tSZ emission coming from CMASS galaxies to cross-correlation of
$\tau$ reconstructed with a QE with galaxies. We call this map as ‘tSZ high’.

##  Appendix G Simulation Methods for TI SNR Forecast

###  G.1 Signal maps

First, we create a $\tau$ map correlated with CMASS galaxy positions. We
choose map dimensions of $10^{\circ}\times 10^{\circ}$ as a balance between
sufficient statistics and computational efficiency. To make the map, we create
a plate carrée (CAR) projection map of the
$200^{\circ}<\text{RA}<210^{\circ}$, $10^{\circ}<\text{Dec}<20^{\circ}$ region
with a pixel scale of 0.5’. We choose this region to avoid any large masked
regions in the CMASS sample footprint. At each galaxy position, we draw a 2D
isotropic Gaussian with FWHM=5’ and amplitude of 0.000245. These values are
chosen based on the fit to the signed estimator stacked profile using
ACT$\times$unWISE data Coulton_2023 .

As shown in Fig. 6, stacking on the unWISE (and potentially LSST) positions
causes an excess signal from the two-halo term. To avoid artificially boosting
the signal for unWISE and LSST during stacking, we use the measured $\tau$
signal from the CMASS map for all forecasts.

Finally, we create a separate signal map for each of the ACT, SPT, SO and
CMB-S4 forecasts by convolving the map with a Gaussian kernel with a FWHM of
1.6’, 1.6’, 1.4’ and 1.0’, respectively, to match the beam of the
corresponding experiment’s noise map (see below).

###  G.2 Noise maps

The noise maps for the forecast consist of Gaussian random field realizations
with power spectra corresponding to each CMB experiment: ACT, SPT, SO, and
CMB-S4. The spectra are plotted in Fig. 4. For ACT, we use the total power
spectrum measured on the Planck+ACT needlet internal linear combination (ILC)
map, which includes all the CMB anisotropies and is convolved with a 1.6’
Gaussian beam Coulton2023_NILCmaps . For SPT, we use the 5-year projected ILC
residual power spectrum curve including the noise and foreground contributions
Raghunathan2023.954.83R ; Sobrin_2022 . For SO, we use forecast ILC power
spectrum with the “goal” sensitivity specifications from SOScience . For
CMB-S4, we use detector noise-only (i.e. no foregrounds) spectrum assuming
$1\,\mu$K$\cdot$arcmin sensitivity. We add an analytic lensed CMB contribution
to the SPT, SO and CMB-S4 power spectra before convolving with a Gaussian beam
of 1.6’, 1.4’ and 1.0’, respectively. The maps have the same dimensions, pixel
scale and CAR projection as the signal maps.

###  G.3 Filtering & Stacking

We give a basic overview of the filtering and stacking algorithm here but
refer the reader to Schaan2021a  for further details. To obtain unbiased
$T_{L}$ and $T_{S}$ estimates, we apply scale-dependent filtering to our maps.
We find the compensated aperture photometry filter used in Schaan2021a
results in a significant noise bias. Instead, we apply low- and high-pass
filtering to the maps to decorrelate the large- and small-scale temperature
estimates:

$f^{\text{low-pass}}_{\ell}=\begin{cases}1,&\text{if $\ell<2000$}\\\
\cos\left(\frac{(\ell-2000)\pi}{300}\right)&\text{otherwise}\\\ 0&\text{if
$\ell>2150$}\end{cases}$ (83) $f^{\text{high-
pass}}_{\ell}=\begin{cases}0&\text{if $\ell<2350$}\\\
\sin\left(\frac{(\ell-2350)\pi}{300}\right)&\text{otherwise}\\\ 1&\text{if
$\ell>2500$}.\end{cases}$ (84)

For the stacked analysis we need catalogs for each galaxy sample. For the
CMASS and unWISE samples, the catalog consists of the galaxies that lie in the
region $200^{\circ}<\text{RA}<210^{\circ}$,
$10^{\circ}<\text{Dec}<20^{\circ}$. As the DC2 footprint does not cover this
region, we simply take a 10 $\times$ 10 deg2 patch from the DC2 footprint
(specifically $55^{\circ}<\text{RA}<65^{\circ}$,
$-40^{\circ}<\text{Dec}<-30^{\circ}$) and shift these positions to the above
region. Using real (or realistically simulated) galaxy positions ensures the
effects of clustering are taken into account in our forecast. The CMASS,
unWISE, and DC2 samples have 7659, 515,095 and 11,904,577 galaxies in this
region, respectively.

After applying the low- and high-pass filtering, we use bilinear interpolation
to make small $\sim$10’$\times$10’ cutouts of the temperature maps centered on
each galaxy position. We estimate $T_{L}$ for a galaxy position by measuring
the mean temperature of the low-pass filtered map in a $\sim$9’ disk centered
on the galaxy. The $T_{S}$ radial profile is estimated using nine annular ring
aperture photometry filters centered on the galaxy with outer radii equally-
spaced between 1’ and 6’. In each aperture $i$, we measure the mean
temperature of the high-pass filtered map, $T_{S,i}$. The estimated mean
$\tau$ in each aperture is calculated as

$\hat{\tau_{i}}=-\frac{T_{L}T_{S,i}}{(T_{L})^{2}}.$ (85)

The optimal stacked estimator is inverse-variance weighted, as shown in Sec.
IV.3. In the forecast, we use uniform weighting for simplicity. Thus, the
stacked $\tau$ profile simplifies to

$\hat{\tau}=-\frac{\sum_{j=1}^{N_{\rm gal}}T_{L,j}T_{S,j}}{\sum_{j=1}^{N_{\rm
gal}}(T_{L,j})^{2}}.$ (86)

###  G.4 Covariance Estimation & SNR Calculation

To estimate the covariance matrix, we simulate 128 independent noise map
realizations, filter and stack in the same manner as above and measure a mean
$\tau$ profile for each realization. We calculate the covariance as the sample
covariance over the 128 mean $\tau$ profiles. The Hartlap correction
Hartlap2007  to this covariance matrix is
$(n_{\mathrm{sims}}-n_{\mathrm{bins}}-2)/(n_{\mathrm{sims}}-1)=1.08$ for 128
simulations and 9 profile bins.

Using simulations to estimate the covariance matrix is superior to bootstrap
resampling as it accounts for the covariance added by the clustering of the
galaxies. We find the difference between the covariance estimated from
simulations versus from 10,000 bootstrap samples is negligible for the CMASS
sample. However, bootstrap resampling underestimates the covariance by $\sim
20$% for the unWISE sample and by more than 50% for the LSST sample, as shown
in Fig. 5.

For a given experiment with noise covariance matrix, $C$, and mean $\tau$
profile, $\mathbf{v}$, we measure the detection SNR of the signal profile over
the $10^{\circ}\times 10^{\circ}$ map as

$\left(\frac{S}{N}\right)^{2}_{\rm map}=\mathbf{v}^{\rm T}C^{-1}\mathbf{v}.$
(87)

We rescale this result to the full survey sky fraction and galaxy number
density as shown in Table 2.

##  Appendix H SNR calculation for the $\tau$ quadratic estimator

For the QE of $\tau$, the SNR for the cross-power spectrum of the $\tau$ map
reconstructed using QE and galaxy map $g$ is calculated analytically as

$\left(\frac{S}{N}\right)^{2}=f_{\rm sky}\sum_{{\ell_{b}}_{\rm
min}}^{{\ell_{b}}_{\rm
max}}\frac{(2{\ell_{b}}+1)\Delta\ell\left(C_{\ell_{b}}^{\tau
g}\right)^{2}}{\left(C_{\ell_{b}}^{\hat{\tau}g}\right)^{2}+C_{\ell_{b}}^{\hat{\tau}\hat{\tau}}\times
C_{\ell_{b}}^{gg}}\,,$ (88)

where

$C_{\ell_{b}}=\frac{1}{\Delta\ell}\sum_{\ell\in[\ell_{1},\ell_{2}]}C_{\ell}\,,$
(89)

and $\Delta\ell$ is the bin width. Here, $C_{\ell_{b}}^{\tau g}$ is the cross-
power spectrum between input $\tau$ map and $g$ map,
$C_{\ell_{b}}^{\hat{\tau}\hat{\tau}}=C_{\ell_{b}}^{\tau\tau}+N_{\ell_{b}}^{\tau\tau}$
where $C_{\ell_{b}}^{\hat{\tau}\hat{\tau}}$ is the power spectrum of the input
$\tau$ map and $N_{\ell_{b}}^{\tau\tau}$ is the reconstruction noise for QE
given by Eq. 12. $C_{\ell_{b}}^{\hat{\tau}\hat{\tau}}$ is therefore, power
spectrum of the reconstructed $\tau$ map. $C_{\ell_{b}}^{gg}$ is the total
galaxy power spectrum including the clustering and shot noise components, and
$f_{\rm sky}$ is the common sky fraction between the CMB map and galaxy
survey. We start with the $\tau$ map and CMASS galaxy map as described in App.
G, and cross-correlate them to calculate the signal $C_{\ell}^{\tau g}$. For
the noise part, we take into account different clustering and shot noise
levels in the CMASS, unWISE, and LSST galaxy samples. For CMASS sample and
LSST galaxies, we calculate the galaxy clustering power spectrum using a halo
model approach where we make sure that our model matches with the clustering
power spectrum of the CMASS maps used in cross-correlations. For unWISE
galaxies, we use the measured power spectra for the ‘blue’ sample
Krolewski_2020 ; Krolewski_2021 . The shot noise levels for each sample are
calculated using the number density values provided in Tab. 2.

For different CMB experiments, the reconstruction noise for the quadratic
estimator for $\tau$ is different and is calculated using Eq. 12.