format: text

text:
# Untitled Document

THE TWO-EXPONENTIAL LIOUVILLE THEORY AND

THE UNIQUENESS OF THE THREE-POINT FUNCTION

L. O’Raifeartaigh1 e-mail: lor@stp.dias.ie, J. M. Pawlowski2 e-mail: jmp@stp.dias.ie, and V. V. Sreedhar3 e-mail: sreedhar@stp.dias.ie

School of Theoretical Physics

Dublin Institute for Advanced Studies

10 Burlington Road, Dublin 4, Ireland

Abstract

It is shown that in the two-exponential version of Liouville theory the coefficients of the three-point functions of vertex operators can be determined uniquely using the translational invariance of the path integral measure and the self-consistency of the two-point functions. The result agrees with that obtained using conformal bootstrap methods. Reflection symmetry and a previously conjectured relationship between the dimensional parameters of the theory and the overall scale are derived.

Introduction: Although the quantisation of two-dimensional Liouville theory with a potential of the form $`V_b=\mu e^{2b\stackrel{~}{\varphi }(x)}`$, where $`\stackrel{~}{\varphi }(x)`$ is a scalar field and $`\mu `$ and $`b`$ are constants, has been widely investigated , it still presents some problems. What is perhaps the most disturbing problem is the following: The three point functions of vertex operators $`\mathrm{exp}[2\alpha _I\stackrel{~}{\varphi }(x_I)]`$ play a central role in the theory, in the sense that all $`N`$-point functions can be obtained from these by integration, and have the following form

$$𝒢_3=𝐂_3\underset{I=1}{\overset{3}{}}|\frac{x_Ix_J}{L}|^{\mathrm{\Delta }_{(IJ)}(\alpha )}\text{where}𝐂_3=\left(g_b(\xi )Z_0^b(\xi )\right)\left(\frac{K(\xi ,\alpha _I)}{k^{}(\xi )}\right)$$
$`(1)`$
Here $`L`$ is a constant scale, the $`\mathrm{\Delta }_{(IJ)}`$ are known combinations of constant conformal weights, $`\xi =q_{I=1}^3\alpha _I`$, and $`q=b+1/b`$, while $`Z_0^b(\xi )`$ is related to the zero-mode integration and can be computed explicitly. $`K`$ and $`k`$ are the functions defined in (26). The problem is that $`𝐂_3`$ can actually be computed only at the points $`\xi =mb`$ for $`mZ_+`$ and the extrapolation from these points to general values of $`\xi `$ leaves a factor $`g_b(\xi )`$ in (1) undetermined. So far it has not been possible to determine this factor from first principles. Instead, what has been done is to make the so-called DOZZ Ansatz that

$$𝐂_3=\frac{K(\xi ,\alpha _I)}{k(\xi )}h_b(\xi )g_bZ_0^b\frac{k}{k^{}}=1$$
$`(2)`$
and check that the resulting three-point function satisfies some reasonable, but extra, conditions such as reflection symmetry and crossing symmetry. The latter check is the decisive one, but uses some special four-point functions .

In a previous paper  it was suggested that the gap in the direct derivation of $`h(\xi )`$ could be closed by using a potential of the form $`V(\varphi )=\mu _be^{2b\stackrel{~}{\varphi }}+\mu _ce^{2c\stackrel{~}{\varphi }}`$, where $`bc=1`$, – this being the most general potential one can use in the path integral context for a conformal field theory of a single real scalar field without derivative interactions  – rather than the standard form $`V_b`$. As in the one-exponential theory, the $`𝐂_3`$ in the two-exponential theory can be calculated only at a discrete, but much larger, set of points, namely $`\xi _{mn}=mb+nc`$ and the equation corresponding to (1) takes the form

$$𝐂_3=\left(g(\xi )Z_0(\xi )\right)\left(\frac{K(\xi ,\alpha _I)}{k^{}(\xi )}\right)$$
$`(3)`$
where $`g(\xi )`$ is a priori undetermined, the functions $`K`$ and $`k`$ are exactly the same as in (1) but $`Z_0`$ is not. It was shown in  that, subject to two conditions, the factor $`g(\xi )`$ could be fixed and thus the DOZZ Ansatz (2) could be derived. The conditions were (a) that the dimensional parameters $`\mu _b`$ and $`\mu _c`$ were related to the overall scale $`L`$ of the system by an equation of the form $`\mathrm{\Omega }(\mu _b,\mu _c,L)=1`$, where $`\mathrm{\Omega }`$ is the function defined in (25), and (b) that $`h(\xi )`$, defined as $`h(\xi )=g(\xi )Z_0(\xi )k(\xi )/k^{}(\xi )`$, had no singular points.

In the present paper we refine these results considerably and extend the analysis of the two-exponential theory. In particular we show that, subject to the mild technical condition given in (37), the DOZZ Ansatz follows from the translational invariance of the path integral measure and the self-consistency of the two-point functions. We also present an extrapolation of the fluctuating part of the path integral from the lattice points $`\xi _{mn}`$, for which the $`Z_0`$ part of the path integral is automatically

$$Z_0(\xi )=\frac{k^{}(\xi )}{k(\xi )}$$
$`(4)`$
This means that the zero mode integral for the two-exponential theory produces exactly the factor $`k^{}/k`$ that was postulated in the DOZZ Ansatz.

As in , we use the path-integral formalism; and as the symmetries of the path-integral and the associated sum rules are of interest in their own right, a secondary purpose of the paper is to present these in a systematic way. The most important symmetries are those connected with the translational invariance of the measure and conformal (Weyl) invariance.

A by-product is an analysis of the two-point function $`𝒢_2`$. Although $`𝒢_2`$ cannot be defined directly because of conformal invariance, it can be defined both as a limit of $`𝒢_3`$ when one of the $`\alpha `$’s becomes zero, and as a volume integral of $`𝒢_3`$ when the corresponding $`\alpha `$ is $`b`$ or $`c`$. The compatibility of the two definitions and the sum rule mentioned above lead to a simple linear homogeneous sum rule for the quantity $`h(\xi )`$. Together with the boundary conditions $`h(\xi _{mn})=\mathrm{\Omega }^{m+n}`$, obtained by direct computation, this sum rule fixes $`h(\xi )=1`$ uniquely. It also fixes $`\mathrm{\Omega }=1`$ which is the relation between the dimensional constants $`\mu _b`$ and $`\mu _c`$ and the overall scale $`L`$ that was conjectured in an earlier paper. The corresponding sum rule and boundary conditions in the one-exponential theory would imply only that $`h(\xi )`$ is periodic. If $`𝒢_2`$ is interpreted as an inner-product of states, then $`h(\xi )=1`$ implies that for each conformal weight a zero-norm state decouples, so that there is only one physical state. The decoupling is equivalent to reflection invariance, which therefore emerges as an output.

The generating functional: The generating functional of the two-exponential Liouville theory is defined as

$$Z[J]=[d\stackrel{~}{\varphi }]e^{{\scriptscriptstyle d^2x\sqrt{g(x)}\left[{\scriptscriptstyle \frac{1}{4\pi }}\stackrel{~}{\varphi }\mathrm{\Delta }\stackrel{~}{\varphi }+{\scriptscriptstyle \frac{q}{4\pi }}\stackrel{~}{\varphi }+\mu _be^{2b\stackrel{~}{\varphi }}+\mu _ce^{2c\stackrel{~}{\varphi }}J(x)\stackrel{~}{\varphi }(x)\right]}}$$
$`(5)`$
where $`\mathrm{\Delta }`$ is the Laplace-Beltrami operator, $``$ is the Ricci scalar, the coupling constants $`\mu _b,\mu _c`$ have dimensions of mass squared, and $`b`$ and $`c`$ are dimensionless constants. For the $`N`$-point functions of vertex operators we have $`e^{{\scriptscriptstyle d^2x\sqrt{g}J\stackrel{~}{\varphi }}}=_{I=1}^N\nu _I\sqrt{g}e^{2\alpha _I\stackrel{~}{\varphi }}`$, but we need not yet specialise to this case. The $`\nu _I`$ are constants which, like the $`\mu _b`$ and $`\mu _c`$, have to be renormalised because of the short distance singularities of the Green’s function defined by the equation $`\mathrm{\Delta }G(x,y)=\frac{\pi }{\sqrt{g}}\delta ^2(xy)`$. As explained in detail in , the Green’s function for the non-coincident and coincident arguments is defined as:

$$G(x,y)=\frac{1}{2}\text{ln}\frac{xy}{L}+𝒪(V^1)\text{and}G(x,x+dx)\frac{1}{2}\mathrm{ln}[\frac{ds}{L}]+\frac{1}{4}\mathrm{ln}\sqrt{g(x)}$$
$`(6)`$
respectively, where $`ds`$ is the infinitesimal geodesic separation, $`L`$ sets the overall scale, and $`V`$ is the volume of the space.

As the Green’s function $`G(x,x)`$ occurs only in the form $`\alpha ^2G(x,x)`$ at each vertex, its infinite part may be absorbed by renormalising $`\mu _b`$, $`\mu _c`$ and $`\nu _I`$. Let us define the renormalisations of $`\mu _b`$, $`\mu _c`$ and $`\nu _I`$ to be

$$\mu _b\mu _b(\mathrm{\Lambda }\frac{ds}{L})^{2(b^2qb)},\mu _c\mu _c(\mathrm{\Lambda }\frac{ds}{L})^{2(c^2qc)},\nu _I\nu _I(\mathrm{\Lambda }\frac{ds}{L})^{2(\alpha _I^2q\alpha _I)}$$
$`(7)`$
where $`\mathrm{\Lambda }`$ is the dimensionless ultra-violet renormalisation scale. Anticipating the result that $`b+c=q,bc=1`$ from Weyl invariance, it is easy to see that $`\mu _b`$ and $`\mu _c`$ have the naive scaling dimensions with respect to $`\mathrm{\Lambda }`$.$`^1`$ There is arbitrariness in the definition of renormalisation associated with the translational invariance of the path integral measure. However, this arbitrariness does not affect equation (8). Clearly the $`b^2`$, $`c^2`$, $`\alpha ^2`$ powers of $`ds`$ in the above equation cancel the infinities coming from the Green’s functions at coincident points. Counting the powers of the remaining factors of $`ds`$ and $`\mathrm{\Lambda }`$, we find, after a little algebra, that the renormalised functional integral is proportional to

$$[\frac{ds}{L}]^{2q^2}\underset{I=1}{\overset{N}{}}\mathrm{\Lambda }^{2(\alpha _I^2q\alpha _I)}$$
$`(8)`$
Since the factor $`[ds/L]^{q^2}`$ may be absorbed in the Polyakov conformal anomaly term which is proportional to $`q^2`$, we see that the $`\mathrm{\Lambda }`$ dependence of the renormalised functional integral cancels except for the contribution coming from the external sources.

The conformal weights $`\mathrm{\Delta }_\alpha `$ of the fields exp \[$`2\alpha _I\stackrel{~}{\varphi }(x_I)`$\] may be read off directly by varying the path integral with respect $`\sqrt{g(x_I)}`$. Taking the $`\sqrt{g}`$ dependence of $`\alpha _I^2G(x_I,x_I)`$ and the cross-term $`q\alpha _Id^2xd^2y\sqrt{g(x)}(x)G(x,y)\delta (yx_I)`$ in the Gaussian integration into account, this yields $`\mathrm{\Delta }_\alpha =\alpha (q\alpha )`$. From (7) it is clear that the scaling dimensions of $`\nu _I`$ are minus the conformal weight of corresponding field. Thus fields of the same conformal weight are renormalised in the same way.

A sum rule for the generating functional: If we now specialise to the case where the external current is of the form $`_{I=1}^N\nu _Ie^{{\scriptscriptstyle d^2x\sqrt{g}J\stackrel{~}{\varphi }}}=_{I=1}^N\nu _I\sqrt{g}e^{2\alpha _I\stackrel{~}{\varphi }}`$, the renormalised functional integral is defined by

