Superclustering with the Atacama Cosmology Telescope and Dark Energy Survey: I. Evidence for thermal energy anisotropy using oriented stacking

The cosmic web contains filamentary structure on a wide range of scales. On the largest scales, superclustering aligns multiple galaxy clusters along inter-cluster 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 multi-scale 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 for 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$\times$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.


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 Mpc, but also including nearcluster-core sizes of a few Mpc as well as some correlated structures which 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.

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  (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;Bond & Myers 1996). Meanwhile, the small-scale details are products of complex local gravitational interactions in the latetime 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,b). 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 towards 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 3-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 which consider local alignments and can be applied to a wide range of scales.

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. Measure-ments 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 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) which 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 superclustering 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 Table 1. Motivation for making localized measurements of gas anisotropy at various scales.

Object class(es)
Approx. long-axis length L Science case(s) Galaxy clusters and local surroundings L 3h −1 Mpc Filaments feeding clusters, cluster assembly processes, baryonic feedback Inter-cluster filaments, 3h −1 Mpc L 40 h −1 Mpc Cosmological model (dark energy, dark matter), small to mid-sized superclusters gas dynamics from filament compression, baryonic feedback Largest superclusters L 40 h −1 Mpc Potential tests of primordial non-Gaussianity 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 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 Kuchner et al. 2022;Kotecha 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 which 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 inter-cluster filaments and cluster pairs (∼ 4 − 12h −1 comoving Mpc, further motivated in Sec 4.2). For this scale, we focus on the redshift range where the available galaxy and cluster data are most abundant (described in Sec 2), 0.2 < z < 0.7.

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 which 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 (Zeldovich & Sunyaev 1969;Sunyaev & Zeldovich 1970, 1972. In this inverse Compton scattering process, the photons gain energy and the observed CMB spectrum shifts towards 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 < 218GHz, and increments in higher frequency maps. The strength of the effect is parameterized by the dimensionless Compton-y parameter, a line-of-sight integral of the number density of electrons n e and the electron temperature T e . σ 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 catalogues 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 c and angular-integrated Compton-y parameter, Y , which is proportional to P e V c (the electron pressure times cluster volume), is fairly consistent with Y ∝ M 5/3 , as occurs in a gas with adiabatic index γ = 5/3 in equilibrium (see Bond & Myers 1996;Mc-Carthy et al. 2003;Giodini et al. 2013, and references therein.) However, for lower-mass halos (M ≤ 10 14 M )such as those that host galaxy groups along inter-cluster 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 P e 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 lowsignal 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 (SNR) than any of individual images by a factor of √ 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 inter-cluster 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. Clusterpair 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 Mpc 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 multi-scale information content of tSZ anisotropy.

Oriented stacking
Oriented stacking is a more general approach to stacking elongated structures which 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 Figure 1. A visualization of unoriented (top) versus 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.
paper presents its first application to a 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 large-scale structure. 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 multi-scale 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 cos-mology, 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 superclustering. 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 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 superclustering 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.cosmology 1 package which implements a Flat ΛCDM model with parameters from the Planck Collaboration et al. (2016c). The model has Ω M = 0.308, Ω k = 0, Ω Λ = 0.691, a single species of massive neutrinos with mass 0.06 eV, and H 0 = 67.7 km s −1 Mpc −1 . All quoted distances are in comoving units.

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 observational results to simulations of large scale structure.

ACT Compton-y map
ACT is a 6-meter off-axis Gregorian telescope, located in the Atacama Desert of northern Chile, at an eleva-tion 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 ACT-Pol receiver was deployed, enabling ACT to observe both temperature and polarization data at 98 GHz 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;. ACT has recently produced maps covering 18,000 square degrees of the sky Aiola 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 square degree 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 night-time 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 < < 300 for the Planck low-frequency instruments (LFI) at 30 GHz and 44 GHz respectively, 20 < < 2000 for the 70 GHz LFI and modes between 20 < < 5800 for the high frequency instruments (HFI) at 100, 143, 217, 353 and 545 GHz. Data from ACT include modes 500 < < 24000. 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 arcminute 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 sub-dominant 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 Compton-y. For our key results in this paper (Sec. 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-sigma errorbars. 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 nony 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 Sec. 2.3 and collected by a different survey (DES).
A complete description of the y map procedure is found in Madhavacheril et al. (2020) 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.

Galaxy cluster data
We stack cutout images from the ACT Compton-y map on locations of galaxy clusters identified in the Dark Energy Survey (DES) Y3 data. DES (The Dark Energy Survey Collaboration 2005) recently completed a six year survey (2013-2019) of 5,000 square degrees of the southern sky in five optical filters (grizY ). The survey was conducted by the 4-meter Blanco Telescope, fitted with the Dark Energy Camera (Flaugher et al. 2015), at the Cerro Tololo Inter-American Observatory (CTIO) in Chile. The redMaPPer algorithm, 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 . The algorithm determines a value for richness, λ, for each cluster by summing the membership probability of each galaxy which has some likelihood of belonging to that cluster, within a defined radius. λ 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 three years of the survey (Sevilla-Noarbe et al. 2021). This redMaPPer catalog (titled v6.4.22+2 Full) extends out to z ∼ 1, includes all λ > 5 clusters, and provides a photometric redshift (photo-z) for 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 degrees 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 degree to the edge of DES galaxy data (discussed in the 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 λ > 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 Comptony signal from the LSS surrounding higher-λ clusters is stronger, and (b) that higher-λ 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 λ < 20 due to their known impurities (e.g. Abbott et al. 2020;To et al. 2021;Costanzi et al. 2021). These low-λ 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 miscentering . However, in our study, a λ > 20 cutoff is not feasible. There are not enough λ > 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 λ threshold, we test λ > 10, λ > 15, and λ > 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 λ > 10. Since this cutoff results in the highest SNR, 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 λ (Zhang et al. 2019) and we will examine signals beyond 1.5R λ . However, the effects of false cluster detections may be nonnegligible, and we would expect them to bias our results lower. We discuss this further in the conclusions.
There are 5,494 clusters in the remaining sample. The median photometric redshift uncertainty is σ 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 λ > 10 sample. The figure also shows three colored lines which represent the different richness thresholds tested in simulations in Section 5.2.

Galaxy data
We use galaxy data from DES to orient each Comptony map cutout with respect to the axis of surrounding elongated structure. The redMaGiC algorithm selects Luminous Red Galaxies from photometric surveys using a matched-filter technique ). It has recently been applied to Y3 DES data over the full DES footprint, as detailed in Pandey et al. (2021). The algo- 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 towards 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 SNR on extra-cluster gas. In section 5.2, we use simulations to test 3 richness cutoff values (10,15,20), shown here in blue, purple, and green.
rithm is designed to minimize errors in galaxy photo-z, resulting in an average scatter of σ z /(1 + z) ∼ 0.013. In this study we use the High Density catalog, which covers redshifts 0.1 to 0.7 with fairly consistent number density (∼ [8.5, 10.5] × 10 4 h 3 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 × 10 13 h −1 M (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.

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×DES results. Buzzard assigns realistic galaxies to dark matter halos in N -body dark-matter only simulations using the ADDGALS method (Wechsler et al. 2022). The implemented cosmology has Ω m = 0.286, h = 0.7, σ 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 massrichness 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. 2022). Despite the lower cluster number density, the number of Buzzard clusters available to use for the noiseless theory cmoparison 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 ∼ 12x 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 from Stein et al. (2020), described in more detail in Sec. 2.5. We apply the model down to halo masses of 10 12 M , 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 arcminute beam and remove modes with < 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 uncer-tainties 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 Sec. 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.

Websky Simulations
We wish to use simulations not only to make direct comparison with ACTxDES, but also to generate puretheory expectations for the superclustering of thermal energy over time. The Buzzard redMaGiC and redMaP-Per 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: Ω m = 0.31, h = 0.68, and σ 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 large-scale structure 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 ∼ 2, which allows us to go well beyond the DES redMaGiC limit of z ∼ 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 arcminute beam and remove all power for modes < 20.
We select cluster-mass halos in a few mass ranges and use mass-weighted halos to provide information about large-scale structure, as further detailed in Section 5.

GALAXY FIELD CHARACTERISTICS
We begin by studying large-scale properties of the projected galaxy number density field. The bulk of this study is implemented with the Cosmology Object Oriented Package (COOP 3 , 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 which satisfy these criteria using projected galaxy overdensity maps in bins of redshift. The overdensity is defined as δ g = n g /n g − 1, where n g is the two-dimensional number density of galaxies andn g is its mean. In practice, each map is created in the Healpix 4 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 δ 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 fullwidth at half-maximum by FWHM=2 √ 2 ln 2R G . We vary the smoothing scale to observe the LSS properties at a range of scales. Specifically, we examine results for Gaussian smoothing with full-width at half-maximum (FWHM) ranging from 6 to 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 ∼ 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 a FWHM of 14 Mpc for demonstration purposes.
The primary property we use to determine superclustering is the field excursion ν, where F is the field value at some position and σ is its root mean square (RMS) after smoothing on some scale R ( Bardeen et al. 1986). Given that the δ g maps have a mean of zero, σ is equivalent to the standard deviation, making ν a measure of signal-to-fluctuation-noise. Points with higher ν 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 Gaussiansmoothed 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 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 −λ i and corresponding eigenvectors. Note that with the negative  sign, λ i are defined to be positive at peaks and negative at troughs. We order the eigenvalues as | λ 1 |>| λ 2 |, such that λ 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: 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 ν and e parameters are visualized in Fig. 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 λ 2 is 0, and consequently e ≤ 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 non-zero. These large-scale gradients point toward supercluster centers rather than toward the individual clusters that make up each supercluster. For some clusters, where the field has two equal-sign eigenvalues, e ≤ 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 towards 0, extends well beyond 0.5 and can become very large when the denominator of equation 4 is small (in other words, when λ 2 ∼ −λ 1 ). We find that points with higher ν, 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: where σ 2 is the root mean squared value of ∇ 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 ∝ (λ 1 + λ 2 ), which is the denominator of e, so points where λ 1 ∼ −λ 2 have small x and large e. Due to the divergence of e at x ∼ 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 ν, |x|e, and x distributions of the galaxy overdensity field with that of a Gaussian random field we have constructed with the same power spectrum. We divide the Buzzard red-MaGiC 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 gg (Pandey et al. 2021), we generate a Gaussian Random Field (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 λ = 15 cluster using the mass-richness rela-tion 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× the mean matter density of the universe at that redshift. λ = 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 discussion on top-hat to Gaussian conversion, see Section 4.2.) We smooth the pseudo-galaxy field with a Gaussian function with FWHM ∼ 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 ν threshold on the small peak scale. We choose a ν 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 (ν, x, |x|e) for the galaxy overdensity fields smoothed on larger scales. Figures 5 (6) show the distribution for ν (x) versus xe 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 ν and x distribution which is more stretched compared to the GRF. The distributions are both skewed towards the high end, displaying more high-ν 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 ν. 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 ν and e threshold for the RedMaPPer cluster sample enhances the supercluster gas signal in stacks. We choose ν > 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 coarsegrained galaxy maps: effectively, regions of strong superclustering. In Section 4.6, we examine the effect that Figure 5. The distribution of |x|e (field elongation) versus ν (field excursion) for ∼ 9, 500 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 Gaussian Random Field with the same galaxy power spectrum (C gg ) as Buzzard. The maps have each been Gaussian smoothed with FWHM=14 Mpc. The entire cluster / peak sample consists of the red plus blue dots. The one, two, and three-sigma contours are drawn in black, and the points which are included after our chosen cuts of ν > 2 and e > 0.3 are colored in blue. The cutoff is sharp at ν = 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.
these constraints have on the non-Gaussianity of oriented stacks.
We also find that the ν and e properties on the 14 Mpc scale are not highly correlated with cluster richness λ, a small-scale property. Figure 7 demonstrates this: despite the weak correlation, the remaining clusters after the cut ν > 2 are still distributed across a wide range in λ. 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 ν value. For e (not shown), there is a weak anti-correlation: clusters in very high-ellipticity regions tend to be 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 Figure 6. The distribution of xe (field elongation) versus x (signal-to-noise of the curvature) for galaxy overdensity maps from DES observations (top), Buzzard mocks (middle), and a Gaussian Random Field (bottom). The chosen cuts which are used throughout much of this paper, ν > 2 and e > 0.3, are shown in blue. The ν cut shifts the x distribution rightwards because rarer overdensities also tend to have higher second derivatives. The e threshold shows up here as a diagonal |x|e threshold. Figure 7. The relation between ν, as measured on the galaxy overdensity field smoothed at 14 Mpc, and cluster richness λ. There is slight correlation, such that the ν > 2 clusters have a higher-skewed λ distribution, but the constrained sample still features a range of λ.
large due to its normalization. There is no correlation in the non-normalized |x|e property.
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.

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 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).
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 photo-z uncertainties are σ z ∼ 0.01 − 0.02 (∼ 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 Sec. 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 λ 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, inter-cluster filaments, and superclusters as the map is increasingly smoothed.
For each cluster, we take a 4 • × 4 • cutout from the ACT Compton-y map centered on the cluster (RA, dec). At the redshift range explored in this work, this corre- 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 signal-to-noise ratio scales as √ 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 large-scale structure information around each cluster determines its orientation. The number density is relatively constant until z ∼ 0.7, the imposed edge of this volume-limited sample. We cut the data at the dropoff. We also cut the lowest-redshift data because 4x4 degree cutouts do not span a large enough physical extent to observe the full effects from superclustering at these redshifts.
sponds to a coverage of ∼ 80 Mpc (closest slice) to ∼ 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 SNR of the final image is proportional to √ 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 5 (Huang 2016).

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 which revealed filament signal from pairs separated by a transverse distance of 6 − 10h −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 ∼ 6 − 10h −1 Mpc. A simple conversion from tophat to Gaussian filter involves equalizing the volume under each function; with a top-hat of height 1 and a Gaussian with amplitude 1, this leads to FWHM∼ 1.7R TH . Thus, for cluster pairs separated by, e.g., ∼ 7h −1 Mpc, Gaussian smoothing for which the Hessian would best encode alignment at that scale has FWHM∼ 12h −1 Mpc or ∼ 18 Mpc.
We therefore choose 18 Mpc as our key scale to study with ACT×DES data. We also select a range of smaller smoothing scales, down to 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 inter-cluster 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 Mpc is held constant across all redshifts.

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 signal-to-noise 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: ∼ 0.25 < z < 0.72, or ∼1000-2600 Mpc.
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. 5 https://www.cita.utoronto.ca/ ∼ zqhuang/work/coop.php 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 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.

Multipole Decomposition
To quantitatively compare the stacked images, we decompose each image into its multipole components. Each image I can be deconstructed as Since symmetry along the x-axis is enforced in the oriented stacking method, the odd components trend towards 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 where X is 2 for m = 0 and 1 for all other m. The multipole decomposition is visualized 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). . Left: a 4x4 • sky patch from the observed (above) and simulated (below) Compton-y maps. At this redshift, the corresponding physical size is ∼ 180Mpc. 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 FWHM=14 Mpc, in the same sky area. The same clusters are circled and black arrows indicate their orientation with respect to the surrounding large-scale structure, as determined by the Hessian matrix on Fg. Right: Oriented stacks made by a combination of multiple distance slices in the range 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 δg before stacking. This builds up signal-to-noise along the horizontal axis, showing the gas signal from superclustering.

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, Σ, of each C m (r) profile. If a profile is binned into radial bins r i , each element Σ ij is the covariance between the signal in 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 to 48 depending on the data set and cluster selection. Cluster samples in all the regions are approximately equal-sized and collected in (RA, Dec) space into patches 30 − 80 deg 2 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 large-scale structure 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 from 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 depends on the data set. For Websky full-sky maps, we 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 Cm(r). The odd m moments trend towards zero with increased numbers of stacked cutouts. For the even moments, we show the symmetrized representation Cm(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.
split the sample into 48 Healpix pixels. For the Buzzard and DES data, we use the kmeans radec algorithm 6 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 a 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 p th region and N cl 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 6 https://github.com/esheldon/kmeans radec then the covariance matrix element between bins r i and r j is The left bracketed term is a normalization by the degrees of freedom 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 p th region, and the result is matrix-multiplied with the transpose of M (generally, (M × M T ) mn = l M ml M nl ). In practice, we calculate Σ with the NumPy cov 7 function. Σ as defined above represents an estimate of the covariance of a single region's y profile which was '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 sub-regions 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 affects the resulting χ 2 and signal to noise 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: 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 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-instrument noise, and highresidual foregrounds. We subtract the tSZ realization from each map so that only the noise contribution remains. We estimate Σ 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 map = 120 independent stacks. Because the number of images per stack is constant, Equations 8 and 9 simplify to: We estimate Σ signal , the covariance matrix of the signal component of the measurement, as being equal to the Buzzard map-splits covariance matrix (calculated by Eq. 9) scaled by N Buzz /N DES (the relative number of clusters in Buzzard versus DES). We compute Σ 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 Σ tot = Σ signal + Σ noise and examine which component contributes more to the total. Σ noise dominates the variances of the signal in each radial bin, as its diagonal is ∼ 3 − 6× larger than the diagonals of Σ 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 Σ 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×DES and Buzzard in Section 6, we choose the binning more carefully. We determine the convergence of Σ 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 Σ 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×DES sky coverage.
We note that the DES photo-z 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: Σ signal includes variance across small samples of different photo-z realizations for galaxies and clusters within Buzzard. Buzzard reproduces the photo-z uncertainties of the DES catalogs well, shown in D19.
The redshift-dependent bias, however, is unaccounted for in the uncertainties. Due to photo-z 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, ∼ 25% of galaxies and ∼ 20% of clusters are interlopers in a given slice. The combination of the peaked distribution and the inclusion of interlopers has a non-trivial 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 σ z of galaxies / clusters in DES. We then compare the m = 2 radial y profiles from a stack using photo-zs 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 ∼ 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 vs. ACT×DES. Due to Buzzard's accurate reproduction of the photo-z 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×DES.

Comparison with a Gaussian Random Field
Can oriented stacking be used to examine non-Gaussian structures in the late-time universe which result from non-linear evolution? In the cold dark matter 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 non-linear 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 non-linear 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 ∼ 1.5× at z = 1 and k = 1 h Mpc −1 , and by ∼ 6× 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 Gaussian random field.
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.
Next, we stack on 'peaks' in the GRF that roughly match the size of λ > 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 a visualization of the spatial distribution of the full cluster sample (red circles) and the clusters constrained by ν and e limits (over-plotted 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 psuedo-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 δ g and Comptony 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 psuedo-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 Figure 11. Top: a comparison of the m = 0, m = 2, and m = 4 moments of oriented stacks on Buzzard clusters (solid lines) and Gaussian random field 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 ν 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 versus 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 FHWM=14 Mpc. Red circles are plotted around all cluster locations. Many clusters are cut by imposing ν > 2, e > 0.3 and the remaining ones have blue circles over-plotted (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 ν and e thresholds select clusters in high-superclustering, highly non-Gaussian regions.
imposing ν > 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.

Websky
We use the Websky simulations to give pure theoretical results which are not subject to the selection effects of DES and the Buzzard mocks. Two key questions are: • How does superclustering depend on redshift?
• 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 sub-samples, which increases noise in the stacks. Due to the full-sky coverage of this simulation, there are enough clusters that the SNR 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 12 , 1 × 10 15 ] M . This range incorporates most Redmagic-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;Sonnenfeld et al. 2019;Maraston et al. 2013). The lower limit is slightly higher than the Websky minimum halo mass (∼ 1.2 × 10 12 M ) because cutting out the lowest-mass halos saves computational time. Nevertheless, the wide halo range paints a near-complete picture of the largescale structure, while the redMaGiC galaxy data only form a subset of the galaxies which trace the underlying matter. We create mass-weighted halo number-density and overdensity maps using a weight of M h /10 12 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 ∼ 0.25 to z ∼ 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, demonstrating that the anisotropic tSZ signal increases from a low level at high redshifts (early times) until z ∼ 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σ error 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 superclustering 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 towards clusters. These mergers are expected to heat up the gas in galaxy groups. In addition, 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- Figure 12. The m = 2 (top) and m = 4 (bottom) radial profiles for Websky oriented stacks in progressive 400 Mpc slices in line-of-sight distance. These stacks of massive clusters (M > 10 14 M ) were oriented using a mass-weighted halo density field smoothed at the scale of FWHM=14 Mpc. The profiles are binned equally in arcminutes, thus the bin sizes and centers in Mpc 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 towards 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.
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 super-clustering 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 richness range (λ) at z = 0.5 assuming the mass-to-richness relation from McClintock et al. (2019). After stacking the cluster cutouts, we combine stacks from 5 slices from 1200-1800 Mpc to get the images shown in Figure 13. Figure 14 shows the m = 0, 2, 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.
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 redMaGiC 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.

Dependence on parameters
We examine the dependence of superclustering y signal on cluster richness λ, field excursion ν, field ellipticity e, and smoothing scale using the Buzzard mocks. Figure 15 shows the 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 oriented y signal, yet a reduction of factor N to the number of stacked images augments the random noise by √ N . For the different parameter values tested, we assess this trade-off by examining the signal-to-noise of the maximum bin of m = 2. We find that the full λ > 10 sample has the highest signal-to-noise. If including the field constraints, a combination of (λ > 10, ν > 2, e > 0.3) is optimal.
Next we examine changes to the smoothing scale that is applied to the δ 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 nearestneighbor clusters) more important to the determination of curvature. With the Buzzard simulations, we perform oriented stacking with a Gaussian-smoothed map of 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 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 1232-1832 Mpc, or z ∼ 0.4. Figure 14. m = 0, 2, 4 radial profiles for the oriented stacks of three mass bins of Peak Patch halos shown in Fig 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.
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 tempera-ture and density of aligned structure at that scale. We will further analyse 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, where C m (r) is defined in Equation 7. For m = 0, this is similar to the unitless angularly-integrated Comptony parameter, with the difference that our Y has units of Mpc 2 . 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 ν 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.

COMPARISON WITH OBSERVATIONAL DATA
We apply the method to ACT×DES data and the Buzzard mocks using a few combinations of smoothing scales and field constraints. We will highlight results for m = 0 m = 2 m = 4 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 emin only has an impact on the anisotropic moments. The strongest effect to 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.
18 Mpc smoothing, a scale which roughly corresponds to inter-cluster filaments (as motivated in Sec. 4.2). We later show all three smoothing scales. Each result figure shows the observational y measurements in 3 radial bins, as motivated in Sec 4.5. Each figure also shows the Buzzard results, which demonstrate what a 'pure' measurement would look like across the full DES footprint, with errorbars due to variance in the large scale structure but not instrumental noise. The Buzzard results are binned into the same 3 radial bins to assess the consistency with observations, while the 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 λ , where R λ is the redMaP-Per cluster radius (Rykoff et al. 2014). R λ is given by (λ/100) 0.2 h −1 physical Mpc, which we calculate at the median richness in our sample, λ = 15. After conversion, this results in a comoving radius of ∼ 2.5 Mpc for the average redshift of the cluster sample. Therefore, 45 Mpc 90 Mpc 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 large-scale structure, 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.
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 longwavelength noise in the y map from residual low-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-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 to 40 Mpc (the tail of each profile) from the full profile. This is similar to performing aperture photom-   Fig 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 maximum anisotropic gas signal from the cluster, but rather a property which depends almost exclusively on the smoothing scale.
etry, 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 ν > 2 constraint is applied. (The e constraint has little to no effect on the isotropic profile, as previously demonstrated in Figure 15). The ν and e constraints reduce the cluster sample by a factor of ∼ 5, causing an increase in noise of ∼ √ 5. We can determine the signal-to-noise of the 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 were separated by > 40 Mpc from the nearest object containing hot gas. We compute the reduced-χ 2 of the observational data vector y obs (the coarsely binned C 0 (r) profile) with respect to null: where Σ −1 is the inverse covariance matrix of the data. Next, we find the probability that a truly null vector could exceed that χ 2 . This Probability to Exceed (PTE) is measured by drawing 1 million random vectors from a normal distribution with mean zero and the covariance matrix from the data. We compute the χ 2 of each vector with respect to null, then find the fraction of the sample for which the χ 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 a detection of non-zero signal in the data. We then relate this to a signal-tonoise (SNR) estimate which is the number of Gaussian sigmas away from null, where erf is the error function. We also assess the goodness-of-fit to the Buzzard results by finding the reduced χ 2 : where Σ o+s is the summed covariance matrix for observational data and simulations, and y sim refers to the coarsely binned Buzzard profile. The bins are highly correlated for m = 0, so the SNR is only 1.7 and 2.3 for the profiles without and with constraints, respectively. Both profiles agree well with Buzzard, with χ 2 red < 1. The low SNR 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 m = 4. In addition, we achieve a higher-SNR 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 λ > 10 clusters in Buzzard peaks at y = 4.5 × 10 −8 . Enforcing (ν > 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.
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 large-scale structure. Table 3 presents the reduced-χ 2 , PTE, and SNR 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 errorbars compensate for the boost in 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), (b), and (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 can not claim evidence for the m = 4 component in 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 ν and e contraints enforced for the respective smoothed maps. The m = 2, 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 χ 2 red column of Table 3.
There is evidence for m = 2 signal at the 3.2σ level for 14 Mpc smoothing and marginal 2.6σ evidence for 10 Mpc. The SNR 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 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-χ 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 Compton-y map stacks. The orientation and ν and e selection were performed using a Gaussian smoothed galaxy field with FWHM=18 Mpc. The background shaded curves surround 1σ regions above and below the continuous mean profile from Buzzard; the triangles correspond to 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. 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.
in Section 4.6, the boost from applying the ν and e thresholds only occurs in realistic fields and not in a Gaussian random field. Therefore, it is likely 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 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, SNR gives the number of Gaussian σ from null, and χ 2 red is the reduced χ 2 value using 3 degrees of freedom (d.o.f.). χ 2 red is shown for the absolute value with respect to 0 as well as the value with respect to Buzzard. The highest signal-to-noise detections of extended signal come from the m = 2 components of the 14 Mpc and 18 Mpc smoothed stacks. Data and simulation are generally in agreement, with χ 2 red (data-sim) consistently near or below 1. 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 SNR=[1.0, 2.0, 2.7] for the [10,14,18] Mpc smoothing scales. Therefore, we find marginal in-  Table 3 for the Buzzard simulations, with the same 3 radial bins. χ 2 red is only shown for the 'absolute' value, as the (data-sim) value is already presented in Table 3. Because the simulation results are high-SNR, rather than calculating a PTE we simply estimate SNR by taking the square root of χ 2 . 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. dications of integrated anisotropic clustering of thermal energy for the two coarser scales. 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 simulated stacks. There is a marginal indication for 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 The Buzzard finely-binned data are shown for visual purposes, but is 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 ν > 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. a distinction from a Gaussian random field. 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 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σ. data, we found marginal-to-significant evidence for 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 m = 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 × 10 −9 . The team also found that bound gas in halos only contributes ∼ 20% of this signal. A WHIM signal at y ∼ 5 × 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 large-scale structure 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 ∼ 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 ∼ 2 to z ∼ 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.

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, λ < 20 redMaPPer clusters are known to contain false detections. False clusters should be less correlated with the surrounding large-scale structure 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 λ > 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 large-scale structure. We found in section 5.2 that cluster richness is indeed correlated with the anisotropic Comptony signal; however, this is subdominant compared to the dependence on ν evaluated at larger scales. Additionally, because the mock Compton-y map is created using only halo information from the Buzzard dark-matteronly 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-ofthe-art: a recent version of Buzzard improves the colordependent clustering which 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 modelling/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 sq. deg. footprint ). 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 m = 2 signal after ν 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 χ 2 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 superclustering 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, expected to begin taking CMB observations in the mid-2020s. With measurements in 6 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 Tele-scope (fomerly 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 RedMaP-Per 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 (DESI, 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 (LSST, 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 signal-to-noise and redshift extent of superclustering measurements.
This paper demonstrates a flexible tool for examining superclustering in the universe. We have shown its effectiveness 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 clus- 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)