© The Authors 2024

1. Introduction

The geometry and structure of cosmic fields form a complex physical system that develops from homogeneity through the interplay of expansion and long-range forces. Its statistical properties should emulate those of such a class of systems. From the perspective of observational cosmology, its evolution mirrors both the history of the universe’s expansion rate and the dynamic growth of embedded substructures (see e.g. Bardeen et al. 1986; Bernardeau et al. 2002, for an account of the emergence of structure due to gravitational instability). While the small-scale (galaxy-sized) structures are determined by the interaction of gravitational effects and complex baryonic physics, the large-scale structure is shaped solely by the development of the gravitational instabilities.

The most standard approach to describe the outcome of this evolution is to consider the matter correlation functions or (equivalently in Fourier space) the spectra (e.g. Peebles & Groth 1975; Fry 1985; Bernardeau 1994; Scoccimarro et al. 1998; Cappi et al. 2015). Those quantities capture however only partially the outcome of the gravitational processes. In particular gravitational instabilities leads to the formation of large scale structures, such as self-gravitating massive halos embedded in an intricate cosmic web made of voids, walls and filaments that are woven together under the influence of gravity (Bond et al. 1996). Inspired by such emerging features, a related alternative to quantify the properties of the field is to explore the topology of the excursion of the cosmic web. It can be done both in real space (e.g. Matsubara 1994; Gay et al. 2012) and redshift space (e.g. Matsubara 1996; Codis et al. 2013), using tools such as the void probability function (VPD, e.g. White 1979; Sheth & van de Weygaert 2004), the Euler-Poincaré characteristic (e.g. Gott et al. 1986; Park & Gott 1991; Appleby et al. 2018), or (more generally) Minkowski functionals (e.g. Mecke et al. 1994; Schmalzing & Buchert 1997), persistent homology (e.g. Sousbie et al. 2011; Pranav et al. 2017), and Betti numbers (e.g. Park et al. 2013; Feldbrugge et al. 2019). Given the duality established by Morse-Smale theory (Forman 2002) between the geometry or topology of excursion, and the loci of null gradients of the underlying density field, a summary statistics is given by the point process of these critical points (e.g. Bardeen et al. 1986; Bond et al. 1996; Gay et al. 2012; Cadiou et al. 2020). For instance, there has recently been some interest in also using the clustering of such points as cosmological probes (Baldauf et al. 2021; Shim et al. 2021), for instance, extracted from Lyman-α tomography (Kraljic et al. 2022). Unfortunately, the theory for capturing their statistics has been limited to the quasi Gaussian limit (e.g. Pogosyan et al. 2009), or slightly beyond (Bernardeau et al. 2015), while relying on the large deviation principle. This puts limitation on its realm of application to larger scales only (or involves relying on calibration over N-body simulations).

A notable cosmologically relevant counterexample is provided by Rayleigh-Levy flights¹ (Mandelbrot 1975; Peebles 1980; Szapudi & Colombi 1996; Zimbardo & Perri 2013; Uchaikin 2019), which are simply defined as a Markov chain point process whose jump probability depends on some power law of the length of the jump alone. In 3D, it acts as a coarse proxy to describing the motion of dark halos (Sefusatti & Scoccimarro 2005; Trotta & Zimbardo 2015), but its definition can be extended to arbitrary dimensions. Rayleigh-Levy flights were first introduced in astrophysics by Holtsmark (1919) in the context of fluctuations in gravitational systems (Litovchenko 2021). They include some degree of connection to first-passage theory (Metzler 2019), underlying Press Schechter theory for halo formation (Press & Schechter 1974), and subsequent mass accretion (Musso et al. 2018). Rayleigh-Levy flights have also attracted lots of attention beyond cosmology, (from anomalous cosmic rays diffusion Wilk & Włodarczyk 1999; Boldyrev & Gwinn 2003, to ISM scintillation) or indeed beyond astronomy (from turbulence, Shlesinger et al. 1987; Sotolongo-Costa et al. 2000, to earthquakes), and even biology (e.g. Reynolds 2018) or risk management (e.g. Bouchaud & Potters 2003), as they captures anomalous diffusion (see e.g. Bouchaud & Georges 1990; Dubkov et al. 2008, and references therein).

Levy flights also provide an appealing playground with genuine scale-independent developed non-Gaussianities (that cannot be mimicked by a mere local nonlinear transformation of the fields), whose amplitudes resemble what is generally expected in gravitational density fields. These remarks would be of limited interest however if we had no way of exploring the statistical properties of such fields. Building on the early results of Peebles (1980) and of Bernardeau & Schaeffer (1999), Bernardeau (2022) derived a large corpus of properties of such fields. They were concerned to a large extend on the joint density probability distribution function (PDFs) at large separations.

In this paper, we focus on the local behaviour of the field, trying to first grasp the expected behaviour of the local density and its derivatives, so as to identify critical points in such a density field. The starting point of those investigations is based on the derivation of the cumulant generating function (CGF) in multiple cells. This allows us to predict their one and two point statistics to arbitrary order in the variance of the field. For clarity, the main text focuses on results, while all the derivations are given in the appendices.

Section 2 recalls the main relevant properties of Rayleigh-Levy flights. Section 3 presents the number count of extrema of flights in 1D, 2D, and 3D. Section 4 compute their clustering. Comparisons to Monte Carlo simulations are provided throughout for validation for the 1D flights. The companion paper will provide them for 2 and 3D flights. Section 5 presents our conclusions.

2. Rayleigh-Levy flight model

The random model, as introduced first in a cosmological context by J. Peebles (Peebles 1980), is defined as a Markov random walk where the PDF of the step length, ℓ,  follows a cumulative distribution function given by:

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

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

where ℓ₀ is a small-scale regularisation parameter. The random walks are dominated by rare, large events rather than the accumulation of many steps.

Calling f₁(r) density of the subsequent point (the first descendant of a given point at position r₀) at position r, we have:

$$\begin{array}{c}\hfill {\displaystyle {\int }_{\ell }^{\infty }}{f}_1(\mathit{r}){\mathcal{U}}_{\mathrm{sph}}(D)D{r}^{D-1}\mathrm{d}r=P(>\ell ),\end{array}$$

where 𝒰_sph(D) is the volume of a D-dimensional sphere of unit radius, 𝒰_sph(D) = π^D/2/Γ(D/2 − 1). For |r − r₀|> ℓ₀, it leads to

$$\begin{array}{c}\hfill f_1(\mathit{r})=\frac{\alpha \mathrm{\Gamma }(D/2+1)}{D{\pi }^{D/2}}\frac{{\ell }_0^{\alpha }}{|\mathit{r}-{\mathit{r}}_0{|}^{D+\alpha }},\mathrm{in}\mathrm{D}\text{-}\mathrm{dimensional}\mathrm{space}.\end{array}$$(2)

We then consider its Fourier transform, ψ(k) as:

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

The two-point correlation functions for Rayleigh-Levy flights between positions r₁ and r₂ is then given by

$$\begin{array}{c}\hfill \xi ({\mathit{r}}_1,{\mathit{r}}_2)=\frac{1}{n}{\displaystyle \int \frac{{\mathrm{d}}^{\mathrm{D}}\mathit{k}}{{(2\pi )}^{\mathrm{D}}}}{e}^{i\mathit{k}·({\mathit{r}}_2-{\mathit{r}}_1)}\frac{2}{1-\psi (k)},\end{array}$$(4)

where n is the number density of points in the sample. When kl₀ ≪ 1, then 1 − ψ(k)∼(ℓ₀k)^α, the large distance correlation function behaves as:

$$\begin{array}{c}\hfill \xi ({\mathit{r}}_1,{\mathit{r}}_2)\sim r^{\alpha -\mathrm{D}}{\ell }_0^{-\alpha }/n.\end{array}$$(5)

The expected behaviour of the large-distance correlation function is shown on Fig. 1 below, which clearly shows that this form is valid as soon as r ≥ ℓ₀.

We note that the numerical system is driven by two parameters, ℓ₀ and n, with ${\xi }_0\sim 1/(n{\mathcal{l}}_0^{\alpha })$. In practice one wants to make sure that the filtering scale is significantly larger than both ℓ₀ and ℓ_n ≡ 1/n^1/D and that the survey size, ℓ_𝒮, is much larger than the filtering scale. When measuring the correlation properties of critical points, we also want to make sure that ℓ_𝒮 is much larger than the scales at which the correlation functions are measured. It is always possible to meet these requirements, and reach any amplitude of ξ, provided there is access to enough resources to build Rayleigh-Levy flights with large number of points.

In practice, in the 1D numerical simulations we exploited, we have the following:

α	=	0.5
n	=	3/pixel
ℓ0	=	0.02 pixel
rG	=	2 pixel
ℓ𝒮	=	105 pixel size
Nrealisations	=	200

providing very accurate determination of all one-point quantities (density and critical points, etc.) and good estimates of their correlation, as shown below.

A key property of Rayleigh-Levy flights involves the structure of the higher order correlation functions. In such a model, not only are they known but they also display the type of behaviour that is generally expected in the context of cosmological models. We know from theory in the perturbative regime (Bernardeau et al. 2002), and to a large extend from observations (Yu & Hou 2022), that the p-order matter correlation functions scale like the power p − 1 of the two-point function. This assumption and its consequences are addressed in more details in Appendices A–C.

Conversely, from a theoretical perspective, the Rayleigh-Levy flights model falls within the class of models whose p-point correlation functions follow the so-called hierarchical Ansatz (Fry 1984), so that they can be expressed as a sum of products of two-point correlations:

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

This process actually embodies a more specific class of models, which correspond to cases where the Q_p(t) values depend only on the vertex composition of the tree, that is,

$$\begin{array}{c}\hfill Q_p(t)={\displaystyle \prod _{\mathrm{vertices}\in t}}{\nu }_q,\end{array}$$(7)

where the vertices values, ν_q, are numbers associated with each vertex and depend only on their connectivity, q being the number of lines a given vertex is connected to (see Appendix B and Bernardeau & Schaeffer 1999). In the case of the flight models we consider here, we have only two non-zero vertices, ν₁ and ν₂, so that the trees we have to consider are, in fact, simply connected lines.

In general, the corresponding cumulant generating function of the underlying density field in cells of arbitrary profile given by 𝒲_i(x) can be written as:

$$\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!}.\end{array}$$(8)

We note that in the Gaussian limit, this expression is restricted to terms with p₁ + … + p_n = 2. We also note that, in general, the full complexity of this expression is hard to exploit and can only be done in a perturbative sense, exploring, for instance, the consequences of terms at cubic order. This complexity is at the heart of the Edgeworth expansion (Kendall & Stuart 1977; Sellentin et al. 2017).

A remarkable property of the tree models in general, and of the Levy flight model in particular, is that the discrete summations that appear in this expression (over p₁, …, p_n), can be done explicitly. This is clearly an appealing feature of these models, as it allows for an exploration of their properties in a regime that can be arbitrarily far from the Gaussian case. This is the motivation for this study.

The details of the derivation of the resulting form for the CGF, φ(λ₁, …, λ_n), is given in Appendix B, which shows that:

$$\begin{array}{c}\hfill \phi ({\lambda }_1,\cdots ,{\lambda }_n)={\displaystyle \sum _j}{\lambda }_j{\displaystyle \int }\mathrm{d}\mathit{x}{\mathcal{W}}_j(\mathit{x})(1+\tau (\mathit{x})/2),\end{array}$$(9)

where τ(x) is the λ_i-dependent implicit solution of the consistency equation

$$\begin{array}{c}\hfill \tau (\mathit{x})={\displaystyle \sum _j}{\lambda }_j{\displaystyle \int }\mathrm{d}\mathit{x}\prime {\mathcal{W}}_j(\mathit{x}\prime )\xi (\mathit{x},\mathit{x}\prime )(1+\tau (\mathit{x}\prime )/2).\end{array}$$(10)

Formally, τ(x) is actually the rooted-Cumulant Generating Function (the generating function of cumulants that originate from location x) and is denoted r-CGF hereafter. It will play an important role throughout this paper. The expression in Eq. (11) is the direct transcription of the tree model described by Eq. (8). It is the result of combinatoric computations, and does not rely on any approximation². The practical implementation of Eq. (10) relies on some approximations, such as the mean field assumption, in which the r-CGF is assumed to remain constant within each cell (as detailed in Appendix B). This proves to be highly accurate for compact, non-compensated, spherically symmetric density profiles. Assuming it holds, we denote τ_i as the value of the r-CGF for each cell, i. We can then average Eq. (11) over the cell, i, with the profile 𝒲_i(x), which leads to the following system

$$\begin{array}{c}\hfill {\tau }_i={\displaystyle \sum _j}{\lambda }_j(1+{\tau }_j/2){\displaystyle \int }\mathrm{d}\mathit{x}{\displaystyle \int }\mathrm{d}\mathit{x}\prime {\mathcal{W}}_i(\mathit{x}){\mathcal{W}}_j(\mathit{x}\prime )\xi (\mathit{x},\mathit{x}\prime ).\end{array}$$(11)

This becomes a set of n equations coupling the different values of τ_i. As can be seen here, in case of the present minimal tree model, this system is linear in τ_i. It can therefore be explicitly inverted³.

For the one-point cumulant generating function, Eq. (12) corresponds to the implicit equation:

$$\begin{array}{c}\hfill {\tau }_1={\lambda }_1{\xi }_0(1-{\tau }_1/2),\end{array}$$(12)

where ξ₀ is the average two-point correlation function within one cell:

$$\begin{array}{c}\hfill {\xi }_0={\displaystyle \int }\mathrm{d}\mathit{x}\mathrm{d}\mathit{x}\prime {\mathcal{W}}_1(\mathit{x}){\mathcal{W}}_1(\mathit{x}\prime )\xi (\mathit{x},\mathit{x}\prime ).\end{array}$$(13)

The resulting CGF can easily be computed (See Appendix C):

$$\begin{array}{c}\hfill {\phi }_1({\lambda }_{\rho };{\xi }_0)=\frac{2{\lambda }_{\rho }}{2-{\xi }_0{\lambda }_{\rho }}.\end{array}$$(14)

Despite its apparent simplicity, this expression fully captures the whole cumulant hierarchy that we would expect for the density PDF. Furthermore, it turns out that the inverse Laplace transform of such an expression can be explicitly derived, see Bernardeau (2022) and Appendix C, and leads to the closed form expression of the density PDF:

$$\begin{array}{cc}\hfill P(\rho )& ={\displaystyle {{\displaystyle \int }}_{-i\infty }^{+i\infty }}\frac{\mathrm{d}{\lambda }_{\rho }}{2\pi i}exp(-{\lambda }_{\rho }\rho +{\phi }_1({\lambda }_{\rho }))\hfill \end{array}$$(15)

$$\begin{array}{cc}& =e^{-2/{\xi }_0}{\delta }_d(\rho )+\frac{4}{{\xi }_0^2}{e}^{-\frac{2(\rho +1)}{{\xi }_0}}{}_0{\stackrel{\sim }{F}}_1(2;\frac{4\rho }{{\xi }_0^2}).\hfill \end{array}$$(16)

The first (singular) term in Eq. (17) reflects the contribution to empty regions. The last terms involves the regularised confluent hypergeometric function, ${}_0{\stackrel{\sim }{F}}_1(a;z)$, given by $x^{\frac{1-a}{2}}{I}_{a-1}\left(2\sqrt{x}\right)$. We note that this model indeed predicts region that are empty even for an arbitrarily large number of points⁴. More specifically, the void probability function (VPDF) is non-zero, even in the continuous limit, and is given by

$$\begin{array}{c}\hfill P_0=e^{-2/{\xi }_0}.\end{array}$$(17)

This feature can be appreciated on Fig. 2, where we can see empty regions that cover a large fraction of the sample. The VPDF can only be non-zero for a compact support filter, which not formally holds for a Gaussian filter. Yet the mean field result seems to accurately predict the VPDF for a top-hat filter. For a Gaussian filter, the behaviour of the density PDF in the low density regime is correct when corrections to the mean field solutions are included.

We note that for the mean-field solution, the dependence on the spatial dimension, the value of α and the filter shape is entirely contained in the value of ξ₀. This is not the case for derivations of the density PDF beyond the mean-field solution, which are expected to depend on the details of the model and filtering schemes. Numerical investigations beyond the mean field show nonetheless that the high density tail of Eq. (17) is very robust. This is illustrated on Fig. 3, which displays the mean-field solution and its corrections in terms of Hermite polynomials (up to sixth order in the expression of the r-CGF, see Appendix C) and compares it to numerical results.

3. Critical point number counts

From a given realisation of the flight, a convolution by, for instance, a Gaussian filter allows us to define a smooth field, ρ, whose critical point can be studied statistically. The number density of critical points at a density, ρ, is then expressed as:

$$\begin{array}{c}\hfill n_{\mathrm{crit}.}(\rho )={\displaystyle \int \prod _{\mathit{ij}}}\mathrm{d}{\rho }_{,ij}\mathrm{Sgn}\left[{\rho }_{,ij}\right]det\left[{\rho }_{,ij}\right]P(\rho ,{\rho }_{,i}=0,{\rho }_{,ij}),\end{array}$$(18)

where Sgn[ρ_,ij] is either 0, 1 or −1, depending on the sign of the eigenvalues of the matrix, ρ_,ij, to reflect the nature of the critical points considered. Hence, the computation of extrema densities requires the joint cumulant generating function of variables dual to the local density, its first and second order derivatives. Appendix D shows that it is indeed possible to derive such a joint CGF in the mean field limit. In short, the calculation is based on writing the mean field solution for a finite number of cells assumed to be infinitely close to one another; the derivatives are then obtained via finite differences.

Once the cumulant generating function for the field and its derivatives is known, the relevant conditional expectations are computed to predict the extrema and critical number counts and their clustering properties. We present below the results for the Euler and critical points in 1D–3D. As Eq. (19) involves up to the second derivative of the field, in principle, we should compute the joint PDF of the field and its derivative up to that order. In practice, For 1D Euler number counts, Appendix D shows explicitly that thanks to the stationarity of the field, only the first derivative is necessary, and we can write:

$$\begin{array}{c}\hfill n_{1\mathrm{D}}^{\mathrm{Euler}}(\rho )=\frac{\partial }{\partial \rho }{\displaystyle {{\displaystyle \int }}_{-\infty }^0}{\rho }_{,x}\mathrm{d}{\rho }_{,x}P(\rho ,{\rho }_{,x}).\end{array}$$(19)

Overall, and after a significant amount of algebra, we derive in Appendix D the closed form Euler counts in ND:

$$\begin{array}{cc}\hfill n_{\mathrm{ND}}^{\mathrm{Euler}}(\rho )& =\frac{1}{{\xi }_0^2}{\left(\frac{2{\xi }_1}{\pi {\xi }_0^2\rho }\right)}^{N/2}{e}^{-\frac{2(\rho +1)}{{\xi }_0}}\hfill \\ \hfill & \times \{P_{N-1}({\xi }_0,\rho ){}_0{\stackrel{\sim }{F}}_1(1;\frac{4\rho }{{\xi }_0^2})+Q_N({\xi }_0,\rho ){}_0{\stackrel{\sim }{F}}_1(2;\frac{4\rho }{{\xi }_0^2})\},\hfill \end{array}$$(20)

where P_m and Q_m are polynomials of order m in ρ and ξ₁ is defined in Eq. (D.1). Specifically we find in Appendix D that

$$\begin{array}{cc}\hfill P_0& =2{\xi }_0,P_1=-{\xi }_0({\xi }_0+8\rho )/2,\hfill \end{array}$$(21a)

$$\begin{array}{cc}\hfill P_2& =(16{\xi }_0\rho (3\rho +1)+3{\xi }_0^3)/8,\hfill \end{array}$$(21b)

$$\begin{array}{cc}\hfill Q_1& =-({\xi }_0+4\rho ),Q_2=(2{\xi }_0\rho +{\xi }_0^2+8\rho (\rho +1))/2,\hfill \end{array}$$(21c)

$$\begin{array}{cc}\hfill Q_3& =-(6{\xi }_0^2\rho +16{\xi }_0\rho +3{\xi }_0^3+32{\rho }^2(\rho +3))/8.\hfill \end{array}$$(21d)

Here, ${\overline{\xi }}_{\mathrm{D}}\propto {\ell }_0^{-\alpha }{r}_s^{\alpha -D}$ where the pre-factor defining ξ₀ is a function of α only, given by Eqs. (A.7)–(A.9) for 1D–3D.

The simplicity of Eq. (21) is quite remarkable. We note, in particular, that for D = 2 to 3, $n_{\mathrm{ND}}^{\mathrm{Euler}}(\rho )$ also does not involve ξ₂ (see the discussion below and Appendix D).

We could speculate that relations, such as Eqs. (3)a–d, hold in higher dimensions, involving polynomials, P_i and Q_i of increasing order, so that Eq. (21) mirrors the hierarchy for Gaussian random field Euler counts (involving Hermite polynomials, Adler & Taylor 2009), which follow from their cumulant generating functions given in Appendix F.

It is of interest to be able to distinguish between the number counts of maxima and minima separately. It turns out that the extrema counts can also be computed via specific (simple or double) numerical integration path in the complex plane as:

$$\begin{array}{cc}\hfill n_{1\mathrm{D}}^{\pm }(\rho )& ={\displaystyle {{\displaystyle \int }}_{-i\infty \pm ϵ}^{+i\infty \pm ϵ}}\frac{\mathrm{d}{\lambda }_{\mathit{xx}}}{2\pi i}\frac{1}{{\lambda }_{\mathit{xx}}^2}{\chi }_{1\mathrm{D}}(\rho ;{\lambda }_{\mathit{xx}}),\hfill \end{array}$$(22)

and

$$\begin{array}{c}\hfill n_{2\mathrm{D}}^{\pm }(\rho )={\displaystyle {{\displaystyle \int }}_{-i\infty \pm ϵ}^{+i\infty \pm ϵ}}\frac{\mathrm{d}{\lambda }_{\kappa }}{2\pi i}{\displaystyle {{\displaystyle \int }}_0^{\infty }}\frac{\mathrm{d}{\lambda }_{\gamma }}{2\pi }\\ \hfill \frac{12\pi {\lambda }_{\gamma }}{{({\lambda }_{\gamma }^2+{\lambda }_{\kappa }^2)}^{5/2}}{\chi }_{2\mathrm{D}}(\rho ;{\lambda }_{\kappa },{\lambda }_{\gamma }),\end{array}$$

where χ_1D(ρ; λ_xx) and χ_2D(ρ; λ_κ, λ_γ) are given by Eqs. (D.34) and (D.51) respectively, while ϵ is a negative real constant for maxima and a positive real constant for minima. The 2D saddle counts follows from ${\displaystyle \sum _{\pm }}{n}_{2\mathrm{D}}^{\pm }(\rho )-n_{2\mathrm{D}}^{\mathrm{Euler}}(\rho )$.

To gain a bit of insight into these quantities, we consider the expression of $n_{1\mathrm{D}}^{\pm }(\rho )$ as a function of $\gamma ={\xi }_1/\sqrt{{\xi }_2{\xi }_0}$ in the limit⁵γ → 0. The integrand in this case is greatly simplified and it leads to the same leading behaviour for the number density of maxima and minima given by:

$$\begin{array}{c}\hfill n_{\mathrm{ext}}(\rho )\mathrm{d}\rho =\frac{1}{2\pi }{\left(\frac{{\xi }_2}{{\xi }_1}\right)}^{1/2}P(\rho )\mathrm{d}\rho ,\end{array}$$(23)

which, as expected, scales like the inverse of typical distance between extrema, (ξ₁/ξ₂)^1/2. The next to leading order is independent of γ, and corresponds to the Euler number density so that

$$\begin{array}{c}\hfill n_{1\mathrm{D}}^{\pm }(\rho )\mathrm{d}\rho =\frac{1}{2\pi }{\left(\frac{{\xi }_2}{{\xi }_1}\right)}^{1/2}P(\rho )\mathrm{d}\rho \pm \frac{1}{2}{n}_{1\mathrm{D}}^{\mathrm{Euler}}(\rho )\mathrm{d}\rho +\cdots ,\end{array}$$(24)

where the dots represents subsequent terms in a γ expansion. Here n_Euler is less sensitive to small scale high frequency features – or noise – compared to extrema number density. This is also the case for Gaussian field and for the Rayleigh-Levy flights in any dimension. Although we cannot demonstrate it, we speculate that this is true for all hierarchical models in which the high-order correlation function behave like products of the two-point correlations. Figures 4 and 5 present the corresponding number counts in 1D and 2D, respectively. In Fig. 4, in particular, we present a detailed comparison of the extrema, minima and Euler number densities as a function of the local density. We take advantage of the fact that we could rely on large number of realisations making the measurements of the critical points rather precise (numerical error bars due to the finite number of realisations are here negligible). For the 1D case, it is also rather simple to determine the position and type of critical points: they are obtained at locations where the local gradients (obtained after convolution on grid points shifted by 1/2 pixel) change sign from one pixel to the next. The type of critical point is simply given by the sign of the difference. Figure 4 shows that the theoretical predictions capture in exquisite details the measured densities of critical points. The only significant departure between the theoretical predictions and the numerical results is for the low density regime (ρ < 1), for the maxima number densities and to a less extent the minima number densities. These discrepancies are thought to be related to the breaking of the mean field approximation in the low density regime, as shown on Fig. 3 for the density PDF. All these conclusions are found to be valid irrespective of the details of the simulations (this is illustrated with a change of the value of α in the Rayleigh-Levy flights).

4. Clustering of critical points

One of the strength of analytical investigations of tree-hierarchical models is that the same re-summation techniques can be used to also infer the large distance correlation functions of the objects for which we can compute the number density. These calculations were pioneered in Bernardeau & Schaeffer (1999) for tree hierarchical models, applied also to the perturbation theory calculations in Bernardeau (1996), and more recently in Codis et al. (2016), Uhlemann et al. (2017), Munshi (2018), Repp & Szapudi (2021). A general thorough presentation of large-scale biasing is also to be found in the review paper by Desjacques et al. (2018).

The properties of large scale biasing were more specifically investigated in Bernardeau (2022) for the Rayleigh-Levy model, which focused on the consequences of this functional form for the covariance properties of density PDF measurements. These models share the same properties: the joint density PDF at large separation are expected to obey the following functional form

$$\begin{array}{c}\hfill P({\rho }_1,{\rho }_2)=P({\rho }_1)P({\rho }_2)[1+\xi (r_1,r_2)b({\rho }_1)b({\rho }_2)],\end{array}$$(25)

where the densities ρ₁ and ρ₂ are taken at position r₁ and r₂ respectively. At this order, the dependence in the densities ρ₁ and ρ₂ factorises, making it possible to define a density bias function. In this paper, we present how such a relation can be extended to the number counts of critical points and derive the corresponding bias functions.

In Eq. (29) the bias function can be seen as the response function of the density PDFs to a change of the global density. This means that although the bias function cannot be derived from the density PDFs alone⁶, we should be able to derive it if an operational method is available to compute the density PDFs for arbitrary large-scale density (following the derivation of halo-bias function as pioneered by Mo & White (1996) in a separate universe approach). In hierarchical tree models, it can be computed making use of the r-CGF. The shape of the join density CGF at for densities at positions r₁ and r₂ written up to first order in ξ(r₁, r₂)/ξ₀ indeed obeys

$$\begin{array}{c}\hfill \phi ({\lambda }_1,{\lambda }_2)={\phi }_1({\lambda }_1){\phi }_1({\lambda }_2)+{\tau }_1({\lambda }_1)\xi (r_1,r_2){\tau }_1({\lambda }_2),\end{array}$$(26)

where τ₁ is the r-CGF introduced in Eq. (11). The computation of the inverse Laplace transform of such a join CGF can be formally done at leading order in ξ(r₁, r₂) and has the functional form given by Eq. (29). More specifically, equation (16) is to be extended to

$$\begin{array}{c}\hfill b(\rho )P(\rho )={\displaystyle {{\displaystyle \int }}_{-i\infty }^{+i\infty }}\frac{\mathrm{d}{\lambda }_{\rho }}{2\pi i}{\tau }_1({\lambda }_{\rho })exp(-{\lambda }_{\rho }\rho +{\phi }_1({\lambda }_{\rho }))\\ \hfill \end{array}$$(27)

for the bias function. For the Rayleigh-Levy flights model, the large scale density bias function can be derived explicitly (Bernardeau 2022) and it is expressed as:

$$\begin{array}{c}\hfill b(\rho )=\frac{{}_0{\stackrel{\sim }{F}}_1(1;\frac{4\rho }{{\xi }_0^2})}{{}_0{\stackrel{\sim }{F}}_1(2;\frac{4\rho }{{\xi }_0^2})}-\frac{2}{{\xi }_0}.\\ \hfill \end{array}$$(28)

For a high density, its limiting behaviour is:

$$\begin{array}{c}\hfill b_{\mathrm{large}}(\rho )=\frac{1}{4}+\frac{2}{{\xi }_0}(\sqrt{\rho }-1),\end{array}$$(29)

which corresponds to a dependence in $\sqrt{\rho }/{\xi }_0$ for high densities. We note that this is not a priori a generic result for hierarchical models and is at variance with the behaviour found in the context of Perturbation Theory calculations in Codis et al. (2016), where the bias at a high density is expected to be proportional to the density⁷.

The objective of this section is to present the correlation properties of the critical points of Rayleigh-Levy flights. Correlation properties of critical points will depend in general on i) the type of critical points we are considering and ii) the distance between the points and how this distance compares to the scale at which these points are defined. Despite the fact that the model we consider is unambiguously defined, it is near impossible to derive the correlation properties of critical points in their full complexity. Results can be derived at large enough separation and using the mean field approach, implying that our findings will be solid for large enough densities only.

The starting point is a generalisation of Equation (30) for the joint CGF of the density and its derivatives:

$$\begin{array}{cc}\hfill \phi ({\lambda }_{\rho },\{{\lambda }_i\},\{{\lambda }_{\mathit{ij}}\};{\mu }_{\rho },\{{\mu }_i\},\{{\mu }_{\mathit{ij}}\})& =\phi ({\lambda }_{\rho },\{{\lambda }_i\},\{{\lambda }_{\mathit{ij}}\})+\phi ({\mu }_{\rho },\{{\mu }_i\},\{{\mu }_{\mathit{ij}}\})\hfill \\ \hfill & +\tau ({\lambda }_{\rho },\{{\lambda }_i\},\{{\lambda }_{\mathit{ij}}\}){\xi }_0(d)\tau ({\mu }_{\rho },\{{\mu }_i\},\{{\mu }_{\mathit{ij}}\}),\hfill \end{array}$$(30)

where λ_ρ, {λ_i},{λ_ij} are the conjugate variables of the density and its derivatives at a given location, μ_ρ, {μ_i},{μ_ij} those at a second location placed at distance, d, from one another. ξ₀(d) is the filtered matter density correlation function at distance, d. This expansion is valid when d is much larger than the filtering scale and derived from the fact that we then expect ξ₀(d)≪ξ₀ and ξ₂(d)/ξ₂ ≪ ξ₀(d)/ξ₀ (see Appendix E). For the Rayleigh-Levy model we further have τ = ϕ. As for the density field this form implies a similar form for the two-point number density of the critical points:

$$\begin{array}{cc}\hfill n_{\#1,\#2}({\rho }_1;{\rho }_2)& =n_{\#1}({\rho }_1)n_{\#2}({\rho }_1)\{1+b_{\#1}({\rho }_1){\xi }_0(d)b_{\#2}({\rho }_2)\},\hfill \end{array}$$(31)

where #i represents the types of critical points we consider and b_#i(ρ) is the associated bias function whose expression is given by Eq. (E.7).

The computation of bias function of the critical points then makes use of the same techniques as for the number density. It does not lead to an explicit form for the extrema or saddle point position both in 1D and 2D. It is however possible to derive the explicit bias functions of the Euler points for the 1D and 2D case:

$$\begin{array}{cc}\hfill b_{\mathrm{Euler}}^{\mathrm{ND}}(\rho )n_{\mathrm{Euler}}^{\mathrm{ND}}(\rho )& =\frac{1}{{\xi }_0^3}{\left(\frac{2{\xi }_1}{\pi {\xi }_0^2\rho }\right)}^{N/2}{e}^{-\frac{2(\rho +1)}{{\xi }_0}}\hfill \\ \hfill & \times \{{\mathcal{P}}_N({\xi }_0,\rho ){}_0{\stackrel{\sim }{F}}_1(1;\frac{4\rho }{{\xi }_0^2})+{\mathcal{Q}}_N({\xi }_0,\rho ){}_0{\stackrel{\sim }{F}}_1(2;\frac{4\rho }{{\xi }_0^2})\}.\hfill \end{array}$$(32)

where the polynomials obey

$$\begin{array}{cc}\hfill {\mathcal{P}}_1& ={\xi }_0({\xi }_0-4(1+\rho )),{\mathcal{Q}}_1=2({\xi }_0+8\rho ),\hfill \end{array}$$(33a)

$$\begin{array}{cc}\hfill {\mathcal{P}}_2& ={\xi }_0({\xi }_0-3{\xi }_0\rho +4\rho (3+\rho )),{\mathcal{Q}}_2=-({\xi }_0^2+8\rho (1+3\rho )).\hfill \end{array}$$(33b)

Again one could conjecture that higher order polynomials exist for bias functions in higher dimensions. We note that bias functions in general are dimensionless quantities. They can therefore be expressed in terms of dimensionless quantities such as ξ₀ and γ. Furthermore the bias parameter of Euler points is found to be also independent of ξ₂ (as was the case for its number density) and therefore of γ.

The corresponding bias functions are shown in Fig. 6 and 7 in 1D and 2D, respectively. These functions are computed for ξ₀ = 1 and for γ = 0.45, which corresponds roughly to what is expected for α = 0.5. The results however do not depend crucially on these specific choices.

The functions exhibit some expected generic behaviours: (i) the high-density asymptotic behaviour of maxima and Euler points are the same. It also matches the asymptotic form of the local density bias: high-density regions tend to trace the locations of the maxima; (ii) for a given density, minima tend to be more correlated than saddle points or maxima: imposing to have a minima within a high-density region indeed requires an even higher density in the surroundings, which in turn translates into larger biases. One important conclusion that can be drawn here is that peak correlation properties derive, to a large extent, for the behaviour of the density bias functions in the sense of Eq. (29). This justifies to use such a proxy for the computation of peak correlation in the high-density regime for more involved models. Such an approach, however, cannot capture proximity or exclusion effects that would be expected for peak correlations.

To explore these effects specifically, we rely here on results from numerical experiments of 1D Rayleigh-Levy flights, for which accurate measurements can be achieved. The corresponding curves are presented on Figs. 8–12. To avoid too much numerical noise and uncertainties, they are expressed in terms of correlation function of cumulative quantities, where the local density is above a given threshold. More precisely the bias factors we will apply are defined as:

$$\begin{array}{c}\hfill b_{\#}(>{\rho }_s)=\frac{{\displaystyle {{\displaystyle \int }}_{{\rho }_s}^{\infty }}\mathrm{d}\rho b_{\#}(\rho )n_{\#}(\rho )}{{\displaystyle {{\displaystyle \int }}_{{\rho }_s}^{\infty }}\mathrm{d}\rho n_{\#}(\rho )}.\end{array}$$(34)

The 1D correlation functions are directly computed in position space by applying a simple shift to the density fields. No binning is applied and they are measured to a separation of 1000 pixels, which is safely smaller than the sample size.

Fig. 8 shows the correlation functions of the thresholded density for different thresholds. The dashed lines are the theoretical predictions, namely, the expectation when one multiplies the measured matter correlation function of the filtered density field (grey solid line) with b(> ρ_s)². One can see that proximity effects tend to be rather insignificant for low density threshold, but rather large for high density threshold: for ρ_s = 4 proximity effects can be detected up to a distance of about 100 pixels. We can also predict the actual shape of the two-point function for Rayleigh-Levy flights models. This is presented in Fig. 1.

Figures 8–11 display their corresponding clipped two-point functions, while Fig. 12 shows the cross correlation of minima and maxima. The consistency of the results with theoretical predictions show that on large scales, the functional form of the peak biases, Eq. (35), is indeed satisfied, and the auto- and cross correlation properties of critical points can be described with such a factorised form.

5. Discussion and conclusions

Rayleigh-Levy flights can serve as numerical models portraying matter distribution within highly nonlinear fields in cosmology. In contrast to standard Markovian processes, Rayleigh-Levy flights display long-range correlations across all orders, which are precisely understood. Thus, these models serve as simplified representations, capturing statistical properties of the cosmic density field after gravitational instabilities reach its full development, although we are well aware they do not precisely replicate the outcomes of these instabilities. Despite this, they exhibit key properties, such as hierarchical structure in the higher correlation functions, making their study a valuable reference point.

We wish to emphasise that Rayleigh-Levy flights differ from simple nonlinear transformations (like a lognormal fields introduced by Coles & Jones 1991, which are now widely used as toy models), as they exhibit hierarchical properties across scales. This makes them the only known example (at least to us) of an actual random process whose outcome displays the expected scaling properties in correlation functions. Although it exhibits such a non-trivial structure, it is yet simple enough to allow for a wide range of explicit predictions.

Utilising the mean field approximation specifically enables the derivation of closed analytical forms for crucial parameters, such as Euler number densities and their correlations across significant separations: Eqs. (21) and (36), along with corresponding quadratures for extreme value counts, represent the primary outcomes of this study. Importantly, our predictions hold irrespective of the density field’s variance. While aligning with expectations for a Gaussian field in the regime of small variances, our numerical experiments confirm their validity across all variance values. Several key insights emerge from our findings:

1. The critical point density depends⁸ on three quantities related to density variance, gradient variance, and variance of the second-order derivative, denoted as ξ₀, ξ₁, and ξ₂, respectively. These can be re-expressed using two dimensionless parameters, ξ₀ and $\gamma ={\xi }_1/\sqrt{{\xi }_0{\xi }_2}$, along with ξ₁, inversely proportional to distance squared. This implies that number densities scale as ${\xi }_1^{N/2}$ in N dimensions. Remarkably, while extrema number counts are correlated with both ξ₀ and γ, Euler number counts remain independent of γ. This hints at the impact of small-scale features on peak counts, creating pairs of peaks and troughs at roughly the same height, ultimately leaving Euler number densities unaffected. We speculate that this property may extend to all hierarchical models.

2. Concerning point correlations, akin to matter fields, we find that critical points of all types exhibit a factorised structure in the large separation limit, as given in Eq. (35). This pattern resembles the separate universe approach and has been validated through our numerical experiments. It enables the derivation of bias parameters for various types of points, especially as given in Eq. (36).

3. Notably, bias parameters of maxima, Euler points, and threshold-limited density fields converge to the same limit in the high-density limit. This suggests that regions with rare high-density values likely host one extremum, sharing identical correlation properties in this limit. This supports simplified approaches where peak correlations mirror those of high-density regions in more complex systems (see e.g. Bernardeau & Schaeffer 1999).

Beyond the paper’s scope, extending Eqs. (21), (26), and (36) to higher dimensions, other Minkowski functionals, and cross-validating them via Monte Carlo techniques, as done in 1D, will prove valuable and will be the topic of upcoming papers. One interesting line of investigation is for instance to use these results as benchmarks for validating quasi-Gaussian expansion schemes. It could also be contrasted with those derived in the large deviation limit (as presented in Uhlemann et al. 2016). Furthermore, these statistics may find application in analyzing the influence of Rayleigh-Levy flights in anomalous diffusive processes (e.g. in supersonic turbulence, Colbrook et al. 2017) or finite star effects in gravitating systems (Chavanis 2009).

Given their common occurrence, the comprehensive findings from this study hold promise not just in astrophysics but also in broader applications beyond the field.

Data availability

The data underlying this article will be shared on reasonable request to the corresponding author. Some of the codes underpinning this paper are available on Github at the following URL: https://github.com/cncpichon/levyflight

---

Acknowledgments

We warmly thank S. Appleby for feedback and numerical checks in 2 and 3D. We also thank D. Pogosyan for comments and for updating his extrema count codes to suit our purpose, T. Abel for interesting feedback and Patrick Peter for typographic advice. We are grateful for the KITP for hosting the workshop https://www.cosmicweb23.org during which this project was advanced. This work is partially supported by the National Science Foundation under Grant No. NSF PHY-1748958. This work has made use of the Horizon cluster hosted by the Institut d’Astrophysique de Paris. We also thank S. Rouberol for running it smoothly.

References

Appendix A: Rayleigh-Levy flight model

The random Rayleigh-Levy flights model was introduced first (in a cosmological context) by J. Peebles (Peebles 1980). This is a Markov random walk where the PDF of the step length ℓ follows a cumulative cumulative distribution function given by Eq. (2), where ℓ₀ is a small-scale regularisation parameter. As a result, the density of the subsequent point (first descendant) at position r is given by:

$$\begin{array}{cc}\hfill f_1(\mathit{r})& =\frac{\alpha }{2}\frac{{\ell }_0^{\alpha }}{|\mathit{r}-{\mathit{r}}_0{|}^{1+\alpha }}\mathrm{in}1\mathrm{D}\mathrm{space};\hfill \end{array}$$(A.1a)

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

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

Defining ψ(k) as the Fourier transform of f₁(r) with Eq. (4), we can easily show that when assuming there are an infinity of points in the flight, the density in r of descendants of the point at position r₀ is given by a series of convolutions whose resummation is:

$$\begin{array}{c}\hfill f({\mathit{r}}_0,\mathit{r})={\displaystyle \int \frac{{\mathrm{d}}^{\mathrm{D}}\mathit{k}}{{(2\pi )}^{\mathrm{D}}}}{e}^{i\mathit{k}.(\mathit{r}-{\mathit{r}}_0)}\frac{1}{1-\psi (k)}.\end{array}$$(A.2)

The two-point density correlation function is then given by two possible configurations: a neighbourg 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}\hfill {\xi }_2({\mathit{r}}_1,{\mathit{r}}_2)& =\frac{1}{n}[f({\mathit{r}}_1,{\mathit{r}}_2)+f({\mathit{r}}_2,{\mathit{r}}_1)],\hfill \end{array}$$(A.3)

where n is the number density of points in the sample which leads to Eq. (5) of the main text. We note that the two-point correlation within the sample is therefore scaled as 1/n. The latter can be associated with a typical length, ℓ_n,

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

At large distance (compared to ℓ₀), we expect the two-point correlation to behave as power laws. They are given by

$$\begin{array}{cc}\hfill {\xi }_2^{1\mathrm{D}}(r)& =\frac{2tan\left(\frac{\pi \alpha }{2}\right)}{\pi }{r}^{\alpha -1}{\ell }_0^{-\alpha }{\ell }_n,\hfill \end{array}$$(A.5a)

$$\begin{array}{cc}\hfill {\xi }_2^{2\mathrm{D}}(r)& =\frac{\alpha }{\pi }{r}^{\alpha -2}{\ell }_0^{-\alpha }{\ell }_n^2,\hfill \end{array}$$(A.5b)

$$\begin{array}{cc}\hfill {\xi }_2^{3\mathrm{D}}(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}$$(A.5c)

For practical purposes, we give here the resulting expression of the average correlation, ξ₀, for a Gaussian window function of width r_s,

$$\begin{array}{cc}\hfill {\overline{\xi }}_G^{1\mathrm{D}}(r_s)& =\frac{2sin\left(\frac{\pi \alpha }{2}\right)\mathrm{\Gamma }(\frac{1}{2}-\frac{\alpha }{2})\mathrm{\Gamma }(\alpha )r_0^{-\alpha }{r}_s^{\alpha -1}{\ell }_n}{{\pi }^2},\hfill \end{array}$$(A.6a)

$$\begin{array}{cc}\hfill {\overline{\xi }}_G^{2\mathrm{D}}(r_s)& =\frac{2^{\alpha -2}\alpha \mathrm{\Gamma }\left(\frac{\alpha }{2}\right)r_0^{-\alpha }{r}_s^{\alpha -2}{\ell }_n^2}{\pi },\hfill \end{array}$$(A.6b)

$$\begin{array}{cc}\hfill {\overline{\xi }}_G^{3\mathrm{D}}(r_s)& =\frac{(\alpha +1)sin\left(\frac{\pi \alpha }{2}\right)\mathrm{\Gamma }(\frac{3}{2}-\frac{\alpha }{2})\mathrm{\Gamma }(\alpha )r_0^{-\alpha }{r}_s^{\alpha -3}{\ell }_n^3}{{\pi }^3}.\hfill \end{array}$$(A.6c)

The main interest of this model however lies in the fact that its higher order correlation functions can also be computed and that they take a simple form. The reason is that 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 }_3({\mathit{r}}_1,{\mathit{r}}_2,{\mathit{r}}_3)=\frac{1}{n^2}[f({\mathit{r}}_1,{\mathit{r}}_2)f({\mathit{r}}_2,{\mathit{r}}_3)+...],\end{array}$$(A.7)

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

$$\begin{array}{cc}\hfill {\xi }_3({\mathit{r}}_1,{\mathit{r}}_2,{\mathit{r}}_3)& =\frac{1}{2}[{\xi }_2({\mathit{r}}_1,{\mathit{r}}_2){\xi }_2({\mathit{r}}_2,{\mathit{r}}_3)+\hfill \\ \hfill & {\xi }_2({\mathit{r}}_2,{\mathit{r}}_3){\xi }_2({\mathit{r}}_3,{\mathit{r}}_1)+{\xi }_2({\mathit{r}}_3,{\mathit{r}}_1){\xi }_2({\mathit{r}}_1,{\mathit{r}}_2)].\hfill \end{array}$$(A.8)

The type of expansion can be pursued to any order. The resulting shapes of the p-point correlation function is the following:

$$\begin{array}{cc}\hfill {\xi }_p(& {\mathit{r}}_1,\cdots {\mathit{r}}_p)=\frac{1}{n^{p-1}}{\displaystyle \sum _{\sigma \in \mathrm{perm}[p]}}f({\mathit{r}}_{{\sigma }_1},{\mathit{r}}_{{\sigma }_2})\cdots f({\mathit{r}}_{{\sigma }_{n-1}},{\mathit{r}}_{{\sigma }_n}),\hfill \end{array}$$

where σ is any permutation of the p indices 1, …, p. It implies that the p-point correlation function can be expressed in terms of the two-point functions

$$\begin{array}{cc}\hfill {\xi }_p(& {\mathit{r}}_1,\cdots {\mathit{r}}_p)={\left(\frac{1}{2}\right)}^{p-2}{\displaystyle \sum _{\sigma \in {\mathrm{perm}}^{\ast }[p]}}{\xi }_2({\mathit{r}}_{{\sigma }_1},{\mathit{r}}_{{\sigma }_2})\cdots {\xi }_2({\mathit{r}}_{{\sigma }_{n-1}},{\mathit{r}}_{{\sigma }_n}),\hfill \end{array}$$

where the exponent refers to the subset of permutations that lead to a unique un-oriented sequence; identifying, for instance, Eqs. (1, 2, 3) and (3, 2, 1)). As we see hereafter it corresponds to a specific hierarchical tree model.

Appendix B: Hierarchical tree models

The Rayleigh-Levy flight model is one representative of a large class of models, the so-called hierarchical tree models, that have been put forward in cosmology as a way to model the density field in the highly nonlinear regime. Such models are presented in detail in Bernardeau & Schaeffer (1999). We recall here how they are defined, together with the basic equation that allows the derivation of their cumulant-generating function. Hierarchical tree models are a general class of non-Gaussian fields whose n-point correlation functions, ξ_p(r₁, …, r_p), follow the so-called hierarchical Ansatz,

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

where ξ(r_i, r_j) is the two-point function, while the sum is made over all possible trees that join the p points (diagram without loops), and the tree value, Q_p, is obtained by the product of a fixed weight (that depends only on the tree’s topology), and the product of the two-point correlation functions, for all pairs that are connected together in the given tree. More specifically it is assumed that:

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

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

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

Such models are thus entirely defined by i) the two-point functions ξ(r) and ii) the vertex-generating function ζ(τ). What the previous section has shown is that the Rayleigh-Levy flight model (in the continuous limit) corresponds, in the continuous limit, to a specific hierarchical tree model⁹ with ν₂ = 1/2, ν_q = 0 for q ≥ 3.

B.1. Expression for cumulant generating function

The exact generating function of multiple cell correlation functions can be built through simple transforms. We hereafter consider a set of n cells i of profile 𝒲_i(x), thus allowing for generic types of profiles, and not only top-hat boxes as was assumed in early papers. These profiles can obviously overlap but they are assumed to be well localised. The joint cumulants we consider in this formalism are those of the average densities in cells 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 }\mathrm{d}{\mathit{x}}_{1,1}{\mathcal{W}}_1({\mathit{x}}_{1,1})\cdots {\displaystyle \int }\mathrm{d}{\mathit{x}}_{1,p_1}{\mathcal{W}}_1({\mathit{x}}_{1,p_1})\cdots \hfill \\ \hfill & \cdots {\displaystyle \int }\mathrm{d}{\mathit{x}}_{n,1}{\mathcal{W}}_n({\mathit{x}}_{n,1})\cdots {\displaystyle \int }\mathrm{d}{\mathit{x}}_{n,p_n}{\mathcal{W}}_n({\mathit{x}}_{n,p_n})\hfill \\ \hfill & \times {\xi }_p({\mathit{x}}_{1,1},\cdots ,{\mathit{x}}_{1,p_1},\cdots {\mathit{x}}_{n,1},\cdots ,{\mathit{x}}_{n,p_n}).\hfill \end{array}$$(B.4)

We then wish to build the cumulant-generating function (CGF) of the densities in each cells. It is given by:

$$\begin{array}{cc}\hfill \phi ({\lambda }_1,\cdots ,{\lambda }_n)& =log⟨e^{\sum {\lambda }_i{\rho }_i}⟩,\hfill \end{array}$$(B.5)

$$\begin{array}{cc}& ={\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}}{p_n!}.\hfill \end{array}$$(B.6)

This function represents the generating function of (averaged) cumulants, where the power of the λ_i counts the number of points in each cells. When the cumulants follow the tree structure described in the previous paragraph, each term that appear in this function then corresponds to a specific tree. Following a method pioneered by Janninck (1987), Bernardeau & Schaeffer (1992) showed that it is actually possible to perform the summation over such a set of trees, leading to a formal expression we re-derive briefly below.

To do so, we first now define a larger class of objects, ϕ[ν_p(x),ξ(x, x′)] (note the different symbol), which represents the sum of trees with vertices of order p having the weight ν_p(x) and lines having the weight ξ(x, x′), and where x are space variables that are subsequently integrated in the whole domain. In the context of tree theory, a point that reaches a one-point vertex ν₁(x) are usually called a ’leaf’. The function we are interested in, φ(λ₁, …, λ_n), is precisely equal to ϕ[ν_p(x),ξ(x, x′)] for the specific choice:

$$\begin{array}{c}\hfill {\nu }_p(\mathit{x})={\nu }_p{\displaystyle \sum _i}{\lambda }_i{\mathcal{W}}_i(\mathit{x}).\end{array}$$(B.7)

We now consider ϕ[ν_p(x),ξ(x, x′)] as a function of ν₁(x) only, with the other quantities it depends on being fixed. We can then define τ(x) as the functional derivative of ϕ with respect to ν₁(x):

$$\begin{array}{c}\hfill \tau (\mathit{x})=\frac{\delta \varphi }{\delta {\nu }_1(\mathit{x})}.\end{array}$$(B.8)

Then τ(x) appears to be the generating function of all trees with (at least) one leaf at position x. This is a sub-part of rooted-trees¹⁰, and we therefore call τ(x) the rooted-Cumulant Generating Function (r-CGF). The key relevant property is then that τ(x) obeys a consistency relation, as it can be built recursively: when a line emerges from position x, it should reach a vertex of order p ≥ 1 at position x′ which is then connected to p − 1 r-CGFs at position x′. The mathematical transcription of the property is that

$$\begin{array}{cc}\hfill \tau (\mathit{x})=& {\displaystyle \int \mathrm{d}}{\mathit{x}}^{\prime }\xi (\mathit{x},{\mathit{x}}^{\prime }){\nu }_1({\mathit{x}}^{\prime })+{\displaystyle \int \mathrm{d}}{\mathit{x}}^{\prime }\xi (\mathit{x},{\mathit{x}}^{\prime })\sum _{p=2}^{\infty }\frac{{\nu }_p({\mathit{x}}^{\prime })}{(p-1)!}{\tau }^{p-1}({\mathit{x}}^{\prime }).\hfill \end{array}$$(B.9)

Once τ(x) is known, solving for ϕ then simply involves integrating Eq. (B.8). This can be done with the help of its Legendre transform with respect to ν₁(x). Defining ψ as a functional of τ(x) via

$$\begin{array}{c}\hfill \psi ={\displaystyle \int \mathrm{d}}\mathit{x}\tau (\mathit{x}){\nu }_1(\mathit{x})-\varphi ,\end{array}$$(B.10)

the relation (B.8) is satisfied when

$$\begin{array}{c}\hfill \frac{\delta \psi [\tau (\mathit{x})]}{\delta \tau (\mathit{x})}={\nu }_1(\mathit{x}).\end{array}$$(B.11)

As ξ(x, x′) is a definite positive operator, the relation (B.9) can be inverted for ν₁(x) as