$$Z[J]=[d\stackrel{~}{\varphi }]e^{{\scriptscriptstyle d^2x\sqrt{g(x)}\left[{\scriptscriptstyle \frac{1}{4\pi }}\stackrel{~}{\varphi }\mathrm{\Delta }\stackrel{~}{\varphi }+{\scriptscriptstyle \frac{q}{4\pi }}\stackrel{~}{\varphi }+\sqrt{g}^{b^2}\mu _be^{2b\stackrel{~}{\varphi }}+\sqrt{g}^{c^2}\mu _ce^{2c\stackrel{~}{\varphi }}\right]}}\underset{I=1}{\overset{N}{}}\sqrt{g}^{\alpha _I^2}\nu _Ie^{2\alpha _I\stackrel{~}{\varphi }(x_I)}$$
$`(9)`$
where we are using the renormalisation prescription in (7) and we have absorbed, for convenience, the ln$`\sqrt{g}`$ part of $`G(x,x)`$ by letting $`e^{2\alpha \varphi }\sqrt{g}^{\alpha ^2}e^{2\alpha \varphi }`$. It is then understood that $`G_R(x,x)=0`$.

The translation invariance of the path integral measure can now be formulated as follows:

$$[d\stackrel{~}{\varphi }]\frac{\delta }{\delta \stackrel{~}{\varphi }(x)}\left[e^{S[\stackrel{~}{\varphi }]+{\scriptscriptstyle d^2x\sqrt{g}J(x)\stackrel{~}{\varphi }(x)}}\right]=0$$
$`(10)`$
where $`\frac{\delta }{\delta \stackrel{~}{\varphi }(x)}`$ is the generator of translations on the space of fields. From (10) we derive the quantum equations of motion for $`Z[J]`$:

$$\frac{1}{4\pi }\mathrm{\Delta }\frac{\delta Z[J]}{\delta J(x)}=\left(b\mu _b\sqrt{g}^{b^2}Z[J_{b,x}]+c\mu _c\sqrt{g}^{c^2}Z[J_{c,x}]\right)+\frac{1}{2}\left(J(x)\frac{q}{4\pi }(x)\right)Z[J]$$
$`(11)`$
where

$$J_{b,x}(y)=J(y)+2b\frac{1}{\sqrt{g(y)}}\delta ^2(yx),J_{c,x}(y)=J(y)+2c\frac{1}{\sqrt{g(y)}}\delta ^2(yx)$$
$`(12)`$
Carrying out the same procedure for the constant, zero-mode, part $`\varphi _0`$ of $`\stackrel{~}{\varphi }`$, we obtain the integrated form of (11), namely

$$d^2x\sqrt{g}\left(b\mu _b\sqrt{g}^{b^2}Z[J_{b,x}]+c\mu _c\sqrt{g}^{c^2}Z[J_{c,x}]\right)=\frac{1}{2}\left(J_0q\chi \right)Z[J]$$
$`(13)`$
where

$$J_0=d^2x\sqrt{g}J(x)\text{and}\chi =\frac{1}{4\pi }d^2x\sqrt{g}(x)$$
$`(14)`$
$`\chi `$ being the Euler number of the underlying manifold. Eq. (13) embodies the sum rule for the generating functional.

Weyl transformations: A local Weyl transformation can be performed by varying the generating functional with respect to $`\sqrt{g(x)}`$ where $`xx_I`$, the external points. For $`Z[J]`$ in (9), we find

$$\frac{\delta Z}{\delta \sqrt{g}}=(1+b^2)\sqrt{g}^{b^2}\mu _bZ[J_{b,x}]+(1+c^2)\sqrt{g}^{c^2}\mu _cZ[J_{c,x}]+\frac{q}{4\pi }\mathrm{\Delta }\frac{\delta Z[J]}{\delta J(x)}$$
$`(15)`$
The third term may be eliminated using (11) to get

$$\left(1+b^2qb\right)\sqrt{g}^{b^2}\mu _bZ[J_{b,x}]+\left(1+c^2qc\right)\sqrt{g}^{c^2}\mu _cZ[J_{c,x}]\left(\frac{q}{4\pi }(x)J(x)\right)\frac{q}{2}Z(J)$$
$`(16)`$
The Weyl condition is that the variation $`\frac{\delta Z}{\delta \sqrt{g}}`$ should be proportional to the external current namely $`\frac{q}{4\pi }J`$. Since this has to be valid for all currents $`J`$, the condition for Weyl invariance is that the first two terms must vanish and we have

$$q=(b+c)\text{and}bc=1,$$
$`(17)`$
The above approach may be contrasted with the one in  where (17) was derived only after the path integral was evaluated.

The $`N`$-point functions: For the computation of general $`N`$-point functions of vertex operators, we let the underlying manifold be a two dimensional sphere. In that case, $`\chi =2`$ and there is only one zero-mode for $`\stackrel{~}{\varphi }`$, namely the constant $`\varphi _0`$. As explained in detail in , the expression for the $`N`$-point function may be simplified by using a Sommerfeld-Watson transform  for the exponential of an integrated vertex operator. With $`\stackrel{~}{\varphi }=\varphi _0+\varphi `$ the resulting expression for the $`N`$-point function takes the form

$$𝒢_N=𝑑\varphi _0\frac{dudv}{\mathrm{\Gamma }(1+iu)\mathrm{\Gamma }(1+iv)}\frac{e^{2(i(bu+cv)\xi _N)\varphi _0}}{\text{sinh}\pi u\text{sinh}\pi v}\times 𝐂_N(iu,iv)$$
$`(18)`$
where

$$𝐂_N(iu,iv)=𝑑\varphi U_b^{iu}U_c^{iv}e^{{\scriptscriptstyle d^2x\sqrt{g}\left[{\scriptscriptstyle \frac{1}{4\pi }}\varphi \mathrm{\Delta }\varphi +{\scriptscriptstyle \frac{q}{4\pi }}R\varphi \right]}}𝚷_N\text{with}\xi _N=q\underset{I=1}{\overset{N}{}}\alpha _I$$
$`(19)`$
and

$$𝚷_N=\underset{I=1}{\overset{N}{}}\sqrt{g}^{\alpha _I^2}\psi _Ie^{2\alpha _I\stackrel{~}{\varphi }(x_I)},U_b(\varphi )=\mu _bd^2x(\sqrt{g})^{qb}e^{2b\varphi }$$
$`(20)`$
and similarly for $`bc`$. The integral $`𝐂_N(iu,iv)`$ is a Gaussian path integral for the fluctuations $`\varphi `$ which, for $`iu=m`$ and $`iv=n`$, $`m`$ and $`n`$ being positive integers, can be done in a straightforward manner and produces an ordinary multiple integral.

In terms of the $`N`$-point functions $`𝒢_N(x_I,\alpha _I)`$, the sum rule (13) takes the form

$$d^2x\sqrt{g}\left[b\mu _b𝒢_{N+1}(x_I,x,\alpha _I,b)+c\mu _c𝒢_{N+1}(x_I,x,\alpha _I,c)\right]=\xi _N𝒢_N(x_I,\alpha _I)$$
$`(21)`$
and thus relates the (integrated) $`N+1`$-point function to the $`N`$-point function. Note that the above equation requires that $`\xi _N0`$ because $`𝒢_N\mathrm{}\text{as}\xi _N0`$, making the right hand side indefinite.

The Three-point function: As is well-known, the three-point function is the lowest $`N`$-point function for which conformal invariance does not require the extraction of an infinite group volume factor. If we choose $`\xi `$ to be pure imaginary and integrate over the zero-mode $`\varphi _0`$ in (18) we obtain a delta function $`\delta (\xi bucv)`$, in which case the coefficient $`𝐂_3`$ may be written as

$$𝐂_3(\xi ,i(u+v))=𝑑\varphi e^{\frac{1}{4\pi }{\scriptscriptstyle d^2x\sqrt{g}[\varphi \mathrm{\Delta }\varphi +qR\varphi ]}}U_b^{iu}U_c^{iv}𝚷_\mathrm{𝟑}$$
$`(22)`$
Apart from the spectator variables $`\alpha _I\alpha _J`$, we see that, due to the delta-function, $`𝐂_3`$ is a function of only two variables, chosen as $`\xi `$ and $`u+v`$ for convenience. This is the great advantage of using the Sommerfeld-Watson transform. We then have in the infinite volume limit,

$$𝒢_3=\frac{dudv}{\mathrm{\Gamma }(1+iu)\mathrm{\Gamma }(1+iv)}\frac{\delta (\xi ibuicv)}{\text{sinh}\pi u\text{sinh}\pi v}𝐂_3(\xi ,i(u+v))$$
$`(23)`$
The problem is that $`𝐂_3`$ can only be computed at the points $`u,v=im,in`$ for $`m,nZ_+`$, where, as shown in , it is given by, for $`\xi _{mn}bm+cn`$,

$$𝐂_3(\xi _{mn},m+n)=(1)^{m+n}m!n!\lambda ^{\xi _{mn}}\mathrm{\Omega }^{m+n}\left(\frac{K(\xi _{mn},\alpha _I)}{k^{}(\xi _{mn})}\right)\underset{I=1}{\overset{3}{}}\frac{x_{IJ}}{L}^{\mathrm{\Delta }_{IJ}(\alpha )}$$
$`(24)`$
where

$$\lambda =\left(\frac{\mu _b\mathrm{\Phi }_b}{\mu _c\mathrm{\Phi }_c}\right)^{\frac{1}{bc}},\mathrm{\Omega }^{cb}=\frac{(\mu _bL^2\mathrm{\Phi }_b)^c}{(\mu _cL^2\mathrm{\Phi }_c)^b},\mathrm{\Phi }_b=\pi \gamma (b^2)(b^2)^{2qb},\gamma (x)=\frac{\mathrm{\Gamma }(x)}{\mathrm{\Gamma }(1x)}$$
$`(25)`$
The functions $`k(\xi )`$ and $`K(\xi ,\alpha _I)`$ are defined by the following equations:

$$\mathrm{ln}k(\xi )=_0^{\mathrm{}}\frac{dt}{t}\left((\frac{q}{2}x)^2e^{2t}\frac{\text{sinh}^2(\frac{q}{2}x)t}{\text{sinh}bt\text{sinh}ct}\right)\text{and}K(\xi ,\alpha _I)=k^{}(0)\underset{I=1}{\overset{3}{}}\frac{k(2\alpha _I)}{k(\xi +2\alpha _I)}$$
$`(26)`$
To extrapolate $`𝐂_3`$ to other values of $`\xi `$ we note that although the arguments of $`k`$’s are constrained by the relation $`_I\alpha _I=q\xi _{mn}`$, they range over the whole real axis for fixed $`\xi _{mn}`$. Hence the only reasonable extrapolation is the obvious one, $`k(\xi _{mn},\alpha _I)k(\xi ,\alpha _I)`$, in which case $`K(\xi _{mn},\alpha _I)K(\xi ,\alpha _I)`$. Unfortunately, this argument is not valid for the rest of $`𝐂_3(\xi _{mn},m+n)`$, which depends only on $`\xi _{mn}`$. Thus the most general extrapolation of (24) is

$$𝐂_3(iu,iv,\xi )=\text{cosh}\pi u\text{cosh}\pi v\mathrm{\Gamma }(1+iu)\mathrm{\Gamma }(1+iv)\lambda ^\xi \mathrm{\Omega }^{i(u+v)}\frac{K(\xi ,\alpha _I)}{k^{}(\xi )}f(\xi ,u+v)\frac{x_{IJ}}{L}^{\mathrm{\Delta }_{(IJ)}}$$
$`(27)`$
where the cosh terms take care of the $`(1)^{m+n}`$ terms and $`f(\xi ,u+v)`$ is an arbitary function with $`f(\xi =mb+nc,m+n)=1`$. Then $`𝒢_3`$ becomes

