THE LARGE-SCALE THREE-POINT CORRELATION FUNCTION OF SLOAN DIGITAL SKY SURVEY LUMINOUS RED GALAXIES

The 3PCF describes the probability of finding three objects in a particular triangle configuration, compared to that of a random sample. The joint probability of finding three objects in three infinitesimal volumes dV₁, dV₂, and dV₃ is given by (Peebles 1980)

$$\begin{eqnarray} P&=&[1\hspace{-.3pt}\,{+}\,\hspace{-.3pt}\xi (r_{12})\hspace{-.3pt}\,{+}\,\hspace{-.3pt}\xi (r_{23})\hspace{-.3pt}\,{+}\,\hspace{-.3pt}\xi (r_{31})\hspace{-.3pt}\,{+}\,\hspace{-.3pt}\zeta (r_{12},r_{23},r_{31})]\nonumber\\ &&{\times}\,\hspace{-.5pt}\bar{n}^3 dV_1 dV_2 dV_3, \end{eqnarray} \tag{ 1 }$$

where $\bar{n}$ is the mean density of objects, ξ is the 2PCF, and ζ is the (connected) 3PCF:

$$\begin{equation} \xi (r_{12}) = \langle \delta (r_1)\delta (r_2)\rangle, \end{equation} \tag{ 2 }$$

$$\begin{equation} \zeta (r_{12},r_{23},r_{31})=\langle \delta (r_1)\delta (r_2)\delta (r_3)\rangle, \end{equation} \tag{ 3 }$$

where δ is the fractional overdensity of objects or the continuous field studied. The triangle sides r_ij are the distances between objects i and j in the triplet. Since the 3PCF depends upon the configuration of the three sides, it is sensitive to the shapes of spatial structures (Sefusatti & Scoccimarro 2005; Gaztañaga & Scoccimarro 2005; Marín et al. 2008). Since the ratio ζ/ξ² is both predicted and found to be close to unity over a large range of length scales even though ξ and ζ each vary by orders of magnitude (Peebles 1980), it is convenient to define the reduced 3PCF,

$$\begin{equation} Q(s,u,\theta) \equiv \frac{\zeta (s,u,\theta)}{\xi (r_{12})\xi (r_{23})+ \xi (r_{23})\xi (r_{31})+\xi (r_{31})\xi (r_{12})}. \end{equation} \tag{ 4 }$$

Here, s ≡ r₁₂ sets the scale size of the triangle, and the shape parameters are given by the ratio of two sides of the triangle, u ≡ r₂₃/r₁₂, and the angle between those two sides, $\theta =\cos ^{-1}(\hat{r}_{12}\cdot \hat{r}_{23})$, where $\hat{r}_{12}$, $\hat{r}_{23}$ are unit vectors in the directions of those sides (see Figure 1). By measuring the shape and scale dependence of Q(s, u, θ) for galaxies and comparing with that predicted for dark matter, we can constrain the galaxy bias.

Zoom In Zoom Out Reset image size

We calculate the 2PCF using the estimator of Landy & Szalay (1993),

$$\begin{equation} \xi = \frac{{\rm DD}-2{\rm DR}+{\rm RR}}{{\rm RR}}. \end{equation} \tag{ 5 }$$

Here, DD is the number of data pairs normalized by N_D × N_D/2, DR is the number of pairs using data and random catalogs normalized by N_DN_R, and RR is the number of random data pairs normalized by N_R × N_R/2, where N_D and N_R are the number of points in the data and in the random catalog, respectively. The 3PCF is calculated using the Szapudi & Szalay (1998) estimator

$$\begin{equation} \zeta =\frac{{\rm DDD}-3{\rm DDR}+3{\rm DRR}-{\rm RRR}}{{\rm RRR}}, \end{equation} \tag{ 6 }$$

where DDD, the number of data triplets, is normalized by N³_D/6, and RRR, the random data triplets, is normalized by N³_R/6. DDR is normalized by N²_DN_R/2, and DRR by N_DN²_R. These estimators optimally correct for edge effects by incorporating pair and triplet counts from a random catalog. To compute the correlation functions, we use the NPT software developed in collaboration with the Auton Lab at Carnegie Mellon University. NPT is a fast implementation of N-point correlation function estimation using multi-resolution kd-trees to compute the number of pairs and triplets in a data set. For more details and information on the algorithm, see Moore et al. (2001), Gray et al. (2003), and Nichol et al. (2006).

Since measurements of two- and three-point correlations involve counts of pairs and triplets of objects, the question arises of how to choose the bins of pair and triplet separation (Matsubara & Suto 1994; Chen & Szapudi 2005; Gaztañaga & Scoccimarro 2005; Nichol et al. 2006; Kulkarni et al. 2007; Marín et al. 2008). Small bins enable measurements with fine resolution in separation and the possible identification of features in the correlation functions but at the cost of large shot noise due to the small numbers of pairs and triplets per bin, resulting in artificially large diagonal errors and in a poorly determined covariance matrix. Large bins reduce the shot noise but tend to wash out features of the 3PCF. The optimal bin size is therefore a compromise and depends on the relative importance of signal-to-noise ratio versus resolution for the question at hand.

For the 3PCF, there is also the question of which separation parameters to bin in, that is, which triangle configurations to include in each bin. Nichol et al. (2006) and Kulkarni et al. (2007) measured the 3PCF in bins of s, u, and θ, that is, they counted triplets of galaxies that satisfy s ∈ s₀ ± Δs, u ∈ u₀ ± Δu, θ ∈ θ₀ ± Δθ for the bin centered on u₀, s₀, θ₀. In their measurement of the LRG 3PCF using the DR3 sample, Kulkarni et al. (2007) used very small resolution in Δs = ±0.1 Mpc (1% at 10 h⁻¹ Mpc) and Δθ = ±π/100 rad, and an intermediate value for Δu = ±0.1. With such small resolution the result is that Q(θ) is measured with low signal-to-noise ratio (since there are few pairs and triplets in each bin) and since Δu is not small (for u = 2, s = 10, the bin size of the second side of the triangle is 12% of the side), there are triplets counted in two or more configurations, or "bins." With this binning scheme it is not possible to measure the LRG 3PCF on larger scales. But just increasing bin size in this way would lead to have very different triangles accepted in the same bin, as shown in the right side of Figure 2.

Zoom In Zoom Out Reset image size

