The Three-Point Correlation Function in Cosmology

With the advent of high-quality surveys in cosmology the full three-point correlation function will be a valuable statistic for describing structure formation models. It contains information on cosmological parameters and detailed halo properties that cannot be extracted from the two-point correlation function. We use the halo clustering model to analytically calculate the three-point correlation function (3PCF) for general cosmological fields. We present detailed results for the configuration dependence of the 3-dimensional mass and galaxy distributions and the 2-dimensional cosmic shear field. We work in real space, where higher order correlation functions on small scales are easier to measure and interpret, but halo model calculations get rapidly intractable. Hence we develop techniques for accurate evaluations of the 1, 2 and 3-halo contributions to the 3PCF. The 3PCF violates the hierarchical ansatz in both its scale and configuration dependence. We study the behavior of the coefficient Q in the hierarchical expansion from large, quasilinear scales down to about 20 kpc. We find that the nonlinear 3PCF is sensitive to the halo profile of massive halos, especially its inner slope. We model the distribution of galaxies in halos and show that the 3PCF of red galaxies has a weaker configuration and scale dependence than the mass, while for blue galaxies it is very sensitive to the parameters of the galaxy formation model. The 3PCF from weak lensing on the other hand shows different scalings due to projection effects and a sensitivity to cosmological parameters.

In this paper, we focus on developing a theoretical model of the 3PCF in real space for 3D and 2D cosmological fields. In practice, on small scales the 3PCF would be easier to measure from observational data over the bispectrum, since it does not require the Fourier transform for the survey data that usually have a complicated geometry of data fields.
Theoretical models of the weakly nonlinear 3PCF have been well studied in the literature based on perturbation theory (Fry 1984b;Gaztañaga & Bernardeau 1998;Frieman & Gaztañaga 1999; Barriga & Gaztañaga 2002). Perturbation theory can describe properties of the dark matter and galaxy clustering on large scales > ∼ 10 h −1 Mpc and predict that the 3PCF depends on the shape of triangle configuration and, as a result, contains information of the primordial power spectrum and the galaxy biasing. Historically, the pioneering measurement of the galaxy 3PCF done by Peebles & Groth (1975) (also Groth & Peebles 1977) proposed the "hierarchical form", ζ(r12, r23, r31) = Q[ξ(r12)ξ(r23) + ξ(r23)ξ(r31) + ξ(r31)ξ(r12)], with the constant Q ≃ 1.3. However, subsequent work has revealed that the measured 3PCF does not obey the hierarchical form rigorously and the large-scale amplitudes can be explained by the perturbation theory results of the cold dark matter (CDM) models, if the biasing relation is correctly taken into account for the analysis Frieman & Gaztañaga 1999).
On the other hand, a quantitative theoretical model of the 3PCF in the strongly nonlinear regime is still lacking except for studies relying on N -body simulations  1 . Simulations provide only limited physical insight into the complex non-linear phenomena involved in gravitational clustering, and are intractable for performing multiple evaluations in parameter space.
Therefore, the main purpose of this paper is to develop an analytical model for predicting the 3PCF applicable to both the linear and nonlinear regimes. For this purpose, we need a model to correctly describe the redshift evolution and statistical properties of gravitational clustering up to the three-point level. We employ the so-called dark matter halo model, where gravitational clustering is described in terms of correlations between and within dark matter halos. Originally, this model was developed to express nonlinear clustering as the real-space convolution of halo density profiles (Neymann & Scott 1952;Peebles 1974;McClelland & Silk 1977; and also see Scherrer & Bertschinger 1991;Sheth & Jain 1997;Yano & Gouda 1999;Ma & Fry 2000b,c). Most recent works have relied on the the Fourier-space formulation, since the forms of the power spectrum and the bispectrum become much simpler (Seljak 2000;Ma & Fry 2000c;Peacock & Smith 2000;Scoccimarro et al. 2001;Cooray & Hu 2001a,b;Berlind & Weinberg 2002;Hamana, Yoshida & Suto 2002;Takada & Jain 2002 hereafter TJ02;Scranton 2002; and also Cooray & Sheth 2002 for a recent review). The halo model appears remarkably successful in that, even though it relies on rather simplified assumptions, it has reproduced results from numerical simulations (Seljak 2000;Ma & Fry 2000c;Scoccimarro et al. 2001;TJ02) and also allowed for interpretations of observational results of galaxy clustering (Seljak 2000;Scoccimarro et al. 2001).
We formulate the 3PCF model so that it can be applied to general 3D and 2D cosmological fields, such as the mass and galaxy distributions and the cosmic shear fields. Our method is built on the real-space formulation for the correlations of three particles in one halo. This is because the real-space approach enables us to compute the one-halo contribution to the 3PCF by a 4-dimensional integration, which is an advantage compared to the Fourier space approach. For the 2-and 3-halo terms, we rely on the Fourier-space approach and the approximations presented in Scoccimarro et al. (2001;also see TJ02). We study the transition from the quasi-linear to nonlinear regimes and the relative contribution of the different terms to the 3PCF. We show the halo model predictions for the 3PCF of the mass and galaxy distribution and the weak lensing convergence field for the currently favored CDM model. To do this, we will focus on the dependences of the 3PCF on the triangle configurations as well as on the properties of the halo profile, its inner slope and concentration.
This paper is organized as follows. In §2 we present the formalism to compute the 3PCF of the 3D density field based on the real-space halo model and review the Fourier-space halo approach and halo model ingredients (halo density profile, mass function and bias model) we will use in this paper. In §3, we present the method to compute the angular 3PCF. In §4 we show the results of the halo model predictions for the 3PCF of the mass and galaxy distributions and the weak lensing convergence field. Finally, §5 is devoted to a summary and discussion. We give some useful approximations for calculating the 3PCF in Appendix A and B. We use the currently favored CDM model (ΛCDM) with Ωm0 = 0.3, Ω λ0 = 0.7, h = 0.7 and σ8 = 0.9. Here Ωm0 and Ω λ0 are the present-day density parameters of matter and cosmological constant, h is the Hubble parameter, and σ8 is the rms mass fluctuation in a sphere of 8h −1 Mpc radius. The choice of σ8 for this model is motivated by the cluster abundance analysis (Eke, Cole & Frenk 1996).

Real-space halo approach
In this section, we briefly review the real-space dark matter halo approach (McClelland & Silk 1977;Scherrer & Bertschinger 1991;Sheth & Jain 1997;Yano & Gouda 1999;Ma & Fry 2000b,c), where the n-point correlation functions are described by the real-space convolution of halo density profiles.
In the halo approach, we assume that all the matter is in halos with the density profile ρ h (x; m) parameterized by a mass m. The halo mass is given by where Vvir is the virial volume of the halo. It is convenient to introduce the normalized halo profile defined as um(x; m) = ρ h (x; m)/m, satisfying the condition V vir d 3 x um(x; m) = 1. (2) In this paper we employ the virial volume as the boundary of a given halo, which can be formally defined, for example, based on the spherical top-hat collapse model. However, since in reality some matter is distributed outside the virial region, it is non-trivial to choose which boundary condition to use for the halo model. We will discuss how this uncertainty affects the halo model predictions, and show that it does not matter in the non-linear regime.
The density field at some point can be given by the superposition of each halo density profile: ρ(x) = i mium i (x; mi). The mean density of the universe is thus obtained from the ensemble average (Scherrer & Bertschinger 1991) where we have replaced the ensemble average by a spatial average and an average over the halo mass function n(m): i δD(m − mi)δ 3 D (x ′ − xi) = n(m). Note that from equation (2) the integral over x ′ is unity, so the above equation is a statement about the normalization of the halo mass function. Throughout this paper we work in the comoving coordinate.
Within the framework of the halo approach, the two-point correlation function (2PCF) of the density field can be expressed as the sum of correlations within a single halo (1-halo term) and between different halos (2-halo term): where ξ(x1 − x2; m1, m2) is the 2PCF of two halos with mass m1 and m2, and we have used r = x − x ′ . Here we have again employed the virial volume for the integration range. In the other words, the integration cutoff reflects the fact that we only account for mass contributions out to the virial region. If we require that ξ(r) is a function of the separation r alone (from statistical homogeneity and isotropy), then the density profile um which we need here is an average over all halos of a given mass so that the resulting 2PCF has no particular direction dependence of r (see similar discussion in Seljak 2000). This holds for a spherically symmetric density profile, which is assumed in this paper for simplicity. If one considers non-spherically symmetric profiles, it requires some extra averages over all possible shapes of halos of a given mass. Equation (4) shows that we can obtain ξ 1h (r) by a 3-dimensional integration for a spherically symmetric density profile, while the 2-halo term needs an 8-dimensional integration at most. In the same spirit as the derivation of equation (4), we can derive expressions for the 1-halo contribution to the n-point correlation function (n ≥ 2): Interestingly, this equation means that we can obtain any n-point correlation function by a 4-dimensional integration, once the density profile and the mass function are given. The 1-halo contribution is expected to be dominant in the strongly non-linear regime (δ ≫ 1), and so the real-space halo approach can be useful for calculations in this regime. In this paper, we focus on the halo model predictions of the 3PCF, which is the lowest order statistic for describing non-Gaussianity of the density field. The 1-halo term can be expressed as Under the assumption of statistical symmetry, the 3PCF can be described by three independent parameters characterizing triangle configurations. Likewise, the 2-and 3-halo terms in 3PCF can be derived within the real-space halo approach (Ma & Fry 2000c); for example, the 2-halo term becomes This equation implies that we have to perform at most an 8-dimensional integration to get ζ 2h for a given halo-halo correlation function, ξ(r; m, m ′ ), while the 3-halo contribution to ζ similarly requires the 12-dimensional integration.

Fourier-space halo approach
Most recent studies of the halo approach has relied on the method build in the Fourier space (Seljak 2000;Ma & Fry 2000c;Peacock & Smith 2000;Scoccimarro et al. 2001;Cooray & Hu 2001a,b;TJ02), since the Fourier-transformed counterparts of the n-point correlation function (power spectrum, bispectrum and etc.) can be expressed simply.
Following the notation introduced in Cooray & Hu (2001a), the power spectrum can be expressed as the sum of the 1-halo and 2-halo contributions with Hereũm(k) is the Fourier transform of the density profile and the following definition has been used recently (e.g., Seljak 2000;Scoccimarro et al. 2001;Cooray & Sheth 2002): where j0 is the zeroth-order spherical Bessel function, j0(x) = sin x/x, and rvir denotes the virial radius. In the following, quantities with tilde symbols denote the corresponding Fourier transformed quantities. Notice that all relevant quantities also depend on redshift, although we often omit it in the argument for simplicity. In contrast to the truncated halo profile above, Ma & Fry (2000c) employed the non-truncated profile whose Fourier transform is defined over an infinity integration range asũm(k) = ∞ 0 4πr 2 drum(r)j0(kr). The different halo boundary leads to different results in the Fourier-space halo model. In particular, we will focus on the difference between the predictions for the mass bispectrum in Ma & Fry (2000b), Scoccimarro et al. (2001) and Cooray & Sheth (2002) in comparison with the real-space halo model for the 3PCF in this paper. In the derivation of equation (8), we have assumed that the power spectrum for the halo-halo 2PCF can be given by P (k; m1, m2) = b(m1)b(m2)P L (k), where P L (k) is the linear power spectrum and b(m) denotes the bias parameter between the distribution of halos and the underlying density field. Through the definition ξ(r) = ∞ 0 k 2 dk/(2π 2 )P (k)j0(kr), we can compute the 2PCF of the density field from the power spectrum (8) within the framework of the Fourier-space halo model. The two derivations of ξ(r) based on the real-space and Fourier-space halo approaches should be equivalent to each other. However, the relation between ξ(r) and P (k) uses a Fourier transform over the infinite volume, while the integrations used in the halo model (see equations (4) and (10)) are confined to the virial volume of halos. This could lead to some discrepancies between the two halo model predictions for ξ(r), which will be carefully investigated below.
Here B P T denotes the perturbation theory bispectrum (e.g., see Jain & Bertschinger 1994) given by Here we have neglected the weak dependences of B P T on the cosmological parameter Ωm0. Using the bispectrum, the 3PCF can be defined as where k123 = k1 + k2 + k3.

Computational ease in real-space versus Fourier-space
Before specifying the details of our method, we briefly compare the real-space and Fourier-space halo approaches for predicting the 3PCF of the mass distribution. Let us first consider the 1-halo contribution to the 3PCF. Equation (6) shows that, once the halo profile and mass function are given, the real-space halo model enables us to obtain the 1-halo term by a 4-dimensional integration for arbitrary triangle shapes. On the other hand, equations (11) and (13) imply that a 7-dimensional integration at most is required to get the 1-halo term based on the Fourier-space halo model, which is intractable with current computational expenses. Therefore, the real-space halo approach has a great advantage in predicting the 1-halo term of the 3PCF.
However, the Fourier-space approach becomes useful for predicting the 2-halo and 3-halo terms, because the real-space model requires at most 8-and 12-dimensional integrations for calculating the 2-and 3-halo terms for given 2PCF and 3PCF of the halo distribution. On the other hand, from equations (11) and (13), one can see that the Fourier-space model allows us to obtain the 2-and 3-halo terms basically by a 7-dimensional integration. However, direct integration is still intractable. Therefore, in this paper we employ the Fourier-space model and develop approximations for calculating the 2-and 3-halo terms. Those approximations can significantly reduce the computational expenses and make the calculations tractable by regular numerical integrations, but still have adequate accuracy for our purpose.

Ingredients of the halo model
To complete the halo model approach, we need suitable models for the three ingredients: the halo density profile, the mass function of halos and the biasing of the halo distribution, each of which depends on halo mass m and redshift z and have been well studied in the literature.
We consider density profiles of the form where ρs is the central density parameter and c is the concentration parameter. In most parts of this paper, we will take the NFW profile with α = 1 (Navarro, Frenk & White 1997). However, since other simulations with higher spatial resolution have indicated α ≈ −1.5 (Fukushige & Makino 1997;Moore et al. 1998;Jing & Suto 2000), we will consider the effects of the variation in α on the 2-and 3-point correlation functions. The parameter ρs can be eliminated by the definition (1) of halo mass as where f = 1/[ln(1+c)−c/(1+c)] and 2F1 denotes the hypergeometric function. We employ the virial radius given by the spherical top-hat collapse model: m = (4πr 3 vir /3)ρ0∆(z), where ∆(z) is the overdensity of collapse given as a function of redshift (e.g., see Nakamura & Suto 1997 and Henry 2000 for a useful fitting formula). We have ∆ ≈ 340 for the ΛCDM model.
To give the halo profile in terms of m and z, we further need to express the concentration parameter c in terms of m and z; however, this still remains somewhat uncertain. We consider the following c parameterized by the normalization and the mass slope: where m * (z = 0) is the nonlinear mass scale defined as δc(z = 0)/D(z = 0)σ(m * ) = 1. Here δc is the threshold overdensity for spherical collapse model (see Nakamura & Suto 1997 and Henry 2000 for a useful fitting function), σ(m) is the present-day rms fluctuations in the matter density top-hat smoothed over a scale Rm ≡ (3m/4πρ0) 1/3 , and D(z) is the growth factor (e.g., see Peebles 1980). The values c0 ∼ 10 and β ∼ O(10 −1 ) are theoretically expected. The redshift dependence (1 + z) −1 is taken based on the numerical simulation results in Bullock et al. (2001). We take c0 = 10 and β = 0.2 as our reference model for the NFW profile, because TJ02 pointed out that the model can reproduce simulation results for the variance, skewness of weak lensing fields and the shear kurtosis. We will consider the influence of variations in the concentration parameter on the halo model predictions of the 2PCF and 3PCF. We will often use the following useful form for the Fourier transform of the NFW profile , which is needed for the Fourier-space halo model calculations: where η ≡ krvir/c, Ci(x) = − ∞ x dt cos t/t is the cosine integral function and Si(x) = x 0 dt sin(t)/t the sine integral. The profile (14) with general α has the asymptotic behaviors ofũm(k) ≈ 1 and ∝ k −3+α for k ≪ c/rvir and k ≫ c/rvir, respectively. It is worth stressing that the equation above is derived from the NFW profile truncated at the virial radius to maintain mass conservation (see equations (1) and (3)). The truncated profile leads to the consequence that the 1-halo term of the power spectrum behaves like shot-noise as P 1h (k) ∝ k 0 at small k. If one employs the non-truncated NFW profile, the resultingũm(k) is logarithmically divergent at k → 0 (see equation (22) in Ma & Fry 2000b) originating from the mass divergence from ∞ 0 4πr 2 drρ h (r), which does not yield the shot-noise power spectrum. For the halo mass function, we adopt an analytical fitting model proposed by Sheth & Tormen (1999), which is more accurate on cluster mass scales than the original Press-Schechter model (Press & Schechter 1974). The number density of halos with mass in the range between m and m + dm is given by where ν is the peak height defined by and the numerical coefficients a and p are empirically fitted from N -body simulations as a = 0.707 and p = 0.3. The coefficient A is set by the normalization condition ∞ 0 dνf (ν) = 1, leading to A ≈ 0.129. Note that the peak hight ν is specified as a function of m at any redshift once the cosmological model is fixed.
Mo & White (1996) developed a useful formula to describe the bias relation between the halo distribution and the underlying density field. This idea has been improved by several authors using N -body numerical simulations (Mo, Jing & White 1997;Jing 1998;Sheth & Lemson 1999;Sheth & Tormen 1999); we will use the fitting formula of Sheth & Tormen (1999) for consistency with the mass function (18) where we have assumed scale-independent bias and neglected the higher order bias functions (b2, b3, · · ·). This bias model is used for calculations of the 2-halo term in the 2PCF and the 2-and 3-halo terms in the 3PCF (see equations (8) and (11)). It should be noted that the requirement that the 2-halo term in the power spectrum gives the linear power spectrum in the limit k → 0 (ũm(k) → 1) imposes the condition dνf (ν)b(ν) = 1 for the integration mass range, since we have assumed that all the matter is in the form of virialized halos.

ANGULAR N -POINT CORRELATION FUNCTION
In this section, we present the angular n-point correlation functions in analogy with the real-space halo approach to the 3D n-point correlation functions in §2.1. The results are applicable for a general projected cosmological field such as the angular galaxy distribution and weak lensing fields.

Projected density field
We first consider the 2-dimensional projection Σ(θ) of the density field along the line of sight: where W (χ) is the weight function, and χ and dA(χ) are the comoving distance and comoving angular diameter distance, respectively. Note that χ is related to redshift z via the relation dχ = dz/H(z). In this paper we take the projected field to be the weak lensing convergence field κ (e.g., see Bartelmann & Schneider 2001) with the weight function given by where we have considered a single source redshift zs, corresponding to the comoving distance χs, for simplicity. The expressions derived below can be applied to the angular clustering of galaxies as well. Setting W (χ) = dN/dχ, the selection function of galaxies, and replacing the density field ρ by ρ galaxy (see §5 below for explicit expressions in terms of the halo occupation number of galaxies) will give the projected density of galaxies Σ galaxy from equation (21), and can be used to replace Σm in the equations below.

