text:
# Observational Constraints on Higher Order Clustering up to 𝑧≃1

## 1 Introduction

The evolution of the spatial distribution of galaxies is intimately related to the physical processes of galaxy formation, to the initial spectrum and subsequent gravitational growth of matter fluctuations in the early universe, and to the global geometry of space-time. Quantifying the galaxy distribution is, thus, fundamental to cosmology and has dominated extragalactic astronomy for the past two decades. The $`n`$point correlation functions provide a statistical toolkit that can be used to characterize the distribution.

The two-point correlation function is the most widely used statistic because it provides the most basic measure of galaxy clustering – the departure from a pure Poisson distribution. It is also popular because its execution is computationally straight forward. The two-point correlation function is defined as the joint moment of the galaxy fluctuation field, $`\delta _g`$, at two different positions

$$\xi _2=\xi =\delta _{g,1}\delta _{g,2},$$
(1)
where $``$ means ensemble average. The two-point correlation function yields a full description of a Gaussian distribution only, for which all higher order connected moments are zero by definition. The galaxy distribution, however, exhibits non-Gaussian behavior on small scales due to non-linear gravitational amplification of mass fluctuations, even if they grew from an initially Gaussian field. On larger scales, where the density field is well represented by linear perturbation theory, non-Gaussian behavior may still be present if the initial perturbation spectrum was similarly non-Gaussian. In addition, the process of galaxy formation is likely to introduce biases between the clustering properties of the dark and luminous matter. In the presence of such realities, higher order moments are required to obtain a full statistical description of the galaxy distribution and to provide discrimination between different biasing scenarios (e.g., , Fry & Gaztañaga 1993, Fry 1994, Jing 1997).

An accurate determination of higher order clustering statistics requires a large number of galaxies and, to date, the most accurate measurements have been derived from wide-area angular surveys, such as the APM (Szapudi et al. 1995; Gaztañaga 1994) and the EDSGC (Szapudi, Meiksin, & Nichol 1996), although recent redshift surveys are now becoming large enough to make interesting constraints (e.g., Hoyle, Szapudi, & Baugh 1999; Szapudi et al. 2000a). These surveys are, by design, limited to the study of the current epoch galaxy distribution. Deep surveys add a further dimension to the exploration of clustering by enabling the study of its evolution. Ideally one would like to have deep, wide redshift surveys available for such analyses but it’s observationally infeasible at the current time. Projected surveys thus still provide a unique way to study the evolution of clustering, especially if photometric redshifts are available. Several deep surveys have been used to study the evolution of low order galaxy clustering (Lilly et al. 1995; Le Fevre et al. 1995; Neuschaffer & Windhorst 1995; Campos et al. 1995; Connolly et al. 1996; Lidman & Peterson 1996; Woods & Fahlman 1997; Connolly, Szalay, & Brunner 1998) but most suffer from insufficient contiguous area and are therefore barely large enough to measure the two-point correlation function. The Deeprange survey (Postman et al. 1998; hereafter paper I) was designed to study the evolution of clustering out $`z1`$. The resulting catalog contains $`710,000`$ galaxies with $`I_{AB}24`$ derived from a KPNO 4m CCD imaging survey of a contiguous $`4^{}\times 4^{}`$ region. The photometric calibration of the catalog is precise enough over the entire survey to limit zeropoint drifts to $`\stackrel{<}{}0.04`$ mag that translates to a systematic error in the angular two-point correlation function, $`\omega (\theta )`$, of $`\stackrel{<}{}0.003`$ on a $`4^{}`$ scale and proportionally less on smaller scales (see paper I for details). Accurate measurements of the $`I`$-band number counts over the range $`12<I_{AB}<24`$ and the two–point angular correlation function up to degree scales are presented in paper I.

The size and quality of the Deeprange catalog are sufficient to enable reliable estimation of the 3–point and 4–point angular correlation functions as well. This paper presents the first attempt to constrain higher order correlation functions down to flux limits of $`I_{AB}=23`$, corresponding to an effective redshift of $`z0.75`$. We briefly review the astrophysics contained within the higher moments in section 2. Section 3 presents a description of the galaxy sample analyzed and the computational methods used to derive the statistics and section 4 summarizes the results. The implications of our calculations are presented in section 5. Technical issues concerning the fit of the three-point correlation are presented in Appendix A.

## 2 A Brief Review of Higher Order Statistics

We begin with a summary of the basic definitions of the higher order statistical methods used in this paper, highlighting their most fundamental properties. For a more in depth review, the reader should consult the references cited.

The observed large-scale structures in the local universe are characterized by a high degree of coherence (e.g. de Lapperent, Geller, & Huchra 1986; Shectman et al. 1996; da Costa et al. 1998) and some features, like the “Great Wall” have undergone asymmetric gravitational collapse (dell’Antonio, Geller, & Bothun 1996). Furthermore, the galaxy distribution traces the underlying dark matter in a non-linear way that may depend on time and scale. It is the non-zero higher order correlation functions that uniquely characterize such phenomena and allow discrimination between the observed galaxy distribution and a Gaussian distribution with the same variance, i.e. two-point correlation function. The $`N`$-point correlation function

$$\xi _N=\delta _1\delta _2\mathrm{}\delta _N$$
(2)
depends on a large number of parameters ($`3N`$ coordinates, minus 3 rotations and 3 translations). The number of parameters decreases if the function is integrated over part of the configuration space (the geometric distribution of the $`N`$-points is often referred to as their configuration, and the $`N`$-dimensional space describing it as the configuration space).

A simple model for the $`N`$-point correlation functions is the clustering hierarchy (e.g., Peebles 1980) defined as

$$\xi _N(r_1,\mathrm{},r_N)=\underset{k=1}{\overset{K(N)}{}}Q_{Nk}\stackrel{B_{Nk}}{}\stackrel{N1}{}\xi (r_{ij}),$$
(3)
where $`\xi (r)\xi _2(r)=(r/r_0)^\gamma `$, and $`Q_{Nk}`$ are structure constants. Their average is

$$Q_N=\frac{{\displaystyle \underset{k=1}{\overset{K(N)}{}}}Q_{Nk}B_{Nk}F_{Nk}}{N^{(N2)}},$$
(4)
where $`F_{Nk}`$ are the form factors associated with the shape of a cell of size unity (see Boschan, Szapudi & Szalay 1994 for details). For the three-point correlation function, our main concern here, the above form factors amount to a few percent only; if neglected then $`Q_3S_3/3`$ (see the definition of $`S_3`$ below).

If the integration domain is a particular cell of volume $`v`$, with the notation $`\overline{f}=_vf/v`$ for cell averaging, then the amplitude of the $`N`$-point correlation function can be expressed as

$$S_N=\overline{\xi }_N/\overline{\xi }^{N1}=\delta ^N_c/\delta ^2^{N1}.$$
(5)
The $`S_N`$’s, with a suitable normalization motivated by leading order perturbation theory, are commonly used to characterize the higher order, i.e. non-Gaussian, properties of the galaxy distribution in real surveys and N-body simulations. In the second half of the above equation the integration over the cell, i.e. smoothing (or filtering) is implicit. While the $`S_N`$’s do not retain all the information encoded in the $`N`$-point correlation functions, in particular their shape dependence, it is an extremely useful measure of clustering. It is directly related to the distribution of counts in cells (in the same cell $`v`$), as it is the cumulant or connected moment 111Ensemble averages over the connected component have the statistical property of additivity for independent processes and hence receive the name “cumulants”. The alternative terminology “connected” comes from the (Feynman) graph representation of the statistical processes. It can be shown that cumulants correspond to graphs that have exactly one connected component. The mathematical definition relies on a logarithmic mapping of the generating function of ordinary moments (e.g. Szapudi & Szalay 1993) thereof. If the shape of the cell $`v`$ is fixed, these quantities depend only on one parameter, the size of the cell. Note the alternative notation $`Q_N=S_N/N^{N2}`$. These two notations differ only in their normalizations: for $`Q_N`$ it follows from the hierarchical assumption (see later), and for $`S_N`$ from perturbation theory. For Gaussian initial conditions with power spectrum of slope $`n`$, leading order (tree-level) perturbation theory of the underlying density field in an expanding universe predicts the $`S_N`$’s. For $`N=3`$ (skewness) and $`N=4`$ (kurtosis) the prediction is (Peebles 1980; Fry 1984; Juszkiewicz, Bouchet, & Colombi 1993; Bernardeau 1994; Bernardeau 1996):

$`S_3`$ $`={\displaystyle \frac{34}{7}}(n+3)`$ (6)
$`S_4`$ $`={\displaystyle \frac{60712}{1323}}{\displaystyle \frac{62(n+3)}{3}}+{\displaystyle \frac{7(n+3)^2}{3}}.`$ (7)
These equations depend on scale through the local slope, $`n`$, of the initial power spectrum that, except for scale invariant initial conditions, varies slowly with scale in all popular cosmological models. According to perturbation theory and simulations, geometric corrections from $`\mathrm{\Omega },\mathrm{\Lambda }`$ are negligible, (e.g., Bouchet et al. 1995). The above results are valid on scales $`7h^1\mathrm{Mpc}`$. On smaller scales accurate measurements from $`N`$-body simulations exist (e.g., Colombi, Szapudi, Jenkins, & Colberg 2000, and references therein) but the measurements are subject to a particular cosmological model.

If the galaxy field, $`\delta _g`$, is a general non-linear function of the mass density field, $`\delta `$, the we can express this function as a Taylor series $`f(\delta )=b_1\delta +b_2\delta ^2/2+\mathrm{}`$. The $`b_N`$ are the non-linear biasing coefficients, of which $`b_1=b`$ is the usual bias factor connecting the two-point correlation function of galaxies with that of the dark matter as

$$\xi _g=b^2\xi .$$
(8)
For higher order correlation functions, and for the $`S_N`$’s, analogous calculations relate the statistics of the galaxy and dark matter density fields (Kaiser 1984; Bardeen et al. 1986; Grinstein & Wise 1986; Matarrese, Lucchin, & Bonometto 1986; Szalay 1988; Szapudi 1994; Matsubara 1995; Fry 1996; Szapudi 1999). For example, the result for third moment is (Fry & Gaztañaga 1993)

$$S_{3,g}=\frac{S_3}{b}+\frac{3b_2}{b^2}.$$
(9)
This formula is expected to hold on the same or larger scales as leading order (weakly non-linear) perturbation theory, i.e. on scales $`10h^1\mathrm{Mpc}`$. On smaller scales, Szapudi, Colombi, Cole, Frenk, & Hatton (2000) found the following phenomenological rule from N-body simulations

$$S_{N,g}=\frac{S_N}{b_{}^{2(N2)}},$$
(10)
where the function $`b_{}=b_{}(b)b`$. For $`b1`$ the typical effect of biasing is that it decreases the $`S_N`$’s. The above theory does not include stochastic effects, when the galaxy density field is a random function of the underlying dark matter field. Stochasticity typically introduces only a slight extra variance on the parameters of the theory and, thus, it will not be considered further in this paper.

## 3 Sample Definition and Analysis Methods

The construction of the galaxy catalog is described in paper I. Briefly, the $`I`$-band survey consists of 256 overlapping CCD images. The field of view of each CCD is 16 arcminutes but the centers are spaced 15 arcminutes apart. The 1 arcminute overlap enables accurate astrometric and photometric calibration to be established over the entire survey. Objects were identified, photometered, and classified using automated software. For consistency, we analyze the higher order clustering properties in the same magnitude slices222Magnitude ranges are given in the Kron-Cousins $`I`$band; the conversion to AB magnitudes is achieved by adding $`0.5`$ mag to the Kron-Cousins values. used to measure the two–point correlation functions in paper I.

### 3.1 Counts in Cells

As our goal is to describe the higher order clustering of galaxies, we adopt statistical measures that are closely related to the moments of the underlying density field. The estimation of the higher order correlation amplitudes follows closely the method described in Szapudi, Meiksin, & Nichol (1996), which can be consulted for more details. Only the most relevant definitions are given next, together with the outline of the technique.

The probability distribution of counts in cells, $`P_N(\theta )`$, is the probability that an angular cell of dimension $`\theta `$ contains $`N`$ galaxies. The factorial moments of this distribution, $`F_k=P_N(N)_k`$ (where $`(N)_k=N(N1)..(Nk+1)`$ is the $`k`$-th falling factorial of $`N`$), are indeed closely related to the moments of the underlying density field, $`N(1+\delta )`$, through $`(1+\delta )^k=F_k/N^k`$ (Szapudi & Szalay 1993).

The most common method employed to relate the discrete nature of the observed galaxy distribution to the underlying continuous density field is known as infinitesimal Poisson sampling, which is effectively a shot noise subtraction technique. This corresponds to the assumption that, for an infinitesimal cell, the number of galaxies follows a Poisson distribution with the mean determined by the underlying field. This assumption must be approximate for galaxies, especially on small scales, because of possible halo interaction or overlap. Nevertheless, on the scales we are studying, Poisson sampling should be a good approximation. Moreover, even on scales where no underlying continuous process exists, factorial moments are the preferred way to deal with the inherent discreteness of galaxies (e.g., Szapudi & Szalay 1997). In what follows, the continuous version of the theory will be used for simplicity, and factorial moments are implicitly assumed wherever continuous moments are used in the spirit of the above.

Next, the mean correlation function, $`\overline{\xi }=\delta ^2`$, and the amplitudes of the $`N`$-point correlation functions, $`S_N=\delta ^N_c/\delta ^2^{N1}`$, are computed 333$`_c`$ refers to ensemble averages over the connected component. The galaxy correlation function measured via counts in cells is actually a smoothed version of this function referred to as $`\overline{\xi }`$. If the galaxy correlation function is a power-law then the average correlation function, $`\overline{\xi }`$, is a power-law as well with the same slope but with a different (increased444The amplitude is higher because $`\overline{\xi }`$ is the average value of $`\xi `$ UP TO a given scale.) amplitude that is determined with Monte-Carlo integration. For square cells and galaxy correlations with the usual slopes ($`\gamma 1.7`$), the ratio $`\overline{\xi }/\xi `$ is typically less than a factor of 2. The calculation of $`\overline{\xi }`$ is an alternative way to measure the two-point correlation function and it is used in the computation of the normalization of the higher order cumulants, i.e the $`S_N`$’s.

In the sequence of calculating the $`S_N`$’s from counts in cells via factorial moments, the most delicate step is the actual estimation of the distribution of counts in cells in the survey. We begin by transforming celestial coordinates into equal area Cartesian coordinates to assure proper handling of the curvature of the survey boundaries on the sky. We then adopt a series of masks (e.g., around bright stars) to denote regions in the survey that should be excluded from analysis. These masks were defined during the object detection phase and help minimize the inclusion of spurious detections. The infinitely oversampling method of (Szapudi 1997) was used to estimate $`P_N`$ in square cells. This method completely eliminates the measurement error due to the finite number of sampling cells. It can be sensitive, however, to edge effects on larger scales. When all the survey masks are used, only small scales up to $`0.16^{}`$ could be studied because larger cells would nearly always contain an excluded region. As a sensible alternative, we eliminated all but the top 5% of the largest masks and repeated the analysis. This procedure extended the dynamic range of the angular scales probed and typically does not significantly alter the results on smaller scales. In addition to shrinking the possible dynamic range available to our measurement, the large number of masks compromise the accuracy of the error estimate as well: the geometry of the survey (as explained in more detail later) could not be taken into account at the level of including the geometry, position, and distribution of the masks.

### 3.2 Three-point correlation function

The moments of counts in cells (i.e., the $`S_N`$’s) are the simplest descriptors of the non-Gaussianity of a spatial distribution, both in terms of measurement and interpretation. However, there are alternative statistical measures, such as the $`N`$-th order correlation functions, that reveal more information about the galaxy distribution by incorporating information about the geometry of higher order correlations. Indeed, the $`S_N`$’s are just higher order correlation functions that have been integrated over a portion of the available configuration space, i.e. over the $`N`$-points within a given cell. Although averaging over the configuration space suppresses noise, it also erases some important information that is, by contrast, kept by the $`N`$-th order correlation functions. Specifically, the validity of the hierarchical assumption can be tested directly (e.g., the first term in equation A2) and ultimately the effects of bias and gravitational instability can be distinguished. It is therefore desirable to compute as many of the $`N`$-th order correlation functions as the survey volume will allow. Three-point correlation functions have been calculated for the Zwicky, Lick, and Jagellonian galaxy samples by Peebles & Groth (1975), Groth & Peebles (1977), and Peebles (1975), for Abell clusters by Tóth, Hollósy, & Szalay (1989, hereafter THS), for the APM galaxy survey by Frieman & Gaztañaga (1999), and for an IRAS-selected galaxy sample by Scoccimarro et al. (2000; in this case the calculation was performed in Fourier space). Only the measurements by Peebles and Groth were performed on scales comparable to us, the rest of the work concentrated on much larger angular scales.

A larger survey is needed to directly measure the three–point correlation function, $`\xi _3=\zeta `$, over a large dynamic range in scale than that required for the measurement of skewness, $`S_3`$. Physically, this can be understood: the total number of triplets available in the survey are divided into finer subsets in the case of $`\xi _3`$ than for $`S_3`$, in order to examine a larger number of parameters. A smaller number of triangles available for each bin of the three-point correlation function naturally causes larger fluctuations and, hence, larger statistical errors.

Szapudi & Szalay (1998, hereafter SS) proposed a set of new estimators for the $`N`$-point correlation functions using either Monte Carlo or grid methods. We denote these series of estimators SSRN and SSWN, respectively. They have as their primary advantage the most efficient edge correction of any estimator developed to date. The two-point correlation function was estimated both with the Landy-Szalay estimator (Landy & Szalay 1993; hereafter LS), $`(DD2DR+RR)/RR`$, and the grid indicator SSW2. We employ both methods primarily to verify that the SSW2 results agree with those derived from the well-tested and widely used LS estimator. The three-point correlation function was estimated by the grid indicator SSW3.

The estimators SSWN ($`N=2,3`$) were implemented as follows. If $`D`$ denotes the data counts in a sufficiently fine grid, and $`W`$ is the characteristic function of the grid, i.e. taking values of $`1`$ inside the valid survey boundary, and $`0`$ elsewhere, the SSW2 estimator for the two-point correlation function is

$$\omega _2=(DD2DW+WW)/WW,$$
(11)
where $`DD`$, $`DW`$, and $`WW`$ denote the data-data, data-window, and window-window correlation functions. This is analogous to the LS estimator. The SSW3 estimator for the three-point correlation function is

