SPT-3G and SPTpol Collaboration

#  First Constraints on the Epoch of Reionization Using the non-Gaussianity of
the Kinematic Sunyaev-Zel’dovich Effect from the South Pole Telescope and
_Herschel_ -SPIRE Observations

S. Raghunathan  Center for AstroPhysical Surveys, National Center for
Supercomputing Applications, Urbana, IL, 61801, USA  P. A. R. Ade  School of
Physics and Astronomy, Cardiff University, Cardiff CF24 3YB, United Kingdom
A. J. Anderson  Fermi National Accelerator Laboratory, MS209, P.O. Box 500,
Batavia, IL, 60510, USA  Kavli Institute for Cosmological Physics, University
of Chicago, 5640 South Ellis Avenue, Chicago, IL, 60637, USA  Department of
Astronomy and Astrophysics, University of Chicago, 5640 South Ellis Avenue,
Chicago, IL, 60637, USA  B. Ansarinejad  School of Physics, University of
Melbourne, Parkville, VIC 3010, Australia  M. Archipley  Department of
Astronomy, University of Illinois Urbana-Champaign, 1002 West Green Street,
Urbana, IL, 61801, USA  Center for AstroPhysical Surveys, National Center for
Supercomputing Applications, Urbana, IL, 61801, USA  J. E. Austermann  NIST
Quantum Devices Group, 325 Broadway Mailcode 817.03, Boulder, CO, 80305, USA
Department of Physics, University of Colorado, Boulder, CO, 80309, USA  L.
Balkenhol  Institut d’Astrophysique de Paris, UMR 7095, CNRS & Sorbonne
Université, 98 bis boulevard Arago, 75014 Paris, France  J. A. Beall  NIST
Quantum Devices Group, 325 Broadway Mailcode 817.03, Boulder, CO, 80305, USA
K. Benabed  Institut d’Astrophysique de Paris, UMR 7095, CNRS & Sorbonne
Université, 98 bis boulevard Arago, 75014 Paris, France  A. N. Bender  High-
Energy Physics Division, Argonne National Laboratory, 9700 South Cass Avenue,
Lemont, IL, 60439, USA  Kavli Institute for Cosmological Physics, University
of Chicago, 5640 South Ellis Avenue, Chicago, IL, 60637, USA  Department of
Astronomy and Astrophysics, University of Chicago, 5640 South Ellis Avenue,
Chicago, IL, 60637, USA  B. A. Benson  Fermi National Accelerator Laboratory,
MS209, P.O. Box 500, Batavia, IL, 60510, USA  Kavli Institute for Cosmological
Physics, University of Chicago, 5640 South Ellis Avenue, Chicago, IL, 60637,
USA  Department of Astronomy and Astrophysics, University of Chicago, 5640
South Ellis Avenue, Chicago, IL, 60637, USA  F. Bianchini  Kavli Institute for
Particle Astrophysics and Cosmology, Stanford University, 452 Lomita Mall,
Stanford, CA, 94305, USA  Department of Physics, Stanford University, 382 Via
Pueblo Mall, Stanford, CA, 94305, USA  SLAC National Accelerator Laboratory,
2575 Sand Hill Road, Menlo Park, CA, 94025, USA  L. E. Bleem  High-Energy
Physics Division, Argonne National Laboratory, 9700 South Cass Avenue.,
Lemont, IL, 60439, USA  Kavli Institute for Cosmological Physics, University
of Chicago, 5640 South Ellis Avenue, Chicago, IL, 60637, USA  Department of
Astronomy and Astrophysics, University of Chicago, 5640 South Ellis Avenue,
Chicago, IL, 60637, USA  J. Bock  California Institute of Technology,
Pasadena, CA 91125, USA  Jet Propulsion Laboratory, California Institute of
Technology, Pasadena, CA 91109, USA  F. R. Bouchet  Institut d’Astrophysique
de Paris, UMR 7095, CNRS & Sorbonne Université, 98 bis boulevard Arago, 75014
Paris, France  L. Bryant  Enrico Fermi Institute, University of Chicago, 5640
South Ellis Avenue, Chicago, IL, 60637, USA  E. Camphuis  Institut
d’Astrophysique de Paris, UMR 7095, CNRS & Sorbonne Université, 98 bis
boulevard Arago, 75014 Paris, France  J. E. Carlstrom  Kavli Institute for
Cosmological Physics, University of Chicago, 5640 South Ellis Avenue, Chicago,
IL, 60637, USA  Enrico Fermi Institute, University of Chicago, 5640 South
Ellis Avenue, Chicago, IL, 60637, USA  Department of Physics, University of
Chicago, 5640 South Ellis Avenue, Chicago, IL, 60637, USA  High-Energy Physics
Division, Argonne National Laboratory, 9700 South Cass Avenue., Lemont, IL,
60439, USA  Department of Astronomy and Astrophysics, University of Chicago,
5640 South Ellis Avenue, Chicago, IL, 60637, USA  T. W. Cecil  High-Energy
Physics Division, Argonne National Laboratory, 9700 South Cass Avenue.,
Lemont, IL, 60439, USA  C. L. Chang  High-Energy Physics Division, Argonne
National Laboratory, 9700 South Cass Avenue., Lemont, IL, 60439, USA  Kavli
Institute for Cosmological Physics, University of Chicago, 5640 South Ellis
Avenue, Chicago, IL, 60637, USA  Department of Astronomy and Astrophysics,
University of Chicago, 5640 South Ellis Avenue, Chicago, IL, 60637, USA  P.
Chaubal  School of Physics, University of Melbourne, Parkville, VIC 3010,
Australia  H. C. Chiang  Department of Physics and McGill Space Institute,
McGill University, 3600 Rue University, Montreal, Quebec H3A 2T8, Canada
School of Mathematics, Statistics & Computer Science, University of KwaZulu-
Natal, Durban, South Africa  P. M. Chichura  Department of Physics, University
of Chicago, 5640 South Ellis Avenue, Chicago, IL, 60637, USA  Kavli Institute
for Cosmological Physics, University of Chicago, 5640 South Ellis Avenue,
Chicago, IL, 60637, USA  T.-L. Chou  Department of Physics, University of
Chicago, 5640 South Ellis Avenue, Chicago, IL, 60637, USA  Kavli Institute for
Cosmological Physics, University of Chicago, 5640 South Ellis Avenue, Chicago,
IL, 60637, USA  R. Citron  University of Chicago, 5640 South Ellis Avenue,
Chicago, IL, 60637, USA  A. Coerver  Department of Physics, University of
California, Berkeley, CA, 94720, USA  T. M. Crawford  Kavli Institute for
Cosmological Physics, University of Chicago, 5640 South Ellis Avenue, Chicago,
IL, 60637, USA  Department of Astronomy and Astrophysics, University of
Chicago, 5640 South Ellis Avenue, Chicago, IL, 60637, USA  A. T. Crites  Kavli
Institute for Cosmological Physics, University of Chicago, 5640 South Ellis
Avenue, Chicago, IL, 60637, USA  Department of Astronomy and Astrophysics,
University of Chicago, 5640 South Ellis Avenue, Chicago, IL, 60637, USA
Dunlap Institute for Astronomy & Astrophysics, University of Toronto, 50 St.
George Street, Toronto, ON, M5S 3H4, Canada  David A. Dunlap Department of
Astronomy & Astrophysics, University of Toronto, 50 St. George Street,
Toronto, ON, M5S 3H4, Canada  A. Cukierman  Kavli Institute for Particle
Astrophysics and Cosmology, Stanford University, 452 Lomita Mall, Stanford,
CA, 94305, USA  SLAC National Accelerator Laboratory, 2575 Sand Hill Road,
Menlo Park, CA, 94025, USA  Department of Physics, Stanford University, 382
Via Pueblo Mall, Stanford, CA, 94305, USA  C. Daley  Department of Astronomy,
University of Illinois Urbana-Champaign, 1002 West Green Street, Urbana, IL,
61801, USA  K. R. Dibert  Department of Astronomy and Astrophysics, University
of Chicago, 5640 South Ellis Avenue, Chicago, IL, 60637, USA  Kavli Institute
for Cosmological Physics, University of Chicago, 5640 South Ellis Avenue,
Chicago, IL, 60637, USA  M. A. Dobbs  Department of Physics and McGill Space
Institute, McGill University, 3600 Rue University, Montreal, Quebec H3A 2T8,
Canada  Canadian Institute for Advanced Research, CIFAR Program in Gravity and
the Extreme Universe, Toronto, ON, M5G 1Z8, Canada  A. Doussot  Institut
d’Astrophysique de Paris, UMR 7095, CNRS & Sorbonne Université, 98 bis
boulevard Arago, 75014 Paris, France  D. Dutcher  Joseph Henry Laboratories of
Physics, Jadwin Hall, Princeton University, Princeton, NJ 08544, USA  W.
Everett  Department of Astrophysical and Planetary Sciences, University of
Colorado, Boulder, CO, 80309, USA  C. Feng  Department of Physics, University
of Illinois Urbana-Champaign, 1110 West Green Street, Urbana, IL, 61801, USA
K. R. Ferguson  Department of Physics and Astronomy, University of California,
Los Angeles, CA, 90095, USA  K. Fichman  Department of Physics, University of
Chicago, 5640 South Ellis Avenue, Chicago, IL, 60637, USA  Kavli Institute for
Cosmological Physics, University of Chicago, 5640 South Ellis Avenue, Chicago,
IL, 60637, USA  A. Foster  Joseph Henry Laboratories of Physics, Jadwin Hall,
Princeton University, Princeton, NJ 08544, USA  S. Galli  Institut
d’Astrophysique de Paris, UMR 7095, CNRS & Sorbonne Université, 98 bis
boulevard Arago, 75014 Paris, France  J. Gallicchio  Kavli Institute for
Cosmological Physics, University of Chicago, 5640 South Ellis Avenue, Chicago,
IL, 60637, USA  Harvey Mudd College, 301 Platt Boulevard., Claremont, CA,
91711, USA  A. E. Gambrel  Kavli Institute for Cosmological Physics,
University of Chicago, 5640 South Ellis Avenue, Chicago, IL, 60637, USA  R. W.
Gardner  Enrico Fermi Institute, University of Chicago, 5640 South Ellis
Avenue, Chicago, IL, 60637, USA  F. Ge  Department of Physics & Astronomy,
University of California, One Shields Avenue, Davis, CA 95616, USA  E. M.
George  European Southern Observatory, Karl-Schwarzschild-Str. 2, 85748
Garching bei München, Germany  Department of Physics, University of
California, Berkeley, CA, 94720, USA  N. Goeckner-Wald  Department of Physics,
Stanford University, 382 Via Pueblo Mall, Stanford, CA, 94305, USA  Kavli
Institute for Particle Astrophysics and Cosmology, Stanford University, 452
Lomita Mall, Stanford, CA, 94305, USA  R. Gualtieri  Department of Physics and
Astronomy, Northwestern University, 633 Clark St, Evanston, IL, 60208, USA  F.
Guidi  Institut d’Astrophysique de Paris, UMR 7095, CNRS & Sorbonne
Université, 98 bis boulevard Arago, 75014 Paris, France  S. Guns  Department
of Physics, University of California, Berkeley, CA, 94720, USA  N. Gupta
CSIRO Space & Astronomy, PO Box 1130, Bentley WA 6102, Australia  T. de Haan
High Energy Accelerator Research Organization (KEK), Tsukuba, Ibaraki
305-0801, Japan  N. W. Halverson  CASA, Department of Astrophysical and
Planetary Sciences, University of Colorado, Boulder, CO, 80309, USA
Department of Physics, University of Colorado, Boulder, CO, 80309, USA  E.
Hivon  Institut d’Astrophysique de Paris, UMR 7095, CNRS & Sorbonne
Université, 98 bis boulevard Arago, 75014 Paris, France  G. P. Holder
Department of Physics, University of Illinois Urbana-Champaign, 1110 West
Green Street, Urbana, IL, 61801, USA  W. L. Holzapfel  Department of Physics,
University of California, Berkeley, CA, 94720, USA  J. C. Hood  Kavli
Institute for Cosmological Physics, University of Chicago, 5640 South Ellis
Avenue, Chicago, IL, 60637, USA  J. D. Hrubes  University of Chicago, 5640
South Ellis Avenue, Chicago, IL, 60637, USA  A. Hryciuk  Department of
Physics, University of Chicago, 5640 South Ellis Avenue, Chicago, IL, 60637,
USA  Kavli Institute for Cosmological Physics, University of Chicago, 5640
South Ellis Avenue, Chicago, IL, 60637, USA  N. Huang  Department of Physics,
University of California, Berkeley, CA, 94720, USA  J. Hubmayr  NIST Quantum
Devices Group, 325 Broadway Mailcode 817.03, Boulder, CO, 80305, USA  K. D.
Irwin  SLAC National Accelerator Laboratory, 2575 Sand Hill Road, Menlo Park,
CA, 94025, USA  Department of Physics, Stanford University, 382 Via Pueblo
Mall, Stanford, CA, 94305, USA  F. Kéruzoré  High-Energy Physics Division,
Argonne National Laboratory, 9700 South Cass Avenue., Lemont, IL, 60439, USA
L. Knox  Department of Physics & Astronomy, University of California, One
Shields Avenue, Davis, CA 95616, USA  M. Korman  Department of Physics, Case
Western Reserve University, Cleveland, OH, 44106, USA  K. Kornoelje
Department of Astronomy and Astrophysics, University of Chicago, 5640 South
Ellis Avenue, Chicago, IL, 60637, USA  Kavli Institute for Cosmological
Physics, University of Chicago, 5640 South Ellis Avenue, Chicago, IL, 60637,
USA  C.-L. Kuo  Kavli Institute for Particle Astrophysics and Cosmology,
Stanford University, 452 Lomita Mall, Stanford, CA, 94305, USA  Department of
Physics, Stanford University, 382 Via Pueblo Mall, Stanford, CA, 94305, USA
SLAC National Accelerator Laboratory, 2575 Sand Hill Road, Menlo Park, CA,
94025, USA  A. T. Lee  Department of Physics, University of California,
Berkeley, CA, 94720, USA  Physics Division, Lawrence Berkeley National
Laboratory, Berkeley, CA, 94720, USA  K. Levy  School of Physics, University
of Melbourne, Parkville, VIC 3010, Australia  D. Li  NIST Quantum Devices
Group, 325 Broadway Mailcode 817.03, Boulder, CO, 80305, USA  SLAC National
Accelerator Laboratory, 2575 Sand Hill Road, Menlo Park, CA, 94025, USA  A. E.
Lowitz  Kavli Institute for Cosmological Physics, University of Chicago, 5640
South Ellis Avenue, Chicago, IL, 60637, USA  C. Lu  Department of Physics,
University of Illinois Urbana-Champaign, 1110 West Green Street, Urbana, IL,
61801, USA  A. Maniyar  Kavli Institute for Particle Astrophysics and
Cosmology, Stanford University, 452 Lomita Mall, Stanford, CA, 94305, USA
Department of Physics, Stanford University, 382 Via Pueblo Mall, Stanford, CA,
94305, USA  SLAC National Accelerator Laboratory, 2575 Sand Hill Road, Menlo
Park, CA, 94025, USA  J. J. McMahon  Kavli Institute for Cosmological Physics,
University of Chicago, 5640 South Ellis Avenue, Chicago, IL, 60637, USA
Department of Physics, University of Chicago, 5640 South Ellis Avenue,
Chicago, IL, 60637, USA  Department of Astronomy and Astrophysics, University
of Chicago, 5640 South Ellis Avenue, Chicago, IL, 60637, USA  F. Menanteau
Department of Astronomy, University of Illinois Urbana-Champaign, 1002 West
Green Street, Urbana, IL, 61801, USA  Center for AstroPhysical Surveys,
National Center for Supercomputing Applications, Urbana, IL, 61801, USA  M.
Millea  Department of Physics, University of California, Berkeley, CA, 94720,
USA  J. Montgomery  Department of Physics and McGill Space Institute, McGill
University, 3600 Rue University, Montreal, Quebec H3A 2T8, Canada  C. Corbett
Moran  Jet Propulsion Laboratory, Pasadena, CA 91109, USA  Y. Nakato
Department of Physics, Stanford University, 382 Via Pueblo Mall, Stanford, CA,
94305, USA  T. Natoli  Kavli Institute for Cosmological Physics, University of
Chicago, 5640 South Ellis Avenue, Chicago, IL, 60637, USA  J. P. Nibarger
NIST Quantum Devices Group, 325 Broadway Mailcode 817.03, Boulder, CO, 80305,
USA  G. I. Noble  Dunlap Institute for Astronomy & Astrophysics, University of
Toronto, 50 St. George Street, Toronto, ON, M5S 3H4, Canada  David A. Dunlap
Department of Astronomy & Astrophysics, University of Toronto, 50 St. George
Street, Toronto, ON, M5S 3H4, Canada  V. Novosad  Materials Sciences Division,
Argonne National Laboratory, 9700 South Cass Avenue, Lemont, IL, 60439, USA
Y. Omori  Department of Astronomy and Astrophysics, University of Chicago,
5640 South Ellis Avenue, Chicago, IL, 60637, USA  Kavli Institute for
Cosmological Physics, University of Chicago, 5640 South Ellis Avenue, Chicago,
IL, 60637, USA  S. Padin  Kavli Institute for Cosmological Physics, University
of Chicago, 5640 South Ellis Avenue, Chicago, IL, 60637, USA  California
Institute of Technology, 1200 East California Boulevard., Pasadena, CA, 91125,
USA  Z. Pan  High-Energy Physics Division, Argonne National Laboratory, 9700
South Cass Avenue., Lemont, IL, 60439, USA  Kavli Institute for Cosmological
Physics, University of Chicago, 5640 South Ellis Avenue, Chicago, IL, 60637,
USA  Department of Physics, University of Chicago, 5640 South Ellis Avenue,
Chicago, IL, 60637, USA  P. Paschos  Enrico Fermi Institute, University of
Chicago, 5640 South Ellis Avenue, Chicago, IL, 60637, USA  S. Patil  School of
Physics, University of Melbourne, Parkville, VIC 3010, Australia  K. A. Phadke
Department of Astronomy, University of Illinois Urbana-Champaign, 1002 West
Green Street, Urbana, IL, 61801, USA  Center for AstroPhysical Surveys,
National Center for Supercomputing Applications, Urbana, IL, 61801, USA  K.
Prabhu  Department of Physics & Astronomy, University of California, One
Shields Avenue, Davis, CA 95616, USA  C. Pryke  School of Physics and
Astronomy, University of Minnesota, 116 Church Street SE Minneapolis, MN,
55455, USA  W. Quan  Department of Physics, University of Chicago, 5640 South
Ellis Avenue, Chicago, IL, 60637, USA  Kavli Institute for Cosmological
Physics, University of Chicago, 5640 South Ellis Avenue, Chicago, IL, 60637,
USA  M. Rahimi  School of Physics, University of Melbourne, Parkville, VIC
3010, Australia  A. Rahlin  Fermi National Accelerator Laboratory, MS209, P.O.
Box 500, Batavia, IL, 60510, USA  Kavli Institute for Cosmological Physics,
University of Chicago, 5640 South Ellis Avenue, Chicago, IL, 60637, USA  C. L.
Reichardt  School of Physics, University of Melbourne, Parkville, VIC 3010,
Australia  M. Rouble  Department of Physics and McGill Space Institute, McGill
University, 3600 Rue University, Montreal, Quebec H3A 2T8, Canada  J. E. Ruhl
Department of Physics, Case Western Reserve University, Cleveland, OH, 44106,
USA  B. R. Saliwanchik  Instrumentation Division, Brookhaven National
Laboratory, Upton, NY, 11973, USA  K. K. Schaffer  Kavli Institute for
Cosmological Physics, University of Chicago, 5640 South Ellis Avenue, Chicago,
IL, 60637, USA  Enrico Fermi Institute, University of Chicago, 5640 South
Ellis Avenue, Chicago, IL, 60637, USA  Liberal Arts Department, School of the
Art Institute of Chicago, 112 South Michigan Avenue, Chicago, IL,60603, USA
E. Schiappucci  School of Physics, University of Melbourne, Parkville, VIC
3010, Australia  C. Sievers  University of Chicago, 5640 South Ellis Avenue,
Chicago, IL, 60637, USA  G. Smecher  Three-Speed Logic, Inc., Victoria, B.C.,
V8S 3Z5, Canada  J. A. Sobrin  Fermi National Accelerator Laboratory, MS209,
P.O. Box 500, Batavia, IL, 60510, USA  Kavli Institute for Cosmological
Physics, University of Chicago, 5640 South Ellis Avenue, Chicago, IL, 60637,
USA  A. A. Stark  Harvard-Smithsonian Center for Astrophysics, 60 Garden
Street, Cambridge, MA, 02138, USA  J. Stephen  Enrico Fermi Institute,
University of Chicago, 5640 South Ellis Avenue, Chicago, IL, 60637, USA  A.
Suzuki  Physics Division, Lawrence Berkeley National Laboratory, Berkeley, CA,
94720, USA  C. Tandoi  Department of Astronomy, University of Illinois Urbana-
Champaign, 1002 West Green Street, Urbana, IL, 61801, USA  K. L. Thompson
Kavli Institute for Particle Astrophysics and Cosmology, Stanford University,
452 Lomita Mall, Stanford, CA, 94305, USA  Department of Physics, Stanford
University, 382 Via Pueblo Mall, Stanford, CA, 94305, USA  SLAC National
Accelerator Laboratory, 2575 Sand Hill Road, Menlo Park, CA, 94025, USA  B.
Thorne  Department of Physics & Astronomy, University of California, One
Shields Avenue, Davis, CA 95616, USA  C. Trendafilova  Center for
AstroPhysical Surveys, National Center for Supercomputing Applications,
Urbana, IL, 61801, USA  C. Tucker  School of Physics and Astronomy, Cardiff
University, Cardiff CF24 3YB, United Kingdom  C. Umilta  Department of
Physics, University of Illinois Urbana-Champaign, 1110 West Green Street,
Urbana, IL, 61801, USA  T. Veach  Space Science and Engineering Division,
Southwest Research Institute, San Antonio, TX 78238  J. D. Vieira  Department
of Astronomy, University of Illinois Urbana-Champaign, 1002 West Green Street,
Urbana, IL, 61801, USA  Department of Physics, University of Illinois Urbana-
Champaign, 1110 West Green Street, Urbana, IL, 61801, USA  Center for
AstroPhysical Surveys, National Center for Supercomputing Applications,
Urbana, IL, 61801, USA  M. P. Viero  California Institute of Technology,
Pasadena, CA 91125, USA  Y. Wan  Department of Astronomy, University of
Illinois Urbana-Champaign, 1002 West Green Street, Urbana, IL, 61801, USA
Center for AstroPhysical Surveys, National Center for Supercomputing
Applications, Urbana, IL, 61801, USA  G. Wang  High-Energy Physics Division,
Argonne National Laboratory, 9700 South Cass Avenue., Lemont, IL, 60439, USA
N. Whitehorn  Department of Physics and Astronomy, Michigan State University,
East Lansing, MI 48824, USA  W. L. K. Wu  Kavli Institute for Particle
Astrophysics and Cosmology, Stanford University, 452 Lomita Mall, Stanford,
CA, 94305, USA  SLAC National Accelerator Laboratory, 2575 Sand Hill Road,
Menlo Park, CA, 94025, USA  V. Yefremenko  High-Energy Physics Division,
Argonne National Laboratory, 9700 South Cass Avenue., Lemont, IL, 60439, USA
M. R. Young  Fermi National Accelerator Laboratory, MS209, P.O. Box 500,
Batavia, IL, 60510, USA  Kavli Institute for Cosmological Physics, University
of Chicago, 5640 South Ellis Avenue, Chicago, IL, 60637, USA  J. A. Zebrowski
Kavli Institute for Cosmological Physics, University of Chicago, 5640 South
Ellis Avenue, Chicago, IL, 60637, USA  Department of Astronomy and
Astrophysics, University of Chicago, 5640 South Ellis Avenue, Chicago, IL,
60637, USA  Fermi National Accelerator Laboratory, MS209, P.O. Box 500,
Batavia, IL, 60510, USA  M. Zemcov  School of Physics and Astronomy, Rochester
Institute of Technology, Rochester, NY 14623, USA  Jet Propulsion Laboratory,
California Institute of Technology, Pasadena, CA 91109, USA

