How do baryonic effects on the cosmic matter distribution vary with scale and local density environment?

In this study, we investigate how the baryonic effects vary with scale and local density environment mainly by utilizing a novel statistic, the environment-dependent wavelet power spectrum (env-WPS). With four state-of-the-art cosmological simulation suites, EAGLE, SIMBA, Illustris, and IllustrisTNG, we compare the env-WPS of the total matter density field between the hydrodynamic and dark matter-only (DMO) runs at $z=0$. We find that the clustering is most strongly suppressed in the emptiest environment of $\rho_\mathrm{m}/\bar\rho_\mathrm{m}<0.1$ with maximum amplitudes $\sim67-89$ per cent on scales $\sim1.86-10.96\ h\mathrm{Mpc}^{-1}$, and less suppressed in higher density environments on small scales (except Illustris). In the environments of $\rho_\mathrm{m}/\bar\rho_\mathrm{m}\geqslant0.316$ ($\geqslant10$ in EAGLE), the feedbacks also lead to enhancement features at intermediate and large scales, which is most pronounced in the densest environment of $\rho_\mathrm{m}/\bar\rho_\mathrm{m}\geqslant100$ and reaches a maximum $\sim 7-15$ per cent on scales $\sim0.87-2.62\ h\mathrm{Mpc}^{-1}$ (except Illustris). The baryon fraction of the local environment decreases with increasing density, denoting the feedback strength, and potentially explaining some differences between simulations. We also measure the volume and mass fractions of local environments, which are affected by $\gtrsim 1$ per cent due to baryon physics. In conclusion, our results show that the baryonic processes can strongly modify the overall cosmic structure on the scales of $k>0.1\ h\mathrm{Mpc}^{-1}$, which encourages further research in this direction.


