# Experimental Cosmic Statistics I : Variance

## 1 Introduction

Measurements of higher order statistics in galaxy catalogs test theories of structure formation, the nature of the initial fluctuations, and the processes of galaxy formation. The power of such measurements to constrain theories, however, depends crucially on the detailed understanding of the errors. Usually it is tacitly assumed that the underlying distribution of events is Gaussian and thus the term “errors” becomes synonymous with the “variance”. Knowledge of the variance is sufficient only when the error distribution is Gaussian.

For statistics related to counts-in-cells a rigorous theory for the cosmic errors was presented in a suite of papers by Szapudi & Colombi 1996, hereafter SC, Colombi, Szapudi & Szalay 1998, and Szapudi, Colombi & Bernardeau 1999a, hereafter SCB. Nevertheless these calculations relied on approximations, for which the domain of validity could not be checked extensively until the arrival of the Virgo Hubble Volume Simulations. Moreover, the regime where the underlying cosmic distribution is Gaussian could not be examined previously. This paper addresses the first problem by studying the statistical errors and cross-correlations numerically, while a companion paper, Szapudi et al. (1999c, hereafter paper II), discusses the underlying distributions of statistics in their full splendour.

Let us consider a statistic $`A`$ measured in a galaxy catalog of volume $`V`$. The corresponding indicator is denoted by $`\stackrel{~}{A}`$. In practice, only one sample of our local universe is accessible. However, a frequentist numerical experiment can be performed in a large numerical simulation if a sufficient number $`C_{}`$ of galaxy catalogs $`_i`$ can be extracted from it. In each of them a value $`\stackrel{~}{A}_i`$, $`1iC_{}`$, can be measured.

For any statistic $`A`$ the cosmic distribution function $`\mathrm{{\rm Y}}(\stackrel{~}{A})`$ is the probability density of measuring the value $`\stackrel{~}{A}`$ in a particular finite realization. This distribution function can be approximately extracted from the $`C_{}`$ subsamples under the ergodic hypothesis. For simplicity, we dispense with the (logical) notation $`\stackrel{~}{\mathrm{{\rm Y}}}`$, and replace it in what follows with $`\mathrm{{\rm Y}}`$. This expresses the fact that we do not wish to enter one more level of complexity by considering the “error on the error” problem (SC) in greater detail. The smoothness and regularity of our measurements suggest that the number of realizations, which represent a two orders of magnitude improvement over any previous work, is large enough to provide an adequate determination of the quantities measured.

While in practice the function $`\mathrm{{\rm Y}}(\stackrel{~}{A})`$ is the fundamental quantity underlying all measurements, this paper concentrates on its first two moments; paper II examines its shape and skewness in detail.

In the following definitions, integrals are to be understood as summations of the estimator over the distribution function. The first moment of $`\mathrm{{\rm Y}}(\stackrel{~}{A})`$ is the spatial average

$$\stackrel{~}{A}\mathrm{{\rm Y}}(\stackrel{~}{A})𝑑\stackrel{~}{A}=\stackrel{~}{A}A,$$
(1)   
where it is assumed that the estimator $`\stackrel{~}{A}`$ is unbiased. The bias is negligible compared to the relative cosmic error in most meaningful cases (SCB) as illustrated later by practical examples. For completeness, however, the definition of the cosmic bias is

$$b_A\frac{\stackrel{~}{A}A}{A}.$$
(2)   
The second (centered) moment of the cosmic distribution is called the cosmic error,

$$(\stackrel{~}{A}A)^2\mathrm{{\rm Y}}(\stackrel{~}{A})𝑑\stackrel{~}{A}=(\stackrel{~}{A}A)^2(\mathrm{\Delta }A)^2.$$
(3)   
For a biased statistic, the variance should be centered around the biased average and not the true value. It can however be shown formally (SCB) that the above definition is valid to second order in $`\mathrm{\Delta }A/A`$ for any biased statistic.<sup>1</sup><sup>1</sup>1More precisely, to first order in $`(\stackrel{~}{x}_ix_i)(\stackrel{~}{x}_jx_j)`$ where $`\stackrel{~}{x}_i`$ denote the unbiased estimators from which $`\stackrel{~}{A}`$ is constructed in a non-linear fashion.

Finally, the cosmic covariance can be defined analogously to the variance as $`(\stackrel{~}{A}A)(\stackrel{~}{B}B)`$.

The theoretical results for the errors and cross-correlations are summarized below. If $`v`$ and $`V`$ are the cell and catalog volumes respectively, the cosmic error can be approximately separated into three components to leading order in $`v/V`$ (SC):

* The discreteness or shot noise error which is due to the finite number of objects $`N_{\mathrm{obj}}`$ in the catalog. It increases towards small scales and with the order of the statistics considered, but becomes negligible when $`N_{\mathrm{obj}}`$ is very large.
* The edge effect error is due to the uneven weight given to galaxies near the edges of the survey compared to those near the centre. It is especially significant on large scales, comparable to the size of the catalog.
* The finite volume error is due to fluctuations of the underlying density field on scales larger than the characteristic size of the catalog.

The next to leading order correction in $`v/V`$ is proportional to the perimeter of the catalog $`V`$. At this level of accuracy there are also correlations between the three sources of error (e.g., Colombi et al. 1999, hereafter CCDFS).

Colombi, Bouchet & Schaeffer (1995, hereafter CBS) investigated in detail the cosmic error on the void probability function. The groundwork for error calculations of statistics related to counts-in-cells is based on SC where the cosmic error for factorial moments<sup>2</sup><sup>2</sup>2e.g., Appendix A for definitions and notations. was evaluated analytically. SCB, extended the work of SC to cross-correlations, including perturbation theory predictions (e.g., Bernardeau 1996). The cosmic errors, biases (see also Hui & Gaztañaga 1998, hereafter HG) and covariances for cumulants $`\overline{\xi }`$ and $`S_N`$ were calculated as well. The main goal of this paper is to compare the analytical predictions of CBS, SC and SCB to measurements made in the VIRGO $`\tau `$CDM Hubble Volume simulation.

The exhaustive nature of the comparison that follows warrants the questions: is it meaningful to thrive for the detailed numerical understanding of the theory? How much of it is practically useful? Can it accurately estimate the errors on measurements in future surveys? While some of these questions were addressed in SCB, a brief account of supporting arguments is given next.

The analytics do take into account all possible theoretical errors, but systematics, such as those resulting from cut out holes, incompleteness from fiber separation, possible magnitude errors in the case of the 2dF, etc., could in principle corrupt the theory and introduce biases. These effects might even require detailed simulation of the survey. In the case of the UKST and Stromlo surveys such simulations were performed and compared with the predictions: the spectacular agreement surprised even the present authors (Hoyle, Szapudi, & Baugh 1999). Thus systematics do not dominate in all surveys; for another example, where cut out holes found to have insignificant effect on the cosmic probability distribution of the two-point correlations function see Kerscher, Szapudi, & Szalay (1999).

Moreover, the wide theoretical framework is flexible enough to incorporate all systematics, which have the effect of altering certain parameters, such as the factorial moments. In such case any bias can be corrected for.

There might be unforeseen systematics which have such complicated non-linear effect that cannot even be modelled by appropriate alteration of a set of parameters. While it would be difficult to anticipate whether these could dominate for a particular survey, it is still instructive to investigate the potential results in an ideal case, especially during design phase of the survey. Error calculations help optimizing geometry, sampling, and other parameters. During the design of the VIRMOS survey such considerations were taken into account (Colombi et al. 1999). These calculations as well as maximum likelihood analyses need to explore such a large region in parameter space, that they would typically be impractical to carry out with simulations.

