# Superclustering with the Atacama Cosmology Telescope and Dark Energy Survey. I. Evidence for Thermal Energy Anisotropy Using Oriented Stacking

M. Lokken, R. Hložek, A. Van Engelen, M. Madhavacheril, E. Baxter, J.

Derose, C. Doux, S. Pandey, E. S. Rykoff, G. Stein, et al.

► To cite this version:

M. Lokken, R. Hložek, A. Van Engelen, M. Madhavacheril, E. Baxter, et al., Superclustering with the Atacama Cosmology Telescope and Dark Energy Survey. I. Evidence for Thermal Energy Anisotropy Using Oriented Stacking. The Astrophysical Journal, 2022, 933, 10.3847/1538-4357/ac7043. insu-03748248

HAL Id: insu-03748248

https://insu.hal.science/insu-03748248

Submitted on 9 Aug 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Superclustering with the Atacama Cosmology Telescope and Dark Energy Survey. I.
Evidence for Thermal Energy Anisotropy Using Oriented Stacking

M. Lokken¹,²,³ ID, R. Hložek¹,³ ID, A. van Engelen⁴ ID, M. Madhavacheril⁵,⁶ ID, E. Baxter⁷ ID, J. DeRose⁸ ID, C. Doux⁹ ID,
S. Pandey⁹ ID, E. S. Rykoff¹⁰ ID, G. Stein⁸,¹¹ ID, C. To¹² ID, T. M. C. Abbott¹³, S. Adhikari¹⁴,¹⁵ ID, M. Aguena¹⁶ ID, S. Allam¹⁷ ID,
F. Andrade-Oliveira¹⁶,¹⁸, J. Annis¹⁷ ID, N. Battaglia¹⁹ ID, G. M. Bernstein⁹ ID, E. Bertin²⁰,²¹ ID, J. R. Bond² ID, D. Brooks²² ID,
E. Calabrese²³, A. Carnero Rosell¹⁶,²⁴,²⁵ ID, M. Carrasco Kind²⁶,²⁷ ID, J. Carretero²⁸ ID, R. Cawthon²⁹ ID, A. Choi³⁰ ID,
M. Costanzi³¹,³²,³³ ID, M. Crocce³⁴,³⁵ ID, L. N. da Costa¹⁶,³⁶ ID, M. E. da Silva Pereira³⁷, J. De Vicente³⁸ ID, S. Desai³⁹ ID,
J. P. Dietrich⁴⁰ ID, P. Doel²², J. Dunkley⁴¹,⁴² ID, S. Everett⁴³, A. E. Evrard³⁷,⁴⁴ ID, S. Ferraro⁸,¹¹ ID, B. Flaugher¹⁷ ID,
P. Fosalba³⁴,³⁵ ID, J. Frieman¹⁵,¹⁷ ID, P. A. Gallardo⁴⁵ ID, J. García-Bellido⁴⁶ ID, E. Gaztanaga³⁴,³⁵ ID, D. W. Gerdes³⁷,⁴⁴ ID,
T. Giannantonio⁴⁷,⁴⁸ ID, D. Gruen⁴⁰ ID, R. A. Gruendl²⁶,²⁷ ID, J. Gschwend¹⁶,³⁶ ID, G. Gutierrez¹⁷ ID, J. C. Hill⁴⁹,⁵⁰ ID,
M. Hilton⁵¹,⁵² ID, A. D. Hincks¹ ID, S. R. Hinton⁵³ ID, D. L. Hollowood⁴³ ID, K. Honscheid³⁰,⁵⁴ ID, B. Hoyle⁴⁰ ID, Z. Huang⁵⁵ ID,
J. P. Hughes⁵⁶ ID, D. Huterer³⁷ ID, B. Jain⁹ ID, D. J. James⁵⁷ ID, T. Jeltema⁴³ ID, K. Kuehn⁵⁸,⁵⁹ ID, M. Lima¹⁶,⁶⁰ ID, M. A. G. Maia¹⁶,³⁶ ID,
J. L. Marshall⁶¹ ID, J. McMahon¹⁴,¹⁵,⁶²,⁶³, P. Melchior⁴² ID, F. Menanteau²⁶,²⁷ ID, R. Miquel²⁸,⁶⁴ ID, J. J. Mohr⁴⁰,⁶⁵,
K. Moodley⁵¹,⁵² ID, R. Morgan²⁹ ID, F. Nati⁶⁶ ID, L. Page⁴¹ ID, R. L. C. Ogando¹⁶,³⁶ ID, A. Palmese¹⁵,¹⁷ ID, F. Paz-Chinchón²⁶,⁴⁷ ID,
A. A. Plazas Malagón⁴² ID, A. Pieres¹⁶,³⁶ ID, A. K. Romer⁶⁷ ID, E. Rozo⁶⁸ ID, E. Sanchez³⁸ ID, V. Scarpine¹⁷, A. Schillaci⁶⁹,
M. Schubnell³⁷ ID, S. Serrano³⁴,³⁵ ID, I. Sevilla-Noarbe³⁸ ID, E. Sheldon⁷⁰ ID, T. Shin⁹ ID, C. Sifón⁷¹ ID, M. Smith⁷² ID,
M. Soares-Santos³⁷ ID, E. Suchyta⁷³ ID, M. E. C. Swanson²⁶ ID, G. Tarle³⁷ ID, D. Thomas¹⁰,⁷⁴,⁷⁵ ID, D. L. Tucker¹⁷ ID,
T. N. Varga⁶⁵,⁷⁶ ID, J. Weller⁶⁵,⁷⁶ ID, R. H. Wechsler¹⁰,¹²,⁷⁵ ID, R. D. Wilkinson⁶⁷ ID, E. J. Wollack⁷⁷ ID, and Z. Xu⁹,⁷⁸ ID

¹ David A. Dunlap Department of Astronomy & Astrophysics, University of Toronto, 50 Saint George Street, Toronto, ON M5S 3H4, Canada
m.lokken@mail.utoronto.ca

² Canadian Institute for Theoretical Astrophysics, University of Toronto, 60 Saint George Street, Toronto, ON M5S 3H4, Canada

³Dunlap Institute of Astronomy & Astrophysics, 50 Saint George Street, Toronto, ON M5S 3H4, Canada

⁴School of Earth and Space Exploration, Arizona State University, Tempe, AZ 85287, USA

⁵Perimeter Institute for Theoretical Physics, 31 Caroline Street N, Waterloo ON N2L 2Y5 Canada

⁶Department of Physics and Astronomy, University of Southern California, Los Angeles, CA, 90007, USA

⁷Institute for Astronomy, University of Hawaiʻi, 2680 Woodlawn Drive, Honolulu, HI 96822, USA

⁸Lawrence Berkeley National Laboratory, 1 Cyclotron Road, Berkeley, CA 94720, USA

⁹Department of Physics and Astronomy, University of Pennsylvania, Philadelphia, PA 19104, USA

¹⁰Kavli Institute for Particle Astrophysics & Cosmology, P.O. Box 2450, Stanford University, Stanford, CA 94305, USA

¹¹Berkeley Center for Cosmological Physics, 341 Campbell Hall, University of California, Berkeley, CA 94720, USA

¹²Department of Physics, Stanford University, 382 Via Pueblo Mall, Stanford, CA 94305, USA

¹³ Cerro Tololo Inter-American Observatory, NSF's National Optical-Infrared Astronomy Research Laboratory, Casilla 603, La Serena, Chile

¹⁴Department of Astronomy and Astrophysics, University of Chicago, 5640 South Ellis Avenue, Chicago, IL 60637, USA

¹⁵Kavli Institute for Cosmological Physics, University of Chicago, Chicago, IL 60637, USA

¹⁶Laboratório Interinstitucional de e-Astronomia—LIneA, Rua Gal. José Cristino 77, Rio de Janeiro, RJ—20921-400, Brazil

¹⁷Fermi National Accelerator Laboratory, P.O. Box 500, Batavia, IL 60510, USA

¹⁸Instituto de Física Teórica, Universidade Estadual Paulista, São Paulo, Brazil

¹⁹Department of Astronomy, Cornell University, Ithaca, NY 14853, USA

²⁰CNRS, UMR 7095, Institut d'Astrophysique de Paris, F-75014, Paris, France

²¹ Sorbonne Universités, UPMC Université Paris 06, UMR 7095, Institut d'Astrophysique de Paris, F-75014, Paris, France

²²Department of Physics & Astronomy, University College London, Gower Street, London, WC1E 6BT, UK

²³School of Physics and Astronomy, Cardiff University, The Parade, Cardiff, CF24 3AA, UK

²⁴Instituto de Astrofísica de Canarias, E-38205 La Laguna, Tenerife, Spain

²⁵Departamento Astrofísica, Universidad de La Laguna, E-38206 La Laguna, Tenerife, Spain

²⁶Center for Astrophysical Surveys, National Center for Supercomputing Applications, 1205 West Clark Street, Urbana, IL 61801, USA

²⁷Department of Astronomy, University of Illinois at Urbana-Champaign, 1002 West Green Street, Urbana, IL 61801, USA

²⁸Institut de Física d'Altes Energies (IFAE), The Barcelona Institute of Science and Technology, Campus UAB, E-08193 Bellaterra (Barcelona) Spain

²⁹Physics Department, 2320 Chamberlin Hall, University of Wisconsin-Madison, 1150 University Avenue, Madison, WI 53706-1390, USA

³⁰Center for Cosmology and Astro-Particle Physics, The Ohio State University, Columbus, OH 43210, USA

³¹Astronomy Unit, Department of Physics, University of Trieste, via Tiepolo 11, I-34131 Trieste, Italy

³²INAF-Osservatorio Astronomico di Trieste, via G.B. Tiepolo 11, I-34143 Trieste, Italy

³³Institute for Fundamental Physics of the Universe, Via Beirut 2, I-34014 Trieste, Italy

³⁴Institut d'Estudis Espacials de Catalunya (IEEC), E-08034 Barcelona, Spain

³⁵Institute of Space Sciences (ICE, CSIC), Campus UAB, Carrer de Can Magrans, s/n, E-08193 Barcelona, Spain

³⁶Observatório Nacional, Rua Gal. José Cristino 77, Rio de Janeiro, RJ—20921-400, Brazil

³⁷Department of Physics, University of Michigan, Ann Arbor, MI 48109, USA

³⁸Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas (CIEMAT), Madrid, Spain

³⁹Department of Physics, IIT Hyderabad, Kandi, Telangana 502285, India

⁴⁰Faculty of Physics, Ludwig-Maximilians-Universität, Scheinerstraße 1, D-81679 Munich, Germany

⁴¹Joseph Henry Laboratories of Physics, Jadwin Hall, Princeton University, Princeton, NJ, 08544, USA

⁴²Department of Astrophysical Sciences, Princeton University, Peyton Hall, Princeton, NJ 08544, USA

⁴³Santa Cruz Institute for Particle Physics, Santa Cruz, CA 95064, USA

⁴⁴Department of Astronomy, University of Michigan, Ann Arbor, MI 48109, USA

⁴⁵Department of Physics, Cornell University, Ithaca, NY 14850, USA

⁴⁶Instituto de Física Teórica UAM/CSIC, Universidad Autónoma de Madrid, E-28049 Madrid, Spain

⁴⁷Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK
⁴⁸ Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK

⁴⁹ Department of Physics, Columbia University, New York, NY, 10027, USA

⁵⁰ Center for Computational Astrophysics, Flatiron Institute, New York, NY, 10010, USA

⁵¹ 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

⁵³ School of Mathematics and Physics, University of Queensland, Brisbane, QLD 4072, Australia

⁵⁴ Department of Physics, The Ohio State University, Columbus, OH 43210, USA

⁵⁵ School of Physics and Astronomy, Sun Yat-sen University, 2 Daxue Road, Tangjia, Zhuhai, 519082, People's Republic of China

⁵⁶ Department of Physics and Astronomy, Rutgers, The State University of New Jersey, 136 Frelinghuysen Road, Piscataway, NJ 08854-8019, USA

⁵⁷ Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA

⁵⁸ Australian Astronomical Optics, Macquarie University, North Ryde, NSW 2113, Australia

⁵⁹ Lowell Observatory, 1400 Mars Hill Road, Flagstaff, AZ 86001, USA

⁶⁰ Departamento de Física Matemática, Instituto de Física, Universidade de São Paulo, CP 66318, São Paulo, SP, 05314-970, Brazil

⁶¹ George P. and Cynthia Woods Mitchell Institute for Fundamental Physics and Astronomy, and Department of Physics and Astronomy, Texas A&M University, College Station, TX 77843, USA

⁶² Department of Physics, University of Chicago, Chicago, IL 60637, USA

⁶³ Enrico Fermi Institute, University of Chicago, Chicago, IL 60637, USA

⁶⁴ Institució Catalana de Recerca i Estudis Avançats, E-08010 Barcelona, Spain

⁶⁵ Max Planck Institute for Extraterrestrial Physics, Giessenbachstrasse, D-85748 Garching, Germany

⁶⁶ Department of Physics, University of Milano-Bicocca, Piazza della Scienza 3, I-20126 Milano (MI), Italy

⁶⁷ Department of Physics and Astronomy, Pevensey Building, University of Sussex, Brighton, BN1 9QH, UK

⁶⁸ Department of Physics, University of Arizona, Tucson, AZ 85721, USA

⁶⁹ Department of Physics, California Institute of Technology, Pasadena, CA 91125, USA

⁷⁰ Brookhaven National Laboratory, Building 510, Upton, NY 11973, USA

⁷¹ Instituto de Física, Pontificia Universidad Católica de Valparaíso, Casilla 4059, Valparaíso, Chile

⁷² School of Physics and Astronomy, University of Southampton, Southampton, SO17 1BJ, UK

⁷³ Computer Science and Mathematics Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831, USA

⁷⁴ Institute of Cosmology and Gravitation, University of Portsmouth, Portsmouth, PO1 3FX, UK

⁷⁵ SLAC National Accelerator Laboratory, Menlo Park, CA 94025, USA

⁷⁶ Universitäts-Sternwarte, Fakultät für Physik, Ludwig-Maximilians Universität München, Scheinerstraße 1, D-81679 München, Germany

⁷⁷ NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA

⁷⁸ MIT Kavli Institute, Massachusetts Institute of Technology, Cambridge, MA, 02139 USA

Received 2021 July 19; revised 2022 May 5; accepted 2022 May 15; published 2022 July 8

## Abstract

The cosmic web contains filamentary structure on a wide range of scales. On the largest scales, superclustering aligns multiple galaxy clusters along intercluster bridges, visible through their thermal Sunyaev–Zel’dovich signal in the cosmic microwave background. We demonstrate a new, flexible method to analyze the hot gas signal from multiscale extended structures. We use a Compton y-map from the Atacama Cosmology Telescope (ACT) stacked on redMaPPer cluster positions from the optical Dark Energy Survey (DES). Cutout images from the y-map are oriented with large-scale structure information from DES galaxy data such that the superclustering signal is aligned before being overlaid. We find evidence of an extended quadrupole moment of the stacked y signal at the $3.5\sigma$ level, demonstrating that the large-scale thermal energy surrounding galaxy clusters is anisotropically distributed. We compare our ACT × DES results with the Buzzard simulations, finding broad agreement. Using simulations, we highlight the promise of this novel technique for constraining the evolution of anisotropic, non-Gaussian structure using future combinations of microwave and optical surveys.

*Unified Astronomy Thesaurus concepts: Cosmic web (330); Large-scale structure of the universe (902); Intergalactic filaments (811); Sunyaev-Zeldovich effect (1654); Cosmology (343); Superclusters (1657); Galaxy clusters (584)*

## 1. Introduction

The anisotropic clustering of galaxies, galaxy clusters, and intergalactic matter provides a unique insight into the development of large-scale structure (LSS) in our universe. Superclusters and filaments, often referred to as discrete objects, are in reality part of a continuous network of matter. Novel statistical methods focused on the anisotropic scale-dependent aligning tendencies of clusters, galaxies, gas, and dark matter in this network are needed to assess the relative amplitudes that the various species contribute to the overall superclustering. We use this term to refer to elongated nonlinear overdense structures that span a wide range of scales

(typically tens of megaparsecs, but also including near-cluster-core sizes of a few megaparsecs as well as some correlated structures that extend beyond 100 Mpc). The formation and evolution of superclustering is highly dependent on the cosmological model, and thus studying large populations of filaments and superclusters may provide key constraining power to discriminate between different cosmologies (Cen 1994; Frisch et al. 1995; Basilakos et al. 2001; Kolokotronis et al. 2002; Bharadwaj & Pandey 2004; Hopkins et al. 2005; Bagchi et al. 2017; Ho et al. 2018).

In this paper, we present a novel way of using the *oriented stacking* method with a combination of multiwavelength data. As a proof of concept, we use the method to assess the anisotropy of the thermal energy distribution surrounding galaxy clusters. We begin with an overview of the theoretical and observational landscape of superclustering.

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
## 1.1. An Overview of Superclustering

Observations, simulations, and analytic theory have converged upon a common model of LSS formation. The early universe was a near-uniform field of Gaussian random density fluctuations, which evolved through gravitational instabilities into the web-like structure that exists today. This continuous network of matter, dubbed the cosmic web by Bond et al. (1996, hereafter BKP), has complex features across a wide range of scales. Its large-scale pattern is predictable from features in the early universe density field, namely, the locations of rare mass peaks and the surrounding large-scale tidal fields (Bardeen et al. 1986; BKP; Bond & Myers 1996). Meanwhile, the small-scale details are products of complex local gravitational interactions in the late-time universe.