$$𝒢_3=\lambda ^\xi h(\xi )\frac{K(\xi ,\alpha _I)}{k(\xi )}\underset{I=1}{\overset{3}{}}\frac{x_{IJ}}{L}^{\mathrm{\Delta }_{(IJ)}(\alpha _I)}$$
$`(28)`$
where

$$h(\xi )=\frac{k(\xi )}{k^{}(\xi )}Z_0(\xi )\text{and}Z_0(\xi )=\frac{dudv\delta (\xi ibuicv)}{\text{tanh}\pi u\text{tanh}\pi v}\mathrm{\Omega }^{i(u+v)}f[\xi ,i(u+v)]$$
$`(29)`$
Since the function $`f`$ is unknown, we cannot proceed from (29). Hence we take an alternative route using the two-point function.

Uniqueness (An Application of the Sum Rule): For two and three point functions the sum rule (21) in the infinite volume limit is

$$b\mu _bd^2x𝒢_3(x_I,x,\alpha _I,b)+c\mu _cd^2x𝒢_3(x_I,x,\alpha _I,c)=\xi _2𝒢_2(x_I,\alpha _I),\xi _20$$
$`(30)`$
This is not useful unless we have an alternative definition for the two-point function. Such a definition may be obtained by regarding it as twice$`^2`$ Actually any constant independent of $`\xi `$ is allowed a priori, but the requirement that the sum rule be satisfied at the points $`\xi =mb+nc`$ fixes the constant to be two. the limit of the three-point function as $`\alpha _30`$ and $`x_3\mathrm{}`$. It is easy to see that the limit is non-zero only if $`\mathrm{\Delta }_1=\mathrm{\Delta }_2`$ i.e if $`\alpha _1=\alpha _2`$ or $`\alpha _1=q\alpha _2`$. Since for scattering states $`\alpha _I=\frac{q}{2}+i\beta _I`$, the quantities $`\alpha _1\alpha _2`$ and $`\xi _2=q\alpha _1\alpha _2`$ are pure imaginary, and we obtain

$$𝒢_2(\alpha _1,\alpha _2;x_{12})=4\pi \lambda ^{\xi _2}\frac{x_{12}}{L}^{(\mathrm{\Delta }_1+\mathrm{\Delta }_2)}h(\xi _2)R(\xi _2)[\delta (\beta _1\beta _2)+\delta (\beta _1+\beta _2)]$$
$`(31)`$
where $`R(\xi )=k(\xi )/k(\xi )`$.

To compute the left hand side of (30) we note that the first integral is

$$\lambda ^{\xi _2b}h(\xi _2b)\frac{K(\alpha _1,\alpha _2,b)}{k(\xi _2+b)}|\frac{x_{12}}{L}|^{2(\mathrm{\Delta }_b\mathrm{\Delta }_1\mathrm{\Delta }_2)}\times \zeta $$
$`(32)`$
where

$$\zeta =d^2x_3|\frac{x_{31}}{L}|^{2(\mathrm{\Delta }_b+D)}|\frac{x_{23}}{L}|^{2(\mathrm{\Delta }_bD)}=2\pi ^2\frac{x_{12}}{L}^{2\mathrm{\Delta }_b}\delta (D)$$
$`(33)`$
and $`\mathrm{\Delta }_b=1`$ and $`D=\mathrm{\Delta }_1\mathrm{\Delta }_2=\beta _1^2\beta _2^2`$. The delta functions in these equations and the condition $`\xi _20`$ in (30) mean that we only need the coefficient for $`\alpha _1=\alpha _2\alpha `$, which is easily computed to be

$$\frac{K(\alpha ,\alpha ,b)}{k(\xi _2+b)}=\xi _2^2k^{}(0)R(\xi _2)\frac{k(2b)}{k^2(b)}=\xi _2^2R(\xi _2)\frac{\mathrm{\Phi }_b}{\pi b},\xi _2=q2\alpha $$
$`(34)`$
Inserting these formulae into the sum rule (30) and using the identity $`\lambda ^b=\mu _b\mathrm{\Phi }_b\mathrm{\Omega }`$, and similarly for $`c`$, we get the sum rule

$$h(\xi _2+b)+h(\xi _2+c)=2\mathrm{\Omega }h(\xi _2)\text{with}h(mb+nc)=\mathrm{\Omega }^{m+n}$$
$`(35)`$
for the unknown function $`h(\xi _2)`$. The second equation is obtained by explicit computation from (24). The question is whether the equations in (35) determine $`h(\xi _2)`$ uniquely. We show that they do subject to a mild technical assumption to be introduced presently. First let us restrict ourselves to the case when $`b`$ is fractional i.e. $`b=r/s`$ where $`r>s`$ are positive integers with no common factor and lowest common multiple $`rs`$. Next we rewrite the sum rule (35) in terms of new variables defined as follows:

$$ye^\xi h(e^{\frac{r}{s}}y)+h(e^{\frac{s}{r}}y)=2\mathrm{\Omega }h(y)$$
$`(36)`$
We now make the assumption3 It is probable that this could be actually proved directly from the sum rule (35) which rules out a whole class of functions.

$$\underset{y0}{lim}\frac{h\left(e^{(\frac{r}{s})^{\pm 1}}y\right)}{h(y)}C_\pm <\mathrm{}$$
$`(37)`$
Considering the restriction

$$y_t=e^{\frac{t}{rs}}h(e^{\frac{tr^2}{rs}})+h(e^{\frac{ts^2}{rs}})=2\mathrm{\Omega }h(e^{\frac{t}{rs}}),\text{where}tZ_+$$
$`(38)`$
we see that the symmetry group we need to implement is the dilatation group. Thus we may expand a solution of (38) as follows:

$$h(y_t)=\underset{p=0}{\overset{N}{}}C_py_t^{\sigma _p}\text{where}e^{r^2\omega \sigma _p}+e^{s^2\omega \sigma _p}2\mathrm{\Omega }=0\text{and}\omega =\frac{1}{rs}$$
$`(39)`$
$`N(0Nr^2)`$ being a finite number. This follows because the sum rule restricts the number of $`C_p`$s to be the dimension of the solution space of the polynomial equation in (39). This is of course true only on the first sheet of the covering of the $`y`$ variable. In the general case we would also have a sum over the infinite number of coverings.

$$h(y)=\underset{p=0}{\overset{N}{}}\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}C_{p_n}y^{\sigma _{p_n}}\text{where}\sigma _{p_n}=\sigma +2\pi inrs$$
$`(40)`$
We now use the second equation in (35) to write, for the special values $`\xi =mb+nct=mr^2+ns^2`$,

$$\underset{p=0}{\overset{N}{}}C_py_t^{\sigma _p}=\mathrm{\Omega }^{m+n}C_N=1,\sigma _N=0,C_{pN}=0\mathrm{\Omega }=1$$
$`(41)`$
The latter follows as the only consistent solution since the finite sum on the left hand side cannot in general match the right hand side which can be made arbitrarily large by letting $`m`$ and $`n`$ tend to infinity. Note that the covering does not play a role in the discussion of (41). It then follows that

$$h(\frac{t}{rs})=1tZ_+$$
$`(42)`$
Thus for $`b=r/s`$ the function $`h(\xi )`$ is unity when restricted to the subset of points $`\xi =t/rs`$.

This is as far as we can go for a fixed $`b=r/s`$. However we now invoke the fact that (29) implies that $`h(\xi )`$ is continuous in $`b`$ and $`\xi `$. We then define $`\stackrel{~}{b}=\stackrel{~}{r}/\stackrel{~}{s}`$ where $`\stackrel{~}{r}=lr+1`$ and $`\stackrel{~}{s}=ls`$. It follows that $`\stackrel{~}{r}`$ and $`\stackrel{~}{s}`$ have no common factor, $`\stackrel{~}{b}b=1/ls`$ and, by applying the above result to the tilde-variables, $`h(t/l^2sr(1+1/lr)=1`$. But this means that in any $`1/l`$ neighbourhood of $`b`$ there is a $`\stackrel{~}{b}`$ for which $`h(\xi )=1`$ at points which are separated by distances of order $`1/l^2`$. As these distances tend to zero as $`l`$ tends to infinity, we see that this is compatible with the continuity of $`h(\xi )`$ in $`b`$ and $`\xi `$ only if $`h(\xi )=1`$ for all $`\xi `$.

Reflection symmetry: Once $`h(\xi )=1`$ it follows that the denominator in the three-point function is invariant under the reflection $`\alpha _Iq\alpha _I`$ for each $`I`$ and thus the three-point function is covariant with respect to reflection symmetry in the sense that

$$𝒢_3(q\alpha _1,\alpha _2,\alpha _3)=R(q2\alpha _1)𝒢_3(\alpha _1,\alpha _2,\alpha _3)$$
$`(43)`$
where the prefactor depends only on the reflected parameter $`\alpha _1`$. Thus in the two-exponential theory, reflection covariance is an output rather than an input.

It is interesting to note how this reflection covariance expresses itself in terms of the two-point function defined. If we interpret the two-point function normalised by the volume factor $`4\pi `$, as the inner product of primary states

$$<\alpha _1,\alpha _2>\underset{x0}{lim}\frac{1}{4\pi }x^{2\mathrm{\Delta }_\alpha }𝒢_2(q\alpha _1,\alpha _2;x)$$
$`(44)`$
we have

$$\left(\begin{array}{cc}<\alpha ,\alpha >& <\alpha ,q\alpha >\\ <q\alpha ,\alpha >& <q\alpha ,q\alpha >\end{array}\right)=\left(\begin{array}{cc}1& R^1(q2\alpha )\\ R(q2\alpha )& 1\end{array}\right)\delta (0)$$
$`(45)`$
It is clear that the matrix in (45) is hermitian and has zero determinant. Hence one linear combination of the states, namely $`|\alpha >R(q2\alpha )|q\alpha >`$, has zero norm and decouples. Thus effectively,

$$|\alpha >=R(q2\alpha )|q\alpha >$$
$`(46)`$
which means that there is actually only one physical state for each conformal weight $`\mathrm{\Delta }_\alpha `$. This may seem surprising but from (43) it is seen to be a manifestation of the reflection covariance.

Comparison of the one and two-exponential path integrals: In order to compare the one and two exponential theories, we begin by recalling that $`𝐂_3`$ in both the theories is defined in terms of the correlation functions of vertex operators in a free field theory. In the two-exponential theory, the relevant integral is given by

$$𝑑\varphi e^{S(\varphi )}U_b^{iu}U_c^{iv}𝚷_3=\mathrm{\Gamma }(1+iu)\mathrm{\Gamma }(1+iv)\frac{K(\xi ,\alpha _I)}{k^{}(\xi )}f[\xi ,i(u+v)]$$
$`(47)`$
where $`S(\varphi )`$ is the Action for the free theory. The corresponding equation, in the one-exponential theory, is obtained by letting $`v0`$ and takes the form

$$𝑑\varphi e^{S(\varphi )}U_b^{iu}𝚷_3=\mathrm{\Gamma }(1+iu)\frac{K(\xi ,\alpha _I)}{k^{}(\xi )}\times f[\xi ,iu]$$
$`(48)`$
The integral corresponding to the zero-mode integral in (29) in the one-exponential theory for an a priori arbitrary $`f[\xi ,iu]`$ can be performed to yield $`f[\xi ,c\xi ]/\text{tanh}\pi c\xi `$. If we assume that the one-exponential theory leads to the DOZZ Ansatz, then the function $`f[\xi ,iu]`$ is determined uniquely to be $`b\text{tan}(\pi u)k^{}(ibu)/k(ibu)`$ for $`\xi =ibu`$. We may now ask, what choice, if any, for the function $`f[\xi ,i(u+v)]`$ in the two-exponential theory will produce the DOZZ Ansatz. It is easy to see that if we choose

