# The Atacama Cosmology Telescope: Detection of Patchy Screening of the Cosmic
Microwave Background

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  Theo Schutt
0000-0002-7187-9628 SLAC National Accelerator Laboratory 2575 Sand Hill Road
Menlo Park, California 94025, USA  Kavli Institute for Particle Astrophysics
and Cosmology, 382 Via Pueblo Mall Stanford, CA 94305-4060, USA  Department of
Physics, Stanford University, 382 Via Pueblo Mall, Stanford, CA 94305, USA
Abhishek S. Maniyar  0000-0002-4617-9320 SLAC National Accelerator Laboratory
2575 Sand Hill Road Menlo Park, California 94025, USA  Kavli Institute for
Particle Astrophysics and Cosmology, 382 Via Pueblo Mall Stanford, CA
94305-4060, USA  Department of Physics, Stanford University, 382 Via Pueblo
Mall, Stanford, CA 94305, USA  Emmanuel Schaan  0000-0002-4619-8927 SLAC
National Accelerator Laboratory 2575 Sand Hill Road Menlo Park, California
94025, USA  Kavli Institute for Particle Astrophysics and Cosmology, 382 Via
Pueblo Mall Stanford, CA 94305-4060, USA  Rui An  Department of Physics and
Astronomy, University of Southern California, Los Angeles, CA 90089, USA
Zachary Atkins  0000-0002-2287-1603 Joseph Henry Laboratories of Physics,
Jadwin Hall, Princeton University, Princeton, NJ, USA 08544  Nicholas
Battaglia  Department of Astronomy, Cornell University, Ithaca, NY 14853, USA
J Richard Bond  Canadian Institute for Theoretical Astrophysics, University of
Toronto, Toronto, ON, Canada M5S 3H8  Erminia Calabrese  0000-0003-0837-0068
School of Physics and Astronomy, Cardiff University, The Parade, Cardiff,
Wales CF24 3AA, UK  Steve K. Choi  0000-0002-9113-7058 Department of Physics,
Cornell University, Ithaca, NY 14853, USA  Department of Astronomy, Cornell
University, Ithaca, NY 14853, USA  Mark J. Devlin  Department of Physics and
Astronomy, University of Pennsylvania, 209 South 33rd Street, Philadelphia,
PA, USA 19104  Adriaan J. Duivenvoorden  0000-0003-2856-2382 Center for
Computational Astrophysics, Flatiron Institute, New York, NY 10010, USA,
Joseph HenrLaboratories of Physics, Jadwin Hall, Princeton University,
Princeton, NJ, USA 08544  Jo Dunkley  Joseph Henry Laboratories of Physics,
Jadwin Hall, Princeton University, Princeton, NJ, USA 08544  Department of
Astrophysical Sciences, Peyton Hall, Princeton University, Princeton, NJ USA
08544  Simone Ferraro  Physics Division, Lawrence Berkeley National
Laboratory, Berkeley, CA, USA  Department of Physics, University of
California, Berkeley, CA, USA 94720  Vera Gluscevic  Department of Physics and
Astronomy, University of Southern California, Los Angeles, CA 90089, USA  J.
Colin Hill  0000-0002-9539-0835 Department of Physics, Columbia University,
New York, NY, USA  Center for Computational Astrophysics, Flatiron Institute,
162 5th Avenue, New York, NY 10010 USA  Matt Hilton  Wits Centre for
Astrophysics, School of Physics, University of the Witwatersrand, Private Bag
3, 2050, Johannesburg, South Africa  School of Mathematics, Statistics &
Computer Science, University of KwaZulu-Natal, Westville Campus, Durban 4041,
South Africa  Adam D. Hincks  0000-0003-1690-6678 David A. Dunlap Department
of Astronomy & Astrophysics, University of Toronto, 50 St. George St., Toronto
ON M5S 3H4, Canada  Specola Vaticana (Vatican Observatory), V-00120 Vatican
City State  Arthur Kosowsky  Department of Physics and Astronomy, University
of Pittsburgh, Pittsburgh PA 15260 USA  Darby Kramer  School of Earth and
Space Exploration, Arizona State University, Tempe, AZ, USA 85287  Aleksandra
Kusiak  Department of Physics, Columbia University, New York, NY, USA  Adrien
La Posta  0000-0002-2613-2445 University of Oxford, Denys Wilkinson Building,
Keble Road, Oxford, OX1 3RH, UK  Thibaut Louis  0000-0002-6849-4217 Université
Paris-Saclay, CNRS/IN2P3, IJCLab, 91405 Orsay, France  Mathew S. Madhavacheril
Department of Physics and Astronomy, University of Pennsylvania, 209 South
33rd Street, Philadelphia, PA, USA 19104  Gabriela A. Marques  Fermi National
Accelerator Laboratory, P. O. Box 500, Batavia, IL 60510, USA  Kavli Institute
for Cosmological Physics, University of Chicago, Chicago, IL 60637, USA  Fiona
McCarthy  DAMTP, Centre for Mathematical Sciences, University of Cambridge,
Wilberforce Road, Cambridge CB3 OWA, UK  Center for Computational
Astrophysics, Flatiron Institute, 162 5th Avenue, New York, NY 10010 USA  Jeff
McMahon  Kavli Institute for Cosmological Physics, University of Chicago, 5640
S. Ellis Ave., Chicago, IL 60637, USA  Department of Astronomy and
Astrophysics, University of Chicago, 5640 S. Ellis Ave., Chicago, IL 60637,
USA  Department of Physics, University of Chicago, Chicago, IL 60637, USA
Enrico Fermi Institute, University of Chicago, Chicago, IL 60637, USA  Kavilan
Moodley  000-0001-6606-7142 Astrophysics Research Centre, University of
KwaZulu-Natal, Westville Campus, Durban 4041, South Africa  School of
Mathematics, Statistics & Computer Science, University of KwaZulu-Natal,
Westville Campus, Durban 4041, South Africa  Sigurd Naess  0000-0002-4478-7111
Institute of Theoretical Astrophysics, University of Oslo, Norway  Lyman A.
Page  Joseph Henry Laboratories of Physics, Jadwin Hall, Princeton University,
Princeton, NJ, USA 08544  Bruce Partridge  Department of Physics and
Astronomy, Haverford College, Haverford, PA, USA 19041  Frank J. Qu  DAMTP,
Centre for Mathematical Sciences, University of Cambridge, Wilberforce Road,
Cambridge CB3 OWA, UK  Kavli Institute for Cosmology Cambridge, Madingley
Road, Cambridge CB3 0HA, UK  Neelima Sehgal  0000-0002-9674-4527 Physics and
Astronomy Department, Stony Brook University, Stony Brook, NY 11794  Blake D.
Sherwin  DAMTP, Centre for Mathematical Sciences, University of Cambridge,
Wilberforce Road, Cambridge CB3 OWA, UK  Kavli Institute for Cosmology,
University of Cambridge, Madingley Road, Cambridge CB3 OHA, UK  Cristóbal
Sifón  0000-0002-8149-1352 Instituto de Física, Pontificia Universidad
Católica de Valparaíso, Casilla 4059, Valparaíso, Chile  David N. Spergel
Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue, New
York, NY 10010 USA  Department of Astrophysical Sciences, Peyton Hall,
Princeton University, Princeton, NJ USA 08544  Suzanne T. Staggs
0000-0002-7020-7301 Joseph Henry Laboratories of Physics, Jadwin Hall,
Princeton University, Princeton, NJ, USA 08544  Alexander Van Engelen  School
of Earth and Space Exploration, Arizona State University, Tempe, AZ, USA 85287
Cristian Vargas  0000-0001-5327-1400 Instituto de Astrofísica and Centro de
Astro-Ingeniería, Facultad de Fìsica, Pontificia Universidad Católica de
Chile, Av. Vicuña Mackenna 4860, 7820436 Macul, Santiago, Chile  Edward J.
Wollack  0000-0002-7567-4451 NASA/Goddard Space Flight Center, Greenbelt, MD,
USA 20771

###### Abstract

Spatial variations in the cosmic electron density after reionization generate
cosmic microwave background anisotropies via Thomson scattering, a process
known as the “patchy screening” effect. In this paper, we propose a new
estimator for the patchy screening effect that is designed to mitigate biases
from the dominant foreground signals. We use it to measure the cross-
correlation between unWISE galaxies and patchy screening, the latter measured
by the Atacama Cosmology Telescope and Planck satellite. We report the first
detection of the patchy screening effect, with the statistical significance of
the cross-correlation exceeding $7\sigma$. This measurement directly probes
the distribution of electrons around these galaxies and provides strong
evidence that gas is more extended than the underlying dark matter. By
comparing our measurements to electron profiles extracted from simulations, we
demonstrate the power of these observations to constrain galaxy evolution
models. Requiring only the 2D positions of objects and no individual redshifts
or velocity estimates, this approach is complementary to existing gas probes,
such as those based on the kinetic Sunyaev-Zeldovich effect.

###### pacs:

Valid PACS appear here

††preprint: APS/123-Q ED

##  I Introduction