INTRODUCTION
The standard ΛCDM cosmological paradigm depicts the geometry and matter composition of the Universe and predicts how structures like galaxies, galaxy clusters, and the cosmic web formed from tiny density fluctuations, which are confirmed by extensive observations (e.g., see Frenk & White 2012, for a review).However, the collisionless cold dark matter-only (DMO) simulations indicate that the standard ΛCDM model faces some challenges at galactic and subgalactic scales, such as the core-cusp problem, too big to fail problem, missing satellites problem, and diversity problem (e.g.Klypin et al. 1999;de Blok 2010;Boylan-Kolchin et al. 2011;Primack 2012;Oman et al. 2015;Weinberg et al. 2015;Bullock & Boylan-Kolchin 2017).The possible ways to solve these small-scale problems can be summarized into three categories: (i) replace DMO N -body simulations with hydrodynamic simulations incorporating baryons and baryonic physics (e.g.Sales et al. 2022), (ii) propose alternatives to cold dark matter (e.g.Polisensky & Ricotti 2011), or (iii) modify the laws of gravity (e.g.Sanders & McGaugh 2002), in which the first solution holds the most promise, and alleviates the tension between theory and observation considerably (e.g.Chan et al. 2015;Simpson et al. 2018;Engler et al. 2021).
In order to faithfully mimic galaxy formation and evolution, ⋆ E-mail: yunw@jlu.edu.cn† E-mail: hep@jlu.edu.cn the cosmological hydrodynamic simulations solve the Euler equations to simulate baryons, and adopt fine-tuned sub-grid (or subresolution) models to implement baryonic processes including, e.g., radiative cooling, star formation, and feedbacks from stellar winds, supernovae and active galactic nucleus (AGN) (see Vogelsberger et al. 2020;Crain & van de Voort 2023, for review).Radiative cooling and star formation cause baryons to flow into the dark matter halo and collapse within the central region, whereas the feedbacks expel mass and energy of baryons from the halo.Moreover, the turbulent heating of the intergalactic medium may also be crucial in preventing the gravitational collapse of baryons (e.g.Zhu et al. 2010;Zhuravleva et al. 2014;Schmidt et al. 2016Schmidt et al. , 2017;;Yang et al. 2020Yang et al. , 2022)).The inclusion of baryons and baryonic physics is anticipated to affect the spatial distribution of matter, partly due to the redistribution of baryons by baryonic processes, and partly due to the redistribution of dark matter by the gravitational coupling of baryons and dark matter, which is dubbed "back-reaction" (van Daalen et al. 2011).Such baryonic effects on the cosmic matter distribution are still less well understood, which will lead to significant systematic uncertainties in the next-generation large-scale surveys (e.g.DESI 1 , EUCLID 2 , LSST 3 , and WFIRST 4 ), and therefore need to be care-fully studied to extract unbiased cosmological information from survey data.
Numerous studies have made concerted efforts to investigate baryonic effects at low redshifts.These studies employed a variety of simulations with different box sizes, resolutions, cosmological parameters, and sub-grid models.Despite this diversity, the majority of these studies yield qualitatively consistent findings.For instance, studies focusing on the halo shapes demonstrate that the presence of baryons leads to rounder halos (e.g.Schaller et al. 2015;Butsky et al. 2016;Henson et al. 2017;Chua et al. 2019Chua et al. , 2022;;Cataldi et al. 2021Cataldi et al. , 2023)), and results on the halo mass function suggest that feedbacks can reduce the halo mass (e.g.Cui et al. 2014;Castro et al. 2021;Beltz-Mohrmann & Berlind 2021;Sorini et al. 2022).In addition, baryonic processes greatly impact the matter power spectrum.Particularly, at redshift zero, the AGN feedback suppresses the power on scales of 0.1hMpc −1 ≲ k ≲ 20hMpc −1 with the suppression reaching a maximum of > 10 per cent at k ∼ 10hMpc −1 , while on scales of k ≳ 20hMpc −1 , the suppression begins to be less effective due to gas cooling (e.g.van Daalen et al. 2011van Daalen et al. , 2020;;Hellwing et al. 2016;Chisari et al. 2018;Springel et al. 2018;Debackere et al. 2020;Aricò et al. 2020).However, we notice that turbulent heating can also suppress the matter power spectrum and alleviate the over-cooling problem (Voit 2005), which is supported by many studies (e.g.Zhu et al. 2010;Zhuravleva et al. 2014;Schmidt et al. 2016Schmidt et al. , 2017;;Yang et al. 2020Yang et al. , 2022)), while not being sufficiently emphasized by the leading simulations.Furthermore, considering that power spectrum cannot carry the full information of the matter density field, some works investigated baryonic effects on the higher order statistic, i.e. the bispectrum, and found that the equilateral bispectrum is enhanced around k ∼ 2 − 3hMpc −1 due to baryons backreacting on the dark matter within massive haloes (Foreman et al. 2020;Takahashi et al. 2020;Aricò et al. 2021).It should be noted that most works mentioned above predict or explain the baryonic effects based on the strong assumption that baryons only alter the matter distribution within halos.By measuring the probability distribution function of the cosmic density field, however, the recent work of Sunseri et al. (2023) found that feedback processes can also have pronounced effects on voids, walls, and filaments, which highlights the importance of modeling baryonic effects in the whole cosmic structure.
Motivated by previous studies, we intend to delve deeper into how the impact of baryonic physics on the matter distribution varies with scale and local environment.To this end, we adopt the environmentdependent wavelet power spectrum (env-WPS) as our fiducial statistic, which was originally proposed by Wang & He (2022).As a bivariate function of scale and local density environment, the env-WPS is built upon the continuous wavelet transform (CWT) and tells how the clustering strength of matter is distributed over different environments and scales.In Wang & He (2022), we demonstrated that the env-WPS is capable of discriminating between the non-Gaussian density field and the homogeneous Gaussian random field, while the traditional Fourier power spectrum (FPS) is helpless in this regard.Accordingly, the env-WPS is a more advanced and powerful tool to characterize the large-scale structure of the Universe.
To ensure the reliability of our results, we utilize four modern cosmological-scale simulation projects, EAGLE5 , SIMBA6 , Illus-tris7 and IllustrisTNG8 .In this work, we focus on analyzing the total matter density field at z = 0, which is highly nonlinear and non-Gaussian.By measuring and comparing the env-WPS of the total matter between the hydrodynamic and the corresponding DMO runs, we will quantify the scale-and environment-dependence of the effects of baryonic physics on the matter spatial distribution.Therefore, we hope that the present study will shed new light on the modeling baryonic effects.
The remainder of this manuscript is organized as follows.In Section 2, we briefly describe the simulations used here.In Section3, we introduce the env-WPS and its computational procedure.We present and discuss the results in Section 4, followed by the summary and outlook in Section 5.
These simulations track the evolution of billions of dark matter particles and baryonic elements over large cosmological volumes with L box ≳ 100 Mpc from the high-redshift (z > 100) to the present (z = 0).To do so, the Illustris and IllustrisTNG are executed with the Voronoi moving-mesh code AREPO (Springel 2010), while the EAGLE and SIMBA are implemented using a modified version of the smoothed particle hydrodynamics code GADGET-3 (Springel 2005) and the meshless finite-mass hydrodynamics code GIZMO (Hopkins 2015), respectively.All the considered simulations assume a concordance ΛCDM cosmology with parameters fixed by WMAP or Planck satellites.They all take into account the radiative cooling of gas, the formation of stars from the cold gas, stellar evolution, black hole formation and accretion, stellar and AGN feedbacks, but with different implementations.For details of these simulations, the reader is referred to the literature (Crain et al. 2015;Schaye et al. 2015;Davé et al. 2019;Vogelsberger et al. 2013;Torrey et al. 2014;Pillepich et al. 2018b;Weinberger et al. 2018).In particular, AGN feedback plays a vital role in the simulations, which can regulate star formation, thereby alleviating the overcooling problem substantially (McCarthy et al. 2010(McCarthy et al. , 2011)).The Illustris simulation employs an over-violent radio-mode AGN feedback model, resulting in the gas fraction in group-sized haloes being lower than observed (Genel et al. 2014).Comparatively, IllustrisTNG, EAGLE, and SIMBA are more realistic simulations, yielding better consistency with the available observations of galaxy formation and evolution.
In this work, we use the highest resolution runs of simulations at z = 0, which are listed in Table 1.Throughout this paper, we refer to Illustris-1, TNG100-1, TNG300-1, RefL0100N1504, and m100n1024 as Illustris, TNG100, TNG300, EAGLE and SIMBA, correspondingly.Comparing results from them enables us to scrutinize the effects of baryonic physics on the matter distribution in a less biased manner.

