Boundless baryons: how diffuse gas contributes to anisotropic tSZ signal around simulated Three Hundred clusters

The distribution and state of baryons in the cosmic web has become an increasingly important puzzle in cosmology and astrophysics. The fractional contribution of baryons to the total energy budget of the Universe today, $\Omega_{\mathrm{b}}$, can be predicted from early-Universe cosmic microwave background (CMB) data assuming a cosmological model, such as the current concordance model $\Lambda$CDM (e.g., Hinshaw et al., 2013; Planck Collaboration et al., 2020; Aiola et al., 2020). A full census of baryons via direct detection in the late-time Universe would be an excellent test of the cosmological model. However, for years, the estimated baryon fraction from low-redshift observations fell far below the prediction from the early Universe, motivating a search for the ‘missing baryons’ (Cen & Ostriker, 1999).

Theory paints the following picture of the history of baryons. As the Universe evolved from its initial nearly-uniform state, baryons followed the dark matter through linear and then non-linear structure growth, populating the filamentary cosmic web and evacuating void regions. Intergalactic hydrogen at $z\gtrsim 2$ can be observed via the Lyman-$\alpha$ forest, and with these observations the predicted baryons are largely accounted for. However, the intergalactic hydrogen became increasingly ionized and shock-heated by structure formation, while also becoming more diffuse due to cosmological expansion. Additionally, active galactic nuclei (AGN) and supernovae (SN) began to eject and re-heat gas that had previously fallen into galaxies and cooled. The combination of low densities ($n_{b}\sim 10^{-6}-10^{-5}$ cm^-3), intermediate temperatures ($10^{5}-10^{7}$ K), and ionization make such gas challenging to observe, hence the moniker of ‘missing’ baryons. By $z=0$, simulations predict that $\sim 30-70\%$ of the baryonic mass exists in this so-called warm-hot intergalactic medium, or WHIM (Cen & Ostriker, 1999; Davé et al., 2001; Cen & Ostriker, 2006; Penton et al., 2004; Shull et al., 2012; Nicastro et al., 2008; McQuinn, 2016; Cui et al., 2019; Sorini et al., 2022).

The wide range in these predictions is due to various evolutions of the WHIM which occur depending on the implemented feedback model (Davé et al., 2001). Constraining the evolution of the WHIM is important not only for the cosmological energy census, but also for its potential as a tracer of low-redshift cosmic web structure, a playground for testing the behavior of dark energy and dark matter. The sensitivity of the WHIM to galactic feedback processes means it can also be used to constrain the physics of AGN and SN, an area of high interest for cosmology–e.g., the redistribution of baryons from astrophysical feedback is considered to be one of the most important unknowns for utilizing the full power of current and future lensing surveys (Harnois-Déraps et al., 2015). While feedback presents a nuisance for cosmology, it is interesting for enhancing our understanding of black holes and stellar evolution (Kormendy & Ho, 2013; Yuan & Narayan, 2014; Somerville & Davé, 2015). The states of the circumgalactic and intergalactic media are of interest in their own right (see Meiksin, 2009; Tumlinson et al., 2017, for reviews), and have important relationships with galaxy evolution (Oppenheimer & Davé, 2008; Davé et al., 2011).

The WHIM induces very weak signals in maps of the CMB through the Sunyaev-Zel’dovich (SZ) effects (Sunyaev & Zeldovich, 1972). The thermal Sunyaev-Zel’dovich effect is sensitive to gas pressure along the line-of-sight; it is parameterized by the Compton-$y$ parameter,

$y$ is a line-of-sight integral of the number density of electrons $n_{\mathrm{e}}$ and the electron temperature $T_{\mathrm{e}}$. $\sigma_{\mathrm{T}}$ is the Thomson cross section and $m_{\mathrm{e}}$ is the electron mass. Due to the dependence on both density and temperature, signals from galaxy clusters are orders-of-magnitude higher than signals from intergalactic filaments, which are well below the noise in current maps. Nevertheless, maps have recently become sensitive enough for significant detection of tSZ signals from an intermediate regime: a cluster pre-merger bridge (Planck Collaboration et al., 2013; Bonjean et al., 2018; Hincks et al., 2022). Meanwhile, the kinematic SZ effect, depending on gas density and bulk velocity, is significantly smaller than tSZ and more difficult to reconstruct. However, recent advances (Hill et al., 2016; Ferraro et al., 2016; Kusiak et al., 2021; Schaan et al., 2021) have forwarded the kSZ as a promising WHIM probe, especially in combination with tSZ measurements. We will focus only on tSZ in this work.