###### Abstract

We report results from an analysis aimed at detecting the trispectrum of the
kinematic Sunyaev-Zel’dovich (kSZ) effect by combining data from the South
Pole Telescope (SPT) and Herschel-SPIRE experiments over a 100 ${\rm deg}^{2}$
field. The SPT observations combine data from the previous and current
surveys, namely SPTpol and SPT-3G, to achieve depths of 4.5, 3, and 16
$\mu{\rm K{\text{-}}arcmin}$ in bands centered at 95, 150, and 220 GHz. For
SPIRE, we include data from the 600 and 857 GHz bands. We reconstruct the
velocity-induced large-scale correlation of the small-scale kSZ signal with a
quadratic estimator that uses two cosmic microwave background (CMB)
temperature maps, constructed by optimally combining data from all the
frequency bands. We reject the null hypothesis of a zero trispectrum at
$10.3\sigma$ level. However, the measured trispectrum contains contributions
from both the kSZ and other undesired components, such as CMB lensing and
astrophysical foregrounds, with kSZ being sub-dominant. We use the Agora
simulations to estimate the expected signal from CMB lensing and astrophysical
foregrounds. After accounting for the contributions from CMB lensing and
foreground signals, we do not detect an excess kSZ-only trispectrum and use
this non-detection to set constraints on reionization. By applying a prior
based on observations of the Gunn-Peterson trough, we obtain an upper limit on
the duration of reionization of $\Delta z_{\rm re,50}<4.5$ (95% C.L). We find
these constraints are fairly robust to foregrounds assumptions. This
trispectrum measurement is independent of, but consistent with, Planck’s
optical depth measurement. This result is the first constraint on the epoch of
reionization using the non-Gaussian nature of the kSZ signal.

