INTEGRATED SACHS–WOLFE IMPRINT OF SUPERSTRUCTURES ON LINEAR SCALES

According to Pápai & Szapudi (2010), if the average density contrast in a Gaussian random field inside a sphere of radius R is given, we expect the density contrast at radius r from its center to be

$$\begin{equation} \langle \delta (r)\rangle _{\delta _{{\rm in}}(R)} = \frac{\langle \delta (r) \delta _{{\rm in}}(R)\rangle }{\langle \delta _{{\rm in}}^2(R)\rangle }\delta _{{\rm in}}(R), \end{equation} \tag{ 1 }$$

where 〈...〉 stands for ensemble averaging, δ(r) is the density contrast at radius r, and δ_in(R) is the average of δ inside radius R. In the Appendix of Dekel (1981) the author derived a similar formula for a single, point-like location with given density, which is the R → 0 limit. Equation (1) is essentially the two-point function multiplied by a normalization constant. A more general calculation for the density profile around local maxima in a Gaussian random field was carried out by Bardeen et al. (1986).

In order to measure $\langle \delta (r)\rangle _{\delta _{{\rm in}}(R)}$, one needs to select regions with δ_in(R) and calculate the average over these. However, according to Equation (1), $\frac{\delta (r)}{\delta _{{\rm in}}(R)}$ gives an unbiased estimator of $\frac{\langle \delta (r) \delta _{{\rm in}}(R)\rangle }{\langle \delta _{{\rm in}}^2(R)\rangle }$, the shape of the profile:

$$\begin{eqnarray} \left\langle\frac{\delta }{\delta _{{\rm in}}}\right\rangle &=& \int\! {{d}\delta {d}\delta _{{\rm in}} \frac{\delta }{\delta _{{\rm in}}}P(\delta |\delta _{{\rm in}})P(\delta _{{\rm in}}) } \nonumber\\ &=& \frac{\langle \delta \delta _{{\rm in}}\rangle }{\langle \delta _{{\rm in}}^2\rangle } \int { {d}\delta _{{\rm in}} P(\delta _{{\rm in}}) } = \frac{\langle \delta \delta _{{\rm in}}\rangle }{\langle \delta _{{\rm in}}^2\rangle }, \end{eqnarray} \tag{ 2 }$$

where we shorten δ(r) and δ_in(R) as δ and δ_in. We also use the fact that by definition $\langle \delta (r)\rangle _{\delta _{{\rm in}}(R)} = \int {{d}\delta \delta P(\delta |\delta _{{\rm in}}) }$ and $\frac{\langle \delta \delta _{{\rm in}}\rangle }{\langle \delta _{{\rm in}}^2\rangle }$ is independent of δ_in. The argument above shows that it is needless to search for spheres with a certain δ_in(R) to verify Equation (1). In the rest of this subsection we deal with the integral of δ(r) instead of δ(r), because the former is measurable directly

$$\begin{equation} M(r) = 4\pi \int _0^r {d}r r^2\delta (r). \end{equation} \tag{ 3 }$$