The env-WPS of the density field
The total matter distribution of the Universe can be represented by the density contrast, defined as where ρm is the background density of the Universe.Considering that the total matter is composed of the dark matter and baryons, δm(x) can be expressed as where f dm = (Ωm − Ω b )/Ωm, and f b = Ω b /Ωm.By convolving the total matter density field δm with the wavelet Ψ, we get its CWT as in which Ψ(w, x) = w 0 Ψ(k)/kdk| < +∞, the density field can be recovered via Then we define the env-WPS as which is the statistical average of the squared modulus of the CWT over coordinates with the same local matter density at the scale of w.By averaging over all densities, the env-WPS will degenerate to the global wavelet power spectrum (global-WPS): where ⟨. ..⟩V denotes the statistical average over the full space.Obviously, the relationship between the global-WPS Pm(w) and the env-WPS Pm(w, δ) is where fV (δ) = V δ /V is the volume fraction of the density environment of δm(x) = δ.For more details of the wavelet analysis methods we developed, we refer the reader to Wang & He (2021, 2022, 2023).
To isolate the impact of baryonic processes on the matter distribution, we can compare the same statistic between the hydrodynamical simulations and their matched DMO runs with the same initial conditions.We primarily employ the env-WPS to accomplish the objective of this study, while treating others as supplementary statistics.

Numerical procedure
In this subsection, we present the numerical procedure for estimating the env-WPS from simulations in four steps:9

(i) Assign mass to the regular grid
To conduct our analysis, the initial step involves assigning the mass distribution of particles to the regular Cartesian grid.The mass assignment schemes usually employed include Cloud-in-Cell (CIC), Triangular-Shaped Cloud (TSC), Piecewise Cubic Spline (PCS), and the scale functions of Daubechies wavelets (see Cui et al. 2008;Sefusatti et al. 2016).However, despite the last scheme is well-suited for accurately measuring Fourier-based statistics, e.g.power spectrum and bispectrum, the fields in real space obtained using it suffer from spurious artifacts.Here, we assign mass to a grid with Ng = 15363 cells using the PCS assignment scheme, which performs better than the CIC and TSC in underdense regions.Readers are referred to appendix A for the tests of different schemes.(ii) Choose a specific wavelet function In the context of the CWT, there is a wide variety of candidate wavelets.The current study requires that we focus on both the real space and the scale domain, which suggests that a wavelet with a good trade-off between spatial and scale resolutions is needed.Thus, consistent with Wang & He (2022), we choose the isotropic cosine-weighted Gaussian-derived wavelet (CW-GDW) as the analyzing wavelet, the formula of which is  .162, 9) [9, 30.623) [30.623, 99) [99, +∞) , 3.162) [3.162, 10) [10, 31.623) [31.623, 100) [100, +∞) 2e/(9 + 55e) is the normalization constant such that ´|Ψ(x)| 2 d 3 x = 1, and its Fourier transform is where k = |k|.For this wavelet, the scale w can be expressed as the pseudo wavenumber kpseu by the following relation (see Meyers et al. 1993;Torrence & Compo 1998;Wang et al. 2022), with cw ≈ 0.3883.The radial view of the isotropic CW-GDW is shown in Fig. 1.As can be seen, it is well localized in both the real and Fourier domains.(iii) Perform the CWT by FFT The cosmic density fields of simulations are assumed to satisfy the periodic boundary conditions, i.e. δm(x) = δm(x + nL box ), where n = (nx, ny, nz) is an integer vector.Therefore, equation (3) can be rearranged as where V box = L 3 box is the volume of the simulation box, and Ψ P (w, x) is the periodized wavelet given by One can see that δm(w, x) is also periodic, which thus can be expressed in terms of Fourier series as follows where δm(k) and w − 3 2 Ψ(k/w) are Fourier transforms of the density field and the rescaled wavelet, respectively.We see that the CWT is just the inverse Fourier transform of the multiplication of δm(k) and w − 3 2 Ψ(k/w), which can be implemented by the Fast Fourier Transform (FFT) efficiently.Here the CWT of each density field is computed at N k = 25 equally logarithmically spaced scales in where k f = 2π/L box , and kNyq = N 1 3 g π/L box .In this scale range, the smearing and aliasing effects, as well as the shot noise are negligible (see appendix A). (iv) Measure the env-WPS In this step, we divide the densities into J = 8 bins: (a) six logarithmic bins equally divided between ρm/ρm = 0.1 and ρm/ρm = 100, (b) one with ρm/ρm < 0.1, and (c) one with ρm/ρm ⩾ 100, which are listed in Table 2.For the j-th environment δj, we compute the env-WPS with equation (5) as where N δ j is the number of cells within the environment δj, and thus the volume fraction is given by fV(δj) = N δ j /Ng.For the rest of the equations and plots, we will use the pseudo wavenumber kpseu as a proxy for the wavelet scale w, and throw away its subscript "pseu" without any ambiguity.We refer to Appendix B for a better understanding of the env-WPS.