$$\omega _3=(DDD3DDW+3DWWWWW)/WWW.$$
(12)
The key difference between the LS and SSW2 estimators is that the latter statistic is analogous to Euler’s method instead of Monte Carlo integration. In fact, SSRN is the direct generalization of the relation for the two-point function, which is LS $`=`$ SSR2. For the grid estimators SSWN the random catalog R is replaced with the characteristic function of the survey estimated on a grid, $`W`$. The difference amounts to a slight perturbation of the radial bins.

The main advantage of using a grid instead of a random catalog is computational speed, which starts to become prohibitive for the three-point function (at least with current technology computers and the present algorithms). The CPU time for calculating two-point correlation functions with the LS estimator is on the order of hours for about $`n=10^5`$ galaxies, and it roughly scales as $`n^2`$. For the three-point correlation function $`n`$-times more CPU power is thus needed, i.e. it would take $`10^5`$ hours on a fast workstation to calculate the SSR3 estimator. The SSWN series of grid estimators scale again as $`n^N`$, but here $`n`$ is the number of grid points. A significant speed up was achieved for the SSW3 estimator by storing precalculated distances of the grid-point pairs in a look up table; population of the look up table required $`n^2`$ operations. This did not change the scaling of the actual three-point estimator, but did minimize CPU time as floating point operations are eliminated from the main calculation. With this method, the three-point function for all magnitude cuts could be calculated in about a CPU day on a $`256^2`$ grid.

For the two-point correlation function, where it was computationally feasible, we confirmed that SSW2 produced results that were consistent with the LS estimator on scales larger than a few grid spacings, as it should. To repeat the same test for the three-point function would have necessitated supercomputer resources (see discussion above). This was not deemed justifiable as it is clear that no qualitative differences are expected for the three-point function from the slight perturbation of the binning introduced by the grid estimators.

We computed $`\omega _3`$ for all magnitude slices up to about 1 degree scales. The star-galaxy classification accuracy degrades at faint magnitudes ($`I\stackrel{>}{}22`$) in ways that are detailed in paper I. In brief, a fraction of the faintest compact galaxies are often classified as stars. The statistical correction for this effect is to include some faint “stars” into the galaxy catalog based upon a selection function that is derived from an extrapolation of the star counts at brighter ($`I<20.5`$) magnitudes where classification is extremely reliable. The effect of the correction is to restore, at a statistical level, the correct galaxy-to-star ratios at faint magnitudes. While this introduces some minor contamination of the galaxy sample by stars, the stars are presumably randomly distributed on the angular scales being considered. Therefore, except for the highly improbable case when a star cluster is included, any stellar contamination dilutes the correlations by the usual factor $`f_s^2`$, where $`f_s`$ is the fraction of stars in the survey. At the magnitudes where the statistical classification corrections are required ($`I>21.5`$), the number of stars introduced into the galaxy catalog as a result of the correction is estimated to be less than 10% of the galaxy population. The effect of the dilution is, thus, negligible.

An alternative approach to dealing with faint object misclassification is to only include in the analysis those objects initially classified as galaxies. If the misclassified galaxies have the same clustering properties as the correctly classified galaxies then the correlation functions will not change significantly (small changes may occur as a result of cosmic variance). However, if the compact objects likely to be misclassified as stars cluster differently, their omission could introduce a bias. In that case, sampling the objects classified as stars statistically according to an estimated misclassification ratio, as discussed above, is the remedy.

Lastly, one can compute the correlation functions using all detected objects regardless of their classification. In this case, one can again estimate the “true” galaxy correlation function by extrapolating the bright stellar number counts and multiplying the correlation amplitude by the factor $`N_{Obj}^2/(N_{Obj}N_{Star})^2`$ where $`N_{Obj}`$ is the total number of objects in a given magnitude bin and $`N_{Star}`$ is an estimate of the number of stars in the same bin.

To determine if misclassification effects are a significant source of bias, we estimated the two-point correlation function using all three strategies in all slices. The four brightest slices ($`17I<21`$) exhibit no discrepancies between any of the methods, indicating that misclassification is negligible in this flux range. For the deepest slice ($`22I<22.5`$), the results based on the statistical inclusion of faint stars according to an estimate of the misclassification rate is in excellent agreement with the corrected correlation function derived using all detected objects. Agreement with the results using only objects classified as galaxies was marginal. The $`21I<22`$ slice is an intermediate case. These results suggests that correction for misclassification is essential in this survey for $`I\stackrel{>}{}21.5`$.

For the deepest slices, we have thus performed the estimation of the three-point correlation functions using the first statistical correction discussed above (results from this method are denoted with the letter “s” in Tables 1 and 3) as well as an estimation based on using only those objects classified initially as galaxies. The difference in the results between the two methods provides an upper estimate of the systematic error introduced by star-galaxy misclassification. For the two deepest slices ($`21I<22,22I<22.5`$), the statistical procedure should yield the correct results.

## 4 Results

### 4.1 Counts in Cells

The $`s_3`$ and $`s_4`$ results are shown in Figure 1. In what follows lower case letters, such as $`s_N`$, denote the observed angular quantities and uppercase characters, such as $`S_N`$, represent the deprojected, spatial statistics. In each plot, open symbols show the exact estimates when all cells containing survey masks are excluded (hereafter referred to as the measurement type “E”) and closed symbols show the values when only cells containing the largest 5% of the masks were excluded (hereafter referred to as measurement type “M”). The magnitude cuts are shown in the lower right hand corner of each plot. The triangular symbols refer to $`s_3`$ while the rectangular ones to $`s_4`$. The $`s_N`$ data ($`N=3,4`$) and their estimated uncertainties are given in Table 1. The results in Table 1 for $`I<21`$ are for the “M” type measurements whereas those for $`I21`$ are for the “E” type measurements. We provide the “E” type measurements for the faintest slices because they provide the best fidelity and accuracy in the presence of any faint end systematic errors in the catalog at the expense of angular scale dynamic range (see discussion below).

For the scales considered, deprojection was performed using the hierarchical Limber’s equation:

$$s_N=R_NS_N.$$
(13)
(see e.g., Szapudi, Meiksin, & Nichol 1996). Table 2 lists the $`R_N`$’s for a typical choice of the luminosity function. The $`I`$band luminosity function is assumed to be the same evolving Schechter function as in paper I with $`M^{}(z)=20.9+5logh\beta z`$, $`\alpha =1.1`$. Our fiducial parameters for the deprojection are $`\beta =1.5`$, and $`\gamma =1.75`$ for the slope of the two-point spatial correlation function. The cosmological parameters adopted are $`h=0.65`$, $`\mathrm{\Omega }_o=0.3`$, and $`\mathrm{\Lambda }=0`$. The results are independent of the normalization of the luminosity function, and robust against variations in these parameters within the accepted range, or even beyond. For instance, the effect of varying all parameters simultaneously over the ranges $`\beta =12,\gamma =1.51.75,h=0.50.75,\mathrm{\Omega }=0.31`$ is a $`1015`$% change in $`R_3`$. A non-zero cosmological constant $`\mathrm{\Lambda }`$ was not considered; the effect of a reasonable non-zero $`\mathrm{\Lambda }`$ on our estimates of the higher order correlation statistics is expected to be bracketed by the changes in $`\mathrm{\Omega }`$ tested above. Table 2 also gives $`R_3`$ for $`\beta =1`$ and for $`\mathrm{\Omega }=1`$ to explicitly demonstrate the relatively weak dependence of the projection coefficients on these parameters.

The problem with the “E” type measurement is that the large number of excluded cells severely limits the maximum scale for which we can derive constraints. Indeed, on scales above 10 arcminutes any cell would intersect a masks with almost unit probability. With “M” type measurements, we can increase the maximum spatial scale by a factor of 2 at the price of contaminating the estimate slightly. The typical effect of such a contamination is a slight decrease in the $`s_N`$’s, although it depends on the exact spatial and size distribution of the masks and, as such, is extremely complicated to correct for or predict. The complex distribution and approximately power law size distribution of the masks prevented the application of counts in cells on larger scales, still relatively small compared to the characteristic size of the survey. This unfortunate phenomenon can be understood in terms of the masks effectively chopping up the survey into a large number of smaller surveys. Nevertheless, the large area of the Deeprange survey is still essential for the calculation of higher order statistics, as otherwise the finite volume errors would render the measurement of the $`s_N`$’s meaningless, even on the present small scales. However, the contamination decreases as the scale increases, since the area-ratio of the ignored mask is decreasing, and, thus, if the “M” and “E” type measurements agree out to a given angular scale, then the “M” measurements on larger scales can be considered relatively free of contamination whereas the “E” measurements may show edge effects. The two faintest slices, however, exhibit some disagreement between the “E” and “M” measurements on the smallest scales (see Figure 1) and therefore it is the “E” measurements that should be considered the most reliable results in these cases.

#### 4.1.1 Error Estimation

The errors of the measurements (columns 4, 6, 9, and 11 in Table 1) are estimated using the full non-linear error calculation by (Szapudi & Colombi 1996, Szapudi, Colombi, & Bernardeau 1999), obtained from the FORCE (FORtran for Cosmic Errors) package. All the parameters needed for such a calculation were estimated from the survey self-consistently. These are the perimeter and area of the survey, the cell area, the average galaxy count and the average two-point correlation function over a cell, and the average two-point correlation function over the survey (estimated from the integral constraint fits of paper I). The higher order $`S_N`$’s up to 8 were needed as well. They were calculated from extended perturbation theory (Colombi et al. 1997; Szapudi, Meiksin, & Nichol 1996). This theory is defined by one parameter, $`n_{eff}`$, which is identical to $`n`$, the local slope of the power spectrum on weakly non-linear scales, and is purely phenomenological on small scales. It can be determined from $`S_3`$ according to Equation 10, and was measured from the average $`S_3`$ estimated from the survey itself. In addition, the package uses a model for the cumulant correlators, the connected joint moments of counts in cells. The model by Szapudi & Szalay 1993 was chosen, as this was shown to provide a good description of projected data (Szapudi & Szalay 1997).

To estimate the uncertainties in the $`S_N`$’s, the FORCE package uses a perturbative procedure that is numerically accurate for small errors. When the fractional errors become $`>1`$, the value computed is meaningless and only indicates, qualitatively, that the errors are large. In those instances we renormalized the plotted error to $`100\%`$ and the corresponding entry in Table 1 is given as “L”.

The accuracy of the fully non-linear error calculation method has been thoroughly checked (see Szapudi, Colombi, & Bernardeau 1999; Colombi et al. 1999). However, one cautionary note is in order: although the calculation takes into account edge effects to second order (the ratio of the survey area to survey perimeter), the use of masks introduces complex geometrical constraints and extra edge effects. This is the main reason why the counts in cells method is limited to fairly small scales; the errors associated with these effects are not accurately modeled by the FORCE.

For our two faintest bins ($`I>21`$), we also provide an estimate of the fractional systematic errors associated with uncertainties in the statistical correction for object misclassification. These are given as the second error value in the parentheses in Table 1. The systematic errors are just the standard deviations in the mean $`s_3`$ and $`s_4`$ computed from 10 realizations of the statistical misclassification correction described in §3.2. The fractional systematic errors are shown by the right shifted error bars in Figure 1. Systematic errors for the brighter bins are negligible ($`\stackrel{<}{}2`$%).

### 4.2 Three-point correlations

Figure 2 shows our measurements of the angular three-point correlation function, $`z(a,b,c)`$, as a function of the hierarchical term $`\omega (a)\omega (b)+\omega (b)\omega (c)+\omega (c)\omega (a)`$, where $`\omega `$ is the two-point correlation function. The magnitude ranges are shown in the lower right hand corner of each subplot. The three-point correlation function was estimated in logarithmic bins for each side of the corresponding triangle. The figure makes no attempt to display shape information, although it is available from our estimator. The interpretation of projected shape space is extremely difficult because several different three-dimensional triangles project onto identical angular configurations. In the weakly non-linear regime, it is possible to project down the firm predictions for triangular configurations from a given theory and compare them with the observed angular three-point correlation function (e.g., Frieman & Gaztañaga 1999). For the small scales we are considering, however, a typical shape would deproject to a mixture of triangular galaxy configurations that are in the highly and mildly non-linear regimes. Therefore, we have chosen a different, entirely phenomenological, approach pioneered by THS. This essentially consists of fitting the three-point correlation function with a class of functions motivated by a general expansion; details of the fitting procedure and error estimation are presented in Appendix A.

Table 3 summarizes the results of the three-point correlation function fitting procedure. Column 2 of this table gives the number of degrees freedom in the fit. We perform fits with the third order terms ($`q^{111}`$ and $`q^{21}`$) free or locked to zero - hence the two different sets of results for each magnitude bin. On scales where the correlations are small ($`\xi _210^4`$) and the relative fluctuations of the estimator are becoming increasingly large, the division by the hierarchical term (see Appendix A) becomes unstable and produces a multitude of outliers that are easily identified in Figure 2. From the measurement of $`S_3`$ from Deeprange, $`q_3`$ can be safely bracketed with $`0.1<q_3<10`$, at least for the brighter magnitude cuts. These are fairly conservative limits and are in agreement with other measurements (e.g., Szapudi, Meiksin, & Nichol 1996) as well. These limits are displayed as long dashes on the figures representing three-point correlation functions, and were used in the fits to eliminate outliers. For a few of the deepest slices the fit could be sensitive to the exact placement of the lower $`q_3`$ outlier cut. All the fits were, thus, reevaluated with the limits $`0.01<q_3<10`$ as well. These extended fits are given in columns 7 through 10 in Table 3. The $`\chi ^2`$ values in columns 3 and 7 of Table 3 are the full (unreduced) values.

Any significant difference between the fits in Table 3 for a given magnitude interval suggests large systematic errors. All the fits were performed by standard computer algebra packages, and the parameters were found to be robust with respect to the initial value assigned to them. The results in Table 3 can be compared with those in Table 1 using the fact that $`s_33q_3`$. This equation is not exact, because the shape of the cell influences the integral performed in the definition of $`s_3`$. The resulting form factors amount to only $`23\%`$, which is the accuracy of the approximation (Boschan, Szapudi & Szalay 1994). Our fitting procedure cannot yield an accurate error on the fitted parameters, since the input errors for the $`\chi ^2`$ were not accurately determined. Nevertheless, the FORCE error bars should give a good indications from the previous tables for the expected statistical uncertainty since the hierarchy is approximately true. Furthermore, any differences in the results introduced either by neglecting the statistical correction for star/galaxy misclassification or by altering the $`q_3`$ fit limits (see Table 3) should give a reliable indication of the systematic errors.

We did not attempt to compute systematic errors for the three-point correlation function because it would have been very time consuming computationally. We note, however, that one can gauge the amplitude of any systematic errors by comparing the results with and without the statistical misclassification correction in Table 3. Furthermore, the realization of the misclassification-corrected galaxy catalog used in the three-point correlation function computation has $`s_3`$ values that are very similar to the values derived from the counts-in-cells method.

## 5 Discussion and Summary

Our constraints on the higher order clustering of the brightest galaxies in our survey are in good agreement with previous work. This can be seen in Figure 1 for our $`17I<18`$ subsample where we also display the results reported by Szapudi & Gaztañaga (1998, see also Gaztañaga 1994) for the APM $`17b_j<20`$ galaxy sample and those reported by Szapudi, Meiksin & Nichol (1996) for the EDSGC with dash-dots, and long dash-dots, respectively. Comparison of our $`17I<18`$ subsample with the above larger surveys is justifiable because the projection coefficients only vary slowly with depth. The APM and EDSGC measurements are the present state-of-the-art for shallow angular surveys, to be superseded only with the Sloan Digital Sky Survey. The APM and EDSGC measurements agree with each other on intermediate-large scales (see Szapudi & Gaztañaga 1998 for detailed comparison of counts in cells measurements of the two surveys). On small scales, however, the EDSGC results are higher than that of the APM. The discrepancy between the APM and EDSGC results can most likely be attributed to differences between the deblending procedures used in the construction of the two catalogs (Szapudi & Gaztañaga 1998). Our measurements, especially for $`S_3`$, appear to lie mostly between the two where they disagree and consistent with them, although perhaps somewhat lower, where they agree. We conclude that the Deeprange counts-in-cells measurements are in agreement with previous estimates in shallow angular surveys, despite the difference in wavelength and the smaller area. Constraints from the three-point correlation function mirror the above results as well – our measurement of $`Q_3=1.57`$ from the $`17I<18`$ subsample is remarkably close to that of Groth and Peebles (1977), $`Q_3=1.29\pm 0.21`$ compiled from the Zwicky, Lick, and Jagellonian samples.

For the collection of magnitude limited cuts, $`S_3`$ and $`S_4`$ are approximately scaling. The $`S_N`$’s appear to decrease with depth, which suggest a small evolution with $`z`$, even though the redshift distribution of each slice is fairly broad. This is in accordance with theories of structure formation where the initial Gaussian fluctuations grow under the influence of gravity (for perturbation theory see e.g., Peebles 1980; Juszkiewicz, Bouchet, & Colombi 1993; Bernardeau 1992; Bernardeau 1994); for $`N`$-body simulations see e.g., Colombi et al. 1995; Baugh, Gaztañaga, & Efstathiou 1995; Szapudi et al. 2000b), and where there is a small bias.

The wide area of the Deeprange survey enables a self-contained study of the time evolution of the bias – there are a statistically sufficient number of low-$`z`$ galaxies contained in the survey that a self-consistent constraint can be derived. Our measurement of the the bias evolution uses a simple model based on the parameterization of the evolution of the two-point correlation function discussed in paper I. Briefly, we parameterize the redshift dependence of the correlation function as

$$\xi (r,z)=(\frac{r}{r_0})^\gamma (1+z)^{3ϵ},$$
(14)
where $`\gamma `$ and $`r_0`$ are the slope and correlation length, and $`ϵ`$ is a phenomenological parameter describing the evolution (e.g. Peebles 1980; Efstathiou et al. 1991; Woods & Fahlman 1997). For $`\mathrm{\Omega }=1`$ and linear evolution, the correlation function at any $`z`$ can then be compared to that at the current epoch through the mapping $`(1+z)^2\xi (r/(1+z),z)`$. Although this mapping is strictly for $`\mathrm{\Omega }=1`$, variations in the mapping due to scenarios with $`\mathrm{\Omega }<1`$ are well within the range of uncertainties due to the exclusion of all details associated with specific galaxy formation processes in this simple model. Under the same mapping the $`S_N`$’s would be invariant. Therefore, we can define the ratio of the mapped correlation function to that at the current epoch as

$$b(z)^2=(1+z)^{\gamma 1ϵ}.$$
(15)
Since our higher order measurements probe galaxies in the highly non-linear regime, it is appropriate to use Equation 10 for biasing, which yields $`S_3(z)=S_3(0)/b(z)^2`$. Figure 3 displays the predicted $`s_3`$, normalized using the measurement from the $`17I<18`$ sample, and the actual measurements at $`\theta =0.04^{}`$ (solid symbols). Redshifts are assigned to the observed data by computing the median $`z`$ for each magnitude slice as predicted by the $`\beta =1.5`$ LF evolution model described in §4.1. The predicted data were computed using $`\gamma =1.75`$, and $`ϵ=1`$. Strictly, $`S_3`$ should be compared at the same comoving scale, but the flat scaling of the $`s_3`$ allows comparison at the same angular scale instead. This model is clearly an oversimplification of the time evolution of the bias, and yet, it captures the trend presented by the data remarkably well. In addition, the time evolution of the bias is constrained, according to Equation 10, to be less than a factor of 2 between the current epoch and $`z0.75`$.

On the other hand, the above observations are in strong contrast with expectations from non-Gaussian scenarios. In these models, primordial non-Gaussianity (e.g., skewness) grows in linear theory. For Gaussian initial conditions, the growth of the skewness is a second order effect. Thus according to theoretical calculations (Fry & Scherrer 1994), and simulations (Colombi 1992), the $`S_N`$’s should have been larger in the past in non-Gaussian models compared to their Gaussian counterparts. As an example, $`S_3`$ is expected to be a factor of 2 larger at $`z1`$ than at $`z0`$ for the typical initial conditions in cosmic string models (Stebbins 1996, private communication; Colombi 1992). The growth of $`S_4`$ with increasing $`z`$ is expected to be even more prominent. These effects would have been detectable in our survey, despite the dilution effects of projection and possible systematic errors at the survey magnitude limit. Formally, the non-Gaussian expectation is about $`78\sigma `$ from our measurement at the highest redshift. Even generously doubling our error bars (which would make the previous naive biasing model based on Gaussian initial conditions a perfectly reasonable fit to the data in Figure 3) would still exclude typical cosmic string non-Gaussian initial conditions at about a $`3.54\sigma `$ level. The Deeprange data, thus, strongly favor Gaussian initial conditions.

While string initial conditions are not favored for many reasons (e.g., Pen, Seljak, & Turok 1997; Albrecht, Battye, & Robinson 1998), our expectation is that the above arguments would hold for a large class of non-Gaussian models. However, it is possible to invoke biasing schemes that are in accord with our observations but that also mask the signature of non-Gaussian initial conditions. For example, strong bias at early times decreases the higher order moments, counteracting the effects of the initial non-Gaussianity. Later the bias naturally decreases, resulting in an increase of the $`S_N`$’s. While perhaps such a scenario is not completely unimaginable physically, e.g. by invoking a strong feedback during galaxy formation, it would require unnatural fine tuning in order to assure that the time evolution of biasing effectively cancels the time evolution of the non-Gaussian initial conditions. Thus, while our results cannot rule out initial non-Gaussianity with high certainty, the most natural explanation is that Gaussian initial conditions of the fluctuation field grew via gravity. Galaxies appear to trace mass quite accurately, and the small evolution of bias predicted by our naive model appears to describe the trends of the data fairly well.

The three-point correlation function appears to be hierarchical down to $`I=22.5`$. Our estimates of the non-hierarchical term, $`Q^{21}`$, are uniformly small in amplitude (see Table 3), although a cubic term $`Q^{111}`$ cannot be excluded with high significance. While the $`\chi ^2`$ analysis is only approximate in nature (the simple error model did not attempt to estimate bin-to-bin cross-correlations), the inclusion of a cubic term does not typically result in a significant change in the goodness of fit. Gravitational instability predicts $`Q^{111}0`$, thus a small cubic term is likely to mean mild bias, as predicted by general bias theory. In the case of a Gaussian field with a completely general bias function $`Q^{111}=Q_3^3`$ is expected under fairly general conditions (Szalay 1988). Since we find $`Q^{111}1`$ and $`Q_31`$, the galaxy density field cannot be a biased version of a Gaussian field but, rather, a mildly biased version of an underlying non-Gaussian field. We emphasize that the non-Gaussianity being referred to here is that which is induced by non-linear gravitational amplification and does not refer to the nature of the initial perturbation spectrum. This interpretation is subject to cosmic variance on $`Q^{111}`$, possible systematic errors, and possible stochasticity of the bias, all of which could contribute to the cubic term.

The best fitting $`S_3=3\times Q_3Q^{11}`$ from the three-point correlation measurements are also displayed in Figure 1 as two horizontal lines: the dotted lines show result of the full fit, while dashed lines display the results of the restricted fit. Despite the fact that the three-point correlation formula extends to larger scales than the counts in cells analysis, the best fitting $`Q^{11}`$ appears to be in excellent agreement with $`S_3`$ obtained from counts in cells, with the exception of the $`19I<20`$ results where there is a factor of two discrepancy. The generally good agreement between the two methods, however, is a further indication of the insignificance of the cubic term $`Q^{111}`$. In Figure 2, the dashed lines show the best fitting $`Q^{11}`$ and the dotted lines display $`Q_3=S_3/3`$ as obtained from the counts in cells analysis. This demonstrates the agreement between the counts in cells analysis and the direct three-point correlation function estimation from a different perspective. The disagreement seen above for the $`19I<20`$ sample is less significant here and, thus, it must stem from the smallest scales. The points display a slight non-hierarchical curvature as well. More accurate measurements from even larger surveys will be required to assess whether this is a significant behavior.

In summary we have measured moments of counts in cells and the three-point correlation function in the Deeprange survey. These constitute the deepest higher order clustering measurements to date. The moments measured on small scales appear to be hierarchical, and the three-point function, extending to larger scales continues this hierarchy. While the cubic term resulting from possible bias could not be excluded with high significance, the hierarchical assumption holds to a good approximation. This argues that gravity is the dominant process in creating galaxy correlations with bias having a detectable but minor role. Qualitatively, models with Gaussian initial conditions and a small amount of biasing, which increases slightly with redshift, are favored. The large area of the Deeprange survey allows us to study the evolution of bias over a relatively broad magnitude range. We find that the bias between $`I`$band selected galaxies and the underlying matter distribution increases slightly with increasing redshift (up to $`z0.8`$) but not by more than a factor of 2.