In this work, we use the parameters s, u, and θ to determine the centers of the triplet bins, but we construct the bins themselves directly in each of the pair separations, with Δr_ij∝r_ij, i.e., we use bins of constant Δlog r_ij (Marín et al. 2008; similar schemes were used also in Matsubara & Suto 1994; Chen & Szapudi 2005; McBride et al. 2011). As Figure 2 shows, this binning scheme groups together triplet configurations that are more similar in shape (both mathematically and colloquially) and size than do bins in s, u, and θ. As a result, we can achieve either higher resolution in triangle shape and size or higher signal-to-noise ratio for fixed resolution. Also to obtain a better signal in the 3PCF, we increase somewhat the size of the bins, which leads to a bigger correlation between the points, and therefore we do not need to measure the 3PCF for a large number of angles θ as Kulkarni et al. (2007) did.

We show in Figure 3 the results of using the different binning schemes: using the DR3 LRG sample (see Section 2.4 for a description of the sample), we compare the resultant 3PCF using the binning scheme used in Kulkarni et al. (2007, green solid lines), with the scheme used in this paper, for Δr_ij = ±0.1r_ij (open squares); it is clear that in our binning scheme Q(θ) fluctuates much less than when using the "old" scheme, it reproduces better the theoretical predictions (note that there is no plateau in the new scheme that happens in the θ > 100° in the s = 10 h⁻¹ Mpc configurations), and has smaller errors (calculated, as in Kulkarni et al. 2007, using jackknife (JK) resampling, which will be discussed in Section 4).

Zoom In Zoom Out Reset image size

We take advantage of this binning scheme to measure the 3PCF over a large range of triangle sizes and shapes: we study central configurations with s ranging from 7 h⁻¹ Mpc (which allows us to compare these measurements to previous works) to 30 h⁻¹ Mpc (on the quasi-linear regime, in order to estimate bias parameters), with u = 2 and 15 equally spaced values of θ. For other studies (outside the scope of this paper) such as finding best-fit galaxy population (HOD) parameters from fits to the 3PCF, it is necessary to carry out measurements on smaller scales. The results are limited by shot noise on small scales (due to the small spatial density of the LRG field) and by finite volume on large scales.

In Figure 4, we illustrate the sensitivity of the reduced 3PCF measurement to the choice of binning resolution. We show measurements with s = 7 and 15 h⁻¹ Mpc for the SDSS DR7-Dim LRG sample (see Section 2.4). The blue dotted, solid red, and dashed magenta curves represent results with Δr_ij/r_ij = ±0.05, 0.1, and 0.2, respectively. We can see significant bin-dependent differences in the reduced 3PCF amplitude, particularly for elongated triangle configurations (θ close to 0° or 180°). The differences and the covariance increase for larger scales. Proceeding from smaller to larger bins, the configuration dependence of Q is smoothed out, and the measurements for different values of θ become more strongly correlated. On the other hand, the signal-to-noise ratio for each measurement increases with bin width.

Zoom In Zoom Out Reset image size

After experimenting with different bin sizes using the LRG data catalogs as well as multiple realizations of artificial data resembling the LRG clustering (i.e., mock catalogs, described in Section 2.5), we found that higher bin sizes (Δr ≳ 0.1r) reduce greatly the configuration dependence of the 3PCF; lower bin sizes (higher resolution, Δr ≲ 0.05r) limited our ability to estimate the LRG 3PCF on the largest scales. Therefore, we opted to use a resolution of 14%, i.e., the bin size per triangle side to be Δr = ±0.07r, for all LRG samples we studied (see Section 2.4).

The HOD model (e.g., Jing et al. 1998; Seljak 2000; Scoccimarro et al. 2001; Berlind & Weinberg 2002; Cooray & Sheth 2002; Zheng et al. 2005) provides a parameterized prescription for the galaxy spatial distribution. The first component is a model for the distribution of dark halos. One specifies the halo mass function, n(M), the spatial correlations of the halos, and the density profiles for halos of mass M, all based on analytical models and fits to N-body simulations (e.g., Press & Schechter 1974; White & Rees 1978; Sheth et al. 2001; Tinker et al. 2008). The second component is the HOD itself, a model specifying how galaxies occupy dark halos as a function of halo mass.

The HOD is a parameterized model of the P(N|M), the probability that a halo of mass M contains N galaxies with specified properties. The mean occupation function, 〈N(M)〉, is well modeled when one separates the contribution from central galaxies, N_cen(M), which is roughly a step function in halo mass (for halos above a certain mass M_min), and the contribution of the satellite galaxies, N_sat(M), which appears to be well characterized by a power law in halo mass (Guzik & Seljak 2002; Kravtsov et al. 2004; Zheng et al. 2005), i.e., N_sat ∼ (M/M₁)^α, where M₁ characterizes the mass of halos that host satellite galaxies, and the exponent α describes the high-mass slope of the satellite occupation number.

There is no unique model for the galaxy HOD; in fact, there are many functional forms of 〈N(M)〉, as well as the number of parameters describing a particular HOD. Also, apart from the mean HOD, one needs to specify the higher moments of P(N|M). The distributions of N_cen(M) and N_sat(M) are found to be consistent with nearest-integer (Bernoulli) and Poisson distributions (Kravtsov et al. 2004; Zheng et al. 2005), and the satellite galaxies are assumed to be distributed within halos according to the Navarro–Frenk–White halo mass profiles (Navarro et al. 1997). The HOD parameters can be found using different statistics of galaxy clustering (e.g., Zehavi et al. 2004, 2005b; Kulkarni et al. 2007; Blake et al. 2008; Reid & Spergel 2009; Zheng et al. 2009; Ross & Brunner 2009).

In the case of LRGs, there is consensus (e.g., Kulkarni et al. 2007; Zheng et al. 2009; Sabiu 2009, and references therein) that these galaxies form mostly in high-mass halos with mass M ∼ a few times 10¹³–10¹⁴ h⁻¹ M_☉, while a few percent of them (in halos with M ≳ 10¹⁴ h⁻¹ M_☉) have to be satellites within halos in order to produce the 2PCF exhibited on smaller scales (Masjedi et al. 2006); the exact values differ depending on the statistic used and models of 〈N(M)〉 that were used. As we will see below in Section 2.5, the HOD model is a powerful way to create mock galaxy catalogs, and in Section 5.2, the HOD is also effective in modeling the nonlinear and the large-scale bias.