The linear size of the HVS is 3000 h⁻¹ Mpc. We chose to place spheres at 25³ evenly spaced coordinates. We counted galaxies in concentric spheres with radii of $10{\,{h^{-1}\;\rm Mpc}}, 20{\,{h^{-1}\;\rm Mpc}}, \ldots, 160{\,{h^{-1}\;\rm Mpc}}$ and subtracted the average to get M(r). We normalized this as $\frac{M(r)}{M(R)}\frac{4\pi }{3}R^3$ and calculated the average of the 25³ instances to get an estimate of the integral of Equation (2). As R changes, technically, by weighting M(r) differently, one can get the profile of linear superstructures of different sizes. A nuisance we have to deal with is that the estimator given by Equation (2) is unstable because δ_in can be arbitrarily small. The ratio of two Gaussian variables with zero means has a Lorentzian distribution, which has an undefined variance. This problem can be avoided if cases when |δ_in| is below a threshold are ignored. This constraint leaves the estimator unchanged as can be seen from Equation (2). We chose the threshold based on the variance of δ_in in the particular catalog, e.g., $2\sqrt{\langle \delta _{{\rm in}}^{2}\rangle }$. The result is not sensitive to the exact choice but the right choice can decrease the variance of the estimator substantially. In Figure 1 we plotted the measured profiles and their theoretical counterparts for R = 60 h⁻¹ Mpc on the upper panel and for R = 100 h⁻¹ Mpc on the lower panel. We obtained the error bars by repeating the measurement on 100 mock catalogs generated with a second-order Lagrangian (2LPT) code (Crocce et al. 2006). The mocks were set up to have the same cosmological parameters, size, and density as the HVS. To calculate the theoretical profiles of Equation (1), we used a matter power spectrum calculated via CAMB (Lewis et al. 2000). The measured profiles, both in 2LPT and the HVS, appear to be in good agreement with theory within the uncertainty up to 400 h⁻¹ Mpc, the largest radius in our measurement.

Zoom In Zoom Out Reset image size

In Pápai & Szapudi (2010) the ISW effect in the direction of the center of a superstructure was calculated. This was done simply by calculating the spherically symmetric potential of the superstructure, then computing the following integral along the path of the CMB photon (Sachs & Wolfe 1967):

$$\begin{equation} \Phi (r) = {-}\frac{3\Omega _m}{8\pi }\bigg (\frac{H_0}{c} \bigg)^2 \int _{r}^{\infty } \frac{M(\tilde{r})}{\tilde{r}^2}{d}\tilde{r}, \end{equation} \tag{ 4 }$$

$$\begin{equation} \frac{\Delta T}{T} = {-}\frac{2}{c^2}\int \! {d}\tau \frac{\partial \Phi (r(\tau),\tau) }{\partial \tau }, \end{equation} \tag{ 5 }$$

where τ denotes conformal time. Similarly, this calculation can be carried out for any direction easily in order to predict the full ISW profile. In linear theory the potential at a particular redshift is just $\Phi (r,z) = \Phi (r,z=0)\frac{D(z)}{a(z)}$. We leave the study of nonlinear corrections for future work. Since these integrals are linear, the linear ISW signal is predicted to be proportional to the average density of the superstructure, δ_in. In this paper, we place the centers of the superstructures at z = 0.52, which is the median redshift of the galaxy catalog we use in Section 3.

To test the ISW profile we traced rays through the HVS along the z-axis. We calculated δ on an 800³ grid and calculated Φ at z = 0 by using the simple formula

$$\begin{equation} \Phi (k) = -\frac{3\Omega _m}{2}\bigg (\frac{H_0}{c} \bigg)^2\frac{\delta (k)}{k^2}. \end{equation} \tag{ 6 }$$

Since we aimed to analyze imprints of superstructures at redshift 0.52, we defined the boundaries of the integral of Equation (5) in a way that put the 25³ preselected locations (see Section 2.1) at this redshift. We used periodic boundary conditions and integrated up to 3000 h⁻¹ Mpc from starting points ensuring that 25² of the locations were at redshift 0.52 each time. This yielded 25 ISW images of the HVS in the x–y plane. Each of these images contains imprints of 25² superstructures. After this, one can follow the procedure discussed for density profiles in Section 2.1. If the ISW profile of a superstructure is normalized by the average density of the superstructure, the result will be an estimator of the shape of the ISW profile:

$$\begin{equation} \left\langle \frac{\Delta T(r_{{\rm 2D}})}{\delta _{{\rm in}}(R)}\right\rangle = f(r_{{\rm 2D}},R). \end{equation} \tag{ 7 }$$

The arguments are r_2D, the radius in the x–y plane, and R, the radius in which we know the galaxy count in the corresponding volume. The interpretation of Equation (7) is similar to that of Equation (2). The density fluctuations are replaced by the temperature on the left-hand side of the equation; while on the right-hand side, we symbolically refer to the normalized ISW profile which can be derived via Equations (3)–(5).