In addition to applications to surveys, the theory can be applied reliably to assess significance of measurements in simulations where multiple runs would be too costly (e.g., Szapudi et al. 1999d). All these present and potential future applications motivate the detailed investigations performed in this article.

The exposition is organized as follows. § 2 describes the $`N`$-body data used for the purpose of our study. § 3 analyses the count-in-cells distribution function $`P_N`$, its cumulants $`\overline{\xi }`$ and $`S_N`$’s, and the scaling function of the void probability distribution $`\sigma \mathrm{ln}(P_0)/F_1`$. These quantities are measured in the full simulation as well as in $`C_{}=4096`$ subsamples. The accuracy of simulation is assessed by comparing the measurements to the non-linear Ansatz of Hamilton et al. (1991) improved by Peacock and Dodds (1996, hereafter PD), and to perturbation theory predictions (hereafter PT). The model of Fosalba & Gaztañaga (1998) and extended perturbation theory (hereafter EPT, see Colombi, Bernardeau, Bouchet & Hernquist 1997) are considered as well. § 4 extends these investigations to the cosmic error and the variance of the cosmic distribution function. A preliminary investigation of the cross-correlations is done for factorial moments and cumulants. The measurements are compared where possible to the theoretical predictions of SC, SCB and CBS, including extended perturbation theory. Finally § 5 recapitulates the results and discusses their implications. In addition, Appendix A gives a summary of the definitions and notations used in this paper for counts-in-cells statistics. It will be useful for the reader unfamiliar with these concepts.

## 2 The $`N`$-body data

The $`\tau `$CDM Hubble volume simulation (e.g., Evrard et al. 1999) was carried out using a parallel $`\mathrm{P}^3\mathrm{M}`$ code described in MacFarland et al. (1998). The code was run on 512 processors of the Cray T3E-600 at the Rechenzentrum in Garching.

Initial conditions were laid down by imposing perturbations on an initially uniform state represented by a “glass” distribution of particles generated by the method of White (1996). Because of the size of the simulation, a glass file of $`10^6`$ particles was tiled 10 times in each direction. As the initial glass file was created with periodic boundary conditions tiling does not create any non-uniformities at the interface between the tiles.

A Gaussian random density field was set up by perturbing the positions of the particles and assigning velocities to them according to the growing mode linear theory solutions, using the algorithm described by Efstathiou et al. (1985). Individual modes were assigned random phases and the power for each mode was selected at random from an exponential distribution with mean power corresponding to the desired power spectrum $`|\delta _k^2|`$. Unlike Efstathiou et al. (1985), however, the initial velocities were set up exactly proportional to the initial displacements, according to the Zel’dovich (1970) approximation. As shown by Scoccimarro (1998) this leads to larger initial transients. To compensate for this the simulation was started at a high redshift, $`z=29`$.

The cosmological model used for the simulation $`\tau `$CDM is described in more detail in Jenkins et al. (1998). The approximation to the linear CDM power spectrum (Bond & Efstathiou 1984) was used

$$|\delta _k^2|=\frac{Ak}{\left[1+\left[aq+(bq)^{3/2}+(cq)^2\right]^\nu \right]^{2/\nu }},$$
(4)   
where $`q=k/\mathrm{\Gamma }`$, $`a=6.4h^1`$ Mpc, $`b=3h^1`$ Mpc, $`c=1.7h^1`$ Mpc and $`\nu =1.13`$. The value of $`\mathrm{\Gamma }`$ was set equal to $`0.21`$. The normalization constant, $`A`$, is chosen by fixing the value of $`\sigma _8^2`$ (the linear variance of the matter distribution in a sphere of radius $`8h^1`$ Mpc at $`z=0`$). A value of $`\sigma _8=0.6`$ was motivated by estimates based on cluster abundances (White, Efstathiou & Frenk 1993; Eke, Cole & Frenk 1996).

The simulation was integrated using a leapfrog scheme as described in Hockney & Eastwood (1981), section 11-4-3. The simulation was completed in 500 equal steps in time. The softening used was 100 kpc/$`h`$ comoving Plummer equivalent - see Jenkins et al. (1998).

## 3 Counts-in-cells analysis : the underlying statistics

The count probability distribution function (CPDF) $`P_N`$ is defined as the probability of finding $`N`$ objects in a cell of volume $`v`$ thrown at random in the catalog. CPDF was measured in the whole simulation $``$ for cubic cells of size $`L_{\mathrm{box}}/512\mathrm{}L_{\mathrm{box}}/8`$, where $`L_{\mathrm{box}}=2000h^1`$ Mpc is the size of the simulation cube (see Table 1). Then the simulation cube was divided into $`16^3`$ contiguous cubic subsamples $`_i`$ of size $`L=125h^1`$ Mpc. $`P_N`$ was evaluated in each of these for $`L/512\mathrm{}L/2`$ (see Table 1). The successive convolution algorithm of Szapudi et al. (1999d, hereafter SQSL) allowed the determination of the CPDF on all scales simultaneously in only a few minutes of CPU on a workstation<sup>3</sup><sup>3</sup>3This estimate does not include the reading in of the file. with $`512^3`$ sampling cells. The accuracy is thus $`P_NP_{\mathrm{min},1}=1/512^37.45\times 10^9`$ for the measurement in $``$ and for each individual $`_i`$; the accuracy increases by averaging over all subsamples: $`P_NP_{\mathrm{min},2}=1/(512\times 16)^31.82\times 10^{12}`$. For $`4\mathrm{}63h^1`$ Mpc the measurements in $``$ and $`_i`$ overlap (Table 1). This is illustrated by Fig. 1, displaying $`P_N`$ as a function of $`N`$: the figure presents the CPDF extracted from both the full cube and averaged over all the sub-cubes. In the overlap region, the difference can be detected as slight irregularities of the high-$`N`$ tail from the full cube measurements. The figure suggests that at least on the smallest scales considered in $``$ (or each $`_i`$), our sampling is probably insufficient by the standards of SC. However, this does not affect significantly the calculations as indicated by the agreement of the moments measured in $``$ and those calculated from averages obtained from the subsamples. Therefore measurement errors will be neglected in what follows, i.e. infinite sampling is assumed. Note that this ideal can be achieved in practice for two-dimensional and small three-dimensional catalogs via the algorithm of Szapudi (1998), however, the present simulation is too large for this method.

The smallest scale considered is only $`2.4`$ times larger than the softening length $`\lambda _ϵ=100h^1`$ kpc. As discussed extensively in Colombi, Bouchet & Hernquist (1996), contamination by softening restricts the validity of the simulation on small scales. For spherical cells of radius $`R`$, at least $`R4\lambda _ϵ`$ should hold. For the cubic cells of the present simulation this condition translates to $`\mathrm{}6.5\lambda _ϵ0.65h^1`$ Mpc. Thus the two smallest cell sizes i.e. the two leftmost points could be contaminated by softening, a fact that should be borne in mind, especially when comparing with theoretical calculations which employ models motivated by dynamics. On the other hand, for statistical purposes the dynamics can be ignored and the simulation can be regarded as a set with prescribed statistics. Then the possible contamination is irrelevant at the level of the approximations taken in the next sections.

Another possible source of contamination could be in principle the anticorrelation introduced by the glass initial positions. The effect of this is, however, extremely small as evidenced by the measurement of $`\overline{\xi }`$ shown below.

Figure 2 displays the average correlation function $`\overline{\xi }`$ as a function of scale. By definition