The distribution of gas throughout the Universe contains a wealth of
information on the astrophysical processes governing the thermodynamics of the
circumgalactic and intracluster media (ICM)[1, 2, 3, 4, 5]. It is also a key
ingredient in cosmological analyses, especially galaxy weak lensing analyses
as highlighted in e.g., Refs. [6, 7]. Whilst the properties of this gas in the
centers of massive clusters have been well characterized by X-ray and thermal
Sunyaev-Zeldovich (tSZ) measurements [e.g., 3, 8, 9, 10, 11, 12, 13], the gas
distribution in cluster outskirts and lower mass systems is less well
understood. Recent measurements of the kinetic Sunyaev-Zeldovich (kSZ) effect
offer one means of studying this gas [14, 15, 16, 17, 18]. In this work we
explore a complementary approach based on the “patchy screening” effect which,
like the “projected fields” kSZ estimator [19, 20, 21], does not require
spectroscopic data.

Cosmic microwave background (CMB) photons are Thomson scattered by electrons
encountered on their path from the surface of last scattering to the observer.
This scattering damps the primary CMB anisotropies. The isotropic component
scales with the optical depth to reionization and is one of the six
$\Lambda$CDM parameters [22, 23, 24, 25]. However, the damping is anisotropic
as lines of sight with more (fewer) electrons will be more (less) suppressed,
and thereby spatially modulates the underlying CMB anisotropies [26, 27, 28,
29, 30, 31, 32, 33]. This effect is known as the “patchy screening” effect
and, on small scales ($\theta\lesssim 1^{\circ}$), the new anisotropies are

$\displaystyle\Delta T^{\text{patchy
screen}}(\mathbf{n})=-\delta\tau(\mathbf{n})\Delta
T^{\mathrm{primary}}(\mathbf{n}),$ (1)

where $\delta\tau(\mathbf{n})$ is the fluctuation in the optical depth to last
scattering in direction $\mathbf{n}$ and $\Delta
T^{\mathrm{primary}}(\mathbf{n})$ is from the primordial CMB anisotropy in
that direction. The physical intuition is that on small scales the primary CMB
has no power, so patchy screening is simply the damping of the larger scale
CMB fluctuations by the small scale variations in $\tau$. See Ref. [33, 32] or
the Supplementary Material for more details.

Patchy screening is important during two epochs of the Universe: 1) during the
reionization of the Universe, when the Universe has large spatial fluctuations
in the cosmic ionization fraction, and 2) in the late-time Universe, when
excesses in the density of electrons occur around gravitationally collapsed
objects [e.g. 34, 35, 36, 37]. The observational signature from the latter
epoch is damping of the primary CMB anisotropies at a scale of a few
arminutes. These anisotropies are well suited to studying the distribution of
gas in the Universe as the damping is linearly proportional to the excess of
electrons (for the typical small values of optical depths encountered). Thus,
the angular size and damping amplitude characterize the projected gas profile.
Patchy screening is complementary to the kinetic Sunyaev-Zeldovich (kSZ)
effect, which arises from Thomson scattering of CMB photons off moving
electrons. Both effects probe the optical depth, but with different
weightings: patchy screening depends on the background primary CMB
anisotropies whilst the kSZ effect depends on the line-of-sight velocity.

Measuring the late-time patchy screening signal requires two components: a
precise measurement of the large scale CMB anisotropies, as a measure of the
CMB before damping, and low-noise measurements of the new small scale
anisotropies. Ground based CMB experiments, such as the Atacama Cosmology
Telescope (ACT), South Pole Telescope, the Simons Observatory, CMB-HD, and
CMB-S4 [38, 39, 40, 41, 42], are designed for the sensitivities and resolution
necessary to provide the latter component. These experiments are complemented
by the precise measures of large scale CMB anisotropies made by the Planck and
WMAP satellites [43, 44]. Exploiting this synergy, we use data from Planck and
ACT to measure patchy screening at small scales.

CMB experiments are only sensitive to the patchy screening signal integrated
along the line of sight. To understand the host mass and redshift dependence
of the signal, we cross-correlate the CMB data with another tracer of the gas
distribution [37]. If the second tracer is measured at high precision, then
cross correlation based methods also provide a means of boosting the
detectability of the patchy screening effect. In this work we cross correlate
the ACT-Planck CMB measurements with galaxies in the unWISE catalog, using a
stacking-based approach [45, 46, 47]. The unWISE galaxy catalog contains over
$500$ million galaxies as obtained from the extended observations of the Wide-
field Infrared Survey Explorer (WISE) satellite [48, 49]. The high source
number density of this catalog makes it ideal for studying the patchy
screening signal.

##  II Cross-Correlation Methodology

A key challenge to detecting patchy screening is the host of other small-scale
millimeter sky signals. These include emission from high redshift star-forming
galaxies and radio galaxies, as well as CMB secondaries such as gravitational
lensing, and the thermal and kinetic Sunyaev-Zeldovich effects [see e.g., 50,
for a review]. These foreground signals are expected to be at least 30 times
larger than the patchy screening signal and are spatially correlated with it.
To avoid potential biases, we combine well-tested component-separation methods
with a new estimator called the “sign estimator.”

The sign estimator uses two parts to isolate the patchy screening effect: the
first part is a map of the large scale CMB anisotropies, which quantifies how
much damping can occur. Current multifrequency measurements, such as those
from the Planck satellite, can be combined to produce a map of the large-scale
anisotropies that is dominated by the primary CMB. The sign operation, a non-
linear function which sets all pixels to $\pm 1$, is applied to this map. This
operation isolates the primary CMB, and eliminates contaminating signals, as
the sign of each pixel of the map is set by the dominant component. The second
part is to multiply this map by a map of the small scale CMB anisotropies.
Except for the patchy screening effect, the sign of the primary CMB
anisotropies is uncorrelated with the small scale sky signals; thus the
foreground contribution will average to zero and leave the patchy screening.
This robustness to foregrounds means that we can apply our estimator to maps
containing bright contaminants without biases. The sign estimator trades
optimality for increased immunity to foreground biases. This is one of several
conservative choices made to ensure a robust detection. In the Supplementary
Material we show how this approach is motivated by the squeezed limit of the
standard quadratic estimator. An in-depth comparison of different patchy
screening estimators, along with forecasts for current and future experiments,
is given in Ref. [51].

We detect patchy screening via its correlation with the distribution of
galaxies. Specifically, we apply the sign estimator to cutouts of the CMB data
extracted at the location of known galaxies. These cutouts are then averaged
together in a 2D stack. A detailed description of our methodology is provided
in the Supplementary Material.

##  III Data Sets

The CMB data used in this work are the component-separated CMB maps obtained
from combining data from the Planck satellite with data from ACT, as detailed
in Ref. [52]. Here we briefly review the key products and refer the reader to
Ref. [52] for more details.

The Planck satellite observed the full sky over four years with an angular
resolution of $7.22^{\prime}$ at $143\,$GHz [43]. The component-separated maps
are constructed from the NPIPE Planck data release [53] and use data between
$30\,$GHz and $545\,$GHz. These maps provide the large-scale CMB measurement;
the small scales come from high resolution measurements of the millimeter sky
by ACT. ACT observed around one third of the sky over the course of 15 years
[54, 55, 56, 57, 58, 59, 60, 61, 62]. This work uses data from Data Release 4
(DR4) and Data Release 6 (DR6), which contains data from 2013-2023. It
comprises observations at three frequencies: 98, 148, and 225 GHz (known as
f090, f150, f220). At $148\,$GHz the ACT beam full-width-half-maximum (FWHM)
is $1.4$′. Typical noise levels in the ACT maps are 13/12/48 $\mu$K arcmin at
f090/f150/f220. Atmospheric noise, primarily from the vibrational and
rotational modes of water vapor, has a red spatial power spectrum and inhibits
ACT from making observations of the large scale CMB anisotropies [63]. Thus
the combination of Planck and ACT data is well suited to studying the patchy
screening signal.

The ACT and Planck multifrequency observations are combined using a needlet
internal linear combination (NILC) to produce a map isolating the CMB
anisotropies. This CMB map is dominated by the primary CMB anisotropies on
angular scales ranging from full-sky to a few arcminutes. The resulting map is
convolved to a 1.6′ FWHM Gaussian beam. The NILC map we use here was optimized
to produce an unbiased CMB map with minimal “noise” variance, with noise
defined to be all signals (sky or instrumental) that are not the component of
interest [64, 65]. As discussed further in the Supplementary Material, the
robustness of our measurement can be tested by making variants of the NILC
maps. These variants are produced by either using only subsets of the data
during the NILC process (if only one frequency is used the NILC method can be
thought of as a means of coadding data) or by applying constraints to ensure
that a known contaminant signal does not contribute to the map “noise.” To
achieve this latter goal, we use the constrained NILC method [66, 67].

Sample | Mean Redshift | Mean Halo Mass | Source density
---|---|---|---
|  | (M⊙/h) | (arcmin-2 )
Blue | 0.6 | $1.4\times 10^{13}$ | 0.58
Green | 1.1 | $1.3\times 10^{13}$ | 0.32
Red | 1.5 | $1.6\times 10^{13}$ | 0.018
Table 1: The key properties of the unWISE galaxies as measured by Refs [46,
68, 69]. The source density is computed after our catalog cuts (see
Supplementary Material).

The galaxy catalogs used in our cross correlation analysis come from the Wide-
field Infrared Survey Explorer (WISE) satellite [48]. The WISE satellite
observed the sky from 2010 to 2011 and then again from 2013 to the present. We
use the “blue,” “green” and “red” galaxy samples from Ref. [46] and refer the
reader to Ref. [46] for details on the selection criteria. These samples,
whose properties are summarized in Table 1, are obtained from the deep 3.4 and
4.6 micron unWISE catalogs [70, 71, 49], with Gaia data used to help remove
contaminant stars [72]. The high source densities and large overlap with ACT
are the key reasons for using the unWISE catalogs – note that the patchy
screening cross-correlation only needs the galaxy positions and not their
redshifts.

In the Supplementary Materials we describe the simulations used to validate
our analysis and to construct the covariance matrices used to interpret our
results. We also detail a data-based test of the covariance matrix.

Figure 1:  High-pass filtered 2D stacked optical depth map, obtained by
running the “sign” estimator on ACT data, at the positions of the unWISE blue
sample. The optical depth map is convolved with the filter of Eq. (13), on top
of the 1.6′ FWHM Gaussian ACT beam convolution. The visible excess in the
center comes from the gas around the unWISE galaxies.

##  IV Results

Figure 2: Stacked 1D profiles of the optical depth, obtained from the sign
estimator, Section B.3, applied to filtered CMB data stacked at the locations
of galaxies. The points show the measurements of the unWISE red, green, blue
and combined blue-green samples. These are obtained by azimuthally averaging
the 2D stacks. After accounting for the covariance of the points, the combined
sample rejects the null hypothesis of no signal at $7.9\,\sigma$ (see Table
2). The dot-dashed line is the signal for a point source and the dotted line
is a Gaussian electron density profile, with parameters fit to the combined
blue-green sample. The point-source and Gaussian profiles are convolved to the
beam (1.6′ FWHM) and filtered in the same manner as the data. The error bars
are computed from Monte-Carlo simulations as described in the Supplementary
Material.  Figure 3: The correlation matrices of the $\tau$ measurements from
the unWISE blue and green samples, computed from 150 independent simulations
that do not include the patchy screening signal. The lower left corner of the
plot is the blue sample, the upper right is the green sample and the diagonals
are the cross terms. Significant off-diagonal correlations arise from the
high-pass filtering. No strong correlations are detected between the two
samples.  Figure 4:  A qualitative comparison of the measured profile from the
blue sample to a set of theoretical models of the electron distribution: the
Illustris (red), Illustris-TNG (green) and the SIMBA (purple) simulations, and
the best-fit Generalized Navarro Frenk White profile from the kSZ analysis of
BOSS galaxies[16] (orange). The dark matter profile from Illustris-TNG is
shown in dot-dashed grey, re-normalized to match the integrated electron
density. Each of the theoretical profiles has been filtered in a manner
equivalent to the data. The black lines are samples from the fit to the
Gaussian model. Perfect agreement between the data and simulations comparison
is not expected as this comparison does not account for the observational
selection function. The qualitative comparison already suggests that the gas
does not follow the dark matter, and that future measurements will contain
useful information about feedback and numerical simulations of galaxy
formation [5].

In Fig. 1 we present the 2D stack of the sign-weighted CMB data on the unWISE
blue galaxies, where the patchy screening signal can be clearly seen. We show
the azimuthally averaged 1D profiles of optical depth, computed with Section
B.3, for the three unWISE samples in Fig. 2. Whilst the highest redshift stack
(red) is consistent with zero, we detect signals in the blue and green samples
at high signal-to-noise ratio (SNR). This difference is expected given that
the red catalog has an order of magnitude fewer members than either of the
other catalogs. Note that the shape of the profiles, especially the “ringing”
structure, arises from the filtering of the CMB map.

We compute the covariance of the data points, including the covariance between
the green and blue samples, using Monte-Carlo simulations, described in the
Supplementary Material. In Fig. 3 we show the correlation matrix for the blue
and green samples. The strong radial correlations induced by the filtering can
clearly be seen. These correlations are not due to the signal – the
simulations do not include the patchy screening signal – and arise from the
primary CMB, noise and other sky components. These signal-free simulations
quantify the expected noise-induced variance in the profiles. The comparison
of our measurement to this noise quantifies the detection significance. The
correlations between the two samples are consistent with noise. We combine the
blue and green data points – with an inverse covariance weighting– to obtain a
higher SNR sample, also shown in Fig. 2.

Table 2 provides $\chi^{2}$ and $p$-values estimated from the data for two
cases: the null hypothesis, and a two-parameter model described below. The
large $\chi^{2}$ values and low $p$-values for the null hypothesis demonstrate
the high significance of the detection. The final column of Table 2 translates
the $p$-values to the equivalent number of Gaussian sigmas, showing that the
profile from the combination of the blue and green samples for the nine radial
bins of Fig. 2 is inconsistent with the null hypothesis at $7.9\sigma$.
However, quantifying the detection significance from just the null hypothesis
is non-trivial as the $\chi^{2}$ values will clearly depend on the maximum
radius used in the measurement. Note that if we just use three radial points
we still reject the null hypothesis by more than $7\sigma$.

We fit the data to a two-parameter model in order to perform a model
comparison test to better assess the detection significance. In this
phenomenological model the optical depth around the unWISE galaxies is modeled
as a 2D symmetric Gaussian distribution. To this model we apply the CMB beam,
filter in an identical manner to the data and compute the 1D azimuthal
average; then we fit for the amplitude and width. The best fit profile is
plotted in Fig. 2, and its $\chi^{2}$ statistics are in Table 2.

We consider two statistical tests to compare the Gaussian model to the null
hypothesis: the Akaike Information Criterion (AIC) and an F-test [73, 74, 75].
Given the small size of data vector (9 elements), the Bayesian Information
Criterion (a third commonly used model comparison metric) is effectively the
same as the AIC discussed here and hence we only consider the AIC. We find AIC
values of $-39$, $-43$ and $-71$ for the green, blue and combined samples
respectively. These values correspond to strong evidence for the Gaussian
model and hence strong evidence of a non zero signal. With the F-test we find
p-values of $p=0.003,0.013$ and $0.0018$, which also shows evidence for a non-
zero signal model.

Rejection of null hypothesis
---
Sample | $\chi_{\text{null}}^{2}$ |  $p$-value | $N_{\sigma}$
Blue | 67 | $6.8\times 10^{-11}$ | 6.4
Green | 53 | $2.2\times 10^{-8}$ | 5.5
Red | 3.5 | 0.93 | -1.5
Combined blue and green | 90 | $1.4\times 10^{-15}$ | 7.9

Goodness of fit for the 2 parameter Gaussian model
---
Sample | $\chi_{\text{best fit}}^{2}$ |  $p$-value
Blue | 19 | 0.007
Green | 11 | 0.16
Combined blue and green | 15 | 0.04

Table 2: The $\chi^{2}$, $p$-values and equivalent Gaussian sigma significance
( $N_{\sigma}$) of our measurements (top), and of fits with a Gaussian profile
(bottom). Note that the statistics above used the full data covariance for the
nine radial bins in each data set.

###  IV.1 Discussion

Performing a quantitative interpretation of the measured gas profiles, in the
vein of Ref. [16, 76], is challenging due to the large measurement errors and
large uncertainties in the properties of the unWISE galaxies. We leave this
challenge to future work, and restrict our analysis to qualitative
comparisons. We contrast our measurements with a theoretical model and three
state-of-the-art cosmological hydrodynamical simulations: Illustris [77],
IllustrisTNG [78] and SIMBA[79]. Each simulation models the unresolved small-
scale physics differently and the range provides some measure on the diversity
of current theoretical predictions. For each simulation we follow the method
of [80, 81] to compute the profiles. We use the unWISE properties from Ref.
[68, 69] to perform the comparisons, with the mean redshifts and halo masses
listed in Table 1. For the simulations, we use the nearest redshift snapshot
and a mass bin of $1.5\times 10^{13}M_{\odot}/h\leq
M_{\mathrm{halo}}<2.4\times 10^{13}M_{\odot}/h$. We do not attempt to
reproduce the unWISE sample selections, nor the impact of satellite galaxies,
and therefore do not expect a perfect agreement.

First, we compare the measured profiles to the profile from a point source.
Fig. 2 shows that our data strongly favours a more extended source for the
signal. The electrons are not expected to be distributed as a point source.
This comparison demonstrates that our measurements can probe the spatial
distribution of the gas and provides a check on potential point source
contamination. We then compare the measurements of the blue and green samples.
From the fits to the samples’ profiles, Fig. 2, we find that the green and
blue sample profiles are fit by Gaussian models of widths
$\theta_{\mathrm{FWHM}}=5.3\pm 0.5$ and $\theta_{\mathrm{FWHM}}=5.7\pm 0.5$
arcmin. Note that without evolution of profiles, the angular size of the green
sample is expected to be smaller due to its higher redshift. These results
highlight the potential of future, higher SNR measurements to characterize the
redshift evolution of the distribution of gas.

Next we compare our measurements to the expected distribution of dark matter.
As can be seen in Fig. 4, the measured gas profiles do not follow the dark
matter (grey). Note that we show the dark matter from the Illustris-TNG
simulation, but similar results are obtained from all three simulations. These
results reproduce those seen in Refs [14, 15, 16, 76], which used kSZ
measurements stacked on BOSS galaxies to demonstrate that the gas is
significantly more extended than the dark matter. The best-fit generalized
Navarro Frenk White (GNFW)[82] profile obtained from Ref. [16] is shown in
orange. Our measurements are more consistent with [16] than the dark matter;
however, our results appear to have a higher gas density and slightly more
extended profile. The unWISE galaxies used in this work are different from the
BOSS galaxies: the unWISE galaxies are infrared selected objects and trace
star forming galaxies, whilst the BOSS galaxies are a constant stellar mass
population and contain both star forming and quiescent galaxies [83]. Of
particular relevance to this analysis is the proportion of galaxies that are
central versus satellite members of their halos. The GNFW and simulation
profiles do not account for the miscentering of the measured profiles that
arises when we stack on satellite galaxies, which are offset from the center
of the halo. The measurements from [84, 68, 69] suggest that the fraction of
unWISE galaxies that are satellites is larger than in the BOSS sample. The
impact of the satellite fraction on the observed profile is complex as it
couples with the selection effects (see Ref. [85]) and can either broaden or
sharpen the profile. These small-scale measurements of the gas around the
unWISE galaxies complement the larger-scale kSZ-based measurements by [17] and
future comparisons between screening and kSZ measurements will be a powerful
astrophysical and cosmological probe.