From the very earliest cosmological simulations, filaments emerged as the dominant characteristics of the cosmic web. Two distinct theories of their formation (the pancake picture of Zel'Dovich 1970 and hierarchical clustering) were synthesized into the BKP model, as reviewed in van de Weygaert & Bond (2008a, 2008b). This model demonstrates how close pairs of clusters are bridged by filaments with a strength determined by the proximity and alignment of the cluster neighbors. However, the alignment of structure is not limited to the cluster-bridge scale. Clusters themselves can align and cluster, forming superclusters. This superclustering of clusters and galaxies is ubiquitous at lower redshifts, and is also played out at higher redshifts in protoclusters.

Supercluster regions consist not only of clusters and filaments; their 3D complexity includes membranes joining filaments and large voids. Despite being the largest nonlinear structures in the universe, they are far from dynamical equilibrium. Hence, there is no simple (e.g., spherical) shape readily usable for analysis. Though supercluster regions always have some alignment, a simple straight filament picture does not adequately capture the way in which filaments arc between clusters as the orientations of the clusters pass from perfect alignment to partial alignment.

In recent years, the advent of larger and larger surveys has driven the cosmology field away from localized measures of LSS and toward global statistical descriptors. The isotropic two-point correlation function and its Fourier transform, the power spectrum, are frequently used to describe the clustering of matter (a review is presented in Coil 2013) and have provided strong constraints on cosmological parameters (see Alam et al. 2017, for the cosmological results from DR12 of BOSS). However, these methods are only sensitive to isotropic clustering, as they are functions of a directionless distance (for the correlation function) or wavenumber magnitude $|k|$ (for the power spectrum). The lowest-order statistic beyond power is the bispectrum (a review is provided in Section 4 of Desjacques et al. 2018). This is the Fourier transform of the three-point correlation function (Peebles 1980). The bispectrum encodes information about non-Gaussianity in the late universe, measuring significant cosmological information beyond the power spectrum, which can be used to probe dark energy (Takada & Jain 2004; Sefusatti et al. 2006). However, it is expensive to measure in its full glory with all possible wavenumber triangle configurations, and reductions of the wavenumber possibilities to specific choices restricts the measurable information. (Desjacques et al. 2018). In the era of large surveys, localized measurements of superclustering are still needed, especially those that consider local alignments and can be applied to a wide range of scales.

## 1.2. Motivation

With the goal of pushing beyond the limitations of correlation functions, which are ensemble averages of localized measurements centered at every point in the universe, this paper presents a new way to explore anisotropic structure through localized measurements at selected centers in real space. We emphasize the distinction between statistical and localized anisotropy; in this paper we refer only to local anisotropies, i.e., the local variations as a field is rotated around a selected point. While the superclustering of matter can be examined through multiple probes, this work focuses on signals from hot gas for a proof of concept. Measurements of the gas anisotropy surrounding galaxy clusters have the potential to address cosmological and astrophysical questions across a wide range of scales; Table 1 provides a brief overview. We focus the discussion on low redshifts ($z < 1$), where galaxy survey data (a key component of our methods) is most abundant.

For a homogeneous and statistically isotropic Gaussian random field (GRF) in the linear phase, all information, including local anisotropic aspects, is encoded in the 3D power spectrum. If there is primordial non-Gaussianity, the power spectrum is insufficient to capture the full information of the LSS fields even on linear scales. Thus higher-order statistics like the bi- and tri-spectra are often used to search for primordial non-Gaussianity, in the hopes of distinguishing between different models of inflation (Giannantonio et al. 2012). However, some types of primordial non-Gaussianity may also be detectable through methods (such as those presented in this work) that measure, in real space, any excess in the anisotropy of different tracers around selected points in the cosmic web. These searches would be best applied at large scales which are less affected by the non-Gaussianity induced by late-time gravitational evolution.

On scales for which the universe has evolved beyond the linear regime, gravitational collapse drives runaway local anisotropy and non-Gaussianity of the dark matter (and baryons). At late times, dark energy changes the super-clustering pattern. Hence, localized probes of extended structure ranging from scales of a single intercluster bridge to a many-cluster superstructure encode information on the nature of the dark components that may be obscured in the global $k$-space compression of the information onto isotropic power. Understanding the gas and galaxy content of filaments is also important for galaxy evolution studies. For example, studies have shown that the position of a galaxy relative to filaments correlates with its spin and spin alignment (Codis et al. 2018a; Krolewski et al. 2019; Welker et al. 2020) as well as its mass, morphology, star formation rate, nuclear activity, and feedback mechanisms (Darragh Ford et al. 2019; Kraljic et al. 2020; Santiago-Bautista et al. 2020). Overall, the relationship between the dark matter, galaxies, and gas in filamentary structures contains a wealth of cosmological and astrophysical information. Our succeeding paper will focus on measuring the relationship between the galaxy number density and gas thermal energy using the methods presented herein.

On smaller scales, the orientation of the gas and galaxies within and on the outskirts of clusters is determined by highly nonlinear dynamics. While the anisotropy is certainly influenced by the properties of dark energy and dark matter, the many unknowns in complex dynamical processes such as mergers and splashback make cosmological information difficult to disentangle. However, the dynamics and assembly
<table><caption>**Table 1**
Motivation for Making Localized Measurements of Gas Anisotropy at Various Scales</caption><thead><tr><th>Object Class(es)</th><th>Approx. Long-axis Length <i>L</i></th><th>Science Case(s)</th></tr></thead><tbody><tr><td>Galaxy clusters and local surroundings</td><td><i>L</i> &le; 3<i>h</i><sup>-1</sup> Mpc</td><td>Filaments feeding clusters, cluster assembly processes, baryonic feedback</td></tr><tr><td>Intercluster filaments, Small to mid-sized superclusters</td><td>3<i>h</i><sup>-1</sup> Mpc &le; <i>L</i> &le; 40<i>h</i><sup>-1</sup> Mpc</td><td>Cosmological model (dark energy, dark matter), gas dynamics from filament compression, baryonic feedback</td></tr><tr><td>Largest superclusters</td><td><i>L</i> &ge; 40<i>h</i><sup>-1</sup> Mpc</td><td>Potential tests of primordial non-Gaussianity</td></tr></tbody></table>

history of such objects is interesting in its own right. Studying how the cluster gas orientation relates to its surrounding filament(s) could provide insight into the process by which the cosmic web feeds cold gas and galaxies into clusters (e.g., as recently studied in The ThreeHundred project in Kotecha et al. 2022; Kuchner et al. 2022).

At all scales, baryonic feedback processes are a confounding factor in studying anisotropies in the gas thermal energy content for cosmological purposes. The complex small-scale processes that move baryons out of dark matter halos and change gas temperatures and pressures are not well understood, and must be constrained in order to glean cosmological information from thermal energy measurements. Understanding the effects of feedback on clusters, groups, and filaments will in turn improve our understanding of the feedback sources, namely, black holes and massive stars.

Thus, localized measurements of thermal energy anisotropy are motivated by various science cases across a wide range of scales. For any scale, the questions of interest can be addressed by comparing the characteristics of the observational signal with simulations run under varying sets of cosmological and astrophysical assumptions.

For the purposes of a proof of concept we focus this work on the scales of intercluster filaments and cluster pairs (~4–12 $h^{-1}$ comoving Mpc, further motivated in Section 4.2). For this scale, we focus on the redshift range where the available galaxy and cluster data are most abundant (described in Section 2), $0.2 < z < 0.7$.

## 1.3. Introduction to Methods

We measure the superclustering of thermal energy through the signals imprinted on cosmic microwave background (CMB) data from hot gas in galaxy clusters, groups, and filaments. We also incorporate data from a large-sky galaxy survey into the method to provide necessary LSS information. The galaxy number density, a biased tracer of the total matter, is also of interest as a signal and not only an intermediate step; our succeeding paper will address this. The mass field can be more directly probed with weak lensing maps (as done for filaments in Yang et al. 2020); we leave this for future exploration.

Hot gas in the universe is visible through the thermal Sunyaev–Zel’dovich (tSZ) effect. The tSZ effect is observed along lines of sight that pass through hot gas, arising because a small percentage of incoming cold photons from the primary CMB scatter off of hot free electrons in the intervening gas (Zel’dovich & Sunyaev 1969; Sunyaev & Zel’dovich 1970, 1972). In this inverse Compton scattering process, the photons gain energy and the observed CMB spectrum shifts toward

higher frequencies. This causes a shortage of CMB photons at frequencies lower than ~218 GHz and an excess of photons at higher frequencies as compared to the usual CMB spectrum. Thus, galaxy clusters appear as decrements in CMB maps with frequency < 218 GHz, and increments in higher frequency maps. The strength of the effect is parameterized by the dimensionless Compton-$y$ parameter,

$$y = \int d\ell \, n_e \frac{k_B T_e}{m_e c^2} \sigma_T, \quad (1)$$

a line-of-sight integral of the number density of electrons $n_e$ and the electron temperature $T_e$. $\sigma_T$ is the Thomson cross section and $m_e$ is the electron mass. Maps of $y$ are typically constructed with component separation techniques from linear combinations of CMB maps at multiple frequencies both higher and lower than 218 GHz (Remazeilles et al. 2011).

Recent improvements in ground-based tSZ surveys have yielded catalogs of galaxy clusters readily apparent in Compton-$y$-maps due to their high temperatures and densities. Details of the interior cluster gas are revealed through $y$-maps, as the $y$ profile of a cluster is a 2D (cylindrical) projection of the cluster ionized gas 3D pressure profile (Carlstrom et al. 2002; Mroczkowski et al. 2019). The relationship between the cluster halo mass $M_e$ and angular-integrated Compton-$y$ parameter, $Y$, which is proportional to $\langle P_e \rangle V_c$ (the electron pressure times cluster volume), is fairly consistent with $Y \propto M^{5/3}$, as occurs in a gas with adiabatic index $\gamma = 5/3$ in equilibrium (see Bond & Myers 1996; McCarthy et al. 2003; Giodini et al. 2013, and references therein).

However, for lower-mass halos ($M \le 10^{14} M_\odot$)—such as those that host galaxy groups along intercluster filaments—this relationship no longer holds. Lim et al. (2018), Hill et al. (2018), and others have shown that these smaller halos have a steeper $Y-M$ relation, primarily because gas is blown out of halos by feedback from active galactic nuclei (AGN), changing the relation from that which is expected by gravitational arguments alone. The extent to which feedback mechanisms from AGN, supernovae, magnetic fields, and other sources determine the $\langle P_e \rangle V_c$ of lower-mass clusters and groups and redistribute gas beyond halo boundaries is still not fully understood. Even for massive clusters, AGN feedback has been a needed addition to relate cluster observations to theory (Sijacki et al. 2007; Puchwein et al. 2008; Sijacki et al. 2008; Battaglia et al. 2010; Gitti et al. 2012).

Therefore, it is difficult to probe superclustering from the tSZ effect due to the wide variation in gas pressure in different parts of the cosmic web. The weak tSZ signals from lower-mass
halos are drowned out by noise in individual images, and the signals from gas outside of halos are even weaker. Researchers frequently employ *stacking* to extract information from these low-signal regions. In stacking, multiple similar images are averaged. The signal is overlaid throughout the stack and the random noise averages down with more images. If all of the *N* component images have roughly the same signal, the final stack has a higher signal-to-noise ratio (S/N) than any of individual images by a factor of $\sqrt{N}$. Stacking has been used to study the relationship between the tSZ signal and other properties of clusters using data from both space- and ground-based instruments (see, e.g., Plagge et al. 2010; Hand et al. 2011; Sehgal et al. 2011; Planck Collaboration et al. 2013; Sehgal et al. 2013; Planck Collaboration et al. 2016a). It has also been used to study small halos through the thermal and kinetic Sunyaev–Zel’dovich effects in recent work by Schaan et al. (2021), Amodeo et al. (2021), who found that the gas profile extends well beyond the virial radius of the dark matter halo.

In recent years, various teams employed stacking to study intercluster filaments and attempt to detect the warm-hot intergalactic medium (WHIM), low-density gas thought to have been blown out from dark matter halos by feedback processes. This gas at $10^5-10^7$ K may make up ~40%–50% of baryons (Cen & Ostriker 2006; Shull et al. 2012). Thus, despite its lower tSZ signal compared to collapsed objects, it is an important contributor to the census of cosmological baryons (Cen & Ostriker 1999). By stacking ~1 million pairs of luminous red galaxies (typically found at the center of clusters), de Graaff et al. (2019) and Tanimura et al. (2019) detected a filament tSZ signal. The signal emerged from a combination of galaxies, groups, and the WHIM. The former group also claimed detection of the WHIM gas itself by accounting for the massive galaxies in each filament. In addition to tSZ evidence, recent studies have found tentative signals from the WHIM through absorption lines from filaments intersecting the lines of sight of quasars (Tejos et al. 2016; Pessa et al. 2018; Bouma et al. 2021). Altogether, these works have been important confirmations of the existence of a large amount of baryons in low-mass halos and intergalactic gas, accounting for most or all of the so-called *missing baryons* (as reviewed in Tumlinson et al. 2017).

Despite these successes, thus far, most studies of filamentary structure have been limited in scale. Cluster-pair stacking is limited by the distance between pairs: in the aforementioned studies, the filaments ranged from 6–10 $h^{-1}$ Mpc. Therefore, this method cannot be used to study very small-scale alignments of galaxies nearby a cluster, nor superclusters tens or even hundreds of megaparsecs long. Some recent work has used alternative methods to probe a wider range of scales (e.g., Tanimura et al. 2020), but more work is needed to fully explore the multiscale information content of tSZ anisotropy.

## 1.4. *Oriented Stacking*

Oriented stacking is a more general approach to stacking elongated structures that does not necessitate the identification of cluster pairs nor distinct filaments. This method seeks to quantify anisotropic superclustering over a range of scales. Oriented stacking was used in Battaglia et al. (2012a) to examine the anisotropy of gas around clusters on relatively small (~5 Mpc) scales in simulations; it was also applied to studies of CMB polarization in Planck Collaboration et al. (2016b). Our paper presents its first application to a

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Illustration of unoriented (top) vs. oriented (bottom) stacking. When stacking cluster cutout images from the Compton-y-map with random orientation, the signal from external structure is averaged out such that the final image, as the number of cutouts increases, approaches perfect isotropy. In oriented stacking, we determine the direction of strongest alignment around each cluster using the curvature of galaxy overdensity maps. We rotate each cutout image before stacking such that the excess y signal from anisotropic superclustering is overlaid.</figcaption></figure>

combination of CMB and galaxy data. This work complements a recent theoretical exploration of anisotropic superclustering in full 3D and in 2D projections given in Regaldo-Saint Blancard et al. (2021).

In Figure 1 we illustrate how oriented stacking is applied to a combination of a CMB y-map and galaxy and cluster surveys to tease out the faint tSZ signals from superclustering. In this method, we stack the tSZ signal from both clusters and the surrounding gas by aligning along a measured axis of LSS. We go beyond previous filament studies by applying characteristics of the galaxy field to select on areas of high superclustering, which we control, augmenting the signal relative to the fluctuation noise. We demonstrate how multiscale selections on superclustering features can probe how the gas signal changes with scale. After performing oriented stacks, we quantify our results using a multipolar decomposition of the results. We compare the observational results with simulations to search for any discrepancies indicating inaccuracies in the cosmology, galaxy, and cluster sample selection, and/or gas prescription in the simulations.

Because the number of stacked objects in this work is 2–3 orders of magnitude smaller than the cluster-pair stacking studies, we do not expect to significantly detect the WHIM signal. Rather, the dominant signals come from clusters and galaxy groups, and thus our methods generally probe the anisotropic clustering of thermal energy. In future work, we will determine whether the methods described herein can be used to find explicit evidence for hot gas outside of halos. Future work will also explore the potential of oriented stacking
to address the various science questions in Table 1, including
how to disentangle the small-scale baryonic physics from
cosmological effects and ultimately search for signs of physics
beyond the standard model.

This paper is organized as follows. In Section 2, we describe
each data product. Section 3 describes the properties of the
galaxy density field we use to find regions of high super-
clustering. Section 4 describes the stacking methods used.
Section 5 presents the expected signal from the Websky (Stein
et al. 2020) and Buzzard (DeRose et al. 2019) simulations with
various cluster populations. In Section 6 we compare the
Atacama Cosmology Telescope × Dark Energy Survey
(ACT × DES) results with the Buzzard simulations and
forecast for future data. Finally, in Section 7 we discuss the
prospects of applying these novel methods to future super-
clustering analysis.

Throughout the paper, when a cosmological model is
assumed for conversions from redshift to comoving distance
and from angular size to transverse comoving distance, we use
the Planck15 cosmology from the astropy.cosmol-
ogy<sup>79</sup> package which implements a flat Λ cold dark matter
(CDM) model with parameters from Planck Collaboration et al.
(2016c). The model has Ω<sub>M</sub> = 0.308, Ω<sub>k</sub> = 0, Ω<sub>Λ</sub> = 0.691, a
single species of massive neutrinos with mass 0.06 eV, and
H<sub>0</sub> = 67.7 km s<sup>-1</sup> Mpc<sup>-1</sup>. All quoted distances are in comoving
units.

## 2. Data

In this work we combine the tSZ Compton-y signal maps
derived from high-resolution ACT measurements with DES
cluster and galaxy catalogs, and then compare our observa-
tional results to simulations of LSS.

### 2.1. ACT Compton-y-map

ACT is a 6 m off-axis Gregorian telescope, located in the
Atacama Desert of northern Chile, at an elevation of 5190 m on
the Cerro Toco stratovolcano (Fowler et al. 2007; Swetz et al.
2011). The telescope has been operating since 2008, first
measuring only temperature fluctuations in the microwave
regime. In 2013, the ACTPol receiver was deployed, enabling
ACT to observe both temperature and polarization data at
98 and 150 GHz (Thornton et al. 2016). The receiver was
subsequently upgraded to Advanced ACTPol, adding three
more frequencies (Henderson et al. 2016; Ho et al. 2017; Choi
et al. 2018; Li et al. 2018). ACT has recently produced maps
covering 18,000 deg² of the sky (Aiola et al. 2020; Naess et al.
2020). This paper focuses only on the region for which a tSZ
map has been made, which also overlaps with part of the DES
footprint. This 456 deg² region is called “D56” and was first
presented in Louis et al. (2017).

We use the Compton-y-map in D56 first presented in
Madhavacheril et al. (2020). This map was reconstructed from
a combination of 2015–2016 nighttime data from the ACTPol
receiver at 98 and 150 GHz as well as multifrequency data from
the Planck satellite (Planck Collaboration et al. 2016a). Data
from Planck were isolated to include only modes between
20 < l < 300 for the Planck low-frequency instruments (LFI)
at 30 GHz and 44 GHz, respectively, 20 < l < 2000 for
the 70 GHz LFI and modes between 20 < l < 5800 for the

high-frequency instruments (HFI) at 100, 143, 217, 353, and
545 GHz. Data from ACT include modes $500 < l < 24,000$.
The maps are combined via an internal linear combination
(ILC) algorithm to isolate the y component in each Fourier
pixel using the estimated covariance between the map arrays at
different frequencies. The resulting y-map has an effective
1'6 beam.

The y-map may contain residual contribution from the
primary CMB and astrophysical foregrounds. For the purposes
of this work, the residual of greatest concern is the cosmic
infrared background (CIB): this emission from dusty galaxies is
highly correlated with the tSZ effect. Radio point sources are a
far subdominant source of contamination because they are
much less correlated with tSZ sources (Sehgal et al. 2010), and
the primary CMB is not expected to bias our results as it is
uncorrelated. To test the impact of the CIB contamination, we
make use of an additional map described in Madhavacheril
et al. (2020), made with a constrained ILC algorithm in which
the frequency-combined map is required to have a null
response to the CIB in addition to having unit response to
the Compton y. For our key results in this paper (Section 6) we
run our pipeline with this CIB-deprojected map and find that
it biases the signal slightly lower compared to the results with
the original map. The signal is reduced by ~6%, which
corresponds to only ~12% of the original 1σ error bars.
Because it is a negligible effect compared to the errors, we
choose to present results using the less noisy map without CIB
deprojection.

In addition to the residual contamination from non-y
components, the y-map contains stripey instrument noise,
which is oriented along the ACT scan direction and
uncorrelated with any structure in the Galactic or extragalactic
sky. We emphasize that this does not affect the oriented
stacking procedure, as orientations are entirely determined by
the galaxy data described in Section 2.3 and collected by a
different survey (DES).

A complete description of the y-map procedure is found in
Madhavacheril et al. (2020) and the maps are publicly available
on the NASA Legacy Archive Microwave Background Data
Analysis (LAMBDA; Miller & LAMBDA Group 2018).<sup>80</sup>
Figure 2 shows the sky area of the map (outlined in black)
overlaid on a map of the positions of galaxies (red points) and
galaxy clusters (black points) from the DES. We discuss these
cluster samples below.

### 2.2. Galaxy Cluster Data

We stack cutout images from the ACT Compton-y-map on
locations of galaxy clusters identified in the DES Y3 data. DES
(The Dark Energy Survey Collaboration 2005) recently
completed a 6 yr survey (2013–2019) of 5000 deg² of the
southern sky in five optical filters (grizY). The survey was
conducted by the 4 m Blanco Telescope, fitted with the Dark
Energy Camera (Flaugher et al. 2015), at the Cerro Tololo
Inter-American Observatory in Chile. The redMaPPer algo-
rithm, originally introduced in Rykoff et al. (2014), identifies
galaxy clusters by searching for overdensities of red galaxies.
The first redMaPPer catalog for DES was published using the
Science Verification Data (Rykoff et al. 2016). The algorithm
determines a value for richness, λ, for each cluster by summing

<footnote n="79">⁷⁹ https://docs.astropy.org/en/stable/cosmology/index.html#module-astropy.cosmology</footnote>

<footnote n="80">⁸⁰ https://lambda.gsfc.nasa.gov/product/act/act_dr4_derived Maps_info.cfm#compsep</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Selected redMaPPer clusters (black points), redMaGiC galaxies (red points), and the outline of the ACT Compton-y-map (black line). Clusters are cut from the edge of the y-map and galaxy footprint to avoid edge effects from the y-map when stacking, and to ensure that there is ample galaxy information surrounding each cluster.</figcaption></figure>

the membership probability of each galaxy that has some likelihood of belonging to that cluster, within a defined radius. $\lambda$ is therefore related to the mass; a detailed study of the mass-richness relation for DES was done in McClintock et al. (2019). Our study uses a catalog generated from Y3 Gold data from the first 3 yr of the survey (Sevilla-Noarbe et al. 2021). This redMaPPer catalog (titled v6.4.22+2 Full) extends out to $z \sim 1$, includes all $\lambda > 5$ clusters, and provides a photometric redshift (photoz) for each cluster with uncertainties $\sigma_z/(1+z) \sim 0.01-0.02$ (McClintock et al. 2019).

From this catalog, we select only clusters which overlap with the ACT D56 region. We further limit the cluster area by enforcing that clusters must be over 2° inside the edge of D56, as shown in Figure 2 (black points versus black D56 outline). This ensures that no edge effects are present in any of the y-map cutouts. Additionally, we remove clusters that are closer than 1° to the edge of DES galaxy data (discussed in Section 2.3). Thus, every cluster is surrounded by ample LSS information in all directions, necessary for the accurate determination of orientation.

We choose to limit the cluster sample to $\lambda > 10$. This threshold is a trade-off between the disadvantage of small number statistics when imposing a stricter cutoff and the advantages of a higher-richness sample. These advantages include (a) that the anisotropic Compton-y signal from the LSS surrounding higher-$\lambda$ clusters is stronger, and (b) that higher-$\lambda$ cluster data are more pure. Point (a) is determined in Section 5.2 with simulations. As for point (b), past cosmology studies with redMaPPer clusters have almost always excluded clusters with $\lambda < 20$ due to their known impurities (e.g., Abbott et al. 2020; Costanzi et al. 2021; To et al. 2021). These low-$\lambda$ clusters are more likely to be false detections from random fluctuations or line-of-sight projections in the galaxy field. They also suffer from more mis-centering (Rykoff et al. 2016). However, in our study, a $\lambda > 20$ cutoff is not feasible. There are not enough $\lambda > 20$ clusters available in the D56 sky region to achieve a detection of anisotropic thermal LSS, especially given that the Compton-y signals beyond the cluster radius are much weaker than the internal signal, which is often the focus of cluster-stacking research. To choose the specific $\lambda$ threshold, we test $\lambda > 10$, $\lambda > 15$, and $\lambda > 20$ on noiseless simulations (described below). The correlation described in (a) does not provide enough signal boost to offset the increase in random noise from limiting the cluster sample beyond $\lambda > 10$. Since this cutoff results in the highest S/N, we apply it to the real data.

Effects from mis-centering are expected to be negligible in our study, as the center offsets are typically a fraction of the redMaPPer cluster radius $R_\lambda$ (Zhang et al. 2019) and we will examine signals beyond $1.5R_\lambda$. However, the effects of false cluster detections may be non-negligible, and we would expect them to bias our results lower. We discuss this further in the conclusions.

There are 5494 clusters in the remaining sample. The median photometric redshift uncertainty is $\sigma_z/(1+z) = 0.009$. As part of the stacking process described in Section 4, we divide the cluster sample into even slices in comoving distance along the line of sight, which are each 200 Mpc thick. Figure 3 shows the distribution of cluster richness for each distance slice of the $\lambda > 10$ sample. The figure also shows three colored lines which represent the different richness thresholds tested in simulations in Section 5.2.

## 2.3. Galaxy Data

We use galaxy data from DES to orient each Compton-y-map cutout with respect to the axis of the surrounding elongated structure. The redMaGiC algorithm selects luminous red galaxies from photometric surveys using a matched-filter technique (Rozo et al. 2016). It has recently been applied to Y3 DES data over the full DES footprint, as detailed in Pandey et al. (2021). The algorithm is designed to minimize errors in galaxy photoz, resulting in an average scatter of $\sigma_z/(1+z) \sim 0.013$. In this study we use the high-density catalog, which covers redshifts 0.1–0.7 with a fairly consistent number density ($\sim$[8.5, 10.5] $\times 10^4 h^{-1}$ Mpc$^{-3}$). The high-density nature of this catalog is important for accurate orientation. The average halo mass of redMaGiC galaxies is quite large, at $1.5 \times 10^{13} h^{-1}$ M$_\odot$ (h = 0.69, Pandey et al. 2021). In the future, it would be interesting to explore stacking on every redMaGiC galaxy rather than only redMaPPer clusters, which would sample more points in the cosmic web and provide improved statistics.

## 2.4. Buzzard Simulations

For comparison with the observational data, we make use of the Buzzard Version 1 simulations (DeRose et al. 2019, hereafter D19). The Buzzard galaxy and cluster catalogs were explicitly designed to provide a comparison with DES by replicating many of the survey's selection effects; therefore this simulation suite is the most straightforward choice for direct comparison with our ACT $\times$ DES results. Buzzard assigns realistic galaxies to dark matter halos in N-body dark-matter only
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Box plots showing the richness distributions of selected $\lambda > 10$ clusters after separating the sample into 200 Mpc thick slices in line-of-sight comoving distance. The central orange line in each box indicates the median richness in the bin, which is consistently ~15. Each box is drawn from quartile 1 (Q1) to 3 (Q3) of the data, and the whiskers are drawn out to 1.5 × (Q1–Q3) beyond the box on either side. Black circles represent outliers. The richness is skewed toward low values but there are a fair number of high-richness outliers. The higher-richness clusters contribute a stronger tSZ signal from both the cluster itself and surrounding structure; however, lower-richness clusters must be included to achieve a sufficient S/N on extra-cluster gas. In Section 5.2, we use simulations to test three richness cutoff values (10, 15, 20), shown here in blue, purple, and green.</figcaption></figure>

simulations using the ADDGALS method (Wechsler et al. 2021). The implemented cosmology has $\Omega_m = 0.286$, $h = 0.7$, and $\sigma_8 = 0.82$. This is slightly different from the Planck Collaboration et al. (2016c) cosmology, which we use for radial and transverse distance calculations; however, the differences contribute at most a ~2% error to the rescaling of cutout images, which is later described in Section 4. The galaxy catalogs are post-processed by the DES pipelines to create mocks of the redMaPPer and redMaGiC catalogs in the same footprint with the same selection effects. We apply the same redshift selections to Buzzard as we do to the observational data.

The redMaGiC mock galaxies approximate the clustering in the real catalog well as shown in D19. However, the simulation struggles to match the redMaPPer observables from DES. Key differences between the mass-richness relation and the cluster abundance in Buzzard and DES are shown in Figures 12 and 13 of D19. In Buzzard, the number of identified redMaPPer clusters is a factor of 3–5 below the number of real redMaPPer clusters for the redshifts and richnesses used in our work. This deficit is likely due to a reduced number of galaxies in the central regions of galaxy clusters and a reduced number of red galaxies in dense regions, both of which result in fewer richness selected clusters (D19, Wechsler et al. 2021). Despite the lower cluster number density, the number of Buzzard clusters available to use for the noiseless theory comparison is ~5 times larger than it is for DES. This is due to the fact that the Buzzard tSZ, cluster, and galaxy data overlap fully in the entire DES footprint, whereas which the ACT × DES overlap is ~12× smaller. Thus, the random noise in stacks will be lower when using Buzzard. We further discuss the effects of the cluster abundance discrepancy in Section 7.2.

To create a mock Compton-y-map, we paste pressure profiles from Battaglia et al. (2012b) on Buzzard halos from the same simulation run, then convert these to Compton y. This follows the approach of Stein et al. (2020), described in more detail in

Section 2.5. We apply the model down to halo masses of $10^{12} M_{\odot}$, although halos at such low masses contribute very little to the overall tSZ signal. The result is a noiseless projected y-map. The map contains only signal from halos; because Buzzard is not a hydrodynamic simulation there is no prescription for ejecting baryonic material out of halos into the WHIM. Therefore, filaments in Buzzard only contain bound gas in halos. We then convolve this map with a 1.6 beam and remove modes with $l < 20$ to match the filtering of the ACT Compton-y-map.

Ideally, we would create many noisy versions of the Buzzard-y map by combining many realizations of the Buzzard simulation with many realizations of simulated ACT noise. Such maps would be useful to assess uncertainties and provide a direct comparison between simulations and data. This is unfeasible because there is only one readily available realization of the Buzzard simulation which has had all the post-processing steps applied to create galaxy and cluster catalogs as well as a y-map. Also, the relevant ACT noise has only been simulated in a small fraction of the Buzzard sky footprint. Instead, we estimate uncertainties by using spatial splits of the single Buzzard realization combined with many ACT noise simulations in D56; details are further described in Section 4.5. We also use Buzzard for a noiseless comparison to ACT results.

The Buzzard algorithm was recently improved in v2.0 for validation of the DES Y3 results (DeRose et al. 2022); our follow-up paper will use the state of the art for DES mock simulations.

## 2.5. Websky Simulations

We wish to use simulations not only to make direct comparison with ACT × DES, but also to generate pure theory expectations for the superclustering of thermal energy over time. The Buzzard redMaGiC and redMaPPer mock catalogs are not ideal for studying pure theory expectations because of the DES observational limits applied to them, i.e., the limited redshift range and sky coverage. The Buzzard suite also includes more extensive catalogs without DES limitations: however, we choose instead to make use of the full-sky Websky Extragalactic CMB Simulations (Stein et al. 2020, hereafter S20). The notable speed of the Websky algorithms will be useful in future work for generating alternative cosmologies to which we will apply oriented stacking.

The halo catalogs for Websky were produced via the Peak Patch algorithm, which rapidly generates halos from an initial density field using an ellipsoidal collapse model (Bond & Myers 1996; Stein et al. 2019), excluding overlapping halos in the final list. Websky extends out to $z = 4.6$ and was run with Planck 2018 cosmology: $\Omega_m = 0.31$, $h = 0.68$, and $\sigma_8 = 0.81$. These parameters are slightly different than those used for Buzzard, but this is unimportant as we will not directly compare the two simulations.

The Peak Patch halo catalogs were transformed to sky maps of various probes of structure, including the tSZ effect. For the tSZ effect, the Websky simulation assigns each halo a thermal pressure as a function of redshift and mass. This pressure profile relationship was determined by the stacking of halos found in the LSS simulations of Battaglia et al. (2012b), which focused on clusters and galaxy groups and included AGN feedback. The exact prescription in Websky is given in Equation (3.12) of S20. The Websky Compton-y power spectra have been validated against multiple observational data sets.
In this paper, we will stack the Websky-y-map on halo positions extending to $z \sim 2$, which allows us to go well beyond the DES redMaGIC limit of $z \sim 0.7$. Because the full Websky halo catalog extends further (to $z = 4.6$), for any stack, there are ample higher-redshift clusters contributing to a realistic uncorrelated tSZ background. As with the Buzzard map, we convolve the Websky-y-map with a 1.6 beam and remove all power for modes $l < 20$.

We select cluster-mass halos in a few mass ranges and use mass-weighted halos to provide information about LSS, as further detailed in Section 5.

### 3. Galaxy Field Characteristics

We begin by studying the large-scale properties of the projected galaxy number density field. The bulk of this study is implemented with the Cosmology Object Oriented Package (COOP)81 Huang 2016). Later, these properties will be used to constrain the redMaPPer cluster sample and thus limit the stacks of y-map cutouts to special locations in the cosmic web. In particular, we are interested in the gas signal from superclusters, where the tSZ contribution from high- and low-mass clusters, lower-mass groups, and shock-heated gas in filaments and other compressing configurations should be stronger than that of average aligned structure. Superclusters correspond with overdense and elongated regions of coarsely smoothed galaxy maps (Oort 1983; Einasto et al. 1997; van de Weygaert & Bond 2008b).