$$f[\xi ,i(u+v)]=\frac{k^{}(\xi )}{k(\xi )}[b\text{tanh}\pi u+c\text{tanh}\pi v]$$
$`(49)`$
the $`Z_0`$ integral in (29) produces $`k^{}/k`$ in accordance with the DOZZ Ansatz. This choice has the virtue that it reduces to the one-exponential result in the limit $`v0`$. However, one can easily convince oneself that there exist other choices for the function $`f[\xi ,i(u+v)]`$. If we choose, for example, $`f[\xi ,i(u+v)]=\text{sgn}[i(u+v)]`$ (the sgn-factor is necessary to preserve the symmetry of the path integral under a change of sign of $`b,c`$ and the $`\alpha `$’s) and integrate over $`u`$ and $`v`$ first, we obtain,

$$Z_0(\xi )=_0^{\mathrm{}}𝑑t\frac{\text{sinh}((q2\xi )t)}{\text{sinh}(bt)\text{sinh}(ct)},t|\varphi _0|$$
$`(50)`$
Up to a regulating term,$`^4`$ The regulating term can be obtained naturally by modifying the path-integral measure $`d\varphi _0`$ to $`lim_{ϵ0}d\varphi _0\left[1\text{exp}(|\varphi _0|/ϵ)\right]`$. this will be recognized, from the definition of $`k(\xi )`$ in (26), as $`k^{}(\xi )/k(\xi )`$. Since the choice for the extrapolation function is not unique we conclude that any result that is based on an extrapolation is not conclusive. It is therefore desirable to have a direct proof of the DOZZ result for the three-point function without any reference to the extrapolation. This is exactly what was achieved in this paper by showing that the sum rule in the two-exponential theory has a unique (constant) solution. In contrast, in the one-exponential theory, if the specific extrapolation leading to the DOZZ Ansatz is not made, the final result can only be obtained up to a periodic function of $`\xi `$.

Acknowledgements: We thank I. Sachs, J. Teschner and P. Watts for useful discussions.

References

1. H. Poincare, J. Math. Pure App. 5 se 4 (1898) 157; N. Seiberg, Notes on Quantum Liouville Theory and Quantum Gravity, in ”Random Surfaces and Quantum Gravity”, ed. O. Alvarez, E. Marinari, and P. Windey, Plenum Press, 1990; A. Polyakov, Phys. Lett. B103 (1981) 207; T. Curtright and C. Thorn, Phys. Rev. Lett. 48 (1982) 1309; E. Braaten, T. Curtright and C. Thorn, Phys. Lett 118 (1982) 115; Ann. Phys. 147 (1983) 365; J.-L. Gervais and A. Neveu, Nuc. Phys. B238 (1984) 125; B238 (1984) 396; 257\[FS14\] (1985) 59; E. D’Hoker and R. Jackiw, Phys. Rev. D26 (1982) 3517.

2. H. Dorn and H.-J. Otto, Phys. Lett. B291 (1992) 39; Nuc. Phys. B429 (1994) 375.

3. A. and Al. Zamolodchikov, Nuc. Phys. B 477 (1996) 577; Low Dimensional Applications of Quantum Field Theory, ed. L. Baulieu, V. Kazakov, M. Picco, and P. Windey, NATO ASI Series B, Physics Volume 362, Plenum Press 1997.

4. J. Teschner, Phys. Lett. B363 (1995) 65.

5. L. O’Raifeartaigh, J. M. Pawlowski, and V. V. Sreedhar, Annals of Physics 277 (1999), 117.

6. See also Vl. S. Dotsenko, Mod. Phys. Lett. A6 (1991), 3601; E. D’Hoker, Mod. Phys. Lett A6 (1991), 745; P. Mansfield, Phys. Lett. B242 (1990), 242.

7. G. N. Watson, Proc. Soc. London 95 (1918), 83; A. Sommerfeld, ”Partial Differential Equations of Physics”, Academic Press, New York 1949.

8. Vl. Dotsenko and V. Fateev, Nuc. Phys. B251 (1985) 691.


text:
# Deriving the Nonlinear Cosmological Power Spectrum and Bispectrum from Analytic Dark Matter Halo Profiles and Mass Functions

## 1 Introduction

Two conceptual pictures of galaxy clustering have been examined in the literature, the continuous hierarchical clustering model and the power-law cluster model (Peebles 1980, §61). In the hierarchical clustering model, which has emerged as the accepted model over the past two decades, galaxy clustering is characterized by power-law correlation functions: the $`N`$-point correlation function $`\xi _N`$ scales with configuration size as $`\xi _Nr^{\gamma _N}\xi _2^{(N1)}`$, where $`\gamma _N=(N1)\gamma `$ and the two-point correlation function goes as $`\xi _2=\xi r^\gamma `$. The hierarchical model is motivated by the observed power-law behavior $`\gamma 1.8`$ of galaxy correlations (Groth & Peebles 1977; Fry & Peebles 1978), with a theoretical basis in a self-similar, scale-invariant solution to the equations of motion (Davis & Peebles 1977).

The alternative power-law cluster model has an even longer history (Neyman & Scott 1952; Peebles 1974, 1980; McClelland & Silk 1977; Scherrer & Bertschinger 1991; Sheth & Jain 1997; Valageas 1998; Yano & Gouda 1999). In this model, galaxies are placed in spherical clumps that are assumed to follow a power-law density profile $`\rho (r)r^ϵ`$, with the centers of the clumps distributed randomly. The resulting two-point correlation function is also a power law with a logarithmic slope $`\gamma =2ϵ3`$. While it is possible to reproduce the observed two-point function by an appropriate choice of the power index $`ϵ=(3+\gamma )/22.4`$, Peebles and Groth (1975) pointed out that this model produces a three-point function that is too steep to be consistent with observations in the Zwicky and Lick catalogs.

In an earlier paper (Ma & Fry 2000a), we have shown that in the nonlinear regime, the three-point correlation function $`\zeta =\xi _3`$ of the cosmological mass density field does not exactly follow the prediction $`\zeta \xi ^2`$ of the hierarchical clustering model. These conclusions are drawn from study of high resolution numerical simulations of a cold dark matter (CDM) model with cosmological constant and of a model with scale-free initial conditions $`P(k)k^n`$ with $`n=2`$. In experiments replacing simulation dark matter halos with power-law density profiles, $`\rho (r)r^ϵ`$, we have demonstrated that the behavior of the correlation functions in the nonlinear regime are determined by the halo profiles, but that it is not possible to match both the two- and three-point correlations with a single slope $`ϵ`$. These results differ from the predictions of both of these two conceptual models.

In this paper, we expand our previous study of the nonlinear two- and three-point correlation functions by investigating a new prescription that takes into account the non-power-law profiles of halos, the distribution of halo masses, and the spatial correlations of halo centers. Each of these ingredients has been well studied in the literature. We find that this halo model provides a good description of the two- and three-point correlation functions in both the $`n=2`$ and CDM simulations over the entire range of scales from the weak clustering, perturbative regime on large length scales, to the strongly nonlinear regime on small length scales. Our result is approximately hierarchical over an intermediate range of scales, thus uniting the two pictures. An independent recent study by Seljak (2000), which appeared during completion of this work, has also examined the two-point power spectrum in a similar construction and has found that this type of approach can reproduce the power spectrum in the CDM model. The analytic model proposed here can be used to compute the two- and three-point correlation functions and their Fourier transforms, the power spectrum and bispectrum, over any range of scale where the input halo properties are valid.

In a subsequent paper (Ma & Fry 2000c), we study the predictions of this analytic halo model for the asymptotic nonlinear behavior of the $`N`$-point correlation functions and the pairwise velocities and examine the conditions required for stable clustering.

The outline of this paper is as follows. In §2 we describe the three input ingredients of the model: halo density profiles, halo mass functions, and halo-halo correlations. In §3 we assemble these ingredients and construct analytic expressions for the two-point correlation function $`\xi (r)`$ and the power spectrum $`P(k)`$. In §4 we do the same for the three-point correlation function $`\zeta (r_1,r_2,r_3)`$ and its Fourier transform, the bispectrum $`B(k_1,k_2,k_3)`$. In §5 we test the validity of this new model by comparing its predictions with results from numerical simulations of an $`n=2`$ scale free model and a low-density CDM model with a cosmological constant ($`\mathrm{\Lambda }`$CDM). We also present results of the synthetic halo replacement technique used to enhance the numerical resolution. In §6 we discuss further the physical meanings and implications of the model. In particular, we elaborate on two important implications of this model: deviations from the common assumptions of stable clustering and hierarchical clustering. Section 7 is a summary.

## 2 Model Ingredients

### 2.1 Halo Mass Density Profile

It has been suggested recently that the mass density profiles of cold dark matter halos have a roughly universal shape, generally independent of cosmological parameters (Navarro, Frenk, & White 1996, 1997)

$$\frac{\rho (r)}{\overline{\rho }}=\overline{\delta }u(r/R_s),$$
(1)
where $`\overline{\delta }`$ is a dimensionless density amplitude, $`R_s`$ is a characteristic radius, and $`\overline{\rho }`$ is the mean background density. We consider two functional forms for the density profiles

$`u_I(x)`$ $`=`$ $`{\displaystyle \frac{1}{x^p(1+x)^{3p}}},`$
$`u_{II}(x)`$ $`=`$ $`{\displaystyle \frac{1}{x^p(1+x^{3p})}}.`$ (2)
Both forms have asymptotic behaviors $`x^p`$ at small $`x`$ and $`x^3`$ at large $`x`$, but they differ in the transition region. The first form $`u_I(x)`$ with $`p=1`$ is found to provide a good fit to simulation halos by Navarro et al. (1996, 1997), whereas the second form $`u_{II}(x)`$ with a steeper inner slope $`p=3/2`$ is favored by Moore et al. (1998, 1999). Some independent simulations have produced halos that are well fit by the shallower $`p=1`$ inner slope (e.g., Hernquist 1990; Dubinski & Carlberg 1991; Huss, Jain, & Steinmetz 1999), and others the steeper $`p>1`$ slope (e.g., Fukushige and Makino 1997). Jing & Suto (2000) have recently reported a mass-dependent inner slope, with $`p1.5`$ for galactic-mass halos and $`p1`$ for cluster-mass halos. Many of these authors find that the outer profile scales as $`r^3`$, but steeper outer profiles have also been suggested (Hernquist 1990; Dubinski & Carlberg 1991). Given these uncertainties, we will consider in this paper both types of profiles in equation (2).

The parameters $`R_s`$ and $`\overline{\delta }`$ in equation (1) are generally functions of the halo mass $`M`$. A concentration parameter,

$$c=\frac{R_{200}}{R_s},$$
(3)
can be used to quantify the central density of a halo (Navarro et al. 1997), where $`R_{200}`$ is the radius within which the average density is 200 times the mean density of the universe. Using $`M=800\pi \overline{\rho }R_{200}^3/3`$, we can relate $`R_s`$ and $`\overline{\delta }`$ to $`M`$ and $`c`$, where the scale radius $`R_s`$ is

$$R_s=\frac{1}{c}\left(\frac{3M}{800\pi \overline{\rho }}\right)^{1/3}=\frac{1.63\times 10^5}{\mathrm{\Omega }_m^{1/3}c}\left(\frac{M}{h^1M_{}}\right)^{1/3}h^1\mathrm{Mpc},$$
(4)
and the density amplitude $`\overline{\delta }`$ is

$`\overline{\delta }_I`$ $`=`$ $`{\displaystyle \frac{200c^3}{3[\mathrm{ln}(1+c)c/(1+c)]}},p=1,`$
$`\overline{\delta }_{II}`$ $`=`$ $`{\displaystyle \frac{100c^3}{\mathrm{ln}(1+c^{3/2})}},p={\displaystyle \frac{3}{2}}.`$ (5)
Typical values of $`c`$ are in the range of a few to ten for type I and perhaps a factor of three smaller for type II. There is a weak dependence on mass, such that less massive halos have a larger central density (e.g., Cole & Lacey 1996; Tormen, Bouchet, & White 1997; Navarro et al. 1996, 1997; Jing & Suto 2000). This is understood in general terms as reflecting the mean density at the redshift $`z_f`$ when the halo initially collapsed, $`\overline{\delta }(1+z_f)^3`$. For $`\mathrm{\Omega }=1`$ this is $`c\sigma (M)`$, or $`cM^{(3+n)/6}`$ in a scale-free model.

### 2.2 Halo Mass Function

The number density of halos with mass $`M`$ within a logarithmic interval is often approximated by the prescription of Press & Schechter (1974),

$$\frac{dn}{d\mathrm{ln}M}=\sqrt{\frac{2}{\pi }}\frac{d\mathrm{ln}\sigma ^1}{d\mathrm{ln}M}\frac{\overline{\rho }}{M}\nu e^{\nu ^2/2},\nu =\frac{\delta _c}{\sigma (M)},$$
(6)
where $`\delta _c`$ is a parameter characterizing the linear overdensity at the onset of gravitational collapse, and $`\sigma `$ is the linear rms mass fluctuations in spheres of radius $`R`$

$$\sigma ^2(M)=_0^{\mathrm{}}\frac{4\pi k^2dk}{(2\pi )^3}P(k)W^2(kR),$$
(7)
where $`W(x)=3(\mathrm{sin}xx\mathrm{cos}x)/x^3`$ is the Fourier transform of a real-space tophat window function. The mass $`M`$ is related to $`R`$ by $`M=4\pi \overline{\rho }R^3/3`$. For scale free models with a power law initial power spectrum $`Pk^n`$, this is $`\sigma =(M/M_{})^{(3+n)/6}`$. The parameter $`M_{}`$ characterizes the mass scale at the onset of nonlinearity, $`\sigma (M_{})=1`$, and is related to the nonlinear wavenumber $`k_{\mathrm{nl}}`$ (defined as \]$`_0^{k_{\mathrm{nl}}}4\pi k^2𝑑kP(k)/(2\pi )^3=1`$) by

$$M_{}=\frac{4\pi \overline{\rho }}{3}\frac{B(n)}{k_{\mathrm{nl}}^3}=\frac{4\pi \overline{\rho }}{3}R_{}^3,$$
(8)
where

$`B(n)`$ $`=`$ $`(k_{\mathrm{nl}}R_{})^3=\left[(n+3){\displaystyle _0^{\mathrm{}}}𝑑xx^{n+2}W^2(x)\right]^{3/(n+3)},`$
$`B^{(3+n)/3}`$ $`=`$ $`\mathrm{sin}\left[(n+2){\displaystyle \frac{\pi }{2}}\right]\mathrm{\Gamma }(n+2){\displaystyle \frac{9(2^n)(3+n)}{(n)(1n)(3n)}}`$ (9)
(defined for $`3n<1`$). Various modifications to the Press-Schechter mass function have been suggested (e.g., Sheth & Tormen 1999; Lee & Shandarin 1999; Jenkins et al. 2000) to improve the accuracy of the original formula.

### 2.3 Halo-Halo Correlations

Dark matter halos do not cluster in the same way as the mass density field. On large scales, a bias parameter $`b`$ is typically used to quantify this difference. Let $`\xi _{\mathrm{halo}}(r;M,M^{})`$ be the two-point correlation function of halos with masses $`M`$ and $`M^{}`$, $`\xi _{\mathrm{lin}}(r)`$ be the linear correlation function for the mass density field, and $`P_{\mathrm{halo}}`$ and $`P_{\mathrm{lin}}`$ be the corresponding power spectra. On large length scales, we assume a linear bias and write

$`\xi _{\mathrm{halo}}(r;M,M^{})`$ $`=`$ $`b(M)b(M^{})\xi _{\mathrm{lin}}(r),`$
$`P_{\mathrm{halo}}(k;M,M^{})`$ $`=`$ $`b(M)b(M^{})P_{\mathrm{lin}}(k).`$ (10)
Based on the peak and the Press-Schechter formalism, Mo & White (1996) developed a model for the linear bias $`b(M)`$, which is later modified by Jing (1998) to be

$$b(M)=\left(1+\frac{\nu ^21}{\delta _c}\right)\left(\frac{1}{2\nu ^4}+1\right)^{0.060.02n},\nu =\frac{\delta _c}{\sigma (M)}.$$
(11)
The original formula for $`b(M)`$ by Mo & White includes only the first factor above; the second factor, dependent on the primordial spectral index $`n`$, is obtained empirically for an improved fit to simulation results at the lower mass end (Jing 1998). In this bias model, $`b(M)`$ is below unity for $`MM_{}`$ (where $`\sigma (M_{})=1`$) and reaches $`0.5`$ for $`M0.01M_{}`$. Small dark matter halos are therefore anti-biased relative to the mass density. For $`MM_{}`$, $`b(M)`$ increases monotonically with the halo mass and reaches $`b10`$ at $`M100M_{}`$. Nonlinear effects on the bias have been studied (Kravtsov & Klypin 1999 and references therein), but they are unimportant in our model because the halo-halo correlation terms contribute significantly only on large length scales in the linear regime (see §3 and 4).

Similarly, we use higher order bias parameters to relate the higher-order correlation functions for halos and mass density. In this paper we examine the three-point correlation function $`\zeta (r_1,r_2,r_3)`$ and its Fourier transform, the bispectrum $`B(k_1,k_2,k_3)`$ (see §4 for a more detailed discussion). On large length scales where the amplitude of $`\delta `$ is small, perturbation theory can be used to relate the lowest order contribution to the bispectrum of the mass density to the linear power spectrum $`P_{\mathrm{lin}}`$ (Fry 1984):

$`B^{(0)}(k_1,k_2,k_3)`$ $`=`$ $`F_{12}P_{\mathrm{lin}}(k_1)P_{\mathrm{lin}}(k_2)+F_{23}P_{\mathrm{lin}}(k_2)P_{\mathrm{lin}}(k_3)+F_{31}P_{\mathrm{lin}}(k_3)P_{\mathrm{lin}}(k_1),`$
$`F_{ij}`$ $`=`$ $`{\displaystyle \frac{10}{7}}+(k_i/k_j+k_j/k_i)(\widehat{\text{k}}_i\widehat{\text{k}}_j)+{\displaystyle \frac{4}{7}}(\widehat{\text{k}}_i\widehat{\text{k}}_j)^2.`$ (12)
Using this perturbative result and the results of Mo, Jing, & White (1997), we can write the halo bispectrum as

$`B_{\mathrm{halo}}(k_1,k_2,k_3;M,M^{},M^{\prime \prime })`$ $`=`$ $`\left[b(M)b(M^{})b(M^{\prime \prime })F_{12}+b(M)b(M^{})b_2(M^{\prime \prime })\right]P_{\mathrm{lin}}(k_1)P_{\mathrm{lin}}(k_2)`$
$`+`$ $`\left[b(M)b(M^{})b(M^{\prime \prime })F_{23}+b(M)b_2(M^{})b(M^{\prime \prime })\right]P_{\mathrm{lin}}(k_2)P_{\mathrm{lin}}(k_3)`$
$`+`$ $`\left[b(M)b(M^{})b(M^{\prime \prime })F_{31}+b_2(M)b(M^{})b(M^{\prime \prime })\right]P_{\mathrm{lin}}(k_3)P_{\mathrm{lin}}(k_1),`$
where $`b(M)`$ is given by equation (11), and the quadratic bias parameter $`b_2(M)`$ is

$$b_2(M)=\frac{8}{21}\frac{(\nu ^21)}{\delta _c}+\left(\frac{\nu }{\delta _c}\right)^2(\nu ^23).$$
(14)
For the special equilateral case of $`k_1=k_2=k_3=k`$, equation (LABEL:Bbias) simplifies to

$`B_{\mathrm{halo}}^{\mathrm{eq}}(k;M,M^{},M^{\prime \prime })`$ $`=`$ $`[{\displaystyle \frac{12}{7}}b(M)b(M^{})b(M^{\prime \prime })+b(M)b(M^{})b_2(M^{\prime \prime })`$ (15)
$`+`$ $`b(M)b_2(M^{})b(M^{\prime \prime })+b_2(M)b(M^{})b(M^{\prime \prime })]P_{\mathrm{lin}}^2(k).`$
In practice, the terms involving $`b_2(M)`$ in equations (LABEL:Bbias) and (15) make only a small net contribution. For simplicity, we will therefore not include this term in the subsequent derivations and calculations.

## 3 Two-Point Statistics: $`\xi (r)`$ and $`P(k)`$

We now construct our analytic halo model for the two-point correlation function $`\xi (r)`$ and the power spectrum $`P(k)`$. The two-point correlation function of the cosmological mass density field $`\delta =\delta \rho /\overline{\rho }`$ is

$$\xi (\text{r})=\delta (\text{x})\delta (\text{x}+\text{r}).$$
(16)
The Fourier transform of $`\xi (r)`$ is the mass power spectrum $`P(k)=d^3re^{i\text{k}\text{r}}\xi (r)`$, which is related to the density field in $`k`$-space by $`\delta (\text{k}_1)\delta (\text{k}_2)=P(k_1)(2\pi )^3\delta _D(\text{k}_1+\text{k}_2)`$, where $`\delta _D`$ is the Dirac delta-function.

The two-point correlation function measures the excess probability above the Poisson distribution of finding a pair of objects with separation $`r`$ (Peebles 1980). The objects can be taken to be dark matter particles, most of which cluster gravitationally in the form of dark matter halos. One should therefore be able to express $`\xi `$ for the density field in terms of properties of dark matter halos. In this picture, we can write the contributions to $`\xi `$ as two separate terms, one from particle pairs in the same halo, and the other from pairs that reside in two different halos. In realistic cosmological models, dark matter halos exhibit a spectrum of masses that can be characterized by a distribution function $`dn/dM`$, and the halo centers are spatially correlated. Taking these factors into consideration, we can write the two-point correlation function for $`\delta `$ in terms of the halo density profile $`u(x)`$, halo mass function $`dn/dM`$, and halo-halo correlation function $`\xi _{\mathrm{halo}}`$ discussed in §2. We write

$$\xi (r)=\xi _{1h}(r)+\xi _{2h}(r),$$
(17)
where the subscripts “$`1h`$” and “$`2h`$” denote contributions from particle pairs in “1-halo” and “2-halos”, respectively, and

$`\xi _{1h}(r)`$ $`=`$ $`{\displaystyle d^3r^{}𝑑M\frac{dn}{dM}\overline{\delta }^2u(r^{}/R_s)u(|\text{r}^{}+\text{r}|/R_s)}`$
$`\xi _{2h}(r)`$ $`=`$ $`{\displaystyle d^3r^{}d^3r^{\prime \prime }𝑑M^{}\frac{dn}{dM^{}}\overline{\delta }^{}u(r^{}/R_s^{})𝑑M^{\prime \prime }\frac{dn}{dM^{\prime \prime }}\overline{\delta }^{\prime \prime }u(r^{\prime \prime }/R_s^{\prime \prime })\xi _{\mathrm{halo}}(|\text{r}^{}\text{r}^{\prime \prime }+\text{r}|)}`$ (18)
$`=`$ $`{\displaystyle d^3r^{}d^3r^{\prime \prime }\left[𝑑M\frac{dn}{dM}\overline{\delta }u(r^{}/R_s)b(M)\right]\left[𝑑M\frac{dn}{dM}\overline{\delta }u(r^{\prime \prime }/R_s)b(M)\right]}`$
$`\times \xi _{\mathrm{lin}}(|\text{r}^{}\text{r}^{\prime \prime }+\text{r}|).`$
These expressions arise from averaging over displacements $`\text{r}^{}`$, $`\text{r}^{\prime \prime }`$ of halo centers from the particle positions $`\text{r}_1`$, $`\text{r}_2`$, where $`r=|\text{r}_1\text{r}_2|`$. In the last expression above, we have used the bias model of equation (10) to relate the halo-halo correlation function $`\xi _{\mathrm{halo}}`$ to the linear correlation function $`\xi _{\mathrm{lin}}`$ of the mass density field.

As we will show in §5, the dominant contribution to the two-point correlation function in the nonlinear regime on small length scales is from the first, 1-halo term $`\xi _{1h}`$ for particle pairs that reside in the same halos. This makes intuitive sense, because closely spaced particle pairs are most likely to be found in the same halo. This term is determined by the convolution of the dimensionless density profile with itself,

$$\lambda (x)=d^3yu(y)u(|\text{x}+\text{y}|).$$
(19)
For many forms of $`u(x)`$, the angular integration in this equation is analytic, and $`\lambda `$ can be reduced to a simple one-dimensional integral over $`y`$. For some special cases, $`\lambda `$ can even be reduced to an analytic expression. We leave the detailed results for $`\lambda `$ to the Appendix.

In $`k`$-space, the convolutions in equation (18) for $`\xi (r)`$ become simple products. Using $`\stackrel{~}{u}(q)`$ to denote the Fourier transform of $`u(x)`$, where $`\stackrel{~}{u}(q)=d^3xu(x)e^{i\text{q}\text{x}}`$, we can readily transform equation (18) into expressions for the mass power spectrum:

$$P(k)=P_{1h}(k)+P_{2h}(k),$$
(20)
where the 1-halo and 2-halo terms are

$`P_{1h}(k)`$ $`=`$ $`{\displaystyle 𝑑M\frac{dn}{dM}[R_s^3\overline{\delta }\stackrel{~}{u}(kR_s)]^2}`$
$`P_{2h}(k)`$ $`=`$ $`{\displaystyle 𝑑M\frac{dn}{dM}R_s^3\overline{\delta }\stackrel{~}{u}(kR_s)𝑑M^{}\frac{dn}{dM^{}}R_s^3\overline{\delta }^{}\stackrel{~}{u}(kR_s^{})P_{\mathrm{halo}}(k)}`$
$`=`$ $`\left[{\displaystyle 𝑑M\frac{dn}{dM}R_s^3\overline{\delta }\stackrel{~}{u}(kR_s)b(M)}\right]^2P_{\mathrm{lin}}(k).`$
To arrive at the last expression above, we have again used the bias model of equation (10). For computational efficiency, we find that the algebraic expressions

$`\stackrel{~}{u}_I(q)`$ $`=`$ $`{\displaystyle \frac{4\pi \{\mathrm{ln}(e+1/q)\mathrm{ln}[\mathrm{ln}(e+1/q)]/3\}}{(1+q^{1.1})^{(2/1.1)}}},p=1`$
$`\stackrel{~}{u}_{II}(q)`$ $`=`$ $`{\displaystyle \frac{4\pi \{\mathrm{ln}(e+1/q)+0.25\mathrm{ln}[\mathrm{ln}(e+1/q)]\}}{1+0.8q^{1.5}}},p={\displaystyle \frac{3}{2}}`$ (22)
provide excellent fits for the profiles of Navarro et al. (1997) and Moore et al. (1999), with less than 4% rms error for form I and less than 1% rms error for form II. The functional form is chosen to reproduce the asymptotic behaviors: $`\stackrel{~}{u}4\pi \mathrm{ln}q`$ at small $`q`$ (with no radial cutoff), and $`\stackrel{~}{u}q^2`$ (type I) and $`\stackrel{~}{u}q^{3/2}`$ (type II) at large $`q`$.

The two-point $`\xi (r)`$ and $`P(k)`$ can now be computed analytically from equations (18) and (3). The inputs are equation (2) or (22) for the halo density profile $`u(x)`$ or $`\stackrel{~}{u}(q)`$, equations (4) and (5) for $`R_s`$ and $`\overline{\delta }`$, equation (6) for the halo mass function $`dn/dM`$, and equation (10) for the halo-halo correlation function. Since the halo density profile appears to have a nearly universal form regardless of background cosmology, $`\xi (r)`$ and $`P(k)`$ depend on cosmological parameters mainly through $`\sigma (M)`$ of equation (7) and the halo concentration $`c(M)`$ or central density $`\overline{\delta }(M)`$. (See Ma & Fry 2000c for a more detailed discussion of $`c(M)`$.)

## 4 Three-Point Statistics: $`\zeta `$ and $`B`$

Here we construct our analytic halo model for the three-point correlation function $`\zeta `$ and the bispectrum $`B`$. The joint probability of finding three objects in volume elements $`dV_1,dV_2`$, and $`dV_3`$ is given by

$$dP=[1+\xi (r_1)+\xi (r_2)+\xi (r_3)+\zeta (r_1,r_2,r_3)]\overline{n}^3dV_1dV_2dV_3,$$
(23)
where $`\xi (r)`$ and $`\zeta (r_1,r_2,r_3)`$ are the two- and three-point correlation functions, respectively, $`\overline{n}`$ is the mean number density of objects, and $`r_1,r_2`$ and $`r_3`$ are the lengths of the sides of the triangle defined by the three objects (Peebles 1980). The Fourier transform of the three-point correlation function $`\zeta (r_1,r_2,r_3)`$ is the bispectrum $`B(k_1,k_2,k_3)`$, which is related to the density field in $`k`$-space by $`\delta (\text{k}_1)\delta (\text{k}_2)\delta (\text{k}_3)=B(k_1,k_2,k_3)(2\pi )^3\delta _D(\text{k}_1+\text{k}_2+\text{k}_3)`$. The bispectrum depends on any three parameters that define a triangle in $`k`$-space. A particular simple configuration to study is the equilateral triangle ($`k_1=k_2=k_3=k`$), and in this case the bispectrum $`B^{\mathrm{eq}}`$ depends only on a single wavenumber.

Similar to the two-point halo model of §3, we can write the contributions to the three-point correlation function $`\zeta `$ of the mass density as three separate terms, each term representing particle triplets that reside in a single halo, two distinct halos, or three distinct halos. Taking into account the halo mass distribution and halo-halo correlations discussed in §2, we obtain

$$\zeta (r_1,r_2,r_3)=\zeta _{1h}(r_1,r_2,r_3)+\zeta _{2h}(r_1,r_2,r_3)+\zeta _{3h}(r_1,r_2,r_3),$$
(24)
where the separate 1-halo, 2-halo, and 3-halo terms are

$`\zeta _{1h}(r_1,r_2,r_3)`$ $`=`$ $`{\displaystyle d^3r𝑑M\frac{dn}{dM}\overline{\delta }^3u(r/R_s)u(|\text{r}+\text{r}_1\text{r}_2|/R_s)u(|\text{r}+\text{r}_1\text{r}_3|/R_s)}`$
$`\zeta _{2h}(r_1,r_2,r_3)`$ $`=`$ $`{\displaystyle d^3rd^3r^{}𝑑M\frac{dn}{dM}\overline{\delta }^2u(r/R_s)u(|\text{r}+\text{r}_1\text{r}_2|/R_s)𝑑M^{}\frac{dn}{dM^{}}\overline{\delta }^{}u(r^{}/R_s^{})}`$ (25)
$`\times \xi _{\mathrm{halo}}(|\text{r}\text{r}^{}+\text{r}_1\text{r}_3|)+\text{sym.(1,2,3)}`$
$`\zeta _{3h}(r_1,r_2,r_3)`$ $`=`$ $`{\displaystyle d^3rd^3r^{}d^3r^{\prime \prime }𝑑M\frac{dn}{dM}\overline{\delta }u(r/R_s)𝑑M^{}\frac{dn}{dM^{}}\overline{\delta }^{}u(r^{}/R_s^{})}`$
$`\times {\displaystyle }dM^{\prime \prime }{\displaystyle \frac{dn}{dM^{\prime \prime }}}\overline{\delta }^{\prime \prime }u(r^{\prime \prime }/R_s^{\prime \prime })\zeta _{\mathrm{halo}}(\text{r}+\text{r}_1,\text{r}^{}+\text{r}_2,\text{r}^{\prime \prime }+\text{r}_3).`$
The dominant contribution to the three-point correlation function in the nonlinear regime is from the first term $`\zeta _{1h}`$, which comes from particle triplets that reside in the same halo. This term is determined by the convolution $`\gamma (\text{x}_1,\text{x}_2)=d^3yu(y)u(|\text{y}+\text{x}_1|)u(|\text{y}+\text{x}_2|)`$ of three factors of the density profile $`u(x)`$, and is analogous to the convolution $`\lambda `$ in equation (19) for the one-halo term $`\xi _{1h}`$ in the two-point correlation function.

The bispectrum of the mass density field $`\delta `$ in $`k`$-space can be obtained by Fourier transforming the equations above. We find

$$B(k_1,k_2,k_3)=B_{1h}(k_1,k_2,k_3)+B_{2h}(k_1,k_2,k_3)+B_{3h}(k_1,k_2,k_3),$$
(26)
where

$`B_{1h}(k_1,k_2,k_3)`$ $`=`$ $`{\displaystyle 𝑑M\frac{dn}{dM}[R_s^3\overline{\delta }\stackrel{~}{u}(k_1R_s)][R_s^3\overline{\delta }\stackrel{~}{u}(k_2R_s)][R_s^3\overline{\delta }\stackrel{~}{u}(k_3R_s)]}`$
$`B_{2h}(k_1,k_2,k_3)`$ $`=`$ $`{\displaystyle 𝑑M\frac{dn}{dM}[R_s^3\overline{\delta }\stackrel{~}{u}(k_1R_s)][R_s^3\overline{\delta }\stackrel{~}{u}(k_2R_s)]}`$ (27)
$`\times {\displaystyle }dM^{}{\displaystyle \frac{dn}{dM^{}}}R_s^3\overline{\delta }^{}\stackrel{~}{u}(k_3R_s^{})P_{\mathrm{halo}}(k_3;M,M^{})+\text{sym.(1,2,3)}`$
$`B_{3h}(k_1,k_2,k_3)`$ $`=`$ $`{\displaystyle 𝑑M\frac{dn}{dM}R_s^3\overline{\delta }\stackrel{~}{u}(k_1R_s)𝑑M^{}\frac{dn}{dM^{}}R_s^3\overline{\delta }^{}\stackrel{~}{u}(k_2R_s^{})}`$
$`\times {\displaystyle }dM^{\prime \prime }{\displaystyle \frac{dn}{dM^{\prime \prime }}}R_s^{\prime \prime 3}\overline{\delta }^{\prime \prime }\stackrel{~}{u}(k_3R_s^{\prime \prime })B_{\mathrm{halo}}(k_1,k_2,k_3;M,M^{},M^{\prime \prime }).`$
The halo-halo power spectrum $`P_{\mathrm{halo}}(k)`$ and bispectrum $`B_{\mathrm{halo}}(k_1,k_2,k_3)`$ are related to the linear mass power spectrum $`P_{\mathrm{lin}}(k)`$ by equations (10) and (LABEL:Bbias).

The expressions for the mass bispectrum above simplify considerably for the equilateral triangle configuration, and

$$B^{\mathrm{eq}}(k)=B_{1h}^{\mathrm{eq}}(k)+B_{2h}^{\mathrm{eq}}(k)+B_{3h}^{\mathrm{eq}}(k),$$
(28)
where

$`B_{1h}^{\mathrm{eq}}(k)`$ $`=`$ $`{\displaystyle 𝑑M\frac{dn}{dM}[R_s^3\overline{\delta }\stackrel{~}{u}(kR_s)]^3}`$
$`B_{2h}^{\mathrm{eq}}(k)`$ $`=`$ $`3\left[{\displaystyle 𝑑M\frac{dn}{dM}[R_s^3\overline{\delta }\stackrel{~}{u}(kR_s)]^2b(M)}\right]\left[{\displaystyle 𝑑M\frac{dn}{dM}R_s^3\overline{\delta }\stackrel{~}{u}(kR_s)b(M)}\right]P_{\mathrm{lin}}(k)`$ (29)
$`B_{3h}^{\mathrm{eq}}(k)`$ $`=`$ $`\left[{\displaystyle 𝑑M\frac{dn}{dM}R_s^3\overline{\delta }\stackrel{~}{u}(kR_s)b(M)}\right]^3{\displaystyle \frac{12}{7}}P_{\mathrm{lin}}^2(k).`$
Here we have written out explicitly the bias factors $`b(M)`$ using equations (10) and (15), and we have neglected terms with $`b_2(M)`$ as discussed in §2.3.

## 5 $`N`$-body Experiments and Numerical Results

In this section we compare the predictions of our analytical model described in §2, 3, and 4 with results from cosmological $`N`$-body simulations. We examine two cosmological models: an $`n=2`$ scale-free model and a low-density $`\mathrm{\Lambda }`$CDM model. These are the same simulations studied in Ma & Fry (2000a). The $`n=2`$ simulation has $`256^3`$ particles and a Plummer force softening length of $`L/5120`$, where $`L`$ is the box length. The $`\mathrm{\Lambda }`$CDM model is spatially flat with matter density $`\mathrm{\Omega }_m=0.3`$ and cosmological constant $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$. This run has $`128^3`$ particles and is performed in a $`(100\mathrm{Mpc})^3`$ comoving box with a comoving force softening length of $`50\mathrm{kpc}`$ for Hubble parameter $`h=0.75`$. The baryon fraction is set to zero for simplicity. The primordial power spectrum has a spectral index of $`n=1`$, and the density fluctuations are drawn from a random Gaussian distribution. The gravitational forces are computed with a particle-particle particle-mesh (P<sup>3</sup>M) code (Ferrell & Bertschinger 1994). We compute the density field $`\delta `$ on a grid from particle positions using the second-order triangular-shaped cloud (TSC) interpolation scheme. A fast Fourier transform is then used to obtain $`\delta `$ in $`k`$-space. The $`k`$-space TSC window function is deconvolved to correct for smearing in real space due to the interpolation, and shot noise terms are subtracted to correct for discreteness effects. We then compute the second and third moments of the density amplitudes in Fourier space.

We show results for the power spectrum as the dimensionless variance $`\mathrm{\Delta }(k)4\pi k^3P(k)/(2\pi )^3`$. A useful dimensionless three-point statistic is the hierarchical three-point amplitude

$$Q(k_1,k_2,k_3)\frac{B(k_1,k_2,k_3)}{P(k_1)P(k_2)+P(k_2)P(k_3)+P(k_3)P(k_1)}.$$
(30)
The three-point amplitude $`Q`$ has the convenient feature that for the lowest nonvanishing result in perturbation theory, $`Q`$ is independent of time and the overall amplitude of $`P`$; for scale-free models with a power-law $`P`$, $`Q`$ is independent of overall scale as well. To lowest order, it follows from equation (12) that the equilateral bispectrum has a particularly simple form, $`B^{(0)}(k)=\frac{12}{7}P_{\mathrm{lin}}^2(k)`$, and we have $`Q^{(0)}(k)=\frac{4}{7}`$, independent of the power spectrum.

### 5.1 Synthetic Halo Replacement

To investigate the numerical effects of limited resolution in the simulations, we have experimented with the distribution of matter in halos identified in the simulations. In these experiments, we keep the locations and masses of the halos unchanged but redistribute the subset of particles which lies within the virial radius $`R_{200}`$ (the radius within which the mean overdensity is 200) of each halo according to a prescribed density profile. We then recompute the two- and three-point statistics $`\mathrm{\Delta }`$ and $`Q`$ from the redistributed particle positions as well as the original non-halo particles, which remain at their original positions. By using density profiles obtained empirically from higher-resolution simulations of individual halos, this recipe allows us to model accurately the inner regions of the halos on scales below the numerical softening length scale while at the same time preserving all the large-scale information available in the large parent simulation. This technique should also be useful for other studies that are sensitive to the inner halo density profiles, for example the ray-tracing method in gravitational lensing.

Ma & Fry (2000a) have used this replacement technique to experiment with synthetic halos that follow a pure power-law profile $`ur^ϵ`$. It is found that $`\mathrm{\Delta }(k)`$ and $`Q(k)`$ at high-$`k`$ indeed obey $`\mathrm{\Delta }(k)k^{2ϵ3}`$ and $`Q(k)k^{3ϵ}`$ as predicted by the simple power-law model of Peebles (1974). The scaling works even in the presence of the full distribution of matter outside the halo cores. Here we extend this replacement technique to more realistic halo profiles of equation (2). Figures 1 and 2 illustrate the effects on the matter power spectrum and bispectrum when the original halos in large cosmological simulations are replaced by synthetic halos with the density profile $`u_{II}=1/(x^{3/2}+x^3)`$ of equation (2). For the $`n=2`$ scale-free model, the concentration parameter is taken to be $`c(M)=3(M_{}/M)^{1/6}`$, which is consistent with Navarro et al. (1997) and has the expected scaling with mass, $`cM^{(3+n)/6}`$, in a scale-free model. For $`\mathrm{\Lambda }`$CDM models, we use $`c(M)=5(M_{}/M)^{1/6}`$ as suggested by Figure 3 of Moore et al. (1999). We note, however, that $`c(M)`$ from various recent simulations has shown a large scatter, and its functional form depends on the exact form of the density profile used. For the $`\mathrm{\Lambda }`$CDM model and form $`u_{II}=1/(x^{3/2}+x^3)`$, for example, a flatter and smaller $`c(M)=3(M_{}/M)^{0.084}`$ appears to be preferred by Jing & Suto (2000) and Navarro et al. (1997). The results of Tormen et al. (1997) and Cole & Lacey (1996) are also only marginally consistent with each other. A more detailed investigation of the different forms of $`c(M)`$ can be found in Ma & Fry (2000c).

In Figures 1 and 2, the agreement at low values of $`k`$ between the original and synthetic halos is excellent, confirming that the correlation functions on larger length scales are insensitive to the spatial distribution of particles in the halo cores. The only significant difference between the simulation and synthetic halos appears at small length scales, where the coarser resolution of the simulation blurs out the structure of the inner halo and results in an inner profile flatter than in equation (2). This effect is manifested in the bending over of the dashed curves for $`P(k)`$ in Figures 1 and 2 at high $`k`$, and is corrected for when the synthetic halos are used.

### 5.2 $`N`$-body Results vs. Analytic Halo Model

We now proceed to compare the predictions of the analytic model of §2 – §4 with the numerical results from cosmological simulations. Figures 3 and 4 show the $`k`$-space density variance $`\mathrm{\Delta }(k)`$ (upper panel) and the three-point amplitude $`Q_{\mathrm{eq}}(k)`$ for equilateral triangles for the $`n=2`$ scale-free model and the $`\mathrm{\Lambda }`$CDM model. The solid black curves are the model predictions computed from equations (3) and (29). The contribution from the single-halo and multiple-halo terms are shown separately as dashed curves. For the density profile, we use the same $`u_{II}=1/(x^{3/2}+x^3)`$ and concentration parameters as in Figures 1 and 2. For the mass function, we use the Press-Schechter formula but reduce its overall amplitude by 25%, which we find necessary in order to match the halo mass functions for our numerical simulations. This overestimation of halo numbers with $`MM_{}`$ by Press-Schechter is a well known result reported in many other studies (see Jenkins et al. 2000 and references therein). The mass limits for the integrals in equations (3) and (29) do not significantly affect the model predictions for the total $`\mathrm{\Delta }`$ or $`Q`$. Raising the lower mass limit does reduce the contribution from lower mass halos and hence lower the high-$`k`$ amplitudes of the multiple halo terms $`\mathrm{\Delta }_{2h}`$, $`Q_{2h}`$, and $`Q_{3h}`$, but these terms make negligible contributions to the total $`\mathrm{\Delta }`$ and $`Q`$.

As discussed in §3 and 4, the nonlinear parts of both the two- and three-point statistics are determined by the dominant 1-halo term because the closely spaced particle pairs and triplets mostly reside in the same halos. The multiple-halo terms are therefore significant only on larger length scales comparable to the separation between halos. Their inclusion, however, is necessary for the transition into the linear regime.

For the $`n=2`$ model in Figure 3, we plot the results against the scaled $`k/k_{\mathrm{nl}}`$, where $`k_{\mathrm{nl}}`$ characterizes the length scale that is becoming nonlinear and is defined by $`_0^{k_{\mathrm{nl}}}d^3kP_{\mathrm{lin}}(a,k)/(2\pi )^3=1`$. Three time outputs are shown, where the expansion factor (1 initially) and $`k_{\mathrm{nl}}`$ (in units of $`2\pi /L`$) are: $`(a,k_{\mathrm{nl}})=(13.45,29),(19.03,14.5)`$, and $`(26.91,7.25)`$ (from left to right). For the two-point $`\mathrm{\Delta }(k)`$, the agreement between the model prediction and the simulations is excellent. The three simulation outputs also overlap well, indicating that self-similarity is obeyed, as reported in Jain & Bertschinger (1998). For the three-point $`Q_{\mathrm{eq}}`$, however, self-similar scaling does not hold as rigorously (Ma & Fry 2000a). It is interesting to note that the analytic prediction agrees most closely with the earliest output $`(a,k_{\mathrm{nl}})=(13.45,29)`$ (green curve). This provides further evidence to the suggestion of Ma & Fry (2000a) that the later outputs of the $`n=2`$ simulation may be affected by the finite volume of the simulation box. For the $`\mathrm{\Lambda }`$CDM model in Figure 4, the analytic model again provides a good match to the $`N`$-body results within the fluctuations among the simulations. We illustrate the numerical effects due to box sizes by showing results from two runs with volume (100 Mpc)<sup>3</sup> and (640 Mpc)<sup>3</sup>. The model predictions extend well beyond the resolution of the simulations.

The real-space two-point correlation function for the $`n=2`$ and $`\mathrm{\Lambda }`$CDM models is shown in Figures 5 and 6. For the halo model predictions, we have chosen to show only the results for the 1-halo term $`\xi _{1h}`$ because this term dominates the interesting nonlinear portion of $`\xi `$. The agreement between the halo model (dashed curves) and the simulations (symbols) is again excellent. For the 2-halo terms $`\xi _{2h}`$, the computation can be done more easily in $`k`$-space as shown in Figures 3 and 4, so we do not include them here.

For comparison, we plot in Figures 3–6 the results from the commonly used fitting formulas for the nonlinear power spectrum (Hamilton et al. 1991; Jain et al. 1995; Peacock & Dodds 1996; Ma 1998; Ma et al. 1999). While the formulas provide a good approximation to $`\mathrm{\Delta }(k)`$ up to $`k/k_{\mathrm{nl}}50`$ for the $`n=2`$ model and $`k20h`$ Mpc<sup>-1</sup> for the $`\mathrm{\Lambda }`$CDM model, the figures show that significant deviations occur at higher $`k`$, and the fitting formula and our current model predict different high-$`k`$ slopes for $`\mathrm{\Delta }(k)`$. Since the high-$`k`$ behavior of the fitting formulas has been constructed to obey the stable clustering prediction, this discrepancy has an important implication for the validity of stable clustering, which we discuss briefly in the next section and at length in Ma & Fry (2000c).

## 6 Discussion

We have constructed a physical model for the correlation functions of the mass density field in which the correlations are derived from properties of dark matter halos. We have described in detail the input, construction, and results of this model in §2 – §5. We now examine more closely its physical meanings and implications in three separate regimes.

On scales larger than the size of the largest halo, the contributions from separate halos dominate, and (by design) the model reproduces the results of perturbation theory. On intermediate scales, $`1/R_{}k1/R_s(M_{})`$, because of the exponential cutoff in the mass function $`dn/dM`$ at the high mass end, the contribution to the volume integrals in equation (18) is dominated by the large-$`r`$ regime where the halo profiles are roughly $`r^3`$. The correlation functions therefore behave approximately as predicted by the power-law model with $`ϵ=3`$, i.e., $`\mathrm{\Delta }k^{2ϵ3}k^3`$ and $`Qk^{3ϵ}`$ constant. This is why $`Q`$ exhibits an approximately flat plateau at intermediate $`k`$ in the bottom panels of Figures 3 and 4.