Several studies have developed new stacking techniques for teasing out the weak tSZ signals from filaments using combinations of CMB and galaxy survey data. Stacking averages together $N$ images of the same or similar objects, which contain signal in the same location but have varying noise, to create a single image in which the signal-to-noise ratio (SNR) is enhanced by $\sqrt{N}$. By stacking a $y$ map from Planck satellite data on locations of pairs of galaxies, de Graaff et al. (2019) and Tanimura et al. (2019) found evidence for tSZ signal from a stacked intercluster filament. The former group posited that $\sim 80\%$ percent of the $y=(0.6\pm 0.2)\times 10^{-8}$ filament signal came from gas beyond dark matter halos. Lokken et al. (2022, hereafter L22) used a distinct oriented stacking method to measure the large-scale anisotropy of tSZ around clusters in Atacama Cosmology Telescope (ACT) data, detecting a signal from the combined emission of neighboring halos and filament gas at the 3.5$\sigma$ level and demonstrating that the signal is stronger in supercluster regions, as defined by galaxy data.

These methods have been applied only to data from $z<1$, where there is sufficient galaxy number density in the current state-of-the-art large galaxy surveys (e.g., the Sloan Digital Sky Survey (SDSS) and Dark Energy Survey (DES), Bolton et al., 2012; Pandey et al., 2021; Porredon et al., 2021). All oriented stacking methods will improve with higher galaxy number density. As upcoming surveys map the cosmic web of galaxies to higher redshifts with greater number density and precision (e.g., the Vera C. Rubin Observatory and the Dark Energy Spectroscopic Instrument, Ivezić et al., 2019; DESI Collaboration et al., 2016), and anisotropic statistical methods improve, it may finally be possible to observationally constrain the WHIM evolution from $z\sim 1$ to $z=0$.

To motivate such studies, in this work we use oriented stacking to examine tSZ signals from the diffuse gas in hydrodynamic simulations at $z=1$. Following the applications to both simulated and observed data in L22, which determined the orientation of filaments and superclusters on scales of $\sim 10-20$ comoving Mpc, we conduct a similar study on the smaller end of this scale range using the Gizmo-Simba run of The Three Hundred  project¹¹1https://the300-project.org/, part of the data are available through this website: https://the300.ft.uam.es/, upon request. (the300 for short, Cui et al., 2018). In brief, these simulations are run with Simba hydrodynamics (Davé et al., 2019), include both SN and AGN feedback, and re-simulate reasonably large (15 $h^{-1}$Mpc comoving) regions surrounding 324 different massive clusters with a competitive resolution ($\sim 10^{8}$ $h^{-1}{\rm{M_{\odot}}}$). Thus, the300 runs provide an excellent option for studying the filamentary structure around clusters due to the sufficient numbers of clusters for stacking (competitive or better than the numbers formed in large cosmological simulations) and the larger extent than typical cluster zoom simulations. Furthermore, their nature as zoom simulations makes it feasible for future work to re-simulate selected regions with different feedback models. We focus our study at $z=1$, the new frontier of galaxy-CMB cross-correlation studies, where the proto-cluster regions are still coalescing along complex filamentary structures. The detailed nature of the simulations enables us to address a key question: how much tSZ signal does diffuse gas contribute to filaments?

We begin by dividing the particles into diffuse and concentrated categories by several methods. The fiducial method applies a threshold at the halo virial radius $R_{200c}$ (where the density enclosed is 200 times the critical density), considering the particles within this radius as halo-associated and those beyond as diffuse. Compton-$y$ maps are then created for the different particle sets, oriented using information from the galaxy field, and stacked. We examine the stacks to determine the fractional contribution to the $y$ signal from halo versus diffuse gas.

