© F. Bernardeau 2022

1. Introduction

In the context of cosmological studies, the concept of counts-in-cells statistics has been put forward for a long time as a unique way to quantify the statistical properties of the cosmological fields (White 1979; Colombi et al. 1995; Balian & Schaeffer 1989; Bernardeau & Schaeffer 1999). It was then shown in particular that counts-in-cells statistics, which represents a discrete representation of the local density probability distribution function (PDF), could be directly related to the correlation hierarchy of the density field.

Interest in these types of observables was recently renewed for several reasons. The size of the surveys makes accurately measuring these quantities more realistic. This is already the case for surveys such as the Dark Energy Survey (DES collaboration; Abbott et al. 2018), the Kilo-Degree Survey (KIDS; Heymans et al. 2021), and the Hyper Suprime Cam (HSC; Hikage et al. 2019). The future promises even larger and more powerful surveys such as Euclid (Laureijs et al. 2011; Amendola et al. 2018) and the Rubin Observatory (Ivezić et al. 2019). Moreover, the theoretical foundations for these constructions (at least in the cosmological context) has been considerably strengthen with the realization that the large-deviation theory (LDT; for a general review, see Touchette 2011) could successfully be invoked, as shown in Bernardeau & Reimberg (2016). It clarifies the applicability of the theory to the cosmological density field and places previous works on a much more solid foundation (Valageas 2002; Bernardeau et al. 2014). The ability of density PDF to constrain cosmology was emphasized in Codis et al. (2016b) and completed in Friedrich et al. (2020) and in Uhlemann et al. (2020), who showed that these observable could efficiently constrain the neutrino mass or primordial non-Gaussianities. Finally, although the matter PDF is not a direct observable, as is matter density, it can be closely approached with the help of luminous tracer statistics (Repp & Szapudi 2020), more convincingly in weak-lensing fields, as advocated in numerous recent papers (Barthelemy et al. 2021; Bernardeau & Valageas 2000), or with combined approaches such as density-split statistics (Gruen et al. 2018; Friedrich et al. 2018; Brouwer et al. 2018), which proved to be particularly promising.

The construction of a full theory of these observable requires a detailed analysis of its global error budget, however, due to finite-size surveys, imperfect tracers, and so on. Some of these aspects have been explored in early studies such as Szapudi & Colombi (1996) and Szapudi et al. (1999), but a full theory is still lacking. The developments presented in this paper are made in this context. More precisely, the purpose of this study is to explore what determines the expression of the covariance of data vectors whose elements are local quantities, such as the density contrast and density profiles, in cosmological contexts, that is, in classical random fields with long-range correlations. Derivations were made furthermore assuming statistical homogeneity and isotropy. The domain of application encompasses both counts-in-cells statistics, basically 2D or 3D counts of density tracers, or proxies to projected densities such as mass maps for weak-lensing tomographic observations.

In order to gain insights into the different contributions and the effects that might contribute to the covariance, we rely on the use of the hierarchical models to derive results we think rather general. The immense advantage of using such models is that they naturally incorporate many of the features expected in density cosmological fields (e.g., the magnitude of the high-order correlation functions), and there are also models for which many exact results can be obtained in particular for counts-in-cells statistics. The goal of these constructions is to eventually infer precisely what the performance of PDF measurements would be on the determination of cosmological parameters, taking advantage of results such as those found in Boyle et al. (2021), which give the response function of these observable to various cosmological parameters

Section 2 is devoted to the presentation of the general framework. The subsequent section explores different contributions, from large-scale effects with the derivation of several bias functions to short-distance contributions. Results are derived in a framework as general as possible, including discrete noise associated with the use of a finite number of tracers. Section 4 presents the general hierarchical models, and more specifically, the Rayleigh-Levy flight model that we use as a toy model to evaluate the performances of approximate schemes. In Sect. 5, simplified models for the covariance matrix are presented and evaluated with the help of a set of numerical experiments. Section 6 summarizes the results that have be found and specifies their expected range of application.

The text is complemented by appendices that contain a large amount of material. They present the hierarchical models, their mathematical description, and the mean-field approximation that is used throughout for explicit computations. Appendix C is more specifically devoted to the minimal tree model and the construction of the exact mean-field covariance matrix.

2. General framework. Construction of covariance matrices

The purpose of this section is to show how the elements of the covariance matrix are related to the joint density PDFs within a given survey. We first formalize this relation in a general framework before we explore its consequence in the context we are interested in. We assume we are interested in the PDF of some local quantity, μ, that can be evaluated within a survey, thus defining a field μ(x) throughout the survey. The a priori typical example of this quantity is the density (see below for a more precise illustration of what this quantity could be). The value of μ is assumed to lie in some ensemble ℳ (that can be simply the real numbers), and the data vector we are interested in consists of the probabilities p_i that μ lie within the subsets ℳ_i (which are a priori nonzero within ℳ). The one-point PDF of μ is then given by

$$\begin{array}{c}\hfill p_i(\mathbf{x})={\displaystyle {\int }_{{\mathcal{M}}_i}}\mathrm{d}\mu \mathcal{P}(\mu ,\mathbf{x}),\end{array}$$(1)

if ℘(μ,x) dμ is the PDF of μ at location x. p_i(x) is then assumed to be independent of x in the context we are interested in, for which statistical homogeneity is assumed. More formally, we can define the characteristic function χ_i(x), which takes the value 1 where μ(x) ∈ ℳ_i and 0 otherwise.

An estimation of p_i would then be given by the volume fraction of the survey where μ(x) ∈ ℳ_i. We note this estimate as P_i ¹,

$$\begin{array}{c}\hfill P_i=\frac{1}{V}{\displaystyle \int \mathrm{d}}\mathbf{x}{\chi }_i(\mathbf{x}),\end{array}$$(2)

which is then itself a random variable, the properties of which we are interested in. More precisely, we would like to derive an operational form for the likelihood function of a set of P_i variable. We limit our investigation here to the construction of the likelihood from the covariance matrix, assuming that the likelihood of P_i is close enough to a Gaussian distribution².

The ensemble average of P_i is

$$\begin{array}{c}\hfill ⟨P_i⟩=\frac{1}{V}{\displaystyle \int \mathrm{d}}\mathbf{x}⟨{\chi }_i(\mathbf{x})⟩=\frac{1}{V}{\displaystyle \int \mathrm{d}}\mathbf{x}{p}_i(\mathbf{x})=p_{i.}\end{array}$$(3)

We can further define a joint PDF of the same field, ℘(μ,x;μ′,x′), which is the joint PDF of μ and μ′ in locations x and x′. Defining p_ij(x, x′) as the joint ensemble average of ℘(μ,x;μ′,x′), we have

$$\begin{array}{c}\hfill p_{\mathit{ij}}(\mathbf{x},{\mathbf{x}}^{\prime })={\displaystyle {\int }_{{\mathcal{M}}_i}}\mathrm{d}\mu {\displaystyle {\int }_{{\mathcal{M}}_j}}\mathrm{d}{\mu }^{\prime }\mathcal{P}(\mu ,\mathbf{x};{\mu }^{\prime },{\mathbf{x}}^{\prime }).\end{array}$$(4)

The elements of the covariance matrix of P_i are then formally

$$\begin{array}{cc}& ⟨P_i{P}_j⟩=\frac{1}{V^2}{\displaystyle \int \mathrm{d}}\mathbf{x}{\displaystyle \int \mathrm{d}}{\mathbf{x}}^{\prime }⟨{\chi }_i(\mathbf{x}){\chi }_j({\mathbf{x}}^{\prime })⟩\hfill \\ \hfill & =\frac{1}{V^2}{\displaystyle \int \mathrm{d}}\mathbf{x}{\displaystyle \int \mathrm{d}}{\mathbf{x}}^{\prime }{p}_{\mathit{ij}}(\mathbf{x},{\mathbf{x}}^{\prime })\equiv {\overline{p}}_{\mathit{ij}}.\hfill \end{array}$$(5)

This gives the relation between the covariances and joint PDF. If p_ij(x, x′) depends only on the relative distance between x and x′, this expression can be recast in terms of the distribution of such distances, P_s(r_d), in the form

$$\begin{array}{c}\hfill ⟨P_i{P}_j⟩={\displaystyle \int \mathrm{d}}{r}_d{P}_s(r_d)p_{\mathit{ij}}(r_d).\end{array}$$(6)

The precise form of P_s(r_d) depends on the detail of the survey. Explicit forms can be given in case of simple regular surveys such as square surveys³. In the context of statistically homogeneous and isotropic random fields, this latter expression is used. In particular, we wish to determine the configurations that contribute most to p̄_ij. They obviously depend on both the random processes we consider and on the definition of ℳ_i and ℳ_j. In order to be more specific, we assume in the following that μ is a local density assigned to be in bins (i) centered on ρ_i and with width dρ_i (assumed a priori to be arbitrarily small), so that

$$\begin{array}{c}\hfill P_i=P({\rho }_i)\mathrm{d}{\rho }_i.\end{array}$$(7)

If necessary, local densities could be obtained after the field μ(x) has been convolved with a window function W_R(x), associated with a scale R that is

$$\begin{array}{c}\hfill \rho (\mathbf{x})={\displaystyle \int \mathrm{d}}{\mathbf{x}}^{\prime }\mu (\mathbf{x}-{\mathbf{x}}^{\prime })W_R({\mathbf{x}}^{\prime }).\end{array}$$(8)

It is then assumed R is small compared to the sample size in order to identify what the leading contributions to the joint PDFs might be. In practice, W_R might also be a simple top-hat window function, but this is not necessarily the case. It could be more elaborated filters, such as compensated filters (of zero average), such as those introduced for cosmic shear analysis (Schneider 1996; Kaiser 1998; Bernardeau & Valageas 2000).

We furthermore allow the estimated densities ρ_i to be defined with respect the overall density of the sample ρ_s,

$$\begin{array}{c}\hfill {\rho }_s=\frac{1}{V}{\displaystyle \int \mathrm{d}}\mathbf{x}\mu (\mathbf{x}).\end{array}$$(9)

For instance, we could be interested in ${\widehat{\rho }}_i\equiv {\rho }_i/{\rho }_s$ or ${\overline{\rho }}_i\equiv {\rho }_i-({\rho }_s-1),$ which are frequently encountered situations in praxis. Then ρ_s is itself a random variable whose correlation with ρ(x) ought to be taken into account. We then need to explore the properties of either P(ρ_i, ρ_j; x, x′) or P(ρ_s, ρ_i, ρ_j; x, x′) from which functions of interest can be built, that is,

$$\begin{array}{cc}\hfill P(\widehat{\rho },{\widehat{\rho }}^{\prime })& ={\displaystyle \int \mathrm{d}}{\rho }_s{\rho }_s^{2}P({\rho }_s,\widehat{\rho }{\rho }_s,{\widehat{\rho }}^{\prime }{\rho }_s;\mathbf{x},{\mathbf{x}}^{\prime })\hfill \end{array}$$(10)

$$\begin{array}{cc}\hfill P(\overline{\rho },{\overline{\rho }}^{\prime })& ={\displaystyle \int \mathrm{d}}{\rho }_{s}P({\rho }_s,{\overline{\rho }}_i-1+{\rho }^{\prime },{\overline{\rho }}^{\prime }-1+{\rho }_s;\mathbf{x},{\mathbf{x}}^{\prime }),\hfill \end{array}$$(11)

from which the covariance elements such as

$$\begin{array}{cc}& \mathrm{Cov}({\rho }_i,{\rho }_j)\mathrm{d}{\rho }_i\mathrm{d}{\rho }_j={\displaystyle \int \mathrm{d}}{r}_d{P}_s(r_d)P({\rho }_i,{\rho }_j;r_d)\mathrm{d}{\rho }_i\mathrm{d}{\rho }_j\hfill \\ \hfill & -P({\rho }_i)P({\rho }_j)\mathrm{d}{\rho }_i\mathrm{d}{\rho }_j\hfill \end{array}$$(12)

can be derived and whose properties we wish to explore. We wish in particular to build a model of the likelihood function from such a covariance, requiring full knowledge of its eigenvalues and eigendirections. In this respect, it is implicit that the number of bins (i) to be used is finite. We nonetheless present at least in this first section the results in the continuous limit for ρ_i. It is finally to be noted that as stated before, we restrict our analysis to covariance matrices, but higher-order correllators might also be considered by generalizing the relation (5) to a higher number of variables.

3. PDF covariances in the context of cosmological models

3.1. Modeling the joint PDF

To make progress, we need to make further assumptions about the mathematical structure of the joint PDF. In the following, we assume in particular that joint PDFs can be obtained from their cumulant generating functions (CGF)⁴, φ(λ_i, λ_j; r_d). The latter is defined as

$$\begin{array}{cc}& exp(\phi ({\lambda }_i,{\lambda }_j;r_d))=⟨exp\left({\lambda }_i{\rho }_i+{\lambda }_j{\rho }_j\right)⟩\hfill \\ \hfill & ={\displaystyle \int \mathrm{d}}{\rho }_i\mathrm{d}{\rho }_{j}P({\rho }_i,{\rho }_j;r_d)exp\left({\lambda }_i{\rho }_i+{\lambda }_j{\rho }_j\right),\hfill \end{array}$$(13)

and it is assumed that this relation can be inverted to give the joint PDF from Laplace inverse transformations,

$$\begin{array}{c}\hfill P({\rho }_i,{\rho }_j;r_d)={\displaystyle \int \frac{\mathrm{d}{\lambda }_i}{2\pi \mathrm{i}}}\frac{\mathrm{d}{\lambda }_j}{2\pi \mathrm{i}}{e}^{-{\lambda }_i{\rho }_i-{\lambda }_j{\rho }_j+\phi ({\lambda }_i,{\lambda }_j;r_d)},\end{array}$$(14)

where the integrations are made a priori along the imaginary axis. The CGFs themselves are closely related to the averaged correlation functions of the underlying field. In the following, we develop models for which these correlation functions can be computed precisely.

3.2. Large-distance contributions to the joint density PDF

We start by assuming that covariances are dominated by long-range correlation and not by proximity effects (e.g., densities taken in nearby regions). Whether this assumption is correct obviously depends on the particular model and setting we consider, as we detail below. There are large sets of models for which general expressions can be given in this regime. They are the so-called hierarchical models, originally introduced in Peebles (1980), discussed in more detail in Fry (1984a,b), Balian & Schaeffer (1989), Bernardeau & Schaeffer (1992), and further formalized in Bernardeau & Schaeffer (1999) as described below; it is also true in the quasilinear regime as originally pointed out in Bernardeau (1996) and derived in more detail in Codis et al. (2016a). In these regimes, we obtain the following functional form (see previous references and the detailed derivation in Appendix B):

$$\begin{array}{cc}& \phi ({\lambda }_s,{\lambda }_i,{\lambda }_j)={\lambda }_s+{\phi }_0({\lambda }_i)+{\phi }_0({\lambda }_j)\hfill \\ \hfill & +\frac{{\lambda }_s^2}{2}{\displaystyle \int \mathrm{d}}{\mathbf{x}}_s\mathrm{d}{\mathbf{x}}_s^{\prime }\xi ({\mathbf{x}}_s,{\mathbf{x}}_s^{\prime })+{\lambda }_s{\displaystyle \int \mathrm{d}}{\mathbf{x}}_s\xi ({\mathbf{x}}_s,{\mathbf{x}}_1){\phi }_1({\lambda }_i)\hfill \\ \hfill & +{\lambda }_s{\displaystyle \int \mathrm{d}}{\mathbf{x}}_s\xi ({\mathbf{x}}_s,{\mathbf{x}}_2){\phi }_1({\lambda }_j)+{\phi }_1({\lambda }_i)\xi ({\mathbf{x}}_1,{\mathbf{x}}_2){\phi }_1({\lambda }_j),\hfill \end{array}$$(15)

where ξ(x, x′) is the two-point correlation function of the density field at positions x and x′, and φ₀(λ) and φ₁(λ) are specific functions of λ that depend on the details of the model.

Then, setting λ_s to zero, we can easily obtain the expression of the joint PDF at leading order in ξ(r_d),

$$\begin{array}{c}\hfill P({\rho }_i,{\rho }_j;r_d)=P({\rho }_i)P({\rho }_j)(1+b({\rho }_i)\xi (r_d)b({\rho }_j)).\end{array}$$(16)

Here P(ρ_i) is the one-point density PDF (i.e., implicitly at scale R), and b(ρ_i) is the density-bias function (to be distinguished from the standard halo-bias function). It also depends on ρ_i (and on the scale R) so that in the previous expression, the dependence on ρ_i, ρ_j, and r_d can be factorized out.

To be more precise, P(ρ_i) is given by the inverse Laplace transform of the CGF (see, e.g., Balian & Schaeffer 1989 and Bernardeau 2013 for a detailed derivation of this inversion),

$$\begin{array}{c}\hfill P({\rho }_i)={\displaystyle \int \frac{\mathrm{d}\lambda }{2\pi \mathrm{i}}}exp(-\mathrm{i}\lambda {\rho }_i+{\phi }_0(\lambda )),\end{array}$$(17)

where φ₀(λ) is the CGF of the density taken at scale R (i.e., for the filter W_R). The function b(ρ_i) is defined through a similar relation,

$$\begin{array}{c}\hfill b({\rho }_i)P({\rho }_i)=\int \frac{\mathrm{d}\lambda }{2\pi \mathrm{i}}{\phi }_1(\lambda )exp(-\mathrm{i}\lambda {\rho }_i+{\phi }_0(\lambda )).\end{array}$$(18)

The function φ₁(λ) can be explicitly computed in the context of perturbation theory calculations (Codis et al. 2016a). This is the case also for models in the so-called hierarchical models (see appendices). In the latter case, these calculations can be extended to higher order, as we describe below, providing ways to better assess the domain of validity of this expansion.

According to the previous relation, this implies that this form translates into the expression of the covariance coefficients for the density PDF. More precisely, we expect

$$\begin{array}{c}\hfill \mathrm{Cov}({\rho }_i,{\rho }_j)=b({\rho }_i)P({\rho }_i){\overline{\xi }}_{s}b({\rho }_j)P({\rho }_j),\end{array}$$(19)

where ${\overline{\xi }}_s$ is the average value of the two-point correlation function ξ(r_d) within the sample.

It is to be noted, however, that this is true if

• The term in ${\overline{\xi }}_s$ is indeed the leading contribution of the expansion (15). This is obviously not the case for samples with periodic boundary conditions, for which ${\overline{\xi }}_s$ vanishes by construction;

• The density is defined regardless of the density of the sample. Its expectation value therefore does not coincide with ρ_s for a given sample.

It can also be noted that in the Gaussian limit, we have b(ρ_i)=δ_i/ξ. Applying the relation (16) to the density within one cell and to the sample density ρ_s = 1 + δ_s leads then to the following expression for the conditional expression of density PDF,

$$\begin{array}{c}\hfill P({\rho }_i|{\rho }_s)=P({\rho }_i)(1+{\delta }_{s}b({\rho }_i)).\end{array}$$(20)

This leads to the interpretation of the function b(ρ_i) as the response function of the density PDF to the sample density. This means that although the density-bias function cannot be derived from the density PDF alone, we should be able to derive it if we are in possession of an operational method to compute the density PDF for arbitrary cosmological parameters (in a way similar to the derivation of halo-bias function as pioneered in Mo & White 1996). Undoubtedly, this result is closely related to the so-called supersample effects (as described for the power spectra covariance in Takada & Hu 2013), that is, the impact of modes of scale comparable to or larger than the sample size. This is not necessarily their only contribution (because subdominant large-scale contributions can also contribute to the covariance), however, but likely to be the most important contribution, as described below.

The density-bias function obeys the following consistency relations:

$$\begin{array}{cc}& {\displaystyle \int \mathrm{d}}\rho b(\rho )P(\rho )=0,\hfill \end{array}$$(21)

$$\begin{array}{cc}& {\displaystyle \int \mathrm{d}}\rho \rho b(\rho )P(\rho )=1,\hfill \end{array}$$(22)

as initially pointed out in Bernardeau & Schaeffer (1992).

3.3. Case of relative densities

The previous formula applies to the local densities, evaluated regardless of the sample density. It does not apply in particular to standard settings (e.g., densities measured out of galaxy counts) where the density is defined with respect to the mean density of the sample. To address this case in particular, we should consider

$$\begin{array}{c}\hfill {\widehat{\rho }}_i=\frac{{\rho }_i}{{\rho }_{\mathrm{s}}}\end{array}$$(23)

as the observable for which the covariance is to be computed. In this case, the formal derivation of the PDFs is presented in the appendix, and it leads to

$$\begin{array}{cc}& P({\widehat{\rho }}_i)=\int \frac{\mathrm{d}{\lambda }_i}{2\pi \mathrm{i}}{\left[\frac{\partial \phi }{\partial {\lambda }_s}\right]}_{{|}_{{\lambda }_s=-{\lambda }_i{\widehat{\rho }}_i}}exp[\phi (-{\lambda }_i{\widehat{\rho }}_i,{\lambda }_i)]\hfill \end{array}$$(24)

$$\begin{array}{cc}& P({\widehat{\rho }}_i,{\widehat{\rho }}_j)=\int \frac{\mathrm{d}{\lambda }_i}{2\pi \mathrm{i}}\frac{\mathrm{d}{\lambda }_j}{2\pi \mathrm{i}}{[{\left(\frac{\partial \phi }{\partial {\lambda }_s}\right)}^2+\frac{{\partial }^2\phi }{\partial {\lambda }_s^2}]}_{{|}_{{\lambda }_s=-{\lambda }_i{\widehat{\rho }}_i-{\lambda }_j{\widehat{\rho }}_j}}\hfill \\ \hfill & \times exp[\phi (-{\lambda }_i{\widehat{\rho }}_i-{\lambda }_j{\widehat{\rho }}_j,{\lambda }_i,{\lambda }_j)].\hfill \end{array}$$(25)

We then use relation (15) to compute the form of this function. At this stage, it is to be noted that the expressions $\int \mathrm{d}{\mathbf{x}}_0\mathrm{d}{\mathbf{x}}_0^{\mathrm{\prime }}\xi ({\mathbf{x}}_0,{\mathbf{x}}_0^{\mathrm{\prime }})$, ∫dx₀ ξ(x₀, x₁) and ξ₁₂ all take the same averaged value when integrated over the sample. We note here this common value as ${\overline{\xi }}_s$ . Inserting the resulting expressions of the CGF in both the expressions of $P({\widehat{\rho }}_i)$ and $P({\widehat{\rho }}_i,{\widehat{\rho }}_j)$ and expanding all terms at linear oder in ${\overline{\xi }}_s$, we obtain

$$\begin{array}{cc}& P({\widehat{\rho }}_i,{\widehat{\rho }}_j)-P({\widehat{\rho }}_i)P({\widehat{\rho }}_j)={\overline{\xi }}_s{\displaystyle \int \frac{\mathrm{d}{\lambda }_i}{2\pi \mathrm{i}}}\frac{\mathrm{d}{\lambda }_j}{2\pi \mathrm{i}}\hfill \\ \hfill & \times (1+{\phi }_1({\lambda }_i)-{\lambda }_i{\widehat{\rho }}_i)(1+{\phi }_1({\lambda }_j)-{\lambda }_j{\widehat{\rho }}_j)\hfill \\ \hfill & \times exp[-{\lambda }_i{\widehat{\rho }}_i-{\lambda }_j{\widehat{\rho }}_j+{\phi }_0({\lambda }_i)+{\phi }_0({\lambda }_j)].\hfill \end{array}$$(26)

This leads to the definition of the first sample-bias function,

$$\begin{array}{c}\hfill b_{\mathrm{s}1}({\widehat{\rho }}_i)=\frac{1}{P({\widehat{\rho }}_i)}{\displaystyle \int \frac{\mathrm{d}\lambda }{2\pi \mathrm{i}}}(1+{\phi }_1(\lambda )-\lambda {\widehat{\rho }}_i)exp[-\lambda {\widehat{\rho }}_i+{\phi }_0(\lambda )],\end{array}$$(27)

which can be re-expressed in terms of the density-bias function defined in Eq. (18) and the derivative of $P({\widehat{\rho }}_i)$ with respect to ${\widehat{\rho }}_i$

$$\begin{array}{c}\hfill b_{\mathrm{s}1}({\widehat{\rho }}_i)=b({\widehat{\rho }}_i)+1+\frac{\partial log(P({\widehat{\rho }}_i))}{\partial log{\widehat{\rho }}_i}.\end{array}$$(28)

In this case, the covariance matrix elements are then expected to be given by

$$\begin{array}{c}\hfill \mathrm{Cov}({\widehat{\rho }}_i,{\widehat{\rho }}_j)=b_{\mathrm{s}1}({\widehat{\rho }}_i)P({\widehat{\rho }}_i){\overline{\xi }}_s{b}_{\mathrm{s}1}({\widehat{\rho }}_j)P({\widehat{\rho }}_j).\end{array}$$(29)

Remarkably, b_s1(ρ) can entirely be expressed in terms of b(ρ). For the sake of completeness, we also consider the case of ${\overline{\rho }}_i={\rho }_i-({\rho }_s-1)$. In this case, it is easy to show that

$$\begin{array}{cc}& P({\overline{\rho }}_i)={\displaystyle \int \frac{\mathrm{d}{\lambda }_i}{2\pi \mathrm{i}}}exp[-{\lambda }_i{\overline{\rho }}_i+\phi (-{\lambda }_i,{\lambda }_i)]\hfill \end{array}$$(30)

$$\begin{array}{cc}& P({\overline{\rho }}_i,{\overline{\rho }}_j)={\displaystyle \int \frac{\mathrm{d}{\lambda }_i}{2\pi \mathrm{i}}}\frac{\mathrm{d}{\lambda }_j}{2\pi \mathrm{i}}\hfill \\ \hfill & \times exp[-{\lambda }_i{\overline{\rho }}_i-{\lambda }_j{\overline{\rho }}_j+\phi (-{\lambda }_i-{\lambda }_j,{\lambda }_i,{\lambda }_j)].\hfill \end{array}$$(31)

Following the same approach as for the previous case, the leading-order expression in ${\overline{\xi }}_s$ of the connected joint PDF is

$$\begin{array}{cc}& P({\overline{\rho }}_i,{\overline{\rho }}_j)-P({\overline{\rho }}_i)P({\overline{\rho }}_j)={\overline{\xi }}_s{\displaystyle \int \frac{\mathrm{d}{\lambda }_i}{2\pi \mathrm{i}}}\frac{\mathrm{d}{\lambda }_j}{2\pi \mathrm{i}}\hfill \\ \hfill & \times ({\phi }_1({\lambda }_i)-{\lambda }_i)({\phi }_1({\lambda }_j)-{\lambda }_j)\hfill \\ \hfill & \times exp[-{\lambda }_i{\overline{\rho }}_i-{\lambda }_j{\overline{\rho }}_j+{\phi }_0({\lambda }_i)+{\phi }_0({\lambda }_j)].\hfill \end{array}$$(32)

It leads to the definition of the second sample-bias function,

$$\begin{array}{c}\hfill b_{\mathrm{s}2}({\overline{\rho }}_i)=b({\overline{\rho }}_i)+\frac{\partial log(P({\overline{\rho }}_i))}{\partial {\overline{\rho }}_i},\end{array}$$(33)

so that

$$\begin{array}{c}\hfill \mathrm{Cov}({\overline{\rho }}_i,{\overline{\rho }}_j)=b_{\mathrm{s}2}({\overline{\rho }}_i)P({\overline{\rho }}_i){\overline{\xi }}_s{b}_{\mathrm{s}2}({\overline{\rho }}_j)P({\overline{\rho }}_j).\end{array}$$(34)

The three bias functions are therefore closely related. Although the density-bias function b(ρ) cannot be derived from the shape of P(ρ) alone, as mentioned before, the relations between b(ρ) and either b_s1(ρ) and b_s2(ρ) depend on the PDF alone. Furthermore, the two relative density bias functions obey the following consistency relations:

$$\begin{array}{cc}& {\displaystyle \int \mathrm{d}}{\rho }_i{b}_{s\#}({\rho }_i)P({\rho }_i)=0\hfill \end{array}$$(35)

$$\begin{array}{cc}& {\displaystyle \int \mathrm{d}}{\rho }_i{\rho }_i{b}_{s\#}({\rho }_i)P({\rho }_i)=0.\hfill \end{array}$$(36)

The second relation is at variance with the corresponding relation (22) for the density-bias function. It indicates that for typical values of ρP(ρ), the sample bias functions, b_s#(ρ), are likely to be smaller than the density-bias function b(ρ).

3.4. Structure of the covariance matrix

The consequences of these formulae on the structure of the covariance matrix are illustrated below with the help of the Rayleigh-Levy flight model. Figure 3 compares the results from exact derivations of the covariance matrix with these prescriptions. The diagonal parts of the covariance matrices are well accounted for by these formulae. The root mean square of the measured local density PDF in particular exhibits the expected density dependence, at least for mild values of the density.

In all the formulae (19), (29), and (34), the expression of the covariance exhibits a simple structure, as it is factorizable in the two densities. This implies, for instance, that the reduced covariance matrix

$$\begin{array}{c}\hfill {\mathrm{Cov}}_{\mathrm{reduced}}({\rho }_i,{\rho }_j)=\frac{\mathrm{Cov}({\rho }_i,{\rho }_j)}{\sqrt{\mathrm{Cov}({\rho }_i,{\rho }_i)\mathrm{Cov}({\rho }_j,{\rho }_j)}}\end{array}$$(37)

has an extremely simple structure: it is given by the sign of the product of the bias functions (i.e., sign[b(ρ_i)b(ρ_j)], sign[b_s1(ρ_i)b_s1(ρ_j)], and sign[b_s2(ρ_i)b_s2(ρ_j)] for the three different measurement strategies). This leads to the butterfly-like structure in the plotted matrices, as illustrated in Fig. 4. This simple form betrays the fact that the density covariance is only poorly known. To be more specific, formulae (19), (29), and (34) give only a single eigendirection of the covariance matrix (namely b(ρ_i)P(ρ_i)) and the amplitude of a single eigenvalue associated with it. The numerical calculations suggest that it is the leading one when ${\overline{\xi }}_s$ does not vanish, as illustrated on Fig. 7. These formulae do not offer any indication of the amplitude of the covariance in orthogonal directions, however. Taken at face value, they imply that the other eigenvalues all vanish, preventing the covariance matrix from being invertible. These formulae therefore cannot be used alone to model the covariance for practical purposes, and complementary contributions have to be derived from other (and a priori subdominant) effects.

3.5. Beyond leading-order effects

In the previous subsection, we identified the long-distance leading contributions. As mentioned before, this leads to only limited information of the covariance structure. This difficulty is even more acute for covariances evaluated in numerical experiments consisting of a collection of independent samples, each of them with periodic boundary conditions (this does not have to be so, but it is often the case in practice). By construction, the mean correlation function within the sample then vanishes, ${\overline{\xi }}_s\to 0,$ making the term we have computed identically zero. All these considerations indicate that further contributions need to be identified. The identification of the next-to-leading order effects in Eq. (16) is difficult to do a priori, however:

• One natural next-to-leading contribution is obtained by taking into account second-order terms in ξ(d) in Eq. (16), that is, by considering doubles lines between cells in a diagrammatic representation. This would induce a term of about ξ(d)², whose average never vanishes⁵. As shown in the appendix, these contributions can be formally derived in the context of the hierarchical models. This leads to correction terms that can be organized in a sum of factorized terms. Therefore, although it can indeed provide corrective terms to the covariance matrix, only a limited number of eigendirections can be generated.

• Other contributions naturally come from proximity effects due to the fact that cells are finite, and could even overlap, which makes the expansion in ξ(d) ineffective. In a diagrammatic point of view, they are due to the fact that many more diagrams contributed when cells are too close. This has dramatic effects for overlapping cells. For hierarchical models, an approximate form can be used to help model these effects, which we use below.

• Finally, effects due to the fact that discrete tracers are used in count-in-cells statistics might also play a role at short distances. They are also tentatively modeled below.

In the following, we propose some modeling of these effects and explore how they depend on the properties of the survey.

3.5.1. Joint PDF at short distances

There are no general forms for the joint PDF at close distance. The hierarchical models suggest the following form (derived from the saddle point approximation, which is valid for moderate values of $\overline{\xi }$ and of the density contrast), however:

$$\begin{array}{cc}& P_{\mathrm{short}\mathrm{dist}.}({\rho }_i,{\rho }_j)\mathrm{d}{\rho }_i\mathrm{d}{\rho }_j=P({\rho }_m)\hfill \\ \hfill & \times exp[-\frac{{\delta }_{\rho }^2}{{\rho }_m^{\alpha }{\mathrm{\Delta }}_{\xi }(d)}]\frac{\mathrm{d}{\rho }_m\mathrm{d}{\delta }_{\rho }}{\sqrt{\pi {\rho }_m^{\alpha }{\mathrm{\Delta }}_{\xi }(d)}},\hfill \end{array}$$(38)

where ρ_m = (ρ_i + ρ_j)/2 and δ_ρ = (ρ_i − ρ_j)/2, and where α is model-dependent parameter. In other words, the PDF of the difference between ρ_i and ρ_j can be described by a simple Gaussian with a known width driven by the expression of ${\mathrm{\Delta }}_{\xi }(d)\equiv \overline{\xi }-{\xi }_{12}(d),$ provided it is small compared to $\overline{\xi }$. We note that Δ_ξ(d) obviously vanishes at d = 0, it then leads to a Kronecker δ function at zero separation as expected, and generically scales like d² at short distances⁶. Interestingly, for the minimal tree model the form (38) is exact for α = 1 (see appendix). In general, this is also the expected form based on the saddle point approximation (valid when $\overline{\xi }$ is small) for generic hierarchical tree models. The value of α can be identified from small-order cumulants,

$$\begin{array}{c}\hfill \alpha =\frac{2}{3}{S}_3,\end{array}$$(39)

where S₃ is the reduced third-order cumulant,

$$\begin{array}{c}\hfill S_3=\frac{⟨{\delta }^3⟩}{{⟨{\delta }^2⟩}^2}.\end{array}$$(40)

This form is probably not very accurate in general. It can be used to model the impact of close distances to the covariance matrix, however, as shown below.

3.5.2. Poisson noise and minimum separation

A further contribution to this joint PDF can come from discrete effects that arise because the density is evaluated from the counting of discrete tracers (as explored in Szapudi & Colombi 1996 or more recently in Repp & Szapudi 2021). In this subsection, we assume that the density corresponds to the density obtained after application of a top-hat filter and that tracers are Poisson realizations of continuous fields (although it is possible to encounter sub- or over-Poissonnian noises; Friedrich et al. 2018). The use of other filters can be explored but would require specific developments that we do not pursue here. Within such hypotheses then, the joint distribution of counts-in-cells N_i is given by the convolution of the joint density PDF, ℘({ρ_i}), in the continuous limit convolved by Poisson counts-in-cells as

$$\begin{array}{c}\hfill \mathcal{P}(\{N_i\})={\displaystyle \int {\mathrm{\Pi }}_i}\mathrm{d}{\rho }_i{P}_{\mathrm{Poisson}}(N_i;{\overline{N}}_i{\rho }_i)\mathcal{P}(\left\{{\rho }_i\right\}),\end{array}$$(41)

where P_Poisson(N;N̄) is more precisely the probability of having N tracers in a cell where the mean density of tracers is N̄. In practice, N̄_i is given by nV_i, where n is the number density of tracers and V_i is the volume of the cell V_i.

Discrete effects would then induce further scatter between the estimated values of ρ_i and ρ_j. The latter are given by Poisson noise induced by the nonoverlapping parts of the cells, as shown in Szapudi & Colombi (1996), further contributing to the scatter. The scatter in the difference in the number of points is

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{Poisson}}^2=\frac{2}{\overline{N}}{\rho }_m{f}_e(d).\end{array}$$(42)

It can be incorporated as a contribution to the variance of the PDF of δ_ρ of the form

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{Poisson}}^2=\frac{1}{2\overline{N}}{\rho }_m{f}_e(d),\end{array}$$(43)

where f_e(d) is the fraction of the volume of each cell that does not overlap with the other as a function of the cell distance. For short distances (i.e., for about d ≲ R), it is in the 2D case given by

$$\begin{array}{c}\hfill f_e(d)=\frac{2d}{\pi R}.\end{array}$$(44)

The expression (43) is then a priori to be added to the variance term that appears in Eq. (38) so that the total variance for the density difference reads

$$\begin{array}{c}\hfill {\sigma }_{{\delta }_{\rho }}^2(d)=\frac{1}{2}{\rho }_i^{\alpha }{\mathrm{\Delta }}_{\xi }(d)+\frac{1}{2\overline{N}}{\rho }_i{f}_e(d).\end{array}$$(45)

Equation (38) then fully encodes the fact that nearby cells are likely to have similar densities. This encodes, for instance, that nearby cells are within the same haloes. This contribution is expected to enhance the covariance terms. It shows that the amount of information is limited at small scales: there is therefore a minimum separation between cells smaller than which no gain in precision is expected of PDF measurements. The minimum distance depends on the bin size: d_min is the distance such that the densities in two cells separated by less than d_min are almost certainly in the sale density bin. d_min therefore depends on the bin width Δ_i. From Eq. (45), it is possible to infer this value. We desire

$$\begin{array}{c}\hfill {\sigma }_{{\delta }_{\rho }}^2(d_{min})\ll {\mathrm{\Delta }}_i^2.\end{array}$$(46)

This suggests that a minimum distance between cells can be adopted, given by

$$\begin{array}{c}\hfill d_{min\mathrm{Poisson}}=\frac{\pi {\mathrm{\Delta }}_i^2}{\overline{N}}R.\end{array}$$(47)

The other upper limit comes on d from the expression of Δ_ξ as a function of d. The latter depends on both the shape of the power spectrum and on the filter that is used. In general (e.g., for a Gaussian filter), Δ_ξ(d) scales like d²/R², where R is the filtering radius, with a coefficient c_ns that depends on the power spectrum index n_s and is proportional to its amplitude. Top-hat filters have different analytical properties. We give here the formal expression of Δ_ξ(d) at 2D for a power-law spectrum of index n_s,

$$\begin{array}{cc}\hfill \frac{{\mathrm{\Delta }}_{\xi }(d)}{\overline{\xi }}& =-\frac{2^{n_s-1}\mathrm{\Gamma }(1-\frac{n_s}{2})\mathrm{\Gamma }(2-\frac{n_s}{2})\mathrm{\Gamma }\left(\frac{1}{2}(n_s-1)\right)}{\sqrt{\pi }\mathrm{\Gamma }(\frac{1}{2}-\frac{n_s}{2})\mathrm{\Gamma }(\frac{3}{2}-\frac{n_s}{2})\mathrm{\Gamma }(\frac{n_s}{2}+1)}{\left(\frac{d}{R}\right)}^{1-n_s}\hfill \\ \hfill & \approx 0.72\frac{d^{3/2}}{R^{3/2}}\mathrm{for}{n}_s=0.5.\hfill \end{array}$$(48)

This is the situation we encounter below in the numerical experiments we perform. This leads to the following form:

$$\begin{array}{c}\hfill d_{min\mathrm{halo}}={\left(\frac{{\mathrm{\Delta }}_i^2}{0.72\overline{\xi }}\right)}^{2/3}R.\end{array}$$(49)

It is to be noted that it can be in practice a rather short distance, shorter than the filtering scale R. For instance, for a bin width of 1/4, a variance of about unity, d_minhalo is about R/5.

Equation (38), together with the expressions of the bias functions described above, is the main results of this paper. We illustrate below how they can be used to compute the covariance matrices.

4. Hierarchical models

In order to illustrate the previous findings, we make use of toy models for which explicit computations can be made.

4.1. General formalism

Hierarchical models are a class of non-Gaussian fields whose correlation functions follow the so-called hierarchical ansatz,

$$\begin{array}{c}\hfill {\xi }_p({\mathbf{r}}_1,\cdots ,{\mathbf{r}}_p)={\displaystyle \sum _{t\in \mathrm{trees}}}{Q}_p(t){\displaystyle \prod _{\mathrm{lines}\in t}}\xi ({\mathbf{r}}_i,{\mathbf{r}}_j),\end{array}$$(50)

where the sum is made over all possible trees that join the p points (diagram without loops), and the tree value is obtained by the product of a fixed weight (that depends only on the tree topology) and the product of the two-point correlation functions for all pairs that are connected together in the given tree. This construction ensures that the average p-point function, $\overline{{\xi }_p}$, scales like the ${\overline{\xi }}^{p-2}$, where $\overline{\xi }$ is the averaged two-point function. More precisely, there are S_p parameters such that

$$\begin{array}{c}\hfill \overline{{\xi }_p}=S_p{\overline{\xi }}^{p-2}.\end{array}$$(51)

The precise value of the S_p parameters depend on the Q_p parameters and on the averages of the product of ξ(r_ij) functions. A very good approximation is to assume that the average of the products of this function is given by the product of these averages. Then the S_p coefficients depend solely on Q_p,

$$\begin{array}{c}\hfill S_p={\displaystyle \sum _t}{Q}_p(t).\end{array}$$(52)

4.2. The (minimal) tree model

The tree models are based on a further assumption on the Q_p parameters. It is basically assumed that tree expressions can be computed locally⁷, that is,

$$\begin{array}{c}\hfill Q(t)={\displaystyle \prod _{\mathrm{vertices}\in t}}{\nu }_p,\end{array}$$(53)

where ν_p is a weight attributed to all vertices with p incoming lines (ν₀ = ν₁ = 1 for completion). In this formalism, the vertex generating function is generally introduced,

$$\begin{array}{c}\hfill \zeta (\tau )={\displaystyle \sum _p}\frac{{\nu }_p}{p!}{\tau }^p.\end{array}$$(54)

The minimal tree model is a model in which ν₂ alone does not vanish. In the minimal model⁸, its value is fixed and is given by ν₂ = 1/2, so that

$$\begin{array}{c}\hfill {\zeta }_{\mathrm{RL}}(\tau )={(1+\tau /2)}^2.\end{array}$$(55)

Together with the Gaussian case (which corresponds to ζ(τ)=1 + τ), this is the only case for which we are sure that it can be effectively built (in the sense that other models may be unphysical).

In this model, it is possible to build the cumulant generating function of the local density. For the one-point case, assuming the mean-field approximation, the CGF is given by

$$\begin{array}{c}\hfill \phi (\lambda )=\lambda [\zeta (\tau )-\frac{1}{2}\tau {\zeta }^{\prime }(\tau )]\end{array}$$(56)

with

$$\begin{array}{c}\hfill \tau =\lambda \overline{\xi }{\zeta }^{\prime }(\tau ).\end{array}$$(57)

This is not the result of large deviation principle calculations, but of mere combinatorics, although it leads to the same formal transformation between the CGF and the vertex-generating function. In case of the minimal model, Eq. (57) takes a simple form that can be easily solved. We finally have

$$\begin{array}{c}\hfill \phi (\lambda )=\frac{\tau (\lambda )}{\overline{\xi }},\tau (\lambda )=\frac{\lambda \overline{\xi }}{1-\lambda \overline{\xi }/2}.\end{array}$$(58)

The one-point PDF of the density can then be computed explicitly (see appendix),

$$\begin{array}{c}\hfill P(\rho )=\frac{4}{{\overline{\xi }}^2}exp[-\frac{2}{\overline{\xi }}(1+\rho )]{}_0{F}_1(2,\frac{4\rho }{\overline{\xi }}),\end{array}$$(59)

as can the density-bias function,

$$\begin{array}{c}\hfill b(\rho )=\frac{{}_0{F}_1(1,\frac{4\rho }{\overline{\xi }})}{{}_0{F}_1(2,\frac{4\rho }{\overline{\xi }})}-\frac{2}{\overline{\xi }},\end{array}$$(60)

where $\overline{\xi }$ is the averaged two-point correlation function within the cell.

4.3. Rayleigh-Levy flight model

The minimal tree model can be implemented with Rayleigh-Levy random walks (or rather Rayleigh-Levy flights, as described in Peebles 1980). This is a Markov random walk where the PDF of the step length ℓ follows a power law,

$$\begin{array}{c}\hfill P(\ell )\sim \frac{1}{{\ell }^{\alpha }},\end{array}$$(61)

with a regularizing cutoff at small separation, and where α satisfies

$$\begin{array}{c}\hfill 0<\alpha <2.\end{array}$$(62)

The sample points are then all the step points reached by the walker.

More precisely, the cumulative distribution function of step of length ℓ is

$$\begin{array}{cc}& P(>{\ell }_0)=1,\hfill \end{array}$$(63)

$$\begin{array}{cc}& P(>\ell )={\left(\frac{{\ell }_0}{\ell }\right)}^{\alpha }\mathrm{for}\ell >{\ell }_0,\hfill \end{array}$$(64)

where ℓ₀ is a small-scale parameter. The two- and higher-order correlation functions can then be explicitly computed. Starting with a first point at position r₀, the density of the subsequent point (first descendant) at position r is given by

$$\begin{array}{cc}& f_1(\mathbf{r})=\frac{\alpha }{2\pi }\frac{{\ell }_0^{\alpha }}{|\mathbf{r}-{\mathbf{r}}_0{|}^{2+\alpha }}\mathrm{in}2\mathrm{D}\mathrm{space};\hfill \end{array}$$(65)

$$\begin{array}{cc}& f_1(\mathbf{r})=\frac{\alpha }{4\pi }\frac{{\ell }_0^{\alpha }}{|\mathbf{r}-{\mathbf{r}}_0{|}^{3+\alpha }}\mathrm{in}3\mathrm{D}\mathrm{space}.\hfill \end{array}$$(66)

In the following, the dimension of space is denoted D. The density of the descendants, assuming there are an infinity of them, of a point at position r₀ is then given by a series of convolutions,

$$\begin{array}{c}\hfill f({\mathbf{r}}_0,\mathbf{r})=f_1(\mathbf{r})+{\displaystyle \int {\mathrm{d}}^D}{\mathbf{r}}_1{f}_1(\mathbf{r}-{\mathbf{r}}_1)f_1({\mathbf{r}}_1-{\mathbf{r}}_0)+\cdots ,\end{array}$$(67)

with subsequent convolutions and where the integral is done in the whole space. Defining the Fourier transform of f₁(r) as ψ(k),

$$\begin{array}{c}\hfill \psi (k)={\displaystyle \int {\mathrm{d}}^D}\mathbf{r}{f}_1(\mathbf{r})e^{-\mathrm{i}\mathbf{k}.\mathbf{r}},\end{array}$$(68)

which is then a function of k only, it is easy to see that

$$\begin{array}{cc}& f({\mathbf{r}}_0,\mathbf{r})={\displaystyle \int \frac{{\mathrm{d}}^D\mathbf{k}}{{(2\pi )}^D}}{e}^{\mathrm{i}\mathbf{k}.(\mathbf{r}-{\mathbf{r}}_0)}[\psi (k)+{\psi }^2(k)+\cdots ]\hfill \\ \hfill & ={\displaystyle \int \frac{{\mathrm{d}}^D\mathbf{k}}{{(2\pi )}^D}}{e}^{\mathrm{i}\mathbf{k}.(\mathbf{r}-{\mathbf{r}}_0)}\frac{1}{1-\psi (k)},\hfill \end{array}$$(69)

where we take advantage of the expression of convolutions in Fourier space and their resummations. The two-point correlation function is then given by two possible configurations: a neighbor can either be an ascendant or a descendant, so that the two-point correlation functions between positions r₁ and r₂ are given by

$$\begin{array}{cc}& {\xi }_2({\mathbf{r}}_1,{\mathbf{r}}_2)=\frac{1}{n}[f({\mathbf{r}}_1,{\mathbf{r}}_2)+f({\mathbf{r}}_2,{\mathbf{r}}_1)]\hfill \\ \hfill & =\frac{1}{n}{\displaystyle \int \frac{{\mathrm{d}}^D\mathbf{k}}{{(2\pi )}^D}}{e}^{\mathrm{i}\mathbf{k}.({\mathbf{r}}_2-{\mathbf{r}}_1)}\frac{2}{1-\psi (k)},\hfill \end{array}$$(70)

where n is the number density of points in the sample that can be associated with a typical length ℓ_n,

$$\begin{array}{c}\hfill n=\frac{1}{{\ell }_n^D}.\end{array}$$(71)

At large scale, this expression causes the power spectra to be power laws. They scale like k^−α, and the resulting two-point correlation function then takes the form in the large separation limit,

$$\begin{array}{cc}& {\xi }_{\alpha ,2D}(r)=\frac{\alpha }{\pi }{r}^{\alpha -2}{\ell }_0^{-\alpha }{\ell }_n^2\hfill \end{array}$$(72)

$$\begin{array}{cc}& {\xi }_{\alpha ,3D}(r)=\frac{(1-{\alpha }^2)tan\left(\frac{\pi \alpha }{2}\right)}{{\pi }^2}{r}^{\alpha -3}{\ell }_0^{-\alpha }{\ell }_n^3.\hfill \end{array}$$(73)

It is to be noted, however, that this expression does not take into account the boundary conditions, in particular if they are assumed to be periodic. This case is examined in some detail in the next paragraph. It is to be noted, however, that in this case, the function ξ(r) has a more complex form. It is in particular no more isotropic.

Higher-order correlation functions can also be computed in this model: n points are correlated when they are embedded in a chronological sequence that can be run in one direction or the other. Thus the three-point function is simply given by

$$\begin{array}{c}\hfill {\xi }_{\alpha }({\mathbf{r}}_1,{\mathbf{r}}_2,{\mathbf{r}}_3)=\frac{1}{n^2}[f({\mathbf{r}}_1,{\mathbf{r}}_2)f({\mathbf{r}}_2,{\mathbf{r}}_3)+\cdots ],\end{array}$$(74)

with five other terms obtained by all combinations of the indices. Expressing the result in terms of the two-point function, we have

$$\begin{array}{cc}& {\xi }_{\alpha }({\mathbf{r}}_1,{\mathbf{r}}_2,{\mathbf{r}}_3)=\frac{1}{2}[{\xi }_{\alpha }({\mathbf{r}}_1,{\mathbf{r}}_2){\xi }_{\alpha }({\mathbf{r}}_2,{\mathbf{r}}_3)+\hfill \\ \hfill & {\xi }_{\alpha }({\mathbf{r}}_2,{\mathbf{r}}_3){\xi }_{\alpha }({\mathbf{r}}_3,{\mathbf{r}}_1)+{\xi }_{\alpha }({\mathbf{r}}_3,{\mathbf{r}}_1){\xi }_{\alpha }({\mathbf{r}}_1,{\mathbf{r}}_2)],\hfill \end{array}$$(75)

corresponding to a tree structure with ν₂ = 1/2.

Higher-order correlation functions can be computed similarly. They follow a tree structure in the sense above, with ν₂ = 1/2 and ν_p = 0 for p ≥ 3.

4.4. Periodic boundary conditions

We briefly explore here the case of periodic boundary conditions. Then the multipoint density field g^PBC(r_i) for periodic boundary conditions can be expressed in terms of the former density field g(r_i) as sums of all copies of the sample, that is,

$$\begin{array}{c}\hfill g^{\mathrm{PBC}}(\{{\mathbf{r}}_i\})={\displaystyle \sum _{{\mathbf{n}}_i}}g(\{{\mathbf{r}}_i+{\mathbf{n}}_{i}L\}),\end{array}$$(76)

where n_i are vectors whose components are integers, ${\mathbf{n}}_i=(n_i^x,n_i^{\mathit{y}},n_i^z)$ and the sums run over all integer values for all i; L is the size of the sample (assumed to be the same in all directions).

When it is applied in this context, we can construct the n-point density function out of the density function f computed previously. Thus the two-point density function is given by

$$\begin{array}{c}\hfill g^{\mathrm{PBC}}({\mathbf{r}}_1,{\mathbf{r}}_2)=n^{\mathrm{PBC}}{\displaystyle \sum _{{\mathbf{n}}_{12}}}f({\mathbf{r}}_1,{\mathbf{r}}_2-{\mathbf{r}}_1+{\mathbf{n}}_{12}L),\end{array}$$(77)

where n₁₂ = n₂ − n₁ and n^PBC is the resulting one-point (and therefore homogeneous) density in the sample. This expression is therefore written in terms of the function

$$\begin{array}{c}\hfill f^{\mathrm{PBC}}({\mathbf{r}}_0,\mathbf{r})={\displaystyle \sum _{\mathbf{n}}}f({\mathbf{r}}_0,\mathbf{r}-{\mathbf{r}}_0+\mathbf{n}L).\end{array}$$(78)

We can now compute its expression in terms of the power spectra, or more specifically, the function ψ(k) defined previously. We have

$$\begin{array}{c}\hfill f^{\mathrm{PBC}}({\mathbf{r}}_0,\mathbf{r})={\displaystyle \int \frac{{\mathrm{d}}^D\mathbf{k}}{{(2\pi )}^D}}{e}^{\mathrm{i}\mathbf{k}.(\mathbf{r}-{\mathbf{r}}_0)}{\displaystyle \sum _{\mathbf{n}}}{e}^{\mathrm{i}\mathbf{n}.\mathbf{k}L}\frac{1}{1-\psi (k)},\end{array}$$(79)

and the latter sum ensures that the contributing wave modes k are only those that are periodic in all three directions, that is, those whose components are multiples of 2π/L so that

$$\begin{array}{c}\hfill f^{\mathrm{PBC}}({\mathbf{r}}_0,\mathbf{r})={\displaystyle \sum _{{\mathbf{n}}^{\ast }}}\frac{1}{L^D}{e}^{2\pi \mathrm{i}\mathbf{n}.(\mathbf{r}-{\mathbf{r}}_0)/L}\frac{1}{1-\psi (k_n)},\end{array}$$(80)

with

$$\begin{array}{c}\hfill k_n={(\mathbf{n}.\mathbf{n})}^{1/2}\frac{2\pi }{L},\end{array}$$(81)

and where the sum is all over possible integer triplets but n = (0, 0, 0). The two-point correlation function is now given by

$$\begin{array}{c}\hfill {\xi }_{\alpha }^{\mathrm{PBC}}({\mathbf{r}}_1,{\mathbf{r}}_2)=\frac{1}{n^{\mathrm{PBC}}}[f^{\mathrm{PBC}}({\mathbf{r}}_1,{\mathbf{r}}_2)+f({\mathbf{r}}_2,{\mathbf{r}}_1)].\end{array}$$(82)

A similar result can be obtained for the three-point correlation function with

$$\begin{array}{c}\hfill {\xi }_{\alpha }^{\mathrm{PBC}}({\mathbf{r}}_1,{\mathbf{r}}_2,{\mathbf{r}}_3)=\frac{1}{{(n^{\mathrm{PBC}})}^2}[f^{\mathrm{PBC}}({\mathbf{r}}_1,{\mathbf{r}}_2)f^{\mathrm{PBC}}({\mathbf{r}}_2,{\mathbf{r}}_3)+\cdots ].\end{array}$$(83)

As a consequence, the functional relation between the three-point correlation function and the two-point correlation function is left unchanged. This is also the case at all orders.

4.5. Covariance matrix of the minimal tree model

Remarkably, in case of the minimal tree model, the derivation of the CGF can also be made for multiple cells, and in particular, for two cells. Its expression is derived in the appendix. It takes the form

$$\begin{array}{c}\hfill \phi ({\lambda }_1,{\lambda }_2)=\frac{{\lambda }_1+{\lambda }_2+({\xi }_{12}-\overline{\xi }){\lambda }_1{\lambda }_2}{1-({\lambda }_1+{\lambda }_2)\overline{\xi }/2-{\lambda }_1{\lambda }_2({\xi }_{12}^2-{\overline{\xi }}^2)/4}.\end{array}$$(84)

In this case, it is then possible to expand its expression in powers of ξ₁₂ for distant cells or in powers of $(\overline{\xi }-{\xi }_{12})$ for close cells, and in both cases, closed forms can be obtained to any order. It leads to the possibility of computing the joint PDF for any configuration (see the appendix for details) and finally to evaluate the covariance matrix directly. This is even possible for any of the thre sets of variables we consider, {ρ_i}, $\{{\widehat{\rho }}_i\}$, or $\{{\overline{\rho }}_i\}$.

We performed these computations for the minimal tree model with a power-law behavior ξ(r)∼r^−1.5 (α = 0.5), a 2D survey with a size of 200² pixels, and a top-hat smoothing radius of 4.25 pixels. The amplitude of the correlation function was fixed to give $\overline{\xi }=1.09$ at the smoothing scale. It precisely corresponds to the setting of the numerical simulations of Rayleigh-Levy flights we also performed, as described in the next section. It allows us to compare the two approaches. These analytic results have two limitations: the results are based on the mean-field approximation for the computation of the two-variable GCF, and the covariance elements are computed ignoring the bin sizes (i.e., by evaluating the expression of the covariance for their central values). Although in most cases this should not be an issue, it still might have a non-negligible impact when the PDF varies rapidly with the density.

5. Simplified models of the covariance matrix

The purpose of this section is then to propose two levels of modeling of the covariance matrix based on the previous results and to compare these propositions with results of either the full analytic results presented before or with the results of numerical experiments based on Rayleigh-Levy flights.

5.1. Modeling the covariance matrix

More specifically, we considered two approximate forms for the full covariance. The first approximation is fully analytic. It makes use of the large-scale contributions and those from the short distance expression (38). It reads as the sum of the two contributions

$$\begin{array}{cc}& \mathrm{Cov}({\rho }_i,{\rho }_j)=b_{\#}({\rho }_i)P({\rho }_i){\overline{\xi }}_s{b}_{\#}({\rho }_j)P({\rho }_j)\hfill \\ \hfill & +{\displaystyle {\int }_0^{r_{max}}}\mathrm{d}{r}_d{P}_s(r_d)P_{\mathrm{short}\mathrm{dist}.}({\rho }_i,{\rho }_j,r_d).\hfill \end{array}$$(85)

In this expression, the only free parameter is r_max. This is indeed a crucial parameter as it determines to a large extent the amplitude of the short-distance effects. In the following, we take r_max = R, that is, the filtering scale. It is found to give a good result for the 2D case and for n_s = −1/2, but this choice is likely to depend on the shape of the power spectrum. In general, this formula is intended to give a good account of the general properties of the covariance matrix, it cannot provide reliable quantitative results a priori.

The other form we propose is intended to be much more precise quantitatively. Is is given by the following expression:

$$\begin{array}{c}\hfill \mathrm{Cov}({\rho }_i,{\rho }_j)=b_{\#}({\rho }_i)P({\rho }_i){\overline{\xi }}_s{b}_{\#}({\rho }_j)P({\rho }_j)+{\mathrm{Cov}}^{\mathrm{PBC}}({\rho }_i,{\rho }_j),\end{array}$$(86)

where Cov^PBC(ρ_i, ρ_j) is the expression of the covariant matrix for periodic boundary conditions. It is obtained here simply by replacing $P({\rho }_i,{\rho }_j,\overline{\xi },{\xi }_{12}(r_d))$ by $P({\rho }_i,{\rho }_j,\overline{\xi }-{\overline{\xi }}_s,{\xi }_{12}(r_d)-{\overline{\xi }}_s)$ before integrating over r_d so that the averaged joint correlations vanish identically. The rationale for this proposition is that Cov^PBC(ρ_i, ρ_j) could be more easily estimated from specific numerical experiments. In both cases, the short-distance contributions are the same for the three types of observables ρ_i, ${\widehat{\rho }}_i$, and ${\overline{\rho }}_i$. These forms are then compared to numerical results.

5.2. Numerical experiments with the Rayleigh-Levy flight model

A series of experiments of 2D walks with a large number of samples were performed as described below. We restricted our analysis to α = 0.5 with l₀ = 0.003 pixel size (the dependence on l₀ was tested as illustrated on Fig. 2, where l₀ = .006 was also used, but the analyses were made for a fixed value of l₀). Figure 1 illustrates how points are distributed in these samples. The point distribution does not show the filamentary structure of realistic cosmological simulations. It exhibits the presence of concentrated halos surrounded by empty regions, however, which are reminiscent of the structure of the largest matter concentrations of the cosmic web.

Two different setting were employed to explore different aspects of the results that were found:

• Set 𝒜: 1600 samples extracted from a single numerical realization (with periodic boundary conditions) with a size of 8000 × 8000 pixels² containing 64 × 10⁶ points. Each sample then has 200 × 200 pixels² containing an average of 200² points each. For this set of samples, the average and covariance of the PDF were extracted following the three procedures mentioned before: either the density was taken with respect to the mean density of the realization, with respect to the density of each sample, or by subtracting the sample density. It therefore corresponds to an evaluation of the mean and covariance of the PDF of ρ_i, ${\widehat{\rho }}_i$, and ${\overline{\rho }}_i$, respectively.

• Set 𝓑: 1600 samples, each with periodic boundary conditions, with a size of 200 × 200 pixels² containing 200² points each. By construction, the average two-point function in the sample, ${\overline{\xi }}_s$, vanishes in this case, and covariance is entirely due to proximity effects.

In each case, the local density was obtained after a filtering procedure. The point positions were first pixelized, that is, each point was attributed to a pixel so that the mean number of points per pixel was one. The field was then filtered by a (quasi) circular top-hat functions. In practice, the number of pixels in the window function was 57. This makes the effective smoothing radius about 4.25 in pixel units. The resulting density was then measured at each pixel location. Their histograms were then computed after density binning. To avoid large undue discrete effects, the bin width was chosen to be a multiple of 1/57, and in order to ensure that the requirement (49) was met at the pixel distance, we chose a bin width of about 1/4, more precisely, of 14/57.

Figure 2 shows the resulting PDF as measured in the simulations and how it compares to the theoretical prediction, Eq. (59), for two different choices of l₀. The expected scaling for $\overline{\xi }$ is recovered. The measured PDF also follows theoretical predictions for a wide range of probabilities remarkably well. It gives us confidence in the whole procedure and in the approaches used to compute PDFs in this model. The detailed comparisons were made for l₀ = 0.003, leading to $\overline{\xi }=1.09,$ and a sample density variance in sets 𝒜 given by ${\overline{\xi }}_s=0.09$.

The measured variance of the density PDF is obtained from 1600 samples in each case. The resulting shapes are presented in Figs. 3–5 for the different cases, density in a supersample realization, and in samples with periodic boundary conditions. The results show the comparison between results obtained in the numerical simulations with yellow symbols, and results derived from the analytic prescriptions as blue dots, based on the mean-field approximation. The agreement between the two is very good. The overall shape of the variance and its dependence on the density is well reproduced. Discrepancies can be observed for densities above 4 or 5, however, where the theoretical predictions are seen to underestimate the results. The reasons for these discrepancies are not clear at this stage. A possible explanation might be the finite number of samples that is used to infer the variances⁹. The variance of the density PDF is also compared with the large-scale contributions (19), (29), and (34) for set 𝒜 depending on the cases (at this order, the covariance vanishes for set 𝒝). It shows that this formula captures some features of the variance (especially at low and moderate densities), but does not account for all. This is also illustrated in Fig. 4, which shows the reduced covariance. The fact that the covariance is determined to a large extent by its leading large-scale contribution leads to values of the reduced covariance close to 1 or −1, leading to these butterfly patterns. Proximity effects, not captured in these forms, also contribute to the covariances at a significant level, however. This is already apparent in Fig. 3.

5.3. Testing models of covariance matrices

Expressions (85) and (86) are precise propositions to show how the large-scale contributions can be completed to account for the full form of the covariance. The comparisons between the predicted form and those obtained from the numerical experiments are explored in detail at different levels and using the following criteria:

• Amplitude of the PDF variance,

• Density dependence of the first eigenvalue of the covariance matrix,

• Amplitude of the eigenvalues of the covariance matrix, and

• Resulting χ² distribution of a set of data vectors drawn from the original covariance.

These comparisons are shown in Figs. 6–8. For model (86), the term Cov^PBC(ρ_i, ρ_j) is taken from the measured covariance of set 𝒝. Figures 6 and 7 show that these two prescriptions give a good account of the leading behavior of the covariance matrix. The conclusion is quite sensitive for the choice of r_max for prescription (85). On the other hand, there is no free parameter that can be adjusted for model (86). Interestingly, Fig. 7 shows that the PDF variance also departs significantly from the large-scale term. The first eigenvector reproduces the functional form of the large-scale density-bias functions very faithfully.

The last two criteria are designed to verify that the reconstructed covariances also capture the subleading behavior of the matrix and can eventually be safely inverted and used as a model of likelihood. To avoid numerical uncertainties and make the comparison tractable, we chose to reduce the binning to six bins (through a rebinning of the histograms and densities ranging from 0.5 to 6.5). The resulting eigenvalues are shown in the top panel of Fig. 8. It shows that the eigenvalues decrease rapidly in amplitude, suggesting that the eigendirections are well sequenced and that the approximate form captures their values rather accurately. Form (86) in particular reproduces all six eigenvalues almost exactly.

Finally, χ² distributions were computed from a set of random values $P_i^{\mathrm{e}\mathrm{x}}$ drawn in each case from a Gaussian likelihood built from the measured covariance (with six bins). The values of ${\chi }^2(P_i^{\mathrm{e}\mathrm{x}})$ were then computed for each data vector, and their histogram was computed from each of the proposed models (including the original model for reference),

$$\begin{array}{c}\hfill {\chi }_{\mathrm{model}}^2(P_i^{\mathrm{ex}})=\frac{1}{2}{\displaystyle \sum _{\mathit{ij}}}{\mathcal{N}}_{\mathrm{model}}^{\mathit{ij}}{P}_i^{\mathrm{ex}}{P}_j^{\mathrm{ex}},\end{array}$$(87)

where ${\mathcal{N}}_{\text{model}}^{\mathit{ij}}$ is the inverse of the covariance matrix, either computed from Eq. (85) or from Eq. (86). For the original model, the expected distribution of the χ² values is then expected to be precisely that of a χ² distribution with six degrees of freedom. This is indeed what is almost exactly obtained for model (86). Results obtained from prescription (85) are not quite as good. This is expected as the short-distance effects are estimated rather crudely in Eq. (85). The performance of this prescription deteriorates when the dimension of the data vector (i.e., the number of bins) increases.

6. Conclusions and lessons

We presented key relations that give the large-scale behavior of the joint PDF, and hence the leading behavior of the covariance matrix of the density PDF. These contributing terms do not give the expression of a covariance matrix that can be used to build a likelihood function, however, as it is not invertible. Further significant contributions are found to be due to small separation effects, and an approximate form is proposed in Eq. (38). The latter is found to encapsulate most of the proximity effects, that is, it informs on the fact that nearby regions are likely to be correlated. They also give an indication on the minimal grid size that can be used the maximum bin size that can be used without information loss for a given bin width.

We then used a toy model for which numerical experiments can easily be performed and for which the exact PDF and large-scale covariance can be derived. It allows us to evaluate the efficiency of approximate schemes precisely. The conclusions of these comparisons are listed below.

• The theoretical forms Eqs. (19), (29), and (34) give the leading-order expression of the covariance elements when supersample effects are taken into account. It gives an accurate prediction of the leading eigenvalue and eigendirection of the covariance matrix.

• Whether subdominant effects can be accounted for by subsequent terms depends on the behavior of the two-point function: if the rms of the two-point function is dominated by large separations, then next-to-leading-order effect need to be taken into account; otherwise, short-distance effects will be the dominant contributor.

• In case short-distance effects dominate, the covariance matrix can be accessed from small simulations provided the relevant dominant large-scale contributions are added.

• This suggests that in realistic situations, the supersample effects, that is, the effects due to modes whose wavelength is larger than the size of the survey, have limited impact on the structure of the covariance matrix and that they can be captured by the only leading large-scale contribution. This is supported by a further analysis of the behavior of ξ(r_d) in realistic cosmological settings. For the standard model of cosmology (as derived from cosmic microwave background observations, Planck Collaboration VI 2020), the behavior of the matter correlation function can be derived. This is illustrated in Fig. 9, which illustrates the scales that are the main contributors to the first two moments of the two-point correlation function. Whether in 2D or in 3D, the first moment is dominated by large-scale contributions, whereas the second moment is dominated by small-scale contributions.

  Fig. 9. Scale dependence of the matter correlation functions for a realistic cosmological model (cosmological parameters derived from Plank, Planck Collaboration VI 2020) for the 3D density and the projected density (for a uniformly sampled survey with a depth of about 800 h−1Mpc between z = 0.75 and z = 1.25). The top panel shows $r_d^3\xi (r_d)$ (solid blue line) and $r_d^3{\xi }^2(r_d)$ (dashed red line) for the 3D density field, and the bottom panel shows ${\theta }_d^2\xi ({\theta }_d)$ and ${\theta }_d^2{\xi }^2({\theta }_d)$ for the projected density. In both cases, the average value of the first moment of the two-point correlation function is dominated by large-distance contributions, whereas short-distance contributions dominate the second moment, assuming survey sizes of about 100 h−1 Mpc or above.

• In the context of this study, we assumed that the measured P_i were Gaussian distributed. Although it is difficult to assess the accuracy of this hypothesis, the structure uncovered in Sect. 3 can be used to make such an attempt. In tree models, higher-order expressions of the joint density PDFs are expected to preserve the tree structure; see Bernardeau & Schaeffer (1999). The connected part of the three-point joint density PDF is then expected to take the form

  $$\begin{array}{cc}& \mathrm{Cov}({\rho }_i,{\rho }_j,{\rho }_k)\hfill \\ \hfill & =b_2({\rho }_i)P({\rho }_i){\overline{\xi }}_{s}b({\rho }_j)P({\rho }_j){\overline{\xi }}_{s}b({\rho }_k)P({\rho }_k)+\mathrm{sym}.,\hfill \end{array}$$(88)

  where b₂(ρ) is the two-line bias function of amplitude similar to b²(ρ). This implies in particular that the third-order cumulant is about $b{(\rho )}^{4}P{(\rho )}^3{{\overline{\xi }}_s}^2$, much smaller than ${\left[b{(\rho )}^{2}P{(\rho )}^2{\overline{\xi }}_s\right]}^{3/2}$, making the distribution of the measured values of P(ρ) (quasi-) Gaussian distributed. There might be some combination of ρ_i and values of ${\overline{\xi }}_s$, however, for which a higher-order term could play a role in the expression of the likelihood function.

For the application of these formulae in practical cases, some limitations have to be noted. We list them below.

• In the proposed form, the fact that in practice, PDFs are generally measured on a grid, that is, on a finite set of locations, is not taken into account. For instance, the exclusion of nonoverlapping cells is not considered. this is expected to introduce additional noise in the PDF estimates. The covariance matrix for these constructions cannot then be derived from general formulae (6) even when the integral in r_d is restricted above a given threshold.

• Relation (38) has been derived in a specific regime (using saddle point approximations) for tree hierarchical models. They are expected to capture the phenomenon at play for “typical” values of the densities, but they may not perform so well in the rare event tails (the exception being the minimum model, for which it is exact). Further checks of the validity of (38) should therefore certainly be done.

• The general formulae (19), (29), and (34) are valid for any type of filtering schemes, even for a compensated filter. This is not the case for relation (38). The proximity effects for compensated filters ought to be considered specifically.

• Prescription (86) is found to give a very precise account of the properties of the covariance matrix. It is based on the proposition that large-scale (supersample) effects can be added separately from the proximity effects and that the latter can be evaluated with small-scale mocks in which supersample effects are absent (with periodic boundary conditions). This is not an exact result, however,. It relies in particular on the fact that the r.m.s. of the ξ_s is dominated by scales much smaller than the sample size.

• Prescription (85) is less solid. It can be used for a quick assessment of the different contributing terms, or to build fully invertible covariance matrices, but it is unlikely to give reliable predictions at the χ² level.

In all cases, prescriptions (85) and (86) can be the starting point of a more precise evaluation of the covariance from specific numerical experiments that can complement its evaluation following the approach presented in Friedrich & Eifler (2018), for instance. The authors also showed that some strategies could be adopted to limit the number of realizations required to reach a specific precision. This point is not discussed here.

---

Acknowledgments

The author of this article is indebted to Cora Uhlemann, Alex Gough, Oliver Friedrich, Sandrine Codis, Aoife Boyle and Alexandre Barthelemy for many comments and careful examination of the preparatory notes of this manuscript.

References

Appendix A: Hierarchical tree models

In hierarchical tree models, the p-point matter correlation functions are assumed to follow tree structures in the sense described in the main text. They are thus entirely defined by the two-point functions ξ(r) and the vertex-generating function ζ(τ). The exact generating function of multiple cell correlation functions can be built through simple transforms. We therefore consider a set of n cells V_i. These cells can overlap. The joint cumulants we consider are those of the average densities in cells V_i that can be expressed in terms of spatial averages¹⁰ of correlation functions,

$$\begin{array}{cc}\hfill {⟨{\rho }_1^{p_1}\cdots {\rho }_n^{p_n}⟩}_c& ={\displaystyle {\int }_{V_1}}\frac{\mathrm{d}{\mathbf{x}}_{1,1}}{V_1}\cdots {\displaystyle {\int }_{V_1}}\frac{\mathrm{d}{\mathbf{x}}_{1,p_1}}{V_1}\cdots \hfill \\ \hfill & ...{\displaystyle {\int }_{V_n}}\frac{\mathrm{d}{\mathbf{x}}_{n,1}}{V_n}\cdots {\displaystyle {\int }_{V_n}}\frac{\mathrm{d}{\mathbf{x}}_{n,p_n}}{V_n}\hfill \\ \hfill & \times {\xi }_p({\mathbf{x}}_{1,1},\cdots ,{\mathbf{x}}_{1,p_1},\cdots {\mathbf{x}}_{n,1},\cdots ,{\mathbf{x}}_{n,p_n}).\hfill \end{array}$$(A.1)

We then wish to build the cumulant-generating function,

$$\begin{array}{c}\hfill \phi ({\lambda }_1,\cdots ,{\lambda }_n)={\displaystyle \sum _{p_1,\cdots ,p_n}}{⟨{\rho }_1^{p_1}\cdots {\rho }_n^{p_n}⟩}_c\frac{{\lambda }_1^{p_1}}{p_1!}\cdots \frac{{\lambda }_n}{p_n!}.\end{array}$$(A.2)

This function represents the generating function of (averaged) tree diagrams where λ_i counts the number of points in each cells. As shown in Bernardeau & Schaeffer (1999), this is obtained with the help of the intermediate function τ(x) solution of the consistency equation¹¹,

$$\begin{array}{c}\hfill \tau (\mathbf{x})={\displaystyle \sum _j}{\lambda }_j{\displaystyle {\int }_{V_j}}\frac{\mathrm{d}\mathbf{x}}{V_j}\xi (\mathbf{x},{\mathbf{x}}^{\prime }){\zeta }^{\prime }(\tau ({\mathbf{x}}^{\prime })),\end{array}$$(A.3)

and then

$$\begin{array}{cc}\hfill \phi ({\lambda }_1,\cdots ,{\lambda }_n)& ={\displaystyle \sum _j}{\lambda }_j{\displaystyle {\int }_{V_j}}\frac{\mathrm{d}\mathbf{x}}{V_j}\zeta (\tau (\mathbf{x}))\hfill \\ \hfill & -\frac{1}{2}{\displaystyle \sum _j}{\lambda }_j{\displaystyle {\int }_{V_j}}\frac{\mathrm{d}\mathbf{x}}{V_j}\tau (\mathbf{x}){\zeta }^{\prime }(\tau (\mathbf{x})).\hfill \end{array}$$(A.4)

This is an exact result based on pure combinatorics.

For cases of interest, it is possible to do a mean-field approximation that consists of assuming that τ(x) is constant within each cell. We then have the system of equations for τ_i,

$$\begin{array}{c}\hfill {\tau }_i={\displaystyle \sum _j}{\lambda }_j{\xi }_{\mathit{ij}}{\zeta }^{\prime }({\tau }_j),\end{array}$$(A.5)

where

$$\begin{array}{c}\hfill {\xi }_{\mathit{ij}}={\displaystyle {\int }_{V_i}}\frac{\mathrm{d}{\mathbf{x}}_i}{V_i}{\displaystyle {\int }_{V_i}}\frac{\mathrm{d}{\mathbf{x}}_i}{V_i}\xi ({\mathbf{x}}_i,{\mathbf{x}}_j)\end{array}$$(A.6)

and

$$\begin{array}{c}\hfill \phi ({\lambda }_1,\cdots ,{\lambda }_n)={\displaystyle \sum _j}{\lambda }_j[\zeta ({\tau }_j)-\frac{1}{2}{\tau }_j{\zeta }^{\prime }({\tau }_j)].\end{array}$$(A.7)

Bernardeau & Schaeffer (1999) found this to be very accurate, and we extensively use this approximation in the following, in particular for the minimal tree model.

Appendix B: Joint PDF, density-bias function in PT, and hierarchical tree models

Here we consider the joint distribution of densities in two cells whose centers are at distance d. The calculation is based on the inverse Laplace transform of the joint cumulant-generating function φ(λ₁, λ₂),

$$\begin{array}{c}\hfill \phi ({\lambda }_1,{\lambda }_2)={\displaystyle \sum _{p,q}}{⟨{\rho }_1^p{\rho }_2^q⟩}_c\frac{{\lambda }_1^p}{p!}\frac{{\lambda }_2^q}{q!},\end{array}$$(B.1)

where ${⟨{\rho }_1^p{\rho }_2^q⟩}_c$ are the cumulants of the local density fields. They depend on the size and distance d between the cells. We assume in particular that the density correlation function between cells ξ(d) is small compared to unity and can serve as a small parameter.

B.1. Leading order in the mean-field approximation

Expanding with respect to ξ(d) then leads to the following form:

$$\begin{array}{c}\hfill \phi ({\lambda }_1,{\lambda }_2)={\phi }_0({\lambda }_1)+{\phi }_0({\lambda }_2)+\xi (d){\phi }_1({\lambda }_1){\phi }_1({\lambda }_2),\end{array}$$(B.2)

that is, a factorization of the linear term in ξ(d). This can explicitly be shown in case of tree models (as described in the main text). This is also the case in perturbation theory.

In case of the tree models, we have

$$\begin{array}{c}\hfill {\phi }_1(\lambda )=\frac{\tau (\lambda )}{\overline{\xi }},\end{array}$$(B.3)

where $\overline{\xi }$ is the averaged correlation function within the cells. This is obtained assuming a mean-field approximation. We expect subleading corrections when d becomes comparable to the size of the cells.

B.2. Extending the previous case to the three variable case

In addition to the two variables ρ₁ and ρ₂, we introduce here the variable ρ_s, which is the density within the sample. We assume that the correlation functions are built with the same model. Here the small parameter is the correlation between two cells V₁ and V₂ (at positions x₁ and x₂) and the correlation function at sample size. It is natural in the context we are interested in to assume that these two quantities are on the same order.

We further assume we can use the mean-field approximation for the two cells V₁ and V₂. This is not a priori the case for the density in the whole sample, however. We therefore derive the results without this approximation. The general expression is then

$$\begin{array}{cc}\hfill \phi ({\lambda }_s,{\lambda }_1,{\lambda }_2)=& {\lambda }_s{\displaystyle \int \mathrm{d}}{\mathbf{x}}_s(\zeta (\tau ({\mathbf{x}}_s))-\frac{1}{2}\tau ({\mathbf{x}}_s){\zeta }^{\prime }(\tau ({\mathbf{x}}_s)))\hfill \\ \hfill & +{\lambda }_1[\zeta ({\tau }_1)-\frac{1}{2}{\tau }_1{\zeta }^{\prime }({\tau }_1)]+{\lambda }_2[\zeta ({\tau }_2)-\frac{1}{2}{\tau }_2{\zeta }^{\prime }({\tau }_2)],\hfill \end{array}$$(B.4)

with the consistency relations

$$\begin{array}{cc}\hfill \tau ({\mathbf{x}}_s)& ={\lambda }_s{\displaystyle \int \mathrm{d}}{\mathbf{x}}_s^{\prime }\xi ({\mathbf{x}}_s,{\mathbf{x}}_s^{\prime }){\zeta }^{\prime }(\tau ({\mathbf{x}}_s^{\prime }))\hfill \\ \hfill & +{\lambda }_1\xi ({\mathbf{x}}_s,{\mathbf{x}}_1){\zeta }^{\prime }({\tau }_1)+{\lambda }_2\xi ({\mathbf{x}}_s,{\mathbf{x}}_2){\zeta }^{\prime }({\tau }_2)\hfill \end{array}$$(B.5)

$$\begin{array}{cc}\hfill {\tau }_1& ={\lambda }_s{\displaystyle \int \mathrm{d}}{\mathbf{x}}_s\xi ({\mathbf{x}}_1,{\mathbf{x}}_s){\zeta }^{\prime }(\tau ({\mathbf{x}}_s))\hfill \\ \hfill & +{\lambda }_1\overline{\xi }{\zeta }^{\prime }({\tau }_1)+{\lambda }_2{\xi }_{12}{\zeta }^{\prime }({\tau }_2)\hfill \end{array}$$(B.6)

$$\begin{array}{cc}\hfill {\tau }_2& ={\lambda }_s{\displaystyle \int \mathrm{d}}{\mathbf{x}}_s\xi ({\mathbf{x}}_2,{\mathbf{x}}_s){\zeta }^{\prime }(\tau ({\mathbf{x}}_s))\hfill \\ \hfill & +{\lambda }_1{\xi }_{12}{\zeta }^{\prime }({\tau }_1)+{\lambda }_2\overline{\xi }{\zeta }^{\prime }({\tau }_2).\hfill \end{array}$$(B.7)

We therefore derive the expression of φ(λ_s, λ₁, λ₂) up to linear order jointly in $\xi ({\mathbf{x}}_s,{\mathbf{x}}_s^{\mathrm{\prime }})$, ξ(x_s, x₁), ξ(x_s, x₁) and ξ₁₂ ≡ ξ(x₁, x₂). At zeroth order, τ(x_s) vanishes and

$$\begin{array}{c}\hfill {\tau }_1^{(0)}={\lambda }_1\overline{\xi }{\zeta }^{\prime }({\tau }_1),\end{array}$$(B.8)

with a similar relation for ${\tau }_2^{(0)}$. At linear order, we have

$$\begin{array}{cc}\hfill {\tau }^{(1)}({\mathbf{x}}_s)& ={\lambda }_s{\displaystyle \int \mathrm{d}}{\mathbf{x}}_s^{\prime }\xi ({\mathbf{x}}_s,{\mathbf{x}}_s^{\prime })\hfill \\ \hfill & +{\lambda }_1\xi ({\mathbf{x}}_s,{\mathbf{x}}_1){\zeta }^{\prime }({\tau }_1^{(0)})+{\lambda }_2\xi ({\mathbf{x}}_s,{\mathbf{x}}_2){\zeta }^{\prime }({\tau }_2^{(0)})\hfill \end{array}$$(B.9)

$$\begin{array}{cc}\hfill \mathcal{F}({\tau }_1^{(0)}){\tau }_1^{(1)}& ={\lambda }_s{\displaystyle \int \mathrm{d}}{\mathbf{x}}_s\xi ({\mathbf{x}}_1,{\mathbf{x}}_s)+{\lambda }_2{\xi }_{12}{\zeta }^{\prime }({\tau }_2^{(0)})\hfill \end{array}$$(B.10)

$$\begin{array}{cc}\hfill \mathcal{F}({\tau }_2^{(0)}){\tau }_2^{(1)}& ={\lambda }_s{\displaystyle \int \mathrm{d}}{\mathbf{x}}_s\xi ({\mathbf{x}}_2,{\mathbf{x}}_s)+{\lambda }_1{\xi }_{12}{\zeta }^{\prime }({\tau }_1^{(0)}),\hfill \end{array}$$(B.11)

where

$$\begin{array}{c}\hfill \mathcal{F}(\tau )\equiv (1-\frac{\tau {\zeta }^{″}(\tau )}{{\zeta }^{\prime }(\tau )}).\end{array}$$(B.12)

The resulting cumulant-generating function reads

$$\begin{array}{cc}\hfill \phi ({\lambda }_s,{\lambda }_1,{\lambda }_2)=& {\lambda }_s+{\lambda }_1[\zeta ({\tau }_1^{(0)})-\frac{1}{2}{\tau }_1^{(0)}{\zeta }^{\prime }({\tau }_1^{(0)})]\hfill \\ \hfill & +{\lambda }_2[\zeta ({\tau }_2^{(0)})-\frac{1}{2}{\tau }_2^{(0)}{\zeta }^{\prime }({\tau }_2^{(0)})]+\frac{{\lambda }_s}{2}{\displaystyle \int \mathrm{d}}{\mathbf{x}}_s{\tau }^{(1)}({\mathbf{x}}_s)\hfill \\ \hfill & +\frac{{\lambda }_1}{2}{\zeta }^{\prime }({\tau }_1^{(0)})\mathcal{F}({\tau }_1^{(0)}){\tau }_1^{(1)}+\frac{{\lambda }_2}{2}{\zeta }^{\prime }({\tau }_2^{(0)})\mathcal{F}({\tau }_2^{(0)}){\tau }_2^{(1)},\hfill \end{array}$$(B.13)

where the first three terms are at zeroth order and the last three are at linear order. Using the previous expression, we obtain

$$\begin{array}{cc}\hfill \phi ({\lambda }_s,{\lambda }_1,{\lambda }_2)=& {\lambda }_s+{\phi }_0({\lambda }_1)+{\phi }_0({\lambda }_2)\hfill \\ \hfill & +\frac{{\lambda }_s^2}{2}{\displaystyle \int \mathrm{d}}{\mathbf{x}}_s\mathrm{d}{\mathbf{x}}_s^{\prime }\xi ({\mathbf{x}}_s,{\mathbf{x}}_s^{\prime })+{\lambda }_s{\displaystyle \int \mathrm{d}}{\mathbf{x}}_s\xi ({\mathbf{x}}_s,{\mathbf{x}}_1){\phi }_1({\lambda }_1)\hfill \\ \hfill & +{\lambda }_s{\displaystyle \int \mathrm{d}}{\mathbf{x}}_s\xi ({\mathbf{x}}_s,{\mathbf{x}}_2){\phi }_1({\lambda }_2)+{\phi }_1({\lambda }_1){\xi }_{12}{\phi }_1({\lambda }_2).\hfill \end{array}$$(B.14)

This relation is used to derive the expression of the sample bias functions in the next section.

B.3. Second order in the mean-field approximation

Results of Sect. B.1 can naturally be extended to any order in the cross-cell correlation function in the context of the tree-hierarchical models (as illustrated on Fig B.1). Up to second order, it takes the form

$$\begin{array}{cc}\hfill \phi ({\lambda }_1,{\lambda }_2)& ={\phi }_0({\lambda }_1)+{\phi }_0({\lambda }_2)+{\phi }_1({\lambda }_1){\xi }_{12}{\phi }_1({\lambda }_2)\hfill \\ \hfill & +\frac{1}{2}{\phi }_1^2({\lambda }_1){\xi }_{12}^2{\phi }_2({\lambda }_2)+\frac{1}{2}{\phi }_2({\lambda }_1){\xi }_{12}^2{\phi }_1^2({\lambda }_2),\hfill \end{array}$$(B.15)

where the function φ₂(λ) takes the form

$$\begin{array}{c}\hfill {\phi }_2(\lambda )=\frac{\lambda {\zeta }^{″}(\tau )}{1-\tau {\zeta }^{″}(\tau )/{\zeta }^{\prime }(\tau )}.\end{array}$$(B.16)

This last expression can be directly obtained through a perturbative expansion as presented in the previous subsection¹².

In case of the minimal tree model, this perturbative expansion can be directly compared with exact mean-field results¹³. This is shown in Fig. B.2. It shows that for a large regime in λ, the relation (B.15) provides a very accurate description of the joint cumulant-generating function down to a distance corresponding to overlapping cells. For overlapping cells, relation (B.15) continues to be accurate except for high values of λ. In general,

$$\begin{array}{c}\hfill \phi ({\lambda }_1,{\lambda }_2)\to {\phi }_0({\lambda }_1+{\lambda }_2)\end{array}$$(B.17)

when d → 0 (more accurate results in case of the minimal tree model are given in Appendix D.).

Using (B.15), we then derive corrective perturbative terms to the joint density PDF. More specifically, we have

$$\begin{array}{cc}\hfill P({\rho }_i,{\rho }_j)& =P({\rho }_i)P({\rho }_j)[1+b_1({\rho }_i){\xi }_{12}{b}_1({\rho }_j)\hfill \\ \hfill & +\frac{1}{2}{b}_2({\rho }_i){\xi }_{12}^2{b}_2({\rho }_j)+\frac{1}{2}{b}_2({\rho }_i){\xi }_{12}^2{q}_1({\rho }_j)\hfill \\ \hfill & +\frac{1}{2}{q}_1({\rho }_i){\xi }_{12}^2{b}_2({\rho }_j)],\hfill \end{array}$$(B.18)

where

$$\begin{array}{cc}\hfill b_2({\rho }_i)P({\rho }_i)& ={\displaystyle \int \frac{\mathrm{d}\lambda }{2\pi \mathrm{i}}}{\left[{\phi }_1(\lambda )\right]}^{2}exp[-\lambda {\rho }_i+\phi (\lambda )]\hfill \end{array}$$(B.19)

$$\begin{array}{cc}\hfill q_1({\rho }_i)P({\rho }_i)& ={\displaystyle \int \frac{\mathrm{d}\lambda }{2\pi \mathrm{i}}}{\phi }_2(\lambda )exp[-\lambda {\rho }_i+\phi (\lambda )].\hfill \end{array}$$(B.20)

For a sample with periodic boundary conditions, the average of ξ₁₂ vanishes, which a priori makes the other terms the leading contributors to the covariance elements. Equation (B.18) can be written as a sum of symmetric factorized terms,

$$\begin{array}{cc}\hfill P({\rho }_i,{\rho }_j)& =P({\rho }_i)P({\rho }_j)[1+b_1({\rho }_i){\xi }_{12}{b}_1({\rho }_j)\hfill \\ \hfill & +\frac{1}{2}(b_2({\rho }_i)+q_1({\rho }_i)){\xi }_{12}^2(b_2({\rho }_j)+q_1({\rho }_j))\hfill \\ \hfill & -\frac{1}{2}{q}_1({\rho }_i){\xi }_{12}^2{q}_1({\rho }_j)],\hfill \end{array}$$(B.21)

showing the eigenstructure of the resulting matrix and showing that it defines three different eigendirections at most.

B.4. Relative density joint PDFs and bias functions

We wish to compute the joint PDF of the density when expressed in terms of the survey average density ρ_s. In order to do so, we consider the joint density P(ρ_s, ρ_i, ρ_j), where ρ_s is the density in the sample and ρ_i and ρ_j are the densities in two cells at distance d.

We wish to compute the joint probability distribution function of ${\widehat{\rho }}_i$ and ${\widehat{\rho }}_j$, defined as

$$\begin{array}{c}\hfill {\widehat{\rho }}_i=\frac{{\rho }_i}{{\rho }_s},\end{array}$$(B.22)

and the joint distribution functions of ${\overline{\rho }}_i$ and ${\overline{\rho }}_j$, defined as

$$\begin{array}{c}\hfill {\overline{\rho }}_i={\rho }_i-{\rho }_s+1.\end{array}$$(B.23)

From these changes of variables, we have

$$\begin{array}{c}\hfill P({\widehat{\rho }}_i)={\displaystyle \int \mathrm{d}}{\rho }_s{\rho }_{s}P({\rho }_s,{\widehat{\rho }}_i{\rho }_s)\end{array}$$(B.24)

and

$$\begin{array}{c}\hfill P({\widehat{\rho }}_i,{\widehat{\rho }}_j)={\displaystyle \int \mathrm{d}}{\rho }_s{\rho }_s^{2}P({\rho }_s,{\widehat{\rho }}_i{\rho }_s,{\widehat{\rho }}_j{\rho }_s).\end{array}$$(B.25)

Similarly, we also have

$$\begin{array}{c}\hfill P({\overline{\rho }}_i)={\displaystyle \int \mathrm{d}}{\rho }_{s}P({\rho }_s,{\overline{\rho }}_i+{\rho }_s-1)\end{array}$$(B.26)

and

$$\begin{array}{c}\hfill P({\overline{\rho }}_i,{\overline{\rho }}_j)={\displaystyle \int \mathrm{d}}{\rho }_{s}P({\rho }_s,{\overline{\rho }}_i+{\rho }_s-1,{\overline{\rho }}_j+{\rho }_s-1).\end{array}$$(B.27)

We continue the calculations by expressing the joint PDF with the help of inverse Laplace transforms,

$$\begin{array}{cc}\hfill P({\rho }_s,{\rho }_i)& ={\displaystyle \int \frac{\mathrm{d}{\lambda }_s}{2\pi \mathrm{i}}}\frac{\mathrm{d}\lambda }{2\pi \mathrm{i}}\hfill \\ \hfill & \times exp[-{\lambda }_s{\rho }_s-\lambda {\rho }_i+\phi ({\lambda }_s,\lambda )]\hfill \end{array}$$(B.28)

$$\begin{array}{cc}\hfill P({\rho }_s,{\rho }_i,{\rho }_j)& ={\displaystyle \int \frac{\mathrm{d}{\lambda }_s}{2\pi \mathrm{i}}}\frac{\mathrm{d}{\lambda }_1}{2\pi \mathrm{i}}\frac{\mathrm{d}{\lambda }_2}{2\pi \mathrm{i}}\hfill \\ \hfill & \times exp[-{\lambda }_s{\rho }_s-{\lambda }_1{\rho }_i-{\lambda }_2{\rho }_j+\phi ({\lambda }_s,{\lambda }_1,{\lambda }_2)].\hfill \end{array}$$(B.29)

As a result,

$$\begin{array}{cc}\hfill P({\widehat{\rho }}_i)& ={\displaystyle \int \frac{\mathrm{d}{\lambda }_s}{2\pi \mathrm{i}}}\frac{\mathrm{d}\lambda }{2\pi \mathrm{i}}\frac{1}{{({\lambda }_s+\lambda {\widehat{\rho }}_i)}^2}exp[\phi ({\lambda }_s,\lambda )]\hfill \end{array}$$(B.30)

$$\begin{array}{cc}\hfill P({\widehat{\rho }}_i,{\widehat{\rho }}_j)& ={\displaystyle \int \frac{\mathrm{d}{\lambda }_s}{2\pi \mathrm{i}}}\frac{\mathrm{d}{\lambda }_1}{2\pi \mathrm{i}}\frac{\mathrm{d}{\lambda }_2}{2\pi \mathrm{i}}\frac{1}{{({\lambda }_s+{\lambda }_1{\widehat{\rho }}_i+{\lambda }_2{\widehat{\rho }}_j)}^3}\hfill \\ \hfill & \times exp[\phi ({\lambda }_s,{\lambda }_1,{\lambda }_2)]\hfill \end{array}$$(B.31)

after integration over ρ_s. The latter expressions can be expressed as

$$\begin{array}{cc}\hfill P({\widehat{\rho }}_i)& ={\displaystyle \int \frac{\mathrm{d}\lambda }{2\pi \mathrm{i}}}{\left[\frac{\partial \phi }{\partial {\lambda }_s}\right]}_{{|}_{{\lambda }_s=-\lambda {\widehat{\rho }}_i}}exp[\phi (-\lambda {\widehat{\rho }}_i,{\lambda }_i)]\hfill \end{array}$$(B.32)

$$\begin{array}{cc}\hfill P({\widehat{\rho }}_i,{\widehat{\rho }}_j)& ={\displaystyle \int \frac{\mathrm{d}{\lambda }_1}{2\pi \mathrm{i}}}\frac{\mathrm{d}{\lambda }_2}{2\pi \mathrm{i}}{[{\left(\frac{\partial \phi }{\partial {\lambda }_s}\right)}^2+\frac{{\partial }^2\phi }{\partial {\lambda }_s^2}]}_{{|}_{{\lambda }_s=-{\lambda }_1{\widehat{\rho }}_i-{\lambda }_2{\widehat{\rho }}_j}}\hfill \\ \hfill & \times exp[\phi (-{\lambda }_1{\widehat{\rho }}_i-{\lambda }_2{\widehat{\rho }}_j,{\lambda }_1,{\lambda }_2)].\hfill \end{array}$$(B.33)

In a similar manner, we can obtain the form of the joint PDF for $\{{\overline{\rho }}_i\}$,

$$\begin{array}{cc}\hfill P({\overline{\rho }}_i)=& {\displaystyle \int \frac{\mathrm{d}\lambda }{2\pi \mathrm{i}}}exp[-\lambda {\overline{\rho }}_i+\phi (-\lambda ,\lambda )]\hfill \end{array}$$(B.34)

$$\begin{array}{cc}\hfill P({\overline{\rho }}_i,{\overline{\rho }}_j)=& {\displaystyle \int \frac{\mathrm{d}{\lambda }_1}{2\pi \mathrm{i}}}\frac{\mathrm{d}{\lambda }_2}{2\pi \mathrm{i}}\hfill \\ \hfill & \times exp[-{\lambda }_1{\overline{\rho }}_i-{\lambda }_2{\overline{\rho }}_j+\phi (-{\lambda }_1-{\lambda }_2,{\lambda }_1,{\lambda }_2)].\hfill \end{array}$$(B.35)

We then use the relation (B.14) to compute the form of these functions.

Noting that the expressions $\int \mathrm{d}{\mathbf{x}}_0\mathrm{d}{\mathbf{x}}_0^{\mathrm{\prime }}\xi ({\mathbf{x}}_0,{\mathbf{x}}_0^{\mathrm{\prime }})$, ∫dx₀ ξ(x₀, x₁) take all the same averaged value when integrated over the sample, which we note ${\overline{\xi }}_s$, then at linear order in ${\overline{\xi }}_s$,

$$\begin{array}{cc}\hfill \frac{\partial \phi }{\partial {\lambda }_s}=& 1+{\overline{\xi }}_s[{\lambda }_s+{\phi }_1({\lambda }_1)+{\phi }_1({\lambda }_2)].\hfill \end{array}$$(B.36)

$$\begin{array}{cc}\hfill \frac{{\partial }^2\phi }{\partial {\lambda }_s^2}& ={\overline{\xi }}_s.\hfill \end{array}$$(B.37)

At the same order, we then have

$$\begin{array}{c}\hfill {\left(\frac{\partial \phi }{\partial {\lambda }_s}\right)}^2+\frac{{\partial }^2\phi }{\partial {\lambda }_s^2}=1+{\overline{\xi }}_s[1+2({\lambda }_s+{\phi }_1({\lambda }_1)+{\phi }_1({\lambda }_2))].\end{array}$$(B.38)

Combining both the expressions of $P({\widehat{\rho }}_i)$ and $P({\widehat{\rho }}_i,{\widehat{\rho }}_j)$ and expanding all terms at linear order in ${\overline{\xi }}_s$, we obtain

$$\begin{array}{cc}\hfill P({\widehat{\rho }}_i,{\widehat{\rho }}_j)-P({\widehat{\rho }}_i)P({\widehat{\rho }}_j)& ={\displaystyle \int \frac{\mathrm{d}{\lambda }_1}{2\pi \mathrm{i}}}\frac{\mathrm{d}{\lambda }_2}{2\pi \mathrm{i}}\hfill \\ \hfill & \times [1+{\overline{\xi }}_s(1+{\phi }_1({\lambda }_1)-{\lambda }_1{\widehat{\rho }}_i)(1+{\phi }_1({\lambda }_2)-{\lambda }_2{\widehat{\rho }}_j)]\hfill \\ \hfill & \times exp[-{\lambda }_1{\widehat{\rho }}_i-{\lambda }_2{\widehat{\rho }}_j+{\phi }_0({\lambda }_1)+{\phi }_0({\lambda }_2)].\hfill \end{array}$$(B.39)

This leads to the definition of the first sample bias function,

$$\begin{array}{c}\hfill b_{\mathrm{s}1}({\widehat{\rho }}_i)=\frac{1}{P({\widehat{\rho }}_i)}{\displaystyle \int \frac{\mathrm{d}\lambda }{2\pi \mathrm{i}}}(1+{\phi }_1(\lambda )-\lambda {\widehat{\rho }}_i)exp[-\lambda {\widehat{\rho }}_i+{\phi }_0(\lambda )],\end{array}$$(B.40)

which can be re-expressed in terms of the density-bias function defined in Eq. (18) and the derivative of $P({\widehat{\rho }}_i)$ with respect to ${\widehat{\rho }}_i$

$$\begin{array}{c}\hfill b_{\mathrm{s}1}({\widehat{\rho }}_i)=b({\widehat{\rho }}_i)+1+\frac{\mathrm{d}log(P({\widehat{\rho }}_i))}{\mathrm{d}log{\widehat{\rho }}_i}.\end{array}$$(B.41)

The second sample-bias function can be obtained in a similar manner. We indeed have

$$\begin{array}{cc}\hfill P({\overline{\rho }}_i,{\overline{\rho }}_j)-P({\overline{\rho }}_i)P({\overline{\rho }}_j)=& {\displaystyle \int \frac{\mathrm{d}{\lambda }_1}{2\pi \mathrm{i}}}\frac{\mathrm{d}{\lambda }_2}{2\pi \mathrm{i}}\hfill \\ \hfill & \times [1+{\overline{\xi }}_s({\phi }_1({\lambda }_1)-{\lambda }_1)({\phi }_1({\lambda }_2)-{\lambda }_2)]\hfill \\ \hfill & \times exp[-{\lambda }_1{\overline{\rho }}_i-{\lambda }_2{\overline{\rho }}_j+{\phi }_0({\lambda }_1)+{\phi }_0({\lambda }_2)],\hfill \end{array}$$(B.42)

which eventually leads to

$$\begin{array}{c}\hfill b_{\mathrm{s}2}({\overline{\rho }}_i)=b({\overline{\rho }}_i)+\frac{\mathrm{d}log(P({\overline{\rho }}_i))}{\mathrm{d}{\overline{\rho }}_i}.\end{array}$$(B.43)

B.5. Response to a change in amplitude in ξ

A close notion related to the density-bias function is how the PDF is changed when the parameters of the simulations are changed. In particular for tree models, the statistical properties are entirely determined by the amplitude of the two-point function, for instance, at cell size. This dependence can be made explicit by writing Eq. (17) as

$$\begin{array}{c}\hfill P({\rho }_i,\overline{\xi })=\frac{1}{\overline{\xi }}{\displaystyle \int \frac{\mathrm{d}\widehat{\lambda }}{2\pi \mathrm{i}}}exp(-\widehat{\lambda }\frac{{\rho }_i}{\overline{\xi }}+\frac{1}{\overline{\xi }}\psi (\widehat{\lambda }))\end{array}$$(B.44)

after the change of variable and function,

$$\begin{array}{c}\hfill \widehat{\lambda }=\overline{\xi }\lambda ,\psi (\widehat{\lambda })=\overline{\xi }\phi (\lambda ),\end{array}$$(B.45)

where then the expression $\psi (\widehat{\lambda })$ does not depend on $\overline{\xi }$ (only on the functional form of ζ). It follows that

$$\begin{array}{cc}\hfill -\overline{\xi }\frac{\partial P({\rho }_i,\overline{\xi })}{\partial \overline{\xi }}& =P({\rho }_i,\overline{\xi })+\rho \frac{\partial P({\rho }_i,\overline{\xi })}{\partial {\rho }_i}\hfill \\ \hfill & +{\displaystyle \int \frac{\mathrm{d}\lambda }{2\pi \mathrm{i}}}\phi (\lambda )exp(-\widehat{\lambda }\frac{{\rho }_i}{\overline{\xi }}+\phi (\lambda )).\hfill \end{array}$$(B.46)

This expression can be used to defined the function b_ξ(ρ_i) as

$$\begin{array}{c}\hfill b_{\xi }({\rho }_i)=-\frac{\partial logP({\rho }_i,\overline{\xi })}{\partial log\overline{\xi }}.\end{array}$$(B.47)

It appears that b_ξ(ρ_i) is very similar to b_s1(ρ_i), but the two are not equal in general.

B.6. Close cell results

B.6.1. CGF for 2 close cells

Saddle point approximation and close cell results. In case of two cells, the general system in the mean-field approximation leads to

$$\begin{array}{cc}\hfill {\tau }_1& ={\lambda }_1\overline{\xi }{\zeta }^{\prime }({\tau }_1)+{\lambda }_2{\xi }_{12}{\zeta }^{\prime }({\tau }_2)\hfill \end{array}$$(B.48)

$$\begin{array}{cc}\hfill {\tau }_2& ={\lambda }_2{\overline{\xi }}_{12}{\zeta }^{\prime }({\tau }_1)+{\lambda }_2\overline{\xi }{\zeta }^{\prime }({\tau }_2)\hfill \end{array}$$(B.49)

and

$$\begin{array}{cc}\hfill \phi ({\lambda }_1,{\lambda }_2)& ={\lambda }_1(\zeta ({\tau }_1)-\frac{1}{2}{\tau }_1{\zeta }^{\prime }({\tau }_1))\hfill \\ \hfill & +{\lambda }_2(\zeta ({\tau }_2)-\frac{1}{2}{\tau }_2{\zeta }^{\prime }({\tau }_2)).\hfill \end{array}$$(B.50)

We are interested here in the behavior of φ(λ₁, λ₂) when the two cells are close, that is, when ${\xi }_{12}\to \overline{\xi }$. When ${\xi }_{12}=\overline{\xi }$, τ₁ and τ₂ are also equal, making φ(λ₁, λ₂) a sole function of λ₁ + λ₂ and therefore forcing the joint PDF to be proportional to δ_Dirac(ρ_i − ρ_j). To be more precise, in this regime, ${\xi }_{12}\to \overline{\xi }$, δ_​ρ ≡ (ρ_i − ρ_j)/2 is expected to be distributed with a width of about ${\mathrm{\Delta }}_{\xi }^{1/2}={(\overline{\xi }-{\xi }_{12})}^{1/2}$. This suggests that in this limit, the difference λ₁ − λ₂ should be treated as a large quantity of about ${(\overline{\xi }-{\xi }_{12})}^{-1/2}$. The limit behavior of the joint CGF can then be explicitly computed in terms of

$$\begin{array}{cc}\hfill \lambda & ={\lambda }_1+{\lambda }_2,\hfill \end{array}$$(B.51)

$$\begin{array}{cc}\hfill \mu & ={\lambda }_1-{\lambda }_2,\hfill \end{array}$$(B.52)

$$\begin{array}{cc}\hfill {\mathrm{\Delta }}_{\xi }& =\overline{\xi }-{\xi }_{12}.\hfill \end{array}$$(B.53)

In this limit, we obtain

$$\begin{array}{cc}\hfill \tau & =\lambda \overline{\xi }{\zeta }^{\prime }(\tau )+\mu \overline{\xi }{\zeta }^{″}(\tau )\delta \tau ,\hfill \end{array}$$(B.54)

$$\begin{array}{cc}\hfill \delta \tau & =\frac{\mu }{2}{\mathrm{\Delta }}_{\xi }{\zeta }^{\prime }(\tau ),\hfill \end{array}$$(B.55)

leading to

$$\begin{array}{c}\hfill \tau =\lambda \overline{\xi }{\zeta }^{\prime }(\tau )+\frac{{\mu }^2\overline{\xi }}{4}{\mathrm{\Delta }}_{\xi }{\left({\zeta }^{\prime }{(\tau )}^2\right)}^{\prime }\end{array}$$(B.56)

and

$$\begin{array}{cc}\hfill \phi ({\lambda }_1,{\lambda }_2)& =\lambda (\zeta (\tau )-\frac{1}{2}\tau {\zeta }^{\prime }(\tau ))\hfill \\ \hfill & +\frac{{\mu }^2}{4}{\mathrm{\Delta }}_{\xi }({\zeta }^{\prime }{(\tau )}^2-\frac{1}{2}\tau {\left({\zeta }^{\prime }{(\tau )}^2\right)}^{\prime }).\hfill \end{array}$$(B.57)

The joint PDF of ρ_m = (ρ_i + ρ_j)/2 and δ_​ρ then reads

$$\begin{array}{c}\hfill P({\rho }_{\mathrm{m}},{\delta }_{\rho })={\displaystyle \int \frac{\mathrm{d}\lambda }{2\pi \mathrm{i}}}\frac{\mathrm{d}\mu }{2\pi \mathrm{i}}exp[-\lambda {\rho }_{\mathrm{m}}-\mu {\delta }_{\rho }+\phi (\lambda ,\mu )],\end{array}$$(B.58)

for which there is in general no closed form. We then need to rely on approximation schemes to complete the calculations.

B.6.2. Saddle point approximation

One of the approximations that can be used to evaluate Eq. (B.58) is to use the saddle point approximation. It has been used in the literature to compute the PDF (see Balian & Schaeffer 1989; Bernardeau 1992; Valageas 2002; Codis et al. 2016a). It is a priori valid when $\overline{\xi }$ is small (and not for too high values of the density). In this approach, the expression under the exponential is approximated by a quadratic form at its minimum. In practice, the latter is obtained from the resolution of the system

$$\begin{array}{cc}\hfill \frac{\partial \phi (\lambda ,\mu )}{\partial \lambda }& ={\rho }_{\mathrm{m}}\hfill \end{array}$$(B.59)

$$\begin{array}{cc}\hfill \frac{\partial \phi (\lambda ,\mu )}{\partial \mu }& ={\delta }_{\rho },\hfill \end{array}$$(B.60)

which leads to the implicit or explicit values of λ, μ, and τ at the saddle point position (we hereafter denote this with the subscript s),

$$\begin{array}{c}\hfill {\rho }_{\mathrm{m}}=\zeta ({\tau }_s),{\delta }_{\rho }={\mu }_s\frac{{\mathrm{\Delta }}_{\xi }}{2}{\left({\zeta }^{\prime }({\tau }_s)\right)}^2.\end{array}$$(B.61)

It is to be noted here that the value of τ_s is independent of δ_​ρ. At the saddle point position, we then have

$$\begin{array}{c}\hfill -{\lambda }_s{{\rho }_{\mathrm{m}}}_s-{\mu }_s{\delta }_{\rho }+\phi ({\lambda }_s,{\mu }_s)=-\frac{{{\delta }_{\rho }}^2}{{\mathrm{\Delta }}_{\xi }{\left[{\zeta }^{\prime }({\tau }_s)\right]}^2}-\frac{{\tau }_s^2}{2}.\end{array}$$(B.62)

This then suggests that the joint PDF is given by the product of the one-point PDF of ρ_m and a Gaussian distribution of δ_​ρ of width ${\mathrm{\Delta }}_{\xi }^{1/2}\zeta \prime ({\tau }_s)/\sqrt{2}$. For usual models, as described above, ζ′(τ_s)² scales like a power of ρ_m so that one suggested form for the joint PDF is the following:

$$\begin{array}{c}\hfill P({\rho }_{\mathrm{m}},{\delta }_{\rho })=P({\rho }_{\mathrm{m}})exp(-\frac{{{\delta }_{\rho }}^2}{{\mathrm{\Delta }}_{\xi }{\rho }_m^{\alpha }})\frac{1}{\sqrt{\pi {\mathrm{\Delta }}_{\xi }{{\rho }_{\mathrm{m}}}^{\alpha }}}.\end{array}$$(B.63)

Interestingly, the value of α can be related to the reduced skewness of the density field from the computation of ⟨δ_​ρ²ρ_m⟩_c, and in the context of tree hierarchical models, it leads to

$$\begin{array}{c}\hfill \alpha =\frac{2}{3}{S}_3.\end{array}$$(B.64)

The validity of this form clearly ought to be checked. Its simplicity nonetheless offers a good grasp of the contribution of close cells to the covariance matrix.

Appendix C: Minimal tree model

In the previous section, general formulae were given. The aim of this section is to account for more precise results obtained in the case of a specific hierarchical model, namely the minimal tree model, as described below. It can then serve as a toy model for the construction of the approximate form for the covariance matrix. We first recall that this model describes the Rayleigh Levy flights model.

C.1. One-point results in the mean-field approximation

The Rayleigh-Levy flight model makes it possible to build synthetic samples whose statistical properties follow the minimal model, that is, a hierarchical model with the following vertex-generating function:

$$\begin{array}{c}\hfill {\zeta }_{\mathrm{RL}}(\tau )=1+\tau +\frac{1}{4}{\tau }^2.\end{array}$$(C.1)

In the one-cell mean-field approximation, we have the equation

$$\begin{array}{c}\hfill \tau =\lambda \overline{\xi }{{\zeta }_{\mathrm{RL}}}^{\prime }(\tau ),\end{array}$$(C.2)

which can be solved in

$$\begin{array}{c}\hfill \tau =\frac{\lambda \xi }{1-\lambda \xi /2},\end{array}$$(C.3)

which leads to the following expression for the CGF:

$$\begin{array}{c}\hfill \phi (\lambda )=\frac{\lambda }{1-\lambda \overline{\xi }/2}.\end{array}$$(C.4)

The one-point PDF of the density can then be computed explicitly. It takes the form in the continuous limit of

$$\begin{array}{cc}\hfill P(\rho )& ={\displaystyle \int \frac{\mathrm{d}\lambda }{2\pi \mathrm{i}}}{e}^{-\lambda \rho +\phi (\lambda )}\hfill \\ \hfill & =e^{-\frac{2}{\overline{\xi }}}{\delta }_{\mathrm{Dirac}}(\rho )+\frac{4}{{\overline{\xi }}^2}{e}^{-\frac{2}{\overline{\xi }}(1+\rho )}{}_0{F}_1(2,\frac{4\rho }{\overline{\xi }}).\hfill \end{array}$$(C.5)

For this particular model, the void probability distribution (VPF) is nonzero even in the continuous limit. We recall here that the general expression of the VPF is given by exp(φ(−N̄)), which for the minimal model leads to

$$\begin{array}{c}\hfill P_0=exp(-2/\overline{\xi })\end{array}$$(C.6)

when N̄ → ∞.

The density-bias function can also be computed explicitly. For the minimal model, we have φ₁(λ)=φ(λ) so that

$$\begin{array}{cc}\hfill b_1(\rho )& =\frac{1}{P(\rho )}{\displaystyle \int \frac{\mathrm{d}\lambda }{2\pi \mathrm{i}}}{\phi }_1(\lambda )e^{-\lambda \rho +\phi (\lambda )}\hfill \\ \hfill & =\frac{{}_0{F}_1(1,\frac{4\rho }{\overline{\xi }})}{{}_0{F}_1(2,\frac{4\rho }{\overline{\xi }})}-\frac{2}{\overline{\xi }}\hfill \end{array}$$(C.7)

for ρ >  0. For this model, the fact that φ₁(λ)=φ(λ) implies that

$$\begin{array}{c}\hfill b_{\mathrm{s}1}(\rho )=b_{\xi }(\rho ).\end{array}$$(C.8)

This means that in the case of the minimal model, the density-bias function can be extracted from the functional form of the one-point PDF as

$$\begin{array}{c}\hfill b(\rho )P(\rho )=[-1-\frac{\mathrm{d}}{\mathrm{d}log\rho }-\frac{\mathrm{d}}{\mathrm{d}log\overline{\xi }}]P(\rho ,\overline{\xi }).\end{array}$$(C.9)

This is a somewhat remarkable identity (which can be extended to higher orders, as shown below.)

In this case, the second-order expansion leads to the form φ₂(λ) given by

$$\begin{array}{c}\hfill {\phi }_2(\lambda )=\frac{1}{2}{\phi }_1(\lambda ),\end{array}$$(C.10)

and we note that φ_p(λ) all vanish for p ≥ 3.

C.2. Two-cell results in the mean-field approximation

These mean-field calculations can be extended to the two-cell case. In this case, we have the system

$$\begin{array}{cc}\hfill {\tau }_1& ={\lambda }_1\overline{\xi }{{\zeta }_{\mathrm{RL}}}^{\prime }({\tau }_1)+{\lambda }_2{\xi }_{12}{{\zeta }_{\mathrm{RL}}}^{\prime }({\tau }_2)\hfill \end{array}$$(C.11)

$$\begin{array}{cc}\hfill {\tau }_2& ={\lambda }_1{\xi }_{12}{{\zeta }_{\mathrm{RL}}}^{\prime }({\tau }_1)+{\lambda }_2\overline{\xi }{{\zeta }_{\mathrm{RL}}}^{\prime }({\tau }_2)\hfill \end{array}$$(C.12)

when the two cells are of the same size. This leads to the following expression for the joint CGF:

$$\begin{array}{c}\hfill \phi ({\lambda }_1,{\lambda }_2)=\frac{{\lambda }_1+{\lambda }_2+({\xi }_{12}-\overline{\xi }){\lambda }_1{\lambda }_2}{1-({\lambda }_1+{\lambda }_2)\overline{\xi }/2-{\lambda }_1{\lambda }_2({\xi }_{12}^2-{\overline{\xi }}^2)/4}.\end{array}$$(C.13)

Remarkably, this expression can be written in terms of the one-cell CGF as

$$\begin{array}{c}\hfill \phi ({\lambda }_1,{\lambda }_2)=\frac{\phi ({\lambda }_1)+\phi ({\lambda }_2)+{\xi }_{12}\phi ({\lambda }_1)\phi ({\lambda }_2)}{1-{\xi }_{12}^2\phi ({\lambda }_1)\phi ({\lambda }_2)/4}.\end{array}$$(C.14)

This opens the possibility of computing the joint PDF to any order of ξ₁₂. The calculation of this expansion is made simple by the following observations: The corrective terms will make intervene functions of the forms

$$\begin{array}{c}\hfill b_n(\rho )P(\rho )\equiv {\displaystyle \int \frac{\mathrm{d}\lambda }{2\pi \mathrm{i}}}{[\phi (\lambda )]}^{n}exp(-\lambda \rho +\phi (\lambda )).\end{array}$$(C.15)

We further note that

$$\begin{array}{c}\hfill \frac{\mathrm{d}}{\mathrm{d}\lambda }\phi (\lambda )={\left(\frac{\overline{\xi }}{2}\right)}^2{[\psi (\lambda )]}^2\end{array}$$(C.16)

with

$$\begin{array}{c}\hfill \psi (\lambda )=\phi (\lambda )+\frac{2}{\overline{\xi }}.\end{array}$$(C.17)

We the define c_n(ρ) as

$$\begin{array}{c}\hfill c_n(\rho )P(\rho )\equiv {\displaystyle \int \frac{\mathrm{d}\lambda }{2\pi \mathrm{i}}}{[\psi (\lambda )]}^{n}exp(-\lambda \rho +\phi (\lambda )).\end{array}$$(C.18)

We then have on one side

$$\begin{array}{c}\hfill b_n(\rho )={\displaystyle \sum _{q=0}^n}{C}_n^q{(-\frac{2}{\overline{\xi }})}^q{c}_{n-q}(\rho ),\end{array}$$(C.19)

and on the other side

$$\begin{array}{c}\hfill \rho c_n(\rho )={\left(\frac{\overline{\xi }}{2}\right)}^2(c_{n+2}(\rho )+n{c}_{n+1}(\rho )),\end{array}$$(C.20)

which derives from the fact that

$$\begin{array}{c}\hfill \rho c_n(\rho )={\displaystyle \int \frac{\mathrm{d}\lambda }{2\pi \mathrm{i}}}exp(-\lambda \rho )\frac{\mathrm{d}}{\mathrm{d}\lambda }\{{[\psi (\lambda )]}^{n}exp(\phi (\lambda ))\}\end{array}$$(C.21)

after integration by parts. As a result, the expression of the join PDF to any order can be written as polynomials making intervene P(ρ₁), b(ρ₁) P(ρ₂) and b(ρ₁) only.

C.3. Perturbative expansion for close cells

Another interesting result is when the cell centers are close (so that cells overlap), as described above. In this case, the limit behavior of the joint CGF is given by

$$\begin{array}{c}\hfill \phi ({\lambda }_1,{\lambda }_2)=\frac{\lambda +{\mu }^2{\mathrm{\Delta }}_{\xi }/4}{1-\lambda \overline{\xi }/2-{\mu }^2{\mathrm{\Delta }}_{\xi }\overline{\xi }/8},\end{array}$$(C.22)

with

$$\begin{array}{cc}\hfill \lambda & ={\lambda }_1+{\lambda }_2,\hfill \end{array}$$(C.23)

$$\begin{array}{cc}\hfill \mu & ={\lambda }_1-{\lambda }_2,\hfill \end{array}$$(C.24)

$$\begin{array}{cc}\hfill {\mathrm{\Delta }}_{\xi }& =\overline{\xi }-{\xi }_{12}.\hfill \end{array}$$(C.25)

It is then remarkable to see the result can be expressed with the sole one-cell CGF,

$$\begin{array}{c}\hfill \phi ({\lambda }_1,{\lambda }_2)=\phi (\lambda +{\mu }^2{\mathrm{\Delta }}_{\xi }/4).\end{array}$$(C.26)

In other words, the GFC of the variables ρ_m = (ρ₁ + ρ₂)/2 and δ_​ρ = (ρ₁ − ρ₂)/2 is given by Eq. (C.26). It is possible to compute the joint PDF,

$$\begin{array}{cc}\hfill P({\rho }_{\mathrm{m}},{\delta }_{\rho })& ={\displaystyle \int \frac{\mathrm{d}\lambda }{2\pi \mathrm{i}}}{\displaystyle \int \frac{\mathrm{d}\mu }{2\pi \mathrm{i}}}\hfill \\ \hfill & \times exp(-\lambda {\rho }_{\mathrm{m}}-\mu {\delta }_{\rho }+\phi (\lambda +{\mu }^2{\mathrm{\Delta }}_{\xi }/4)),\hfill \end{array}$$(C.27)

with the change of variable

$$\begin{array}{c}\hfill \stackrel{\sim }{\lambda }=\lambda +{\mu }^2\frac{{\mathrm{\Delta }}_{\xi }}{4}.\end{array}$$(C.28)

The integral in $\stackrel{\sim }{\lambda }$ leads to the one-cell PDF of the density ρ_m, whereas the integral in μ can be done explicitly as it is a quadratic form in μ, leading to a Gaussian distribution in δ_​ρ. The final PDF is given by

$$\begin{array}{c}\hfill P({\rho }_{\mathrm{m}},\delta \rho )=P({\rho }_{\mathrm{m}})\frac{1}{{\left[\pi {\mathrm{\Delta }}_{\xi }{\rho }_{\mathrm{m}}\right]}^{1/2}}exp(-\frac{{{\delta }_{\rho }}^2}{{\mathrm{\Delta }}_{\xi }{\rho }_{\mathrm{m}}}).\end{array}$$(C.29)

This shows that the joint PDF peaks for ρ₁ ∼ ρ₂ with a width that depends on the distance between the cells through the difference $\overline{\xi }-{\xi }_{12}$. Moreover, this form has no overlapping regime with the previous expansions of the joint PDF. It captures different aspects of the covariance calculations as listed below.

• The previous expression says that close cells contribute more specifically to the covariance when ρ₁ and ρ₂ are close. This suggests that Eq. (C.29) contributes mostly to the near diagonal terms, whereas off diagonal terms could still be well described by perturbative expansions, as described before.

• As noted before, perturbative expansions are closely related to supersample effects. They encode the way in which the local densities are jointly correlated with long-wavelength modes. This is not the case in Eq. (C.29). It rather captures how a rare event, such as a peak, can contribute to the covariance elements: if there is a peak somewhere, nearby cells are likely to have a similar density up to distances for which $\overline{\xi }-{\xi }_{12}$ remains small enough.

The above development can be pursued to any order in ${\mathrm{\Delta }}_{\xi }^{1/2}$ provided the following recipe is applied:

$$\begin{array}{c}\hfill {\delta }_{\rho }\sim {\mathrm{\Delta }}_{\xi }^{1/2}\mathrm{and}\mu \sim \frac{1}{{\mathrm{\Delta }}_{\xi }^{1/2}}.\end{array}$$(C.30)

Then the joint density can be computed to any order in ${\mathrm{\Delta }}_{\xi }^{1/2}$, making use of the very same expressions b_n(ρ).

The next-to-leading order in ${\mathrm{\Delta }}_{\xi }^{1/2}$ is thus given by

$$\begin{array}{cc}& P^{(2)}({\rho }_{\mathrm{m}},\delta \rho )=P^{(0)}({\rho }_{\mathrm{m}},\delta \rho )\hfill \\ \hfill & \times \{\frac{{\delta }_{\rho }^2}{2\overline{\xi }{\rho }_m}-\frac{{\rho }_m{\mathrm{\Delta }}_{\xi }}{{\overline{\xi }}^2}-\frac{{\mathrm{\Delta }}_{\xi }}{4\overline{\xi }}-\frac{{\mathrm{\Delta }}_{\xi }}{{\overline{\xi }}^2}-\frac{{\delta }_{\rho }^4}{4{\rho }_m^3{\mathrm{\Delta }}_{\xi }}+\frac{3{\delta }_{\rho }^2}{4{\rho }_m^2}-\frac{3{\mathrm{\Delta }}_{\xi }}{16{\rho }_m}\hfill \\ \hfill & +\frac{(\overline{\xi }{\rho }_m{\mathrm{\Delta }}_{\xi }-2\overline{\xi }{\delta }_{\rho }^2+8{\rho }_m^2{\mathrm{\Delta }}_{\xi })}{8\overline{\xi }{\rho }_m^2}\frac{{}_0{F}_1(1,\frac{4{\rho }_m}{{\overline{\xi }}^2})}{{}_0{F}_1(2,\frac{4{\rho }_m}{{\overline{\xi }}^2})}\},\hfill \end{array}$$(C.31)

and the expansion can be extended in any (even) order in ${\mathrm{\Delta }}_{\xi }^{1/2}$. Fig. C.1 illustrates the convergence properties of these expansions. Depending on ${\xi }_{12}/\overline{\xi }$, either the expansion in ξ₁₂ or that in ${\mathrm{\Delta }}_{\xi }^{1/2}$ gives a very accurate estimate of the joint PDF. It opens the way to computing the covariance matrix starting in the two-cell mean-field approximation (C.13).

C.4. Construction of the theoretical covariance matrix for the minimal tree model

The previous form can be used to compute the covariance matrix for the minimal tree model in simple implementations. It relies on analytic forms for both the two-point cell correlation functions, which can formally be written as

$$\begin{array}{c}\hfill {\xi }_{12}(r)={\displaystyle \int {\mathrm{d}}^2}\mathbf{k}{W}_{2D}^2(kR)J_0(kr)P(k)\end{array}$$(C.32)

for a given power spectrum. We also make use of the form P_s(r_d) given in footnote 2 to derive the PDF of cell distances. We then have all the required ingredients to compute the elements of the covariance matrix in the mean-field approximation,

$$\begin{array}{c}\hfill {\mathrm{Cov}}_{\mathrm{MF}}({\rho }_1,{\rho }_2)={\displaystyle \int \mathrm{d}}r{P}_d(r)P_{\mathrm{MF}}^{(8)}({\rho }_1,{\rho }_2;\overline{\xi },{\xi }_{12}(r{L}_{\mathrm{sample}})).\end{array}$$(C.33)

In practice, $P_{\mathrm{MF}}^{(8)}({\rho }_1,{\rho }_2;\overline{\xi },{\xi }_{12})$ is computed from the eighth-order expansion either in ξ₁₂ when ${\xi }_{12}/\overline{\xi }<0.4$ or in Δ_ξ when ${\xi }_{12}/\overline{\xi }>0.4$. This is used to explore the detailed properties of the covariance matrix and the validity of approximate schemes.

C.5. Joined PDF for relative densities

The minimal model allows us also to pursue the computation of the joint PDF for the variables $\{{\widehat{\rho }}_i\}$ or $\{{\overline{\rho }}_i\}$ in all regimes. The first step is to extend Eq. (B.14) to a regime in which ξ₁₂ is not assumed to be small. We find that

$$\begin{array}{cc}\hfill \phi ({\lambda }_s,{\lambda }_1,{\lambda }_2)& ={\lambda }_s+\phi ({\lambda }_1,{\lambda }_2)+\frac{{\lambda }_s^2}{2}{\displaystyle \int \mathrm{d}}{\mathbf{x}}_s\mathrm{d}{\mathbf{x}}_s^{\prime }\xi ({\mathbf{x}}_s,{\mathbf{x}}_s^{\prime })\hfill \\ \hfill & +{\lambda }_s{\displaystyle \int \mathrm{d}}{\mathbf{x}}_s\xi ({\mathbf{x}}_s,{\mathbf{x}}_1){\phi }_c({\lambda }_1,{\lambda }_2)\hfill \\ \hfill & +{\lambda }_s{\displaystyle \int \mathrm{d}}{\mathbf{x}}_s\xi ({\mathbf{x}}_s,{\mathbf{x}}_2){\phi }_c({\lambda }_2,{\lambda }_1),\hfill \end{array}$$(C.34)

where φ_c(λ₁, λ₂) is given by

$$\begin{array}{c}\hfill {\phi }_c({\lambda }_1,{\lambda }_2)=\frac{{\lambda }_1+({\xi }_{12}-\overline{\xi }){\lambda }_1{\lambda }_2/2}{1-({\lambda }_1+{\lambda }_2)\overline{\xi }/2-{\lambda }_1{\lambda }_2({\xi }_{12}^2-{\overline{\xi }}^2)/4},\end{array}$$(C.35)

and we can note that

$$\begin{array}{c}\hfill {\phi }_c({\lambda }_1,{\lambda }_2)+{\phi }_c({\lambda }_2,{\lambda }_1)=\phi ({\lambda }_1,{\lambda }_2).\end{array}$$(C.36)

At leading order in ξ_s, that is, when we assume that the density fluctuations at sample size are much smaller than at smoothing scale, this expression then reduces to

$$\begin{array}{c}\hfill \phi ({\lambda }_s,{\lambda }_1,{\lambda }_2)={\lambda }_s+\frac{{\lambda }_s^2}{2}{\xi }_s+\phi ({\lambda }_1,{\lambda }_2)(1+{\lambda }_s{\xi }_s).\end{array}$$(C.37)

We can then exploit this relation to compute the $P({\widehat{\rho }}_i,{\widehat{\rho }}_j)$ and $P({\overline{\rho }}_i,{\overline{\rho }}_j)$ from Eqs. (25) and (31), respectively. We then have at leading order in ξ_s

$$\begin{array}{cc}\hfill P_{s1}({\widehat{\rho }}_i,{\widehat{\rho }}_j)& =P({\widehat{\rho }}_i,{\widehat{\rho }}_j)+{\xi }_s{\displaystyle \int \frac{\mathrm{d}{\lambda }_1}{2\pi \mathrm{i}}}\frac{\mathrm{d}{\lambda }_2}{2\pi \mathrm{i}}\hfill \\ \hfill & \times (1+2{\lambda }_s+\frac{1}{2}{\lambda }_s^2+(2+{\lambda }_s)\phi ({\lambda }_1,{\lambda }_2))\hfill \\ \hfill & \times exp{[{\lambda }_s+\phi ({\lambda }_1,{\lambda }_2)]}_{{|}_{{\lambda }_s=-{\widehat{\rho }}_i{\lambda }_1-{\widehat{\rho }}_j{\lambda }_2}}\hfill \end{array}$$(C.38)

and

$$\begin{array}{cc}\hfill P_{s2}({\overline{\rho }}_i,{\overline{\rho }}_j)& =P({\overline{\rho }}_i,{\overline{\rho }}_j)\hfill \\ \hfill & +{\xi }_s{\displaystyle \int \frac{\mathrm{d}{\lambda }_1}{2\pi \mathrm{i}}}\frac{\mathrm{d}{\lambda }_2}{2\pi \mathrm{i}}{(\frac{1}{2}{\lambda }_s^2+{\lambda }_s\phi ({\lambda }_1,{\lambda }_2))}_{{|}_{{\lambda }_s=-{\lambda }_1-{\lambda }_2}}\hfill \\ \hfill & \times exp[-{\lambda }_1{\overline{\rho }}_i-{\lambda }_2{\overline{\rho }}_j+\phi ({\lambda }_1,{\lambda }_2)].\hfill \end{array}$$(C.39)

To complete the formal calculation of these expressions, we introduce the function

$$\begin{array}{cc}\hfill P_b({\rho }_i,{\rho }_j,\overline{\xi },{\xi }_{12})=& {\displaystyle \int \frac{\mathrm{d}{\lambda }_1}{2\pi \mathrm{i}}}\frac{\mathrm{d}{\lambda }_2}{2\pi \mathrm{i}}\phi ({\lambda }_1,{\lambda }_2)\hfill \\ \hfill & \times exp[-{\lambda }_1{\rho }_i-{\lambda }_2{\rho }_j+\phi ({\lambda }_1,{\lambda }_2)].\hfill \end{array}$$(C.40)

We can first note that Eq. (C.9) can be extended to

$$\begin{array}{cc}\hfill P_b({\rho }_i,{\rho }_j,\overline{\xi },{\xi }_{12})=& (-2-\frac{\mathrm{d}}{\mathrm{d}log{\rho }_i}\hfill \\ \hfill & -\frac{\mathrm{d}}{\mathrm{d}log{\rho }_j}-\frac{\mathrm{d}}{\mathrm{d}log\overline{\xi }}-\frac{\mathrm{d}}{\mathrm{d}log{\xi }_{12}})P({\rho }_i,{\rho }_j,\overline{\xi },{\xi }_{12}).\hfill \end{array}$$(C.41)

This comes from the observation that

$$\begin{array}{cc}\hfill P({\rho }_i,{\rho }_j,\overline{\xi },{\xi }_{12},\eta )& \equiv {\displaystyle \int \frac{\mathrm{d}{\lambda }_1}{2\pi \mathrm{i}}}\frac{\mathrm{d}{\lambda }_2}{2\pi \mathrm{i}}\hfill \\ \hfill & \times exp[-{\lambda }_1{\rho }_i-{\lambda }_2{\rho }_j+\eta \phi ({\lambda }_1,{\lambda }_2,\overline{\xi },{\xi }_{12})]\hfill \end{array}$$

can also be written

$$\begin{array}{cc}\hfill P({\rho }_i,{\rho }_j,\overline{\xi },{\xi }_{12},\eta )& =\frac{1}{{\eta }^2}{\displaystyle \int \frac{\mathrm{d}{\widehat{\lambda }}_1}{2\pi \mathrm{i}}}\frac{\mathrm{d}{\widehat{\lambda }}_2}{2\pi \mathrm{i}}\hfill \\ \hfill & \times exp[-{\widehat{\lambda }}_1\frac{{\rho }_i}{\eta }-{\widehat{\lambda }}_2\frac{{\rho }_j}{\eta }+\phi ({\widehat{\lambda }}_1,{\widehat{\lambda }}_2,\frac{\overline{\xi }}{\eta },\frac{{\xi }_{12}}{\eta })]\hfill \\ \hfill & =\frac{1}{{\eta }^2}P(\frac{{\rho }_i}{\eta },\frac{{\rho }_j}{\eta },\frac{\overline{\xi }}{\eta },\frac{{\xi }_{12}}{\eta },1)\hfill \end{array}$$(C.42)

and that

$$\begin{array}{c}\hfill P_b({\rho }_i,{\rho }_j,\overline{\xi },{\xi }_{12})={\frac{\partial }{\partial \eta }}_{{|}_{\eta =1}}P({\rho }_i,{\rho }_j,\overline{\xi },{\xi }_{12},\eta ).\end{array}$$(C.43)

The final expression of the PDF of the relative densities can then be obtained by noting that applying a multiplicative factor λ_i to the moment-generating function is equivalent to the application of the operator ∂/∂ρ_i to the final expression this finally leads to the following forms:

$$\begin{array}{cc}& P_{s1}({\widehat{\rho }}_i)=P({\widehat{\rho }}_i)+{\xi }_s[({\widehat{\rho }}_i\frac{\partial }{\partial {\widehat{\rho }}_i}+\frac{{\widehat{\rho }}_i^2}{2}\frac{{\partial }^2}{\partial {\widehat{\rho }}_i^2})P({\widehat{\rho }}_i)\hfill \\ \hfill & +(1+{\widehat{\rho }}_i\frac{\partial }{\partial {\widehat{\rho }}_i})(b({\widehat{\rho }}_i)P({\widehat{\rho }}_i))]\hfill \end{array}$$(C.44)

$$\begin{array}{cc}& P_{s1}({\widehat{\rho }}_i,{\widehat{\rho }}_j)=P({\widehat{\rho }}_i,{\widehat{\rho }}_j)\hfill \\ \hfill & +{\xi }_s[(1+2{\displaystyle \sum _i}{\widehat{\rho }}_i\frac{\partial }{\partial {\widehat{\rho }}_i}+\frac{1}{2}{\displaystyle \sum _i}{\widehat{\rho }}_i^2\frac{{\partial }^2}{\partial {\widehat{\rho }}_i^2})P({\widehat{\rho }}_i,{\widehat{\rho }}_j)\hfill \\ \hfill & +(2+{\displaystyle \sum _i}{\widehat{\rho }}_i\frac{\partial }{\partial {\widehat{\rho }}_i})P_b({\widehat{\rho }}_i,{\widehat{\rho }}_j)]\hfill \end{array}$$(C.45)

for the ${\widehat{\rho }}_i={\rho }_i/{\rho }_s$ and

$$\begin{array}{cc}\hfill P_{s2}({\overline{\rho }}_i)& =P({\overline{\rho }}_i)+{\xi }_s[\frac{1}{2}\frac{{\partial }^2}{\partial {\overline{\rho }}_i^2}P({\overline{\rho }}_i)\hfill \\ \hfill & +\frac{\partial }{\partial {\overline{\rho }}_i}(b({\overline{\rho }}_i)P({\overline{\rho }}_i))]\hfill \end{array}$$(C.46)

$$\begin{array}{cc}{P}_{s2}({\overline{\rho }}_i,{\overline{\rho }}_j)& =P({\overline{\rho }}_i,{\overline{\rho }}_j)+{\xi }_s[\frac{1}{2}\sum _i\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{\overline{\rho }}_i^2}P({\overline{\rho }}_i,{\overline{\rho }}_j)\\ & +\sum _i\frac{\mathrm{\partial }}{\mathrm{\partial }{\overline{\rho }}_i}{P}_b({\hat{\rho }}_i,{\hat{\rho }}_j)]\end{array}$$(C.47)

for ${\overline{\rho }}_i={\rho }_i-{\rho }_s+1$. These relations can then be applied to the expressions of the joint density such as $P_{\mathrm{MF}}^{(8)}({\rho }_i,{\rho }_j;\overline{\xi },{\xi }_{12})$ found in the previous subsection.

All Figures