The SDSS LRGs, as selected by the algorithm developed by Eisenstein et al. (2001), are an excellent tracer of dark matter on large scales due to their high luminosity, which allow us to map them up to z ∼ 0.5 and comoving volumes up to 1.5 (h⁻¹ Gpc)³. The selection algorithm and the high success rate of the spectroscopy make (average sky completeness of 98%) the LRG sample nearly volume-limited up to z ∼ 0.35 and flux limited up to z ∼ 0.47.

We use galaxies drawn from a sample comprising 105,831 spectroscopically selected LRGs, based on SDSS Data Release 7 (DR7; Abazajian et al. 2009). The LRG galaxy selection algorithm is described in Eisenstein et al. (2001), and the DR7 LRG samples we use are drawn from the "DR7-Full" catalog described in Kazin et al. (2010), which is publicly available.¹ The DR7-Full sample spans the redshift range 0.16 < z < 0.47 and includes a range of luminosities in SDSS g band of −23.2 < M_g < −21.2, where the magnitudes have been K-corrected and passively evolved to z = 0.3, near the median redshift of the sample. Although the sample is not volume-limited, the number density is close to uniform for redshifts z ≲ 0.36 and drops off due to the flux limit at higher redshift. The sample covers a sky area of 7908 deg² (close to 20% of the sky) and includes data from both the northern and southern Galactic hemispheres.

From this main catalog, two important samples are extracted: one that includes all LRGs, but is strictly volume-limited, which we call "DR7-Dim," which has a maximum redshift of z = 0.36 and median spatial density of ∼10⁻⁴ (h⁻¹ Mpc)⁻³, and a sample of the bright LRGs, where M_g < −21.8, which we call "DR7-Bright," which reaches redshifts up to z = 0.44 and has a small spatial density ∼2.5 × 10⁻⁵ (h⁻¹ Mpc)⁻³. Both of these samples are the same as the ones presented in Kazin et al. (2010), and use data from the northern cap only, in order to have more mock catalogs when estimating the errors and covariance of our measurements (see Section 2.5 for more details).

To calculate the correlation functions, we use the same DR7-Full random catalogs as in Kazin et al. (2010), which contain 10 times the number of LRGs as the data. These random catalogs have been built having the same redshift selection function as the data, as well as the same sky completeness. On the smallest scales we ensure that we have enough random triplets in all configurations. For the purposes of calculating the 3PCF in a reasonable time, on the largest scales we lower the number of random points; when measuring the 3PCF on our largest scales, its random catalog has four times the number of data points. This is a conservative choice to ensure we have enough random points from the binning chosen, and, as was shown in Kulkarni et al. (2007), the Poisson errors on large scales are much smaller than the JK errors expected on these distances.

For comparison, we have also analyzed the earlier SDSS DR3 LRG sample used by Kulkarni et al. (2007). That sample covers 3816 deg² and contains 50,967 LRGs spanning the redshift range 0.15 < z < 0.55, similar to the DR7-Full sample except that it has a smaller sky area coverage. In this case, we use the data and random catalogs of Eisenstein et al. (2005), kindly provided by the author.

The information on these samples is summarized in Table 1.

The LRGs, i.e., giant red ellipticals, are associated with massive halos, and they are highly clustered (see, for instance, Zehavi et al. 2005a; Masjedi et al. 2006). This is illustrated in Figure 5, where we show the redshift-space 2PCF of the DR7-Dim and DR7-Bright LRGs, along with the dark matter redshift 2PCF from N-body simulations with a cosmological parameters based on the WMAP5 results (Komatsu et al. 2009, described in Section 2.5). Both LRG samples have a higher 2PCF than the dark matter, and on large scales the offset seems constant, and due to this fact we can say that the LRGs are a biased tracer of dark matter on large scales. The DR7-Bright LRGs have a higher 2PCF than the DR7-Dim LRGs, suggesting that they are associated with more massive halos (Zheng et al. 2009).

Zoom In Zoom Out Reset image size

To help analyze and interpret the 3PCF results, we use mock catalogs drawn from N-body simulations. The mock catalogs serve two functions: (1) to measure the 3PCF of dark matter, with which the galaxy 3PCF measurements can be compared to infer the bias and (2) to determine the expected error covariance matrix for the 3PCF from the variance among a large number of mock catalogs. We use different simulations tailored to these different functions.

First, to study the dark matter two- and three-point correlation functions, we use six N-body simulations carried out by C. Sabiu et al. (2010, in preparation), kindly provided by the author, using the GADGET code. The cosmological model is a spatially flat ΛCDM (Λ-cold dark matter) universe, with matter density parameter Ω_m = 0.27, baryon density Ω_b = 0.045, power spectrum amplitude σ₈ = 0.8, reduced Hubble parameter h = 0.7, and initial power spectrum index n_s = 0.95. The simulation box volume is 1 (h⁻¹ Gpc)³ and contains 512³ dark matter particles with a particle mass of 5 × 10¹¹ h⁻¹ M_☉, with a spatial softening of 60 h⁻¹ kpc, and we use the output at z = 0.3, the mean redshift of the LRGs, to compare their correlation functions.

Since using all the dark matter particles for the 3PCF calculations would take a prohibitively long time, we randomly select particles from the simulation, to a degree where we are sure that the dilution will not affect the results. For the 3PCF measurement, this depends on scale, and we use between 5% and 0.5% of the particles (a higher fraction for smaller scales); this fraction is enough to obtain reliable measurements of the 2PCF and 3PCF on the scales studied in this work. To obtain the spatial coordinates in the redshift-space dark matter catalog, we use the long-distance approximation in the z-direction (see, e.g., Bernardeau et al. 2002) to account for the effect of peculiar velocities on the measurement of the correlation functions.

We use the dark matter halo catalogs obtained from these simulations to study redshift distortions in the estimation of the bias parameters (see Section 5.2 below), found using a friends-of-friends algorithm (Davis et al. 1985) with a linking length of b = 0.2 times the mean interparticle separation. To populate halos with mock LRGs, we use the best-fit HOD parameters from fits to the projected 2PCF by Zheng et al. (2009) in their five-parameter model. We assign mock galaxy positions, based on halo profiles, and peculiar velocities (based on a Gaussian distribution, which depends on the mass of the host halo) of galaxies inside halos using a code kindly provided by Jeremy Tinker.

Second, to estimate the errors in the 3PCF estimates, we use the covariance among public mock catalogs from the LasDamas simulations (McBride et al. 2011).² As will be discussed in Section 4, we prefer this error estimation method to the use of JK resampling. For the LasDamas simulations, the cosmological model is similar but not identical to that above: flat, ΛCDM with Ω_m = 0.25, Ω_b = 0.04, σ₈ = 0.8, h = 0.7, and n_s = 1. For this set of parameters, the LasDamas team has carried out N-body simulations in boxes of different sizes to model galaxies of different luminosities; for the LRGs, they use 40 boxes with volume of 2400³(h⁻¹ Mpc)³, each one containing 1280³ dark matter particles with a particle mass of 4.57 × 10¹¹ h⁻¹ M_☉.

These simulations are used to construct mock galaxy catalogs by placing galaxies in dark matter halos using an HOD model (Berlind & Weinberg 2002), with HOD parameters fit from the observed SDSS LRG two-point clustering. From this main set they create two types of galaxy mock catalogs, one that resembles the LRG DR7-Dim catalog, with z_max = 0.36, and one set that resembles the Bright LRG sample with z_max = 0.44. The mocks reproduce the SDSS DR7 geometry, and they include redshift-space distortions from peculiar velocities. These galaxy mocks have been shown to reproduce well the form of the DR7 LRG 2PCF on large scales (Kazin et al. 2010), and we show below that they are in general agreement with the 3PCF measurements as well. Since our DR7-Dim and DR7-Bright samples use only the northern cap of the SDSS, in each LasDamas box it is possible to extract four mock LRG catalogs; in this way, we have 160 mock catalogs for each sample.

Figure 6 shows the reduced 3PCF of the DR7 LRG samples on intermediate scales, for s = 7, 10, and 15 h⁻¹ Mpc. For clarity, we show diagonal error bars only for the DR7-Dim sample, calculating using the 3PCF variance between the LasDamas mocks. For comparison, we also show the 3PCF for the DR3 LRG sample. The solid black curves show the 3PCF of the dark matter from the N-body simulations.

The reduced 3PCF shows the general shape dependence expected from gravitational instability theory (Bernardeau et al. 2002): the amplitude is higher for elongated triangle configurations (θ ∼ 0°, 180°), reflecting anisotropic velocity flows along density gradients. The configuration dependence is less pronounced on smaller scales, where quasi-virialized, more isotropic flows and structures dominate. Since these measurements are in redshift space, the configuration dependence also reflects the effects of redshift distortions due to peculiar velocities.

Since LRGs are biased tracers of the dark matter, the 3PCF of the LRGs differs from that predicted for the dark matter. On small scales, the reduced 3PCF amplitude of the LRGs is smaller by a nearly constant factor from that of the dark matter. On larger scales, the nature of the 3PCF bias changes: there is a shape dependence in the offset between the LRG and dark matter amplitudes, with the galaxy 3PCF showing less configuration dependence than that of the dark matter.

Comparing the DR7-Full results to the DR3 measurements, which cover half the volume of the DR7 sample, we find general consistency in the 3PCF amplitude on these intermediate scales. For s = 15 h⁻¹ Mpc, the DR3 amplitude is slightly lower than that of the DR7 sample for all configurations; this is a consequence of the smaller volume sampled to estimate the correlation functions: as observed in other studies (see Marín et al. 2008, and references therein), the 3PCF is more sensitive to the sample volume than the 2PCF.

In general, there are no important differences between the DR7-Full and the volume-limited DR7-Dim sample results. There are differences between these and the DR7-Bright sample: the reduced 3PCF for the bright LRGs fluctuates much more and it deviates from the DR7-Dim sample results for the elongated triangles (large θ) for the s = 7 and 15 h⁻¹ Mpc configurations, and in the collapsed triangles (small θ) for the s = 10 h⁻¹ Mpc configurations. Since we have used a binning more tailored to measure the DR7-Dim correlations on the DR7-Bright (a much less dense sample) measurements, we expect more fluctuations due to poorer statistics and effects of large structures (Nichol et al. 2006); we will explore the variance of bright LRGs in Section 3.3.

In Figure 7, we show the reduced 3PCF for LRGs on large scales, for s = 20 and 30 h⁻¹ Mpc, with u = 2.0. As expected from second-order perturbation theory, the reduced 3PCF shows a stronger configuration dependence on these scales compared to smaller scales, with a dramatic increase in the difference in amplitude between elongated and rectangular configurations at s ⩾ 20 h⁻¹ Mpc. That behavior is seen for both the dark matter and the LRG samples. For s = 30 h⁻¹ Mpc, the LRG 3PCF dependence in θ is more asymmetrical between collapsed and elongated configurations, as it happens with the dark matter 3PCF. The errors are greatly increased as well, and there are strong anticorrelations, i.e., negative values of Q(θ), for the rectangular configurations; this is again a consequence of the filamentary shape of the large-scale structure.

The large-scale 3PCF of the DR7-Bright galaxies is similar to the DR7-Dim 3PCF, with less significant fluctuations on larger scales. On these scales the finite-volume effects can be observed more clearly in the 3PCF: comparing the DR7-Full and DR7-Dim LRG samples, we can see that in general they agree within the error bars, and the small differences are due to a combination of finite-volume effects and biasing effects (the DR7-Full sample includes brighter galaxies at large redshifts, which have a different clustering bias than the DR7-Dim sample). Comparing the DR3-Full results compared to DR7-Full and even DR7-Dim, we can see that the DR3-Full 3PCF deviates from the DR7 samples on the elongated configurations (larger scales).

In Figure 8, we compare the 3PCF measurements for the DR7 LRGs, as shown previously in Sections 3.1 and 3.2 to the ones obtained from the LasDamas mocks, with dashed lines for redshift-space measurements, solid lines for real-space measurements. Here we present a representative set of values, with s = 7, 15, and 30 h⁻¹ Mpc. In the top panels we show the comparison for DR7-Dim galaxies; in general the agreement is good, and usually within the 1σ error bars.

Zoom In Zoom Out Reset image size

In the lower panels, we compare DR7-Bright galaxies and the mock 3PCFs. It can be seen that the variance of the bright mocks is larger compared to the dim ones, and the agreement between the mocks and the data is not as good, especially on smaller scales for the elongated configurations. Although due to their low density less signal-to-noise ratio and more fluctuations in Q(θ) are expected, the variance from the mocks should reflect this, and that is not the case for the smaller scales. Since the mocks are built in order to fit the measured two-point statistics, it can happen that the mock misses the higher-order statistics, and for future modeling these statistics have to be considered.

Note also the differences between the real- and redshift-space measurements in the LasDamas mocks. Both have a dependence on shape, but on small scales the real-space 3PCF has a more pronounced shape dependence. On bigger scales, especially around s = 15 h⁻¹ Mpc, they are very similar within the errors. Comparing these differences to the ones found by Marín et al. (2008), there are two causes of this: first, the LRGs are galaxies that inhabit massive halos (e.g., Kulkarni et al. 2007; Zheng et al. 2009, among others), which are in the center of the gravitational potential wells of the large matter structures, and therefore they are less affected peculiar velocities. Second, the 3PCF on larger scales is much more affected by the morphology of the structures than from the peculiar velocities, and therefore on the largest scales the distortions of the reduced 3PCF will be smaller than on the small scales. There is also a resolution (or binning) component to this; in very thin bins we might be able to detect better the distortions on large scales. We will use this similarity between the real- and redshift-space 3PCFs below in Section 5.2 to estimate real-space bias parameters from the 3PCF.

For the JK and variance methods, the diagonal errors are given by

$$\begin{equation} \sigma _Q^2 = \frac{E_{{\rm method}}}{N}\sum _{i=0}^N(Q_i-\bar{Q})^2. \end{equation} \tag{ 7 }$$

For the JK method, N = N_JK, the number of subvolumes, and E_JK = N_JK − 1; for the variance among mocks, N is the number of independent realizations, and E_var = 1.

To have a qualitative estimate of the differences in the diagonal errors between the two methods, we calculate, for each sample, the average of the reduced 3PCF for all angles in a particular configuration, i.e., we calculate $\bar{Q}(s=15,u=2,\theta)$, the mean of Q(θ) for triangles with s = 15.0  h⁻¹ Mpc, u = 2.0, as a function of the number N of samples. We do this for three values of θ, one corresponding to an almost collapsed triangle (θ = 18°), another for a rectangular configuration (θ = 90°), and one for an almost completely extended configuration (θ = 160°). We chose this scale since the 3PCF here is well behaved and does not have big diagonal errors, but also is large enough to make visible the effects of the JK resampling in the 3PCF (i.e., to test the "independent regions" hypothesis).

The results are shown in Figure 9. In the left panel, we compare $\bar{Q}(s=15,u=2,\theta)$ for mocks and for JK subsamples. The average values converge very fast with N, and the fact that the means of the JK samples and the mocks are close to each other (differences at the level of 5% on these scales) is a good signal that we can use the mocks for the covariance measurements. In the right panel, we show the variance (standard deviation) of Q(s, u, θ) as a function of N, for the different methods. We notice that while the mock's variance seems to converge at large N, the JK errors tend to increase when we increase the number of subsamples, with differences reaching close to 50% on the elongated configurations.

Zoom In Zoom Out Reset image size

For N-point correlation functions, measurements for different configurations are strongly correlated, therefore the off-diagonal errors must be included in model fitting and parameter estimation. Having as many mock realizations or subsamples as possible help us define better the covariance matrix. Since we cannot have an infinite number of subsamples, if our number of data points is larger than the number of mocks used, when calculating the inverse of the covariance matrix for constraining bias parameters, we need to take into account the singular and numerically ill-defined matrix modes. We describe our approach later in Section 5.

Figure 10 shows the correlation matrices of the 3PCF for configurations with s = 15 h⁻¹ Mpc, u = 2, and varying θ, comparing the JK (top row) and mock-covariance (middle row) methods, for different numbers of JK subvolumes and mock realizations. For both methods, the off-diagonal correlation coefficients increase with increasing N, which goes from 10 to 160 subsamples/realization.

Zoom In Zoom Out Reset image size

For these configurations, the off-diagonal terms tend to converge when N ⩾ 40. As can be seen in the bottom row of Figure 10, where we show the ratio between the off-diagonal terms of the covariance matrices, |C^mocks_ij(N)/C^JK_ij(N)|, both methods yield similar results except for the terms relating very collapsed and very elongated configurations, where the JK configurations have higher correlations. This enhanced covariance for the JK methods is due to the fact that this method uses only one "realization"; while on small scales the regions can be effectively independent, on large scales (i.e., when examining large triangles) these regions are correlated, and then the JK assumption breaks down.

Even though for large N the covariance matrix converges for the JK method, and we would be able to extract more information from the covariance matrix, we would be introducing artificial correlations since the bigger N is, the smaller and "less independent" the volumes of the subsamples are. From this analysis we can conclude that the JK method leads to significant biases in the estimation of the covariance matrix for the 3PCF, and therefore using adequate mock catalogs is preferable. In our analysis of large-scale clustering (for bias and cosmological parameters) we will use the 3PCF error covariance matrices from the variation among the 160 mock catalogs.

We compare the dark matter (from the N-body simulations described in Section 2.5) and LRG 2PCF and 3PCF in order to constrain bias parameters using relations (11) and (12); in this case assuming these relations are valid in redshift space, i.e., ξ^r-space → ξ^z-space and Q^r-space → Q^z-space. We use all configurations with s ⩾ 10 h⁻¹ Mpc, and s ⩾ 15 h⁻¹ Mpc with u = 2 to measure the joint likelihood of the bias parameters and σ₈.

As described by Gaztañaga & Scoccimarro (2005) and used in Gaztañaga et al. (2005) and Marín et al. (2008), we minimize

$$\begin{equation} \chi ^2=\sum _{i=1}^{i=2N_b}\sum _{j=1}^{j=2N_b}\Delta _iC_{ij,{\rm SVD}}^{-1}\Delta _j, \end{equation} \tag{ 13 }$$

where N_b is the number of configurations used (N_b = 60 for triangles with s ⩾ 10 h⁻¹ Mpc and N_b = 45 for triangles with s ⩾ 15 h⁻¹ Mpc). Since we use both two- and three-point correlations we have

$$\begin{equation} \Delta _i = (\xi (r_3)^{{\rm obs}}_i - \xi (r_3)^{{\rm model}}_i)/\sigma _{\xi (r_3)i}, \hbox{ for } i\le N_b, \end{equation} \tag{ 14 }$$

$$\begin{equation} \Delta _i = (Q^{{\rm obs}}_i-Q^{{\rm model}}_i)/\sigma _{Q(i)}, \hbox{ for } i>N_b, \end{equation} \tag{ 15 }$$