The 1-halo contribution to angular n-point correlation functions
Analogous to the 3D real-space halo approach, we can express the 1-halo contribution to the angular 2-point correlation function of the convergence field as where d 2 V /dχdΩ is the comoving differential volume given by d 2 V /dχdΩ = d 2 A (χ) for a flat universe and we have assumed a circularly symmetric profile Σm for halo of a given mass. We have also employed the flat-sky approximation, the Limber approximation (Kaiser 1992), which means that we account for contributions to w 1h from lens structures at the same redshift, and the Born approximation. Finally, as in equation (4), we have restricted the integration range to the circular area Ωvir enclosed by the projected virial radius. Equation (23) implies that we can obtain w 1h by a 4-dimensional integration if we have an analytical expression for Σm. In the halo model, Σm is the normalized column density field for a halo with mass m defined as The factor d 2 V /dχdΩ = d 2 A (χ) in equation (23) arises from the substitution d 2 x ⊥ = d 2 A d 2 ϕ, if one wishes to relate equation (23) to equation (4) for the 3-dimensional case. It is worth noting that, in analogy with equation (2), Σm should be defined so that it satisfies the mass conservation condition inside the virial radius: For the NFW profile (14), Σm(θ) can be analytically derived using formulae (2.266) and (2.269.2) in Gradshteyn & Ryzhik (2000) as with where θvir is the angular size of the virial radius of a halo at a given distance χ: θvir = rvir/dA(χ). Note that our expression for Σm differs from the result of Bartelmann (1996) because he took the projection over r = [−∞, ∞] (his expression includes typos, so, e.g., see Wright & Brainerd 2000 for the correct expression). The difference turns out to be pronounced, e.g. using the expression for c → ∞ in the above equation leads to discrepancies between the real-and Fourier-space halo model calculations of w 1h (θ), as shown below. It is useful to replace the projected halo profile Σm in equation (23) by the convergence field, κm(θ), for a halo of mass m: where we have used (3/2)Ωm0H 2 0 /ρ0 = 4πG. Hence, equation (23) can be rewritten in a more physically transparent form This equation means that the 2PCF of the weak lensing convergence can be expressed as the line-of-sight integration of the lensing contribution of halos at different redshift, weighted with the halo number density. An interesting application of this formulation is the study of the 2PCF of the reduced shear field, g = γ/(1 − κ) (γ is the shear due to lensing) which is a direct observable of the cosmic shear measurement. Equation (29) no longer relies on the power spectrum to compute the 2PCF. Therefore, replacing κm in the above equation with g m for a given halo will yield predictions of the reduced shear 2PCF. So far cosmological interpretations of cosmic shear data have been made by comparing the data with the theoretical model of the shear 2PCF computed from a model of the 3D power spectrum (e.g., see Van Waerbeke et al. 2001). The non-linear correction in using the reduced shear rather than γ could be important on sub-arcminute scales, where κ is of O(10 −1 ) for massive halos. This study will be presented elsewhere.
Performing the 2-dimensional Fourier transform of equation (29) yields the 1-halo contribution to the angular power spectrum of the convergence field, Cκ(l): where we have assumed that the Fourier transform confined to the virial region is well approximated as the 2D Fourier transform over an infinite integration range andκm(l) is the 2D Fourier transform of κm defined as Here J0(x) is the zeroth-order Bessel function. The expression (30) coincides with the form in Cooray et al. (2000) (see also the original formulation in Cole & Kaiser 1988).
To derive Cκ(l), the usual approach is to model the 3D power spectrum and then employ Limber's equation (e.g., see TJ02 and see references therein) to get Using the expression (8) of P (k) based on the Fourier-space halo model leads to the another expression of w 1h (θ): where k = l/dA(χ). The difference between the real-space and Fourier-space approaches in equations (23) and (33) is only in the order of the projection. Therefore, the two approaches should be equivalent to each other (see also discussions in Cooray et al. 2000;Cooray & Hu 2001a), if the weight function is a smooth function with respect to redshift as required for Limber's approximation. We will check this by comparing the predictions from the two methods In a similar manner, we can derive the 1-halo term for the angular n-point correlation function based on the real-space halo approach as This equation implies that we can obtain the 1-halo term of any n-point angular correlation functions by a 4-dimensional integration. This expression will be useful in investigating statistical properties of the projected field at non-linear angular scales. In this paper, we focus on the 3-point correlation function of the weak lensing convergence, given by From statistical symmetry, Z can be given as a function of three independent parameters characterizing the triangle shape on the sky. On the other hand, in the conventional approach to the 3PCF, Z is expressed in terms of the 3-dimensional bispectrum as Figure 1. The 2-point correlation function of the mass density field, ξ(r), at z = 0 for ΛCDM. The solid curve is the prediction from the Peacock & Dodds fitting formula. The thick and thin dashed curves show the 1-halo contribution to ξ(r), calculated using the real-space halo approach with different boundary conditions. The thick and thin dotted curves are the 1-halo and 2-halo terms from the Fourier-space halo approach. The dot-dashed curve is the total prediction of the 1-halo plus 2-halo terms, where we have used the real-space (thick dashed curve) and Fourier-space (thin dotted curve) models for the 1-halo and 2-halo predictions, respectively.
where ki = li/dA(χ). We use this formulation for evaluating the 2-halo and 3-halo contributions to Z, using the expressions (11) for the 2and 3-halo bispectrum. The resulting 6-dimensional integrations are intractable, so we develop approximations, presented in Appendix A, to evaluate them.

RESULTS
In this section, we present halo model predictions for the three-point correlation function for the mass and galaxy distribution and for the weak lensing field. We employ a scale invariant spectrum of primordial fluctuations with the BBKS transfer function (Bardeen et al. 1986).

2-point correlation function
Before considering the 3PCF, we demonstrate the validity of the real-space halo model by comparing our results for the two-point correlation function with those of the Fourier-space halo model well studied in the literature. For the NFW profile, the 1-halo term in equation (4) can be further analytically simplified as with where we have taken the virial volume for the integration range of d 3 s in equation (4) and the boundary condition um(x) = 0 for x > rvir.
The predictions for ξ 1h on nonlinear scales are not sensitive to the boundary condition, as shown below. Note that, although the integration over ds can be further analytically calculated as done in Ma & Fry (2000c) and in Sheth et al. (2001), the resulting expression for ξ 1h is lengthy for the boundary condition and so we stopped at the above expression. Figure 1 plots the mass two-point correlation function ξ(r) at z = 0 for the ΛCDM model. The halo model prediction is shown by the dot-dashed curve, which matches the solid curve given by the fitting formula of Peacock & Dodds (1996; hereafter PD). We consider  (14)). (c) The right panel shows the dependence on the concentration parameter. The result for our fiducial model (c 0 , β) = (10, 0.2) is shown by the thick solid curve. the NFW profile with the concentration specified by (c0, β) = (10, 0.2) in equation (16) as the fiducial model. During the preparation of this paper, Smith et al. (2002) have proposed a new fitting formula of the non-linear power spectrum that can better match high-resolution N -body simulations. For the ΛCDM model, the effect is small: the power spectrum differs from the PD prediction by ∼ 10% at most at k < ∼ 30h/Mpc. The thick dashed curve shows the real-space halo model prediction for the 1-halo contribution to ξ computed from equation (37), while the thick dotted curve is the prediction from the Fourier-space halo model. The two halo model predictions are in remarkable agreement for small scales r < ∼ 2 h −1 Mpc, but deviate slightly at the transition scale of ∼ 2 h −1 Mpc between the non-linear and linear regimes. The discrepancy between the real-space and Fourier-space halo models could be ascribed to the boundary condition used in the integration. The thin dashed curve shows the result obtained without using the boundary condition (um(x) = 0 for x > rvir) in equation (37), implying that the change alters the predictions only on scales r > ∼ 1 h −1 Mpc. However on these scales the 2-halo term dominates as shown by the thin dotted curve. Hence, to summarize the results in this figure, the real-space halo model can be used to predict the statistical properties of the non-linear density field with the same accuracy as the Fourier-space halo model whose validity has been carefully investigated in the literature (Seljak 2000;Ma & Fry 2000c;Scoccimarro et al. 2001;TJ02).
The left panel of Figure 2 shows the dependence on the maximum mass cutoff used in the calculation on ξ(r). From top to bottom, the solid curves show the results for the maximum cutoff of 10 16 , 10 15 , 10 14 , 10 13 and 10 12 M⊙. Massive halos with m > ∼ 10 14 M⊙ yield the dominant contribution (greater than ∼ 80%) to ξ 1h (r) for r > ∼ 0.1 h −1 Mpc, while less massive halos become relevant at smaller scales. The dashed curve denotes the full halo model prediction for ξ(r) (1h+2h) as shown in Figure 1.
The middle and right panels of Figure 2 show the dependences of ξ(r) on the variations in the halo profile. From the middle panel, one can see that increasing α for the halo profile of equation (14) leads to a higher amplitude and steeper slope for ξ(r) on small scales (see also Jain & Sheth 1997;Ma & Fry 2000c). For generic α, the angular integration in d 3 s in equation (4) is not analytic, so we performed the 3-dimensional integration to get ξ 1h (r) numerically. The right panel shows the dependences on the concentration parameter parameterized in terms of the normalization c0 and the slope β as c(m) = c0(m/m * ) −β . The result of our fiducial model (c0, β) = (10, 0.2) is shown by the thick solid curve. It is apparent that increasing c0 with fixed β or decreasing β with fixed c0 increases ξ at small scales. This is because these variations lead to more concentrated density profiles for halos more massive than the nonlinear mass scale m * . Since the massive halos dominate the contribution to ξ at the scales considered here, this has the effect of increasing ξ. The results in Figure 2 imply that we can adjust the inner slope α and the concentration parameter, both of which are somewhat uncertain theoretically and observationally, so that the halo model can reproduce the PD result for a given mass function (see Seljak 2000 for a similar discussion). We will argue that the 3PCF can be used in combination with the 2PCF to constrain properties of the halo profile. Figure 3 shows the 2PCF of the weak lensing convergence field. The real-space halo model prediction (dashed curve) for the 1-halo contribution is compared with the result of the Fourier-space halo model (thick dotted curve), computed from equations (29) and (33), respectively. Here we have considered the NFW profile and source redshift zs = 1 for simplicity. The different halo model predictions perfectly agree over the angular scales we have considered. This agreement is encouraging, because the real-space halo approach can then be used to predict the 3PCF with the same accuracy expected as for the Fourier-space model well studied in the literature (TJ02 and references therein). For comparison, the solid curve denotes the PD result and the thin dotted curve is the Fourier-space halo model prediction of the 2-halo term. It clarifies that the 2-halo term is important on θ > ∼ 3 ′ for the ΛCDM model. The dot-dashed curve is the total contribution to the 2PCF from the 1-and 2-halo terms, which slightly deviates from the PD result on the small scales θ < ∼ 1 ′ . As discussed below equation (26), if we use the expression in Bartelmann (1996) for the NFW convergence field in the calculation of the 1-halo term, one finds discrepancies  (1h+2h), if the NFW convergence field in our equation (29) is replaced by Bartelmann's (1996)  between the real-space and Fourier-space halo models. The broken curve shows the result (1-halo + 2-halo terms), which overestimates wκ by 10-25 % over the scales considered.

Dark matter correlation function
We first consider the 3PCF of the mass density field. To obtain the 1-halo contribution, where all the three matter particles reside in one halo, we perform the following 4-dimensional numerical integration; where |s + r| = (s 2 + r 2 + 2sr cos θ) 1/2 and |s + q| = [s 2 + q 2 + 2sq(sin ψ sin θ cos ϕ + cos ψ cos θ)] 1/2 . Note that we employ the boundary condition um(r) = 0 for r > rvir, but this does not affect the final results on the non-linear scales, as explained in Figure 1 and explicitly demonstrated in Figure 8. In the above equation, we have used the fact that the 3PCF can be expressed as a function of the three independent parameters r, q and ψ specifying the triangle configuration (see Figure 4).
To complete our halo model predictions, we develop approximations for calculating the 2-and 3-halo terms in Appendix A. The 2-halo term is relevant only for small range of scales in the transition between the quasi-linear and non-linear regimes (r ∼ 1 h −1 Mpc), while the 3-halo term is dominant on larger scales. As in the literature, we consider the hierarchical 3PCF amplitude defined as The results for the NFW profile at z = 0 are plotted for the ΛCDM model. The triangle shape is parameterized by three parameters r, q and ψ as illustrated in Figure 4. The four panels show the results for Q vs. ψ, with r = 0.02, 0.1, 1.0 and 10 h −1 Mpc as indicated in the panels. In each panel, the solid, dashed, dot-dashed, dotted, and broken curves denote the Q parameter for r/q = 1, 2, 3, 4 and 5, respectively.
where ζ = ζ 1h + ζ 2h + ζ 3h . In the following halo model predictions for Q, we will maintain self-consistency by using the halo model for calculating the 2PCF in the denominator of Q. Figure 5 plots the dependences of Q(r, q, ψ) on the size and shape of triangle configurations. In this figure we have used the NFW profile. The four panels show the results for r = 0.02, 0.1, 1 and 10 h −1 Mpc. In each panel, the solid, dashed, dot-dashed, dotted and broken curves plot the Q values for q/r = 1, 2, 3, 4 and 5 as a function of the interior angle ψ. It is apparent that on the highly non-linear scales r < ∼ 0.1 h −1 Mpc the 3PCF has a weak dependence on the triangle configuration; even so, for plausible halo model parameters the hierarchical form Q ≈ const. does not hold in that Q depends both on the size and shape of the triangle configuration (also see discussion in Ma & Fry 2000b). The curves with q/r = 1 show that Q decreases with decreasing ψ at ψ < ∼ 0.2π because of the strong suppression due to ξ 2 in the denominator of Q. On larger scales r > ∼ 1 Mpc, an oscillatory feature in Q appears as predicted by perturbation theory (e.g., Fry 1984b; Barriga & Gaztañaga 2002). Thus, the configuration dependence of the 3PCF is more prominent on quasi-linear scales than on strongly non-linear scales.
The features observed in Q can be understood using Figure 6, which shows the 1-, 2-and 3-halo contributions to Q separately for the triangle configurations in Figure 5. Note that these separate contributions are shown for the 3PCF, the numerator of Q, while the 2PCF factors in the denominator include the full contribution from the 1-plus 2-halo terms. This figure clarifies that the 1-halo term yields the dominant contribution to the 3PCF on small scales r < ∼ 1 Mpc. The 2-halo term becomes relevant over the transition scales 1 < ∼ r < ∼ 5 Mpc between the non-linear and linear regimes, while the 3-halo term dominates for r > ∼ 10 Mpc. Figure 7 displays the halo model predictions for Q for equilateral triangles as a function of the side length r. The three dashed curves plot the 1-, 2-and 3-halo contributions separately. For comparison, the perturbation theory result is shown by the dotted curve. The figure again clarifies that, at the small scales r < ∼ 0.5 h −1 Mpc, the 1-halo term dominates the total contribution, while the 3-halo term captures the large-scale correlations in the quasi-linear regime. A comparison between the perturbation theory result and the 3-halo term reveals that our approximation (A4) reproduces the perturbation theory result at large scales r > ∼ 5 h −1 Mpc. Thus the way the halo model provides a separate description of the 1-, 2-and 3-halo terms can clarify how gravitational clustering transits from the linear regime to the strongly non-linear regime as one goes to smaller scales.
The halo model predictions match N -body simulation results on small scales r < ∼ 1 h −1 Mpc (R. Scoccimarro; private communication). However, the Q parameter from the halo model has a bump feature at r ≃ 2 h −1 Mpc which does not seem to be present in simulation data. The corresponding bump feature in the bispectrum has been found at k ∼ 1 hMpc −1 in previous work (see Figure 2 in Scoccimarro et al. 2001 and Figure 13a in Cooray & Sheth 2002). Hence, it is unlikely to be due to inaccuracies in our approximations used for the 2-and 3-halo term calculations, although the approximations tend to overestimate the true amplitudes to some extent as carefully investigated in TJ02. It probably does not reflect real properties of dark matter clustering either (as discussed below for Figure 10). Rather, it can be ascribed to an inaccuracy in the standard implementation of the halo model which appears only on transition scales of order a few Mpc, in contrast to the success of the halo model in the non-linear and linear regimes.
One possible origin of the inaccuracy is the sharp cutoff at the virial radius for the integrals for the 2PCF and 3PCF, as have been used for the real-space halo model in this paper as well as for the Fourier-space model in Scoccimarro et al (2001) and Cooray & Sheth (2002). As shown in Figure 1, modifications of the integration range could alter the halo model prediction at the transition scales. In the upper panel of Figure 8, the thick dashed curve shows the result obtained without using the boundary condition um(x) = 0 for x > rvir for the calculations of the 1-halo terms of the 2PCF and 3PCF, which enter the denominator and numerator of Q, respectively. One can see that the bump feature is weakened through complex dependences -the 2PCF amplitude is enhanced by the modified boundary condition more strongly than the 3PCF. In fact, this trend is verified by the results in Figure 3 and 4 in Ma & Fry (2000b), which show that the Fourier-space halo model yields bispectra with no bump feature. Ma & Fry employed the Fourier transform of the non-truncated profile,ũm(k), calculated over an infinite integration range, which includes mass contribution outside the virial region. However, in this paper we have used the virial boundary condition since it preserves mass conservation. We have also clarified in Figure 3 that, if we include mass contributions beyond the virial radius, the halo model prediction for the lensing 2PCF overestimates the amplitude (see the dot-dashed and broken curves in Figure 3). In reality, the mass distribution outside the virial region eventually merges into quasi-linear structures, such as filamentary structures in the universe, which are relevant for the transition scale (δ > ∼ 1), and are unlikely to follow the halo profile far outside the virial region. Therefore, a careful investigation of which boundary conditions to use for the halo model integrals will be needed to achieve more accurate predictions at the transition scale.
Another effect that could be manifested at these scales is the exclusion effect of different halos for the 2-and 3-halo terms. We could impose the condition that different halos are separated by scales larger than the sum of their virial radii. The halo model does not explicitly account for this effect. Since the transition scale ∼ 1 Mpc is of the order of the virial radius of massive halos that make a significant contribution to the integrals, the halo model used here might overestimate the contribution from the 2-and 3-halo terms at these scales. The dot-dashed curve in the lower panel of Figure 8 shows an estimate of this effect. The estimate is obtained by imposing the condition rvir(m) ≤ r/2, an approximate prescription to exclude pairs and triplets of halos contributing to the 2-and 3-halo terms at scales smaller than the sum of their virial radii. The plot shows that including the exclusion effect suppresses the bump feature in Q. The dotted curve in   the lower panel of Figure 8 shows the result of applying both alterations (boundary and exclusion effects) to the calculations. A systematic resolution of these problems is beyond the scope of our paper, and will be considered elsewhere. Finally, the figure implies that, except for the problematic bump feature, a transition between the non-linear amplitude (Q > ∼ 3) and the quasi-linear amplitude (Q < ∼ 1) occurs at scales of a few h −1 Mpc, corresponding to a flattening of the 2PCF ξ(r) (see Figure 1). Figure 9 shows the dependences of Q on the cosmological model. We consider three models that differ from our fiducial model. One is the Λ CDM model normalized with σ8 = 1.1, motivated by the possible detection of the Sunyaev-Zel'dovich effect in the CMB (Bond et al. 2002). The second is a flat model without cosmological constant (SCDM), with parameters Ωm0 = 1, h = 0.5 and σ8 = 0.6. To clarify the dependence on the shape of the matter power spectrum, we also consider the τ CDM model, where the shape parameter Γ = 0.21, but the other parameters are same as in the SCDM model. The figure reveals that the Q parameter is not very different for the different models but does vary slightly with changes in the various model parameters. Thus accurate measurements must be interpreted with precise predictions that explore all parameter dependences.
For completeness we present our results with the triangles specified using a set of alternative parameters which have been used in the literature (e.g., Peebles 1980): with the condition r12 ≤ r23 ≤ r31 which imposes the constraints u ≥ 1 and 0 ≤ v ≤ 1. The 3PCF for any triangle can then be expressed as a function of three parameters r12, u and v, where u and v characterize the shape and r12 the size of a triangle. Figure 10 shows the dependence of Q on r12, u and v. The values of r12, u and v are chosen such that our model predictions can be compared with the N -body simulation results shown in Figure 2 in . The data with the error bars represent the simulation results, which were kindly made available to us by Yipeng Jing. Although the cosmological parameters in  are slightly different from our model, this fact does not matter because of the weak dependence of Q on the cosmological model as explained in Figure 9. It is also noted that the halo model prediction does not account for the effect of averaging over bin widths in u and r12 which is taken for the Q measurement from the simulations. This figure reveals that our halo model agrees well with the simulation results for the configuration dependence of Q and its amplitude at the non-linear scales r12 = 0.2 and 0.32 h −1 Mpc, as well as the quasi-linear scale r12 = 3.25 h −1 Mpc. However, as discussed above, for the case in which one side length of the triangle is comparable to 1 h −1 Mpc, the halo model tends to overestimate the simulation result. More detailed comparisons with recent higher resolution simulations will be presented elsewhere. Figure 11 shows how the 1-halo contribution to the mass 3PCF depends on the maximum mass cutoff used in the calculation against the side length of equilateral triangles. This can be compared to the dependence of the 2PCF on mass cutoff shown in the left panel of Figure  2. For the 3PCF more massive halos dominate the contribution at a given length scale: for example, at r = 0.5 h −1 Mpc, over half the contribution to the 3PCF is from halos with m > 10 15 M⊙, whereas for the 2PCF the mass range is m > 10 14 M⊙.
In Figure 12 and 13, we present dependences of the Q parameter on possible variations in the halo inner slope and the concentration parameter, as in the middle and right panels of Figure 2. To do this, we use the parameters r, q and ψ for the triangle configuration as in Figure 5. Note that the results shown are computed from only the 1-halo contributions to the 3PCF and 2PCF, which enter the numerator and denominator of Q, respectively, since the variations of the halo profile are relevant only at the non-linear scales. Increasing the inner slope parameter α for the profile (14) leads to higher Q (see similar discussions in Ma & Fry 2000a for the bispectrum). A comparison of Figures  12 and 13 shows that Q is more sensitive to α for triangles with smaller size. Compared with the effect on the 2PCF shown in Figure 2, this result also means that the 3PCF is more sensitive to the inner slope of halo profile than the 2PCF. We find the following fitting formula of the  (41)). The values of r 12 , u and v are the same as those used in Figure 2 in , where the Q parameters were computed from N -body simulations. The data with error bars denote the simulation results. Figure 11. The dependence of the mass 3PCF on the maximum mass cutoff used in the halo model calculations, plotted against the side length of the equilateral triangle. From top to bottom, the maximum mass is 10 16 , 10 15 , 10 14 , 10 13 and 10 12 M ⊙ . The dashed curve is the full 3PCF.
for r = 0.02 h −1 Mpc, while ξ 2 (r) = [ξ 2 NF W ] 1+0.048(α−1)+0.034(α 2 −1) , ζ(r) = ζNF W 1+0.06(α−1)+0.033(α 2 −1) , Q = QNF W 1+0.38(α−1)+0.057(α 2 −1) for r = 0.1 h −1 Mpc. Here ξNF W , ζNF W and QNF W denote the use of the NFW profile (α = 1); we have used equilateral triangle configurations for the 3PCF calculation. Figure 14 explicitly plots the α dependence of Q. The right panels in Figure 12 and 13 illustrate the sensitivity of Q to possible variations in the halo concentration parameter (16). One can see that the change in concentration affects the Q amplitudes, but does not strongly alter the configuration dependence. Increasing c0 with fixed β or decreasing β with fixed c0 leads to a higher amplitude for Q, which is the same trend as for the 2PCF shown in Figure 2. This implies that the 3PCF is more sensitive to the halo concentration than the 2PCF. These results in Figure 12 and 13 show that measuring the 3PCF could constrain halo profile properties complementary to the measurements of the 2PCF, in analogy with determinations of the galaxy bias parameter from measurements at quasi-linear scales (Feldmann et al. 2001;Verde et al. 2002).
In Figure 12-14, we have demonstrated the sensitivity of the mass 3PCF to the halo profile properties. We have compared our halo model predictions for Q with the asymptotic shape derived in Ma & Fry (2000b) based on the halo model. The asymptotic Q amplitude depends only on the mass slope β of the concentration parameter (see equation (16)), the slope α ′ in the low mass end of the mass dependence of the mass function, n(m) ∝ ν α ′ , and a primordial spectral index n. From equation (7) in Ma & Fry (2000b), the asymptotic Q for equilateral triangles is given by where β ′ ≈ 0.8β and we have considered equilateral triangle configuration. This equation was derived under three assumptions: the contribution from the exponential cutoff of the mass function in the high mass end was ignored, a scale-free primordial power spectrum was employed and the k-dependence of the Fourier-transformed halo profile,ũm(k) was ignored. The third assumption leads to the consequence that the asymptotic Q has no dependence on the inner slope of the halo profile. Since for the profile in equation (14)ũm(k) does depend on the inner slope asũm(k) ≈ 1 and ∝ k −3+α for k ≪ c/rvir and k ≫ c/rvir, the asymptotic Q should depend on the inner slope or more generally the halo profile shape. By comparing our predictions of Q with the analytical result of Ma & Fry (2000b), we find that our halo model predicts a steeper slope for Q(r) than expected from the asymptotic shape in equation (44), for plausible values of β and α ′ . This discrepancy is probably due to the assumptions used in analytically deriving the asymptotic shape. The r-dependence of Q arises from a complex superposition of contributions from the halo profile, the power spectrum shape and the shape of the mass function. Indeed,  as an example of this possibility, Taruya, Hamana & Kayo (2002) explicitly showed, based on the halo model, that the non-Gaussian tail of the probability distribution function of the density field arises from such a superposition for small smoothing scales.
An important implication pointed out in Ma & Fry (2000b,c) is that plausible halo model properties are unlikely to follow the stable clustering hypothesis. This hypothesis is useful because it allows us to analytically predict the behavior of non-linear gravitational clustering (e.g. Peebles 1980;Jain 1997). In fact, the popular PD fitting formula widely used in the literature was derived based on this hypothesis (also see Hamilton et al. 1991;Jain, Mo & White 1995). Very recently, Smith et al. (2002) showed that the non-linear power spectra measured from the high-resolution N -body simulations indeed showed a weak violation of the stable clustering hypothesis. Figure 9 in Smith et al. (2002) compares the k-slope of the measured non-linear power spectrum with the asymptotic prediction from Ma & Fry (2000b), which also shows some discrepancies. Hence, a careful investigation will be needed to clarify how non-linear gravitational clustering can be described by the halo model ingredients in connection with the stable clustering hypothesis and the hierarchical ansatz (see also the next subsection). Figure 15. The hierarchical amplitude Qn for the mass n-point correlation function in the strongly non-linear regime is shown versus the order n. We consider the configuration with n equal sides of length r = 0.1h −1 Mpc and equal interior angles. The solid curve shows the halo model prediction. The dashed curve shows the scaling for Qn as a function of Q 3 and n proposed by Fry (1984a), while the dotted curve shows the formula proposed in Hamilton (1988) (we have used the halo model result for Q 3 in these formulae). The dot-dashed curve is the prediction from hyper-extended perturbation theory (Scoccimarro & Frieman 1999).