Finally, we compare to electron density profiles from three state-of-the-art
simulations. This qualitative comparison of the profiles, Fig. 4, shows that
our measurements are broadly similar to all of the simulation suites. Better
agreement is found with the SIMBA and Illustris-TNG simulations, which are
known to have less heating and redistribution of the ICM than Illustris.
Future measurements, with well-characterized galaxy samples, will enable
quantitative comparisons and be highly informative for the next generation of
hydrodynamical simulations.

In summary, CMB anisotropy data have reached the precision necessary to study
the patchy screening effect. In this paper and Ref. [51], we introduce a new
estimator that is robust to foreground biases, enabling the use of arcminue-
scale temperature anisotropies with significant foreground contributions. A
key feature of this analysis technique is that it relies only on the sky
positions of galaxies, and thus complements existing methods of probing the
gas distribution which also require galaxy distance estimates [14, 15, 76].
Our results already show that patchy screening measurements provide a powerful
method to characterize the physics governing the cosmic gas distribution. In
Ref. [51], we show that future CMB observations with the Simons Observatory
[40] and CMB-S4 [42] will allow the precision of these measurements to rapidly
increase and greatly improve differentiation among galaxy formation models.

###### Acknowledgements.

We are very grateful to Leander Thiele for providing the code used to measure
the gas profiles in the Illustris, Illustris-TNG and SIMBA simulations.
Support for ACT was through the U.S. National Science Foundation through
awards AST-0408698, AST-0965625, and AST-1440226 for the ACT project, as well
as awards PHY-0355328, PHY-0855887 and PHY-1214379. Funding was also provided
by Princeton University, the University of Pennsylvania, and a Canada
Foundation for Innovation (CFI) award to UBC. ACT operated in the Parque
Astronómico Atacama in northern Chile under the auspices of the Agencia
Nacional de Investigación y Desarrollo (ANID). The development of multichroic
detectors and lenses was supported by NASA grants NNX13AE56G and NNX14AB58G.
Detector research at NIST was supported by the NIST Innovations in Measurement
Science program. Computing for ACT was performed using the Princeton Research
Computing resources at Princeton University, the National Energy Research
Scientific Computing Center (NERSC), and the Niagara supercomputer at the
SciNet HPC Consortium. SciNet is funded by the CFI under the auspices of
Compute Canada, the Government of Ontario, the Ontario Research Fund–Research
Excellence, and the University of Toronto. We thank the Republic of Chile for
hosting ACT in the northern Atacama, and the local indigenous Licanantay
communities whom we follow in observing and learning from the night sky. JCH
acknowledges support from NSF grant AST-2108536, NASA grants 21-ATP21-0129 and
22-ADAP22-0145, the Sloan Foundation, and the Simons Foundation. KM and MHi
acknowledge support from the National Research Foundation of South Africa. AK
acknowledge support from NSF grant AST-2108536 SN acknowledges support from a
grant from the Simons Foundation (CCA 918271, PBL). GAM is part of the Fermi
Research Alliance, LLC under Contract No. DE-AC02-07CH11359 with the U.S.
Department of Energy, Office of Science, Office of High Energy Physics The
Flatiron Institute is supported by the Simons Foundation. CS and CV
acknowledges support from the Agencia Nacional de Investigación y Desarrollo
(ANID) through BASAL project FB210003. MM acknowledges support from NSF grants
AST-2307727, AST-2153201 and NASA grant 21-ATP21-0145. NS acknowledges support
from DOE award number DE-SC0020441. ZA and JD acknowledge support from NSF
grant AST-2108126. ADH acknowledges support from the Sutton Family Chair in
Science, Christianity and Cultures, from the Faculty of Arts and Science,
University of Toronto, and from the Natural Sciences and Engineering Research
Council of Canada (NSERC) [RGPIN-2023-05014, DGECR-2023- 00180]. 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.