$$\overline{\xi }\frac{1}{v^2}_vd^3r_1d^3r_2\xi (|r_1r_2|),$$
(5)   
where $`\xi (r)`$ is the two-point correlation function. In practice, it is obtained as the variance of the counts-in-cells, corrected for discreteness effects automatically via the use of factorial moments (e.g., see SQSL and Appendix A for the detailed description of the method used in this paper to obtain the cumulants including the variance from counts-in-cells). The measured $`\overline{\xi }`$ is compared with linear theory (dots) and with the non-linear Ansatz of Hamilton et al. (1991) improved by PD (dashes). As expected, the agreement with linear theory in the regime $`\overline{\xi }1`$ is excellent, even on the largest scales where the anticorrelations introduced by the glass initial condition could cause contamination. The two leftmost points are slightly below the dashes, because of softening effects as discussed above, otherwise the results are in perfect accord with theory.

Figure 3 plots the extracted cumulants $`S_N`$’s against $`\overline{\xi }`$. They are compared with predictions of various models, including perturbation theory (PT, dots). By definition (e.g., Balian & Schaeffer 1989a)

$$S_N=N^{N2}Q_N\overline{\xi }_N/\overline{\xi }^{N1},$$
(6)   
where $`\overline{\xi }_N`$ is the $`N`$-point correlation function averaged over a cell:

$$\overline{\xi }_N=\frac{1}{v^N}_vd^3r_1\mathrm{}d^3r_N\xi _N(r_1,\mathrm{},r_N).$$
(7)   
Perturbation theory predictions have been calculated for spherical cells by Juszkiewicz, Bouchet & Colombi (1993) for $`S_3`$ and extended to arbitrary order by Bernardeau (1994):

$$S_N(\mathrm{})=f_N(\gamma _1,\mathrm{},\gamma _{N2}),$$
(8)   
$$\gamma _i\frac{d^i\mathrm{log}\overline{\xi }}{(d\mathrm{log}\mathrm{})^i}.$$
(9)   
For example

$$S_3=\frac{34}{7}+\gamma _1,$$
(10)   
$$S_4=\frac{60712}{1323}+\frac{62}{3}\gamma _1+\frac{7}{3}\gamma _1^2\frac{2}{3}\gamma _2.$$
(11)   
The dots on Fig. 3 assume $`\gamma _i=0`$, $`i2`$. While this is incorrect in principle for a scale dependent spectrum such as $`\tau `$CDM, the long dashes on the left-hand panels prove that the contribution of $`\gamma _2`$ is insignificant. Higher order $`\gamma _i`$ terms, as discussed also by Baugh, Gaztañaga & Efstathiou (1995), have an even smaller effect and can be rightly neglected.

PT predictions are accurately fulfilled in the weakly non-linear regime. This confirms again numerous earlier works (see, e.g. Juszkiewicz, Bouchet & Colombi 1993; Bernardeau 1994; Juszkiewicz et al. 1995; Gaztañaga & Baugh 1995; Baugh, Gaztañaga & Efstathiou 1995; SQSL). In fact the textbook quality agreement with PT demonstrates the accuracy of the $`\tau `$CDM Hubble Volume simulation.

The dashes give the predictions obtained from extended perturbation theory (EPT, Colombi et al. 1997; see also Szapudi, Meiksin & Nichol 1996 for EPT applied to galaxy data, and Scoccimarro & Frieman, 1998 for “hyperextended” perturbation theory). EPT assumes that the same forms of the higher order moments are preserved in the highly non-linear regime. There $`\gamma _1`$ above is simply an adjustable parameter without any particular meaning, i.e.

$$\gamma _{1,\mathrm{eff}}=\gamma _1(S_3)=S_3\frac{34}{7},$$
(12)   
where $`S_3`$ is the measured one. With this value of $`\gamma _1`$ the $`S_N`$’s, $`N4`$, can be computed using equation (8) (with $`\gamma _i=0`$, $`i2`$). The dashed curves matches the measurements quite well even in the highly non-linear regime thereby reconfirming the efficiency of EPT (see also SQSL). The agreement is not expected to be absolutely perfect from this Ansatz: on Fig. 3, EPT tends to underestimate slightly the measured values of $`S_N`$ when $`1\overline{\xi }10`$.

The dynamic range in the upper left panel of Fig. 3 is narrower than in the lower left panel: on large scales the agreement between PT and measurement becomes less accurate for the $`S_N`$’s, especially if $`N`$ is large. This might be related to transients due to the initial setup of the particles on a glass perturbed by using the Zel’dovich approximation. On the one hand, the transients related to pure Zel’dovich should decrease the value of the $`S_N`$’s (e.g., Juszkiewicz et al. 1993 and Scoccimarro 1998) while, on the other hand, the anticorrelations due to the glass could have the opposite effect by decreasing $`\overline{\xi }^{N1}`$ more than $`\overline{\xi }_N`$. Although this problem was not examined in detail, the glass contamination on $`\overline{\xi }`$ appears to be inconsequential. Alternatively, finite volume effects can degrade the high-$`N`$ tail of the CPDF (e.g., Colombi, Bouchet & Schaeffer 1994; CBS; Colombi et al. 1996). In addition, it is worth reemphasising that two rightmost points are prone to errors caused by softening as discussed earlier.

The right-hand panels of Fig. 3 zoom in on the transition between the weakly and highly nonlinear regime. For comparison, PT (with $`\gamma _i=0`$, $`i2`$, dots), EPT (dashes), and the one loop perturbation theory of Fosalba & Gaztañaga (1998) (dots-long dashes) are displayed. The last model yields agreement with the extracted values of $`S_N`$ for $`\overline{\xi }1`$, or even larger when the order $`N`$ is high enough (see upper right panel). This affirms the success of one-loop perturbation theory (see also Lokas et al. 1996; Scoccimarro et al. 1998). Interestingly, EPT produces almost identical results to the spherical model when $`\overline{\xi }1`$.

Finally, figure 4 shows $`\sigma =\mathrm{ln}(P_0)/\overline{N}`$ as a function of scale, compared with EPT predictions. By definition (White 1979; Balian & Schaeffer 1989a; see also Appendix A)

$$\sigma =\underset{N=1}{\overset{\mathrm{}}{}}(1)^{N1}\frac{S_N}{N!}\left(\overline{N}\overline{\xi }\right)^{N1},$$
(13)   
where $`\overline{N}`$ is the average count in a cell. This function is thus sensitive to low order statistics when $`N_\mathrm{c}\overline{N}\overline{\xi }1`$, and to high order statistics when $`N_\mathrm{c}1`$. According to Fig. 4, EPT is an accurate Ansatz on small scales where $`\sigma `$ is close to unity and is dominated by low order $`S_N`$. It is a less precise approximation on the largest scales probed, as expected. Indeed, the rightmost point of Fig. 4 corresponds to where $`\overline{\xi }1`$ in Fig. 3. There EPT increasingly underestimates the $`S_N`$’s when $`N`$ is high. Note the remarkable power-law behavior of $`\sigma \mathrm{}^{D_0}`$, $`D_00.25`$, in agreement with the predictions of the scaling model of Balian & Schaeffer (1989a). This reflects a non-trivial (multi)fractal particle distribution (Balian & Schaeffer 1989b) with a Hausdorff dimension $`D_0`$. Such behavior was found in a standard CDM model by Bouchet, Schaeffer & Davis (1991). Subsequently, the fractal distribution with $`D_00.5`$ was established by Colombi, Bouchet & Schaeffer (1992).

## 4 The Cosmic Error

