text:
# Manuscript Title: with Forced Linebreak

## I Introduction

In the absence of a reionization, the photons of the Cosmic Microwave Background (CMB) are not scattered from decoupling to present time, except inside galaxy clusters (Sunyaev-Zel’dovich effect); nevertheless, the propagation direction of these photons changes due to the gravitational action of the cosmological inhomogeneities (lensing).

The so-called primary CMB anisotropies were produced by linear structures at high redshifts. The corresponding temperature contrast, $`\mathrm{\Delta }__P(\mathrm{\Delta }T/T)__P`$, is a statistical field which can be expanded in spherical harmonics as follows: $`\mathrm{\Delta }__P=_{m=\mathrm{}}^+\mathrm{}a_\mathrm{}mY_\mathrm{}m`$, where the $`a_\mathrm{}m`$ quantities are statistically independent variables depending on the sky realization we are expanding. If many sky realizations are averaged, the resulting $`a_\mathrm{}m`$ coefficients have zero means and variances $`a_\mathrm{}m^2=C_{\mathrm{}}`$. Since $`\mathrm{\Delta }__P`$ is a homogeneous and isotropic statistical field, quantities $`a_\mathrm{}m^2`$ do not depend on $`m`$. In the model under consideration, scalar energy density fluctuations are initially Gaussian, and they remain Gaussian during linear evolution, by this reason, the primary CMB anisotropies –which were produced by linear inhomogeneities– also are Gaussian.

Gaussian primary anisotropy appears to be superimposed to secondary anisotropies generated well after decoupling, and also to the anisotropy of the microwave radiation emitted by both our galaxy and some extragalactic sources (contaminant foregrounds). These foregrounds are sub-dominant non Gaussian components contributing to the total observable signal at microwave wavelengths. In the model we are considering here, there are various secondary anisotropies. A few comments about them are worthwhile: the late integrated Sachs-Wolfe anisotropy was produced by quasi-linear inhomogeneities ful00 (1) and, consequently, deviations with respect to Gaussianity are not expected to be very important; another secondary gravitational anisotropy was produced by strongly nonlinear structures as galaxy clusters and substructures (Rees-Sciama effect), this component is sub-dominant and non Gaussian, and its deviations with respect to Gaussianity were studied by various authors mol95 (2, 3, 4); finally, the Sunyaev-Zeldovich anisotropy was produced by hot plasma in galaxy clusters, and it is also non Gaussian coo01 (5).

The main goal of this paper is the analysis of the non Gaussianity caused by lensing, and such an analysis can be performed using the operative ray-tracing procedure designed in next sections. In the absence of lensing, there are directional frequency shifts which cause an anisotropic temperature distribution in the CMB. The lens effect does no cause frequency shifts, but it produces angular deviations of the propagation directions. This means that the lens effect does not produce any alteration of the CMB temperature. It changes the propagation direction of the CMB photons and, then, the angular distribution of the CMB temperature changes accordingly. In other words, the lens effect deforms the maps of the CMB temperature which appear as a result of pure frequency shifts. Any numerical estimate of the alterations of the CMB maps produced by lensing (including non Gausianity generation) involves three steps: (i) maps of the CMB temperature distribution are built up for the model under consideration (without any lensing), (ii) these maps are deformed taking into account the deviations of the propagation directions produced by lensing, and (iii) the resulting deformations of unlensed maps (hereafter called either lens deformations or lens effect) are analyzed. Let us now consider each of these steps in more detail.

In order to built up unlensed maps of the CMB, we assume that, in the absence of lensing, the Gaussian primary anisotropy dominates and, consequently, the maps to be deformed can be built up using the angular power spectrum of this primary anisotropy, which has been calculated using CMBFAST seza96 (6) for the $`\mathrm{\Lambda }`$CDM model under consideration. As it is discussed below, only small Gaussian maps having sizes of a few degrees and resolutions of the order of one arc-minute are necessary and they are built up using the Fourier Transform (FT) and the mentioned spectrum ($`C_{\mathrm{}}`$ quantities). See sae96 (7) and references cited therein for details. The resulting maps have spots with different sizes and amplitudes (angular structure). Spots with angular size $`\pi /\mathrm{}`$ have a mean amplitude proportional to $`C_{\mathrm{}}^{1/2}`$; hence, spots with angular size close to one degree ($`\mathrm{}200`$) have large amplitudes, whereas spots corresponding to one arcminute ($`\mathrm{}10000`$) have negligible amplitudes as a result of the smallness of $`C_{10000}`$. For the spectrum under consideration, significant spots have sizes of various arcminutes.

Cosmological inhomogeneities produce lensing and this effect is associated to the deviation field $`\stackrel{}{\delta }`$. The anisotropy observed in the $`\stackrel{}{n}`$ direction is the primary anisotropy corresponding to the direction $`\stackrel{}{n}_0`$ in the absence of lensing, where $`\stackrel{}{n}=\stackrel{}{n}_0\stackrel{}{\delta }`$; hence,

$$\mathrm{\Delta }(\stackrel{}{n})=\mathrm{\Delta }__P(\stackrel{}{n}_0)=\mathrm{\Delta }__P(\stackrel{}{n}+\stackrel{}{\delta }).$$
(1)
The unit vectors $`\stackrel{}{n}`$ and $`\stackrel{}{n}_0`$ point towards two points of the last scattering surface, and the deviation field $`\stackrel{}{\delta }`$ gives the angular excursion on this surface due to lensing. We must estimate the deviation field $`\stackrel{}{\delta }`$ (see below) to get the deformed observable map $`\mathrm{\Delta }`$. If the $`\stackrel{}{\delta }`$ deviations are smaller than the size of the smallest significant spots (various minutes), Eq. (1) can be expanded in $`\stackrel{}{\delta }`$ powers to get

$$\mathrm{\Delta }(\stackrel{}{n})\mathrm{\Delta }__P(\stackrel{}{n})+\frac{\mathrm{\Delta }__P}{\stackrel{}{n}}\stackrel{}{\delta }$$
(2)
We see that map deformations, namely, the differences in the temperature contrasts corresponding to associate directions (second term of the r.h.s.) are the product of two statistical fields: $`\mathrm{\Delta }__P`$ and $`\stackrel{}{\delta }`$ and, even if these field are Gaussian and statistically independent, the resulting product would be non Gaussian.

Several authors have studied angular correlations in lensed CMB maps. In the most theoretical papers on this subject, some correlations have been estimated without using simulations; the approach used in these papers only applies for those correlations which can be written in terms of the time varying form of the matter power spectrum. In the most basic of these papers sel96 (8), the mentioned spectrum was modelled beyond the linear regime and the second order correlations of lensed CMB maps were estimated. A similar method was used by Bernardeau ber97 (9) to estimate fourth order correlations corresponding to various sets of four directions (with distinct configurations); furthermore, this last author pointed out that third order correlations vanish. Other methods are necessary to estimate correlations of order $`n`$ for $`n>4`$, and also to simulate CMB maps with the true statistics induced by lensing, namely, maps having the true correlations at any order (not only at second and fourth orders). The most powerful technique to study lens distortions in CMB maps is the use of ray-tracing through N-body simulations. Two methods to apply this tecnique are described in references jai00 (10) and whi01 (11) (other relevant papers in this field are cited in these two references). Another new method is proposed here. We agree with White & Hu whi01 (11), who stated that: ….. On subdegree scales, a full description of week lensing therefore requires numerical simulations, the most natural being N-body simulations ……; in fact, using ray-tracing through N-body simulations, CMB photons are deviated by fully nonlinear structures (as galaxy clusters) which move in the simulation box. Nowadays, only N-body simulations lead to a proper description of a distribution of strongly nonlinear structures. Unfortunately, these simulations constraint us to work in a periodic universe, and relevant problems associated to periodicity must be solved. Each method –designed to use ray-tracing through N-body simulations– corresponds to a different way for preventing periodicity effects. In this paper, one of these methods is described and tested.