_Introduction_. — The kinematic Sunyaev-Zel’dovich (kSZ) effect originates
when electrons with bulk motion Compton-scatter cosmic microwave background
(CMB) photons (Sunyaev and Zeldovich, 1980). Detecting the kSZ signal can
provide crucial insights on both structure formation (Mueller  _et al._ ,
2015; Bianchini and Silvestri, 2016) and the physics of reionization (Knox
_et al._ , 1998; McQuinn  _et al._ , 2005; Zahn  _et al._ , 2012; Reichardt,
2016; Choudhury, 2022; Jain  _et al._ , 2023). This is because the source of
the kSZ signal can be decomposed into two main categories: a low-redshift
($z\lesssim 3$) component referred to as homogeneous- or post-reionization kSZ
and a high-redshift ($z\gtrsim 6$) component referred to as inhomogeneous or
reionization-kSZ. The post-reionization kSZ signal is due to the bulk flow of
halos with free electrons in the local Universe. The reionization kSZ, on the
other hand, is due to motion of the ionized bubbles containing free electrons
during the epoch of reionization (EoR). Several observations (e.g, Becker  _et
al._ , 2001; Fan  _et al._ , 2006; Robertson  _et al._ , 2010) suggest that
the energetic ultraviolet light from the first stars and galaxies at $z\gtrsim
6$ was responsible for ionizing the neutral hydrogen in the early
Universe—although more data is required to precisely understand the process,
timing, and duration of the EoR (Robertson, 2022). From the CMB data, besides
kSZ, the EoR can also be probed using the large-scale bump in the CMB $EE/TE$
power spectra (Zaldarriaga, 1997). While measuring the large-scale $EE/TE$
bump is hard using ground-based experiments, there is also interest in other
EoR probes like kSZ, because the low-$\ell$ $EE/TE$ bump is sensitive to the
cumulative optical depth and cannot distinguish between different reionization
histories (Planck Collaboration  _et al._ , 2016a).

The kSZ signals are sub-dominant compared to other signals in maps of total
intensity at typical CMB observing frequencies and detecting them has proved
to be challenging. In the past, the kSZ signals have been detected through
cross-correlation of CMB maps with galaxy surveys (Hand  _et al._ , 2012; Hill
_et al._ , 2016; Calafut  _et al._ , 2021; Schiappucci  _et al._ , 2023) but
these measurements only probe the post-reionization kSZ signal. Although
cross-correlations with high redshift galaxy catalogs are in principle
possible for the reionization kSZ, the expected signal-to-noise ($S/N$) is
small ($2-3\sigma$) even for future Stage-4 experiments because of the
difficulty in obtaining galaxy catalogues at $z\gtrsim 6$ (La Plante  _et al._
, 2022). Hence, forecasts for kSZ constraints on the EoR have typically relied
on kSZ power spectrum (2-pt function) measurements or reconstructions of the
optical depth (Dvorkin and Smith, 2009). The kSZ 2-pt function, however,
receives contributions from both of the kSZ components and disentangling the
two is difficult due to their similar shapes and amplitudes (Shaw  _et al._ ,
2012; Battaglia  _et al._ , 2013). The presence of astrophysical foregrounds
in the CMB maps, especially thermal SZ (tSZ) and cosmic infrared background
(CIB) signals, complicates the interpretation of the kSZ power spectrum
further as these foreground signals are much brighter than kSZ (Raghunathan
and Omori, 2023, hereafter RO23). As a result, EoR constraints from recent
measurements of the total kSZ power spectrum (Reichardt  _et al._ , 2021;
Gorce  _et al._ , 2022) have been limited by knowledge of the foreground and
post-reionization kSZ signal.

Smith and Ferraro (2017, hereafter SF17) proposed a novel method of using the
kSZ trispectrum (4-pt function) to probe the physics of reionization (also see
Ferraro and Smith, 2018; Alvarez  _et al._ , 2020). Since the reionization kSZ
signal depends on both the free electron density and the velocity, the small-
scale fluctuations in electron density due to inhomogeneous reionization will
get modulated on larger scales by the velocity field, leading to position-
dependent non-Gaussianities in the CMB maps. Thus, the kSZ non-Gaussianity
arises due to the large scale correlations of small scale clustering of halos,
similar to CMB lensing (SF17). However, in our case the large scale
correlations are due to the bulk velocity flow in the Universe rather than due
to gravitational lensing (Okamoto and Hu, 2003). Although the post-
reionization kSZ also starts out as a non-Gaussian signal at different epochs,
the comoving line-of-sight distance over which it gets integrated in the local
Universe is much larger and the signal ends up being Gaussian based on the
central limit theorem. This characteristic allows the two kSZ components to be
easily distinguished using the trispectrum (see Fig. 2 of SF17, ). Despite
this advantage, the constraints on reionization from kSZ 4-pt alone are not
expected to be as competitive as the optical depth measurements from Planck
even for future CMB surveys, due to the degeneracy between the parameters that
govern reionization. However, as demonstrated by Alvarez  _et al._ (2020), the
joint constraints from kSZ 4-pt, Planck primary CMB, and kSZ 2-pt can
effectively break that degeneracy, resulting in a significant improvement on
reionization compared to what can be achieved individually by any of these
probes.

In this work, we present results from an analysis aimed at detecting the kSZ
trispectrum using CMB temperature maps obtained by combining South Pole
Telescope (SPT) and Herschel-SPIRE datasets. The observed trispectrum receives
a contribution from reionization kSZ but is dominated by CMB lensing and
astrophysical foreground signals. We build a template for the latter using the
Agora simulations (Omori, 2022). Given the difficulties in correctly modeling
the foreground signals and the lower amplitude of the kSZ signals compared to
the other undesired signals, we adopt different strategies to handle the
foregrounds. In the baseline case, we marginalize over the amplitude of the
CMB lensing and foreground signals. We also take the approach of fixing the
CMB lensing and foreground signals. In neither case do we observe an excess
kSZ trispectrum. We use this non-detection along with a prior based on the
measurements of the Gunn-Peterson (GP) trough to set upper limits (95% C.L.)
on the duration of reionization $\Delta z_{\rm re,50}$ corresponding to the
difference in redshifts at which the Universe has been 25% and 75% reionized.
We show that our results are consistent with Planck’s optical depth
measurement. We do not combine our results with the kSZ 2-pt measurements from
the literature (Reichardt  _et al._ , 2021; Gorce  _et al._ , 2022) owing to
assumptions made about foregrounds and the post-reionization kSZ signals made
in those works.

_CMB maps_. — This work uses data from two different experiments: SPT (Padin
_et al._ , 2008; Carlstrom  _et al._ , 2011) and Herschel-SPIRE (Pilbratt  _et
al._ , 2010; Griffin  _et al._ , 2010). For SPT, we use data from two surveys:
SPTpol (Austermann  _et al._ , 2012) and SPT-3G (Benson  _et al._ , 2014;
Bender  _et al._ , 2018; Sobrin  _et al._ , 2022). The SPTpol observations
were carried out between 2012 and 2016, and in this work we only use the 150
GHz observations (Henning  _et al._ , 2018) since the noise level of 95 GHz
SPTpol is roughly $\times 3$ higher than the equivalent SPT-3G data. The
SPT-3G observations used in this work were carried out between 2019 and 2020,
and we include data from all the three bands: 95, 150, and 220 GHz. After
combining SPTpol and SPT-3G, the map depths for the three bands are: 4.5, 3,
and 16 $\mu$K-arcmin respectively. For Herschel-SPIRE, we use the data from
600 GHz (500 $\mu$m) and 857 GHz (350 $\mu$m) bands. Since Herschel-SPIRE is
primarily used for CIB mitigation and since the CIB-decorrelation between SPT
bands and Herschel-SPIRE’s 1200 GHz (250 $\mu m$) band is high (Viero  _et
al._ , 2019), we do not use the 1200 GHz band in this work. We limit the SPT
footprint in this work to the region that has overlap with Herschel-SPIRE,
which is a roughly 100 ${\rm deg}^{2}$ region centered at (RA, Decl.) =
(23h30m, -55∘). We provide details about the data processing in the Appendix.
In short, the raw SPT data is filtered and binned into maps with a pixel
resolution of $0.^{\prime}5$. The effect of filtering is accounted for by
using the transfer function (TF) calculated from simulations. The individual
frequency maps are then calibrated by cross-correlating with Planck. The SPIRE
and calibrated SPT maps are combined to produce a minimum-variance (MV) map
which has an unbiased response to CMB temperature. This is done using a scale-
dependent linear combination technique (Cardoso  _et al._ , 2008), and we
refer to the product as the MV-LC.

_Quadratic estimator for $\hat{K}({\bf\hat{n}})$_. — We follow the work of
SF17 and develop a quadratic estimator (QE) to reconstruct
$\hat{K}({\bf\hat{n}})$ which captures the degree-scale correlations of small-
scale clustering of the kSZ signal. The QE is similar to the one used for CMB
lensing reconstruction but here the estimator looks for the velocity-induced,
rather than lensing-induced, correlation between the otherwise independent
$\ell$ modes.

The QE works by computing the product of two Wiener-filtered CMB temperature
maps $X(\hat{\bf n})$ and $Y(\hat{\bf n})$ to construct
$\hat{K}({\bf\hat{n}})=X(\hat{\bf n})Y(\hat{\bf n})$ or equally in harmonic
space as $\hat{K}_{L}=\int d^{2}{\ell}X_{\ell}Y^{\ast}_{\ell-L}$. The Wiener
filter $W_{\ell}$ is employed to down-weight the modes contaminated by sources
such as CMB, astrophysical foregrounds, and instrumental noise. As mentioned
previously, we use the MV-LC map for both $X(\hat{\bf n})$ and $Y(\hat{\bf
n})$. Given that the CMB is much brighter than the kSZ on large scales, we
explicitly zero modes in our maps at $\ell<\ell_{\rm min}$, following SF17.
Similarly, we also remove the contribution from small scales $\ell>\ell_{\rm
max}$ as they are dominated by foregrounds and instrumental noise. The values
used for $\ell_{\rm min}$ and $\ell_{\rm max}$ are different from SF17,
though, and are based on a thorough investigation of the impact of foreground
signals as discussed below. The combination of these bandpass filters and the
Wiener filter for the MV-LC map yields

$\displaystyle W_{\ell}$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{l l}\dfrac{C_{\ell}^{\rm
kSZ}}{(C_{\ell}^{\rm kSZ}+C_{\ell}^{\rm CMB}+N_{\ell}^{\rm MV-
LC})}&,\ell\in[\ell_{\rm min},\ell_{\rm max}]\\\ 0&,~{}{\rm
otherwise}\end{array}\right.$ (3)

where $C_{\ell}^{\rm kSZ}$ is the expected total kSZ power spectrum which we
assume to be flat in $D_{\ell}^{\rm kSZ}\equiv\ell(\ell+1)C_{\ell}^{\rm
kSZ}/2\pi$ with an amplitude of $3\mu K^{2}$ (Battaglia  _et al._ , 2013;
Reichardt  _et al._ , 2021), $C_{\ell}^{\rm CMB}$ corresponds to the lensed
CMB power spectrum calculated using CAMB (Lewis  _et al._ , 2000), and
$N_{\ell}^{\rm MV-LC}$ is the residual noise and foreground power spectrum in
the MV-LC map. We set $\ell_{\rm min}=3300$ and $\ell_{\rm max}=4300$ as the
fiducial values and discuss more about this choice in the Appendix.

The desired kSZ trispectrum
$\hat{C}_{L}^{KK}=\delta(L-L^{\prime})\hat{K}_{L}\hat{K}_{L^{\prime}}^{\ast}$
is the power spectrum of the reconstructed $\hat{K}({\bf\hat{n}})$ map (SF17)
after removing an estimate of the mean-field that arises due to masking. We
debias $\hat{C}_{L}^{KK}$ to account for the above filtering $W_{\ell}$, the
effective beam $B_{\ell,{\rm eff}}$ of the MV-LC map, TF $F_{\ell}$ and the
sky fraction lost ($f_{\rm sky,lost}$) due to masking by dividing the measured
$\hat{C}_{L}^{KK}$ by $\bar{F}^{2}$ where $\bar{F}=f_{\rm sky,final}^{-4}\int
W^{2}_{\ell}B^{2}_{\ell,{\rm eff}}\ F_{\ell}\ d{\rm ln}\ell$ with $f_{\rm
sky,final}=\frac{100}{4\pi}\left(\frac{\pi}{180}\right)^{2}(1-f_{\rm
sky,lost})$. As mentioned before, besides kSZ, the measured trispectrum
$\hat{C}_{L}^{KK}$ contains contributions from CMB lensing ($C_{L,\kappa_{\rm
CMB}}^{KK}$) and foreground ($C_{L,{\rm FG}}^{KK}$) signals. There is also
contribution from the Gaussian disconnected piece arising due to chance
correlations of the two CMB maps, which we represent as $\hat{N}^{(0,KK)}_{L}$
to continue with the CMB lensing analogy. Together, the total measured
$\hat{C}_{L}^{KK}$ can be defined as

$\displaystyle\hat{C}_{L}^{KK}$ $\displaystyle=$ $\displaystyle
C_{L}^{KK}+\hat{C}_{L,{\rm sys}}^{KK}+\hat{N}^{(0,KK)}_{L}.$ (4)

The second term on the right is the systematic $\hat{C}_{L,{\rm
sys}}^{KK}=C_{L,\kappa_{\rm CMB}}^{KK}+C_{L,{\rm FG}}^{KK}+C_{L,\kappa_{\rm
CMB}-FG}^{KK}$ where $C_{L,\kappa_{\rm CMB}}^{KK}$ is the lensing power
spectrum, $C_{L,{\rm FG}}^{KK}$ is the foreground trispectrum which consists
of both unlensed and lensed foreground signals, and the $C_{L,\kappa_{\rm
CMB}-FG}^{KK}$ is the CMB lensing-foreground bispectrum. The foregrounds
correspond to auto- and cross-trispectra from CIB, tSZ, and radio galaxies. We
discuss these contributions in the Appendix. In this work, we ignore higher-
order bias terms like $N^{(1,KK)}_{L},N^{(2,KK)}_{L},N^{(3/2,KK)}_{L}$, some
of which are important for CMB lensing reconstruction with current low-noise
CMB datasets (for definitions of each term, see Madhavacheril  _et al._ ,
2020). We calculate the bias terms $\hat{N}^{(0,KK)}_{L}$ and $\hat{C}_{L,{\rm
sys}}^{KK}$ using simulations described below.

_Simulations_. — As mentioned before, we use the Agora (Omori, 2022)
simulations to model the effects of astrophysical foregrounds (both lensed and
unlensed) and CMB lensing. Agora is a simulation suite containing correlated
multi-component extragalactic sky maps, namely lensed CMB, CIB, kSZ, radio,
and tSZ. The components are also correlated across the different bands. It has
been generated using the dark matter particles and halo catalogs from
MultiDark Planck 2 (MDPL2, Klypin  _et al._ 2016). Note that the kSZ map in
Agora corresponds to the post-reionization kSZ signal due to halos in the
local Universe and not the reionization kSZ.