In the previous section we demonstrated that good agreement was obtained comparing measurements made on the $`\tau `$CDM Hubble Volume dataset with previous work regarding higher order clustering statistics. Having established the accuracy of the dataset this section concentrates on the the determination of cosmic errors and their comparison to the available theoretical predictions, where possible. In § 4.1 we summarise analytic calculations of the cosmic errors and their cross-correlations. From this follows a systematic study of the experimental cosmic error of low-order statistics, i.e. factorial moments $`F_k`$, $`1k4`$ (§ 4.2), and cumulants $`\overline{\xi }`$, $`S_3`$ and $`S_4`$ (§ 4.3) together with a thorough comparison with the theoretical predictions. Also in § 4.3 we discuss the cosmic bias of the cumulants. Then the void probability and its scaling function $`\sigma `$ are explored (§ 4.4) followed by the cosmic error on the CPDF itself (§ 4.5). Finally, in § 4.6, there is a preliminary investigation of the cosmic cross-correlations of factorial moments and cumulants.

In all subsequent figures, except for the cross-correlations, there are errorbars plotted on the symbols corresponding to measurements due to the finite number of realizations $`C_{}=4096`$. These measurement errors, proportional to $`1/\sqrt{C_{}}`$ (SC), are negligible for our simulation, and the errorbars are smaller than the size of the symbols in most cases. As discussed in the Introduction, we neglect the cosmic error on the determination of the cosmic error (which is due to the finite size of the Hubble Volume itself) because in practice it is insignificant.

### 4.1 Cosmic Error: Theoretical Predictions

Before making any comparison with the analytic predictions, we outline the main ideas in CBS, SC, and SCB – more details can be found in these papers. Spherical cells of radius $`\mathrm{}`$ are assumed throughout for simplicity.

The bivariate CPDF $`P_{N,M}(\mathrm{},r)`$ is the probability of finding $`N`$ and $`M`$ points in two cells of size $`\mathrm{}`$ at distance $`r=|r_1r_2|`$ from each other. According to SC the cosmic error is computed via a double integral of $`P_{N,M}(\mathrm{},r)`$ over $`r_1`$, and $`r_2`$, conveniently split according to whether the cells overlap or not:

give rise to the discreteness and edge effect errors (see Introduction). The locally Poissonian assumption (CBS, SC) enables the approximate representation of the generating function $`P(x,y)`$ for overlapping cells by using only the monovariate generating function $`P(x)`$, i.e. the calculation depends on $`\overline{\xi }`$, $`S_N`$, $`N3`$ and the average count $`\overline{N}`$.

generate the finite volume error (see Introduction). To simplify the writing of $`P_{N,M}(\mathrm{},r)`$, the distance $`r`$ is assumed to be large enough compared to the cell size such that the bivariate CPDF can be Taylor expanded (to first order) in terms of $`\xi (r)/\overline{\xi }`$. This approximation is surprisingly accurate even when the cells touch each other (Szapudi, Szalay, & Boschán 1992, Bernardeau 1996, hereafter B96). Three models are used: two particular but still quite general forms of the hierarchical model, SS and BeS, introduced by Szapudi & Szalay (1993a, hereafter SSa, 1993b) and by Bernardeau & Schaeffer (1992), respectively, and perturbation theory, hereafter PT (B96). See SC and SCB for more details. The former two models depend only on monovariate statistics, i.e. on $`\overline{\xi }`$ and $`S_N`$, $`N3`$ and $`\overline{N}`$. PT on the other hand is expressed in terms of $`\gamma _i`$, $`\overline{\xi }`$ and $`\overline{N}`$ (B96). In principle, PT is accurate only in the weakly non-linear regime, for which it was originally designed, but it can be extended to the nonlinear regime as well: for monovariate distributions, EPT was proposed by Colombi et al. (1997), as discussed and tested versus measurements in § 3. This Ansatz can actually be naturally generalized to the bivariate CPDF (Szapudi & Szalay 1997, SCB). Our version, denoted by E<sup>2</sup>PT, takes the measured (non-linear) value for $`\overline{\xi }`$, $`\gamma _{1,\mathrm{eff}}`$ from equation (12) and it assumes, as EPT, $`\gamma _i=0`$ for $`i2`$.

Except for the error on the void probability and its scaling function $`\sigma `$ detailed in CBS, the theoretical results shown in this section were computed to leading order in $`v/V`$, where $`v`$ is the cell volume and $`V=L^3`$ is the sample volume.

The calculation of the error on a statistics of order $`k`$ depends on $`\overline{N}F_1`$, $`\overline{\xi }`$, $`\overline{\xi }(\widehat{L})`$, the average of the correlation function over the survey (see below), and $`S_N`$, $`3N2k`$. PT is determined by $`\gamma _i`$, $`i2k2`$ (§ 3) and E<sup>2</sup>PT by $`\gamma _{1,\mathrm{eff}}`$ as explained above. In all cases, we use the measured value of $`\overline{N}`$. Other parameters are chosen as follows:

1. linear theory is employed to compute $`\overline{\xi }`$ and $`\overline{\xi }(\widehat{L})`$ \[the catalog is assumed to be spherical to simplify the calculation of integral (16) below\] while higher order statistics are evaluated according to eq. (8) with $`\gamma _i=0`$, $`i2`$.
2. : the experimental $`\overline{\xi }`$ is used (open symbols on Fig. 2). The quantity $`\overline{\xi }(\widehat{L})`$ is computed numerically with the non-linear Ansatz of PD discussed in § 3 (assuming that the catalog is spherical). For the $`S_N`$’s, the measurements (open symbols on left panels of Fig. 3) are used for $`\mathrm{}15h^1`$ Mpc. On larger scales, EPT is more appropriate to determine $`S_N`$, $`N4`$: the increasing inaccuracy of the $`S_N`$’s on large scales and for large $`N`$ require this procedure. It is justified all the more since when $`\overline{\xi }0.27`$ EPT matches quite well to the PT predictions (see Fig. 3).

There is a subtlety worth mentioning which concerns the finite volume error, proportional to the integral

$$\overline{\xi }(\widehat{L})=\frac{1}{\widehat{V}}_{r_{12}2\mathrm{}}d^3r_1d^3r_2\xi (|r_1r_2|).$$
(14)   
To leading order in $`v/V`$, this integral reads (CCDFS)

$$\overline{\xi }(\widehat{L})=\overline{\xi }_0(\widehat{L})\frac{8v}{\widehat{V}}\overline{\xi }_1(2\mathrm{}),$$
(15)   
with

$$\overline{\xi }_0(\widehat{L})=\frac{1}{\widehat{V}^2}_{r_1,r_2\widehat{V}}d^3r_1d^3r_2\xi (|r_1r_2|),$$
(16)   
$$\overline{\xi }_1(\mathrm{})\frac{1}{v}_r\mathrm{}4\pi r^2\xi (r)𝑑r.$$
(17)   
In the above equations, $`\widehat{V}`$ corresponds to the volume covered by cells of volume $`v`$ included in the catalog.

The next to leading order correction, $`\overline{\xi }_1`$, can be identified as a negligible correction to the edge effects for most practical purposes. Although it did not make a significant difference, we included this correction nonetheless.

### 4.2 Cosmic Error: Factorial Moments

Figure 5 presents the cosmic error measured for the factorial moments $`F_k`$, $`1k4`$. By definition

$$F_k(N)_kN(N1)\mathrm{}(Nk+1)=\underset{N}{}(N)_kP_N.$$
(18)   
The factorial moments directly estimate the moments of the underlying continuous density field: $`F_k=\overline{N}^k\rho ^k`$ where $`\overline{N}=F_1`$ is the average count (e.g., SSa). On Fig. 5, the dotted, dash, long dash and dotted-long dash curves correspond to SS, BeS, E<sup>2</sup>PT and PT.

