Estimating the galaxy two-point correlation function using a split random catalog

TL;DR

This paper proposes a split random catalog method to efficiently estimate the galaxy two-point correlation function, significantly reducing computation time while maintaining accuracy, crucial for large upcoming galaxy surveys.

Contribution

Introducing a split random catalog approach that optimizes the calculation of the galaxy two-point correlation function by reducing computational costs without sacrificing precision.

Findings

Splitting the random catalog into subcatalogs improves efficiency.

Excluding cross-subcatalog pairs maintains estimator accuracy.

Reduces computation time by over tenfold for large catalogs.

Abstract

The two-point correlation function of the galaxy distribution is a key cosmological observable that allows us to constrain the dynamical and geometrical state of our Universe. To measure the correlation function we need to know both the galaxy positions and the expected galaxy density field. The expected field is commonly specified using a Monte-Carlo sampling of the volume covered by the survey and, to minimize additional sampling errors, this random catalog has to be much larger than the data catalog. Correlation function estimators compare data-data pair counts to data-random and random-random pair counts, where random-random pairs usually dominate the computational cost. Future redshift surveys will deliver spectroscopic catalogs of tens of millions of galaxies. Given the large number of random objects required to guarantee sub-percent accuracy, it is of paramount importance to improve the efficiency of the algorithm without degrading its precision. We show both analytically and numerically that splitting the random catalog into a number of subcatalogs of the same size as the data catalog when calculating random-random pairs, and excluding pairs across different subcatalogs provides the optimal error at fixed computational cost. For a random catalog fifty times larger than the data catalog, this reduces the computation time by a factor of more than ten without affecting estimator variance or bias.

Tables

Table 1. Table 1: Mean computation time over the 300 runs and the mean variance over four different ranges of r 𝑟 r bins (given in units of h − 1   Mpc times h 1 megaparsec {\mathrm{\mathnormal{h}}}^{-1}\text{\,}\mathrm{Mpc} ) for each method. The first three are standard LS. The variance cannot be measured accurately enough from 300 realizations to correctly show all the differences between methods. Thus the table shows some apparent improvements (going from M r = 50 subscript 𝑀 r 50 M_{\rm r}=50 to M r = 25 subscript 𝑀 r 25 M_{\rm r}=25 and from M s = 4 subscript 𝑀 s 4 M_{\rm s}=4 to M s = 16 subscript 𝑀 s 16 M_{\rm s}=16 ), which are not to be taken as real. See Fig.  6 for the fifth vs. second column with error bars.

method	time	mean variance
	[s]	$0$–$2$	$5$–$15$	$80$–$120$	$150$–$200$
		[$\times {10}^{-2}$]	[$\times {10}^{-5}$]	[$\times {10}^{-6}$]	[$\times {10}^{-7}$]
$M_r=50$	7889	1.01	4.42	1.23	7.01
$M_r=25$	2306	1.01	4.39	1.23	7.04
$M_r=1$	16.5	5.96	5.53	1.46	8.05
$d=0.5$	2239	1.01	4.41	1.25	7.11
$d=0.25$	812	1.02	4.41	1.25	7.42
$d=0.14$	487	1.08	4.41	1.28	7.72
$M_s=4$	1854	1.01	4.42	1.23	7.02
$M_s=16$	763	1.00	4.42	1.23	7.01
$M_s=50$	593	1.09	4.42	1.23	7.02

Equations

$${\hat{\xi }}_{\mathrm{LS}}(\mathbf{r}):=\frac{N_{\mathrm{r}}(N_{\mathrm{r}}-1)}{N_{\mathrm{d}}(N_{\mathrm{d}}-1)}\frac{DD(\mathbf{r})}{RR(\mathbf{r})}-\frac{N_{\mathrm{r}}-1}{N_{\mathrm{d}}}\frac{DR(\mathbf{r})}{RR(\mathbf{r})}+1,$$

$$DD(\mathbf{r})$$

$$DR(\mathbf{r})$$

$$RR(\mathbf{r})$$

$$⟨DD(\mathbf{r})⟩$$

$$⟨DR(\mathbf{r})⟩$$

$$⟨RR(\mathbf{r})⟩$$

$$G^{\mathrm{p}}(\mathbf{r}):=\frac{2}{K^2}\sum _{i<j}^K{\Theta }^{ij}(\mathbf{r})$$

$$⟨{\hat{\xi }}_{\mathrm{LS}}⟩$$

$$\begin{array}{rl}& \text{mbox}Cov\left[{\hat{\xi }}_{\mathrm{LS}}({\mathbf{r}}_1),{\hat{\xi }}_{\mathrm{LS}}({\mathbf{r}}_2)\right]\equiv ⟨{\hat{\xi }}_{\mathrm{LS}}({\mathbf{r}}_1){\hat{\xi }}_{\mathrm{LS}}({\mathbf{r}}_2)⟩-⟨{\hat{\xi }}_{\mathrm{LS}}({\mathbf{r}}_1)⟩⟨{\hat{\xi }}_{\mathrm{LS}}({\mathbf{r}}_2)⟩\\ & \eqsim (1+{\xi }_1)(1+{\xi }_2)⟨{\alpha }_1{\alpha }_2⟩+4⟨{\beta }_1{\beta }_2⟩+(1-{\xi }_1)(1-{\xi }_2)⟨{\gamma }_1{\gamma }_2⟩\\ & -2(1+{\xi }_1)⟨{\alpha }_1{\beta }_2⟩-2(1+{\xi }_2)⟨{\beta }_1{\alpha }_2⟩\\ & -2(1-{\xi }_1)⟨{\gamma }_1{\beta }_2⟩-2(1-{\xi }_2)⟨{\beta }_1{\gamma }_2⟩.\end{array}$$

$$\begin{array}{rl}& \Delta \text{mbox}Cov\left[{\hat{\xi }}_{\mathrm{LS}}({\mathbf{r}}_1),{\hat{\xi }}_{\mathrm{LS}}({\mathbf{r}}_2)\right]\eqsim 4\Delta ⟨{\beta }_1{\beta }_2⟩+(1-{\xi }_1)(1-{\xi }_2)⟨{\gamma }_1{\gamma }_2⟩\\ & -2(1-{\xi }_1)⟨{\gamma }_1{\beta }_2⟩-2(1-{\xi }_2)⟨{\beta }_1{\gamma }_2⟩.\end{array}$$

$$⟨D{D}_{1}D{D}_2⟩=⟨D{D}_1⟩⟨D{D}_2⟩\left(1+⟨{\alpha }_1{\alpha }_2⟩\right),$$

$$⟨{\alpha }_1{\alpha }_2⟩$$

$$⟨{\beta }_1{\beta }_2⟩$$

$$⟨{\gamma }_1{\gamma }_2⟩$$

$$⟨{\alpha }_1{\beta }_2⟩$$

$$⟨{\beta }_1{\gamma }_2⟩$$

$$⟨{\alpha }_1{\gamma }_2⟩$$

$$G_{12}^{\mathrm{t}}:=G^{\mathrm{t}}({\mathbf{r}}_1,{\mathbf{r}}_2):=\frac{1}{K^3}\sum _{ijk}^{\ast }{\Theta }_1^{ik}{\Theta }_2^{jk}$$

$$⟨n_i{n}_j{n}_k⟩\propto 1+{\xi }_{ij}+{\xi }_{jk}+{\xi }_{ik},$$

$$⟨n_i{n}_j{n}_k{n}_l⟩\propto 1+{\xi }_{ij}+{\xi }_{jk}+{\xi }_{ik}+{\xi }_{il}+{\xi }_{jl}+{\xi }_{kl}+{\xi }_{ij}{\xi }_{kl}+{\xi }_{ik}{\xi }_{jl}+{\xi }_{il}{\xi }_{jk}.$$

$$(1+{\xi }_1)(1+{\xi }_2)⟨{\alpha }_1{\alpha }_2⟩\simeq \frac{4}{N_{\mathrm{d}}}(1+{\xi }_1)(1+{\xi }_2)\left[\frac{G_{12}^{\mathrm{t}}}{G_1^{\mathrm{p}}{G}_2^{\mathrm{p}}}-1\right]$$

$$+\frac{2(1+{\xi }_1)}{N_{\mathrm{d}}(N_{\mathrm{d}}-1)}\left[\frac{{\delta }_{12}}{G_1^{\mathrm{p}}}-2(1+{\xi }_2)\frac{G_{12}^{\mathrm{t}}}{G_1^{\mathrm{p}}{G}_2^{\mathrm{p}}}+(1+{\xi }_2)\right]$$

$$+\frac{4(N_{\mathrm{d}}-2)}{N_{\mathrm{d}}(N_{\mathrm{d}}-1)}({\xi }_{12}-{\xi }_1{\xi }_2)\frac{G_{12}^{\mathrm{t}}}{G_1^{\mathrm{p}}{G}_2^{\mathrm{p}}},$$

$$\displaystyle\langle\beta_{1}\beta_{2}\rangle\ \simeq\ \frac{1}{N_{\rm d}N_{\rm r}}\Biggl{\{}N_{\rm r}\left[\frac{G^{\rm t}_{12}}{G^{\rm p}_{1}G^{\rm p}_{2}}-1\right]+N_{\rm d}\left[\frac{G^{\rm t}_{12}}{G^{\rm p}_{1}G^{\rm p}_{2}}-1\right]$$

$$\displaystyle\qquad+1-\frac{2G^{\rm t}_{12}}{G^{\rm p}_{1}G^{\rm p}_{2}}+\frac{\delta_{12}}{G^{\rm p}_{1}}\Biggr{\}}+\frac{N_{\rm d}-1}{N_{\rm d}N_{\rm r}}\xi_{12}\frac{G^{\rm t}_{12}}{G^{\rm p}_{1}G^{\rm p}_{2}}\,,$$

$$⟨{\gamma }_1{\gamma }_2⟩=\frac{2}{N_{\mathrm{r}}(N_{\mathrm{r}}-1)}\left\{2(N_{\mathrm{r}}-2)\left[\frac{G_{12}^{\mathrm{t}}}{G_1^{\mathrm{p}}{G}_2^{\mathrm{p}}}-1\right]+\frac{{\delta }_{12}}{G_1^{\mathrm{p}}}-1\right\},$$

$$⟨{\alpha }_1{\beta }_2⟩\simeq \frac{2}{N_{\mathrm{d}}}\left[\frac{G_{12}^{\mathrm{t}}}{G_1^{\mathrm{p}}{G}_2^{\mathrm{p}}}-1\right],$$

$$⟨{\beta }_1{\gamma }_2⟩=\frac{2}{N_{\mathrm{r}}}\left[\frac{G_{12}^{\mathrm{t}}}{G_1^{\mathrm{p}}{G}_2^{\mathrm{p}}}-1\right],$$

$$⟨{\alpha }_1{\gamma }_2⟩=0,$$

$$t_{12}$$

$$t_{12}^{\mathrm{r}}$$

$$p_{12}$$

$$p_{12}^{\mathrm{c}}$$

$$p_{12}^{\mathrm{r}}$$

$$q_{12}$$

$$q_{12}^{\mathrm{r}}$$

$$(1+{\xi }_1)(1+{\xi }_2)⟨{\alpha }_1{\alpha }_2⟩\simeq (1+{\xi }_1)(1+{\xi }_2)(4t_{12}+p_{12})$$

$$+4\frac{N_{\mathrm{d}}-2}{(N_{\mathrm{d}}-1)}({\xi }_{12}-{\xi }_1{\xi }_2)q_{12},$$

$$⟨{\beta }_1{\beta }_2⟩\simeq t_{12}+t_{12}^{\mathrm{r}}+p_{12}^{\mathrm{c}}+\frac{N_{\mathrm{d}}-1}{N_{\mathrm{d}}}{\xi }_{12}{q}_{12}^{\mathrm{r}},$$

$$⟨{\gamma }_1{\gamma }_2⟩=4t_{12}^{\mathrm{r}}+p_{12}^{\mathrm{r}},$$

$$⟨{\alpha }_1{\beta }_2⟩\simeq 2t_{12},⟨{\beta }_1{\gamma }_2⟩=2t_{12}^{\mathrm{r}},\text{mbox}and⟨{\alpha }_1{\gamma }_2⟩=0.$$

$$\Delta ⟨{\beta }_1{\beta }_2⟩\simeq t_{12}^{\mathrm{r}}+p_{12}^{\mathrm{c}}+\frac{N_{\mathrm{d}}-1}{N_{\mathrm{d}}}{\xi }_{12}{q}_{12}^{\mathrm{r}}.$$

$$⟨{\hat{\xi }}_{\mathrm{LS}}⟩=\xi +\left(\xi -1\right)(4t_{\mathrm{r}}+p_{\mathrm{r}})+4t_{\mathrm{r}}.$$

Full text

1 Introduction

The spatial distribution of luminous matter in the Universe is a key diagnostic for studying cosmological models and the physical processes involved in the assembly of structure. In particular, light from galaxies is a robust tracer of the overall matter distribution, whose statistical properties can be predicted by cosmological models. Two-point correlation statistics are very effective tools for compressing the cosmological information encoded in the spatial distribution of the mass in the Universe. In particular, the two-point correlation function in configuration space has emerged as one of the most popular cosmological probes. Its success stems from the presence of characterized features that can be identified, measured, and effectively compared to theoretical models to extract clean cosmological information.

One such feature is baryon acoustic oscillations (BAOs), which imprint a characteristic scale in the two-point correlation function that can be used as a standard ruler. After the first detection in the two-point correlation function of SDSS DR3 and 2dFGRS galaxy catalogs (Eisenstein et al. 2005; Cole et al. 2005), the BAO signal was identified, with different degrees of statistical significance, and has since been used to constrain the expansion history of the Universe in many spectroscopic galaxy samples (see e.g., Percival et al. (2010) Blake et al. (2011), Beutler et al. (2011), Anderson et al. (2012), Anderson et al. (2014), Ross et al. (2015), Alam et al. (2017), Ross et al. (2017), Vargas-Magaña et al. (2018), Bautista et al. (2018), and Ata et al. (2018)). Several of these studies did not focus on the BAO feature only but also analyzed the anisotropies in the two-point correlation function induced by the peculiar velocities (Kaiser 1987), the so-called redshift space distortions (RSD), and by assigning cosmology-dependent distances to the observed redshifts (the Alcock & Paczyński (1979) test). For RSD analyses, see also, for example, Peacock et al. (2001), Guzzo et al. (2008), Beutler et al. (2012), Reid et al. (2012), de la Torre et al. (2017), Pezzotta et al. (2017), Zarrouk et al. (2018), Hou et al. (2018), and Ruggeri et al. (2019).

Methods to estimate the galaxy two-point correlation function (2PCF) $\xi({\bf r})$ from survey data are based on its definition as the excess probability of finding a galaxy pair. One counts from the data ($D$) catalog the number $DD({\bf r})$ of pairs of galaxies with separation ${\bf x}_{2}-{\bf x}_{1}\in{\bf r}$, where ${\bf r}$ is a bin of separation vectors, and compares it to the number of pairs $RR({\bf r})$ in a corresponding randomly generated ($R$) catalog and to the number of data-random pairs $DR({\bf r})$. The bin may be a 1D ($r\pm{\textstyle\frac{1}{2}}\Delta r$), 2D, or a 3D bin. In the 1D case, $r$ is the length of the separation vector and $\Delta r$ is the width of the bin. From here on, ‘separation ${\bf r}$’ indicates that the separation falls in this bin.

Several estimators of the 2PCF have been proposed by Hewett (1982), Davis & Peebles (1983), Hamilton (1993), and Landy & Szalay (1993), building on the original Peebles & Hauser (1974) proposal. These correspond to different combinations of the $DD$, $DR$, and $RR$ counts to obtain a 2PCF estimate $\hat{\xi}({\bf r})$; see Kerscher (1999) and Kerscher (2000) for more estimators. The Landy–Szalay (Landy & Szalay 1993) estimator

[TABLE]

(we call this method ‘standard LS’ in the following) is the most commonly used, since it provides the minimum variance when $|\xi|\ll 1$ and is unbiased in the limit $N_{\rm r}\rightarrow\infty$. Here $N_{\rm d}$ is the size (number of objects) of the data catalog and $N_{\rm r}$ is the size of the random catalog. We define $M_{\rm r}:=N_{\rm r}/N_{\rm d}$. To minimize random error from the random catalog, $M_{\rm r}\gg 1$ should be used (for a different approach, see Demina et al. (2018)).

One is usually interested in $\xi({\bf r})$ only up to some $r_{\mathrm{max}}\ll L_{\mathrm{max}}$ (the maximum separation in the survey), and therefore pairs with larger separations can be skipped. Efficient implementations of the LS estimator involve pre-ordering of the catalogs through kd-tree, chain-mesh, or other algorithms (e.g., Moore et al. 2000; Alonso 2012; Jarvis 2015; Marulli et al. 2016) to facilitate this. The computational cost is then roughly proportional to the actual number of pairs with separation $r\leq r_{\mathrm{max}}$.

The correlation function is small for large separations, and in cosmological surveys $r_{\mathrm{max}}$ is large enough so that for most pairs $|\xi({\bf r})|\ll 1$. The fraction $f$ of $DD$ pairs with $r\leq r_{\mathrm{max}}$ is therefore not very different from the fraction of $DR$ or $RR$ pairs with $r\leq r_{\mathrm{max}}$. The computational cost is dominated by the part proportional to the total number of pairs needed, ${\textstyle\frac{1}{2}}fN_{\rm d}(N_{\rm d}\!-\!1)+fN_{\rm d}N_{\rm r}+{\textstyle\frac{1}{2}}fN_{\rm r}(N_{\rm r}\!-\!1)\approx{\textstyle\frac{1}{2}}fN_{\rm d}^{2}(1+2M_{\rm r}+M_{\rm r}^{2})$, which in turn is dominated by the $RR$ pairs as $M_{\rm r}\gg 1$. The smaller number of $DR$ pairs contribute much more to the error of the estimate than the large number of $RR$ pairs, whereas the cost is dominated by $RR$. Thus, a significant saving of computation time with an insignificant loss of accuracy may be achieved by counting only a subset of $RR$ pairs, while still counting the full set (up to $r_{\mathrm{max}}$) of $DR$ pairs.

A good way to achieve this is to use many small (i.e, low-density) $R$ catalogs instead of one large (high-density) catalog (Landy & Szalay 1993; Wall & Jenkins 2012; Slepian & Eisenstein 2015), or, equivalently, to split an already generated large $R$ catalog into $M_{\rm s}$ small ones for the calculation of $RR$ pairs while using the full $R$ catalog for the $DR$ counts. This method has been used by some practitioners (e.g., Zehavi et al. (2011)111I. Zehavi, private communication), but this is usually not documented in the literature. One might also consider obtaining a similar cost saving by diluting (subsampling) the $R$ catalog for $RR$ counts, but, as we show below, this is not a good idea. We refer to these two cost-saving methods as ‘split’ and ‘dilution’.

In this work we theoretically** **derive the additional covariance and bias due to the size and treatment of the $R$ catalog; test these predictions numerically with mock catalogs representative of next-generation datasets, such as the spectroscopic galaxy samples that will be obtained by the future Euclid satellite mission (Laureijs et al. 2011); and show that the ‘split’ method, while reducing the computational cost by a large factor, retains the advantages of the LS estimator.

We follow the approach of Landy & Szalay (1993; hereafter, LS93), but generalize it in a number of ways: In particular, since we focus on the effect of the random catalog, we do not work in the limit $M_{\rm r}\rightarrow\infty$. Also, we calculate covariances, not just variances, and make fewer approximations (see Sect. 2.2).

The layout of the paper is as follows. In Sect. 2 we derive theoretical results for bias and covariance. In Sect. 3 we focus on the split LS estimator and its optimization. In Sect. 4 we test the different estimators with mock catalogs. Finally, we discuss the results and present our conclusions in Sect. 5.

2 Theoretical results: bias and covariance

2.1 General derivation

We follow the derivation and notations in LS93 but extend to the case that includes random counts covariance. We consider the survey volume as divided into $K$ microcells (very small subvolumes) and work in the limit $K\rightarrow\infty$, which means that no two objects will ever be located within the same microcell.

Here, $\alpha$, $\beta$, and $\gamma$ represent the relative deviation of the $DD({\bf r})$, $DR({\bf r})$, and $RR({\bf r})$ counts from their expectation values (mean values over an infinite number of independent realizations):

[TABLE]

By definition $\langle\alpha\rangle=\langle\beta\rangle=\langle\gamma\rangle=0$. The factors $\alpha$, $\beta$, and $\gamma$ represent fluctuations in the pair counts, which arise as a result of a Poisson process. As long as the mean pair counts per bin are large ($\gg 1$) the relative fluctuations will be small. We calculate up to second order in $\alpha$, $\beta$, and $\gamma$, and ignore the higher-order terms (in the limit $M_{r}\rightarrow\infty$, $\gamma\rightarrow 0$, so LS93 set $\gamma=0$ at this point).

The expectation values for the pair counts are:

[TABLE]

where $\xi({\bf r})$ is the correlation function normalized to the actual number density of galaxies in the survey and

[TABLE]

is the fraction of microcell pairs with separation ${\bf r}$. Here $\Theta^{ij}({\bf r}):=1$ if ${\bf x}_{i}-{\bf x}_{j}$ falls in the ${\bf r}$-bin, otherwise it is equal to zero.

The expectation value of the LS estimator (1) is

[TABLE]

A finite $R$ catalog thus introduces a (small) bias. (In LS93, $\gamma=0$, so the estimator is unbiased in this limit.) This expression is calculated to second order in $\alpha$, $\beta$, and $\gamma$ (we denote equality to second order by ‘$\eqsim$’). Calculation to higher order is beyond the scope of this work. Since data and random catalogs are independent, $\langle\alpha\gamma\rangle=0$.

We introduce shorthand notations $\langle\alpha_{1}\alpha_{2}\rangle$ for $\langle\alpha({\bf r}_{1})\alpha({\bf r}_{2})\rangle$, $\langle DD_{1}DD_{2}\rangle$ for $\langle DD({\bf r}_{1})DD({\bf r}_{2})\rangle$, and similarly for other terms.

For the covariance we get

[TABLE]

Terms with $\gamma$ represent additional variance due to finite $N_{\rm r}$, and are new compared to those of LS93. Also, $\langle\beta_{1}\beta_{2}\rangle$ collects an additional contribution, which we denote by $\Delta\langle\beta_{1}\beta_{2}\rangle$, from variations in the random field (see Sect. 2.3). The cross terms $\alpha_{1}\beta_{2}$ and $\alpha_{2}\beta_{1}$, instead depend linearly on the random field, and average to the $N_{\rm r}\rightarrow\infty$ result. The additional contribution due to finite $N_{\rm r}$ is thus

[TABLE]

From (2),

[TABLE]

and so on, so that the covariances of the deviations are obtained from

[TABLE]

2.2 Quadruplets, triplets, and approximations

We use

[TABLE]

to denote the fraction of ordered microcell triplets, where ${\bf x}_{i}-{\bf x}_{k}\in{\bf r}_{1}$ and ${\bf x}_{j}-{\bf x}_{k}\in{\bf r}_{2}$. The notation $\sum^{\ast}$ means that only terms where all indices (microcells) are different are included. Here $G^{\rm t}_{12}$ is of the same magnitude as $G^{\rm p}_{1}G^{\rm p}_{2}$ but is larger.

Appendix A gives examples of how the $\langle DD_{1}\ DD_{2}\rangle$ and so on in (9) are calculated. These covariances involve expectation values $\langle n_{i}n_{j}n_{l}n_{k}\rangle$, where $n_{i}$ is the number of objects (0 or 1) in microcell $i$ and so on, and only cases where the four microcells are separated pairwise by ${\bf r}_{1}$ and ${\bf r}_{2}$ are included. If all four microcells $i$, $j$, $k$, and $l$ are different, we call this case a quadruplet; it consists of two pairs with separations ${\bf r}_{1}$ and ${\bf r}_{2}$. If two of the indices, that is, microcells, are equal, we have a triplet with a center cell (the equal indices) and two end cells separated from the center by ${\bf r}_{1}$ and ${\bf r}_{2}$.

We make the following three approximations:

For microcell quadruplets, the correlations between unconnected cells are approximated by zero on average. 2. 2.

Three-point correlations vanish. 3. 3.

The part of four-point correlations that does not arise from the two-point correlations vanishes.

With approximations (2) and (3), we have for the expectation value of a galaxy triplet

[TABLE]

where $\xi_{ij}:=\xi({\bf x}_{j}-{\bf x}_{i})$, and for a quadruplet

[TABLE]

We use ‘$\simeq$’ to denote results based on these three approximations. Approximation (1) is good as long as the survey size is large compared to $r_{\mathrm{max}}$. It allows us to drop terms other than $1+\xi_{ij}+\xi_{kl}+\xi_{ij}\xi_{kl}$ in (12). Approximations (2) and (3) hold for Gaussian density fluctuations, but in the realistic cosmological situation they are not good: the presence of the higher-order correlations makes the estimation of the covariance of $\xi({\bf r})$ estimators a difficult problem. However, this difficulty applies only to the contribution of the data to the covariance, that is, to the part that does not depend on the size and treatment of the random catalog. The key point in this work is that while our theoretical result for the total covariance does not hold in a realistic situation (it is an underestimate), our results for the difference in estimator covariance due to different treatments of the random catalog hold well.

In addition to working in the limit $N_{\rm r}\rightarrow\infty$ ($\gamma=0$), LS93 considered only 1D bins and the case where $r_{1}=r_{2}\equiv r$ (i.e., variances, not covariances) and made also a fourth approximation: for triplets (which in this case have legs of equal length) they approximated the correlation between the end cells (whose separation in this case varies between [math] and $2r$) by $\xi(r)$. We use $\xi_{12}$ to denote the mean value of the correlation between triplet end cells (separated from the triplet center by ${\bf r}_{1}$ and ${\bf r}_{2}$). (For our plots in Sect. 4 we make a similar approximation of $\xi_{12}$ as Landy&Szalay, see Sect. 4.2.) Also, LS93 only calculated to first order in $\xi$, whereas we do not make this approximation. Bernstein (1994) also considered covariances, and included the effect of three-point and four-point correlations, but worked in the limit $N_{\rm r}\rightarrow\infty$ ($\gamma=0$).

2.3 Poisson, edge, and $q$ terms

After calculating all the $\langle DD_{1}\ DD_{2}\rangle$ and so on (see Appendix A), (9) becomes

[TABLE]

for the standard LS estimator.

Following the definition of $t$ and $p$ in LS93, we define

[TABLE]

For their diagonals (${\bf r}_{1}={\bf r}_{2}$), we write $t$, $t_{\rm r}$, $p$, $p_{\rm c}$, $p_{\rm r}$, $q$, and $q_{\rm r}$. Thus, $t\equiv t_{11}\equiv t_{22}$, $t_{\rm r}\equiv t^{\rm r}_{11}\equiv t^{\rm r}_{22}$ and so on. (We use superscripts for the matrices, e.g., $t^{\rm r}({\bf r}_{1},{\bf r}_{2}),$ and subscripts for their diagonals, e.g., $t_{\rm r}({\bf r})$.)

Using these definitions, (13) becomes

[TABLE]

The new part in $\langle\beta_{1}\beta_{2}\rangle$ due to finite size of the random catalog is

[TABLE]

Thus only $\langle\alpha_{1}\alpha_{2}\rangle$ and $\langle\beta_{1}\beta_{2}\rangle$ are affected by $\xi({\bf r})$ (in our approximation its effect cancels in $\langle\alpha_{1}\beta_{2}\rangle$). The results for $\langle\gamma_{1}\gamma_{2}\rangle$, $\langle\beta_{1}\gamma_{2}\rangle$, and $\langle\alpha_{1}\gamma_{2}\rangle$ are exact. The result for $\langle\alpha_{1}\alpha_{2}\rangle$ involves all three approximations mentioned above, $\langle\alpha_{1}\beta_{2}\rangle$ involves approximations (1) and (2), and $\langle\beta_{1}\beta_{2}\rangle$ involves approximation (1).

We refer to $p$, $p^{\rm c}$, and $p^{\rm r}$ as ‘Poisson’ terms and $t$ and $t^{\rm r}$ as ‘edge’ terms (the difference between $G^{\rm t}_{12}$ and $G^{\rm p}_{1}G^{\rm p}_{2}$ is due to edge effects). While the Poisson terms are strongly diagonal dominated, the edge terms are not. Since $N_{\rm d}t_{12}=N_{\rm r}t^{\rm r}_{12}\ll 1$, the $q$ terms are much larger than the edge terms, but they get multiplied by $\xi_{12}-\xi_{1}\xi_{2}$ or $\xi_{12}$. In the limit $N_{\rm r}\rightarrow\infty$: $\langle\beta_{1}\gamma_{2}\rangle\rightarrow 0$, $\langle\gamma_{1}\gamma_{2}\rangle\rightarrow 0$, $\langle\beta_{1}\beta_{2}\rangle\rightarrow t_{12}$; $\langle\alpha_{1}\alpha_{2}\rangle,$ and also $\langle\alpha_{1}\beta_{2}\rangle$ are unaffected.

We see that $DD$-$DR$ and $DR$-$RR$ correlations arise from edge effects. If we increase the density of data or random objects, the Poisson terms decrease as $N^{-2}$ but the edge terms decrease only as $N^{-1}$ so the edge effects are more important for a higher density of objects.

Doubling the bin size (combining neighboring bins) doubles $G^{\rm p}({\bf r})$ but makes $G^{\rm t}({\bf r}_{1},{\bf r}_{2})$ four times as large, since triplets where one leg was in one of the original smaller bins and the other leg was in the other bin are now also included. Thus, the ratio $G^{\rm t}_{12}/(G^{\rm p}_{1}G^{\rm p}_{2})$ and $t$ are not affected, but the dominant term in $p$, $1/(1+\xi)G^{\rm p}$ is halved. Edge effects are thus more important for larger bins.

2.4 Results for the standard Landy–Szalay estimator

Inserting the results for $\langle\alpha_{1}\alpha_{2}\rangle$ and so on into Eqs. (5) and (6), we get that the expectation value of the standard LS estimator (1) is

[TABLE]

This holds also for large $\xi$ and in the presence of three-point and four-point correlations. A finite $R$ catalog thus introduces a bias $\left(\xi-1\right)(4t_{\rm r}+p_{\rm r})+4t_{\rm r}=-p_{\rm r}+(4t_{\rm r}+p_{\rm r})\xi$; the edge ($t_{\rm r}$) part of the bias cancels in the $\xi\rightarrow 0$ limit.

For the covariance we get

[TABLE]

Because of the approximations made, this result for the covariance does not apply to the realistic cosmological case; not even for large separations $r$, where $\xi$ is small, since large correlations at small $r$ increase the covariance also at large $r$. However, this concerns only $\langle\alpha_{1}\alpha_{2}\rangle$ and $\langle\alpha_{1}\beta_{2}\rangle$. Our focus here is on the additional covariance due to the size and handling of the random catalog, which for standard LS is

[TABLE]

To zeroth order in $\xi$ the covariance is given by the Poisson terms and the edge terms cancel to first order in $\xi$. This is the property for which the standard LS estimator was designed. To first order in $\xi$, the $q$ terms contribute. This $q$ contribution involves the triplet correlation $\xi_{12}$, which, depending on the form of $\xi({\bf r})$, may be larger than $\xi_{1}$ or $\xi_{2}$.

If we try to save cost by using a diluted random catalog with $N^{\prime}_{\rm r}\ll N_{\rm r}$ for $RR$ pairs, $\langle\gamma_{1}\gamma_{2}\rangle$ is replaced by $\langle\gamma^{\prime}_{1}\gamma^{\prime}_{2}\rangle=4t^{{\rm r}^{\prime}}_{12}+p^{{\rm r}^{\prime\prime}}_{12}$ with $N^{\prime}_{\rm r}$ in place of $N_{\rm r}$, but $\langle\beta_{1}\gamma_{2}^{\prime}\rangle=\langle\beta_{1}\gamma_{2}\rangle$ and $\langle\beta_{1}\beta_{2}\rangle$ are unaffected, so that the edge terms involving randoms no longer cancel. In Sect. 4 we see that this is a large effect. Therefore, one should not use dilution.

3 Split random catalog

3.1 Bias and covariance for the split method

In the split method one has $M_{\rm s}$ independent smaller $R^{\mu}$ catalogs of size $N^{\prime}_{\rm r}$ instead of one large random catalog $R$. Their union, $R,$ has a size of $N_{\rm r}=M_{\rm s}N^{\prime}_{\rm r}$. The pair counts $DR({\bf r})$ and $RR^{\prime}({\bf r})$ are calculated as

[TABLE]

that is, pairs across different $R^{\mu}$ catalogs are not included in $RR^{\prime}$. The total number of pairs in $RR^{\prime}$ is ${\textstyle\frac{1}{2}}M_{\rm s}N^{\prime}_{\rm r}(N^{\prime}_{\rm r}-1)={\textstyle\frac{1}{2}}N_{\rm r}(N^{\prime}_{\rm r}-1)$. Here, $DR$ is equal to its value in standard LS.

The split Landy–Szalay estimator is

[TABLE]

Compared to standard LS, $\langle\alpha_{1}\alpha_{2}\rangle$, $\langle\beta_{1}\beta_{2}\rangle$, and $\langle\alpha_{1}\beta_{2}\rangle$ are unaffected. We construct $\langle RR^{\prime}\rangle$, $\langle RR^{\prime}\cdot RR^{\prime}\rangle$, and $\langle RR^{\prime}\cdot DR\rangle$ from the standard LS results, bearing in mind that the random catalog is a union of independent catalogs, arriving at

[TABLE]

where

[TABLE]

The first is the same as in standard LS and dilution, but the second differs both from standard LS and from dilution, since it involves both $N_{\rm r}$ and $N^{\prime}_{\rm r}$.

For the expectation value we get

[TABLE]

so that the bias is $(\xi-1)(4t_{\rm r}+p^{\prime}_{\rm r})+4t_{\rm r}\ =\ -p^{\prime}_{\rm r}+(4t_{\rm r}+p^{\prime}_{\rm r})\xi$. In the limit $\xi\rightarrow 0$ the edge part cancels, leaving only the Poisson term.

The covariance is

[TABLE]

The change in the covariance compared to the standard LS method is

[TABLE]

which again applies in the realistic cosmological situation. Our main result is that in the split method the edge effects cancel and the bias and covariance are the same as for standard LS, except that the Poisson term $p^{\rm r}$ from $RR$ is replaced with the larger $p^{{\rm r}^{\prime}}$.

3.2 Optimizing computational cost and variance of the split method

The bias is small compared to variance in our application (see Fig. 1 for the theoretical result and Fig. 5 for an attempted bias measurement), and therefore we focus on variance as the figure of merit. The computational cost should be roughly proportional to

[TABLE]

and the additional variance due to finite $R$ catalog in the $\xi\rightarrow 0$ limit becomes

[TABLE]

Here, $N_{\rm d}$ and $p$ are fixed by the survey and the requested ${\bf r}$ binning, but we can vary $M_{\rm r}$ and $M_{\rm s}$ in the search for the optimal computational method. In the above we defined the ‘cost’ and ‘variance’ factors $c$ and $v$.

We may ask two questions:

For a fixed level of variance $v$, which combination of $M_{\rm r}$ and $M_{\rm s}$ minimizes computational cost $c$? 2. 2.

For a fixed computational cost $c$, which combination of $M_{\rm r}$ and $M_{\rm s}$ minimizes the variance $v$?

The answer to both questions is (Slepian & Eisenstein 2015)

[TABLE]

Thus, the optimal version of the split method is the natural one where $N^{\prime}_{\rm r}=N_{\rm d}$. In this case the additional variance in the $\xi\rightarrow 0$ limit becomes

[TABLE]

and the computational cost factor $N_{\rm d}^{2}+2N_{\rm d}N_{\rm r}+N_{\rm r}N^{\prime}_{\rm r}$ becomes

[TABLE]

meaning that $DR$ pairs contribute twice as much as $RR$ pairs to the variance and also twice as much computational cost is invested in them. The memory requirement for the random catalog is then the same as for the data catalog. The cost saving estimate above is optimistic, since the computation involves some overhead not proportional to the number of pairs.

For small scales, where $\xi\gg 1$, the situation is different. The greater density of $DD$ pairs due to the correlation requires a greater density of the $R$ catalog so that the additional variance from it is not greater. From Eq. (19) we see that the balance of the $DR$ and the $RR$ contributions is different for large $\xi$ (the $p_{\rm c}$ term vs. the other terms). We may consider recomputing $\hat{\xi}$ for the small scales using a smaller $r_{\mathrm{max}}$ and a larger $R$ catalog. Considering just the Poisson terms ($p_{\rm c}$ and $p_{\rm r}$ or $p^{\prime}_{\rm r}$) with a ‘representative’ $\xi$ value, (27) and (28) become $c=1+\xi+2M_{\rm r}+M_{\rm r}^{2}/M_{\rm s}$ and $v=2/M_{\rm r}+(1-\xi)^{2}M_{\rm s}/M_{\rm r}^{2}$ which modifies the above result (Eq. 29) for the optimal choice of $M_{\rm s}$ and $M_{\rm r}$ to

[TABLE]

This result is only indicative, since it assumes a constant $\xi$ for $r<r_{\mathrm{max}}$. In particular, it does not apply for $\xi\approx 1$, because then the approximation of ignoring the $q^{\rm r}$ and $t^{\rm r}$ terms in (19) is not good.

4 Tests on mock catalogs

4.1 Minerva simulations and methodology

The Minerva mocks are a set of 300 cosmological mocks produced with $N$-body simulations (Grieb et al. 2016; Lippich et al. 2019), stored at five output redshifts $z\in\{2.0,1.0,0.57,0.3,0\}$. The cosmology is flat $\Lambda$CDM with $\Omega_{m}=0.285$, and we use the $z=1$ outputs. The mocks have $N_{\rm d}\approx 4\times 10^{6}$ objects (“halos” found by a friends-of-friend algorithm) in a box of $1500${\mathrm{\mathnormal{h}}}^{-1}\text{,}\mathrm{Mpc}$$ cubed.

To model the survey geometry of a redshift bin with $\Delta z\approx 0.1$ at $z\sim 1$, we placed the observer at comoving distance $2284.63${\mathrm{\mathnormal{h}}}^{-1}\text{,}\mathrm{Mpc} from the center of the cube and selected from the cube a shell $2201.34$–$2367.92${\mathrm{\mathnormal{h}}}^{-1}\text{\,}\mathrm{Mpc} from the observer. The comoving thickness of the shell is $166.58${\mathrm{\mathnormal{h}}}^{-1}\text{,}\mathrm{Mpc}$$. The resulting mock sub-catalogs have $N_{\rm d}\approx 4.5\times 10^{5}$ and are representative of the galaxy number density of the future Euclid spectroscopic galaxy catalog.

We ignore peculiar velocities, that is, we perform our analysis in real space. Therefore, we consider results for the 1D 2PCF $\xi(r)$. We estimated $\xi(r)$ up to $r_{\mathrm{max}}=200${\mathrm{\mathnormal{h}}}^{-1}\text{,}\mathrm{Mpc} using $\Delta r=1${\mathrm{\mathnormal{h}}}^{-1}\text{\,}\mathrm{Mpc} bins.

We chose standard LS with $M_{\rm r}=50$ as the reference method. In the following, LS without further qualification refers to this. The random catalog was generated separately for each shell mock to measure their contribution to the variance. For one of the random catalogs we calculated also triplets to obtain the edge effect quantity $N_{\rm d}t_{12}=G^{\rm t}_{12}/G^{\rm p}_{1}G^{\rm p}_{2}-1$.

While dilution can already be discarded on theoretical grounds, we show results obtained using dilution, since these results provide the scale for edge effects demonstrating the importance of eliminating them with a careful choice of method. For the dilution and split methods we also used $M_{\rm r}=50$ , and tried out dilution fractions $d:=N^{\prime}_{\rm r}/N_{\rm r}=0.5,0.25,0.14$ and split factors $M_{\rm s}=4,16,50$ (chosen to have similar pairwise computational costs). In addition, we considered standard LS with $M_{\rm r}=25$, which has the same number of $RR$ pairs as $d=0.5$ and $M_{\rm s}=4$, but only half the number of $DR$ pairs; and standard LS with $M_{\rm r}=1$ to demonstrate the effect of a small $N_{\rm r}$.

The code used to estimate the 2PCF implements a highly optimized pair-counting method, specifically designed for the search of object pairs in a given range of separations. In particular, the code provides two alternative counting methods, the chain-mesh and the kd-tree. Both methods measure the exact number of object pairs in separation bins, without any approximation. However, since they implement different algorithms to search for pairs, they perform differently at different scales, both in terms of CPU time and memory usage. Overall, the efficiency of the two methods depends on the ratio between the scale range of the searching region and the maximum separation between the objects in the catalog.

The kd-tree method first constructs a space-partitioning data structure that is filled with catalog objects. The minimum and maximum separations probed by the objects are kept in the data structure and are used to prune object pairs with separations outside the range of interest. The tree pair search is performed through the dual-tree method in which cross-pairs between two dual trees are explored. This is an improvement in terms of exploration time over the single-tree method.

On the other hand, in the chain-mesh method the catalog is divided in cubic cells of equal size, and the indexes of the objects in each cell are stored in vectors. To avoid counting object pairs with separations outside the interest range, the counting is performed only on the cells in a selected range of distances from each object. The chain-mesh algorithm has been imported from the CosmoBolognaLib, a large set of free software C++/python libraries for cosmological calculations (Marulli et al. 2016).

For our test runs we used the chain-mesh method.

4.2 Variance and bias

In Fig. 1 we show the mean (over the 300 mock shells) estimated correlation function and the scatter (square root of the variance) of the estimates using the LS, split, and dilution methods; our theoretical approximate result for the scatter for LS; and our theoretical result for bias for the different methods.

The theoretical result for the scatter is shown with and without the $q$ terms, which include the triplet correlation $\xi_{12}$, for which we used here the approximation $\xi_{12}\approx\xi(\max(r_{1},r_{2}))$. This behaves as expected, that is, it underestimates the variance, since we neglected the higher-order correlations in the $D$ catalog. Nevertheless, it (see the dash-dotted line in Fig. 1) has similar features to the measured variance (dashed lines).

In Fig. 2 we plot the diagonals of the $p$, $t$, and $q$ quantities. This shows how their relative importance changes with separation scale. It also confirms that our initial assumption on small relative fluctuations is valid in this simulation case.

Consider now the variance differences (from standard LS with $M_{\rm r}=50$), for which our theoretical results should be accurate. Figure 3 compares the measured variance difference to the theoretical result. For the diluted estimators and LS with $M_{\rm r}=1$ the measured result agrees with theory, although clearly the measurement with just 300 mocks is rather noisy. For the split estimators and LS with $M_{\rm r}=25$ the difference is too small to be appreciated with 300 mocks, but at least the measurement does not disagree with the theoretical result.

In Fig. 4 we show the relative theoretical increase in scatter compared to the best possible case, which is LS in the limit $M_{\rm r}\rightarrow\infty$. Since we do not have a valid theoretical result for the total scatter, we estimate it by subtracting the theoretical difference from LS with $M_{\rm r}=50$ from the measured variance of the latter.

At scales $r\lesssim 10${\mathrm{\mathnormal{h}}}^{-1}\text{,}\mathrm{Mpc} the theoretical prediction is about the same for dilution and split and neither method looks promising for $r\ll 10${\mathrm{\mathnormal{h}}}^{-1}\text{\,}\mathrm{Mpc} where $\xi\gg 1$. This suggests that for optimizing cost and accuracy, a different method should be used for smaller scales than that used for large scales. The number of $RR$ pairs with small separations is much less. Therefore, for the small-scale computation there is no need to restrict the computation to a subset of $RR$ pairs, or alternatively, one can afford to increase $M_{\rm r}$. For the small scales, we may consider the split method with increased $M_{\rm r}$ as an alternative to standard LS. We have the same number of pairs to compute as in the reference LS case, if we use $M_{\rm r}=866$ and $M_{\rm s}=866$. We added this case to Fig. 4. It seems to perform better than LS at intermediate scales, but for the smallest scales LS has the smaller variance. This is in line with our conclusion in Sect. 3.2, which is that when $\xi\gg 1$, it is not optimal to split the $R$ catalog into small subsets.

We also compared the differences in the mean estimate from the different estimators to our theoretical results on the bias differences (see Fig. 5), but the theoretical bias differences are much smaller than the expected error of the mean from 300 mocks; and we simply confirm that the differences we see are consistent with the error of the mean and thus consistent with the true bias being much smaller. We also performed tests with completely random ($\xi=0$) mocks, and with a large number ($10\,000$) of mocks confirmed the theoretical bias result for the different estimators in this case. Since the bias is too small to be interesting we do not report these results in more detail here.

However, we note that for the estimation of the 2D 2PCF and its multipoles, the 2D bins will contain a smaller number of objects than the 1D bins of these test runs and therefore the bias is larger. Using the theoretical results (17) or (24) the bias can be removed afterwards with accuracy depending on how well we know the true $\xi$.

4.3 Computation time and variance

The test runs were made using a single full 24-core node for each run. Table 1 shows the mean computation time and mean estimator variance for different $r$ ranges for the different cases we tested. Of these $r$ ranges, the $r=80$–$120${\mathrm{\mathnormal{h}}}^{-1}\text{,}\mathrm{Mpc}$$ is perhaps the most interesting, since it contains the important BAO scale. Therefore, we plot the mean variance at this range versus mean computation time in Fig. 6, together with our theoretical predictions. The theoretical estimate for the computation time for other dilution fractions and split factors is

[TABLE]

assuming $M_{\rm r}=50$. For standard LS with other random catalog sizes, the computation time estimate is

[TABLE]

5 Conclusions

The computational time of the standard Landy–Szalay estimator is dominated by the $RR$ pairs. However, except at small scales where correlations are large, these make a negligible contribution to the expected error compared to the contribution from the $DD$ and $DR$ pairs. Therefore, a substantial saving of computation time with an insignificant loss of accuracy can be achieved by counting a smaller subset of $RR$ pairs.

We considered two ways to reduce the number of $RR$ pairs: dilution and split. In dilution, only a subset of the $R$ catalog is used for $RR$ pairs. In split, the $R$ catalog is split into a number of smaller subcatalogs, and only pairs within each subcatalog are counted. We derived theoretical results for the additional estimator covariance and bias due to the finite size of the random catalog for these different variants of the LS estimator, extending in many ways the original results by Landy & Szalay (1993), who worked in the limit of an infinite random catalog. We tested our results using 300 mock data catalogs, representative of the $z=0.95$–$1.05$ redshift range of the future Euclid survey. The split method maintains the property the Landy–Szalay estimator was designed for, namely cancelation of edge effects in bias and variance (for $\xi=0$), whereas dilution loses this cancellation and therefore should not be used.

For small scales, where correlations are large, one should not reduce $RR$ counts as much. The natural dividing line is the scale $r$ where $\xi(r)=1$. Interestingly, the difference in bias and covariance between the different estimators (split, dilution, and LS) vanishes when $\xi=1$. We recommend the natural version of the split method, $M_{\rm s}=M_{\rm r}$, for large scales where $|\xi|<1$. This leads to a saving in computation time by more than a factor of ten (assuming $M_{\rm r}=50$) with a negligible effect on variance and bias. For small scales, where $\xi>1$, one should consider using a larger random catalog and one can use either the standard LS method or the split method with a more modest split factor. Because the number of pairs with these small separations is much smaller, the computation time is not a similar issue as for large separations.

The results of our analysis will have an impact also on the computationally more demanding task of covariance matrix estimation. However, since in that case the exact computational cost is determined by the balance of data catalogs and random catalogs, which does not need to be the same as for the individual two-point correlation estimate, we postpone a quantitative analysis to a future, dedicated study. The same kind of methods can be applied to higher-order statistics (three-point and four-point correlation functions) to speed up their estimation (Slepian & Eisenstein 2015).

Acknowledgements.

We thank Will Percival and Cristiano Porciani for useful discussions. The 2PCF computations were done at the Euclid Science Data Center Finland (SDC-FI, urn:nbn:fi:research-infras-2016072529), for whose computational resources we thank CSC – IT Center for Science, the Finnish Grid and Cloud Computing Infrastructure (FGCI, urn:nbn:fi:research-infras-2016072533), and the Academy of Finland grant 292882. This work was supported by the Academy of Finland grant 295113. VL was supported by the Jenny and Antti Wihuri Foundation, AV by the Väisälä Foundation, AS by the Magnus Ehrnrooth Foundation, and JV by the Finnish Cultural Foundation. We also acknowledge travel support from the Jenny and Antti Wihuri Foundation. VA acknowledges funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No 749348. FM acknowledges the grants ASI n.I/023/12/0 “Attivit‘a relative alla fase B2/C per la missione Euclid” and PRIN MIUR 2015 “Cosmology and Fundamental Physics: illuminating the Dark Universe with Euclid”.

Appendix A Derivation examples

As examples of how the variances of the different deviations, $\langle\alpha\beta\rangle$ and so on, are derived, following the method presented in LS93, we give three of the cases here: $\langle\alpha_{1}\alpha_{2}\rangle$ (common for all the variants of LS), and $\langle\beta_{1}\gamma^{\prime}_{2}\rangle$ and $\langle\gamma^{\prime}_{1}\gamma^{\prime}_{2}\rangle$ for the split method. The rest are calculated in a similar manner.

To derive the $\langle DD_{1}\ DD_{2}\rangle$ appearing in the $\langle\alpha_{1}\alpha_{2}\rangle$ of (15) we start from

[TABLE]

where both $i,j,$ and $k,l$ sum over all microcell pairs; $n_{i}=1\mbox{\ or\ }0$ is the number of galaxies in microcell $i$. There are three different cases for the terms $\langle n_{i}n_{j}n_{k}n_{l}\rangle$ depending on how many indices are equal ($i\neq j$ and $k\neq l$ for all of them).

The first case (quadruplets, i.e., two pairs of microcells) is when $i,j,k,l$ are all different. We use $\sum^{*}_{ijkl}$ to denote this part of the sum. There are ${\textstyle\frac{1}{2}}K(K-1)\times{\textstyle\frac{1}{2}}(K-2)(K-3)={\textstyle\frac{1}{4}}K^{4}$ (we work in the limit $K\rightarrow\infty$) such terms and they have

[TABLE]

where $\delta_{i}:=\delta({\bf x}_{i}$) is the density perturbation (relative to the actual mean density of galaxies in the survey), $\zeta$ is the three-point correlation, and $\eta$ is the connected (i.e., the part that does not arise from the two-point correlation) four-point correlation. The fraction of microcell quadruplets (pairs of pairs) that satisfy ${\bf r}_{ij}\equiv{\bf x}_{j}-{\bf x}_{i}\in{\bf r}_{1}$ and ${\bf r}_{kl}\in{\bf r}_{2}$ is $G^{\rm p}({\bf r}_{1})G^{\rm p}({\bf r}_{2})=:G^{\rm p}_{1}G^{\rm p}_{2}$. In the limit of large $K$ the number of other index quadruplets is negligible compared to those where all indices have different values, so we have

[TABLE]

For the connected pairs $(i,j)$ and $(k,l)$ we have $\xi({\bf r}_{ij})=\xi({\bf r}_{1})\equiv\xi_{1}$ and $\xi({\bf r}_{kl})=\xi({\bf r}_{2})\equiv\xi_{2}$.

The second case (triplets of microcells) is when $k$ or $l$ is equal to $i$ or $j$. We denote this part of the sum by

[TABLE]

It turns out that it goes over all ordered combinations of $\{i,j,k\}$, where $i,j,k$ are all different, exactly once, so there are $K(K-1)(K-2))=K^{3}$ such terms (triplets). Here

[TABLE]

and

[TABLE]

The third case (pairs of microcells) is when $i=k$ and $j=l$. This part of the sum becomes

[TABLE]

where

[TABLE]

and

[TABLE]

that is, the sum vanishes unless the two bins are the same, ${\bf r}_{1}={\bf r}_{2}$.

We now apply the three approximations listed in Sect. 2.2:

1. $\xi({\bf r}_{ik})=\xi({\bf r}_{il})=\xi({\bf r}_{jk})=\xi({\bf r}_{jl})=0$ in (36);
2. $\zeta=0$ in (36) and (39);
3. $\eta=0$ in (36). We obtain

[TABLE]

and using $\langle DD\rangle$ from (3) we arrive at the $\langle\alpha_{1}\alpha_{2}\rangle$ result given in Eq. (15).

For the split method, we give the calculation for

[TABLE]

below.

First the $\langle\gamma^{\prime}_{1}\beta_{2}\rangle$:

[TABLE]

where

[TABLE]

$s_{i}$ is the number (0 or 1) of $R^{\mu}$ objects in microcell $i$, and $r_{l}$ is the number of $R$ objects in microcell $l$.

There are ${\textstyle\frac{1}{2}}K^{4}$ quadruplet terms, for which (see (37)) $\sum^{\ast}\Theta\Theta={\textstyle\frac{1}{2}}K^{4}G^{\rm p}_{1}G^{\rm p}_{2}$, and

[TABLE]

Triplets where $i$ or $j$ is equal to $k$ have $\langle s_{i}s_{k}n_{k}r_{l}\rangle=0$, since $s_{k}n_{k}=0$ always (we cannot have two objects, from $R^{\mu}$ and $D$, in the same cell). There are $K^{3}$ triplet terms where $i$ or $j$ is equal to $l$: for them $\sum^{\ast}\Theta\Theta=K^{3}G^{\rm t}_{12}$, and

[TABLE]

where if $s_{l}=1$ then also $r_{l}=1$ since $R^{\mu}\subset R$.

Pairs where $(i,j)=(k,l)$ or $(i,j)=(l,k)$ have $\langle s_{k}s_{l}n_{k}r_{l}\rangle=0$, since again we cannot have two different objects in the same cell. Thus

[TABLE]

and we obtain that

[TABLE]

which is equal to the $\langle\gamma_{1}\beta_{2}\rangle$ of standard LS.

Then the $\langle\gamma^{\prime}_{1}\gamma^{\prime}_{2}\rangle$:

[TABLE]

where there are $M_{\rm s}(M_{\rm s}-1)$ terms with $\mu\neq\nu$ giving

[TABLE]

and $M_{\rm s}$ terms with $\mu=\nu$ giving

[TABLE]

Adding these up gives

[TABLE]

and

[TABLE]

This differs both from standard LS and from dilution, since it involves both $N_{\rm r}$ and $N^{\prime}_{\rm r}$.