where ξ(r₃)^model and Q^model are given by Equations (9) and (12) for the galaxies. The matrix C_{ij, SVD} is the normalized covariance matrix. This covariance matrix is calculated from the LasDamas LRG mock catalogs; $\sigma _{\xi (r_3)i}$ and σ_Q(i) are the uncertainties in ξ(r₃)_i and Q(i), respectively, and we assign them the error bars from the mocks. As mentioned in Section 4.2, we need to use a number of modes smaller than the number of mocks to avoid adding or subtracting noise from places where the inverse of covariance matrix is singular. Using N_mocks = 160, in principle we do not have that problem, but still we are subject to numerical rounding errors when inverting the covariance matrix. Therefore, we follow the suggestion of Gaztañaga & Scoccimarro (2005): when using a limited number of mocks, we recalculate the inverse of the covariance matrix with the highest modes obtained from a singular value decomposition (SVD), and use the modes where its SVD eigenvalues are $\lambda < \sqrt{2/N_{{\rm mocks}}}$. This in practice lowers the number of degrees of freedom when estimating χ²: they will be only equal to the number of eigenmodes used minus the parameters constrained.

The results are shown in Figures 11 and 12, where we use the covariance matrix from the 160 LasDamas mocks for each sample to constrain c₁, c₂, and σ₈. In Figure 11 we show the results for the DR7-Dim sample, and in Figure 12, the results for the DR7-Bright sample. In the top panels of each figure we show constraints on the bias parameters and σ₈, for two sets of configurations: solid lines use triangles with s ⩾ 10 h⁻¹ Mpc and dotted lines use only the largest triangles s ⩾ 15 h⁻¹ Mpc. The contours correspond to Δχ² = 1.0, 2.3, and 6.2 (corresponding to 1σ limits for one-parameter fit, 1σ limit for two-parameter fit, and 2σ limit for two parameters, respectively).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We can see that in both samples a zero nonlinear bias is excluded at 2σ significance. We can see that, as expected, DR7-Dim galaxies have a lower linear bias c₁ than the DR7-Bright, and there is a bigger difference in the nonlinear bias parameter c₂. Using different sets of triangles changes not only the area of the confidence contours, but also the best-fit results; going to smaller scales put us in the direction of the nonlinear regime, where the relations in Equations (11) and (12) are less valid. Observe that the sizes of the confidence intervals are not small: this is an effect of using a fraction of the total number of eigenmodes due to the limited number of mock catalogs, and also that our errors increase on large scales.

With respect to the constraints on σ₈, in general they agree with the WMAP5 results (Komatsu et al. 2009), but they are more dependent on which configurations are used.

The best-fit values and 1σ errors are also presented in Table 2. We can see that the values of χ²/dof are robust for the DR7-Dim sample and show that the model is a good fit for the data. For the highly biased DR7-Bright galaxies the χ²/dof values are somewhat worse; we believe this is primarily due to the higher uncertainties and systematics when measuring the 3PCF for the brighter galaxies rather than an inadequacy of the bias model used.

Table 2. Best-fit Bias Parameters in Redshift Space

Sample	smina	Fit	c1	c2	σ8	$\frac{\chi ^2}{\hbox{dof}}$/dofb
DR7-Dim	10	2PCF and 3PCF, three-parameter	1.92+0.2− 0.1	0.38+0.02− 0.08	0.82+0.05− 0.13	1.31/(70 − 3)
DR7-Dim	15	2PCF and 3PCF, three-parameter	2.02+0.2− 0.1	0.58+0.12− 0.08	0.76+0.05− 0.07	0.95/(57 − 3)
DR7-Bright	10	2PCF and 3PCF, three-parameter	1.92+0.2− 0.1	0.55+0.12− 0.15	0.84+0.04− 0.05	1.76/(68 − 3)
DR7-Bright	15	2PCF and 3PCF, three-parameter	1.95+0.2− 0.1	0.27+0.15− 0.08	0.87+0.06− 0.07	1.44/(55 − 3)
DR7-Dim	15	2PCF and 3PCF, σ8 = 0.9	1.72+0.02− 0.02	0.31+0.05− 0.05	⋅⋅⋅	0.79/(46 − 2)
DR7-Dim	15	2PCF and 3PCF, σ8 = 0.8	1.93+0.02− 0.02	0.45+0.05− 0.05	⋅⋅⋅	0.73/(47 − 2)
DR7-Dim	15	3PCF, two-parameter	1.83+0.12− 0.1	0.32+0.12− 0.05	⋅⋅⋅	0.85/(29 − 2)
DR7-Bright	15	2PCF and 3PCF, σ8 = 0.9	1.95+0.02− 0.02	0.32+0.1− 0.1	⋅⋅⋅	1.32/(54 − 2)
DR7-Bright	15	2PCF and 3PCF, σ8 = 0.8	2.18+0.02− 0.02	0.51+0.1− 0.1	⋅⋅⋅	1.28/(52 − 2)
DR7-Bright	15	3PCF, two-parameter	2.05+0.1− 0.1	0.48+0.1− 0.1	⋅⋅⋅	1.01/(31 − 2)
DR7-Dim	⋅⋅⋅	Z08 HOD fit, σ8 = 0.8	2.22	0.55	⋅⋅⋅	⋅⋅⋅
DR7-Bright	⋅⋅⋅	Z08 HOD FoF fit, σ8 = 0.8	2.36	0.79	⋅⋅⋅	⋅⋅⋅

Notes. ^aAll triangles used with u = 2.0 where the first side s ⩾ s_min, in h⁻¹ Mpc. ^bDegrees of freedom (dof) are equal to the number of modes used minus parameters constrained.

Download table as:  ASCIITypeset image

In the lower two panels we show the LRG 3PCF along with the dark-matter-biased 3PCF, using the best-fit bias values mentioned above. For the DR7-Dim sample the fits are quite good, and there are no significant differences between the fits from using different triangles (solid line corresponds to fit using s ⩾ 10 h⁻¹ Mpc triangles whereas dotted lines use only s ⩾ 15 h⁻¹ Mpc triangles). For the brighter sample, the fits are less good for configurations with θ > 100°.

We explore how the constraints on the bias parameters change if σ₈ is fixed (instead of being marginalized), and compare results to the analytical bias prediction from the HOD models. In Figure 13 we show these results, using only s ⩾ 15 h⁻¹ Mpc configurations. On the left side, the vertical ellipses correspond to constraints on the bias parameters using both the 2PCF and 3PCF, using a model where σ₈ is fixed. As expected, a lower σ₈ means a higher c₁. The constraints on c₁ are tight, due to the lower errors in ξ, leaving mostly the uncertainty on c₂ = b₂/b₁.