## References

  * [1] J. P. Ostriker, P. Bode, and A. Babul, Astrophys. J. 634, 964 (2005).
  * [2] D. Nagai, A. V. Kravtsov, and A. Vikhlinin, Astrophys. J. 668, 1 (2007).
  * [3] B. R. McNamara and P. E. J. Nulsen, Annual Review of Astronomy and Astrophysics 45, 117 (2007).
  * [4] N. Battaglia et al., Astrophys. J. 725, 91 (2010).
  * [5] E. Moser et al., Astrophys. J. 933, 133 (2022).
  * [6] A. Amon and G. Efstathiou, Mon. Not. R. Astron. Soc.516, 5355 (2022).
  * [7] G. Aricò et al., arXiv e-prints arXiv:2303.05537 (2023).
  * [8] J. J. Ruan, M. McQuinn, and S. F. Anderson, Astrophys. J. 802, 135 (2015).
  * [9] D. Crichton et al., Mon. Not. R. Astron. Soc.458, 1478 (2016).
  * [10] A. de Graaff, Y.-C. Cai, C. Heymans, and J. A. Peacock, Astron. Astrophys.624, A48 (2019).
  * [11] D. Eckert et al., Astron. Astrophys.621, A40 (2019).
  * [12] A. Zenteno et al., Mon. Not. R. Astron. Soc.495, 705 (2020).
  * [13] S. Pandey et al., Phys. Rev. D105, 123526 (2022).
  * [14] E. Schaan et al., Phys. Rev. D93, 082002 (2016).
  * [15] E. Schaan et al., Phys. Rev. D103, 063513 (2021).
  * [16] S. Amodeo et al., Phys. Rev. D103, 063514 (2021).
  * [17] A. Kusiak et al., Phys. Rev. D104, 043518 (2021).
  * [18] E. Schiappucci et al., Phys. Rev. D107, 042004 (2023).
  * [19] O. Doré, J. F. Hennawi, and D. N. Spergel, Astrophys. J. 606, 46 (2004).
  * [20] J. C. Hill et al., Phys. Rev. Lett. 117, 051301 (2016).
  * [21] S. Ferraro et al., Phys. Rev. D94, 123526 (2016).
  * [22] Planck Collaboration Int. XLVI, Astron. Astrophys.596, A107 (2016).
  * [23] L. Pagano et al., Astron. Astrophys.635, A99 (2020).
  * [24] Planck Collaboration Int. LVII, Astron. Astrophys.643, 42 (2020).
  * [25] R. de Belsunce, S. Gratton, W. Coulton, and G. Efstathiou, Mon. Not. R. Astron. Soc.507, 1072 (2021).
  * [26] J. P. Ostriker and E. T. Vishniac, Astrophys. J. Lett.306, L51 (1986).
  * [27] E. T. Vishniac, Astrophys. J. 322, 597 (1987).
  * [28] W. Hu, D. Scott, and J. Silk, Phys. Rev. D49, 648 (1994).
  * [29] S. Dodelson and J. M. Jubas, Astrophys. J. 439, 503 (1995).
  * [30] F. M. Persi, D. N. Spergel, R. Cen, and J. P. Ostriker, Astrophys. J. 442, 1 (1995).
  * [31] A. Gruzinov and W. Hu, Astrophys. J. 508, 435 (1998).
  * [32] C. Dvorkin, W. Hu, and K. M. Smith, Phys. Rev. D79, 107302 (2009).
  * [33] C. Dvorkin and K. M. Smith, Phys. Rev. D79, 043003 (2009).
  * [34] A. Roy, A. Lapi, D. Spergel, and C. Baccigalupi, J. Cosm. Astropart. Phys.2018, 014 (2018).
  * [35] C. Feng and G. Holder, Phys. Rev. D99, 123502 (2019).
  * [36] T. Namikawa et al., Phys. Rev. D104, 063514 (2021).
  * [37] A. Roy, A. van Engelen, V. Gluscevic, and N. Battaglia, arXiv e-prints arXiv:2201.05076 (2022).
  * [38] B. A. Benson et al., in Millimeter, Submillimeter, and Far-Infrared Detectors and Instrumentation for Astronomy VII, Vol. 9153 of Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, edited by W. S. Holland and J. Zmuidzinas (PUBLISHER, ADDRESS, 2014), p. 91531P.
  * [39] S. W. Henderson et al., Journal of Low Temperature Physics 184, 772 (2016).
  * [40] P. Ade et al., J. Cosm. Astropart. Phys.2019, 056 (2019).
  * [41] The CMB-HD Collaboration et al., arXiv e-prints arXiv:2203.05728 (2022).
  * [42] K. N. Abazajian et al., arXiv e-prints arXiv:1610.02743 (2016).
  * [43] Planck Collaboration I, Astron. Astrophys.641, A1 (2020).
  * [44] C. L. Bennett et al., Astrophys. J. Suppl. Ser.208, 20 (2013).
  * [45] E. F. Schlafly, A. M. Meisner, and G. M. Green, Astrophys. J. Suppl. Ser.240, 30 (2019).
  * [46] A. Krolewski, S. Ferraro, E. F. Schlafly, and M. White, J. Cosm. Astropart. Phys.2020, 047 (2020).
  * [47] A. Krolewski, S. Ferraro, and M. White, J. Cosm. Astropart. Phys.2021, 028 (2021).
  * [48] E. L. Wright et al., Astron. J.140, 1868 (2010).
  * [49] A. M. Meisner, D. Lang, and D. J. Schlegel, Research Notes of the American Astronomical Society 2, 1 (2018).
  * [50] N. Aghanim, S. Majumdar, and J. Silk, Reports on Progress in Physics 71, 066902 (2008).
  * [51] T. Schutt et al., (in prep.).
  * [52] W. R. Coulton et al., arXiv e-prints arXiv:2307.01258 (2023).
  * [53] Planck Collaboration et al., Astron. Astrophys.643, A42 (2020).
  * [54] R. Dünner et al., Astrophys. J. 762, 10 (2013).
  * [55] J. Dunkley et al., Astrophys. J. 739, 52 (2011).
  * [56] S. Das et al., Astrophys. J. 729, 62 (2011).
  * [57] J. L. Sievers et al., J. Cosm. Astropart. Phys.2013, 060 (2013).
  * [58] M. B. Gralla et al., Astrophys. J. 893, 104 (2020).
  * [59] T. Louis et al., J. Cosm. Astropart. Phys.2017, 031 (2017).
  * [60] S. Naess et al., J. Cosm. Astropart. Phys.2014, 007 (2014).
  * [61] S. Aiola et al., J. Cosm. Astropart. Phys.2020, 047 (2020).
  * [62] S. Naess et al., J. Cosm. Astropart. Phys.2020, 046 (2020).
  * [63] T. W. Morris et al., Phys. Rev. D105, 042004 (2022).
  * [64] C. L. Bennett et al., Astrophys. J. Lett.396, L7 (1992).
  * [65] C. L. Bennett et al., Astrophys. J. Suppl. Ser.148, 97 (2003).
  * [66] X. Chen and E. L. Wright, Astrophys. J. 694, 222 (2009).
  * [67] M. Remazeilles, J. Delabrouille, and J.-F. Cardoso, Mon. Not. R. Astron. Soc.410, 2481 (2011).
  * [68] A. Kusiak, B. Bolliet, A. Krolewski, and J. C. Hill, Phys. Rev. D106, 123517 (2022).
  * [69] A. Kusiak, K. M. Surrao, and J. C. Hill, arXiv e-prints arXiv:2303.08121 (2023).
  * [70] A. M. Meisner, D. Lang, and D. J. Schlegel, Astron. J.153, 38 (2017).
  * [71] A. M. Meisner, D. Lang, and D. J. Schlegel, Astron. J.154, 161 (2017).
  * [72] Gaia Collaboration et al., Astron. Astrophys.616, A1 (2018).
  * [73] H. Akaike, in Selected Papers of Hirotugu Akaike, edited by E. Parzen, K. Tanabe, and G. Kitagawa (Springer New York, New York, NY, 1998), pp. 199–213.
  * [74] H. Akaike, IEEE Transactions on Automatic Control 19, 716 (1974).
  * [75] D. Hahs-Vaughn and R. Lomax, Statistical Concepts - A Second Course (Routledge, London, 2012).
  * [76] M. Mallaby-Kay et al., Phys. Rev. D108, 023516 (2023).
  * [77] D. Nelson et al., Astronomy and Computing 13, 12 (2015).
  * [78] D. Nelson et al., Computational Astrophysics and Cosmology 6, 2 (2019).
  * [79] R. Davé et al., Mon. Not. R. Astron. Soc.486, 2827 (2019).
  * [80] L. Thiele et al., Machine Learning: Science and Technology 3, 035002 (2022).
  * [81] B. K. K. Lee, W. R. Coulton, L. Thiele, and S. Ho, Mon. Not. R. Astron. Soc.517, 420 (2022).
  * [82] J. F. Navarro, C. S. Frenk, and S. D. M. White, Astrophys. J. 462, 563 (1996).
  * [83] K. S. Dawson et al., Astron. J.145, 10 (2013).
  * [84] B. A. Reid et al., Mon. Not. R. Astron. Soc.444, 476 (2014).
  * [85] B. Hadzhiyska et al., arXiv e-prints arXiv:2305.00992 (2023).
  * [86] M. J. Mortonson and W. Hu, Astrophys. J. 657, 1 (2007).
  * [87] O. Doré et al., Phys. Rev. D76, 043002 (2007).
  * [88] M. Su et al., arXiv e-prints arXiv:1106.4313 (2011).
  * [89] A. Roy et al., J. Cosm. Astropart. Phys.2021, 003 (2021).
  * [90] F. Bianchini and M. Millea, Phys. Rev. D107, 043521 (2023).
  * [91] T. Namikawa, D. Hanson, and R. Takahashi, Mon. Not. R. Astron. Soc.431, 609 (2013).
  * [92] M. S. Madhavacheril and J. C. Hill, Phys. Rev. D98, 023534 (2018).
  * [93] E. Schaan and S. Ferraro, Phys. Rev. Lett. 122, 181301 (2019).
  * [94] N. Sailer et al., Phys. Rev. D104, 123514 (2021).
  * [95] O. Darwish et al., arXiv e-prints arXiv:2111.00462 (2021).
  * [96] N. MacCrann et al., arXiv e-prints arXiv:2304.05196 (2023).
  * [97] B. Horowitz, S. Ferraro, and B. D. Sherwin, Mon. Not. R. Astron. Soc.485, 3919 (2019).
  * [98] B. Hadzhiyska, B. D. Sherwin, M. Madhavacheril, and S. Ferraro, Phys. Rev. D100, 023547 (2019).
  * [99] F. J. Qu, A. Challinor, and B. D. Sherwin, Phys. Rev. D108, 063518 (2023).
  * [100] M. Hilton et al., Astrophys. J. Suppl. Ser.253, 3 (2021).
  * [101] F. J. Qu et al., arXiv e-prints arXiv:2304.05202 (2023).
  * [102] Our simulations showed that this source of bias is generally negligible, however we retain this cut out of an abundance of caution.
  * [103] G. Stein et al., Journal of Cosmology and Astroparticle Physics 2020, 012–012 (2020).
  * [104] N. Sehgal et al., Astrophys. J. 709, 920 (2010).
  * [105] Z. Atkins et al., arXiv e-prints arXiv:2303.04180 (2023).
  * [106] R. Datta et al., Journal of Low Temperature Physics 184, 568 (2016).
  * [107] M. B. Gralla et al., Astrophys. J. 893, 104 (2020).
  * [108] M. Hilton et al., Astrophys. J. Suppl. Ser.235, 20 (2018).
  * [109] W. Coulton, A. Miranthis, and A. Challinor, Mon. Not. R. Astron. Soc.523, 825 (2023).
  * [110] B. Stölzner, A. Cuoco, J. Lesgourgues, and M. Bilicki, Phys. Rev. D97, 063506 (2018).
  * [111] L. Amendola et al., Living Reviews in Relativity 21, 2 (2018).
  * [112] Planck Collaboration XXI, Astron. Astrophys.594, A21 (2016).

##  Appendix A Supplementary Material

##  Appendix B I. Overview of Patchy Screening Estimators

In this Section we first review the standard estimator as derived in Ref. [33]
before presenting our new method. We refer the reader to Ref.[51] for a
detailed analysis of the relative merits of the different estimators.

In the screening effect, free electrons along a given line of sight (LOS)
Thomson scatter a fraction $1-e^{-\tau}$ of the CMB photons from this LOS away
from the observer. At the same time, Thomson scattering by the same free
electrons deflects other CMB photons (which would not otherwise reach the
observer) into the observer’s LOS. Because Thomson scattering is a time-
reversible process, the probability for a LOS photon to be deflected out of
the LOS is the same as the probability for any photon to be deflected into the
LOS. Thus screening is simply replacing a fraction $1-e^{-\tau}$ of LOS photon
flux at $\tilde{T}(\hat{n})$ with an equal fraction $1-e^{-\tau}$ of all other
photons incident on the electron, whose average temperature is
$\bar{T}_{\text{CMB}}$[33]:

$T(\hat{n})=\tilde{T}(\hat{n})e^{-\tau}+\bar{T}_{\text{CMB}}\left[1-e^{-\tau(\hat{n})}\right].$
(2)

This gives the observed temperature fluctuation:

$\displaystyle\Delta T(\hat{n})$ $\displaystyle\equiv
T(\hat{n})-\bar{T}_{\text{CMB}}$ (3)
$\displaystyle=\Delta\tilde{T}(\hat{n})e^{-\tau(\hat{n})}\quad\text{with
}\Delta\tilde{T}(\hat{n})\equiv\tilde{T}(\hat{n})-\bar{T}_{\text{CMB}},$
$\displaystyle=\Delta\tilde{T}(\hat{n})e^{-\bar{\tau}}e^{-\delta\tau(\hat{n})}\quad\text{with
}\tau(\hat{n})=\bar{\tau}+\delta\tau(\hat{n}),$ $\displaystyle=\Delta
T^{\text{primary}}(\hat{n})e^{-\delta\tau(\hat{n})}\quad\text{with }\Delta
T^{\text{primary}}(\hat{n})\equiv\Delta\tilde{T}(\hat{n})e^{-\bar{\tau}},$
$\displaystyle\simeq\Delta
T^{\text{primary}}(\hat{n})\left[1-\delta\tau(\hat{n})\right].$

Here we have introduced the average optical depth $\bar{\tau}$ to Thomson
scattering, through which the primary fluctuations are observed ($\Delta
T^{\text{primary}}(\hat{n})$ versus the temperature fluctuations at the
surface of last scattering $\Delta\tilde{T}(\hat{n})$). While this mean
optical depth is not negligible ($\bar{\tau}\sim 5-6\%$ from Planck [22, 23,
24, 25]) the fluctuations $\delta\tau(\hat{n})$ in the optical depth are much
smaller, of order $10^{-5}-10^{-3}$ for galaxy/cluster halos, justifying the
linear expansion in the last line above.

Screening of the CMB also occurs in polarization. The equations above remain
exact when replacing the temperature $T$ by the Stokes parameters $Q$ and $U$.
In terms of $E$ and $B$ instead, the multiplication of $Q$ and $U$ with a
scalar function $1-\delta\tau(\hat{n})$ mixes $E$ and $B$ modes, thus leading
to the generation of $B$-modes from $E$-modes [86, 87, 33, 32, 88, 34, 89].

###  B.1 The quadratic estimator

As described in Ref. [33], we can construct a quadratic estimator for
$\delta\tau$ as

$\displaystyle\delta\tilde{\tau}(\mathbf{L})=\int\frac{\mathrm{d}^{2}\bm{\ell}_{1}}{\left(2\pi\right)^{2}}\frac{\mathrm{d}^{2}\bm{\ell}_{2}}{\left(2\pi\right)^{2}}g(\bm{\ell}_{1},\bm{\ell}_{2})T(\bm{\ell}_{1})T(\bm{\ell}_{2})(2\pi)^{2}\delta^{(2)}(\mathbf{L}-\bm{\ell}_{1}-\bm{\ell}_{2})$
(4)

where $T(\bm{\ell})$ are the 2D Fourier modes of the observed CMB temperature
fluctuation map, as the intuition is clearer in the flat-sky approximation,
and $g(\bm{\ell}_{1},\bm{\ell}_{2})$ are a set of weights. The weights are
found by searching for the minimum variance unbiased estimator giving

$\displaystyle\delta\tilde{\tau}(\mathbf{L})=-\frac{1}{N_{L}}\int\frac{\mathrm{d}^{2}\bm{\ell}_{1}}{\left(2\pi\right)^{2}}\frac{\mathrm{d}^{2}\bm{\ell}_{2}}{\left(2\pi\right)^{2}}(2\pi)^{2}\delta^{(2)}(\mathbf{L}-\bm{\ell}_{1}-\bm{\ell}_{2})(C^{TT}_{\bm{\ell}_{1}})\frac{T(\bm{\ell}_{1})T(\bm{\ell}_{2})}{C^{\mathrm{tot}}_{\bm{\ell}_{1}}C^{\mathrm{tot}}_{\bm{\ell}_{2}}},$
(5)

where $C^{TT}_{\ell}$ is the power spectrum of the
$\tilde{T}(\mathbf{n})e^{-\bar{\tau}}$ field, i.e. the standard CMB power
spectrum, $C^{\mathrm{tot}}_{\ell}$ is the total power spectrum including
foregrounds and instrumental noise and $N_{L}$ is a normalization constant
given by

$\displaystyle
N_{L}=2\int\frac{\mathrm{d}^{2}\bm{\ell}_{1}}{\left(2\pi\right)^{2}}\frac{\mathrm{d}^{2}\bm{\ell}_{2}}{\left(2\pi\right)^{2}}(2\pi)^{2}\delta^{(2)}(\mathbf{L}-\bm{\ell}_{1}-\bm{\ell}_{2})\frac{\left.C^{TT}_{\bm{\ell}_{1}}\right.^{2}}{C^{\mathrm{tot}}_{\bm{\ell}_{1}}C^{\mathrm{tot}}_{\bm{\ell}_{2}}}.$
(6)

###  B.2 General small-scale $\tau$ estimator

For small scale reconstructions, $L\gtrsim$ 2000, on data sets with noise and
foregrounds, this can be approximated as

$\displaystyle\delta\tilde{\tau}(\mathbf{L})\approx$
$\displaystyle-\int\mathrm{d}^{2}\mathbf{n}e^{i\mathbf{n}\cdot\mathbf{L}}T_{\mathrm{large-
scales}}(\mathbf{n})\left[C^{-1}T_{\mathrm{small-
scales}}\right](\mathbf{n})C^{\mathrm{tot}}_{L}\left[\int_{\ell_{\mathrm{small-
scales}}}\frac{\mathrm{d}^{2}\bm{\ell}}{(2\pi)^{2}}C^{TT}_{\ell}\right]^{-1},$
(7)

where $T_{\mathrm{large-scales}}(\mathbf{n})$ is a map that contains only the
large-scale modes of the map, $\left[C^{-1}T_{\mathrm{small-
scales}}\right](\mathbf{n})$ is an inverse variance filtered map of small-
scale modes and $\ell_{\mathrm{small-scales}}$ are the modes that contribute
to the small-scale map. This estimator arises from two parts: first the filter
$C^{TT}_{\ell}/C^{\mathrm{tot}}_{\ell}$ in Eq. 5 selects only modes where the
CMB dominates. This filtering results in a map containing only large-scale
modes. Second, the Dirac delta function in Eq. 5 forces the second part of the
estimator to only select small-scale modes. Finally, we have approximated
$C^{tot}_{\bm{\ell}+\mathbf{L}}\approx C^{tot}_{\mathbf{L}}$, which is valid
on small scales where the power spectrum varies slowly. For the small-scale
regime, these approximations will not significantly lose information compared
to the full quadratic estimator or Bayesian estimators [90].

The key challenge on small scales is separating out this signal from other
non-Gaussian sky components such as gravitational lensing, the thermal and
kinetic Sunyaev-Zeldovich effects and the cosmic infrared background (CIB).
These signals all produce quadratic couplings that will contribute to this
estimator and can bias measurements of $\delta\tau$. In a manner analogous to
lensing these can be removed by bias hardening the estimator [91, 36] and by
using component separation techniques to remove signals with different
spectral signatures (e.g. the tSZ and CIB)[69]. This approach will suppress
these biases at the cost of reduced signal to noise. The difficulty for patchy
screening is that the biases are much larger than the signal and so need to be
suppressed to a very high degree. Fully removing similar biases to lensing
estimators is very challenging [92, 93, 94, 95, 96] and it has yet to be
demonstrated that such bias mitigation methods are sufficient for patchy
screening analyses. This motivates the consideration of alternative estimators
that are naturally robust to foreground biases.

By examining the structure of the small scale patchy screening effect

$\displaystyle T(\mathbf{n})\approx-T_{\mathrm{large-
scale}}(\mathbf{n})\delta\tau(\mathbf{n}),$ (8)

it can be seen that an estimator for patchy screening can be obtained by
replacing the $T_{\mathrm{large-scale}}$ term in Eq. 7 with just the sign of
the large scale mode and changing the normalization of the estimator. This
gives the following estimator

$\displaystyle\delta\tilde{\tau}(\mathbf{L})\approx$
$\displaystyle-\int\mathrm{d}^{2}\mathbf{n}e^{i\mathbf{n}\cdot\mathbf{L}}\,\mathrm{Sign}\left[T_{\mathrm{large-
scales}}(\mathbf{n})\right]\left[C^{-1}T_{\mathrm{small-
scales}}\right](\mathbf{n})C^{\mathrm{tot}}_{L}\left[\langle|T_{\mathrm{large-
scales}}|\rangle\right]^{-1},$ (9)

where $\langle|T_{\mathrm{large-scales}}|\rangle$ is the expectation of the
modulus of the large scale field. It is straightforward to show that this
estimator provides an unbiased estimator of the filtered $\tau$ field. Modes
with wavelengths significantly larger than those in the small-scale field will
not be recovered. Further, for small scale reconstructions the variance of
this new estimator is only a factor of $\langle T_{\mathrm{large-
scales}}^{2}\rangle/\langle|T_{\mathrm{large-scales}}|\rangle^{2}=\pi/2$
larger. This result is obtained by evaluating the two expectations using a
Gaussian distribution and using the relation

$\displaystyle\langle|T|\rangle=\int\limits_{-\infty}^{\infty}\mathrm{d}T\,|T|p(T)=2\int\limits_{0}^{\infty}\mathrm{d}T\,Tp(T),$
(10)

where $p(T)$ is the probability distribution of the temperature anisotropies.
At the cost of this extra noise, this estimator suppresses all foreground
biases. To understand why, note that the large scales of the CMB are dominated
by the primary CMB anisotropies. By taking the sign of this fluctuation we
isolate the primary CMB from the other contributions to the large scale CMB.
This simple method removes almost all possible biases: biases such as the tSZ,
CIB and kSZ are not correlated with the primary CMB so average to zero. A
subtle bias could arise from correlations with the integrated Sachs-Wolfe
(ISW) effect, which give a non-trivial contribution to the large scale modes.
These can be removed by also removing the largest scale CMB modes
($\ell\lesssim 20$) or subtracting a map of the ISW field from the large scale
leg. The latter approach is discussed in Section IV of the Supplementary
Material.

Lensing is correlated with the large scale primary CMB; however this term also
does not contribute to the estimator. To see this we first examine the effect
of lensing on the small scales. This is given by

$\displaystyle T(\mathbf{n})$
$\displaystyle=\tilde{T}\left(\mathbf{n}+\nabla\phi(\mathbf{n})\right)\approx\tilde{T}\left(\mathbf{n}\right)+\nabla\phi(\mathbf{n})\nabla\tilde{T}\left(\mathbf{n}\right),$
(11)

where $\phi(\mathbf{n})$ is the lensing potential. As in the patchy screening
case, the small scale lensing is given by the modulation of a large scale
primary CMB term ($\nabla\tilde{T}$) by the small scale term (for lensing
$\nabla\phi$). However lensing depends on the gradient of the large-scale CMB
anisotropies, which is not correlated with the sign of the CMB anisotropy.
Thus lensing also does not contribute to this estimator.

There are similarities between our estimator and the ‘gradient inversion’
lensing estimators [e.g. 97, 98]. However the gradient-inversion estimator
generally increases the SNR of lensing reconstruction when lensing is
dominant. This regime generally does not exist for the patchy screening
estimator so we do not consider further dividing out the specific large scale
mode, which would lead to the equivalent of the ’gradient inversion’ estimator
for patchy screening. Rather our estimator can be thought of as an analogy to
the shear-only lensing estimators [93, 99]. We consider only part of the
signal at the cost of signal to noise, but gain robustness to foreground
emission. A detailed comparison of these different estimators is available in
Ref. [51].

###  B.3 Implementation of the cross correlation estimator

To apply Eq. 9, we first need to create two filtered maps: one of the large-
scale CMB anisotropies, and one of the small-scale CMB. For the low-pass maps
we apply the following harmonic space filter, $f^{\mathrm{low-pass}}_{\ell}$,