BARYONIC EFFECTS
Here, we explore the baryonic effects on the total matter distribution by mainly examining the relative difference between the total matter env-WPS in hydrodynamic simulations and their corresponding DMO runs, which is given explicitly below10 10 −1 10 0 10 1 Illustris SIMBA EAGLE TNG100 TNG300 Wavelet Fourier Then for the total matter global-WPS, with equation ( 7), we have where rV(δ) + 1 = fV(δ)/f DMO V (δ), and QDMO(k, δ) = f DMO V (δ) PDMO(k, δ)/ PDMO(k) with δ QDMO(k, δ) = 1.In the above equation, the second and third equalities tell us the contribution to the relative global-WPS Pm(k)/ PDMO(k) from individual density environment.
Note that due to the limited volume of the simulation box and the fact that there is only one single realization of initial conditions for each simulation, the estimated global-and env-WPSs, are inevitably affected by cosmic variance.However, since the hydrodynamical run shares the same initial condition with its matched DMO run, some of the variance will cancel out in the relative differences, Rm(k) and Rm(k, δ).To estimate the residual cosmic variance, we divide the simulation volume into eight equal-sized sub-volumes and measure the relative differences in each sub-volume separately, as is similar to the approaches in Chisari et al. (2018) and Foreman et al. (2020).Furthermore, we quantify the spread of them by the 1-σ dispersion around the median values (see Appendix D).

Effects on the FPS and global-WPS
The effects of baryonic physics on the FPS have been extensively investigated and discussed (e.g.van Daalen et al. 2011van Daalen et al. , 2020;;Hellwing et al. 2016;Chisari et al. 2018;Springel et al. 2018;Debackere et al. 2020;Aricò et al. 2020).The consensus of these investigations is that the AGN feedback can suppress the FPS relative to the DMO simulations out to large scales of k < 1 hMpc −1 at z = 0.It is necessary to test the validity of our global-WPS measurements against those of the FPS, the results of which are shown in Fig. 2. For each simulation, the global-WPS is suppressed in the same fashion as the The "Suppression" row lists the maximum suppression and the scale where it occurs in the emptiest environment, δ0, for each simulation, while the "Enhancement" row shows the greatest enhancement in the densest environment, δ7.
FPS, albeit with subtle variations.This is expected since the global-WPS is a wavelet-weighted average of the FPS over the wavenumber space (see equation (B11) in Wang et al. 2022).Given this consistency, we do not discuss the baryonic effects on the global-WPS in detail here.