To model the reionization kSZ, we use the Abundance Matching Box for the Epoch
of Reionization (AMBER111https://github.com/hytrac/amber/tree/main) simulation
package (Trac  _et al._ , 2022; Chen  _et al._ , 2023). AMBER has the
capability to produce the reionization kSZ signal as a function of the
following parameters: midpoint of reionization $z_{\rm re}^{\rm mid}$;
duration222In AMBER, the default parameterization of the duration of
reionization $\Delta z_{\rm re,90}$ is the period corresponding to the
difference in redshifts at which the Universe has been 5% and 95% reionized.
$\Delta z_{\rm re,90}$; asymmetry of reionization $A_{z}$ with respect to
$z_{\rm re}^{\rm mid}$; minimum halo mass $M_{\rm min}^{h}$ used for
reionization; and the opacity of the ionized bubbles parameterized using mean
free path $\lambda_{\rm MFP}$. We find that the changes in the kSZ 4-pt
function are negligible when we tweak the parameters $\theta\in[A_{z},{\rm
log}(M_{\rm min}^{h}),\lambda_{\rm MFP}]$ and fix them to $A_{z}=3$, ${\rm
log}(M_{\rm min}^{h})=8$, and $\lambda_{\rm MFP}=3$ (Chen  _et al._ , 2023).
We set the boxsize to $1{\rm Gpc}/h$ with $512$ particles corresponding to a
resolution of $\sim 2{\rm Mpc}$. We do not find a significant change in the
kSZ 4-pt function signal when we increase the boxsize by a factor of $\times
1.5$ at fixed resolution.

Upon running AMBER for different values of $z_{\rm re}^{\rm mid}$ and $\Delta
z_{\rm re,90}$, and averaging over 10 realizations, we find that the average
kSZ 4-pt function in $L\in[50,300]$ can be described by

$\displaystyle LC_{L}^{KK}=A\left[\dfrac{\Delta z_{\rm re,90}}{\Delta z_{\rm
re,90}^{\ast}}\right]^{\alpha}\left[\dfrac{z_{\rm re}^{\rm mid}}{z_{\rm
re}^{\rm mid^{\ast}}}\right]^{\beta}$ (5)

with $A=2.73\times 10^{-5}$, $\alpha=1.74$, and $\beta=2.52$ for pivot points
$\Delta z_{\rm re,90}^{\ast}=4$ and $z_{\rm re}^{\rm mid^{\ast}}=8$. We
convert $\Delta z_{\rm re,90}$ to $\Delta z_{\rm re,50}$, the standard
definition in literature, using the relation $\Delta z_{\rm re,50}=0.387\Delta
z_{\rm re,90}+0.009$ for $A_{z}=3$ (Chen  _et al._ , 2023).

To test our pipeline we require multiple realizations of the reionization kSZ
and the astrophysical foreground signals. While the former is easier to
generate using AMBER, generating multiple realizations of correlated non-
Gaussian foregrounds is computationally expensive. Instead we add the AMBER-
kSZ signal to the single available Agora full-sky map and extract 100
different sky patches corresponding to the size of our field. Given that our
field is only 100 ${\rm deg}^{2}$, the 100 patches have no overlap, and we
treat each of these patches as a separate sky realization.

We also generate $250$ correlated Gaussian realizations of lensed CMB and
foregrounds using the power spectra $C_{\ell}$ estimated from Agora
simulations. These realizations are used to estimate the mean-field,
$\hat{N}^{(0,KK)}_{L}$, and the covariance. We also include the reionization
kSZ with the above simulations. This is done by creating Gaussian realizations
using the power spectrum of the AMBER kSZ map with $\Delta z_{\rm re,50}=4$
and $z_{\rm re}^{\rm mid}=8$. Modifying the values of $\Delta z_{\rm re,50}$
and $z_{\rm re}^{\rm mid}$ has negligible impact on the above terms.

We filter the simulations using the 2D TF calculated from mock observations.
The filtered simulations are then convolved with the beam corresponding to
each band. Next we add noise to the simulations obtained using the sign-flip
realizations. The simulations are processed in the same way as the data (see
Appendix) and then we construct the MV-LC maps from both the non-Gaussian and
Gaussian versions.

_$\hat{C}_{L}^{KK}$ estimation and Likelihood_. — We pass the MV-LC maps from
both data and the two sets of simulations through the QE pipeline to
reconstruct the $\hat{K}({\bf\hat{n}})$ maps and estimate $\hat{C}_{L}^{KK}$.
For the likelihood calculation, we set a bin width of $\Delta_{L}=50$ and use
$\hat{C}_{L}^{KK}$ in the range $L\in[50,300]$, since beyond that the kSZ 4-pt
function is expected to be sub-dominant compared to other components (SF17).
We assume a Gaussian likelihood of the form

$-2\ln{}\mathcal{L}(\mathbf{d}|\theta)=\sum_{LL^{\prime}}\left(\hat{C}_{L}^{KK}-C_{L}^{KK}\right)\mathbf{\hat{\Sigma}}_{LL^{\prime}}^{KK^{-1}}\left(\hat{C}_{L^{\prime}}^{KK}-C_{L^{\prime}}^{KK}\right)\,,$
(6)

where $\hat{C}_{L}^{KK}$ is the reconstructed 4-pt containing the input kSZ,
CMB lensing and foregrounds after removing $\hat{N}^{(0,KK)}_{L}$;
$C_{L}^{KK}(\theta)\equiv C_{L}^{KK}$ is the model vector as a function of
parameters $\theta\in[z_{\rm re}^{\rm mid},\Delta z_{\rm re,50}]$; and
$\mathbf{\hat{\Sigma}}_{LL^{\prime}}^{KK}$ is the covariance matrix.

_Mean-field, $\hat{N}^{(0,KK)}_{L}$, $\hat{C}_{L,{\rm sys}}^{KK}$ and
$\mathbf{\hat{\Sigma}}_{LL^{\prime}}^{KK}$_. — Similar to CMB lensing
reconstruction, the presence of point source masks, boundary masks, and
spatially varying signals will result in a non-zero “mean-field” bias. Note
that we do not use point source masks as we inpaint the locations of detected
point sources and clusters, and our noise is uniform across the 100 ${\rm
deg}^{2}$ field. As a result, our mean-field should be dominated by our
boundary mask. We derive the mean-field by averaging the
$\hat{K}({\bf\hat{n}})$ maps from 250 Gaussian realizations. This estimate is
subtracted from all the $\hat{K}({\bf\hat{n}})$ before computing
$\hat{C}_{L}^{KK}$. The $\hat{C}_{L}^{KK}$ from $250$ Gaussian realizations is
averaged to estimate the disconnected Gaussian piece $\hat{N}^{(0,KK)}_{L}$.

We use the non-Gaussian realizations to build a template for $\hat{C}_{L,{\rm
sys}}^{KK}$. Specifically, we estimate it as the mean reconstructed 4-pt
function of all the 100 simulations. Given the imperfect knowledge of
foregrounds, we also estimate a systematic error on this template by scaling
the tSZ/CIB signals in the input Agora simulations by 20%. As we show later,
these uncertainties do not have a significant impact on our results.

The covariance matrix used for likelihoods is computed using
$\hat{C}_{L}^{KK}$ measurements from the non-Gaussian realizations. This takes
into account the scatter from both the Gaussian disconnected terms and the
non-Gaussian signals. We also compute the covariance from the Gaussian
realizations $\mathbf{\hat{\Sigma}}_{\rm Gau}^{KK}$ which we use to report the
raw detection significance of the trispectrum measurement.

_$\hat{K}({\bf\hat{n}})$ and $\hat{C}_{L}^{KK}$ measurements_. — The left
panel of Fig. 1 shows the $\hat{K}({\bf\hat{n}})$ reconstructed from a single
non-Gaussian simulation run while the right panel is for data. We note that
the statistical properties look qualitatively similar between the panels. The
simulation includes CMB, foregrounds, and noise but does not contain the
reionization kSZ signal. The similarity between the simulation and data
suggests that our reconstructed $\hat{K}({\bf\hat{n}})$ should be dominated by
$\hat{N}^{(0,KK)}_{L}$ and foregrounds, and the kSZ signal must be sub-
dominant.

Figure 1: Mean-field subtracted $\hat{K}({\bf\hat{n}})$ maps from a single
simulation run in the left panel and data in the right panel. Both the panels
have been smoothed using a Gaussian beam with $\theta_{\rm FWHM}=30^{\prime}$.
The figure illustrates that our data and simulations have qualitatively
similar-looking features. The simulation only contains CMB, foregrounds, and
noise. It does not include the reionization kSZ signal. The similarity of the
simulation with the data map suggests that the reconstructed
$\hat{K}({\bf\hat{n}})$ is dominated by $\hat{N}^{(0,KK)}_{L}$ and
foregrounds.  Figure 2: Reconstructed $\hat{C}_{L}^{KK}$ signals: The violins
represent the distribution of the the measurements from the 100 non-Gaussian
reionization kSZ-free simulations in each $L$ bin and green data points
represent data. The mean of all the simulations ($\hat{C}_{L,{\rm sys}}^{KK}$)
is shown in black. The band around the black curve indicates the systematic
uncertainty in the template obtained by scaling the tSZ signal in the input
simulations by $\pm 20\%$. For each curve, the estimate of
$\hat{N}^{(0,KK)}_{L}$ (blue dash-dotted curve) has been removed. For
reference, we also show the expected kSZ signal from AMBER simulations for
$z_{\rm re}^{\rm mid}=7.69$ and different values of $\Delta z_{\rm re,50}$ in
shades of orange. From the figure, it is evident that the data is consistent
with the non-Gaussian simulations, which do not contain any reionization kSZ
signal, in agreement with the $\hat{K}({\bf\hat{n}})$ maps in Fig. 1.

In Fig. 2 we present the $\hat{C}_{L}^{KK}$ measurements from simulations and
data. These correspond to the results with our fiducial filter choices
$(\ell_{\rm min},\ell_{\rm max})=(3300,4300)$ in Eq.(3). The result from data
(green circles) is consistent with the distribution of the simulations
represented by the violins. With our data, we reject the null hypothesis of a
zero trispectrum at $10.3\sigma$. This raw $S/N$ has been calculated just with
the Gaussian covariance $\mathbf{\hat{\Sigma}}_{\rm Gau}^{KK}$ in Eq.(6). The
mean of all the simulations, which we use as the estimate of $\hat{C}_{L,{\rm
sys}}^{KK}$ is the black solid curve. The semi-transparent band around the
mean is the systematic uncertainty in $\hat{C}_{L,{\rm sys}}^{KK}$ obtained by
scaling the tSZ signal in the input simulations by $\pm 20\%$ roughly
consistent with the uncertainty in the hydrostatic mass bias parameter (Planck
Collaboration  _et al._ , 2016b). The change in the results is only marginal
if we scaled the CIB instead of the tSZ. The blue dash-dotted curve is
$\hat{N}^{(0,KK)}_{L}$ calculated using 250 simulations, and it has been
removed from all the other curves in the figure. The error bars include
contribution both from the scatter in the Gaussian $\hat{N}^{(0,KK)}_{L}$ and
the non-Gaussian signals. For reference, we also show the expected kSZ 4-pt
function signal in shades of orange. These assume a fixed midpoint $z_{\rm
re}^{\rm mid}=7.69$ and different values of duration $\Delta z_{\rm
re,50}\in[1,15]$.

The probability to exceed (PTE or $p$-value), obtained by comparing the
individual simulations (distributions shown using the violins) and data
(green) with the $\hat{C}_{L,{\rm sys}}^{KK}$ template (black) and computing
$p=\chi^{2}_{\rm sims}\geq\chi^{2}_{\rm data}$ is $p=0.23$ indicating the
consistency between the data and the simulations. Note that, similar to Fig.
1, the simulations do not include the reionization kSZ signal. This suggests
that the reionization kSZ signal must be sub-dominant compared to
$\hat{N}^{(0,KK)}_{L}$ and 4-pt function contributions from CMB lensing and
foregrounds $\hat{C}_{L,{\rm sys}}^{KK}$. Indeed, when we remove the
$\hat{C}_{L,{\rm sys}}^{KK}$ estimate from data, the residual
$\hat{C}_{L}^{KK}$ measurement is consistent with a null signal with $p=0.39$.

_Reionization constraints_. — We compare the $\hat{C}_{L}^{KK}$ measurements
obtained above to the expected kSZ signal from AMBER to place constraints on
the EoR parameters. We fit for three parameters $[z_{\rm re}^{\rm mid}$,
$\Delta z_{\rm re,50}$, and $A_{\rm CMB-FG}]$ where $A_{\rm CMB-FG}$ is the
amplitude term for the Agora $\hat{C}_{L,{\rm sys}}^{KK}$ template. In the
Appendix, we use simulations to show that this approach is robust to the
assumptions about the $\hat{C}_{L,{\rm sys}}^{KK}$ template and returns
unbiased results on simulated data.

Figure 3:  Constraints on EoR parameters ($\Delta z_{\rm re,50}$ and $z_{\rm
re}^{\rm mid}$) for different dataset combinations after marginalising over
$A_{\rm CMB-FG}$: We present both 2D and 1D posteriors. The two parameters are
highly degenerate for kSZ 4-pt alone, as shown in the 2D posterior plot. The
black line represents the 68% C.L. joint upper limit to the two parameters in
the 2D plot and the marginalized 1-parameter posteriors in the 1D plots. The
green lines show the analogous curves after applying a prior based on GP
trough measurements which removes regions where reionization ends later than
$z=6$. We adopt these as our baseline constraints, and for this case we set a
95% C.L. upper limit of $\Delta z_{\rm re,50}<4.5$. We also show 68% C.L.
allowed 2D regions and 1D posterior curves derived from Planck’s low-$\ell$
optical depth measurement (blue dash-dotted curve) and from combining Planck
and kSZ 4-pt with (without) the GP-based prior in red (pink). The band around
the combinations that include kSZ 4-pt data represents the uncertainties in
$\hat{C}_{L,{\rm sys}}^{KK}$ template. To avoid cluttering, we show the band
only in the 1d posteriors.

In Fig. 3, we present the constraints on $\Delta z_{\rm re,50}$ and $z_{\rm
re}^{\rm mid}$ for different dataset combinations. As evident from the figure,
the two parameters are degenerate for the kSZ 4-pt-only case shown in the 2D
posterior. On the other hand, evidence based on the measurements of GP trough
(Becker  _et al._ , 2001; Fan  _et al._ , 2006) in the spectra of high
redshift quasars suggest that hydrogen in the Universe must be fully ionized
by $z\sim 5-6$. We use this information to set a binary GP-based prior on the
reionization histories $X_{\rm e}(z)$ from AMBER. To this end, we remove
regions in the $(z_{\rm re}^{\rm mid}-\Delta z_{\rm re,50})$ plane where
reionization ends later than $z<6$ by setting their prior to zero. We define
the end of reionization as the redshift at which $X_{\rm e}(z)>0.95$. The kSZ
4-pt result combined with the GP-based prior is shown in green, and we adopt
this as our baseline result. For this case, we are able to set an upper limit
on $\Delta z_{\rm re,50}$ of $<4.5$ (95% C.L.). Modifying the GP-based prior
to have the reionization end at $z<5$ rather than $z<6$, slightly increases
the upper limit to $\Delta z_{\rm re,50}<5$ (95% C.L.).

We also show constraints based on the Planck optical depth measurement (blue
dash-dotted curve) as well as the kSZ 4-pt measurement + Planck with and
without the GP-based prior (pink and red curves respectively) in Fig. 3. As
expected, Planck (blue) tightly constrains the upper tail of the midpoint of
reionization but is not sensitive to the duration with an upper limit of
$\Delta z_{\rm re,50}$ of $<13.0$ (95% C.L.). Adding kSZ 4-pt to Planck
without the GP-based prior (pink) leads to marginal ($<10\%$) changes in the
results. On the other hand, by including the GP-based prior to kSZ 4-pt +
Planck (red), we set an upper limit on $\Delta z_{\rm re,50}$ of $<2.6$ (95%
C.L.) and obtain $z_{\rm re}^{\rm mid}=7.4\pm 0.6$. This combination is
dominated by Planck and the GP-based prior and kSZ 4-pt only adds marginal
improvement, however, and we do not quote it as our main result.

For the amplitude of the CMB lensing and foreground template, we obtain a
best-fit amplitude of $A_{\rm CMB-FG}=1.25\pm 0.32$. The detection
significance for $A_{\rm CMB-FG}$ is consistent with what we expect from
simulations shown in the inset plots of Fig. 4. As shown in Fig. 5, we do not
observe a strong degeneracy between $A_{\rm CMB-FG}$ and the EoR parameters.

To summarise, we quote the constraint from kSZ 4-pt with the GP-based prior
(green curves in Fig. 3) as our baseline result and we set an upper limit on
the duration of reionization of $\Delta z_{\rm re,50}<4.5$ (95% C.L.). These
are the first constraints on EoR parameters using the non-Gaussianity of the
kSZ signal.

_Foreground uncertainties_. — In Fig. 3, the band around the dataset
combinations that include kSZ 4-pt data represents the uncertainties in
$\hat{C}_{L,{\rm sys}}^{KK}$ propagated to uncertainties in the parameter
constraints. These are obtained by scaling the tSZ signal in the simulations
by 20% (semi-transparent black band around $\hat{C}_{L,{\rm sys}}^{KK}$ in
Fig. 2). As evident from the figure, the uncertainties in $\hat{C}_{L,{\rm
sys}}^{KK}$ have negligible effect and thus our constraints appear to be
relatively robust to assumptions about the CMB+FG template. The above results
are for our fiducial values for $(\ell_{\rm min},\ell_{\rm max})=(3300,4300)$
with $\Delta_{\ell}=1000$. We find these values to be optimal both in terms of
$S/N$ and the foreground biases. In the Appendix, we justify this by
discussing the changes to $S/N$ and the impact of foregrounds when the values
of $\ell_{\rm min},\ell_{\rm max},\Delta_{\ell}$ are modified. We also show
that the systematics due to uncertainties in beam and TF to be negligible.

_Conclusion_. — In this work, we reported results from an analysis aimed at
detecting the non-Gaussian nature of the kSZ signal from reionization. We
combined data from SPTpol, SPT-3G, and Herschel-SPIRE surveys in a 100 ${\rm
deg}^{2}$ field. Using MV-LC map and a QE, we detected the total trispectrum
$\hat{C}_{L}^{KK}$ at $10.3\sigma$. The measurement is dominated by
$\hat{N}^{(0,KK)}_{L}$ bias from the disconnected term, and the contributions
from CMB lensing and foreground signals. After accounting for these undesired
signals, we do not measure an excess kSZ trispectrum and our results are
consistent with a null signal ($p=0.39$). The results from data are consistent
with the expectations from simulations ($p=0.23$). We quantified the biases
due to uncertainties in instrumental beam and TF, and found them to be
negligible. We also thoroughly checked the biases due to mismatch between the
Agora foreground template and data, and found that our results are robust to
uncertainties in both amplitude and shape of the template.

We found that the constraints on $z_{\rm re}^{\rm mid}$ and $\Delta z_{\rm
re,50}$ are highly degenerate from kSZ 4-pt alone. Hence, we applied a loose
prior based on GP trough measurements on the reionization histories from AMBER
to remove $z_{\rm re}^{\rm mid}$ and $\Delta z_{\rm re,50}$ from the parameter
space where the reionization ends later than $z_{\rm end}=6$. With this prior,
we set a upper limit (95% C.L.) of $\Delta z_{\rm re,50}<4.5$. This result is
independent of, but consistent with, the optical depth measurement from
Planck. This work represents the first constraints on EoR parameters using the
trispectrum of the kSZ signal.

Besides the kSZ measurements, the high redshift measurements of quasars
(recently, Dayal  _et al._ , 2024) from JWST (for a review, see Robertson,
2022), the upcoming 21cm measurements from experiments like HERA and SKA (Liu
and Shaw, 2020; The HERA Collaboration  _et al._ , 2022), and the cross-
correlation between multiple probes (for example, Namikawa  _et al._ , 2021;
La Plante  _et al._ , 2023; Georgiev  _et al._ , 2023) are all expected to
significantly enhance our understanding of the physics of EoR in the next
decade. The kSZ results are also expected to improve with the upcoming low-
noise multi-frequency CMB datasets; namely from the full-depth SPT-3G and
SPT-3G+ (Benson  _et al._ , 2014; Bender  _et al._ , 2018; Anderson  _et al._
, 2022), Simons Observatory (Simons Observatory Collaboration, 2019) and
CMB-S4 (CMB-S4 Collaboration  _et al._ , 2019) surveys, which all cover much
wider sky area compared to this work. The combination of kSZ 4-pt function
measurements from these upcoming surveys along with high-significance
measurements of the kSZ power spectrum (RO23) and optical depth measurements
from LiteBIRD (LiteBIRD Collaboration  _et al._ , 2023) can help to break
degeneracies between $z_{\rm re}^{\rm mid}$ and $\Delta z_{\rm re,50}$
(Alvarez  _et al._ , 2020) forming a powerful probe of the epoch of
reionization. This work forms a key first step towards such future high
precision measurements.

The data files and plotting scripts used in this work are available from this
link${}^{\text{\faGithub}}$.

## Acknowledgments

SR acknowledges support by the Illinois Survey Science Fellowship from the
Center for AstroPhysical Surveys at the National Center for Supercomputing
Applications.

This work made use of the following computing resources: Illinois Campus
Cluster, a computing resource that is operated by the Illinois Campus Cluster
Program (ICCP) in conjunction with the National Center for Supercomputing
Applications (NCSA) and which is supported by funds from the University of
Illinois at Urbana-Champaign; the computational and storage services
associated with the Hoffman2 Shared Cluster provided by UCLA Institute for
Digital Research and Education’s Research Technology Group; and the computing
resources provided on Crossover, a high-performance computing cluster operated
by the Laboratory Computing Resource Center at Argonne National Laboratory.

The South Pole Telescope program is supported by the National Science
Foundation (NSF) through award OPP-1852617. Partial support is also provided
by the Kavli Institute of Cosmological Physics at the University of Chicago.

Work at Argonne National Lab is supported by UChicago Argonne LLC, Operator of
Argonne National Laboratory (Argonne). Argonne, a U.S. Department of Energy
Office of Science Laboratory, is operated under contract no. DE-
AC02-06CH11357.

##  Appendix A APPENDIX

##  Appendix B A. Data processing

###  B.1 A.1. Map making

The mapmaking procedure for both the SPTpol and SPT-3G surveys is similar to
previous SPT works (Schaffer  _et al._ , 2011; Henning  _et al._ , 2018;
Dutcher  _et al._ , 2021) and we direct the readers to those works for more
details. Briefly, the time-ordered data (TOD) from each detector are binned
using a flat-sky approximation in the Sanson-Flamsteed projection (Calabretta
and Greisen, 2002; Schaffer  _et al._ , 2011) with a pixel resolution of
$0.^{\prime}5$. We employ the following filtering schemes to the TOD before
binning them into maps. First is the removal of the common-mode, which
corresponds to the mean signal from detectors in a given band. Next, to remove
excess noise in the scan direction $\ell_{x}$, we fit for and remove a
combination of Legendre polynomials (up to $7^{th}$ order) and sines and
cosines (up to an effective high-pass cutoff of $\ell_{x}=300$. Finally, to
prevent the aliasing in the step of binning into map pixels, we low-pass
filter the TOD at a frequency equivalent to $\ell_{x}=13000$. The Herschel-
SPIRE maps used in this work are the same as those used in previous SPT works
(Holder  _et al._ , 2013; Viero  _et al._ , 2019) and we refer the reader to
those works for more details. We calculate the effect of the above TOD
filtering using end-to-end simulations. We make Gaussian realizations of an
underlying power spectrum and mock-observe them with our map-making pipeline.
The ratio of the power spectrum of the output maps over the input maps –
filter transfer function (TF) – captures the effect of filtering. To account
for anisotropic filtering in SPT data, we compute the TF in the 2D as a
function of $\ell_{x}$ and $\ell_{y}$. We denote the azimuthal average of the
2D TF as $F_{\ell}$. We use 250 mock observations of a white noise power
spectrum with $\Delta_{T}=10\mu{\rm K{\text{-}}arcmin}$ to estimate the TF. We
do not note significant differences in TF when we replace the white noise
mocks with CMB and foregrounds realizations. The TF is slightly different at
low $\ell$ for SPTpol and SPT-3G because of the size of the focal-plane, but
the $F_{\ell}$ is close to 95% and roughly flat for all the SPT bands in the
range $\ell\in[3000,5000]$, which is the $\ell$ range most relevant to this
analysis. For Herschel-SPIRE data, a high-pass filter along both
$|\ell_{x}|\lesssim 200$ and $|\ell_{y}|\lesssim 200$ has been applied to
remove excess large-scale noise. This, however, has negligible impact on
$F_{\ell}$ which is unity in our desired $\ell$ range (See Fig. 2 of Viero
_et al._ , 2013). The SPT beam window functions $B_{\ell}$ are measured by
combining dedicated observations of planets and point source signals from CMB
field data (Henning  _et al._ , 2018; Dutcher  _et al._ , 2021). The real-
space beams roughly correspond to Gaussians with $\theta_{\rm
FWHM}=1.^{\prime}7$, $1.^{\prime}2$, $1^{\prime}$ in the 95, 150, and 220 GHz
bands, respectively, though we use the measured $B_{\ell}$ in this analysis.
For Herschel-SPIRE, we approximate the beams as Gaussian with $\theta_{\rm
FWHM}=36.^{\prime\prime}6$ and $25.^{\prime\prime}2$ for the 600 and 857 GHz
bands (Viero  _et al._ , 2019).

###  B.2 A.2. Post-map processing

We obtain the calibration factor $T_{\rm cal}$ for SPT maps by cross-
correlating them with Planck individual frequency maps. Since the SPT
filtering will affect this cross-correlation estimate, we mock-observe the
Planck maps using our map-making pipeline. The value of $T_{\rm cal}$ depends
on the frequency band and the SPT survey. It is obtained as $T_{\rm
cal}=\dfrac{C_{\ell}^{{\rm SPT}\times{\it Planck}}}{C_{\ell}^{\rm SPT-H1\times
SPT-H2}}$ averaged in the multipole range $\ell\in[100,1500]$. Note that we
use SPT observations from the two halves (H1 and H2) to remove the effect of
noise bias.

We optimally combine the calibrated maps from SPT and SPIRE observations using
a harmonic space linear combination technique $S_{\ell}=\sum_{i=1}^{\rm N_{\rm
bands}}w_{\ell}^{i}M_{\ell,m}^{i}$ (Cardoso  _et al._ , 2008). In this study,
we work with a minimum variance MV-LC map for which the frequency-dependent
weights correspond to $w_{\ell}^{\rm MV}=\dfrac{{\bf
C}_{\rm\ell}^{-1}A_{S}}{A_{S}^{\dagger}{\bf C}_{\rm\ell}^{-1}A_{S}}$, where
${\bf C}_{\ell}$ is the covariance matrix containing the covariance between
maps in multiple frequencies at a given $\ell$ and has dimension $\rm N_{\rm
bands}\times\rm N_{\rm bands}$ and $A_{S}=[1...1]$ is a $1\times\rm N_{\rm
bands}$ vector containing the frequency response vector of the CMB. Since the
SPT noise is anisotropic, with strong features along the scan direction, we
perform this combination using a 2D-LC and the weights are a function of both
$\ell_{x}$ and $\ell_{y}$ modes. We deconvolve the beams and the TF from
individual bands before the linear combination step and apply an effective
beam ($B_{\ell,{\rm eff}}$) corresponding to the beam of the 150 GHz channel
of the SPT-3G survey ($B_{\ell}^{\rm 150,SPT-3G}$) to the MV-LC map.

The covariance matrix ${\bf C}_{\ell}={\bf C}_{\ell}^{\rm astro}+{\bf
N}_{\ell}$ receives contribution from both astrophysical signals and the
experimental noise terms. For the signal portion ${\bf C}_{\ell}^{\rm astro}$,
we make use of the Agora simulations. We do not find a significant difference
in the residuals when we switch the Agora-based covariance to data-based (see
Appendix C and Fig. C4 of RO23, ). The noise power spectra ${\bf N}_{\ell}$
for different bands are obtained using sign-flip realizations (Dutcher  _et
al._ , 2021).

Given that the weights are applied in Fourier space, the presence of bright
point sources and clusters in our maps can introduce artefacts. Moreover, like
in CMB lensing, source masking can also lead to significant mean-field bias
(Benoit-Lévy  _et al._ , 2013). To mitigate such effects, we use the
inpainting333 https://github.com/sriniraghunathan/inpainting technique
(Benoit-Lévy  _et al._ , 2013; Raghunathan  _et al._ , 2019) and reconstruct
the primary CMB at the location of sources and clusters. The radius used for
inpainting depends on the flux level of the source and the size of the
clusters. We inpaint locations of point sources detected in our maps with
$S/N\gtrsim 5$ which corresponds to a flux level of $S_{150}\geq 2\ {\rm mJy}$
and clusters detected with $S/N\geq 4.5$. While we lose roughly 18% ($f_{\rm
sky,lost}=0.18$) of the sky due to this stringent masking scheme, it helps in
reducing the level of the foreground signals.

###  B.3 B. Pipeline validation

We use simulations to validate our pipeline and the $\hat{C}_{L,{\rm
sys}}^{KK}$ template. The distribution of $\hat{C}_{L}^{KK}$ from the
simulations are shown as the violins in Fig. 2. The simulations contain noise,
and the the non-Gaussian CMB and foreground signals. But for these tests we
also inject a reionization kSZ signal assuming $z_{\rm re}^{\rm mid}=7.69$ and
$\Delta z_{\rm re,50}=4$ to the above simulations. We report the results using
constraints on the $\Delta z_{\rm re,50}$ parameter. Since $z_{\rm re}^{\rm
mid}$ and $\Delta z_{\rm re,50}$ are highly degenerate, for these tests we
also apply a Gaussian prior on $z_{\rm re}^{\rm mid}$ centred at $z_{\rm
re}^{\rm mid}=7.69\pm 0.715$. This prior roughly corresponds to Planck’s
measurement based on the low-$\ell$ $E$-mode reionization bump. We do note
that Planck’s optical depth measurement error cannot be directly translated
into a flat prior on $z_{\rm re}^{\rm mid}$ since that conversion depends on
the reionization history $X_{\rm e}(z)$. For the simulations, however, the
signal is injected assuming a fixed $z_{\rm re}^{\rm mid}=7.69$ and not based
on $\tau_{\rm re}$, and as a result this flat prior on $z_{\rm re}^{\rm mid}$
should not bias the results. We follow two approaches to account for the
potential mismatch between the $\hat{C}_{L,{\rm sys}}^{KK}$ from data and
Agora.

  * •

Case (1): In this approach, we remove the $\hat{C}_{L,{\rm sys}}^{KK}$
template computed using Agora from the measured $\hat{C}_{L}^{KK}$ and only
fit for the two EoR parameters $\theta\in[z_{\rm re}^{\rm mid},\Delta z_{\rm
re,50}]$. This is optimistic in terms of $S/N$.

  * •

Case (2): In the second approach, we use the Agora $\hat{C}_{L,{\rm
sys}}^{KK}$ template but rather than removing it, we fit for an additional
amplitude term $A_{\rm CMB-FG}$. In this case, we fit for three parameters
$\theta\in[z_{\rm re}^{\rm mid},\Delta z_{\rm re,50},A_{\rm CMB-FG}]$. We
assume a flat prior on $A_{\rm CMB-FG}\in[0,5]$. This is our baseline
approach.

Figure 4: Marginalized posterior distributions from simulations for $\Delta
z_{\rm re,50}$ for Case (1) in the left panels and the baseline Case (2) in
the right panels. The gray curves represent the individual simulations and the
orange curve is the combined posterior from 100 simulations. In the top
panels, the template matches the input simulations (i.e:) $\hat{C}_{L,{\rm
sys}}^{KK}$ is the mean of all the 100 simulations in Fig. 2. The bottom
panels are for the case when the input simulations have the tSZ signal scaled
by $\times 1.2$ but the templates used are the same as in the top panels. In
both the top panels, when the templates match, the orange curve recovers the
input $\Delta z_{\rm re,50}=4$ (green dash-dotted curve) in an unbiased way.
We also note that the width of the orange curve in the right panel is wider
when $A_{\rm CMB-FG}$ is marginalized compared to when it was fixed in the
left panel. In the bottom left panel, however, when the templates are mis-
modeled, the orange curve is significantly biased ($5.5-6\sigma$) when $A_{\rm
CMB-FG}$ is fixed. But in the bottom right panel, we again recover unbiased
results. Also shown in the inset plots are the marginalized posterior for
$A_{\rm CMB-FG}$ parameter. We note in the bottom right, the orange curve
shifts slightly to higher values compared to the top right to account for the
enhanced tSZ signal in the input simulations.

The results are presented in Fig. 4. The individual simulations are in grey
and the combined likelihood from all the $100$ simulations is in orange. In
the figure, left and right panels correspond to Case (1) and Case (2)
described above. The simulations used in the top are the fiducial 100 non-
Gaussian simulations from Fig. 2 while in the bottom panel we scale the tSZ
portion of the simulations by $\times 1.2$. The templates used are always the
same – the black curve from Fig. 2 which is the mean of all the 100
simulations –.

In the top panels, the orange curve recovers the input in both Case (1) and
Case (2) with values $\Delta z_{\rm re,50}=3.9\pm 0.5$ and $\Delta z_{\rm
re,50}=4\pm 1.2$, respectively. The error on $\sigma(\Delta z_{\rm re,50})$ is
higher for Case (2) because of the additional degree of freedom, as expected.
In the bottom panel, however, the orange curve in the left panel is
significantly biased ($5.5-6\sigma$) while the right panel, in which we fit
for the amplitude of $A_{\rm CMB-FG}$, returns an unbiased result. The inset
plots in the right panels show the marginalised posteriors for $A_{\rm CMB-
FG}$. The combined posterior shifts slightly to higher values in the bottom
panel compared to the top panel to account for the enhanced tSZ signal in the
simulations. We also note that the joint likelihood from the 100 simulations
in the top panel results in $A_{\rm CMB-FG}$ that is smaller than 1 and this
could be due to the mismatch in the template shapes between different
simulations. This parameter, however, is not of interest and the shift should
be much smaller for a single simulation run. Moreover, this does not cause any
impact on the recovered $\Delta z_{\rm re,50}$ and hence we do not investigate
this further.

Given these results, we choose Case (2) as our baseline approach. For both
Cases (1) and (2), the detection significance for data or equivalently a
single simulation run is roughly $7.8\sigma/\sqrt{100}$ and
$3.4\sigma/\sqrt{100}$ which are both $<1\sigma$. Subsequently, we do not
expect a detection of $\Delta z_{\rm re,50}$ and only place upper limits (95%
C.L) in this work.

Below, we discuss more about how the biases change when the CIB and tSZ
signals are further modified, and also when $\ell$ ranges used for the filter
in Eq.(3) are modified.

##  Appendix C C. Importance of $\ell$ ranges to maximize $S/N$ and mitigate
biases

The choice of $\ell_{\rm min}$ and $\ell_{\rm max}$ in Eq.(3) depends on two
things: the measurement $S/N$ and the systematic biases arising from CMB
lensing and foregrounds. The value for $\ell_{\rm min}$ is primarily motivated
to reduce the contribution from CMB to the $\hat{N}^{(0,KK)}_{L}$ and
$\hat{C}_{L,{\rm sys}}^{KK}$ estimates. Similarly, the value for $\ell_{\rm
max}$ is chosen to reduce the contribution of foregrounds, which tend to
increase on small scales. Like CMB, foregrounds also impact both the
$\hat{N}^{(0,KK)}_{L}$ and $\hat{C}_{L,{\rm sys}}^{KK}$. While the bias from
the CMB lensing term is relatively easy to model, thanks to high-significance
measurement of CMB lensing from SPTpol in our 100 ${\rm deg}^{2}$ field (Wu
_et al._ , 2019; Millea  _et al._ , 2021), modeling the foreground 4-pt
function signal is hard as it involves understanding each of the foregrounds
(tSZ, CIB, and radio) individually and also the correlation between them.

As described above, we use the MV-LC map without explicitly projecting out any
foreground components using, e.g., constrained-ILC (cILC, Remazeilles  _et
al._ 2011). While removing the tSZ signal can be helpful for tSZ-induced
biases, the process can enhance the CIB residuals significantly (See Fig. 3 of
RO23, ). Similarly, nulling the CIB signal enhances the tSZ residuals. The
process of foreground removal also generally increases the noise in the final
ILC maps. Another approach is to plug in different cILC in the two legs of the
QE as proposed by RO23, although this operation also increases the noise in
the reconstructed $\hat{K}({\bf\hat{n}})$ map. Instead, we use different
$\ell$ cuts on the MV-LC map to minimize the foreground contamination.

Besides the CMB lensing and foreground 4-pt function, which are both
positive444Note that some of the foreground cross-correlations terms can be
negative but the total foreground trispectrum is positive., there is also a
bispectrum bias arising from the correlation of CMB lensing and foregrounds
which is negative on large scales (See Fig.14 of van Engelen  _et al._ ,
2014). If this negative bispectrum bias dominates our measurement, a non-
detection of the kSZ 4-pt function cannot be used to set upper limits on the
EoR parameters without a proper model for the bispectrum bias. In the absence
of the negative bispectrum bias, however, foreground modeling is not necessary
to set place upper limits on EoR using the measured 4-pt function signal.

Our choice of the fiducial $\ell_{\rm min}=3300$ and $\ell_{\rm max}=4300$ is
based on the above arguments. Reducing $\ell_{\rm min}$ increases the
contribution due to CMB lensing. Although that can be modeled as stated above,
it also significantly increases the $\hat{N}^{(0,KK)}_{L}$ and the scatter due
to sample variance of the CMB from these large scales. Similarly, increasing
$\ell_{\rm max}$ increases the contribution from foregrounds. In the case
where the CMB lensing and foregrounds can be perfectly modeled in the
simulations, we do not expect any bias. However, given that the foregrounds
modeled using Agora can be different from data, we estimate the biases due to
imperfect modeling of the simulations as a function of different $\ell_{\rm
min}$ and $\ell_{\rm max}$ values. Below, we quantify these biases and the
impact on $S/N$ when the filter $\ell$ ranges are modified.

###  C.1 C.1. Impact on $S/N$

Here we discuss the changes to $S/N$ when we modify our fiducial values for
$(\ell_{\rm min},\ell_{\rm max})=(3300,4300)$ and $\Delta_{\ell}=1000$. We
limit this test to Case (1) where the CMB+FG is fixed. When we reduce
$\Delta_{\ell}$ from 1000 to 700 (500) keeping the $\ell_{\rm min}$ and
$\ell_{\rm max}$ values fixed, the best-fit $\Delta z_{\rm re,50}$ from the
combined $100$ simulations is $3.8\pm 0.6$ ($3.8\pm 0.7$) compared to $3.9\pm
0.5$ using $\Delta_{\ell}=1000$, implying a reduction in the $S/N$ by roughly
$20\%$ ($40\%$) for the reduced $\Delta_{\ell}$ values. This is not surprising
given the reduced number of modes used for $\hat{K}({\bf\hat{n}})$
reconstruction. Similarly, we also checked the effects of sliding the
($\ell_{\rm min},\ell_{\rm max}$) values lower or higher keeping the
$\Delta_{\ell}=1000$ fixed. Moving them to lower values (large scales) with
$(\ell_{\rm min},\ell_{\rm max})=(3000,4000)$, reduces the $S/N$ by $14\%$ due
to the increased sample variance of CMB. Sliding the filters to higher values
with $(\ell_{\rm min},\ell_{\rm max})=$ (3700, 4700) and (3500, 4500) also
reduces the $S/N$ by 5% and 14% respectively due to the increased scatter from
foregrounds and noise on small scales.

###  C.2 C.2. Bias due to mismatch in $\hat{C}_{L,{\rm sys}}^{KK}$ template

Now, we discuss the bias in the recovered $\Delta z_{\rm re,50}$ from the
mismatch between the true CMB+FG 4-pt function $\hat{C}_{L,{\rm sys}}^{KK}$
template and the one estimated using Agora. We perform this test for both
Case(1) and Case (2) using simulations where we modify the tSZ or the CIB
signals in the input simulations by $\times 1.2$. The $\hat{C}_{L,{\rm
sys}}^{KK}$ template used for the fitting process, however, does not include
the tSZ or CIB scaling.

For the combined 100 simulations in Case (1), we see a $\sim 5.5-6\sigma$
shift due to the mismatch in the tSZ signal used for input and fitting (see
bottom left panel of Fig. 4). While this is large, note that the bias in a
single simulation run is much smaller at $\leq 0.5\sigma$. Sliding the filters
to larger scales $(\ell_{\rm min},\ell_{\rm max})=(3000,4000)$ results in a
slightly smaller (2%) level of bias. Moving the $(\ell_{\rm min},\ell_{\rm
max})$ values of the filters to higher values increases the bias with roughly
$0.73\sigma$ for $(\ell_{\rm min},\ell_{\rm max})=(3700,4700)$. The biases are
$\sim 10\%$ higher when the CIB signal is scaled by $\times 1.2$ rather than
the tSZ. When we reduce the scaling to $\times 1.1$, the bias goes down, as
expected.

For Case (2), when we fit for the amplitude of $A_{\rm CMB-FG}$, however, we
do not observe any bias (see bottom right panel of Fig. 4). When we use the
baseline $\hat{C}_{L,{\rm sys}}^{KK}$ template without tSZ scaling to analyse
the results from simulations where the tSZ is enhanced by $\times 1.2$, we
recover $\Delta z_{\rm re,50}=4\pm 1.3$. At this juncture, we would like to
emphasise that we do not simply mismatch the $\hat{C}_{L,{\rm sys}}^{KK}$
template by a simple amplitude scaling but instead modify the tSZ or CIB
signal to generate $\hat{C}_{L,{\rm sys}}^{KK}$. This process also changes the
shape of the $\hat{C}_{L,{\rm sys}}^{KK}$ template.

Based on these results, we claim our baseline approach Case (2) is not
affected by the mismatch in the shape of the template used.

###  C.3 C.3. Bispectrum bias from CMB lensing and foregrounds

While the bias due to tSZ and CIB is 2% smaller when sliding the filters
towards lower $\ell$ for $(\ell_{\rm min},\ell_{\rm max})=(3000,4000)$
compared to the fiducial $(3300,4300)$, as described above, the $S/N$ goes
down for the former by $14\%$. Another important factor to be considered when
sliding the filters to larger scales is the bispectrum bias due to the
correlation between CMB lensing and foregrounds (van Engelen  _et al._ ,
2014). We estimate this by running two sets of non-Gaussian simulations. The
first case is our baseline 100 simulations from Agora where the lensed CMB and
foreground signals are correlated. For the second case, we randomize the
lensed CMB signals, thus removing the correlation between CMB lensing and
foregrounds. Next we compare the $\hat{C}_{L}^{KK}$ measured from the mean of
the 100 simulations using the two sets. The CMB lensing and foreground 4-pt
function signals should be the same in both the sets but the second set should
not have the bispectrum bias. Subsequently, when we subtract the latter from
the former, we should obtain negative signals if the bispectrum bias dominates
the measurement. For all our filter choices $\ell_{\rm
min}\in[3000,3300,3500,3700]$ with $\Delta_{\ell}=1000$, the difference
between the two sets are consistent with a null signal indicating that the
bispectrum bias is negligible.

Based on the above arguments, our fiducial choice of $(\ell_{\rm
min},\ell_{\rm max})=(3300,4300)$ is an optimal choice in terms of both the
$S/N$ and biases. While it is possible to set $\Delta_{\ell}>1000$ retaining
the fiducial $\ell_{\rm min}=3300$, we do not do so because of the potential
issues due to foregrounds on smaller scales.

##  Appendix D D. Other systematic checks

We now estimate the systematics due to other sources namely the uncertainties
in beam and TF. As mentioned previously, the SPT-3G beams are obtained by
combining measurements of planets and point sources in our field. While
planets have a higher $S/N$, we find issues due to detector nonlinearities
near the peak response and we use stacked observations of point sources in
those regions. Modifying the radius where these two observations are combined
results in differences in beam $B_{\ell}$ which is close to $2\%$ for 90 and
150 GHz bands in our fiducial $(\ell_{\rm min},\ell_{\rm max})=(3300,4300)$
range. Similarly, for the simulations used in this work, we have approximated
the TF for all SPT bands to the 150 GHz SPT-3G TF. In the desired $\ell$
range, we find that there are $\sim 2-5\%$ differences in TF (in power units).
We check the impact of these using simulations and find negligible shifts
($0.05-0.2\sigma$) in the recovered $\Delta z_{\rm re,50}$ values.

##  Appendix E E. Constraints on EoR and CMB+FG template from kSZ 4-pt

Figure 5:  Constraints on EoR ($\Delta z_{\rm re,50}$ and $z_{\rm re}^{\rm
mid}$) and $A_{\rm CMB-FG}$ parameters from kSZ 4-pt only. Similar to Fig. 3,
we present both the 1d and 2D posteriors, and the curves in the 2D plots
correspond to the 68% C.L. allowed regions. We do not observe a strong
degeneracy between $A_{\rm CMB-FG}$ and EoR parameters. The green curves
correspond to the results with the fiducial $\hat{C}_{L,{\rm sys}}^{KK}$ while
blue (yellow) correspond to low (high) $\hat{C}_{L,{\rm sys}}^{KK}$ after
scaling the tSZ signal down (up) by 20% in the Agora simulations.

Here we present the constraints on both EoR and $A_{\rm CMB-FG}$ parameters,
in particular to argue that they are not degenerate. Given this is the case,
we do not report the constraints on $A_{\rm CMB-FG}$ after including the GP-
based prior. This is also evident from the differences in shapes between the
AMBER kSZ 4-pt (in orange shade) and the $\hat{C}_{L,{\rm sys}}^{KK}$ (black)
in Fig. 2. While it is true that we could better constraints on $A_{\rm CMB-
FG}$ by including higher multipoles $L$, given that the degeneracy between
$z_{\rm re}^{\rm mid},\Delta z_{\rm re,50}$ and $A_{\rm CMB-FG}$ is not strong
and the fact that kSZ 4-pt is small compared to foregrounds, we do not extend
the multipole range for parameter constraints. We present constraints for
three cases. Green curves for the fiducial $\hat{C}_{L,{\rm sys}}^{KK}$
template and in this case we obtain $A_{\rm CMB-FG}=1.25\pm 0.32$. The blue
(yellow) curves are when we scale the tSZ signal low (high) by 20% in the
Agora simulations during the template construction. As expected, switching to
a low (high) $\hat{C}_{L,{\rm sys}}^{KK}$ results in higher (lower) values of
$A_{\rm CMB-FG}=1.49\pm 0.32\ (0.88\pm 0.31)$.

## References

  * Sunyaev and Zeldovich (1980) R. A. Sunyaev and Y. B. Zeldovich, MNRAS 190, 413 (1980).
  * Mueller  _et al._ (2015) E.-M. Mueller, F. de Bernardis, R. Bean, and M. D. Niemack, Astrophys. J. 808, 47 (2015), arXiv:1408.6248 [astro-ph.CO] .
  * Bianchini and Silvestri (2016) F. Bianchini and A. Silvestri, Physical Review D 93 (2016), 10.1103/physrevd.93.064026.
  * Knox  _et al._ (1998) L. Knox, R. Scoccimarro, and S. Dodelson, Physical Review Letters 81, 2004 (1998), arXiv:astro-ph/9805012 .
  * McQuinn  _et al._ (2005) M. McQuinn, S. R. Furlanetto, L. Hernquist, O. Zahn, and M. Zaldarriaga, Astrophys. J. 630, 643 (2005), astro-ph/0504189 .
  * Zahn  _et al._ (2012) O. Zahn, C. L. Reichardt, L. Shaw, _et al._ , Astrophys. J. 756, 65 (2012), arXiv:1111.6386 [astro-ph.CO] .
  * Reichardt (2016) C. L. Reichardt, in  _Understanding the Epoch of Cosmic Reionization: Challenges and Progress_, Astrophysics and Space Science Library, Vol. 423, edited by A. Mesinger (2016) p. 227, arXiv:1511.01117 [astro-ph.CO] .
  * Choudhury (2022) T. R. Choudhury, General Relativity and Gravitation 54, 102 (2022), arXiv:2209.08558 [astro-ph.CO] .
  * Jain  _et al._ (2023) D. Jain, T. R. Choudhury, S. Raghunathan, and S. Mukherjee, arXiv e-prints , arXiv:2311.00315 (2023), arXiv:2311.00315 [astro-ph.CO] .
  * Becker  _et al._ (2001) R. H. Becker, X. Fan, R. L. White, _et al._ , Astrophys. J. , submitted (2001).
  * Fan  _et al._ (2006) X. Fan, C. L. Carilli, and B. Keating, ARA&A 44, 415 (2006), astro-ph/0602375 .
  * Robertson  _et al._ (2010) B. E. Robertson, R. S. Ellis, J. S. Dunlop, R. J. McLure, and D. P. Stark, Nature (London) 468, 49 (2010), arXiv:1011.0727 [astro-ph.CO] .
  * Robertson (2022) B. E. Robertson, ARA&A 60, 121 (2022), arXiv:2110.13160 [astro-ph.CO] .
  * Zaldarriaga (1997) M. Zaldarriaga, Phys. Rev. D 55, 1822 (1997), arXiv:astro-ph/9608050 [astro-ph] .
  * Planck Collaboration  _et al._ (2016a) Planck Collaboration, R. Adam, N. Aghanim, M. Ashdown, _et al._ , AAP 596, A108 (2016a), arXiv:1605.03507 [astro-ph.CO] .
  * Hand  _et al._ (2012) N. Hand, G. E. Addison, E. Aubourg, _et al._ , Physical Review Letters 109, 041101 (2012), arXiv:1203.4219 [astro-ph.CO] .
  * Hill  _et al._ (2016) J. C. Hill, S. Ferraro, N. Battaglia, J. Liu, and D. N. Spergel, Physical Review Letters 117 (2016), 10.1103/physrevlett.117.051301.
  * Calafut  _et al._ (2021) V. Calafut, P. Gallardo, E. Vavagiakis, _et al._ , Physical Review D 104 (2021), 10.1103/physrevd.104.043502.
  * Schiappucci  _et al._ (2023) E. Schiappucci, F. Bianchini, M. Aguena, _et al._ , Phys. Rev. D 107, 042004 (2023), arXiv:2207.11937 [astro-ph.CO] .
  * La Plante  _et al._ (2022) P. La Plante, J. Sipple, and A. Lidz, Astrophys. J. 928, 162 (2022), arXiv:2111.13717 [astro-ph.CO] .
  * Dvorkin and Smith (2009) C. Dvorkin and K. M. Smith, Phys. Rev. D 79, 043003 (2009), arXiv:0812.1566 .
  * Shaw  _et al._ (2012) L. D. Shaw, D. H. Rudd, and D. Nagai, Astrophys. J. 756, 15 (2012), arXiv:1109.0553 [astro-ph.CO] .
  * Battaglia  _et al._ (2013) N. Battaglia, A. Natarajan, H. Trac, R. Cen, and A. Loeb, Astrophys. J. 776, 83 (2013), arXiv:1211.2832 [astro-ph.CO] .
  * Raghunathan and Omori (2023) S. Raghunathan and Y. Omori, Astrophys. J. 954, 83 (2023), arXiv:2304.09166 [astro-ph.CO] .
  * Reichardt  _et al._ (2021) C. L. Reichardt, S. Patil, P. A. R. Ade, _et al._ , Astrophys. J. 908, 199 (2021), arXiv:2002.06197 [astro-ph.CO] .
  * Gorce  _et al._ (2022) A. Gorce, M. Douspis, and L. Salvati, AAP 662, A122 (2022), arXiv:2202.08698 [astro-ph.CO] .
  * Smith and Ferraro (2017) K. M. Smith and S. Ferraro, Physical Review Letters 119, 021301 (2017).
  * Ferraro and Smith (2018) S. Ferraro and K. M. Smith, Phys. Rev. D 98, 123519 (2018), arXiv:1803.07036 [astro-ph.CO] .
  * Alvarez  _et al._ (2020) M. A. Alvarez, S. Ferraro, J. C. Hill, R. Hložek, and M. Ikape, arXiv e-prints , arXiv:2006.06594 (2020), arXiv:2006.06594 [astro-ph.CO] .
  * Okamoto and Hu (2003) T. Okamoto and W. Hu, Phys. Rev. D 67, 083002 (2003), arXiv:astro-ph/0301031 .
  * Omori (2022) Y. Omori, arXiv e-prints , arXiv:2212.07420 (2022), arXiv:2212.07420 [astro-ph.CO] .
  * Padin  _et al._ (2008) S. Padin, Z. Staniszewski, R. Keisler, _et al._ , Appl. Opt. 47, 4418 (2008).
  * Carlstrom  _et al._ (2011) J. E. Carlstrom, P. A. R. Ade, K. A. Aird, _et al._ , PASP 123, 568 (2011), arXiv:0907.4445 .
  * Pilbratt  _et al._ (2010) G. L. Pilbratt, J. R. Riedinger, T. Passvogel, _et al._ , AAP 518, L1 (2010), arXiv:1005.5331 [astro-ph.IM] .
  * Griffin  _et al._ (2010) M. J. Griffin, A. Abergel, A. Abreu, _et al._ , AAP 518, L3 (2010), arXiv:1005.5123 [astro-ph.IM] .
  * Austermann  _et al._ (2012) J. E. Austermann, K. A. Aird, J. A. Beall, _et al._ , in  _Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series_, SPIE, Vol. 8452 (2012) arXiv:1210.4970 [astro-ph.IM] .
  * Benson  _et al._ (2014) B. A. Benson, P. A. R. Ade, Z. Ahmed, _et al._ , in  _Millimeter, Submillimeter, and Far-Infrared Detectors and Instrumentation for Astronomy VII_, SPIE, Vol. 9153 (2014) p. 91531P, arXiv:1407.2973 [astro-ph.IM] .
  * Bender  _et al._ (2018) A. N. Bender, P. A. R. Ade, Z. Ahmed, _et al._ , in  _SPIE_, SPIE, Vol. 10708 (2018) p. 1070803, arXiv:1809.00036 [astro-ph.IM] .
  * Sobrin  _et al._ (2022) J. A. Sobrin, A. J. Anderson, A. N. Bender, _et al._ , ApJS 258, 42 (2022), arXiv:2106.11202 [astro-ph.IM] .
  * Henning  _et al._ (2018) J. W. Henning, J. T. Sayre, C. L. Reichardt, _et al._ , Astrophys. J. 852, 97 (2018), arXiv:1707.09353 .
  * Viero  _et al._ (2019) M. P. Viero, C. L. Reichardt, B. A. Benson, _et al._ , Astrophys. J. 881, 96 (2019), arXiv:1810.10643 [astro-ph.CO] .
  * Cardoso  _et al._ (2008) J.-F. Cardoso, M. Le Jeune, J. Delabrouille, M. Betoule, and G. Patanchon, IEEE Journal of Selected Topics in Signal Processing 2, 735 (2008).
  * Lewis  _et al._ (2000) A. Lewis, A. Challinor, and A. Lasenby, Astrophys. J. 538, 473 (2000).
  * Madhavacheril  _et al._ (2020) M. S. Madhavacheril, K. M. Smith, B. D. Sherwin, and S. Naess, arXiv e-prints , arXiv:2011.02475 (2020), arXiv:2011.02475 [astro-ph.CO] .
  * Klypin  _et al._ (2016) A. Klypin, G. Yepes, S. Gottlöber, F. Prada, and S. Heß, MNRAS 457, 4340 (2016), arXiv:1411.4001 [astro-ph.CO] .
  * Trac  _et al._ (2022) H. Trac, N. Chen, I. Holst, M. A. Alvarez, and R. Cen, Astrophys. J. 927, 186 (2022), arXiv:2109.10375 [astro-ph.CO] .
  * Chen  _et al._ (2023) N. Chen, H. Trac, S. Mukherjee, and R. Cen, Astrophys. J. 943, 138 (2023), arXiv:2203.04337 [astro-ph.CO] .
  * Planck Collaboration  _et al._ (2016b) Planck Collaboration, P. A. R. Ade, N. Aghanim, M. Arnaud, _et al._ , AAP 594, A24 (2016b), arXiv:1502.01597 .
  * Dayal  _et al._ (2024) P. Dayal, M. Volonteri, J. E. Greene, _et al._ , arXiv e-prints , arXiv:2401.11242 (2024), arXiv:2401.11242 [astro-ph.GA] .
  * Liu and Shaw (2020) A. Liu and J. R. Shaw, PASP 132, 062001 (2020), arXiv:1907.08211 [astro-ph.IM] .
  * The HERA Collaboration  _et al._ (2022) The HERA Collaboration, Z. Abdurashidova, T. Adams, J. E. Aguirre, _et al._ , arXiv e-prints , arXiv:2210.04912 (2022), arXiv:2210.04912 [astro-ph.CO] .
  * Namikawa  _et al._ (2021) T. Namikawa, A. Roy, B. D. Sherwin, N. Battaglia, and D. N. Spergel, Phys. Rev. D 104, 063514 (2021), arXiv:2102.00975 [astro-ph.CO] .
  * La Plante  _et al._ (2023) P. La Plante, J. Mirocha, A. Gorce, A. Lidz, and A. Parsons, Astrophys. J. 944, 59 (2023), arXiv:2205.09770 [astro-ph.CO] .
  * Georgiev  _et al._ (2023) I. Georgiev, A. Gorce, and G. Mellema, arXiv e-prints , arXiv:2312.04259 (2023), arXiv:2312.04259 [astro-ph.CO] .
  * Anderson  _et al._ (2022) A. J. Anderson, P. Barry, A. N. Bender, _et al._ , in  _Millimeter, Submillimeter, and Far-Infrared Detectors and Instrumentation for Astronomy XI_, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 12190, edited by J. Zmuidzinas and J.-R. Gao (2022) p. 1219003, arXiv:2208.08559 [astro-ph.IM] .
  * Simons Observatory Collaboration (2019) Simons Observatory Collaboration, JCAP 2019, 056 (2019), arXiv:1808.07445 [astro-ph.CO] .
  * CMB-S4 Collaboration  _et al._ (2019) CMB-S4 Collaboration, K. Abazajian, G. Addison, P. Adshead, _et al._ , arXiv e-prints , arXiv:1907.04473 (2019), arXiv:1907.04473 [astro-ph.IM] .
  * LiteBIRD Collaboration  _et al._ (2023) LiteBIRD Collaboration, E. Allys, K. Arnold, J. Aumont, _et al._ , Progress of Theoretical and Experimental Physics 2023, 042F01 (2023), arXiv:2202.02773 [astro-ph.IM] .
  * Schaffer  _et al._ (2011) K. K. Schaffer, T. M. Crawford, K. A. Aird, _et al._ , Astrophys. J. 743, 90 (2011), arXiv:1111.7245 [astro-ph.CO] .
  * Dutcher  _et al._ (2021) D. Dutcher, L. Balkenhol, P. A. R. Ade, _et al._ , Phys. Rev. D 104, 022003 (2021), arXiv:2101.01684 [astro-ph.CO] .
  * Calabretta and Greisen (2002) M. R. Calabretta and E. W. Greisen, AAP 395, 1077 (2002), arXiv:astro-ph/0207413 .
  * Holder  _et al._ (2013) G. P. Holder, M. P. Viero, O. Zahn, _et al._ , ApJL 771, L16 (2013), arXiv:1303.5048 [astro-ph.CO] .
  * Viero  _et al._ (2013) M. P. Viero, L. Wang, M. Zemcov, _et al._ , Astrophys. J. 772, 77 (2013), arXiv:1208.5049 [astro-ph.CO] .
  * Benoit-Lévy  _et al._ (2013) A. Benoit-Lévy, T. Déchelette, K. Benabed, J.-F. Cardoso, D. Hanson, and S. Prunet, AAP 555, A37 (2013), arXiv:1301.4145 [astro-ph.CO] .
  * Raghunathan  _et al._ (2019) S. Raghunathan, G. P. Holder, J. G. Bartlett, S. Patil, C. L. Reichardt, and N. Whitehorn, JCAP 2019, 037 (2019), arXiv:1904.13392 [astro-ph.CO] .
  * Wu  _et al._ (2019) W. L. K. Wu, L. M. Mocanu, P. A. R. Ade, _et al._ , Astrophys. J. 884, 70 (2019), arXiv:1905.05777 [astro-ph.CO] .
  * Millea  _et al._ (2021) M. Millea, C. M. Daley, T. L. Chou, _et al._ , Astrophys. J. 922, 259 (2021), arXiv:2012.01709 [astro-ph.CO] .
  * Remazeilles  _et al._ (2011) M. Remazeilles, J. Delabrouille, and J.-F. Cardoso, MNRAS 410, 2481 (2011), arXiv:1006.5599 [astro-ph.CO] .
  * van Engelen  _et al._ (2014) A. van Engelen, S. Bhattacharya, N. Sehgal, G. P. Holder, O. Zahn, and D. Nagai, Astrophys. J. 786, 13 (2014), arXiv:1310.7023 .