There are two main motivations for dividing the halo from the diffuse gas in the context of oriented SZ studies. First, the tSZ effect has the potential to be a probe of cosmological structure beyond 2-point statistics. With advances in oriented stacking, the tSZ can start to provide information about the cosmic web shape and structure beyond the densest, highest-pressure nodes. However, a cosmological probe is only as good as its respective modelling. To be able to simulate the effects of varied cosmologies, the best approach is to use analytic theory, N-body, or rapid predictive simulations for the dark matter, and post-process to apply a prescription for the baryons (e.g., as done in Stein et al., 2019; Pandey et al., 2019; Raghunathan, 2022). The alternative – many large-volume hydrodynamic simulations with varied cosmologies – is currently computationally infeasible, although it is being tested in smaller volumes by the Camels project (Villaescusa-Navarro et al., 2021). The leading prescriptions for tSZ are based on the halo model, assigning baryons to halo locations with a pressure profile which is well-understood to some radius. However, it is unclear whether this prescription sufficiently captures the tSZ signal, especially for locally-anisotropic studies such as oriented stacking, due to its lack of description for the diffuse gas beyond halos. Our study asks: is it possible to capture the majority of the extended tSZ signal by using halo modeling? If so how far, and at what halo masses, need to be accurately modeled? This will guide future oriented SZ studies: if the unassociated gas is unimportant, then the focus should remain on perfecting the halo pressure profile models in order to be able to use the tSZ to probe beyond-two-point statistics. If it is important, more attention should be devoted to modelling the field component of the tSZ.

The second major motivation is to guide studies which seek to find the ‘missing’ baryons. For the observational analyses that aim to do this with tSZ, the contribution from known baryons must be subtracted: e.g., the $y$ map should be masked at halo positions out to some radius, or the contribution from halos should be modelled (as done in de Graaff et al. 2019; Tanimura et al. 2019). Thus it is useful to know what fraction of $y$ is captured by halo models and what residual is expected in the current best hydrodynamical simulations.

The answers to these questions depend on the implemented feedback model and thus will be specific to the300 feedback mechanisms. Analyzing how the diffuse $y$ signal changes for different feedback prescriptions is an interesting and important question for future work. Nevertheless, it is still useful to explore the questions given the particular setup of the300, which is one of the leading simulations currently available with enough volume for studies of extended tSZ. Even an approximate estimate for the importance of diffuse gas in tSZ studies will help motivate analysis with upcoming tSZ and galaxy survey data.

Throughout the paper, we use the cosmological parameters from the Planck Collaboration et al. (2016), also used in the300 simulations. All quoted distances are in comoving Mpc. Sec. 2 describes details of the simulations. Sec. 3 explains the methods including selection of a central halo for the stacks, reduction of the snapshot sample, splitting of the halo-associated and diffuse gas, creation of $y$ maps, and orientation of the maps using galaxy data. In Sec. 4, we present the stacked maps and analyse their radial profiles in a multipole decomposition. We also discuss the association of gas with halos and compare the filament to off-filament axes. In Sec. 5, $y$ maps are sub-selected based on measures of density and elongation, and the resulting differences in signal are analyzed. Sec. 6 presents discussions and conclusions.

The cluster samples used in this study are coming from the300 Gizmo-Simba simulations. the300 selects the most massive galaxy clusters from the MultiDark simulation (MDPL2, Klypin et al., 2016), which is a series of simulations that cover a range of masses from $10^{10}M_{\odot}-10^{15}M_{\odot}$ and volumes up to $50~{}\mathrm{Gpc}^{3}$. The simulation re-generates the initial conditions using the Ginnungagap code²²2Code available online: https://github.com/ginnungagapgroup/ginnungagap with a zoomed-in technique, which allows us to only simulate the cluster regions with different hydrodynamic codes. The zoom regions are extended to 15 $h^{-1}$Mpc from the centres of these clusters and cover over $5\times R_{200c}$, which allows us to investigate the environments around the clusters. Many studies, such as Haggar et al. (2020); Kuchner et al. (2020); Rost et al. (2021); Kuchner et al. (2021), benefit from this large re-simulation region. Since MDPL2 adopted the Planck cosmology parameters (Planck Collaboration et al., 2016), these re-simulations follow the same cosmology. These cluster regions have been simulated with a range of simulation prescriptions including Gadget-MUSIC (Sembolini et al., 2013), Gadget-X (Rasia et al., 2015; Beck et al., 2016) and Gizmo-Simba (Davé et al., 2019; Cui et al., 2022). In this study, we will focus on the Gizmo-Simba result and leave a cross-model comparison for later work.

The Gizmo-Simba run uses similar input physics to the Simba simulation (Davé et al., 2019), which is an advanced version of Mufasa (Davé et al., 2016) based on the Gizmo cosmological hydrodynamics code (Hopkins, 2015) with its meshless finite mass hydro-solver. Because of the different simulation resolutions ($m_{\mathrm{gas}}\approx 2\times 10^{7}$ M_⊙ for Simba vs $\approx 2\times 10^{8}$ M_⊙ for the300) and different objects (cosmological run for Simba and galaxy clusters for the300), the parameters for the baryon model are re-calibrated; see Cui et al. (2022) for detailed changes. Here, we only briefly list the baryon modelling included in Gizmo-Simba. Radiative cooling and photon-heating/ionization processes of gas are implemented using the Grackle-3.1 library (Smith et al., 2017) with a spatially-uniform ultraviolet background model (Haardt & Madau, 2012) and the self-shielding prescription following that in Rahmati et al. (2013), and an H₂-based star formation model is included from Mufasa (Davé et al., 2016). The decoupled two-phase model SN feedback is also adopted from Mufasa, but with the mass loading factor from Anglés-Alcázar et al. (2017b), while the chemical enrichment model tracks eleven elements with metals from supernovae Type Ia and Type II and Asymptotic Giant Branch stars. More importantly, the AGN model adopts two BH accretion descriptions – a torque-limited accretion model for the cold gas (Anglés-Alcázar et al., 2015, 2017a) and the Bondi accretion (Bondi, 1952) for hot gas – and three feedback modes: (1) ‘radiative mode’ or ‘quasar mode’ with bipolar outflows with velocities $\sim 500-1,500\,{\rm km}\,{\rm s}^{-1}$; (2) ‘jet mode’ with bipolar ejection up to $15,000\,{\rm km}\,{\rm s}^{-1}$; and (3) X-ray feedback mode following Choi et al. (2012). More details of the baryon models can be found in Davé et al. (2016); Davé et al. (2019).

The AGN jet mode in the Simba model is particularly effective at ejecting and heating gas to large radii from halo centres, more so than in comparable cosmological hydrodynamic simulations (for comparisons, see Christiansen et al., 2020; Yang et al., 2022; Ayromlou et al., 2022; Cui et al., 2022). It should be noted that the jet-mode feedback is implemented even more strongly in the300 run than the Simba cosmological volumes: jets are allowed to reach over twice-as-high velocities (maximally $15,000\,{\rm km}\,{\rm s}^{-1}$) than in Simba (maximally $7,000\,{\rm km}\,{\rm s}^{-1}$). This change was calibrated on the most massive clusters, necessary to achieve realistic quenching of cluster galaxies by $z=0$ (Cui et al., 2022) – however, the velocity formulae were applied identically in the lower-mass halos without a separate calibration. It is challenging to assess exactly how this affects the feedback with respect to the Simba cosmological volumes, as the jet speed varies rapidly but snapshots are only saved at coarsely-spaced time intervals. The overall picture is certainly that the300 Gizmo-Simba runs live on the extreme end of WHIM heating and gas ejection from lower-mass halos compared to other recent simulations like Illustris TNG. However, differences from other simulations does not signify inaccuracy. At the time of writing, there exists scant observational evidence for the gas properties of halos with $M<10^{14}$M_⊙, especially at $z=1$, meaning the simulations assessed in this paper are viable models.

These simulations are analysed with two different algorithms to create halo catalogues: Caesar³³3https://github.com/dnarayanan/caesar and The Amiga Halo Finder⁴⁴4http://popia.ft.uam.es/AHF/ (AHF, Gill et al., 2004; Knollmann & Knebe, 2009). The Caesar package identifies halos by a Friends-of-Friends (FOF) algorithm with linking length $b=0.2$. Such a linking length approximately leads to the identification of halos with an overdensity of 200 times the critical density (White, 2001), the same overdensity that is also frequently used to describe the virial radius (Navarro et al., 1997, e.g.,). The FOF halos are identified by firstly linking dark matter particles; then, all the other particles within the halo are identified by attaching each to the nearest dark matter particle. For Caesar FOF halos, the mass is given by the sum of all particles belonging to the halo. Galaxies in Caesar are identified through 6D FOF linking.

Meanwhile, AHF defines halos based on the spherical overdensity method. Each halo has a position given by the halo’s highest density peak, and a radius $R_{200c}$ within which the average density equals 200 times the critical density $\rho_{c}$. Thus, the halo mass is $M_{200c}$. Furthermore, the subhalos are identified through an unbinding process. However, we do not examine subhalos in this work. For larger halos hosting subhalos, the subhalo mass is included within the host halo mass.

As previously stated, we are mainly interested in results at $z=1$. The closest snapshot redshift is $z=0.987$; only simulation snapshots at that redshift are included in the main analysis and are referred to as $z=1$ for simplicity. A single brief comparison at $z=2$ is clearly marked.