Effects on the env-WPS
In this subsection, we will consider how the baryonic effects on the total matter distribution may vary with scale and local density environment.Fig. 3 presents the relative difference of the total matter env-WPS in hydrodynamic runs with respect to that in DMO runs.By visual inspection, we can see that the cosmic variance has a sizable effect on large scales.A closer look shows that this effect is the minimal in TNG300, which has the largest cosmic volume among the simulations used here.However, for extremely overdense environments, the cosmic variance appears to be more severe in Illustris and SIMBA than in other simulations.Particularly, the relative differences in the environments of δ j⩾6 measured from the full volume of SIMBA, fall outside the 1-σ uncertainty of those from eight subvolumes, on scales of k ≲ 1 hMpc −1 .This phenomenon reflects the inhomogeneity of the AGNs in SIMBA, as shown in Fig. D1.
Clearly, the relative difference Rm(k, δ) displays a prominent scale-and environment-dependence for total matter density field.Regarding the general features, it can be seen that the env-WPS is significantly suppressed on small scales for all density environments, while is enhanced at intermediate and large scales for some environments.To be specific, the scale (labeled k th for convenience) between the suppression and enhancement regimes is equal to 2 hMpc −1 in Illustris, 3.3 hMpc −1 in SIMBA, 8 hMpc −1 in EAGLE, 7 hMpc −1 in TNG100, and 5 hMpc −1 in TNG300, respectively.On small scales of k > k th , the env-WPS is suppressed more strongly in lower density environments, with the exception of Illustris, in which the suppression is not monotonic with the local density.We also note that in the overdense regions of δ j⩾6 , the env-WPS is enhanced by 1 − 10 per cent on the galactic scales of k > 10 hMpc −1 for Illustris.This phenomenon may be related to more efficient star formation in Illustris (Nelson et al. 2015;Beltz-Mohrmann & Berlind 2021).Nevertheless, all simulations agree that the most significant suppression happens in the emptiest environment δ0, over almost all scales, which becomes greater than 1 per cent on k ∼ 0.1 hMpc −1 , and then reaches a maximum of ∼ 67−89 per cent on scales of 1.86 ≲ k ≲ 10.96 hMpc −1 (see Table 3 for exact values).This result essentially confirms the finding of Sunseri et al. (2023), which revealed that the probability distribution function of the density field is suppressed by nearly 100 per cent in the emptiest regions.On the scales of k < k th , we see that the env-WPS is enhanced in overdense environments of δ j⩾5 (for EAGLE), . The relative difference between the total matter env-WPS in hydrodynamic simulations and that in DMO runs at z = 0, measured from Illustris, SIMBA, EAGLE, TNG100, and TNG300, respectively.As labeled by the colorbar, the colors range from dark blue to dark red, corresponding to local densities from the lowest to the highest.The light-colored areas represent the cosmic variances of Rm(k, δ), quantified by the 1-σ dispersion of the measurements in sub-volumes.The matter clustering is suppressed in all environments when the scales are smaller than that indicated by the vertical dotted line, which is equal to 2 hMpc −1 in Illustris, 3.3 hMpc −1 in SIMBA, 8 hMpc −1 in EAGLE, 7 hMpc −1 in TNG100, and 5 hMpc −1 in TNG300.as well as in several lower density environments of δ j⩾2 .Among these environments, the enhancement is the most significant in the densest region at scales of k ≳ 1 hMpc −1 (except Illustris), and which peaks at scales of 0.87 − 2.62 hMpc −1 with amplitudes of 7 − 15 per cent (see Table 3).This result is more explicit for the EA-GLE, TNG100, and TNG300 simulations, in which measurements of Rm(k, δ7) are less affected by cosmic variance.
According to van Daalen et al. ( 2020), the FPS suppression due to baryonic physics can be predicted by the baryon fraction of haloes, which is indicative of the feedback strength.In our context, it is fascinating to investigate the baryon fraction of local density environments: in which ∆V cell = V box /Ng is the cell volume of the simulation box.For comparison between simulations with different cosmolo- 17), which demonstrates the absolute contribution of each density environment to the impact of baryons on the global-WPS.In these panels, the gray shaded regions denote the scale range where the global-WPS suppression is less than 1 per cent (see Fig. 2).
gies, we re-normalize the fraction f b (δj) as the measurements of which are displayed in Fig. 4. It is impossible here to tightly constrain the correlation between the baryon fraction and baryonic effects on the env-WPS, since we only have five simulations.Even so, these results can offer us some guidance.We find that the baryon fraction is above the mean fraction Ω b /Ωm at the low-density end, and then decreases to below it with increasing density.This finding is attributed to the gas ejection from the very overdense regions driven by feedback processes, in which the AGN feedback plays a vital role at the redshift of z = 0 (e.g.Weinberger et al. 2018;Sorini et al. 2022).Nonetheless, we also notice that Yang et al. ( 2020) reported a consistent trend by examining the scatter plot of ρ b /ρm vs. ρ dm (see Fig. 7 therein), in which they used the WIGEON simulation (Feng et al. 2004) without incorporating star formation and any feedback.Due to the excellent performance of WIGEON in capturing shock wave and complex structures, they concluded that the gas falling into the potential well of haloes is prevented by turbulence generated from the non-linear structure formation.Turbulent heating could be a crucial mechanism that suppresses the matter power spectrum and mitigates the over-cooling problem, thereby warranting further investigation.For Illustris and SIMBA, the baryon fraction is lower in overdense environments of δ j⩾4 , whereas higher in underdense environments of δ j⩽2 than that for other simulations.This means that AGN feed-back is more efficient in Illustris and SIMBA, which can explain the heavy global-WPS suppression in these two simulations (see Fig. 2).Furthermore, by referring to Table 3 and Figs. 3 and 4, it can be observed that in the emptiest environment, a higher baryon fraction is generally linked to a more pronounced suppression.A possible reason for this could be that the more effective AGN feedback injects more gas into the underdense regions, thereby raising gas pressure and impeding the gravitational clustering significantly.On the other hand, we see that for Illustris and SIMBA with lower baryon fractions in overdense environments of δ j⩾4 , the enhancement peaks in the densest environment are shifted to larger scales, and the enhancement amplitudes in the lower densities are more significant on scales of k ≲ 1 hMpc −1 , compared to other simulations.Noteworthy, previous works show that the equilateral bispectrum is also enhanced on scales of k ∼ 2 − 3 hMpc −1 in the hydrodynamic runs relative to that in the DMO runs (except Illustris), which is attributed to the gas reaccretion into massive haloes caused by the decreased AGN feedback strength at late-time (see Foreman et al. 2020;Takahashi et al. 2020;Aricò et al. 2021).The env-WPS enhancement in the environment of δ7 probably comes from the same source, as this environment largely overlaps with massive haloes.