In Durham, IS was supported by the PPARC rolling grant for Extragalactic Astronomy and Cosmology. The FORCE (FORtran for Cosmic Errors) package can be obtained from its authors, S. Colombi, and IS (http://www.cita.utoronto.ca/szapudi/istvan.html). IS would like to thank Alex Szalay for stimulating discussions.

## Appendix A Fit for the Three-point Correlation Function

The three-point correlation function can be expanded in terms of powers of the two-point correlation function. This is called a Meyer clustering expansion and it is commonly used in the field of statistical physics as well. The second (leading) order term in the expansion corresponds to the hierarchical assumption, which is predicted by gravitational hierarchy (Peebles 1980). The higher order terms, which are determined by the nature of the galaxy–matter biasing (Szalay 1988), represent corrections to the hierarchy that result in non-trivial shape dependencies.

The Meyer clustering expansion, just like any spatial expansion, can be represented by a Feynman-like graphical representation, shown on Figure 4. This figure also illuminates the terminology of connected components, since they correspond to connected components of the graph in the pictorial representation. By definition, three-point function is a connected third order moment: it is the extra probability of a triangle above those predicted by Poisson and Gaussian (i.e. two-point) terms. Consequently, it makes sense to use a connected Meyer expansion for it.

Our aim is to project a full third order connected Meyer clustering expansion of the three dimensional three-point function to its angular analog. Specifically, the spatial three-point correlation function can be expressed as a third order connected expansion of the two-point correlation function, $`\xi `$, as

$`\xi _3=\zeta (1,2,3)`$ $`=`$ $`Q^{11}\left(\xi (1)\xi (2)+\xi (2)\xi (3)+\xi (3)\xi (1)\right)+`$ (A2)
$`Q^{111}\left(\xi (1)\xi (2)\xi (1)\right)+Q^{21}\left(\xi (1)^2\xi (2)+sym\right).`$
This equation is directly translated to the graph on Figure 4. This expansion contains all the possible connected terms; THS included all possible terms, including disconnected ones. Their notation is modified here slightly to avoid confusion with the cumulant correlator notation introduced since THS.

The projected three-point correlation function can be written as

$$z(a,b,c)=\underset{i=11,111,21}{}q^i\omega _i(a,b,c),$$
(A3)
where $`a,b,`$ and $`c`$ are the angles of a triangle on the sky. The $`q^i`$ terms are related to the three-dimensional $`Q^i`$ via $`q^i=R_iQ^i`$. For details see THS. The projection coefficient for the hierarchical term $`R_{11}R_3`$ is identical to that of $`S_3`$. Since we have found that all other terms are consistent with zero (i.e. the hierarchy is a good approximation) we only deproject the hierarchical term. The individual terms of equation A2 project down to two dimensions as (THS)

$`\omega _{11}(a,b,c)`$ $`=\omega (a)\omega (b)+\omega (b)\omega (c)+\omega (c)\omega (a)`$ (A4)
$`\omega _{111}(a,b,c)`$ $`=\left({\displaystyle \frac{\omega (a)\omega (b)\omega (c)}{a+b+c}}\right){\displaystyle \frac{\pi ^2}{H(\gamma )^3}}\left({\displaystyle \frac{180}{\pi }}\right)`$ (A5)
$`\omega _{21}(a,b,c)`$ $`=\left({\displaystyle \frac{\omega (a)^2\omega (b)+\omega (a)\omega (b)^2}{a}}+sym\right){\displaystyle \frac{H(2\gamma )}{H(\gamma )^2}}\left({\displaystyle \frac{180}{\pi }}\right),`$ (A6)
where $`H(\gamma )=_{\mathrm{}}^{\mathrm{}}𝑑x(1+x^2)^{\gamma /2}`$. The second equation is only an approximation, being exact for $`\gamma =2`$. Using the above terms, a minimum $`\chi ^2`$ fit to the angular three-point correlation function was performed with the following choice

$$\chi ^2=\underset{a,b,c}{}\frac{(z(a,b,c)z_{mod}(a,b,c))^2}{\sigma ^2},$$
(A7)
where $`z_{mod}`$ is equation A3. The variance, from Szapudi & Szalay (1998), is

$$\sigma ^2=Var(z)=\frac{6S_3}{S^2\lambda ^3}\frac{6}{DDD}.$$
(A8)
Here $`S_3=\mathrm{\Phi }(a,b,c)^2`$, a configuration integral over the definition of the angular bin, described by the function $`\mathrm{\Phi }(a,b,c)=1`$ when a triangle is in the bin, and zero otherwise.

The simplification arises since $`S_3=S=\mathrm{\Phi }`$, and $`S\lambda ^3DDD`$ for our choice of the characteristic function describing the bin ($`\mathrm{\Phi }=\mathrm{\Phi }^2`$) used in the three-point estimator (see Szapudi & Szalay (1998) for details). $`DDD`$ is the number of data triplets in the bin.

### A.1 Error Estimation

The variance defined in the above equation was derived for a Poisson distribution and therefore accounts only for discreteness effects. Two additional error contributions, edge and finite volume effects (Szapudi & Colombi 1996; Szapudi, Colombi, & Bernardeau 1999), arising from the uneven weights given to data points and from the fluctuations of the universe on scales larger than the survey size, respectively, are not accounted for. However, edge effects are expected to be small for these edge corrected estimators, and it can be shown that at least partial correction for finite volume effects are contained in $`DDD`$. The $`\chi ^2`$’s suggest that this ansatz for the variance is within factor of 2 of the truth. This is the limit of this simple model ignoring cross correlations of different bins. Without complicated treatment of the full correlation matrix, it would not make sense to improve the above formula with the inclusion of an additional phenomenological term for finite volume effects to tune $`\chi ^2`$ to the number of degrees of freedom. While the results should not sensitively depend on the exact choice of the error model, it should be kept in mind that our errors for $`z(a,b,c)`$ could be off by as much as a factor of 2. The most likely sense of the offset is that our errors may be overestimated, as judged from the $`\chi ^2`$’s.