The hierarchy of higher-order correlation functions
As stated in §2.1, the real-space halo model can be a useful tool for predicting the n-point correlation functions in the strongly non-linear regime, where correlations within one halo dominate the contribution to the n-point functions. The exact expressions for the the hierarchy are provided by the BBGKY equations (Peebles 1980), which govern collisionless gravitational system. However, their complex form makes it intractable to obtain solutions for the higher order correlations in the non-linear regime. Here we will use the halo model to shed some light on a fundamental question about gravity that has been discussed in the literature over the last decades: what is the asymptotic small-scale behavior of the hierarchy of the n-point functions under gravitational clustering (e.g., Peebles 1980;Fry 1984a,b;Hamilton 1988;Ma & Fry 2000b;Bernardeau et al. 2002a)?
The reduced n-point correlation function, Qn, is defined as where the denominator is given by all topologically distinct tree diagrams (the different n n−2 ways of drawing (n − 1) links that connect n points). This form is motivated by the expectation ξ (n) ∝ ξ n−1 , indicated observationally by the pioneering measurements of the two-and three-point functions from galaxy surveys (e.g., Groth & Peebles 1977;Peebles 1980), and theoretically by perturbation theory (Fry 1984b). The amplitudes Qn are a natural set of statistics to describe the non-Gaussianity that results from gravitational clustering. The solid curve in Figure 15 shows the halo model predictions for the hierarchical amplitudes Qn for 3 ≤ n ≤ 6. We consider configurations with n equal sides, each of length r = 0.1h −1 Mpc, and equal interior angles for evaluations of the n-point functions. We take the halo profile to have the NFW form. For comparison, the dashed curve shows the scaling law for Qn proposed in Fry (1984a): Qn = (4Q3/n) (n−2) n/(2n − 2). The dotted curve shows the scaling suggested by Hamilton (1988): Qn = (Q3/n) n−2 n!/2. The dotdashed curve is the prediction from hyper-extended perturbation theory for 3 ≤ n ≤ 5 (Scoccimarro & Frieman 1999). The halo model predicts stronger clustering with increasing n than the formulae of Fry and of Hamilton, but not as strong as hyper-extended perturbation theory. Note that higher-order correlations are more sensitive to variations in the inner slope and halo concentrations of massive halos.
The results in Figure 15 might be compared with Figure 3 in Szapudi et al. (1999), which show the cumulants Sn (3 ≤ n ≤ 10) measured from the N -body simulations. From equation (45) we can expect following rough relation between the cumulants and the reduced n-point correlation functions: [ξ (2) (r)] (n−1) ∼ n n−2 Qn(r).  Figure 15. Note that the simulation results were read from the dashed curves in Figure 3 in Szapudi et al. (1999) for the SCDM model (the cosmological parameters are slightly different from the model we have considered), and we expect similar values for the ΛCDM because of the weak dependence of Sn on cosmological models. One can see that the rapid increase in the amplitudes of Sn with n seen in the simulation result is better described by the halo model predictions than those of Fry and of Hamilton. However, as explicitly pointed out Table 1. Comparison of the theoretical predictions for the cumulants, Sn, at r = 0.1h −1 Mpc with the simulation result. The analytical predictions are estimated from the reduced n-point functions, Qn, shown in Figure 15 using equation (46). The simulation results correspond to the values denoted by the dashed curves in Figure 3 in Szapudi et al (1999), which are measured from N -body simulations for the SCDM model.