$\displaystyle f^{\mathrm{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}$ (12)

For the high-pass filter we use the following filter, $f^{\mathrm{high-
pass}}_{\ell}$,

$\displaystyle f^{\mathrm{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}$ (13)

The structure of the high- and low-pass filters was chosen to satisfy the
following criteria. First, the low-pass map needs to contain enough modes to
preserve the large scale sign. The fact that the CMB power spectrum is red
means that this is satisfied providing the transition is above $\ell\sim
1000$. The location of the high-pass filter was chosen such that the power
spectrum of the maps is approximately white (i.e. flat). When the power
spectrum is white, Eq. 9 further simplifies as the harmonic-space, inverse-
variance filter is just a constant and factorizes out. The smooth tapering of
the filters prevents excess ringing from the harmonic space filtering. A
separation between the two filters is used to prevent correlations between the
high- and low-pass filtered maps. Moving the location, and width, of these
filters by $\pm 250$ was found to have minimal effect on the analysis. It is
left to future work to optimize location of the filters, and widths of the
transitions, which could be further optimized to maximize signal to noise.

Note that we do not deconvolve the beam from the scale-scale temperature map
and instead include the beam and the impact of the filtering in our theory
models. However, we do deconvolve the pixel window function, which would be
difficult to account for in the theory model.

After the cuts described below, we extract a 20$\prime$ x 20$\prime$ cutout
from the small scale map at the location of each unWISE galaxy. To account for
sky curvature we use the tangent plane projection and interpolation methods
described in Ref. [15]. We stack these cutouts, weighted by the sign of the
large scale mode, as

$\displaystyle\widehat{\delta\tau}^{\mathrm{stacked}}(\mathbf{n})=$
$\displaystyle\frac{-\sum\mathrm{Sign}[T^{\mathrm{large-
scale}}(\mathbf{n_{o}}+\mathbf{n})]T^{\mathrm{small-
scale}}(\mathbf{n_{o}}+\mathbf{n})}{N_{\mathrm{objects}}\langle|T^{\mathrm{large-
scale}}|\rangle},$ (14)

where the sum is over objects at locations $\mathbf{n}_{o}$ and
$N_{\mathrm{objects}}$ is the number of objects in the stack. We stack on tens
of millions of galaxies for the blue and green samples; see Table 1 for source
number densities. The sign of the large-scale field is consistent over $\sim
1/2\,\mathrm{deg}^{2}$ patches. A 1D azimuthally averaged profile is computed
from the stacks.

The non-linearity introduced by the sign operation suppresses foreground
biases in the large scale map as the sign is determined by the primary CMB
anisotropies. The immunity to foregrounds provided by the sign operation means
we can include small scale temperature measurements, which are dominated by
foregrounds, without becoming biased.

We do not stack on all the objects in the input catalog. Instead we perform a
series of cuts to ensure that our results are robust. Whilst the foregrounds
will average to zero over many objects, due to the sign weighting, for a
finite number there can be a non-trivial residual [15]. To minimize the impact
of especially bright objects skewing the average, we exclude regions from our
maps where there are known point sources or clusters. Specifically we avoid
stacking on objects that are within 6′ of a detected SZ cluster (those in the
updated Ref. [100] cluster catalog ), 6′ of a subtracted point source (those
detected in any of the input maps with SNR$>5$), 10$\prime$ of the edge of the
map and 12′ of sources that are inpainted. The latter are a set of especially
bright objects that are inpainted during the construction of the NILC maps and
are described in Ref. [101].

Next, we ensure that in the stack there are exactly equal numbers of
positively and negatively weighted regions. This helps further suppress the
foreground emission and is implemented by selecting a random subsample. A
heuristic of this can be understood by splitting the sky into two regions:
areas where $\mathrm{Sign}[T^{\mathrm{large-scale}}(\mathbf{n_{o}})]$ is
positive and areas where it is negative. Within each region foregrounds will
not average down, and so contribute a bias given by the size of the average
source present in the map, e.g., the average amplitude of a dusty galaxy or
thermal SZ cluster. The average is obtained from many objects ($\sim 10^{6}$)
and so will be similar in the positive and negative regions of the sky. Thus,
the degree to which these foregrounds cancel in the final estimate is set by
the difference in number of sources in the positive and negative regions of
the sky. The large scale CMB has a large coherence scale ($\sim 1^{\circ}$) so
in any single sky realization there will be O(10000) independent positive and
negative regions. The Poissionian-like scatter means we expect $\sim$1 $\%$
more sources in one region than the other. Without explicitly balancing the
number of sources in each region, the degree of cancellation will leave a
$\sim$1 $\%$ residual of the mean foreground signal. If we instead manually
force the same number of sources in each region, the cancellation will only
depend on the difference in the average residual in each region. This will be
approximately $1/\sqrt{N_{\mathrm{objects}}}$ and in our case this is expected
to reduce the residual by more than an order of magnitude.

Finally we can apply cuts based on the amplitude of the large-scale field. The
argument that the sign operation removes contaminants works only when the
large scale is dominated by the CMB. Whilst this is generally true, there are
regions of the sky where the primary CMB anisotropies are small and these
patches could be dominated by foreground emission. To mitigate this as a
source of potential contamination we exclude the regions of the sky where the
large-scale temperature fluctuation is $\leq 40\,\mu$K, where tSZ and CIB
emissions could potentially change the sign of the large-scale temperature.
This threshold was based on the typical large scale variance of the tSZ and
CIB in simulated maps.111Our simulations showed that this source of bias is
generally negligible, however we retain this cut out of an abundance of
caution.

##  Appendix C II. Validation with Simulations

To validate this approach we apply the method to non-Gaussian CMB sky
simulations. The goal is to verify both that the other sky signals do not bias
the estimator and that the estimator is able to recover the input signal.
These tests are structured to approximate the sky as observed by ACT and
Planck satellite so that they serve as a validation of the data analysis
presented in the main text.

###  C.1 Simulation Properties

We consider two non-Gaussian sky simulations: the Websky [103] and the Sehgal
et al. [104] simulations. Each simulation contains a range of correlated and
non-Gaussian CMB secondary anisotropies including: the tSZ, kSZ, CIB, lensing
and radio galaxies. Both simulations are dark matter only simulations with the
secondary anisotropies then ‘painted’ on in post processing. As the two
simulations use distinct methods and models for this painting procedure, they
provide some measure of the theoretical uncertainty on the properties of the
CMB secondaries. We use the rescaled Sehgal et al. simulations presented in
Ref. [40] that have been adjusted to be consistent with recent measurements.

Using these simulations of the sky, we construct simulated observations that
match the properties of the ACT and Planck data sets used in this work (see
the main text for further details). These simulations match the instrumental
properties including map noise, beams, and observation frequencies. We refer
the reader to Ref. [105] and Appendix B of Ref. [52] for a thorough
description of these simulations. We use the Needlet Internal Linear
Combination (NILC) component separation pipeline, described in Ref. [52], to
combine simulations of the individual observations into a map of the CMB
anisotropies. Through this procedure we obtain maps with the appropriate
levels of residual foregrounds that are present in the observations.

The Websky and Sehgal et al. simulations do not include any patchy screening
signal. To validate that our method provides an unbiased measurement of non-
zero signals, we created an additional version of the simulations that contain
an approximate patchy screening signal. We first construct an approximate
$\tau$ map, by rescaling the Compton-y map from each simulation by a constant
factor to “convert” from Compton-y to optical depth. This factor is an
approximation to the mean gas temperature of the halos and is chosen to give
an average $\tau$ similar to the expected value ($\sim 5\times 10^{-4}$). With
this approximate $\tau$ map we create the patchy screening effect by using it
to screen the primary CMB anisotropies, via Eq. 1.

The catalogs used in the stacking analysis are the halo catalogs of the two
simulations. For the Websky simulations we use all halos with mass
$M_{h}>1.8\times 10^{13}$ and in Sehgal et al. simulations we stack on all the
available halos. Like the processing of the data, we avoid stacking on objects
that are close to point sources, clusters or the edge of the mask. For the
data the sources are identified in the maps using point source and cluster
finders [e.g., 106, 107, 108, 100], however for the simulations we instead use
the true source catalogs. The cut thresholds are set to approximately match
the point source and cluster thresholds of the ACT data: CIB sources with
amplitudes $\gtrsim 350\mu$ K ($>30$ mJy) at $217\,$GHz, radio sources with
amplitudes $\gtrsim 150\mu$ K ($>9$ mJy) at $90\,$ GHz and clusters with
$M_{200}\gtrsim 5\times 10^{14}$ M⊙/h.

###  C.2 Analysis of the Simulations

(a) Websky simulations

(b) Sehgal et al. simulations

Figure 5: An application of our stacked estimator to mock ACT & Planck
observations constructed from the Websky and Sehgal et al. simulations. The
cross and circle points correspond to the two realizations. We apply the
estimator to simulations that contain foregrounds but no patchy screening
signal (blue lines). These are consistent with zero (see Table 3) and
demonstrate that our method is immune to foreground biases. We also show that
the estimator is unbiased, as the measurements on simulations with a $\tau$
signal (orange) match the input signal when filtered in the same manner
(green).

In Fig. 5a and Fig. 5b we apply the method to the Websky and Sehgal et al.
simulations. As can be seen when applied to the default simulations, those
without any patchy screening signal, the estimator measurements are consistent
with zero. This demonstrates that our method is not biased by the foregrounds
present in these simulations. When applied to the simulations with a mock
screening signal, we recover measurements that are consistent with the input
$\tau$ signal. The oscillatory structure seen in the signal arises from the
high-pass filtering performed on the small-scale CMB map.

These two tests provide evidence that our pipeline is able to provide robust
and unbiased measurements of a patchy screening signal.

##  Appendix D III. Uncertainty quantification

To compute the errors we generate a set of Gaussian simulations of the sky
that match the properties of the sky. As the signal is highly subdominant to
all the sky signals, we can simply treat the NILC map as a map of “noise.”
This is used as the input to the tiled noise model, described in Ref. [105].
With this formalism we can accurately capture the variation in the “noise”
properties in both scale and location across the map. From these we generate
an ensemble of 150 simulations. For each simulation we perform the same
filtering as the data and then apply the stacked sign estimator at the
locations of the unWISE galaxies. The weighting is given by the sign of a
simulated large-scale CMB map. This is repeated over the ensemble of
simulations to obtain an estimate of the covariance matrix.

As a cross check of this approach, we compute an approximate covariance matrix
from the data themselves by computing the covariance matrix of the individual
$\tau$ measurements from each galaxy. If each galaxy’s profile were
independent the covariance matrix of the stack would be
$1/\sqrt{N_{\mathrm{objects}}}$ times the single galaxy covariance matrix. In
practice the galaxies are not completely independent and thus this covariance
matrix would slightly underestimate the true variance. We found that this
approximate covariance matrix was in reasonable agreement with our simulated
covariance matrix ($\sim 15\%$ agreement between the diagonal elements of the
two matrices), providing a mild validation of our simulations.

##  Appendix E IV. Robustness Tests

Figure 6: A summary of the null and consistency tests performed to validate our result. The two null tests, stacking on random location and randomizing the sign used in the stack, are both consistent with zero and provide a data based check of our method (see Table 3). The two foreground tests, which are labeled “Deproj CIB and tSZ” and “ISW removed”, explicitly remove contamination from the tSZ, CIB and ISW effect. These demonstrate that our result likely contains no significant biases from these effects. The final tests, using only data from 93GHz or 148GHz, are consistent with anisotropies that have a blackbody spectrum, as expected for the patchy screening effect.  |  |  | Equivalent
---|---|---|---
| $\chi^{2}$ |  $p$-value | Gaussian $\sigma$
Random positions | 12 | 0.19 | 0.89
Random signs | 13 | 0.16 | 0.98
Table 3: The $\chi^{2}$, $p$-values and equivalent Gaussian $\sigma$ for two
of the null tests.

To validate our results we performed a series of data based robustness and
null tests.

As a test of potential foreground biases, and of the computation of the
covariance matrices, we performed two null tests. First we repeated our
stacking analysis on random locations. Our estimator, Section B.3, is not
sensitive to the mean optical depth and so it is expected that if we stack on
random locations we should see a signal consistent with zero. The number of
locations is the same as the number in the original catalog, for which we use
the largest sample, the blue sample. As seen in Fig. 6, when random locations
are used, we find no signal. The $\chi^{2}$ of this null test compared to
zero, reported in Table 3, is consistent with no signal. This tests both our
estimator and the covariance matrices.

Second we repeat our stacking analysis, but randomize the sign used to weight
each profile. This is a test of residual foreground components. Foregrounds
are expected to average to zero in our standard estimator. However, if there
are any anomalously large values they may average down very slowly and bias
our result. This “random sign” test should identify such cases as randomizing
the sign will remove the patchy screening signal. As seen in Fig. 6 and Table
3, this null test is also consistent with no signal.

Next we perform tests aimed at identifying contamination in our measurement
from two known sources: the thermal Sunyaev-Zeldovich effect (tSZ) and the
cosmic infrared background (CIB). To test for biases from the tSZ and CIB, we
make a constrained NILC map that explicitly removes the tSZ and CIB effects as
described in Ref. [66, 67]. This map is referred to as the CIB and tSZ
deprojected map and we refer the reader to [52] for a detailed description of
our implementation. This operation results in a significant increase in the
noise in the component-separated CMB map. We could deproject the tSZ and CIB
separately. However as shown in [94, 109] removing only one component often
increases the other component and so it is hard to predict how our results
would change. However, on large scales the map is still dominated by CMB
modes. To remove any biases from our measurement we note that we only need to
ensure that the contamination is removed from one of the maps used in our
estimator. In this case we can use the tSZ and CIB deprojected map as the low-
pass filtered map for the sign operation. We make this choice as the scales
relevant for the large-scale sign are largely unaffected by the increased
noise in the deprojected map. The small scale noise is increased to such an
extent that the signal to noise is too low for a useful consistency test. The
resulting measurement, shown in Fig. 6, is consistent with our baseline
measurement. This provides further support that our measurement is not
contaminated by the tSZ or CIB.

We also check for contamination from the ISW effect, which leads to a bias
that is not removed from our estimator. This is because the ISW contributes a
significant amount of power to the CMB map and cannot be removed by
multifrequency cleaning. This means that there may be regions where the sign
is set by the ISW effect, not the CMB. The ISW effect traces the large scale
potential field and so is correlated with other cosmic structures, such as the
tSZ, CIB and unWISE galaxies, and thus could produce a signal in our analysis.
It is expected that any bias from the ISW effect will be small as, even with
optimal estimators, it is difficult to detect the ISW effect at high
significance [110, 111]. To test for this we subtract the Planck
collaboration’s map of the ISW effect (see Ref. [112] for details of the
construction of this map) from our component-separated CMB map. We then repeat
our stacking analysis. By explicitly removing the ISW effect we should
mitigate any potential bias. As shown in Fig. 6, our measurement is largely
unchanged. This suggests that our measurement is not dominated by bias from
the ISW effect. This null test has two limitations: first, it is difficult to
produce a map of the ISW anisotropies. The map produced in Ref. [112] is noisy
and not a perfect tracer of the true ISW effect. Thus our subtraction will not
perfectly and noiselessly remove any biases. The ISW effect is only likely to
bias our measurement on the largest scales, where it may impact the measured
sign, and these are the scales where the ISW effect is best measured. Second,
the ISW map has a slightly different mask compared to our fiducial analysis
(see Fig. 20. of Ref. [112]). This will introduce some noise and scatter in
our measurements. Given the very small impact of subtracting the ISW map, we
expect both of these effects to be minor. An alternative method would be to
filter out the largest scales. We do not consider that approach here as the
filtering adds large amounts of noise to the sign measurement.

The final test performed is aimed at characterizing the spectral response of
the observed signature signal. The frequency dependence of the patchy
screening signal is expected to be the same as the primary CMB anisotropies.
In the differential thermodynamic units used in this analysis, the
anisotropies are expected to be the same size at all frequencies. We test this
by generating maps that contain only observations at $93\,$GHz or $148\,$GHz.
These maps are then high-pass filtered and used in the estimator. These two
measurements will be consistent if the observed signal has no frequency
dependence. We choose these two frequencies as they have sufficiently low
noise and high resolution to obtain a high SNR measurement. Fig. 6 shows that
these two measurements are consistent and thereby provides evidence that our
observed signal has the expected frequency dependence.

##  Appendix F V. Gaussian Profile Parameter Posteriors

In the main text, we use a Gaussian model to aid the quantification of the
detection significance. We parameterized the Gaussian via the central
amplitude $\tau$ and the FWHM. The Gaussian is then convolved with the beam
and filtered in a manner identical to the data. In Fig. 7 we show the 2D
contour plot for the best fit parameters of the Gaussian model fit to the
green and blue samples. Whilst the exact degeneracy direction depends on the
parameterization (instead of central optical depth, we could have considered
the integrated optical depth), we find that both the amplitude and width are
constrained. Note that the filtering leads to a complex interplay between the
amplitude and width parameters. Filtering reduces the observed amplitude more
significantly for wider profiles - a profile that is wider than the filtering
scale is completely removed. This is discussed in detail in [51].

Figure 7: Constraints on the central optical depth and Full-Width-Half-Maximum
of the Gaussian profiles for the blue and green samples.