$$\begin{array}{c}\hfill {\nu }_1(\mathit{x})={\displaystyle \int \mathrm{d}}{\mathit{x}}^{\prime }{\xi }^{-1}(\mathit{x},{\mathit{x}}^{\prime }){\tau }_1({\mathit{x}}^{\prime })-\sum _{p=2}^{\infty }\frac{{\nu }_p(\mathit{x})}{(p-1)!}{\tau }^{p-1}(\mathit{x}),\end{array}$$(B.12)

and formally integrated via Eq. (B.11) into

$$\begin{array}{c}\hfill \psi =\frac{1}{2}{\displaystyle \int \mathrm{d}}\mathit{x}\mathrm{d}{\mathit{x}}^{\prime }\tau (\mathit{x}){\xi }^{-1}(\mathit{x},{\mathit{x}}^{\prime })\tau ({\mathit{x}}^{\prime })-{\displaystyle \int \mathrm{d}}\mathit{x}{\displaystyle \sum _{p=2}^{\infty }}\frac{{\nu }_p(\mathit{x})}{p!}{\tau }^p(\mathit{x}),\end{array}$$

with the correct boundary conditions (it should vanish when ν₁(x) does). From Eq. (B.10), the function ϕ is then expressed as:

$$\begin{array}{c}\hfill \varphi ={\displaystyle \int \mathrm{d}}\mathit{x}{\displaystyle \sum _{p=1}^{\infty }}\frac{{\nu }_p(\mathit{x})}{p!}\tau {(\mathit{x})}^p-\frac{1}{2}{\displaystyle \int \mathrm{d}}\mathit{x}\mathrm{d}{\mathit{x}}^{\prime }\tau (\mathit{x}){\xi }^{-1}(\mathit{x},{\mathit{x}}^{\prime })\tau ({\mathit{x}}^{\prime }),\end{array}$$(B.13)

which solves the formal calculation. We can then apply this result for our specific setting (using Eq. B.7) to derive the CGF we are interested in. It is then convenient to re-express τ(x) in Eq. (B.13) in terms of the ζ function (Eq. B.3), the profiles, 𝒲_i(x), and the λ_i variables. It eventually yields the following expression:

$$\begin{array}{cc}\hfill \phi ({\lambda }_1,\cdots ,{\lambda }_n)& ={\displaystyle \sum _j}{\lambda }_j{\displaystyle \int }\mathrm{d}\mathit{x}{\mathcal{W}}_j(\mathit{x})\zeta (\tau (\mathit{x}))\hfill \\ \hfill & -\frac{1}{2}{\displaystyle \sum _j}{\lambda }_j{\displaystyle \int }\mathrm{d}\mathit{x}{\mathcal{W}}_j(\mathit{x})\tau (\mathit{x}){\zeta }^{\prime }(\tau (\mathit{x})),\hfill \end{array}$$(B.14)

and the equation for the r-CGF function reads:

$$\begin{array}{c}\hfill \tau (\mathit{x})={\displaystyle \sum _j}{\lambda }_j{\displaystyle \int }\mathrm{d}{\mathit{x}}^{\prime }{\mathcal{W}}_j({\mathit{x}}^{\prime })\xi (\mathit{x},{\mathit{x}}^{\prime }){\zeta }^{\prime }(\tau ({\mathit{x}}^{\prime })).\end{array}$$(B.15)

The expressions presented in the main text, Eqs. (10) and (11), are obtained for the specific expression of ζ corresponding to the Rayleigh-Levy flight statistics; namely, ζ(τ) = (1 + τ/2)². The generic expression for the multiple CGF of hierarchical tree models was presented in Bernardeau & Schaeffer (1999).

B.2. Mean field equation

To make the resolution of the implicit Eqs. (11) more tractable, it is possible to assume that within each cell, τ(x) can be approximated by a constant equal to τ_i. The consistency equations for the τ_is, Eq. (11) is then expressed as:

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

where ξ_ij is the average density correlation between cells i and j,

$$\begin{array}{c}\hfill {\xi }_{\mathit{ij}}={\displaystyle \int }\mathrm{d}\mathit{x}{\mathcal{W}}_i(\mathit{x})\mathrm{d}{\mathit{x}}^{\prime }{\mathcal{W}}_j({\mathit{x}}^{\prime })\xi (\mathit{x},{\mathit{x}}^{\prime }).\end{array}$$(B.17)

From Eq. (10), we then have

$$\begin{array}{c}\hfill {\phi }_n(\{{\lambda }_i\})={\displaystyle \sum _i}{\lambda }_i(\zeta ({\tau }_i)-\frac{1}{2}{\tau }_i{\zeta }^{\prime }({\tau }_i)).\end{array}$$(B.18)

This is the form we will mostly use in the main text.

Appendix C: Minimal tree model

The minimal tree model (MTM) is defined as a tree model for which only vertices with two outcoming lines exist. It is therefore associated with a vertex generating function of the form:

$$\begin{array}{c}\hfill \zeta (\tau )=1+\tau +\frac{{\nu }_2}{2}{\tau }^2.\end{array}$$(C.1)

It has been shown (in Balian & Schaeffer 1989, from the behaviour of the void probability function) that the only possible value for ν₂ is ν₂ = 1/2. This is precisely the case of the Rayleigh Levy flight model. From Eq. (C.1), this model is thus characterised by

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

In this case the stationary equations for a set of cells, Eq. (B.16) is expressed as:

$$\begin{array}{c}\hfill {\tau }_i={\displaystyle \sum _j}{\lambda }_j{\xi }_{\mathit{ij}}(1+{\tau }_j/2),\end{array}$$(C.3)

and Eq. (B.18) becomes

$$\begin{array}{c}\hfill {\phi }_n(\{{\lambda }_i\})={\displaystyle \sum _i}{\lambda }_i(1+\frac{{\tau }_i}{2}).\end{array}$$(C.4)

As the right hand side of Eq. (C.3) is linear in τ, this system can be solved by the simple inversion of a matrix. In practice it is therefore relatively easy to derive the expression of τ_i (and therefore of the cumulant generating function) for a finite number of cells.

C.1. Mean field solution

From Eqs. (C.3)-(C.4) when only one cell is considered, the one-point mean field solution turns out to be

$$\begin{array}{c}\hfill {\phi }_1({\lambda }_{\rho };{\xi }_0)=\frac{2{\lambda }_{\rho }}{2-{\xi }_0{\lambda }_{\rho }}.\end{array}$$(C.5)

It gives the expression of the CGF for a single variable ρ of mean unity and mean square ξ₀. The successive reduced cumulants of such a quantity can be easily computed by Taylor expansion.

We note that if the random variable ρ is rescaled to have a mean of $\overline{\rho }$ (instead of 1), its CGF is ${\phi }_1(\overline{\rho }\lambda ;{\xi }_0)$, which, given Eq. (C.5), is equal to $\overline{\rho }{\phi }_1(\lambda ;\overline{\rho }{\xi }_0)$.

C.2. MTM composition law and scale convolution

We generically define the cumulant generating function, φ_n({λ_i};{ξ_ij}), of a set of n cells, of variables λ_i for cells whose two-point correlations are ξ_ij. Then the following identity is satisfied by φ_n:

$$\begin{array}{c}\hfill {\phi }_n(\{{\lambda }_i\};\{{\xi }_{\mathit{ij}}+\xi \})={\phi }_1({\phi }_n(\{{\lambda }_i\};\{{\xi }_{\mathit{ij}}\});\xi ).\end{array}$$(C.6)

This important scaling, which we use continuously throughout, may seem at first view as an awkward relation, as it states that the ξ in φ_n can be shifted by successive application of the function φ₁. This result seems specific to Rayleigh-Levy flights.

We first briefly demonstrate the property. To avoid confusion in the notation, we define τ_i and ${\widehat{\tau }}_i$ as the solutions of the following systems:

$$\begin{array}{cc}\hfill {\tau }_i& ={\displaystyle \sum _j}{\lambda }_j{\xi }_{\mathit{ij}}(1+{\tau }_j/2),\hfill \end{array}$$(C.7a)

$$\begin{array}{cc}\hfill {\widehat{\tau }}_i& ={\displaystyle \sum _j}{\lambda }_j({\xi }_{\mathit{ij}}+\xi )(1+{\widehat{\tau }}_j/2).\hfill \end{array}$$(C.7b)

Then we also define μ_i and ${\widehat{\mu }}_i$ as μ_i = 1 + τ_i/2 and ${\widehat{\mu }}_i=1+{\widehat{\tau }}_i/2$ so that

$$\begin{array}{cc}\hfill {\phi }_n& \equiv {\phi }_n(\{{\lambda }_i\};\{{\xi }_{\mathit{ij}}\})={\displaystyle \sum _i}{\lambda }_i{\mu }_i,\hfill \end{array}$$(C.8a)

$$\begin{array}{cc}\hfill {\widehat{\phi }}_n& \equiv {\phi }_n(\{{\lambda }_i\};\{{\xi }_{\mathit{ij}}+\xi \})={\displaystyle \sum _i}{\lambda }_i{\widehat{\mu }}_i.\hfill \end{array}$$(C.8b)

Then ${\widehat{\mu }}_j$ is the only solution of the system

$$\begin{array}{c}\hfill {\displaystyle \sum _j}{M}_{\mathit{ij}}{\widehat{\mu }}_j=1,\end{array}$$(C.9)

where the matrix M_ij is defined as

$$\begin{array}{c}\hfill M_{\mathit{ij}}=({\delta }_{\mathit{ij}}-\frac{1}{2}{\displaystyle \sum _j}{\lambda }_j({\xi }_{\mathit{ij}}+\xi )).\end{array}$$(C.10)

Now, since

$$\begin{array}{c}\hfill {\displaystyle \sum _j}{M}_{\mathit{ij}}{\mu }_j=1-\frac{\xi }{2}{\displaystyle \sum _j}{\lambda }_j{\mu }_j,\end{array}$$(C.11)

this implies that ${\mu }_i/\left(1-\frac{\xi }{2}{\displaystyle \sum _j}{\lambda }_j{\mu }_j\right)$ is the solution of Eq. (C.11), so that by identification we have:

$$\begin{array}{c}\hfill {\widehat{\mu }}_i=\frac{{\mu }_i}{1-\frac{\xi }{2}{\sum }_j{\lambda }_j{\mu }_j}.\end{array}$$(C.12)

Hence,

$$\begin{array}{c}\hfill {\widehat{\phi }}_n={\displaystyle \sum _i}{\lambda }_i\frac{{\mu }_i}{1-\frac{{\displaystyle \xi }}{2}{\sum }_j{\lambda }_j{\mu }_j}=\frac{{\phi }_n}{1-\frac{\xi }{2}{\phi }_n}={\phi }_1({\phi }_n;\xi ),\end{array}$$(C.13)

which establishes the relation.

Equation (C.6) reflects some interesting physical properties it is related to: scale composition when cell correlations are built as a two-steps procedure. Indeed, we consider a set of random Rayleigh-Levy flights experiments, in which instead of having of fixed number of points in each sample, the sample density is itself drawn from the one-point PDF derived from an MTM process. We denote ρ_s the sample density (whose average is set to unity). Its CGF is then φ₁(λ_s; ξ_s), where ξ_s is the variance of the sample density.

We then note that for each realisation the cell densities, ρ_i scale like 1/ρ_s (by definition ρ_i is defined with respect to the mean of the survey) and that its correlation functions ξ_ij scale like 1/ρ_s (this is a consequence of the expression of the scaling of the two-point function in the Rayleigh Levy model). To be more precise, we define ${\overline{\xi }}_{\mathit{ij}}$ as the cell correlations when the sample density is equal to unity. We then have

$$\begin{array}{c}\hfill {\xi }_{\mathit{ij}}={\overline{\xi }}_{\mathit{ij}}/{\rho }_s.\end{array}$$(C.14)

We then define the densities ${\widehat{\rho }}_i$ which represent the ‘true’ density in the sample (that is when the sample density is taken into account) as:

$$\begin{array}{c}\hfill {\widehat{\rho }}_i={\rho }_s{\rho }_i,\end{array}$$(C.15)

and aim at building the joint PDF of ${\widehat{\rho }}_i$. This PDF is formally given by

$$\begin{array}{cc}\hfill P& (\{{\widehat{\rho }}_i\})={\displaystyle \int \frac{\mathrm{d}{\rho }_s}{{\rho }_s^n}}P({\rho }_s)P(\{{\rho }_i/{\rho }_s\}),\hfill \\ \hfill & ={\displaystyle \int \frac{\mathrm{d}{\rho }_s}{{\rho }_s^n}}P({\rho }_s)\int {\mathrm{\Pi }}_i\frac{\mathrm{d}{\lambda }_i}{2\pi i}exp(-{\displaystyle \sum _i}{\lambda }_i{\widehat{\rho }}_i/{\rho }_s+{\phi }_n(\{{\lambda }_i\};\{{\xi }_{\mathit{ij}}\})),\hfill \end{array}$$(C.16)

where the dependence in ρ_s is also present in the expression of ξ_ij. Then, making the change of variable

$$\begin{array}{c}\hfill {\widehat{\lambda }}_i={\lambda }_i/{\rho }_s,\end{array}$$(C.17)

and noting that:

$$\begin{array}{c}\hfill {\phi }_n(\{{\lambda }_i\};\{{\xi }_{\mathit{ij}}\})={\rho }_s{\phi }_n(\{{\widehat{\lambda }}_i\};\{{\overline{\xi }}_{\mathit{ij}}\}),\end{array}$$(C.18)

we have

$$\begin{array}{cc}\hfill P(\{{\widehat{\rho }}_i\})& ={\displaystyle \int }\mathrm{d}{\rho }_s{\displaystyle \int \frac{\mathrm{d}\mu }{2\pi i}}\int {\mathrm{\Pi }}_i\frac{\mathrm{d}{\lambda }_i}{2\pi i}\hfill \\ \hfill & \times exp(-{\displaystyle \sum _i}{\widehat{\lambda }}_i{\widehat{\rho }}_i-\mu {\rho }_s+{\rho }_s{\phi }_n(\{{\widehat{\lambda }}_i\};\{{\overline{\xi }}_{\mathit{ij}}\})+{\phi }_1(\mu ,{\xi }_s)).\hfill \end{array}$$(C.19)

Now the integral over ρ_s leads to a Dirac delta function in $\mu -{\phi }_n(\{{\widehat{\lambda }}_i\};\{{\overline{\xi }}_{\mathit{ij}}\})$ leading to

$$\begin{array}{c}\hfill P(\{{\widehat{\rho }}_i\})=\int {\mathrm{\Pi }}_i\frac{\mathrm{d}{\lambda }_i}{2\pi i}exp(-{\displaystyle \sum _i}{\widehat{\lambda }}_i{\widehat{\rho }}_i+{\phi }_1({\phi }_n(\{{\widehat{\lambda }}_i\};\{{\overline{\xi }}_{\mathit{ij}}\}),{\xi }_s)).\end{array}$$(C.20)

This expression means formally that the CGF of ${\widehat{\rho }}_i$ is this two-step construction given by ${\phi }_1({\phi }_n(\{{\widehat{\lambda }}_i\};\{{\overline{\xi }}_{\mathit{ij}}\}),{\xi }_s)$. The relation (C.6) ensures that it is also given by ${\phi }_n(\{{\widehat{\lambda }}_i\};\{{\overline{\xi }}_{\mathit{ij}}+{\xi }_s\})$, which states that survey density fluctuations can, in this model, be taken into account via a simple shift in the cell correlation amplitudes.

A useful practical consequences of property (C.6) is that we could set any peculiar element of ξ_ij to zero. For instance, the CGF of two identical cells of density variance $\overline{\xi }$ and of cross-correlation ξ₁₂ is given by:

$$\begin{array}{c}\hfill {\phi }_2({\lambda }_1,{\lambda }_2;\overline{\xi },{\xi }_{12})={\phi }_1({\phi }_1({\lambda }_1;\overline{\xi }-{\xi }_{12})+{\phi }_1({\lambda }_2;\overline{\xi }-{\xi }_{12});{\xi }_{12}).\end{array}$$

It is possible to extend such a construction for a larger number of cells¹¹, but for specific configurations only. Note finally that when one needs to construct the density CGF for a large number of cells it is convenient to set $\overline{\xi }=0$ for all cells as it makes the system in τ_i sparser.

C.3. Computing PDF in the minimal tree model

The basic quantity one wishes to have access to is the PDF, defined as the inverse Laplace transform of the cumulant generating function φ₁(λ_ρ), that is,

$$\begin{array}{c}\hfill P(\rho )={\displaystyle {{\displaystyle \int }}_{-i\infty }^{+i\infty }}\frac{\mathrm{d}{\lambda }_{\rho }}{2\pi i}exp(-{\lambda }_{\rho }\rho +{\phi }_1({\lambda }_{\rho })),\end{array}$$(C.21)

where the integral runs formally along the imaginary axis, but can be moved in the complex plane as long as no poles or singularities are encountered along the path.

The following relation will be exploited throughout this paper

$$\begin{array}{c}\hfill {\displaystyle \int \frac{\mathrm{d}{\lambda }_{\rho }}{2\pi i}}exp\left(-{\lambda }_{\rho }\rho +\frac{a_c}{{\lambda }_c-{\lambda }_{\rho }}\right)={\delta }_d(\rho )+a_c{e}^{-{\lambda }_c\rho }{}_0{\stackrel{\sim }{F}}_1(2;a_c\rho ),\end{array}$$(C.22)

where ${\delta }_{\mathrm{D}}$ is the 1D Dirac distribution and ${}_0{\stackrel{\sim }{F}}_1(a;z)$ is the regularised confluent hypergeometric function  ₀F₁(a; z)/Γ(a). This formula can be applied directly to φ₁(λ_ρ; ξ₀) to derive P(ρ). It can also be applied after successive derivatives with respect to a_c from which on can establish that