The 3-point correlation function of galaxies
It is straightforward to extend the dark matter halo approach to calculate the 3PCF of galaxies which can be measured from galaxy surveys.
To do this, we need to know the mean number of galaxies per halo of given mass, known as the halo occupation number, N gal (m), and the second and third moments of the galaxy distribution ( N gal (N gal − 1) and N gal (N gal − 1)(N gal − 2) ). Following the method in Scoccimarro et al. (2001), the real-space halo model expressions for the 1-halo contributions to the galaxy 2PCF and 3PCF are given by whereng denotes the mean number density of galaxies defined as ng = dm n(m) N gal (m).
In a similar manner, we can derive halo model predictions for the 2-halo contribution to the galaxy 2PCF, and the 2-and 3-halo contributions to the galaxy 3PCF. For simplicity we assume that the distribution of galaxies within a given halo follows the dark matter profile um(r). Note that in the large-scale limit, the galaxy bias parameter is defined as The halo occupation number has been used in the literature to explain galaxy clustering properties (Jing, Mo & Börner 1998;Seljak 2000;Peacock & Smith 2000;Scoccimarro et al. 2001;Sheth et al. 2001;Berlind & Weinberg 2002;Moustakas & Somerville 2001;Scranton 2002;Cooray 2002). In this paper, we employ the model of Scranton (2002), which approximately reproduces the results for N gal (m) from semi-analytic models in Kauffmann et al. (1999). The model gives N gal for red and blue galaxies 2 where the parameters γR = 1.1, mR = 1.8 × 10 13 h −1 M⊙ and mR0 = 4.9 × 10 12 h −1 M⊙ for red galaxies, while γB = 0.93, mB = 2.34 × 10 13 h −1 M⊙, A = 0.65, A0 = 6.6 and mB s = 11.73 for blue galaxies. In the following, we employ a lower mass cutoff of m ≥ 10 11 h −1 M⊙ for the calculations, since in small halos effects such as supernova winds can blow away the gas from halos and thus suppresses further star formation. The model above yields bias parameters b g,red = 1.51 and b g,blue = 0.82 for red and blue galaxies, respectively. Figure 16 plots the mean number of galaxies, N gal (m), as a function of parent halo mass. The figure shows that N gal ≤ 1 at m < ∼ 10 13 h −1 M⊙ for both types of galaxies; a significant fraction of low mass halos are thus likely to contain at most one galaxy. For such low mass halos, we need to carefully model the higher order moments of the galaxy distribution within them, as discussed below. It is also apparent that N gal for blue galaxies has a bump feature for halos of m < ∼ 10 12 h −1 M⊙. The bump is due to the prescriptions used in the semi-analytical model. Although stars are formed from the cold gas in halos within a dynamical time scale, it is assumed that the gas cooling in halos with circular velocity Vc > 350 km s −1 does not form visible stars so that the output galaxy catalog can fit the observed Tully-Fisher relation (see Kauffmann et al. 1999 for details). This sharp cutoff of the star formation prevents further formation of galaxies in halos with Vc > 350 km s −1 unless the halos experience merging, which is taken as the trigger to induce the starburst. Note that the circular velocity Figure 16. Mean number of red and blue galaxies as a function of parent halo mass using the forms in equation (50). The thin dashed curve shows the model for blue galaxies without the bump feature, which we will use to demonstrate the effect on the halo model prediction for the 3PCF in Figure 18. Vc = 350km s −1 roughly corresponds to halo of m ∼ 10 13 M⊙ at the present epoch; the velocity is larger for the same mass halo at earlier epochs. In fact, we will show that the bump feature drastically affects the configuration dependence and amplitude of the reduced 3PCF of blue galaxies.
As mentioned above, we also need the second and third moments of the galaxy distribution within parent halos, N gal (N gal − 1) and N gal (N gal − 1)(N gal − 2) . We follow Scranton (2002), who expressed them as N gal This model implies that galaxies follow a Poisson distribution for massive halos m > ∼ 10 13 h −1 M⊙ and a sub-Poisson distribution for smaller halos. Other treatments of the sub-Poisson distribution have been proposed by Sheth et al. (2001) and Scoccimarro et al. (20001).
We expect that the models of the 2nd and 3rd moments of N gal (m) will affect Q gal at small scales which are dominated by the 1halo term. Large scales > ∼ 5 h −1 Mpc, in the quasi-linear regime, are dominated by the 3-halo term. The behavior of Q gal on the large scales can be roughly expressed in terms of the bias parameter (49) as Q gal ∼ Qmass/bg. This follows from Q gal ∼ ζ gal /(3ξ 2 gal ) and that: ζ gal ∼ b 3 g ζmass and ξ gal ∼ b 2 g ξmass. As a result, the anti-biasing and bias of blue and red galaxies should lead to larger or smaller amplitudes of Q gal compared to Qmass, respectively. Figure 17 shows Q gal for the galaxy 3PCF with triangle configurations parameterized by r, q and ψ as in Figure 5. For red galaxies, the halo model combined with the halo occupation number leads to a weaker configuration dependence and smaller amplitude for Q gal than the mass 3PCF, as suggested from the measurements in . On the other hand, the 3PCF of blue galaxies displays complex features. This is mainly due to the bump feature in N gal (m) at 10 11 ≤ m < ∼ 10 13 h −1 M⊙, as shown below. Hence, a detailed knowledge of the configuration dependence of Q gal could be used to quantitatively constrain the occupation number, or more generally, the galaxy formation scenario as a function of observed galaxy properties (morphology, luminosity and so on) within a halo of given mass.
The sensitivity of the input model of the halo occupation number to the galaxy 3PCF can be more clearly understood from Figure 18, which shows the Q gal parameters for equilateral triangles. The 1-, 2-and 3-halo contributions to the galaxy 3PCF are explicitly shown. This figure may be compared with the result for the galaxy bispectrum shown in Figure 6 of Scoccimarro et al. (2001). It is again apparent that the amplitude of Q gal for red galaxies is suppressed compared to the mass 3PCF, while the result for blue galaxies exhibits complex features. Figure 19 shows halo model predictions for the 2PCF and 3PCF of mass and blue and red galaxies separately. It is worth noting that the 2PCF and 3PCF of blue galaxies are well approximated by power laws at scales r > ∼ 0.2 h −1 Mpc, unlike the case for the mass or for red galaxies. This is due to the fact that galaxy formation is inefficient in massive halos and thus suppresses the 1-halo term for blue galaxies compared to the mass or to red galaxies which are formed by mergers. Also, Q gal for blue galaxies has contributions from the 2-and 3-halo terms over a wider range of scales compared to that of red galaxies. This is mainly due to the bump feature in N gal at low mass scales in Figure 16, and partly due to the shallower mass slope in N gal ∝ m γ .
The lower panel of Figure 18 explores the origin of the complex features in Q gal for blue galaxies at 0.1 < ∼ r < ∼ 2 h −1 Mpc. The dotdashed curve shows how the halo model prediction for Q gal for blue galaxies changes if we suppress the bump feature in N gal at low mass scales by setting A = 0 in equation (50). The model N gal used is shown by the thin dashed curve in Figure 16. One can see that the complex features in Q gal disappear, except for the problematic bump feature at r ∼ 2Mpc (see the discussion around Figure 8), and its amplitude becomes smaller than the mass 3PCF because of the change of the bias parameter to bg = 1.24 from bg = 0.82, which results in behavior similar to the red galaxies. Thus, the input model of N gal has a drastic effect on Q gal . We also investigate how the halo model predictions are affected by altering the sub-Poisson distribution for the galaxy third-order moment, which is relevant for the range N gal ≤ 1 or equivalently m < ∼ 10 13 h −1 M⊙. Scoccimarro et al. (2001) proposed the binomial distribution, where the third-order moments is given by . Note that this model leads to negative values for the third moment for αg ≪ 1. The thin dotted curve in Figure 18 shows the result of using this prescription -it changes Q gal only for scales r < ∼ 0.2 h −1 Mpc. Figure 17. Halo model predictions for the galaxy 3PCF as a function of triangle shapes parameterized by r, q and ψ (see Figure 4). The solid, dashed and dotted curves show the results for the mass and for blue and red galaxies, respectively. We have employed the model of equation (50) for the the halo occupation number of red and blue galaxies.