Zoom In Zoom Out Reset image size

We compare these bias parameters with that expected from the analytical predictions of the best-fit HOD from Zheng et al. (2009, using projected 2PCF only) for both linear and nonlinear biases in the last row of Table 2, and we have good agreement. Our error bars, in the case of σ₈ = 0.8, confirm that a zero nonlinear bias parameter is rejected at the 2σ level. On the right side, we show the bias constraints using the 3PCF only. As can be seen from the results shown in Table 2, there is a slight disagreement in the best-fit values using only the 3PCF versus using both the 2PCF and 3PCF.

We explore here the validity of using redshift-space correlations to calculate galaxy bias parameters. When mapping the galaxies in the sky, we use the redshift of the galaxy as an indicator of distance. But since galaxies (and dark matter) are subject to gravitational forces from the surrounding structures, peculiar velocities are added to the cosmological redshift, distorting our mapping of the large-scale structure and, as a consequence, the correlation functions are affected as well (see Hamilton 1998, for a review).

In the case of the 3PCF, there are differences between the real- and redshift-space 3PCF of galaxies and dark matter (e.g., Gaztañaga & Scoccimarro 2005; Marín et al. 2008), and they will affect bias parameters as well. Marín et al. (2008), using mock galaxy catalogs, found good agreement, but not perfect, on the bias parameters in the L_* range, between using real- or redshift-space correlations, using s ∼ 10 h⁻¹ Mpc configurations, and poorer agreement for brighter galaxies. That was caused by a combination of studying the bias using smaller scales (in the mildly nonlinear regime) and by the small size of the box (only 200 h⁻¹ Mpc on the side) for studying these galaxies which have a low spatial density.

As shown in Figure 8, on large scales (with s ⩾ 15 h⁻¹ Mpc), for the mocks, the difference between the redshift- and real-space 3PCFs is small and falls within the errors. In the case of the dark matter, there are differences, but because of the large bias, the mocks indicate that the LRG 3PCF will not show significant differences between the redshift- and real-space measurements on large scales.

In Figure 14, we show the differences in the correlation functions of dark matter and mock LRG galaxies from the N-body simulations by C. Sabiu et al. (2010, in preparation), and how the bias parameters and σ₈ constraints change from using real- and redshift-space correlations. Note that on small scales there are significant differences between redshift- and real-space 3PCFs, but on larger scales, they are very similar. We have to keep in mind that these similarities occur when using this particular binning scheme, and might not be necessarily true with smaller binning 3PCF resolutions. For the 2PCF, the differences are small, but significant. For the bias constraints, note that real-space measurements have a higher linear and nonlinear bias, and usually the best-fit value of σ₈ is lower in the real-space fits.

Zoom In Zoom Out Reset image size

There have been analytical attempts to relate the real- and redshift-space three-point function in Fourier space (i.e., bispectrum, Scoccimarro et al. 1999; Smith et al. 2008), but there are practically no attempts of this in the configuration space, since nonlinearities and mixing of large- and small-scale modes make the calculations very complicated unless one takes very simplified cases or does not take full advantage of the shape dependence of the three-point statistics (e.g., Pan & Szapudi 2005; Sefusatti et al. 2006). In our case, we take an empirical approach.

Using the similarity between the mock LRG 3PCF in real and redshift space for the large scales, we estimate constraints on the bias and σ₈ parameters using the real-space dark matter correlations, and identifying the real-space 3PCF with the redshift-space 3PCF. In the case of the 2PCF, since the differences between real- and redshift-space correlations are significant, we will use the first-order linear approximation in the case of redshift distortion (Kaiser 1987). Following Pan & Szapudi (2005),

$$\begin{equation} \xi ^{z\scriptsize {\hbox{-space}}}_{{\rm gal}}=\left(1+\frac{2}{3}f+\frac{1}{5}f^2\right)c_1^2\left(\frac{\sigma _8}{0.8}\right)^2\xi _{{\rm dm}}^{r\scriptsize {\hbox{-space}}}, \end{equation} \tag{ 16 }$$

where f = Ω^0.6_m/c₁; we tested this prescription for our HOD mocks and it works well on the largest scales r > 10 h⁻¹ Mpc. The results are shown in Figure 15, for both the LRG-Dim and LRG-Bright samples, using triangles with s ⩾ 15 h⁻¹ Mpc. The best-fit values are summarized in Table 3. In general, we see the same trend that was seen in Figure 14: both c₁ and c₂ are higher in real space, and we get reasonable constraints on σ₈. We notice that here the constraints are in good agreement with the values obtained from the HOD fit of by Zheng et al. (2009; shown in Table 2).

Zoom In Zoom Out Reset image size

Table 3. Best-fit Bias Parameters in Real Space

Sample	smina	Fit	c1	c2	σ8	$\frac{\chi ^2}{\hbox{dof}}$/dofb
DR7-Dim	10	2PCF and 3PCF, three-parameter	2.51+0.15− 0.1	0.82+0.12− 0.04	0.75+0.03− 0.02	1.17/(70 − 3)
DR7-Dim	10	2PCF and 3PCF, three-parameter	2.38+0.2− 0.1	0.92+0.1− 0.1	0.77+0.02− 0.05	0.79/(57 − 3)
DR7-Bright	15	2PCF and 3PCF, three-parameter	2.22+0.2− 0.1	0.75+0.2− 0.1	0.91+0.1− 0.1	1.65/(68 − 3)
DR7-Bright	15	2PCF and 3PCF, three-parameter	2.45+0.45− 0.25	0.79+0.4− 0.1	0.85+0.1− 0.1	1.22/(55 − 3)

Notes. ^aAll triangles used with u = 2.0 where the first side s ⩾ s_min, in h⁻¹ Mpc. ^bDegrees of freedom (dof) are equal to the number of modes used minus parameters constrained.

Download table as:  ASCIITypeset image

It is possible to use the bias parameter constraints to estimate best-fit values on the LRG HOD parameters. In the Appendix, we show how it is possible to extract the large-scale bias parameters from the HOD, where only the average of the HOD function 〈N(M)〉 is needed to obtain the bias parameters b^h₁ and b^h₂. Given the uncertainties in c₁, c₂ from the 3PCF, and for a fixed σ₈, it is possible in principle to obtain constraints on HOD parameters (Sefusatti & Scoccimarro 2005) from large-scale bias parameters. We exemplify this using the LRG DR7-Dim sample.

There are many models for the mean 〈N(M)〉, with different number of parameters. It is worth mentioning that even though the meaning of common HOD parameters such as M_min, M₁, and α is similar for different functional forms of 〈N(M)〉, the best-fit parameters can change significantly depending on the function used (Zheng et al. 2005). For the sake of simplicity, and to avoid marginalizing over many parameters, we use a simple model that uses three parameters, with a soft function for N_cen(M) (Sefusatti et al. 2006; Kulkarni et al. 2007):

$$\begin{equation} \langle N(M) \rangle =N_{{\rm cen}}(M)+N_{{\rm sat}}(M), \end{equation} \tag{ 17 }$$

$$\begin{equation} N_{{\rm cen}} = \exp (-M_{{\rm min}}/M), \end{equation} \tag{ 18 }$$

$$\begin{equation} N_{{\rm sat}} = \exp (-M_{{\rm min}}/M)\left(\frac{M}{M_1}\right)^\alpha. \end{equation} \tag{ 19 }$$

We effectively constrain the HOD parameters M₁ and α, since we fix the spatial density of model galaxies to the mean density of the sample, and in practice, this constrains M_min.

We create a grid with different values of M₁ and α, where we calculate the bias parameters b^h₁ = c₁ and b^h₂ = c₁c₂ (see the Appendix), and then we relate and use the results for the constraints on the bias parameters c₁ and c₂ of the DR7-Dim sample when we fix σ₈ = 0.8. This gives us a likelihood on the values of M₁ and α HOD parameters.

In Figure 16, we show the constraints on M₁ and α using the constraints on c₁ and c₂ for the DR7-Dim LRG sample, using s ⩾ 15 h⁻¹ Mpc configurations, for fixed σ₈ = 0.8. We notice the constraints are not very strong, especially on α. The best-fit HOD parameters, M₁ = 3.2 × 10¹⁴ h⁻¹ M_☉, and α = 2.35, and the constraints on M₁ are consistent with what is expected from fits to two-point correlations, and although our best-fit value for α is high, it is consistent with the findings of high α for the LRG HOD (Blake et al. 2008; Zheng et al. 2009), and higher from the estimates from 3PCF fits by Kulkarni et al. (2007). Note that the strongest constraints are on M₁, and comparing the shape of these contours with the contours of Figure 13, we can conclude that M₁, the "average" mass of a halo containing at least two LRGs (a central and a satellite) is mostly controlled by the linear bias c₁, and α is more sensitive to the c₂ constraints.

Zoom In Zoom Out Reset image size

The distribution of dark matter halos is biased with respect to the dark matter (Press & Schechter 1974), therefore the galaxies, formed inside these halos, are biased tracers as well. As mentioned above, in the HOD model we proceed in two steps: first, we described how dark matter halos are distributed in the universe, and second, how galaxies populate these halos.

To describe the mass function of dark matter halos (i.e., the spatial density of halos as a function of mass) we use the Sheth & Tormen mass function (Sheth et al. 2001; Sheth & Tormen 2002):

$$\begin{equation} n_h(M) = \frac{\bar{\rho }_m}{M^2}\left|\frac{d\ln \sigma }{d\ln M}\right| f(\nu), \end{equation} \tag{ A1 }$$

where $\bar{\rho }_m\equiv \Omega _M\rho _c$ is the average matter density in the universe at the epoch of observation, ν = δ_c/σ(M, z), with δ_c = 1.69 the threshold for growth of linear fluctuations, and σ, the rms variance of a sphere of radius R(M) at redshift z:

$$\begin{equation} \sigma ^2(R,z)= \int \frac{d^3k}{(2\pi)^3} \left| W(k,R) \right|^2P(k,z). \end{equation} \tag{ A2 }$$

Here, P(k, z) is the linear matter power spectrum and W(k, R) is the Fourier transform of a top-hat window of radius R. The function f(ν) is motivated by the ellipsoidal collapse model (Sheth & Tormen 2002) and has the form

$$\begin{equation} f(\nu) = A\sqrt{\frac{2a\nu ^2}{\pi }}[1+(a\nu ^2)^{-p}]e^{-a\nu ^2/2}, \end{equation} \tag{ A3 }$$

where for normalization A = 0.322, p = 0.3, and a = 0.707.

Knowing the probability function P(N|M) of how galaxies populate halos as a function of mass, we can infer the clustering at any scale. The HOD model allow us, knowing the mean occupation number <N(M) >, in the large-scale limit to estimate the bias parameters:

$$\begin{equation} b_i = \frac{1}{n_g}\int _{M_{{\rm min}}}^{\infty } {\it dM} n_h(M)b_i^h(M)\langle N_g(M)\rangle, \end{equation} \tag{ A4 }$$

where i = 1, 2 represent the linear and first nonlinear bias term, n_h(M) the dark matter halo mass function, b^h_i(M) is the halo bias and <N_g(M) > is the mean number of galaxies for a halo of mass M, and n_g is the galaxy number density, calculated as

$$\begin{equation} n_g = \int _{M_{{\rm min}}}^{\infty } {\it dM} n_h(M)\left<N(M)\right>. \end{equation} \tag{ A5 }$$

The dark matter halo bias parameters are given by (Scoccimarro et al. 2001)

$$\begin{equation} b_1^h(M) = 1+\epsilon _1+E_1, \end{equation} \tag{ A6 }$$

$$\begin{equation} b_2^h(M) = \frac{8}{21}\left(\epsilon _1+E_1\right)+\epsilon _2+E_2, \end{equation} \tag{ A7 }$$

where the coefficients are given by

$$\begin{equation} \epsilon _1 = \frac{a\nu ^2-1}{\delta _c}, \end{equation} \tag{ A8 }$$

$$\begin{equation} \epsilon _2 = \frac{a\nu ^2(a\nu ^2-3)}{\delta _c^2}, \end{equation} \tag{ A9 }$$

$$\begin{equation} E_1 = \frac{2p}{\delta _c(1+(a\nu ^2)^p)}, \end{equation} \tag{ A10 }$$

$$\begin{equation} \frac{E_2}{E_1} = \frac{1+2p}{\delta _c}+2\epsilon _1. \end{equation} \tag{ A11 }$$

It is worth emphasizing that the bias relations mentioned in Equation (A4) are only valid on large scales, in the quasi-linear regime, and they are scale independent. On smaller scales the bias becomes very scale dependent and higher moments of the P(N|M) function need to be known.