Gaussianity of Cosmic Velocity Fields and Linearity of the Velocity-Gravity Relation

We present a numerical study of the relation between the cosmic peculiar velocity field and the gravitational acceleration field. We show that on mildly non-linear scales (4-10 Mpc Gaussian smoothing), the distribution of the Cartesian coordinates of each of these fields is well approximated by a Gaussian. In particular, their kurtoses and negentropies are small compared to those of the velocity divergence and density fields. We find that at these scales the relation between the velocity and gravity field follows linear theory to good accuracy. Specifically, the systematic errors in velocity-velocity comparisons due to assuming the linear model do not exceed 6% in beta. To correct for them, we test various nonlinear estimators of velocity from density. We show that a slight modification of the alpha-formula proposed by Kudlicki et al. yields an estimator which is essentially unbiased and has a small variance.


INTRODUCTION
According to the gravitational instability paradigm, structures in the universe formed by the growth of small inhomogeneities present in the early Universe. Gravitational instability gives rise to a coupling between the density and peculiar velocity fields on scales larger than the size of clusters of galaxies, the largest bound objects in the Universe. On very large, linear scales, the relation between the density contrast δ and the peculiar velocity v in co-moving coordinates can be expressed in differential form, or in integral form, v(r) = g(r) .
Here,  (Lahav et al. 1991). Hence, comparing the observed density and velocity fields of galaxies allows one to constrain Ωm, or the degenerate combination β ≡ Ω 0.6 m /b in the presence of galaxy biasing (e.g., Strauss & Willick 1995 for a review). This comparison is done by extracting the density field from full-sky redshift surveys (such as the PSCz; Saunders et al. 2000), and comparing it to the observed velocity field from peculiar velocity surveys. The methods for doing this fall into two broad categories. One can use equation (2) to calculate the predicted velocity field from a redshift survey, and compare the result with the measured peculiar ve-locity field; this is referred to as a velocity-velocity comparison. Alternatively, one can use the differential form, equation (1), and calculate the divergence of the observed velocity field to compare directly with the density field from a redshift survey; this is called a density-density comparison. Nonlinear extensions of equation (1) have been developed by a number of workers (Nusser et al. 1991, Bernardeau 1992, Gramann 1993, Mancinelli et al. 1994, Mancinelli & Yahil 1995, Chodorowski 1997, Chodorowski & Lokas 1997, Chodorowski et al. 1998, Bernardeau et al. 1999, Dekel et al. 1999, Kudlicki et al. 2000; see also the discussion below). However, very little work has been done to test on what scales the integral relation, equation (2) holds, and how it might be extended into the mildly nonlinear regime; thus the motivation for this paper. Attempts have been made to carry out velocity-velocity comparisons with very large smoothing lengths (e.g., the inverse Tully-Fisher method of Davis, Nusser, & Willick 1996), and very small smoothing lengths (the VELMOD method of Willick et al. 1997;Willick & Strauss 1998). Davis et al. (1991), and more recently Berlind et al. (2000) discuss the systematic errors caused by mismatch of smoothing scales between the velocity and density fields. In this paper, we concentrate on the velocity-velocity comparison after smoothing on scales of 4 h −1 Mpc or larger: any smaller would be affected by strongly nonlinear effects, while larger smoothing would reduce the number of independent volumes over which the comparison could be made.
The amplitude of the velocity field smoothed on a given scale R with the window W depends on the density field power spectrum, P (k), as while for the density field the relation is as follows: (see the discussion in chapter 2 of Strauss & Willick 1995).
Here, W is the Fourier transform of the smoothing window. The absence of the k 2 term in equation (5) means that the velocity field is more heavily weighted by modes with low values of the wavenumber k, i.e., large scales which are fully in the linear regime. Therefore, we expect the relation between the velocity, v, and the gravity, g, to be closer to linear than that between the density and velocity divergence. KCPR have shown that the relation between (see equation 1) and δ (proportional to ∇ · g(r)) is nonlinear on small scales, and have proposed a semi-empirical formula accurately describing the dependence of θ δ ≡ θ|δ on δ: Here, the constant ǫ is introduced in order to keep θ = 0 as it must; it is approximately where σ 2 δ ≡ δ 2 denotes the variance of the density field. Thus, this approximation to ǫ gets progressively worse as the fluctuations become more nonlinear and in the following we use the exact value of ǫ, obtained numerically. KCPR have found that α = 1.9 is a good fit over a large range of smoothing scales. Solving Equation 7 for v in the case of an irrotational flow gives: where equation (8) gives an expression for θ δ .
In this paper, we investigate the nonlinearities of the relationship between the velocity and gravity fields using a set of grid-based simulations. The simulations are described in §2. In §3, we ask how well the probability distribution functions (PDFs) of the Cartesian components of v and g are fitted with a Gaussian. Given that the source fields for the velocity and gravity fields (i.e. the velocity divergence and density contrast respectively) have mildly non-Gaussian distributions, it is not a priori obvious what the distribution of the integral quantities should be. In §4 we directly measure the relation between v and g on various scales, and test the extent to which linear theory, or non-linear extensions to it, may hold. This is important in determining whether the existing velocity-velocity comparisons which use linear theory give biased results. We present our conclusions in §5. Two appendices contain derivations of results used in the text.