On the smallest and most nonlinear scales, the correlation functions probe the innermost regions of the halos. Intriguingly, the halo model predicts on these scales a behavior that is different from either the frequently-assumed stable clustering result of $`\mathrm{\Delta }(k)k^\gamma `$ with $`\gamma =(9+3n)/(5+n)`$ (Davis & Peebles 1977), or the power-law profile result of $`\gamma =2ϵ3`$. The implication of departure from stable clustering is significant because all the fitting formulas for the nonlinear $`P(k)`$ in the literature (see §5.2) have been constructed to approach the stable clustering limit at high $`k`$. A more detailed study on the criteria for stable clustering in this model is given in a separate paper (Ma & Fry 2000c).

The origin of the deviation from stable clustering in the model at high-$`k`$ can be understood as follows. For the two-point function, as $`k`$ becomes large, the one-halo integral $`P_{1h}(k)`$ in equation (3) converges before the exponential cutoff, and is dominated by contributions near the mass scale for which $`kR_s=1`$. The behavior now depends on the mass distribution function. The various mass functions discussed in §2.2 have the same general behavior of $`dn/dMM^2\nu ^\alpha e^{\nu ^2/2}`$, where $`\nu =\delta _c/\sigma `$. The Press-Schechter form assumes $`\alpha =1`$ (see eq. ), while others (e.g., Sheth & Tormen 1999; Jenkins et al. 2000) suggest a flatter slope of $`\alpha 0.4`$ for the lower mass halos. Since the scale radius $`R_s`$ depends on mass as $`R_s=R_{200}/cM^{1/3}/M^{(3+n)/6}M^{(5+n)/6}`$, and $`R_s^3\overline{\delta }M`$ (up to logarithmic factors), we find from equation (3) that the power spectrum at high $`k`$ goes as

$$\mathrm{\Delta }(k)\mathrm{\Delta }_{1h}(k)k^3𝑑M\nu ^\alpha \stackrel{~}{u}^2(kR_s).$$
(31)
Changing variables to $`y=kR_sk(M/M_{})^{(5+n)/6}`$, we see that

$$\mathrm{\Delta }(k)k^\gamma ,\gamma =\left(\frac{9+3n}{5+n}\right)\alpha \left(\frac{3+n}{5+n}\right),$$
(32)
where the first term in $`\gamma `$ is the prediction of stable clustering. The departure arises from the factor $`\nu ^\alpha `$ in the mass function, and would vanish only if $`\alpha =0`$ or $`n=3`$. This is the origin of the difference in $`\mathrm{\Delta }(k)`$ at high $`k`$ between the model prediction (solid curves) and the fitting formula (dotted curves) shown in Figures 3 and 4.

For the three-point function, the one-halo integral $`B_{1h}^{\mathrm{eq}}`$ in equation (29) converges (barely, for $`p=\frac{3}{2}`$ and $`n=2`$), giving

$`B^{\mathrm{eq}}(k)`$ $``$ $`k^{\gamma _36},\gamma _3=2\left({\displaystyle \frac{9+3n}{5+n}}\right)\alpha \left({\displaystyle \frac{3+n}{5+n}}\right)`$
$`Q^{\mathrm{eq}}(k)`$ $``$ $`k^{\alpha (3+n)/(5+n)}`$ (33)
This again disagrees with the prediction of stable clustering that $`Q`$ is constant, but it appears to be consistent with numerical simulations as shown in Figures 3 and 4.

For yet higher order correlations, details of the halo profile begin to matter. For $`p=1`$, the pattern of equations (32) and (33) persists to all orders, but for $`p=\frac{3}{2}`$ they apply only for the two- and three-point functions; for four-point and higher functions the nonlinear scale $`M_{}`$ and $`\gamma _n=np3`$ for $`n4`$. Thus there seems to be some potentially interesting behavior that is tested only in the four-point function and higher.

## 7 Summary

We have presented an analytic model for the two- and three-point correlation functions $`\xi (r)`$ and $`\zeta (r_1,r_2,r_3)`$ of the cosmological mass density field and their Fourier transforms, the mass power spectrum $`P(k)`$ and the bispectrum $`B(k_1,k_2,k_3)`$. In this model, the clustering statistics of the density field are derived from a superposition of dark matter halos with a given set of input halo properties. These input ingredients include realistic halo density profiles of equation (2), halo mass distribution of equation (6), and halo-halo spatial correlations of equations (10) and (15). The main results of the model are given by equations (18) and (3) for the two-point statistics $`\xi `$ and $`P`$, and by equations (25) and (29) for the three-point statistics $`\zeta `$ and $`B`$. This model provides a rapid way to compute the correlation functions over all length scales where the model inputs are valid; it also gives a physical interpretation of the clustering process of matter in the universe.

We have tested the validity of this model by comparing its predictions with results from cosmological simulations of an $`n=2`$ scale-free model and a $`\mathrm{\Lambda }`$CDM model. As Figures 3 – 6 illustrate, the model describes well the simulation results spanning the entire range of behavior from the perturbative regime on large scales to the strongly nonlinear regime on small scales. To probe the critical high-$`k`$ range in the deeply nonlinear regime, we have used a halo replacement technique to increase the resolution of the large parent simulations. As Figures 1 and 2 illustrate, this method of replacing the original halos that suffer from numerically softened cores with synthetic halos of analytic profiles is a reasonable way to improve the resolution of numerical simulations. By using density profiles obtained empirically from higher-resolution simulations of individual halos, this recipe allows us to model accurately the inner regions of the halos on scales below the numerical softening length scale, while at the same time preserving all the large-scale information available in the large parent simulation. This technique should also be useful for other studies that depend on the inner halo density profiles, for example, the ray-tracing method in gravitational lensing.

Given that dark matter halos in simulations (and presumably in nature) are not perfectly spherical, cleanly delineated objects, it is intriguing that the model constructed in this paper works as well as it does at matching the simulation results. Nevertheless, this analytic model provides a good qualitative and quantitative description over the entire range of scales covered by the simulation, and it can be used to make predictions beyond these scales. This is the first model prescription that successfully reproduces both two- and three-point mass correlations. We believe that it will prove to be a generally useful framework.

We have enjoyed stimulating discussions with John Peacock and David Weinberg. We thank Edmund Bertschinger for valuable comments and for providing the $`n=2`$ scale-free simulation. Computing time for this work is provided by the National Scalable Cluster Project and the Intel Eniac2000 Project at the University of Pennsylvania. C.-P. M. acknowledges support of an Alfred P. Sloan Foundation Fellowship, a Cottrell Scholars Award from the Research Corporation, a Penn Research Foundation Award, and NSF grant AST 9973461.

## Appendix A Appendix

In this Appendix we display analytic forms for the convolution of the dimensionless profile shape

$$\lambda (x)=d^3yu(y)u(|\text{x}+\text{y}|)$$
(A1)
discussed in §3. These analytic expressions are useful for computing the nonlinear two-point correlation function $`\xi `$ of the mass density field, which is dominated by the 1-halo term $`\xi _{1h}`$ in equation (18) and is related to $`\lambda `$ by

$$\xi (r)\xi _{1h}(r)=𝑑M\frac{dn}{dM}\overline{\delta }^2R_s^3\lambda (r/R_s),\mathrm{for}\xi 1.$$
(A2)
For the type-I profile $`u_I`$ of equation (2), the angular integration in equation (A1) is analytic, and $`\lambda `$ is reduced to a simple integral

$$\lambda _I(x)=\frac{2\pi }{(2p)x}_0^{\mathrm{}}\frac{ydy}{y^p(1+y)^{3p}}\left[\frac{(x+y)^{2p}}{(1+x+y)^{2p}}\frac{|xy|^{2p}}{(1+|xy|)^{2p}}\right].$$
(A3)
For the special case $`p=1`$, this integral can be further reduced to the analytical form

$$\lambda _I(x)=\frac{8\pi }{x^2(x+2)}\left[\frac{(x^2+2x+2)\mathrm{ln}(1+x)}{x(x+2)}1\right],p=1.$$
(A4)
For $`u_{II}`$ of equation (2), we are able to simplify $`\lambda `$ to

$$\lambda _{II}(x)=\frac{2\pi }{x}_0^{\mathrm{}}\frac{ydy}{y^p(1+y^{3p})}F_p(x,y),$$
(A5)
where the function $`F_p(x,y)`$ represents the angular part of the integration in equation (A1) and

$$F_p(x,y)=_{|xy|}^{x+y}\frac{zdz}{z^p(1+z^{3p})}.$$
(A6)
The integral in $`F_p`$ can be reduced to analytic forms for special values of $`p`$. Here we display the six cases $`p=0`$, $`1/2`$, 1, $`3/2`$, 2, and $`5/2`$:

$`F_0`$ $`=`$ $`{\displaystyle \frac{1}{6}}\left\{2\sqrt{3}\mathrm{tan}^1\left[{\displaystyle \frac{1+2(x+y)}{\sqrt{3}}}\right]+\mathrm{ln}\left[{\displaystyle \frac{1(x+y)+(x+y)^2}{1+2(x+y)+(x+y)^2}}\right]\right\}`$ (A7)
$`{\displaystyle \frac{1}{6}}\left\{\text{replace }(x+y)\text{ above with }|xy|\right\}`$
$`F_{1/2}`$ $`=`$ $`{\displaystyle \frac{1}{10}}\{2\sqrt{10+2\sqrt{5}}\mathrm{tan}^1\left({\displaystyle \frac{1+\sqrt{5}4\sqrt{x+y}}{\sqrt{102\sqrt{5})}}}\right)`$ (A8)
$`2\sqrt{102\sqrt{5}}\mathrm{tan}^1\left({\displaystyle \frac{1+\sqrt{5}+4\sqrt{x+y}}{\sqrt{10+2\sqrt{5})}}}\right)`$
$`+4\mathrm{ln}\left(1+\sqrt{x+y}\right)(1+\sqrt{5})\mathrm{ln}\left[1+{\displaystyle \frac{1}{2}}(1+\sqrt{5})\sqrt{x+y}+x+y\right]`$
$`(1\sqrt{5})\mathrm{ln}[1{\displaystyle \frac{1}{2}}(1+\sqrt{5})\sqrt{x+y}+x+y]\}`$
$`{\displaystyle \frac{1}{10}}\left\{\text{replace }(x+y)\text{ above with }|xy|\right\}`$
$`F_1`$ $`=`$ $`\mathrm{tan}^1(x+y)\mathrm{tan}^1(|xy|)`$ (A9)
$`F_{3/2}`$ $`=`$ $`{\displaystyle \frac{1}{3}}\left\{2\sqrt{3}\mathrm{tan}^1\left[{\displaystyle \frac{1+2\sqrt{x+y}}{\sqrt{3}}}\right]+\mathrm{ln}\left[{\displaystyle \frac{1+2\sqrt{x+y}+x+y}{1\sqrt{x+y}+x+y}}\right]\right\}`$ (A10)
$`{\displaystyle \frac{1}{3}}\left\{\text{replace }(x+y)\text{ above with }|xy|\right\}`$
$`F_2`$ $`=`$ $`\mathrm{ln}\left[{\displaystyle \frac{x+y}{1+x+y}}\right]\mathrm{ln}\left[{\displaystyle \frac{|xy|}{1+|xy|}}\right]`$ (A11)
$`F_{5/2}`$ $`=`$ $`{\displaystyle \frac{2}{\sqrt{|xy|}}}{\displaystyle \frac{2}{\sqrt{x+y}}}+\mathrm{ln}\left[{\displaystyle \frac{(1+2\sqrt{x+y}+x+y)|xy|}{(1+2\sqrt{|xy|}+|xy|)(x+y)}}\right]`$ (A12)