Continuous Wavelet Analysis of Matter Clustering Using the Gaussian-derived Wavelet

The standard cosmological model (Lambda cold dark matter) states that the hierarchical matter clustering evolved from the random primordial fluctuations via gravitational instability, as revealed by the observed pattern of spatial distribution for luminous objects (Tegmark et al. 2004). However, luminous objects are made up of baryons, whose clustering is not completely consistent with that of the dark matter (DM) component because of a series of physical processes such as radiative cooling, star formation, and feedback processes. With hybrid N-body/hydrodynamic simulations, the baryonic effects on the large-scale clustering of matter have been extensively studied in Fourier space, and these studies show that the deviation between baryons and DM in mass distribution depends on the scale and redshift (e.g., Chisari et al. 2019; van Daalen et al. 2020; Yang et al. 2020). Moreover, many studies indicate that the clustering of both DM and baryonic matter is environmentally dependent (e.g., Abbas & Sheth 2005; Peng et al. 2010; Wang et al. 2018; Man et al. 2019). However, usual statistical schemes based on the ordinary Fourier transform, e.g., power spectrum and bispectrum, cannot be used to measure the environmental and scale dependence of the matter clustering at the same time, since these schemes do not contain any positional information on the signal. A straightforward solution to overcome this problem is to use the technique of the windowed Fourier transform (WFT), which takes harmonic waves multiplied by a window of fixed width as basis functions, hence enabling us to gain the positional and scale information simultaneously (Gabor 1946).

However, the fixed window size of the WFT leads to a poor space-scale resolution (e.g., Kaiser & Hudgins 1994; Romeo et al. 2004; Gao & Yan 2011). If we set the window size to be small enough, then the spatial resolution certainly matches the signal, but the scale/frequency localization is too poor to resolve the large scales/low frequencies. If we choose a wider window in order to have a finer scale/frequency resolution, the spatial resolution becomes too coarse to detect small-scale structures. Therefore, the WFT is not a suitable tool for analyzing signals with a wide range of scales, such as spatial distributions of the DM and galaxies.

An alternative way to achieve simultaneous analysis of space and scale is the wavelet transform (WT), which provides an adaptive space-scale resolution (e.g., Daubechies 1992; Kaiser & Hudgins 1994; Chui 1997; Torrence & Compo 1998; Van den Berg 2004; Mallat 2009; Addison 2017). WT analysis decomposes a signal into separate scale components using a set of scaled and shifted wavelets, which are well localized in both real and Fourier space. Consequently, the local features of the signal at different scales are revealed by this decomposition. WT methods can be classified into the discrete wavelet transform (DWT) and the continuous wavelet transform (CWT). The DWT, using orthogonal wavelet bases, operates over scales and positions based on the integer power of 2, hence giving the most compact representation of the signal. This leads to the effectiveness and ease of implementation of the DWT, which hence is particularly useful for information on compression (e.g., Khalifa et al. 2008; Abdulazeez et al. 2020). Due to its advantages, the DWT has also been applied to study the large-scale structure of the universe (e.g., Fujiwara & Soda 1996; Pando & Fang 1996; Pando et al. 1998, 2004; Fang & Feng 2000; Romeo et al. 2004; Liu & Fang 2008; Lu et al. 2010). However, there are mainly two drawbacks in the DWT caused by dyadic scales. First, the DWT provides a poor scale resolution so that some meaningful features of a signal cannot be detected. Second, it lacks translational invariance. The so-called translational invariance means that if a signal is translated, then its wavelet coefficients are translated by the same amount without other modifications at every scale (Addison 2017). Obviously, this is not the case for the DWT. A small translation of the signal can make the discrete wavelet coefficients vary substantially on different scales, thereby the total energy in the wavelet domain is not conserved after the signal was shifted. These drawbacks suggest that the DWT is not suitable for analyzing signals with extremely high complexity.

In contrast to the DWT, the CWT allows the scale and translation parameters to continuously change, which makes it translationally invariant and redundant. The redundancy guarantees that the CWT can provide high-resolution results, which are much easier to interpret than those obtained with the DWT (Aguiar-Conraria & Soares 2014; Addison 2018). As a result, the CWT is becoming more and more popular across different disciplines, including geophysics, biomedicine, economics, astrophysics, fluid mechanics, and so on. For instance in astrophysics and cosmology, the CWT is used for the detection of 21 cm signals (e.g., Gu et al. 2013), the detection of baryonic acoustic oscillations (e.g., Tian et al. 2011; Arnalte-Mur et al. 2012; Labatie et al. 2012), the detection of substructures in two-dimensional (2D) mass maps (e.g., Flin & Krywult 2006; Schwinn et al. 2018), the analysis of turbulence evolution in the intracluster medium (e.g., Shi et al. 2018; Roh et al. 2019), the correlation analysis of galactic images (e.g., Frick et al. 2001; Tabatabaei et al. 2013), the analysis of the multifractal character of the galaxy distribution (e.g., Martínez et al. 1993; Rozgacheva et al. 2012), and the analysis of cosmic microwave background maps (e.g., Cayón et al. 2001; Starck et al. 2004; González-Nuevo et al. 2006; Curto et al. 2011).

The main problem faced by the CWT, e.g., when analyzing one-dimensional (1D) signals, is that its classical inverse formula is a double integration, which results in a heavy computational effort in recovering the original signal. Although there is an alternative inverse transform formula in the form of a single integral for the complex-valued wavelets, the real wavelets have generally been considered to have no such simple inverse transformation (Delprat et al. 1992; Aguiar-Conraria & Soares 2014). In our previous work (Wang & He 2021), we proposed a novel scheme of constructing continuous wavelets to overcome this problem, in which the wavelet functions are obtained by taking the first derivative of smoothing functions with regard to the positively defined scale parameter. With this scheme, the original signal is recovered easily by integrating the wavelet coefficients with respect to the scale parameter. As an inspired example, we took the Gaussian function as a smoothing function to derive the wavelet dubbed the Gaussian-derived wavelet (hereafter GDW) and briefly discussed its preliminary application to the matter power spectrum, demonstrating the success of this scheme.