Here we would like to remind the reader that 〈...〉 refers to ensemble averaging. For density, averaging over a relatively large volume proved to be a sufficient substitute in Section 2.1. Despite the fact that the cosmic variance is large on the largest scales, these modes contribute only a little to the density distribution since their amplitudes are small. (The power spectrum goes to zero with k.) For the potential, this picture changes due to the 1/k² factor. As a result, the ISW effect has a large cosmic variance regardless of the large size of the simulation. To demonstrate this, we projected the density on the x–y plane. We split this projection in the middle of the x-axis to create two 1500 × 3000 Mpc² h⁻² areas. In Figure 2 in the top row we plotted the histograms of these. We did the same with one of the 25 ISW images, which is essentially the projection of the potential. Those histograms are shown in the middle row. The cosmic variance is significantly larger for the ISW map. To demonstrate that the difference is due to modes on the largest scales we removed the low k-modes from the potential. After experimenting we found that removing modes with $|k|<2\pi \frac{5}{L}$, where L = 3000 h⁻¹ Mpc for the HVS, gave a visually compelling result (Figure 2, bottom row). After applying this high-pass filter to a two-dimensional projection, the effect of the cosmic variance on the modes with the lowest wavenumber will be around 18%. This argument is valid for any survey or simulation. Generally, modes can be combined optimally with the proper covariance matrix. In this paper, we simply use the sum of the modes as an estimator of the ISW signal. In other words, our covariance matrix is diagonal with a truncation at 1/5 of the linear size. This subjective choice, which we consistently carry along through the paper, is based purely on the histogram of the potential (Figure 2).

Zoom In Zoom Out Reset image size

Having removed the low k-modes from both the density and the potential, we computed the average in Equation (7). As before, we omitted cases when the denominator was too close to zero. (See Section 2.1.) The results for R = 60 h⁻¹ Mpc and R = 100 h⁻¹ Mpc are shown in Figure 3. The error bars show the uncertainty calculated from the HVS itself due to the long CPU time required to calculate the ISW images. Predictions of the linear model are calculated by removing the low k-modes from the power spectrum when doing calculation according to Equation (1).

Zoom In Zoom Out Reset image size

This approach has to be modified slightly when the CMB is given and late time anisotropies are affected by every mode of the galaxy distribution. The previous exercise still proves that low k-modes are to be ignored if we want to study the imprints of superstructures statistically. This can be achieved by filtering out these modes from the CMB before performing such analysis. In addition, one has to remove the same components from the theoretical ISW profile before fitting the data. This is shown in Figure 4. We filtered the ISW map and the ISW profile, which is implicitly expressed as f(r_2D, R), on the right-hand side of Equation (7), so that only modes with $|k|>2\pi \frac{5}{L}$ remained. After these manipulations the measurement remained consistent with linear theory.

Zoom In Zoom Out Reset image size

The success of these tests proves the validity of the Gaussian model of Pápai & Szapudi (2010). We studied Euclidean ISW projections because or their simplicity but there is no reason to think that our findings are restricted to the Euclidean case. The only major difference is that on a sphere one deals with spherical harmonics and instead of low k-modes we speak of low l-modes analogously.

Large regions with significantly high overdensities or low underdensities are called here superclusters or supervoids. In a given a data set these regions can be defined in many different ways. Void finder algorithms were compared in Colberg et al. (2008). Each implementation is based on a certain notion of a void. Voids can have spherical or irregular shapes, finders can be parameter free or there can be a hard coded density threshold or radius to describe what a void is like. It is important to emphasize that despite the similarities, any analysis using the output of a void finder will be sensitive to its particulars. This renders our study qualitative in nature.