All the models converge and agree quite well with the measurements on large scales $`\mathrm{}\mathrm{}_07.1h^1`$ Mpc, as expected, since PT predictions should be valid. In contrast, on small scales $`\mathrm{}<\mathrm{}_0`$ the models overestimate slightly the numerically obtained error, E<sup>2</sup>PT being the most accurate. It is worth remembering that the leftmost two points may be contaminated by smoothing effects and should not be over-interpreted. Nevertheless, the decrease of precision on small scales suggests that our assumptions (i) or (ii) in § 4.1 are becoming more and more approximate in the non-linear regime, i.e. either the local Poisson assumption or the particular hierarchical decompositions loose their accuracy. To test this idea the contribution of overlapping cells (edge $`+`$ discreteness effects) were separated from the contribution of disjoint cells (finite volume effects), as shown respectively as solid and dash-long dash curves on Figure 6,

which concentrates on E<sup>2</sup>PT (long dashes). Note that the solid curve represents the SS and BeS models as well. Finite volume effects appear to dominate on small scales because our subsamples are dense enough to suppress discreteness error as expected (SC). This pinpoints assumption (ii) as the source of inaccuracy. Note that naively one would suspect additional loss of precision in the Taylor expansion of the bivariate CPDF. However, the finite volume error is a double integral over all the cells included in the catalog and separated by more than $`2\mathrm{}`$. The contribution of close cells is small, especially when $`\mathrm{}/L`$ is small. Thus E<sup>2</sup>PT itself appears to break down in the non-linear regime (SS and BeS are even less accurate), at least for the particular experiment we are analysing. Despite that EPT itself fares quite well (Fig. 3), its simplest natural extension to bivariate distributions, E<sup>2</sup>PT, is less accurate, as noticed earlier by Szapudi & Szalay (1997) in connection with the cumulant correlators of the APM galaxy catalog. However, the accuracy of the calculation based on E<sup>2</sup>PT should be adequate for most practical uses, and future work on the representation of the bivariate distribution in the highly non-linear regime will result in increased precision.

The solid curves in Fig. 6 represent the main contribution of the cosmic error on large scales. Here, as expected (SC), the cosmic error is dominated by edge effects. Despite the fact that theoretical predictions were determined to leading order in $`v/V`$ and the largest scale considered is $`\mathrm{}=L/2`$, i.e. $`v/V=1/8`$, the agreement between theory and measurement is surprisingly good. CCDFS have computed the next to leading order contribution proportional to the perimeter $`V`$ of the survey. With this correction, which increases the cosmic error especially on the largest scales, next to leading order theory would be inferior to the leading order one. The reason is that the calculation of CCDFS assumes a perimetric curvature radius much larger than the cell size. This assumption, which is useful for deep galaxy surveys with small sky coverage, obviously fails for a compact catalogue such as this one where the cell size $`\mathrm{}`$ becomes comparable to $`L`$.

### 4.3 Cosmic Error and Cosmic Bias: Variance and Cumulants

So far only the full moments $`F_k`$ have been examined. The cumulants $`\overline{\xi }`$ and $`S_N`$, however, are the more physically motivated quantities. But the statistics of these is complicated by the fact that they are ratios. For example (see Appendix A)

$$\overline{\xi }=F_2/F_1^21.$$
(19)   
As is well known in statistics (e.g., HG, SCB) $`A/BA/B`$. In other words, the estimator

$$\stackrel{~}{\overline{\xi }}=\stackrel{~}{F}_2/\stackrel{~}{F}_1^21$$
(20)   
is biased. Note that this is a general feature for any statistic constructed from unbiased estimators in a non-linear fashion (e.g. SCB). However, SCB showed theoretically that the cosmic bias defined in the Introduction, given here by

$$b_{\overline{\xi }}(\stackrel{~}{\overline{\xi }}\overline{\xi })/\overline{\xi },$$
(21)   
is of same order of $`(\mathrm{\Delta }\overline{\xi }/\overline{\xi })^2`$ in the regime $`\mathrm{\Delta }\overline{\xi }/\overline{\xi }1`$. Similar reasoning applies to the $`S_N`$’s. Thus leading order theoretical calculations neglect the bias. This can be done safely in the domain of validity of the perturbative approach used to expand a non-linear combination of biased estimators. A reasonable criterion proposed by SCB for this domain is that the cosmic bias be small compared to the relative cosmic error which itself should be small compared to unity. For an arbitrary (possibly biased) statistic $`A`$ this reads

$$b_A\mathrm{\Delta }A/A1.$$
(22)   
Left panels of Figure 7 are analogous to Fig. 5

and show the measured cosmic error as a function of scale for the biased estimators of $`\overline{\xi }`$, $`S_3`$ and $`S_4`$. The middle panels show the absolute value of the cosmic bias (open symbols) compared to the cosmic error (filled symbols). For additional clarity, the cosmic bias is plotted in linear coordinates as well in the right-hand panels.

It is interesting first to compare the cosmic error for factorial moments and cumulants of same order. The discreteness error is negligible for the scaling regime and the statistics considered here. The cumulants fare better/worse than the factorial moments in the non-linear/weakly nonlinear regimes, respectively. The finite volume error, dominating on small scales, is the limiting factor for factorial moments, while the edge effect error, dominating on large scales, drives the errors of the cumulants. This is in full accord with the predictions of SCB which can be consulted for more details. 

The theoretical models on Fig. 7 use the analytic calculations of SCB and are computed analogously to Fig. 5, as explained in § 4.1. E<sup>2</sup>PT only is presented in the middle and right-hand panels. Again, it is worth remembering that the leftmost points are dangerously close to the limit of possible contamination from artificial smoothing effects introduced by the force softening.

For the variance $`\overline{\xi }`$, the theory systematically overestimates the errors and the cosmic bias, except for the latter on large scales. This is not at all unexpected in light of the previous findings on small scales, where the three models SS, BeS and E<sup>2</sup>PT loose precision. In the weakly non-linear regime, $`\mathrm{}>\mathrm{}_0=7.1h^1`$ Mpc, where perturbation theory is valid, this is somewhat disappointing. However, the dynamic range is limited by the criterion (22), which is hardly, if at all, fulfilled here. Hence the leading order perturbative approach is likely to be insufficient.

For higher order statistics $`S_3`$ and $`S_4`$, the theory again tends to overestimate the amplitude of the measured cosmic bias on small scales. On large scales, where the predicted $`|b_{S_k}|`$ presents a sudden turn-up, condition (22) breaks down, thus the theory is inapplicable. The measured cosmic errors, on the other hand, are in accord with the theory within the range of its validity. The agreement on small scales is even better for $`\mathrm{\Delta }S_k/S_k`$ than for $`\mathrm{\Delta }F_k/F_k`$, $`k=3,4`$. This, however, should not be over-interpreted, as it is probably a coincidence due to cancellation effects of the ratios $`S_3=\overline{\xi }_3/\overline{\xi }^2`$ and $`S_4=\overline{\xi }_4/\overline{\xi }^3`$.

The cosmic bias is always negative (right-hand panels of Fig. 7), i.e. the biased estimators tend to underestimate real values (SCB; HG). In this particular experiment, the measured cosmic bias is always dominated by the measured cosmic error as predicted by the perturbative approach, except for the largest scales. Here the cosmic bias can become of same order as the cosmic error. HG suggested that the cosmic bias should be corrected for when measuring cumulants. Whether this makes sense depends on the magnitude of the cosmic skewness, i.e., the skewness of the cosmic distribution function itself. This will be discussed in more detail by paper II. However, it is worth noting that function $`\mathrm{{\rm Y}}(\stackrel{~}{A})`$ is positively skewed and that its maximum corresponds to the most likely measurement. This is in general smaller than the average, $`\stackrel{~}{A}`$. Thus, as pointed out already by SC, the measured value $`\stackrel{~}{A}`$ in a finite sample is likely to underestimate the real value $`A`$ even if $`\stackrel{~}{A}`$ is unbiased. If the cosmic skewness and/or the cosmic variance are large compared to the cosmic bias, it is pointless to correct for the cosmic bias. Either of the above is true for most surveys, including the upcoming wide-field surveys such as the 2dF and SDSS, thus bias-corrected estimators are unlikely to be useful in the future.