$$\begin{array}{cc}\hfill {\displaystyle \int \frac{\mathrm{d}{\lambda }_{\rho }}{2\pi i}}{\left(\frac{1}{{\lambda }_c-{\lambda }_{\rho }}\right)}^n& exp(-{\lambda }_{\rho }\rho +\frac{a_c}{{\lambda }_c-{\lambda }_{\rho }})=\hfill \\ \hfill & e^{-{\lambda }_c\rho }{\left(\frac{\mathrm{d}}{\mathrm{d}{a}_c}\right)}^n\left\{a_c{}_0{\stackrel{\sim }{F}}_1(2;a_c\rho )\right\},\hfill \end{array}$$(C.23)

for n ≥ 1. Then, the application of Eq. (C.24) to Eq. (C.5) leads to:

$$\begin{array}{c}\hfill P(\rho )=e^{-2/{\xi }_0}{\delta }_d(\rho )+\frac{4}{{\xi }_0^2}{e}^{-\frac{2(\rho +1)}{{\xi }_0}}{}_0{\stackrel{\sim }{F}}_1(2;\frac{4\rho }{{\xi }_0^2}),\end{array}$$(C.24)

which yields the PDF in the mean field approximation. One can see that it involves the sum of two terms, a Dirac term term for ρ = 0, and a continuous contribution. This highlights one of the key feature of this model, which is that the probability of having empty regions of finite size remains finite even in the continuous limit.

C.4. Beyond the mean field approximation

We explore the possibility of solving the consistency relations beyond the mean field approximation, which provides means to explore the validity of this approximation. The calculations will be limited to the 1D case and to a Gaussian filter and its derivatives. We therefore assume that

$$\begin{array}{c}\hfill {\mathcal{W}}_0(x)=\frac{1}{\sqrt{2\pi }}exp(-\frac{x^2}{2}).\end{array}$$(C.25)

The idea is then to expand τ(x) in, for instance, Eq. (11) on its natural ortho-normal basis, namely the basis of the Hermite polynomials. More precisely, τ(x)𝒲₀(x) being bounded at large distance, it is possible to expand it as

$$\begin{array}{c}\hfill \tau (x)={\displaystyle \sum _{n=0}^{\infty }}{\tau }_n{h}_n(x),\mathrm{with}{h}_n(x)=\frac{1}{2^{n/2}\sqrt{n!}}{H}_n\left(\frac{x}{\sqrt{2}}\right),\end{array}$$(C.26)

noting that

$$\begin{array}{c}\hfill {\displaystyle \int \mathrm{d}}x{h}_n(x)h_{n^{\prime }}(x){\mathcal{W}}_0(x)={\delta }_{n{n}^{\prime }}.\end{array}$$(C.27)

We are specifically interested in the joint CGF, φ({λ₀, λ₁, λ₂}) involving fields dual to the density, the first and second derivatives. The window functions for the latter can be expressed in terms of the Hermite polynomials as

$$\begin{array}{c}\hfill {\mathcal{W}}_q(x)={\left(\frac{\mathrm{d}}{\mathrm{d}x}\right)}^q{\mathcal{W}}_0(x)=\sqrt{q!}{h}_q(x){\mathcal{W}}_0(x),\end{array}$$(C.28)

for q = 1 and q = 2, and we have

$$\begin{array}{c}\hfill {\displaystyle \int \mathrm{d}}x{\mathcal{W}}_q(x)h_n(x)=\sqrt{q!}{\delta }_{\mathit{qn}}.\end{array}$$(C.29)

As a result, the consistency relation in Eq. (11) for τ(x) is now expressed as:

$$\begin{array}{c}\hfill \tau (x)=\int \mathrm{d}{x}^{\prime }\xi (x,x^{\prime })(1+\frac{\tau (x^{\prime })}{2}){\mathcal{W}}_0(x^{\prime }){\displaystyle \sum _{q=0}^2}{\lambda }_q\sqrt{q!}{h}_q(x^{\prime }),\end{array}$$(C.30)

which can then be transformed into a system in τ_p after integration with a weight of h_p(x)C₀(x):

$$\begin{array}{c}\hfill {\tau }_p={\displaystyle \sum _{q=0}^2}{\lambda }_q[{\mathrm{\Xi }}_{p0}^{(q)}+{\displaystyle \sum _{p^{\prime }=0}^{\infty }}\frac{{\tau }_{p^{\prime }}}{2}{\mathrm{\Xi }}_{p{p}^{\prime }}^{(q)}],\end{array}$$(C.31)

with

$$\begin{array}{c}\hfill {\mathrm{\Xi }}_{p{p}^{\prime }}^{(q)}={\displaystyle \int \mathrm{d}}x\mathrm{d}{x}^{\prime }{\mathcal{W}}_0(x)h_p(x)\xi (x,x^{\prime })h_{p^{\prime }}(x^{\prime }){\mathcal{W}}_q(x^{\prime }).\end{array}$$(C.32)

The expression of CGF for the density and its gradients derives from Eq. (10), and is simply given by

$$\begin{array}{c}\hfill \phi (\{{\lambda }_q\})={\lambda }_0+{\displaystyle \sum _{q=0}^2}{\lambda }_q\frac{{\tau }_q}{2}\sqrt{q!},\end{array}$$(C.33)

taking advantage of the orthogonality relations between 𝒲_q(x) and h_p(x). We note that those quantities depend on the actual shape and amplitude of the correlation function. In the following, we simply assume that ξ(x, x′) ∼ |x − x′|^−1/2 corresponding to α = 1/2 for a 1D Rayleigh-Levy flight. It implies that ${\Xi }_{p{p}^{\prime }}^{(q)}/{\xi }_0$ are all fixed quantities and that the result can therefore be expressed in terms of ξ₀ only.

We present hereafter the result of such derivations for the density, when the Hermite expansion in τ is truncated at increasing order, from 0 (where the r-CGF is assumed to be constant) to 4 (in practice we have been able to derive the CGF up to 10th order as illustrated in the figures below). We thus have

$$\begin{array}{cc}\hfill {\phi }_0(\lambda )& =\frac{2\lambda }{2-\lambda {\xi }_0},\hfill \end{array}$$(C.34a)

$$\begin{array}{cc}\hfill {\phi }_2(\lambda )& =-\frac{\lambda (5\lambda {\xi }_0-64)}{2{\lambda }^2{\xi }_0^2-37\lambda {\xi }_0+64},\hfill \end{array}$$(C.34b)

$$\begin{array}{cc}\hfill {\phi }_4(\lambda )& =-\frac{3\lambda (75{\lambda }^2{\xi }_0^2-8240\lambda {\xi }_0+65536)}{80{\lambda }^3{\xi }_0^3-10849{\lambda }^2{\xi }_0^2+123024\lambda {\xi }_0-196608},\hfill \end{array}$$(C.34c)

for subsequent truncation orders. The first expression reproduces the one-cell mean field approximation. The others correspond to corrections to it.

These expressions exhibit a number of worthwhile properties. They predict values of reduced cumulants that slightly evolve with the order of the truncation. This is illustrated on Figure C.1 for the skewness. It is actually possible to compute its exact expression for the Rayleigh-Levy flight model with a given slope for the two-point function (or, equivalently, the index of the power spectrum) and a given filter shape (here, a Gaussian filter). Figure C.1 shows that the prediction from the mean field expansion rapidly converges to the exact value. The rapid convergence can also be observed for higher order cumulants. This suggests that the high density tail of the PDF is well captured by mean field approach.

On the other hand the low density part turns out to be much more difficult to capture with a Gaussian filter. As illustrated in Figure 2, the density fields exhibit genuine large empty regions. This leads to a non-zero void probability density function for filters that have a compact support. However, for filters with extended radial tails, the probability of finding the filtered density to be exactly zero can only vanish, leading however to an excess of probability at low densities. This is illustrated in Figure 3. Unfortunately the mean field extensions as described here always lead to finite VPDF to all finite orders. This effect is illustrated in Fig. C.2, which also shows that the convergence of the PDF in the low density regions is slow. This is also illustrated at the level of the density PDF in Fig. 3. To conclude, i) the mean field solution is very efficient at predicting the high density regions but fails in the low density regions; ii) solutions beyond the mean field, Eqs. (C.4), can account for the behaviour of the PDFs in the low density regions, but the convergence is slow. Extension beyond the mean field in higher dimensions will be the topic of future work.

Appendix D: Extrema and Euler densities

D.1. General formalism

The computation of extrema densities is relies on the knowledge of the joint PDF of the local density, its first and second order derivatives (Bardeen et al. 1986). The latter will be derived from the CGF of those quantities.

The number of relevant variables depend on the spatial dimension. To fix the notation we define ξ₀, ξ₁, and ξ₂ as the variance of respectively the local density, the one-component density gradient and second derivatives:

$$\begin{array}{c}\hfill {\xi }_0=\langle {\rho }^2\rangle -1;{\xi }_1=\langle {\rho }_{,x}^2\rangle ;{\xi }_2=\langle {\rho }_{,xx}^2\rangle ,\end{array}$$(D.1)

assuming 〈ρ〉 = 1. We further note the consistency relations, which can be derived by integration via:

$$\begin{array}{c}\hfill \langle \rho {\rho }_{,xx}\rangle =-{\xi }_1;\langle {\rho }_{,xx}{\rho }_{,yy}\rangle =\frac{1}{3}{\xi }_2;\langle {\rho }_{,xy}^2\rangle =\frac{1}{3}{\xi }_2.\end{array}$$(D.2)

We then denote λ_ρ, λ_i, and λ_ij the conjugate variables to respectively ρ, ρ_,i and ρ_,ij (keeping in mind that ρ_,ij and ρ_,ji are identical).

In general the number density of extrema is given by Eq. (19). For instance, maxima are obtained when the integral is restricted to the regions where all eigenvalues are negative; minima when all eigenvalues are positive and in 2D, saddle points are obtained when the integral is restricted to the regions where the sign of the two eigenvalues is different. We finally note that Euler number densities are obtained with χ^extr.[ρ_,ij] = 1. Its advantage is that it preserves the analytical structure of the operator making in general its derivation easier. Specifically, the formulae we obtain for the 1D and 2D cases are:

$$\begin{array}{cc}\hfill n_{1\mathrm{D}}^{\mathrm{Euler}}(\rho )& ={\displaystyle \int }\mathrm{d}{\rho }_{,xx}{\rho }_{,xx}P(\rho ,{\rho }_{,x}=0,{\rho }_{,xx}),\hfill \\ \hfill n_{2\mathrm{D}}^{\mathrm{Euler}}(\rho )& ={\displaystyle \int }\mathrm{d}{\rho }_{,xx},\mathrm{d}{\rho }_{,yy},\mathrm{d}{\rho }_{,xy}({\rho }_{,xx}{\rho }_{,yy}-{\rho }_{,xy}^2)P(\rho ,{\rho }_{,i}=0,{\rho }_{,ij}),\hfill \end{array}$$(D.3)

and a similar (if more involved) expression for the 3D case. For the Euler number density, taking advantage of its analyticity, it is possible to simply re-express it in terms of the CGF through inverse Laplace transforms. More precisely we have in 1D:

$$\begin{array}{cc}\hfill n_{1\mathrm{D}}^{\mathrm{Euler}}(\rho )& ={\displaystyle \int \frac{\mathrm{d}{\lambda }_{\rho }}{2\pi i}}\frac{\mathrm{d}{\lambda }_x}{2\pi i}{\frac{\partial }{\partial {\lambda }_{\mathit{xx}}}}_{{|}_{{\lambda }_{\mathit{xx}}=0}}exp(-{\lambda }_{\rho }\rho +\phi ({\lambda }_{\rho },{\lambda }_x,{\lambda }_{\mathit{xx}})).\hfill \end{array}$$

It shows that the Euler number density, unlike the extrema density in general, depends only on limited information from the joint CGF. For the 2D case, the calculations are similar but a bit more involved, as:

$$\begin{array}{cc}\hfill n_{2\mathrm{D}}^{\mathrm{Euler}}(\rho )& ={\displaystyle \int \frac{\mathrm{d}{\lambda }_{\rho }}{2\pi i}}\frac{\mathrm{d}{\lambda }_x}{2\pi i}\frac{\mathrm{d}{\lambda }_y}{2\pi i}{(\frac{{\partial }^2}{\partial {\lambda }_{\mathit{xx}}\partial {\lambda }_{\mathit{yy}}}-\frac{{\partial }^2}{\partial {\lambda }_{\mathit{xy}}^2})}_{{|}_{{\lambda }_{\mathit{ij}}=0}}\hfill \\ \hfill & \times exp(-{\lambda }_{\rho }\rho +\phi ({\lambda }_{\rho },{\lambda }_i,{\lambda }_{\mathit{ij}})),\hfill \end{array}$$(D.5)

with a similar construction for the 3D case. Explicit derivation of the Euler number densities for the Rayleigh-Levy flight model are presented hereafter. We start with general consideration regarding the symmetry properties of both the CGF and the joint PDFs.

D.2. Spatial homogeneity & isotropy

In the cosmological context, density fields are statistically homogeneous and isotropic. For the hierarchical tree models, it is ensured by the fact that the two-point correlation function ξ₂(x, x′) depends on |x − x′| only.

A number of properties follows. The first immediate one, for the 1D case, is that the expectation value of any gradient should vanish. That implies in particular that the expectation value of ${\rho }_{,x}^2+\rho {\rho }_{,xx}$ vanishes, which implies Eq. (D.2). More generally a consequence of this invariance is that, for any integer p and q,

$$\begin{array}{c}\hfill p{\langle {\rho }^{p-1}{\rho }_{,x}^{q+1}\rangle }_c+q{\langle {\rho }^p{\rho }_{,x}^{q-1}{\rho }_{,xx}\rangle }_c=0,\end{array}$$(D.6)

which in terms, after summing over p and q, implies

$$\begin{array}{c}\hfill {\lambda }_{\rho }\frac{\partial }{\partial {\lambda }_x}\phi (\lambda ,{\lambda }_x)+{\lambda }_x{\frac{\partial }{\partial {\lambda }_{\mathit{xx}}}|}_{{\lambda }_{\mathit{xx}}=0}\phi (\lambda ,{\lambda }_x,{\lambda }_{\mathit{xx}})=0.\end{array}$$(D.7)

An alternative formulation at the level of the PDF reads

$$\begin{array}{c}\hfill {\displaystyle \int \mathrm{d}}{\rho }_{,xx}\{{\rho }_{,x}\frac{\partial P(\rho ,{\rho }_{,x},{\rho }_{,xx})}{\partial \rho }+{\rho }_{,xx}\frac{\partial P(\rho ,{\rho }_{,x},{\rho }_{,xx})}{\partial {\rho }_{,x}}\}=0.\end{array}$$(D.8)

Integrating the second term of this equation over ρ_,x up to ρ_,x = 0 leads to the expression (D.3) of $n_{1\mathrm{D}}^{\mathrm{Euler}}(\rho )$ so that the latter can eventually be written in terms of the joint probability of the density and its gradient as

$$\begin{array}{c}\hfill n_{1\mathrm{D}}^{\mathrm{Euler}}(\rho )=\frac{\partial }{\partial \rho }{\displaystyle {{\displaystyle \int }}_{-\infty }^0}{\rho }_{,x}\mathrm{d}{\rho }_{,x}P(\rho ,{\rho }_{,x}).\end{array}$$(D.9)

It shows that the Euler number density does not depend actually on the way ρ_,xx is correlated to the density and its gradient. This is not so for the extrema counts. This points to an interesting feature of Euler number density: it is significantly more insensitive to small scale fluctuations compared to other observables.

Unfortunately this simplification does not directly extend to higher dimensions. There are however a number of consistency relations that derive from statistical invariance under translation, parity change and rotation. They are

1. Translation invariance¹²: as in 1D, for any direction, i, we expect

   $$\begin{array}{c}\hfill {\lambda }_{\rho }\frac{\partial }{\partial {\lambda }_i}\phi (\lambda ,\{{\lambda }_i\})+{\displaystyle \sum _j}{\lambda }_j{\frac{\partial }{\partial {\lambda }_{\mathit{ij}}}}_{{|}_{{\lambda }_{\mathit{ij}}=0}}\phi (\lambda ,\{{\lambda }_i\},\{{\lambda }_{\mathit{ij}}\})=0.\end{array}$$(D.10)

2. Parity invariance: For dimensions above 2D, there are other combinations that vanish such that the expectation values of T_,xU_,y − T_,yU_,x (as can be shown by successive integrations by parts). For the 2D case, this is the transcription of parity invariance: the expectation of any rotational is expected to vanish. It implies that

   $$\begin{array}{c}\hfill {\lambda }_y{(\frac{{\partial }^2\phi }{\partial {\lambda }_{\mathit{xx}}\partial {\lambda }_{\mathit{yy}}}-\frac{{\partial }^2\phi }{\partial {\lambda }_{\mathit{xy}}^2})}_{{|}_{{\lambda }_{\mathit{ij}}=0}}={\lambda }_{\rho }{(\frac{{\partial }^2\phi }{\partial {\lambda }_x\partial {\lambda }_{\mathit{xy}}}-\frac{{\partial }^2\phi }{\partial {\lambda }_y\partial {\lambda }_{\mathit{xx}}})}_{{|}_{{\lambda }_{\mathit{ij}}=0}}\end{array}$$

   and a similar equation after substitution y ↔ x. We note that the left-hand-side of this equation is the operator for the computation of the Hessian in 2D. It can be reduced to first order derivatives in λ_ij. Contributions from joint cumulants between the density and its first and second order derivatives still however contribute to the Euler number density.

3. Rotation invariance: for a sake a completeness one can also derive the consequences of rotational invariance. We do it here for the 2D case: forcing λ_xρ_x + λ_yρ_y + λ_xxρ_xx + λ_xyρ_xy + λ_yyρ_yy to be rotation invariant, a rotation of the coordinate system by infinitesimal angle θ leads to

   $$\begin{array}{cc}\hfill \delta {\lambda }_x& =-\theta {\lambda }_y,\delta {\lambda }_y=\theta {\lambda }_x,\hfill \end{array}$$(D.10a)

   $$\begin{array}{cc}\hfill \delta {\lambda }_{\mathit{xx}}& =-\theta {\lambda }_{\mathit{xy}},\delta {\lambda }_{\mathit{yy}}=\theta {\lambda }_{\mathit{xy}},\hfill \end{array}$$(D.10b)

   $$\begin{array}{cc}\hfill \delta {\lambda }_{\mathit{xy}}& =-2\theta ({\lambda }_{\mathit{xx}}-{\lambda }_{\mathit{xx}}),\hfill \end{array}$$(D.10c)

   so that

   $$\begin{array}{c}\hfill \mathcal{R}[\phi ({\lambda }_{\rho },{\lambda }_i,{\lambda }_{\mathit{ij}})]=0,\end{array}$$(D.11)

   with

   $$\begin{array}{cc}\hfill {\mathcal{R}}_{\theta }& =\frac{{\lambda }_x\partial }{\partial {\lambda }_y}-\frac{{\lambda }_y\partial }{\partial {\lambda }_x}+(\frac{{\lambda }_{\mathit{xy}}\partial }{\partial {\lambda }_{\mathit{yy}}}-\frac{{\lambda }_{\mathit{xy}}\partial }{\partial {\lambda }_{\mathit{xx}}})-2({\lambda }_{\mathit{yy}}-{\lambda }_{\mathit{xx}})\frac{{\lambda }_y\partial }{\partial {\lambda }_{\mathit{xy}}}.\hfill \end{array}$$

We now turn by exploring scaling properties of the cumulants. Hierarchical models are by nature such that

$$\begin{array}{c}\hfill {\langle {\rho }^n\rangle }_c\sim {\xi }_0^{n-1}.\end{array}$$(D.12)

In general, the cumulants of the product of powers of the local density, first and second order derivatives behave like ${\xi }_0^s{\xi }_1^t{\xi }_2^u$ where s + t + u = n − 1, if the total power is n. We note that for dimensional reasons, if only the local density and first derivative is involved then

$$\begin{array}{c}\hfill {\langle {\rho }^p{\rho }_x^q\rangle }_c\sim {\xi }_0^{p+q/2-1}{\xi }_1^{q/2}.\end{array}$$(D.13)

It is not possible to give rules when both density, first and second order derivatives are present. We can however note that

$$\begin{array}{c}\hfill {\langle {\rho }^p{\rho }_{,xx}\rangle }_c=-p{\langle {\rho }^{p-1}{\rho }_{,x}^2\rangle }_c\sim {\xi }_0^{p-1}{\xi }_1.\end{array}$$(D.14)

Finally we conjecture¹³ that the following scaling is expected,

$$\begin{array}{c}\hfill {\langle {\rho }^p{\rho }_{,i}^q{\rho }_{,ij}\rangle }_c\sim {\xi }_0^{p+q/2-1}{\xi }_1^{1+q/2}.\end{array}$$(D.15)

It is indeed verified in the expressions derived below.

D.3. Derivation of join CGF in mean field limit

We now turn to writing the joint CGF for a given set of cells. The calculation is based on the mean field solution for a finite number of cells assumed to be infinitely close to one another; the derivatives are then obtained via finite differences. More precisely if two cells are at distance d, we then assume that the average correlation between the 2 cells, ξ(d), can be expanded in powers of d. Its expression can be derived using the correlation between the density and its gradient (again via integration by part)

$$\begin{array}{c}\hfill \langle \rho {\rho }_{,xx}\rangle =-{\xi }_1,\langle \rho {\rho }_{,xxxx}\rangle ={\xi }_2,\end{array}$$(D.16)

so that

$$\begin{array}{c}\hfill \xi (d)\approx {\xi }_0-\frac{1}{2}{\xi }_1{d}^2+\frac{1}{24}{\xi }_2{d}^4,\end{array}$$(D.17)

where the coefficients ξ_i are precisely those defined in Eq. (D.1).

The density gradient is then represented by a finite difference (ρ_i − ρ_j)/d is the two cells are at a distance, d. The conjugate variable associated with the density gradient, λ_x, would then contribute to both λ_i and λ_j with a weight, respectively, of 1/d and −1/d. In general we then build the multi-values CGF as a function of such quantities. To make the construction more explicit, we consider the 2D case. To do so we need at least 6 cells taken from a 3x3 grid. Each λ_α where α stands for ρ or its derivatives contribute through the following pattern

$$\begin{array}{cc}\hfill L_{\rho }& =\left(\begin{array}{ccc}0& 0& 0\\ 0& 1& 0\\ 0& 0& 0\end{array}\right),L_x=\frac{1}{2d}\left(\begin{array}{ccc}0& 0& 0\\ -1& 0& 1\\ 0& 0& 0\end{array}\right),\hfill \end{array}$$(D.18a)

$$\begin{array}{cc}\hfill L_{\mathit{xx}}& =\frac{1}{d^2}\left(\begin{array}{ccc}0& 0& 0\\ 1& -2& 1\\ 0& 0& 0\end{array}\right),L_{\mathit{xy}}=\frac{1}{d^2}\left(\begin{array}{ccc}0& -1& 1\\ 0& 1& -1\\ 0& 0& 0\end{array}\right),\hfill \end{array}$$(D.18b)

with a similar definition for L_y and L_yy. Defining the 3x3 matrix

$$\begin{array}{c}\hfill L_{\lambda }={\lambda }_{\rho }{L}_{\rho }+{\lambda }_x{L}_x+{\lambda }_y{L}_y+{\lambda }_{\mathit{xy}}{L}_{\mathit{xy}}+{\lambda }_{\mathit{xx}}{L}_{\mathit{xx}}+{\lambda }_{\mathit{yy}}{L}_{\mathit{yy}},\end{array}$$(D.19)

the stationary equations become

$$\begin{array}{c}\hfill {\tau }_i={\displaystyle \sum _{j=1}^9}{\xi }_{\mathit{ij}}(d){(L_{\lambda })}_j(1+{\tau }_j/2),\end{array}$$(D.20)

where (L_λ)_j is the jth element of the matrix L_λ. The CGF is then expressed as:

$$\begin{array}{c}\hfill \phi ({\lambda }_{\rho },{\lambda }_i,{\lambda }_{\mathit{ij}})={\displaystyle \sum _{i=1}^9}{(L_{\lambda })}_i(1+{\tau }_i/2).\end{array}$$(D.21)

The result is a function of λ_α and of d, ξ₀, ξ₁ and ξ₂ through the ξ_ij dependence. The set of nine Eqs. (D.23) can be solved via MATHEMATICA, yielding a rather complex expression.

D.4. 1D Euler and extrema counts

In 1D, in the low d limit and for ξ₀ = 0, the generating 1D function, φ_1D(λ_ρ, λ_x, λ_xx; 0, ξ₁, ξ₂), takes the remarkably simple form of:

$$\begin{array}{c}\hfill -\frac{2({\xi }_1{\lambda }_x^2+{\lambda }_{\rho }(2-2{\xi }_1{\lambda }_{\mathit{xx}})+{\xi }_2{\lambda }_{\mathit{xx}}^2(1-{\xi }_1{\lambda }_{\mathit{xx}}))}{({\xi }_1{\lambda }_{\mathit{xx}}-1){({\xi }_1{\lambda }_{\mathit{xx}}+2)}^2}.\end{array}$$(D.22)

Taking now advantage of the shifting relation (C.6), we can write, for a general ξ₀,

$$\begin{array}{cc}& {\phi }_{1\mathrm{D}}({\lambda }_{\rho },{\lambda }_{\mathrm{x}},{\lambda }_{\mathrm{xx}};{\xi }_0,{\xi }_1,{\xi }_2)=\hfill \\ \hfill & {\phi }_1[-\frac{2({\xi }_1{\lambda }_x^2+{\lambda }_{\rho }(2-2{\xi }_1{\lambda }_{\mathit{xx}})+{\xi }_2{\lambda }_{\mathit{xx}}^2(1-{\xi }_1{\lambda }_{\mathit{xx}}))}{({\xi }_1{\lambda }_{\mathit{xx}}-1){({\xi }_1{\lambda }_{\mathit{xx}}+2)}^2};{\xi }_0].\hfill \end{array}$$(D.23)

Equation (D.26) displays advantageous properties to complete the computation of number density of extrema. We indeed have:

$$\begin{array}{cc}\hfill {\phi }_{1\mathrm{D}}({\lambda }_{\rho },{\lambda }_x,{\lambda }_{\mathit{xx}};{\xi }_0,{\xi }_1,{\xi }_2)& =\hfill \\ \hfill & {\phi }_{1\mathrm{D}}({\lambda }_{\rho }+\frac{{\xi }_1{\lambda }_x^2}{2(1-{\xi }_1{\lambda }_{\mathit{xx}})},0,{\lambda }_{\mathit{xx}};{\xi }_0,{\xi }_1,{\xi }_2),\hfill \end{array}$$(D.24)

which leads naturally to the change of variable

$$\begin{array}{c}\hfill {\lambda }_e={\lambda }_{\rho }+\frac{{\xi }_1{\lambda }_x^2}{2(1-{\xi }_1{\lambda }_{\mathit{xx}})},\end{array}$$(D.25)

when computing the inverse Laplace transformation over λ_x, necessary for the joint PDF, Eq. (C.22). We are left with

$$\begin{array}{cc}& P(\rho ,{\rho }_x=0,{\rho }_{,xx})={\displaystyle \int \frac{\mathrm{d}{\lambda }_e}{2\pi i}}\frac{\mathrm{d}{\lambda }_{\mathit{xx}}}{2\pi i}\frac{{(1-{\xi }_1{\lambda }_{\mathit{xx}})}^{1/2}}{\sqrt{2\pi {\xi }_1\rho }}\hfill \\ \hfill & \times exp(-{\lambda }_e\rho -{\lambda }_{\mathit{xx}}{\rho }_{\mathit{xx}}+{\phi }_{1\mathrm{D}}(\{{\lambda }_e,0,{\lambda }_{\mathit{xx}}\},{\xi }_0,{\xi }_1,{\xi }_2)).\hfill \end{array}$$(D.26)

We can further note that the dependence in λ_e in φ_1D({λ_e,0,λ_xx};ξ₀,ξ₁,ξ₂) now takes a very simple form

$$\begin{array}{c}\hfill {\phi }_{1\mathrm{D}}(\{{\lambda }_e,0,{\lambda }_{\mathit{xx}}\};{\xi }_0,{\xi }_1,{\xi }_2)=-\frac{2}{{\xi }_0}+\frac{a_c({\lambda }_{\mathit{xx}})}{{\lambda }_c({\lambda }_{\mathit{xx}})-{\lambda }_e},\end{array}$$(D.27)

with

$$\begin{array}{cc}\hfill a_c({\lambda }_{\mathit{xx}})& =\frac{{({\xi }_1{\lambda }_{\mathit{xx}}+2)}^2}{{\xi }_0^2},\hfill \end{array}$$(D.28)

$$\begin{array}{cc}\hfill {\lambda }_c({\lambda }_{\mathit{xx}})& =\frac{1}{2{\xi }_0}{({\xi }_1{\lambda }_{\mathit{xx}}+2)}^2-\frac{{\xi }_2}{2}{\lambda }_{\mathit{xx}}^2.\hfill \end{array}$$(D.29)

The integral over λ_e in Eq. (D.29) can then be done with the use of Eq. (C.24). It leads to the following expression (for the continuous, non singular, component)

$$\begin{array}{c}\hfill P(\rho ,{\rho }_x=0,{\rho }_{\mathit{xx}})={\displaystyle \int \frac{\mathrm{d}{\lambda }_{\mathit{xx}}}{2\pi i}}{e}^{-{\lambda }_{\mathit{xx}}{\rho }_{,xx}}{\chi }_{1\mathrm{D}}(\rho ;{\lambda }_{\mathit{xx}}),\end{array}$$(D.30)

where

$$\begin{array}{cc}& {\chi }_{1\mathrm{D}}(\rho ;{\lambda }_{\mathit{xx}})=\frac{2\sqrt{2}}{{\xi }_0^2\sqrt{\pi {\xi }_1\rho }}{(1-{\xi }_1{\lambda }_{\mathit{xx}})}^{1/2}{(1+{\xi }_1{\lambda }_{\mathit{xx}}/2)}^2\times \hfill \\ \hfill & exp(-\frac{2}{{\xi }_0}+\frac{\rho {\xi }_2}{2}{\lambda }_{\mathit{xx}}-\frac{2\rho }{{\xi }_0}{(1+{\xi }_1{\lambda }_{\mathit{xx}}/2)}^2){}_0{\stackrel{\sim }{F}}_1(2;\frac{4\rho }{{\xi }_0^2}{(1+{\xi }_1{\lambda }_{\mathit{xx}}/2)}^2).\hfill \end{array}$$(D.31)

We are now in position to compute the number densities of extrema of the various kinds. The Euler number density is obtained from the computation of Eq. (D.3) which can be obtained directly via the derivative of χ_1D(ρ, λ_xx) defined in (D.34) with respect to λ_xx at position λ_xx = 0,

$$\begin{array}{c}\hfill n_{1\mathrm{D}}^{\mathrm{Euler}}(\rho )={\frac{\partial }{\partial {\lambda }_{\mathit{xx}}}}_{{|}_{{\lambda }_{\mathit{xx}}=0}}{\chi }_{1\mathrm{D}}(\rho ;{\lambda }_{\mathit{xx}}).\end{array}$$(D.32)

This finally leads to the following remarkable closed form expression for the 1D Euler number density of Levy flights

$$\begin{array}{cc}\hfill n_{1\mathrm{D}}^{\mathrm{Euler}}(\rho )& =\sqrt{\frac{2{\xi }_1}{\pi \rho }}\frac{e^{-\frac{2(\rho +1)}{{\xi }_0}}}{{\xi }_0^3}\hfill \\ \hfill & \times \{2{\xi }_0{}_0{\stackrel{\sim }{F}}_1(1;\frac{4\rho }{{\xi }_0^2})-({\xi }_0+4\rho ){}_0{\stackrel{\sim }{F}}_1(2;\frac{4\rho }{{\xi }_0^2})\}.\hfill \end{array}$$(D.33)

This is to be contrasted to the number density of Euler points in the Gaussian limit, which reads

$$\begin{array}{c}\hfill n_{1\mathrm{D}}^{\mathrm{Gaussian}}(\rho )=\frac{\sqrt{{\xi }_1}(1-\rho )e^{-\frac{{(\rho -1)}^2}{2{\xi }_0}}}{2\pi {\xi }_0^{3/2}},\end{array}$$(D.34)

which can also be obtained from more direct calculations (see, e.g. Bardeen et al. 1986). In general, it is of interest to be able to distinguish between the number counts of maxima and minima separately. This can be achieved via specific integration path in the complex plane as follows

$$\begin{array}{c}\hfill n_{1\mathrm{D}}^{\pm }(\rho )={\displaystyle {{\displaystyle \int }}_{-i\infty \pm ϵ}^{+i\infty \pm ϵ}}\frac{\mathrm{d}{\lambda }_{\mathit{xx}}}{2\pi i}\frac{1}{{\lambda }_{\mathit{xx}}^2}{\chi }_{1\mathrm{D}}(\rho ;{\lambda }_{\mathit{xx}}),\end{array}$$(D.35)

where χ_1D(ρ; λ_xx) is given by Eq. (D.34), while ϵ is a negative real constant for maxima and a positive real constant for minima¹⁴. This quadrature, together with its Euler counterpart, Eq. (D.36), is one of the key result of this paper.

D.5. 2D Euler and extrema counts

In 2D, it is best to re-organise the field variables. Indeed it is well known that while ρ_x and ρ_y form a vector field, κ = (ρ_xx + ρ_yy)/2 is a scalar field and γ₁ = (ρ_xx − ρ_yy)/2 and γ₂ = ρ_xy form a spin-2 field. It is then worth defining the CGF for those quantities and so to write it as a function of the conjugate variables of those components, that is, respectively λ_ρ, λ_g1, λ_g2, λ_κ, λ_γ1, and λ_γ2. Furthermore, because the CGF must be rotationally invariant, it can be expressed in terms of the following scalar quantities:

$$\begin{array}{cc}\hfill {\mathrm{\Lambda }}_1& ={\lambda }_{g_1}^2+{\lambda }_{g_2}^2,\hfill \end{array}$$(D.36a)

$$\begin{array}{cc}\hfill {\mathrm{\Lambda }}_2& ={\lambda }_{\kappa }={\lambda }_{\mathit{xx}}+{\lambda }_{\mathit{yy}},\hfill \end{array}$$(D.36b)

$$\begin{array}{cc}\hfill {\mathrm{\Lambda }}_3& ={\lambda }_{{\gamma }_1}^2+{\lambda }_{{\gamma }_2}^2={({\lambda }_{\mathit{xx}}-{\lambda }_{\mathit{yy}})}^2+{\lambda }_{\mathit{xy}}^2,\hfill \end{array}$$(D.36c)

$$\begin{array}{cc}\hfill {\mathrm{\Lambda }}_4& ={\lambda }_{{\gamma }_1}({\lambda }_{g_1}^2-{\lambda }_{g_2}^2)+2{\lambda }_{{\gamma }_2}{\lambda }_{g_1}{\lambda }_{g_2}.\hfill \end{array}$$(D.36d)

Following the same steps as in 1D, using MATHEMATICA to compute the small d limit of the solution of the system of Eqs. (D.23), relying on the shifting relation (Eq. C.6) and re-expressing the resulting generating function in terms of the Λ_is, the final result for the generating function is expressed as:

$$\begin{array}{cc}\hfill {\phi }_{2\mathrm{D}}& ({\lambda }_{\rho },{\mathrm{\Lambda }}_i;{\xi }_0,{\xi }_1,{\xi }_2)={\phi }_1[\frac{{\lambda }_{\rho }+\frac{1}{6}(2{\mathrm{\Lambda }}_1^2+{\mathrm{\Lambda }}_3){\xi }_2}{{(\frac{{\mathrm{\Lambda }}_1{\xi }_1}{2}+1)}^2}\hfill \\ \hfill & -\frac{{\xi }_1({\mathrm{\Lambda }}_2({\mathrm{\Lambda }}_1{\xi }_1-2)-{\mathrm{\Lambda }}_4{\xi }_1)}{{(\frac{{\mathrm{\Lambda }}_1{\xi }_1}{2}+1)}^2(({\mathrm{\Lambda }}_1^2-{\mathrm{\Lambda }}_3){\xi }_1^2-4{\mathrm{\Lambda }}_1{\xi }_1+4)};{\xi }_0].\hfill \end{array}$$(D.37)

The calculation of the number density of extrema also follows the same steps as for the 1D case. We first remark that

$$\begin{array}{c}\hfill {\phi }_{2\mathrm{D}}({\lambda }_{\rho },{\mathrm{\Lambda }}_i)={\phi }_{2\mathrm{D}}({\lambda }_e,{\mathrm{\Lambda }}_1=0,{\mathrm{\Lambda }}_2,{\mathrm{\Lambda }}_3,{\mathrm{\Lambda }}_4=0),\end{array}$$(D.38)

with

$$\begin{array}{c}\hfill {\lambda }_e={\lambda }_{\rho }+\frac{{\xi }_1(-2{\mathrm{\Lambda }}_1+{\xi }_1{\mathrm{\Lambda }}_1{\mathrm{\Lambda }}_2-{\xi }_1{\mathrm{\Lambda }}_4)}{-4+4{\xi }_1{\mathrm{\Lambda }}_2+{\xi }_1^2({\mathrm{\Lambda }}_3-{\mathrm{\Lambda }}_2^2)},\end{array}$$(D.39)

which allows us to integrate over λ_g1 and λ_g2. This leads to the joint PDF,

$$\begin{array}{cc}\hfill P(\rho ,& {\rho }_i=0,{\rho }_{\mathit{ij}})=\int \frac{\mathrm{d}{\lambda }_e}{2\pi i}{\displaystyle \prod _{\mathit{ij}}}\frac{\mathrm{d}{\lambda }_{\mathit{ij}}}{2\pi i}\frac{{(4-4{\xi }_1{\mathrm{\Lambda }}_2+{\xi }_1^2({\mathrm{\Lambda }}_2^2-{\mathrm{\Lambda }}_3))}^{1/2}}{4\pi {\xi }_1\rho }\hfill \\ \hfill & \times exp(-{\lambda }_e\rho -{\displaystyle \sum _{\mathit{ij}}}{\lambda }_{\mathit{ij}}{\rho }_{\mathit{ij}}+{\phi }_{2\mathrm{D}}({\lambda }_e,0,{\mathrm{\Lambda }}_2,{\mathrm{\Lambda }}_3,0)).\hfill \end{array}$$(D.40)

Then, defining:

$$\begin{array}{cc}\hfill a_c({\mathrm{\Lambda }}_2)& =\frac{{({\xi }_1{\mathrm{\Lambda }}_2+2)}^2}{{\xi }_0^2},\hfill \end{array}$$(D.41)

$$\begin{array}{cc}\hfill {\lambda }_c({\mathrm{\Lambda }}_2,{\mathrm{\Lambda }}_3)& =\frac{1}{2{\xi }_0}{({\xi }_1{\mathrm{\Lambda }}_2+2)}^2-\frac{{\xi }_2}{6}(2{\mathrm{\Lambda }}_2^2+{\mathrm{\Lambda }}_3),\hfill \end{array}$$(D.42)

allows us to rewrite φ_2D in Eq. (D.43) as:

$$\begin{array}{c}\hfill {\phi }_{2\mathrm{D}}({\lambda }_e,0,{\mathrm{\Lambda }}_2,{\mathrm{\Lambda }}_3,0)=-\frac{2}{{\xi }_0}+\frac{a_c({\mathrm{\Lambda }}_2)}{{\lambda }_c({\mathrm{\Lambda }}_2,{\mathrm{\Lambda }}_3)-{\lambda }_e}.\end{array}$$(D.43)

Then we can carry out the integration over λ_e in Eq. (D.46) with the help of Eq. (C.24). This leads to the following expression for the joint PDF

$$\begin{array}{c}\hfill P(\rho ,{\rho }_i=0,{\rho }_{,ij})=\int {\displaystyle \prod _{\mathit{ij}}}\frac{\mathrm{d}{\lambda }_{\mathit{ij}}}{2\pi i}{e}^{-{\lambda }_{\mathit{ij}}{\rho }_{,ij}}{\chi }_{2\mathrm{D}}(\rho ;{\mathrm{\Lambda }}_2,{\mathrm{\Lambda }}_3),\end{array}$$(D.44)

where

$$\begin{array}{cc}\hfill {\chi }_{2\mathrm{D}}& (\rho ;{\mathrm{\Lambda }}_2,{\mathrm{\Lambda }}_3)=\frac{{({\mathrm{\Lambda }}_2{\xi }_1+2)}^2\sqrt{({\mathrm{\Lambda }}_2^2-{\mathrm{\Lambda }}_3){\xi }_1^2-4{\mathrm{\Lambda }}_2{\xi }_1+4}}{4\pi {\xi }_0^2{\xi }_1\rho }\hfill \\ \hfill & \times exp(-\frac{2}{{\xi }_0}(\rho \{{(\frac{{\mathrm{\Lambda }}_2{\xi }_1}{2}+1)}^2-\frac{1}{12}(2{\mathrm{\Lambda }}_2^2+{\mathrm{\Lambda }}_3){\xi }_0{\xi }_2\}+1))\hfill \\ \hfill & \times {}_0{\stackrel{\sim }{F}}_1(2;\frac{\rho {({\mathrm{\Lambda }}_2{\xi }_1+2)}^2}{{\xi }_0^2}),\hfill \end{array}$$(D.45)

that has to be integrated over λ_κ, λ_γ1 and λ_γ2. In the following we note γ the two component vector (γ₁, γ₂) and γ_n its norm, ${\gamma }_n=\sqrt{{\gamma }_1^2+{\gamma }_2^2}$. We also note ${\lambda }_{\gamma }\equiv \sqrt{{\lambda }_{{\gamma }_1}^2+{\lambda }_{{\gamma }_2}^2}$. The resulting number density of Euler points follows through differentiation

$$\begin{array}{cc}\hfill n& {}_{2\mathrm{D}}^{\mathrm{Euler}}(\rho )={\displaystyle \int \mathrm{d}}\kappa {\mathrm{d}}^2\mathit{\gamma }({\kappa }^2-{\gamma }_n^2)P(\rho ,{\rho }_i=0,{\rho }_{\mathit{ij}})\hfill \\ \hfill & ={\{{\left(\frac{\partial }{\partial {\lambda }_{\kappa }}\right)}^2-{\left(\frac{\partial }{\partial {\lambda }_{{\gamma }_1}}\right)}^2-{\left(\frac{\partial }{\partial {\lambda }_{{\gamma }_2}}\right)}^2\}}_{{|}_{{\lambda }_{\rho }={\lambda }_{{\gamma }_1}={\lambda }_{{\gamma }_2}=0}}{\chi }_{2\mathrm{D}}(\rho ;{\mathrm{\Lambda }}_2,{\mathrm{\Lambda }}_3),\hfill \end{array}$$(D.46)

which after some algebra yields the following final closed form for the 2D Euler count of Rayleigh-Levy flights

$$\begin{array}{cc}\hfill n_{2\mathrm{D}}^{\mathrm{Euler}}(\rho )=& -\frac{{\xi }_1{e}^{-\frac{2(\rho +1)}{{\xi }_0}}}{\pi {\xi }_0^4\rho }\{{\xi }_0({\xi }_0+8\rho ){}_0{\stackrel{\sim }{F}}_1(1;\frac{4\rho }{{\xi }_0^2})\hfill \\ \hfill & -(2{\xi }_0\rho +{\xi }_0^2+8\rho (\rho +1)){}_0{\stackrel{\sim }{F}}_1(2;\frac{4\rho }{{\xi }_0^2})\},\hfill \end{array}$$(D.47)

to be once again contrasted to the Gaussian counts

$$\begin{array}{c}\hfill n_{2\mathrm{D}}^{\mathrm{Gaussian}}(\rho )=-\frac{{\xi }_1{e}^{-\frac{{(\rho -1)}^2}{2{\xi }_0}}({\xi }_0-{(\rho -1)}^2)}{2\sqrt{2}{\pi }^{3/2}{\xi }_0^{5/2}}.\end{array}$$(D.48)

The computation of extrema number densities in 2D require us to distinguish the signs of the eigenvalues of the Hessian. In practice that means one should be able to split the integral over κ into 3 domains, κ < −γ_n, −γ_n < κ < γ_n and κ < γ_n. We first note that

$$\begin{array}{c}\hfill {\displaystyle {{\displaystyle \int }}_{-{\gamma }_n}^{{\gamma }_n}}\mathrm{d}\kappa ({\kappa }^2-{\gamma }_n^2)e^{i\kappa {\lambda }_{\kappa }}=-\frac{4}{{\lambda }_{\kappa }^3}\sqrt{\frac{\pi }{2}}{\left({\gamma }_n{\lambda }_{\kappa }\right)}^{3/2}{J}_{3/2}\left(\gamma {\lambda }_{\kappa }\right),\end{array}$$(D.49)

which implies that

$$\begin{array}{cc}& {\displaystyle {\int }_{{\mathbb{R}}^2}}{\mathrm{d}}^2\mathit{\gamma }{\displaystyle {{\displaystyle \int }}_{-{\gamma }_n}^{{\gamma }_n}}\mathrm{d}\kappa ({\kappa }^2-{\gamma }_n^2)e^{i\kappa {\lambda }_{\kappa }+i{\gamma }_1{\lambda }_{{\gamma }_1}+i{\gamma }_2{\lambda }_{{\gamma }_2}}\hfill \\ \hfill & =-8\pi \sqrt{\frac{\pi }{2{\lambda }_{\kappa }^3}}{\displaystyle {\int }_0^{\infty }}\mathrm{d}\gamma {\gamma }^{5/2}{J}_0\left(\gamma {\lambda }_{\gamma }\right)J_{3/2}\left(\gamma {\lambda }_{\kappa }\right),\hfill \\ \hfill & =-12\pi {({\lambda }_{\gamma }^2-{\lambda }_{\kappa }^2)}^{-5/2},\hfill \end{array}$$(D.50)

which leads to the following quadrature for the extrema counts

$$\begin{array}{c}\hfill n_{2\mathrm{D}}^{\mathrm{extr}.}(\rho )={\displaystyle {{\displaystyle \int }}_{-i\infty +ϵ}^{+i\infty +ϵ}}\frac{\mathrm{d}{\lambda }_{\kappa }}{2\pi i}{\displaystyle {{\displaystyle \int }}_0^{\infty }}\frac{\mathrm{d}{\lambda }_{\gamma }}{2\pi }\frac{12\pi {\lambda }_{\gamma }}{{({\lambda }_{\gamma }^2+{\lambda }_{\kappa }^2)}^{5/2}}{\chi }_{2\mathrm{D}}(\rho ;{\mathrm{\Lambda }}_2,{\mathrm{\Lambda }}_3),\end{array}$$

where ϵ is a negative real constant for maxima and a positive real constant for minima. This double quadrature, together with its Euler counterpart (D.53), is also one of the key result of this paper.

D.6. The 3D Euler and extrema counts

The 3D Euler and extrema counts can be similarly derived, although the calculation is somewhat more sophisticated. To sort out the results, we strive to take full advantage of the rotation invariance of the result and to be a bit more precise. The quantities, ρ, ρ_,i and ρ_,ij form respectively a tensor of rank 0, 1, and 2 with respect to coordinate changes ${\mathcal{R}}_i^j$. Therefore we would like to derive the conjugate expression ${\widehat{\lambda }}_i$ and ${\widehat{\lambda }}_{\mathit{ij}}$ of the quantities obtained after a change of coordinates,

$$\begin{array}{c}\hfill {\widehat{\rho }}_{,i}={\mathcal{R}}_i^{i^{\prime }}{\rho }_{,i^{\prime }},{\widehat{\rho }}_{,ij}={\mathcal{R}}_i^{i^{\prime }}{\mathcal{R}}_j^{j^{\prime }}{\rho }_{,i^{\prime }{j}^{\prime }}.\end{array}$$(D.51)

We should have

$$\begin{array}{c}\hfill {\displaystyle \sum _i}{\lambda }_i{\rho }_{,i}={\displaystyle \sum _i}{\widehat{\lambda }}_i{\widehat{\rho }}_{,i},{\displaystyle \sum _{i\ge j}}{\lambda }_{\mathit{ij}}{\rho }_{,ij}={\displaystyle \sum _{i\ge j}}{\widehat{\lambda }}_{\mathit{ij}}{\widehat{\rho }}_{,ij}.\end{array}$$(D.52)

Replacing ${\widehat{\rho }}_{,i}$ by its expression in the right hand side of Eq. (D.58) ensures that

$$\begin{array}{c}\hfill {\lambda }_i={\mathcal{R}}_i^{i^{\prime }}{\widehat{\lambda }}_{,i^{\prime }},\end{array}$$(D.53)

making clear that λ_i is a rank-1 co-tensor. This is not the case for λ_ij as the sum in Eq. (D.58) excludes repeated symmetric quantities. We are then led to define

$$\begin{array}{c}\hfill {\lambda }_{\mathit{ij}}^R=\left(\begin{array}{ccc}{\lambda }_{\mathit{xx}}& {\lambda }_{\mathit{xy}}/2& {\lambda }_{\mathit{xz}}/2\\ {\lambda }_{\mathit{xy}}/2& {\lambda }_{\mathit{yy}}& {\lambda }_{\mathit{yz}}/2\\ {\lambda }_{\mathit{xz}}/2& {\lambda }_{\mathit{yz}}/2& {\lambda }_{\mathit{zz}}\end{array}\right),\end{array}$$(D.54)

so that

$$\begin{array}{c}\hfill {\displaystyle \sum _{i\ge j}}{\lambda }_{\mathit{ij}}{\rho }_{,ij}={\displaystyle \sum _{\mathit{ij}}}{\lambda }_{\mathit{ij}}^R{\rho }_{,ij}.\end{array}$$(D.55)

It is then clear that ${\lambda }_{ij}^R$ transforms as a rank-2 co-tensor. From the definition of Eq. (B.5), the CGF itself is also clearly invariant under rotations. We then expect it to depend on combinations of λ_i and ${\lambda }_{ij}^R$ that are scalar invariant. Eventually we found the expression of the CGF to depend on the following six invariant quantities

$$\begin{array}{cc}\hfill {\mathrm{\Lambda }}_1& ={\displaystyle \sum _i}{\lambda }_i{\lambda }_i,{\mathrm{\Lambda }}_2={\displaystyle \sum _i}{\lambda }_{\mathit{ii}}^R,\hfill \end{array}$$(D.56a)

$$\begin{array}{cc}\hfill {\mathrm{\Lambda }}_3& ={\displaystyle \sum _{\mathit{ij}}}{\lambda }_{\mathit{ij}}^R{\lambda }_{\mathit{ij}}^R,{\mathrm{\Lambda }}_4=det\left[{\lambda }_{\mathit{ij}}^R\right],\hfill \end{array}$$(D.56b)

$$\begin{array}{cc}\hfill {\mathrm{\Lambda }}_5& ={\displaystyle \sum _{\mathit{ij}}}{\lambda }_{\mathit{ij}}^R{\lambda }_i{\lambda }_j,{\mathrm{\Lambda }}_6={\displaystyle \sum _{\mathit{ijk}}}{\lambda }_{\mathit{ij}}^R{\lambda }_i{\lambda }_{\mathit{jk}}^R{\lambda }_k,\hfill \end{array}$$(D.56c)

with

$$\begin{array}{cc}\hfill \phi & {}_{3\mathrm{D}}({\lambda }_{\rho },{\mathrm{\Lambda }}_i;0,{\xi }_1,{\xi }_2)=\{2{\xi }_1^3(12{\mathrm{\Lambda }}_4{\lambda }_{\rho }+{\mathrm{\Lambda }}_2^2(2{\mathrm{\Lambda }}_4{\xi }_2-3{\mathrm{\Lambda }}_1)+4{\mathrm{\Lambda }}_3{\mathrm{\Lambda }}_4{\xi }_2)\hfill \\ \hfill & +2{\xi }_1^3(6{\mathrm{\Lambda }}_5{\mathrm{\Lambda }}_2+3{\mathrm{\Lambda }}_1{\mathrm{\Lambda }}_3-6{\mathrm{\Lambda }}_6)-2{\xi }_1^2({\mathrm{\Lambda }}_2^4{\xi }_2-6{\mathrm{\Lambda }}_1{\mathrm{\Lambda }}_2+6{\mathrm{\Lambda }}_5)\}\hfill \\ \hfill & -2{\xi }_1^2({\mathrm{\Lambda }}_2^2(6{\lambda }_{\rho }+{\mathrm{\Lambda }}_3{\xi }_2)-2{\mathrm{\Lambda }}_3(3{\lambda }_{\rho }+{\mathrm{\Lambda }}_3{\xi }_2))\hfill \\ \hfill & +4{\xi }_1(6{\mathrm{\Lambda }}_2{\lambda }_{\rho }+{\mathrm{\Lambda }}_2^3{\xi }_2+2{\mathrm{\Lambda }}_3{\mathrm{\Lambda }}_2{\xi }_2-3{\mathrm{\Lambda }}_1)-4(6{\lambda }_{\rho }+{\mathrm{\Lambda }}_2^2{\xi }_2+2{\mathrm{\Lambda }}_3{\xi }_2)\}\hfill \\ \hfill & /\{6{\mathrm{\Lambda }}_2^2{\mathrm{\Lambda }}_4{\xi }_1^5+3{\mathrm{\Lambda }}_2(-{\mathrm{\Lambda }}_2^3+{\mathrm{\Lambda }}_3{\mathrm{\Lambda }}_2+8{\mathrm{\Lambda }}_4){\xi }_1^4-24\hfill \\ \hfill & -(3(2{\mathrm{\Lambda }}_2^3-4{\mathrm{\Lambda }}_2{\mathrm{\Lambda }}_3)-24{\mathrm{\Lambda }}_4){\xi }_1^3+(6{\mathrm{\Lambda }}_2^2+12{\mathrm{\Lambda }}_3){\xi }_1^2\},\hfill \end{array}$$(D.57)

and

$$\begin{array}{c}\hfill {\phi }_{3\mathrm{D}}({\lambda }_{\rho },{\mathrm{\Lambda }}_i;{\xi }_0,{\xi }_1,{\xi }_2)={\phi }_1({\phi }_{3\mathrm{D}}({\lambda }_{\rho },{\mathrm{\Lambda }}_i;0,{\xi }_1,{\xi }_2);{\xi }_0).\end{array}$$(D.58)

As in 1 and 2D, the integral over λ_i and λ_ρ can be carried out explicitly. The computation of the Euler number density from Eq. (19) relies on the following rules for the integration over λ_ij: det[ρ_,ij] → det[∂/∂λ_ij], identifying λ_ij and λ_ji. After some significant algebra, we eventually get

$$\begin{array}{cc}\hfill n_{3\mathrm{D}}^{\mathrm{Euler}}(\rho )& =-\frac{e^{-\frac{2(\rho +1)}{{\xi }_0}}{\left({\xi }_1\rho \right)}^{3/2}}{2\sqrt{2}{\pi }^{3/2}{\xi }_0^5{\rho }^3}\hfill \\ \hfill & \times \{(6{\xi }_0^2\rho +16{\xi }_0\rho +3{\xi }_0^3+32{\rho }^2(\rho +3)){}_0{\stackrel{\sim }{F}}_1(2;\frac{4\rho }{{\xi }_0^2})\hfill \\ \hfill & -(16{\xi }_0\rho (3\rho +1)+3{\xi }_0^3){}_0{\stackrel{\sim }{F}}_1(1;\frac{4\rho }{{\xi }_0^2})\},\hfill \end{array}$$(D.59)

while the corresponding Gaussian limit is expressed as:

$$\begin{array}{c}\hfill n_{3\mathrm{D}}^{\mathrm{Gaussian}}(\rho )=-\frac{{\xi }_1^{3/2}(\rho -1)e^{-\frac{{(\rho -1)}^2}{2{\xi }_0}}({(\rho -1)}^2-3{\xi }_0)}{4{\pi }^2{\xi }_0^{7/2}}.\end{array}$$(D.60)

This completes our results regarding Euler number densities. We note that the 1D, 2D, and 3D results are strikingly similar. Their Gaussian limit is consistent with the results presented in (Codis et al. 2013), noting that ξ₁ = ⟨|∇ρ|²⟩/3 and that the Euler characteristic and the genus differ by one integration.

Appendix E: Extrema correlation

The computation of the expected correlations of extrema is based on the derivation of the joint CGF for the density and its gradients at finite distance,

$$\begin{array}{c}\hfill \phi ({\lambda }_{\rho },\{{\lambda }_i\},\{{\lambda }_{\mathit{ij}}\};{\mu }_{\rho },\{{\mu }_i\},\{{\mu }_{\mathit{ij}}\}).\end{array}$$(E.1)

What makes its estimation complicated is that it depends not only on the distance, d, but also on the ratio of d with the smoothing scale, R. It is however possible to get a simple expression in the large separation limit. To carry out the calculation, we define ξ₁(d) and ξ₂(d) as, respectively, the first- and second-order derivative of the filtered densities at distance, d, and Taylor expand w.r.t. d as:

$$\begin{array}{c}\hfill \xi (d+\delta l)\approx {\xi }_0(d)+\delta l{\xi }_1(d)+\frac{\delta l}{2}{\xi }_2(d).\end{array}$$(E.2)

Using the same approach as before we can then derive joint CGF. An intermediate result is given by the expression of the following quantity

$$\begin{array}{cc}& {\phi }^c({\lambda }_{\rho },{\lambda }_x,{\lambda }_{\mathit{xx}})\equiv {\frac{\partial }{\partial {\mu }_{\rho }}}_{{|}_{{\mu }_{\rho }=0}}\phi ({\lambda }_{\rho },{\lambda }_x,{\lambda }_{\mathit{xx}};{\mu }_{\rho })\hfill \\ \hfill & =\{-2{\xi }_0(d){\lambda }_{\rho }+{\xi }_1(-{\xi }_0(d)({\lambda }_x^2+{\xi }_2{\lambda }_{\mathrm{xx}}^3-2{\lambda }_{\rho }{\lambda }_{\mathit{xx}}))\hfill \\ \hfill & +{\xi }_1({\xi }_2(d){\lambda }_{\mathit{xx}}^2+{\xi }_0({\lambda }_x^2+{\xi }_2{\lambda }_{\mathit{xx}}^3-2{\lambda }_{\rho }{\lambda }_{\mathit{xx}}))\hfill \\ \hfill & -{\xi }_1(d){\lambda }_x({\xi }_1{\lambda }_{\mathit{xx}}+2)+{\xi }_1^2{\lambda }_{\mathit{xx}}^2({\xi }_2(d){\lambda }_{\mathit{xx}}+3)+{\xi }_0(d){\xi }_2{\lambda }_{\mathit{xx}}^2\hfill \\ \hfill & {-2{\xi }_2(d){\lambda }_{\mathit{xx}}+2{\xi }_0{\lambda }_{\rho }+{\xi }_1^3{\lambda }_{\mathit{xx}}^3-{\xi }_0{\xi }_2{\lambda }_{\mathit{xx}}^2-4\}}^2\hfill \\ \hfill & /\{2{\xi }_0{\lambda }_{\rho }+{\xi }_0{\xi }_1({\lambda }_x^2+{\xi }_2{\lambda }_{\mathit{xx}}^3-2{\lambda }_{\rho }{\lambda }_{\mathit{xx}})+{\xi }_1^3{\lambda }_{\mathit{xx}}^3+3{\xi }_1^2{\lambda }_{\mathit{xx}}^2\hfill \\ \hfill & {-{\xi }_0{\xi }_2{\lambda }_{\mathit{xx}}^2-4\}}^2,\hfill \end{array}$$(E.3)

which in principle can be used to build cross matter-extrema correlations at any separation¹⁵. It is worth investigating this expression in the large separation limit. We assume that ξ₀(d), ξ₁(d) and ξ₂(d) are all small quantities. We can further note that ξ₂(d)ξ₀ ≪ ξ₀(d)ξ₂ if d is much larger that the smoothing scale. We can then observe the following property in the large separation limit,

$$\begin{array}{c}\hfill {\phi }^c({\lambda }_{\rho },{\lambda }_x,{\lambda }_{\mathit{xx}})=1+\phi ({\lambda }_{\rho },{\lambda }_x,{\lambda }_{\mathit{xx}})+{\xi }_0(d)\phi ({\lambda }_{\rho },{\lambda }_x,{\lambda }_{\mathit{xx}}).\end{array}$$(E.4)

And more generally we can explicitly verify that

$$\begin{array}{cc}\hfill \phi ({\lambda }_{\rho },{\lambda }_x,{\lambda }_{\mathit{xx}};{\mu }_{\rho })& =\phi ({\lambda }_{\rho },{\lambda }_x,{\lambda }_{\mathit{xx}})+{\phi }_1({\mu }_{\rho })\hfill \\ \hfill & +\phi ({\lambda }_{\rho },{\lambda }_x,{\lambda }_{\mathit{xx}}){\xi }_0(d){\phi }_1({\mu }_{\rho }).\hfill \end{array}$$(E.5)

This is a property which is generic to MTM (Bernardeau 2022). More generally we expect that Eq. (34) holds. We note that the correction fully factorises in variables of each type.

In the large separation limit, the two-point number densities of critical points has correspondingly the following functional form:

$$\begin{array}{cc}\hfill n_{\mathrm{crit}.}({\rho }_1;{\rho }_2)& =n_{\mathrm{crit}.}({\rho }_1)n_{\mathrm{crit}.}({\rho }_1)\{1+b_{\mathrm{crit}.}({\rho }_1){\xi }_0(d)b_{\mathrm{crit}.}({\rho }_2)\},\hfill \end{array}$$(E.6)

where b_crit.(ρ) is defined by

$$\begin{array}{cc}\hfill b_{\mathrm{crit}.}(\rho )n_{\mathrm{crit}.}(\rho )=& {\displaystyle \int }{\displaystyle \prod _{\mathit{ij}}}\mathrm{d}{\rho }_{,ij}\mathrm{Sgn}\left[{\rho }_{,ij}\right]det\left[{\rho }_{,ij}\right]P_b(\rho ,{\rho }_{,i}=0,{\rho }_{,ij}),\hfill \end{array}$$(E.7)

where Sgn[ρ_,ij] is given in Eq. (19) and we have

$$\begin{array}{cc}\hfill P_b& (\rho ,{\rho }_{,i}=0,{\rho }_{,ij})={\displaystyle \int \frac{\mathrm{d}{\lambda }_{\rho }}{2\pi i}}\int \frac{{\mathrm{d}}^{\mathrm{D}}{\lambda }_i}{{(2\pi i)}^{\mathrm{D}}}\int {\displaystyle \prod _{\mathit{ij}}}\frac{\mathrm{d}{\lambda }_{\mathit{ij}}}{2\pi i}\phi ({\lambda }_{\rho },\{{\lambda }_i\},\{{\lambda }_{\mathit{ij}}\})\hfill \\ \hfill & \times exp(-i{\lambda }_{\rho }\rho -i{\displaystyle \sum _{\mathit{ij}}}{\lambda }_{\mathit{ij}}{\rho }_{,ij}+\phi ({\lambda }_{\rho },\{{\lambda }_i\},\{{\lambda }_{\mathit{ij}}\})).\hfill \end{array}$$(E.8)

The calculation then follows the same articulation as before.

E.1. 1D bias asymptotic

For the 1D case we are led to

$$\begin{array}{c}\hfill P_b(\rho ,{\rho }_{,x}=0,{\rho }_{,xx})={\displaystyle \int \frac{\mathrm{d}{\lambda }_{\mathit{xx}}}{2\pi i}}{e}^{-{\lambda }_{\mathit{xx}}{\rho }_{,xx}}{\chi }_b^{1\mathrm{D}}(\rho ;{\lambda }_{\mathit{xx}}),\end{array}$$(E.9)

where

$$\begin{array}{cc}\hfill {\chi }_b^{1\mathrm{D}}(\rho ;{\lambda }_{\mathit{xx}})& =\frac{2\sqrt{2}}{{\xi }_0^3\sqrt{\pi {\xi }_1\rho }}{(1-{\xi }_1{\lambda }_{\mathit{xx}})}^{1/2}{(1+{\xi }_1{\lambda }_{\mathit{xx}}/2)}^2\hfill \\ \hfill & \times exp(-\frac{2}{{\xi }_0}+\frac{{\xi }_2}{2}\rho {\lambda }_{\mathit{xx}}-\frac{2}{{\xi }_0}\rho {(1+{\xi }_1{\lambda }_{\mathit{xx}}/2)}^2)\hfill \\ \hfill & \times \{{\xi }_0{}_0{\stackrel{\sim }{F}}_1(1;\frac{4\rho }{{\xi }_0^2}{(1+{\xi }_1{\lambda }_{\mathit{xx}}/2)}^2)-2{}_0{\stackrel{\sim }{F}}_1(2;\frac{4\rho }{{\xi }_0^2}{(1+{\xi }_1{\lambda }_{\mathit{xx}}/2)}^2)\}\hfill \end{array}$$(E.10)

that derives from the use of Eqs. (C.24) and (C.25). It leads to the expression of the biased number density of Euler numbers and of maxima. We thus have

$$\begin{array}{cc}\hfill b& {}_{\mathrm{Euler}}^{1\mathrm{D}}(\rho )n_{\mathrm{Euler}}^{1\mathrm{D}}(\rho )=\sqrt{\frac{2{\xi }_1}{\pi \rho }}\frac{e^{-\frac{2(\rho +1)}{{\xi }_0}}}{{\xi }_0^4}\times \hfill \\ \hfill & \{{\xi }_0({\xi }_0-4(1+\rho )){}_0{\stackrel{\sim }{F}}_1(1;\frac{4\rho }{{\xi }_0^2})+2({\xi }_0+8\rho ){}_0{\stackrel{\sim }{F}}_1(2;\frac{4\rho }{{\xi }_0^2})\}.\hfill \end{array}$$(E.11)

which gives the behaviour of the large scale bias factor.

E.2. 2D bias asymptotic

For the 2D case the corresponding form is:

$$\begin{array}{cc}& {\chi }_b^{2\mathrm{D}}(\rho ;{\mathrm{\Lambda }}_2,{\mathrm{\Lambda }}_3)=\frac{{({\mathrm{\Lambda }}_2{\xi }_1+2)}^2\sqrt{{\mathrm{\Lambda }}_2^2{\xi }_1^2-{\mathrm{\Lambda }}_3{\xi }_1^2-4{\mathrm{\Lambda }}_2{\xi }_1+4}}{4\pi {\xi }_0^3{\xi }_1\rho }\hfill \\ \hfill & \times exp(-\frac{3\rho {({\mathrm{\Lambda }}_2{\xi }_1+2)}^2-(2{\mathrm{\Lambda }}_2^2+{\mathrm{\Lambda }}_3){\xi }_0{\xi }_2\rho +12}{6{\xi }_0})\hfill \\ \hfill & \times \{{\xi }_0{}_0{\stackrel{\sim }{F}}_1(1;\frac{\rho {({\mathrm{\Lambda }}_2{\xi }_1+2)}^2}{{\xi }_0^2})-2{}_0{\stackrel{\sim }{F}}_1(2;\frac{\rho {({\mathrm{\Lambda }}_2{\xi }_1+2)}^2}{{\xi }_0^2})\},\hfill \end{array}$$(E.12)

which, after some algebra, leads to the bias factor

$$\begin{array}{cc}\hfill b_{\mathrm{Euler}}^{2\mathrm{D}}& (\rho )n_{\mathrm{Euler}}^{2\mathrm{D}}(\rho )=\frac{2{\xi }_1{e}^{-\frac{2(\rho +1)}{{\xi }_0}}}{\pi {\xi }_0^5\rho }\{{\xi }_0({\xi }_0-3{\xi }_0\rho +4\rho (3+\rho )){}_0{\stackrel{\sim }{F}}_1(1;\frac{4\rho }{{\xi }_0^2})\hfill \\ \hfill & -({\xi }_0^2+8\rho (1+3\rho )){}_0{\stackrel{\sim }{F}}_1(2;\frac{4\rho }{{\xi }_0^2})\}.\hfill \end{array}$$(E.13)

The expressions (D.34), (D.51), (E.10) and (E.12) are the central building blocks for the construction of the extrema correlations for the 1D and 2D cases.

Appendix F: Gaussian CGFs

The cumulant generating functions in their Gaussian limits the take the following form:

$$\begin{array}{cc}\hfill \phi & {}_{1\mathrm{D}}^{\mathrm{Gaussian}}({\lambda }_{\rho },{\lambda }_x,{\lambda }_{\mathit{xx}})={\lambda }_{\rho }+\frac{1}{2}{\xi }_0{\lambda }_{\rho }^2+\frac{1}{2}{\xi }_1{\lambda }_x^2\hfill \\ \hfill & -{\xi }_1{\lambda }_{\rho }{\lambda }_{\mathit{xx}}+\frac{1}{2}{\xi }_2{\lambda }_{\mathit{xx}}^2,\hfill \end{array}$$(F.1a)

$$\begin{array}{cc}\hfill \phi & {}_{2\mathrm{D}}^{\mathrm{Gaussian}}({\lambda }_{\rho },{\lambda }_i,{\lambda }_{\mathit{ij}})={\lambda }_{\rho }+\frac{1}{2}{\xi }_0{\lambda }_{\rho }^2+\frac{1}{2}{\xi }_1{\displaystyle \sum _i}{\lambda }_i^2\hfill \\ \hfill & -{\xi }_1{\displaystyle \sum _i}{\lambda }_{\rho }{\lambda }_{\mathit{ii}}+\frac{1}{2}{\xi }_2{\displaystyle \sum _i}{\lambda }_{\mathit{ii}}^2+\frac{1}{6}{\xi }_2{\lambda }_{\mathit{xy}}^2+\frac{1}{3}{\xi }_2{\lambda }_{\mathit{xx}}{\lambda }_{\mathit{yy}},\hfill \end{array}$$(F.1b)

$$\begin{array}{cc}\hfill \phi & {}_{3\mathrm{D}}^{\mathrm{Gaussian}}({\lambda }_{\rho },{\lambda }_i,{\lambda }_{\mathit{ij}})={\lambda }_{\rho }+\frac{1}{2}{\xi }_0{\lambda }_{\rho }^2+\frac{1}{2}{\xi }_1{\displaystyle \sum _i}{\lambda }_i^2\hfill \\ \hfill & -{\xi }_1{\displaystyle \sum _i}{\lambda }_{\rho }{\lambda }_{\mathit{ii}}+\frac{1}{2}{\xi }_2{\displaystyle \sum _i}{\lambda }_{\mathit{ii}}^2+\frac{1}{6}{\xi }_2({\lambda }_{\mathit{xy}}^2+{\lambda }_{\mathit{yz}}^2+{\lambda }_{\mathit{xz}}^2)\hfill \\ \hfill & +\frac{1}{3}{\xi }_2({\lambda }_{\mathit{xx}}{\lambda }_{\mathit{yy}}+{\lambda }_{\mathit{yy}}{\lambda }_{\mathit{zz}}+{\lambda }_{\mathit{xx}}{\lambda }_{\mathit{zz}}).\hfill \end{array}$$(F.1c)

All Figures