text:
# Non-Gaussian Tails of Cosmological Density Distribution Function from Dark Halo Approach

## 1 Introduction

The probability distribution function (PDF) of the cosmic matter density is one of most fundamental quantities characterizing statistical properties of large-scale structures. In most standard scenarios of the structure formation, the large-scale structures are thought to form by gravitational instabilities from small Gaussian density fluctuations. A deviation in the PDF from the Gaussian distribution is thus most likely caused by the nonlinear gravitational evolution and is a key to understand the structure formation history. The higher order moments of the density field such like the skewness and kurtosis have been used to study the non-Gaussian nature in cosmic density field (e.g., Peebles 1980; Juszkiewicz, Bouchet & Colombi 1993; Bernardeau 1994a). The one-point PDF contains more statistical information than them.

The galaxy PDFs have been measured from various catalogs (e.g., Hamilton 1985; Bouchet et al. 1993; Kofman et al. 1994) and shows significant non-Gaussianity. Indeed it was pointed out that the galaxy PDF is well fitted by the log-normal distribution (e.g., Coles & Jones 1991; Kofman et al. 1994). The dark matter PDFs measured from $`N`$-body simulations are also fitted by the log-normal distribution very well (e.g., Bernardeau & Kofman 1995; Taylor & Watts 2000). Although the galaxy distribution cannot be regarded as a direct tracer of the dark matter distribution because of unknown biasing relation between them, the similarity in their PDF suggests that there exists a nature common to their distributions.

The origin of the phenomenological log-normal PDF has not yet been clarified. Bernardeau (1994b) and Benardeau & Kofman (1995) argued that the perturbative construction of one-point PDF shows a strong dependence on the primordial power spectrum and the PDF approaches the log-normal form only when the initial spectrum is proportional to $`P(k)k^n`$ with $`n=1`$. On the other hand, using a series of high-resolution N-body simulations, Kayo, Taruya & Suto (2001) found that the one-point density PDFs with Gaussian initial conditions are approximately described by the log-normal distribution in the strongly non-linear regime. Remarkably, the accuracy of log-normal model turns out to become almost insensitive to the underlying power spectrum of density fluctuation, in contrast to the prediction from perturbation theory. While the results of Kayo et al. (2001) show a clear dependence on the primordial initial spectrum in a weakly non-linear regime, consistent with the perturbation theory, there remains crucial to clarify the reason why the PDFs in the simulations tend to become insensitive to the initial power spectrum and approach to the log-normal form in a strongly non-linear regime.

The primary purpose of this paper is to discuss this issue using a simple analytic model. That is, employing the so-called dark halo approach as advocated recently, we construct a model for the non-Gaussian tails of the one-point PDF, which includes the essential ingredient for gravitational clustering and is, at least, consistent with N-body simulations. The basic assumption of our model is that the high density regions after smoothed over with a filter function are mainly attributed to a single massive halo. Clearly, this assumption is violated at the low density regions and we thus restrict our analysis to the highly non-Gaussian tails of mass density PDF. Then, the dark halo approach allows us to derive a simple analytic expression for the one-point PDF in cases with the scale-free initial spectra and power-law profiles of dark halo, from which we explain the qualitative behavior of tails of PDF and its dependences on the halo density profile and on the initial power spectrum. Further, adopting the realistic halo profile and the mass function, the results of our analytic model are quantitatively compared with N-body simulations with various initial power spectra.

In section 2, the basic ingredients for the dark halo approach and an analytic model of tails of one-point PDF are described. Section 3 presents the detailed investigation of the dependence on the halo profile and the sensitivity to the initial power spectra on the non-Gaussian tails of one-point PDF. Assuming the scale-free initial spectrum and the power-law profile of dark halos, we analytically derive an approximate expression for non-Gaussian tails of PDF. Using this result, we attempt to clarify the sensitivity of the one-point PDF to the initial spectra and the halo profiles. In section 4, the analytic model of one-point PDF is quantitatively compared with those obtained from cosmological N-body simulations. Finally, section 5 is devoted to the conclusions and the discussion.

## 2 Mass density PDF from the dark halo approach

### 2.1 Basic ingredients

The basic idea of the dark halo approach which has become very popular recently(e.g., Seljak 2000; Ma & Fry 2000a; Scoccimarro et al. 2001) is that (i) all the dark matter in the universe belong to one of the virialized clumps (dark halos), and (ii) the matter distribution in the universe is approximated by a sum of the matter distribution within each halo. To construct an analytic model of mass density PDF from the dark halo approach, we need to specify the halo mass function and the density profiles of halos. In the following, we describe these two ingredients.

For the halo mass function, the approximate prescription for counting the number of dark halos has been first presented by Press & Schechter (1974) and recently an accurate analytic fitting model has been proposed by Sheth & Tormen (1998). Both models are simply expressed in a unified manner as follows:

$$n(M)dM=A\left\{1+\left(\frac{\sigma _\mathrm{M}}{\sqrt{a}\delta _{\mathrm{crit}}}\right)^{2p}\right\}\sqrt{\frac{a}{2\pi }}\frac{\overline{\rho }}{M}\frac{\delta _{\mathrm{crit}}}{\sigma _\mathrm{M}^3}\mathrm{exp}\left\{\frac{a}{2}\left(\frac{\delta _{\mathrm{crit}}}{\sigma _\mathrm{M}}\right)^2\right\}\left|\frac{d\sigma _\mathrm{M}^2}{dM}\right|dM,$$
(1)
where $`\overline{\rho }`$ is the mean mass density and $`\delta _{\mathrm{crit}}`$ is the critical over-density of the spherical collapse model(see Nakamura & Suto 1997 for useful fitting formula). The quantity $`\sigma _\mathrm{M}`$ denotes the linear rms fluctuations of the matter density in a sphere of radius $`R_M=(3M/4\pi \overline{\rho })^{1/3}`$. The numerical coefficients $`a`$ and $`p`$ are set to $`a=1`$, $`p=0`$ for the Press–Schechter mass function and $`a=0.707`$, $`p=0.3`$ for the mass function by Sheth & Tormen. The normalization factor $`A`$ is determined by $`n(M)M𝑑M=\overline{\rho }`$. Throughout the paper, we extensively use the mass function by Sheth & Tormen (1998).

The density profile of dark halos is defined by a spherical average for each halo and is given by a function of radius and halo mass. The numerical simulation by Navarro et al. (1996,1997; hereafter NFW) revealed that the outer slope of profile becomes $`3`$, while the inner part asymptotically approaches $`1`$. In this paper, we adopt their fitting form:

$$\frac{\rho _{\mathrm{NFW}}(r;M)}{\overline{\rho }}=\frac{\delta _c}{(r/r_s)(1+r/r_s)^2},$$
(2)
where $`\delta _c`$ and $`r_s`$ respectively denote the characteristic over-density and the scale radius. These parameters are not independent and are characterized in terms of concentration parameter $`cr_{\mathrm{vir}}/r_s`$, with $`r_{\mathrm{vir}}`$ being virial radius of the halo. Using this quantity, the characteristic over-density $`\delta _c`$ is given by $`\delta _c=(\mathrm{\Delta }_{\mathrm{vir}}/3)c^3/[\mathrm{ln}(1+c)c/(1+c)]`$, where $`\mathrm{\Delta }_{\mathrm{vir}}`$ denotes the virial over-density of spherical collapse model (Nakamura & Suto 1997). For the relation between the concentration parameter and the halo mass, we adopt the following fitting form according to the simulation by Bullock et al. (2001),

$$c(M)=c_{}\left(\frac{M}{M_{}}\right)^\beta ,$$
(3)
where $`M_{}`$ corresponds to the characteristic mass scale satisfying $`\sigma _\mathrm{M}(M_{})=\delta _{\mathrm{crit}}`$. Below in comparing the PDF with the N-body simulation (Sec.4), the parameters $`c_{}`$ and $`\beta `$ are chosen as $`(c_{},\beta )=(10,0.15)`$ for cold dark matter(CDM) model and $`(10,(n+3)/6)`$ for scale-free model with spectral index $`n`$. The systematic influences of the uncertainty of $`c(M)`$ on the evaluation of PDF are discussed in detail in Appendix B.

Indeed, recent high-resolution N-body study reveals that the inner slope of halo profile might not converge to $`1`$, but it rather approaches $`1.5`$ (Moore et al. 1999; Fukushige & Makino 2001a,b; see also Jing & Suto 2000). Hence, it becomes important to clarify the profile dependence on the non-Gaussian tails of PDF. In this paper, in addition to the NFW profile, we further treat the halo models characterized by a single power-law profile:

$$\frac{\rho _{\mathrm{power}}(r;M)}{\overline{\rho }}=\delta _c\left(\frac{r}{r_s}\right)^\alpha ,(0<\alpha <3)$$
(4)
which enables us to derive the analytic expressions of PDFs (see Sec.3). Here, the characteristic density and scale radius are determined by the relation between virial mass and virial radius:

$$\delta _cr_s^\alpha =\frac{3\alpha }{3}\mathrm{\Delta }_{\mathrm{vir}}r_{\mathrm{vir}}^\alpha =\frac{3\alpha }{3}\mathrm{\Delta }_{\mathrm{vir}}\left(\frac{3M}{4\pi \overline{\rho }\mathrm{\Delta }_{\mathrm{vir}}}\right)^{\alpha /3}.$$
(5)
### 2.2 Mass density PDFs

Provided the basic ingredients, we now describe a simple analytic model of mass density PDF.

Let us first describe the cumulative distribution function $`P(>\delta ;R)`$, which gives a probability that a given density field smoothed over the radius $`R`$ is larger than $`\delta `$. Defining the smoothed halo profile:

$$\rho _{\mathrm{smooth}}(r;R,M)=d^3𝒓^{}W(|𝒓𝒓^{}|;R)\rho (r^{};M),$$
(6)
with the function $`W(x;R)`$ being the spherically symmetric smoothing kernel. Since we only consider simple spherical symmetric density distributions of halos, one can calculate the characteristic radius $`r(\delta ;R,M)`$ by solving the following equation:

$$\frac{\rho _{\mathrm{smooth}}[r(\delta ;R,M);R,M]}{\overline{\rho }}=\delta +1\delta .$$
(7)
Now we introduce the volume factor $`V(>\delta ;R,M)`$ which represents, for one halo with mass $`M`$, the comoving volume in which the smoothed mass density field is always greater than $`\delta `$. Since the smoothed halo profile is a monotonic decreasing function of the radius, the volume factor is simply given by,

$$V(>\delta ;R,M)=\frac{4\pi }{3}\left\{r(\delta ;R,M)\right\}^3.$$
(8)
The volume factor (8) explicitly depends on the halo mass and thereby the volume fraction of the entire region greater than $`\delta `$ is evaluated from the integration of equation (8) over the different halo mass weighting the comoving number density of dark halos. Indeed, this volume fraction is equivalent to the cumulative distribution function $`P(>\delta ;R)`$. We thus obtain:

$$P(>\delta ;R)=_0^{\mathrm{}}dMV(>\delta ;R,M)n(M).$$
(9)
Once provided the cumulative distribution, it is easy to calculate the one-point PDF from the definition $`P(>\delta ;R)=_\delta ^+\mathrm{}d\delta P(\delta ;R)`$:

$$P(\delta ;R)=\left|\frac{d}{d\delta }P(>\delta ;R)\right|.$$
(10)
Equation (10) with equation (9) is the heart of our subsequent analysis. It provides a simple analytic estimate for the degree of non-Gaussian tails, which is, in principle applicable to both weakly non-linear and strongly non-linear regimes. Note that the expression (9) heavily relies on the assumption that any of the dark matter particles reside in each dark halo. This would be true for high density parts of the mass density PDF, but might not be the case for the low density regions, $`1<\delta \stackrel{<}{}1`$. We shall therefore focus on the high density tails, $`\delta 1`$. As shown in the following sections, this approach is very useful to understand the essential physical mechanisms that determine the behavior of the high density tail of the PDF.

While the above model prescription is general and does not restrict the choice of the smoothing filter, hereafter, we specifically treat the density field smoothed over the top-hat filter function:

$$W_{\mathrm{th}}(𝒓;R_{\mathrm{th}})=\frac{3}{4\pi R_{\mathrm{th}}^3}\mathrm{\Theta }(R_{\mathrm{th}}|𝒓|),$$
(11)
where the quantity $`\mathrm{\Theta }(x)`$ denotes Heaviside step function.