### 4.4 Cosmic Error and Cosmic Bias: Void Probability and Scaling Function

Upper panel of Fig. 8 shows $`\mathrm{\Delta }P_0/P_0`$ as a function of

scale compared to the prediction of CBS (long dashes), with the finite volume error contribution (dashes-long dashes) and with the edge + discreteness contribution (solid curve). The agreement between theory and prediction is excellent.

The lower panel of Figure 8 corresponds to the scaling function $`\sigma `$. As for $`\overline{\xi }`$ and $`S_N`$, the indicator $`\stackrel{~}{\sigma }=\mathrm{ln}(\stackrel{~}{P}_0)/\stackrel{~}{\overline{N}}`$ is biased. This bias (open symbols) is of order $`(\mathrm{\Delta }\sigma /\sigma )^2`$ and can be neglected.<sup>4</sup><sup>4</sup>4The theoretical and measured errors displayed on bottom part of Fig. 8 correspond to the biased indicator. The agreement between theory and measurement is less impressive than for $`P_0`$, but this is mostly due to the difference of dynamic range covered by the error in upper and lower panels of Fig. 8. Moreover, the calculation of $`\mathrm{\Delta }\sigma /\sigma `$ by CBS is only approximate and could certainly be improved (see the discussion in CBS).

The errorbars about $`\sigma `$ are quite small: nearly an order of magnitude smaller than in Figs. 5 and 7. According to equation (13), $`\sigma `$ reflects the low order statistics when $`N_\mathrm{c}=\overline{N}\overline{\xi }1`$ ($`\sigma 1`$ in Fig. 4) and the high order statistics when $`N_\mathrm{c}1`$ ($`\sigma <1`$). From the point of view of the errors, function $`\sigma `$ is an excellent higher order indicator (as discussed earlier by CBS); it is better than the low order factorial moments or cumulants, at least in the non-linear regime $`\mathrm{}\mathrm{}_0`$. This fact alone unfortunately does not guarantee the usefulness of this statistic as various models of large scale structure formation could be degenerate with respect to the void probability. The thorough work of Little & Weinberg (1994) suggests that this is indeed the case. It is tempting although dangerous to extrapolate the results of their analysis to the function $`\sigma `$.

### 4.5 Cosmic Error: Counts-in-Cells

Upper panel of Fig. 9 shows the cosmic error in the

CPDF as a function of $`N`$ for the various scales considered in $`_i`$. The scale increases with the $`x`$-coordinate of the upper right part of each curve. In the lower panel $`\mathrm{\Delta }P_N/P_N`$ is represented in a similar manner as a function of $`N/N_{\mathrm{max}}`$ where $`N_{\mathrm{max}}`$ is the value of $`N`$ for which $`P_N`$ is a maximum. \[We did not display the (small) scales corresponding to $`N_{\mathrm{max}}=0`$ or $`N_{\mathrm{max}}=1`$\]. In agreement with intuition, the cosmic error reaches its minimum in the vicinity of $`NN_{\mathrm{max}}`$ and becomes increasingly large in the tails. Thus the shape of the CPDF near its maximum has the most power to constrain in terms of errors. Kim & Strauss (1998) have measured the cumulants $`S_3`$ and $`S_4`$ by fitting an Edgeworth expansion convolved with a Poisson distribution to the measured CPDF in the 1.2 Jy IRAS galaxy catalog. According to their recipe, the best determined part of the CPDF near the maximum was kept for the fit. Their maximum likelihood approach uses a simple model for the cosmic error, but their method is promising. Its main weakness is the necessity to make a strong prior assumption for the shape of the CPDF. A natural consequence is that the estimated errorbars on the measured cumulants are considerably smaller than with the standard methods.

### 4.6 Cosmic Correlations

So far this section has dealt only with the second moment of the cosmic distribution function, i.e. with the cosmic errors. For a full description in the Gaussian limit, however, the moments of the joint distribution function are needed. These moments form the cosmic (cross-correlation) matrix (SCB). It is defined as $`(\stackrel{~}{A}\stackrel{~}{A})(\stackrel{~}{B}\stackrel{~}{B})`$, where $`\stackrel{~}{A}`$ and $`\stackrel{~}{B}`$ are any counts-in-cells related indicators, for example $`A=F_k(\mathrm{})`$ and $`B=F_k^{}(\mathrm{}^{})`$, or $`A=\overline{\xi }(\mathrm{})`$ and $`B=S_N(\mathrm{}^{})`$, etc. A detailed theoretical investigation can be found in SCB (for $`\mathrm{}=\mathrm{}^{}`$). By definition, for two statistics $`A`$ and $`B`$, the correlation coefficient $`1\rho 1`$ reads as

$$\rho \frac{\delta \stackrel{~}{A}\delta \stackrel{~}{B}}{\mathrm{\Delta }A\mathrm{\Delta }B}\frac{(\stackrel{~}{A}A)(\stackrel{~}{B}B)}{\mathrm{\Delta }A\mathrm{\Delta }B}.$$
(23)   
The cosmic cross-correlation coefficient together with the errors form the full correlation matrix. The inverse of this is the central quantity for the joint probability distribution function in the Gaussian limit. As a preliminary numerical analysis, Figs. 10 and 11 present the correlation coefficients as functions of scale ($`\mathrm{}^{}=\mathrm{}`$) for factorial moments and cumulants, respectively. As in Fig. 5, the dots, dashes and long dashes show the theoretical predictions given by the SS, BeS and E<sup>2</sup>PT models, respectively, as computed by SCB. The computation of $`\delta \stackrel{~}{A}\delta \stackrel{~}{B}`$ in eq. (23) is analogous to that of the cosmic error (see SCB for more details). \[For $`\mathrm{\Delta }A`$ and $`\mathrm{\Delta }B`$, and to have completely self-consistent calculations, we take the theoretical results as well in eq. (23)\].

The agreement between theory and measurement is less convincing for the cosmic cross-correlations than for the cosmic error. This appearance is due partly to the linear coordinates of the figures which emphasize deviations, but nonetheless real.

On Fig. 10 there is a significant discord between theory and measurements for factorial moments in the middle top, middle bottom, and top right panels. On small scales, this result is quite natural: it is probably due to the inaccuracy of the models SS, BeS and E<sup>2</sup>PT employed to describe the underlying bivariate distributions (§ 4.2). In the weakly nonlinear regime, this discrepancy is apparently puzzling, since the predicted cosmic error matches perfectly the measurements (Fig. 5). The disagreement increases with $`|kl|`$, where $`k`$ and $`l`$ are the corresponding orders. On large scales, the cross-correlations are dominated by edge effects leading to the suspicion that the local Poisson approximation (SC, § 4.1) is becoming increasingly inaccurate with $`|kl|`$.<sup>5</sup><sup>5</sup>5This is not surprising: this approximation neglects local correlations. This is all the more inaccurate as the difference between the “weights” given to two overlapping cells, i.e. $`(N)_k`$ and $`(N)_l`$ for factorial moments, increases. Another although less likely possibility, is that the leading order approach in $`v/V`$ is insufficient and higher order corrections are necessary to calculate cross-correlations. It would go beyond the scope of this paper to analyse in detail these effects which are left for future research.

For the cumulants, in addition to the above arguments, our perturbative approach to compute cross-correlation allows only a narrow dynamic range for analytic predictions, defined by criterion (22). In Fig. 11, this condition is chosen for practical purpose to be