Contributions from different environments to the global-WPS ratio
By examining the summand fV(δ) Pm(k, δ)/ PDMO(k) in equation (17), we can learn about the absolute contributions of indi- Illustris SIMBA EAGLE TNG100 TNG300 vidual environments to the global-WPS ratio Rm(k) at different scales, which are illustrated in Fig. 5. Interestingly, the variation of fV(δ) Pm(k, δ)/ PDMO(k) with scale and local density looks like a "bow-tie", which holds for all smiluations.The knots of the bow-tie are located at scales of k ∼ 0.7, 0.9, 0.9, and 0.75 hMpc −1 in the SIMBA, EAGLE, TNG100, and TNG300 simulations, respectively.However, in Illustris, the knot is much more diffuse, which covers the scales of 0.5 ≲ k ≲ 1 hMpc −1 .At the bow-tie knot, all environments contribute equally to the global-WPS ratio.On the right side of the knot, the contribution comes mainly from the densest environment and decreases with decreasing density.This result is consistent with the previous findings, which point out that the FPS and bispectrum ratios can be interpreted well by the one-halo term on scales of k ≳ 2 hMpc −1 (e.g.Foreman et al. 2020).Conversely, the densest environment contributes the least to the global-WPS ratio on the left side of the knot.For the Illustris and SIMBA simulations with more efficient AGN feedback, we note that on scales of ∼ 0.1 hMpc −1 , the few per cent suppression on the global-WPS is primarily due to underdense environments.

Volume and mass fractions of local density environments
When we consider the relationship between the env-WPS and global-WPS, the volume fraction of the local density environment fV(δ) is an ingredient that cannot be neglected.In Fig. 6, we compare the volume fractions between the hydrodynamic and DMO runs.From the top panel of this figure, we see that the environment volume fractions in the different hydrodynamic simulations are approximately identical to each other, and they all decrease with increasing density.Compared to the DMO case, all the volumes of the environments are affected by ≳ 1 per cent, as shown by the bottom panel.Specifically, volumes of the extreme environments, δ0, δ6, and δ7, are reduced, while intermediate environments of δ 1⩽j⩽5 are expanded.
For completeness, we also measure the mass fractions of environments by fM(δj) = δm(x)∈δ j [δm(x) + 1]/Ng, the results of which are shown in Fig. 7.It can be seen that in hydrodynamic simulations, the environment mass fraction grows as the local density increases.Obviously, the densest environment, which occupies only ∼ 0.1 percent of the volume, contains ∼ 40 per cent of the mass.In contrast, with ∼ 50 per cent by volume, the emptiest environment makes a mere ∼ 2 per cent of the mass of the Universe.By examining the relative difference of the environment mass fraction with respect to the DMO runs, rM(δj) = fM(δj)/f DMO M (δj) − 1, the bottom panel of Fig. 7 shows that overdense environments of δ6 and δ7 lose ∼ 1 − 10 per cent of their mass, while lower density environments gain a substantial mass, which is caused by the AGN feedback transporting large amounts of gas from inner halo to regions outside the halo.Cui et al. (2018) and Sunseri et al. (2023) also made comparisons of the mass and volume fractions of cosmic structures between the hydrodynamic and DMO runs.Using the TNG300 simulation, the latter found that haloes lose ∼ 10 percent mass to filaments, walls, and voids, while the structures change less in volume.This is in rough agreement with our results for TNG300.Surprisingly, the former supports that baryonic processes have almost no impact on large-scale structures.The sources of discrepancies may be twofold.First, the spatial and mass resolutions of the simulations used by Cui et al. (2018) are lower than ours.Second, differences in the cosmic web classification schemes also induce differences in the results.Specifically, Cui et al. ( 2018) used a velocity-shear-tensor code (Hoffman et al. 2012) and a tidal-tensor code (Forero-Romero et al. 2009), Sunseri et al. (2023) used the NEXUS algorithm (Cautun et al. 2013), and we divide the density field into different environments according to its face values.
Finally, taking Figs.5-7 together, we find that the summand fV(δ) Pm(k, δ)/ PDMO(k) is positively correlated with fV(δ) and anti-correlated with fM(δ) on scales larger than the bow-tie knot.Conversely, the situation is reversed on scales smaller than the bowtie knot.This indicates that the global-WPS suppression is mainly contributed by the large-volume structures on large scales, but dominated by the massive objects on small scales.

SUMMARY AND OUTLOOK
The purpose of this study is to understand how the baryonic effects on the cosmic matter distribution vary with scale and local density environment.For this, we primarily use the env-WPS, defined in equation ( 5), to interpret the present-day density fields of the total matter, which come from the SIMBA, EAGLE, Illustris, TNG100, and TNG300 simulations with volumes of ≳ (100 Mpc) 3 .Through a comparative analysis of the calculations between the hydrodynamic and DMO runs, we obtain some meaningful findings as follows: (i) The matter clustering is suppressed more severely in lower density environments on small scales (except Illustris).In particular, the suppression is most severe in the emptiest environment of ρm/ρm < 0.1 over all scales, with a maximum of 67 − 89 per cent on 1.86 ≲ k ≲ 10.96 hMpc −1 .(ii) On intermediate and large scales, the clustering is enhanced in the environments of ρm/ρm ⩾ 0.316 (ρm/ρm ⩾ 10 in EA-GLE).The most pronounced enhancement occurs in the densest environment of ρm/ρm ⩾ 100, which peaks on scales of 0.87 ≲ k ≲ 2.62 hMpc −1 with levels of 7 − 15 per cent (except Illustris).(iii) The baryon fraction of the local density environment decreases with increasing density, due to the feedback processes ejecting significant amounts of gas from the highly overdense regions.
Comparisons between different simulations indicate that the AGN feedback is most effective in Illustris as well as SIMBA, which could explain the main differences in the impact of baryons between them and other simulations.(iv) The suppression on the global-WPS is almost the same as on the Fourier power spectrum.The contributions of local density environments to the global-WPS suppression exhibit a "bow-tie" pattern, where its knot is located around 1 hMpc −1 .On scales larger than this, the suppression mainly comes from the large-volume underdense regions, whereas from the massive objects on smaller scales.(v) The mass and volume fractions of the local density environment are also affected by ≳ 1 per cent due to baryon physics.
In detail, the moderate-density environments of 0.1 ⩽ ρm/ρm < 31.623 are expanded, while the extreme density environments of ρm/ρm < 0.1 and ρm/ρm ⩾ 31.623 are contracted.There is a total of ∼ 10 per cent mass transferred from the highly overdense environment of ρm/ρm ⩾ 31.623 to lower density environments.
In summary, these findings support that the baryonic effects are substantial across all density environments, not limited to just the extremely overdense regions (haloes).Therefore, this study provides insights for further understanding and modeling the effects of bary-onic physics.However, several issues still remain to be clarified and require further investigation, which are: (i) The env-WPS is affected by the cosmic variance visibly on large scales and in extremely overdense environments, which suggests that the simulations with much larger volumes should be used.For example, the latest hydrodynamic simulations, MillenniumTNG (Pakmor et al. 2023) and FLAMINGO (Schaye et al. 2023), are proper choices, the former with a volume of ∼ (740 Mpc) 3 and the latter with volumes of ∼ (1 − 2.8 Gpc) 3 .(ii) We just qualitatively discuss the possible physical reasons underlying the results.For one thing, a more quantitative explanation requires us to employ data from multiple epochs, as well as to examine the redistribution of the baryons and the back-reaction on the dark matter.

APPENDIX A: TESTS OF DIFFERENT MASS ASSIGNMENT SCHEMES
To minimize the smearing and aliasing effects on the FPS or other Fourier-based statistics, Cui et al. (2008) found that the scale functions of Daubechies wavelet transformations are the optimal mass assignment window functions.However, does this scheme yield a physical density field in real space?The answer is no, as shown in Fig. A1, in which we examine the probability distribution functions of the total matter density fields obtained by using three types of Daubechies scale functions (db6, db12, and db20, see Daubechies 1992;Hand et al. 2018).It can be seen that Daubechies windows lead to too many unphysical values of ρm/ρm < 0.
As is well known, we can get the physical density field with ρm/ρm ⩾ 0 by using traditional mass assignment window functions, e.g. the CIC (2nd-order), TSC (3rd-order), and PCS (4thorder) (Sefusatti et al. 2016).Then we test different choices of window functions and investigate how the relative differences Rm(k) and Rm(k, δ) vary between them, the results of which are presented in Fig. A2.Most of the variation is seen passing from CIC to TSC, whereas the results converge for higher order schemes.We also see that the effects of windows on Rm(k) are marginal, while the effects on Rm(k, δ) are a bit more noticeable in the extremely underdense environments on small scales.Consequently, we use the PCS scheme to compute the density field, and consider only the measurements on scales of k ⩽ 0.4kNyq to avoid numerical effects, e.g.smearing, aliasing, and shot noise.

APPENDIX D: THE RESIDUAL COSMIC VARIANCE
Here we describe in more detail how the residual cosmic variance effects are estimated.First, we perform the CWT of the density field at N k scales.Then we divide the CWT field into N sub = 8 equal-sized subfields at each scale, and calculate the relative differences within each subfield, i.e.Rm,i (k, δ) = Pm,i (k, δ) / PDMO,i (k, δ) − 1 and Rm,i (k) = Pm,i (k) / PDMO,i (k) − 1, where "i" denotes the i-th subfield.Finally, if we use Rm,med to denote the median of relative differences between the different subfields, then the cosmic variances are indicated by the 1-σ (68 percentile) dispersion around the median, Rm,med, as shown in Figs. 2 and 3. We find that all relative differences measured from the total simulation volume are within the 1-σ uncertainty, except for those at k ≲ 1 hMpc −1 and in δ j⩾6 for SIMBA.Given this scale and environment, this result may reflect the inhomogeneous distribution of AGNs in SIMBA.To check this, we count the number fraction of black hole particles in each subvolume, the results of which are shown in Fig. D1.Compared to the other simulations, it can be seen that the dispersion of the number fractions across sub-volumes is much larger in SIMBA, demonstrating the inhomogeneity of the AGNs.

Figure 1 .
Figure1.The isotropic CW-GDW (top) and its Fourier transform (bottom), at the scales of w = cw, 2cw, and 3cw.For presentational convenience, all variables are dimensionless.

Figure 2 .
Figure2.A comparison of baryonic effects on the global-WPSs (solid lines) and that on the FPS (dashed lines) at z = 0, measured from different simulations.The color bands show the cosmic variances of the relative difference Rm(k), estimated by the 1-σ dispersion of the measurements in subvolumes.

Figure 4 .
Figure4.The re-normalized baryon fraction fb (δ) as a function of the local density at z = 0, measured from different simulations as labeled.We see that the baryon fraction is above the mean fraction Ω b /Ωm at the low-density end, while below the mean at the high-density end.

Figure 6 .
Figure 6.Top panel: the volume fraction occupied by each local density environment in hydrodynamic simulations.Bottom panel: the relative difference of the environment volume fraction between the hydrodynamic and DMO runs.

Figure 7 .
Figure 7. Top panel: the mass fraction of each local density environment in hydrodynamic simulations.Bottom panel: the relative difference of the environment mass fraction between the hydrodynamic and DMO runs.

Figure D1 .
Figure D1.The number fraction of black hole particles contained in the subvolume, f BH = N BH,sub /N BH , where N BH,sub and N BH are the number of black hole particles in the sub-volume and in the full volume, respectively.

Table 1 .
Simulation data employed in this study.From left to right, we present the simulation name, code used to conduct simulations, choice of cosmological parameters, redshift range, box size, run name, number of dark matter particles (N dm ), initial number of gas cells (Ngas ), dark matter mass resolution (m dm ), and baryon mass resolution (mgas).

Table 2 .
The local density environments.

Table 3 .
The maximum effects of baryonic physics observed in the extreme density environments.
For another, different baryonic processes impact the matter distribution in different ways, which can be quantified by adjusting or disabling a certain process in the simulation, such as the star formation, black hole accretion, stellar feedback or AGN feedback.(iii) Predicting the baryonic effects on the env-WPS from the baryon fraction of the local density environment would require a large number of simulations.It is possible to do this with the CAMELS project (Villaescusa-Navarro et al. 2021), which is a suite of 4233 cosmological simulations with varying cosmological and feedback parameters.(iv) In this work, we classify the local environments only by the face value of the density field.A comparison between different classification schemes (e.g.Hahn et al. 2007; Forero-Romero et al. 2009; Hoffman et al. 2012; Cautun et al. 2013) is necessary in order to obtain more robust conclusions.