Hereafter, quantity h is the reduced Hubble constant $`h=10^2H_0`$, where $`H_0`$ is the Hubble constant in units of $`Kms^1Mpc^1`$, the density parameters corresponding to baryons, dark matter, and vacuum, are $`\mathrm{\Omega }_b`$, $`\mathrm{\Omega }_d`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }`$, respectively, the total density parameter is $`\mathrm{\Omega }=\mathrm{\Omega }_b+\mathrm{\Omega }_d+\mathrm{\Omega }_\mathrm{\Lambda }`$, and the matter density parameter is $`\mathrm{\Omega }_m=\mathrm{\Omega }_b+\mathrm{\Omega }_d`$. In the flat inflationary universe under consideration, the above parameters take on the following values: $`h=0.65`$, $`\mathrm{\Omega }_b=0.05`$, $`\mathrm{\Omega }_d=0.25`$, and $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$ and, then, according to Eke et al. eke96 (12), the power spectrum of scalar energy density perturbations must be normalized with $`\sigma _8=0.93`$ in order to have cluster abundances compatible with observations. Units are chosen in such a way that $`c=8\pi G=1`$, where $`c`$ is the speed of light and $`G`$ the gravitation constant. Whatever quantity ”$`A`$” may be, $`A__L`$ and $`A_0`$ stand for the $`A`$ values on the last scattering surface and at present time, respectively. The scale factor is $`a(t)`$, where $`t`$ is the cosmological time, and its present value, $`a_0`$, is assumed to be unity, which is always possible in flat universes.

## II Formalism

The CMB photons move on null geodesics and the line element is:

$$ds^2=(1+2\varphi )dt^2+(12\varphi )a^2\delta _{ij}dx^idx^j$$
(3)
where function $`\varphi `$ satisfies the equation:

$$\mathrm{\Delta }\varphi =\frac{1}{2}a^2(\rho \rho _{m0}),$$
(4)
and $`\rho _{m0}=\mathrm{\Omega }_{m0}\rho _{crit}`$ is the background energy density for matter. On account of this last equation, function $`\varphi `$ can be interpreted as the peculiar gravitational potential created by the cosmological structures. Equations (3) and (4) are valid for linear inhomogeneities located well inside the horizon and also for nonlinear structures (potential approximation), see mar90 (13).

The deviation field $`\stackrel{}{\delta }`$ is given by the following integral sel96 (8):

$$\stackrel{}{\delta }=2_{\lambda __L}^{\lambda _0}W(\lambda )\stackrel{}{}_{}\varphi 𝑑\lambda ,$$
(5)
where $`\stackrel{}{}_{}\varphi =\stackrel{}{n}\stackrel{}{n}\stackrel{}{}\varphi `$ is the transverse gradient of the potential, and $`W(\lambda )=(\lambda __L\lambda )/\lambda __L`$. The variable $`\lambda `$ is

$$\lambda (a)=H_0^1_a^1\frac{db}{(\mathrm{\Omega }_{m0}b+\mathrm{\Omega }_\mathrm{\Lambda }b^4)^{1/2}}.$$
(6)
The integral in the r.h.s. of Eq. (5) is to be evaluated along the background null geodesics. In our flat background, the equations of the null geodesics passing by point $`x__P^i`$ are:

$$x^i=x__P^i+\lambda (a)n^i,$$
(7)
In next section, spherical clusters with Navarro-Frenk-White (NFW) density profiles nav96 (14) are considered for qualitative analysis. These profiles appear in N-body simulations of dark matter halos and they have the form:

$$\rho (r)=\frac{\delta _c\rho _{m0}}{(r/r_s)[1+(r/r_s)]^2},$$
(8)
where $`\rho _{m0}=\rho _{crit}\mathrm{\Omega }_{m0}`$, and $`\rho _{crit}=3H^2`$ is the critical density of the universe. The parameters $`\delta _c`$ and $`r_s`$ are related to both the concentration parameter $`c_{_{NFW}}`$ and the mass $`M_{200}`$. This mass is that contained inside a sphere whose mean density is $`200`$. The radius of this sphere is denoted $`r_{200}`$. The following relations hold:

$$r_s=(3M_{200}/800\pi \rho _{m0})^{1/3}/c_{_{NFW}}$$
(9)
$$\delta _c=\frac{200c_{_{NFW}}^3}{ln(1+c_{_{NFW}})c_{_{NFW}}/(1+c_{_{NFW}})}.$$
(10)
The profile (8) remains almost unchanged from virialization ($`z1`$) to present time. Some authors have tried to assume the NFW profile with variable values of the involved parameters to go beyond virialization –until $`z5`$ in reference bul01 (15) –. Other authors use generalized profiles wyi01 (16). Fortunately, our qualitative analysis –presented in next section– does not require any detailed z-dependent profile.

In Sec. III, three clusters –which are hereafter called $`C__I`$, $`C_{_{II}}`$, and $`C_{_{III}}`$– are considered. Cluster $`C__I`$ is a rich one having $`M_{200}=1000M_g`$ (where $`M_g=10^{12}M_{}`$) and $`c_{_{NFW}}=6.8`$, cluster $`C_{_{II}}`$ is a standard cluster having a few cents of galaxies with $`M_{200}=250M_g`$ and $`c_{_{NFW}}=7.7`$, and cluster $`C_{_{III}}`$ is a galaxy group having a few tens of galaxies with $`M_{200}=25M_g`$ and $`c_{_{NFW}}=12.5`$. Since the NFW profile diverges as physical radius $`r`$ tends to zero, an uniform core of radius $`r_c`$ has been assumed and, consequently, the density profile has the form (8) for $`rr_c`$ and a constant value for $`rr_c`$. Both profiles match continuously at $`r=r_c`$. For this mass distribution, function $`d\varphi /dr`$ has the following form:

$$\frac{d\varphi }{dr}=\frac{1}{6}\rho (r_c)r$$
(11)
for $`rr_c`$, where $`\rho (r_c)`$ is given by Eq. (8), and

$$\frac{d\varphi }{dr}=\frac{1}{r^2}[A+Dr_s^3\{E(r)E(r_c)\}]$$
(12)
for $`rr_c`$, where $`A=(1/6)r_c^3\rho (r_c)`$, $`D=\delta _c\rho _{m0}/2`$, and $`E(r)=ln[1+(r/r_s)]+[1+(r/r_s)]^1`$. Core radius of $`0.2Mpc`$, $`0.1Mpc`$, and $`0.05Mpc`$ are assigned to clusters $`C__I`$, $`C_{_{II}}`$, and $`C_{_{III}}`$, respectively. In next section, our conclusions are probed to be highly independent on this assignation.

Spherical linear inhomogeneities evolving inside the effective horizon are also considered in Sec. III for qualitative discussion. These spherical regions are assumed to be uniform and they are characterized by their comoving radius $`R__C`$. This radius fixes the associated comoving scale $`k__C=\pi /2R__C`$ and, then, the present density contrast is fixed using the power spectrum of the energy density perturbations –in the model under consideration– which gives the average amplitude of the Fourier mode $`k__C`$. The density contrast evolves proportional to the growing mode D(t), which is given in pee80 (17). Three of these linear inhomogeneities will be considered, they correspond to $`R__C=250Mpc`$, $`R__C=50Mpc`$, and $`R__C=12.5Mpc8h^1Mpc`$, and they are denoted $`L__I`$, $`L_{_{II}}`$, and $`L_{_{III}}`$, respectively. The derivative of the gravitational potential with respect to the physical radius $`r`$ is

$$\frac{d\varphi }{dr}=\frac{1}{6}(\rho \rho _{m0})r$$
(13)
for $`raR__C`$, and

$$\frac{d\varphi }{dr}=\frac{a^2R__C^2}{r^2}\frac{d\varphi }{dr}(aR__C)$$
(14)
for $`r>aR__C`$.

## III Angular scales, boxes and resolution

Recently, various groups have used simulations of structure formation to estimate the effects of lensing on the CMB. An important problem is that there are many scales producing lensing and that, according to White & Hu whi01 (11), simulating the full range of scales implied is currently a practical impossibility. These authors proposed the tiling method to circumvent this problem and also the problems with periodicity (see Sec. I); this method employs many independent simulations with different boxes and resolutions to tile the photon trajectories. White & Hu whi01 (11) explained the advantages of their method with respect to the most traditional one based on plane projections, which has been extensively used in the literature (see jai00 (10) and references cited therein). The method proposed here does not require many independent simulations –but only one– and, consequently, it is simpler and numerically faster than the tiling one. Furthermore, it does not use plane projections whose potential problems were pointed out in whi01 (11). In this section, we are concerned with the scales relevant for CMB lensing, whereas periodicty is considered in next sections.

We are interested in the deviations with respect to Gaussianity produced by lensing jai00 (10). These deviations can be measured by angular correlations corresponding to three and more directions. Given an angular scale, the correlations are produced by a set of linear and nonlinear cosmological structures and, a detailed study of the type of structures contributing to these correlations is necessary in order to design our N-body simulations, which should include –in the same box– all the structures responsible for the effects we are looking for. In order to estimate scales, we can consider the evolving linear inhomogeneities $`L__I`$, $`L_{_{II}}`$, and $`L_{_{III}}`$, and the nonevolving NFW density profiles $`C__I`$, $`C_{_{II}}`$, and $`C_{_{III}}`$. Cluster evolution would not affect our qualitative estimates for small redshifts $`z<5`$, although it could be crucial in other contexts. In short, our study is based on six appropriate structures, three linear ones and three galaxy clusters.

First of all, Eq. (5) and the potentials given in Sec. II have been used to calculate the lens deviations produced by the six selected structures. Each linear structure is placed at redshifts 2, 10 and 100, whereas each cluster is located at redshifts 0.5, 2 and 5 (cluster $`C_{II}`$ is also placed at $`z=0.2`$ and $`z=0.1`$). The deviation $`\delta `$ is calculated for each observation direction, which is characterized by the angle $`\theta `$ formed by the line of sight and the direction pointing towards the structure center. Results are presented in Fig. 1, where we see that: (i) the deviations grow as the redshift decreases, and the growing is small for $`z0.5`$ (compare the curves of the left middle panel of Fig. 1 corresponding to redshifts 0.5, 0.2 and 0.1), (ii) the deviations increase as the inhomogeneity size increases (this size is fixed by $`M_{200}`$ for clusters and by $`R__C`$ for linear structures), (iii) small linear structures ($`L_{_{III}}`$) produces very small deviations whatever their location may be, (iv) small clusters (or groups, $`C_{_{III}}`$) produce a maximum deviation which is about $`10`$ % of the maximum deviation produced by cluster $`C__I`$ (located at the same redshift), and (v) the most important deviations are produced by massive clusters and by very extended linear inhomogeneities.

We are interested in the angular correlations of some CMB maps. Given a map of the variable $`\zeta `$ and $`m`$ directions, these correlations are defined as follows:

$$C_m=\zeta (\stackrel{}{n}_1)\zeta (\stackrel{}{n}_2)\mathrm{}\zeta (\stackrel{}{n}_m),$$
(15)
where the average is over many realizations of the CMB sky. Some of these averages are estimated below for two, three and four directions in maps of both primary anisotropy and lens deformations. For $`m=4`$, the chosen directions draw the vertices of a square on the last scattering surface with $`\stackrel{}{n}_1\stackrel{}{n}_2=\mathrm{cos}\alpha `$, and for $`m=3`$, they point towards the vertices of an isosceles rectangle triangle with $`\stackrel{}{n}_1\stackrel{}{n}_2=\mathrm{cos}\alpha `$ and $`\stackrel{}{n}_1\stackrel{}{n}_3=\mathrm{cos}\alpha `$. The relevance of the deviations presented in the panels of Fig. 1 depends on the angular scale, $`\alpha `$, chosen for correlation computations. Suppose the scale $`\alpha `$ ($`\mathrm{}=\pi /\alpha `$) and $`m=2`$. From the primary temperatures corresponding to many pairs of directions ($`\stackrel{}{n}_1`$, $`\stackrel{}{n}_2`$) forming angle $`\alpha `$, the average $`\mathrm{\Delta }__P(\stackrel{}{n_1})\mathrm{\Delta }__P(\stackrel{}{n_2})`$ can be calculated and the result measures correlations in the absence of lensing. According to Eq. (1), after lensing, the correlations are $`\mathrm{\Delta }(\stackrel{}{n_1})\mathrm{\Delta }(\stackrel{}{n_2})=\mathrm{\Delta }__P(\stackrel{}{n_{01}})\mathrm{\Delta }__P(\stackrel{}{n_{02}})`$ and, consequently, the correlations would change –due to lensing– if, as a result of deviations, the directions $`\stackrel{}{n}_{01}`$ and $`\stackrel{}{n}_{02}`$ form an angle $`\alpha _0`$ different enough from $`\alpha `$; therefore, in order to see if a certain spherical inhomogeneity located at redshift $`z`$ can contribute to the correlation at the angular scale $`\alpha `$, we may consider two directions pointing towards $`R`$ and $`S`$, two points of a certain diameter of the spherical structure; the first (second) of these directions forms angle $`\theta `$ ($`\theta +\alpha `$ ) with the direction of the symmetry center, and the angle $`\mathrm{\Theta }=\theta +\alpha /2`$ corresponds to the middle point of the segment $`RS`$; then, the differences $`d(\mathrm{\Theta })=\delta (\theta +\alpha )\delta (\theta )`$ measure the deformations of the angle formed by the two chosen directions. The top panels of Fig. 2 give these differences for the structures $`C__I`$, $`C_{_{II}}`$, and $`C_{_{III}}`$ (at z=0.5), and the bottom panels for $`L__I`$ and $`L_{_{II}}`$ (at z=2). The correlation scale is $`\alpha =13.5^{}`$ ($`\mathrm{}=800`$) in the left panels and $`\alpha =6.75^{}`$ ($`\mathrm{}=1600`$) in the right ones. From the analysis of this Figure it follows that: (a) big close clusters of type $`C__I`$ produce the largest deformations reaching the maximum value $`1.68^{}`$ ($`1.98^{}`$) for $`\alpha =13.5^{}`$ ($`\alpha =6.75^{}`$), (b) for $`\alpha =13.5^{}`$, inhomogeneities like $`C_{_{II}}`$, $`C_{_{III}}`$, $`L__I`$, and $`L_{_{II}}`$ produce maximum deformations which are about $`37.5\%`$, $`7.7\%`$, $`33\%`$ and $`6.8\%`$ of the maximum ones corresponding to $`C__I`$ clusters ($`1.68^{}`$) and, (c) for $`\alpha =6.75^{}`$, the maximum deformations produced by $`C_{_{II}}`$, $`C_{_{III}}`$, $`L__I`$, and $`L_{_{II}}`$ structures are about $`36.9\%`$, $`7.1\%`$, $`13.7\%`$ and $`2.9\%`$ of the maximum deformation appearing in case $`C__I`$ ($`1.98^{}`$). Points (b) and (c) lead to the conclusion that the deformations produced by big linear structures are more important for $`\alpha =13.5^{}`$, but they are not negligible in the case $`\alpha =6.75^{}`$. The deformations corresponding to structures like $`L_{_{III}}`$ are too small and they have not been presented in Fig. 2.

From the above discussion, it follows that –for the scales under consideration– the most relevant structures are clusters of type $`C__I`$ and $`C_{_{II}}`$ at low redshifts, and big linear inhomogeneities with $`R__C>50Mpc`$ located at small redshifts. This conclusion is important. It tell us that simulations which do not resolve galaxy groups and small clusters can lead to good estimations, and also that the box size can be chosen to prevent the existence of extended linear inhomogeneities contributing to lensing. Boxes with sizes of $`128Mpc`$ should be appropriated, sizes of $`256Mpc`$ could be acceptable (although a part of the lensing could be due to linear structures with $`R__C>50Mpc`$, see Sec. V) and, finally, boxes of $`512Mpc`$ and greater could contain an extended part of a $`L__I`$ structure giving important contributions to lensing. The idea is that the lens effect produced by big linear inhomogeneities can be studied without simulations, see references in jai00 (10), whereas the effect of any other significant structure ($`C__I`$ and $`C_{_{II}}`$ clusters) can be included in simulations for appropriate boxes (between $`128Mpc`$ and $`256Mpc`$ size). Let us now analyze if results from Figs. 1 and 2 are robust against variations of the core radius $`r_c`$.

Imagine one of the above three clusters located at redshift z. It produces a significant deviation of the CMB photons if and only if these photons cross the region where the function to be integrated in Eq. (5) (proportional to $`\stackrel{}{}_{}\varphi `$) is not negligible; hence, if the dependence of this function on the core radius is weak, lens deviations (integral in Eq. (5)) also depend weakly on $`r_c`$, and previous results based on Figs. 1 and 2 are robust. In order to study the mentioned function for clusters $`C__I`$, $`C_{_{II}}`$, and $`C_{_{III}}`$, these structures are located at redshift 0.5 (where cluster evolution is not expected to be important) and, then, the function $`d\delta (\theta )/d\lambda =2W(\lambda )\stackrel{}{}_{}\varphi `$ involved in Eq. (5) is calculated at points of the null geodesic corresponding to the angle $`\theta =2.5^{}`$. Results are presented in Fig. 3, where function $`d\delta (2.5)/d\lambda `$ is given (with arbitrary normalization) in terms of the comoving distance to the cluster center $`\lambda \lambda _c`$, where $`\lambda `$ is the comoving distance from the observer to an arbitrary point of the null geodesic, and $`\lambda _c`$ is the $`\lambda `$ value corresponding to the cluster center. Continuous lines display the results obtained for clusters $`C__I`$, $`C_{_{II}}`$, and $`C_{_{III}}`$ with the radius cores we have previously fixed, whereas the associated dotted lines correspond to the same clusters with smaller cores, the new core radius being $`0.1Mpc`$, $`0.05Mpc`$, and $`0.025Mpc`$ in cases $`C__I`$, $`C_{_{II}}`$, and $`C_{_{III}}`$, respectively (the initially chosen radius have been reduced by a factor $`1/2`$).

The continuous and dotted associated lines of Fig. 3 are very similar for any cluster, which implies that, for $`\theta =2.5^{}`$, the lens deviation $`\delta `$ has a very weak dependence on the core radius $`r_c`$. The direction $`\theta =2.5^{}`$ has been arbitrarily chosen in the $`\theta `$ interval where the lens deformations are relevant (see Fig. 1); nevertheless, other angles covering this interval have been also studied. For example, we have investigated the directions corresponding to the maximum deviation in each of the cases reported in Fig. 1, and lens deviations also depends weakly on the core radius. The same occurs for any $`\theta `$ value, except for very small angles, which are not relevant for the angular scales under consideration (see below). In Fig. 3, we can also see that the region where $`\stackrel{}{}_{}\varphi `$ is significantly contributing to the integral of Eq. (5) has a comoving size of various Megaparsecs for any cluster of interest. It has been checked that the same conclusion is also valid for any relevant $`\theta `$ values.

A scale for the computation of angular correlations has been chosen. In order to do that, the CMB angular power spectrum ($`C_{\mathrm{}}`$ coefficients) due to lensing by nonlinear structures has been estimated using CMBFAST. Results lead to the conclusion that $`C_2(\alpha )`$ correlations from this lensing are important for $`\alpha 13.5^{}`$ ($`\mathrm{}800`$); hence we must work with some scale $`\mathrm{}>800`$. The scale $`\mathrm{}=1600`$ ($`\alpha =6.75^{}`$) has been chosen to calculate correlations because this scale has the following properties: (a) for $`\mathrm{}=1600`$, the effect of linear inhomogeneities –which should be reduced as much as possible– is smaller than that of the case $`ł=800`$ (compare the bottom panels of Fig. 2), (b) according to CMBFAST estimations, lensing from nonlinear structures gives a $`C_{1600}`$ coefficient which is smaller than the maximum $`C_{\mathrm{}}`$ only by a factor $`0.6`$, and (c) the value $`\mathrm{}=1600`$ is in the range of the multipoles to be observed by PLANCK (the most sensitive projected experiment for CMB anisotropy measurements).

According to previous comments, a box of $`128Mpc`$ have been assumed as the basic one and, then N-body simulations have been performed inside this box with a standard PM code hoc88 (18), which was used, tested, and described in detail in qui98 (19). The resolution has been succesively increased starting from a very poor one of $`2Mpc`$. Of course, such a low resolution leads to non-linear structures whose amplitudes (sizes) are smaller (larger) than those of the true clusters. These simulations only traces a good spatial distribution of extended structures with total masses as those of the clusters (which would become clusters for appropriate resolutions). Resolutions of $`1Mpc`$ and $`0.5Mpc`$ have been also considered, but no higher resolutions have been used by the reasons given in Secs. V and VI.

## IV Departures from Gaussianity

Suppose an observer which is located at a certain point with comoving spatial coordinates $`x__P^i`$. The main question is: which is the deviation, $`\stackrel{}{\delta }(\stackrel{}{x}__P,\stackrel{}{n})`$, observed from point $`\stackrel{}{x}__P`$ in the $`\stackrel{}{n}`$ direction?

If $`\varphi _\stackrel{}{k}(t)`$ is the FT of the potential $`\varphi `$, namely,

$$\varphi (\stackrel{}{x},t)=\frac{1}{(2\pi )^{3/2}}d^3ke^{i\stackrel{}{k}\stackrel{}{x}}\varphi _\stackrel{}{k}(t),$$
(16)
and $`\delta _\stackrel{}{k}`$ is the FT of the density contrast $`(\rho \rho _{m0})/\rho _{m0}`$, Equation (4) leads to the following relation in Fourier space: $`\varphi _\stackrel{}{k}=B\delta _\stackrel{}{k}/k^2`$, where $`B=\rho _{m0}/a`$.

Using elemental Fourier algebra and Eqs. (5), and (7), the following basic equations –giving the required deviation– are easily obtained:

$$\stackrel{}{\delta }(\stackrel{}{x}__P,\stackrel{}{n})=\frac{2i}{(2\pi )^{3/2}}\frac{\stackrel{}{k}_{}}{k^2}F_\stackrel{}{k}(\stackrel{}{n})e^{i\stackrel{}{k}\stackrel{}{x}__P}d^3k,$$
(17)
where

$$F_\stackrel{}{k}(\stackrel{}{n})=_{\lambda __L}^{\lambda _0}W(\lambda )B(\lambda )e^{i\lambda \stackrel{}{k}\stackrel{}{n}}\delta _\stackrel{}{k}(\lambda )𝑑\lambda ,$$
(18)
and $`\stackrel{}{k}_{}=\stackrel{}{k}(\stackrel{}{n}\stackrel{}{k})\stackrel{}{n}`$.

According to Eq. (17), each component of $`\stackrel{}{\delta }(\stackrel{}{x}__P,\stackrel{}{n})`$ is the FT –extended to all the space– of a component of vector $`\stackrel{}{k}_{}F_\stackrel{}{k}(\stackrel{}{n})/k^2`$. The integral of the r.h.s. of Eq. (18) must be estimated to get the function $`F_\stackrel{}{k}(\stackrel{}{n})`$ involved in the FT. The integration variable $`\lambda `$ can be seen as a generalized time coordinate. Equations (17) and (18) are general and they describe how the photons deviate –following null geodesics– in a realization of the full universe; nevertheless, such a realization is not available in practice, and we are constrained to use these equations in a fictitious periodic universe, let us now discuss this fact and its consequences in detail.

Structure evolution is simulated with a certain N-body code (see Sec. III), which involves a box and a certain resolution (a network in position space). This code has a certain time step (or $`\lambda `$ step) and, consequently, it gives data at a set of time values; among these data, we are particularly interested in the function $`\delta _\stackrel{}{k}(\lambda )`$ which appears in Eq. (18). This function is easily evaluated at each node, $`\stackrel{}{k}`$, in Fourier space, and at each time step $`t_i`$, but not at arbitrary times. The question is: Can we calculate the integral (18) using only the data corresponding to the time discretisation of the N-body code? Fortunately, the answer is positive. We have studied the function $`\delta _\stackrel{}{k}(\lambda )`$ for many $`\stackrel{}{k}`$ values and, thus, we have verified that the term $`I(\stackrel{}{k},\lambda )=W(\lambda )B(\lambda )\delta _\stackrel{}{k}(\lambda )`$ –which appears in Eq. (18)– can be very well approximated by a straight line between $`\lambda _i`$ and $`\lambda _{i+1}`$, namely, between two successive time steps of the N-body simulation. The equation of this line is of the form:

$$I_i(\stackrel{}{k},\lambda )=A_i(\stackrel{}{k})\lambda +B_i(\stackrel{}{k})$$
(19)
where quantities $`A_i`$ and $`B_i`$ are obtained from $`\delta _\stackrel{}{k}(\lambda _i)`$ and $`\delta _\stackrel{}{k}(\lambda _{i+1})`$; namely, from quantities given by the N-body in two successive steps. If Eq. (19) is substituted into Eq. (18), the integral of the r.h.s. can be analytically calculated in the interval ($`t_i,t_{i+1}`$). The addition of these integrals for all the N-body intervals gives function $`F_\stackrel{}{k}(\stackrel{}{n})`$. After this integration, the FT (17) can be performed and, thus, the deviation field in the $`\stackrel{}{n}`$ direction is simultaneously calculated for many observers located at the nodes $`\stackrel{}{x}__P`$ of the spatial Fourier box. The importance of this fact is discussed below.

Discretisation implies that the CMB photons do not move in a true realization of the universe, but in a periodic universe, which is, at each time, an ensemble of boxes identical to that of the N-body simulation. This periodicity should lead to errors in our estimations. In order to understand the problem and its solution, let us establish a two-dimensional (2D) analogy. The squares of Fig. 4 (boxes in the 3D case) have a size of $`128Mpc`$. These squares contain small circles (clusters), and let us imagine an observer, S, at the center of the right square. If a photon propagates along the pointed line of Fig. 4, it crosses all the squares in the same way, namely, it passes at the same distances from the same clusters and undergoes the same systematic repeated deviations. Only evolution induces small deformations from square to square, but these deformations would be only relevant after a large enough number of squares have been crossed; hence, a certain wrong accumulative effect should appear. However, if the photon moves following the dashed line of Fig. 4, the situation is very different because the photon enters each square around a very different place and, then, it can be affected by distinct independent clusters. The key point is the selection of the angle $`\mathrm{\Phi }`$ and its associate distance $`L`$. Two facts are relevant to make this choice: (i) clusters significantly deviate the CMB photons when their impact parameter is smaller than a few Megaparsecs (see Fig. 3), and (ii) correlations in the cluster distribution (superclusters) extend near $`10Mpc`$; hence, if we choose $`L`$ between $`15Mpc`$ and $`20Mpc`$, photons move from square to square, but they are always influenced by different almost uncorrelated clusters and, consequently, the effect of periodicity must be minimum. All these considerations are easily generalized to the 3D case, in which there are two angles defining the line of sight. Hereafter, these angles are identified with the spherical coordinates $`\theta `$ and $`\varphi `$ (with respect to the Cartesian axis defined by the edges of the box). These angles can be easily chosen to maximize the length of the path where the photons undergo the action of distinct uncorrelated clusters; namely, to minimize the effects of periodicity. In the 3D case, with boxes of $`128Mpc`$, it has been found that for $`\theta =10.75^{}`$ and $`\varphi =41.19^{}`$, the photons would cross 56 boxes before arriving to the initial region of the box. The CMB photons moving along this direction (hereafter denoted $`\stackrel{}{N}_0`$ and called the preferred direction) do not feel periodicity during their travel from $`r7400Mpc`$ to $`r=0`$, namely, from $`z=3.62`$ to $`z=0`$ and, taking into account that the lens effect we are evaluating is produced at these low redshifts by galaxy clusters, it follows that, for direction $`\stackrel{}{N}_0`$ and close ones, periodicity must be rather irrelevant (see Sec. V for more discussion). Furthermore, the photon enters in two successive boxes through points that –when placed in a unique box– are located at a distance $`L=15.9Mpc`$. For the chosen values of $`\theta `$ and $`\varphi `$ and a box of $`256Mpc`$, we have $`L=31.8Mpc`$, and photons does not feel periodicity at all from decoupling to present time.

On the angular scales under consideration, the cluster distribution producing lensing is statistically independent on the structure distribution causing the primary anisotropy; in fact, for $`\alpha `$ values between $`6.75^{}`$ and $`13.5^{}`$ ($`1600\mathrm{}800`$), primary anisotropy is due to temperature fluctuations previous to decoupling, which appeared as a result of the coupling between matter an radiation in the presence of energy density fluctuations; hence, the scales of temperature fluctuations and those of energy perturbations coincide. Inhomogeneities with a comoving size of $`41Mpc`$ ($`20.5Mpc`$) located at $`z=1100`$ subtend an angle of $`13.5^{}`$ ($`6.75^{}`$); hence, taking into account that –in our $`\mathrm{\Lambda }`$CDM model– the comoving scale of an Abell cluster with $`10^{15}M_{}`$ is $`25.6Mpc`$, we can conclude that, for $`800\mathrm{}1600`$, the primary anisotropy is produced by comoving scales similar to those of clusters. These small structures located close to our last scattering surface ($`z1100`$) and the clusters producing lensing (small $`z`$ values) are very far, different, and independent structure distributions; therefore, we can consider that the distribution of clusters is producing deformations on independent primary CMB maps. In order to calculate the averages of Eq. (15) (correlations), various full realizations of the cluster distribution should be considered to produce deviation fields and, then, each of these fields should be used to deform a large enough number of independent primary maps. In practice (computational limitations), only a few pairs of directions (close to $`\stackrel{}{N}_0`$ to minimize periodicity effects) are studied for each cluster distribution (N-body simulation), and these directions should be used to deform different primary maps. Let us now describe in detail the method implemented in this paper to calculate correlations.

For each N-body simulation, nine directions $`\stackrel{}{n}_i`$ are considered. These directions are chosen in such a way that they depict the vertices of four neighboring squares with size $`\alpha `$ on the last scattering surface. These squares form another bigger one with size $`2\alpha `$ with vector $`\stackrel{}{N}_0`$ pointing towards its center. It is important that, for each of the chosen directions, the deviations $`\stackrel{}{\delta }(\stackrel{}{x}__P,\stackrel{}{n}_i)`$ are calculated for the $`N^3`$ observers, $`\stackrel{}{x}__P`$, at the same time (see above). Let us discuss the importance of this fact from the statistical point of view. Put a $`4.25^{}\times 4.25^{}`$ map, $`M__P`$, of primary anisotropy in such a way that, from the box center (point $`\stackrel{}{x}__C`$), the direction $`\stackrel{}{N}_0`$ points towards the map center $`C`$. Consider now one of the observation direction $`\stackrel{}{n}_i`$. From the box center, the chosen direction points towards a certain point $`C_i`$; nevertheless, from another point of the box $`\stackrel{}{x}__Q`$, direction $`\stackrel{}{n}_i`$ does not point towards point $`C_i`$ but towards another point $`Q_i`$ of the map $`M__P`$. Two parallel lines –with direction $`\stackrel{}{n}_i`$– starting from points $`\stackrel{}{x}__C`$ and $`\stackrel{}{x}__Q`$ of the box would intersect the map $`M__P`$ in points $`C_i`$ and $`Q_i`$, respectively; hence, the $`Q_i`$ coordinates in the map $`M__P`$ can be calculated as follows: project the vector $`\stackrel{}{x}__Q\stackrel{}{x}__C`$ on the plane normal to $`\stackrel{}{n}_i`$, translate this projection parallely from the hypersurface $`t=0`$ to the last scattering surface at $`z=1100`$, calculate the angle subtended by the final vector (as it is seen by the observers) and, then, the angular coordinates of point $`Q_i`$ in the map $`M__P`$ can be trivially obtained. Of course, the deviation, $`\stackrel{}{\delta }(\stackrel{}{x}__Q,\stackrel{}{n}_i)`$, corresponding to the observer $`\stackrel{}{x}__Q`$ in the direction $`\stackrel{}{n}_i`$, must be applied to point $`Q_i`$ in the map $`M__P`$. Hence, for each direction $`\stackrel{}{n}_i`$, we have not a unique point in $`M__P`$, but many points, one for each $`\stackrel{}{x}__Q`$ observer, and these points cover a certain region on map $`M__P`$. In order to estimate the size of the covered region, we can consider a fixed direction, for example $`\stackrel{}{n}=(0,0,1)`$, and a box of $`128Mpc`$ ($`256Mpc`$) per edge; thus, given two observers $`\stackrel{}{x}_{_{Q1}}`$ and $`\stackrel{}{x}_{_{Q2}}`$ separated by a distance of $`128Mpc`$ ($`256Mpc`$), the angular distance between points $`Q_1`$ and $`Q_2`$ in the map $`M__P`$ appears to be $`\mathrm{\Xi }=0.5^{}`$ ($`\mathrm{\Xi }=1^{}`$), which means that, for $`\alpha =13.5^{}`$, a squared region with a size of $`43.5^{}`$ ($`73.5^{}`$) is covered, and for $`\alpha =6.75^{}`$, the size of the corresponding region is $`36.75^{}`$ ($`66.75^{}`$); hence, we see that, for boxes of $`128Mpc`$ ($`256Mpc`$) and $`\alpha =6.75^{}`$, the size of this region is larger than $`\alpha `$ by a factor close to five (ten) and, consequently, the averages in the $`\stackrel{}{x}__Q`$ observers –a novel aspect of our ray-tracing procedure– play a crucial statistical role. The greater the boxes and the smaller the angular scales $`\alpha `$, the larger the above factor and the greater the statistical significance of the $`\stackrel{}{x}__Q`$ averages.

On account of the above considerations, correlations are calculated as follows: (i) twelve different pairs of directions, sixteen triads, and four tetrads are defined using the above nine directions, (ii) given an observer $`\stackrel{}{x}__Q`$, a map $`M__P`$ (primary anisotropy), and a N-body simulation (for $`\stackrel{}{\delta }`$ computation), the deviations $`\stackrel{}{\delta }(\stackrel{}{x}__Q,\stackrel{}{n}_i)`$ and Eq. (1) (or Eq. (2)) are used to calculate the lens effect $`\mathrm{\Delta }[\stackrel{}{n}_i+\stackrel{}{\delta }(\stackrel{}{x}__Q,\stackrel{}{n}_i)]\mathrm{\Delta }[\stackrel{}{n}_i]`$ (or $`\stackrel{}{\delta }(\stackrel{}{x}__Q,\stackrel{}{n}_i)[\mathrm{\Delta }__P/\stackrel{}{n}]_{\stackrel{}{n}_i}`$) for each of the $`\stackrel{}{n}_i`$ directions and, then, the resulting data are used to compute the averages (15) for the above pairs ($`m=2`$), triads ($`m=3`$), and tetrads ($`m=4`$), (iii) calculations in (ii) are repeated for each observer $`\stackrel{}{x}__Q`$, and results are averaged again, (iv) all the process is repeated for a certain number of $`M__P`$ maps and a new average is performed and, (v) finally, the above calculations are repeated for a certain number of N-body simulations, and the final averages are taken as our estimations of the correlations (15). The question is: How many primary maps and how many cluster realizations would be necessary to calculate the required correlations?

In order to answer this question, the following method has been implemented: In a first step, the process (i)-(iv) has been used to analyze pure primary maps in the absence of lens deviations. The primary maps are not deformed and, consequently, the resulting averages should approach the true two, three, and four direction correlations of these maps, which are hereafter called primary correlations. These correlations can be estimated by a direct analysis of the simulated primary maps (see next section) and the results of such an analysis can be compared with the correlations estimated by the process (i)-(iv) in the absence of deviations. Since no deviations have been considered, this study measures the capability of our method (based on nine direction, $`N^3`$ observers, and various primary maps) to create pairs, triads and tetrads covering a significant part of the primary maps; namely, allowing a good statistical analysis of maps. In the second and final step, the number of primary maps suggested by our previous study (in the absence of deviations) is taken and, then, deviations from more and more N-body simulations are considered in order to compute the required correlations as it has been described above. When the averages reach almost stable values, the process is stopped and no more N-body realizations are performed. Results are presented in next section.

## V Results

First of all, the methods applied –in this paper– to calculate correlations of the lens effect are tested. In order to do that, these methods are used to analyze well known maps of primary anisotropy. The angular scale for correlation estimates is $`\alpha =6.75^{}`$. A set of two hundred maps of primary anisotropy has been used. Results are presented in Fig. 5, where variable $`N_p`$ (horizontal axis) is the number of primary maps used to get the correlation value appearing in the vertical axis. Left, central, and right panels correspond to the two ($`m=2`$), three ($`m=3`$) and four ($`m=4`$) direction correlations (15), respectively. In all the panels, the solid line shows the correlations obtained by means of an exhaustive coverage of the $`M__P`$ maps using many pairs, triads and tetrads of directions. This is our best estimate of the required correlations; results from this method –which cannot be applied to analyze the lens component– approach the theoretical values of the correlations we have used to simulate the $`M__P`$ maps. Other methods must give similar correlations to be acceptable. As $`N_p`$ increases, the solid lines of the top and bottom panels seem to tend to a certain value, whereas the continuous line of the middle panel seems to be compatible with a vanishing $`C_3`$ correlation. The long dashed lines display the correlations obtained using nine directions and only one observer (only nine directions in each map). These lines slowly approaches the solid one as $`N_p`$ increases and the level of approximation seems to be good for $`N_p=200`$. The dotted (dashed) lines give the correlations obtained with nine directions and $`128^3`$ ($`256^3`$) observers. As it follows from Fig. 5, dashed and dotted curves approach the solid ones faster than the long dashed lines and, consequently, the use of many observers located in the nodes of a network (see Sec. IV) is statistically significant; the greater the box size, the faster the statistical approximation to the solid lines. The best situation corresponds to $`256^3`$ observers in a $`256Mpc`$ box. Two hundred primary maps suffice for a good enough estimate of $`C_2(6.75)`$ and $`C_4(6.75)`$; particularly, when $`N^3`$ observers are considered.

After computing the correlations $`C_2(6.75)`$, $`C_3(6.75)`$, and $`C_4(6.75)`$ for lens deformations, it would be worthwhile a comparison between the correlation level of these deformations and that of primary maps. The ratios $`r_2=C_2(6.75)/C_2(0)`$, $`r_3=C_3(6.75)/C_2^{3/2}(0)`$, and $`r_4=C_4(6.75)/C_2^2(0)`$ are apropriate to do this comparison; these quantities will be calculated for lens deformations and primary maps and, then, results will be compared. We begin with the case of primary maps. Using two hundred of these maps and exhaustive coverage, we have found: $`C_1=\mathrm{\Delta }__P(\stackrel{}{n})=2.12\times 10^6`$, $`C_2(0)=\mathrm{\Delta }__P(\stackrel{}{n})\mathrm{\Delta }__P(\stackrel{}{n})=7.29\times 10^8`$, and $`r_1=C_1/C_2^{1/2}(0)=7.85\times 10^3`$. We know that there are many positive and negative temperature contrasts –in the maps– with absolute values close to $`C_2^{1/2}(0)`$ and, consequently, only an strong cancellation of these contrasts can explain the resulting value of the mean contrast $`C_1(0)`$, which is much smaller than $`C_2^{1/2}(0)`$ (small ratio $`r_17.9\times 10^3`$). Strong cancellation leading to a small $`r_1`$ value indicates that, as the number of maps increases, quantity $`\mathrm{\Delta }__P(\stackrel{}{n})`$ approaches either a small value or zero. Similarly, in order to calculate the correlation $`C_2(6.75)=\mathrm{\Delta }__P(\stackrel{}{n_1})\mathrm{\Delta }__P(\stackrel{}{n_2})`$, positive and negative numbers with absolute values close to $`C_2(0)`$ are averaged. The resulting ratio $`r_20.85`$ indicates that these numbers have not cancelled among them and a nonvanishing $`C_2(6.75)`$ correlation exists. For $`m=3`$ (m=4), the relevant ratio is $`r_38.8\times 10^3`$ ($`r_42`$) and, this small (large) ratio $`r_3`$ ($`r_4`$) strongly suggests a vanishing $`C_3`$ (significant $`C_4`$) correlation.

Since the simulated primary maps are Gaussian, the relations $`C_3(6.75)=0`$ and $`C_4(6.75)=3C_2^2(6.75)`$ must be satisfied. The correlation $`C_3`$ obtained by exhaustive coverage (solid lines of Fig. 5) is $`C_31.76\times 10^{13}`$, and it has been proved to be small by computing $`r_3`$; furthermore, we have used the same method to get the correlations $`C_26.24\times 10^8`$ and $`C_41.09\times 10^{14}`$, which lead to the ratio $`3C_2^2/C_41.07`$; hence, we can say that the correlations extracted from the maps are compatible with Gaussianity, as it should be for the primary maps we are analyzing. If we use other methods to analyze the $`M_p`$ maps; for example, the methods used to built up the dotted, dashed, and long dashed curves of Fig. 5 (see above), the extracted correlations appear to be similar to those reported above and so they are compatible with Gaussianity; this means that all the methods used to analyze maps work very well for primary maps, and the same should occur in our applications to the analysis of lens deformations.

Correlations of lens deformations have been calculated in various cases using distinct PM simulations; in case 1, the resolution is $`0.5Mpc`$ (the best one), the box size is $`128Mpc`$, directions are close to the preferred one, and the angular scale is $`\alpha =6.75^{}`$. Correlations obtained in this case are presented in the first row of table 1, and also in the continuous lines of the left panels of Fig. 6. All the data of table 1 are mean correlations with $`2\sigma `$ errors calculated from the correlations of forty simulations. These errors are statistical ones assigned in case 1, but there are other sources of uncertainty as periodicity and simulation features, which are not included in the errors reported in table 1. The lines mentioned above display the variation of the estimated correlations as the number N of simulations (used to perform averages) increases; for $`N>20`$, correlations $`C_2`$ and $`C_4`$ seem to coverge towards a certain value, whereas the correlation $`C_3`$ slowly decreases in the $`N`$-interval (20,40) (see continuous lines in the middle panels of Fig. 6, which are two different representations of the same function). We assign to $`C_2`$, $`C_3`$, and $`C_4`$ the values and errors given in the first row of table 1; the tabulated values of $`C_2`$ and $`C_4`$ are estimations of these correlations, whereas that of $`C_3`$ can be better considered as an upper limit due to the slow decreasing we have mentioned above. The values assigned to $`C_2`$, $`C_3`$, and $`C_4`$ have been used to calculate the ratios $`r_20.081`$, $`r_31.17\times 10^4`$, and $`r_40.014`$, which are smaller than the corresponding ratios of the primary maps (see above). The small $`r_3`$ value indicates that lensing (by nonlinear structures) does not introduce any appreciable $`C_3`$ correlation, whereas the small values of $`C_2`$ and $`C_4`$ suggest small levels of these correlations. Finally, the ratio $`3C_2^2/C_41.43`$ confirms a deviation with respect to Gaussianity.

Since the method proposed in previous sections is based on N-body simulations, results could depend on both resolution and box size. These dependences are now studied. In order to analyze the importance of resolution, cases 2 and 3 have been considered, they are identical to the case 1 described above except for resolution, which is $`0.5Mpc`$ in case 1, $`1.Mpc`$ in case 2, and $`2.Mpc`$ in case 3. Results of case 2 (3) are displayed in the second (third) row of table 1, and also in the dashed (dotted) lines of the left panels of Fig. 6, where we see that the correlations $`C_2`$, $`C_3`$, and $`C_4`$ of case 3 ($`2.Mpc`$ resolution) are smaller than those of case 2 ($`1.Mpc`$) and case 1 ($`0.5Mpc`$), which are very similar between them. In order to quantify these comparisons between results corresponding to pairs of cases, the relative variations $`\mathrm{\Delta }C_m/C_m`$ for $`m=2,3,4`$ are obtained for the pairs 1-2 and 2-3, and the comparation percentages $`100\mathrm{\Delta }C_m/C_m`$ for these pairs are given in the first and second rows of table 2, where we see that the comparation percentages of the pair 1-2 are much smaller than those of 2-3. These considerations strongly suggests that a $`2.Mpc`$ resolution is poor, whereas a $`1.Mpc`$ resolution is good because a better one ($`0.5Mpc`$) does not lead to significantly better results; particularly, for the correlations $`C_2`$ ($`100\mathrm{\Delta }C_2/C_2=0.56`$ for cases 1-2) and $`C_4`$ ($`100\mathrm{\Delta }C_4/C_4=0.37`$ for cases 1-2). Taking into account that cluster description can be improved using higher resolutions, these results have two possible interpretations: (1) the calculation of the lens correlations does not require a more detailed description of the clusters (for the chosen angular scale $`\alpha =6.75^{}`$) and, (2) results from resolutions of $`1Mpc`$ and $`0.5Mpc`$ are similar because of a lack of resolution in both cases, but higher resolutions would lead to –different– better results. We have not found any theoretical reason supporting one of these alternatives. In this situation, the most effective procedure is the use of new simulations with higher resolutions, which should support one of the above possibilities (see Sec. VI for more discussion).

After concluding that a resolution of $`1.Mpc`$ seems to be good (at least it is equivalent to a better resolution of $`0.5Mpc`$), and before analyzing the importance of the box size, a study of the relevance of periodicity is worthwhile. In order to perform such a study, results of case 2 are compared with those of case 4. In both cases, the same forty simulations are used, but the propagation directions of the CMB photons are different. In cases 2 and 4, directions form small angles with the preferred direction and with an edge of the simulation box, respectively. Results of case 4 are displayed in the forth row of table 1, and also in the solid line of the right panels of Fig. 6 (where dashed lines correspond to case 2). Comparing the solid and dashed lines of these panels, one easily concludes that periodicity strongly (weakly) affects the correlations $`C_2`$ and $`C_4`$ ($`C_3`$). Accordingly, the comparation percentages of the correlations $`C_2`$ and $`C_4`$ for the pair of cases 2-4 (third row of table 2) are much greater than that of the $`C_3`$ correlation, which appears to be as small as $`2.09\%`$. This dependence of $`C_2`$ and $`C_4`$ on the direction constraints us to work with directions close to the preferred one, which strongly minimizes periodicity effects (see Sec. IV).

Finally, results of case 2 are compared with those of case 5 to analyze the importance of the box size. Case 5 is identical to case 2 except for this size, which takes on the values $`256Mpc`$ and $`128Mpc`$ in cases 5 and 2, respectively. Possible differences between cases 2 and 5 can be explained in two different ways: (i) boxes of $`256Mpc`$ contain linear inhomogeneities more extended than those involved in boxes of $`128Mpc`$ and, according to our discussion of Sec. III, linear inhomogeneities with different sizes produce distinct lensing and, (ii) there is a residual effect of periodicity, which would depend on box size for the fixed preferred direction used in our calculations (initially associated to $`128Mpc`$ boxes, see Sec. IV). Results of case 5 are shown in the fifth row of table 1, and also in the solid lines of Fig. 7 (where dashed lines correspond to case 2). Comparing the solid and dashed lines of Fig. 7, and taking into account that the comparation percentages of the pair 2-5 (forth row of table 2) are $`13.73\%`$ for $`C_2`$, $`46\%`$ for $`C_3`$, and $`23\%`$ for $`C_4`$, we conclude that box size is important. Is that due to residual periodicity? The answer to this question is negative because the comparation percentage in $`C_3`$ due to periodicity has been proved to be very small ($`2.09\%`$), whereas the corresponding percentage has appeared to be $`46\%`$ when the box size has been varied. Hence, the box size is important as a result of the presence of large scale inhomogeneities and, consequently, our best estimate of the correlations induced by nonlinear lensing corresponds to the smallest box size ($`128Mpc`$) and the best resolution ($`0.5Mpc`$) we have used (case 1). The comparation percentages of the pair 2-5 are not to be confused with the error percentages in the correlations of case 2, which would be expected to be smaller, perhaps about a few per cent. The same has occurred for resolution, where the comparation percentages of the pair 2-3 (second row of table 2) are much larger than the error percentages of case 2 (comparable to the quantities in the first row of table 2).

## VI Discussion

A ray-tracing method to calculate correlations in lens deformations of CMB maps has been designed, tested, and applied. It is based on the use of PM N-body simulations to evolve nonlinear cosmological structures. Although correlations are calculated from the deformations associated to a reduced number of directions (no complete maps of deformations are created), numerical calculations are time consuming because a large enough number of N-body simulations must be performed to calculate good statistical averages (correlations). The method is similar to that proposed by Aliaga et al. ali02 (4) to study the Rees-Sciama effect, and it can be easily implemented starting from an N-body code.

Firts of all, the main structures producing lens deformations in the CMB sky have been identified. A qualitative estimation of the lens deformations produced by a representative set of linear and nonlinear cosmological structures has been used to achieve this identification. Various reasons have motivated the choice of the angular scale $`\alpha =6.75^{}`$ to compute correlations (see Sec. III). For this scale, it has been proved that lens deformations are mainly produced by big and standard clusters, although there is a moderate contribution from linear inhomogeneities with diameters greater than $`100Mpc`$. This information has been crucial to choose the N-body simulations in such a way that the effect of the linear inhomogeneities contained in the simulation box is much smaller than that of the nonlinear structures. Using these simulations, the deviation field $`\stackrel{}{\delta }`$ produced by nonlinear structures can be numerically estimated and, independently, the field $`\stackrel{}{\delta }`$ due to large scale linear inhomogeneities can be analytically calculated. Qualitative arguments suggest that box sizes between $`128Mpc`$ and $`256Mpc`$ could be appropriate to achieve the above requirements (Sec. III) and, then, experiments with simulations (see Sec. V) have been used to reach the following conclusions: (i) a size of $`128Mpc`$ is preferable against larger sizes including extended linear structures (which would produce a significant lensing) and, (ii) resolutions of $`0.5Mpc`$ and $`1.Mpc`$ lead to similar correlations; hence, results converge as resolution increases. Is this convergence real or apparent? As it has been pointed out in Sec. V, the resolution of our simulations must be increased to answer this question. How much should we increase the resolution? We should consider –at least– a resolution of $`0.25Mpc`$ and, perhaps another of $`0.125Mpc`$. These resolutions lead to severe computational problems because many simulations are necessary for statistical purposes. Higher-resolution codes –for example tree or $`P^3M`$ codes– could be more appropriate than our PM one. Fortunately, our simulations lead to boxes containing nonlinear structures which allow us to apply and test our ray-tracing procedure and, furthermore, our best resolution ($`0.5Mpc`$) should lead to a rather acceptable description of clusters, because the central region, where the density is more peaked (distances $`1Mpc`$ from the center), is covered by about $`2^3`$ pairs of cells. We can state that our main goal (description of a new ray-tracing procedure and estimations of nongaussianity from nonlinear structures) has been reached, whereas results from higher resolutions are out of the scope of this work, and they will be developed in due time and presented elsewhere.

In our model, photons move in a periodic universe which is the repetition of the simulation box. We have verified that periodicity is important, but photons do not feel periodicity when they move along directions close to a preferred one (see Sec. IV); this fact constraints us to calculate deformations for this type of directions. Fortunately, these directions suffice to find the required correlations. It is worthwhile to emphasize some novel aspects of our calculations: the basic equations are (17) and (18); the r.h.s of Eq. (17) has the form of a Fourier transform, and it gives –at the same time– the deviation field $`\stackrel{}{\delta }(\stackrel{}{x}__P,\stackrel{}{n})`$ for all the observers located in the nodes ( $`\stackrel{}{x}__P`$ ) of the simulation grid. This fact has been probed to be relevant when averages are performed. The r.h.s. of Eq. (18) is a time integral which can be analytically calculated after using the interpolation formula (19), this procedure simplifies the calculations. The use of directions close to the preferred one is the last novel aspect deserving attention. For these directions and fixed boxes, photons enter the n-th box across a certain region and the (n+1)-th box through a different independent zone; hence, roto-traslations are not necessary to avoid periodicity effects; furthermore, it is worthwhile to emphasize that –without roto-traslations– function $`\stackrel{}{}_{}\varphi `$ is continuous at the points where photons cross boundary boxes (periodicity). However, if the boxes are moved (rototralations), discontinuities at these points are unavoidable and, consequently, the function to be integrated in Eq. (5) has artificial finite discontinuities which influence deviation calculations. Perhaps these discontinuities do not produce any relevant total effect when many boxes are crossed, but potential problems with discontinuities are clearly surmounted by our method in the absence of roto-traslations, which seem not to be either necessary or appropriate for us.

Since periodicity effects are prevented along any direction forming an angle of a few degrees with the preferred one, our method can be used: (1) to simulate squared maps with sizes of a few degrees, for example, $`4^{}\times 4^{}`$ maps with the preferred direction pointing towards their centers, and (2) to compute angular correlations $`C_n`$ for $`n>4`$. In previous sections, correlations $`C_3`$ and $`C_4`$ have been computed using sets of three and four directions close to the preferred one and, evidently, using the same methods and other sets of these directions, the correlation $`C_6`$ and higher ones can be computed. On account of these considerations, future applications of the proposed method are certainly promising.

Hereafter, we only make reference to the lens effect (or lens deformations) produced by nonlinear structures. As it was predicted by Bernardeau ber97 (9), the lens effect does not create any significant $`C_3`$ correlation ($`r_31.17\times 10^4`$). Correlations $`C_2`$ and $`C_4`$ appear to be weaker than those of the primary maps (low correlation level, see Sec. V). The resulting ratio $`3C_2^2/C_41.43`$ points out a deviation from Gaussianity. In the observational maps, the lens deformations are superimposed to both the dominant primary anisotropies and other contributions described in the introduction (foregrounds and nonlinear effects). Let us analyze the superposition of primary maps and lens deformations taking into account that: (1) the correlation $`C_2`$ ($`C_4`$) of primary maps has appeared to be greater than that of lens deformations by a factor close to $`1700`$ ($`4\times 10^6`$), (2) the correlation $`C_3`$ vanishes in primary maps and it has been found to be negligible for deformations, and (3) the two components under consideration are statistically independent. Simple calculations give the correlations of the superposition of components, which are almost identical to those of the primary maps. Deviations with respect to Gaussianity are almost negligible in this superposition. In the observational maps, there are other deviations due to the presence of other components, which would contribute to hide the small deviations associated to lens deformations.

In future, our calculations could be repeated using other N-body codes and resolutions, the lens effect due to linear structures contained in boxes with different sizes could be studied in more detail; the present study could be extended to other components of the observational CMB maps, in particular to the Rees Sciama component ali02 (4); maps of lens deformations and other components could be simulated –with our method– and these maps (which would have the true statistical properties at any order of correlation) could be analyzed using various estimators of nonGaussianity proposed in the literature tak01 (20, 21); our simulations could be also used to study the correlation between galaxy shear and CMB temperature distortions van00 (22) and, finally, using other N-body codes and greater resolutions, it would be worthwhile the study of very small (large) angular scales ($`\mathrm{}`$ values), for which, the $`C_{\mathrm{}}`$ coefficients of the primary anisotropy are very small. Some of these studies and other possible ones could be useful to improve on the method proposed in this paper; nevertheless, a good version has been already implemented and numerical codes are operative.

###### Acknowledgements.

This work has been partially supported by the Spanish MCyT (project AYA2000-2045). VQ is Ramón y Cajal Fellow of the Spanish Ministry of Science and Technology, and PC thanks the same Ministry for a fellowship. Calculations were carried out on a SGI Origin 2000s at the Centro de Informática de la Universitat de València.


format: text

text:
# SEARCHING FOR NEW CONDITIONS FOR FERMION 𝑁-REPRESENTABILITY

## 1. Introduction

The $`N`$-representability problem  appears in non-relativistic quantum mechanics of $`N`$-fermion systems. The Hamiltonian $`H^N`$ for systems of $`N`$-fermions contain only operators with 1- and 2-body interactions:

$$H^N=\underset{i=1}{\overset{N}{}}H^1(i)+\underset{1i<jN}{}H^2(i,j)=\underset{1i<jN}{}h^2(i,j).$$

The ground state energy of the system can be determined variationally by minimizing the $`N`$-particle functional

$$E=\begin{array}{c}inf\\ D^N𝒫^N\end{array}Tr(H^ND^N)$$

over the set $`𝒫^N`$ consisting of all fermion $`N`$-particle density operators. Instead, because of the appearance of at most 2-body interactions between the particles, the ground state energy of the system of $`N`$-fermions could be, in principle, determined variationally by minimizing the 2-particle functional

$$E=\begin{array}{c}inf\\ D^2𝒫_N^2\end{array}\left(\genfrac{}{}{0pt}{}{N}{2}\right)Tr(h^2D^2)$$

over the set $`𝒫_N^2`$ consisting of all fermion 2-particle reduced density operators, i.e. such 2-particle density operators $`D^2`$ which possess an $`N`$-particle fermion preimage $`D^N`$:

$$D^2=L_N^2D^N=Tr_{3,\mathrm{},N}D^N.$$

It is known  that the set $`𝒫_N^2`$ is a proper convex subset of $`𝒫^2`$, the set of all fermion 2-density operators. However, the complete characterization of $`𝒫_N^2`$ as a convex proper subset of $`𝒫^2`$ is not yet known. It has been shown  that a knowledge of all exposed points (which are extreme) of $`𝒫_N^2`$ is sufficient to characterize the closure of $`𝒫_N^2`$; only some of them are known. The dual characterization of $`𝒫_N^2`$ involves a determination of the dual (polar) cone $`\stackrel{~}{𝒫}_N^2`$ consisting of 2-particle self-adjoint operators $`b^2`$ for which $`Tr(b^2D^2)0`$, $`D^2𝒫_N^2`$, that is equivalent to the positive-semidefiniteness of the $`N`$-particle operator $`\mathrm{\Gamma }_2^Nb^2=b^2I^{(N2)}=A^Nb^2I^{(N2)}A^N0`$. The dual cone $`\stackrel{~}{𝒫}_N^2`$ is a convex one. Any element of the dual cone $`\stackrel{~}{𝒫}_N^2`$ provides an $`N`$-representability condition. Those coming from the extreme elements of $`\stackrel{~}{𝒫}_N^2`$ are the strongest ones. They give the hyperplane characterization of $`𝒫_N^2`$ and thus the solution of the $`N`$-representability problem. Several necessary conditions for $`N`$-representability have been derived and some of their structural features and mutual interrelations are established .

In this paper we find new elements of the dual cone $`\stackrel{~}{𝒫}_N^2`$, and thus define new necessary conditions for $`N`$-representability of a trial 2-particle density matrix $`D^2`$. In the case N=3, in which the new elements of the dual cone are given explicitly, they lead to necessary conditions for 3-representability of the trial $`D^2`$ which are stronger in comparison with the known $`B`$ and $`C`$ conditions . The results show that in the spectral decomposition of a 2-particle fermion density operator $`D^2`$ there exists an intimate relation between the eigenvalue and the corresponding eigenfunction which has to be satisfied in order that $`D^2`$ can be 3-representable ( in general $`N`$-representable). More precisely, the condition obtained for 3-representability says that the upper bound of the eigenvalue of $`D^2`$ is a functional of the corresponding eigenfunction. It is worthwhile to remember that for $`N`$-representability of a 1-particle density matrix $`D^1`$ the upper bound on the eigenvalues does not depend on the eigenfunctions. This explains why the $`N`$-representability problem for $`D^2`$ is so much harder than for $`D^1`$.

## 2. The dual P-condition.

In this paper the underlying 1-particle Hilbert space $`^1`$ is finite dimensional $`dim^1=n`$. $`^2`$ is a 2-particle antisymmetric Hilbert space, the Grassmann product $`^1^1=A^2^1^1`$. $`P_g^2=g^2\overline{g}^2`$ is the projection operator onto a 2-particle antisymmetric function $`g^2^2`$. The operator $`\left(\genfrac{}{}{0pt}{}{N}{2}\right)P_g^2I^{(N2)}=\left(\genfrac{}{}{0pt}{}{N}{2}\right)A^NP_g^2I^{(N2)}A^N=A^N_{1i<jN}P_g^2(i,j)I^{N2}(1,\mathrm{},i1,i+1,\mathrm{},j1,j+1,\mathrm{},N)A^N`$ is the simplest (elementary) antisymmetric operator with ”2-body interactions” acting on $`^N`$. Because $`\left(\genfrac{}{}{0pt}{}{N}{2}\right)P_g^2I^{(N2)}`$ is positive-semidefinite, $`P_g^2`$ belongs to the dual cone $`\stackrel{~}{𝒫}_N^2`$, and therefore $`Tr(D^2P_g^2)0`$ is a necessary condition for $`N`$-representability of $`D^2`$ (the P-condition). Besides, there exists another element of $`\stackrel{~}{𝒫}_N^2`$ that is generated by the 2-particle antisymmetric function $`g^2^2`$. Let $`\mathrm{\Lambda }_{max}(g^2)`$ denote the maximal eigenvalue of $`\left(\genfrac{}{}{0pt}{}{N}{2}\right)P_g^2I^{(N2)}`$, which in general depends on $`g^2`$, then the operator $`\mathrm{\Lambda }_{max}(g^2)I^N\left(\genfrac{}{}{0pt}{}{N}{2}\right)P_g^2I^{(N2)}0`$ (is positive-semidefinite), and therefore the 2-particle operator $`\frac{\mathrm{\Lambda }_{max}(g^2)}{\left(\genfrac{}{}{0pt}{}{N}{2}\right)}I^2P_g^2`$ belongs to the dual cone $`\stackrel{~}{𝒫}_N^2`$ . Thus, we have obtained

###### Theorem 2.1:

For any $`g^2^2`$ the operator

$$\frac{\mathrm{\Lambda }_{max}(g^2)}{\left(\genfrac{}{}{0pt}{}{N}{2}\right)}I^2P_g^2\stackrel{~}{𝒫}_N^2,$$
(2.1)
and gives the following necessary condition for $`N`$\- representability of a 2-particle fermion density matrix $`D^2`$: $`Tr[(\frac{\mathrm{\Lambda }_{max}(g^2)}{\left(\genfrac{}{}{0pt}{}{N}{2}\right)}I^2P_g^2)D^2]0`$, i.e.

$`Tr(D^2P_g^2){\displaystyle \frac{\mathrm{\Lambda }_{max}(g^2)}{\left(\genfrac{}{}{0pt}{}{N}{2}\right)}},`$ $`g^2^2.`$ (2.2)
Here, $`\mathrm{\Lambda }_{max}(g^2)`$ is the maximal eigenvalue of the operator $`\left(\genfrac{}{}{0pt}{}{N}{2}\right)P_g^2I^{(N2)}`$.

We propose to call the new condition ”the dual P-condition”.

In particular, from the above theorem follows

###### Theorem 2.2:

If $`D^2=_{i=1}^{\left(\genfrac{}{}{0pt}{}{n}{2}\right)}\lambda _iP_{g_i}^2`$ is the spectral decomposition of $`D^2`$, then it is a necessary condition for $`N`$-representability that the eigenvalues $`\lambda _i`$ must satisfy the inequalities $`\lambda _i\frac{\mathrm{\Lambda }_{max}(g_i^2)}{\left(\genfrac{}{}{0pt}{}{N}{2}\right)}`$. Here, $`\mathrm{\Lambda }_{max}(g_i^2)`$ is the maximal eigenvalue of the operator $`\left(\genfrac{}{}{0pt}{}{N}{2}\right)P_{g_i}^2I^{(N2)}`$, where $`g_i^2`$ is the eigenfunction corresponding to the eigenvalue $`\lambda _i`$.

This theorem shows that the bound on the eigenvalue is a functional of the corresponding eigenfunction. For the $`N`$-fermion 1-body elementary operator $`NP_g^1I^{(N1)}`$, $`g^1^1`$, the maximal eigenvalue $`\mathrm{\Lambda }_{max}(g^1)=1`$ for any $`g^1`$, because $`NP_g^1I^{(N1)}`$ is a projection operator. Hence, the equivalent of Theorem 2.1 says that $`Tr(D^1P_g^1)\frac{1}{N}`$ for arbitrary $`g^1^1`$, which in particular means that the bound on an eigenvalue of an $`N`$-representable 1-density operator does not depend on the corresponding eigenfunction, contrary to the 2-density matrix case.

As seen from Theorem 2.1, the new necessary condition for $`N`$-representability of a 2-particle density matrix $`D^2`$ requires knowledge of the maximal eigenvalue $`\mathrm{\Lambda }_{max}(g^2)`$ of the $`N`$\- particle operator $`\left(\genfrac{}{}{0pt}{}{N}{2}\right)P_g^2I^{(N2)}`$. It is rather hopeless to find $`\mathrm{\Lambda }_{max}(g^2)`$ for arbitrary $`N`$ and any $`g^2^2`$. But even the solution for some ”simple” $`g^2`$ would contribute to knowledge of the structure of both the dual cone $`\stackrel{~}{𝒫}_N^2`$ and the set $`𝒫_N^2`$ of fermion $`N`$-representable 2-density operators. So far, we have results for $`N=3`$, in which case it was possible to find the spectral decomposition of the elementary operator $`3P_g^2I^1`$ for arbitrary $`g^2^2`$. The details and proofs concerning this spectral decomposition and the reduced 2- and 1-particle density operators corresponding to the eigenstates of the operator $`3P_g^2I^1`$ will be published in a separate paper . In this paper we use only the following result which we formulate as

###### Theorem 2.3:

Let $`^1`$ be a finite dimensional Hilbert space (dim $`^1=n`$), and $`^2=^1^1`$ denotes the 2-particle antisymmetric space generated by $`^1`$ (the 2-fold Grassmann product of $`^1`$). Let $`P_g^2`$ denote the 1-dimensional projection operator onto a 2-particle antisymmetric function $`g^2^2`$ of 1-rank $`r=2s`$ possessing the canonical decomposition $`g^2=_{i=1}^s\xi _i2i1,2i`$ with $`_{i=1}^s\xi _i^2=1`$, where $`2i1,2i=\sqrt{2}\varphi _{2i1}^1\varphi _{2i}^1=\frac{1}{\sqrt{2}}det(\varphi _{2i1}^1,\varphi _{2i}^1)`$ is the 2-particle normalized Slater detrminant. Let the identity operator $`I^1`$ on $`^1`$ possess the decomposition $`I^1=_{i=1}^{r=2s}P_i^1+_{i=r+1}^nP_i^1`$, where $`P_i^1=ii=\varphi _i^1\overline{\varphi }_i^1(i=1,\mathrm{},n)`$ are 1-dimensional mutually orthogonal projection operators onto the functions $`i=\varphi _i^1`$. Then, the 3-particle operator $`3P_g^2I^1`$ possesses the following spectral decomposition

$$3P_g^2I^1=\underset{k=1}{\overset{s=r/2}{}}(1\xi _k^2)(P_{g_{2k1}}^3+P_{g_{2k}}^3)+\underset{l=r+1}{\overset{n}{}}P_{g_l}^3+0Ker(3P_g^2I^1).$$

Here, $`P_{g_{2k1}}^3,P_{g_{2k}}^3(k=1,\mathrm{},s=r/2),P_{g_l}^3(l=r+1,\mathrm{},n)`$ are 1-dim projectors onto the following functions:

$$g_{2k1}^3=\frac{1}{\sqrt{1\xi _k^2}}\underset{\begin{array}{c}i=1\\ (ik)\end{array}}{\overset{s}{}}\xi _i2i1,2i,2k1=\sqrt{\frac{3}{1\xi _k^2}}g^22k1,$$

$$g_{2k}^3=\frac{1}{\sqrt{1\xi _k^2}}\underset{\begin{array}{c}i=1\\ (ik)\end{array}}{\overset{s}{}}\xi _i2i1,2i,2k=\sqrt{\frac{3}{1\xi _k^2}}g^22k,$$

$$g_l^3=\underset{i=1}{\overset{s}{}}\xi _i2i1,2i,l=\sqrt{3}g^2l,$$

$`Ker(3P_g^2I^1)`$ denotes the projection operator onto the null-space of the operator $`3P_g^2I^1`$, which is of dimension $`\left(\genfrac{}{}{0pt}{}{n}{3}\right)n`$. The symbols of the type $`2i1,2i,2k1`$ denote the appropriate 3-particle Slater determinants.

Now, we can apply results of Theorem 2.3 to find new necessary conditions for 3-representability according to the formulae of Theorems 2.1 and 2.2 . The required $`\mathrm{\Lambda }_{max}(g^2)`$ is either $`1\xi _{min}^2`$, if the 1-rank of $`g^2`$ is equal to the dimension $`n`$ of the 1-particle Hilbert space $`^1`$ ($`r=2s=n`$), or $`\mathrm{\Lambda }_{max}(g^2)=1`$, if $`r`$ is less than $`n`$ ($`r<n`$). Hence, we have

###### Theorem 2.4:

If $`D^2`$ ($`TrD^2=1`$) is a 2-fermion density matrix, then for 3-representability it must satisfy for any $`g^2^2`$, with the canonical decomposition $`g^2=_{i=1}^{r/2}\xi _i2i1,2i`$, the following inequality:

$$Tr(D^2P_g^2)\frac{1}{3}(1\xi _{min}^2),$$
(2.3)
if the 1-rank of $`g^2`$ equals $`n=dim^1(r=n),`$

$$Tr(D^2P_g^2)\frac{1}{3},$$
(2.4)
if the 1-rank $`r`$ of $`g^2`$ is less than $`dim^1(r<n).`$

In particular, if we take as $`g^2`$ the eigenfunctions $`g_i^2`$ of $`D^2`$, we get

###### Theorem 2.5:

If $`D^2`$ has the spectral decomposition $`D^2=_{i=1}^{\left(\genfrac{}{}{0pt}{}{n}{2}\right)}\lambda _iP_{g_i}^2`$ with $`g_i^2=_{j=1}^{s_i}\xi _{ij}2j1,2j`$, then $`\lambda _i\frac{1}{3}(1\xi _{ij}_{\begin{array}{c}min\\ j\end{array}}^2)`$, if the eigenfunction $`g_i^2`$ has 1-rank $`r_i=2s_i=n`$, and $`\lambda _i\frac{1}{3}`$, if the eigenfunction belonging to $`\lambda _i`$ has 1-rank $`r_i<n=dim^1`$.

It is worthwhile to observe that the eigenfunctions belonging to the maximal eigenvalue $`\lambda =\frac{1}{3}`$ cannot have full 1-rank ( $`r=n=dim^1`$) if $`D^2`$ is 3-representable, and on the other hand, the eigenvalues corresponding to the eigenfunctions with full 1-rank must be strictly less than $`\frac{1}{3}`$.

## 3. Strength of the dual P-condition for fermion 3-representability

It is important to compare the new condition with the already known conditions for 3-representability because the effort to find $`\mathrm{\Lambda }_{max}(g^2)`$ for at least some $`g^2`$ in the general case (arbitrary $`N`$) might not pay off. Since the new condition is an estimation from the above on the expectation value of $`D^2`$ in any state $`g^2^2`$, it can be compared with the $`B`$\- and $`C`$\- conditions for $`N`$-representability of $`D^2`$ :

$`Tr[B_N^2(g^2)D^2]0,`$ $`Tr[C_N^2(g^2)D^2]0,`$ $`g^2^2,`$ (3.1)
where

$$B_N^2(g^2)=I^2(N2)L_2^1P_g^2I^1(N1)P_g^2\stackrel{~}{𝒫}_N^2,$$

$$C_N^2(g^2)=(nN+2)L_2^1P_g^2I^1(N1)P_g^2\stackrel{~}{𝒫}_N^2.$$

In the case $`N=3`$, conditions (3.1) take the form:

$`Tr(D^2P_g^2){\displaystyle \frac{1}{2}}[1Tr(L_2^1D^2L_2^1P_g^2)],`$ $`g^2^2`$ (3.2)
$`Tr(D^2P_g^2){\displaystyle \frac{n1}{2}}Tr(L_2^1D^2L_2^1P_g^2),`$ $`g^2^2.`$ (3.3)
In order to compare the inequalities (2.2), (3.2), (3.3), we choose as $`g^2`$ the so called ”extreme geminal”  $`g_{extr}^2`$ of 1-rank $`r=2s=n=dim^1`$ possessing the following canonical decomposition:

$$g_{extr}^2=\underset{i=1}{\overset{n/2}{}}\sqrt{\frac{2}{n}}2i1,2i.$$

The required by inequality (2.3) $`\xi _{min}^2=\frac{2}{n}`$, and therefore the new condition gives

$$Tr(D^2P_{g_{extr}}^2)\frac{1}{3}(1\frac{2}{n}).$$
(3.4)
On the other hand, since $`L_2^1P_{g_{extr}}^2=_{i=1}^n\frac{1}{n}P_i^1`$ , $`P_i^1=ii`$, $`_{i=1}^nP_i^1=I^1`$, we have from both (3.2) and (3.3) the same result

$$Tr(D^2P_{g_{extr}}^2)\frac{1}{2}(1\frac{1}{n}).$$
(3.5)
Comparing (3.4) with (3.5), we see that the new condition is stronger than the B and C ones. So far, we were unable to establish the relation between the new condition (2.3) and the G-condition  for 3-representability. This would be important because it is known  that the G-condition implies the B- and C- ones.

## 4. The strengthened B-condition

While comparing the new condition with the B-condition we found out that the B-condition could in principle be improved for N odd. $`B_N^2(g^2)\stackrel{~}{𝒫}_N^2`$, $`g^2^2`$ means that the operator $`\left(\genfrac{}{}{0pt}{}{N}{2}\right)B_N^2(g^2)I^{(N2)}0`$, $`g^2^2`$. If $`\mathrm{\Lambda }_{min}^{B(g)}`$ is the minimal eigenvalue of the operator $`\left(\genfrac{}{}{0pt}{}{N}{2}\right)B_N^2(g^2)I^{(N2)}`$, then the operator $`\left(\genfrac{}{}{0pt}{}{N}{2}\right)B_N^2(g^2)I^{(N2)}\mathrm{\Lambda }_{min}^{B(g)}I^N`$ is positive-semidefinite and therefore

$$\left(\genfrac{}{}{0pt}{}{N}{2}\right)B_N^2(g^2)\mathrm{\Lambda }_{min}^{B(g)}I^2\stackrel{~}{𝒫}_N^2,$$
(4.1)
giving a necessary condition for $`N`$-representability (the strengthened B-condition). It is known  that for $`N`$ even $`\mathrm{\Lambda }_{min}^{B(g)}=0`$, $`g^2^2`$. For $`N=3`$ and $`g_{extr}^2`$, the eigenvalue $`\mathrm{\Lambda }_{min}^{B(g_{extr})}`$ is greater than zero, and can be found explicitly. The operator

$$3B_3^2(g_{extr}^2)I^1=3(1\frac{1}{n})I^32(3P_{g_{extr}}^2I^1)$$

possesses the following minimal eigenvalue

$$\mathrm{\Lambda }_{min}^{B(g_{extr})}=3(1\frac{1}{n})2(1\frac{2}{n})=1+\frac{1}{n},$$

where Theorem 2.3 has been used. Hence, the strengthened B-condition in this case is

$$3B_3^2(g_{extr}^2)\mathrm{\Lambda }_{min}^{B(g_{extr})}I^2=\frac{2(n2)}{n}I^22(3P_{g_{extr}}^2)\stackrel{~}{𝒫}_N^2,$$
(4.2)
and therefore

$$Tr(D^2P_{g_{extr}}^2)\frac{1}{3}(1\frac{2}{n}).$$
(4.3)
Comparing inequalities (3.4) and (4.3) we see that both the new conditions for 3-representability in the case under consideration ($`g^2=g_{extr}^2`$) give the same bound from above on the expectation value $`Tr(D^2P_{g_{extr}}^2)`$.

It seems that the above results suggest, it would be worthwhile to make an effort to extend the new conditions to arbitrary $`N`$ for at least some $`g^2`$ for which the maximum eigenvalue $`\mathrm{\Lambda }_{max}(g^2)`$ of the operator $`\left(\genfrac{}{}{0pt}{}{N}{2}\right)P_g^2I^{(N2)}`$ can be found, e.g. the extreme geminal $`g_{extr}^2`$. So far, we have obtained only partial information about the spectral decomposition of the operator $`\left(\genfrac{}{}{0pt}{}{N}{2}\right)P_{g_{extr}}^2I^{(N2)}`$ for arbitrary $`N`$. However, at least in this case it is realistic to succeed. The other known conditions for $`N`$-representability could be treated in a similar way as the P- and B- conditions considered in this paper, provided the maximal and minimal eigenvalues of the appropriate operators were known. The obtained results concerning arbitrary $`N`$ will be published later.

Acknowledgements

This work was supported in part by the Polish Committee of Scientific Research (KBN) under grant No. 5 P03B 083 20.


text:
# Untitled Document

Journal of Statistical Physics, Vol. 96, Nos. 5/6, 1999

A Remark on the Kramers Problem

Alex A. Samoletov111Laboratoire de Physique, Ecole Normale Supérieure de Lyon, 69364 Lyon cedex 07, France; Permanent address: Institute for Physics and Technology, National Academy of Sciences of the Ukraine, 83114 Donetsk, Ukraine. E-mail: samolet@kinetic.ac.donetsk.ua

We present new point of view on the old problem, the Kramers problem. The passage from the Fokker-Planck equation to the Smoluchowski equation, including corrections to the Smoluchowski current, is treated through an asymptotic expansion of the solution of the stochastic dynamical equations. The case of an extremely weak force of friction is also discussed.   KEY WORDS: Kramers problem; Klyatskin-Novikov theory; stochastic equation; current.

1. Evolution of a physical system can be ordered in multi-time-scales. Details of evolution on short-time-scale do not need for description in a closed form of a system evolution on long-time-scale and appears on this scale only in an average form. The prototype of such kind physical systems is dissipative Brownian motion of a particle in an external potential field. In this problem, with the exception of extremely short characteristic time scales of random forces, there are two time scales: (1) time scales of a particle motion in an external field; (2) time scales of relaxation (rate of dissipation) of Brownian particle in a media. It is intuitively absolutely clear that, if the friction force is strong enough (time of free motion is extremely short), then probability distribution of a particle velocity to be rapidly relaxed to the Maxwell distribution and on this background a particle position will be undergoing to slow process of diffusion. In the following we deal with the consideration of approximate reduction of the Fokker-Planck equation for phase-space probability density to the Smoluchowski equation which deals with probability density of a particle position only. In the opposite case of extremely weak force of friction we have energy (or action) variable as evident slow one.

2. The Kramers problem consist in mathematical description of approximate reductions of the Fokker-Planck equation for dissipative Brownian motion of a particle in an external field to the two limiting cases: (1) to the Smoluchowski equation (extremely strong force of friction) or (2) to equation for probability density of energy (or action) variable (extremely weak force of friction)(1). These reduction procedures are the prototypes of all adiabatic elimination procedures or procedures of separation on slow and fast variables(2).

The Kramers model(1), firstly formulated for kinetics of chemical reactions, consists of a particle of mass $`m`$ moving in an one-dimensional potential field $`U(x)`$ under influence of a random force $`f(t)`$ and a linear friction force with a constant dissipation rate $`\lambda `$. The corresponding set of Langevin equations has the following form

$$\dot{x}=u,m\dot{u}=U^{}(x)\lambda mu+f(t),$$
(1)
where the random force $`f(t)`$ is generalized Gaussian $`\delta `$-correlated stochastic processes (white noise) with the following properties (including the fluctuation-dissipation relation)

$$f(t)=0,f(t)f(t^{})=2\lambda mk_BT\delta (tt^{});$$
(2)
$`\mathrm{}`$ denotes average over all realizations of random force.

The Langevin (1)-(2) dynamics is stochastically equivalent to the Fokker-Planck equation for the rate of change of probability density $`P(u,x;t)`$ which has the form (e.g.(2))

$$_tP(u,x;t)=u_xP+\frac{1}{m}U^{}(x)_uP+\frac{\lambda }{m}_u[uP+\frac{k_BT}{m}_uP],$$
(3)
where $`_t=/t,`$ $`_x=/x,`$ and $`_u=/u.`$

Keeping in mind (1)-(3) we formulate the problem in the following manner: in the case of extremely strong force of friction beginning with (3) or equivalently (1)-(2) to derive the approximate reduction to an equation for the rate of change of probability density $`P(x;t)`$ in the form of asymptotic expansion by the parameter $`\lambda ^1`$:

$$_tP(x;t)=_x[\lambda ^1J_S+o(\lambda ^1)],$$
(4)
where $`J_S`$ is the Smoluchowski current,

$$J_S=[U^{}(x)P(x;t)+k_BT_xP(x;t)]/m$$
(5)
In other words it is asymptotics of strong force of friction on time scales $`\lambda t1`$.

This problem has long history starting since 1940, the date of publication of the Kramers famous work(1). For relevant references including reviews of the problem see refs. 2–5. The first treatment of the problem has been down in ref. 6 and the first correct solution has been down in ref. 7 and then in refs. 8–10. The works(10-14) are of importance for the following in respect of treatment of the corrections older then $`\lambda ^3`$ which break the Fokker-Planck structure of (4). Most general treatment of the problem has been down in ref. 15.

3. All of the cited works deal with (3) as the input equation for a solution of the problem. The purpose of this paper is to take notice of the fact that (1)-(2) are indeed convenient input equations for an answer to the problem. With respect to solutions of (1)-(2) we use the method of asymptotic expansion by the parameter $`\lambda ^1`$. In the way, the Fokker-Planck type equations (4) to be derived from approximate stochastic dynamical equations in each order of $`\lambda ^1`$. Moreover, Fokker-Planck equation is an approximate equation and in any case must be derived from input dynamical equations. Convenient and powerful method of derivation, in particular, of the Fokker-Planck type equations immediately from stochastic equations has been initiated by Novikov(16) and then it has been sufficiently developed by Klyatskin(17). In the following we use this method systematically. In this connection it should be pointed out that the Klyatskin-Novikov theory, generally, interprets a stochastic differential equation in the sense of Stratonovich. However, in the case under consideration it does not matter. We refer the reader to refs. 16 and 17 for further information.

The probability density $`P(x;t)`$ can be written in the form(17)

$$P(x;t)=\delta (xx(t)),$$

where $`x(t)`$ is a stochastic process and $`\delta (\mathrm{})`$ is $`\delta `$-function. Differentiating this definition by time we obtain the equation

$$_tP(x;t)=_x\dot{x}(t)\delta (xx(t))_xJ(x;t),$$
(6)
which has the form of a conservation law and is the proforma for an equation of the type (4). If $`x(t)`$ is defined by (1), our problem is in calculation of asymptotic expansion of $`\dot{x}(t)`$ by $`\lambda ^1`$ and then the corresponding average in (6). Further insight is gained by making the following construction. Rewrite (1) for $`\lambda t1`$ in the form

$$\dot{x}(t)=\frac{1}{m}\mathrm{\Lambda }^1U^{}(x)+\frac{1}{m}\xi (t),$$
(7)
where operator $`\mathrm{\Lambda }`$ has the form $`\mathrm{\Lambda }=d/dt+\lambda ,`$ and the Ornstein-Uhlenbeck process is introduced:

$$\xi (t)=\mathrm{\Lambda }^1f(t)=\mathrm{exp}(\lambda t)_0^t\mathrm{exp}(\lambda t^{})f(t^{})𝑑t^{}.$$

Formal expansion of $`\mathrm{\Lambda }^1`$ by $`\lambda ^1`$ has the form

$$\mathrm{\Lambda }^1U^{}(x)=\frac{1}{\lambda }\underset{n=0}{\overset{N}{}}\frac{(1)^n}{\lambda ^n}\left(\frac{d^n}{dt^n}\right)U^{}(x)+\mathrm{}.$$

Hence, (7) can be written in the form

$`\dot{x}(t)`$ $`\left[{\displaystyle \frac{1}{m\lambda }}U^{}(x)+{\displaystyle \frac{1}{m}}\xi (t)\right]+{\displaystyle \frac{1}{m\lambda ^2}}\left[U^{\prime \prime }(x)\dot{x}(t)\right]`$
$`{\displaystyle \frac{1}{m\lambda ^3}}\left[U^{\prime \prime \prime }(x)(\dot{x}(t))^2+U^{\prime \prime }(x)\ddot{x}(t)\right]+\mathrm{}.`$
The $`\dot{x}(t)`$, $`\ddot{x}(t)`$, and so on, can be excluded from right hand side of last equation repeatedly using iterations of this equation and its time derivatives. Then, and it is important, we are in need of expansion by $`\lambda ^1`$ of the stochastic process $`\xi (t)`$ or, more precisely, of expansion by $`\lambda ^1`$ of the average in (6) which involves $`\xi (t)`$. First of all we must remark that the derivatives $`\dot{f}(t)`$ and so on, have sense only as derivatives of the generalized stochastic process $`f(t)`$(18) and break the simple Fokker-Planck structure of (6) as of an second order partial differential equation. It is evident in the frame of Klyatskin-Novikov theory(17). Namely in the process of calculation of the corresponding averages according to ref. 17 we easy detect a complex form of (6) including memory as well as an integral-operator structure. Further, in respect of $`\lambda ^1f(t)`$ it is easy to verify(17) that corresponding averages in (6) have factor $`\lambda ^1`$ because the noise (2) intensity has the order $`\lambda `$.

Taking into account what has been outlined above the first terms of expansion of $`\dot{x}(t)`$ by $`\lambda ^1`$ that lead to the current $`J(x;t)`$ expansion up to order $`\lambda ^3`$ (the maximum-order of saving of the simple Fokker-Planck structure of (6)) can be written in the form

$$\dot{x}(t)\left(1+\frac{1}{m\lambda ^2}U^{\prime \prime }(x)\right)\left[\frac{1}{m\lambda }U^{}(x)+\frac{1}{m\lambda }f(t)\right]+\mathrm{}.$$

Substituting last expression into the current $`J(x;t)`$ (6) and performing averaging exactly follow Klyatskin-Novikov theory(17) we obtain

$`J(x;t)`$ $`=\left(1+{\displaystyle \frac{1}{m\lambda ^2}}U^{\prime \prime }(x)\right)\left[{\displaystyle \frac{1}{m\lambda }}U^{}(x)+{\displaystyle \frac{1}{m\lambda }}f(t)\right]\delta (xx(t))+o(\lambda ^3)`$
$`=\left(1+{\displaystyle \frac{1}{m\lambda ^2}}U^{\prime \prime }(x)\right)J_S(x;t)+o(\lambda ^3).`$ (8)
where $`J_S`$ is the Smoluchowski current (5). Asymptotic expansions of $`\dot{x}(t)`$ and $`J(x;t)`$(8) together with (2) lead to conventional conclusion: the Smoluchowski equation is valid if: $`\lambda t1`$, $`l|U^{}(x)|k_BT`$, $`l^2|U^{\prime \prime }(x)|k_BT`$, - where a length scale $`l=\sqrt{k_BT/m\lambda ^2}`$ is introduced. (8) contains lowest order correction to the Smoluchowski equation(8-10). Higher order corrections in $`\lambda ^1`$, involving in averaging time derivatives of the generalized stochastic process $`f(t)`$, lead to break of simple structure of the Smoluchowski equation as a second order partial differential equation. It was pointed out also in traditional approach(10-14,5).

4. Consider now the case of extremely weak force of friction. This case more complicated then previous but not so interesting in calculation.

Let $`m=1`$ in (1). In the case of extremely weak force of friction and on the time-scale $`\lambda t1`$ the energy $`E=u^2/2+U(x)`$ of unperturbed system is evident candidate for slow variable. But previously $`E`$ must be averaged over period of relatively rapid dynamical oscillations. It is more convenient, however , to consider the action variable $`J`$ instead of $`E`$, $`J=J(E)`$.(1) Let $`J`$ is action variable averaged over period of rapid dynamical oscillations. Then an equation for the rate of change of probability density $`P(J;t)`$ can be written in the form (see (6))

$$_tP(J;t)=_J\dot{J}(t)\delta (JJ(t)),P(J;t)=\delta (JJ(t)),$$
(9)
In usual way(1) and taking into account the change of time-scale of the white noise(2,7) we obtain the equations of motion for slow variables

$$\dot{J}(t)=\lambda J+\frac{V}{\omega }f(t),\dot{V}=\lambda \frac{V}{\omega }J+f(t),$$
(10)
where $`\omega =\omega (J)=dE/dJ`$ is frequency and $`V`$ is velocity averaged over period of dynamical motion. Substituting (10) in (9) we obtain

$$_tP(J;t)=_J\left[\lambda JP+\frac{1}{\omega (J)}f(t)V(t)\delta (JJ(t))\right].$$
(11)
For calculation of the average in right hand side of (11) we can use the Klyatskin-Novikov procedure again. Using (9) and the causality condition(17) we obtain for functional derivatives

$$\frac{\delta V(t)}{\delta f(t)}=1;\frac{\delta J(t)}{\delta f(t)}=\frac{V}{\omega }.$$

Taking also into account that $`V^2/\omega =J`$ if $`J`$=const, according to ref. 17 we finally obtain

$`_tP(J;t)`$ $`=_J\left[\lambda JP{\displaystyle \frac{\lambda k_BT}{\omega (J)}}{\displaystyle \frac{\delta V(t)}{\delta f(t)}}\delta (JJ(t))V(t)_J\delta (JJ(t)){\displaystyle \frac{\delta J(t)}{\delta f(t)}}\right]`$
$`=_J\left[\lambda JP{\displaystyle \frac{\lambda k_BT}{\omega (J)}}P+{\displaystyle \frac{\lambda k_BT}{\omega (J)}}_J\left(JP\right)\right]`$
$`=_J\left[\lambda J+\lambda k_BT{\displaystyle \frac{J}{\omega (J)}}_J\right]P(J;t).`$
It is exactly the Kramers equation. We can derive corrections to this equation but it is slightly more difficult task then in the case of extremely strong force of friction and does not take special interest in the context of this paper.

5. In conclusion, we have presented in a simplest framework a unique approach to the kinetic equations for slow variables by taking stochastic dynamical equations as the input instead of the Fokker-Planck equation. We hope that this approach is general enough.

ACKNOWLEDGMENTS

I would like to thank Michel Peyrard for his financial support and kind hospitality at the ENS-Lyon. This work was also supported in part by National Foundation for Basic Research (Grant No. F4/310-97),

REFERENCES

1. H. A. Kramers, Physica 7:284 (1940).
2. C. W. Gardiner, Handbook of Stochastic Methods (2nd ed.) (Springer, Berlin, 1997).
3. P. Hänggi, P. Talkner and M. Berkovec, Rev. Mod. Phys. 62:251 (1990).
4. N. G. van Kampen, Stochastic Processes in Physics and Chemistry (North-Holland,
Amsterdam, 1984).
5. N. G. van Kampen, Phys. Rep. 124:69 (1985).
6. H. C. Brinkman, Physica 22:29 (1956).
7. R. L. Stratonovich, Topics in the Theory of Random Noise. Vol.1.(Gordon and Breach,
New York, 1963).
8. G. Wilemski, J. Stat. Phys. 14:153 (1976).
9. U. M. Titulaer, Physica A91:321 (1978).
10. U. M. Titulaer, Physica  A100:251 (1980).
11. H. Risken and H.D. Vollmer, Z. Phys. B33:297 (1979).
12. H. Risken and H.D. Vollmer, Z. Phys. B35:177 (1979).
13. H. Risken, H.D. Vollmer and H. Denk, Phys. Lett A78:22 (1980).
14. H. Risken, H.D. Vollmer and M. Mörsch, Z. Phys. B40:343 (1980).
15. V. I. Mel’nikov and S. V. Meshkov, J. Chem. Phys. 85:1018 (1986).
16. E. A. Novikov, ZhETF 47:1919 (1964) \[ Sov.Phys.-JETP 20:1990 (1965) \]. See also: M. D. Donsker, Proc. Conf. on the Theory and Applications of Analysis in Function Space (MIT, Cambridge, 1964).-P.17-30; K. Furutsu, J. Res. NBS 67:303 (1963); V. I. Klyatskin and V. I. Tatarskii, Teor. Mat. Fiz. 17:273 (1973).
17. V. I. Klyatskin, Statisticheskoe Opisanie Dinamicheskikh Sistem s Fluktuiruyushchimi Parametrami (Nauka, Moscow, 1975), in Russian; V. I. Klyatskin, Stokhasticheskie Uravneniya i Volny v Sluchaino-Neodnorodnykh Sredakh (Nauka, Moscow, 1980), in Russian; V. I. Klyatskin and V. I. Tatarskii,  Usp. Fiz. Nauk (Sov.Phys.- Uspekhi) 110:499 (1973).
18. I. M. Gel’fand and N. Ya. Vilenkin, Generalized Functions. Vol.4. Applications of Harmonic Analysis (Acad. Press, New York, 1964).