$$|b_A|\mathrm{\Delta }A/A1.$$
(24)   
This is necessary but not sufficient: the theory appears to disagree significantly with the measurements on large scales at the top left, lower left and lower middle panels of Fig. 11.

Despite some of the discrepancies, the general features of the cross-correlations are well described by the theoretical predictions. For instance the cross-correlation between two statistics $`A_k`$ and $`A_l`$ decreases with the difference between the orders $`|kl|`$ as predicted (SCB). In our particular experiment $`\overline{N}`$ is significantly correlated with $`\overline{\xi }`$, but only weakly (anticorrelated) with $`S_k`$, $`k=3,4`$. Similarly, $`\overline{\xi }`$ and $`S_3`$ are weakly, but $`S_3`$ and $`S_4`$ are strongly correlated. A detailed discussion on these effects can be found in SCB.

## 5 Summary and Discussion

In this paper we have studied experimentally the properties of the moments of the cosmic distribution function of measurements $`\mathrm{{\rm Y}}(\stackrel{~}{A})`$, where $`\stackrel{~}{A}`$ is an indicator of a counts-in-cells statistic. For a thorough examination of $`\mathrm{{\rm Y}}(\stackrel{~}{A})`$ itself the reader is referred to paper II also in this volume.

We examined the factorial moments $`F_k`$, the cumulants $`\overline{\xi }`$ and $`S_N`$’s, the void probability $`P_0`$, its scaling function, $`\sigma \mathrm{ln}(P_0)/F_1`$, and the count-in-cells themselves $`P_N`$. $`\mathrm{{\rm Y}}(\stackrel{~}{A})`$ was measured in the largest available $`\tau `$CDM simulation divided into 4096 cubical subsamples. In each of these many subsamples, $`\stackrel{~}{A}`$ was extracted, and its probability distribution function $`\mathrm{{\rm Y}}`$ was estimated with great accuracy. The main results of our analysis are the following:

1. The measured count-in-cells in the whole simulation, in particular the cumulants $`S_N`$, are in excellent agreement with perturbation theory predictions in the weakly nonlinear regime. This confirms the results of numerous previous investigations in an unprecedented dynamic range. The textbook quality agreement demonstrates the state of the art accuracy of the simulation. Similarly, the measurements confirm extended perturbation theory (EPT) in the full available dynamic range $`0.05\overline{\xi }50`$, for $`S_N`$, $`N10`$. In addition one loop perturbation theory predictions based on the spherical model (Fosalba & Gaztañaga 1998) were found to be an excellent description of the measured $`S_N`$ up to a $`\overline{\xi }1`$.
2. The variance of $`\mathrm{{\rm Y}}`$ is the square of the expected cosmic error, $`\mathrm{\Delta }A`$, in the measurement of $`A`$ in a subsample, identified with a realization of the local observed universe. The measurement of $`\mathrm{\Delta }A/A`$, for $`A=P_0`$, $`\sigma `$, $`F_k`$ and $`S_N`$ appears to be globally in good accord with the theoretical predictions of Colombi, Bouchet & Schaeffer (1995), Szapudi & Colombi (1996, SC) and Szapudi, Bernardeau & Colombi (1998a, SCB).

In the highly non-linear regime, the theoretical predictions of SC and SCB tend to overestimate the cosmic error slightly, except for the ratios $`S_3=\overline{\xi }_3/\overline{\xi }^2`$ and $`S_4=\overline{\xi }_4/\overline{\xi }^3`$. In the latter case, there are some cancellations and the agreement between theory and measurement is good even on small scales, but this is probably a coincidence. It appears thus that none of the three variants of the hierarchical model in SC and SCB, can give an accurate enough account of the non-linear behavior of gravitational dynamics for the bivariate distribution functions.<sup>6</sup><sup>6</sup>6As discussed in § 4.1., the analysis of the cosmic error indirectly probes the bivariate probability distribution function $`P_{N,M}(r,\mathrm{})`$ of having $`N`$ and $`M`$ galaxies respectively in two cells of size $`\mathrm{}`$ separated by distance $`r`$ (see, e.g. SC).

In the weakly non-linear regime, agreement between theory and predictions is excellent for the factorial moments, but less good for the cumulants, due to the limitations of the perturbative approach used to expand such ratios.

Nonetheless EPT yields the most precise overall agreement with theory for our particular experiment. On small scales $`1h^1`$ Mpc $`\mathrm{}4h^1`$ Mpc, EPT overestimates the errors perhaps by a factor of two in the worst case.
3. In addition to the cosmic errors, the cosmic bias, $`b_A`$, was studied in detail as well. An estimator is biased when its ensemble average is different from the real value: $`b_A\stackrel{~}{A}/A10`$. This is always the case when unbiased estimators are combined in a non-linear fashion to form a new estimator (SCB, Hui & Gaztañaga 1998, HG), such as the cumulants.

In agreement with SCB, the measured cosmic bias is of order $`(\mathrm{\Delta }A/A)^2`$ and thus negligible when the cosmic error is small. However, as for the errors, the theory tends to overestimate the bias in the non-linear regime. On large scales, where the cosmic bias becomes significant because of edge effects, the perturbative approach used by SBC to compute theoretical predictions is then outside of its domain of validity.

Note that in the regime where the cosmic bias is significant, the cosmic error is likely to be large. For instance, in the particular numerical experiment used in this paper, the cosmic bias was always smaller than the cosmic errors and in most cases negligible. Moreover, in the regime where the bias could be significant, the cosmic distribution function $`\mathrm{{\rm Y}}(\stackrel{~}{A})`$ is significantly positively skewed (paper II). This implies that the measured $`\stackrel{~}{A}`$ is likely to underestimate the true value even for an unbiased estimator. The result is an effective cosmic bias, at most of order $`\mathrm{\Delta }A/A`$. As already shown by SC, this effective bias can contaminate even unbiased estimators such as $`\stackrel{~}{F}_k`$ and $`\stackrel{~}{P}_N`$. As a consequence, it is pointless correcting for the cosmic bias, in contrast with the proposition of HG, unless it is done in the framework of a maximum likelihood approach which takes into account fully the effects of the shape of the cosmic distribution function.
4. To complete the analysis of second moments, a preliminary investigation of the cosmic correlation coefficients for factorial moments and cumulants was conducted. Together with the cosmic errors, these coefficients form the cosmic cross-correlation matrix which underlies maximum likelihood analysis in the Gaussian limit.

Theoretical predictions of SBC give good qualitative account of the measured correlation coefficients, although they become increasingly approximate with the difference between the corresponding orders. This is likely to be a consequence of the local Poisson assumption (SC) employed for analytic predictions.

Provided that the Gaussian limit is reached in terms of the error distribution, the formalism of SC and SBC allows for a maximum likelihood analysis of the CPDF measured in three-dimensional galaxy catalogs. Two preliminary investigations are currently being undertaken. Szapudi, Colombi & Bernardeau (1999b) reanalyse already existing joint measurements of $`F_1`$ and $`\overline{\xi }`$, and Bouchet, Colombi & Szapudi (1999) perform a likelihood analysis of the count-in-cells measured in the 1.2Jy IRAS survey (Bouchet et al. 1993).

Paper II probes the domain of the Gaussian approximation for the cosmic distribution function, together with preliminary investigations for the bivariate cosmic distributions $`\mathrm{{\rm Y}}(\stackrel{~}{A},\stackrel{~}{B})`$. As shown there, the Gaussian limit is reached when the relative cosmic error is small compared to unity. This is expected to hold for a large dynamic range in future large galaxy surveys such as the 2dF and the SDSS (Colombi et al. 1998).

Statistical analyses of weak lensing surveys are similar to counts-in-cells measurements (e.g., Bernardeau, Van Waerbeke & Mellier 1997; Mellier 1999; Jain, Seljak & White 1999). As a result, slight modification of the formalism of SC and SCB is fruitful to compute theoretical cosmic errors and cross-correlations (Bernardeau, Colombi & Szapudi, 1999).

Finally, it is worth to mention a few questions which were not addressed so far by the investigations presented in this paper. As light might not trace mass, the distribution of galaxies may be biased (not to be confused with the cosmic bias), and also realistic galaxy surveys are subject to redshift distortion. While the above results were obtained for the mass, note that the theory which served as a basis of comparison is quite general and was formulated to describe phenomenologically either the mass or the galaxies. It appears that there should be no qualitative changes introduced by biasing or redshift distortions (e.g., Szapudi et al., 1999e), thus the same theory can be used for the galaxies as for the mass, except perhaps with slightly different parameters, or underlying statistical models. In fact, two of the models (SS, BeS) were entirely motivated by the galaxy and not by the mass distribution; they are expected to be more accurate for realistic catalogs if used in a self-consistent fashion. The scaling properties underlying these models is even more accurate in redshift space, as is well known. EPT, on the other hand, was originally motivated by theoretical considerations of the mass distribution and numerical simulations (Colombi et al. 1997), and therefore it is no wonder that it is the most successful model for the mass (but see also Scoccimarro & Frieman 1998). Nonetheless, even EPT was found to be a fairly good model for the galaxy distribution, at least in the EDSGC survey (Szapudi, Meiksin & Nichol 1996), a possible indication that galaxies approximately trace mass after all. In addition, it is worth mentioning that biasing models can be nondeterministic, i.e. stochastic in nature, but this again does not introduce anything new qualitatively which could not be handled in the framework of the theory of SCB. Finally, the theory outlined in this paper was contrasted against measurements in a $`\tau `$CDM simulation. However, the analytical framework is general enough to accommodate any cosmological model, and there are no qualitative differences in this respect between different cosmologies with Gaussian initial conditions and hierarchical clustering. Thus repeating the same analysis for a different CDM-like cosmogony would be superfluous and inconsequential.

## Acknowledgments

The FORTRAN routine for computing $`S_N`$, $`3N10`$, using one loop perturbation theory predictions based on the spherical model was provided by P. Fosalba (see the right-hand panels of Fig. 3). We thank F. Bernardeau, P. Fosalba, C. Frenk, R. Scoccimarro, A. Szalay and S. White for useful discussions. It is a pleasure to acknowledge support for visits by IS and SC to the MPA, Garching and by SC to the dept Physics, Durham, during which part of this work was completed. IS and AJ were supported by the PPARC rolling grant for Extragalactic Astronomy and Cosmology at Durham.

The Hubble volume simulation data was made available by the Virgo Supercomputing Consortium (http://star-www.dur.ac.uk/frazerp/virgo/virgo.html). The simulation was performed on the T3E at the Computing Centre of the Max-Planck Society in Garching. We would like to give our thanks to the many staff at the Rechenzentrum who have helped us to bring this project to fruition. The FORCE package (FORtran for Cosmic Errors) used for the error-calculations in this paper is available on request from its authors SC and IS.

## Appendix A Definitions and Notations

The count probability distribution function (CPDF) $`P_N`$, gives the probability of finding $`N`$ objects in a cell of volume $`v`$ thrown at random in the catalog.

Factorial moments, $`F_k`$, are defined as follows

$$F_k(N)_kN(N1)\mathrm{}(Nk+1)=\underset{N}{}(N)_kP_N,$$
(25)   
where the falling factorial $`(N)_k`$ is defined in the first part of the equation. The $`F_k`$ are proportional to the moments of the underlying density field $`\rho `$ smoothed over the cell of volume $`v`$: $`F_k=\overline{N}^k\rho ^k`$ (SSa; assuming the normalization $`\rho =1`$), where $`\overline{N}`$ is the average count in a cell:

$$\overline{N}N=F_1.$$
(26)   
Counts-in-cells are related to quantities of dynamical interest, such as the (connected) $`N`$-point correlation functions, $`\xi _N`$ (e.g., Peebles 1980). The averaged $`N`$-point correlation function over a cell is given by

$$\overline{\xi }_N\frac{1}{v^N}_vd^3r_1\mathrm{}d^3r_N\xi (r_1,\mathrm{},r_N).$$
(27)   
This is the connected moment of the smoothed density field, $`\overline{\xi }_N=\delta ^N_\mathrm{c}`$ (with $`\delta \rho 1`$). The connected moments or cumulants of a Gaussian field are identically zero for $`N3`$. In this paper, normalized cumulants are defined as

$$S_N\frac{\overline{\xi }_N}{\overline{\xi }^{N1}},$$
(28)   
with the short-hand notation $`\overline{\xi }\overline{\xi }_2`$. By definition, $`S_1S_21`$, thus for second order $`\overline{\xi }`$ is used.

The quantities $`S_3`$ and $`S_4`$ are often called skewness and kurtosis in the astrophysical literature, although their definition differs slightly from the original usage in statistics. The reason for normalization in eq. (28) is dynamical. The $`S_N`$’s exhibit a weak scale dependence only due to the scale-free nature of gravity. In the highly nonlinear regime stable clustering is expected to set in, (e.g., Peebles 1980) and in the weakly nonlinear regime perturbation theory predicts approximate scaling depending on the initial fluctuation spectrum (e.g., Juszkiewicz, Bouchet & Colombi 1993; Bernardeau 1994).

The counts-in-cells generating function,

$$P(x)\underset{N=0}{\overset{\mathrm{}}{}}x^NP_N,$$
(29)   
writes (White 1979; Balian & Schaeffer 1989a; SSa)

$$P(x)=\mathrm{exp}\{\overline{N}(1x)\sigma [N_\mathrm{c}(1x)]\},$$
(30)   
where

$$N_\mathrm{c}\overline{N}\overline{\xi }$$
(31)   
is the typical number of objects in an overdense cell (e.g., Balian & Schaeffer 1989a), and

$$\sigma (N_\mathrm{c})=\underset{N=1}{\overset{\mathrm{}}{}}(1)^{N1}\frac{S_N}{N!}N_\mathrm{c}^{N1}.$$
(32)   
It is worth noticing that (White 1979; Balian & Schaeffer 1989a; SSa)

$$P(x)=P_0[\overline{N}(1x)],$$
(33)   
if the void probability is expressed in terms of average counts $`\overline{N}`$. The measurement of $`P_0`$ is particularly interesting since it probes directly the count probability generating function:

$$\sigma (N_\mathrm{c})=\mathrm{ln}(P_0)/\overline{N}.$$
(34)   
The exponential generating function for factorial moments,

$$F(x)=\underset{k0}{}F_k\frac{x^k}{k!}$$
(35)   
is directly related to $`P(x)`$ (SSa) through

$$F(x)=P(x+1).$$
(36)   
Combining eqs. (30), (32) and (36), one can obtain a useful relation between cumulants and factorial moments (SSa)

$$S_N=\frac{\overline{\xi }F_N}{N_\mathrm{c}^N}\frac{1}{N}\underset{k=1}{\overset{N1}{}}\left(\begin{array}{c}N\\ k\end{array}\right)\frac{(Nk)S_{Nk}F_k}{N_\mathrm{c}^k}.$$
(37)   
The state of the art practical recipe consists of measuring the CDPF with high oversampling (Sect. 3), computing the factorial moments from eq. (25), and finally calculating the cumulants from the above recursion eq. (37). This procedure eliminates the need for explicit discreteness correction.