Our work is based on supervoids and superclusters found in a subsample of the SDSS DR6 LRG catalog by ZOBOV and VOBOZ, parameter free void and cluster finder algorithms (Neyrinck 2008; Neyrinck et al. 2005). ZOBOV and VOBOZ calculate the Voronoi tessellation of the data in order to approximate the density field with a simple function, a function that is constant inside polyhedra. Voids are found by using the concept of drainage basins from geography. Clusters are identified analogously.

A table with the properties of 50 supervoids and 50 superclusters can be found in Granett et al. (2008b). Their sizes and densities cannot be quantified simply due to their shapes and survey boundary effects. The bias, redshift error, and the shot noise are also sources of uncertainty for both size and density. This is why we used only the directions of the centers of these superstructures in creating ISW maps in Section 3.3. The error in the angular coordinates is relatively small.

We used the Internal Linear Combination (ILC) Map from the WMAP 7 year release (Gold et al. 2011) with the KQ75y7 galactic foreground and point-source mask. We employed Healpix pixelization at resolution n_side = 64 (Górski et al. 2005). In general, we strove to carry out the measurement with up-to-date maps and masks, while remaining consistent with Granett et al. (2009).

Section 2 concluded that a Gaussian model gave a good approximation to the expected ISW signal from a large spherical region in the dark matter field with given density contrast. It also became apparent that in an analysis of statistical ISW maps low l-modes were to be filtered out due to the large cosmic variance of the potential affecting these modes. This also lowers the cosmic variance. The lowest mode to be taken into account is determined by the sky coverage of the ISW map. Experimenting with the HVS gave us $k_{{\rm min}}\approx 2\pi \frac{5}{L}$ for a Euclidean ISW map, where L is the linear size of the map. This translates into $l_{{\rm lmin}}\approx \frac{2\pi 5r}{\lambda }$ via the Limber approximation, where λ and r are the angular size and the median distance of the SDSS DR6 LRG sample, respectively. A rough estimate yields l_min ≈ 12.

Due to the large uncertainty in their size and density, we assume that each supervoid is a realization of the same statistical entity when we create ISW maps. We consider superclusters similar in nature but with an opposite density contrast.

An ISW map is made simply by placing the expected ISW profiles of superstructures at the locations given in Granett et al. (2008b). We set δ_in to 1 for clusters and −1 for voids in formulae of Section 2. For every radius there is a corresponding ISW map. We calculated the a_lm coefficients of the spherical expansion of both the ISW maps and the CMB map, and zeroed the modes below l_min = 12 before transforming the a_lm's back. These manipulations were carried out with routines from the Healpix package (Górski et al. 2005). Figure 5 shows how the ISW profile changes when modes are removed. This figure also tells about the decrease of the signal. In the case of l_min = 12, it is around 50% at the center and grows outward. Our heuristic method aims to find an acceptable compromise between the loss of signal and small sample variance.

Zoom In Zoom Out Reset image size

We continue with measuring the amplitude of the ISW maps created in Section 3.2 in the WMAP ILC map. To do this we used the same matched filter technique as Granett et al. (2009). A likelihood function is defined as

$$\begin{equation} L(\lambda) = -\frac{1}{2}(T_{{\rm cmb}}-\lambda T_{{\rm ISW}})C^{-1}(T_{{\rm cmb}}-\lambda T_{{\rm ISW}}), \end{equation} \tag{ 8 }$$

where T_cmb is a foreground cleaned CMB map, in this case the WMAP ILC map, and T_ISW refers to the ISW maps. The pixel covariance matrix, C, is computed from the best-fitting ΛCDM power spectrum. This likelihood is based on the assumption that the T_ISW and the primary CMB are the realizations of two very different Gaussian random fields and the latter is given by the minimum of Equation (8):

$$\begin{equation} \lambda ^* = \frac{T_{{\rm ISW}}C^{-1}T_{{\rm cmb}}}{T_{{\rm ISW}}C^{-1}T_{{\rm ISW}}}, \end{equation} \tag{ 9 }$$