In this work, we use local statistics established on wavelet coefficients to further investigate the potential of the GDW for analyzing the large-scale clustering of matter. Specifically, the wavelet power spectrum (WPS), the wavelet cross correlation (WCC), and the wavelet bicoherence (WBC) are employed here. The WPS and WCC were first introduced by Hudgins et al. (1993) to examine atmospheric turbulence. The WPS measures the variance of a signal at various scales within some local region, and the WCC is used to quantify similarities between two signals. The WBC was originally proposed by van Milligen et al. (1995a, 1995b) to detect the short-lived structures induced by phase coupling in turbulence. Moreover, the cosmic baryonic density distribution at late times is similar to a fully developed turbulence (Shandarin & Zel'dovich 1989; He et al. 2006), which enlightens us to apply these tools to the context of the structure formation of the universe. To illustrate our approach in a more intuitive way, we shall use the 1D matter distributions, since a great deal of work on the large-scale structure of the universe is accomplished using 1D cosmology models (Gouda & Nakamura 1989; Fujiwara & Soda 1996; Tatekawa & Maeda 2001; Miller & Rouet 2010; Manfredi et al. 2016). The Zel'dovich approximation is a simple model that provides a good approximated solution for the nonlinear evolution of collisionless matter (Zel'dovich1970). In the 1D case, the Zel'dovich approximation is proved to be an exact solution in the fully nonlinear regime until the first singularity appears (Soda & Suto 1992). It is straightforward and efficient to calculate the 1D Zel'dovich solution, and so we discuss our analysis method with it. First, we decompose the density fields obtained from the 1D Zel'dovich approximation into wavelet components at different positions and scales. Then, by measuring the WPS and WBC of matter densities based on these components, we investigate the effects of the density environment on the matter clustering.

To further demonstrate the ability of wavelet statistics in characterizing the environmental dependence of the matter clustering, we generalize them to the three-dimensional (3D) versions and then apply them to the 3D density fields for the different matter components in the IllustrisTNG simulation (Springel et al. 2018; Naiman et al. 2018; Nelson et al. 2018; Pillepich et al. 2018; Marinacci et al. 2018). By measuring the WPS and WCC, we can investigate (1) the variation of the clustering strength with the density environment at given scales, (2) the variation of the cross correlations between different matter components with the density environment at given scales, and (3) the baryonic effects on the matter clustering in different density environments.

This paper is organized as follows. In Section 2, we describe the characteristics of the GDW, the theory of the CWT, and wavelet-based statistical tools. In Section 3, we briefly introduce the data we used, including 1D Zel'dovich approximation and density fields in IllustrisTNG. In Section 4, we give the numerical results for the matter clustering. Finally, in Section 5, we discuss and summarize our results.

In this section, we briefly review the definition of the CWT based on the GDW in 1D and then introduce the wavelet-based statistical tools, i.e., the WPS, the WCC function, and the WBC. To assess their significance, statistical errors are then discussed. Finally, we give the correspondence between the wavelet scale and the Fourier wavenumber for the GDW. The 3D isotropic CWT and the statistics formulated on it are given in Appendix B.

2.1. GDW and CWT

The GDW used in this work is defined as the first derivative of the Gaussian smoothing function with respect to the scale parameter w, multiplied by a factor $\sqrt{w}$,

$$\begin{eqnarray}\begin{array}{rcl}\psi (w,x) & \equiv & \sqrt{w}\displaystyle \frac{\partial G(w,x)}{\partial w}\\ & = & \displaystyle \frac{\sqrt{w}}{4\sqrt{\pi }}(2-{w}^{2}{x}^{2})\exp \left(-\displaystyle \frac{{w}^{2}{x}^{2}}{4}\right),\end{array}\end{eqnarray} \tag{ 1 }$$

where $G(w,x)=w\exp \left(-{w}^{2}{x}^{2}/4\right)/(2\sqrt{\pi })$ denotes the Gaussian smoothing function with the scale parameter w greater than zero. The prefactor $\sqrt{w}$ on the right-hand side of Equation (1) ensures that the energy of the wavelet is unaffected by the scale parameter. ³ The Fourier transform of the GDW is given by

$$\begin{eqnarray}&&\hat{\psi }(w,k)=\displaystyle \frac{2}{\sqrt{w}}{\left(\displaystyle \frac{k}{w}\right)}^{2}\exp \left[-{\left(\displaystyle \frac{k}{w}\right)}^{2}\right],\end{eqnarray} \tag{ 2 }$$

through which the scale parameter w is defined as the peak wavenumber of $\hat{\psi }(w,k)$, i.e., the wavenumber at which $\hat{\psi }(w,k)$ takes the maximum value, as demonstrated in Figure 1. It can be seen from Equations (1) and (2) that the GDW satisfies the conditions of admissibility, similarity, and regularity. In fact, the GDW is the same as the Mexican hat wavelet in 1D, except for their constant factors. This is illustrated by substituting $w=\sqrt{2}/a$ into Equation (1), where a is the scale parameter for the traditional CWT. ⁴ However, as pointed out by Wang & He (2021), the 3D GDW is an anisotropic separable wavelet function, which is not a 3D Mexican hat wavelet. At present, we are only concerned with the 1D case.

Zoom In Zoom Out Reset image size

The wavelet function ψ(w, x) is then used to perform the WT of a 1D signal f(x) as follows:

$$\begin{eqnarray}&&{W}_{f}(w,x)={\int }_{-\infty }^{+\infty }f(u)\psi (w,x-u){\rm{d}}u,\end{eqnarray} \tag{ 3 }$$

Note that there is no complex conjugate in Equation (3) since both signal and wavelet are real. The WT, Equation (3), is nothing but a convolution of the signal with wavelets at different scales, which can be implemented efficiently by the fast Fourier transform (FFT) technique. As a function of scale w and space x, from the WT we can clearly see how different scale features are localized in space. However, it is not possible to achieve arbitrarily good resolution in space and scale simultaneously (Chui 1997; Addison 2017). As illustrated in Figure 1, a narrower (wider) wavelet provides better (poorer) spatial resolution accompanied by poorer (better) frequency resolution. This fact is quantified in terms of the uncertainty principle ΔxΔk ≳ 1/2 for the GDW, where Δx is the standard deviation of the wavelet in real space and Δk is the standard deviation in frequency (Chui 1997).

The greatest convenience brought to us by the definition of GDW is that the signal can be reconstructed through a single integral of the wavelet coefficients. Combining Equations (1) and (3), we have

$$\begin{eqnarray}&&{W}_{f}(w,x)=\sqrt{w}\displaystyle \frac{\partial {f}_{{\rm{s}}}(w,x)}{\partial w},\end{eqnarray} \tag{ 4 }$$

where f_s(w, x) refers to the smoothed field under the scale w and is given by

$$\begin{eqnarray}&&{f}_{{\rm{s}}}(w,x)={\int }_{-\infty }^{+\infty }f(u)G(w,x-u){\rm{d}}u.\end{eqnarray} \tag{ 5 }$$

By integrating Equation (4) with respect to w, we get the single-integral inverse transform as shown below,

$$\begin{eqnarray}&&f(x)=C+{\int }_{0}^{+\infty }\displaystyle \frac{{W}_{f}(w,x)}{\sqrt{w}}{\rm{d}}w,\end{eqnarray} \tag{ 6 }$$

where the integral constant C is equal to f_s(w = 0, x), and generally C = 0 for, say, the density contrast field of the universe. Compared with the usual reconstruction formula for the CWT found in most wavelet literature, e.g., Addison (2017), our reconstruction formula defined by Equation (6) is much easier to compute numerically and generalize to integrations of higher dimensions.

2.2. Parseval's Theorem for the CWT

Parseval's theorem is an important result in Fourier transform, which states that the inner product between signals is preserved in going from time to the frequency domain. Similarly, there is an analog of Parseval's theorem for the WT (Hudgins et al. 1993), the form of which is expressed as

$$\begin{eqnarray}&&\begin{array}{l}{\displaystyle \int }_{0}^{+\infty }{\displaystyle \int }_{-\infty }^{+\infty }{W}_{f}(w,x){W}_{g}(w,x){\rm{d}}x{\rm{d}}w\\ \,=\,{C}_{\psi }{\displaystyle \int }_{-\infty }^{+\infty }f(u)g(u){\rm{d}}u,\end{array}\end{eqnarray} \tag{ 7 }$$

where W_f (w, x) and W_g (w, x) are WTs for f(x) and g(x), respectively, ${C}_{\psi }={\int }_{0}^{+\infty }(\hat{\psi }{(1,k)}^{2}/k){\rm{d}}k$ is the admissibility constant, and equals 1/2 for the GDW. If f and g are the same, then we have

$$\begin{eqnarray}&&\begin{array}{l}{\displaystyle \int }_{0}^{+\infty }{\displaystyle \int }_{-\infty }^{+\infty }| {W}_{f}(w,x){| }^{2}{\rm{d}}x{\rm{d}}w\\ =\,{C}_{\psi }{\displaystyle \int }_{-\infty }^{+\infty }| f(u){| }^{2}{\rm{d}}u.\end{array}\end{eqnarray} \tag{ 8 }$$

2.3. WPS and WCC

We assume that the signal f(x) satisfies the periodic boundary condition with period L_b, which is a usual choice for studying a typical region of the universe. Then for a subregion of the signal L ≤ L_b, the local WPS is defined as

$$\begin{eqnarray}&&{P}_{f,L}^{W}(w)=\displaystyle \frac{1}{L}{\int }_{L}| {W}_{f}(w,x){| }^{2}{\rm{d}}x.\end{eqnarray} \tag{ 9 }$$

We can see that Equation (9) refers to the variance at the scale w within the spatial region L, since ∣W_f (w, x)∣² is just the variance per area at the space-scale plane, as referred to by the Parseval's theorem, Equation (8). According to the Fourier Parseval's theorem for periodical signals, we can obtain the relationship between the global wavelet and Fourier power spectrum, which is given by

$$\begin{eqnarray}\begin{array}{rcl}{P}_{f,{L}_{b}}^{W}(w) & = & \displaystyle \frac{1}{{L}_{b}}{\displaystyle \int }_{0}^{{L}_{b}}| {W}_{f}(w,x){| }^{2}{\rm{d}}x\\ & = & \displaystyle \frac{1}{{L}_{b}}\displaystyle \sum _{k}\displaystyle \frac{| {\hat{W}}_{f}(w,k){| }^{2}}{{L}_{b}}\\ & = & \displaystyle \sum _{k}{P}_{\psi }^{F}(w,k){P}_{f}^{F}(k),\end{array}\end{eqnarray} \tag{ 10 }$$

where ${\hat{W}}_{f}(w,k)=\hat{\psi }(w,k)\hat{f}(k)$ is the Fourier transform of W_f (w, x), ${P}_{f}^{F}(k)=| \hat{f}(k){| }^{2}/{L}_{{\rm{b}}}$ is the Fourier power spectrum of the signal, and ${P}_{\psi }^{F}(w,k)=| \hat{\psi }(w,k){| }^{2}/{L}_{{\rm{b}}}$ denotes the Fourier power spectrum of the wavelet. Obviously, the global WPS of a signal is the average of its Fourier power spectrum weighted by the Fourier power spectrum of the corresponding wavelet function over all wavenumbers.

Then given two signals f(x) and g(x) with WTs W_f (w, x) and W_g (w, x), we can define the cross-wavelet transform (XWT) as

$$\begin{eqnarray}&&{{\text{XWT}}}_{f-g}(w,x)={W}_{f}(w,x){W}_{g}(w,x),\end{eqnarray} \tag{ 11 }$$

i.e., it measures the local covariance at each spatial position and scale, as revealed by Equation (7). By integrating the XWT over a finite spatial region, we get the normalized local WCC as follows:

$$\begin{eqnarray}&&{C}_{f-g,L}^{W}(w)=\displaystyle \frac{1}{L}\displaystyle \frac{{\int }_{L}{{\text{XWT}}}_{f-g}(w,x){\rm{d}}x}{{\left({P}_{f,L}^{W}{P}_{g,L}^{W}\right)}^{1/2}},\end{eqnarray} \tag{ 12 }$$

which can take on values between −1 (perfect anticorrelation) and 1 (perfect correlation).

2.4. WBC and Error Estimation

The main content of this subsection is based on van Milligen et al. (1995a, 1995b), and we refer the interested readers to these two references for the details.

The Fourier bispectrum is the lowest order statistic that measures the amount of phase coupling of harmonic modes within a signal. By analogy, the WBS is given as

$$\begin{eqnarray}\begin{array}{rcl}{B}_{f,L}^{W}({w}_{1},{w}_{2}) & = & \displaystyle \frac{1}{L}{\displaystyle \int }_{L}{W}_{f}({w}_{1},x){W}_{f}({w}_{2},x)\\ & & \times \,{W}_{f}(w,x){\rm{d}}x,\end{array}\end{eqnarray} \tag{ 13 }$$

where w₁ + w₂ = w (frequency sum rule). The WBS measures the nonlinear interplay within the local region L between scale components w₁, w₂, and w such that the sum rule is satisfied. In the case of completely random phases of the signal, ${B}_{f,L}^{W}({w}_{1},{w}_{2})$ is statistically zero. However, once a coherent structure is formed by the phase coupling, ${B}_{f,L}^{W}({w}_{1},{w}_{2})$ will take significant nonzero values. The WBS usually is normalized in the following way:

$$\begin{eqnarray}\begin{array}{cc} & {b}_{f,L}^{W}({w}_{1},{w}_{2})\,=\\ \, & \,\,{\left[\displaystyle \frac{\left|{B}_{f,L}^{W}({w}_{1},{w}_{2}){| }^{2}\right.}{\tfrac{1}{L}{\displaystyle \int }_{L}{\left|\left.{W}_{f}({w}_{1},x){W}_{f}({w}_{2},x\right)\right|}^{2}{\rm{d}}x{P}_{f,L}^{W}}\right]}^{1/2},\end{array}\end{eqnarray} \tag{ 14 }$$

which is called the WBC, attaining values between 0 and 1. Throughout this paper we will use it to measure the nonlinear behaviors of matter clustering instead of the WBS. In addition, it is convenient to introduce the summed WBC defined as

$$\begin{eqnarray}&&{b}_{f,L}^{W}(w)={\left[\displaystyle \frac{1}{s(w)}\sum _{{w}_{1},{w}_{2}}{[{b}_{f,L}^{W}({w}_{1},{w}_{2})]}^{2}\right]}^{1/2},\end{eqnarray} \tag{ 15 }$$

where the summation is taken over all w₁ and w₂ such that w₁ + w₂ equals w and s(w) is the number of summands in the summation. In addition, as a measure of the total nonlinearity in the chosen region of the signal, the total WBC is defined by averaging the squared WBC over all points in the scale-scale plane as

$$\begin{eqnarray}&&{b}_{f,L}^{W}={\left[\displaystyle \frac{1}{S}\sum _{{w}_{1},{w}_{2}}{[{b}_{f,L}^{W}({w}_{1},{w}_{2})]}^{2}\right]}^{1/2},\end{eqnarray} \tag{ 16 }$$

where S is the total number of points (w₁, w₂) in the scale-scale plane.

In practice for discrete sampled signals, integrations over the interval L involved in calculating the statistics mentioned above are carried out with summation over N sample points. By the law of large numbers, an integration over L suffers a relative statistical error of $1/\sqrt{N}$. In addition, the fact that CWTs are non-orthogonal leads to the wavelet coefficients not all being statistically independent (van Milligen et al. 1995a, 1995b, 1997). If continuous wavelets ψ(x) and ψ^*(x + d₀) satisfy ∫ψ(x)ψ^*(x +d₀)dx = 0, then they can be regarded as approximately orthogonal and the corresponding wavelet coefficients are statistically independent. In the case of the GDW, independent wavelet coefficients are separated by distance d₀/w at each scale, where ${d}_{0}=2\sqrt{3-\sqrt{6}}$. Then the number of statistical independent points on the interval L is $N^{\prime} =L/({d}_{0}/w)\,=2\pi {Nw}/({d}_{0}{k}_{\mathrm{samp}})$, where k_samp = 2π/Δx is the sampling frequency. Thus for the WPS, its statistical error is estimated by

$$\begin{eqnarray}&&\begin{array}{l}\epsilon \left[{P}_{f,L}^{W}(w)\right]\approx \displaystyle \frac{{P}_{f,L}^{W}(w)}{\sqrt{N^{\prime} }}\\ =\,{P}_{f,L}^{W}(w){\left(\displaystyle \frac{{d}_{0}{k}_{\mathrm{samp}}}{2\pi {Nw}}\right)}^{1/2}.\end{array}\end{eqnarray} \tag{ 17 }$$

By applying similar estimates for all integral terms in Equation (12) and according to error propagation, we obtain the statistical noise level for the WCC,

$$\begin{eqnarray}&&\epsilon \left[| {C}_{f-g,L}^{W}(w)| \right]\approx {\left(\displaystyle \frac{3{d}_{0}{k}_{\mathrm{samp}}}{4\pi {Nw}}\right)}^{1/2}.\end{eqnarray} \tag{ 18 }$$

Equation (18) is called the statistical noise level because it is the cross-correlation value that can be achieved by Gaussian noise and is caused by using a limited number of values in the integration. Notice that this noise level is scale dependent, suggesting that the interpretation of the signal becomes increasingly significant as the scale decreases. Just like the approach above, we can obtain the noise level of the WBC as shown below,

$$\begin{eqnarray}&&\epsilon \left[{b}_{f,L}^{W}({w}_{1},{w}_{2})\right]\approx {\left(\displaystyle \frac{{d}_{0}{k}_{\mathrm{samp}}}{2\pi N\mathrm{Min}({w}_{1},{w}_{2},{w}_{1}+{w}_{2})}\right)}^{1/2}.\end{eqnarray} \tag{ 19 }$$

2.5. Relationship between Scale and Wavenumber

To facilitate the comparison between wavelet and Fourier spectra, we need to ascertain the relationship between the scale parameter w and the equivalent Fourier wavenumber. Meyers et al. (1993) and Torrence & Compo (1998) suggest that the relationship between them can be derived analytically for a particular wavelet function by performing a WT of a cosine wave with a known period, such as $\cos ({k}_{F}x)$, and then computing the scale w at which the scalogram reaches its maximum. Following their method, the scalogram of $\cos ({k}_{F}x)$ is first computed with the GDW, and the result is

$$\begin{eqnarray}&&| {W}_{\cos }(w,x){| }^{2}=\displaystyle \frac{4{k}_{F}^{4}{\cos }^{2}({k}_{F}x)\exp (-\tfrac{2{k}_{F}^{2}}{{w}^{2}})}{{w}^{5}},\end{eqnarray} \tag{ 20 }$$

which is depicted in Figure 2. The scale parameter that makes $| {W}_{\cos }(w,x){| }^{2}$ take the maximum value should be equivalent to the wavenumber k_F , since the Fourier power spectrum of a cosine wave is an impulse at k_F . Therefore, by solving $\partial | {W}_{\cos }(w,x){| }^{2}/\partial w=0$, we get the correspondence between the wavelet scale and Fourier wavenumber,

$$\begin{eqnarray}&&w=\displaystyle \frac{2}{\sqrt{5}}{k}_{F}\approx 0.89{k}_{F}.\end{eqnarray} \tag{ 21 }$$

The relation of Equation (21) shows that the wavelet scale is proportional to the Fourier wavenumber. ⁵ In the next sections, we will present our results in terms of the equivalent Fourier wavenumber k_F instead of w, and we drop the subscript "F" for simplicity of notation.

Zoom In Zoom Out Reset image size

For the convenience of the readers, in Table 1, we list all the notations used in our paper, with their meanings and corresponding acronyms.

Table 1. Notations Used in This Paper with Their Meanings and Acronyms ^a

Notation	Meaning	Acronym
ψ(w, x)	Gaussian-derived wavelet	GDW
$\hat{\psi }(w,k)$	Fourier transform of the GDW
Wf (w, x)	Wavelet transform	WT
${\hat{W}}_{f}(w,k)$	Fourier transform of the WT
${P}_{f,L}^{W}(w)$	Wavelet power spectrum	WPS
${P}_{f}^{F}(k)$	Fourier power spectrum
XWTf−g (w, x)	Cross-wavelet transform	XWT
${C}_{f-g,L}^{W}(w)$	Wavelet cross correlation	WCC
${B}_{f,L}^{W}({w}_{1},{w}_{2})$	Wavelet bispectrum	WBS
${b}_{f,L}^{W}({w}_{1},{w}_{2})$	Wavelet bicoherence	WBC
${b}_{f,L}^{W}(w)$	Summed WBC
${b}_{f,L}^{W}$	Total WBC

Note.

^a In this table, we only list the symbols associated with the 1D CWT. The same notation convention is also used for the quantities based on the 3D isotropic CWT, except that the scalars x and k are replaced with vectors r and k, and the length L is replaced with volume V.

Download table as:  ASCIITypeset image

In this work, we use two types of data, i.e., the 1D density fields obtained from a Zel'dovich approximation and the 3D density fields of the IllustrisTNG simulations. In the 1D case, the Zel'dovich approximation provides the exact nonlinear solution for the perturbative equations of collisionless matter up to the first appearance of orbit-crossing singularities (Soda & Suto 1992). Due to easy implementation and fast computation of the 1D Zel'dovich exact solution, we will demonstrate the usefulness of our method with it before further analyzing the cosmological simulation IllustrisTNG.

3.1. Zel'dovich Approximation in One Dimension

The fundamental idea of the Zel'dovich approximation is the transformation between Eulerian and Lagrangian coordinates, i.e.,

$$\begin{eqnarray}&&x(t,q)=q-\theta (t)f(q),\end{eqnarray} \tag{ 22 }$$

where x and q are the Eulerian and Lagrangian coordinates, respectively. Then by applying the mass conservation to Equation (22), the density contrast is given explicitly by

$$\begin{eqnarray}&&\delta \left(\eta ,t\right)+1=\displaystyle \frac{1}{1-\theta \left(t\right)F\left(\eta \right)},\end{eqnarray} \tag{ 23 }$$

where $F(\eta )=\tfrac{{\rm{d}}f(\eta {L}_{b})}{{\rm{d}}\eta }/{L}_{b}$ with η = q/L_b being the dimensionless coordinates divided by the length size of the density field, and θ(t) is the growth factor. For simplicity, we normalize θ(t) to unity at the initial time and use it as a time variable instead of t. It is easy to see from Equation (23) that F(η) is the initial density contrast if small enough. For more details on the 1D Zel'dovich approximation, we refer the readers to Soda & Suto (1992) and Fujiwara & Soda (1996). The initial condition is set to be

$$\begin{eqnarray}&&F(\eta )=2\sum _{k\gt 0}\sqrt{{P}_{{\rm{i}}}(k)}[{B}_{k}\cos (k\eta )+{C}_{k}\sin (k\eta )],\end{eqnarray} \tag{ 24 }$$

where B_k and C_k are drawn from a Gaussian with a standard deviation of 1. We impose a periodic boundary condition on the interval η ∈ [0, 1], then divide it into 1024 equally spaced segments. So the wavenumber, as an integer multiple of 2π, has a maximum value of 512 × 2π, i.e., the Nyquist frequency. In this work, we assume that the spectral index of the initial power spectrum _Pi(k) = Akⁿ in Equation (24) is equal to n = −2. The amplitude A is chosen to be 2.5 × 10⁻⁶ such that the initial density perturbation is between −0.01 and 0.01. Hence, the evolution of the density field is totally determined by Equations (23) and (24).

As shown in Figure 3, we select the density fields at θ = 10, 20, 40, 60, 80, and 100 to examine their evolution from linear to nonlinear stages.

Zoom In Zoom Out Reset image size

3.2. IllustrisTNG Data

4.1. Wavelet Analysis of 1D Zel'dovich Density Fields

From Figure 3 we can see that the characteristics of the 1D density growth in our case are similar to the 3D simulations. With the evolution of this simple system from linear (δ ≪ 1) to the highly nonlinear regime (δ ≫ 1) due to gravitational effects, structures become increasingly significant. For example, there is an obvious underdense region roughly ranging from η ∼ 0.3 to ∼ 0.5, as well as a large overdense region next to it. Since these two distinct local environments in the universe are thought to have different effects on matter clustering (Abbas & Sheth 2005), we hope that wavelet methods are able to distinguish between them in our 1D toy model. In this section, we present the results for 1D matter clustering based on the CWT analysis.

4.1.1. Space-scale Decomposition of the Density Fields

As stated by Equation (3), the CWT is defined for the infinite input signal. Hence, the finite signal must be padded with some values before the transform is performed. The usual padding schemes include zero padding, decay padding, periodic padding, and symmetric padding (Addison 2017). In the present work, the density contrasts are periodically padded with themselves since they satisfy the periodic boundary condition. Then by employing the CWT, different features are picked out by GDW at each scale while retaining positional information, as illustrated in Figure 4. From this figure, we can see that the wavelet coefficients W_θ (k, η) take values from negative (blue) to positive (red), reflecting that the density field is anti-correlated and correlated with the GDW. Positive wavelet coefficients correspond to the location of the density peaks, and those negative coefficients are distributed between the density peaks. This is due to the fact that the shape of the GDW matches well with the density peaks, which is more evident in Figure 15.

Zoom In Zoom Out Reset image size

Let us focus on ∣W_θ (k, η)∣. The characteristics of matter clustering will be seen qualitatively from the space-scale plane. From θ = 10–20, ∣W_θ (k, η)∣ evolves little with time and is dominated by large-scale components with a relatively random spatial distribution. This indicates that the density field is almost homogeneous. At θ = 40, the small-scale components start to become apparent owing to the gravitational interactions. Afterward, ∣W_θ (k, η)∣ progressively increases with time. As a consequence, strongly structured patterns are formed in the space-scale plane at θ = 100, suggesting that the density field is highly nonhomogeneous at this time. It is noteworthy that all scale components grow very significantly in the region from η ∼ 0.5 to ∼ 0.8. In contrast, there are almost no small-scale components generated between η ∼ 0.3 and ∼ 0.5, while there is a moderate increase at small scales in other regions. Based on these simple analyses, we find that the CWT of the density contrast reproduces its local evolutionary features very well.

4.1.2. WPS Measurements

By averaging squared wavelet coefficients over a local spatial region, we thus obtain the local WPS. For a better understanding of the WPS, we first compare the global wavelet and Fourier power spectrum of 1D Zel'dovich fields. A main feature visible in Figure 5 is that the WPSs are smoother than their Fourier counterparts, as expected by Equation (10). These two types of power spectra have approximately the same trends with comparable amplitudes. In fact, the global WPS is able to reproduce the correct exponent of the power-law Fourier power spectrum (Wang & He 2021), although that is not obvious for our 1D noisy data. From the WPS shown in Figure 5, we can see that the strength of clustering on small scales increases with time. Furthermore, by measuring the change in the WPS at each time relative to the initial WPS, it should help better demonstrate the nonlinear evolution of clustering. Since δ(η, t) approximately equals θ(t)F(η) in the linear regime, there is a relation between the WPS of δ(η, t) and the initial one given by

$$\begin{eqnarray}&&{P}_{\theta ,{L}_{b}}^{W}(k)\approx {\theta }^{2}{P}_{{\rm{i}},{L}_{b}}^{W}(k).\end{eqnarray} \tag{ 25 }$$

Hence, if the ratio of these two WPSs, ${P}_{\theta ,{L}_{b}}^{W}(k)/{P}_{{\rm{i}},{L}_{b}}^{W}(k)$, is not scale dependent, then the growth of the density fluctuations is linear. We show this ratio for times θ = 10, 20, 40, 60, 80, and 100 in Figure 6. Obviously, ${P}_{\theta ,{L}_{b}}^{W}(k)/{P}_{{\rm{i}},{L}_{b}}^{W}(k)$ is almost constant over the entire scale range at θ = 10 and 20, suggesting that the density contrast is linear. However, this ratio becomes scale dependent at the time θ = 40 and there is a mild bump around k ≈ 300, which means that there are a few newly generated scale components and the density contrast is quasi-nonlinear. Since then, this bump has been enhanced by nonlinear effects with time to θ = 80, while the ratio remains approximately flat on large scales of k ≲ 20. At the time θ = 100, the ratio increases more significantly on small scales thus leading to a plateau at scales of k ≳ 300, implying that the density field has developed into the highly nonlinear regime.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Next, we turn our attention to the environmental dependence of matter clustering. To do this, we split the density field into four consecutive segments labeled I–IV. The spatial extent spanned by each segment and the mean density within that local space at each time are listed in Table 2. The evolutionary trends of these four local mean densities are more clearly depicted in Figure 7. We see that Segment III represents the highly overdense environment, and both Segments I and IV represent the slightly overdense environment, while Segment II represents the underdense environment. These different density environments are expected to exhibit distinctly different clustering characteristics from each other. Therefore, it is instructive to examine the local WPS for each environment, as shown in Figure 8. Each subplot shows the WPS for different regions at the same time. At linear stages (θ = 10 and 20), all of these local power spectra roughly converge to the global WPS. Starting from θ = 40, however, the divergence between these local power spectra becomes more and more significant with time. For the highly overdense environment, i.e., Segment III, its WPS is more enhanced than all other regions over the entire scale range. In contrast, the WPS for the underdense environment (Segment II) is more suppressed than all other regions on scales of k ≳ 40. For Segments I and IV, as slightly overdense regions, the amplitudes of their WPS are very close to each other and fall between Segments II and III on scales of k ≳ 40. Notice that as the scale becomes larger, the statistical error in the power spectrum becomes increasingly larger to the extent that it does not give a significant description of matter clustering in the range k ≲ 40, due to the wavelet coefficients being non-orthogonal and the length of the local region being too short. Even so, the WPS is still able to give a meaningful interpretation over a relatively large-scale range.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2. Density Field at Each Time Split into Four Consecutive Segments ^a

Segment	Spatial Range	Mean Density Contrast ${\overline{\delta }}_{L}$
		θ = 10	20	40	60	80	100
I	0.00 ≲ η ≲ 0.27	0.012	0.025	0.056	0.095	0.144	0.206
II	0.27 ≲ η ≲ 0.54	−0.045	−0.085	−0.155	−0.214	−0.265	−0.309
III	0.54 ≲ η ≲ 0.78	0.044	0.094	0.216	0.381	0.629	1.102
IV	0.78 ≲ η ≲ 1.00	0.000	0.001	0.007	0.018	0.037	0.066

Note.

^a The spatial range spanned by each segment and the corresponding local mean density contrast at times of θ = 10, 20, 40, 60, 80, and 100 are listed.

Download table as:  ASCIITypeset image

In Figure 9, we also measure the local WPS at each time relative to its corresponding initial WPS within the scale range k ≳ 40 where statistical errors are smaller. For all segments, we see that ${P}_{\theta ,L}^{W}(k)/{P}_{{\rm{i}},L}^{W}(k)$ approximates to θ² at stages of θ = 10 and 20. However, at late stages, the evolutionary trends of Segments II and III are completely opposite. For Segment II, the growth of small-scale components is getting slower, but for Segment III, the small-scale components are growing faster and faster. Although the WPSs of Segments I and IV are very similar, there are visible differences between them. Specifically, the former has a slightly upward tilt on small scales, while the latter remains roughly horizontal.

Zoom In Zoom Out Reset image size

Based on these facts, we find that WPS is fully capable of detecting the effects of density environments on matter clustering. However, as pointed out above, the density field evolves into the nonlinear regime at late times, indicating that there are couplings between different scale components. The WPS, as two-point statistics, cannot determine such scale coupling since it does not contain phase information. To detect scale coupling, we measure the WBC, which is the lowest order statistics sensitive to nonlinear couplings between different scales of matter density field.

4.1.3. WBC Measurements

Figure 10 shows the time evolution of matter WBC with statistical noise subtracted, i.e., ${b}_{\theta ,{L}_{b}}^{W}({k}_{1},{k}_{2})-\epsilon [{b}_{\theta ,{L}_{b}}^{W}({k}_{1},{k}_{2})]$. Note that ${b}_{\theta ,{L}_{b}}^{W}({k}_{1},{k}_{2})$ is symmetric about the diagonal where k₁ = k₂, as indicated by Equations (13) and (14). So we will only address the part below (or above) this diagonal. For visual convenience, values of bicoherence less than the noise level are not considered because they are less significant physically. To better understand the statistical noise level, as suggested by van Milligen et al. (1995a, 1995b), we perform FFT on the highly nonlinear density field at time θ = 100, and give each Fourier component a random phase while maintaining its amplitude, then we perform inverse FFT to get a new set of data. Such a new phase-randomized density field is expected to have no structures induced by scale coupling, while its power spectrum is identical to that of the raw density field, in which structures are well formed, as is illustrated in Figure 11. It can be seen that the summed WBC of the density field at θ = 100 is much higher than that of the phase-randomized version, which falls below the noise level along with the initial density field. Thus, the statistical noise level provides us a criterion to discriminate matter distributions with structures formed from those without. Accordingly, looking at Figures 10 and 11 together, there is a very weak coupling between large scales around k ∼ 10 and intermediate scales around k ∼ 300 at times θ = 10 and 20. Gradually, the scale coupling becomes more and more significant with time and eventually spreads to the entire scale range, which implies that the matter distribution becomes more and more structured.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 12, we consider the evolution of the total WBC with time in the consecutive Segments I–IV. An interesting phenomenon is that the degree of nonlinearity in Segment II, the underdense region, is much higher than the noise level at the initial time, and this nonlinearity is less evolved. On the other hand, the nonlinearity is more evolved in overdense regions, although it is very weak at the initial time. In particular, Segment III, as a highly overdense region, has the strongest nonlinearity at late times. This implies that structure formation occurs mainly in Segment III. Furthermore, we examine the scale coupling for these segments by measuring their WBC and summed WBC at different times and corresponding results are shown in Figures 13 and 14. Consistent with the results in Figure 12, the estimated WBC and summed WBC in Segment II exhibit significant values between widely separated scales at the initial time because the matter distribution within such an underdense region is left-skewed. Both the WBC and summed WBC in Segment II barely evolve with time, indicating that there is almost no structure formation in this region. Also, Segment III shows a notable coupling between scale bands 20 ≲ k ≲ 100 and 100 ≲ k ≲ 500 at the initial stages, due to the matter distribution in such a highly overdense region being right-skewed. Unlike Segment II, the WBC and summed WBC of Segment III evolve the most dramatically, which indicates that the matter becomes increasingly clustered. For Segments I and IV, the bicoherences are below or close to the statistical noise level at the initial stages and evolve slowly with time, which means that the matter experiences moderate structure formation in these two slightly overdense regions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.2. Wavelet Analysis of the Density Fields of the TNG100 Simulation

In previous sections, we focused on the wavelet analysis of the nonlinear evolution of the 1D Zel'dovich density field, revealing the considerable potential of the CWT for characterizing the matter clustering. Therefore, we further apply the 3D continuous wavelet method to the analysis of the TNG100 simulation. For simplicity, we consider the 3D isotropic CWT methods based on the GDW described in Appendix B, which is the simplest extension of the 1D case. We first perform the isotropic CWT of the 3D density fields of the total matter, DM, gas, and stars ⁶ at redshift z = 0. In order to avoid the aliasing, we take the maximum scale as ${k}_{\max }={k}_{\mathrm{Nyq}}/2$, where k_Nyq =1024π/L_box is the Nyquist frequency of the density grid. For the purpose of presentation, WTs of density slices at three scales are shown in Figure 15. By visual inspection, We see that for all matter components, only large voids and filaments are captured at the large scale of k ≈ 0.35 hMpc⁻¹, whereas as the scale becomes smaller, the WT is increasingly accurate in tracking the finer web-like structure of the universe. More quantitative statements can be made by using statistics formulated on the WT.

Zoom In Zoom Out Reset image size

In Figure 16, we compare the global wavelet and Fourier power spectra for different matter distributions. It can be seen that these two types of spectra for each matter have a similar shape. The global WPS of the total matter falls slightly below that of the DM at small scales. The gas becomes significantly less clustered than DM at intermediate and small scales, while the stellar mass shows a very strong clustering. However, the magnitude of the WPS is different from the Fourier power spectrum. More precisely, the deviations between these two spectra are measured by the ratio between them, which appears to be roughly a constant determined by Equation B13.

Zoom In Zoom Out Reset image size

To study the impact of a large-scale environment on the matter clustering, we divide the whole simulation volume ${V}_{b}={L}_{b}^{3}$ into ${N}_{{\rm{sub}}}^{3}$ sub-volumes of $V={({L}_{b}/{N}_{{\rm{sub}}})}^{3}$, where N_sub is set to be 8 and hence the side length L_b /N_sub is 9.375 h⁻¹Mpc. In Figure 17, we show the local WPSs of the total matter spatial distribution measured in the sub-volumes. These local WPSs are separated into six categories, according to the mean densities of sub-volumes as follows:

• 1.  
  ${\bar{\delta }}_{V}\gt 3$,
• 2.  
  $1\lt {\bar{\delta }}_{V}\lt 3$,
• 3.  
  $0\lt {\bar{\delta }}_{V}\lt 1$,
• 4.  
  $-0.5\lt {\bar{\delta }}_{V}\lt 0$,
• 5.  
  $-0.75\lt {\bar{\delta }}_{V}\lt -0.5$,
• 6.  
  $-1\lt {\bar{\delta }}_{V}\lt -0.75$,

where ${\bar{\delta }}_{V}=\langle {\delta }_{\mathrm{tot}}({\bf{r}}){\rangle }_{V}$ denotes the local mean density of the total matter within a sub-volume. We see that the local WPSs are greater for sub-volumes with higher local mean density for scales up to a few times 0.1 hMpc⁻¹. To make it more explicit, we consider the mean of the local WPSs over all sub-volumes within each density interval, which is shown below:

$$\begin{eqnarray}&&{P}_{\mathrm{tot},{\bar{\delta }}_{V}}^{W}(k)=\langle {P}_{\mathrm{tot},V}^{W}(k){\rangle }_{{\rm{\Delta }}{\bar{\delta }}_{V}},\end{eqnarray} \tag{ 26 }$$

where the density interval ${\rm{\Delta }}{\bar{\delta }}_{V}$ is specified by the abovementioned list (i)–(vi). These results are presented in Figure 18. Obviously, at the largest scales, ${P}_{\mathrm{tot},{\bar{\delta }}_{V}}^{W}(k)$ in all density intervals converge to the global WPS, while as scales become smaller, discrepancies between ${P}_{\mathrm{tot},{\bar{\delta }}_{V}}^{W}(k)$ in different density intervals as well as the global WPS become more pronounced. We noticed that the mean local WPS in ${\bar{\delta }}_{V}\gt 3$ is the most enhanced among all density intervals, particularly at scales around 3 hMpc⁻¹, where they are dominated by the most massive halos, as pointed out by van Daalen & Schaye (2015). As the local mean density ${\bar{\delta }}_{V}$ gets smaller, the mean local WPS gets smaller. In ${\bar{\delta }}_{V}\lt 1$, the mean local WPSs are most severely suppressed on scales of k ∼ 3 hMpc⁻¹. Therefore, our results suggest that most massive halos tend to reside in the densest environments, while those less dense environments show a deficit of massive halos, which are qualitatively consistent with the results of Fisher & Faltenbacher (2018) and Zhu et al. (2022).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Note that the total matter density field is formed via

$$\begin{eqnarray}&&{\bar{\rho }}_{\mathrm{tot}}{\delta }_{\mathrm{tot}}({\bf{r}})={\bar{\rho }}_{\mathrm{dm}}{\delta }_{\mathrm{dm}}({\bf{r}})+{\bar{\rho }}_{\mathrm{gas}}{\delta }_{\mathrm{gas}}({\bf{r}})+{\bar{\rho }}_{\star }{\delta }_{\star }({\bf{r}}).\end{eqnarray} \tag{ 27 }$$

Combining this equation and Equation (B6), the local WPS of the total matter in a sub-volume can be separated into

$$\begin{eqnarray}\begin{array}{rcl}{\bar{\rho }}_{\mathrm{tot}}^{2}{P}_{\mathrm{tot},V}^{W}(k) & = & {\bar{\rho }}_{\mathrm{dm}}^{2}{P}_{\mathrm{dm},V}^{W}(k)+2{\bar{\rho }}_{\mathrm{dm}}{\bar{\rho }}_{\mathrm{gas}}{P}_{\mathrm{dm}-\mathrm{gas},V}^{W}(k)\\ & & +{\bar{\rho }}_{\mathrm{gas}}^{2}{P}_{\mathrm{gas},V}^{W}(k)+2{\bar{\rho }}_{\mathrm{gas}}{\bar{\rho }}_{\star }{P}_{\mathrm{gas}-\star ,V}^{W}(k)\\ & & +{\bar{\rho }}_{\star }^{2}{P}_{\star ,V}^{W}(k)+2{\bar{\rho }}_{\mathrm{dm}}{\bar{\rho }}_{\star }{P}_{\mathrm{dm}-\star ,V}^{W}(k),\end{array}\end{eqnarray} \tag{ 28 }$$

in which the cross terms ${P}_{f-g,V}^{W}(k)$ encode the coherence between different matter species. This can be examined by the local WCC ${C}_{f-g,V}^{W}(k)$, which is defined as Equation (B7). The mean values of local WCCs in different density intervals ${C}_{f-g,{\bar{\delta }}_{V}}^{W}(k)=\langle {C}_{f-g,V}^{W}(k){\rangle }_{{\rm{\Delta }}{\bar{\delta }}_{V}}$ are shown in Figure 19. On large scales, ${C}_{f-g,{\bar{\delta }}_{V}}^{W}(k)$ is close to one for all pairs of matter components in all density environments, whereas it deviates more from one on smaller scales, and also shows more distinct discrepancies between different density environments. Specifically, in ${\bar{\delta }}_{V}\gt -0.5$, the correlation between the DM and stars remains higher than other pairs on scales of k ≳ 3 hMpc⁻¹. However, this correlation decreases with decreasing density. This could be explained by the fact that stars are mostly formed in the center of halos, which are very few in less dense environments. For the DM and gas pair, the correlation between them is very close to 1 and hardly varies with the density environment on the scales of k ≲ 2 hMpc⁻¹. On scales smaller than 2 hMpc⁻¹, this correlation first decreases with decreasing density until reaching the minimum in $-0.5\lt {\bar{\delta }}_{V}\lt 0$ and then increases with decreasing density, so that it becomes even higher than that between the DM and stars in ${\bar{\delta }}_{V}\lt -0.75$. For the gas and star pair, the correlation between them shows a density dependence similar to that between the DM and gas, but it is still lower than other pairs of matter in all density environments. Within the scale range we investigate, the variation of gas-DM and gas-star correlations with scale and density may be related to the feedback processes, mainly active galactic nuclei (AGNs) feedback at z = 0, which can expel gas from halos to large distances and hence reshape the matter distribution.

Zoom In Zoom Out Reset image size

In order to investigate the impact of the AGN feedback on the matter clustering in different density environments, we consider the ratio between the local WPS from TNG100-1 and that from the corresponding DM-only simulation TNG100-1-Dark, i.e., ${P}_{f,V}^{W}(k)/{P}_{\mathrm{DM}-\mathrm{only},V}^{W}(k)$, where f refers to the total matter, DM, and total baryons (gas+stars), respectively. The mean values of ${P}_{f,V}^{W}(k)/{P}_{\mathrm{DM}-\mathrm{only},V}^{W}(k)$ in different density intervals are plotted in Figure 20. For the total matter, we see that its global WPS is reduced on scales k ≳ 1 hMpc⁻¹ relative to the DM-only's global WPS, with a 14% suppression at k ∼ 20 hMpc⁻¹. This effect is also seen in the local WPS, but with some variation in different density environments. The change in the total matter clustering compared to that in the DM-only simulation comes from two aspects: (1) the redistribution of baryons by nongravitational physics, and (2) the change in the DM distribution resulted from the gravitational coupling of baryons and DM, which is called the back reaction (van Daalen et al. 2020). We see from Figure 20 that the former effect is the determinant of the change in the clustering of total matter for all density intervals. In particular, the clustering strength of baryons is the least suppressed in the density ranges of ${\bar{\delta }}_{V}\gt 3$ and ${\bar{\delta }}_{V}\lt -0.75$, revealing that the AGN feedback is the weakest in the densest and most underdense environments. Considering the fact that most massive halos prefer to reside in the densest environments, our results are consistent with the findings of Springel et al. (2018), showing that AGN feedback has a minimal impact on the most massive halos. As for the lowest AGN feedback in the most underdense environments, it may be related that these environments contain few or no halos.

Zoom In Zoom Out Reset image size

The clustering of matter is an intricate process. In addition to redshift and scale dependence, the density environment plays a very important role in the clustering. When the usual Fourier statistics are utilized to examine the scale dependence, we are unable to take into account the effect of the environment, since the local density information contained in the physical space is smeared out. Therefore, we need better tools to consider the dependence of both scale and environment simultaneously, and the continuous WT is such a tool with good performance.

In this work, we introduce some statistical quantities formulated on the specific wavelet function, GDW, which is a real valued continuous wavelet designed by taking the first derivative of the Gaussian smoothing function with respect to the scale parameter. With such a wavelet construction scheme, the original signal can be recovered from wavelet coefficients easily. These wavelet statistics we used include the WPS, the WCC, and the WBC.

To reveal the usefulness of wavelet statistics in analyzing matter clustering, we analyze the time evolution of the 1D density distributions obtained from the Zel'dovich approximation by measuring their wavelet power spectra and bicoherences. Measurements show that the global WPS on small scales increases more significantly with time than that on large scales, which is generally in agreement with the Fourier case.

To manifest the capability of wavelet statistics to perform local spectral analysis, we divide the 1D Zel'dovich density field at each time into four consecutive segments. All the local wavelet power spectra for these segments almost converge to the global WPS at linear stages. However, the difference between the local power spectra and the global WPS becomes progressively larger with time due to nonlinear effects. In particular, the WPS of Segment III (highly overdense environment) is significantly greater than those of the other segments on scales k ≳ 40, where statistical errors are smaller, which implies that structures on these scales are generated in this region at nonlinear stages. Another striking feature is that the growth of the WPS in Segment II (underdense environment) is severely suppressed on scales k ≳ 40, meaning that there are very few structures on these scales generated in this region at later times. Moreover, measurements of the WBCs show that the scale coupling occurs mainly in Segment III, whereas there is almost no scale coupling in Segment II at late epochs. Both the WPS and the WBC in Segments I and IV exhibit similar behaviors to and fall between those in Segments II and III, probably because Segments I and IV are slightly overdense environments. These results demonstrate the great potential of wavelet statistics for taking into account the effects of both environment and scale on matter clustering.

Furthermore, we generalize the 1D CWT methods to the 3D isotropic version and apply them to the 3D density fields of the TNG100-1 simulation at redshift z = 0. Due to space limitations, WBCs of the 3D density fields are left for future works. In order to examine the impact of the large-scale environment on the matter clustering, we split the whole simulation volume evenly into 512 small cubes and then measure the local WPS and WCC for different matter species in each sub-volume. To clearly see the dependence of the matter clustering on the density environment, we consider the mean of local WPSs and WCCs over all sub-volumes within each density interval. Our main findings are summarized as follows:

• 1.  
  On the largest scales, the clustering of the total matter shows no dependence on the density environment, whereas it increases monotonically with increasing density on scales up to several times 0.1 hMpc⁻¹. Compared to the global WPS of the total matter, its mean local WPS is most enhanced in the density of ${\bar{\delta }}_{V}\gt 3$ around scales of 3 hMpc⁻¹ where it is dominated by the most massive halos. As the density ${\bar{\delta }}_{V}$ becomes smaller and smaller, the mean local WPS is more and more suppressed around scales of 3 hMpc⁻¹, suggesting that lower density environments contain less massive halos.
• 2.  
  For the pairs of the DM gas, DM stars, and gas stars, the cross correlations in all density environments are close to one on large scales, whereas they deviate from one significantly on small scales. Specifically, we find that the DM and stars are less correlated in less dense environments on small scales. The correlation between the DM and gas remains close to one in all environments on the scales of k ≲ 3 hMpc⁻¹. On smaller scales, this correlation first decreases and then increases with the increasing density. The correlation between the stars and gas also shows a similar density dependence, but is lower than other pairs.
• 3.  
  As revealed by the ratio between the local WPS from the TNG100-1 and that from its companion DM-only simulation, the impact of the AGN feedback on the matter clustering also varies with the density environment. Particularly, AGN feedback has less impact on the matter clustering in the densest environments ${\bar{\delta }}_{V}\gt 3$ and the most underdense environments ${\bar{\delta }}_{V}\lt -0.75$. Compared to these extreme environments, the matter clustering strength at small scales is more suppressed by the AGN feedback in densities of $-0.75\lt {\bar{\delta }}_{V}\lt 3$.

These results are qualitatively consistent with previous research on matter clustering (e.g., van Daalen & Schaye 2015; Fisher & Faltenbacher 2018; Springel et al. 2018; van Daalen et al. 2020; Zhu et al. 2022). This encourages us that we can further apply statistics based on the continuous WT to various cosmological simulations and characterize the environmental dependence of the matter clustering and the baryonic effects on it in more detail.

Y.W. especially thanks Dr. van Milligen for helpful discussions on the statistical error estimation for wavelet statistical quantities. P.H. acknowledges the support of the Natural Science Foundation of Jilin Province, China (No. 20180101228JC), and by the National Science Foundation of China (Nos. 12047569, 12147217). In this work, we used the data from IllustrisTNG simulations. The IllustrisTNG simulations were undertaken with computing time awarded by the Gauss Centre for Supercomputing (GCS) under GCS Large-Scale Projects GCS-ILLU and GCS-DWAR on the GCS share of the supercomputer Hazel Hen at the High Performance Computing Center Stuttgart (HLRS), as well as on the machines of the Max Planck Computing and Data Facility (MPCDF) in Garching, Germany.

In principle, there is a diverse set of wavelets available in the CWT. For instance, the most common ones are the Poisson wavelet, the Mexican hat wavelet, the spline wavelet, the Meyer wavelet, and the Morlet wavelet (see, e.g., Akujuobi 2022 for a review of wavelets). The choice of an appropriate wavelet depends on the goal of the performed analysis. In this study, we chose to use the GDW for two reasons. First, its shape matches well that of the density peak. Second, it provides a good trade-off between spatial and frequency resolution. Subsequently, we will highlight the superiority of the GDW by comparing it with the Poisson and Morlet wavelets.

The Poisson wavelet described by (Holschneider 2000) is

$$\begin{eqnarray}&&{\psi }^{\mathrm{Poi}}(x)={C}^{\mathrm{Poi}}\displaystyle \frac{1-{x}^{2}}{{\left(1+{x}^{2}\right)}^{2}}\end{eqnarray} \tag{ A1 }$$

with Fourier transform

$$\begin{eqnarray}&&{\hat{\psi }}^{\mathrm{Poi}}(k)={C}^{\mathrm{Poi}}| k| {e}^{-| k| },\end{eqnarray} \tag{ A2 }$$

where the constant C^Poi is used to make it have the same energy as the GDW. As can be seen from Figure A1, the Poisson wavelet is more localized in the spatial domain but less localized in the frequency domain than the GDW is. The Morlet wavelet (also called the Gabor wavelet) given by Hong & Kim (2004) is written as

$$\begin{eqnarray}&&{\psi }^{\mathrm{Mor}}(x)={C}^{\mathrm{Mor}}{e}^{-{x}^{2}/(2{\sigma }^{2})}{e}^{i\eta x},\end{eqnarray} \tag{ A3 }$$

which is a harmonic function with frequency η modulated by a Gaussian envelope with variance σ², and its Fourier transform is

$$\begin{eqnarray}&&{\hat{\psi }}^{\mathrm{Mor}}(k)=\sqrt{2\pi }{C}^{\mathrm{Mor}}\sigma {e}^{-\displaystyle \frac{1}{2}{\sigma }^{2}{\left(\eta +k\right)}^{2}},\end{eqnarray} \tag{ A4 }$$

where the constant C^Mor is used to make it have the same energy as the GDW. In principle, the Morlet wavelet is not a true wavelet function, since it does not satisfy the admissibility condition, i.e., ${\int }_{0}^{\infty }| {\hat{\psi }}^{\mathrm{Mor}}(k){| }^{2}/k{\rm{d}}k$ is infinite. This means that the Morlet WT has no inverse transform (Daubechies 1992). Only when ∣σ η∣ is large enough, the admissibility condition is approximately satisfied. A reasonable choice would be ∣σ η∣ ≥ 5 (Hong & Kim 2004), and here we set σ = 5, η = −1. In this case, compared to the GDW and the Poisson wavelet, the Morlet wavelet is the most extended in the spatial domain but the narrowest in the frequency domain, as shown in Figure A1.

Zoom In Zoom Out Reset image size

In order to discern the applicability of the discussed wavelets, we use them individually to calculate the squared CWT of the 1D Zel'dovich density field at θ = 100, which is shown in Figure A2. We observe that in the Poisson case, structures of the density field are resolved quite well along the spatial axis, whereas they are very poor along the scale axis. The Morlet wavelet yields the best separation of scales but the worst separation of spatial features among these wavelets. The GDW performance is intermediate between the Poisson and the Morlet wavelets, giving moderate spatial and scale resolution. Therefore, using GDW to perform local spectral analysis may provide more reliable results.

Zoom In Zoom Out Reset image size

In general, the 3D CWT uses a scale vector containing three scale parameters, i.e., w = (w_x , w_y , w_z ), which measures the scale along each of the three axes (x, y, and z), as pointed out by Wang & He (2021). For simplicity, we consider the isotropic case, in which the same scale parameter w is used for all directions. Then the isotropic CWT of the 3D field is written as

$$\begin{eqnarray}&&{W}_{f}(w,{\bf{r}})={\int }_{{{\mathbb{R}}}^{3}}f({\bf{u}}){\rm{\Psi }}(w,{\bf{r}}-{\bf{u}}){{\rm{d}}}^{3}{\bf{u}},\end{eqnarray} \tag{ B1 }$$

where ${{\mathbb{R}}}^{3}$ represents the 3D real space, and Ψ(w, r) is the isotropic GDW, which is defined by

$$\begin{eqnarray}&&\begin{array}{l}{\rm{\Psi }}(w,{\bf{r}})\equiv {w}^{\kappa }\displaystyle \frac{\partial G(w,{\bf{r}})}{\partial w}\\ =\,\displaystyle \frac{{w}^{\kappa +2}}{16{\pi }^{3/2}}(6-{w}^{2}| {\bf{r}}{| }^{2})\exp \left(-\displaystyle \frac{{w}^{2}| {\bf{r}}{| }^{2}}{4}\right),\end{array}\end{eqnarray} \tag{ B2 }$$

where $G{(w,{\bf{r}})={w}^{3}\exp (-{w}^{2}| {\bf{r}}{| }^{2}/4)/(2\sqrt{\pi })}^{3}$ is the isotropic Gaussian smoothing function, and the index κ is set to be −1/2 so that the integral ∫∣Ψ(w, r)∣²d³ r is constant. The Fourier transform of the isotropic GDW is

$$\begin{eqnarray}&&\hat{{\rm{\Psi }}}(w,{\bf{k}})=2{w}^{\kappa -1}\displaystyle \frac{| {\bf{k}}{| }^{2}}{{w}^{2}}\exp \left(-\displaystyle \frac{| {\bf{k}}{| }^{2}}{{w}^{2}}\right).\end{eqnarray} \tag{ B3 }$$

In line with the method used to obtain Equation (6), combining Equations (B1) and (B2) yields the inverse CWT,

$$\begin{eqnarray}&&f({\bf{r}})=C+{\int }_{0}^{+\infty }{w}^{-\kappa }{W}_{f}(w,{\bf{r}}){\rm{d}}w.\end{eqnarray} \tag{ B4 }$$

Recalling that the 1D wavelet-based statistics are defined through the integral over a length L, the 3D wavelet-based statistics can be defined through the integral over a volume V = L³. The local WPS defined in this way is shown below:

$$\begin{eqnarray}&&{P}_{f,V}^{W}(w)=\displaystyle \frac{1}{V}{\int }_{V}| {W}_{f}(w,{\bf{r}}){| }^{2}{{\rm{d}}}^{3}{\bf{r}},\end{eqnarray} \tag{ B5 }$$

which is the volume average of squared wavelet coefficients within the local volume V, and hence can be written in a more compact form as

$$\begin{eqnarray}&&{P}_{f,V}^{W}(w)=\langle | {W}_{f}(w,{\bf{r}}){| }^{2}{\rangle }_{V}.\end{eqnarray} \tag{ B6 }$$

Then the local WCC between two fields, f(r) and g(r), is given by

$$\begin{eqnarray}&&{C}_{f-g,V}^{W}(w)=\displaystyle \frac{{P}_{f-g,V}^{W}(w)}{\sqrt{{P}_{f,V}^{W}(w){P}_{g,V}^{W}(w)}},\end{eqnarray} \tag{ B7 }$$

where ${P}_{f-g,V}^{W}(w)\equiv \langle {W}_{f}(w,{\bf{r}}){W}_{g}(w,{\bf{r}}){\rangle }_{V}$.

For simulated cosmic fields with a periodic box of volume ${V}_{b}={L}_{b}^{3}$, the global WPS is obtained by replacing the local volume V with the total volume V_b ,

$$\begin{eqnarray}&&{P}_{f}^{W}(w)=\langle | {W}_{f}(w,{\bf{r}}){| }^{2}{\rangle }_{{V}_{b}},\end{eqnarray} \tag{ B8 }$$

and the global WCC is

$$\begin{eqnarray}&&{C}_{f-g}^{W}(w)=\displaystyle \frac{{P}_{f-g}^{W}(w)}{\sqrt{{P}_{f}^{W}(w){P}_{g}^{W}(w)}},\end{eqnarray} \tag{ B9 }$$

where ${P}_{f-g}^{W}(w)\equiv \langle {W}_{f}(w,{\bf{r}}){W}_{g}(w,{\bf{r}}){\rangle }_{{V}_{b}}$. In the 1D case, the global WPS is a smoothed version of the Fourier power spectrum, as stated in Section 2.3. Here we go further and give a tighter relationship between the wavelet and the Fourier power spectra. According to Parseval's theorem, we have

$$\begin{eqnarray}\begin{array}{rcl}{P}_{f}^{W}(w) & = & \langle | {W}_{f}(w,{\bf{r}}){| }^{2}{\rangle }_{{V}_{b}}\\ & = & \displaystyle \frac{1}{{V}_{b}}{\displaystyle \int }_{{V}_{b}}| {W}_{f}(w,{\bf{r}}){| }^{2}{{\rm{d}}}^{3}{\bf{r}}\\ & = & \displaystyle \frac{1}{{\left(2\pi \right)}^{3}{V}_{b}}\displaystyle \int | \hat{f}({\bf{k}}){| }^{2}| \hat{{\rm{\Psi }}}(w,k){| }^{2}{k}^{2}{\rm{d}}k\sin \theta {\rm{d}}\theta {\rm{d}}\phi ,\end{array}\end{eqnarray} \tag{ B10 }$$

in which we exploit the fact that GDW is isotropic, namely, $\hat{{\rm{\Psi }}}(w,k)=\hat{{\rm{\Psi }}}(w,{\bf{k}})$. If $\hat{f}({\bf{k}})$ is also isotropic, then the above equation has a more simple form,

$$\begin{eqnarray}&&{P}_{f}^{W}(w)=\displaystyle \frac{1}{2{\pi }^{2}}{\int }_{0}^{+\infty }{P}_{f}^{F}(k)| \hat{{\rm{\Psi }}}(w,k){| }^{2}{k}^{2}{\rm{d}}k,\end{eqnarray} \tag{ B11 }$$

where ${P}_{f}^{F}(k)=\tfrac{| \hat{f}(k){| }^{2}}{V}$ is the Fourier power spectrum. Integrating the scale factor w at both ends of this equation yields

$$\begin{eqnarray}&&\begin{array}{l}{\displaystyle \int }_{0}^{\infty }{P}_{f}^{W}(w){\rm{d}}w=\displaystyle \frac{1}{2{\pi }^{2}}\displaystyle \int {P}_{f}^{F}(k)\left({\displaystyle \int }_{0}^{\infty }| \hat{{\rm{\Psi }}}(w,k){| }^{2}{\rm{d}}w\right){k}^{2}{\rm{d}}k\\ \,=\,\displaystyle \frac{1}{4{\pi }^{2}}{\displaystyle \int }_{0}^{+\infty }{P}_{f}^{F}(k){\rm{d}}k,\end{array}\end{eqnarray} \tag{ B12 }$$

in which ${\int }_{0}^{\infty }| \hat{{\rm{\Psi }}}(w,k){| }^{2}{\rm{d}}w=1/2{k}^{2}$. Finally, with the help of the relation $w=2k/\sqrt{7}$ obtained by following the method in Section 2.5, we find that the global wavelet and Fourier power spectra differ by only one constant factor, as shown below:

$$\begin{eqnarray}&&\displaystyle \frac{8{{\pi }}^{2}}{\sqrt{7}}{P}_{f}^{W}(k)={P}_{f}^{F}(k).\end{eqnarray} \tag{ B13 }$$