We search for regions that satisfy these criteria using projected galaxy overdensity maps in bins of redshift. The overdensity is defined as $\delta_g = n_g/\bar{n}_g - 1$, where $n_g$ is the two-dimensional number density of galaxies and $\bar{n}_g$ is its mean. In practice, each map is created in the Healpix 82 pixelization scheme through the Python package Healpy (Górski et al. 2005; Zonca et al. 2019). Each $n_g$ map is created with NSIDE = 4096 by adding 1 to the appropriate pixel for every galaxy within the redshift bin. After transforming to $\delta_g$, we smooth the map by convolving with a 2D Gaussian function, creating a smoothed map $F_g$. Various choices of smoothing kernel would be valid, and the top-hat function is another that is frequently used in the LSS literature; we choose a Gaussian as it is most conveniently implemented in COOP. The Gaussian filter scale $R_G$ is related to the FWHM by $FWHM = 2\sqrt{2\ln 2}R_G$. We vary the smoothing scale to observe the LSS properties at a range of scales. Specifically, we examine results for Gaussian smoothing with the FWHM ranging from 6–18 Mpc. The approximate equivalent range in top-hat radius $R_{TH}$, if the maps were smoothed with a top-hat filter to produce similar field properties, is $R_{TH} \sim 5–11$ Mpc. The conversion is done by enforcing an equal volume under the top-hat function with radius $R_{TH}$ and height 1, and a Gaussian function with amplitude 1. The range is chosen to examine highly nonlinear structure beyond a typical cluster radius; smoothing scales are discussed further in Section 4. In the current section, we use an FWHM of 14 Mpc for demonstration purposes.

The primary property we use to determine superclustering is the field excursion $\nu$,

$$ \nu = F/\sigma(R), \qquad (2) $$

where $F$ is the field value at some position and $\sigma$ is its rms after smoothing on some scale $R$ (Bardeen et al. 1986). Given that the $\delta_g$ maps have a mean of zero, $\sigma$ is equivalent to the standard deviation, making $\nu$ a measure of signal-to-fluctuation noise. Points with higher $\nu$ correspond to rarer overdensities.

We also use the asymmetry of the field as a metric of superclustering. To measure the alignment and elongation at any point in the field, it is natural to consider using the tidal field, because the tidal shear in the early universe is key to producing filaments (van de Weygaert & Bond 2008b). In addition, it has the same power spectrum as the density. The derivative of a Gaussian-smoothed tidal field with respect to scale $R$ is the Hessian matrix of the Gaussian-smoothed density, which has frequently been used to characterize cosmic web phenomenology (as reviewed in Libeskind et al. 2018). The Hessian is defined in 2D as

$$ H = \begin{bmatrix} \frac{\partial^2 F}{\partial x^2} & \frac{\partial^2 F}{\partial x \partial y} \\ \frac{\partial^2 F}{\partial y \partial x} & \frac{\partial^2 F}{\partial y^2} \end{bmatrix} \qquad (3) $$

for a field $F$, evaluated at some point. We use the Hessian to determine the asymmetry and alignment of the smoothed projected galaxy overdensity field at cluster positions. The choice of smoothing scale defines a characteristic radius from each selected field point at which the Hessian encompasses maximal information; this can also be thought of as a shell region of the tidal field.

We adopt the notation and conventions of Bond & Efstathiou (1987, hereafter BE87), in defining dimensionless eigenvalues of the Hessian. BE87 first applied the study of these field properties to the CMB. At any field point, the Hessian has eigenvalues $-\lambda_i$ and corresponding eigenvectors. Note that with the negative sign, $\lambda_i$ are defined to be positive at peaks and negative at troughs. We order the eigenvalues as $|\lambda_1| > |\lambda_2|$, such that $\lambda_1$ corresponds to the eigenvector along which curvature is changing most rapidly: the short axis of curvature. Using the eigenvalues, we can define the ellipticity $e$:

$$ e = \frac{\lambda_1 - \lambda_2}{2(\lambda_1 + \lambda_2)}. \qquad (4) $$

This follows the definition in BE87, which succeeded the 3D representations for galaxy fields in Bardeen et al. (1986, hereafter BBKS). The numerator describes how elongated the field is at a certain point by the difference in eigenvalues there. This is normalized by the trace of the Hessian to provide an equitable comparison between different-amplitude peaks/ troughs. The $\nu$ and $e$ parameters are illustrated in Figure 4.

In BBKS, the ellipticity was applied in the context of cluster-scale smoothing to the 3D galaxy density field. In this context, clusters lie at peaks in the field. At a location with negative curvature in both directions, such as a peak, the minimum value for $\lambda_2$ is 0, and consequently $e \le 0.5$. However, for the larger scale (e.g., 14 Mpc) smoothing used in this work, most clusters do not lie at peaks in the projected, smoothed galaxy overdensity maps. Many clusters exist where the gradient of $F_g$ is nonzero. These large-scale gradients point toward supercluster centers rather than toward the individual clusters that make up each supercluster. For some clusters, where the

<footnote n="81">81 https://www.cita.utoronto.ca/~zqhuang/work/coop.php</footnote>

<footnote n="82">82 https://healpix.sourceforge.io/</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Illustration of $\nu$ and $e$ at cluster locations overlaid on a smoothed, projected redMaGC galaxy map $F_g$. A $3.5^\circ \times 3.5^\circ$ patch of $F_g$ is shown, corresponding to $\sim 100$ Mpc at the center of the galaxy redshift bin ($0.35 < z < 0.5$). The smoothing was Gaussian with an FWHM = $28'$ ($\sim 14$ Mpc). redMaPPer cluster locations are identified by white points. $\nu$ and $e$ are calculated at the cluster locations using $F_g$; because $F_g$ has been coarsely smoothed, they are properties of the LSS in which each cluster is embedded. Black circles have an area scaled by the $\nu$ value at the cluster point. Headless vectors have a length corresponding to $e$ and are oriented along the eigenvector $\nu_2$ of the Hessian (Section 4). We define supercluster regions as points where $\nu$ and $e$ are both large.</figcaption></figure>

field has two equal-sign eigenvalues, $e \le 0.5$. Many clusters also lie at regions where the eigenvalues have opposite signs (which are saddle points if the background gradient is 0), where $e > 0.5$. Thus, the distribution of $e$ values, while concentrated toward 0, extends well beyond 0.5 and can become very large when the denominator of Equation (4) is small (in other words, when $\lambda_2 \sim -\lambda_1$). We find that points with higher $\nu$, at rarer overdensities, are more likely to have $e < 0.5$. Later, we choose to apply a minimum $e$ threshold to select for locations of the field that are highly elongated, which is a feature of superclustering. However, we do not limit our sample by applying a maximum $e$ threshold; $e$ is allowed to be arbitrarily large.

We also consider the field curvature excursion, $x$, related to the trace of the Hessian:

$$x = \nabla^2 F / \sigma_2, \qquad (5)$$

where $\sigma_2$ is the rms value of $\nabla^2 F$. This property was defined in Bond & Efstathiou (1987); our definition differs by an absolute value sign such that our $x$ is allowed to be negative. Points with high $x$ have large curvature in one or both directions. $x \propto (\lambda_1 + \lambda_2)$, which is the denominator of $e$, so points where $\lambda_1 \sim -\lambda_2$ have small $x$ and large $e$. Due to the divergence of $e$ at $x \sim 0$, we choose to show the non-normalized elongation $|x|e$ in upcoming figures for visual purposes.

To better understand these characteristics and how they may be used to find regions of strong superclustering in the universe, we compare the $\nu$, $|x|e$, and $x$ distributions of the galaxy overdensity field with that of a GRF we have constructed with the same power spectrum. We divide the Buzzard redMaGC galaxy catalog into 3 redshift bins: $0.15 < z < 0.35$, $0.35 < z < 0.5$, and $0.5 < z < 0.65$. Using

measurements of the Buzzard galaxy autospectrum, $C_l^{gg}$ (Pandey et al. 2021), we generate a GRF realization from the power spectrum by using the Synfast function from Healpy. We will refer to this map as the pseudo-galaxy map.

As we are interested in properties of the Buzzard field at the locations of mock redMaPPer clusters, it is necessary to identify redMaPPer-like peaks in the GRFs. Accordingly, we identify peaks in the field of a similar angular size as redMaPPer clusters for each $z$ bin. We determine the size by finding the approximate mass $M_{200m}$ of a $\lambda = 15$ cluster using the mass-richness relation from McClintock et al. (2019). $M_{200m}$ refers to the mass enclosed within a sphere of radius $R_{200m}$ within which the density is, on average, $200 \times$ the mean matter density of the universe at that redshift. $\lambda = 15$ was chosen because it is the median richness of our sample. Next, we convert $M_{200m}$ to $R_{200m}$ and find the Gaussian filter equivalent of a top-hat filter with that radius. (For a discussion on top-hat to Gaussian conversion, see Section 4.2.) We smooth the pseudo-galaxy field with a Gaussian function with an FWHM $\sim 1.6$ Mpc. The equivalent angular size varies for each redshift bin. For the smoothed field, we find all peaks and sub-select them by applying a $\nu$ threshold on the small peak scale. We choose a $\nu$ threshold for each $z$ bin which results in the same number of GRF peaks as Buzzard mock redMaPPer clusters. In summary, we perform approximate abundance matching to find GRF peaks which are similar in size and amplitude to redMaPPer clusters.

Next, we examine the distributions of $(\nu, x, |x|e)$ for the galaxy overdensity fields smoothed on larger scales. Figures 5 and 6 show the distribution for $|x|e$ versus $\nu$ and $|x|e$ versus $x$, respectively, for the galaxy/GRF field smoothed at a 14 Mpc scale, at the chosen cluster/peak positions in the middle $z$ bin.

In the real and mock galaxy fields, the clusters (red + blue points in the top two panels) display a $\nu$ and $x$ distribution, which is more stretched compared to the GRF. The distributions are both skewed toward the high end, displaying more high-$\nu$ and high-$x$ points than the GRF, which by definition has no skew. Additionally, the field elongation extends to larger values in the real and mock data than in the GRF. We examine these distributions for a few galaxy field smoothing scales. For finer-grained smoothing, the distribution of GRF peaks is roughly the same but the skew of the real and simulated galaxy fields increases, especially in $\nu$. This is expected because with finer smoothing, the galaxy field is more non-Gaussian due to nonlinear structure formation. As the smoothing scale becomes coarser, the galaxy fields approach the GRF result.

We later demonstrate, in Section 5.2, that enforcing a minimum $\nu$ and $e$ threshold for the redMaPPer cluster sample enhances the supercluster gas signal in stacks. We choose $\nu > 2$, $e > 0.3$ as the optimal cuts (justified in Section 5.2). Clusters satisfying this constraint are shown in blue in Figures 5 and 6. This selection furthers the distinction between the GRF and galaxy overdensity fields, as remaining points are more concentrated in the GRF in all properties. The field constraints select for clusters in highly overdense, elongated regions of coarse-grained galaxy maps: effectively, regions of strong superclustering. In Section 4.6, we examine the effect that these constraints have on the non-Gaussianity of oriented stacks.

We also find that the $\nu$ and $e$ properties on the 14 Mpc scale are not highly correlated with cluster richness $\lambda$, a small-scale
<figure><img src="image_8.png" /></figure>

**Figure 5.** Distribution of $|x|e$ (field elongation) vs. $v$ (field excursion) for ~9500 points in the redshift bin 0.35 < z < 0.5. From top to bottom, the distribution is shown for real redMaPPer cluster positions in the redMaGiC galaxy overdensity field, the same for Buzzard mocks of the DES data products, and the GRF with the same galaxy power spectrum ($C_{\ell}^{gg}$) as Buzzard. The maps have each been Gaussian smoothed with an FWHM = 14 Mpc. The entire cluster/peak sample consists of the red plus blue dots. The $1\sigma$, $2\sigma$, and $3\sigma$ contours are drawn in black, and the points which are included after our chosen cuts of $v > 2$ and $e > 0.3$ are colored in blue. The cutoff is sharp at $v = 2$ but less visually apparent for $|x|e$ because the threshold is in $e$ alone. These cuts select for high-superclustering regions in the real and mock galaxy fields.

property. Figure 7 demonstrates this: despite the weak correlation, the remaining clusters after the cut $v > 2$ are still distributed across a wide range in $\lambda$. A higher-richness cluster is more likely to be in a higher large-scale overdensity, but low-richness clusters may also belong to such regions. For example, a small cluster may lie on the edge of a supercluster and thus have a high $v$ value. For $e$ (not shown), there is a weak anticorrelation: clusters in very high-ellipticity regions tend to be of lower richness. This suggests that these small clusters are more likely to lie at saddle points in the field, such as between two massive overdensities, where $e$ is large due to its normalization. There is no correlation in the non-normalized $|x|e$ property.

# 4. Stacking Methods

The stacking procedures were originally developed and implemented as part of the Cosmology Object Oriented Package (Huang 2016), which we use in our stacking pipeline. This work is its first application to low-redshift objects (galaxies and clusters), whereas the program had been previously used for primary CMB analysis.

## 4.1. Oriented Stacking with ACT and DES

In individual cutouts of clusters in the ACT-tSZ data, only the most massive clusters are easily detectable by eye because the map is noise dominated. The tSZ signal of lower-mass

<figure><img src="image_9.png" /></figure>

**Figure 6.** Distribution of $xe$ (field elongation) vs. $x$ (S/N of the curvature) for galaxy overdensity maps from DES observations (top), Buzzard mocks (middle), and a GRF (bottom). The chosen cuts that are used throughout much of this paper, $v > 2$ and $e > 0.3$, are shown in blue. The $v$ cut shifts the $x$ distribution rightward because rarer overdensities also tend to have higher second derivatives. The $e$ threshold shows up here as a diagonal $|x|e$ threshold.

<figure><img src="image_10.png" /></figure>

**Figure 7.** The relation between $v$, as measured on the galaxy overdensity field smoothed at 14 Mpc, and cluster richness $\lambda$. There is slight correlation, such that the $v > 2$ clusters have a higher-skewed $\lambda$ distribution, but the constrained sample still features a range of $\lambda$.

clusters, groups, and galaxies lies well below the noise and only emerges through stacking. In addition, stacking averages over the diversity of distributions and shapes of thermal energy along cosmic filaments. While the physics of individual superclusters is interesting in its own right, our goal is to measure the ensemble average of the anisotropic clustering of thermal energy around galaxy clusters.

Therefore, for each selected cluster sample, we stack cutout images from each Compton-y-map with orientation. Oriented stacking aligns and combines the gas signal from the most massive extended structures surrounding each cluster while driving down the noise (Figure 1).
<figure><img src="image_11.png" /><figcaption>**Figure 8.** The number of redMaPPer clusters (top) and projected comoving number density of redMaGiC galaxies (bottom) in each 200 Mpc wide bin. Gray-shaded regions cover data not used in this work. The total cluster number per bin is relevant because the S/N scales as $\sqrt{N_{cl}}$ for a stack. $N_{cl}$ generally rises with redshift because the same sky area covers a progressively larger physical area. For galaxies, the comoving number density is the relevant quantity, as the available LSS information around each cluster determines its orientation. The number density is relatively constant until $z \sim 0.7$, the imposed edge of this volume-limited sample. We cut the data at the drop-off. We also cut the lowest-redshift data because $4^\circ \times 4^\circ$ cutouts do not span a large enough physical extent to observe the full effects from superclustering at these redshifts.</figcaption></figure>

Applying the method to our observational data begins by dividing the selected redMaPPer clusters and all redMaGiC galaxies into equal-sized slices in comoving distance along the line-of-sight direction. Each slice is 200 Mpc thick. The thickness is chosen to minimize projection effects from uncorrelated structure, yet ensure that most clusters and galaxies are placed in the correct redshift slice given that the DES photoz uncertainties are $\sigma_z \sim 0.01 - 0.02$ ($\sim$30–60 Mpc at the redshifts of interest). The distributions of the cluster and galaxy data within the slices are shown in Figure 8.

In each slice, we create a smoothed projected galaxy overdensity map as described in Section 3. The local Hessian matrix (Equation (3)) of the galaxy overdensity at the position of each cluster provides information on the strongest axis of anisotropic clustering. The eigenvector $v_2$ with eigenvalue $\lambda_2$ points along the axis of slowest change, which we define to be the superclustering axis. The choice of smoothing scale defines the characteristic radius at which this axis is identified; thus, for the same galaxy map smoothed at different levels, the

eigenvector for a given cluster can rotate. The vector will point along near-cluster structure, intercluster filaments, and superclusters as the map is increasingly smoothed.

For each cluster, we take a $4^\circ \times 4^\circ$ cutout from the ACT Compton-y-map centered on the cluster (R.A., decl.). At the redshift range explored in this work, this corresponds to a coverage of $\sim$80 Mpc (closest slice) to $\sim$120 Mpc (furthest slice) on each side of the square cutout. Thus, if the central cluster is a member of a supercluster, the cutout should almost always contain the entire structure (see Borgani 1995 for a detailed discussion of supercluster scales). The square cutout is oriented along the superclustering axis, and each cutout is rotated so as to align the superclustering axis of all cutouts along the horizontal axis of the stacked image. After alignment, the final stack is the average of all the cutouts. The S/N of the final image is proportional to $\sqrt{N}$, where N is the total number of stacked cutouts. The cutout, orientation, and stacking processes are implemented with the *GetPeaks* and *Stack* programs from the Cosmology Object Oriented Package<sup>83</sup> (Huang 2016).

## 4.2. Choice of Smoothing Scale

Multiple smoothing scales are explored in Sections 5.2 and 6; they span a range of Gaussian FWHM from 6–18 Mpc. The upper end of this range is motivated by recent cluster-pair stacking studies that revealed filament signal from pairs separated by a transverse distance of 6–10$h^{-1}$ Mpc (Tanimura et al. 2019). Considering that our stacking method centers on single clusters, the simplest way to enforce orientation at these typical scales would be to apply top-hat smoothing to the galaxy overdensity map with a radius of $R_{TH} \sim 6–10h^{-1}$ Mpc. A simple conversion from top-hat to Gaussian filter involves equalizing the volume under each function; with a top hat of height 1 and a Gaussian with an amplitude of 1, this leads to an FWHM $\sim 1.7R_{TH}$. Thus, for cluster pairs separated by, e.g., $\sim 7h^{-1}$ Mpc, Gaussian smoothing for which the Hessian would best encode alignment at that scale has an FWHM $\sim 12h^{-1}$ Mpc or $\sim 18$ Mpc.

We therefore choose 18 Mpc as our key scale to study with ACT $\times$ DES data. We also select a range of smaller smoothing scales, down to an FWHM = 6 Mpc, to apply to both observational and simulated data in order to study aligned superclustering on highly nonlinear scales. Finer-grained smoothing causes the orientation to be determined primarily using intercluster filament galaxies and groups.

For each chosen comoving smoothing scale, the corresponding angular size is determined at the center of each slice. The angular size varies from map to map such that the FWHM in comoving megaparsec is held constant across all redshifts.

## 4.3. Redshift Slice Combination

Ideally, it would be most interesting to observe the change in the average superclustering signal with cosmic time by comparing multiple redshift slices. We explore this in simulations in Section 5; however, we find that the S/N in the real data analyzed here is too low for each individual 200 Mpc slice. For the observed data we therefore combine stacks over nearly the full range of redshifts available with the DES data: $\sim 0.25 < z < 0.72$, or $\sim 1000–2600$ Mpc.

<footnote n="83">83 https://www.cita.utoronto.ca/~zqhuang/work/coop.php</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 9.** The ingredients for making an oriented stack in the slice [1632, 1832] Mpc ($z = [0.41, 0.47]$). Left: a 4°–4° sky patch from the observed (above) and simulated (below) Compton-y-maps. At this redshift, the corresponding physical size is ~180 Mpc. redMaPPer cluster locations within the slice are circled in white. The ACT y-map is noise dominated, whereas the Buzzard simulated map is noiseless. Both y-maps display the projected tSZ signal from all contributing sources along each line of sight. Because of this projection, as well as the dependence of tSZ strength on cluster mass, not all circled clusters are easily detectable by eye. Middle: observed (above) and simulated (below) galaxy overdensity maps, Gaussian smoothed with an FWHM = 14 Mpc, in the same sky area. The same clusters are circled and black arrows indicate their orientation with respect to the surrounding LSS, as determined by the Hessian matrix on $F_g$. Right: oriented stacks made by a combination of multiple distance slices in the range of 1032–2632 Mpc. Each stack is made by taking cutouts from the y-map centered on cluster locations and rotating it with the information from $\delta_g$ before stacking. This builds up S/N along the horizontal axis, showing the gas signal from superclustering.</figcaption></figure>

Each stack is the same angular size (4° × 4°) by construction, and therefore stacks on structure in more distant slices span a larger transverse comoving distance. To properly combine multiple slices at different redshifts, we must first adjust all images to the same comoving size. We begin by calculating the transverse comoving size of the nearest image at the midpoint of the slice in the line-of-sight direction (~1100 Mpc). At this distance, 4° spans ~80 Mpc. For more distant slices, the stacks are cropped to the angular size, which spans 80 Mpc at each slice midpoint and is rescaled to the same pixelization via interpolation. Finally, all images are averaged together. By using the slice midpoints for the conversion from angular to transverse comoving size, we make the approximation that the slices are very thin, while in reality the conversion varies across the 200 Mpc slice thickness.

Figure 9 displays both the observational and simulated ingredients and outputs of the stacking process.

## 4.4. Multipole Decomposition

To quantitatively compare the stacked images, we decompose each image into its multipole components. Each image $I$ can be deconstructed as

$$I(\theta, r) = \sum_{m} (C_{m}(r)\cos(m\theta) + S_{m}(r)\sin(m\theta)). \quad (6)$$

Since symmetry along the x-axis is enforced in the oriented stacking method, the odd components trend toward zero as the number of component images grows. The odd moments are consistent with zero in our results. For even $m$, due to the alignment along the x-axis, we are interested in only the cosine

component. The sine term is like a noise term and fluctuates around zero when the number of stacked images is large. The radial profile of the cosine component is taken by

$$C_{m}(r) = \frac{1}{X\pi} \int_{0}^{2\pi} d\theta F(\theta, r)\cos(m\theta), \quad (7)$$

where $X$ is 2 for $m=0$ and 1 for all other $m$. The multipole decomposition is illustrated in Figure 10. Throughout the rest of the paper, we focus on a comparison of $C_2(r)$ and $C_4(r)$.

All $m=2$ and $m=4$ profiles shown in this paper feature a rise, peak, and fall. This feature is almost entirely dependent on the smoothing scale, as we later demonstrate in Section 5.2. In short, determining the orientation of each cluster with a galaxy overdensity map smoothed at a chosen scale enforces a radius at which structure is maximally aligned between all stacked images, translating to a peak at that radius.

Because the $m=4$ moment sums signal not only from the horizontal image axis, which contains signal by construction, but also from the vertical axis, which contains only noise, $C_4(r)$ is expected to be noisier than $C_2(r)$.

## 4.5. Uncertainties

Uncertainties in the stacked profiles are expected to stem from a variety of sources, including random noise, contamination of the y-map from other components such as dust and the primary CMB, and photometric redshift uncertainties in the DES data.

We characterize the uncertainties by estimating the covariance matrix, $\Sigma$, of each $C_m(r)$ profile. If a profile is binned into radial bins $r_i$, each element $\Sigma_{ij}$ is the covariance between the signal in
<figure><img src="image_13.png" /><figcaption>Figure 10. A representation of the three lowest even multipolar moments of an oriented stack from the Buzzard-y-map (left). The stack image is passed as I through Equation (7) to get C_m(r). The odd m moments trend toward zero with increased numbers of stacked cutouts. For the even moments, we show the symmetrized representation C_m(r)cos(mθ). Power in the m = 4 moment comes from the horizontal wings of the stacked image, but gets evenly redistributed to four poles in this representation.</figcaption></figure>

the ith and jth radial bins, and the elements along the diagonal of the matrix are the variances of the signal in each bin. We attempt two different methods to estimate Σ. Method 1 is sufficient for the simulated data but insufficient for the real data. In method 1, we split the clusters into N_reg separate regions on the map where N_reg ranges from 12–48, depending on the data set and cluster selection. Cluster samples in all the regions are approximately equally sized and collected in (R.A., decl.) space into patches 30–80 deg² in area, depending on N_reg. The clusters in each region are stacked with orientation, where the orientation is given by the full F_g map such that information beyond the region edges can be incorporated. Splitting the data into spatially separate regions is motivated by the fact that our measurement contains both spatially correlated noise from LSS fluctuations (data and simulations) and long-wavelength noise from residual low-ℓ primary CMB contamination (data only). With large-area regions, each sample is reasonably independent of the rest, save for some inevitable LSS and long-wavelength noise overlap between neighboring regions.

The number of regions and the angular area per region depend on the data set. For Websky full-sky maps, we split the sample into 48 Healpix pixels. For the Buzzard and DES data, we use the kmeans_radec algorithm<fnref n="84" /> to split the clusters into approximately equal-sized regions on the sky. The full λ > 10 cluster sample is split into 16 regions in DES and 48 regions in Buzzard. Adding in the field constraints ν > 2 and e > 0.3 shrinks the number of clusters, so when stacking this constrained sample we reduce the number of splits to 12 for DES and 24 for Buzzard.

After stacking the clusters in each region, we have N_reg stacks. We decompose each stack into multipole moments m and measure the Compton-y radial profiles C_m(r). Each profile is binned in radius to make a data vector C_m = C_m(r_i), where i ∈ (1, N_bin). We combine the N_reg profiles into an N_bin × N_reg matrix, called X. In calculating the covariance matrix, we weight each region by w_p = N_{cl,p} / N_{cl}, where N_{cl,p} is the number of clusters in the pth region and is the average over all regions. If X becomes the modified matrix M by subtracting the weighted average of each bin across all regions (i.e., the row mean), so that element M_{ip} is given by

$$M_{ip} = X_{ip} - \frac{\sum_{p'=1}^{N_{\text{reg}}} X_{ip'} w_{p'}}{\sum_{p'=1}^{N_{\text{reg}}} w_{p'}}, \quad (8)$$

then the covariance matrix element between bins r_i and r_j is

$$\Sigma_{ij} = \left[ \frac{\sum_{p'=1}^{N_{\text{reg}}} w_{p'}}{\left(\sum_{p'=1}^{N_{\text{reg}}} w_{p'}\right)^2 - \sum_{p'=1}^{N_{\text{reg}}} w_{p'}^2} \right] \sum_{p'=1}^{N_{\text{reg}}} (w_{p'} M_{ip'} M_{jp'}). \quad (9)$$

The left bracketed term is a normalization by the degrees of freedom (dof) that would simply be 1/(N_reg - 1) if there were no observation weights w_p. Expressing the rightmost unbracketed sum in words, the weight w_p is applied to every radial bin in the zero-mean profile of the pth region, and the resulting matrix is multiplied by the transpose of M (generally, (M × M^T)_{mn} = ∑_l M_{ml} M_{nl}). In practice, we calculate Σ with the NumPy cov<fnref n="85" /> function. Σ as defined above represents an estimate of the covariance of a single region's y profile if it had been observed N_reg times, so we further divide Σ by N_reg to achieve our estimate of the covariance of the full map data.

The map-split method works sufficiently well for the Buzzard and Websky simulations, which both can be split into a larger number of regions and have more clusters per region than in the real data. However, we find that there are not enough possible subregions of the smaller-footprint ACT y-map to achieve convergence of the covariance matrix. With only 12 regions and 3–5 radial bins, we find that the position and number of bins significantly affect the resulting χ² and S/N estimates.

Due to the lack of convergence, we apply a different method to the observed data. In method 2, we assume that each covariance matrix for the final ACT y profiles can be decomposed into the sum of components:

$$\Sigma_{\text{tot}} = \Sigma_{\text{noise}} + \Sigma_{\text{signal}}, \quad (10)$$

where Σ/noise refers to the covariance matrix from all non-signal components in the y-maps, and Σ/signal refers to the covariance matrix, which would emerge from a noiseless-y-map.

To estimate the former, and dominant, source of error (Σ/noise), we use 120 simulated ACT Compton-y-maps described in Section VII of Madhavacheril et al. (2020). Each map contains an independent realization of a Gaussian y signal generated with a tSZ power spectrum as well as independent realizations of the estimated noise contribution from all known contaminants. Each map covers the D56 footprint. The contaminants include the low-ℓ primary CMB, high-ℓ

<footnote n="84">⁸⁴ https://github.com/esheldon/kmeans_radec</footnote>

<footnote n="85">⁸⁵ www numpy.org/doc/stable/reference/generated/numpy cov.html</footnote>
instrument noise, and high-ℓ residual foregrounds. We subtract the tSZ realization from each map so that only the noise contribution remains. We estimate $\Sigma_{\text{noise}}$ by generating a stack for each different map. The noise in the simulated maps is uncorrelated with true cluster and galaxy positions, so for each map we stack on a sample of random points, which is the same size as the cluster sample and apply random orientation. We bin the data to sample the radial profile in positions of interest, and the resulting covariance matrix emerges from $N_{\text{map}} = 120$ independent stacks. Because the number of images per stack is constant, Equations (8) and (9) simplify to

$$M_{ip} = X_{ip} - \frac{\sum_{p'=1}^{N_{\text{map}}} X_{ip'}}{N_{\text{map}}}; \qquad (11)$$

$$\Sigma_{\text{noise,ij}} = \frac{\sum_{p'=1}^{N_{\text{map}}} M_{ip'} M_{jp'}}{N_{\text{map}} - 1}. \qquad (12)$$

We estimate $\Sigma_{\text{signal}}$, the covariance matrix of the signal component of the measurement, as being equal to the Buzzard map-splits covariance matrix (calculated by Equation (9)) scaled by $N_{\text{Buzz}}/N_{\text{DES}}$ (the relative number of clusters in Buzzard versus DES). We compute $\Sigma_{\text{signal}}$ for each separate smoothing scale that we apply to the galaxy maps; therefore, it includes the correlation induced between radial bins by smoothing. This is an imperfect estimate of the true signal variances and covariances due to the known inaccuracies in Buzzard's mass-richness relation and our approximate, prescriptive approach to adding gas to the simulation. We compute $\Sigma_{\text{tot}} = \Sigma_{\text{signal}} + \Sigma_{\text{noise}}$ and examine which component contributes more to the total. $\Sigma_{\text{noise}}$ dominates the variances of the signal in each radial bin, as its diagonal is $\sim 3-6 \times$ larger than the diagonals of $\Sigma_{\text{signal}}$. However, the off-diagonal covariances are similar in magnitude. Therefore, although the noise contribution is more important overall, any inaccuracies in the Buzzard simulation may have a non-negligible impact on $\Sigma_{\text{tot}}$. Future work using expanded ACT sky coverage will avoid this concern by estimating uncertainties through data only.

When working with simulations to determine general theory expectations in Section 5, we use arbitrarily fine binning for the radial profiles $C_m(r)$. However, when making robust comparisons between ACT $\times$ DES and Buzzard in Section 6, we choose the binning more carefully. We determine the convergence of $\Sigma_{\text{noise}}$ for $n$ bins by examining the matrix and its inverse as the number of contributing maps increases from 60 to 100 to 120. We also examine the condition number (the ratio of the largest to smallest eigenvalue) of $\Sigma_{\text{tot}}$. We find instability and a high condition number for 5 or more bins when using the constrained cluster sample because the Buzzard covariance matrix was calculated using only 24 separate map regions. Therefore, we select 3 bins for stable results and useful placement along the radial profile. It may be possible to achieve a more stable finely binned covariance matrix with alternative uncertainty methods, but we leave that for future work with larger ACT $\times$ DES sky coverage.

We note that the DES photoz uncertainties are expected not only to contribute to the overall uncertainty of the anisotropic stacked y signal, but also to bias the signal in a redshift-dependent manner. The contribution to uncertainty is accounted for within method 2: $\Sigma_{\text{signal}}$ includes variance across small samples of different photoz realizations for galaxies and clusters within Buzzard. Buzzard reproduces the photoz uncertainties of the DES catalogs well, shown in D19.

The redshift-dependent bias, however, is unaccounted for in the uncertainties. Due to photoz error, for any redshift slice, some objects near the edge are not included and some interlopers just outside the slice are included. This results in a distribution of the true redshifts of objects, which is concentrated toward the middle of each slice with tails extending beyond the slice edges. On average, $\sim 25\%$ of galaxies and $\sim 20\%$ of clusters are interlopers in a given slice. The combination of the peaked distribution and the inclusion of interlopers has a nontrivial effect on the determination of orientation, as it both increases the correlation of objects within the slice but also adds in distant, uncorrelated structure. We test the effect this has on the anisotropic stacked y signal using the Websky simulations. For each redshift slice, we transform the true Websky halo redshifts to photometric redshifts by drawing from a Gaussian-z distribution with the average $\sigma_z$ of galaxies/clusters in DES. We then compare the $m=2$ radial y profiles from a stack using photozs versus a stack using true zs. We find that the profile can be boosted, decreased, and/or shifted in radial units due to the use of photozs, depending on the redshift slice. The most significant effect is in the most distant slice, where there is a decrement of $\sim 25\%$ in the maximum height of the y profile. To make physical inferences from the values of our observational results, these effects must be fully characterized and corrected for, which will be one of the goals of our succeeding paper. However, they do not impact the conclusions of this paper as we limit our comparison to only Buzzard versus ACT $\times$ DES. Due to Buzzard's accurate reproduction of the photoz uncertainties from DES, oriented stacks using Buzzard will be biased in the same manner. We do not attempt to compare Websky directly with Buzzard or ACT $\times$ DES.

## 4.6. Comparison with a GRF

Can oriented stacking be used to examine non-Gaussian structures in the late-time universe that result from nonlinear evolution? In the CDM model, gravitational instabilities cause the primordial matter field, which is Gaussian or very nearly Gaussian, to cluster. At early times, the amplitudes of the matter fluctuations grow linearly at almost all scales. As the universe evolves, small-scale overdensities begin to exceed a critical threshold for collapse. The increase in clustering on small scales causes non-Gaussianity in the matter field and nonlinear deviations in the matter power spectrum. Over time, as the overall power spectrum grows, increasingly larger scales pass the threshold, leading to the collapse of the first stars, then galaxies, then clusters. These nonlinear effects are important at the scales and redshifts studied in our paper, e.g., nonlinearity increases the amplitude of the matter power spectrum $P(k)$ by $\sim 1.5 \times$ at $z=1$ and $k=1 h \text{ Mpc}^{-1}$, and by $\sim 6 \times$ at the same scale at $z=0$ (Springel et al. 2018). Accordingly, we search for a characteristic signal of late-time non-Gaussianity in the stacked maps by comparing the Buzzard simulations to a purely GRF.

We generate GRFs using the measurements of the galaxy autospectrum, galaxy-y cross spectrum, and y autospectrum in the Buzzard simulations (Pandey et al. 2021). These measurements were made for three wide redshift bins in the galaxy simulation, rather than thin 200 Mpc slices, which are used elsewhere. For each bin, the three power spectra contain all necessary information to generate two correlated GRFs, which represent the galaxy field and y field without non-Gaussianities.
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Top: a comparison of the *m* = 0, *m* = 2, and *m* = 4 moments of oriented stacks on Buzzard clusters (solid lines) and GRF peaks (dashed lines); note the three different y-axes. Red indicates that all clusters were used in the stack, whereas blue is the sample that remains after a constraint on *v* and *e* at large scales. The *m* = 0 moment is not well matched near the cluster interior, but at larger radii the Buzzard and GRF profiles are consistent both before and after large-scale field constraints are imposed. For *m* = 2, the red profiles are similar to each other, showing that the quadrupole does not strikingly display signs of non-Gaussianity when the stack includes all clusters. However, the field constraints induce a much larger boost in the Buzzard signal than the GRF (blue dashed vs. continuous line). The *m* = 4 moment is notably consistent with zero for the GRF, whereas in Buzzard the profile peaks near 10 Mpc, which also is boosted with the field constraints. We interpret this sign of non-Gaussianity to be related to the flatness of filaments, arising from nonlinear growth of structure at late times. Bottom: a cutout of the galaxy (Buzzard, left) and pseudo-galaxy (GRF, right) fields that were used for the orientation and field constraints. Both fields are Gaussian smoothed with an FHWM = 14 Mpc. Red circles are plotted around all cluster locations. Many clusters are cut by imposing *v* > 2, *e* > 0.3 and the remaining ones have blue circles overplotted (in other words, blue circles also belong to the red sample). The remaining clusters are themselves clustered. Clumps of blue clusters in the left field are more aligned than in the right. This demonstrates how the *v* and *e* thresholds select clusters in high-superclustering, highly non-Gaussian regions.</figcaption></figure>

Next, we stack on *peaks* in the GRF that roughly match the size of $\lambda > 10$ galaxy clusters in redMaPPer, as described in Section 3. We pass the maps and peak positions through the same oriented stacking pipeline such that cutouts from the pseudo-y-map are rotated based on information from the pseudo-galaxy map, then stacked.

The comparisons between the resulting GRF stacks versus the Buzzard stacks are shown in Figure 11. The lower plots show the smoothed projected galaxy overdensity field for both Buzzard (left) and the GRF (right). These figures provide an illustration of the spatial distribution of the full cluster sample (red circles) and the clusters constrained by $v$ and $e$ limits (overplotted blue circles; all blue circles belong to the red sample as well). The radial profiles for stacks on these points are shown above, for three multipoles, with the same color scheme for the unconstrained and constrained cluster samples. First, we examine the isotropic signal, $m=0$, which is equivalent to what it would be in an unoriented stack. It is consistent beyond ~5 Mpc for both the full (red) and constrained (blue) cluster sample, demonstrating that large-scale differences between the GRF pseudo-y and Buzzard y-maps are indistinguishable from unoriented stacking alone. For $r<5$ Mpc, the Buzzard profiles are significantly higher than the GRF profiles. The discrepancy at this small scale is expected because the realistic Buzzard $\delta_g$ and Compton-y fields have high skew at small scales from the collapse of massive, rare peaks in the density field and the strong Y-M relation for massive clusters. A GRF, by definition, has no skew at any

scale. Thus, by generating the GRF, we redistribute Buzzard power evenly to both low and high values in the pseudo-galaxy and pseudo-y-maps. This redistribution of power causes there to be fewer high-y peaks in the GRF. Some additional discrepancy may come from the imperfect peak finding in the GRF, namely, the selection of peaks at only one scale.

However, in the anisotropic components $m=2$ and $m=4$, significant distinctions appear. Generally, oriented stacking with the GRF results in some anisotropic signal because, as in the real universe, the field around any peak has a preferred alignment. The imposed correlation between the pseudo-galaxy and pseudo-y fields results in an aligned pseudo-y signal. We find that the $m=2$ moment is similar in shape and peak height for the full cluster sample between the GRF and Buzzard, although not perfectly in agreement. Notably, though, imposing the field constraints makes a significant difference. The Buzzard profile rises by ~2.5 × while the constrained GRF stack remains nearly the same as before the constraints were applied. This demonstrates that imposing $v>2$ and $e>0.3$ on the Buzzard galaxy field selects for cluster locations in the y-map where the local anisotropy is highly non-Gaussian. Finally, the $m=4$ moment demonstrates the most striking signal of non-Gaussianity. The profile of the GRF stack for both peak samples is consistent with zero, while the Buzzard profile shows a significant peak, which increases with the imposed field constraints. We interpret the $m=4$ moment as a sign of filament structure, flattened from small-scale gravitational effects in the late-time universe.
To support this claim, in the following section we examine
the *m* = 4 signal in oriented stacks from the Websky
simulation.

# 5. Expectations from Theory

## 5.1. Websky

We use the Websky simulations to give pure theoretical
results that are not subject to the selection effects of DES and
the Buzzard mocks. Two key questions are:

1. How does superclustering depend on redshift?

2. How does the average superclustering signal from hot gas depend on the mass of the clusters being stacked?

Websky is useful for addressing the former question because
the simulations extend to redshifts of *z* = 4.6. In addition,
answering either question requires splitting data into smaller
subsamples, which increases noise in the stacks. Due to the
full-sky coverage of this simulation, there are enough clusters
that the S/N of subsample stacks is sufficient to distinguish
between the different results.

We create the analog to a galaxy overdensity map by using
all Websky dark matter halos in the mass range [1.5 × 10¹²,
1 × 10¹⁵] M☉. This range incorporates most redMaGC-galaxy-
hosting halos (Clampitt et al. 2017; Pandey et al. 2021), and
also includes most halos that would host galaxies from the
Sloan Digital Sky Survey BOSS-CMASS sample (Dawson
et al. 2013; Maraston et al. 2013; Sonnenfeld et al. 2019). The
lower limit is slightly higher than the Websky minimum halo
mass (~1.2 × 10¹² M☉) because cutting out the lowest-mass
halos saves computational time. Nevertheless, the wide halo
range paints a near-complete picture of the LSS, while the
redMaGC galaxy data only form a subset of the galaxies that
trace the underlying matter. We create mass-weighted halo
number-density and overdensity maps using a weight of
M_h/10¹² M☉. These maps are roughly proportional to galaxy
overdensity maps because the number of galaxies in a halo is
approximately linearly proportional to the mass (Kravtsov et al.
2004). With smoothing scales far larger than a cluster, the
missing details of the subhalo distribution are unimportant.

We test how the extended gas around massive clusters varies
with redshift by stacking on halos with $M > 10^{14}$ M☉ from
$z \sim 0.25$ to $z \sim 2$. We rescale the stacked images for pairs of
consecutive 200 Mpc slices to the same physical size and
combine them. The results are shown in Figure 12, demonstrat-
ing that the anisotropic tSZ signal increases from a low level at
high redshifts (early times) until $z \sim 0.75$, after which it
becomes fairly stable. As time evolves past $z = 0.75$, the peak
exhibits a slight increase and subsequent decrease. However, as
most of these low-$z$ profiles are consistent within the $1\sigma$
regions, we leave a detailed study of the low redshifts to
future work.

These results can be interpreted as a combination of the halo
mass function evolution as well as the evolution of super-
clustering in the late-time universe. As the simulation evolves,
the cosmic web structure in which clusters are embedded
becomes more pronounced. Halos merge to form larger and
more massive halos, which are prescribed larger tSZ profiles in
post-processing. This prescription mimics the theoretical and
observational understanding that, over time, galaxies merge
along filaments and flow toward clusters. These mergers are
expected to heat up the gas in galaxy groups. In addition,

<figure><img src="image_15.png" /></figure>

**Figure 12.** The *m* = 2 (top) and *m* = 4 (bottom) radial profiles for oriented stacks of Websky y in progressive 400 Mpc slices in the line-of-sight distance. These stacks of massive clusters (*M* > 10¹⁴ M☉) were oriented using a mass-weighted halo density field smoothed at the scale of an FWHM = 14 Mpc. The profiles are binned equally in arcminutes, thus the bin sizes and centers in megaparsec vary across the different redshifts. We remove the first *m* = 4 bin to account for uncertainty in the angular integration near the center of the stacked image. In both moments, the peak location is fairly consistent with redshift and there is a trend toward lower overall signal at higher redshifts. The profiles at lower redshifts have smaller uncertainties because of the progressively larger number of massive halos at later times. In the range available in DES (0.25 < *z* < 0.7), the profiles are roughly consistent within the shaded 1σ error regions. This suggests that the evolution of superclustering on these scales slows down during these redshifts, due in part to the onset of Λ domination at *z* ~ 1.

mergers and AGN feedback can also heat up intergalactic gas. The effect of feedback and mergers are de facto included in the Websky simulations, encoded in the gas response to the presence of the evolving halo population, but there are many possibilities that are not accounted for. When galaxies fall into clusters, they undergo ram pressure stripping and the clusters gain mass and heat up through shocks. Overall, the tSZ signal from the contributions of clusters, groups, and intergalactic gas should grow over time and thus the anisotropic superclustering signal should be stronger at lower redshifts. Our results match this expectation.

We next explore how the strength of the gas signal from
superclustering varies with the mass of the stacked clusters. We
divide the data into three bins, the ranges for which are shown
in Table 2. The table also shows the equivalent DES cluster
<figure><img src="image_16.png" /><figcaption>**Figure 13.** Stacks of the Websky Compton-y-map on dark matter halo locations. From left to right, stacks are on the low-, mid-, and high-mass bins of Websky halos described in Table 2. Each stack combines halos in the distance range of 1232–1832 Mpc, or $z \sim 0.4$.</figcaption></figure>

<table><caption>**Table 2**
Websky Halo Mass Bins</caption><thead><tr><th>Bin Title</th><th>Mass Range (M<sub>200m</sub> (M<sub>c</sub>))</th><th>Equiv. λ</th><th>N</th></tr></thead><tbody><tr><td>Low mass</td><td>[1 × 10<sup>13</sup>, 5 × 10<sup>13</sup>]</td><td>[3, 11]</td><td>3,017,917</td></tr><tr><td>Mid mass</td><td>[5 × 10<sup>13</sup>, 10<sup>14</sup>]</td><td>[11, 18]</td><td>314,389</td></tr><tr><td>High mass</td><td>[10<sup>14</sup>, 2.6 × 10<sup>15</sup>]</td><td>[18, 200]</td><td>169,316</td></tr></tbody></table>

**Note.** The equivalent cluster/group richness in the third column is calculated from Equation (52) in McClintock et al. (2019) for a redshift of 0.5, with the caveat that the mass-richness relation should not be extrapolated to masses as small as the low-mass bin. The final column gives the number of halos in a 200 Mpc slice at $z \sim 0.4$.

richness range ($\lambda$) at $z=0.5$ assuming the mass-to-richness relation from McClintock et al. (2019). After stacking the cluster cutouts, we combine stacks from five slices from 1200–1800 Mpc to get the images shown in Figure 13. Figure 14 shows the $m=0, 2$, and 4 moments of the stacks in each mass bin. Larger halos have higher isotropic tSZ profiles, as the y signal is prescribed in the simulations to be proportional to $M^{5/3}$. Our results demonstrate that additionally, the anisotropic y signal is stronger for structures surrounding higher-mass halos.

This stronger signal could result from a combination of effects. Because the Websky simulations are hydrodynamical only in the limited-spatial-range response to the presence of halos, albeit of all masses, we emphasize that this dependence is *not* from temperature differences in intergalactic filament gas, since these are not included in these simulations. The contributing tSZ sources are clusters and groups along the alignment axis. The increase in the $m=2, 4$ and signal around more massive clusters is due to their being embedded in overall denser filaments on average. These regions would have more halos along the alignment axis contributing tSZ signal. In addition, massive clusters are known to have higher connectivity to the cosmic web (Aragón-Calvo et al. 2010; Codis et al. 2018b). In other words, more massive clusters are connected to a higher number of filaments. This could contribute to some of the correlation between mass and tSZ signal in $m=2$ and 4; massive clusters may have multiple filaments partially overlapping along the line of sight of the alignment axis, boosting the signal. Both effects could be contributing simultaneously.

<figure><img src="image_17.png" /><figcaption>**Figure 14.** $m=0, 2, 4$ radial profiles for the oriented stacks of three mass bins of Peak Patch halos shown in Figure 13. The highest-mass bin has the largest signal in all components, demonstrating that not only is the isotropic gas signal higher for larger clusters (as prescribed), but the surrounding supercluster structure also has a stronger signal. Notably, the $m=2$ and $m=4$ components rise above the $m=0$ component at their peaks, demonstrating that outside of the central stacked cluster, the anisotropic structure contributes a larger tSZ signal than the isotropic component.</figcaption></figure>

We attempt to address the same questions with the observational data. However, we find that splitting the small amount of available cluster data into richness bins increases the noise in the stacks enough that the results for different bins are all consistent within the errors. In addition, the redMaGC high-density galaxy catalog extends only to $z=0.7$. Within this range, the Websky results indicate that the superclustering signal is expected to remain fairly constant. We determine that the error bars in ACT × DES and even the Buzzard mocks are too large to observe a redshift evolution in this range, and any apparent evolution may be only a function of the data selection functions.

## 5.2. Dependence on Parameters

We examine the dependence of superclustering y signal on cluster richness $\lambda$, field excursion $\nu$, field ellipticity $e$, and smoothing scale using the Buzzard mocks. Figure 15 shows the
<figure><img src="image_18.png" /><figcaption>Figure 15. Results of changing the lower threshold for cluster richness (λ, left column), field excursion on a 14 Mpc scale (ν, middle column), and field elongation on a 14 Mpc scale (e, right column). For ν changes, the λ threshold is fixed at 10; for e changes, the ν threshold is fixed at 2. From top to bottom, the plots are of m = 0, 2, 4. 1σ error bars are taken from the diagonal of the covariance matrices calculated through cluster sample splits. Increases in λmin and νmin augment the signal in all three moments, while e min only has an impact on the anisotropic moments. The strongest effect on m = 2 comes from changes to νmin. We adopt (λ > 10, ν > 2, e > 0.3) as the combination of constraints to apply to observed ACT × DES data, as they increase the anisotropic signal without overly depleting the cluster sample.</figcaption></figure>

effects of imposing minimum limits on λ, ν, and e on the stacked Compton-y signal. The ν property has the strongest impact, demonstrating that clusters embedded in large-scale highly overdense regions tend to be members of more massive filaments and superclusters.

Certain cuts to the cluster sample significantly boost the oriented y signal, yet a reduction of factor N to the number of stacked images augments the random noise by $\sqrt{N}$. For the different parameter values tested, we assess this trade-off by examining the S/N of the maximum bin of m = 2. We find that the full $\lambda > 10$ sample has the highest S/N. If including the field constraints, a combination of ($\lambda > 10$, $\nu > 2$, $e > 0.3$) is optimal.

Next, we examine changes to the smoothing scale that is applied to the $\delta_g$ maps. As Figure 16 visually demonstrates, finer-grained smoothing causes the orientation for each cluster to be more dependent on local features such as nearby filament galaxies, whereas coarser smoothing makes larger-scale features (like the nearest-neighbor clusters) more important to the determination of curvature. With the Buzzard simulations, we perform oriented stacking with a Gaussian-smoothed map of an FWHM = [6, 10, 14, 18] Mpc. The results are shown in Figure 17. We only show m = 2 for brevity, but m = 4 displays similar effects.

We find that the location of the peak of the m = 2 profile scales nearly linearly with the smoothing scale. This is because the smoothing scale sets the radius from the cluster at which orientation is determined, and thus by construction, anisotropic structure is maximally aligned along the horizontal image axis of the stack at that radius. We will call this the *radius of maximal alignment*. Figure 18 demonstrates this concept by rescaling the x-axis by the Gaussian FWHM, bringing nearly all peaks into alignment. m = 4, not shown, behaves similarly. We emphasize that the peak location therefore contains little to no physical information. The *height* of the profile at the radius of maximal alignment, however, is physical, as it is determined by the average temperature and density of aligned structure at that scale. We will further analyze the meaning of peak height and its relationship to physical properties in a subsequent paper.

Another useful quantity is the total integrated Compton-y signal for each moment,

$$Y_m = \int_0^R C_m(r) r dr, \qquad (13)$$

where $C_m(r)$ is defined in Equation (7). For m = 0, this is similar to the unitless angularly integrated Compton-y parameter, with the difference that our Y is in units of square
<figure><img src="image_19.png" /><figcaption>**Figure 16.** Overdensity maps smoothed at three scales and the respective orientations determined for two clusters. For the left cluster, smaller-scale surrounding structure is fairly aligned with larger-scale structure so the orientation is only slightly rotated between the different maps. The right cluster lies at a kink in the surrounding LSS, so the orientation angle depends more strongly on the smoothing scale. All three maps are smoothed at scales larger than a cluster, so the small-scale peaks corresponding to the two clusters are not visible in the maps.</figcaption></figure>

megaparsecs. $Y_2$ is shown in the lower panel of Figure 17. $Y_2$ depends not only on the gas properties at the radius of maximal alignment, but also on how coherent the structure is within and beyond that radius. It also depends on the galaxy field constraints $v$ and $e$, which limit the cluster sample in a smoothing-dependent manner. As the smoothing becomes coarser, we observe that the $m=2$ $y$ profiles broaden and $Y_2$ increases. This may suggest that LSS is more coherently aligned along the direction of orientation determined for larger scales than for smaller scales. However, because of the nontrivial impacts of smoothing on the galaxy field constraints, we leave a more robust physical interpretation for the succeeding paper.

The flexibility of our method to smoothing demonstrates that it can be applied to scales as small as individual clusters, to study the alignment of galaxies and cluster gas, and as large as the longest superclusters.

**6. Comparison with Observational Data**

We apply the method to ACT $\times$ DES data and the Buzzard mocks using a few combinations of smoothing scales and field constraints. We will highlight results for 18 Mpc smoothing, a scale that roughly corresponds to intercluster filaments (as motivated in Section 4.2). We later show all three smoothing

<figure><img src="image_20.png" /><figcaption>**Figure 17.** Effect of changing the smoothing scale applied to $\delta_g$. Top: changes to the $m=2$ radial profile; bottom: changes to the integrated $m=2$ Compton-$y$ signal out to 40 Mpc. All profiles are constrained by ($\lambda > 10$, $v > 2$, $e > 0.3$), where the latter two parameters are calculated post-smoothing. The smoothing scale thus impacts both the selection and orientation of clusters. The location $r$ of the maximum highly depends on the FWHM and thus contains little physical information.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 18.** The same figure as Figure 17 (upper), but with $r$ scaled by the FWHM of the Gaussian to create a unitless radial quantity. This scaling brings nearly all curves into alignment, demonstrating that the location of the peak of $m=2$ should not be interpreted as the true position of the maximum anisotropic gas signal from the cluster, but rather a property that depends almost exclusively on the smoothing scale.</figcaption></figure>

scales. Each result figure shows the observational $y$ measurements in 3 radial bins, as motivated in Section 4.5. Each figure also shows the Buzzard results, which demonstrate what a pure measurement would look like across the full DES footprint, with error bars due to variance in the LSS but not instrumental noise. The Buzzard results are binned into the same 3 radial bins to assess the consistency with observations, while the
<figure><img src="image_22.png" /><figcaption>Figure 19. Comparison of the binned $m = 0$ component of a Buzzard oriented stack (triangles) vs. ACT × DES (circles), with continuous Buzzard profiles also shown in the background for reference. Red indicates the stack using the full $\lambda > 10$ cluster sample, and blue indicates the sample constrained by $\nu$ and $e$ thresholds for 18 Mpc smoothing. The data are binned evenly, including only data beyond the cluster radius (black vertical dashed line). Each profile has been adjusted by subtracting its tail, the average value from 33–40 Mpc. The binned profiles represent the average tSZ signal from LSS beyond the central stacked cluster. $\nu$ thresholding boosts this signal within the first ~20 Mpc.</figcaption></figure>

figures also display the continuous Buzzard profiles for visual purposes.

To begin, we briefly examine the isotropic ($m = 0$) component of the stacked images in Figure 19. This is identical for an oriented and unoriented stack and will not provide specific information about filaments. Nevertheless, comparing the $m = 0$ Compton-y signal beyond the stacked cluster radius between simulations and observations is informative of how well the simulations reproduce the large-scale clustering and gas content of halos.

Because our focus lies beyond the cluster interior, we exclude the majority of central-cluster tSZ contributions to the radial y profile by choosing an inner cutoff radius. Our choice is $R_c = 1.5R_\lambda$, where $R_\lambda$ is the redMaPPer cluster radius (Rykoff et al. 2014). $R_\lambda$ is given by $(\lambda/100)^{0.2}h^{-1}$ physical Mpc, which we calculate at the median richness in our sample, $\lambda = 15$. After conversion, this results in a comoving radius of ~2.5 Mpc for the average redshift of the cluster sample. Therefore, in both the simulations and real data, we begin binning the signal beyond 2.5 Mpc.

The raw ACT × DES profiles each have a constant positive offset with respect to the Buzzard profiles (not shown in any figure). This is primarily due to long-wavelength noise in the y-map from residual low-$\ell$ primary CMB contamination, which does not average down with more stacked clusters. Stacks on random points in simulations of the ACT y-map described in Section 4.5 each result in a different constant offset depending on different realizations of the low-$\ell$ primary CMB. To account for the offset in the real ACT y-map and each simulation, we subtract the average value of $C_0(r)$ from 33–40 Mpc (the tail of each profile) from the full profile. This is similar to performing aperture photometry, a technique typically applied to unoriented stacks to subtract the noise calculated in an annulus from the signal within some inner radius.

The adjusted $m = 0$ profiles for 18 Mpc smoothing, after subtraction of the tail, are shown in Figure 19. There is a nonzero signal in the first two bins, which is boosted when the $\nu > 2$ constraint is applied. (The $e$ constraint has little to no effect on the isotropic profile, as previously demonstrated in Figure 15.) The $\nu$ and $e$ constraints reduce the cluster sample by a factor of ~5, causing an increase in noise of ~$\sqrt{5}$. We can determine the S/N of the observational results by comparing them to a null profile. This null result corresponds to the unphysical hypothesis that there is no tSZ signal outside of the average cluster radius. This would occur in a universe where every cluster was separated by > 40 Mpc from the nearest object containing hot gas. We compute the reduced $\chi^2$ of the data vector $y_{\text{obs}}$ (the coarsely binned $C_0(r)$ profile) with respect to null:

$$ \chi_{\text{red}}^2 = y_{\text{obs}}^{\text{T}} \Sigma^{-1} y_{\text{obs}} / N_{\text{bins}} \quad (14) $$

where $\Sigma^{-1}$ is the inverse covariance matrix of the data. Next, we find the probability that a truly null vector could exceed that $\chi^2$. This probability to exceed (PTE) is measured by drawing 10⁶ random vectors from a normal distribution with mean zero and the covariance matrix from the data. We compute the $\chi^2$ of each vector with respect to null, then find the fraction of the sample for which the $\chi^2$ exceeds that of the real data vector. This fraction is the PTE; lower values indicate that the random vector is unlikely to exceed the true data vector, providing stronger evidence for the detection of a nonzero signal in the data. We then relate this to an S/N estimate, which is the number of Gaussian sigmas away from the null,

$$ (1 - \text{PTE}) = \text{erf}(S/N / \sqrt{2}), \quad (15) $$

where erf is the error function. We also assess the goodness of fit to the Buzzard results by finding the reduced $\chi^2$:

$$ \chi_{\text{red}}^2 = (y_{\text{obs}} - y_{\text{sim}})^{\text{T}} \Sigma_{o+s}^{-1} (y_{\text{obs}} - y_{\text{sim}}) / N_{\text{bins}} \quad (16) $$

where $\Sigma_{o+s}$ is the summed covariance matrix for observational data and simulations, and $y_{\text{sim}}$ refers to the coarsely binned Buzzard profile.

The bins are highly correlated for $m = 0$, so the S/N is only 1.7 and 2.3 for the profiles without and with constraints, respectively. Both profiles agree well with Buzzard, with $\chi_{\text{red}}^2 < 1$. The low S/N of this extended isotropic signal, and the requirement to subtract a constant offset, are both arguments that extended structure is not probed well by unoriented stacking. Notably, the same constant offset does not appear in $m = 2$ or 4. In addition, we achieve a higher-S/N detection with $m = 2$.

The $m = 2$ moment of each stacked image, for the same 18 Mpc smoothing scale, is presented in Figure 20. We again remove the inner cluster region when coarsely binning the data to avoid cluster mis-centering effects. The binning is slightly uneven, designed to place the first bin at the location of the Buzzard continuous-profile peak. The signal for all $\lambda > 10$ clusters in Buzzard peaks at $y = 4.5 \times 10^{-8}$. Enforcing ($\nu > 2, e > 0.3$) raises the signal at the peak by a factor of ~2.3. The Buzzard results are in strong agreement with ACT × DES. We also show a null test, the $m = 2$ of an unoriented stack in which each cluster cutout was randomly rotated before stacking. This profile is consistent with zero as expected. Additionally, the figure shows the covariance matrices for both ACT × DES data vectors.
<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /><figcaption>**Figure 20.** Top: comparison of the *m* = 2 component of the Buzzard (triangles) vs. ACT × DES (circles) oriented Compton-y-map stacks. The orientation and *ν* and *e* selection were performed using a Gaussian-smoothed galaxy field with an FWHM=18 Mpc. The background shaded curves surround 1σ regions above and below the continuous mean profile from Buzzard; the triangles correspond to the coarse binning of this profile for direct comparison with real data. The first triangle falls below the curve because of the wide bin size which includes lower values on either side of the peak. Points are artificially spaced in *r* for visual distinction. The attached lower panel shows a null test for comparison: the quadrupole profile for an ACT × DES unoriented stack using the same selected superclustering points. Bottom: the covariance matrices for each ACT data vector. The discrepancy between the visual appearance of the significance of the real data and the values reported in Table 3 is due to the correlations between bins.</figcaption></figure>

We compare the *m* = 2 signal to a completely null profile—corresponding to a perfectly isotropic stack—to determine its significance. Having no signal in *m* = 2 would indicate that either (a) the average thermal energy distribution surrounding clusters is highly isotropic, (b) the galaxy distribution is uncorrelated with the gas distribution, or (c) our oriented stacking method does not effectively align LSS. Table 3 presents the reduced-χ², PTE, and S/N values. Without constraints on the galaxy field, the ACT × DES oriented stacks have 3.5σ level evidence for signal in the *m* = 2 moment. Introducing the field constraints does not change the significance, as the larger error bars compensate for the boost in

<table><caption>**Table 3**
Summary of the ACT × DES Results</caption><thead><tr><th rowspan="2">m</th><th rowspan="2">Cuts (Y/N)</th><th rowspan="2">N</th><th rowspan="2">PTE</th><th rowspan="2">S/N</th><th colspan="2">χ<sup>2</sup><sub>red</sub> (dof = 3)</th></tr><tr><th>Absolute</th><th>Data - Sim</th></tr></thead><tbody><tr><td colspan="7" style="text-align: center;">FWHM of smoothing = 10 Mpc</td></tr><tr><td>2</td><td>Y</td><td>1190</td><td>0.06</td><td>2.6</td><td>1.8</td><td>0.9</td></tr><tr><td>4</td><td>Y</td><td>1190</td><td>0.3</td><td>2.4</td><td>3.9</td><td>1.2</td></tr><tr><td colspan="7" style="text-align: center;">FWHM of smoothing = 14 Mpc</td></tr><tr><td>2</td><td>Y</td><td>1103</td><td>0.004</td><td>3.2</td><td>5.2</td><td>0.5</td></tr><tr><td>4</td><td>Y</td><td>1103</td><td>0.2</td><td>1.5</td><td>1.8</td><td>1.1</td></tr><tr><td colspan="7" style="text-align: center;">FWHM of smoothing = 18 Mpc</td></tr><tr><td>2</td><td>N</td><td>5494</td><td>0.0004</td><td>3.5</td><td>6.0</td><td>0.1</td></tr><tr><td>4</td><td>N</td><td>5494</td><td>0.4</td><td>0.9</td><td>1.0</td><td>0.2</td></tr><tr><td>2</td><td>Y</td><td>975</td><td>0.0005</td><td>3.5</td><td>5.9</td><td>0.07</td></tr><tr><td>4</td><td>Y</td><td>975</td><td>0.1</td><td>1.5</td><td>1.8</td><td>0.3</td></tr></tbody></table>

**Note.** From left to right, *m* is the multipole moment, “Cuts” refers to whether or not the cluster sample has been constrained by *ν* > 2, *e* > 0.3; *N* is the number of stacked clusters, PTE is the probability to exceed, S/N gives the number of Gaussian *σ* from null, and χ²*red* is the reduced χ² value using 3 dof. χ²*red* is shown for the absolute value with respect to 0 as well as the value with respect to Buzzard. The highest S/N detections of extended signal come from the *m* = 2 components of the 14 and 18 Mpc smoothed stacks. Data and simulation are generally in agreement, with χ²*red* (data-sim) consistently near or below 1.

<table><caption>**Table 4**
Summary of the Buzzard Results</caption><thead><tr><th>m</th><th>Cuts (Y/N)</th><th>N</th><th>S/N</th><th>χ<sup>2</sup><sub>red</sub> (dof = 3)</th></tr></thead><tbody><tr><td colspan="5" style="text-align: center;">FWHM of smoothing = 10 Mpc</td></tr><tr><td>2</td><td>Y</td><td>6165</td><td>10</td><td>31</td></tr><tr><td>4</td><td>Y</td><td>6165</td><td>4.8</td><td>7.8</td></tr><tr><td colspan="5" style="text-align: center;">FWHM of smoothing = 14 Mpc</td></tr><tr><td>2</td><td>Y</td><td>5401</td><td>11</td><td>37</td></tr><tr><td>4</td><td>Y</td><td>5401</td><td>5.3</td><td>9.2</td></tr><tr><td colspan="5" style="text-align: center;">FWHM of smoothing = 18 Mpc</td></tr><tr><td>2</td><td>N</td><td>24,922</td><td>12</td><td>51</td></tr><tr><td>4</td><td>N</td><td>24,922</td><td>3.9</td><td>5.1</td></tr><tr><td>2</td><td>Y</td><td>4,836</td><td>9.5</td><td>30</td></tr><tr><td>4</td><td>Y</td><td>4,836</td><td>2.8</td><td>2.6</td></tr></tbody></table>

**Note.** The table repeats the calculations in Table 3 for the Buzzard simulations, with the same 3 radial bins. χ²*red* is only shown for the value with respect to null, as the (data-sim) value is already presented in Table 3. Because the simulation results are high S/N, rather than calculating a PTE we simply estimate the S/N by taking the square root of χ². Generally, the detection significance is much higher for Buzzard than for ACT × DES because the simulation is noiseless. *m* = 2 is detected at a much higher significance than *m* = 4.

signal. Table 4 presents similar summary statistics for the Buzzard results, which have much higher significance (12σ and 10σ) because of the lack of instrumental noise and larger cluster sample.

Finally, the *m* = 4 results are shown in Figure 21. The same bins are applied to *m* = 4 as *m* = 2. The Buzzard profiles peak at ~4 times smaller *y* values than *m* = 2. For this component, a null profile could correspond to any of the (a)–(c) possibilities listed for *m* = 2 or the signal from a perfectly Gaussian random field. As shown in Tables 4 and 3, the Buzzard *m* = 4 component from the constrained cluster sample is at 2.8σ, but the corresponding ACT × DES result is only 1.5σ. Therefore, we cannot claim evidence for the *m* = 4 component in
<figure><img src="image_25.png" /><figcaption>**Figure 21.** Top: comparison of the *m* = 4 component of Buzzard vs. ACT × DES for smoothing with an FWHM = 18 Mpc. The attached lower panel shows the *m* = 4 profile for an unoriented stack as a null test reference. The second red point and third blue point are significantly (1.3σ and 1.1σ, respectively) higher in the real data than in the simulation. Because the tension is not extreme, this may be due to random fluctuations in the ACT-y-map at these cluster locations; larger data sets will address whether there is a true physical difference in the average profile. Overall, however, the Buzzard and ACT × DES profiles are in statistical agreement. Bottom: the covariance matrices for each ACT data vector.</figcaption></figure>

observations. The *m* = 4 profiles are in agreement with Buzzard.

We repeat the stacking procedure for two smaller smoothing scales, FWHM = 10 Mpc and FWHM = 14 Mpc, with the *v* and *e* constraints enforced for the respective smoothed maps. The *m* = 2 and 4 plots are shown in Figure 22. The binning is adjusted for each scale to align the first bin with the peak location. The continuous profiles in the figures are shown for visual purposes but cannot be directly compared to the binned real data; instead, this comparison is best captured in the $\chi^2_{\text{red}}$ column of Table 3.

There is evidence of the *m* = 2 signal at the 3.2σ level for 14 Mpc smoothing and marginal 2.6σ evidence for 10 Mpc. The S/N of *m* = 4 is smaller in all cases, although there is marginal 2.4σ for an *m* = 4 component in the smallest smoothing scale. This is consistent with the expectation that

<figure><img src="image_26.png" /><figcaption>**Figure 22.** Comparison of the *m* = 2 component (top) and *m* = 4 component (bottom) of Buzzard vs. ACT × DES for three smoothing scales: FWHM = [10, 14, 18] Mpc. The Buzzard finely binned data are shown for visual purposes, but are not directly comparable to the ACT × DES points because of the difference in binning (see Table 3 for a statistical comparison with the same binning). Points are artificially spaced in *r* for visual purposes. The cluster sample varies slightly for each scale due to the scale-dependent *v* > 2, *e* > 0.3 galaxy field constraints. The differences in profiles are mostly due to the dependence of orientation on smoothing, as discussed in Section 5.2. Measurements for the three scales are highly correlated with each other due to the alignment of cosmic web structure across a wide range of scales, so they should not be interpreted as independent. All data points at *r* ~ 33 Mpc are significantly higher than the respective binned Buzzard point, but as this tension is at the 1.5σ level at most, more data are necessary to determine whether this is a random fluctuation or true difference.</figcaption></figure>

non-Gaussianity is more pronounced on smaller scales. However, because the evidence does not meet the 3σ threshold, we leave further analysis to future work with larger data sets. Generally, the observed data are in very good agreement with the simulations, reaching at most a reduced $\chi^2$ of 1.2.

Due to the lack of detection in *m* = 4, we cannot claim to find evidence for non-Gaussian structure with the currently available data. For *m* = 2, as discussed in Section 4.6, the boost from applying the *v* and *e* thresholds only occurs in realistic fields and not in a Gaussian random field. Therefore, it is likely
<figure><img src="image_27.png" /><figcaption>**Figure 23.** Comparison of the integrated *m* = 2 profile to *R* = 40 Mpc between Buzzard and ACT × DES for three smoothing scales. At all scales, the simulated and observed results are in agreement within 1σ.</figcaption></figure>

that the strength of *m* = 2 in stacks on the constrained cluster sample is indicative of non-Gaussianity. However, to robustly address this, we would need to repeat the study in Section 4.6 with the exact methods applied to the final results (using thin redshift slicing, rescaling, and combining multiple slices). We leave this detailed comparison for future work.

As discussed in Section 5.2, the shape of each profile is highly dependent on smoothing scale and thus the integrated profile provides a useful single-value quantity for comparison between data and simulations. Figure 23 shows the integrated *Y* signal (computed with a Riemann sum) over the Buzzard (red) and ACT × DES (blue) coarsely binned *m* = 2 profiles for three smoothing scales. Errors are propagated with the covariance matrix for each scale. At all scales, the Buzzard and ACT integrals are within 1σ of each other. The ACT integrals have an S/N = [1.0, 2.0, 2.7] for the [10, 14, 18] Mpc smoothing scales. Therefore, we find marginal indications of integrated anisotropic clustering of thermal energy for the two coarser scales.

# 7. Conclusions

## 7.1. Overview

This paper introduced a new real-space method for probing anisotropies in gas signal in the cosmic web. By combining millimeter-wavelength data from the ACT CMB survey with optical data from the DES galaxy survey, we measured the average superclustering signal from hot gas surrounding DES redMaPPer clusters. We showed that there is a significant Compton-y signal from non-Gaussianity in the late-time universe visible in the *m* = 4 moment of the simulated stacks. There is a marginal indication of the *m* = 4 signal in the observational data. Using characteristics of the galaxy field, we identified highly overdense, elongated regions. Selecting for clusters in these high-superclustering areas enhances the gas signal in the *m* = 2 moment of oriented stacks, causing a distinction from a GRF. We visually demonstrated this enhancement in the observational data, leaving a more rigorous proof of non-Gaussianity to future work. Generally, with observed ACT × DES data, we found marginal-to-significant evidence for an extended signal in the *m* = 2 moments; the

significance depended on the chosen galaxy field smoothing. This evidence demonstrates that the average thermal energy distribution around clusters is not isotropic. This, of course, is expected because clusters are embedded in filamentary structures in the cosmic web spanning a wide range of scales. Comparing to the Buzzard mocks showed broad agreement in the *m* = 2 and 4 radial profiles. For the integrated Compton y from *m* = 2, the results at all scales are in agreement.

We did not attempt to constrain the contribution of the WHIM, due to gas outside individual halos, to the anisotropic y signal. In a previous filament-stacking study, de Graaff et al. (2019) found a filament signal with $y = 6 \times 10^{-9}$. The team also found that bound gas in halos only contributes ~20% of this signal. A WHIM signal at $y \sim 5 \times 10^{-9}$, then, is ~5% of the peak $m = 2$ signal measured in our work. If the de Graaff et al. (2019) study was correct, this indicates that as expected, the dominant tSZ signal from the extended LSS surrounding clusters comes from bound gas in halos. The previously measured WHIM signal is smaller than any of the 1σ error bars in our work and thus undetectable with the current method and data. A detection may be possible with the currently available data by stacking on all galaxies instead of clusters, as well as by applying a more sophisticated approach to the combination of smoothing scales and field constraints. A WHIM detection may also be possible with the future, expanded version of the ACT y-map and greater overlap with DES cluster data. We leave a detailed study of how oriented stacking can be used to characterize the WHIM for future work.

With the Websky full-sky simulations, we demonstrated the potential of the oriented stacking method to probe the evolution of superclustering out to $z \sim 2$ and the relationship between cluster mass and the surrounding environment. These theory results demonstrated that the anisotropic gas signal from superclustering is expected to grow with time from $z \sim 2$ to ~0.7, then stabilize. For fixed redshift, the strength of superclustering correlates with cluster mass, indicating that massive clusters are embedded in more massive and/or denser filaments. The limited number of available clusters in the ACT y-map overlapping with DES prohibited the same analyses on real data.

## 7.2. Systematics and Future Outlook

We identify a few sources of systematic uncertainty, which were negligible in this study but will become important as the ACT data expands and improves in the next few years. First, as noted in Section 2, $\lambda < 20$ redMaPPer clusters are known to contain false detections. False clusters should be less correlated with the surrounding LSS than real clusters, and therefore their inclusion should bias our results low. Nevertheless, we found evidence for superclustering in ACT × DES, indicating that the impact of cluster sample impurities was not enough to drown out the signal. Estimating how much higher the signal would be with a completely clean $\lambda > 10$ sample is beyond the scope of this work, as these low-richness impurities have not yet been well characterized.

A possible source of systematic error in the comparison with simulations, as noted in Section 2, is the failure of Buzzard to reproduce the redMaPPer mass-richness relation. As the errors on ACT × DES decrease with future data, this may become a source of tension between the simulation and data. However, there are reasons to suspect this has a minor effect on our measurements compared to other factors such as the imperfect
gas prescription applied to Buzzard. First, because this work
does not focus on the cluster interior, discrepancies in the
richness distribution are only important to the extent that they
correlate with the surrounding LSS. We found in Section 5.2
that cluster richness is indeed correlated with the anisotropic
Compton-y signal; however, this is subdominant compared to
the dependence on v evaluated at larger scales. Additionally,
because the mock Compton-y-map is created using only halo
information from the Buzzard dark-matter-only simulation, the
y-map signal does not depend on the realistic identification of
redMaPPer clusters. The most important steps in our pipeline
are the selection and orientation of clusters using the redMaGiC
mock galaxy data, where Buzzard does well in matching DES.
It may be possible to entirely avoid the mass-richness
systematic in future work by using the Buzzard state of the
art: a recent version of Buzzard improves the color-dependent
clustering that heavily affects redMaPPer selection (DeRose
et al. 2022). In general, to fully understand the factors
influencing the oriented tSZ signal, we must disentangle the
effects of cosmology, galaxy and cluster modeling/selection,
and gas prescription. A succeeding paper will address this by
applying oriented stacking to the galaxy distribution alone and
comparing the DES results to Buzzard.

Moving forward, future Compton-y-maps from ACT will
cover a sky area with complete overlap with DES over the 5000
deg² footprint (Naess et al. 2020). They will additionally
include high-resolution 30, 40, and 230 GHz data collected
using Advanced ACT arrays, allowing for better removal of
contamination from Galactic synchrotron sources, extragalactic
radio sources, the cosmic infrared background, and dust
(Madhavacheril et al. 2020). Our study will be repeated with
the full cluster sample from DES, which is approximately 14
times larger than the sample used in this study. The systematic
error contribution is at the few percent level for the current y-
map, so random errors will still dominate. Therefore, the error
bars on the data are expected to shrink by a factor of
√14 = 3.7. If the binned y values in the m = 2 profile for
18 Mpc smoothing remained the same with 3.7 × smaller error
bars, the detection of the m = 2 signal after v and e cuts will be
at the 13σ level. The m = 4 component would be detectable at
the 5σ level for the same scale. Slight tension would emerge
with Buzzard, with a χ² red of 2.5. With the next iteration of
ACT data, oriented stacking will be able to address questions of
non-Gaussianity and assess the accuracy of the tSZ pasting
prescriptions in simulations in more detail.

Additionally, the larger cluster sample will allow us to
remove low-λ clusters and repeat the study with a more pure,
higher-richness sample. This is expected to boost the super-
clustering gas signal. In a succeeding paper, we will use this
expanded data to study the impact of the WHIM on the
superclustering signal from hot gas.

Looking further to the future, The Simons Observatory (SO;
Ade et al. 2019) is currently being built in Chile, and is
expected to begin taking CMB observations in the mid-2020s.
With measurements in six frequencies, the SO Compton-y-map
will attain lower noise levels and higher resolution than that of
ACT. In addition, upcoming telescopes like the Fred Young
Sub-millimeter Telescope (formerly CCAT-p; Terry et al.
2019) and the future CMB S4 mission (Carlstrom et al. 2019)
will further advance tSZ measurements by mapping the
millimeter sky with improved sensitivity and frequency
coverage.

To optimize the science returns from oriented stacking, these
higher-quality CMB maps will need to be coupled with
augmented cluster and galaxy data. The SO team predicts that
the number of clusters detected through the tSZ effect will rise
by an order of magnitude from current levels (Ade et al. 2019).
The sample will extend to further redshifts than the DES red-
sequence-detected sample because the tSZ effect is redshift
independent, whereas received flux from optical galaxy
observations (as in DES) diminishes with redshift. SO forecasts
~200 detected clusters in a bin of width Δz = 0.1 at z ~ 1.5,
whereas the DES redMaPPer catalog extends only to z ~ 0.9.
Higher-redshift clusters will enable us to measure the evolution
of superclustering further into the universe's past. In addition,
the Dark Energy Spectroscopic Instrument (Dey et al. 2019)
will provide many new spectroscopic redshifts for galaxies and
clusters, enabling improved tomographic slicing to study
redshift evolution of the LSS (Zou et al. 2021). The Vera C.
Rubin Observatory Legacy Survey of Space and Time (Ivezić
et al. 2019) will make unprecedented contributions to the
optical galaxy and cluster data with deep observations of the
southern sky resulting in an estimated 20 billion galaxies.
Looking to other wavelengths, eROSITA is expected to deliver
hundreds of superclusters detected in the X-ray: additional
objects on which to stack the tSZ signal (Ghirardini et al.
2021). Applying the methods detailed in this paper to
combinations of these varied data sets, where they overlap,
will greatly increase the S/N and redshift extent of super-
clustering measurements.

This paper demonstrates a flexible tool for examining
superclustering in the universe. We have shown its effective-
ness for a combination of galaxy and gas data, but it can also be
used to combine galaxy and weak lensing data or to probe the
anisotropies of the galaxy field alone. By combining different
probes, this method can be implemented to study the
anisotropic bias of galaxies and gas with respect to dark
matter. A follow-up paper will use oriented stacking on galaxy
number density maps and compare the results to extended gas
signal to probe information about the anisotropic bias, the
baryonic content of halos, and possibly the baryons outside of
halos. This may provide more insight into the census of
baryonic content in filaments. In addition, future work will
determine whether this method of measuring anisotropic
clustering can complement two-point clustering statistics in
constraining viable models of dark energy.

Canadian co-authors acknowledge support from the Natural
Sciences and Engineering Research Council of Canada.
Computations were performed in part on the Niagara super-
computer at the SciNet HPC Consortium. SciNet is funded by
the Canada Foundation for Innovation; the Government of
Ontario; Ontario Research Fund—Research Excellence; and
the University of Toronto.

The authors thank Bruce Partridge for useful and extensive
comments on the draft.

J.R.B. was funded by the Natural Sciences and Engineering
Research Council of Canada Discovery Grant Program and a
fellowship from the Canadian Institute for Advanced Research
(CIFAR) Gravity and Extreme Universe program.

A.D.H. acknowledges support from the Sutton Family Chair
in Science, Christianity and Cultures.

R.H. is a CIFAR Azrieli Global Scholar, Gravity and the
Extreme Universe Program, 2019, and a 2020 Alfred P. Sloan
Research Fellow. R.H. is supported by the Natural Sciences and Engineering Research Council of Canada Discovery Grant Program and the Connaught Fund. The Dunlap Institute is funded through an endowment established by the David Dunlap family and the University of Toronto.

J.P.H. acknowledges funding for SZ cluster studies from NSF AAG No. AST-1615657.

M.L. acknowledges the support of the National Sciences and Engineering Research Council of Canada (NSERC) [PGSD - 559296 - 2021] and the Queen Elizabeth II / Graduate Scholarships in Science and Technology (QEII-GSST).

K.M. acknowledges support from the National Research Foundation of South Africa.

This work was supported by 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 operates in the Parque Astronómico Atacama in northern Chile under the auspices of the Comisión Nacional de Investigación (CONICYT). 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.

Funding for the DES Projects has been provided by the U.S. Department of Energy, the U.S. National Science Foundation, the Ministry of Science and Education of Spain, the Science and Technology Facilities Council of the United Kingdom, the Higher Education Funding Council for England, the National Center for Supercomputing Applications at the University of Illinois at Urbana-Champaign, the Kavli Institute of Cosmological Physics at the University of Chicago, the Center for Cosmology and Astro-Particle Physics at the Ohio State University, the Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M University, Financiadora de Estudos e Projetos, Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro, Conselho Nacional de Desenvolvimento Científico e Tecnológico and the Ministério da Ciência, Tecnologia e Inovação, the Deutsche Forschungsgemeinschaft and the Collaborating Institutions in the Dark Energy Survey.

The Collaborating Institutions are Argonne National Laboratory, the University of California at Santa Cruz, the University of Cambridge, Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas-Madrid, the University of Chicago, University College London, the DES-Brazil Consortium, the University of Edinburgh, the Eidgenössische Technische Hochschule (ETH) Zürich, Fermi National Accelerator Laboratory, the University of Illinois at Urbana-Champaign, the Institut de Ciències de l'Espai (IEEC/CSIC), the Institut de Física d'Altes Energies, Lawrence Berkeley National Laboratory, the Ludwig-Maximilians Universität München and the associated Excellence Cluster Universe, the University of Michigan, NFS' NOIRLab, the University of Nottingham, The Ohio State University, the University of Pennsylvania, the University of Portsmouth, SLAC National Accelerator Laboratory, Stanford University, the University of Sussex, Texas A&M University, and the OzDES Membership Consortium.

Based in part on observations at Cerro Tololo Inter-American Observatory at NSF's NOIRLab (NOIRLab Prop. ID 2012B-0001; PI: J. Frieman), which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation.

The DES data management system is supported by the National Science Foundation under grant Nos. AST-1138766 and AST-1536171. The DES participants from Spanish institutions are partially supported by MICINN under grant Nos. ESP2017-89838, PGC2018-094773, PGC2018-102021, SEV-2016-0588, SEV-2016-0597, and MDM-2015-0509, some of which include ERDF funds from the European Union. IFAE is partially funded by the CERCA program of the Generalitat de Catalunya. Research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Program (FP7/2007-2013), including ERC grant agreements 240672, 291329, and 306478. We acknowledge support from the Brazilian Instituto Nacional de Ciência e Tecnologia (INCT) do e-Universo (CNPq grant 465376/2014-2). We thank the anonymous referee for providing valuable comments which improved the quality of the manuscript.

This manuscript has been authored by 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.

This work received support from the U.S. Department of Energy under contract No. DE-AC02-76SF00515 at SLAC National Accelerator Laboratory. This research used computing resources at SLAC National Accelerator Laboratory and at the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility located at Lawrence Berkeley National Laboratory, operated under contract No. DE-AC02-05CH11231.

**Software:** Astropy (Astropy Collaboration et al. 2013, 2018), COOP (Huang 2016), healpy (Zonca et al. 2019), HEALPix (Górski et al. 2005), NumPy (Harris et al. 2020).

## ORCID iDs

M. Lokken https://orcid.org/0000-0001-5917-955X

R. Hložek https://orcid.org/0000-0002-0965-7864

A. van Engelen https://orcid.org/0000-0002-3495-158X

M. Madhavacheril https://orcid.org/0000-0001-6740-5350

E. Baxter https://orcid.org/0000-0002-6836-3196

J. DeRose https://orcid.org/0000-0002-0728-0960

C. Doux https://orcid.org/0000-0003-4480-0096

S. Pandey https://orcid.org/0000-0001-5780-637X

E. S. Rykoff https://orcid.org/0000-0001-9376-3135

G. Stein https://orcid.org/0000-0002-5193-516X

C. To https://orcid.org/0000-0001-7836-2261

S. Adhikari https://orcid.org/0000-0002-0298-4432

M. Aguena https://orcid.org/0000-0001-5679-6747

S. Allam https://orcid.org/0000-0002-7069-7857

J. Annis https://orcid.org/0000-0002-0609-3987

N. Battaglia https://orcid.org/0000-0001-5846-0411

G. M. Bernstein https://orcid.org/0000-0002-8613-8259

E. Bertin https://orcid.org/0000-0002-3602-3664

J. R. Bond https://orcid.org/0000-0003-2358-9949

D. Brooks https://orcid.org/0000-0002-8458-5047

A. Carnero Rosell https://orcid.org/0000-0003-3044-5150

M. Carrasco Kind https://orcid.org/0000-0002-4802-3194
J. Carretero https://orcid.org/0000-0002-3130-0204

R. Cawthon https://orcid.org/0000-0003-2965-6786

A. Choi https://orcid.org/0000-0002-5636-233X

M. Costanzi https://orcid.org/0000-0001-8158-1449

M. Crocce https://orcid.org/0000-0002-9745-6228

L. N. da Costa https://orcid.org/0000-0002-7731-277X

J. De Vicente https://orcid.org/0000-0001-8318-6813

S. Desai https://orcid.org/0000-0002-0466-3288

J. P. Dietrich https://orcid.org/0000-0002-8134-9591

J. Dunkley https://orcid.org/0000-0002-7450-2586

A. E. Evrard https://orcid.org/0000-0002-4876-956X

S. Ferraro https://orcid.org/0000-0003-4992-7854

B. Flaugher https://orcid.org/0000-0002-2367-5049

P. Fosalba https://orcid.org/0000-0002-1510-5214

J. Frieman https://orcid.org/0000-0003-4079-3263

P. A. Gallardo https://orcid.org/0000-0001-9731-3617

J. García-Bellido https://orcid.org/0000-0002-9370-8360

E. Gaztanaga https://orcid.org/0000-0001-9632-0815

D. W. Gerdes https://orcid.org/0000-0001-6942-2736

T. Giannantonio https://orcid.org/0000-0002-9865-0436

D. Gruen https://orcid.org/0000-0003-3270-7644

R. A. Gruendl https://orcid.org/0000-0002-4588-6517

J. Gschwend https://orcid.org/0000-0003-3023-8362

G. Gutierrez https://orcid.org/0000-0003-0825-0517

J. C. Hill https://orcid.org/0000-0002-9539-0835

M. Hilton https://orcid.org/0000-0002-8490-8117

A. D. Hincks https://orcid.org/0000-0003-1690-6678

S. R. Hinton https://orcid.org/0000-0003-2071-9349

D. L. Hollowood https://orcid.org/0000-0002-9369-4157

K. Honscheid https://orcid.org/0000-0002-6550-2023

B. Hoyle https://orcid.org/0000-0002-2571-1357

Z. Huang https://orcid.org/0000-0002-1506-1063

J. P. Hughes https://orcid.org/0000-0002-8816-6800

D. Huterer https://orcid.org/0000-0001-6558-0112

D. J. James https://orcid.org/0000-0001-5160-4486

T. Jeltema https://orcid.org/0000-0001-6089-0365

K. Kuehn https://orcid.org/0000-0003-0120-0808

M. A. G. Maia https://orcid.org/0000-0001-9856-9307

J. L. Marshall https://orcid.org/0000-0003-0710-9474

P. Melchior https://orcid.org/0000-0002-8873-5065

F. Menanteau https://orcid.org/0000-0002-1372-2534

R. Miquel https://orcid.org/0000-0002-6610-4836

K. Moodley https://orcid.org/0000-0001-6606-7142

R. Morgan https://orcid.org/0000-0002-7016-5471

F. Nati https://orcid.org/0000-0002-8307-5088

L. Page https://orcid.org/0000-0002-9828-3525

R. L. C. Ogando https://orcid.org/0000-0003-2120-1154

A. Palmese https://orcid.org/0000-0002-6011-0530

F. Paz-Chinchón https://orcid.org/0000-0003-1339-2683

A. A. Plazas Malagón https://orcid.org/0000-0002-2598-0514

A. Pieres https://orcid.org/0000-0001-9186-6042

A. K. Romer https://orcid.org/0000-0002-9328-879X

E. Rozo https://orcid.org/0000-0002-1666-6275

E. Sanchez https://orcid.org/0000-0002-9646-8198

M. Schubnell https://orcid.org/0000-0001-9504-2059

S. Serrano https://orcid.org/0000-0002-0211-2861

I. Sevilla-Noarbe https://orcid.org/0000-0002-1831-1953

E. Sheldon https://orcid.org/0000-0001-9194-0441

C. Sifón https://orcid.org/0000-0002-8149-1352

M. Smith https://orcid.org/0000-0002-3321-1432

M. Soares-Santos https://orcid.org/0000-0001-6082-8529

E. Suchyta https://orcid.org/0000-0002-7047-9358

M. E. C. Swanson https://orcid.org/0000-0002-1488-8552

G. Tarle https://orcid.org/0000-0003-1704-0781

D. Thomas https://orcid.org/0000-0002-6325-5671

D. L. Tucker https://orcid.org/0000-0001-7211-5729

J. Weller https://orcid.org/0000-0002-8282-2010

R. H. Wechsler https://orcid.org/0000-0003-2229-011X

E. J. Wollack https://orcid.org/0000-0002-7567-4451

Z. Xu https://orcid.org/0000-0001-5112-2567

## References

Abbott, T. M. C., Aguena, M., Alarcon, A., et al. 2020, *PhRvD*, 102, 023509

Ade, P., Aguirre, J., Ahmed, Z., et al. 2019, *JCAP*, 2019, 056

Aiola, S., Calabrese, E., Maurin, L., et al. 2020, *JCAP*, 2020, 047

Alam, S., Ata, M., Bailey, S., et al. 2017, *MNRAS*, 470, 2617

Amodeo, S., Battaglia, N., Schaan, E., et al. 2021, *PhRvD*, 103, 063514

Aragón-Calvo, M. A., van de Weygaert, R., & Jones, B. J. T. 2010, *MNRAS*, 408, 2163

Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, *AJ*, 156, 123

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, *A&A*, 558, A33

Bagchi, J., Sankhyayan, S., Sarkar, P., et al. 2017, *ApJ*, 844, 25

Bardeen, J. M., Bond, J. R., Kaiser, N., & Szalay, A. S. 1986, *ApJ*, 304, 15

Basilakos, S., Plionis, M., & Rowan-Robinson, M. 2001, *MNRAS*, 323, 47

Battaglia, N., Bond, J. R., Pfrommer, C., & Sievers, J. L. 2012a, *ApJ*, 758, 74

Battaglia, N., Bond, J. R., Pfrommer, C., & Sievers, J. L. 2012b, *ApJ*, 758, 75

Battaglia, N., Bond, J. R., Pfrommer, C., Sievers, J. L., & Sijacki, D. 2010, *ApJ*, 725, 91

Bharadwaj, S., & Pandey, B. 2004, *ApJ*, 615, 1

Bond, J. R., & Efstathiou, G. 1987, *MNRAS*, 226, 655

Bond, J. R., Kofman, L., & Pogosyan, D. 1996, *Natur*, 380, 603

Bond, J. R., & Myers, S. T. 1996, *ApJS*, 103, 1

Borgani, S. 1995, *PhR*, 251, 1

Bouma, S. J. D., Richter, P., & Wendt, M. 2021, *A&A*, 647, A166

Carlstrom, J., Abazajian, K., Addison, G., et al. 2019, *BAAS*, 51, 209

Carlstrom, J. E., Holder, G. P., & Reese, E. D. 2002, *ARA&A*, 40, 643

Cen, R. 1994, *ApJ*, 424, 22

Cen, R., & Ostriker, J. P. 1999, *ApJ*, 514, 1

Cen, R., & Ostriker, J. P. 2006, *ApJ*, 650, 560

Choi, S. K., Austermann, J., Beall, J. A., et al. 2018, *JLTP*, 193, 267

Clampitt, J., Sánchez, C., Kwan, J., et al. 2017, *MNRAS*, 465, 4204

Codis, S., Jindal, A., Chisari, N. E., et al. 2018a, *MNRAS*, 481, 4753

Codis, S., Pogosyan, D., & Pichon, C. 2018b, *MNRAS*, 479, 973

Coil, A. L. 2013, in *Planets, Stars and Stellar Systems*, ed. T. D. Oswalt & W. C. Keel, Vol. 6 (Dordrecht: Springer), 387

Costanzi, M., Saro, A., Bocquet, S., et al. 2021, *PhRvD*, 103, 043522

The Dark Energy Survey Collaboration 2005, *arXiv:0510346*

Darragh Ford, E., Laigle, C., Gozaliasl, G., et al. 2019, *MNRAS*, 489, 5695

Dawson, K. S., Schlegel, D. J., Ahn, C. P., et al. 2013, *AJ*, 145, 10

de Graaff, A., Cai, Y.-C., Heymans, C., & Peacock, J. A. 2019, *A&A*, 624, A48

DeRose, J., Wechsler, R. H., Becker, M. R., et al. 2019, *arXiv:1901.02401*

DeRose, J., Wechsler, R. H., Becker, M. R., et al. 2022, *PhRvD*, 105, 123520

Desjacques, V., Jeong, D., & Schmidt, F. 2018, *PhR*, 733, 1

Dey, A., Schlegel, D. J., Lang, D., et al. 2019, *AJ*, 157, 168

Einasto, M., Tago, E., Jaaniste, J., Einasto, J., & Andernach, H. 1997, *A&AS*, 123, 119

Flaugher, B., Diehl, H. T., Honscheid, K., et al. 2015, *AJ*, 150, 150

Fowler, J. W., Niemack, M. D., Dicker, S. R., et al. 2007, *ApOpt*, 46, 3444

Frisch, P., Einasto, J., Einasto, M., et al. 1995, *A&A*, 296, 611

Ghirardini, V., Bulbul, E., Hoang, D. N., et al. 2021, *A&A*, 647, A4

Giannantonio, T., Porciani, C., Carron, J., Amara, A., & Pillepich, A. 2012, *MNRAS*, 422, 2854

Giodini, S., Lovisari, L., Pointecouteau, E., et al. 2013, *SSRv*, 177, 247

Gitti, M., Brighenti, F., & McNamara, B. R. 2012, *AdAst*, 2012, 950641

Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, *ApJ*, 622, 759

Hand, N., Appel, J. W., Battaglia, N., et al. 2011, *ApJ*, 736, 39

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, *Natur*, 585, 357

Henderson, S. W., Allison, R., Austermann, J., et al. 2016, *JLTP*, 184, 772

Hill, J. C., Baxter, E. J., Lidz, A., Greco, J. P., & Jain, B. 2018, *PhRvD*, 97, 083501
Ho, A., Gronke, M., Falck, B., & Mota, D. F. 2018, A&A, 619, A122

Ho, S.-P. P., Austermann, J., Beall, J. A., et al. 2017, Proc. SPIE, 9914, 301

Hopkins, P. F., Bahcall, N. A., & Bode, P. 2005, ApJ, 618, 1

Huang, Z. 2016, PhRvD, 93, 043538

Ivezić, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, ApJ, 873, 111

Kolokotronis, V., Basilakos, S., & Plionis, M. 2002, MNRAS, 331, 1020

Kotecha, S., Welker, C., Zhou, Z., et al. 2022, MNRAS, 512, 926

Kraljic, K., Pichon, C., Codis, S., et al. 2020, MNRAS, 491, 4294

Kravtsov, A. V., Berlind, A. A., Wechsler, R. H., et al. 2004, ApJ, 609, 35

Krolewski, A., Ho, S., Chen, Y.-C., et al. 2019, ApJ, 876, 52

Kuchner, U., Haggar, R., Aragón-Salamanca, A., et al. 2022, MNRAS, 510, 581

Li, Y., Austermann, J. E., Beall, J. A., et al. 2018, Proc. SPIE, 10708, 107080A

Libeskind, N. I., van de Weygaert, R., Cautun, M., et al. 2018, MNRAS, 473, 1195

Lim, S. H., Mo, H. J., Li, R., et al. 2018, ApJ, 854, 181

Louis, T., Grace, E., Hasselfield, M., et al. 2017, JCAP, 2017, 031

Madhavacheril, M. S., Hill, J. C., Næss, J. C., et al. 2020, PhRvD, 102, 023534

Maraston, C., Pforr, J., Henriques, B. M., et al. 2013, MNRAS, 435, 2764

McCarthy, I. G., Babul, A., Holder, G. P., & Balogh, M. L. 2003, ApJ, 591, 515

McClintock, T., Varga, T. N., Gruen, D., et al. 2019, MNRAS, 482, 1352

Miller, N. & LAMBDA Group 2018, AAS Meeting Abstracts, 231, 430.05

Mroczkowski, T., Nagai, D., Basu, K., et al. 2019, SSRv, 215, 17

Naess, S., Aiola, S., Austermann, J. E., et al. 2020, JCAP, 2020, 046

Oort, J. H. 1983, ARA&A, 21, 373

Pandey, S., Krause, E., DeRose, J., et al. 2021, arXiv:2105.13545

Pessa, I., Tejos, N., Barrientos, L. F., et al. 2018, MNRAS, 477, 2991

Peebles, P. J. E. 1980, The Large-scale Structure of the Universe (Princeton, NJ: Princeton Univ. Press)

Plagge, T., Benson, B. A., Ade, P. A. R., et al. 2010, ApJ, 716, 1118

Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2013, A&A, 550, A131

Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016b, A&A, 594, A16

Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016c, A&A, 594, A13

Planck Collaboration, Aghanim, N., Arnaud, M., et al. 2016a, A&A, 594, A22

Puchwein, E., Sijacki, D., & Springel, V. 2008, ApJL, 687, L53

Regaldo-Saint Blancard, B., Codis, S., Bond, J. R., & Stein, G. 2021, MNRAS, 504, 1694

Remazeilles, M., Delabrouille, J., & Cardoso, J.-F. 2011, MNRAS, 410, 2481

Rozo, E., Rykoff, E. S., Abate, A., et al. 2016, MNRAS, 461, 1431

Rykoff, E. S., Rozo, E., Busha, M. T., et al. 2014, ApJ, 785, 104

Rykoff, E. S., Rozo, E., Hollowood, D., et al. 2016, ApJS, 224, 1

Santiago-Bautista, I., Caretta, C. A., Bravo-Alfaro, H., Pointecouteau, E., & Andernach, H. 2020, A&A, 637, A31

Schaan, E., Ferraro, S., Amodeo, S., et al. 2021, PhRvD, 103, 063513

Sefusatti, E., Crocce, M., Pueblas, S., & Scoccimarro, R. 2006, PhRvD, 74, 023522

Sehgal, N., Addison, G., Battaglia, N., et al. 2013, ApJ, 767, 38

Sehgal, N., Bode, P., Das, S., et al. 2010, ApJ, 709, 920

Sehgal, N., Trac, H., Acquaviva, V., et al. 2011, ApJ, 732, 44

Sevilla-Noarbe, I., Bechtol, K., Carrasco Kind, M., et al. 2021, ApJS, 254, 24

Shull, J. M., Smith, B. D., & Danforth, C. W. 2012, ApJ, 759, 23

Sijacki, D., Pfrommer, C., Springel, V., & Enßlin, T. A. 2008, MNRAS, 387, 1403

Sijacki, D., Springel, V., Di Matteo, T., & Hernquist, L. 2007, MNRAS, 380, 877

Sonnenfeld, A., Wang, W., & Bahcall, N. 2019, A&A, 622, A30

Springel, V., Pakmor, R., Pillepich, A., et al. 2018, MNRAS, 475, 676

Stein, G., Alvarez, M. A., & Bond, J. R. 2019, MNRAS, 483, 2236

Stein, G., Alvarez, M. A., Bond, J. R., van Engelen, A., & Battaglia, N. 2020, JCAP, 2020, 012

Sunyaev, R. A., & Zeldovich, Y. B. 1970, CoASP, 2, 66

Sunyaev, R. A., & Zeldovich, Y. B. 1972, CoASP, 4, 173

Swetz, D. S., Ade, P. A. R., Amiri, M., et al. 2011, ApJS, 194, 41

Takada, M., & Jain, B. 2004, MNRAS, 348, 897

Tanimura, H., Aghanim, N., Bonjean, V., Malavasi, N., & Douspis, M. 2020, A&A, 637, A41

Tanimura, H., Hinshaw, G., McCarthy, I. G., et al. 2019, MNRAS, 483, 223

Tejos, N., Prochaska, J. X., Crighton, N. H. M., et al. 2016, MNRAS, 455, 2662

Terry, H., Battaglia, N., Basu, K., et al. 2019, BAAS, 51, 213

Thornton, R. J., Ade, P. A. R., Aiola, S., et al. 2016, ApJS, 227, 21

To, C., Krause, E., Rozo, E., et al. 2021, PhRvL, 126, 141301

Tumlinson, J., Peeples, M. S., & Werk, J. K. 2017, ARA&A, 55, 389

van de Weygaert, R., & Bond, J. R. 2008a, in Observations and Morphology of the Cosmic Web, ed. M. Plionis, O. López-Cruz, & D. Hughes, Vol. 740 (Dordrecht: Springer), 24

van de Weygaert, R., & Bond, J. R. 2008b, in Clusters and the Theory of the Cosmic Web, ed. M. Plionis, O. López-Cruz, & D. Hughes, Vol. 740 (Dordrecht: Springer), 335

Wechsler, R. H., DeRose, J., Busha, M. T., et al. 2021, ApJ, 931, 145

Welker, C., Bland-Hawthorn, J., Van de Sande, J., et al. 2020, MNRAS, 491, 2864

Yang, T., Hudson, M. J., & Afshordi, N. 2020, MNRAS, 498, 3158

Zel'dovich, Y. B. 1970, A&A, 5, 84

Zel'dovich, Y. B., & Sunyaev, R. A. 1969, Ap&SS, 4, 301

Zhang, Y., Jeltema, T., Hollowood, D. L., et al. 2019, MNRAS, 487, 2578

Zonca, A., Singer, L., Lenz, D., et al. 2019, JOSS, 4, 1298

Zou, H., Gao, J., Xu, X., et al. 2021, ApJS, 253, 56