THE SIMULATIONS
We performed our simulations using the CPPA (Cosmological Pressureless Parabolic Advection) code (see Kudlicki et al. 1996, KCPR). Matter in this code is represented as a non-relativistic pressureless fluid, and its equations of motion are solved on a uniform grid fixed in comoving coordinates. Periodic boundary conditions are applied. Parabolic interpolation of the hydrodynamical state (like in the Piecewise Parabolic Method scheme, see Colella & Woodward 1984) assures low internal diffusion of the code and accurate treatment of high density and velocity gradients.
We chose to use a grid-based code rather than an Nbody code, because it produces a volume-weighted velocity field directly. This is important because equation (10) is a solution to equation (7) only when v is a potential (curlfree) field, and the mass-weighted velocity field exhibits curl even in the linear regime. 2 Moreover, the a priori unknown galaxy bias does not allow one to treat observational data as purely mass-weighted anyway.
All our simulations assume an Einstein-de Sitter universe. The relation between the velocity and the (scaled) gravity in the mildly nonlinear regime is insensitive to the Cosmological Density Parameter and Cosmological Constant, as demonstrated both analytically (Gramann 1993, Chodorowski 1997, Nusser & Colberg 1998; see also App. B.3 Table 1. Parameters of the studied models. Hereafter, the model with the cell size equal to 1.56 h −1 Mpc will be called the 'highresolution model', while the remaining models will be called 'lowresolution models'. of Scoccimarro et al. 1998), and by means of N-body simulations (Bernardeau et al. 1999), thus our results should be valid for any cosmology. Our simulations start from Gaussian density fluctuations with the linear APM power spectrum (Baugh & Efstathiou 1993, Baugh & Gaztañaga 1996. The initial velocity field is obtained from equation (2). The initial fields are evolved until the linear variance of the density in spheres of radius 8 h −1 Mpc, σ8, is unity. Then for subsequent analysis the output is selected for which nonlinear σ8 = 0.87 (the cluster normalization in the currently preferred model Ωm = 0.3, ΩΛ = 0.7; Eke, Cole & Frenk 1996).
In order to test the dependence of the results on the spatial resolution and box size we have studied three numerical models with parameters given in Table 1. To improve the statistics, we have performed six realizations of each of the models, with different random phases of the initial density field. We have also performed a few additional simulations with the timestep reduced five times; the results did not change noticeably.

THE MARGINAL DISTRIBUTIONS OF V AND G
The typical correlation length of the density field is of the order of 5 h −1 Mpc, while g is influenced by density fluctuations in a much larger region. For instance, the gravitational acceleration of the Local Group receives considerable contributions from shells up to at least 150 h −1 Mpc (see e.g. Rowan-Robinson et al. 2000). Thus, g, and similarly v, come from integration over an effective domain containing a large number of essentially independent regions. Hence, the Central Limit Theorem suggests that they should have close to Gaussian distributions. We test the Gaussianity hypothesis with our simulations. We plot the measured distribution functions for individual Cartesian components of the peculiar velocity and gravity fields, which we label v and g. It turns out that on mildly nonlinear scales, the distributions are well-fitted by Gaussians. Figures 1 and 2 show the distributions for 4 h −1 Mpc Gaussian smoothing. (As there is no preferred direction in space, the distributions must be even functions. Therefore, without the loss of information we plot them as functions of the absolute value of a Cartesian component.) The closeness to a Gaussian distribution is remarkable, especially for the velocity field.
A standard way of quantifying modest departures from Gaussianity is to decompose the PDF with an Edgeworth expansion: a leading Gaussian, plus correction terms with amplitudes proportional to the higher-order connected moments of the field (Longuet-Higgins 1963, Juszkiewicz et al. 1995, Bernardeau & Kofman 1995. As there is no preferred direction in space, the skewness of the PDF of the Cartesian coordinates of the velocity and gravity fields must equal zero. Thus the first non-vanishing connected moment of a Cartesian component of the velocity field is its kurtosis, Kv = ( v 4 − 3 v 2 2 )/ v 2 2 , and similarly for the gravity field. Therefore, the kurtosis measures the leading-order departure from Gaussianity of the fields (Kofman et al. 1994, Catelan & Moscardini 1994. Specifically, the Edgeworth expansion in our case reads: Here, P (µ) is the PDF for the variable v/σv or g/σg, where σ 2 v = v 2 and σ 2 g = g 2 are the variances of the velocity and gravity fields respectively, and φ(µ) = (2π) −1/2 exp(−µ 2 /2) is the standardized normal distribution. The symbol H4 denotes the fourth order Hermite polynomial. The quantity S4 is related to the velocity, or gravity, kurtosis and the variance of the density field in the following way: The quantities S4v and S4g are called hierarchical amplitudes. For a given smoothing scale, during the weakly non- linear evolution they are constant, independent of the normalization of the power spectrum.
We plot Kv and Kg as a function of the Gaussian smoothing radius, R, in Figure 3. We find that Kv is close to zero for smoothing between 3 and 15 h −1 Mpc. In contrast, the gravity field develops a detectable kurtosis. Still, this is substantially less kurtosis than the kurtosis of the gravity's source field, density. Measured from the high-resolution simulations (i.e., with the grid size 128 3 and the box size of 200 h −1 Mpc), for the 4 h −1 Mpc smoothing the density kurtosis is K δ = 22.2 ± 5.0, thus more than an order of magnitude bigger than the gravity kurtosis. (For the same smoothing, the kurtosis of the velocity divergence is K θ = 4.1 ± 0.6). Therefore, qualitatively the picture is clear: our simulations reveal substantially less kurtosis of the velocity and gravity fields than those of the velocity divergence and density contrast respectively; consequently the distributions of v and g are much closer to Gaussian. Our findings are inconsistent with the results of perturbative calculations of Catelan & Moscardini (1994), and consistent with the results of N-body simulations of the velocity field of Kofman et al. (1994). 3 Quantitatively, the agreement of the evolved velocity and gravity PDFs with the perturbatively motivated equation (11) is less than perfect. In the lower panel of Figure 2 we plot the difference between the measured distribution for a Cartesian component of the gravity field and the standardized Gaussian. As the thick solid line we plot the difference predicted from equation (11), with the value Kg = 0.83 (the mean value for all models). The predicted difference has somewhat too high amplitude.
Expansion (11) is an expansion in the density variance. In our simulations, for 4 h −1 Mpc smoothing σ 2 δ = 0.65. For such a big variance, higher-order contributions to formula (11) may simply not be negligible. We have checked this conjecture by analyzing the gravity field at an earlier output time. In Figure (4), we plot the gravity PDF for σ 2 δ = 0.25. Formula (11) fits the simulated distribution better.
We have also checked the perturbative scaling of the gravity (and velocity) kurtosis with the variance of the density field, i.e. equation (12). Although for a given smoothing scale, during weakly non-linear evolution the hierarchical amplitude S4g remains constant, it depends moderately on the smoothing scale via an effective spectral index of the power spectrum at this scale. However, we have performed an additional simulation for a power-law initial spectrum (n = −1), and found qualitatively similar results; a change of sign in the kurtosis at several megaparsecs.
Intrigued by this feature, we have studied the values of the kurtoses for even larger, apparently linear, scales. The results for the gravity kurtosis are shown in Figure 5. (The results for the velocity kurtosis are similar.) The values of the kurtosis clearly depend on the box size used in the simulations. For the box size equal to 200 h −1 Mpc, they tend, instead to zero, to an asymptotic value −1.5. We interpret this as a purely numerical effect. Namely, on large scales the velocity field of the simulation is dominated by a small number of Fourier waves, and thus the Central Limit Theorem does not force the distribution to be Gaussian, even in the initial conditions. Equations (5) and (6) show that this effect should be much stronger for g than for δ. On the largest scales, approaching the simulation box size, Kg converges to −1.5, i.e. the value expected for a single mode, as shown in Appendix A. 4 On scales shown in Figure 5, simulations with the 400 h −1 Mpc box size are much less affected by this effect, as expected. Why is it that this seems to affect scales down to ∼ 20 megaparsecs? Aren't there enough modes with scales ranging from ∼ 20 megaparsecs to the box size, to warrant Gaussianity? The answer to this question is provided by equation (5). To the variance of the velocity field on a given scale contribute all larger modes with amplitudes proportional to P (k), the density power spectrum. The spectrum of the APM galaxies, employed by us, has a maximum at ∼ 300 h −1 Mpc. This is larger than the 200 h −1 Mpc box size of some of our simulations; thus there are relatively few modes on the scale of the box size that dominate the velocity field on large scales.
To settle this issue definitely would require additional simulations with box size comparable to the Hubble radius (like the Hubble Volume simulations of the virgo consortium, Evrard et al. 2002). 5 However, the ultimate goal of this paper is to study the relation between velocity and gravity for 4 h −1 Mpc smoothing. Figure 3 shows that on such a small scale the velocity and gravity kurtoses are positive, thus most likely induced by nonlinear dynamics rather than due to finite volume of our simulations. However, given the limitations of the Edgeworth expansion, 6 we would like to have an alternative measure of non-Gaussianity. As such a measure we borrow the concept of negentropy from information theory (e.g. Cover & Thomas 1991, Papoulis 1991.
We define the entropy S[f ] of a probability distribution f : It can be shown that, for a given variance, the entropy is maximal for Gaussian fields. Hence the difference between the entropy of a given field and the entropy of a Gaussian of the same variance, the negentropy N [f ], can be used as a measure of departure from Gaussianity: The negentropy was calculated by numerically integrating the integral of equation (13), using a PDF binned over the range from −5σ to 5σ. We have found that this technique applied to a Gaussian with the same range and binning yielded a negentropy of less than 10 −5 , which assures us that our results are not affected by numerical effects. We plot the results as a function of Gaussian smoothing scale in Figure 6. The negentropy of the density field is compared to the negentropy of the gravity in the upper panel. The values for velocity and its divergence are compared in the lower panel. The negentropy of the velocity field is practically zero. The negentropy of the gravity field, even on the smallest scales, is extremely small compared to the negentropy of the density. The contrast between the negentropy of the gravity and the negentropy of the density is even greater than between the corresponding kurtoses. It is so because the density field, unlike the gravity field, has significant skewness, also contributing to the value of the negentropy. In other words, the density field is more non-Gaussian than the gravity field not only because it has bigger kurtosis, but also because it has non-vanishing skewness.
We conclude that on mildly non-linear scales, the non-Gaussianity of v and g is completely negligible compared to the non-Gaussianity of θ and δ. As stated earlier, this finding is consistent with the results of N-body simulations of the velocity field of Kofman et al. (1994). It is also consistent with the results of Gooding et al. (1992) and Scherrer (1992), where non-Gaussianity of the density field was (at least on large scales) due to initial conditions rather than nonlinear gravitational evolution.

THE ONE-COMPONENT V-G RELATION
In grid simulations, the shortest Fourier modes correspond to the Nyquist wavelength, of two cells. Smaller structures are not well resolved by the code. In low-resolution simulations (cell size 3.13 h −1 Mpc, see Table 1), the Nyquist wavelength is larger than the scale at which we would like to study the velocity-gravity relation (4 h −1 Mpc). In high-resolution simulations (cell size 1.56 h −1 Mpc), it is smaller. Therefore, to model the velocity-gravity relation at 4 h −1 Mpc scale, we will use only high-resolution simulations.
On scales of a few h −1 Mpc, the relation between δ and θ is non-linear. Nevertheless, since the probability distributions of the Cartesian components of both v and g are nearly Gaussian, we hypothesize that their joint probability distribution is a (bivariate) Gaussian as well. In Figure 7 we present the simulated joint PDF for v and g, measured on the 4 h −1 Mpc scale. It is indeed quite close to a bivariate Gaussian. The only substantial deviation is for the probability contour of 99%, along the gravity coordinate. There, positive kurtosis of the gravity field broadens the simulated isocontour with respect to the Gaussian one. Figure 7 also shows that v and g are strongly correlated. In App. B we show that in the bivariate Gaussian case, the relationship between v and g is purely linear; we now show that this is approximately the case in the simulations.
In the linear regime v = g and thus each of the Cartesian components of these quantities, which we denote v and g, respectively, are also equal. In practice, the relationship between these two quantities has some finite scatter (Fig. 7), thus we will characterize the relation between the mean v at constant g, v|g , and the converse, g|v . Fitting a straight line to v|g as a function of g gives a slope of 0.94. Therefore assuming pure linear theory on this scale would give a systematic 6% error in β. We now consider going beyond linear theory.
The relationship between v|g and g should be invariant to coordinate inversions, thus it must be an odd function; similarly for the relationship between g|v and v.
Hence, we shall adopt third-order polynomials as the simplest odd non-linear model, and use the unitless parameters c3 and d3 as a measure of the non-linearity of the relation. As the deviations from linear theory are small, we expect c1 and d1 to be close to unity, and c3 and d3 to be small. Since in velocity-velocity comparisons one reconstructs velocities based on the gravity field, the quantity v|g is more relevant, and we shall concentrate on it in this paper. Note however that relation (15) can be used to transform velocities in the first step of density-density comparisons such as POTENT (Dekel et al. 1999).
The parameters c1 and c3 as functions of the smoothing scale are shown in Figure 8, while Figure 9 shows d1 and d3. We have fitted them independently for the three coordinates in each of the six realizations of every numerical model; the lines are the averages over these 18 datasets and errorbars indicate their standard deviations. As expected, c3 and d3 are much smaller than unity. On large scales, c1 and d1 tend to the linear theory value with accuracy better than half per cent. Moreover, they remain close to unity even on small scales, which implies that the systematic bias in estimating β based on the linear theory approximation (2) is small.
To quantify this bias, in Figure 10 we show the parameters c1 and d1 when c3 and d3 are set to zero. For R > 20 h −1 Mpc, c1 and d1 deviate from unity by less than 2%. Therefore, at these scales one can apply linear theory with such good accuracy.
Note in Figures 8 and 9 that at small scales, highresolution simulations yield the highest values of c3 and 1−c1 (similarly for |d3| and 1 − d1), but overall the effects of resolution and of the simulation box size are not large.  We plot the mean predicted velocity as a function of the true velocity for a 4 h −1 Mpc Gaussian filter in Figure 11. The velocity predicted with the linear-theory model (i.e., the gravity) is plotted as a dashed line. Combined data from six realizations of our high-resolution model (three components each) are binned with respect to the predicted velocity and averaged in the bins. All velocities are in kilometers per second (the one-dimensional rms true velocity σv = 290 km· s −1 on this scale). A long-dashed line shows the velocity predicted by the polynomial formula (16). We see that linear theory, v = g, is a good fit for velocities up to about 2σv. The polynomial approximation works well up to about 3σv, but fails in the tails too. Therefore, we decided to seek a better estimator of velocity from density than a simple function of gravity.   KCPR have found that formula (8) is a good description of the tails of the relation between the density and the velocity divergence, so it should also work well on the integral level (eq. 10). It turned out, however, that at small scales a small modification is needed. Formula (8) implies that the coefficient of the linear term in density in the densityvelocity divergence relation is unity, θ δ = δ−r2(δ 2 −σ 2 δ )+. . ., where r2 = (α − 1)/2α, regardless the value of α, while we have found here that at small scales it is not strictly true. To cure this problem, we introduced an additional parameter, γ, so that For given α and γ the predicted velocity, v pred , was calculated from equation (10). The best values of α and γ were found by minimizing the quantity where vtrue is the true velocity and the sum is over all grid points (N 3 in total) and all Cartesian components. The resulting values of the parameters, for the smoothing scale 4 h −1 Mpc and different values of σ8 at different output times, are shown in Figure 12.
In the limit of linear theory, α = γ = 1. The offset ǫ, fully determined by α and γ, is then equal to zero. In general, on the integral level the value of ǫ is irrelevant, since its contribution to velocity in the integral in equation (10) averages out to zero. We see that for small σ8 the parameter γ tends to unity, while α decreases only weakly and remains well away from this value. This is consistent with the findings of KCPR, that α-formula (eq. 8, or eq. 17 with γ = 1) describes well the density-velocity relation in the weakly nonlinear regime (i.e., for σ8 smaller than, say, 0.3). Though the best-fit value of γ is close to unity even for large σ8, its inclusion markedly decreases X 2 .
In Figure 13 we plot the quantity χ(v) ≡ ∆X 2 /∆N , that is the square root of the contribution to X 2 from a given bin in velocity, divided by the number of points in that bin. (χ has units of km · s −1 .) We have X 2 = χ 2 (v) P (v) dv, where P (v) is the probability distribution function of a Cartesian component of the velocity field. Thus, X 2 is a number-weighted average of χ 2 . Dashed line shows χ for the velocity predicted according to the linear theory, longdashed -for the velocity approximated by a polynomial in the gravity, and solid -for the velocity predicted by the αγformula (inserting eq. 17 in 10). The linear-theory estimator of velocity yields the largest values of χ. The polynomial formula results in smaller χ in the tails, but in equal χ for 'typical' values of velocity (−σv ≤ v ≤ σv). The αγ-formula yields the smallest χ in the whole range of velocity. For typical values of velocity, i.e. where most of the data comes from, χ is slightly smaller than for the previous estimators. Furthermore, it grows only moderately in the tails.
In Figure 11, a solid line shows the mean velocity predicted by the αγ-formula as a function of the true velocity. The smoothing scale is that of velmod, i.e. 4 h −1 Mpc, and σ8 = 0.87. For this value of σ8, the best-fitted values of α and γ are α = 1.56 and γ = 1.054. We see that the αγformula yields a practically unbiased estimator of velocity. The difference between the predicted and true velocity in the tails is very small. Figure 14 demonstrates that the estimator of velocity from density provided by the αγ-formula is not only essentially unbiased but also has a small variance (as suggested already by Figure 13). The figure shows a joint PDF for the velocity predicted by the αγ-formula and the true velocity, for 4 h −1 Mpc smoothing. It is instructive to compare it with Figure 7, where v pred = g. We have found that for larger values of the smoothing scale the correlation between the predicted and true velocity is even higher.
In Figure 15 we show the dependence of α and γ on the smoothing scale, for σ8 = 0.87. For large R the parameter γ tends to unity, as expected. The parameter α grows from the value around 1.55 for 3 h −1 Mpc, to around 2.05 for 20 h −1 Mpc. This is in contrast with the results of KCPR, Figure 14. The joint probability distribution of a Cartesian component of the velocity predicted by the αγ-formula and the true velocity, both smoothed with a 4 h −1 Mpc Gaussian filter. Solid contours represent combined numerical results from six high-resolution simulations, dotted ones -a bivariate Gaussian of the same correlation coefficient and variances. The contours correspond to the probability levels of 68%, 95%, and 99%. where α ≃ 1.9 was a good fit for a large range of scales. We can see at least three likely sources of this discrepancy. First, the power spectrum used by KCPR was a pure power law, P (k) ∝ k −1 , while ours is not and the effective spectral index of our APM spectrum depends (slightly) on the smoothing scale. Secondly, KCPR studied the α-formula, equation 8, which is equivalent to our αγ-formula only when γ = 1. Finally, KCPR's fit was obtained for all bins in velocity having equal weight, while we weight the bins by the number of points. As a result, our fit is much less affected by the tails of the velocity distribution.
To summarize, we have tested various estimators of velocity from density. Nonlinear effects can be accounted for either at the integral level, employing equation (16), or at the differential level, employing equation (17). We have shown that the differential method works much better.

SUMMARY AND CONCLUSIONS
In this paper we have studied the cosmic velocity-gravity relation. First, we have measured the non-Gaussianities of the cosmic velocity and gravity fields, evolved from Gaussian initial conditions, by computing their kurtoses and negentropies. We have shown that, on scales of a few h −1 Mpc, the non-Gaussianities of the cosmic velocity and gravity fields are small compared to the non-Gaussianities of velocity divergence and density. A similar finding for the velocity field was reported by Kofman et al. (1994). Guided by this result, we have shown that the relation between v and g is nearly linear. Moreover, its proportionality coefficient is close to that predicted by linear theory. Specifically, we have shown that the systematic errors in velocity-velocity comparisons due to assuming linear theory do not exceed 6% in β. (Strictly speaking, if σ8 < 0.9 and the smoothing scale is not smaller than 4 h −1 Mpc). To correct for this small bias, we have tested various nonlinear estimators of velocity from density. We have shown that the αγ-formula (a slight modification of the α-formula proposed by KCPR) yields an estimator which is essentially unbiased and of small variance.
The smoothing of observed peculiar velocity data, with its sparse and noisy coverage of the velocity, is technically difficult. Thus velocity-velocity analyses like velmod compare unsmoothed peculiar velocities with minimally smoothed predicted velocity fields from redshift surveys. The finite resolution of a grid code does not allow us to test this effect; it would be worthwhile to repeat the calculations with a high resolution N -body code, or with CPPA enhanced with Adaptive Mesh Refinement. The former has been done by Berlind et al. (2000), but they did not separate the effects of non-linear evolution from the effects of different smoothing of the velocity and gravity fields.
An alternative method of estimating β from cosmic velocities is by comparing the cosmic microwave background dipole with the density dipole inferred from a galaxy redshift survey. In maximum-likelihood approaches to this problem (Strauss et al. 1992, Schmoldt et al. 1999, Chodorowski & Ciecielag 2002) the relevant quantity is the joint distribution for the velocity and gravity, which is commonly approximated as a multivariate Gaussian. We have shown that this indeed holds to fairly good accuracy (Fig. 7). Small nonlinear effects can be corrected for by a more thorough modelling of the joint distribution. There is, however, a better approach: from a redshift survey, instead of calculating the acceleration on the Local Group (i.e., the density dipole, or gravity), one can calculate the predicted velocity from the αγ-formula. The joint distribution for the predicted velocity and the true velocity is perfectly Gaussian (see Fig. 14). Moreover, the relation between these two variables is tighter than the relation between the velocity and the gravity. This implies that in the proposed method, the estimated value of β will be unbiased and its inferred errors will be smaller.

APPENDIX A: CALCULATION OF THE KURTOSIS OF VELOCITY FOR A SINGLE FOURIER MODE
Let us assume a single-mode distribution of the density field. All variables will depend on one dimension only. In the linear regime the density field is δ(x, y, z) = δ(x) = sin(x) .
(A1) (The amplitude of the wave is irrelevant here, since it will cancel out in the calculation of the dimensionless kurtosis.) Then, in the linear regime, [θ denotes the scaled velocity divergence, so in equation (A2) there is no f (Ωm, ΩΛ) term.] Hence, and v(x) = cos(x) .
Let P (v) denote the volume-weighted probability distribution function of a single component of the velocity field: Hence, and finally P (v) = π −1 1 − v 2 −1/2 .
Now let us compute the four first moments of this distribution: We can calculate the kurtosis of this field now: