A&A 449, 891-902 (2006)
DOI: 10.1051/0004-6361:20053939

Acoustic oscillations in the SDSS DR4 luminous red galaxy sample power spectrum[*]

G. Hütsi1,2


1 - Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str. 1, 86740 Garching bei München, Germany
2 - Tartu Observatory, Tõravere 61602, Estonia

Received 29 July 2005 / Accepted 1 December 2005

Abstract
We calculate the redshift-space power spectrum of the Sloan Digital Sky Survey (SDSS) Data Release 4 (DR4) Luminous Red Galaxy (LRG) sample, finding evidence for a full series of acoustic features down to the scales of $\sim $ $0.2~h~{\rm Mpc}^{-1}$. This corresponds to the 7th peak in the CMB angular power spectrum. The acoustic scale derived, $(105.4 \pm 2.3)~h^{-1}~{\rm Mpc}$, agrees very well with the "concordance'' model prediction and also with the one determined via the analysis of the spatial two-point correlation function by Eisenstein et al. (2005, ApJ, 633, 560). The models with baryonic features are favored by $3.3 \sigma$ over their "smoothed'' counterparts without any oscillatory behavior. This is not only an independent confirmation of results by Eisenstein et al. (2005), made with different methods and software but also, to our knowledge, is the first determination of the power spectrum of the SDSS LRG sample.

Key words: large-scale structure of Universe

1 Introduction

At the beginning of the 1970's it was already realized that acoustic waves in the tightly coupled baryon-photon fluid prior to the epoch of recombination will lead to characteristic maxima and minima in the post-recombination matter power spectrum. The same mechanism is also responsible for the prominent peak structure in the CMB angular power spectrum (Peebles & Yu 1970; Doroshkevich et al. 1978; Sunyaev & Zeldovich 1970). The scale of these features reflects the size of the sound horizon, which itself is fully determined given the physical densities $\Omega_{\rm b} h^2$ and $\Omega_{\rm m} h^2$. The acoustic horizon can be calibrated using the CMB data, thus turning it into a standard ruler which can be used to carry out classical cosmological tests. For example, if we are able to measure the redshift and angular intervals corresponding to the physically known acoustic scale in the matter power spectrum at a range of redshifts, we can immediately find the angular diameter distance $d_{\rm A}$ and Hubble parameter H as a function of redshift. Having good knowledge of these dependencies allows us to place constraints on the properties of the dark energy. To carry out this project one needs a tracer population of objects whose clustering properties with respect to the underlying matter distribution are reasonably well understood. There have been several works discussing the use of galaxies (Hu & Haiman 2003; Blake & Glazebrook 2003; Linder 2003; Seo & Eisenstein 2003) and clusters of galaxies (Hu & Haiman 2003; Hütsi 2005; Majumdar & Mohr 2004) for this purpose. What is most important is that already currently existing galaxy redshift surveys have lead to the detection of acoustic features in the spatial distribution of galaxies, thus providing clear support for the standard gravitational instability picture of cosmic structure formation. In Eisenstein et al. (2005) the detection of the acoustic "bump'' in the two-point redshift-space correlation function of the SDSS[*] LRG sample was announced. The discovery of similar features in the power spectrum of 2dF[*] galaxies is presented in Cole et al. (2005). These results clearly demonstrate the promise of future dedicated galaxy redshift surveys like KAOS[*] Similarly, useful measurements of the acoustic scale may be obtained by the planned SZ cluster surveys like those of the PLANCK Surveyor[*] spacecraft and SPT[*] (Hütsi 2005) and also with future large photometric redshift surveys (Blake & Bridle 2005). For the SZ surveys one needs an additional optical follow-up to obtain estimates for the cluster redshifts. In this paper we calculate the redshift-space power spectrum of the SDSS LRG sample finding evidence for acoustic oscillations down to the scales of $\sim $ $0.2~h~{\rm Mpc}^{-1}$, which effectively correspond to the 7th peak in the CMB angular power spectrum. These scales in the CMB are very strongly damped due to the finite width of the last-scattering surface and also due to the Silk damping (Silk 1968). This can be seen in Fig. 1[*] where the CMB data is plotted in a somewhat unusual way to enhance the acoustic features at the high wavenumber damping tail. Also, at those scales the secondary CMB anisotropies (mostly the thermal Sunyaev-Zeldovich effect; Sunyaev & Zeldovich 1980,1972) start to dominate over the primary signal. On the other hand, features in the matter power spectrum, although small ($\sim $$5\%$ fluctuations), are preserved by the linear evolution and so allow probe of acoustic phenomena at scales smaller than the ones accessible by CMB studies.

The paper is structured as follows. In Sect. 2 we describe the dataset to be analyzed. Section 3 presents the method of the power spectrum calculation. In Sect. 4 we determine power spectrum errors and the covariance matrix. Section 5 discusses the convolution effect of the survey window. Analytical model spectra are presented in Sect. 6. The results of the measurement of the acoustic scale are given in Sect. 7. A correlation function analysis is carried out in Sect. 8. In Sect. 9 we compare the measured power spectrum with the published results for the 2dF and SDSS main sample, and conclude in Sect. 10.


  \begin{figure}
\par\includegraphics[width=7.7cm,clip]{3939_f01.eps}
\end{figure} Figure 1: Acoustic oscillations in the CMB ( upper panel) and linear matter power spectrum ( lower panel) for the "concordance'' cosmological model. Here, as we have plotted the spectra against spatial wavenumber k, we have changed the standard notation of $C_\ell $ to Ck. Due to the k3 factor the first CMB acoustic peak is barely visible. Density fluctuations in matter at smaller scales, being mostly induced by the velocity fields, are out of phase with respect to the fluctuations in the CMB component. The fluctuation period is also twice as large.
Open with DEXTER


  \begin{figure}
\par\includegraphics[width=7.5cm,clip]{3939_f02.eps}
\end{figure} Figure 2: Comoving number density of galaxies as a function of comoving distance. Smooth solid line shows a cubic spline fit to the number density estimated for 50 discrete radial bins.
Open with DEXTER

2 Data

We analyze the publicly available data from the SDSS DR4 (Adelman-McCarthy et al. 2005). Specifically, we carry out our power spectrum measurements using the subset of the SDSS spectroscopic sample known as the Luminous Red Galaxy (LRG) sample. The LRG selection algorithm (Eisenstein et al. 2001) selects $\sim $12 galaxies per square degree meeting specific color and magnitude criteria[*]. The resulting set of galaxies consists mostly of early types populating dense cluster environments and as such are significantly biased (bias factor $b \sim 2$) with respect to the underlying matter distribution. The selection method is very effective in producing a galaxy sample with reasonably high density up to a redshift of $z \sim 0.5$.

Since the selection criteria are very complicated, involving both cuts in magnitude and in color, and also due to the steepness of the luminosity function, the usual method of using only the luminosity function to determine radial selection function does not work here (Zehavi et al. 2005). Rather we build the radial selection function as a smooth spline fit to the number density profiles shown in Fig. 2. To calculate distances we use the cosmological parameters as given by the WMAP[*] "concordance'' model (Spergel et al. 2003). Unfortunately the coverage masks of the SDSS DR4 spectroscopic sample are not available in a readily accessible format and so we chose to build the angular survey masks using the galaxy data itself[*]. As the number density of galaxies in the sample is rather high, one can determine relatively accurately the beginning, end and also possible gaps in the scan stripes. We have built angular masks using both the whole DR4 galaxy sample and LRGs only. The measured power spectra are practically identical with only some minor differences at smaller scales (see Fig. 6). This can be seen as an indication that our power spectrum measurements are stable against small uncertainties in the survey geometry. The angular distribution of the galaxies and also the boundaries of the survey mask built in the above-mentioned way (using all the galaxies) is shown in Fig. 3. Here the angular positions are plotted using the survey coordinate system of the SDSS[*].


  \begin{figure}
\par\includegraphics[width=7.5cm,clip]{3939_f03.eps}
\end{figure} Figure 3: Angular distribution of galaxies given in the SDSS survey coordinates  $(\lambda ,\eta )$. The survey mask is shown with solid lines. The vertical dashed lines show the division of the sample into 22 separate regions each containing $\sim $2350 galaxies. This division will be exploited in the "jackknife'' error analysis of the correlation function.
Open with DEXTER

We have selected all the objects that have a spectrum classified as galaxy (i.e. SpecClass=2) and are additionally flagged as GALAXY_RED or GALAXY_RED_II (i.e. PrimTarget bit mask set as 0x20 or 0x4000000, respectively). Only galaxies for which the redshift confidence parameter, zConf, is greater than 0.95 were used. We apply lower and upper redshift cutoffs of 0.16 and 0.47 as also done in Eisenstein et al. (2005). The lower cutoff is needed since the color cuts that define the LRG sample break down for redshifts below $\sim $0.2 (Eisenstein et al. 2001). For the analysis presented in this paper we have excluded the three southern stripes since these just increase the sidelobes of the survey window without adding much extra volume. We have also removed some minor parts of the sample to obtain a more continuous and smoother chunk of volume. In total the analyzed galaxy sample covers $\sim $ $0.75 ~h^{-3}~{\rm Gpc}^3$ over $\sim $3850 square degrees on the sky and contains 51 763 galaxies.

