text:
# Evolution of Power Spectrum in Non-Gaussian Models

## 1 INTRODUCTION

The power spectrum (or equivalently, the two-point correlation function as its Fourier transform) is the most fundamental measure to quantify matter clustering. Observational determination of it will bring us basic information of our universe. For example we can make constraints on cosmological parameters or on statistical aspects of the initial matter fluctuations that are the seeds of structure formation. Evolution of power spectrum by gravitational instability is also an interesting problem of nonlinear physics.

It is often assumed that the initial density fluctuations obey random Gaussian distributions. Though this is the simplest assumption from statistical point of views, its origin is explained by the simplest inflation scenario. Many analyses of the large-scale structure (LSS) in the universe are performed under this assumption. But other models (including some inflationary models and defect models) predict non-Gaussian initial fluctuations (Vilenkin 1985, Vachaspati 1986, Srednicki 1993, Peebles 1997, Linde & Mukahanov 1997) and possibilities of these models should be observationally investigated and have recently called much attentions (e.g. Peebles 1999a, 1999b, Ferreira, Górski & Magueijo 1998, Robinson, Gawiser & Silk 1998, 2000, Koyama, Soda & Taruya 1999, Willick 2000). These recent analyses involve the large-scale structure, properties (e.g. time evolution, spatial clustering) of clusters and CMB anisotropies. It now becomes clear that defect models cannot generate the observed CMB data by their own account but inflation models can explain them well (e.g. Albrecht 2000). Recent CMB data by BOOMERanG (Lange et al. 2000) and MAXIMA (Balbi et al. 2000) seem to suggest a weaker secondary peak than expected in a simple inflation scenario. Bouchet et al. (2000) discussed that this might be resolved by a hybrid scenario (inflation+cosmic string) in which some component of perturbations come from topological defects.

As for the quantitative measurement of the primordial non-Gaussianity, Contaldi et al. (2000) studied the third-order moment of the temperature anisotropies using 4yr COBE Differential Microwave Radiometer (DMR) data and found their most conservative estimate as $`(\mathrm{\Delta }T)^3=6.5\pm 8.7[10^4\mu \mathrm{K}^3]`$ (95% CL). Even though the magnitude of the error bar is large due to the limited sensitivity and angular resolution, traditional Gaussian model ($`(\mathrm{\Delta }T)^3=0`$) is consistent with their result. Feldman et al. (2000) analyzed the bispectrum of the IRAS PSCz catalog (Saunders et al. 2000) and obtained a constraint for the primordial dimensionless skewness as $`S\delta ^2/\delta ^2^{3/2}<0.52`$ (95% CL) for $`\chi _N^2`$-models. Physics of CMB are much simpler than that of the galaxy distribution. The signal of the primordial non-Gaussianity could be masked by effects of biasing or nonlinear gravitational evolution (Verde et al. 2000). Observed data of both CMB and LSS will be significantly improved from MAP and Plank satellites, SDSS and 2dF surveys. Verde et al. (2000) insisted that CMB would be a better probe of the primordial non-Gaussianity for several non-Gaussian models and LSS data would be useful to constrain the biasing if primordial Gaussianity is suggested by CMB map.

So far nonlinear effects of matter clustering from non-Gaussian initial conditions are studies by numerical simulations (White 1999, Robinson & Baker 2000 and references therein) or by perturbative methods (e.g. Fry & Scherrer 1994, Chodorowski & Bouchet 1996, Scoccimarro 2000, Durrer et al. 2000). In this article we study evolution of the matter power spectrum for non-Gaussian models using the dark-matter-halo approach that can be traced back to Peebles (1974) and McClelland and Silk (1977) (see also Scherrer & Bertschinger 1991). For Gaussian models this approach is confirmed to reproduce numerical results very well from linear to nonlinear scales and regarded as an excellent tool to explore the matter clustering. It seems interesting to see how the nonlinear spectrum can be different within the currently allowed region of the primordial non-Gaussianity.

Many works have been recently performed with this approach for Gaussian models. For example nonlinear evolution of the power spectrum or the bispectrum of density field is studies by Ma & Fry (2000a, 2000b, 2000c), Seljak (2000), White (2000), Peacock & Smith (2000), Scoccimarro et al. (2000) and weak lensing effects by Cooray, Hu & Miralda-Escude (2000) and Cooray & Hu (2000). Komatsu & Kitayama (1999) investigated the Sunyaev & Zeldovich effect using a similar method. Here we simply extend this approach to non-Gaussian models. This approach is made from three basic ingredients (i) the mass function of dark-matter-halos, (ii) the bias parameter of halos relative to the large-scale density fluctuations and (iii) the density profile of halos. Numerical analysis (Robinson & Baker 2000) supports that a straightforward extension of the Press & Schechter (1974) mass function is effective for non-Gaussian models. But validity of simple extensions of two other elements have not been checked numerically. Therefore we make a careful analysis and limit our investigation in the range where results are expected not to depend largely on the details of these two elements. Our goals in this article are set to the following two points. Firstly, we intend to roughly understand the effects of primordial non-Gaussianity on the weakly nonlinear evolution of the power spectrum. Secondary, we estimate the difference of the evolved spectra within the current observational constraint of the primordial non-Gaussianity.

This article is organized as follows. In §2 we describe basic properties of the dark-matter-halo approach with its basic three ingredients, namely, the mass function of halos (§2.1), the bias parameter (§2.2) and the density profiles of halos (§2.3). In §3 we analyze toy models to clarify the general effects of non-Gaussianity on the weakly nonlinear evolution of power spectrum. In §4 we investigate realistic models. We study a flat $`\mathrm{\Lambda }`$-dominated model ($`\mathrm{\Omega }_0=0.3,\lambda _0=0.7`$ and $`h=0.7`$) with a fixed linear CDM spectrum and evaluate evolution of power spectra for log-normal probability distribution functions that are motivated by work of Robinson & Baker (2000). §5 is devoted to a brief summary. In appendix A we discuss the scale dependence of the skewness parameter for $`\chi _N^2`$-models.

## 2 DARK MATTER HALO APPROACH

In the dark-matter-halo approach the nonlinear density field is given by superposition of halos with various masses. The density profile of each halo is assumed to be determined by its mass. The nonlinear power spectrum $`P_{NL}(k)`$ is constituted by two terms as

$$P_{NL}(k)=P_{2h}(k)+P_{1h}(k).$$
(1)
The two-halo term $`P_{2h}(k)`$ counts contributions of two points coming from two different halos and is given by

$$P_{2h}(k)\left[\rho _b^1𝑑mN(m)mb(m)u(k,m)\right]^2P_L(k),$$
(2)
where $`P_L(k)`$ is the linear power spectrum. The one-halo term $`P_{1h}(k)`$ counts two particles within the same halo

$$P_{1h}(k)\rho _b^2𝑑mN(m)m^2u(k,m)^2.$$
(3)
In the above equations $`\rho _b`$ is the background density of the universe and $`N(m)dm`$ is the number density (mass function) of halos with mass from $`m`$ to $`m+dm`$. The factor $`b(m)`$ is the bias parameter and represents distribution of halos relative to the large-scale density fluctuations. The function $`u(k,m)`$ is the Fourier transform of the density profile of a dark-matter-halo with mass $`m`$. The original form of the dark-matter halo approach was proposed by Peebles (1974). This approach has recently called much attention as our understandings of its ingredients have largely developed. In the following three subsections we study these basic ingredients.

### 2.1 Mass Function

We use the Press & Schechter (1974) formalism for the mass function $`N(m)`$ of collapsed objects. This formalism was first given for Gaussian random fields but can be straightforwardly extended to general models (e.g. Lucchin & Matarrese 1989, Chiu, Ostriker & Strauss 1998, Koyama, Soda & Taruya 1999, Robinson, Gawsier & Silk 2000, Willick 2000, Matarrese, Verde & Jimenez 2000). For simplicities we assume that the linear one-point probability distribution function (hereafter PDF) $`p(\delta ,\sigma )`$ for the smoothed density field is written as

$$p(\delta ,\sigma )d\delta =\sigma ^1p(\delta /\sigma )d\delta ,$$
(4)
which means that the shape of the PDF is scale invariant. This ansatz is correct for the random Gaussian fluctuations, but not guaranteed for general non-Gaussian models. We comment on this scale dependence in §4 and appendix A.

The volume fraction of points that belong to collapsed objects with mass larger than $`m`$ is given as follows

$$F^1_{\delta _c/\sigma (m)}^{\mathrm{}}p(\nu )𝑑\nu ,$$
(5)
where $`\delta _c3(12\pi )^{2/3}/20=1.69`$ is the critical linear density contrast for the spherical top-hat collapse and $`\sigma (m)`$ is the root-man-square mass fluctuation in a sphere of mass scale $`m`$. $`F`$ is the normalization factor to account for all masses in the universe

$$F=_0^{\mathrm{}}p(\nu )𝑑\nu .$$
(6)
For a symmetric profile $`p(\nu )=p(\nu )`$ we have $`F=1/2`$ (see also Bond et al. 1991, Lacey & Cole 1993). From equation (5) the number density of halos with mass from $`m`$ to $`m+dm`$ is given as

$$N(m)dm=F^1\rho _bm^1\frac{}{m}_{\delta _c/\sigma (m)}^{\mathrm{}}p(\nu )𝑑\nu =F^1\rho _bp(\nu )m^1\nu \frac{d\mathrm{ln}\sigma (m)}{dm},$$
(7)
where we have denoted the normalized density contrast $`\nu \delta _c/\sigma (m)`$. We define the nonlinear mass $`m_{}`$ by equation $`\sigma (m_{})=\delta _c`$. For power-law models with $`P_L(k)k^n`$ the linear mass fluctuation $`\sigma (m)`$ is given as

$$\sigma (m)=\delta _c\left(\frac{m}{m_{}}\right)^{(n+3)/6}.$$
(8)
In this subsection we have followed the basic analysis of the Press & Schechter method. For random Gaussian fluctuations some refinements have been proposed to reproduce N-body data better (e.g. Sheth & Tormen 1999, Jenkins et al. 2000). If we write down the mass function in the following form

$$N(m)=A\sqrt{\frac{2}{\pi }\alpha }\frac{d\mathrm{ln}\sigma }{d\mathrm{ln}m}\frac{\rho _b}{m^2}\nu e^{\alpha \nu ^2/2}(1+(\alpha \nu ^2)^p).$$
(9)
The original Press & Schechter function corresponds to $`(A,p,\alpha )=(0.5,0,1)`$. The mass function given by Sheth & Tormen is given by parameters $`(A,p,\alpha )=(0.322,0.3,0.707)`$.

Robinson & Baker (2000) have recently confirmed that the simple extension of the Press & Schecheter formula (7) can accurately fit the mass function of clusters for several non-Gaussian models with typical error $`25`$ percent (see also Avelino & Viana 2000). This result encourages our analytical study.

### 2.2 Bias Parameter

Next we study the bias parameter $`b(m)`$ using (i) the peak-background splitting method (Kaiser 1984) and (ii) a simple extension of the Press & Schechter method (Mo & White 1997, Mo, Jing & White 1997, Sheth & Lemson 1999, Sheth & Tormen 1999). When a large-scale perturbation $`\delta _{LS}`$ is added to the local (small scale) density field $`\delta _{SS}(𝒙)`$, the density contrast for $`\delta _{SS}(𝒙)`$ that is critical to the spherical top-hat collapse becomes $`\delta _c\delta _{LS}`$. Then the Eulerian bias parameter for halos with mass $`m`$ becomes (Mo & White 1997, Robinson, Gawiser & Silk 1998, Koyama, Soda & Taruya 1999)

$$b(m)=1\frac{d\mathrm{ln}N(m)}{d\delta _c}=1\frac{1}{\delta _c}\left(\frac{p^{}(\nu )}{p(\nu )}\nu +1\right).$$
(10)
From equations (7) and (10) we can easily confirm the following equation 111This equation holds for general non-scaling PDF in the form $`p(\delta ,\sigma )`$.

$$\rho _b^1_0^{\mathrm{}}b(m)N(m)m𝑑m=1.$$
(11)
As we see in §4 this equation is useful to know the behavior of the nonlinear power spectrum $`P_{NL}(k)`$ at large scale (small $`k`$). There are several analytical works about the bias parameter for non-Gaussian models (Robinson, Gawiser & Silk 1998, Koyama, Soda & Taruya 1999). But validity of the formulas (10) for a non-Gaussian PDF have not been checked with numerical simulations. But we expect that it works well, as this formula is simply derived from the Press & Schechter method that has been checked numerically for some non-Gaussian models as commented in the previous subsection.

For the mass function (9) given for random Gaussian fluctuations we have

$$b(\nu )=1+\frac{\alpha \nu ^21}{\delta _c}+\frac{2p}{\delta _c(1+(\alpha \nu ^2)^p)}.$$
(12)
### 2.3 Density Profile of Dark Matter Halo

Recent numerical simulations suggest that the average density profile of cold-dark-mater (CDM) halos has nearly universal shape $`\rho _m(r)`$ characterized by its mass $`m`$ (Navarro, Frenk & White 1996, 1997, Moore et al. 1998, 1999, Kravsov et al. 1998, Fukushige & Makino 1998, Jing & Suto 2000). These numerical simulations are performed with Gaussian initial conditions, but here we assume that the density profile depends only weakly on the linear PDF as the halo profile would strongly reflect nonlinear gravitational dynamics rather than the initial conditions. This assumption may be a potential weakness of our dark-matter-halo approach especially at strongly nonlinear regimes. However as we see in §4, the power spectrum $`P_{NL}(k)`$ at weakly nonlinear regimes does not depend details of the halo profile. In this article we adopt the universal profile of dark matter halo given in Navarro, Frenk & White (1997) as follows

$`u_m(r)`$ $``$ $`{\displaystyle \frac{\rho _m(r)}{m}}`$ (13)
$`=`$ $`{\displaystyle \frac{fc^3}{4\pi R_{vir}^3}}{\displaystyle \frac{1}{(cr/R_{vir})^\alpha (1+cr/R_{vir})^{3+\alpha }}},`$ (14)
where $`R_{vir}`$ is the virial radius of the halo. The halo mass $`m`$ is written by $`R_{vir}`$ as $`m=4\pi R_{vir}^3\rho _b\mathrm{\Delta }/3`$ with the average density contrast $`\mathrm{\Delta }=200`$ (Einstein de-Sitter model) and $`\mathrm{\Delta }=340`$ (flat model with $`\mathrm{\Omega }_0=0.3`$ and $`\lambda _0=0.7`$). The concentration parameter $`c`$ determines the transition of inner and outer regions of a halo and depends on its mass $`m`$, roughly reflecting the formation epoch of the halo (Ma & Fry 2000b). We have $`c=\beta (n)\sigma (m)`$ with a coefficient $`\beta (n)`$ that depends on the linear power spectrum. The normalization factor $`f`$ in equation (14) is given by

$$f=(\mathrm{ln}(1+c)c/(1+c))^1.$$
(15)
The density profile (14) in central region becomes $`\rho (r)r^\alpha `$. The original result in Navarro, Frenk & White (1996) was $`\alpha =1`$. But some other simulations show different results. For example, Moore et al. (1998, 1999) found $`\alpha =1.5`$ and Jing & Suto (2000) argued that the index $`\alpha `$ is not a universal value but depends on the mass scale of halos. As we see later, the nonlinear evolution of power spectrum in weakly nonlinear regime is not sensitive to the details of the halo profile (see also Seljak 2000). Thus we use $`\alpha =1`$ in this article.

Finally the Fourier transformed function $`u(k,m)`$ is defined in terms of the real space density profile $`u_m(r)`$ as

$`u(k,m)`$ $`=`$ $`{\displaystyle d^3x\mathrm{exp}(i𝒙𝒌)u_m(r)}`$ (16)
$`=`$ $`4\pi {\displaystyle \frac{\mathrm{sin}kr}{kr}r^2u_m(r)𝑑r}.`$ (17)
For the profile (14) an explicit analytic formula of $`u(k,m)`$ was given by Scoccimarro et al. (2000). At small wave number $`k`$ the above Fourier transformation becomes a simple volume integral and we have $`u(0,m)=1`$ due to the definition of $`u_m(r)`$ given in equation (13).

## 3 TOY MODELS

In this section we discuss effects of initial non-Gaussianity on the evolution of the power spectrum in Einstein de-Sitter background. We intend to extract out their general relations. From this standpoint we do not try to mimic realistic PDFs so seriously but study toy models that are different from Gaussian distribution. We also limit our analysis to scale-free initial spectra. With results obtained in this section we study models relevant for observational cosmology in §4.

### 3.1 Density Distribution Functions

We use the following four PDFs for the linear smoothed density fields. As commented in §2.1 we assume that the PDFs are scale invariant. The first example is the standard Gaussian model (G-model),

$$p_G(\nu )=\frac{1}{\sqrt{2\pi }}\mathrm{exp}\left(\frac{\nu ^2}{2}\right)\mathrm{}<\nu <\mathrm{}.$$
(18)
This PDF is symmetric around the mean density $`\nu =0`$ and we have $`F=0.5`$. Here the factor $`F`$ is defined in equation (6).

The second example is the Chi-square ($`\chi ^2`$)-model (C-model) with two degrees of freedom. A $`\chi ^2`$-model is often adopted as a typical non-Gaussian model and its origin in cosmological context was recently proposed by Peebles (1999a) using an inflation model. The PDF for the unsmoothed field is

$$p_C(\nu )=\{\begin{array}{cc}\mathrm{exp}(\nu 1)\hfill & (1<\nu <\mathrm{}),\hfill \\ 0\hfill & (\nu 1),\hfill \end{array}$$
(19)
with $`F=1/e`$. We should notice that the PDF for the smoothed density field is same as the above unsmoothed one. This point is discussed in the next section and appendix A.

Next example is a log-normal distribution. This distribution is characterized by a parameter $`\alpha (0)`$ as follows

$$p_{LN}(\nu ,\alpha )=\{\begin{array}{cc}\frac{\alpha }{(1+\nu \alpha )\sqrt{2\pi \mathrm{ln}(1+\alpha ^2)}}\mathrm{exp}\left[\frac{\left\{\mathrm{ln}(\sqrt{1+\alpha ^2}(1+\nu \alpha ))\right\}^2}{2\mathrm{ln}(1+\alpha ^2)}\right]\hfill & \alpha ^1<\nu <\mathrm{}\text{,}\hfill \\ 0\hfill & \nu \alpha ^1\text{.}\hfill \end{array}$$
(20)
At $`\alpha =0`$ we recover the Gaussian distribution $`p_G(\nu )=p_{LN}(\nu ,0)`$. The normalization factor $`F`$ for this distribution is given in terms of the error function as

$$F=\frac{1}{2}\mathrm{erfc}\left[\frac{1}{2}\sqrt{\frac{\mathrm{ln}(1+\alpha ^2)}{2}}\right],$$
(21)
where the error function is defined as $`\mathrm{erfc}(x)=\sqrt{2/\pi }_x^{\mathrm{}}\mathrm{exp}(y^2/2)𝑑y`$.

Here we add two models generated from equation (20). We define the LN1-model by

$$p_{LN1}(\nu )=p_{LN}(\nu ,1)$$
(22)
This profile has larger positive tail than the G-model. We have $`F=0.3386`$. The final example (LN2-model) is obtained by a simple replacement of the previous model as

$$p_{LN2}(\nu )p_{LN1}(\nu )=p_{LN}(\nu ,1).$$
(23)
Note that we have $`p_{LN2}(\nu )=0`$ for $`\nu >1`$. Existence of halos with mass larger than $`m_{}`$ is prohibited in this model. As we see in the following section, the shape of the positive tail is very important to study nonlinear evolution of the power spectrum. The probability for $`\nu 3`$ is $`1.34\times 10^3`$ (G-model), $`1.83\times 10^2`$ (C-model), $`1.87\times 10^2`$ (LN1-model) and 0 (LN2-model). For $`\nu 5`$ we have $`2.87\times 10^7`$ (G-model), $`2.48\times 10^3`$ (C-model), $`5.11\times 10^3`$ (LN1-model) and 0 (LN2-model). We can easily confirm that all of these models have vanishing means

$$\nu =_{\mathrm{}}^{\mathrm{}}\nu p(\nu )𝑑\nu =0,$$
(24)
and are normalized properly

$$_{\mathrm{}}^{\mathrm{}}p(\nu )𝑑\nu =1,\nu ^2=_{\mathrm{}}^{\mathrm{}}\nu ^2p(\nu )𝑑\nu =1.$$
(25)
The skwness parameter $`S\nu ^3/\nu ^2^{3/2}`$ is a fundamental measure to characterize PDFs. 222Definition of skewness parameter is different from that usually used for Gaussian initial condition (Bouchet et al. 1992).  We have $`S=0`$ (G-model), 2 (C-model), 4 (LN1-model) and $`4`$ (LN2-model). For the general log-normal distribution $`p_{LN}(\nu ,\alpha )`$ we have

$$S_{LN}=\alpha (3+\alpha ^2).$$
(26)
We use this relation in section 4.

### 3.2 Results

Using ingredients described in §2 and various linear PDFs given in §3 we calculate the nonlinear power spectrum $`P_{NL}(k)`$ with equations (1), (2) and (3). We use power-law models for initial spectra

$$P_L(k)=Ak^n,$$
(27)
where $`A`$ is the normalization factor but irrelevant for present analysis. Time evolution of power spectrum can be reduced to rescaling of spatial length due to our scale-free initial conditions. We introduce the nonlinear wave-number $`k_{NL}`$ that is defined by

$$\frac{1}{2\pi ^2}_0^{k_{NL}}P_L(k)k^2𝑑k=1.$$
(28)
This wave-number is the most fundamental scale for the present analysis. In this article we call the weakly nonlinear scale for the wave number $`k`$

$$0<kk_{NL}.$$
(29)
We first calculate the nonlinear power spectra $`P_{NL}(k)`$ for random Gaussian fields and compare our analytical predictions with fitting formula that is obtained from N-body simulations with the same Gaussian initial conditions (Peacock & Dodds 1996, see also Hamilton et al. 1991, Jain, Mo & White 1995). We fix the coefficient $`\beta (n)=3`$ for the concentration parameter $`c(m)`$, but our results are insensitive to this parameter. It is found that our model cannot reproduce the numerical results for spectral index with $`n=0`$ in contrast to the case with $`n\stackrel{<}{}\text{ }1`$. We can see this tendency from figures 1.a to 1.c. It is not clear why this dependence appears. Note that the evolved spectrum (thick solid line) for $`n=1`$ obtained from the fitting formula is smaller than the linear spectrum (dash-dotted line) as predicted by the second-order perturbation theory (e.g. Makino, Sasaki & Suto 1992). In the framework of the dark-matter-halo approach this fact might be related to the possibility that the total halo mass fraction does not converge to unity (Jenkins et al. 2000) and the constraint (11) does not hold. As we see below, without the constraint (11) the two-halo term $`P_{2h}(k)`$ can be smaller than the linear spectrum $`P_L(k)`$ with nonnegligible amount at weakly nonlinear scale.

Here we limit our analysis for spectral indexes with $`n1`$ where our approach is expected to work well. These are close to the effective slope of a typical CDM power spectrum at weakly nonlinear regime $`\sigma \stackrel{<}{}\text{ }1`$. In figures 1.a to 1.c we plot the nonlinear power spectrum $`P_{NL}(k)`$ in the form $`\mathrm{\Delta }(k)`$ defined by

$$\mathrm{\Delta }(k)\frac{P_{NL}(k)k^3}{2\pi ^2},$$
(30)
for spectral indexes $`n=2,1.5`$ and $`1`$. In these figures contributions from the one-halo term $`P_{1h}(k)`$ and the two-halo term $`P_{2h}(k)`$ are presented separately. Note that our analytical prediction for the Gaussian models reproduce the results from N-body simulations well in the range $`kk_{NL}`$ for $`n1.5`$ This fact indicates validity of the dark-matter-halo approach. For $`n=2`$ model we also evaluate the nonlinear spectrum with using generalized formulas (9) (12) as in Seljak (2000). Our results shows good agreement with this results. in figures 1.a to 1.c w have used the fitting formula of Peacock & Dodds (1996). We find that the fitting formula of Jain et al. (1995) agrees with our analytical result better for $`n=2`$ model but worse for $`n=1`$ model.

Next we discuss behavior of the nonlinear power spectrum $`P_{NL}(k)`$ for various PDFs given in §3.1. We present numerical results in figures 1.a to 1.c. These are obtained by fully using the dark-matter-halo approach. In this calculation we evaluate all the three elements, the mass function $`N(M)`$, the bias parameter $`b(M)`$ and the halo profile $`u(k,M)`$. But information of the mass function is most important at weakly nonlinear regime, as we see below.

For all models the two-halo term in the range $`k\stackrel{<}{}\text{ }k_{NL}`$ is almost same as the linear power spectrum $`P_L(k)`$, as shown on figures for the two-halo term. In this regime we have $`u(k,m)1`$ in the dominant contribution of the integral in equation (2) (see also figure 2). Then we obtain the following relation from equation (11)

$$_0^{\mathrm{}}𝑑mN(m)mb(m)u(k,m)_0^{\mathrm{}}𝑑mN(m)mb(m)=\rho _b,$$
(31)
and

$$P_{2h}(k)P_L(k),$$
(32)
(see also Sheth et al. (2000) for a same kind of analysis for the two-point correlation function).

Therefore we can roughly understand the weakly nonlinear evolution of power spectrum only by analyzing the one-halo term

$$P_{1h}(k)=\rho _b^2m^2N(m)u(k,m)^2𝑑m.$$
(33)
As shown in figures 1.a to 1.c, the contribution of $`\mathrm{\Delta }(k)`$ from one-halo term $`P_{1h}(k)k^3/(2\pi ^2)`$ is nearly on a straight line proportional to $`k^3`$. This means $`P_{1h}(k)const`$ and indicates that $`k`$-dependence caused by the halo profile $`u(k,m)`$ is weak. Thus we approximate the one-halo term at weakly nonlinear regimes $`k\stackrel{<}{}\text{ }k_{NL}`$ by setting $`u(k,m)=1`$

$$P_{1h}(k)P_{1h}(0)=\rho _b^2𝑑mN(m)m^2.$$
(34)
We have a constraint for the function $`mN(m)`$ due to the conservation of mass as

$$_0^{\mathrm{}}𝑑mN(m)m=\rho _b,$$
(35)
and $`P_{1h}(0)`$ in equation(34) becomes larger for a mass function $`N(m)`$ with more abundant massive objects. Thus we can roughly understand the weakly nonlinear correction of power spectrum using information of mass function in the form given in $`P_{1h}(0)`$.

Note that the present analysis at weakly nonlinear regime $`kk_{NL}`$ based on equations (32) and (33) does not depend on details of the halo density profile. Therefore in this regime we should not be troubled seriously with potential problems about the reliability of the universal profile for non-Gaussian models. This also means that our analysis is not sensitive to the bias parameter $`b(m)`$ (eq.). In the region where $`u(k,m)=1`$ is a good approximation, we need the bias parameter only in the form of equation (11). But this equation represent a simple requirement that the large-scale fluctuation is analyzed perturbatively $`P_{NL}(k)P_L(k)`$. Here we should notice that we have not taken into account the nonlinear biasing effects or the exclusion effects of halo-halo clustering. These effects might be important in weakly nonlinear regime (Scoccimarro et al. 2000), but we neglect them in this article. At present there is no simple prescription to deal with these two opposite effects. This is a generic problem of the dark-matter-halo approach.

Using results for $`n=2`$ model we discuss effects of the linear PDF on the weakly nonlinear evolution of the power spectrum. As shown in figure 1.a the nonlinear correction for the LN2-model is very weak. In this model there are no halos with $`mm_{}`$ (see figure 2) and the integral $`N(m)m^2𝑑m`$ is much smaller than other models. The weakly nonlinear effect is strongest for the LN1-model as expected from its mass function in figure 2. Even at $`kk_{NL}`$ magnitude of the nonlinear power spectra can differ by a factor of 10 depending on the initial PDFs.

Next we compare the behavior of nonlinear power spectrum $`P_{NL}(k)`$ for different spectral indexes $`n`$. As $`n`$ becomes larger, the nonlinear correction becomes smaller (figures 1a-1c). All nonlinear power spectra $`P_{NL}(k)`$ for $`n=1`$ model are very close to the linear spectrum $`P_L(k)`$, but these for $`n=2`$ model depend largely on the adopted PDF. This $`n`$-dependence can be also understood by using the integral $`P_{1h}(0)`$ given in equation (34). We have the following simple relation for the root-mean-square mass fluctuations $`\sigma (m)`$

$$\frac{\delta _c}{\sigma (m)}=\nu =\left(\frac{m}{m_{}}\right)^{(n+3)/6}.$$
(36)
As the index $`n`$ increases, abundance of large mass objects becomes smaller. Then the magnitude of the one-halo term $`P_{1h}(0)`$ decreases and we have smaller nonlinear corrections.

## 4 CDM SPECTRUM

In this section we study nonlinear evolution of the power spectrum for realistic models. We fix the cosmological parameters at $`h=0.7`$, $`\mathrm{\Omega }_0=0.3`$ and $`\lambda _0=0.7`$ that are compatible with CMB (e.g. Lange et al. 2000) and high-redshift SNe data (Riess et al. 1998, Perlmutter et al. 1999). We use the primordially scale invariant ($`n=1`$) spectrum with CDM transfer function from Bardeen et al. (1986). The power spectrum is COBE-normalized and we have $`\sigma _80.9`$. For these parameters the nonlinear wavenumber $`k_{NL}`$ defined in equation (28) becomes $`k_{NL}=0.19`$Mpc.

Robinson & Baker (2000) studied evolution of the cluster abundance in representative non-Gaussian models, such as, Peebles isocurvature cold-dark-matter (ICDM) model or cosmic string model. This ICDM model corresponds to a $`\chi ^2`$-model with one degree of freedom (see appendix A). Peebles (1999b) show that the power spectrum for this model with parameters $`\mathrm{\Omega }_0=0.2`$, $`\lambda =0.8`$ and $`n_\varphi =2.4`$ (see eq. \[A14\] for definition of $`n_\varphi `$) is in good agreement with observed spectrum of CMB and the large-scale structure.

Robinson & Baker (2000) also analyzed PDFs of smoothed linear density fields for these non-Gaussian models and found that the PDFs are well fitted by the log-normal distributions that are nearly scale independent at weakly nonlinear scale. For example, they did not detect scale dependence for the ICDM model and the string+HDM model (their table 1). Thus, motivated by their work, we investigate the log-normal PDFs given in equation (20) and neglect the scale dependence. We characterize the non-Gaussianity of a log-normal PDF by the skewness parameter $`S\delta ^2/\delta ^2^{3/2}`$ (see eq.). For models in Robinson & Baker (2000) we can estimate the parameter $`S`$ from their table 1. We have $`S2.2`$ (somewhat smaller than result $`2.5`$ given in Peebles 1999b) for the ICDM model with $`n_\varphi =2.4`$ and $`S0.5`$ for the string+HDM model. The parameter $`S`$ for the ICDM model is larger than the observed constraint $`S<0.52`$ (95% CL) given in Feldman et al. (2000). A simple improvement is to increase the number $`N`$ of scalar fields. This number $`N`$ corresponds to the degree of freedom in $`\chi _N^2`$-model. As discussed in appendix A, the parameter $`S`$ scales are $`SN^{1/2}`$ and the required number of freedom $`N`$ is $`N\stackrel{>}{}\text{ }(2.2/0.52)^220`$ (Feldman et al. 2000).

In figure 3.a we show the nonlinear spectrum at the wave number $`k=k_{NL}`$ as a function of the parameter $`S`$. The concentration parameter is fixed at $`c(m)=9(m/m_{})^{0.13}`$ (Bullock et al. 1999). We first calculate the ratio

$$\frac{P_{NL}(k)}{P_L(k)}\frac{P_{2h}(k)+P_{1h}(k)}{P_L(k)},$$
(37)
in the full dark-matter-halo approach. We also evaluate the simple approximation with putting $`u(k,M)=1`$, namely

$$\frac{P_{NL}(k)}{P_L(k)}1+\frac{P_{1h}(0)}{P_L(k)}$$
(38)
to check the effects of the bias parameter $`b(M)`$ and the halo inner density profile $`u(k,M)`$. As commented before, validity of these two elements for non-Gaussian models have not been clarified numerically so far. The simple approximation (38) shows reasonable agreement with the full result (37) and we can expect that our analysis is not seriously affected by uncertainty of two elements in the framework of the dark-matter-halo approach. Figure 3.a shows that nonlinear spectrum depends strongly on the primordial non-Gaussianity measured by the parameter $`S`$. For the Gaussian model ($`S=0`$) we have $`P_{NL}/P_L1.3`$ but $`P_{NL}/P_L2`$ for $`S1`$. For the value $`S=0.52`$ (on 95% CL of Feldman et al. 2000) the correction becomes $`P_{NL}/P_L\stackrel{<}{}\text{ }1.55`$ and about 20% larger than the Gaussian model. We also calculate the ratio $`P_{NL}/P_L`$ for the Gaussian model using the fitting formula of Peacock & Dodds (1996). The result is $`1.2`$ and close to our result $`1.3`$.

In figure 3.b we show the nonlinear spectrum at the wavenumber $`k=0.66k_{NL}`$ where the linear CDM spectrum becomes $`P_L(0.66k_{NL})=2P_L(k_{NL})`$. Apparently the correction becomes smaller than at $`k=k_{NL}`$. We can easily understand this fact. At weakly nonlinear scale we have

$$P_{NL}(k)P_L(k)P_{1h}(0)$$
(39)
and the nonlinear correction behaves as $`P_{NL}/P_L1+P_{1h}(0)/P_L(k)`$ in the dark-matter-halo approach.

## 5 SUMMARY

In this article we have studied evolution of the matter power spectrum in non-Gaussian models of structure formation. We have used the dark-matter-halo approach that is an excellent method for analyzing matter clustering from linear to nonlinear scales (Seljak 2000, Ma & Fry 2000b). This approach contains three basic ingredients (i) the mass function $`N(m)`$ of dark-matter-halos, (ii) the bias factor $`b(m)`$ of halos relative to the large-scale density fluctuations and (iii) the density profiles $`u(k,m)`$ of halos. Recent numerical simulations (Robinson & Baker 2000) have suggested that a simple extension of the Press & Schechter formula provides good fits to mass function of clusters also in non-Gaussian models. For the bias factor $`b(m)`$ and the density profiles $`u(k,m)`$ it has not checked numerically whether we can simply extend the formulas of Gaussian models to non-Gaussian models. Therefore we take a moderate position to limit our investigation up to weakly nonlinear regimes $`k\stackrel{<}{}\text{ }k_{NL}`$. In these regions our analysis supports that the evolved power spectrum does not depend on details of the bias factor $`b(m)`$ or the density profile $`u(k,m)`$.

The nonlinear power spectrum predicted by the dark-matter-halo approach is constituted by two terms. The two-halo term counts two points in different halos and the one-halo term represents contributions of two points within same halos. We have shown that the two-halo term is almost same as the linear power spectrum in weakly nonlinear regimes and nonlinear correction mainly comes from the one-halo term. This term is roughly approximated by the integral of mass function $`N(m)`$ as $`N(m)m^2𝑑m`$. Using this integral we can understand various aspects of the weakly nonlinear effect. For example, nonlinear correction is larger than a PDF with more positive tail or smaller spectral index $`n`$.

We also study a COBE-normalized CDM spectrum in a flat model ($`\mathrm{\Omega }_0=0.3,\lambda _0=0.7`$ and $`h=0.7`$) with log-normal PDFs. This shape of PDF is motivated by work of Robinson & Baker (2000) who found numerically that PDFs of many representative non-Gaussian models are fitted well by log-normal distributions. We show that the nonlinear spectrum $`P_{NL}`$ at $`k=k_{NL}0.2\mathrm{Mpc}^1`$ becomes $`P_{NL}/P_L=1.3`$ for the Gaussian model but $`1.55`$ for the log-normal model with $`S=0.52`$ that is on 95% CL of the observed constraint by Feldman et al. (2000). These results suggest that the ratios $`P_L/P_{NL}`$ for $`\mathrm{\Lambda }`$CDM model at $`k<k_{NL}0.2\mathrm{Mpc}^1`$ could not largely (say factor of 2) vary within currently allowed region of the primordial non-Gaussianity.

The author would like to thank an anonymous referee for valuable comments to improve this manuscript. My work is supported by Japanese Grant-in-Aid No. 0001416.

## Appendix A SCALE DEPENDENCE OF THE SKEWNESS IN $`\chi _N^2`$-MODELS

The scale dependence of a non-Gaussian PDF becomes highly complicated as it is generally very difficult to deal with infinite degrees of freedom of nonlocal quantities. Here we study the scale dependence of the PDF using the skewness parameter for $`\chi _N^2`$ model.

Let us discuss $`\chi _N^2`$ distribution with $`N`$ degrees of freedom. The unsmoothed density field $`\rho (𝒙)`$ is written by $`N`$ independent Gaussian random fields $`\varphi _i(𝒙)`$ ($`i=1,\mathrm{},N`$)

$$\rho (𝒙)=\underset{i=1}{\overset{N}{}}\varphi _i(𝒙)^2.$$
(A1)
We assume that these $`N`$ fields have same statistical character. Moments of the unsmoothed field $`\rho (𝒙)`$ is given in terms of the variance $`\varphi _i^2`$ of the basic Gaussian fields $`\varphi _i`$ (the bracket $``$ represents ensemble average). We can easily calculate their moments as

$`\rho (𝒙)`$ $`=`$ $`N\varphi _i^2,\delta \rho (𝒙)^2=2N\varphi _i^2^2,`$
$`\delta \rho (𝒙)^3`$ $`=`$ $`8N\varphi _i^2^3,\delta \rho (𝒙)^4=(48N+12N^2)\varphi _i^2^4.`$ (A2)
where we have defined the overdensity $`\delta \rho (𝒙)`$ by $`\delta \rho (𝒙)\rho (𝒙)\rho (𝒙)`$. Furthermore we can derive the explicit formula for the PDF of $`\rho `$ as

$$p(\nu )d\nu =\frac{(1+\nu \sqrt{2N^1})^{N/21}}{(2N^1)^{(N1)/2}\mathrm{\Gamma }(N/2)}\mathrm{exp}\left(\frac{N}{2}\left(1+\sqrt{\frac{2}{N}}\nu \right)\right)d\nu ,$$
(A3)
where we have defined the normalized overdensity $`\nu `$ by $`\nu \delta \rho /\delta \rho ^2^{1/2}`$ and $`\mathrm{\Gamma }()`$ is the Gamma function. As the number $`N`$ increases, the PDF becomes closer to the Gaussian distribution. This seems reasonable considering the central limit theorem. For smaller degrees of freedom $`N`$, the PDFs show strong non-Gaussianity.

The skewness parameter $`S`$ is a fundamental quantity to characterize the non-Gaussian profile. For $`\chi _N^2`$-model we have the following expression from equation (A2) for moments

$$S\frac{\delta \rho (𝒙)^3}{\delta \rho (𝒙)^2^{3/2}}=\sqrt{\frac{8}{N}}.$$
(A4)
Note that using the skewness parameter $`S`$, we can specify the $`\chi _N^2`$-distribution that is completely characterized by the degrees of freedom $`N`$.

The simplicities of relations given so far are mainly due to local nature of the density field $`\rho (𝒙)`$. Every statistical characters of the field are reduced to simple Gaussian statistics of the basic fields $`\varphi _i(𝒙)`$. When we discuss smoothing effects on the density field $`\rho (𝒙)`$ that is an nonlinear combination of Gaussian fields $`\varphi _i`$, the situation becomes highly complicated. First we define the smoothed density field $`\rho _R(𝒙)`$ by

$$\rho _R(𝒙)=d𝒙^{}\rho (𝒙^{})w(𝒙^{}𝒙:R),$$
(A5)
where $`w(𝒙;R)`$ is a filter function with smoothing radius $`R`$. In this appendix we mainly use the Gaussian filter

$$w(𝒙:R)=\frac{1}{\sqrt{(2\pi R^2)^3}}\mathrm{exp}(\frac{𝒙^2}{2R^2}).$$

We expand the fields $`\varphi _i(𝒙)`$ in terms of the Fourier modes $`\varphi _i(𝒌)`$

$$\varphi _i(𝒙)=\frac{1}{(2\pi )^3}𝑑𝒌\varphi _i(𝒌)\mathrm{exp}(i𝒌𝒙).$$
(A6)
As the basic fields $`\varphi _i`$ obey independent random-Gaussian distributions, their statistical nature is completely determined by their power spectrum $`P_\varphi (k)`$ defined by

$$\varphi _i(𝒌)\varphi _j(𝒍)=P_\varphi (k)\delta _{ij}\delta _D(𝒌+𝒍),$$
(A7)
where $`\delta _{ij}`$ is the Kronecker’s delta, $`\delta _D()`$ is the Dirac’s delta function and we have assumed that fluctuations are isotropic. Similarly the smoothed density field $`\rho _R(𝒙)`$ is written by $`\varphi _i(𝒌)`$ as

$$\rho _R(𝒙)=\underset{i=1}{\overset{N}{}}\frac{1}{(2\pi )^6}𝑑𝒌𝑑𝒍\varphi _i(𝒌)\varphi _i(𝒍)\mathrm{exp}(i(𝒌+𝒍)𝒙)W(|𝒌+𝒍|R),$$
(A8)
where $`W(kR)`$ is the Fourier transformed filter function of $`w(𝒙:R)`$. For the Gaussian filter we have $`W(kR)=\mathrm{exp}(k^2R^2/2)`$. With expression (A8) we can write down the second- and third-order moments for the smoothed density field in terms of the power spectrum $`P_\varphi (k)`$. After some algebra we arrive at

$`\delta \rho _R^2`$ $`=`$ $`{\displaystyle \frac{2N}{(2\pi )^6}}{\displaystyle 𝑑𝒌𝑑𝒍P_\varphi (k)P_\varphi (l)W(|𝒌+𝒍|R)^2},`$ (A9)
$`\delta \rho _R^3`$ $`=`$ $`{\displaystyle \frac{8N}{(2\pi )^9}}{\displaystyle 𝑑𝒌𝑑𝒍𝑑𝒎P_\varphi (k)P_\varphi (l)P_\varphi (m)W(|𝒌+𝒍|R)W(|𝒍+𝒎|R)W(|𝒎+𝒌|R)}.`$ (A10)
Using these equations we evaluate the skewness parameter $`S_R`$ for the smoothed density field. We can easily confirm that the unsmoothed value $`S=\sqrt{8/N}`$ is recovered by setting $`R=0`$ in these equations. The smoothing effect changes the skewness parameter $`S_R`$ from the unsmoothed value $`S=\sqrt{8/N}`$. To discuss quantitative effects of smoothing we define a parameter $`f_R`$ defined by

$$S_R\frac{\delta \rho _R(𝒙)^3}{\delta \rho _R(𝒙)^2^{3/2}}=\sqrt{\frac{8}{f_RN}}.$$
(A11)
Note that the parameter $`f_R`$ does not depend on the degrees of freedom $`N`$. The effective degrees of freedom $`N^{}`$ determined form the smoothed skewness parameter $`S_R`$ becomes $`f_R`$ times original value.

Next let us evaluate the smoothed skewness $`S_R`$ numerically and the scaling parameter $`f_R`$ for a sequence of matter fluctuations. We adopt the following from for the power spectrum of the Gaussian fields $`\varphi _i`$

$$P_\varphi (k)=Ak^{n_\varphi }\mathrm{exp}(k^2r^2).$$
(A12)
As the normalization factor $`A`$ is irrelevant for our linear analysis we simply put $`A=1`$. Due to the Gaussian cut-off $`\mathrm{exp}(k^2r^2)`$, small scale power at wave-number $`k\stackrel{>}{}\text{ }r^1`$ is strongly suppressed and the fields $`\varphi _i(𝒙)`$ and $`\rho (𝒙)`$ are very smoothed at spatial scale smaller than $`r`$. This means that the smoothed density field $`\rho _R(𝒙)`$ is close to the unsmoothed field $`\rho (𝒙)`$ for smoothing radius $`R`$ with $`R\stackrel{<}{}\text{ }r`$. Therefore radius $`r`$ represents the spatial scale where the PDF of the smoothed field $`\rho _R`$ would become close to the original $`\chi _N^2`$ distribution and the scaling factor would be $`f_R1`$.

Using equation (A9) for variance of the smoothed density field, we obtain the following result for the Gaussian filter

$$\delta \rho _R^2=\frac{NR^{62n_\varphi }}{4\pi ^4}_0^{\mathrm{}}𝑑s_0^{\mathrm{}}𝑑t(st)^{n_\varphi +1}\mathrm{exp}\left[\left(1+\frac{r^2}{R^2}\right)(s^2+t^2)\right]\mathrm{sinh}(2st).$$
(A13)
At large smoothing radius $`R`$ we have $`(\delta \rho )^2R^{62n_\varphi }`$ and the matter power spectrum at small wave-number would behave as (Peebles 1999b, White 1999, Scoccimarro 2000)

$$P_\rho (k)k^{2n_\varphi +3}=k^{n_\rho },$$
(A14)
where we have defined the spectral index $`n_\rho `$ of density field by $`n_\rho 2n_\varphi +3`$.

The expression (A10) for the third-order moment is a nine-dimensional integral. But using symmetry with respect to variables $`𝒌,𝒍`$ and $`𝒎`$, we can simplify the expression. After some algebra we arrive at

$`\delta \rho _R^3`$ $`=`$ $`{\displaystyle \frac{NR^{93n_\varphi }}{4\pi ^6}}{\displaystyle _0^{\mathrm{}}}𝑑s{\displaystyle _0^{\mathrm{}}}𝑑t{\displaystyle _0^{\mathrm{}}}𝑑u{\displaystyle _0^\pi }𝑑a\mathrm{sin}a{\displaystyle _0^\pi }𝑑b\mathrm{sin}b(stu)^{2+n_\varphi }I_0(tu\mathrm{sin}a\mathrm{sin}b)`$ (A15)
$`\times \mathrm{exp}\left[\left(1+{\displaystyle \frac{r^2}{R^2}}\right)(s^2+t^2+u^2)\left(st\mathrm{cos}a+tu\mathrm{cos}a\mathrm{cos}b+us\mathrm{cos}b\right)\right],`$
where $`I_0(x)`$ is the $`0`$-th modified Bessel function. From equations (A13) and (A15), it is apparent that the factor $`f_R`$ is written in the form

$$f_R=f(R/r).$$
(A16)
In the large scale limit $`R\mathrm{}`$, the small scale Gaussian cut-off in the power spectrum $`P_\varphi (k)`$ (eq.\[A12\]) would have no effects on the second- and third-order moments as expected from equations (A13) and (A15). In this limit we can put $`r=0`$ and our results should coincide with results for pure power-law models $`P_\varphi (k)k^{n_\varphi }`$ (Peebles 1999b, Scoccimarro 2000). Simple analytic expressions for the matter power spectrum $`P_\rho `$ and bispectrum $`B_\rho `$ for $`n_\varphi =2`$ model are given by Scoccimarro (2000) as follows

$$P_\rho (k)=\frac{2\pi ^3N}{k},B_\rho (𝒌,𝒍,𝒎)=\frac{8\pi ^3N}{klm}.$$
(A17)
With the latter expression we can numerically study the third-order moment $`(\delta \rho _R)^3`$ at $`R\mathrm{}`$. In this case we need only two-dimensional numerical integration that is much simpler than the five-dimensional integration (A15). We find the numerical result 333For the top-hat filter we have similar result $`f(\mathrm{})=2.6`$.

$$f(\mathrm{})=2.6(n_\varphi =2).$$
(A18)
This asymptotic value can be used to check our five-dimensional numerical integration required to evaluate the smoothed skewness at an intermediate scale.

In figure 4 we show scale dependence of the factor $`f(R/r)`$ for power spectrum (A12) with various spectrum indices, $`n_\varphi =1`$, $`n_\varphi =1.5`$, $`n_\varphi =2`$ and $`n_\varphi =2.4`$. These indices correspond to $`n_\rho =1`$, $`n_\rho =0`$, $`n_\rho =1`$, and $`n_\rho =1.8`$, respectively. The last choice $`n_\varphi =2.4`$ is same as Peebles (1999b). As expected, we have $`f(R/r)1`$ for $`R\stackrel{<}{}\text{ }r`$ reflecting the original profile of the unsmoothed PDF. For a larger smoothing radius $`R\stackrel{>}{}\text{ }20r`$ the parameter $`f(R/r)`$ relaxes to an asymptotic value $`f(\mathrm{})`$ that does not depend on the small scale cut-off as explained earlier. The asymptotic value for $`n_\varphi =2`$ model is close to the value given in equation (A18) and show the validity of our numerical integration. Peebles (1999b) calculated the skewness parameter smoothed by the top-hat filter for a model with one degree of freedom $`N=1`$ and spectral index at $`n_\varphi =2.4`$, using Monte Carlo integration. His result is $`S_R=2.46`$ and corresponds to $`f_R1.32`$. This value is close to our asymptotic value $`f(\mathrm{})1.37`$ obtained for the Gaussian filter.

As shown in figure 4, the factor $`f(R/r)`$ largely depends on the spectral indices $`n_\varphi `$. We have $`f(\mathrm{})200`$ for $`n_\varphi =1(n_\rho =1)`$ model. This large value indicates that the PDF would become close to the Gaussian profile. The factor keeps $`f(R/r)\stackrel{<}{}\text{ }3`$ for $`n_\varphi =2(n_\rho =1)`$ model. White (1999) investigated this model using N-body simulations and found that the PDF is skewed for Gaussian smoothing spanning a factor 5 in scale. This behavior is consistent with our result. For $`n_\varphi =2.4`$ model, the factor $`f(R/r)`$ changes only $`40\%`$ in the range $`0R<\mathrm{}`$. These dependence on the spectral indices $`n_\varphi `$ seems reasonable considering that as $`n_\varphi `$ increases, number of statistically independent region would also increases for a given smoothing volume and the PDF would become a more Gaussian like profile.


format: text

text:
# Application and explicit solution of recurrence relations with respect to space-time dimension

## 1 Introduction

The increasing precision of current experiments and the expected precision of future experiments reveal a need to evaluate the two-loop and higher order radiative corrections . An evaluation of mass dependent two-loop corrections for many physical parameters such as the $`\rho `$ parameter  and $`Z\overline{b}b`$ decay width  is greatly desirable. Of special interest is the problem of calculating one-loop corrections to the process $`e^+e^{}4`$ fermions. To solve this problem efficient algorithms for evaluating multi leg one-loop tensor integrals are needed. Some progress was recently achieved in this direction in .

In the present paper we consider the algorithms based on recurrence relations w.r.t. space-time dimension or $`d`$-shift recurrence relations which we hope will be useful for evaluating multi loop and multi leg integrals. In the first part we describe the method for the reduction of tensor integrals to a set of basic integrals proposed in . In the second part we present a solution of recurrence relations w.r.t. $`d`$ for one-loop $`n`$-point scalar integrals with arbitrary momenta and masses. This is a short report on work done in collaboration with J.Fleischer and F.Jegerlehner .

## 2 Reduction of tensor integrals

The reduction of one-loop tensor integrals to scalar integrals with shifted dimensions was given in . The reduction of multi loop tensor integrals was considered in . In the approach of  any tensor $`L`$ loop integral is represented as

$`{\displaystyle \underset{i=1}{\overset{L}{}}}{\displaystyle \frac{d^dk_i}{\pi ^{dL/2}}\underset{j=1}{\overset{N}{}}P_{\overline{k}_j,m_j}^{\nu _j}\underset{r=1}{\overset{n_1}{}}k_{1\mu _r}\mathrm{}\underset{s=1}{\overset{n_L}{}}k_{L\lambda _s}}=`$
$`T_{\mu _1\mathrm{}\lambda _{n_L}}(\{p_i\},\{_j\},𝐝^+){\displaystyle \underset{i=1}{\overset{L}{}}}{\displaystyle \frac{d^dk_i}{\pi ^{dL/2}}\underset{j=1}{\overset{N}{}}P_{\overline{k}_j,m_j}^{\nu _j}},`$ (1)
where $`\overline{k}_j`$ is the momentum of the $`j`$-th line, $`p_i`$ are external momenta, $`T_{\mu _1\mathrm{}\lambda _{n_L}}`$ is a tensor operator and

$`_j{\displaystyle \frac{}{m_j^2}},𝐝^+G^{(d)}=G^{(d+2)},`$
$`P_{k,m}={\displaystyle \frac{1}{k^2m^2+iϵ}}.`$
It is assumed that the scalar integral on the r.h.s. of (1) has arbitrary masses and only after differentiation w.r.t. $`m_i^2`$ are these set to their concrete values.

To derive the explicit formula for the operator $`T`$, $`L`$ independent auxiliary vectors $`a_i`$ were introduced to represent the tensor in the integrand as

$`k_{1\mu _1}\mathrm{}k_{L\lambda _{n_L}}={\displaystyle \frac{1}{i^{n_1+\mathrm{}+n_L}}}{\displaystyle \frac{}{a_{1\mu _1}}}\mathrm{}{\displaystyle \frac{}{a_{L\lambda _{n_L}}}}`$
$`\times \mathrm{exp}\left[i(a_1k_1+\mathrm{}+a_Lk_L)\right]|_{a_i=0},`$
and after that the integral was transformed into an $`\alpha `$\- parametric representation. From the parametric representation it is easy to deduce that

$`T_{\mu _1\mathrm{}\lambda _{n_L}}(\{p_i\},\{_j\},𝐝^+)={\displaystyle \frac{e^{iQ(\{\{p\},\alpha ,\{0\})\rho }}{i^{n_1+\mathrm{}+n_N}}}`$
$`\times {\displaystyle \frac{}{a_{1\mu _1}}}\mathrm{}{\displaystyle \frac{e^{i[Q(\{p\},\alpha ,\{a\})]\rho }}{a_{L\lambda _{n_L}}}}|_{\genfrac{}{}{0pt}{}{a_j=0}{\genfrac{}{}{0pt}{}{\alpha _j=i_j}{\rho =i^L𝐝^+}}},`$
where $`Q`$ is a polynomial in $`\alpha ,a`$ and $`p_ip_j`$. Applying $`T`$ gives a sum made of external momenta, $`g_{\mu \nu }`$ multiplied by scalar integrals with shifted $`d`$. Advantages of the proposed tensor reduction are:

* no contractions with external momenta and metric tensor and no solution of a linear system of equations are needed
* it is easy to select a particular tensor structure
* representation in terms of integrals with shifted $`d`$ is compact

Thus the problem is reduced to the evaluation of scalar integrals with different powers of scalar propagators and shifts in $`d`$. They can be evaluated by applying generalized recurrence relations . To obtain these relations we start with the identity:

$$\underset{i=1}{\overset{L}{}}d^dk_i\frac{}{k_{r\mu }}\left\{R_{\{\mu \}}(\{k\},\{p\})\underset{j=1}{\overset{N}{}}P_{\overline{k}_j,m_j}^{\nu _j}\right\}0,$$

where $`R`$ is an arbitrary tensor polynomial. Such kind of identities for finding relations between Feynman integrals was used in ,. A systematic method was proposed in . Upon differentiating two representations for scalar products can be used:

$`\mathrm{a})integrationbypartsmethod\text{[13]}.`$
$`k_ip_l={\displaystyle \frac{1}{2}}(k_i^2+p_l^2(k_ip_l)^2),`$
$`\mathrm{b})generalizedrecurrencerelations\text{[8]}:`$
$`p_{l\mu }{\displaystyle k_{i\mu }\underset{r=1}{\overset{L}{}}d^dk_r\underset{j=1}{\overset{N}{}}P_{\overline{k}_j,m_j}^{\nu _j}}=`$
$`p_{l\mu }T_\mu (\{p\},\{\},𝐝^+){\displaystyle \underset{r=1}{\overset{L}{}}d^dk_r\underset{j=1}{\overset{N}{}}P_{\overline{k}_j,m_j}^{\nu _j}}.`$
Combining these methods gives many recurrence relations connecting integrals with different indices $`\nu _j`$ and shifts of $`d`$.

### 2.1 One-loop $`n`$-point integrals

The method for the reduction of one-loop tensor integrals

$$I_{n,r}^{(d)}=\frac{d^dq}{i\pi ^{d/2}}\underset{j=1}{\overset{n}{}}\frac{q_{\mu _1}\mathrm{}q_{\mu _r}}{c_j^{\nu _j}},$$

with $`c_j=(qp_j)^2m_j^2+iϵ,`$ was presented in . The integral $`I_{n,r}^{(d)}`$ can be written as a combination of scalar integrals by using the relation

$$I_{n,r}^{(d)}=T_{\mu _1\mathrm{}\mu _r}(\{p_s\},\{_j\},𝐝^+)I_n^{(d)},$$

where

$$I_n^{(d)}=\frac{d^dq}{i\pi ^{d/2}}\underset{j=1}{\overset{n}{}}\frac{1}{c_j^{\nu _j}},$$
(2)
and

$`T_{\mu _1\mathrm{}\mu _r}(\{p_j\},\{_j\},𝐝^+)={\displaystyle \frac{1}{i^r}}`$
$`\times {\displaystyle \underset{j=1}{\overset{r}{}}}{\displaystyle \frac{}{a_{\mu _j}}}\mathrm{exp}[i({\displaystyle \underset{k=1}{\overset{n}{}}}(ap_k)\alpha _k{\displaystyle \frac{a^2}{4}})\rho ]|{}_{\genfrac{}{}{0pt}{}{a_j=0}{\genfrac{}{}{0pt}{}{\alpha _j=i_j}{\rho =i𝐝^+}}}{}^{}.`$
Two recurrence relations are needed to reduce any scalar integral $`I_n^{(d)}`$ to basic ones:

$`2\mathrm{\Delta }_n\nu _j𝐣^+I_n^{(d)}={\displaystyle \underset{k=1}{\overset{n}{}}}(1+\delta _{jk})`$
$`\times \left({\displaystyle \frac{\mathrm{\Delta }_n}{Y_{jk}}}\right)\left[d{\displaystyle \underset{i=1}{\overset{n}{}}}\nu _i(𝐤^{}𝐢^++1)\right]I_n^{(d)},`$ (3)
$`(d{\displaystyle \underset{i=1}{\overset{n}{}}}\nu _i+1)I_n^{(d+2)}=`$
$`\left[{\displaystyle \frac{2\mathrm{\Delta }_n}{G_{n1}}}+{\displaystyle \underset{k=1}{\overset{n}{}}}{\displaystyle \frac{(_k\mathrm{\Delta }_n)}{G_{n1}}}𝐤^{}\right]I_n^{(d)},`$ (4)
where $`𝐣^\pm `$ etc. shift the indices $`\nu _j\nu _j\pm 1`$ and

$`\mathrm{\Delta }_n\mathrm{\Delta }_n(\{p_1,m_1\},\mathrm{}\{p_n,m_n\})=`$
$`\left|\begin{array}{cccc}Y_{11}& Y_{12}& \mathrm{}& Y_{1n}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ Y_{1n}& Y_{2n}& \mathrm{}& Y_{nn}\end{array}\right|,`$ (8)
$`G_{n1}G_{n1}(p_1,\mathrm{},p_n)=2^n`$
$`\times \left|\begin{array}{ccc}p_{1n}p_{1n}& \mathrm{}& p_{1n}p_{n1n}\\ p_{1n}p_{2n}& \mathrm{}& p_{2n}p_{n1n}\\ \mathrm{}& \mathrm{}& \mathrm{}\\ p_{1n}p_{n1n}& \mathrm{}& p_{n1n}p_{n1n}\end{array}\right|,`$ (13)
$$Y_{ij}=p_{ij}^2+m_i^2+m_j^2,p_{ij}=p_ip_j.$$

Here $`p_i`$ is the external momentum flowing through the $`i`$-th line and $`m_j`$ is the mass of the propagator corresponding to the $`j`$-th line.

Relations similar to (4) for the integrals with zero Gram determinants were given in , .

In the next sections we will also use an index notation for $`\mathrm{\Delta }_n`$ and $`G_n`$

$`\lambda _{i_1i_2\mathrm{}i_n}=\mathrm{\Delta }_n(\{p_{i_1},m_{i_1}\},\mathrm{},\{p_{i_n},m_{i_n}\}),`$
$`g_{i_1\mathrm{}i_n}=G_{n1}(p_{i_1},\mathrm{},p_{i_n}).`$ (14)
We shall use the index notation for integrals obtained from $`I_n^{(d)}`$ by contracting some lines. The characteristic variables occurring in our derivations will be ratios of Gram determinants and therefore it is worthwhile to introduce the shorthand notation

$$r_{ij\mathrm{}k}=\frac{\lambda _{ij\mathrm{}k}}{g_{ij\mathrm{}k}}.$$
(15)
### 2.2  Two-loop propagator integrals

The reduction of the two-loop propagator type tensor integrals to basic scalar integrals can be done by using the relation

$`{\displaystyle \frac{d^dk_1d^dk_2}{c_1^{\nu _1}c_2^{\nu _2}c_3^{\nu _3}c_4^{\nu _4}c_5^{\nu _5}}k_{1\mu _1}\mathrm{}k_{1\mu _r}k_{2\lambda _1}\mathrm{}k_{2\lambda _s}}`$
$`=T_{\mu _1\mathrm{}\lambda _s}(q,\{_j\},𝐝^+){\displaystyle \frac{d^dk_1d^dk_2}{c_1^{\nu _1}c_2^{\nu _2}c_3^{\nu _3}c_4^{\nu _4}c_5^{\nu _5}}},`$ (16)
$`c_1=k_1^2m_1^2+iϵ,c_3=(k_1q)^2m_3^2+iϵ,`$
$`c_2=k_2^2m_2^2+iϵ,c_4=(k_2q)^2m_4^2+iϵ,`$
$`c_5=(k_1k_2)^2m_5^2+iϵ.`$
An explicit form for the operator $`T`$ and all necessary recurrence relations are given in . The recursive procedure allows one to transform any diagram into a sum over 30 basic integrals:

$$I^{(d)}(q^2)=\underset{j=1}{\overset{30}{}}R_j(q^2,\{m_i^2\},d)I_j^{(d)}(q^2).$$

This algorithm was implemented in  as a Mathematica package.

The method of tensor reduction  was used for the evaluation of radiative corrections to several important physical quantities. For example, for the evaluation of the 2-loop correction to the static potential in QCD  and the evaluation of the 2-loop correction in the gauge Higgs system in 3-dimensions . Another important application may be the evaluation of 2-loop correction to the Bhabha scattering , .

## 3 Explicit solution of $`d`$ shift recurrence relations

The solution of the $`d`$ shift relations for one-loop propagator type integral was first discussed in . Here we present the solution of the relation (4) for the integral $`I_n^{(d)}`$ with arbitrary momenta and masses and first powers of scalar propagators. In this case

$$(dn+1)I_n^{(d+2)}=[\frac{2\mathrm{\Delta }_n}{G_{n1}}+\underset{k=1}{\overset{n}{}}\frac{_k\mathrm{\Delta }_n}{G_{n1}}𝐤^{}]I_n^{(d)},$$
(17)
where the operator $`𝐤^{}`$ removes the $`k`$-th line from $`I_n^{(d)}`$. If we assume that $`n1`$ point functions are already known then relation (17) is an inhomogeneous first order difference equation w.r.t. $`d`$. Methods of solution of this kind of equations are well described in the mathematical literature . The redefinition

$$I_n^{(d)}=\frac{1}{\mathrm{\Gamma }\left(\frac{dn+1}{2}\right)}\left(\frac{\mathrm{\Delta }_n}{G_{n1}}\right)^{\frac{d}{2}}\overline{I}_n^{(d)}$$
(18)
leads to a simpler equation

$`\overline{I}_n^{(d+2)}=\overline{I}_n^{(d)}+{\displaystyle \frac{\mathrm{\Gamma }\left(\frac{dn+1}{2}\right)}{2\mathrm{\Delta }_n}}`$
$`\times \left({\displaystyle \frac{G_{n1}}{\mathrm{\Delta }_n}}\right)^{\frac{d}{2}}{\displaystyle \underset{k=1}{\overset{n}{}}}(_k\mathrm{\Delta }_n)𝐤^{}I_n^{(d)}.`$ (19)
Without loss of generality $`d`$ can be parametrized as

$$d=2l2\epsilon ,$$

with $`l`$ being an integer number and $`\epsilon `$ an arbitrary small parameter. The solution of equation (19) then reads

$`\overline{I}_n^{(2l2\epsilon )}=\stackrel{~}{b}_n(\epsilon )+{\displaystyle \underset{r=0}{\overset{l}{}}}{\displaystyle \frac{\mathrm{\Gamma }\left(r1\epsilon \frac{n1}{2}\right)}{2\mathrm{\Delta }_n}}`$
$`\times \left({\displaystyle \frac{G_{n1}}{\mathrm{\Delta }_n}}\right)^{r1\epsilon }{\displaystyle \underset{k=1}{\overset{n}{}}}(_k\mathrm{\Delta }_n)𝐤^{}I_n^{(2r22\epsilon )},`$ (20)
where $`\stackrel{~}{b}_n(\epsilon )`$ is an $`l`$ independent constant. By changing $`\stackrel{~}{b}_n(\epsilon )`$ the solution (20) can be rewritten in another form

$`\overline{I}_n^{(2l2\epsilon )}=\overline{b}_n(\epsilon ){\displaystyle \underset{r=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\mathrm{\Gamma }\left(r+\frac{dn+1}{2}\right)}{2\mathrm{\Delta }_n}}`$
$`\times \left({\displaystyle \frac{G_{n1}}{\mathrm{\Delta }_n}}\right)^{r+\frac{d}{2}}{\displaystyle \underset{k=1}{\overset{n}{}}}(_k\mathrm{\Delta }_n)𝐤^{}I_n^{(d+2r)}.`$ (21)
The result for $`I_n^{(d)}`$ then reads

$`I_n^{(d)}={\displaystyle \frac{1}{\mathrm{\Gamma }(\frac{dn+1}{2})}}({\displaystyle \frac{\mathrm{\Delta }_n}{G_{n1}}})^{\frac{d}{2}}\overline{b}_n{\displaystyle \underset{k=1}{\overset{n}{}}}({\displaystyle \frac{_k\mathrm{\Delta }_n}{2\mathrm{\Delta }_n}})`$
$`\times {\displaystyle \underset{r=0}{\overset{\mathrm{}}{}}}({\displaystyle \frac{dn+1}{2}})_r({\displaystyle \frac{G_{n1}}{\mathrm{\Delta }_n}})^r𝐤^{}I_n^{(d+2r)},`$ (22)
where $`\overline{b}_n`$ can be determined from the asymptotic value of $`I_n^{(d)}`$ in the limit $`d\mathrm{}`$.

### 3.1 Asymptotic form of $`I_n^{(d)}`$ at $`d\mathrm{}`$

The asymptotic form of $`I_n^{(d)}`$ as $`d\mathrm{}`$ can be derived from its parametric integral representation. $`I_n^{(d)}`$ can be written in a parametric form by employing, for example, the following formula

$`{\displaystyle \frac{1}{c_1c_2\mathrm{}c_n}}=\mathrm{\Gamma }(n){\displaystyle _0^1}\mathrm{}{\displaystyle _0^1}𝑑x_1\mathrm{}𝑑x_{n1}`$
$`\times {\displaystyle \frac{x_1^{n2}x_2^{n3}\mathrm{}x_{n2}}{D^n}},`$ (23)
with

$`D=[c_1x_1\mathrm{}x_{n1}`$
$`+c_2x_1\mathrm{}x_{n2}(1x_{n1})+\mathrm{}+c_n(1x_1)],`$ (24)
then shifting integration variable $`q`$ and applying the formula

$$\frac{d^dq}{[i\pi ^{d/2}]}\frac{1}{(q^2m_i^2)^\alpha }=\frac{(1)^\alpha \mathrm{\Gamma }\left(\alpha \frac{d}{2}\right)}{\mathrm{\Gamma }(\alpha )(m_i^2)^{\alpha \frac{d}{2}}}.$$
(25)
Any $`n`$-point function will be represented by a multiple parametric integral of the form:

$`I_n^{(d)}=\mathrm{\Gamma }\left(n{\displaystyle \frac{d}{2}}\right){\displaystyle _0^1}\mathrm{}{\displaystyle _0^1}𝑑x_1\mathrm{}𝑑x_{n1}`$
$`\times f(\{x\})(h_n(\{p_j,m_s\},\{x\}))^{\frac{d}{2}n},`$ (26)
where $`\{x\}=\{x_1,\mathrm{},x_{n1}\}`$ and $`h_n`$ is a polynomial. In our opinion this particular parametric integral representation is very convenient for performing asymptotic expansions. The behavior of this integral, as $`d\mathrm{}`$, can be established by using asymptotic methods . The main contribution to the integral will come from the maximum of $`h_n`$. The maximum can be either on the boundary or in the interior of the integration region. In case the solution $`\{x\}=\{\overline{x}\}`$ of the system of equations

$$\frac{h_n}{x_i}=0,(1<i<n1)$$
(27)
lie in the interior of integration region then it is possible to show that

$$\overline{x}_i=\frac{_{k=1}^{ni}_k\mathrm{\Delta }_n}{_{j=1}^{ni+1}_j\mathrm{\Delta }_n},$$
(28)
and at this point

$$h_n(\{p_j,m_s\},\{\overline{x}\})=r_{1\mathrm{}n}.$$
(29)
In the domain of analyticity $`h_n(\{p_j,m_s\},\{x\})0`$ and $`I_n^{(d)}`$ is an integral of Laplace type. In the domain of non-analyticity there are subdomains of the integration region where $`h_n(\{p_j,m_s\},\{\overline{x}\})<0`$ and so the integral $`I_n^{(d)}`$ has an imaginary part. For this kinematic configuration $`I_n^{(d)}`$ may be represented as

$`I_n^{(d)}=\mathrm{\Gamma }(n{\displaystyle \frac{d}{2}})\{{\displaystyle }{\displaystyle \frac{\{dx\}\theta (h_n)f(\{x\})}{(h_n(\{p_j,m_s\},\{x\}))^{n\frac{d}{2}}}}`$
$`+\left(\mathrm{cos}{\displaystyle \frac{\pi }{2}}(d2n)+i\mathrm{sin}{\displaystyle \frac{\pi }{2}}(d2n)\right)`$
$`\times {\displaystyle \underset{j}{}}{\displaystyle _{\mathrm{\Omega }_j^{}}}{\displaystyle \frac{\{dx\}f(\{x\})}{\left|h_n(\{p_j,m_s\},\{x\})\right|^{n\frac{d}{2}}}}\}.`$ (30)
Here $`\mathrm{\Omega }_j^{}`$ are subdomains of the integration region where $`h_n<0`$. The boundaries of these subdomains can be found from the solutions of the equation

$$h_n(\{p_j,m_s\},\{\overline{x}\})=0.$$
(31)
In (30) we took the main branch of $`h_n`$ and assumed that the momenta and masses are real parameters. A concrete example of representation (30) will be given in the next section.

The imaginary part of $`I_n^{(d)}`$

$`\mathrm{Im}I_n^{(d)}={\displaystyle \frac{\pi }{\mathrm{\Gamma }\left(\frac{d}{2}n+1\right)}}`$
$`\times {\displaystyle \underset{j}{}}{\displaystyle _{\mathrm{\Omega }_j^{}}}{\displaystyle \frac{\{dx\}f(\{x\})}{\left|h_n(\{p_j,m_s\},\{x\})\right|^{n\frac{d}{2}}}}`$ (32)
is a sum of integrals of Laplace type and therefore its asymptotic form at $`d\mathrm{}`$ can be established.

It should be noted that representations like (30) are valid for multi loop integrals including vacuum ones and also they can be used as an alternative to the method for evaluating imaginary part of the diagram with the help of Cutkosky’s rules.

The integrals (26) are multiple Laplace type integrals of the form

$$F(\lambda )=_\mathrm{\Omega }f(x)\mathrm{exp}[\lambda S(x)]𝑑x,$$
(33)
where $`\mathrm{\Omega }`$ is a bounded simply connected domain in $`k`$ dimensional Euclidean space, $`x=(x_1,\mathrm{}x_k)`$, $`\lambda `$ is a parameter, $`S(x)`$ is a real function. Supposing that the absolute maximum of $`S(x)`$ in $`\mathrm{\Omega }`$ is achieved only at the interior point $`x=\overline{x}`$, then at $`\lambda \mathrm{}`$

$$F(\lambda )\mathrm{exp}[\lambda S(\overline{x})]\underset{r=0}{\overset{\mathrm{}}{}}a_k\lambda ^{r\frac{k}{2}}.$$
(34)
This expansion may be differentiated w.r.t. $`\lambda `$ any number of times. The leading term of the expansion is

$$F(\lambda )=\mathrm{exp}[\lambda S(\overline{x})](2\pi /\lambda )^{\frac{k}{2}}\frac{f(\overline{x})+O(\lambda ^1)}{\sqrt{\left|\mathrm{det}S_{xx}(\overline{x})\right|}},$$
(35)
where $`\left|\mathrm{det}S_{xx}(\overline{x})\right|`$ is a Hessian determinant. The leading contribution from the integral at $`d\mathrm{}`$ comes from the region where $`\left|h_n\right|`$ has maximal value. If an extremum of $`h_n`$ occurs inside the integration region then we must compare $`|h_n|`$ at the extremum point (28) with its value on the boundary of the integration region and take the largest one. If an extremum of $`h_n`$ lies outside the integration region then the maximal value of $`\left|h_n\right|`$ will be reached on the boundary. Since $`h_n`$ in the integral is multiplied by the factor $`x_1^{n2}x_2^{n3}\mathrm{}x_{n2},`$ the non-zero contribution to the leading asymptotic comes either from the boundary region $`x_{n1}=0`$ or from the boundary region where one of

$$x_i=1,(1in1).$$

At these $`n`$ boundary edges $`h_n`$ is equal to $`h_{n1}`$ of one of the $`n`$ integrals which one can obtain by contracting one of the lines in $`I_n`$. Similar to $`h_n`$ the maximum value of these $`h_{n1}`$ is equal either to their extremum value

$$h_{n1}=\left|r_{i_1\mathrm{}i_{n1}}\right|$$
(36)
or to the value on the boundary. Again at the boundary $`h_{n1}`$ will be equal to the $`h_{n2}`$ corresponding to integrals one can obtain from $`I_n^{(d)}`$ by contracting two lines and so on. Splitting the integral for $`h_n<0`$ into real and imaginary parts (30) does not influence our considerations because splitting adds a boundary where $`h_n`$ vanishes and therefore there will be no contribution to the asymptotic value of the integral.

According to the previous discussion and from Eq.(22) it follows that in order to find $`\overline{b}_n`$ the contribution of $`I_n^{(d)}`$ at $`d\mathrm{}`$ is needed only if it is proportional to $`r_{1\mathrm{}n}^{d/2}`$. Such an asymptotic value will occur for kinematic configurations when all $`\overline{x}1`$ and $`h_n`$ has either maximum or negative minimum. In the latter case $`I_n^{(d)}`$ must be split into real and imaginary parts according to (30) and the asymptotic value of the term giving a contribution proportional to $`r_{1\mathrm{}n}^{d/2}`$ are to be determined. It should be noted that the contribution of order $`r_{1\mathrm{}n}^{d/2}`$ may come from the analytic continuation of the sums in (22). As will be seen in the next sections these sums represent hypergeometric functions.

From the formula for the leading asymptotic of the multiple integral (35) we obtain:

$$\overline{b}_n=(2\pi )^{\frac{n}{2}}\frac{\mathrm{\Gamma }\left(1\frac{d}{2}\right)\mathrm{\Gamma }\left(\frac{d}{2}\right)}{\sqrt{\pi |G_{n1}|}}\frac{i^d}{r_n^{\frac{n+1}{2}}},$$
(37)
and therefore the general solution of Eq.(17) for the kinematic configuration satisfying the conditions $`0\overline{x}_i1`$ will be

$`I_n^{(d)}=(2\pi )^{\frac{n}{2}}{\displaystyle \frac{\mathrm{\Gamma }(1\frac{d}{2})\mathrm{\Gamma }(\frac{d}{2})}{\mathrm{\Gamma }\left(\frac{dn+1}{2}\right)}}{\displaystyle \frac{r_n^{\frac{dn1}{2}}}{\sqrt{\pi |G_{n1}|}}}`$
$`{\displaystyle \underset{k=1}{\overset{n}{}}}({\displaystyle \frac{_k\mathrm{\Delta }_n}{2\mathrm{\Delta }_n}}){\displaystyle \underset{r=0}{\overset{\mathrm{}}{}}}({\displaystyle \frac{dn+1}{2}})_r`$
$`\times ({\displaystyle \frac{G_{n1}}{\mathrm{\Delta }_n}})^r𝐤^{}I_n^{(d+2r)},`$ (38)
where $`(a)_r\mathrm{\Gamma }(r+a)/\mathrm{\Gamma }(a)`$ is the Pochhammer symbol. A representation of one-loop integrals in terms of Lauricella functions was given in .

## 4 2-point function

Expression (22) for $`I_2^{(d)}`$ includes two one-fold infinite sums over tadpole integrals $`I_1^{(d)}`$ (see formula (25)):

$$I_1^{(d)}=\mathrm{\Gamma }\left(1\frac{d}{2}\right)(m_i^2)^{\frac{d}{2}1}.$$
(39)
It is convenient to label the lines of $`I_2^{(d)}`$ by $`i,j`$ because the integrals $`I_2^{(d)}`$ will be encountered in calculating $`I_3^{(d)}`$, … as a result of contraction of different lines.

At $`n=2`$, substituting (39) into (22) gives

$`{\displaystyle \frac{2\lambda _{ij}I_2^{(d)}}{\mathrm{\Gamma }\left(1\frac{d}{2}\right)}}=b_2+{\displaystyle \frac{_i\lambda _{ij}}{(m_j^2)^{1\frac{d}{2}}}}{\displaystyle \underset{r=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(\frac{d1}{2})_r}{\left(\frac{d}{2}\right)_r}}\left({\displaystyle \frac{m_j^2}{r_{ij}}}\right)^r`$
$`+{\displaystyle \frac{_j\lambda _{ij}}{(m_i^2)^{1\frac{d}{2}}}}{\displaystyle \underset{r=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\left(\frac{d1}{2}\right)_r}{\left(\frac{d}{2}\right)_r}}({\displaystyle \frac{m_i^2}{r_{ij}}})^r.`$ (40)
One can easily recognize that each of the two infinite sums can be written in terms of Gauss’ hypergeometric function $`{}_{2}{}^{}F_{1}^{}`$. The parametric integral for $`I_2^{(2)}`$ is

$$I_2^{(d)}=\mathrm{\Gamma }\left(2\frac{d}{2}\right)_0^1𝑑x_1h_2^{\frac{d}{2}2},$$
(41)
with

$$h_2=p_{ij}^2x_1^2x_1(p_{ij}^2m_i^2+m_j^2)+m_j^2.$$

An extremum of $`h_2`$ is at $`\overline{x}_1=(p_{ij}^2m_i^2+m_j^2)/(2p_{ij}^2).`$ The maximum of $`h_2`$ exists if

$$\frac{^2h_2}{x_1^2}=2p_{ij}^2<0,$$
(42)
i.e. only for Euclidean momentum.

The formula (32) gives the imaginary part of $`I_2^{(d)}`$ on the cut

$$\mathrm{Im}I_2^{(d)}=\pi _{x_1^{}}^{x_1^+}𝑑x_1\frac{((x_1x_1^{})(x_1^+x_1))^{\frac{d}{2}2}}{p_{ij}^{4d}\mathrm{\Gamma }\left(\frac{d2}{2}\right)},$$
(43)
where

$$x_1^\pm =\frac{p_{ij}^2m_i^2+m_j^2\pm \sqrt{\lambda _{ij}}}{2p_{ij}^2}$$

is the solution of the equation

$$h_2(\{p_{ij},m_i\},x_1)=0.$$
(44)
The change of the integration variable $`x_1=(x_1^+x_1^{})z+x_1^{}`$ makes the integration trivial and gives

$$\mathrm{Im}I_2^{(d)}=\frac{\pi \mathrm{\Gamma }\left(\frac{d2}{2}\right)}{p_{ij}^{4d}\mathrm{\Gamma }\left(d2\right)}\left(\frac{\sqrt{\lambda _{ij}}}{p_{ij}^2}\right)^{d3}.$$
(45)
The real part of $`I_2^{(d)}`$ coming from the region where $`h_2<0`$ is

$$\mathrm{Re}I_2^{(d)}=\mathrm{Im}I_2^{(d)}\mathrm{ctg}\frac{\pi }{2}(d4).$$

At large $`d`$ the asymptotic value of $`\mathrm{Im}I_2^{(d)}`$ and $`\mathrm{Re}I_2^{(d)}`$ can be easily found. At $`p_{ij}^2>(m_1+m_2)^2`$ both $`0x_1^\pm 1`$.

Now several remarks concerning the evaluation of $`b_2`$ are in order. From (22) we see that at large dimension $`b_2r_{ij}^{d/2}`$. Such a contribution may come either from the asymptotic value of $`I_2^{(d)}`$ (this may happen only if $`0<\overline{x}_1<1`$) at large $`d`$ or from the analytic continuation of infinite sums when their expansion parameter exceeds 1. To proceed further let us assume that $`m_i^2>m_j^2`$. At $`p_{ij}^2<m_j^2m_i^2`$ and $`p_{ij}^2>(m_i+m_j)^2`$ the value of $`b_2`$ is determined from the asymptotic of $`I_2^{(d)}r_{ij}^{d/2}`$. At $`m_j^2m_i^2<p_{ij}^2<(m_im_j)^2`$ both $`I_2^{(d)}`$ and the infinite sums have no contributions $`r_{ij}^{d/2}`$ and therefore in this region $`b_2=0`$. In the region $`(m_im_j)^2<p_{ij}^2<(m_1+m_2)^2`$ the value of $`b_2`$ is determined from the asymptotic value of infinite sums at large $`d`$.

Combining all values of $`b_2`$ and writing infinite sums in terms of hypergeometric functions we obtained the following result

$`{\displaystyle \frac{2\lambda _{ij}I_2^{(d)}}{\mathrm{\Gamma }\left(1\frac{d}{2}\right)}}=`$
$`{\displaystyle \frac{\sqrt{\pi }\mathrm{\Gamma }\left(\frac{d}{2}\right)}{\mathrm{\Gamma }\left(\frac{d1}{2}\right)}}r_{ij}^{\frac{d2}{2}}\left[{\displaystyle \frac{_i\lambda _{ij}}{\sqrt{1\frac{m_j^2}{r_{ij}}}}}+{\displaystyle \frac{_j\lambda _{ij}}{\sqrt{1\frac{m_i^2}{r_{ij}}}}}\right]`$
$`+{\displaystyle \frac{_i\lambda _{ij}}{(m_j^2)^{1\frac{d}{2}}}}{}_{2}{}^{}F_{1}^{}\left[\begin{array}{c}1,\frac{d1}{2};\\ \frac{d}{2};\end{array}{\displaystyle \frac{m_j^2}{r_{ij}}}\right]`$ (48)
$`+{\displaystyle \frac{_j\lambda _{ij}}{(m_i^2)^{1\frac{d}{2}}}}{}_{2}{}^{}F_{1}^{}\left[\begin{array}{c}1,\frac{d1}{2};\\ \frac{d}{2};\end{array}{\displaystyle \frac{m_i^2}{r_{ij}}}\right].`$ (51)
The analytic continuation of $`{}_{2}{}^{}F_{1}^{}`$ leads to a compact form for $`I_2^{(d)}`$ (see also )

$`{\displaystyle \frac{g_{ij}I_2^{(d)}}{\mathrm{\Gamma }(2\frac{d}{2})}}={\displaystyle \frac{_i\lambda _{ij}}{(m_j^2)^{\frac{4d}{2}}}}{}_{2}{}^{}F_{1}^{}\left[\begin{array}{c}1,\frac{4d}{2};\\ \frac{3}{2};\end{array}1{\displaystyle \frac{r_{ij}}{m_j^2}}\right]`$ (54)
$`+{\displaystyle \frac{_j\lambda _{ij}}{(m_i^2)^{\frac{4d}{2}}}}{}_{2}{}^{}F_{1}^{}\left[\begin{array}{c}1,\frac{4d}{2};\\ \frac{3}{2};\end{array}1{\displaystyle \frac{r_{ij}}{m_i^2}}\right].`$ (57)
The expansion of $`I_2^{(d)}`$ to all orders in $`\epsilon =(4d)/2`$ is given in .

## 5  3-point function

The expression for $`I_3^{(d)}`$ includes summation over $`{}_{2}{}^{}F_{1}^{}`$ coming from $`I_2^{(d)}`$. To derive an explicit formula we first use the relation

$$(1z)^a{}_{2}{}^{}F_{1}^{}\left[\begin{array}{c}a,b;\\ c;\end{array}z\right]={}_{2}{}^{}F_{1}^{}\left[\begin{array}{c}a,cb;\\ c;\end{array}\frac{z}{z1}\right],$$
(58)
in order to remove the factor $`d`$ from the second parameter of $`{}_{2}{}^{}F_{1}^{}`$. Summing over $`I_2^{(d)}`$ gives the Appell hypergeometric functions
$`F_3((d2)/2,1,1,1/2,d/2;x,y)`$ defined as

$`F_3(\alpha ,\alpha ^{},\beta ,\beta ^{},\gamma ;x,y)=`$
$`{\displaystyle \underset{m,n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(\alpha )_m(\alpha ^{})_n(\beta )_m(\beta ^{})_n}{(\gamma )_{m+n}m!n!}}x^my^n,`$ (59)
which may be reduced to $`F_1`$ by means of

$`F_3(\alpha ,\alpha ^{},\beta ,\beta ^{},\alpha +\alpha ^{};x,y)=`$
$`(1y)^\beta ^{}F_1(\alpha ,\beta ,\beta ^{},\alpha +\alpha ^{};x,{\displaystyle \frac{y}{y1}}).`$ (60)
The Appell hypergeometric function $`F_1`$ defined as

$`F_1(\alpha ,\beta ,\beta ^{},\gamma ;x,y)=`$
$`{\displaystyle \underset{m,n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(\alpha )_{m+n}(\beta )_m(\beta ^{})_n}{(\gamma )_{m+n}m!n!}}x^my^n.`$ (61)
Similar to $`I_2^{(d)}`$ we will label the lines of $`I_3^{(d)}`$ by $`i,j,k`$. The result for $`I_3^{(d)}`$ valid when $`h_3`$ has maximum inside the integration region reads

$`{\displaystyle \frac{\lambda _{ijk}}{\mathrm{\Gamma }\left(2\frac{d}{2}\right)}}I_3^{(d)}=2^{\frac{3}{2}}\pi \sqrt{g_{ijk}}r_{ijk}^{\frac{d2}{2}}`$
$`+\theta _{ijk}_k\lambda _{ijk}+\theta _{kij}_j\lambda _{ijk}+\theta _{jki}_i\lambda _{ijk},`$ (62)
where

$`\lambda _{ij}\theta _{ijk}=`$
$`{\displaystyle \frac{(m_i^2)^{\frac{d2}{2}}}{2(d2)}}{\displaystyle \frac{_j\lambda _{ij}}{\sqrt{1\frac{m_i^2}{r_{ij}}}}}F_1({\displaystyle \frac{d2}{2}},1,{\displaystyle \frac{1}{2}},{\displaystyle \frac{d}{2}};{\displaystyle \frac{m_i^2}{r_{ijk}}},{\displaystyle \frac{m_i^2}{r_{ij}}})`$
$`+{\displaystyle \frac{(m_j^2)^{\frac{d2}{2}}}{2(d2)}}{\displaystyle \frac{_i\lambda _{ij}}{\sqrt{1\frac{m_j^2}{r_{ij}}}}}F_1({\displaystyle \frac{d2}{2}},1,{\displaystyle \frac{1}{2}},{\displaystyle \frac{d}{2}};{\displaystyle \frac{m_j^2}{r_{ijk}}},{\displaystyle \frac{m_j^2}{r_{ij}}})`$
$`\left[{\displaystyle \frac{_i\lambda _{ij}}{\sqrt{1\frac{m_j^2}{r_{ij}}}}}+{\displaystyle \frac{_j\lambda _{ij}}{\sqrt{1\frac{m_i^2}{r_{ij}}}}}\right]r_{ij}^{\frac{d2}{2}}{\displaystyle \frac{\sqrt{\pi }\mathrm{\Gamma }(\frac{d2}{2})}{4\mathrm{\Gamma }\left(\frac{d1}{2}\right)}}`$
$`\times {}_{2}{}^{}F_{1}^{}\left[\begin{array}{c}1,\frac{d2}{2};\\ \frac{d1}{2};\end{array}{\displaystyle \frac{r_{ij}}{r_{ijk}}}\right].`$ (65)
The function $`F_1`$ has a simple integral representation

$`F_1({\displaystyle \frac{d2}{2}},1,{\displaystyle \frac{1}{2}},{\displaystyle \frac{d}{2}};x,y)={\displaystyle \frac{d2}{2}}`$
$`\times {\displaystyle _0^1}{\displaystyle \frac{u^{\frac{d4}{2}}}{(1xu)\sqrt{1yu}}}du.`$ (66)
The polynomial

$`h_3=x_1x_2(1x_1)p_{13}^2x_1^2x_2(1x_2)p_{12}^2`$
$`x_1x_2m_1^2x_1(1x_1)(1x_2)p_{23}^2`$
$`+x_1(1x_2)m_2^2+(1x_1)m_3^2.`$ (67)
has a maximum at

$$x_1=\frac{_1\mathrm{\Delta }_3+_2\mathrm{\Delta }_3}{G_2},x_2=\frac{_1\mathrm{\Delta }_3}{_1\mathrm{\Delta }_3+_2\mathrm{\Delta }_3}$$

if

$`{\displaystyle \frac{^2h_3}{x_2x_2}}=2x_1^2p_{12}^2<0,`$ (68)
$`{\displaystyle \frac{^2h_3}{x_1x_1}}{\displaystyle \frac{^2h_3}{x_2x_2}}\left({\displaystyle \frac{^2h_3}{x_1x_2}}\right)^2`$
$`={\displaystyle \frac{(_1\mathrm{\Delta }_3+_2\mathrm{\Delta }_3)^2}{2G_2}}>0.`$ (69)
From (68),(5) it follows that the maximum will be achieved if

$$p_{12}^2<0,\mathrm{and}G_2<0.$$
(70)
The behavior of the integrand for $`I_3^{(d)}`$ having the maximum inside the integration region is demonstrated in Fig.1.

Fig. 1: Function $`x_1h_3^{d/23}`$ at $`p_{13}^2=10,p_{23}^2=5,p_{12}^2=9`$, $`m_i^2=1`$, $`d=36`$.

To evaluate $`I_3^{(d)}`$ in different kinematic regions an analytic continuation of $`{}_{2}{}^{}F_{1}^{}`$ and $`F_1`$ is needed. The corresponding formulae are given in . The representations of $`I_3^{(d)}`$ in terms of other hypergeometric functions can be found in .

### 5.1  4-point function

The most problematic part of the computation of $`I_4^{(d)}`$ is the summation of terms from $`I_3^{(d)}`$ with $`F_1`$ function. To simplify the summation we remove the factor $`d`$ in the first argument of $`F_1`$ by applying the formula

$`(1x)^\beta (1y)^\beta ^{}F_1(\alpha ,\beta ,\beta ^{},\gamma ;x,y)=`$
$`F_1(\gamma \alpha ,\beta ,\beta ^{},\gamma ;{\displaystyle \frac{x}{x1}},{\displaystyle \frac{y}{y1}}).`$ (71)
The result for the $`4`$-point function then reads

$`{\displaystyle \frac{\lambda _{ijkl}}{\mathrm{\Gamma }\left(2\frac{d}{2}\right)}}I_4^{(d)}=4\pi ^{\frac{3}{2}}\sqrt{g_{ijkl}}r_{ijkl}^{\frac{d3}{2}}{\displaystyle \frac{\mathrm{\Gamma }\left(\frac{d2}{2}\right)}{\mathrm{\Gamma }\left(\frac{d3}{2}\right)}}`$
$`+\varphi _{ijkl}_l\lambda _{ijkl}+\varphi _{lijk}_k\lambda _{ijkl}`$
$`+\varphi _{klij}_j\lambda _{ijkl}+\varphi _{jkli}_i\lambda _{ijkl},`$ (72)
where

$`\varphi _{ijkl}={\displaystyle \frac{\pi \sqrt{2}}{\sqrt{g_{ijk}}}}r_{ijkl}^{\frac{d4}{2}}{}_{2}{}^{}F_{1}^{}\left[\begin{array}{c}1,\frac{d3}{2};\\ \frac{d2}{2};\end{array}{\displaystyle \frac{r_{ijk}}{r_{ijkl}}}\right]`$ (75)
$`+\mathrm{terms}\mathrm{with}F_1\mathrm{and}F_S`$ (76)
and $`F_S`$ is the Lauricella-Saran function  of three variables. In our case it has a simple integral representation:

$`F_S={\displaystyle \frac{\mathrm{\Gamma }\left(\frac{d}{2}\right)(yz)^{\frac{1}{2}}}{\mathrm{\Gamma }\left(\frac{d3}{2}\right)\mathrm{\Gamma }\left(\frac{3}{2}\right)}}`$
$`\times {\displaystyle _0^1}{\displaystyle \frac{\mathrm{arcsin}\sqrt{\frac{(yz)t}{1tz}}(1t)^{\frac{d5}{2}}}{(1x+tx)\sqrt{1ty}}}dt.`$ (77)
All other details one can find in .

## 6 Conclusions

The recurrence relations with respect to the space-time dimension turn out to be useful for the reduction of tensor integrals. The relations for master integrals in different dimensions are rather simple even for integrals with several momenta and masses. Preliminary investigation of recurrence relations w.r.t. $`d`$ for the two-loop propagator type integrals revealed that scalar integral in (16) with all $`\nu _j=1`$ can be written as multiple hypergeometric series. Asymptotic expansion at large $`d`$ may be used as yet another tool for approximate evaluation of Feynman integrals. As is known such an expansion works well in quantum mechanics and in quantum field theory in lattice calculations.

Acknowledgements. I wish to thank organizers of the conference “Loops and Legs 2000” Johannes Blümlein and Tord Riemann for the useful and well organized conference. I am thankful to C. Ford and V. Ravindran for carefully reading the manuscript and useful remarks. I would like also to thank the DFG for financial support.