$$\begin{equation} T_{{\rm prim}} = T_{{\rm cmb}} - \lambda ^* T_{{\rm ISW}}. \end{equation} \tag{ 10 }$$

As explained in Section 3.2, we filtered out modes up to l_min = 12 from both T_ISW and T_cmb, and before building the covariance matrix from C_l's we set C₀ = C₁ = ⋅⋅⋅ = C₁₁ = 0.

Since the exact inverse of the covariance matrix is unstable for partial sky maps, we regularized C by calculating a pseudo-inverse omitting the noisiest modes. The likelihood function above defines an error on λ* as

$$\begin{equation} \Delta \lambda ^2 = \frac{1}{T_{{\rm ISW}}C^{-1}T_{{\rm ISW}}}. \end{equation} \tag{ 11 }$$

We calculated λ* according to Equation (9). Since we have ISW maps for different void/cluster radii, the result is a function, λ(R). This is presented in Figure 6 in the bottom row on the left panel. The error bars are given by Equation (11). From now on every time we refer to a λ we mean the value at the minimum of the likelihood function, so we drop the superscript *. If we marginalize λ(R)/Δλ(R) over R, we get a 3.24σ detection. Similarly, the marginalized value of R is 55 ± 28 h⁻¹ Mpc.

Zoom In Zoom Out Reset image size

λ(R) can be interpreted as the average density contrast in 100 superstructures inside radius R, where the supervoids are taken into account with a minus sign:

$$\begin{equation} \lambda (R) \approx \frac{1}{2}\big (\delta _{{\rm in}}^{50c}(R) - \delta _{{\rm in}}^{50v}(R)\big), \end{equation} \tag{ 12 }$$

where superscripts 50c and 50v refer to the average of 50 superclusters and 50 supervoids. This interpretation originates from the linear relationship between the ISW signal and the density fluctuations, and our choice of δ_in = ±1 when building ISW maps in Section 3.2. In addition to the data curves, in Figure 6 we plotted the prediction of linear theory for the average of the 50 largest superclusters. We assumed a Gaussian probability distribution function (PDF) which is reasonable in light of the findings of Pápai & Szapudi (2010). The variance of the PDF was calculated by CAMB (Lewis et al. 2000) using the latest WMAP cosmological parameters (Jarosik et al. 2011) and it was scaled with the growth function to the median redshift of the superstructures (z = 0.52). The number of independent clusters was estimated as V_survey/V_cluster, where V_survey is the volume of the LRG sample used by Granett et al. (2009). λ is consistently several times larger than the predicted δ_in. As a consistency test we split the voids and clusters into two groups. Group 1 consists of the first 25 voids and 25 clusters with the highest significance in the supplementary tables of Granett et al. (2008b). The other 25 voids and 25 clusters make up group 2. We created two ISW maps, one from each of the groups. After repeating the whole procedure for these, we got the center and the right panels in the bottom row of Figure 6. As expected, superstructures associated with lower significance produce lower signal on the CMB and allow larger deviation from the average λ. This supports the hypothesis that the signal is due to the ISW effect.

We also checked the sensitivity of the results to the choice of l_min, the lowest spherical harmonics included in the analysis. The result for different choices of l_min is shown in the middle and top rows of Figure 6. The top row shows λ(R, l_min), while the middle row shows σ(R, l_min) = λ(R, l_min)/Δλ(R, l_min) on contour plots. The result seems robust for a wide range of l_mins.

We performed another consistency check to acquire a better understanding of the uncertainty of λ and its relationship with δ_in. We created ISW maps from single superstructures. We used the matched filter technique to determine a λ and its error as before. The average λ for voids and for clusters separately are plotted in Figure 7. In this case we defined the uncertainty of λ as the standard deviation of its mean. We also plotted the previous results for the previous combined ISW map. All of this is for l_min = 12. The three measured curves are consistent with each other, the error bars show the 2σ deviation. The single void and cluster fit allow wider range of λ than the combined fit.

Zoom In Zoom Out Reset image size