The 3-point correlation function of weak lensing fields
Next we consider the halo model predictions for the angular 3PCF of the weak lensing convergence field. It should be noted that the following method can be easily extended to the predict the angular 3PCF of the galaxy distribution, as discussed in §3. The 1-halo term for the 3PCF of κ is given by the following 4-dimensional integral: where |θ + r| = (θ 2 + r 2 + 2θr cos ψ) 1/2 and |θ + q| = (θ 2 + q 2 + 2θq cos(ϕ − ψ)) 1/2 . For the NFW profile, κm is analytically given by equations (26) and (28). Note that r and q denote angular scales, although we use the same notations as for the 3D case for simplicity. In analogy with equation (40), the Q parameter for the angular 3PCF is defined as Qκ(r, q, ψ) = Z(r, q, ψ) wκ(r)wκ(q) + wκ(r)wκ(|r − q|) + wκ(q)wκ(|r − q|) , where Z = Z 1h + Z 2h + Z 3h . Note that the 2-and 3-halo terms, Z 2h and Z 3h , are computed from the approximations of equations (A6) and (A7), respectively. Figure 20 shows halo model predictions for Qκ with triangle configurations parameterized as in Figure 5. The halos are taken to have NFW profiles. The hierarchical ansatz, Qκ = const., does not hold rigorously over the scales we have considered, as in the 3D case. However, the configuration dependence of Qκ for r = 0.1 ′ is much weaker than for r = 1 ′ and 10 ′ . The strong configuration dependence in the two panels on the right can be understood as follows. First, the results in Figure 5 show that the Q parameter for the 3D mass distribution has a stronger configuration dependence for r > ∼ 1 h −1 Mpc than on strongly non-linear scales. The weak lensing convergence is a weighted projection of the mass distribution. As a result of the lensing projection, a range of 3-dimensional length scales contributes to Qκ at a given θ, which is therefore more sensitive to the 2-and 3-halo terms than the 3D case. This is explicitly shown in Figure 21, which plots the 1and 2-and 3-halo contributions to Qκ for the triangle shapes in Figure 20. For the smallest scale, r = 0.1 ′ , the 1-halo term indeed yields the dominant contribution to Qκ, but the 2-halo term still has a non-negligible contribution (∼ 10%). For r = 1 ′ , the 2-halo term becomes important ( > ∼ 20% contribution). The 3-halo term becomes dominant at for r > ∼ 10 ′ , at which Qκ displays a characteristic, oscillatory feature as predicted by perturbation theory.
In Figure 22 we show Qκ for equilateral triangles as a function of side length. One can readily see that projection effects lead to greater Note that the bias parameters are b g,red = 1.51 and b g,blue = 0.82 for the red and blue galaxies, respectively. The lower panel explores the dependence of Q gal for blue galaxies on the model ingredients. The dot-dashed curve shows the result for the Q gal of blue galaxies, if we suppress the bump feature in N gal (m) at 10 11 ≤ m < ∼ 10 13 h −1 M ⊙ for the calculation. The thin dotted curve shows how the model of Scoccimarro et al. (2001) for the sub-Poisson distribution of blue galaxies changes the prediction at r < ∼ 0.2 h −1 Mpc. Figure 19. The 3PCF (left) and 2PCF (right) of mass, red and blue galaxies.  contributions from the 2-and 3-halo terms, in contrast to the 3D case shown in Figure 7. One caveat that should be noted, as discussed for Figure 7 and 8, is that the halo model might overestimate Q at the transition scale between the quasi-linear and non-linear regimes. This would be reflected in Qκ over the angular scales 0.1 ′ < ∼ r < ∼ 20 ′ , where the 2-halo term is relevant. This uncertainty needs further investigation by comparison with ray-tracing simulation results. The figure also shows that, on very small scales < ∼ 1 ′ , the 1-halo term of Q decreases slightly with decreasing angular scale. The small scale slope of the 1-halo term is determined by the combined effects of the NFW profile, the mass function and the halo concentration in the halo model calculation. For lensing, the small scale slope of the 2PCF and 3PCF is not transparent due to projection effects.
In the lower panel of Figure 22, we plot Qκ for different cosmological models as shown in Figure 9. Qκ is sensitive to cosmological parameters, in particular to Ωm0 parameter, in analogy with the sensitivity of the skewness parameters of weak lensing to Ωm0 (see e.g. Bernardeau et al. 2003;TJ02). Hence, it is expected that measuring Qκ can be used to break the degeneracies in the determinations Ωm0 and σ8 so far from measurements of the shear 2PCF (e.g., Van Waerbeke et al. 2001).

Performance of our approximation for calculating the 1-halo term of the 3PCF
In Appendix B we present a useful approximation for calculating the one halo contribution to the mass 3PCF, which uses the Fourier-space halo model combined with the approximation developed in TJ02 for calculating the skewness and kurtosis parameters. Figure 23 demonstrates the performance of our approximation by comparing the results with the exact values computed from the real-space halo model. We have used the approximation for the 3PCF calculation in the numerator of Q, and the 2PCF calculation in the denominator includes full contributions from the 1-and 2-halo terms. We consider some of the triangle shapes shown in Figure 5. The left and right panels show the results for r = 0.02 and 0.1 h −1 Mpc. In each panel, the solid, dashed and dotted curves are for q/r = 1, 3 and 5, respectively. Note that all the triangles considered are in the strongly non-linear regime since we wish to check the validity of our approximation for the 1-halo term. The thin curves are the exact predictions from the real-space halo model computed from equation (39), while the thick curves denote the results of the approximation (B4). The comparison demonstrates that the approximation works remarkably well with an accuracy of < ∼ 5% for q/r = 1 and 3, and < ∼ 10% even for elongated triangles with q/r = 5. The three dot-dashed curves show the results of each term in the approximation (B4) for q/r = 3. Interestingly, although each term does not work well, the approximation obtained by summing them is close to the exact value. We have confirmed that the approximation of equation (B5) for the 1-halo term of the weak lensing 3PCF also works to about the same accuracy as the 3D case shown here.