## 3 Dependences of power spectrum index and halo profile

Before addressing the comparison with N-body simulations, it is instructive to consider how the halo density profile and initial power spectrum alter the non-Gaussian tails of PDF in a simple analytical manner. In this section, simply assuming the power-law model of halo density profile (4) and scale-free initial spectra $`P(k)k^n`$, we here derive an approximate expression for the one-point PDF $`P(\delta ;R_{\mathrm{th}})`$.

To compute the PDFs, one must first evaluate the volume factor (8), which is given in terms of the characteristic radius $`r(\delta ;R_{\mathrm{th}},M)`$. In Appendix A, the exact expressions for top-hat smoothed density profile is presented in the case of power-law profiles (eqs.). The resultant smoothed profiles become constant at an inner region $`rR_{\mathrm{th}}`$ and show the power-law behavior at $`rR_{\mathrm{th}}`$ (see left panel of Fig.6). Based on these exact results, in Appendix A, we further obtain the approximate expression for the smoothed density profiles, which enable us to calculate the quantity $`r(\delta ;R_{\mathrm{th}},M)`$ analytically (see eq.):

$$\frac{\rho _{\mathrm{smooth}}(r;R_{\mathrm{th}},M)}{\overline{\rho }}\{\begin{array}{ccc}\delta _0(R_{\mathrm{th}};M)\hfill & ;& \hfill r<\left(\frac{3\alpha }{3}\right)^{1/\alpha }R_{\mathrm{th}}\\ & & \\ \left(\frac{3\alpha }{3}\right)\left(\frac{r}{R_{\mathrm{th}}}\right)^\alpha \delta _0(R_{\mathrm{th}};M)\hfill & ;& \hfill r\left(\frac{3\alpha }{3}\right)^{1/\alpha }R_{\mathrm{th}}\end{array},$$
(12)
where the function $`\delta _0(R_{\mathrm{th}},M)`$ represents the maximum value of the smoothed density profile, $`\rho _{\mathrm{smooth}}(0;R_{\mathrm{th}},M)/\overline{\rho }`$, and is given by a function of halo mass:

$$\delta _0(R_{\mathrm{th}},M)\mathrm{\Delta }_{\mathrm{vir}}\left(\frac{R_{\mathrm{th}}}{r_{\mathrm{vir}}}\right)^\alpha =\mathrm{\Delta }_{\mathrm{vir}}\left(\frac{M_R}{M}\right)^{\alpha /3},$$
(13)
The quantity $`M_R`$ means the effective halo mass inside the smoothing radius $`R_{\mathrm{th}}`$, $`M_R(4\pi /3)\overline{\rho }\mathrm{\Delta }_{\mathrm{vir}}R_{\mathrm{th}}^3`$. For a given smoothing length $`R_{\mathrm{th}}`$, equation (13) relates the mass of a halo with its characteristic smoothed density. Adopting the above approximation (12), one can easily solve (7), which yields

$$r(\delta ;R_{\mathrm{th}},M)=\{\begin{array}{ccc}0\hfill & ;& \hfill \delta >\delta _0(R_{\mathrm{th}};M)\\ & & \\ \left\{\frac{3\alpha }{3}\frac{\delta _0(R_{\mathrm{th}},M)}{\delta }\right\}^{1/\alpha }R_{\mathrm{th}}\hfill & ;& \hfill \delta \delta _0(R_{\mathrm{th}};M)\end{array}.$$
(14)
Then, with a help of (13), substitution of equation (14) into (8) leads to

$$V(>\delta ;R_{\mathrm{th}},M)=\{\begin{array}{ccc}0\hfill & ;& \hfill \delta >\delta _0(R_{\mathrm{th}};M)\\ & & \\ \left(\frac{3\alpha }{3}\frac{\mathrm{\Delta }_{\mathrm{vir}}}{\delta }\right)^{3/\alpha }\frac{M}{\overline{\rho }\mathrm{\Delta }_{\mathrm{vir}}}\hfill & ;& \hfill \delta \delta _0(R_{\mathrm{th}};M)\end{array}.$$
(15)
Clearly, the volume factor decreases as $`\delta ^{3/\alpha }`$ and is cut off at the maximum density $`\delta _0(R_{\mathrm{th}};M)`$. This indicates that for a certain value of $`\delta `$, the dominant contribution to the integral in cumulative distribution $`P(>\delta ;R_{\mathrm{th}})`$ is mainly attributed to the halos whose maximum density is larger than that.

Now, turn to focus on the evaluation of the integral in (9). For the cosmological models with initial scale-free spectra, the rms fluctuation $`\sigma _\mathrm{M}`$ in (1) can be expressed as

$$\sigma _\mathrm{M}=\delta _{\mathrm{crit}}\left(\frac{M}{M_{}}\right)^{(n+3)/6}.$$
(16)
Substituting the mass function (1) and the volume factor (15) into (9), the cumulative distribution function is analytically evaluated as follows:

$`P(>\delta ;R_{\mathrm{th}})`$ $`=`$ $`\left({\displaystyle \frac{3\alpha }{3}}{\displaystyle \frac{\mathrm{\Delta }_{\mathrm{vir}}}{\delta }}\right)^{3/\alpha }{\displaystyle \frac{1}{\overline{\rho }\mathrm{\Delta }_{\mathrm{vir}}}}{\displaystyle _{(\delta /\mathrm{\Delta }_{\mathrm{vir}})^{3/\alpha }M_R}^+\mathrm{}}n(M)M𝑑M`$ (17)
$`=`$ $`\left({\displaystyle \frac{3\alpha }{3}}{\displaystyle \frac{\mathrm{\Delta }_{\mathrm{vir}}}{\delta }}\right)^{3/\alpha }G(\delta ;\alpha ,n),`$
where

$`G(\delta ;\alpha ,n)`$ $``$ $`{\displaystyle \frac{1}{\overline{\rho }\mathrm{\Delta }_{\mathrm{vir}}}}{\displaystyle _{(\delta /\mathrm{\Delta }_{\mathrm{vir}})^{3/\alpha }M_R}^+\mathrm{}}n(M)M𝑑M`$ (18)
$`=`$ $`{\displaystyle \frac{A}{\mathrm{\Delta }_{\mathrm{vir}}\sqrt{\pi }}}{\displaystyle _{y_\mathrm{R}(\delta )/2}^+\mathrm{}}\left\{1+{\displaystyle \frac{1}{(2x)^p}}\right\}x^{1/2}e^x{\displaystyle \frac{dx}{x}}`$
$`=`$ $`{\displaystyle \frac{A}{\mathrm{\Delta }_{\mathrm{vir}}\sqrt{\pi }}}\left\{\mathrm{\Gamma }({\displaystyle \frac{1}{2}},{\displaystyle \frac{y_\mathrm{R}(\delta )}{2}})+2^p\mathrm{\Gamma }({\displaystyle \frac{1}{2}}p,{\displaystyle \frac{y_\mathrm{R}(\delta )}{2}})\right\}`$
with the function $`\mathrm{\Gamma }(a,x)`$ being the incomplete Gamma function, $`\mathrm{\Gamma }(a,x)_x^{\mathrm{}}𝑑tt^{a1}e^t`$. Here, the function $`y_\mathrm{R}(\delta )`$ is defined by

$$y_\mathrm{R}(\delta )a\left(\frac{M_R}{M_{}}\right)^{(n+3)/3}\left(\frac{\delta }{\mathrm{\Delta }_{\mathrm{vir}}}\right)^{(n+3)/\alpha }=a\left(\frac{\delta }{\delta _{}}\right)^{(n+3)/\alpha },$$
(19)
where $`\delta _{}\delta _0(R_{\mathrm{th}},M_{})`$ and denotes the characteristic smoothed density of the halo with mass $`M_{}`$, which play a key role to reveal the sensitivity of the non-Gaussian tails.

Equation (17) indicates that a shape of the PDF is determined by the combination of two contributions: One is the halo profile. The scaling $`\delta ^{3/\alpha }`$ in equation (17) reflects the slope of the halo profile. The flatter the halo profile is, the more rapidly the volume factor varies, thus the steeper the PDF becomes. Another contribution is the shape of the mass function, since the factor $`M/\overline{\rho }\mathrm{\Delta }_{\mathrm{vir}}`$ gives the characteristic volume of the halo with the mass $`M`$, and the lower limit of the integration in equation (18) is, roughly speaking, determined by the halo mass whose characteristic smoothed density (defined by eq.) just coincides with $`\delta `$. Therefore, the function $`G(\delta ;\alpha ,n)`$ gives the volume fraction of halos (i.e., a sum of characteristic volume of halos) whose mass is larger than such lower limit, and directly reflects the shape of the halo mass function. Note that the dependence of the spectral index $`n`$ on the PDF enters only through the function $`G(\delta ;\alpha ,n)`$.

Now from (17), the one-point PDF $`P(\delta ;R_{\mathrm{th}})`$ is analytically expressed as follows:

$`P(\delta ;R_{\mathrm{th}})`$ $`=`$ $`{\displaystyle \frac{3}{\alpha }}\left({\displaystyle \frac{3\alpha }{3}}\right)^{3/\alpha }{\displaystyle \frac{1}{\mathrm{\Delta }_{\mathrm{vir}}}}\left({\displaystyle \frac{\delta }{\mathrm{\Delta }_{\mathrm{vir}}}}\right)^{(3/\alpha +1)}F(\delta ;\alpha ,n),`$ (20)
with

$`F(\delta ;\alpha ,n)`$ $`=`$ $`G(\delta ;\alpha ,n)+{\displaystyle \frac{\alpha }{3}}\mathrm{\Delta }_{\mathrm{vir}}\left|{\displaystyle \frac{d}{d\delta }}G(\delta ;\alpha ,n)\right|`$ (21)
$`=`$ $`{\displaystyle \frac{A}{\mathrm{\Delta }_{\mathrm{vir}}\sqrt{\pi }}}\left[\mathrm{\Gamma }({\displaystyle \frac{1}{2}},{\displaystyle \frac{y_\mathrm{R}(\delta )}{2}})+2^p\mathrm{\Gamma }({\displaystyle \frac{1}{2}}p,{\displaystyle \frac{y_\mathrm{R}(\delta )}{2}})+{\displaystyle \frac{n+3}{3\sqrt{2}}}\{y_\mathrm{R}(\delta )\}^{1/2}\left(1+{\displaystyle \frac{1}{\{y_\mathrm{R}(\delta )\}^p}}\right)e^{y_\mathrm{R}(\delta )/2}\right].`$
The expression (20) with (19) and (21) is a main analytical result in this paper. Qualitatively, it reveals the strong dependence of the halo density profiles and the sensitivity to the initial spectra.

Figure 1 shows the behavior of the function $`F(\delta ;\alpha ,n)`$ for various values of $`n`$ (left panel) and for $`\alpha `$ (right panel). Note again that the dependence of the power index $`n`$ on the PDF enters only through the function $`F(\delta ;\alpha ,n)`$. It is clear from both panels that the slope of the function $`F`$ becomes sensitive to both $`n`$ and $`\alpha `$ in the range $`\delta /\delta _{}>1`$, corresponding to the mass range $`MM_{}`$. The sensitivity to the power index $`n`$ simply results from the steep exponential tails of the halo mass function, which becomes very sensitive to the spectral index in such mass range. And the sensitivity to the slope of the halo profile $`\alpha `$ is also explained from the exponential tails of mass function, that is, the variation in the relation between characteristic smoothed density $`\delta _0`$ (eq.) and halo mass $`M`$ due to a small change of $`\alpha `$ is significantly enhanced through such exponential dependences. On the other hand, the slope of the function $`F`$, in turn, becomes insensitive to both $`n`$ and $`\alpha `$ at the region $`\delta /\delta _{}<1`$. This is because in such range (corresponding to the mass range of $`M<M_{}`$) the mass function only weakly depends on the spectral index $`n`$ and scales as $`n(M)M^{(n9)/6}`$. Therefore, within the range $`\delta /\delta _{}<1`$, the PDF is almost solely determined by the contribution from the halo profile and varies as $`\delta ^{(3/\alpha +1)}`$, while at the region $`\delta /\delta _{}>1`$, the function $`F`$ varies much faster than $`\delta ^{(3/\alpha +1)}`$, thus the shape of PDF is determined by the contribution from the exponential tails of mass function.

Figure 2 shows the approximate expression of one-point PDF (20) for various slope of the halo profiles (left) and for various initial spectra (right). For definiteness, we assume the Einstein-de Sitter universe (density parameter $`\mathrm{\Omega }_0=1`$, cosmological constant $`\lambda _0=0`$) and the non-linear mass $`M_{}`$ in each model is determined by introducing the artificial length scale $`L_{\mathrm{box}}`$ and normalizing the rms fluctuation $`\sigma _\mathrm{M}`$ to unity at $`R_M=0.1L_{\mathrm{box}}`$. The smoothing radii are specifically chosen as both the large $`\delta _{}`$ ($`3090`$) ($`R_{\mathrm{th}}=0.02L_{\mathrm{box}}`$) and the small $`\delta `$ ($`0.31.5`$) ($`R_{\mathrm{th}}=0.15L_{\mathrm{box}}`$). Table 1 summarizes the values of $`\delta _{}`$ for the models shown in Figure 2.

In left panel, the spectral index are fixed to $`n=1`$ but the slope of halo profile $`\alpha `$ varies from $`1`$ to $`2.5`$, while in right panel, $`\alpha `$ is fixed to $`2`$ but $`n`$ varies from $`2`$ to $`+1`$. It is clear from the cases of $`R_{\mathrm{th}}=0.02L_{\mathrm{box}}`$ in both panels that at the region $`\delta \delta _{}`$, the slope of PDF hardly depends on $`n`$ but sensitively depends on $`\alpha `$. On the other hand, at $`\delta \delta _{}`$ (see the cases of $`R_{\mathrm{th}}=0.02L_{\mathrm{box}}`$ in the right panel), the shape of PDF strongly depends on $`n`$. It follows from these results that if one takes smaller smoothing scale than the virial radius of the halo mass $`M_{}`$ (thus $`\delta _{}100`$), the slope of the PDF almost solely determined by the slope of the halo profile at the intermediate range of $`5<\delta <50`$, while the shape of PDF strongly depends both on $`\alpha `$ and $`n`$ if we conversely take a larger smoothing scale (so that $`\delta _{}10`$). Figure 3 shows this point clearly. In the figure, the PDF in the case of $`n=1`$ and $`\alpha =2`$ is shown restricting the contribution of halo mass to some small mass ranges. It turns out that at $`\delta >\delta _{}`$, the major contribution to the PDF at a certain value of $`\delta `$ comes from halos whose characteristic mass scale $`\delta _0`$ roughly coincides with $`\delta `$, thus the shape of PDF reflects shape of the halo mass function. At $`\delta <\delta _{}`$, however, halos with $`\delta _0=\delta `$ no longer lead to a major contribution, but halos whose characteristic density $`\delta _0`$ is larger than $`\delta `$ instead contribute to the PDF. In this case, PDF reflects the profile of the halos.

## 4 Comparison with N-body simulations

We are now in a position to compare the one-point PDFs from the dark halo approach with those obtained from N-body simulations. To investigate this issue, the analytic results in the previous section turn out to become inadequate because of the non power-law nature of the realistic halo profiles. Hence, all the model predictions presented below are based on the NFW profile (2). In computing the analytic PDF (10), the top-hat smoothed halo $`\rho _{\mathrm{smooth}}`$ is numerically evaluated in the Fourier space, in which we use the following fitting function $`y(k;M)`$ for the Fourier transform of the NFW profile $`\rho _{\mathrm{NFW}}/\overline{\rho }`$ (Ma & Fry 2000):

$$y(k;M)=\frac{4\pi }{3}\delta _cr_s^3\frac{\mathrm{ln}(e+1/q)\mathrm{ln}[\mathrm{ln}(e+1/q)]}{\left(1+q^{1.1}\right)^{(2/1.1)}},$$
(22)
with $`q=kr_s`$.

For the present purpose, we specifically treat the N-body simulation data for the scale-free models with initial power spectra $`P(k)k^n`$ $`(n=1,0,1,2)`$(Jing 1998) and the CDM model with cosmological constant(Lambda CDM, hereafter LCDM; Jing & Suto 1998). All the models employ $`N=256^3`$ dark matter particles in a periodic comoving cube, $`L_{\mathrm{box}}^3`$, where the box size of the LCDM model is specifically chosen as $`L_{\mathrm{box}}=100h^1`$Mpc. The gravitational force calculation is based on the P<sup>3</sup>M algorithm. While the scale-free models assume an Einstein-de Sitter universe, cosmological parameters of the LCDM model are chosen as $`(\mathrm{\Omega }_0,\lambda _0,h,\sigma _8)=(0.3,0.7,0.7,1.0)`$, where the normalization $`\sigma _8`$ means the linear rms fluctuation at $`R_{\mathrm{th}}=8h^1`$Mpc. As for the scale-free models, the normalization of the density fluctuation is determined by setting the linear rms fluctuation to unity at $`R_{\mathrm{th}}=0.1L_{\mathrm{box}}`$.

Figure 4 shows the N-body results of the one-point PDF for the top-hat smoothed density field in various smoothing radii (crosses). The left panel shows the results in LCDM model, while the right panel plots the PDFs in the scale-free model with the initial spectrum $`n=1`$. The error-bars in the left panel indicate the 1-$`\sigma `$ errors among three different realizations. The smoothing radii in left(right) panel are chosen as $`R_{\mathrm{th}}=2,6,18h^1`$Mpc ($`R_{\mathrm{th}}=0.02,0.05,0.15L_{\mathrm{box}}`$),which are typically scaled from non-linear to weakly non-linear regime. In Table 2, computed values of the rms amplitude $`\sigma `$ are presented(labeled by $`\sigma _{\mathrm{sim}}`$). Obviously, the degree of non-linearity becomes significant as decreasing $`R_{\mathrm{th}}`$ and the resultant PDFs show much longer non-Gaussian tails.

In Figure 4, the one-point PDFs from the dark halo approach are shown in solid lines. Within the validity range $`\delta \stackrel{>}{}1`$, in both panels, the dark halo approach reproduces the non-Gaussian tails of simulated PDF reasonably well. The predictions at large smoothing radius $`R_{\mathrm{th}}=18h^1`$Mpc($`0.15L_{\mathrm{box}}`$) somewhat over-predict the simulations, however, as have been pointed out by Kayo et al.(2001), this is most likely due to the finite volume effect of the N-body simulations.

For comparison, the dashed lines in Figure 4 show the empirical log-normal distribution. The one-point PDF of the log-normal model is analytically expressed as

$$P_{\mathrm{LN}}(\delta )=\frac{1}{\sqrt{2\pi \sigma _{\mathrm{LN}}^2}}\mathrm{exp}\left\{\frac{[\mathrm{log}(1+\delta )+\sigma _{\mathrm{LN}}^2/2]^2}{2\sigma _{\mathrm{LN}}^2}\right\}\frac{1}{1+\delta },$$
(23)
where the single parameter $`\sigma _{\mathrm{LN}}`$ is related to the rms amplitude of density fluctuation at a given smoothing radius $`R_{\mathrm{th}}`$, $`\sigma (R_{\mathrm{th}})`$:

$$\sigma _{\mathrm{LN}}^2=\mathrm{log}[1+\sigma ^2(R_{\mathrm{th}})].$$
(24)
For a given set of cosmological parameters, the quantity $`\sigma (R_{\mathrm{th}})`$ is computed using a fitting formula of non-linear power spectrum by Peacock & Dodds (1996) or even by the dark halo approach (Seljak 2000; Ma & Fry 2000; Scoccimarro et al. 2001; Hamana, Yoshida & Suto 2001). The resultant values are summarized in Table 2 (labeled by $`\sigma _{\mathrm{PD}}`$) and adopting those, the log-normal predictions (23) are plotted. In this sense, both empirical and analytic models in Figure 4 provide a fairly self-consistent prediction without invoking any information of simulation data.

As mentioned in Kayo et al. (2001), the rms amplitude $`\sigma _{\mathrm{PD}}`$ takes a slightly larger value than that of the simulation, $`\sigma _{\mathrm{sim}}`$, which should be ascribed to the finite box-size of the simulations. Apart from this small discrepancy, the log-normal PDFs approximate the simulated PDF fairly accurately in both low-density and high-density regions. The most striking evidence is that the tails of analytic PDF almost resembles those of the log-normal PDF, which is also consistent with the prediction from perturbation theory(Bernardeau 1994b; Bernardeau & Kofman 1995). According to the Table 2, the theoretical estimate of characteristic smoothed density $`\delta _{}`$ roughly varies from $`\delta _{}0.03`$ to $`40`$ as decreasing the smoothing radius. Thus, the tails of PDF at smoothing radius $`R_{\mathrm{th}}=2h^1`$Mpc or $`0.02L_{\mathrm{box}}`$ basically reflect the shape of the halo profiles and they scale as $`\delta ^{(\alpha /3+1)}`$. Comparing the shape of PDFs with the slope of the dependence $`\delta ^{(\alpha /3+1)}`$ superimposed in the right panel of Figure 4, the non-Gaussian tails at smaller smoothing radii can be regarded as a consequence of outer slope of the NFW profile, $`\alpha =3`$.

Figure 5 shows weak and strong sensitivity of the non-Gaussian tails to the initial spectrum in scale-free models. Specifically choosing the smoothing radius as $`R_{\mathrm{th}}=0.15L_{\mathrm{box}}`$(left) and $`R_{\mathrm{th}}=0.02L_{\mathrm{box}}`$(right), the results from the dark halo approach are plotted in right-side of each panel, while the left-side shows the simulation results. For comparison, the log-normal PDFs are also plotted in right-side of the right panel(thin-lines), where rms amplitude $`\sigma `$ is calculated from simulations. In cases with large smoothing radius $`R_{\mathrm{th}}=0.15L_{\mathrm{box}}`$, the simulated PDFs rapidly fall off and they eventually vanish, while the analytic model of PDF shows a slightly longer non-Gaussian tail. This discrepancy is basically ascribed to the finite volume size of the simulations (Kayo et al. 2001), however, the tendency seen in both panels is qualitatively similar. The analytic model reproduces the systematic behavior for the variation of spectral index $`n`$ (see also the right-panel of Fig.2). On the other hand, for smaller smoothing radius with $`R_{\mathrm{th}}=0.02L_{\mathrm{box}}`$, except for the noisy data of scale-free model with $`n=+1`$, both left- and right-side of the panel show not only qualitatively but also quantitatively similar dependence on the initial power spectra. Again, the asymptotic behavior of the PDF is quite similar to the results with steeper power-law profile $`\alpha 3`$ (compare with the slope $`\delta ^{(\alpha /3+1)}`$ superimposed in the figure). As a consequence, the non-Gaussian tails of PDF seem to resemble the log-normal PDF insensitive to the initial spectrum, as has been reported by Kayo et al. (2001). At a closer look at tails of PDF, however, the simulated PDFs still exhibit the initial spectrum dependence and they slightly deviate from log-normal PDFs (thin-lines). In contrast, the analytic PDF from the dark halo approach successfully reproduces this weak dependence. Note that the theoretical estimate of the characteristic smoothed density yields $`\delta _{}10`$$`70`$ at $`R_{\mathrm{th}}=0.02L_{\mathrm{box}}`$, which suggests that the the simulated PDF will show a strong sensitivity to the spectral index at the high density tail $`\delta /\delta _{}1`$, and thereby a large discrepancy between log-normal prediction and simulation is expected. Interestingly, there is some indication at the region $`\delta \stackrel{>}{}200`$ in both simulation and halo approach, where the systematic variation for $`n`$ apparently ceases.

Therefore, we conclude that the analytic model based on the dark halo approach reasonably predicts the general tendency in both weakly and strongly non-linear regime and it quantitatively provides a useful approximation to the non-Gaussian tails of simulated PDF in the strongly non-linear regime.

## 5 Discussion & Conclusion

In this paper, utilizing the dark halo approach, we have presented a simple analytic model characterizing the non-Gaussian tails of one-point PDF and investigated the sensitivity to both the initial power spectrum and the halo density profiles in the non-linear regime. Assuming the power-law profile of dark halos and the scale-free initial spectrum, the approximate expression for the tails of PDF is analytically derived (eq. with  ). Using this result, we explained the sensitivity to the initial spectrum seen in the simulation and the significant dependences on the halo profiles. Further, quantitative comparison with N-body simulation was made, which leads to the conclusion that the analytic model reasonably reproduces the simulated PDF qualitatively and quantitatively in a strongly non-linear regime. Then, the tails of simulated PDF turn out to resemble the log-normal distribution (23) at the region $`\delta \stackrel{<}{}100`$ in a strongly non-linear regime.

In our model, the key quantity is the characteristic smoothed density $`\delta _{}=\rho _{\mathrm{smooth}}(0;R,M_{})/\overline{\rho }`$, which discriminates the weak/strong sensitivity to the halo profile and initial power spectrum (see eq.). At $`\delta /\delta _{}\stackrel{<}{}1`$, the tails of PDF basically follow the halo profile and varies as $`\delta ^{(3/\alpha +1)}`$, while the PDF at $`\delta /\delta _{}\stackrel{>}{}1`$ exhibits a strong sensitivity to the initial spectrum, which can be ascribed to the steep exponential tails of halo mass function. In CDM models at present, the characteristic smoothed density typically takes the value $`\delta _{}1`$ for a small smoothing scale $`R_{\mathrm{th}}1`$Mpc$`/h`$ and conversely $`\delta _{}1`$ for a large smoothing scale $`R_{\mathrm{th}}>10`$Mpc$`/h`$. Therefore, the log-normal behavior seen in the non-linear regime can be regarded as a consequence of the universality of the halo density profile. Further considering the fact that the slope of the simulated PDF is similar to that of the analytic results with $`\alpha 3`$, the weak sensitivity to the initial spectrum is especially due to the universal outer profile of $`\rho r^3`$.

The qualitative feature in the dark matter PDFs might be also the case for the PDF of galaxies, which has been known to be well approximated by the log-normal PDF (e.g., Hamilton 1985; Bouchet et al. 1993; Kofman et al. 1994) as the density PDF is. The number density of galaxies within the clusters is known to be fitted well by the King model (King 1962), which has a radial profile of $`\rho _g[1+(r/r_c)^2]^{3/2}`$, where $`r_c`$ is a core radius of order 0.1$`h^1`$Mpc (Wu & Hammer 1993 and references therein). Thus, on suitably larger smoothing scale $`R_{\mathrm{th}}\stackrel{>}{}r_c`$, smoothed density profile becomes $`\rho _gr^3`$, at an outer part. Therefore, it is very likely that the similarity in the shape of the PDF tails between the dark matter and the galaxy is due to simply their similar distribution within halos.

Further, the good agreement in the non-Gaussian tails of density PDF between the halo approach prediction and N-body data indicates that the halo model is not only used to qualitative study but also allows us to make a quantitative prediction. Although the dark matter PDF is hard to measure, it may be applied to related measurements such like the weak lensing convergence which is a line-of-sight projection of the matter density weighted by the lensing efficiency and is obtained from a coherent distortion in distant galaxy images (see Mellier 1999; Bartelmann & Schneider 2000 for reviews). Taruya et al. (2002) showed that the lensing convergence PDF is well approximated by the log-normal model and pointed out that the lensing convergence field can be regarded as a simple projection of the density field. This implies that one can safely apply the halo approach to the lensing convergence PDF. Actually, Kruse & Schneider (2000) employed the halo approach to predict the non-Gaussian tail in PDF of lensing convergence arguing it may be used to discriminate the cosmological model. Their prediction was tested against numerical experiments and a good agreement was found (Reblinsky et al. 1999). Since the approach adopted by Kruse & Schneider (2000) is essentially the same as described in this paper, one naturally expects that the tails of lensing convergence PDF basically follows the outer slope of the halo profiles at smaller angular scale $`\theta 1^{}`$-$`5^{}`$, and the log-normal behavior found by Taruya et al. (2002) can be explained by similar reason to the mass density PDF.

Therefore, as far as the CDM model with a Gaussian initial condition is concerned, theoretical prediction based on the dark halo approach consistently explains the log-normal nature of observable clustering statistics as a common feature in the hierarchical clustering universe. This conversely indicates that the statistics related to the non-Gaussian tails of PDF provide a powerful cosmological probe to the nature of dark matter clustering, which enables us to constrain the number density and the density profile of dark halos. Probably, one promising application is the weak lensing statistics, as has been demonstrated by Kruse & Schneider (2000). Recently, from the $`2.1`$ deg<sup>2</sup> image taken with a wide field camera on the prime focus of Subaru Telescope, non-Gaussian signature of the lensing convergence was detected evaluating the peak distribution function (Miyazaki et al. 2002). The number count of high positive peaks, which are thought to come from the dark halos, was compared with the theoretical prediction using halo approach and showed a consistent result with NFW profiles under the LCDM cosmology. Note, however, that the analytic models of PDF from the dark halo approach currently limit the applicability to highly non-Gaussian tails in a strongly nonlinear regime. To obtain a more stringent constraint on the cosmological parameter and/or the number density of dark halos, we should further improve the applicability to the low-density region up to the observable scale of the non-linear regime. In the light of this, the perturbative construction of PDF by Bernardeau & Valageas (2000) or the non-perturbative modeling from the hierarchical ansatz (Munshi & Jain 2000; Valageas 2000), as well as the empirical log-normal model (Taruya et al.2002) could be useful and be regarded as a complementary approach. Combining the present work with these approaches would be possible to improve the prediction dramatically and will be discussed elsewhere.

## Acknowledgments

We thank Y.P.Jing for kindly providing us his N-body data, Y.Suto for reading of the manuscript and critical comments, M.Takada and K.Yoshikawa for valuable discussion. I.K is supported by Takenaka-Ikueikai Fellowship. T.H acknowledges supports from Japan Society for Promotion of Science (JSPS) Research Fellowships. A.T is supported by the grand-in-aid for Scientific Research of JSPS (No.$`1470157`$).

## Appendix A: Smoothed density profile for power-law halos

In this appendix, adopting the top-hat filter function (11), we derive the exact expression for the smoothed halo profiles. Substituting the definition (11) into (6), the expressions for smoothed density profile can be reduced to

$`\rho _{\mathrm{smooth}}(r;R_{\mathrm{th}},M)`$ $`=`$ $`{\displaystyle \frac{3}{4\pi R_{\mathrm{th}}^3}}{\displaystyle d^3𝒓^{}\mathrm{\Theta }(R_{\mathrm{th}}|𝒓^{}|)\rho (|𝒓𝒓^{}|;M)}`$ (25)
$`=`$ $`{\displaystyle \frac{3}{2R_{\mathrm{th}}^3}}{\displaystyle _0^{R_{\mathrm{th}}}}𝑑ss^2{\displaystyle _1^1}𝑑x\rho (\sqrt{r^2+s^22rsx};M).`$
For the power-law profiles (4), all the integrals in the expression (25) are analytically evaluated, which lead to

$`\rho _{\mathrm{smooth}}(r;R_{\mathrm{th}},M)`$ $`=`$ $`{\displaystyle \frac{3\delta _c\overline{\rho }r_s^\alpha }{2R_{\mathrm{th}}^3}}{\displaystyle _0^{R_{\mathrm{th}}}}𝑑ss^2{\displaystyle \frac{|rs|^{2\alpha }(r+s)^{2\alpha }}{(\alpha 2)rs}}`$ (29)
$`=`$ $`{\displaystyle \frac{3}{2(\alpha 2)}}\rho _{\mathrm{power}}(R_{\mathrm{th}};M)`$
$`\times \{\begin{array}{ccc}\frac{1}{3\alpha }\left\{(1u)^{3\alpha }+(1+u)^{3\alpha }\right\}+\frac{u^1}{4\alpha }\left\{(1u)^{4\alpha }(1+u)^{4\alpha }\right\}\hfill & ;& (0<rR_{\mathrm{th}})\\ & & \\ \frac{1}{3\alpha }\left\{(1+u)^{3\alpha }(u1)^{3\alpha }\right\}+\frac{u^1}{4\alpha }\left\{(u1)^{4\alpha }(1+u)^{4\alpha }\right\}\hfill & ;& (R_{\mathrm{th}}<r)\end{array}`$
Here, the variable $`u`$ denotes $`r/R_{\mathrm{th}}`$. Notice that the above expressions are valid for the profiles with $`\alpha 2`$. For the power-law profile with $`\alpha =2`$ (singular isothermal sphere), the logarithmic dependence of the radius $`r`$ appears, which can be separately evaluated as follows:

$`\rho _{\mathrm{smooth}}(r;R_{\mathrm{th}},M)={\displaystyle \frac{3\delta _c\overline{\rho }r_s^2}{2R_{\mathrm{th}}^3}}{\displaystyle _0^{R_{\mathrm{th}}}}𝑑s{\displaystyle \frac{s}{r}}\mathrm{log}\left|{\displaystyle \frac{r+s}{rs}}\right|={\displaystyle \frac{3}{2}}\rho (R_{\mathrm{th}};M)\left\{1+{\displaystyle \frac{1}{2}}(u^1u)\mathrm{log}\left|{\displaystyle \frac{u+1}{u1}}\right|\right\}.`$ (30)
While the above results are exact, it is convenient and useful to derive the simple analytic expressions for the one-point PDF $`P(\delta ;R_{\mathrm{th}})`$ (see Sec.3). From the asymptotic behavior of the smoothed halo profiles, we approximate the expressions (29) and (30) as follows:

$$\rho _{\mathrm{smooth}}(r;R_{\mathrm{th}},M)\{\begin{array}{ccc}\left(\frac{3}{3\alpha }\right)\rho _{\mathrm{power}}(R_{\mathrm{th}};M)\hfill & ;& \hfill r<\left(\frac{3\alpha }{3}\right)^{1/\alpha }R_{\mathrm{th}}\\ & & \\ \left(\frac{r}{R_{\mathrm{th}}}\right)^\alpha \rho _{\mathrm{power}}(R_{\mathrm{th}};M)\hfill & ;& \hfill r\left(\frac{3\alpha }{3}\right)^{1/\alpha }R_{\mathrm{th}}\end{array}$$
(31)
The left panel of Figure 6 compares the approximation (31) with the exact results (29) and (30). Apart from the intermediate scales $`r/R_{\mathrm{th}}1`$, the approximation recovers the exact results reasonably well. On the other hand, deviation at the intermediate scale slightly increases as increasing the slope of profile $`\alpha `$, which leads to the under-prediction of volume factor $`V(>\delta ;R_{\mathrm{th}},M)`$. Hence, for steeper halo profile $`\alpha 3`$, the analytic PDF using (31) tends to under-predict the amplitude of PDF as shown in the right panel.

To check the validity of the approximation (20), together with the analytic PDF in the left panel of Figure 2, the PDF computed from the exact results without using approximation (12) is plotted and the results are presented in the right panel of Figure 6 (thin lines). As is expected from the left panel, the approximate expression (20) for steeper halo profiles with $`\alpha =2.5`$ tends to lower its amplitude, but the qualitative behavior still remains correct.

## Appendix B: Effects on uncertainty of the model prediction

The analytic model of one-point PDF presented in section 2 suffers from uncertainties in the model assumptions (mass function of dark halos and halo density profile). Among these, the concentration parameter of the NFW halo profiles, $`c(M)`$ might be crucial in comparing the analytic prediction with simulations (Sec.3), since the sensitivity to $`c(M)`$ has been reported in previous study using dark halo approach (e.g., Seljak 2000). Although we adopt the fitting form of the concentration parameter by Bullock et al. (20001) with a particular set of the parameters $`(c_{},\beta )`$ (see eq.  below), we examine the influence of varying the parameters $`(c_{},\beta )`$.

Figure 7 shows the dependence of the non-Gaussian tails on the variation in concentration parameter assuming the LCDM universe. Upper(lower) panels represent the results fixing $`\beta `$($`c_{}`$). Although the influence of the variation $`(c_{},\beta )`$ turns out to be small, the systematic behavior in the tails of PDF seems to vary in a opposite sense: decreasing the parameter $`\beta `$ or increasing $`c_{}`$ suppresses the non-Gaussian tails.

The small influence on the tails of PDF can be easily deduced from the Fourier transform of the NFW profile (22), in which the concentration parameter appears only through the factor $`\delta _cr_s^3`$, leading to the logarithmic dependence of $`c(M)`$. On the other hand, apparently opposite behavior for the variation of $`(c_{},\beta )`$ is explained as follows. First note that decrease of $`\beta `$ or increase of $`c_{}`$ yields more concentrated density field for halos more massive than the nonlinear mass $`M_{}`$. This implies that while the number of halos in which the central density exceeding a specific value of $`\delta `$ increases, the volume of each density profile exceeding $`\delta `$ conversely decreases in the massive halos. Since the contribution of massive halo is relatively dominant in the tails of PDF, the latter effect finally wins and accordingly the amplitude of one-point PDF decreases.

Except for the details, the systematic uncertainty of the concentration parameter does not affect the final conclusions in section 4 and the dispersion of $`0.2`$ in $`\mathrm{ln}c`$ (Jing 2000; Bullock et al. 2001) can be safely neglected in our calculation. Note, however, that in the case of the lensing convergence PDF, the dependence of $`(c_{},\beta )`$ becomes slightly larger than that in the mass density PDF and the systematic behavior of non-Gaussian tails appears in a opposite sense to $`P(\delta )`$. These differences simply result from the projection of the three-dimensional density field to the two-dimensional field (e.g., Takada & Jain 2002).


format: text

text:
# Darboux Transformation and Variable Separation Approach: the Nizhnik-Novikov-Veselov equation

## Abstract.

Darboux transformation is developed to systematically find variable separation solutions for the Nizhnik-Novikov-Veselov equation. Starting from a seed solution with some arbitrary functions, the once Darboux transformation yields the variable separable solutions which can be obtained from the truncated Painlevé analysis and the twice Darboux transformation leads to some new variable separable solutions which are the generalization of the known results obtained by using a guess ansatz to solve the generalized trilinear equation.

*PACS.*02.30.Jr, 02.30.Ik, 05.45.Yv.

To find some exact solutions of integrable systems has focused many mathematicians and physicists’ attention since the soliton theory came into being at the 1960’s. There are many important methods to obtain the special solutions of a given soliton equation. Some of the most important methods are the inverse scattering transformation (IST) approach, the bilinear form, symmetry reduction, Bäcklund transformation and Darboux transformation  etc. In comparison with the linear case, it is known that IST is an extension of the Fourier transformation in the nonlinear case. In addition to the Fourier transformation, there is another powerful tool called the variable separation method  in the linear case. Recently, some kinds of variable separation approaches had been developed to find new exact solutions of nonlinear models, say, the classical method, the differential Stäckel matrix approach , the geometrical method , the ansatz-based method , functional variable separation approach, the derivative dependent functional variable separation approach, the formal variable separation approach (nonlinearization of the Lax pairs or symmetry constraints)  and the informal variable separation methods .

Especially, for various (2+1)-dimensional nonlinear physics models, a quite universal formula

$$U\frac{2(a_1a_2a_0a_3)q_yp_x}{(a_0+a_1p+a_2q+a_3pq)^2},$$
(1)
where $`a_i,(i=0,1,2,3)`$ are arbitrary constants and $`p=p(x,t)`$ and $`q=q(y,t)`$ are arbitrary constants of the indicated variables, is found by using the informal variable separation approach. Starting from the universal formula (1), abundant localized excitations like the dromions, lumps, ring solitons, breathers, instantons, solitoffs, fractal and chaotic patterns are found. Now a very important question is can we find the universal formula from other well known powerful methods like the IST approach, dressing method, Darboux transformation (DT) etc?

DT is one of the most powerful methods to construct a broad class of considerable physical interest and important nonlinear evolution equations such as the well-known Korteweg-de Vries (KdV) equation, the Kadomtsev-Petviashvili (KP) equation, the Davey-Stewartson (DS) equation, the sine-Gordon (SG) equation  and so on. In this letter, we use the Darboux transformation to study the variable separable solutions for the (2+1) dimensional Nizhnik-Novikov-Veselov system.

$$u_t=u_{xxx}+u_{yyy}+3(vu)_x+3(uw)_y,$$
(2)
$$u_x=v_y,u_y=w_x.$$
(3)
The (2+1)-dimensional NNV equation is an only known isotropic Lax integrable extension of the well-known (1+1)-dimensional KdV equation. Many authors have studied the solutions of the NNV equation. For example, Boiti et al solved the NNV equation via the IST; Tagami and Hu and Li obtained the soliton-like solutions of the NNV equation by means of Bäcklund transformation; Hu also gave out the nonlinear superposition formula of the NNV equation; Some special types of multi-dromion solutions were found by Radha and Lakshmanan; The generalized localized excitations expressed by (1) were given in  and  and the special binary Darboux transformation was given in .

It is known that the NNV equation system (2) and (3) can be represented as a compatibility condition of the linear system

$$\mathrm{\Phi }_{xy}+u\mathrm{\Phi }=0,$$
(4)
$$\mathrm{\Phi }_t=\mathrm{\Phi }_{xxx}+\mathrm{\Phi }_{yyy}+3v\mathrm{\Phi }_x+3w\mathrm{\Phi }_y.$$
(5)
In general, in order to construct the solutions of a given equation by means of Darboux transformation, one may use a fixed solution for a special spectrum parameter of the Lax pair. Then one can get a new solution from an old one with Darboux transformation. Without spectral parameter in the Lax pair of the NNV equation, we have to construct a binary Darboux transformation for the NNV equation. At the same time, with the variable separated approach, two arbitrary functions can be entered into the solution.

It is straightforward to see that the NNV system (2) and (3) possesses the following trivial solution

$$u=0,v=v_0(x,t),w=w_0(y,t),$$
(6)
where $`v_0(x,t)`$ and $`w_0(y,t)`$ are arbitrary functions of $`\{x,t\}`$ and $`\{y,t\}`$ respectively.

In order to find some new solutions via Darboux transformation and the seed solution (6), the key step is to find the fixed solution of the Lax pair (4) and (5) with the seed (6).

It is evident that Eq.(4) have a variable separation solution in the form

$$\mathrm{\Phi }_0=p+q$$
(7)
where $`p=p(x,t)`$, $`q=q(y,t)`$ are two arbitrary functions when $`u=0`$. Substituting Eq. (7) and Eq. (6) into Eq. (5) yields

$$p_t+q_t=p_{xxx}+q_{yyy}+3v_0p_x+3w_0q_y.$$
(8)
It is clear that Eq. (8) can be solve by the usual variable separable approach and the result reads

$$p_t=p_{xxx}+3v_0p_x+c(t),$$
(9)
$$q_t=q_{yyy}+3w_0q_yc(t)$$
(10)
where $`c(t)`$ is an arbitrary function of $`t`$. For given functions $`v_0,w_0`$ and $`c(t)`$, it is still difficult to solve (9) and (10). However, because of the arbitrariness of the functions $`v_0`$ and $`w_0`$, we can solve the problem in an alternative way. If we consider the functions $`p`$ and $`q`$ as arbitrary functions, then $`v_0`$ and $`w_0`$ can be solved out from (9) and (10)

$$v_0=(3p_x)^1(p_tp_{xxx}c(t))$$
(11)
$$w_0=(3q_y)^1(q_tq_{yyy}+c(t)).$$
(12)
As usual, the binary Darboux transformation for the NNV equation can be constructed by introducing the closed 1-form

$`\omega (\mathrm{\Phi },\mathrm{\Phi }_0)`$ $`=`$ $`(\mathrm{\Phi }\mathrm{\Phi }_{0x}\mathrm{\Phi }_x\mathrm{\Phi }_0)dx(\mathrm{\Phi }\mathrm{\Phi }_{0y}\mathrm{\Phi }_0\mathrm{\Phi }_y)dy+\{\mathrm{\Phi }(\mathrm{\Phi }_{0xxx}\mathrm{\Phi }_{0yyy})`$ (13)
$`+\mathrm{\Phi }_0(\mathrm{\Phi }_{yyy}\mathrm{\Phi }_{xxx})+2(\mathrm{\Phi }_{xx}\mathrm{\Phi }_{0x}\mathrm{\Phi }_x\mathrm{\Phi }_{0xx}+\mathrm{\Phi }_{0yy}\mathrm{\Phi }_y\mathrm{\Phi }_{0y}\mathrm{\Phi }_{yy})`$
$`+3v(\mathrm{\Phi }\mathrm{\Phi }_{0x}\mathrm{\Phi }_x\mathrm{\Phi }_0)+3w(\mathrm{\Phi }_0\mathrm{\Phi }_y\mathrm{\Phi }_{0y}\mathrm{\Phi })\}dt.`$
Then the new wave function is

$`\mathrm{\Phi }[1]=\mathrm{\Phi }_0^1{\displaystyle \omega }`$ (14)
and the equation system Eq.(4) and (5) is covariant with respect to the transformation (14) and the transformed coefficients $`u[1]`$, $`v[1]`$ and $`w[1]`$ are

$`u[1]`$ $`=`$ $`2(\mathrm{ln}\mathrm{\Phi }_0)_{xy}={\displaystyle \frac{2p_xq_y}{(p+q)^2}},`$ (15)
$`v[1]`$ $`=`$ $`v_0+2(\mathrm{ln}\mathrm{\Phi }_0)_{xx}=v_0+{\displaystyle \frac{2p_{xx}}{p+q}}{\displaystyle \frac{p_x^2}{(p+q)^2}},`$ (16)
$`w[1]`$ $`=`$ $`w_0+2(\mathrm{ln}\mathrm{\Phi }_0)_{yy}=w_0+{\displaystyle \frac{2q_{yy}}{p+q}}{\displaystyle \frac{q_y^2}{(p+q)^2}},`$ (17)
where $`v_0`$ and $`w_0`$ determined by Eqs.(11) and (12). From (15), we see that the solution $`u[1]`$ obtained by the once Darboux transformation is only a special case of the universal expression (1) with $`a_3=0`$. The solution (15)–(17) can also be obtained by the truncated Painlevé expansion using the similar method as for the AKNS system.

By means of the iteration of the Darboux transformation, we can construct the second Darboux transformation from the first Darboux transformation. We also take two special solutions of the Lax pair (4) and (5) as the variable separable ones, $`\mathrm{\Phi }_1=p_1(x,t)+q_1(y,t)`$ and $`\mathrm{\Phi }_2=p_2(x,t)+q_2(y,t)`$, then the transformed wave function is

$$\mathrm{\Phi }[2]=\omega (\mathrm{\Phi }_1,\mathrm{\Phi }_2)$$
(18)
which leads to the formula

$`u[2]`$ $`=`$ $`2(\mathrm{ln}{\displaystyle \omega (\mathrm{\Phi }_1,\mathrm{\Phi }_2)})_{xy},`$ (19)
$`v[2]`$ $`=`$ $`v_0+2(\mathrm{ln}{\displaystyle \omega (\mathrm{\Phi }_1,\mathrm{\Phi }_2)})_{xx},`$ (20)
$`w[2]`$ $`=`$ $`w_0+2(\mathrm{ln}{\displaystyle \omega (\mathrm{\Phi }_1,\mathrm{\Phi }_2)})_{yy}`$ (21)
According to the selection for the arbitrary functions, the potential functions $`v_0`$ and $`w_0`$ are related to the functions $`p_1,p_2,q_1`$ and $`q_2`$ by

$`p_{1t}p_{1xxx}c_1(t)v_0p_{1x}=0,`$ (22)
$`q_{1t}q_{1yyy}+c_1(t)w_0q_{1y}=0,`$ (23)
$`p_{2t}p_{2xxx}c_2(t)v_0p_{2x}=0,`$ (24)
$`q_{2t}q_{2yyy}+c_2(t)w_0q_{2y}=0.`$ (25)
Now we can consider one of the functions $`p_1=p_1(x,t)`$ and $`p_2=p_2(x,t)`$ as an arbitrary function of $`\{x,t\}`$ and one of the functions $`q_1=q_1(y,t)`$ and $`q_2=q_2(y,t)`$ as an arbitrary function of $`\{y,t\}`$ while the remained functions of $`v_0,w_0,p_1,p_2,q_1`$ and $`q_2`$ should be determined by (22), (23), (24) and (25).

It is not easy to solve the equation system (22), (23), (24) and (25). However, if we choose

$`q_1=0,p_2=0,,`$ (26)
then the Eqs. (19)-(21) can be simplified into

$`u[2]`$ $`=`$ $`2\mathrm{ln}(c+p_1q_2)_{xy},`$ (27)
$`v[2]`$ $`=`$ $`v_0+2\mathrm{ln}(c+p_1q_2)_{xx},`$ (28)
$`w[2]`$ $`=`$ $`w_0+2\mathrm{ln}(c+q_1q_2)_{yy},`$ (29)
where $`c`$ is an arbitrary integral constant, $`p_1`$ and $`q_2`$ are arbitrary functions of $`\{x,t\}`$ and $`\{y,t\}`$ while $`v_0`$ and $`w_0`$ are fixed by (22) and (25). It is easy to prove that the results (27)–(29) are equivalent to those obtained via the usual variable separation of the multi-linear equation. Actually, by re-writing $`c,p_1`$ and $`q_2`$ as

$`c=a_0{\displaystyle \frac{a_1a_2}{a_3}},p_1=a_2+a_3p,q_2=q+{\displaystyle \frac{a_1}{a_3}},`$ (30)
the negative value of the right hand side of (27) is transformed to the universal quantity expressed by (1).

Various types of coherent localized structures for the physical field $`u`$ of the NNV system had been described in  and  thanks to the arbitrariness of the functions $`p`$ and $`q`$. We do not repeated these known localized excitations but write down only one new type of localized solutions for the field $`u`$ expressed by (27) with (30) to complement the results of .

In addition to the continuous localized excitations in (1+1)-dimensional nonlinear systems, some types of significant weak solutions like the peakons and compactons  have been attract much attention of both mathematicians and physicists. In  and , the possible (2+1)-dimensional localized peakons had been given. The (1+1)-dimensional compacton solutions which describes the typical (1+1)-dimensional soliton solutions with finite wavelength when the nonlinear dispersion effects were firstly given by Rosenau et al. and may have many interesting properties and possible physical applications. For instance, the compacton equations may be used to study the motion of ion-acoustic waves and a flow of a two layer liquid. In , the Painlevé integrability of two sets of Korteweg-de Vries (KdV) type and modified KdV type compacton equations are proved. Because of the entrance of arbitrary functions in the (2+1)-dimensional nonlinear physics models, it is easy to find some types of multiple compacton solutions by selecting the arbitrary functions appropriately. For instance, if we fixed the functions $`p`$ and $`q`$ as

$`p={\displaystyle \underset{i=1}{\overset{N}{}}}\{\begin{array}{cc}0\hfill & x+c_itx_{0i}\frac{\pi }{2k_i}\hfill \\ b_i\mathrm{sin}(k_i(x+c_itx_{0i}))+b_i\hfill & x_{0i}\frac{\pi }{2k_i}<x+c_itx_{0i}+\frac{\pi }{2k_i}\hfill \\ 2b_i\hfill & x+c_it>x_{0i}+\frac{\pi }{2k_i}\hfill \end{array},`$ (34)
and

$`q={\displaystyle \underset{j=1}{\overset{M}{}}}\{\begin{array}{cc}0\hfill & yy_{0j}\frac{\pi }{2l_j}\hfill \\ d_j\mathrm{sin}(l_i(yy_{0j}))+d_j\hfill & y_{0j}\frac{\pi }{2l_j}<yy_{0j}+\frac{\pi }{2l_j}\hfill \\ 2d_j\hfill & y>y_{0j}+\frac{\pi }{2l_j}\hfill \end{array}.`$ (38)
where $`b_i,k_i,x_{0i},d_j,l_j`$ and $`y_{0j}`$ are all arbitrary constants, then the solution (27) with (30) becomes a multi-compacton solution.

From (34) and (38), one can see that the piecewise functions $`p`$ and $`q`$ of the compacton solutions are once differentiable

$`p_x={\displaystyle \underset{i=1}{\overset{N}{}}}\{\begin{array}{cc}0\hfill & x+c_itx_{0i}\frac{\pi }{2k_i}\hfill \\ b_ik_i\mathrm{cos}(k_i(x+c_itx_{0i}))\hfill & x_{0i}\frac{\pi }{2k_i}<x+c_itx_{0i}+\frac{\pi }{2k_i}\hfill \\ 0\hfill & x+c_it>x_{0i}+\frac{\pi }{2k_i}\hfill \end{array},`$ (42)
and

$`q_y={\displaystyle \underset{j=1}{\overset{M}{}}}\{\begin{array}{cc}0\hfill & yy_{0j}\frac{\pi }{2l_j}\hfill \\ d_jl_j\mathrm{cos}(l_j(yy_{0j}))\hfill & y_{0j}\frac{\pi }{2l_j}<yy_{0j}+\frac{\pi }{2l_j}\hfill \\ 0\hfill & y>y_{0j}+\frac{\pi }{2l_j}\hfill \end{array},`$ (46)
So the multi-compacton solutions expressed by (1) (i.e., (27) with (30)) are still continuous any where though there are some isolated non-differentiable lines.

In (1+1)-dimensions, a non-differentiable solution, $`u=u_0`$, like the the compacton (and/or peakon) is called a weak solution of a nonlinear (1+1)-dimensional PDE (partial differential equation)

$`F(u,u_x,u_t,u_{xx},\mathrm{})F(u)=0`$ (47)
under the meaning that though the compacton (and/or peakon) solution is non- differentiable at some points $`x=x_i(t)`$, the distribution, $`F(u_0)=f(\delta (xx_i(t))`$ (where $`\delta (xx_i(t))`$ is a Dirac $`\delta `$ function and $`f`$ is a function of the Dirac $`\delta `$ function and its derivatives) is really a zero distribution that means $`_{\mathrm{}}^+\mathrm{}f\psi (x,t)=0`$ for arbitrary $`\psi `$.

However, in (2+1)-dimensions, the (2+1)-dimensional compactons (and the peakons reported in ) are exact solutions of some (2+1)-dimensional nonlinear equations are guaranteed by the variable separation procedure. When we substituting the piecewise solutions into the nonlinear PDEs, the Dirac delta function(s) in $`x`$ ($`y`$) directions will be vanished by the differential operator in other direction, $`_y`$ ($`_x`$). In other words, we can take (have taken) $`_y\psi (x,t)=0`$ no matter the function $`\psi `$ is a continuous function of $`\{x,t\}`$ or a generalized distribution functions (with some Dirac delta functions) of $`\{x,t\}`$.

Fig.1 is the evolution plot of a three compacton solution (27) with (30), (34), (38) and

$`N=3,M=1,a_0=20,a_1=a_2=25a_3=1,b_1=b_3=c_2=c_3=2,`$
$`b_2=c_1=d_1=k_1=k_2=k_3=l_1=1,x_{01}=x_{02}=x_{03}=y_{01}=0.`$ (48)
In , we have pointed out that (i) the interaction between two travelling ring shape soliton solutions is completely elastic and (ii) the interaction between two travelling peakons is not completely elastic, two peakons may completely exchange their shapes. From Fig. 1, we see that the interaction between two compactons exhibits a new phenomenon. The interaction is nonelastic but two compactons do not completely exchange their shapes.

Because of the arbitrariness of the functions $`p`$ and $`q`$ may have also quite rich structures. For instance, (2+1)-dimensional compactons may be only compacted at one direction. For conveniece later we call this type of compactons the partial compactons. Fig. 2 is a plot of a special partial compacton solution (27) with (30), (34), (42),

$`N=1,a_0=20,a_1=a_2=25a_3=c_1=k_1=1,x_{01}=0.`$ (49)
and

$`q=10\mathrm{tanh}(yt),`$ (50)
at $`t=0`$. For the partial compacton shown by Fig. 2, the quantity $`G`$ possesses the following quite simple form ($`u2u[2]`$)

$`u2=\{\begin{array}{cc}0\hfill & |xt|\frac{1}{2}\pi ,\hfill \\ {\displaystyle \frac{100\mathrm{cos}(xt)\mathrm{sech}^2(yt)}{[(9010\mathrm{sin}(xt))+(464\mathrm{sin}(xt))\mathrm{tanh}(yt)]^2}}\hfill & |xt|<\frac{1}{2}\pi .\hfill \end{array}`$ (53)
Similar to the first kind of compacton (full compacted) solution as shown in Fig. 1, the detailed study of the interaction between partial compactons is also non-completely elastic and they will partially exchange their shapes.

If the arbitrary functions $`q_1`$ and $`p_2`$ are not taken as in (26) but are solved out from (23) and (24), then we can obtain some further new types of exact solutions which can not be obtained by the usual variable separation approach. Actually, when the arbitrary functions $`p_1`$ and $`q_2`$ are fixed by (22) and (25) at the same time, then the corresponding solution(s) for the functions $`q_1`$ and $`p_2`$ can be found by solving the linear equations (24) and (23).

In principle, the more kinds of exact solutions can be obtained from further Darboux transformations starting from the seed solution {(18), (19), (20), (21)}. However, we do not discuss these types of solutions further because of their complexity.

In summary, the variable separation solutions can be obtained not only by the truncated Painlevé expansion and the generalized multi-linear equations, but also by other well known approaches especially by the Darboux transformation. In this short paper, the Darboux transformation is successfully used to find the variable separable solutions of the (2+1)-dimensional NNV equation. The recursive Darboux transformation may yield further new types of variable separable solutions while the new variable separable function should satisfy a further constrained condition, say, (24) with $`v_0`$ being given by (22).

By selecting the arbitrary functions appropriately, one may obtain abundant localized excitations like the dromions, lumps, ring solitons, breathers, instantons, solitoffs, peakons, fractal and chaotic patterns. In addition to these types of localized excitations, a further type of the localized excitations, compactons, is given in this paper. The (2+1)-dimensional compactons discussed here possess the different type of interaction properties as that of the ring solitons and peakons. The interactions among compactons are not completely elastic and do not exchange their shapes completely.

The authors are in debt to thank the helpful discussions with the professors X. B. Hu and Y. S. Li and Doctor. C. L. Chen. The work was supported by the National Outstanding Youth Foundation of China (No.19925522), the Research Fund for the Doctoral Program of Higher Education of China (Grant. No. 2000024832) and the Natural Science Foundation of Zhejiang Province of China.


text:
# Tomography of nuclear matter: Comparing Drell-Yan with deep inelastic scattering data

## Drell-Yan process

In the absence of nuclear effects, the leading-order (LO) DY production cross section in hadron-nucleus reactions is given by

$`{\displaystyle \frac{\mathrm{d}\sigma (hA)}{\mathrm{d}x_1}}={\displaystyle \frac{8\pi \alpha ^2}{9x_1s}}{\displaystyle \underset{q}{}}e_q^2{\displaystyle }{\displaystyle \frac{\mathrm{d}M}{M}}[`$ $`Z\left(f_q^h(x_1)f_{\overline{q}}^p(x_2)+f_{\overline{q}}^h(x_1)f_q^p(x_2)\right)`$ (1)
$`+(AZ)(f_q^h(x_1)f_{\overline{q}}^n(x_2)+f_{\overline{q}}^h(x_1)f_q^n(x_2))]`$
where $`x_1`$ (resp. $`x_2`$) is the momentum fraction carried by the beam (resp. target) parton, $`\sqrt{s}`$ the center-of-mass energy of the hadronic collision, and $`M=x_1x_2s`$ the invariant mass. As long as the parton densities do not show any isospin dependence ($`f_i^p(x_2)=f_i^n(x_2)`$), the DY cross section (1) is proportional to the atomic mass number $`A`$. Hence the measured ratio in a heavy over a light nucleus

$$R^h(A/B,x_1)=\frac{B}{A}\left(\frac{\mathrm{d}\sigma (hA)}{\mathrm{d}x_1}\right)\times \left(\frac{\mathrm{d}\sigma (hB)}{\mathrm{d}x_1}\right)^1$$
(2)
is equal to one. However, the multiple scatterings encountered by the incoming (anti-) quark reduce its momentum fraction from $`x_1+ϵ/E_h`$ to $`x_1`$ at the point of fusion, $`ϵ`$ being the energy loss and $`E_h`$ the beam energy in the nucleus rest frame. As a consequence, the projectile parton distribution function should be evaluated at $`(x_1+\mathrm{\Delta }x_1)`$ in the nuclear cross section (1). Because the parton densities dramatically drop at large $`x_1`$, even a tiny shift $`\mathrm{\Delta }x_1`$ may substantially suppress the DY yield in large nuclei. In the general case, the DY cross section reads

$`{\displaystyle \frac{\mathrm{d}\sigma (hA)}{\mathrm{d}x_1}}`$ $`=`$ $`{\displaystyle \frac{8\pi \alpha ^2}{9x_1s}}{\displaystyle \underset{q}{}}e_q^2{\displaystyle }{\displaystyle \frac{\mathrm{d}M}{M}}{\displaystyle }\mathrm{d}ϵD^{in}(ϵ,x_1E_h)\times `$ (3)
$`[Zf_q^h(x_1+\mathrm{\Delta }x_1)f_{\overline{q}}^{p/A}(x_2)+(AZ)f_q^h(x_1+\mathrm{\Delta }x_1)f_{\overline{q}}^{n/A}(x_2)`$
$`+Zf_{\overline{q}}^h(x_1+\mathrm{\Delta }x_1)f_q^{p/A}(x_2)+(AZ)f_{\overline{q}}^h(x_1+\mathrm{\Delta }x_1)f_q^{n/A}(x_2)],`$
where $`D^{in}(ϵ,E_q)`$ is the probability that an incoming quark (with energy $`E_q`$) suffers an energy loss $`ϵ`$ . Furthermore, the effects of shadowing are taken into account through the nuclear parton densities of the target parton, $`f_i^{p/A}(x_2)f_i^p(x_2)`$. The cross section (3) has been fitted to NA3 data in $`E_\pi ^{}=150`$ and 280 GeV pion induced reactions on hydrogen and platinum targets, with the mean energy loss per unit length (hence, the transport coefficient) being the only free parameter in the calculation. We found

$$\frac{dE}{dz}[\mathrm{GeV}/\mathrm{fm}]=\frac{\omega _c}{9L}=\left(0.4\pm 0.3\right)\times \left(\frac{L}{10\mathrm{fm}}\right)$$
(4)
which corresponds to a transport coefficient111Note that in  we used the BDMPS mean energy loss $`ϵ=\omega _c/3`$ (valid for outgoing quarks), while we take here the correct Eq. (48) of  for incoming quarks, hence the discrepancy between the transport coefficient. The energy loss per unit length, fitted to the data, remaining of course unchanged. $`\widehat{q}=0.72\pm 0.54`$ GeV/fm2. Moreover, it is worth emphasizing that this fitted value does not depend much on shadowing corrections . Taking this estimate, we now investigate the effect of quark energy loss on the hadron spectra measured in deep inelastic scattering.

## Hadron production in deep inelastic scattering

The HERMES collaboration at DESY recently reported on hadron yields measured in electron-nucleus collisions . They measured the production ratio

$$R_A^h(z,\nu )=\frac{1}{N_A^e(\nu )}\frac{N_A^h(z,\nu )}{d\nu dz}/\frac{1}{N_D^e(\nu )}\frac{N_D^h(z,\nu )}{d\nu dz}$$
(5)
in a “heavy” (N and Kr) over a light (D) nucleus for a given hadron species $`h`$. Here, $`\nu `$ denotes the virtual photon energy in the lab frame, $`z`$ the momentum fraction carried by the produced hadron, and where the multiplicity of produced electrons $`N_A^e`$ normalizes the hadron yield $`N_A^h`$. The hadron multiplicity in (5) can be computed in perturbation theory. Assuming for simplicity only the valence up quark to contribute when $`x`$ is not too small, Eq. (5) will approximately be given by the ratio of the $`uh`$ fragmentation functions

$$R_A^h(z,\nu )D_u^h(z,Q^2,A)/D_u^h(z,Q^2,D).$$
(6)
Therefore, the nuclear dependence of the fragmentation functions might be revealed through the measure of $`R^h`$. Because of the shift in the quark energy from $`E_q\nu `$ to $`E_q\nu ϵ`$ at the time of the hadronization, the nuclear fragmentation functions may be modeled according to

$$zD_f^h(z,Q^2,A)=_0^{\nu E_h}𝑑ϵD^{\mathrm{o}ut}(ϵ,\nu ϵ)z^{}D_f^h(z^{},Q^2).$$
(7)
with $`z^{}=z/(1ϵ/\nu )`$. Taking the transport coefficient previously extracted from DY data, the $`\nu `$ dependence of $`R^h(\nu )`$ in a krypton over a deuterium target is computed222We use the GRV98 LO parton densities  with the Kretzer parameterization for the fragmentation functions . (Figure 2). The trend is well reproduced for all hadron species, although the calculation for the pions ($`\pi ^++\pi ^{})`$ somehow underpredicts the effect.

We showed that both the DY process as well as hadron production in DIS turn out to be very sensitive to the energy loss of fast quarks traveling through nuclear matter. Using the transport coefficient previously extracted from NA3 DY measurements, the quenching of hadron spectra proved in good agreement with HERMES preliminary data.


text:
# Spin current in the Kondo lattice model

(January 4, 2004)

## Abstract

By using the projection operator technique it is observed that the strong Hund’s rule coupling and s-d interaction in transition metal elements may lead to an effective coupling between the spin current and spin chirality. As a result, the spin chirality can be regarded as a driving force to produce a spin current. The spin current may give rise to a novel type of field acting on the spins. A spin battery is designed based on the interactions between the spin current and spin chirality.

Spin-dependent effects, such as giant or colossal magnetoresistance and magnetization switching, arise from the interaction between spin of charge carrier and an applied external magnetic field. In spintronics, instead of charge, electron spin carries information and it is possible that capability and performance can be enhanced in spintronics devices.Wolf01 ; Ziese01 ; Awschalom02  Possible applications and fundamental science involved make the study of spin-dependent transport an intensive field in condensed matter physics. Recently, several ways to produce a pure spin current are proposed, such as anomalous Hall effect in ferromagnetic metalsHirsch99 ; Zhang00 , ferromagnetic resonance Berger99 ; Brataas02 , and the spin current in a spin spiral state.Shen97 ; Konig01  In this paper, we study spin current in the Kondo lattice model, a simplified model containing essential physics for many ferromagnetic metals or semiconductors. We find a new type of interaction between spin current and spin chirality (defined later) in the strong Hund coupling limit. This indicates that the spin chirality can be regarded as a driving force to produce spin current, similar as an electric field does to the electric current. Furthermore, the spin current induces a new type of an effective magnetic field. Based on the interaction between spin current and spin chirality, we propose a scheme to design a spin battery. As a demonstration, we study the spin transport in the well known spiral state to bring out the above mentioned spin-current effects.

The Hund’s rule coupling and s-d exchange interaction are crucial in transition metal elements. It plays a key role in the formation of metallic ferromagnetism in colossal magnetoresistance materials Zener51 ; Anderson55 ; Kubo72  and diluted magnetic semiconductors.Ohno98 ; Dietl00  In this paper we focus on doped transition metal oxides and/or ferromagnetic semiconductors. The minimal model for these materials is the Kondo lattice model,

$$H=t\underset{n,m,\sigma }{}c_{n,\sigma }^{}c_{m,\sigma }\frac{J}{2}\underset{n,\sigma \sigma ^{}}{}𝐒_n\tau _{\sigma \sigma ^{}}c_{n,\sigma }^{}c_{n,\sigma ^{}},$$
(1)
where $`c_{n,\sigma }^{}`$ and $`c_{n,\sigma }`$ are the creation and annihilation operators for conduction electrons, $`𝐒_n`$ is the spin operators for the localized spin, and $`\tau `$ are the Pauli matrices and $`\tau _0`$ is the $`2\times 2`$ identity matrix. The summation over $`n,m`$ runs for the nearest neighbor pairs of lattice sites. In this paper we will set $`\mathrm{}=c=1`$ and the lattice space is $`a`$ and is set to unit.Scalapino93  For the Hund’s rule coupling, $`J`$ is always positive. The model has been studied extensively.Tsunetsugu97  The spin current operator along $`\alpha `$ axis in the tight binding approximation is defined as $`t\mathrm{}𝐈_{nn+\alpha }^s`$ and

$$𝐈_{nn+\alpha }^s=\frac{i}{2}\underset{\sigma ,\sigma ^{}}{}\left(c_{n,\sigma }^{}\tau _{\sigma \sigma ^{}}c_{n+\alpha ,\sigma ^{}}c_{n+\alpha ,\sigma }^{}\tau _{\sigma \sigma ^{}}c_{n,\sigma ^{}}\right),$$
(2)
where $`\alpha (=x,y,z)`$ are vectors pointing to the nearest neighbor sites, which have one longitudinal and two transverse components. The spin polarized currents and electric current become equal in a fully polarized ferromagnetic system.

Before we discuss the spin transport of Eq.(1) in an arbitrary dimension, let us first to illustrate the physics by considering the well known spin spiral state in one-dimensional system. In general, the spin spiral state may not be the ground state of Eq.(1), however, it can be stabilized by external environments, such as by applying external fields along different directions at the ends of the system.Spiral  Our purpose here is to demonstrate the consequence of the spin spiral state. We parameterize $`𝐒_n`$ by the polar angles $`\theta _n`$ and $`\phi _n,`$ i.e., making the classical spin approximation. The spin spiral state is defined by $`𝐒_n/S=(\mathrm{sin}\theta \mathrm{cos}n\phi ,\mathrm{sin}\theta \mathrm{sin}n\phi ,\mathrm{cos}\theta ).`$ The spin chirality defined as the cross times between the two neighboring spins: $`𝐒_n\times 𝐒_{n+1}_z=S^2\mathrm{sin}^2\theta \mathrm{sin}\phi .`$ In the k- space, we introduce a spinor operator $`\mathrm{\Phi }^{}(k)=(c_{ka+\phi /2,}^{},c_{ka\phi /2,}^{}).`$ The mean field Hamiltonian is written as

$$H=\underset{k}{}\mathrm{\Phi }^{}(k)H(k)\mathrm{\Phi }(k),$$
(3)
where $`H(k)=g_0\tau _0g_x\tau _xg_z\tau _z`$ with $`g_0=\mu +2t\mathrm{cos}\frac{\phi }{2}\mathrm{cos}ka;`$ $`g_x=\frac{J}{2}\mathrm{sin}\theta ;`$ $`g_z=\frac{J}{2}\mathrm{cos}\theta 2t\mathrm{sin}\frac{\phi }{2}\mathrm{sin}ka.`$ $`\mu `$ is the chemical potential. The single particle Green function in a $`2\times 2`$ matrix form is given by

$`G(k,i\omega _n)`$ $`=(i\omega _n\tau _0H(k))^1`$
$`={\displaystyle \frac{\left[(i\omega _n+g_0)\tau _0g_x\tau _xg_z\tau _z\right]}{\left(i\omega _n\omega _++\mu \right)\left(i\omega _n\omega _{}+\mu \right)}},`$ (4)
where $`\omega _n=(2n+1)\pi kT`$ with an integer $`n`$ and $`T`$ is the temperature. The two branches of spectra are

$$\begin{array}{c}\omega _\pm (k)=2t\mathrm{cos}\frac{\phi }{2}\mathrm{cos}ka\hfill \\ \pm \left[\left(\frac{J}{2}\mathrm{cos}\theta 2t\mathrm{sin}\frac{\phi }{2}\mathrm{sin}ka\right)^2+\frac{J^2}{4}\mathrm{sin}^2\theta \right]^{1/2}.\hfill \end{array}$$
(5)
In the k-space, the expectation value of the longitudinal spin current operator is written in terms of the Green’s function,

$$I_z^s=\frac{1}{N}\underset{k}{}\left[A_+n_F(\omega _+\mu )A_{}n_F(\omega _{}\mu )\right],$$
(6)
where

$`A_\pm (k)`$ $`=\pm \mathrm{sin}{\displaystyle \frac{\phi }{2}}\mathrm{cos}ka`$
$`+{\displaystyle \frac{t\mathrm{sin}\phi \mathrm{sin}^2ka\frac{J}{2}\mathrm{cos}\theta \mathrm{cos}\frac{\phi }{2}\mathrm{sin}ka}{\sqrt{\left(\frac{J}{2}\mathrm{cos}\theta 2t\mathrm{sin}\frac{\phi }{2}\mathrm{sin}ka\right)^2+\frac{J^2}{4}\mathrm{sin}^2\theta }}}.`$
We also calculated the electric current $`I_z^c`$ and found that $`I_z^c=0`$. Thus in the spiral state there exists a spin current without a charge current. In Fig. 1 we plot spin current as a function of electron density in the state with $`\theta =\pi /2,`$ and $`\phi =0.2\pi `$ for different $`J`$. We see that a larger $`J`$ produces a higher spin current. For a finite $`J`$ the current may change direction with change of density. In Fig. 2 we plot the relation between spin current and spin chirality at the density of $`0.2`$. The spin current increases with the spin chirality.

In the large $`J`$ limit, the problem can be simplified. An equivalent Hamiltonian of Eq.(1) or so-called double exchange model in the large $`J`$ limit is given byHartman96

$$H_{de}=\underset{n,m}{}t_{nm}\alpha _n^+\alpha _m,$$
(7)
with $`t_{nm}=`$ t$`\left[\mathrm{cos}\frac{\theta _n}{2}\mathrm{cos}\frac{\theta _m}{2}+\mathrm{sin}\frac{\theta _n}{2}\mathrm{sin}\frac{\theta _m}{2}e^{i(\phi _n\phi _m)}\right]`$, and $`\alpha _n=\mathrm{cos}(\theta _n/2)(1n_n)c_n+\mathrm{exp}[i\phi _n]\mathrm{sin}(\theta _n/2)(1n_n)c_n`$ and satisfies the anticommutation relation, $`[\alpha _k,\alpha _k^{}^{}]=\delta _{k,k^{}}.`$ The energy excitations and their eigenstates can be obtained from the diagonalized Hamiltonian,

$$H_{de}=2t\left[\mathrm{cos}^2\frac{\theta }{2}\mathrm{cos}ka+\mathrm{sin}^2\frac{\theta }{2}\mathrm{cos}(ka+\phi )\right]\alpha _k^{}\alpha _k,$$
(8)
where $`\alpha _k`$ is the Fourier transform of $`\alpha _n.`$ In the case when the spins of conduction electrons are frozen by the localized spins, and the quasiparticles become spinless fermions. The two Fermi momenta are $`k_\pm ^Fa=k_0\pm \rho \pi `$ where $`k_0=\mathrm{arctan}[\mathrm{sin}^2\theta \mathrm{sin}\phi /(\mathrm{cos}^2\theta +\mathrm{sin}^2\theta \mathrm{cos}\phi )]`$ and $`\rho `$ is the ratio of the number of conduction electrons to the number of total lattice sites and $`\rho /a`$ is the density of conduction electrons. $`k_0`$ is nonzero when the order parameter for the spin spiralling is nonzero. Its lowest energy state is $`|GS=\underset{k_F^{}kk_F^+}{}\alpha _k^{}|0.`$ The spin current in the state can be calculated from $`𝐈_z^s=_nGS\left|𝐈_{nn+1}^s\right|GS/N.`$ We have an analytical expression for spin current,

$$I_z^s=t\mathrm{}\frac{\mathrm{sin}\rho \pi }{\pi }\frac{\mathrm{sin}^2\theta \mathrm{sin}\phi }{\sqrt{1\mathrm{sin}^2\theta \mathrm{sin}^2(\phi /2)}}.$$
(9)
The result is equal to what we obtain by taking the large limit in Eq.(6). In the point of view of the Berry phase the spin spiral state acquires a constant Berry phase in the spin-dependent renormalized factor in the Hamiltonian (Eq. 7), i.e., $`t_{nn+1}=\left|t_{nn+1}\right|e^{i\delta \phi }.`$ The non-zero phase $`\delta \phi `$ drives a spin current.

By analogy with electric field and resistance, we introduce the concept of spin current resistance or its reciprocal, spin conductance. The role of $`𝐒_n\times 𝐒_m`$ here is analogous to the local electric field. In the spin spiral state the spin chirality is non-zero, $`𝐒_n\times 𝐒_{n+1}_z=S^2\mathrm{sin}^2\theta \mathrm{sin}\phi .`$ In Eq.(9) the current is proportional to the spin chirality $`𝐒_n\times 𝐒_{n+1}_z,`$ and can be written in a compact form $`I_z^s=c𝐒_n\times 𝐒_{n+1}_z/S^2`$ by introducing a spin conductance

$$c=\mathrm{sin}\rho \pi /\left[\pi \sqrt{1\mathrm{sin}^2\theta \mathrm{sin}^2(\phi /2)}\right]$$
(10)
in a large J limit. The spin conductance is determined by the filling of electrons and the magnetic structure. Just like an electric conductance, the spin conductance reflects intrinsic properties of spin transport in materials. It depends on the spin-dependent scattering mechanism of electrons. In Eq. (10), for the case of $`\rho 0`$ or $`1,`$ c is proportional to the density of electrons $`\rho `$ or the density of holes $`1\rho `$. For a small spin chirality $`\mathrm{sin}^2\theta \mathrm{sin}^2(\phi /2)0,`$ the spin conductance aprroaches to $`[\mathrm{sin}\rho \pi ]/\pi `$, which is proportional to the averaging kinetical energy of the charge carriers in this tight binding theory. From the numerical results in Figs. 1 and 2, we deduce that the spin conductance also depends on $`J`$. It becomes more complicated for the finite $`J`$ cases. Similarly, we can also introduce a spin voltage, which is also proportional to the length along the spin spiraling path, $`V_{spin}=𝐒_n\times 𝐒_{n+1}𝑑𝐥/S^2`$ where $`d𝐥=𝐫_{n+1}𝐫_n`$ is the difference of two position vectors of neighboring spins. $`V_{spin}`$ is an important parameter characterizing the ability to produce a spin current.

Thus far we have introduced several important concepts, such as spin chirality or spin current resistance, through the example of the classical spiral state. To determine the microscopic mechanism to produce a spin current, we go back to the Kondo lattice model (Eq.(1)) in the large Hund’s rule coupling limit and address the quantum spin case. The following discussions are not limited to one-dimensional systems and the conclusion is independent of dimensionality. According to the Hund’s rule coupling, the energy of empty or double occupancy of conduction electrons is zero. For a single occupancy the conduction electron and localized spin can form a spin $`S+1/2`$ state with $`JS/2`$ and a spin $`S1/2`$ state with $`+J(S+1)/2`$. In the large $`J`$ limit the spin $`S1/2`$ state and double occupancy should be excluded as they have much higher energy than the spin $`S+1/2`$ state. This process can be realized with the help of the projection technique.Fulde97  For that, we introduce the spinor $`\varphi _n^+=((1n_{n,})c_{n,}^{},(1n_{n,})c_{n,}^{}).`$ The dressed operators $`(1n_{n,\sigma })c_{n,\sigma }^{}`$, instead of $`c_{n,\sigma }^{}`$, are used here to avoid the double occupancy due to the strong coupling. The projection operator for the spin $`S+1/2`$ state is defined as

$$𝒫_n^+=\varphi _n^+𝐏_n^+\varphi _n,$$
(11)
where $`𝐏_n^+=[𝐒_n\tau +(S+1)\tau _0]/(2S+1)`$. The operator $`𝒫_n^+`$ forces the electron spin at site $`n`$ to be fully polarized with the localized spin at the same site. The projection operator for non-occupancy is $`𝒫_n^e=(1n_{n,})(1n_{n,}).`$ The total projection operator to exclude the spin S-1/2 single occupancy and double occupancy is $`𝒫=\underset{n}{}\left(𝒫_n^e+𝒫_n^+\right).`$ Therefore the equivalent Hamiltonian in the strong Hund’s rule coupling limit is expressed asKubo72 ; Shen98

$$H_{eff}=𝒫H𝒫=t\underset{nm}{}[\varphi _n^+𝐏_n^+𝐏_m^+\varphi _m+h.c.],$$
(12)
where a constant is omitted in the last step. In the large S approximation the Hamiltonian of (Eq.12) becomes the double exchange model in Eq.(7). Using a mathematical identity

$$\left(𝐒_n\tau \right)\left(𝐒_m\tau \right)=𝐒_n𝐒_m\tau _0+i\left(𝐒_n\times 𝐒_m\right)\tau ,$$
(13)
the Hamiltonian can be rewritten asNote

$`H`$ $`={\displaystyle \frac{\left(S+1\right)^2t}{\left(2S+1\right)^2}}{\displaystyle \underset{nm}{}}k_{nm}`$
$`{\displaystyle \frac{t}{\left(2S+1\right)^2}}{\displaystyle \underset{nm}{}}k_{nm}𝐒_n𝐒_m`$
$`{\displaystyle \frac{\left(S+1\right)t}{\left(2S+1\right)^2}}{\displaystyle \underset{nm}{}}\left(𝐒_n+𝐒_m\right)𝐌_{nm}`$
$`+{\displaystyle \frac{t}{(2S+1)^2}}{\displaystyle \underset{nm}{}}\left(𝐒_n\times 𝐒_m\right)𝐈_{nm}^s,`$ (14)
where $`k_{nm}=\varphi _n^+\tau _0\varphi _m+\varphi _m^+\tau _0\varphi _n`$ and $`𝐌_{nm}=\varphi _n^+\tau \varphi _m+\varphi _m^+\tau \varphi _n.`$ The first term is the kinetic part for conduction electrons excluding double occupancy, like the $`t`$-term in the $`tJ`$ model. The second term contains an effective Heisenberg-type spin-spin exchange interaction. Approximately, the effective coupling is $`J_{ferro}=tk_{nm}/(2S+1)^2`$, proportional to the kinetic energy in the sense of the mean field approximation. The kinetic energy $`tk_{nm}`$ in the tight binding approximation is always negative if the electron filling $`\rho <1`$ such that the effective coupling is ferromagnetic. This is consistent with the double exchange picture.Furukawa95 ; Millis95  The third term is the interaction between local spin $`𝐒_n`$ and the vector $`𝐌_{nm}`$. For a paramagnetic state, $`𝐌_{nm}`$ is equal to zero, but in a fully polarized state, its value is equal to $`k_{nm},`$ i.e. proportional to the kinetic energy. Among the four terms the most important observation is from the last term: it contains an interaction between the spin current operator $`𝐈_{nm}^s`$ and the spin chirality $`𝐒_n\times 𝐒_m`$. By analogy with the coupling between the electric current and electric field it is straightforward to understand why a non-zero spin chirality induces a pure spin current. Although we cannot derive an analytic expression for a finite $`J`$, it is expected that the mechanism to produce a spin current remains the same.

A new effective field and spin battery. To clarify the physical meaning of the last term, we come to calculate the rate of change of a local spin. From classical mechanics or from the Heisenberg equation of motion the rate of change of angular momentum Sn is equal to the torque $`\mu _n\times 𝐡_{eff}^n`$ which acts on the spin,

$$d𝐒_n/dt=g\mu _B𝐒_n\times 𝐡_{eff}^n,$$
(15)
where the effective field experienced by the spin is

$`𝐡_{eff}^n`$ $`={\displaystyle \frac{t}{g\mu _B(2S+1)^2}}{\displaystyle \underset{\delta }{}}`$
$`\left[k_{nn+\delta }𝐒_{n+\delta }+(S+1)𝐌_{nn+\delta }+𝐈_{nn+\delta }^s\times 𝐒_{n+\delta }\right].`$
The first two terms, again are related to particle’s orbital motions. The third term is the torque between the spin current and the local spin. A new type of force or effective field is caused due to the interaction between the spin current and local spin or local magnetization. The force may be responsible for magnetization precession by the spin injection into a magnetic conductorWeber01  and the anomalous Hall effect in ferromagnetic materials.Hirsch99 ; Niu01  Further discussion along this line will be published elsewhere. Our observation will lead to potential applications. For example, by means of the spin spiraling properties we can design a spin battery to produce a pure spin current. The battery contains two parts as shown in Fig.3. The main part is a spin spiraling material. The spins spirals about the axis which is perpendicular to the interface with a ferromagnetic metal film such that the spin voltage along the axis is non-zero. The spin chirality in the material will produce the spin accumulation at two terminals. The magnetization of the film is also perpendicular to the interface. The film is a source of spins which will flow out of the battery driven by the spin chirality. The key element in this battery is the spin spiraling materials. Possible realization of spiral states have been discussed by several authors,Kubo82 ; Jensen91 ; Inoue95 ; Korling96  and the state may also be stabilized with help of the external magnetic field and exchange coupling between localized spins.

In conclusion, we have shown that the strong coupling between the spins of charge carriers and local spins leads to an effective interaction between spin current and spin chirality in the system. Non-zero spin chirality may give rise to a pure spin current and in return a spin current can cause a novel type of force on the spins. We confirm these findings by studying in detail the spin transport in the spin spiral state.

The work is supported by grants from the Research Grants Council of Hong Kong, China (Project HKU7088/01P) and US-DOE.