3 Power spectrum calculation

We calculate the power spectrum using a direct Fourier method as described in Feldman et al. (1994) (FKP). Power spectra determined this way are the so-called pseudospectra, meaning that the estimates derived are convolved with a survey window. Since in the case of the analyzed LRG sample the volume covered is very large, reaching  $0.75 ~h^{-3}~{\rm Gpc}^3$, and also the survey volume has relatively large dimensions along all perpendicular directions, the correlations in the Fourier space are rather compact. On intermediate scales and if the power spectrum binning is chosen wide enough, the FKP estimator gives a good approximation to the true underlying power.

The FKP estimate for a 3D pseudospectrum is:

 \begin{displaymath}
\tilde{P}(\vec{k}) = \vert F(\vec{k})\vert^2 - P_{\rm shot} ,
\end{displaymath} (1)

where

\begin{displaymath}F(\vec{k}) = \int {\rm d}^3r ~ F(\vec{r})\exp({\rm i}\vec{k}\cdot\vec{r}) .
\end{displaymath} (2)

Here $F(\vec{r})$ is the weighted density contrast field:

 \begin{displaymath}
F(\vec{r}) = w(\vec{r})\left[n_{\rm g}(\vec{r})-\alpha n_{\rm s}(\vec{r})\right] .
\end{displaymath} (3)

$n_{\rm g}(\vec{r})$ and $n_{\rm s}(\vec{r})$ denote the number densities of the analyzed galaxy catalog and a synthetic random catalog with the same selection criteria, respectively. Since we are dealing with discrete point processes, densities can be given as:
$\displaystyle n_{\rm g}(\vec{r}) = \sum \limits_i \delta^D(\vec{r}-\vec{r}^{\rm g}_i) ,$     (4)
$\displaystyle n_{\rm s}(\vec{r}) = \sum \limits_i \delta^D(\vec{r}-\vec{r}^{\rm s}_i) ,$     (5)

where $\vec{r}^{\rm g}_i$ and $\vec{r}^{\rm s}_i$ denote the location of the i-th point in the real and synthetic catalog, respectively, and $\delta^D$ is the 3D Dirac delta function. $\alpha$ in Eq. (3) is the ratio of the number of galaxies to the number of random points in the synthetic catalog i.e. $\alpha=\frac{N_{\rm g}}{N_{\rm s}}$. In our calculations we have $N_{\rm s}=10^7$ and thus $\alpha \simeq 0.0052$. For the weight function $w(\vec{r})$ there are three choices often used in the literature:

\begin{displaymath}w(\vec{r})\propto \left\{
\begin{array}{lll}
\frac{1}{\bar{n...
... {\quad for\ an\ optimal\ FKP\ weighting.}
\end{array} \right.
\end{displaymath} (6)

Here $\bar{n}(\vec{r})$ is the average number density of galaxies at the comoving location $\vec{r}$ i.e. the radial selection function of the survey (see Fig. 2) times the angular mask (Fig. 3). In our calculations we use an optimal FKP weighting scheme, although pure volume weighting would give practically the same results, especially on the larger scales ( $k \lesssim 0.09~h~{\rm Mpc}^{-1}$), since then for the majority of the sample $\bar{n}(\vec{r})\tilde{P} \sim 3$ [*]. The weights in Eq. (3) are normalized such that

\begin{displaymath}\int {\rm d}^3r ~\bar{n}^2(\vec{r}) w^2(\vec{r}) = 1 ,
\end{displaymath} (7)

which can be approximated as the following sum over the synthetic catalog[*]:

\begin{displaymath}\alpha \sum \limits_i \bar{n}(r^{\rm s}_i)w^2(r^{\rm s}_i) = 1 .
\end{displaymath} (8)

The last term in Eq. (1) represents the Poissonian discreteness noise and can be expressed as:

\begin{displaymath}P_{\rm shot} = (1+\alpha) \int {\rm d}^3r ~\bar{n}(\vec{r}) w...
...c{r}) \simeq \alpha(1+\alpha)\sum \limits_i w^2(r^{\rm s}_i) .
\end{displaymath} (9)

Since we are using Fast Fourier Transforms (FFTs) to speed up the calculation of the Fourier sums, we have to deal with some extra complications. As the density field is now restricted on a regular grid with a finite cell size, we have to correct for the smoothing effect this has caused. Also, if our underlying density field contains spatial modes with higher frequency than our grid's Nyquist frequency, $k_{{\rm Ny}}$, then these will be "folded back'' into the frequency interval the grid can support, increasing power close to $k_{{\rm Ny}}$ - the so-called aliasing effect. The relation between the spectra calculated using direct summation and the ones found using FFT techniques was worked out by Jing (2005). It can be expressed as follows:
 
$\displaystyle \vert F(\vec{k})\vert _{{\rm FFT}}^2=\sum \limits_{\vec{n} \in \m...
...\vec{n} \in \mathbb{Z}} \vert\mathcal{W}(\vec{k}+2k_{{\rm Ny}}\vec{n})\vert^2 ,$     (10)

where $\mathcal{W}(\vec{k})$ is the mass assignment function used to build a density grid out of the point set. We use the Triangular Shaped Cloud (TSC) assignment method (Hockney & Eastwood 1988). Since the TSC filter can be obtained by convolving the uniform cube (the Nearest Grid Point filter) two times with itself, the Fourier representation of it follows immediately:

\begin{displaymath}\mathcal{W}(\vec{k}) = \left[\frac{\prod \limits_{i=1}^3 \sin...
...{{\rm Ny}}}\right)}\right]^3, \quad \vec{k} = (k_1,k_2,k_3) .
\end{displaymath} (11)

Here the sum that represents the contribution from aliases runs over all the integer vectors $\vec{n}$. Equation (10) is the direct analog of the previous Eq. (1). The convolution with the mass assignment filter has introduced  $\mathcal{W}^2(\vec{k})$ factors both in the power spectrum and the shot noise term. The sum in the last term of Eq. (10) can be performed analytically for the TSC filter to yield the result (Jing 2005)[*]:
$\displaystyle \sum \limits_{\vec{n} \in \mathbb{Z}} \vert\mathcal{W}(\vec{k}+2k...
...ght)
+ \frac{2}{15}\sin^4\left(\frac{\pi k_i}{2k_{{\rm Ny}}}\right)\right]\cdot$     (12)

To recover the angle averaged pseudospectrum $\tilde{P}(k)$ from Eq. (10) we use an iterative scheme as described in Jing (2005) with a slight modification: we do not approximate the small scale spectrum by a simple power law, but also allow for the possible running of the spectral index i.e. the parametric shape of the power spectrum is taken to be a parabola in log-log. Since on small scales the power spectrum drops fast, the sum over $\vec{n}$ in Eq. (10) converges rapidly. In calculations we use only integer vectors with $\vert\vec{n}\vert \leq 5$. The angular average is taken over all the vectors $\vec{k}$ lying in the same k-space shell with width $\Delta k$. The resulting $\tilde{P}$ is taken to be an estimate for the pseudospectrum at the wavenumber  $k_{{\rm eff}}$ that corresponds to the average length of the k-vectors in that shell.

Thus, our power spectrum calculation consists of the following steps:

1.
Determination of the survey selection function i.e. mean underlying number density $\bar{n}(\vec{r})$ (including the survey geometry);
2.
Calculation of the overdensity field on a grid using a TSC mass assignment scheme;
3.
Fourier transformation of the gridded density field;
4.
Calculation of the raw 3D power spectrum $\vert F(\vec{k})\vert _{{\rm FFT}}^2$;
5.
Subtraction of the shot noise component from the raw spectrum;
6.
Recovery of the angle averaged pseudospectrum $\tilde{P}(k)$ using the iterative method of Jing (2005).
We have applied the above described power spectrum calculation method to a multitude of test problems, the results of which can be found in Hütsi (2005). In Appendix A we show one example where we successfully recover the underlying power spectrum of galaxy clusters from the VIRGO Hubble Volume simulations[*], after applying the selection criteria given in Figs. 2 and 3.

4 Power spectrum errors and covariance matrix

We determine power spectrum errors by three different methods:

1.
Prescription given by FKP that assumes the underlying density field to be Gaussian. This method also does not treat redshift space distortions. Under those simplifying assumptions the power spectrum variance can be expressed as:
 
                             $\displaystyle \sigma_{\tilde{P}}^2(k)$ = $\displaystyle \frac{2}{N_k^2}\sum \limits_{\vec{k'}}\sum \limits_{\vec{k''}}\vert\tilde{P}(k)Q(\vec{k'}-\vec{k''})+S(\vec{k'}-\vec{k''})\vert^2~,$ (13)
$\displaystyle Q(\vec{k})$ = $\displaystyle \int {\rm d}^3r~\bar{n}^2(\vec{r})w^2(\vec{r})\exp({\rm i}\vec{k}\cdot\vec{r})$  
  $\textstyle \simeq$ $\displaystyle \alpha \sum \limits_j \bar{n}(r^{\rm s}_j)w^2(r^{\rm s}_j)\exp({\rm i}\vec{k}\cdot\vec{r}^{\rm s}_j)~,$ (14)
$\displaystyle S(\vec{k})$ = $\displaystyle (1+\alpha)\int {\rm d}^3r~\bar{n}(\vec{r})w^2(\vec{r})\exp({\rm i}\vec{k}\cdot\vec{r})$  
  $\textstyle \simeq$ $\displaystyle \alpha(1+\alpha)\sum \limits_j w^2(r^{\rm s}_j)\exp({\rm i}\vec{k}\cdot\vec{r}^{\rm s}_j) .$ (15)

Here the sum is over all the wavevectors $\vec{k'}$ and $\vec{k''}$ populating the same k-space shell with radius k and thickness $\Delta k$, and Nk denotes the total number of modes in that shell. Since the direct summation over all the vector pairs $\vec{k'}$ and $\vec{k''}$ is very slow for the wide k-space shells and 5123 grid we use, a Monte Carlo sum is performed instead. Thus we calculate the average of the quantity $\vert\tilde{P}(k)Q(\vec{k'}-\vec{k''})+S(\vec{k'}-\vec{k''})\vert^2$ over the random pairs of vectors $\vec{k'}$ and $\vec{k''}$ from the same shell. For the result to converge properly we need on average $\sim $107 random pairs.
2.
The second method is a simple analytical approximation to the first one, also due to FKP (see also Tegmark et al. 1998). Here the variance is given as:

\begin{displaymath}\sigma_{\tilde{P}}^2(k) = \frac{2\tilde{P}^2(k)}{V_{\rm eff}V_{k}} ,
\end{displaymath} (16)

where $V_{k}=4\pi k^2\Delta k/(2\pi)^3$ is the volume of the k-space shell and $V_{\rm eff}$ is the effective volume given by:

\begin{displaymath}V_{\rm eff} = \frac{\left[\int {\rm d}^3r ~\bar{n}^2(\vec{r})...
...\left[1+ \frac{1}{\bar{n}(\vec{r})\tilde{P}(k)}\right]^2}\cdot
\end{displaymath} (17)

3.
The third method is a Monte Carlo approach that uses 1000 mock catalogs generated in the way described in Appendix B. Here, as we use the 2nd order Lagrangian perturbation theory, we get a good approximation for the mode-mode couplings that are induced during the quasi-nonlinear regime of the evolution of the density fluctuations. Also the large-scale redshift distortions are properly accounted for. In terms of the Halo Model (see Appendix C) the halo-halo clustering term is relatively well approximated. Contributions from the one-halo term can be added later, as these allow an analytic treatment.
The results of the power spectra for the SDSS DR4 LRG sample are shown in Figs. 4 and 5. In Fig. 4 the bin width $\Delta k \simeq 0.005 ~h~{\rm Mpc}^{-1}$, while in Fig. 5 $\Delta k \simeq 0.02 ~h~{\rm Mpc}^{-1}$. With different lines we have shown various model spectra, which will be discussed in Sect. 6.


  \begin{figure}
\par\includegraphics[width=7.5cm,clip]{3939_f04.eps}
\end{figure} Figure 4: Power spectrum of the SDSS LRG sample with the bin width $\Delta k \simeq 0.005 ~h~{\rm Mpc}^{-1}$. The upper solid line shows the best fitting model spectrum and the lower one corresponds to the linearly evolved matter power spectrum of the "concordance'' cosmological model multiplied by the square of the bias parameter b=1.95. Both of the spectra are convolved with a survey window. The dashed lines represent the corresponding unconvolved spectra.
Open with DEXTER


  \begin{figure}
\par\includegraphics[width=7.5cm,clip]{3939_f05.eps}
\end{figure} Figure 5: Power spectrum of the SDSS LRG sample with the bin width $\Delta k \simeq 0.02 ~h~{\rm Mpc}^{-1}$. The upper solid line shows the best fitting model spectrum and the lower one corresponds to the linearly evolved matter power spectrum of the "concordance'' cosmological model multiplied by the square of the bias parameter b=1.95. Both of the spectra are convolved with a survey window. The dashed lines represent the "smoothed'' versions of the above model spectra. The dotted line is the cubic spline fit to the data points.
Open with DEXTER


  \begin{figure}
\par\includegraphics[width=7.5cm,clip]{3939_f06.eps}\end{figure} Figure 6: The comparison of the different power spectrum error estimates. For clarity slight relative shifts of the data points have been applied. The errorbars resulting from the 1st method are the rightmost ones and the ones from the 3rd method are displayed in the middle. The lines show cubic spline fits to the data points. The solid line corresponds to the case when all the available galaxy data is used to find the angular mask of the survey, while the dashed line represents the case when LRGs only are used for this purpose.
Open with DEXTER

The comparison of the power spectrum errorbars calculated in the different ways is provided in Fig. 6. We see that the various error estimates are in very good agreement. In the following we will use only the errorbars given by the 3rd method.

So far we have only found the diagonal terms of the covariance matrix. In order to answer the question of how strongly different power spectrum bins are correlated, we must estimate the full covariance matrix.


  \begin{figure}
\par\includegraphics[width=16.5cm,clip]{3939_f07.eps}
\end{figure} Figure 7: Covariance ( left column) and correlation matrices ( right column). Top row represents the results from the FKP prescription (see Eq. (18)) and the middle row the ones from 1000 mock catalogs. The last row displays the nonlinear contribution due to the 1-halo term.
Open with DEXTER

The FKP result for the full covariance matrix, Cij, is a simple generalization of the Eq. (13):

 
$\displaystyle C_{ij}=\frac{2}{N_{k'}N_{k''}}
\sum \limits_{\vec{k'}}\sum \limit...
...ac{k_i+k_j}{2}\right)Q(\vec{k'}-\vec{k''})+S(\vec{k'}-\vec{k''})\right\vert^2~,$     (18)

where the k-vectors $\vec{k'}$ and $\vec{k''}$ lie in shells with width $\Delta k$ and radii ki and kj, respectively. The FKP approach, as mentioned above, does not treat mode couplings arising from nonlinear evolution and also from the redshift space distortions. Linear redshift space distortions can be, in principle, included in the FKP estimate for the covariance matrix. One can generalize the results presented in the Appendix of Zaroubi & Hoffman (1996), where the covariance matrix for the Fourier modes has been found. Since linear redshift distortions applied on a Gaussian field do not change the Gaussian nature, one can still use the result from the Appendix B of FKP that relates the power spectrum covariance matrix to the covariance matrix of the Fourier modes. Also one has to add the shot noise terms to the result of Zaroubi & Hoffman (1996). We have carried out this exercise, leading to the high dimensional integrals (up to 12 dim.) that turn out to be too time consuming to solve in practice. As from the mock catalogs we can hopefully obtain more realistic estimate for the covariance matrix[*] we have not followed this path any further.

The results for the covariance matrix calculation are given in Fig. 7. Here the left hand column shows the covariance and the right hand column the respective correlation matrices:

\begin{displaymath}r_{ij}=\frac{C_{ij}}{\sqrt{C_{ii}C_{jj}}}\cdot
\end{displaymath} (19)

The power spectrum binning is the same as shown in Fig. 5 i.e. $\Delta k \simeq 0.02 ~h~{\rm Mpc}^{-1}$. The top row represents the results from Eq. (18), while the middle row the ones from mock catalogs. Although the diagonal terms of the covariance matrices in the 1st and 2nd row are in very good agreement (see Fig. 6), the off-diagonal components differ strongly. This can be explained as the result of the extra mode-mode couplings that are not accounted for by the FKP approach. We see that even well separated power spectrum bins can be correlated at a $30{-}40\%$ level. The bottommost row in Fig. 7 represents the nonlinear contribution to the covariance matrix arising from the 1-halo term (see Appendix C). We see that this contribution is subdominant at the scales of interest to us[*].

In the following calculations we mostly use the covariance matrix given in the middle row of Fig. 7[*].

5 Relation to the true spectrum

Since masking in real space is equivalent to convolution in Fourier space, our measured power spectrum $\tilde{P}$ is actually a convolution of the real spectrum P with a survey window (see e.g. FKP):

 \begin{displaymath}
\tilde{P}(\vec{k})= \int \frac{{\rm d}^3k'}{(2\pi)^3}P(\vec{k})\vert W(\vec{k}-\vec{k'})\vert^2 ,
\end{displaymath} (20)

where

\begin{displaymath}W(\vec{k})=\int {\rm d}^3r ~ \bar{n}(\vec{r}) w(\vec{r}) \exp({\rm i}\vec{k}\cdot\vec{r}) ,
\end{displaymath} (21)

and the survey window $\vert W(\vec{k})\vert^2$ is normalized as follows:

\begin{displaymath}\int \frac{{\rm d}^3k}{(2\pi)^3}\vert W(\vec{k})\vert^2=1 .
\end{displaymath} (22)

The angle averaged survey window |W(k)|2 is plotted in Fig. 8. Here the core part of the window is well approximated by the functional form:

\begin{displaymath}\vert W(k)\vert^2 = \frac{1}{1+\left(\frac{k}{a}\right)^2+\le...
...frac{k}{b}\right)^4}, \quad (a\simeq 0.0030, b\simeq 0.0028) ,
\end{displaymath} (23)

and asymptotic wings are close to the power law with spectral index -4. These approximations are shown with dashed lines in Fig. 8. With the gray shaded stripe we have marked the scales where |W(k)|2 is above $1\%$ of its maximum value. This stripe serves as a rough guide to the effective width of the survey window and it is also shown in many of the following figures.

Since the survey geometry of the analyzed SDSS LRG sample is far from spherically symmetric, an isotropized window in Fig. 8 gives only a poor representation of the true 3D window, which is displayed as an isosurface corresponding to the isovalue of 0.01 in Fig. 9.


  \begin{figure}
\par\includegraphics[width=7.5cm,clip]{3939_f08.eps}
\end{figure} Figure 8: Isotropized survey window. Here the normalization is taken such that |W(0)|=1. The light gray stripe marks the region where the window is above $1\%$ of its maximum value of 1. Dashed lines show approximations discussed in the text.
Open with DEXTER


  \begin{figure}
\par\includegraphics[angle=90,width=7.5cm,clip]{3939_f09.ps}\end{figure} Figure 9: 3D survey window embedded in a box with a side length of $0.04~h~{\rm Mpc}^{-1}$. Here the isosurface corresponding to $1\%$ of the maximum value of the window is shown. Note the symmetry of the window, $\vert W(\vec{k})\vert^2=\vert W(-\vec{k})\vert^2$, as expected when taking a modulus of the Fourier transform of a real 3D scalar function.
Open with DEXTER

In order to compare theoretical models to the measured power spectrum we have to take into account the smearing effects caused by the survey window. Using Eq. (20) we can express an isotropized power spectrum as:

 \begin{displaymath}
\tilde{P}(k) = \int \frac{{\rm d}\Omega_{\vec{k}}}{4\pi}\tilde{P}(\vec{k}) = \int {\rm d}k'~k'^2P(k')K(k',k) ~,
\end{displaymath} (24)

where the coupling kernels[*]:

\begin{displaymath}K(k',k) = K(k,k') = \frac{1}{2\pi^2}\int \frac{{\rm d}\Omega_...
...rm d}\Omega_{\vec{k'}}}{4\pi}\vert W(\vec{k}-\vec{k'})\vert^2.
\end{displaymath} (25)

Numerically evaluated coupling kernels along with the analytical approximations (see Appendix D) for the analyzed galaxy sample are presented in Fig. 10. Here the solid lines correspond to the numerical results and the dashed ones represent an analytical approximation. We have used the notation $K_i(k)\equiv K(k_i,k)$ where ki denote the central values of the power spectrum bins shown in Fig. 5.


  \begin{figure}
\par\includegraphics[width=8cm,clip]{3939_f10.eps}
\end{figure} Figure 10: Coupling kernels $K_i(k)\equiv K(k_i,k)$ for the power spectrum bins ki shown in Fig. 5. Numerically evaluated kernels are shown with solid lines. The dashed lines correspond to the fitting functions given in Appendix D.
Open with DEXTER

  
6 Model spectra

It is well known that redshift space distortions and nonlinear effects modify simple linear spectra. In order to treat these effects we make use of the very successful analytical model - the Halo Model. For a review see Cooray & Sheth (2002) (see also Seljak 2000). The details of the model we use are presented in Appendix C. The model introduces four free parameters: $M_{{\rm low}}$, $\alpha$, M0 and $\gamma$. Here $M_{{\rm low}}$ is the lower cutoff of the halo mass i.e. below that mass halos are assumed to be "dark''. $\alpha$ and M0 are the parameters of the mean of the halo occupation distribution  $\langle N\vert M \rangle$, which gives the average number of galaxies per halo with mass M. We take $\langle N\vert M \rangle$ to be a simple power law:

\begin{displaymath}\langle N\vert M \rangle = \left ( \frac{M}{M_0} \right )^\alpha\cdot
\end{displaymath} (26)

The last parameter, $\gamma$, is the amplitude factor for the virial velocities of galaxies inside dark matter halos. A one dimensional velocity dispersion of the galaxies inside a halo with mass M is taken to follow the scaling of the isothermal sphere model:

\begin{displaymath}\sigma = \gamma \sqrt {\frac{GM}{2R_{\rm vir}}} ,
\end{displaymath} (27)

where $R_{\rm vir}$ is the virial radius of the halo.

For the model fitting we have used the Levenberg-Marquardt method as described in Press et al. (1992) with modifications (described in Appendix E) that allow us to incorporate correlations between the data points. As the input data we take the power spectrum estimates given in Fig. 5. The covariance matrix used is the one shown in the middle row of Fig. 7. We also perform fits where we use one additional power spectrum bin on a larger scale (not shown in Fig. 5). All of this data are given in a tabular form in Appendix G. The transfer functions needed for the linear spectra are taken from Eisenstein & Hu (1998). There the authors also provide transfer function fits where the baryonic acoustic oscillations have been removed. We use these "smoothed out'' transfer functions to assess the significance of the oscillatory features we see in the data. Throughout this paper we have kept the cosmology fixed to the best-fit WMAP "concordance'' model (Spergel et al. 2003). The implications for the cosmology, and especially for the dark energy equation of state parameter, are planned to be worked out in a future paper.

  \begin{figure}
\par\includegraphics[width=7cm,clip]{3939_f11.eps}
\end{figure} Figure 11: $1\sigma $ error contours for the free model parameters. Best fit parameter values are marked with crosses.
Open with DEXTER

As the cosmology is kept fixed, we have only four free parameters. In order to eliminate some of the degeneracies between the parameters we have imposed one additional constraint. We have demanded that the resulting number of galaxies should agree with the one that is observed with a relative error of $1\%$ i.e. $(51,763 \pm 518)$[*]. The resulting $1\sigma $ error "ellipses'' for the free parameters are shown in Fig. 11. The "ellipses'' appear deformed since instead of $M_{{\rm low}}$ and M0 we have fitted $\log(M_{{\rm low}})$ and $\log(M_0)$. We have marked with crosses the best fit values: $M_{{\rm low}} \simeq 3\times 10^{12} h^{-1}~M_{\odot}$, $\alpha \simeq 0.9$, $M_0 \simeq 1.4 \times 10^{14} h^{-1}~M_{\odot}$ and $\gamma \simeq 0.7$. The model spectra corresponding to these best fit parameters are shown in Figs. 4 and 5. In both figures we have also given the simple linear spectra multiplied by the square of the bias parameter b=1.95. In Fig. 4 we have additionally demonstrated the effect of the window convolution. There the dashed lines correspond to the unconvolved case. In Fig. 5, along with the "wiggly'' spectra we have shown their "smoothed'' counterparts. Using all the 16 power spectrum bins (the 1st not shown in Fig. 5) plus an additional constraint on the total number of galaxies, resulting in 17-4 = 13 independent degrees of freedom, we obtain $\chi^2$ values of 8.8 and 19.9 for the "wiggly'' and "smoothed''[*] models, respectively. So the models with oscillations are favored by  $3.3 \sigma$ over their "smoothed'' counterparts[*]. Since both models have the same number of free parameters, and if additionally the assumption of Gaussianity is valid, the Bayesian approach should give similar results. Actually, Bayesian results should favor "wiggly'' models even more, since prior weight for these should probably be taken higher (assuming knowledge of the other experimental results).

7 Determination of the acoustic scale

To measure the scale of the acoustic oscillations we divide the spectrum shown in Fig. 5 by the best fitting "smoothed'' spectrum. The result of this procedure is given in the upper panel of Fig. 12. There the solid line shows a cubic spline fit to the data points and the long-dashed line corresponds to the best fitting model spectrum also shown in Fig. 5. The above data is fitted with a parametric form:

 \begin{displaymath}
f(x)=1+c_1\cdot\sin(c_2\cdot x)\exp\left[\left(-\frac{x}{c_3}\right)^{c_4}\right]\cdot
\end{displaymath} (28)

Again we use the Levenberg-Marquardt method with the data covariance matrix obtained from mock catalogs. After marginalizing over the other parameters we find the best fitting value of $(105.4 \pm 2.3)~h^{-1}~{\rm Mpc}$ for the parameter c2[*]. The best fitting member of the parametric family in Eq. (28) is shown with short-dashed lines in the upper panel of Fig. 12. Using the FKP covariance matrix instead gives an acoustic scale of  $(105.4 \pm 2.8)~h^{-1}~{\rm Mpc}$.


  \begin{figure}
\par\includegraphics[width=7.5cm,clip]{3939_f12.eps}
\end{figure} Figure 12: Upper panel: power spectrum from Fig. 5 divided by the best fitting "smoothed'' spectrum. Solid line shows a cubic spline fit to the data points and long-dashed line corresponds to the best "wiggly'' model. The short-dashed line represents the most favorable fit from the parametric family of Eq. (28). Lower panel: various input power spectra used to calculate the two-point correlation function. The dashed line is the cubic spline fit from the upper panel. The solid lines represent a transition sequence from the best fitting "wiggly'' model to the best "smoothed'' model. In each step we have erased more oscillatory features. For clarity slight vertical shifts have been introduced.
Open with DEXTER

The sinusoidal modulation in the power spectrum is a consequence of the adiabatic initial conditions. By relaxing this assumption and fitting with a more general functional form:

\begin{displaymath}f(x)=1+c_1\cdot\sin(c_2\cdot x + c_3)\exp\left[\left(-\frac{x}{c_4}\right)^{c_5}\right]
\end{displaymath} (29)

instead, we get the following value for the acoustic scale: $(103.0 \pm 7.6)~h^{-1}~{\rm Mpc}$. In the case of the FKP covariance matrix the corresponding value is $(103.1 \pm 9.1)~h^{-1}~{\rm Mpc}$.

Eisenstein et al. (2005) determined various distance scales (like Dv, which is a specific combination of the comoving distances along and perpendicular to the line of sight (see their Eq. (2))) and their ratios, using SDSS LRGs and the constraints from other cosmological sources. The typical relative accuracy of these measurements is $\sim $$4 \%$, which seems worse than the accuracy of the acoustic scale measurement, $(105.4 \pm 2.3)~h^{-1}~{\rm Mpc}$ i.e. $\sim $$2\%$, presented in this paper. This apparent inconsistency can be attributed to the fact that in our analysis, as stated above, we have kept the cosmology fixed to the WMAP "concordance'' model, whereas the Eisenstein et al. (2005) estimates include the extra uncertainties due to the imperfect knowledge of the various cosmological parameters. Of course, the given length of the acoustic scale, $(105.4 \pm 2.3)~h^{-1}~{\rm Mpc}$, can be easily transformed in order to accommodate other preferences for the background cosmology. We also note that the use of the parametric form in Eq. (28) might be too restrictive, since the acoustic modulation in the case of adiabatic models can be only approximately described as a damped sinusoidal wave (Eisenstein & Hu 1998). For this reason the given sound horizon constraint should not be used in cosmological parameter studies. Instead one should directly use the measured power spectrum in combination with the parametrized models that are physically well motivated.

8 Correlation function analysis

We determine the two-point correlation function of the SDSS LRGs using the edge-corrected estimator given by Landy & Szalay (1993):

 \begin{displaymath}
\xi (r) = \frac{DD - 2DR + RR}{RR} ,
\end{displaymath} (30)

which has minimal variance for a Poisson process. Here DD, DR and RR represent the respective normalized data-data, data-random and random-random pair counts in a given distance range. Random catalogs were generated with 25 times the number of objects in the main catalog. We calculated correlation function for 10 bins ( $r_i,~i=1\ldots10$) in the pair distance range of  $60 \ldots 160~h^{-1}~{\rm Mpc}$. The errors were estimated by a "jackknife'' technique. For this purpose we divided the full sample into 22 separate regions each containing $\sim $2350 galaxies (see Fig. 3). The two-point function was calculated 22 times, each time omitting one of the regions. Denoting the resulting estimates as $\xi_j(r_i),~(j=1\ldots22)$, the "jackknife'' estimate for the variance is (see e.g. Lupton 1993):
$\displaystyle \sigma^2_{\xi}(r_i) = \frac{N-1}{N}\sum \limits_{j=1}^N \left[\xi_j(r_i) -\bar{\xi}(r_i)\right]^2 ,$     (31)
$\displaystyle \bar{\xi}(r_i) = \frac{1}{N}\sum \limits_{j=1}^N \xi_j(r_i) ,$     (32)

where in our case N=22. The results of this calculation are presented in the left panel of Fig. 13. With the crosses and dashed-line errorbars we show the two-point function as determined by Eisenstein et al. (2005). In general our results agree reasonably well with their calculations.
  \begin{figure}
\par\includegraphics[width=7.5cm,clip]{3939_f13.eps}
\end{figure} Figure 13: Left panel: two-point correlation functions as determined in this paper (circles with solid lines) and by Eisenstein et al. (2005). Right panel: correlation functions corresponding to the models shown in the lower panel of Fig. 12 in comparison to the one obtained directly from the data. Here all the data points have been lowered by 0.0035.
Open with DEXTER

It would be interesting to study how the oscillations in the observed power spectrum transform into the peak in the two-point correlation function seen at the scale of $\sim $ $110 ~h^{-1}~{\rm Mpc}$. For this purpose we use the cubic spline fit shown in Fig. 12 and extend it outside of the observed range by smoothly joining it to the power spectrum of the best fitting "smoothed'' model. The correlation function is now simply calculated as the Fourier transform of the power spectrum[*]. The resulting correlation function is plotted with a dashed line in the right panel of Fig. 13. To study the significance of the oscillatory features in the power spectrum in relation to the observed peak in the correlation function, we have calculated correlation functions for several models that have oscillations "switched off'' at various scales. The spectra of these models are shown with solid lines in the lower panel of Fig. 12, where for clarity we have introduced slight vertical shifts between the curves, so that the scales where the transition to the featureless spectrum takes place are easily visible. The corresponding correlation functions are given with solid lines in the right hand panel of Fig. 13. As expected, we see the peak in the correlation function becomes broader and also decreases in amplitude as we successively erase features in the power spectrum. This clearly demonstrates the importance of many of the fluctuations in the power spectrum to produce a relatively sharp feature in the two-point correlation function.

In order to achieve good agreement we have lowered all the data points by 0.0035 in the right hand panel of Fig. 13. Similar shifts were also suggested in Eisenstein et al. (2005) in order to get a better match to the theoretical models. A 0.0035 shift in $\xi$ translates to the $0.175 \%$ shift in the mean density. Thus, if one wishes to determine the amplitude of the correlation function correctly at those large scales, one has to determine the survey selection function with a very high precision, which in practice is very difficult to achieve. By using model spectra that have more large scale power than the "concordance'' cosmology predicts (as might be suggested by Fig. 4), we are able to match the amplitude of the correlation function without any additional vertical shifts. Here we try to avoid making any definite conclusions. The behavior of the power spectrum on the largest scales is a whole interesting topic on its own and there exist much better methods than the direct Fourier approach to investigate these issues (see e.g. Tegmark et al. 1998).

9 Comparison with other surveys

In this section we compare our power spectrum measurements with the ones obtained by Percival et al. (2001) and Cole et al. (2005) for the 2dF redshift survey and by Tegmark et al. (2004) for the SDSS main galaxy sample. The results of this comparison are provided in Figs. 14 and 15. For clarity we have given a variant of Fig. 14 where we have omitted the errorbars. The amplitudes of the SDSS main and 2dF spectra have been freely adjusted to match the clustering strength of the SDSS LRGs. The corresponding bias parameters with respect to the SDSS LRGs are 0.53, 0.61 and 0.50 for the 2dF sample analyzed by Percival et al. (2001), for the one analyzed by Cole et al. (2005), and for the SDSS main sample, respectively. Percival et al. (2001) also provide power spectrum measurements for $k \ga 0.15~h~{\rm Mpc}^{-1}$ but without errorbars. These small-scale measurements along with our SDSS LRG results are shown with solid lines in Fig. 14.

In general the shapes of the spectra agree remarkably well. With the only exception of the Tegmark et al. (2004) results, the power spectrum bins are significantly correlated. Also the Tegmark et al. (2004) measurements are corrected for the redshift space distortions.

10 Discussion and conclusions

In this paper we have calculated the redshift-space power spectrum of the SDSS DR4 LRG sample, finding evidence for a series of acoustic features down to the scales of $\sim $ $0.2~h~{\rm Mpc}^{-1}$. Models with baryonic oscillations are favored by $3.3 \sigma$ over their "smoothed'' counterparts without any oscillatory behavior. Using the obtained power spectrum we predict the shape of the spatial two-point correlation function, which agrees very well with the one obtained directly from the data. Also, the directly calculated correlation function is consistent with the results obtained by Eisenstein et al. (2005). We have made no attempts to put constraints on the cosmological parameters, rather we have assumed the "concordance'' cosmological model in our analysis. The derived acoustic scale $(105.4 \pm 2.3)~h^{-1}~{\rm Mpc}$ agrees well with the best-fit WMAP "concordance'' model prediction of $\simeq$ $106.5 ~h^{-1}~{\rm Mpc}$.


  \begin{figure}
\par\includegraphics[width=7.5cm,clip]{3939_f14.eps}
\end{figure} Figure 14: The comparison of spectra from different surveys.
Open with DEXTER


  \begin{figure}
\par\includegraphics[width=7.5cm,clip]{3939_f15.eps}
\end{figure} Figure 15: The same as Fig. 14 with the errorbars omitted.
Open with DEXTER

The existence of baryonic features in the galaxy power spectrum is very important, allowing one (in principle) to obtain the Hubble parameter H and angular diameter distance dA as a function of redshift, this way opening up a possibility to constrain properties of the dark energy (Hu & Haiman 2003). The currently existing largest redshift surveys, which are still quite shallow, do not yet provide enough information to do this properly. On the other hand, it is extremely encouraging that even with the current generation of redshift surveys we are already able to see traces of acoustic oscillations in the galaxy power spectrum, showing the promise for dedicated future surveys like KAOS. We have seen that acoustic features seem to survive at mildly nonlinear scales ( $k \ga 0.1~h~{\rm Mpc}^{-1}$), which is in agreement with the results of the recent N-body simulations (Springel et al. 2005; Seo & Eisenstein 2005). In order to fully exploit the available information one needs a complete understanding of how nonlinear effects influence these features. Nonlinear bias and redshift space distortions also add complications. In general, redshift-space distortions, biasing and nonlinear evolution do not create any oscillatory modulation in the power spectrum and so acoustic features should be readily observable. So far there have been only a few works studying these important issues (e.g. Springel et al. 2005; White 2005; Seo & Eisenstein 2005) and currently we do not have a full theoretical description of them. In our paper we have modeled the above-mentioned effects using the results from the 2nd order Lagrangian perturbation theory in combination with the Halo Model. Although these models are very successful in capturing many important aspects of structure formation, they are still approximations.

The existence of baryonic oscillations in the galaxy power spectrum gives us important information about the underlying cosmological model and the mechanism of structure formation. First, it confirms the generic picture of the gravitational instability theory where the structure in the Universe is believed to have been formed by the gravitational amplification of the small perturbations layed down in the early Universe. Under linear gravitational evolution all the density fluctuation modes evolve independently i.e. all the features in the power spectrum will be preserved. We are now indeed able to identify features in the low redshift galaxy power spectrum that correspond to the fluctuations seen in the CMB angular power spectrum (which probes redshifts $z \sim 1100$), providing strong support for this standard picture of structure formation. We can even probe scales that are inaccessible to CMB studies due to the strong damping effects and steeply rising influence of the secondary anisotropies, effectively reaching wavenumbers that correspond to the 6th-7th peak in the CMB angular power spectrum. Second, the ability to observe baryonic features in the low redshift galaxy power spectrum demands rather high baryonic to total matter density ratio. In Blanchard et al. (2003) it has been shown that it is possible to fit a large body of observational data with an Einstein-de Sitter type model if one adopts a low value for the Hubble parameter and relaxes the usual assumptions about the single power law initial spectrum. In the light of our results these models are disfavored due to the fact that the high dark matter density completely damps the baryonic features. Finally, purely baryonic models are also ruled out since for them the expected acoustic scale would be roughly two times larger than observed here[*]. Thus the data seems to demand a weakly interacting nonrelativistic matter component and all the models that try to replace this dark matter component with "something else'' e.g. modifying the laws of gravity, might have severe difficulties in fitting these new observational constraints.

Acknowledgements
I thank Rashid Sunyaev and the referee for helpful comments and suggestions. I acknowledge the International Max-Planck Research School on Astrophysics for a graduate fellowship and the support provided through Estonian Ministry of Education and Recearch project TO 0062465S03 and and ESF grant 5347.

Funding for the creation and distribution of the SDSS Archive has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Aeronautics and Space Administration, the National Science Foundation, the US Department of Energy, the Japanese Monbukagakusho, and the Max Planck Society. The SDSS Web site is http://www.sdss.org/

The SDSS is managed by the Astrophysical Research Consortium (ARC) for the Participating Institutions. The Participating Institutions are The University of Chicago, Fermilab, the Institute for Advanced Study, the Japan Participation Group, The Johns Hopkins University, the Korean Scientist Group, Los Alamos National Laboratory, the Max-Planck-Institute for Astronomy (MPIA), the Max-Planck-Institute for Astrophysics (MPA), New Mexico State University, University of Pittsburgh, University of Portsmouth, Princeton University, the United States Naval Observatory, and the University of Washington.

References

 

  
Online Material

  
Appendix A: Test problem

Here we present a test of our power spectrum calculation software[*]. As the input we use the z=0 cluster catalog of the VIRGO Hubble Volume simulations[*], which covers the comoving volume of $3000~h^{-3}~{\rm Mpc}^3$ and contains 1 560 995 clusters above the mass limit of $6.75\times 10^{13} h^{-1}~M_{\odot}$. The average bias parameter of this catalog is b = 1.9, which is comparable to the SDSS LRG value of b = 1.95. The power spectrum of the full sample is shown in Fig. A.1 with a solid line. Here for clarity we have not shown the errorbars, which are rather small for a sample of that size. Out of the full sample we generate 50 mock catalogs that have the same radial and angular selection functions as the SDSS LRG sample analyzed in this paper (see Figs. 2 and 3). The mean number of objects in the resulting catalogs is $\sim $18 500 i.e. the number density is roughly one third of the spatial density of the SDSS LRGs. The observer's location and pointing angles are taken randomly for each of the catalogs. The mean recovered power spectrum with $1\sigma $ errorbars is shown in Fig. A.1. We see that the power spectrum of the underlying sample is recovered very well. On the largest scales there are some deviations, which can be explained as being caused by the smearing effect of the survey window. This is demonstrated by the dotted lines, where the lower/upper curve corresponds to the model spectrum with/without survey window convolution applied.


  \begin{figure}
\par\includegraphics[width=7cm,clip]{3939_f16.eps}
\end{figure} Figure A.1: Power spectra of galaxy clusters from the z=0 Hubble Volume simulation box. The solid line represents the spectrum for the full sample of 1 560 995 clusters. The circles with errorbars denote the recovered spectrum from the 50 mock catalogs having similar selection effects as the analyzed SDSS LRG sample. The dotted lines demonstrate the convolution effect of the survey window on the best fitting model spectrum.

  
Appendix B: Mock catalogs

We build mock catalogs for the SDSS LRG by a 3 step procedure:
1.
Generation of the density field using an optimized 2nd order Lagrangian perturbation calculation (2LPT).
2.
Poisson sampling of the generated density field with the intensity of the process adjusted to obtain a galaxy sample that has a clustering strength enhanced by a factor b2 with respect to the underlying field, and a number density equal to the observed LRG sample density at the minimal used redshift of 0.16 (see Fig. 2).
3.
Extraction of the final catalog by applying the radial and angular selection function as given in Figs. 2 and 3, respectively.
In contrast to the Eulerian perturbation theory, where one does a perturbative expansion of the density contrast field, the Lagrangian approach considers an expansion of the particle trajectories (see e.g. Sahni & Coles 1995; Bernardeau et al. 2002; Bouchet et al. 1995; Buchert & Ehlers 1993). Here the central quantity is the displacement field $\vec{\Psi} (\vec{q})$, which relates particle's initial comoving position (Lagrangian position) $\vec{q}$ to its final Eulerian position $\vec{x}$:

 \begin{displaymath}
\vec{x} = \vec{q} + \vec{\Psi} (\vec{q}) .
\end{displaymath} (B.1)

Due to the decay of the rotational perturbation modes in the expanding Universe each order of the perturbation theory displacement field separates into a time-dependent and a Lagrangian position dependent factor (Ehlers & Buchert 1997). The position dependent part, due to its irrotational nature can be given as a gradient of a scalar potential. As a result, one can expand the displacement field as follows:

 \begin{displaymath}
\vec{\Psi} (\vec{q}) = D_1 \nabla_q \phi^{(1)} + D_2 \nabla_q \phi^{(2)} .
\end{displaymath} (B.2)

Here the 1st term describes the classical Zeldovich approximation (Zel'Dovich 1970). The time independent potentials $\phi^{(1)}$ and $\phi^{(2)}$ are found from the Poisson equations:

\begin{displaymath}\Delta \phi^{(1)}(\vec{q}) = -\delta(\vec{q})
\end{displaymath} (B.3)

and

\begin{displaymath}\Delta \phi^{(1)}(\vec{q}) = \frac{1}{2}\sum \limits_{i} \sum...
...- \phi^{(1)}_{,ij}(\vec{q})\phi^{(1)}_{,ji}(\vec{q})\right ) ,
\end{displaymath} (B.4)

where ,i denotes the partial derivative with respect to the Lagrangian coordinate qi. $\delta(\vec{q})$ is the initial density contrast. We generate $\delta(\vec{q})$ using the standard Zeldovich approximation on a regular cubical grid.


  \begin{figure}
\par\includegraphics[width=11cm,clip]{3939_f17.ps}
\end{figure} Figure B.1: A $25~h^{-1}~{\rm Mpc}$ thick slice through a $1280~h^{-1}~{\rm Mpc}$ computational box. The gray scale image represents the underlying density field obtained by the optimized 2LPT approach. White dots mark the positions of the "galaxies'' generated by the Poisson sampler.

D1 in Eq. (B.2) is the linear growth factor. The second-order growth factor D2 for flat models with a cosmological constant is to a good precision approximated as (Bouchet et al. 1995):

\begin{displaymath}D_2 \simeq -\frac{3}{7}\Omega_{\rm m}^{-1/143}D_1^2 .
\end{displaymath} (B.5)

According to Eqs. (B.1) and (B.2) the peculiar velocity field is given as:

 \begin{displaymath}
\vec{v} = D_1 f_1 H \nabla_q \phi^{(1)} + D_2 f_2 H \nabla_q \phi^{(2)} .
\end{displaymath} (B.6)

Here $H \equiv \frac{\dot{a}}{a}$ and $f_i \equiv \frac{{\rm d} \ln D_i}{{\rm d} \ln a}$. For flat models with a cosmological constant logarithmic derivatives of the growth factors can be approximated as (Bouchet et al. 1995):

\begin{displaymath}f_1 \simeq \Omega_{\rm m}^{5/9},\quad f_2 \simeq 2\Omega_{\rm m}^{6/11} .
\end{displaymath} (B.7)

The Lagrangian perturbative approach works well up to the 1st shell-crossing. After that the formed sharp caustic structures will start to disappear, since the particles keep moving without experiencing the gravitational pull of the dense sheets/filaments. It is possible to resolve this problem significantly by filtering out the small-scale Fourier modes. This is what is meant by the "optimization''. The method applied to the 1st order Lagrangian perturbation calculation is known as the truncated Zeldovich approximation (e.g. Coles et al. 1993; Melott et al. 1994; Weiss et al. 1996). Weiss et al. (1996) suggest to remove the small-scale power by applying a Gaussian k-space filter with a characteristic smoothing scale $k_{\rm gs}$ to the initial density field. Thus the power spectrum of the filtered field is given by:

\begin{displaymath}P_{{\rm optimized}}(k) = P(k)\exp \left (-\frac{k^2}{k^2_{\rm gs}} \right )\cdot
\end{displaymath} (B.8)

They recommend the value $k_{\rm gs}\simeq 1.2 k_{\rm nl}$, where the nonlinearity scale $k_{\rm nl}$ is defined as:

\begin{displaymath}\frac{D_1^2}{(2\pi)^3} \int \limits_0^{k_{\rm nl}}{\rm d}^3k~P(k) = 1 .
\end{displaymath} (B.9)

Although they studied only models with $\Omega_{\rm m}=1$, it was later shown by Hamana (1998) that this "recipe'' performs well for arbitrary Friedmann-Lemaître-Robertson-Walker models.

In our calculations we assume the WMAP "concordance'' cosmology (Spergel et al. 2003). The linear power spectrum is taken from Eisenstein & Hu (1998). We build a 2LPT density field on a 2563-grid with $5~h^{-1}~{\rm Mpc}$ cell size using the same number of particles as the number of grid cells[*]. Four copies of this box are combined to form a large $2560 \times 2560 \times 1280 ~ h^{-3}~{\rm Mpc}^3$ volume. Out of that big box a sample of "galaxies'' is selected with a radial number density as given in Fig. 2 and with an angular mask presented in Fig. 3. The parameters of the Poisson sampler[*] are tuned to give a sample with a bias parameter $b \simeq 2$ in agreement with the observed value for the SDSS LRG sample. The redshift-space catalog is built by altering the radial distances of the "galaxies'' by $\vec{v}_r/H_0$, where $\vec{v}_r$ is the radial component of the peculiar velocity field (see Eq. (B.6)) and $H_0 = 100~h~{\rm km/s/Mpc}$.

In Fig. B.1 we show a $25~h^{-1}~{\rm Mpc}$ thick slice through a box with $1280~h^{-1}~{\rm Mpc}$ side length. The underlying density field is presented as a gray scale image with white dots marking the positions of the "galaxies''. The power spectrum of the sample of $\sim $350 000 "galaxies'' is shown in Fig. B.2[*]. We see that the shape of the spectrum is in good agreement with the linearly evolved power spectrum up to the scales of $k \sim 0.5 ~h~{\rm Mpc}^{-1}$.

This approach gives us a "galaxy'' sample that has realistic large-scale clustering properties. In terms of the Halo Model (see Appendix C) the halo-halo clustering term is properly accounted for. 2LPT also gives reasonably accurate higher order correlations on quasi-nonlinear scales (e.g. Scoccimarro & Sheth 2002; Bouchet et al. 1995).


  \begin{figure}
\par\includegraphics[width=7.5cm,clip]{3939_f18.eps}\end{figure} Figure B.2: The power spectrum of $\sim $350 000 "galaxies'' from the simulation box shown in Fig. B.1. The solid line shows the linearly evolved input spectrum multiplied by the square of the bias parameter b=2.0.

  
Appendix C: Power spectrum from the halo model

The halo model description of the spatial clustering of galaxies is a development of the original idea by Neyman & Scott (1952), where one describes the correlations of the total point set as arising from two separate terms: (i) the 1-halo term, that describes the correlations of galaxies populating the same halo, (ii) the 2-halo term, which takes into account correlations of the galaxies occupying different halos. For a thorough review see Cooray & Sheth (2002). Here we briefly give the results relevant to the current paper (see Cooray 2004; Seljak 2001).

The power spectrum of galaxies in redshift space can be given as:

P(k) = P1h(k) + P2h(k) , (C.1)

where the 1-halo term:
$\displaystyle P^{1h}(k) = \int {\rm d}M ~n(M)\frac{\langle N(N-1)\vert M \rangle}{\bar{n}^2}\mathcal{R}_p(k\sigma)\vert u_{\rm g}(k\vert M)\vert^p ,$     (C.2)
$\displaystyle p = \left\{
\begin{array}{ll}
1 & {\rm\quad if \quad} \langle N(N...
...\
2 & {\rm\quad if \quad} \langle N(N-1)\vert M \rangle > 1
\end{array}\right.$     (C.3)

and the 2-halo term:

\begin{displaymath}P^{2h}(k) = \left ( \mathcal{F}_{\rm g}^2 + \frac{2}{3}\mathc...
...{\rm g} + \frac{1}{5}\mathcal{F}_v^2 \right )P_{\rm lin}(k) .
\end{displaymath} (C.4)

Here:
  
    $\displaystyle \mathcal{R}_p \left ( \alpha = k \sigma \sqrt {\frac{p}{2}} \right ) = \frac{\sqrt{\pi}}{2}\frac{{\rm erf}(\alpha)}{\alpha}~,$ (C.5)
    $\displaystyle \mathcal{F}_{\rm g} = \int {\rm d}M ~n(M) b(M) \frac{\langle N\vert M \rangle}{\bar{n}}
\mathcal{R}_1(k\sigma) u_{\rm g}(k\vert M)~,$ (C.6)
    $\displaystyle \mathcal{F}_v = f \cdot \int {\rm d}M~n(M) b(M) \mathcal{R}_1(k\sigma) u(k\vert M) .$ (C.7)

In the above expressions n(M) is the mass function and b(M) halo bias parameter. We calculate them using the prescription by Sheth & Tormen (1999) and Sheth et al. (2001). $\bar{n}$ represents the mean number density of galaxies:

\begin{displaymath}\bar{n} = \int {\rm d}M~ n(M) \langle N\vert M \rangle .
\end{displaymath} (C.8)

We take the mean of the halo occupation distribution in the following form:

 \begin{displaymath}
\langle N\vert M \rangle = \left ( \frac{M}{M_0} \right )^\alpha ,
\end{displaymath} (C.9)

where M0 and $\alpha$ are free parameters. The second moment is chosen as (see Cooray 2004):
 
    $\displaystyle \langle N(N-1)\vert M \rangle = \beta^2(M) \langle N\vert M \rangle ^2 ,$ (C.10)
    $\displaystyle \beta(M) = \left \{
\begin{array}{ll}
\frac{1}{2} \log \left ( \f...
...0^{13} \ h^{-1}~M_{\odot} \\  [2mm]
1 & {\rm\quad otherwise.}\end{array}\right.$ (C.11)

f in Eq. (C.7) denotes the logarithmic derivative of the linear growth factor: $f \equiv \frac{{\rm d} \ln D_1}{{\rm d} \ln a}$. u(k|M) and $u_{\rm g}(k\vert M)$ are the normalized Fourier transforms of the dark matter and galaxy density distributions within a halo of mass M. In our calculations we take both of these distributions given by the NFW profile (Navarro et al. 1997) and the concentration parameter c(M) is taken from Bullock et al. (2001). The one dimensional velocity dispersion of the galaxies inside a halo with mass M is taken to follow the scaling of the isothermal sphere model:

 \begin{displaymath}\sigma = \gamma \sqrt {\frac{GM}{2R_{\rm vir}}} ,
\end{displaymath} (C.12)

where $R_{\rm vir}$ is the virial radius of the halo and $\gamma$ is a free parameter.

After specifying the background cosmology the above described model has four free parameters: M0, $\alpha$ (Eq. (C.9)), $\sigma$ (Eq. (C.12)) and $M_{{\rm low}}$. The last parameter  $M_{{\rm low}}$ represents the lower boundary of the mass integration i.e. halos with masses below  $M_{{\rm low}}$ are assumed to be "dark''.

One can also use the halo model to estimate nonlinear contributions to the power spectrum covariance matrix. The additional term to the covariance matrix CijNL (i,j denote power spectrum bins) arising from the parallelogram configurations of the trispectrum[*] is given by (Cooray 2004):

 
    $\displaystyle C_{ij}^{NL} = \frac{T_{ij}}{V}=\frac{1}{V}\int \limits_i \frac{{\rm d}^3 k~}{V_i} \int \limits_j \frac{{\rm d}^3 k~}{V_j} \int {\rm d}M~n(M) \cdot$  
    $\displaystyle \frac{\langle N(N-1)(N-2)(N-3)\vert M \rangle}{\bar{n}^4} \vert u_{\rm g}(k_i\vert M)\vert^2 \vert u_{\rm g}(k_j\vert M)\vert^2 .$ (C.13)

$\int \limits_i$ denotes an integral over a k-space shell centered at wavenumber ki with a volume $V_i=4\pi k_i^2 \Delta k$. The 4th moment of the halo occupation distribution is taken as:
$\displaystyle \langle N(N-1)(N-2)(N-3)\vert M \rangle = \beta^2(M)\left[2\beta^2(M)-1\right]\cdot
\left[3\beta^2(M)-2\right]\langle N\vert M \rangle^4 ,$     (C.14)

where $\beta(M)$ and $\langle N\vert M \rangle$ are given in Eqs. (C.11) and (C.9) above. Performing calculations in redshift space a factor of  $\mathcal{R}_p^2(k\sigma)$ (see Eq. (C.5)) must also be included in Eq. (C.13).

  
Appendix D: Fitting formulae for the coupling kernels

In this appendix we provide analytical fitting formulae for the coupling kernels K(k,k') in Eq. (24)[*]. The analytic form is motivated by the fact that the angle averaged survey window |W(k)|2 (see Fig. 8) can be reasonably well approximated by the analytical form:

 \begin{displaymath}
\vert W(k)\vert^2 \equiv f(k) = \frac{1}{1 + \left(\frac{k}{a}\right)^2 + \left(\frac{k}{b}\right)^4}\cdot
\end{displaymath} (D.1)

Now assuming that $\vert W(\vec{k})\vert^2$ is isotropic (which certainly is not the case as seen from Fig. 9), we can find the coupling kernels K(k,k') as:
 
                         K(k,k') = $\displaystyle C\cdot\int{\rm d}\Omega_{\vec{k}}\int{\rm d}\Omega_{\vec{k'}} \vert W(\vec{k}-\vec{k'})\vert^2$  
  = $\displaystyle C\cdot\frac{8\pi^2}{kk'}\int \limits_{\vert k-k'\vert}^{k+k'}f(x)x{\rm d}x .$ (D.2)

For f(k) given by Eq. (D.1) the integral in Eq. (D.2) and the normalization constant C can be found analytically. The kernels are normalized such that

\begin{displaymath}\int K(k,k')k'^2{\rm d}k' = 1
\end{displaymath} (D.3)

is satisfied.

Depending on the values of a and b there are two different solutions.

1.
b4>4a4:

\begin{displaymath}K(k,k') = \frac{C}{kk'}\ln \left[ \frac{g(k-k')}{g(k+k')} \right] ,
\end{displaymath} (D.4)

where
$\displaystyle g(x) = \frac{\mu + b^4 + 2a^2x^2}{\mu - b^4 - 2a^2x^2} ,$     (D.5)
$\displaystyle \mu = b^2\sqrt{b^4 - 4a^4}$     (D.6)

and the normalization constant:

\begin{displaymath}C = \frac{a}{\pi\sqrt{2}\left(\sqrt{b^4+\mu}-\sqrt{b^4-\mu}\right)}\cdot
\end{displaymath} (D.7)

2.
b4<4a4:


\begin{displaymath}K(k,k') = \frac{C}{kk'} \left[ g(k+k') - g(k-k') \right] ,
\end{displaymath} (D.8)

where

\begin{displaymath}g(x) = \arctan\left(\frac{b^4 + 2a^2x^2}{b^2\sqrt{4a^4-b^4}}\right)
\end{displaymath} (D.9)

and the normalization constant:

\begin{displaymath}C = \frac{1}{\pi b \sqrt{2-\left(\frac{b}{a}\right)^2}}\cdot
\end{displaymath} (D.10)

Although the isotropy assumption is certainly not correct, the above parametric family provides a very good fit to the numerically evaluated kernels as seen in Fig. 10. The best fitting a and b for the analyzed SDSS LRG sample are 0.00457 and 0.00475, respectively.

  
Appendix E: Nonlinear model fitting. Correlated data

We find the best fitting parameters for the nonlinear model by minimizing $\chi^2$, which in the case of Gaussian errors is equivalent to finding the maximum likelihood solution. For this purpose we use Levenberg-Marquardt method as described in Press et al. (1992), where it was assumed that data values are uncorrelated. Since we are interested in the case with correlated errors, we have to make slight modifications to their implementation of the algorithm.

Using their notation, $\chi^2$ is now calculated as:

 \begin{displaymath}
\chi ^2(\vec{a})=\sum \limits_{i=1}^{N} \sum \limits_{j=1}^{...
...ot C^{-1}_{ij} \cdot \left [ {y_j - y(x_j;\vec{a})} \right ] ,
\end{displaymath} (E.1)

and the quantities $\beta_k$ and $\alpha_{kl}$ as follows:

\begin{displaymath}\beta_k=\sum \limits_{i=1}^{N} \sum \limits_{j=1}^{N} \left [...
...{-1}_{ij} \cdot \frac{\partial y(x_j;\vec{a})}{\partial a_k} ,
\end{displaymath} (E.2)


\begin{displaymath}\alpha_{kl}=\sum \limits_{i=1}^{N} \sum \limits_{j=1}^{N}
\f...
...}_{ij} \cdot \frac{\partial y(x_j;\vec{a})}{\partial a_l}\cdot
\end{displaymath} (E.3)

In the above relations Cij represents the data covariance matrix.

  
Appendix F: Goodness of fit. Correlated Gaussian data

Under the assumption that statistical fluctuations $\Delta y_i = y_i - y(x_i;\vec{a})$ $(i=1\ldots N)$ in Eq. (E.1) are Gaussian distributed, with covariance matrix Cij, one can easily derive probability density function (pdf) for the quantity $\chi^2$, and thus open up a way to estimate the goodness of fit. The $\chi^2$ goodness-of-fit estimator is usually exploited in the case of independent Gaussian variables. Here we show that calculating $\chi^2$ for the correlated Gaussian data as given in Eq. (E.1), one obtains the same result that is well known for the independently distributed case.

According to our assumption $\vec{\Delta y}$ is Gaussian distributed:

\begin{displaymath}f_{\vec{\Delta y}}(\vec{\Delta y})=\frac{1}{\sqrt{(2\pi)^N \d...
...\vec{\Delta y}^T\cdot {\bf C}^{-1} \cdot\vec{\Delta y}\right).
\end{displaymath} (F.1)

The conditional pdf of $\chi^2$ given $\vec{\Delta y}$:

\begin{displaymath}f_{\chi^2\vert\vec{\Delta y}}(\chi^2\vert\vec{\Delta y})=\del...
...- \vec{\Delta y}^T\cdot{\bf C}^{-1}\cdot\vec{\Delta y}\right),
\end{displaymath} (F.2)

and so the pdf for $\chi^2$ can be written as:
$\displaystyle f_{\chi^2}(\chi^2)=\frac{1}{\sqrt{(2\pi)^N \det {\bf C}}}\int{\rm...
...elta\left(\chi^2 - \vec{\Delta y}^T\cdot{\bf C}^{-1}\cdot\vec{\Delta y}\right).$     (F.3)

Now we define a new set of variables:

\begin{displaymath}\vec{\Delta y}' = {\bf L}^T\cdot\vec{\Delta y},
\end{displaymath} (F.4)

where ${\bf L}$ is the lower triangular matrix appearing in the Cholesky decomposition of ${\bf C}^{-1}$:

\begin{displaymath}{\bf C}^{-1} = {\bf L}\cdot{\bf L}^T.
\end{displaymath} (F.5)

Since ${\bf C}^{-1}$ can be seen as the metric tensor, we can write for the transformation of the volume elements [*]:

\begin{displaymath}\sqrt{\det {\bf C}^{-1}}{\rm d}^N\!\Delta y = \frac{{\rm d}^N\!\Delta y}{\sqrt{\det {\bf C}}}={\rm d}^N\!\Delta y'.
\end{displaymath} (F.6)

In the new frame, after changing to the spherical coordinates and integration over the angles:
$\displaystyle f_{\chi^2}(\chi^2)=(2\pi)^{-\frac{N}{2}}\cdot\frac{\Omega_N}{2}\c...
...left(-\frac{\Delta y'^2}{2}\right)\cdot\delta\left(\chi^2 - \Delta y'^2\right),$     (F.7)

where the total N-dimensional solid angle:

\begin{displaymath}\Omega_N=\frac{2\pi^{\frac{N}{2}}}{\Gamma\left(\frac{N}{2}\right)}\cdot
\end{displaymath} (F.8)

Thus the final result reads as:

\begin{displaymath}f_{\chi^2}(\chi^2)=\frac{1}{2^{\frac{N}{2}}\Gamma\left(\frac{...
...hi^2\right)^{\frac{N}{2}-1}\exp\left(-\frac{\chi^2}{2}\right),
\end{displaymath} (F.9)

which is a chi-square distribution with N degrees of freedom. Fitting P parameters (equivalent to adding P constraints), the effective number of degrees of freedom drops to $N_{{\rm eff}}=N-P$, as usual.

Now it is straightforward to calculate p-values describing the goodness of fit.

  
Appendix G: SDSS LRG power spectrum and covariance matrix

Table G.1: Measured SDSS LRG power spectrum and covariance matrix from 1000 mock catalogs.



Copyright ESO 2006