DISCUSSION AND CONCLUSION
In this paper, we have used the halo clustering model to compute the 3-point correlation function (3PCF) of cosmological fields. We have shown results for the three-dimensional mass and galaxy distributions and the two-dimensional weak lensing convergence field. The halo model enables us to separately consider the contributions to the 3PCF arising from triplets in a single halo or two or three different halos. Thus we can understand how gravitational clustering transitions from the quasi-linear regime to the strongly non-linear regime as one goes to smaller scales. We found that the single halo contribution, which is dominant on small scales, can be computed using the real-space formulation of the halo model far more easily than the Fourier-space approach used in the literature. We also developed approximations for computing the 2-and 3-halo contributions to the 3PCF. Since measuring the real space 3PCF on small scales is likely to be easier than the bispectrum, our model predictions will allow for the extraction of cosmological information from forthcoming galaxy surveys and weak lensing surveys.
We obtained the following results by applying our halo model to the concordance CDM model (ΛCDM). The 3PCF on small scales r < ∼ 1 h −1 Mpc is dominated by the 1-halo term. Hence it probes properties of massive halos and can be used to constrain halo profiles as discussed below. The quasi-linear 3PCF predicted by perturbation theory can be reproduced by the 3-halo contribution for r > ∼ 3 h −1 Mpc. For plausible halo model parameters, the hierarchical ansatz for the reduced 3PCF parameter, Q = constant, does not hold over the range of scales we have considered. However the Q parameter does have a weaker dependence on triangle configurations in the non-linear regime than in the quasi-linear regime (see Figures 5 and 10). These results qualitatively verify the results in Ma & Fry (2000b,c). Ma & Fry also pointed out that the halo model for plausible model parameters violates the stable clustering hypothesis which has been widely used to develop analytical prediction of the non-linear gravitational clustering, as done in the popular PD fitting formula. In fact, very recently Smith et al. (2002) showed a weak violation of the stable clustering hypothesis using high-resolution N -body simulations. However, our halo model predictions for the 3PCF do not match the asymptotic shape proposed in Ma & Fry (2000b) because of the assumptions employed in the asymptotic formula. Therefore, it it merits further study to carefully investigate how small-scale gravitational clustering can be described by the halo model ingredients. Such a study can be carried out with the halo model methods developed in this paper.
We also found that the non-linear 3PCF is most sensitive to the inner slope of the halo profile and to the halo concentration parameter for given cosmological parameters (see Figures 12, 13 and 14). Combinations of the inner slope and the concentration can be adjusted so that the halo model matches the 2PCF result (see Figure 2). However only one combination can match both the 2-and 3PCFs. Hence in combination with the 2PCF, the non-linear 3PCF could be used to constrain the properties of dark matter halos. We suggest that the use of higher order correlations is a useful way of measuring the properties of massive halos. For example, while lensing surveys of clusters can be used to measure cluster halo profiles directly, this involves identifying cluster centers and assigning masses to clusters to measure averaged profiles. In contrast, by measuring the 2-and 3PCFs, parameters of the halo mass function and profile can be fitted for. While this approach is more challenging computationally, and requires some theoretical assumptions, it treats the data more objectively.
Our halo model predictions match the simulation results of  in the non-linear regime as well as the quasi-linear regime, as shown in Figure 10. However, the halo model seems to be less accurate on the intermediate transition scale ∼ 1 h −1 Mpc. Figure 7 shows that the predicted Q parameter for equilateral triangles has a bump at this scale, which corresponds to the bump found for the reduced bispectrum at k ∼ 1 h Mpc −1 by Scoccimarro et al. (2001). This is unlikely to reflect real properties of dark matter clustering. Rather, we argued that the bump feature is sensitive to the sharp cutoff of halos at the virial radius and the exclusion effect for the 2-and 3-halo terms. We explored alternative prescriptions as shown in Figure 8. To resolve this, detailed calculations and comparisons of the halo model predictions with simulation results will be needed.
We extended the halo model to predict the 3PCF of the galaxy distribution. Once we model how galaxies populate a halo of given mass, the halo occupation number, we can straightforwardly predict the galaxy 3PCF based on the halo model developed here. For the halo occupation number for red galaxies we have used, the galaxy 3PCF has a smaller amplitude and weaker dependence on triangle configuration compared to the mass 3PCF (Figures 17 and 18). This trend is indeed consistent with the actual measurements in . On the other hand, the 3PCF of blue galaxies displays complex features reflecting the input model of the halo occupation number. Thus, the galaxy 3PCF can be used to constrain the galaxy formation scenario as a function of host halo properties and galaxy type. Further work is needed to model the expected properties of galaxies by type for specific surveys and compute the resulting 3PCF.
As an example of the angular 3PCF, we have computed the 3PCF of the weak lensing convergence field. In particular, we employed the real-space halo model to compute the single halo contribution, as for the 3D case, which enabled us to compute the 1-halo term by a 4-dimensional integration. We verified that the real-space halo model is equivalent to the Fourier-space model well studied in the literature (see Figure 3). We also developed approximations for calculating the 2-and 3-halo terms. Because of projection effects, the 2-and 3-halo terms contribute to the 3PCF over a wider range of scales for the angular 3PCF compared to the 3D 3PCF. The resulting 3PCF does not obey the hierarchical ansatz over the angular scales we have considered. The lensing 3PCF is sensitive to cosmological parameters, in particular Ωm0. Comparing measurements with model predictions of the 3PCF can be a useful tool to break degeneracies between the power spectrum and Ωm0. We intend to compute 3-point functions of the shear field following Schneider & Lombardi (2002) and Zaldarriaga & Scoccimarro (2002) as these are easier to measure from the data. In Takada & Jain (2003), we presented a brief investigation of the 3PCFs of shear fields, where we found good agreement between the halo model prediction and the results measured from the ray-tracing simulations. It is also expected that the n-point correlations of weak lensing on sub-arcminute scales contain a wealth of information on properties of massive halos, beyond their dependence on cosmological parameters. The real-space halo model we have developed will be a useful analytical tool for such calculations. The halo model presented in this paper allows for several interesting applications. First, the model can be extended to investigate the effect of triaxial halo shapes on the 3PCF, since it is the lowest-order statistical probe of non-sphericity. So far halo model applications assume a single, spherically symmetric profile. Recently, based on high-resolution simulations, Jing & Suto (2002) showed that halos are more accurately described by triaxial halo profiles than the spherically symmetric NFW profile. They claimed that the axis ratios typically have vales of ∼ 0.6 for the smallest and largest axis. The non-sphericity of the halo profile could lead to a characteristic configuration dependence of the 3PCF. This effect should also affect the interpretation of cosmic shear measurements (Jing 2002). Likewise, the formulation developed in Sheth & Jain (2002) can be used to discuss the effects of substructure within halos on the 2PCF and 3PCF.
The real-space halo model can be directly used to compute the covariance of the 2PCF. As discussed in Cooray & Hu (2001b), the mass distribution on small scales displays pronounced non-Gaussian features induced by non-linear gravitational clustering. Hence one needs to take into account the non-Gaussian errors arising from the 4-point correlation function (4PCF). On the scales of interest the single halo contribution should dominate the error. The real-space halo model allows us to analytically compute the error contribution arising from the 4PCF with no additional computational effort than for the 3PCF (see §2.1).
In Appendix B, we constructed an approximation for calculating the 1-halo term of the 3PCF based on the Fourier-space halo model. We showed that the the approximation accurately describes the amplitude and configuration dependence of the 3PCF (see Figure 23). We will employ this approximation to perform an analytical study of the pairwise peculiar velocity dispersion (PVD) within the BBGKY hierarchy formalism (Peebles 1980). The PVD can be measured through the redshift distortions inferred from galaxy surveys (e.g., Zehavi et al. 2002). The BBGKY picture tells us that the PVD arises mainly from the 3PCF on small scales. However, there has been no comprehensive analytical model to describe the non-linear PVD. This is because of the complex form of the BBGKY hierarchy equations. Peebles (1980) (see Mo, Jing & Börner 1997; for a detailed study) assumed the hierarchical form for the non-linear 3PCF, although it turns out to be violated for CDM models. The isothermal assumption for the velocity distribution within a given halo was employed to analytically obtain the PVD based on the halo model . It is important to clarify whether the BBGKY hierarchy leads to a self-consistent PVD for the CDM model. Moreover, we can easily combine the halo model prediction with models of the halo occupation number of galaxies to predict the PVD for different galaxy types, following the approach in §4.2.3.
We would like to thank R. Sheth, R. Scoccimarro, I. Szapudi, A. Taruya and L. Hui for several valuable discussions. We thank Y. P. Jing for kindly providing us with his simulation data results. Helpful comments from the referee, Chung-Pei Ma, led to improvements in the paper. This work is supported by NASA grants NAG5-10923, NAG5-10924 and a Keck foundation grant.

APPENDIX A: APPROXIMATION FOR THE 2-AND 3-HALO TERMS OF THE 3PCF
In this appendix, we present approximations used for the predictions of the 2-and 3-halo contributions to the 3PCF of the density field. The approximations are based on the Fourier-space halo model, combined with an approximation for the configuration dependence of the bispectrum from Scoccimarro et al. (2001; see also TJ02).

A1 3D 3-point correlation function
Following the Fourier-space halo model, as discussed in §2.2, the 2-halo term for the mass 3PCF can be expressed from equations (11) and (13) as ζ 2h (r1, r2, r3) = dmn(m) m ρ0 where the bias parameter b(m) is given by equation (20). As discussed in Scoccimarro et al. (2001), the physical meaning of the 2-halo term tells us that the main contribution to the first term on the r.h.s arises from modes with k1 ≫ k2, leading to the approximation k12 ≈ k1. Replacingũm(k12) withũm(k1) allows us to perform analytically the angular integration in d 3 ki in the first term above. Using similar procedures for the second and third terms on the r.h.s. yields the following approximation: ζ app 2h (r1, r2, r3) = dm n(m) m ρ0 2 b(m) k 2 1 dk1 2π 2 (ũm(k1)) 2 j0(k1r12) where we have used the formula dΩ k exp[ik · r] = j0(kr). The large square bracket is meant to show explicitly that the integrations inside each bracket can be done separately from the other. Thus one needs to perform only 2-dimensional integrations to get the 2-halo term. It was shown that a similar approximation for the skewness calculation turns out to be accurate TJ02). Hence, we expect that this also holds for the 3PCF calculation. The 3-halo contribution to the 3PCF can be expressed as (2π) 3 um(k1)u m ′ (k2)u m ′′ (k12)e ik 1 ·r 13 +ik 2 ·r 23 P L (k1)P L (k2) × 10 7 + 1 k 2 1 + 1 k 2 2 k1 · k2 + 4 7 (k1 · k2) 2 k 2 1 k 2 2 + cyc. .

(A5)
We find that this approximation reproduces the perturbation theory result at scales > ∼ 5 h −1 Mpc for the ΛCDM model (see Figure 7). Note that, for the actual predictions of the 3PCF shown in this paper, we have used the permutation (1) ↔ (2) for terms such as Ξ (2) Ξ (1) in the above equation, so that it satisfies statistical symmetry for permutations between r1, r2 and r3.
By replacing (m/ρ0) in equations (A2) and (A4) with the halo occupation number of galaxies N gal (m)/n gal , as discussed in §4.2.3, we obtain the corresponding approximations for computing the 2-and 3-halo contributions to the 3PCF of galaxies.

A2 2D 3-point correlation function
It is straightforward to extend the approximations for the 3-dimensional 3PCF to the 3PCF of weak lensing fields. As discussed for equation (36), combining equation (A2) with Limber's equation yields the following approximation for the 2-halo term of the weak lensing convergence: Z app 2h (r1, r2, r3) = where ki = li/dA(χ). Note that the redshift dependence of P L (k; χ) is given by the linear growth factor D(z). Likewise, the approximation for the 3-halo term can be obtained by the following replacements in equation (A4) and the additional integration over the redshift: where µ1323 ≡ (r13 · r23)/(r13r23) and

APPENDIX B: APPROXIMATIONS FOR THE 1-HALO TERM OF THE 3PCF
In this section, we present approximations for calculating the 1-halo contribution to the 3D or 2D 3PCF with the Fourier-space halo model. These are combined with the approximation in TJ02 for calculating the skewness parameter of the density field.
The above approximation reduces the 7-dimensional integral in equation (B1) to a 3-dimensional one. However, the resulting 3PCF depends only on the two side lengths r13 and r23, and thus has no dependence on the interior angles. If we had started by eliminating k2 or k1 in equation (B1), the resulting 3PCF would have dependences only on the two side lengths r12 and r23 or r12 and r13, respectively. Thus, the approximate result for the 3PCF has different forms even for the same triangle shape. To resolve this, we propose the following symmetrized form as the approximation for calculating ζ 1h : The dependences of this form on the three parameters of the triangle shape, r12, r23 and r31, are like the hierarchical form ζ = Q[ξ(r12)ξ(r31) + ξ(r12)ξ(r23) + ξ(r13)ξ(r22)]. The factor 1/3 is simply chosen so that it agrees with the third-order moment given in TJ02 for the limit r, q → 0. As shown in §4.4, the approximation (B4) is accurate compared with the exact value computed from equation (6), even though each term on the r.h.s is far from the exact value (see Figure 23).