The oriented stacking method detailed in L22, which we will follow here, first creates many small cutouts from a Compton-$y$ map. Each cutout is centred on the location of a galaxy cluster. The large-scale structure surrounding the galaxy cluster location is then examined via an external galaxy catalogue. Gaussian smoothing is applied to the projected galaxies in a thin redshift bin surrounding the cluster redshift, and the smoothing scale is chosen based on the structures of interest. The curvature of the smoothed galaxy map provides the necessary information to determine the axis of strongest filamentary structure around each cluster (detailed in Sec. 3.4). The $y$ map cutouts are each individually rotated such that this filament axis is aligned across the entire sample. Finally, the sample of cutouts is stacked. The stacked image has a high SNR where structure is overlaid – both in the centre, which displays the average galaxy cluster tSZ profile, and along the filament or supercluster axis, which averages over all contributions to the large-scale extended gas structure. The contributions include those from the intracluster medium in neighboring clusters, the intragroup medium in intermediate-mass halos, the circumgalactic medium around individual galaxies, and the WHIM.

Oriented stacks of $y$ for splits of halo and diffuse gas using $R_{200c}$, averaging over 98 snapshots, are shown in Fig. 5 in logarithmic scale. The $y$ contribution from the warm-hot and hot gas is also shown. The cool-gas $y$ map is not shown as it the stacked contribution from particles with $T<10^{5}$ K is $y\sim 2\times 10^{-9}$ at maximum, negligible in comparison to all other categories. The halo gas map shows a strong central signal from the stacked central cluster with $M>10^{14}$$h^{-1}{\rm{M_{\odot}}}$; this is also seen in the y map of only massive clusters. The diffuse gas map, made up of an entirely distinct set of particles as the halo gas map, also has a non-zero central signal surrounded by a ring of higher $y$. This morphology comes from averaging over shells of high-pressure gas just beyond $R_{200c}$ of the central stacked clusters. Furthermore, careful examination reveals many such shells of $y$ further afield in the diffuse map, whose locations correspond to halos in the massive and low-mass cluster maps. The association between halos and diffuse gas is further explored in Sec. 4.2.

The stronger signal along the horizontal filament axis is visible in each stack, however it appears more concentrated for the halo gas and hot gas maps than for the diffuse and warm-hot maps. This indicates that the diffuse gas is less compressed along the filament axis defined by the galaxies. The puffier appearance relates to the fact that energy injected by AGN sphericalises at large distances from the galaxies (Davé et al., 2019; Cui et al., 2022).

To quantify these observations, we also decompose each stack $I$ into multipole moments $m$:

where $r$ and $\theta$ are the polar coordinates. We will focus only on the first two even cosine amplitudes, $C_{0}$ and $C_{2}$. The $C_{0}$ monopole term quantifies the isotropic profile of the gas, equivalent to the circularly-averaged profile of an unoriented stack. The quadrupole $C_{2}$ quantifies most of the anisotropy of the gas, summing the signal from an opening angle around the $x$ axis. $C_{2}$ is useful especially when comparing observational oriented stacks to expectations from simulations, as it depends on not only the gas processes in the simulation but also the shapes and extents of large-scale structures like filaments, which hold cosmological information. Higher $C_{m}$ moments are cosmologically interesting (see L22), but have lower SNR. The key questions of this study can be addressed without these higher order moments. $S_{m}(r)$ quantifies the contribution along the vertical misaligned axis. It may provide information about underdense regions, but we leave this to future work.

The radial profile of the cosine components are taken by

where $X=2$ when $m=0$ and $X=1$ when $m=2$. To quantify the variance in $C_{m}(r)$ among the 98 snapshots, we perform the decomposition for each individual rotated $y$ map from each snapshot. The standard error on the mean (SEM) for each bin in each multipole profile is the standard deviation of the bin across the 98 independent maps divided by $\sqrt{98}$. In the resulting figures, error bars show $\pm$1 SEM. We examine results only out to 8$h^{-1}$Mpc from the central cluster to avoid contamination and projection effects in the snapshot edges.

In the following sections, the discussion focuses on addressing the key questions for modelling the anisotropic tSZ signal in filaments and superclusters. The first addresses extent–to what radius should a halo model go to capture a sufficient fraction of the $y$ signal?–and the second addresses mass range–is it more important to understand the pressure profiles for halos in a particular mass range than others?

Lastly, we sort the 98 cluster regions by measures of their density and elongation, as determined by the galaxy field, to assess whether these measures are effective at increasing the signal from diffuse gas. The density and elongation of the large-scale galaxy environment are indicators of the overall superclustering of matter in a given region, and are both good predictors of the anisotropic $y$ signal in observations (L22). At the position of the massive halos selected in Sec. 3.1, we examine the Hessian of the smoothed galaxy number density $\stackrel{{\scriptstyle\sim}}{{\smash{n_{g}}\rule{0.0pt}{4.73611pt}}}$. Using the eigenvalues of the Hessian defined in Eq. 2, the ellipticity $e$ is defined as:

We measure this large-scale ellipticity $e$ and the large-scale density $\stackrel{{\scriptstyle\sim}}{{\smash{n_{g}}\rule{0.0pt}{4.73611pt}}}$ at the position of the most massive cluster and combine them to make a normalized measure of superclustering, $s$:

where the means ($\langle...\rangle$) are taken over the sample of 98 snapshots. No further normalization or rescaling (e.g., dividing by the background field rms $\sigma$), is needed because the snapshots all take place at the same redshift.

It is expected that $s$ will enhance the filament signal from halo gas because of the strong correlation between galaxy number and gas pressure for massive halos. However, it is less obvious whether $s$ can also boost the signal from the most diffuse, unassociated gas. Figures 5 and 9 indicate that the warm-hot gas is very diffuse. As this gas is the object of missing baryon searches, we examine whether using $s$ to select regions of high superclustering is helpful for increasing the signal from the warm-hot gas.

We rank-order the snapshots by $s$, then stack the rotated $y$ maps from the top 20 snapshots and bottom 20 snapshots. We also test samples of more or less than 20 to check for robustness. The results are shown in Fig. 11. The effect of $s$ ranking has a clear impact in the all-gas signal, but the impact on the warm-hot gas is less clear. Fig. 12 quantifies the effects with a multipole decomposition. As expected, the $s$ measure has a higher effect on the total $y$ signal than it does on the warm-hot signal. Nevertheless, far beyond the central stacked cluster ($r\gtrsim 4$$h^{-1}$Mpc, the high-$s$ warm-hot stacks show a subtle boost in $m=0$ (the warm-hot gas is $\sim 10-30\%$ stronger) and a more distinct boost ($\sim 30-70\%$) in $m=2$. Interestingly, in the near-outskirts of the central stacked cluster ($r\sim 2$$h^{-1}$Mpc), the effect in $m=2$ is opposite. It is unclear why the low-$s$ stack exhibits stronger anisotropy close to the central cluster.

In summary, this brief investigation demonstrates the potential of using large-scale environmental measures from galaxy catalogues to select for regions with higher diffuse gas content. In future, it would be interesting to investigate whether measuring the large-scale environment with different populations of galaxies (e.g., field galaxies versus cluster galaxies, divided by color) changes the correlation with diffuse gas. Such a study would aid in the observational search for the remainder of the missing baryons.

In this work, we have applied oriented stacking to simulated maps of Compton-$y$ from the300 Gizmo-Simba  runs. We have applied several different cuts to the gas particles to separate the particles associated with halos from the diffuse or unresolved particles. We focused our discussion on the separation that uses a sharp cutoff radius $R_{200c}$ or 2$R_{200c}$, because this is most readily applied observationally, either through masking or through halo prescription modelling.

We find that the contribution to the stacked filament tSZ signal from beyond $R_{200c}$ of all halos is approximately equal in magnitude to the contribution from gas within halos. However, much of this signal comes from the shells of gas between $R_{200c}$ and 2$R_{200c}$, indicating that the gas is still associated with halos despite being diffuse by this basic definition. The gas beyond this radius for massive and low-mass clusters is frequently hotter than $10^{7}$K, placing it outside of the warm-hot definition often used in simulations ($10^{5}\mathrm{K}<T<10^{7}\mathrm{K}$). This heating is likely due to both shocks on infall and AGN heating of outgoing gas, since gas velocities in the simulation show both ingoing and outgoing motions.

Thus, precise modelling of the gas from 1 to 2$R_{200c}$ is very important in order to use extended, stacked tSZ signals for cosmology. This conclusion is consistent with the recent observational study of galaxy stacking at slightly lower redshifts by Schaan et al. (2021), which showed that a significant fraction of the baryons are beyond the virial radius. Meanwhile, Baxter et al. (2021) studied shocks in the300 simulations at $z\sim 0.2$, finding a shock feature in the $y$ signal of relaxed clusters at several virial radii. Further studies of such shocks in simulations at varying redshifts, and observational studies of the outskirts of halos with coming data, will be important for more accurate modelling of the diffuse tSZ signal over a range of redshifts.

In our study, going out to $2R_{200c}$ captures 75% of the signal; however, this is also insufficient for using anisotropic stacked tSZ signals for precision cosmology. Modelling the halo signal further out is a possible solution. For example, there appears to be some association of the gas pressure with halo locations even at $4R_{200c}$ and beyond. However, it is challenging to develop a prescription which extends so far due to the extensive overlap between neighboring profiles. In many overlapping cases, it is unclear which halo the gas should be assigned to. Furthermore, it is challenging to study such large halo extents in the300 simulations due to their limited volume. It would be useful for future work to examine the oriented tSZ signal beyond $4R_{200c}$ in larger-volume simulations, or attempt to separate the halo-associated gas with the field gas at all radii using non-spatial properties, to give insight into whether gas pressure modelling for the field is a necessary addition to halo modelling.

A caveat to this work is that in all definitions applied to separate halo from non-halo gas, halos smaller than $10^{12}$ $h^{-1}{\rm{M_{\odot}}}$ were treated as members of the diffuse category due to their low resolution. It is possible that in higher-resolution simulations, such halos (numerous as they are) would be associated with a significant fraction of the tSZ signal. An additional limitation is that our analysis was restricted to a single suite of simulations with one feedback prescription, and the results will likely vary depending on the feedback model. The small volume of the300 zoom simulations as compared to cosmological volumes make the300 well suited for re-runs with subtle changes to the feedback mechanisms, and future work will explore this.

We also demonstrated that environmental measures of superclustering, as determined by the galaxy field, augment the signals from the warm-hot diffuse gas. Such environmental measures can be applied to help detect this gas in existing data like that from ACT and Planck, but claimed detections of unbound or diffuse baryons should clearly state the modeling that went in to account for halo gas.

Several previous works have claimed a detection of the WHIM with oriented tSZ data. To compare our simulated results with those observations, it is important to again emphasize the role that environment plays in the diffuse gas pressure. Even before imposing the additional constraints for superclustering, the full sample of simulations examined in this work are already constrained to rare environments. They host a rare massive cluster at $z=0$; at $z=1$, we are viewing a snapshot of dynamic, coalescing proto-cluster regions that have higher density and more activity than average. The diffuse $y$ signals in similar regions in the real Universe are likely to be higher than the average intergalactic filament.

To compare with observations, we examine the $y$ profiles from the warm-hot stacks along the filament axis (total signal; roughly equal to the sum of $y_{m=0}+y_{m=2}$ in the red curves in Fig. 9). Various groups have studied the low-$z$ cluster pre-merger bridge Abell 399-401, most recently Hincks et al. (2022) in ACT data, finding a $y$ signal from the bridge of $10^{-5}$. With the bridge being 4.4 $h^{-1}$Mpc in projection, the centre is 2.2 $h^{-1}$Mpc from each cluster; the warm-hot $y$ signal in our work is nearly two orders of magnitude smaller at the same distance. This is likely both because the A399-401 system is in a rare state where the inter-cluster gas has been heated prior to its impending merger, and possibly also related to the difference in redshift. On the other hand, both de Graaff et al. (2019) and Tanimura et al. (2019) found a much smaller residual $y$ signal of $\sim 10^{-8}$ between stacked pairs of massive SDSS galaxies from $z\sim 0.4-0.75$, after subtracting the estimated contributions from the stacked halos, and the former group estimated that $y\sim(5-6)\times 10^{-9}$ comes from diffuse gas. Based on the pair separation in those studies, the signal examined was at a distance of $3-7$ $h^{-1}$Mpc from the halos; the300 warm-hot signal is $\sim 10\times$ larger at that same distance. This may also be due to environmental differences–while the pair condition selects for overdense regions, it does not guarantee them to be in such dense proto-cluster environments as examined in this study. However, it may also be a hint that the feedback in the300 is too powerful in heating diffuse gas.

The setup in this work is most similar to that in L22, which did oriented stacking on DES optically-selected red-sequence clusters with an average halo mass of roughly $5\times 10^{13}$ M_⊙. That work did not attempt to separate the halo gas from non-halo gas. The orientations in that work were determined on similar scales as the current work ($\sim 12$ $h^{-1}$Mpc compared to $10$ $h^{-1}$Mpc), and the corresponding aligned $y_{m=2}$ signal was found to be $\sim 2.7\times 10^{-8}$ at 7 $h^{-1}$Mpc from the cluster. The $m=2$ component of the all-gas $y$ signal found for The Three Hundred  is similar: only $\sim 3\times$ higher. Again, this difference may be due to environment. In L22 it was found that when the DES sample is limited to only those clusters in large scale high-superclustering regions, the signal grows to values very similar to those found in the present work. Some differences may also come from observational limitations: galaxy photometric redshifts, a larger redshift bin size, and galaxy magnitude limits all contribute to orientation inaccuracy in real data which can bias the signal lower compared to the idealized simulation result.

Given the challenges in exactly replicating the environmental selections in existing filament-gas observational studies with the300, it is not possible to comment decisively on whether these simulations produce diffuse gas in accurate amounts and thermal states. The range of values analyzed in this work fall within the ranges seen in observations thus far, and a qualitative description in which environment plays an important role in properties of the WHIM does well at encompassing the results. Future simulations encompassing a larger range of environments would be helpful to complement upcoming advances in observational data, so that anisotropic studies of the tSZ can be used to make quantitative constraints on feedback models.

Beyond the massive clusters, the level of $y$ signals in this work are not detectable in an individual tSZ image in any current data, or even future data, due to the noise levels. Appendix A provides some examples of what $y$ maps from the simulations would look like given forecasted noise for the upcoming Simons Observatory (SO, Ade et al., 2019) and CMB-Stage 4 (CMB-S4, Abazajian et al., 2022). Even with the improvements those surveys will bring, stacking will be necessary to increase the SNR sufficiently to measure the diffuse gas in all but the rarest systems. With current data, the diffuse gas is on the cusp of detectability with the methods employed in this work. At a level of $y\sim 10^{-7}$, the signal is of similar magnitude as the $1\sigma$ uncertainties in the oriented stacking analysis of ACT data (L22), which stacked $y$ map cutouts surrounding 5,500 DES clusters. To achieve a 3$\sigma$ measurement of a $y=10^{-7}$ signal with ACT data, 9$\times$ more regions—with signals similar to the300 snapshots assessed in this work—are needed. Despite recent large increases in the overlap between the ACT and DES footprints, enabling future use of the full DES Y3 cluster catalog ($\sim 10\times$ more clusters), most of these clusters are at lower masses than the300 clusters and thus their filamentary environments tend to produce smaller average tSZ signals. Therefore it is unclear whether the diffuse signal will be readily detectable with this method in incoming ACT data. However, with improvements to the methodology and/or use of upcoming SO data, the detection should be possible given accurate subtraction of halo gas signals, a major challenge to any missing baryons search.

With the coming advances in galaxy and CMB data, improvements in modelling gas both within and beyond halos are needed to fully realize the potential of the tSZ as a cosmological probe. The analysis presented in this work provides guidance to future observation and modelling efforts regarding the elusive boundless baryons.

The authors thank the anonymous referee for a thorough and helpful review of the manuscript.

This work has been made possible by the ‘The Three Hundred’ collaboration⁸⁸8https://www.the300-project.org, which has received financial support from the European Union’s Horizon 2020 Research and Innovation programme under the Marie Sklodowskaw-Curie grant agreement number 734374, i.e. the LACEGAL project. The simulations used in this paper have been performed in the MareNostrum Supercomputer at the Barcelona Supercomputing centre, thanks to CPU time granted by the Red Espa$\stackrel{{\scriptstyle\sim}}{{\smash{\rm n}\rule{0.0pt}{4.73611pt}}}$ola de Supercomputaci$\acute{\rm o}$n. The CosmoSim database used in this paper is a service by the Leibniz-Institute for Astrophysics Potsdam (AIP). The MultiDark database was developed in cooperation with the Spanish MultiDark Consolider Project CSD2009-00064.

ML acknowledges support from the Natural Sciences and Engineering Research Council of Canada Postgraduate Scholarships - Doctoral. ML also thanks Karen Scora for coming up with the alliterative title.

WC is supported by the STFC AGP Grant ST/V000594/1 and the Atracción de Talento Contract no. 2020-T1/TIC-19882 granted by the Comunidad de Madrid in Spain. He also thanks the Ministerio de Ciencia e Innovación (Spain) for financial support under Project grant PID2021-122603NB-C21. He further acknowledges the science research grants from the China Manned Space Project with NO. CMS-CSST-2021-A01 and CMS-CSST-2021-B01.

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

R. H. is a CIFAR Azrieli Global Scholar, Gravity & the Extreme Universe Program, 2019, and a 2020 Alfred P. Sloan Research Fellowship. RH is supported by Natural Sciences and Engineering Research Council of Canada Discovery Grant Program and the Connaught Fund.

The data underlying this work have been provided by The Three Hundred collaboration. The data may be shared on reasonable request to the corresponding author, with the permission of the collaboration.