High-precision Monte-Carlo modelling of galaxy distribution

TL;DR

This paper presents an improved Monte-Carlo method for simulating galaxy distributions with high precision, accounting for non-Gaussian fields, aliasing effects, and light cone evolution, validated through spectral and covariance analysis.

Contribution

It introduces a novel approach that eliminates the need for filtering by managing aliasing and sampling, achieving sub-percent spectral accuracy up to the Nyquist frequency.

Findings

Achieves sub-percent spectral accuracy up to Nyquist frequency.

Demonstrates effective generation of galaxy fields with specified one-point functions.

Provides insights into covariance structures relevant for large surveys.

Abstract

We revisit the case of fast Monte-Carlo simulations of galaxy positions for a non-gaussian field. More precisely we address the question of generating a 3D field with a given one-point function (as a log-normal one, but not only) and some power-spectrum fixed by cosmology. We highlight and investigate a problem that occurs when the field is filtered and identify, for the log-normal case, a regime where it can still be used. However we show that the filtering is unnecessary if one takes into account aliasing effects and finely controls the discrete sampling step. In this way we demonstrate a sub-percent precision of all our spectra up to the Nyquist frequency. We extend the method to generate a full light cone evolution comparing two methods for doing it and validate our method with a tomographic analysis. We investigate analytically and numerically the structure of the covariance matrices obtained with such simulations which may be useful for future large and deep surveys.

Equations

$$\xi (r)=\frac{1}{2{\pi }^2}{\int }_0^{\infty }dk\text{leavevmode}{k}^{2}P(k){\sin }_{c}kr,$$

$${\sigma }^2=\xi (0)=\frac{1}{2{\pi }^2}{\int }_0^{\infty }dk\text{leavevmode}{k}^{2}P(k)$$

$${\delta }_{LN}(x)$$

$${\xi }_{LN}(r)=e^{{\xi }^g(r)}-1.$$

$${\xi }_{\nu }(r)=\ln \left(1+{\xi }_{LN}(r)\right).$$

$$kP(k)=\frac{A}{k}{e}^{-k^2{R}_f^2}.$$

$$r{\xi }_{LN}(r)$$

$$=\frac{1}{2{\pi }^2}{\int }_0^{\infty }\frac{A}{k}{e}^{-k^2{R}_f^2}\sin (kr)\text{leavevmode}dk$$

$$=\frac{A}{4\pi }\mathrm{erf}\left(\frac{r}{2R_f}\right).$$

$$k{P}_{\nu }(k)=$$

$$\xi (r_{ij})={\xi }_{i,j}=\mathbb{E}\left[{\delta }_i{\delta }_j\right]$$

$$a=L/N_s,$$

$$k_N=\pi /a.$$

$$k_{+}=\frac{\pi }{a}=\frac{\pi }{2R_F}.$$

$$\delta =\mathcal{L}({\delta }_{\nu }).$$

$$⟨|{\delta }_{\text{@vec}k}{|}^2⟩=\frac{\mathcal{P}(\text{@vec}k)}{k_F^3},$$

$${\delta }_{\text{@vec}k}=\sqrt{-\mathcal{P}(\text{@vec}k)/k_f^3\ln (1-{ϵ}_1)}{e}^{2\pi {ϵ}_2},$$

$$\hat{\mathcal{P}}(\@vec{k})=\sum_{\@vec{n}}\mathcal{P}\left(|\@vec{k}-2\@vec{n}k_{N}|\right)\$$

$$\hat{C}=\frac{{\Delta }^T\Delta }{N-1},$$

$$\mathbb{V}\left[C_{ij}\right]=(C_{ij}^2+C_{ii}{C}_{jj})/(N-1).$$

$$⟨{\delta }_{\text{@vec}{k}_1}{\delta }_{\text{@vec}{k}_2}{\delta }_{\text{@vec}{k}_3}{\delta }_{\text{@vec}{k}_4}{⟩}_c={\delta }^D(\text{@vec}{k}_1+\text{@vec}{k}_2+\text{@vec}{k}_3+\text{@vec}{k}_4)T(\text{@vec}{k}_1,\text{@vec}{k}_2,\text{@vec}{k}_3,\text{@vec}{k}_4)$$

$$C_{ij}=\frac{\mathcal{P}(k_i{)}^2}{M_{k_i}}{\delta }_{ij}^D+k_F^3\bar{T}(k_i,k_j),$$

$$\bar{T}(k_i,k_j)={\int }_{k_i}{\int }_{k_j}T(\text{@vec}{k}_1,-\text{@vec}{k}_1,\text{@vec}{k}_2,-\text{@vec}{k}_2)\frac{{\mathrm{d}}^3\text{@vec}{k}_1}{V_{k_i}}\frac{{\mathrm{d}}^3\text{@vec}{k}_2}{V_{k_j}},$$

$$\begin{gathered} \bar{T}(k_i,k_i)\sim 8c_1^2\left\{4c_2^2+3c_3{c}_1\right\}{\mathcal{P}}^3(k_i)+ \\ +24\left\{3c_1^2{c}_3^2+4c_1{c}_2^2{c}_3+12c_1^2{c}_2{c}_4\right\}{\mathcal{P}}^2(k_i){\mathcal{P}}^{(2)}(k_i)+ \\ +144c_1^2{c}_3^2{\mathcal{P}}^{(2)}(0){\mathcal{P}}^2(k_i), \end{gathered}$$

$$c_n=\frac{1}{n!}{\int }_{-\infty }^{\infty }\mathcal{L}(\nu )H_n(\nu )\frac{e^{\frac{-{\nu }^2}{2}}}{\sqrt{2\pi }}d\nu ,$$

$${\Lambda }_i={\int }_{v_i}{\rho }^{(i)}(\text{@vec}x){\mathrm{d}}^3\text{@vec}x,$$

$$\hat{\mathcal{P}}(\text{@vec}k)=|W(\text{@vec}k){|}^2\sum _{\text{@vec}n}\mathcal{P}\left(|\text{@vec}k-2\text{@vec}n{k}_N|\right)+\frac{1/(2\pi {)}^3}{{\rho }_0},$$

$${\delta }_{LN}(x,z)$$

$${\delta }_{\ell }^m(r)={\int }_S\delta (r,\theta ,\varphi ){Y_{\ell }^m}^{\ast }(\theta ,\varphi )d^2\Omega ,$$

$$C_{\ell }(r,r^{\prime })=(4\pi {)}^2{\int }_0^{\infty }{k}^2\mathcal{P}(k)j_{\ell }(kr)j_{\ell }(k{r}^{\prime })dk$$

$${\tilde{\delta }}_1(\theta ,\varphi )=\int W(z)\delta (r(z),\theta ,\varphi )\mathrm{d}z.$$

Full text

††offprints: J. Bel,

Introduction

Fast Monte-Carlo methods are essential tools to design analyses over large datasets. Widely used in the Cosmological Microwave Background (CMB) community, thanks to the high quality ealpix software \citep{ealpix:2005 they are less frequently used in galactic surveys where analyses often rely on mock catalogues following a complicated and heavy process chain. The reason is that the problem is more complex since w.r.t to the CMB

• •

the galaxy distribution follows a 3D stochastic point process,

• •

the underlying continuous field is not Gaussian.

The first point, that leads to shot noise, can be accommodated although a Monte Carlo tool cannot provide universal ”window functions”, for correcting voxels effects since data do not lie on a sphere but on some complicated 3D domain. The matter distribution field cannot be Gaussian. A very simple way to see it is to note that even in the so-called ”linear” regime, i.e for scales above $\simeq 8\,h^{-1}\leavevmode\nobreak\ \rm{Mpc}$ one measures $\sigma_{8}\simeq 0.8$ which represents the standard-deviation of the matter density contrast $\delta={\rho}/{\bar{\rho}}-1$. Would the one point distribution $P(\delta)$ follow a Gaussian with such a standard deviation, the energy density $\rho=\bar{\rho}(1+\delta)$ would become negative in about $P(\delta\leq-1)=10\%$ of the cases! This very obvious argument demonstrates that even in what is called the ”linear regime”, the field is not Gaussian and follows some more evolved distribution. This is a serious problem because non Gaussian fields are difficult to characterize (Adler, 1981) and shooting samples following them is a Herculean task. Cosmologists focussed essentially on the subset of fields obtained by applying some transformation to a Gaussian one (Coles & Barrow, 1987). Remarkably, in some (rare) cases the the auto-correlation function of the transformed field has can be expressed analytically from the Gaussian one. This happens for the log-normal (LN) field, obtained essentially by taking the exponential of a Gaussian field (Coles & Jones, 1991) which largely explains the reason for its success in cosmology. Since Hubble conjectured it in 1934 (Hubble, 1934) it still describes surprisingly well the 1-point distribution of galaxies in the $\sigma<1$ regime (Clerkin et al., 2017), given that it has no theoretical foundations. A closer look, based on numerous N-body simulations, reveals it is not perfect in particular for higher variances, thus extensions with more freedom such as the skewed log-normal (Colombi, 1994) or the Gamma expansion (Gaztañaga et al., 2000). More recently, Klypin et al. (2018) proposed some more refined parameterisations. One may prefer a more physical description as the one based on a large deviation principle and spherical infall model (Uhlemann et al., 2016) that provides a fully deterministic formula for the p.d.f in the mildly non linear regime (Codis et al., 2016). Boltzmann codes as LASS \citep{Blas:2011rf}, by solving numerically the perturbation equations in the linear regime and adding some contribution from models for small scales, predict the \texttt{matter power spectrum} for a given cosmology. For any field, this quantity is always defined as the Fourier Transform of the auto-correlation function. Only in the Gaussian case does it contain all the available information. Then if we want to study cosmological parameters we need to provide realisations that follow some given spectrum. In the following, we present a method for generating a matter field (and subsequent catalogs) following any one-point function and some target power-spectrum. Although it is similar to standard methods for generating a LN field \citep[e.g][]{hiang:2013,Greiner:2015,Agrawal:2017 it is more general and solves an important issue. The way of generating a LN distribution with a target power-spectrum by transforming a Gaussian field is ill-defined when the field is smoothed, since it requires an input ”power spectrum” with some negative parts. We will focus on that problem and show its origin in Sect.1. Then we will show in Sect.2 how this problem can be cured by adjusting the discretization step and including aliasing effects. We will use the Mehler transform to show how any form of the probability distribution function (p.d.f) can be achieved still keeping a positive input power-spectrum. We then give an analytical expression of the general tri-spectrum and compare it to the output of the simulations. In Sect.2 we consider the production of a discrete catalog and how the cell window function affects the result. We discuss a linear interpolation scheme that reduces discontinuities between cells. We also consider two methods to account for the redshift evolution, one with the full light-cone reconstruction and the other evolving the perturbation, which will be compared. To qualify our catalogues we then apply in Sect.4 a tomographic analysis to compare the simulated results to the expected theoretical one and focus on the covariance matrix.

Appendices gives more technical details about some properties of the LN distribution, the Mehler transform and the tri-spectrum computation with it. Throughout the paper we target a sub-percent precision of all our spectra up to the Nyquist frequency.

1 The log-normal problem for filtered fields

The autocorrelation of the matter (over-) density field is the Fourier transform of its power-spectrum which, for an isotropic 3D field, reads

[TABLE]

where $\sin_{c}(x)=\dfrac{\sin x}{x}$ and the power-spectrum $P(k)$ can be computed with a Boltzman solver as CLASS. Technically such an integral can be computed efficiently with an FFTLog algorithm (Hamilton, 2000), which is the approach we use in the following, or simply with an FFT by noticing that $(r\xi(r),kP(k))$ form a Fourier (sine) pair. The variance of the field is by definition

[TABLE]

but looking at the shape of the spectrum (Fig. 1), one sees that the variance will increase dramatically with the wavelength and actually even not converge for a non-linear spectrum if no cutoff is introduced. This is why fields are filtered in cosmology (e.g Coles & Lucchin, 2003). This, in practice, always happen for a finite size experiment, but one may want to introduce explicitely a smoothing window as a Gaussian one, modifying $P(k)\to P(k)e^{-k^{2}R_{F}^{2}}$ which band-limits the spectrum below $k\lesssim\tfrac{2}{R_{F}}$.

Let us now recall some properties of a log-normal (LN) field, and refer the reader to Appendix A for their demonstration. Let $\delta^{g}(x)$ represent the Gaussian random field of the density contrast, then its log-normal transform in cosmology is defined as

[TABLE]

where the different factors are here to ensure the mean to be zero. Remarkably, the Gaussian auto-correlation function $\xi_{G}(r)$ transforms as

[TABLE]

This suggests a straightforward way to generate a LN field with some target power-spectrum: just log-transform (Eq. 3) a random Gaussian field with a spectrum $P_{\nu}(k)$ corresponding to

[TABLE]

that we will call an inverse-log transform. From Eq.4 the LN field should then have the desired spectrum. When performing that operation (Fig.2) a problem appears: at large $k$ the corresponding spectrum gets very small and becomes negative. This is clearly unphysical and one cannot use such a ”spectrum” to generate any Gaussian field.

A question one might ask is whether the very small negative values are due to numerical issues (as machine precision or discretization, boundaries, zero-padding in the FFT transform) or to more fundamental reasons. We can get some insight about this question by using the following model. For a smoothing radius of $R_{f}\simeq 4\,h^{-1}\leavevmode\nobreak\ \rm{Mpc}$ the Gaussian window cuts the cosmological spectrum in a region where it falls almost quadratically $P(k)\simeq\dfrac{A}{k^{2}}$ (with $A\simeq 100$ in our case). Then its filtered shape falls like

[TABLE]

The associated auto-correlation (Fourier sine transform) is then analytical. Indeed remembering that the $i/k$ operation corresponds to integration in real space, or using Gradshteyn & Ryzhik (2007, Eq. 3.954)

[TABLE]

Expanding $\ln\left(1+\xi_{LN}(r)\right)$ in series 111note that convergence of the series is not an issue because one can always normalize the Log-Normal field such that its variance is unity meaning that $|\xi_{LN}|\leq 1$ , we can compute the spectrum of the corresponding Gaussian field as

[TABLE]

The first term ($n=1$) corresponds to the Gaussian spectrum (see Eq. 6 ). The higher order terms are smaller due to the $\tfrac{1}{r^{n-1}}$ power suppression, unless the Gaussian term gets very close to zero where they can predominate and lead to a negative total contribution. We have checked this by computing the integrals numerically, keeping up to 30 terms in the sum to reach convergence, and verified that the values for $k\in[0.4,0.6]$ are indeed negative. This completely independent method proves that the negativity of the spectrum on Fig. 2 is not a numerical artefact due to some subbtle FFT effect. So how can a ”power-spectrum” be negative? The answer is actually simple: any function cannot represent an auto-correlation function, it must be positive definite (for an extensive discussion see Yoglom, 1986). If we take the autocorrelation values over a discrete set of points, we end up from its very definition

[TABLE]

to a covariance matrix that must be positive definite. The simplest way to check if the auto-correlation function is positive definite is to study the sign of its Fourier transform. Since the inverse-log transform has a Fourier transform which is not guaranteed to be positive, it cannot represent any genuine auto-correlation function. This is very different from the direct log-transform (Eq. 4) which is constructed to represent the auto-correlation of a transformed field. We must conclude that the very idea of constructing a LN field with a given filtered spectrum, by transforming a Gaussian one, is mathematically ill-defined.

2 Sampling a field with a target p.d.f and spectrum

Let us now consider the sampling of an isotropic field over a regular cubic grid of step size

[TABLE]

$N_{s}$ being the number of sampling points per dimension222From numerical considerations concerning FFT, a power of 2 is generally preferred and $L$ the comoving box size fixed by the cosmology and the maximum wanted redshift. Our goal is to obtain a proper power-spectrum up to the maximal accessible frequency which is the Nyquist one

[TABLE]

2.1 Sampling a filtered field

We fist consider the case of a filtered field, as the solid line shown on Fig.1. As we have seen in Sect.1 there is a mathematical problem for the LN field when the spectrum becomes small leading to a negative contribution to the required input one (Fig.2). One may think that the effect is so small that we can simply clip all negative values to 0 to recover a valid spectrum. Postponing the details of our full pipeline to Sect. 2.2, we have use the clipped field to generate the LN one and reconstructed the power spectrum that is compared to the target one on Fig.3. The result is unsatisfactory well before the Nyquist frequency.

Let us now try to see when the problem happens. To this purpose we counted the fraction of modes with negative values with respect to the smoothing size $R_{f}$ for 3 $N_{s}$ samplings (and fixed box size,so $a$ values). The outcome of this test is presented in Figure 4, which shows that in each case modes with negative power appear when $R_{f}\gtrsim a/2$. Put it in another way, this means that we can reconstruct a proper LN field as far as $R_{f}>a/2$. This is only partially satisfactory since we may adjust the step size (with $N_{s}$) to the smoothing radius, but we will only be able to reconstruct the spectrum up to the Nyquist frequency

[TABLE]

As shown on Fig. 1 as the dashed vertical line we will not sample the spectrum entirely and miss some (tiny) power. This affects the variance of the field although not much (about 2% in this case).

There is no fully satisfactory solution since the procedure itself is mathematically wrong whenever the spectrum reaches small values (the $k_{+}$ cutoff avoids it). The full power LN sampling of a filtered field cannot be achieved by transforming a Gaussian one. But why filter? In practice we use the simulation to generate a discrete set of galaxies and the process introduces some filtering. If we control exactly this filtering (what we will demonstrate in Sect.3 we do not need to perform it explicitly at the very start. Then we can work with an unfiltered field (as the dashed one on Fig.1) that is well behaved for the transforms. However as is clear from the figure there will always be some extra-power above the Nyquist frequency, so the key point is to handle properly aliasing.

2.2 The full pipeline

We then start from an unfiltered field. Although our method lies on a classical ground (e.g Chiang et al., 2013; Greiner & Enßlin, 2015; Agrawal et al., 2017), we introduce two new aspects:

we generalize the p.d.f to any distribution, 2. 2.

we take into account aliasing to deal with the residual power.

The idea to obtain any p.d.f (for the density contrast $\delta$) is to go into configuration space and apply a non-linear local transform to the Gaussian field ($\delta_{\nu}$)

[TABLE]

The ${\mathcal{L}}$ function can be found easily by applying standard probability transformation rules (Bel et al., 2016) and may need to be computed numerically. In the LN case, the transformation is analytical and was given in Eq.3. From now on, be $\delta(\@vec{x})$ a real, $L$ periodic and translational invariant field with null expectation value, let us define $\delta_{\@vec{k}}$ as its Fourier transform. On one hand, the translational invariance imposes that the covariance between wave modes is diagonal, i.e. $\langle\delta_{\@vec{k}}\delta_{\@vec{k}^{\prime}}\rangle=\delta^{D}(\@vec{k}+\@vec{k}^{\prime}){\mathcal{P}}({\@vec{k}})$, on the other hand the periodicity implies that the Fourier transform $\delta_{{\@vec{k}}}$ is non-zero only for ${\@vec{k}}=\@vec{n}k_{F}$, where $k_{F}=2\pi/L$ is the fundamental frequency of the field and $\@vec{n}$ is an integer vector. Adding the fact that the field is real, it follows that the expectation value of the square modulus of the Fourier transform is directly related to the power spectrum

[TABLE]

while the covariance between modes remains null. This property allows to set up a Gaussian field in Fourier space by generating as two uncorrelated centered gaussian random variables (the real and the imaginary part of the Fourier density field $\delta_{\@vec{k}}$). They must have the same variance which should be equal to half the value of the power spectrum evaluated at the considered $k$-mode. This is equivalent to generating the square of the modulus of $\delta_{\@vec{k}}$ following an exponential distribution with parameter ${\mathcal{P}}(\@vec{k})/k_{f}^{3}$ and a random phase peaked from a uniform distribution between [math] and $2\pi$. Thus, in practice the Fourier transform of a Gaussian field can be generated on a Fourier grid as

[TABLE]

with $\epsilon_{1}$ and $\epsilon_{2}$ being two uncorrelated uniformly distributed between [math] and $1$ random variables. In addition to the appealing property of having a null correlation between different modes, generating a Gaussian field in Fourier space allows to take avantage of the Fast Fourier Transforms (3D FFT) algorithm. The novel ingredient is to consider that since we are using here a ”raw” (unfiltered) cosmological spectrum, more power is leaking around the Nyquist frequency.

Then to be coherent with our process of generating the Gaussian field with the required input power spectrum we must take care of adding the aliased power as an input in Fourier space (see Hockney & Eastwood, 1988, for a detailed review) which can be performed by summing the power spectrum aliases

[TABLE]

where $\@vec{n}$ is running over the 3D Fourier wavenumbers. Note that, since the aliasing effect is mixing modes which are uncorrelated the phases remain uniformly distributed in Fourier space while the effective amplitude of the power spectrum is changing according to equation 17. In our analysis we have found that using only the first $125$ contributions from $\@vec{n}=(-2,-2,-2)$ to $\@vec{n}=(2,2,2)$ is enough to reach a percent level accuracy on the power spectrum of the catalogue at the Nyquist Frequency. A more computationally efficient choice would be to take only the first $27$ allias but it would lead to an accuracy of around 5-6%. In turn if one decide to discard all alias contributions then nearly $2$% of the modes would be required to have a negative variance. Clipping those pathological modes to 0 power, would lead to a significant deviation of the power spectrum of the generated field with respect to the expected one, even below the Nyquist frequency. In the following will consider the first $125$ alias contributions. As detailed in Bel et al. (2016) any local transform ${\mathcal{L}}$ applied to a centered Gaussian field $\nu(\@vec{x})$ corresponds to a one to one mapping $\lambda$ of its $2$-point correlation function $\xi_{\nu}(\@vec{r})\equiv\langle\nu(\@vec{x})\nu(\@vec{x}+\@vec{r})\rangle$ such that $\delta(\@vec{x})={\mathcal{L}}[\nu(\@vec{x})]$ and $\xi_{\delta}(\@vec{r})=\lambda[\xi_{\nu}(\@vec{r})]$. The $\lambda$ function is given explicitly in the Appendix (Eqs.56 and 57). As a result, using an inverse Fourier transform we can find the 3D $2$-point correlation of the target non-Gaussian field $\delta(\@vec{x})$, from wich, using the inverse mapping $\lambda^{-1}$ we are able to compute the corresponding $2$-point correlation of the Gaussian field $\nu(\@vec{x})$. In the end one only needs to Fourier transform back in order to get the input power spectrum that will be indeed characterising the input Gaussian field. In the following it will be referred to as virtual: $\mathcal{P}^{virt}(\@vec{k})$ and one can notice that being obtained from the Fourier transform of a regularly (grid) sampled $2$-point correlation it already contains aliasing effects. Thus, the input virtual Gaussian field can be generated on the corresponding Fourier grid from equation 16. A summary of the steps involved in the process computation of the virtual power spectrum are shown in figure 5. Finally, we can inverse Fourier transform the realisation of the Gaussian field and apply to it a local transformation which will automatically turn both the p.d.f and the power spectrum into the expected ones. It is clear that if the input power spectrum was not aliased (as it naturally is) then the corresponding inverse Fourier transform could not be interpreted as a regularly sampled Gaussian field, thus the process would not be self consistent. Figure (6) illustrates the power spectra involved in the generation of the non-Gaussian density field. One can see the raw input power spectrum obtained from the CLASS code and its corresponding aliased version (${\mathcal{P}}^{\delta}(k)$). Note that the aliasing needs to be applied on the $3$-dimensional Fourier grid while we represent only the averaged power spectrum in each Fourier shell of size $k_{f}$. In addition one can see on the same figure the corresponding power spectrum of the Gaussian field that we use to generate Monte-Carlo realisations. We can notice an excess of power at large $k$ which corresponds to the aliasing contribution.

In order to verify the coherency of the method, we generate $1000$ realisations of Log-Normal non gaussian fields in a periodic box of size $L=1200h^{-1}$Mpc with two different spatial resolutions corresponding to a number of sampling per side of $N_{s}=256$ and $512$. From the definition of the power spectrum (Eq. 15), we estimate the power spectrum on the $3$-dimensional Fourier grid by computing the ensemble average of the $1000$ realisations and compare it to the true expected power spectrum (${\mathcal{P}}^{\delta}(\@vec{k})$). In figure 7, we represent the $k$-shell averaged relative difference between the estimated and expected power spectrum for each individual wave modes. One can safely conclude that the accuracy of the proposed method is better than $0.1$% for wave modes close to the Nyquist frequency. In addition, no significant bias can be detected at the sub-percent level on the whole range of wave modes present in the density field independently from the choice of the spatial resolution.

2.3 Covariance matrix

A possible interest of being able to generate non-Gaussian fields with a Monte-Carlo method is related to generating a large number of realisations in order to estimate the covariance matrix (and its inverse) of a cosmological observable with a high level of statistical precision. In the following we define the shell averaged power spectrum as our observable and we estimate its covariance matrix between two shells respectively centered around wave numbers $k_{i}$ and $k_{j}$ as

[TABLE]

where $\Delta$ is a matrix formed by the residual between the estimated power spectrum in each realisation and the estimated ensemble average of the power spectrum in each $k$-shell, $\Delta_{ij}=P_{j}(k_{i})-\bar{P}(k_{i})$ and $\bar{P}(k_{i})=\frac{1}{N}\sum_{j=1}^{N}P_{j}(k_{i})$, the $j$ index refers to the $j$-th realisation. If the deviation elements $\Delta_{ij}$ follow a Gaussian distribution then one can show that the estimated covariance matrix elements $C_{ij}$ follow a Wishart distribution. As a result, the estimator 18 is unbiased and the variance of the covariance matrix elements (see Anderson, 1984) is given by

[TABLE]

In the following we show that the statistical behaviour of the variance of the estimator of the power spectrum is in agreement with what we expect. Having under-control both the target p.d.f and the power spectrum we can predict to some extent the expected covariance matrix $C$ of our power spectrum estimator. Since the density field generated with a Monte-Carlo process is non-Gaussian, the covariance matrix of the estimator of the power spectrum involves contribution of the Fourier space $4$-point correlation function. For a translational invariant density field it reduces to

[TABLE]

where $T$ is defined as the tri-spectrum which is the Fourier transform of the $4$-point correlation function in configuration space. As shown by Scoccimarro et al. (1999), the covariance matrix elements of the power spectrum estimator can be expressed as

[TABLE]

where $M_{k_{i}}$ is the number of independent modes in shell $i$ and

[TABLE]

where the integral is made over two shells of thickness $k_{F}$ centered and encapsulating respectively $k_{i}$ and $k_{j}$. The volume (in Fourier space) of each shell containing independent modes is denoted as $V_{k_{i}}$ and $V_{k_{j}}$, in the limit of thin shells we have that $V_{k}=2\pi k^{2}k_{F}$ thus $M_{k}=2\pi k^{2}/k_{f}^{2}$.

In appendix C, we show how to predict in a perturbative way the tri-spectrum of the generated non-gaussian density field for any local transform. In particular, for the contribution of the tri-spectrum to the diagonal elements of the covariance matrix we obtain the expression

[TABLE]

where $P^{(2)}(k_{i})\equiv{\mathcal{F}}[\xi^{2}]=\int\mathcal{P}(q)\mathcal{P}(|\@vec{q}+\@vec{k}_{i}|){\mathrm{d}}^{3}\@vec{q}$ and the $c_{n}$ are the coefficients of the Hermite transform of the function ${\mathcal{L}}$

[TABLE]

where $H_{n}$ denotes the probabilistic Hermite polynomial of order $n$. From equation 23 one can see that the covariance matrix elements are expected to depend on both the chosen target power spectrum and the probability density distribution of density fluctuations. We generate $7375$ realisations of a LN density field characterised by a $\Lambda$CDM power spectrum at redshift $z=0$, the $c_{n}$ are thus given analytically. We can evaluate the covariance matrix elements of the power spectrum estimator as a simple matrix product (see Eq. 18). In figure 8 we show the diagonal elements, namely the variance at each wave mode compared to the Gaussian contribution and the expected non-Gaussian contribution coming from equation 23. It confirms that for intermediate wave modes the non-Gaussian correction starts being relevant. While it fails to reproduce the full $k$-dependance due to the fact that expression 23 has been obtained in a perturbative way. In Figure 9 we show some combinations of modes $k_{i}$ and $k_{j}$ of the covariance matrix elements, they exhibit a clear dependance in such combination showing that due to the non-Gaussian nature of the created density field long and short wave modes are correlated in our power spectrum estimator. In the following we will extend the case of the continuous sampled density field to the creation of a catalogue of discrete objects, which could be galaxies, clusters haloes or simply dark matter particles.

3 Production of a catalogue

3.1 Poisson sampling

Simulating a galaxy catalogue implies transforming the sampled continuous density field $\delta(\@vec{x})$ into a point-like distribution. The density field must therefore be translated into a number of objects (galaxies, haloes or dark matter particles) per cell imposing an average number density $\rho_{0}$ in the comoving volume such that $\rho(\@vec{x})=\rho_{0}[1+\delta(\@vec{x})]$ and performing a Poisson sampling (Layzer, 1956). To do so, one must choose an interpolation scheme in order to be able to define a continuous density field $\rho^{(i)}(\@vec{x})$ between the sampling nodes $\@vec{x}_{j}$ surrounding a cell centered on position $\@vec{x}_{i}$. This way for each cell $i$ one is able to compute the expected number of object $\Lambda_{i}$ as

[TABLE]

where in practice the integration domain $v_{i}$ corresponds to the volume of a cell. Finally we assign to the cell the corresponding number of galaxies $N_{i}$ such that the probability of observing $N$ objects given the value of the underlying field $\Lambda$ is given by a Poisson distribution $P_{N}=\frac{\Lambda^{N}}{N!}e^{-\Lambda}$. This way one can distribute the right number $N_{i}$ of objects in each cell volume with a spacial probability distribution function proportional to the interpolated density field $\rho^{(i)}(\@vec{x})$ within the cell. The most straightforward interpolation scheme consists in populating cells uniformly with the corresponding number of objects, which is called the Top-Hat scheme. One can guess that on scales comparable to the size of the randomly populated cells the power spectrum of the Poisson sample won’t match the expected power spectrum. Be $\tilde{\delta}(\@vec{x})$ the sampled density contrast field (the true density contrast field multiplied by a Dirac comb) the corresponding interpolated density contrast within the cell $\hat{\delta}(\@vec{x})$ is obtained by convolving the sampled density field with a window function $W(\@vec{x})$ leading in Fourier space to the power spectrum relevant for the Poisson process as $\hat{P}(\@vec{k})=\tilde{P}(\@vec{k})|W(\@vec{k})|^{2}$ where $\tilde{P}(\@vec{k})$ is the power spectrum of the sampled density field, namely the aliased power spectrum. One can thus finally obtain that the expected power spectrum of the created catalogue will be

[TABLE]

where the additional term on the right corresponds to the shot noise contribution due to the auto-correlation of particles with themselves. Note that the Fourier transform of the chosen convolution function $W$ is cutting the power on small scale which is equivalent to smoothing the density field on the size of the cells. The interpolation scheme for the number density within each cell defines the form of the smoothing kernel $W$, in the following we consider two different interpolation scheme. The first one (the first order) is the Top-Hat which consists in defining cells around each node of the grid and assigning the corresponding density within the cell. The second (the second order) is a natural extension which consists in defining a cell as the volume within $8$ grid nodes and adopting a tri-linear interpolation scheme between the nodes. In any of the two cases the window function (see Sefusatti et al., 2016, for higher order smoothing functions) takes the general form $W^{(n)}(\@vec{k})=\left[j_{0}(k_{x}a/2)j_{0}(k_{y}a/2)j_{0}(k_{z}a/2)\right]^{n}$, where $j_{0}$ is the spherical Bessel function of order [math] and the index $n$ corresponds to the order of the interpolation scheme. We estimate the power spectrum of the catalogues with the method described by Sefusatti et al. (2016) employing a particle assignment scheme of order four (Piecewise Cubic Spline) and the interlacing technic to reduce aliasing effects. Note that these choices are intrinsic to the way we estimate the power spectrum of the distribution of generated objects and has nothing to do with the way we generate the catalogues. In figure 10, we compare the power spectra of the catalogues of objects in case of the two interpolation schemes described above. We see that as expected the linear interpolation scheme reduces more the extra power (due to aliasing) on small scales. In the same figure we also show the expected power spectra computed with equation 26 and corresponding to the two mentioned interpolation schemes. We demonstrate in both cases that we control precisely the smoothing of the spectrum up to the Nyquist frequency and even above.

3.2 Light cone

In the following we describe how we build a light cone from our catalogue and compare two methods. Shell method: The first idea is to glue a series of comoving volume at constant time in order to reconstruct the past light cone shell by shell (Fosalba et al., 2013; Crocce et al., 2013). We first select a redshift interval $\Delta z$ labeled by $z_{{\rm min}}$ and $z_{{\rm max}}$ and a number $N_{{\rm shl}}$ of shells within it. For each of these shell, we generate a point-like distribution in a comoving volume at constant cosmic time. Obviously, we perform the poisson sampling only of the cells contributing to the considered redshift shell. In addition, we keep only the objects belonging to the comoving volume spanned by the redshift shell, defined by $[R(z_{i}-dz/2),R(z_{i}+dz/2)]$ where $z_{i}$ corresponds to the redshift of the comoving volume, $dz=\Delta z/N_{shl}$ and $R(z)$ is the radial comoving distance. The light cone covers $4\pi$ steradians of the sky. In the next section we will show the effect of the choice of the number of shells $N_{{\rm shl}}$ used to build the light cone, on the angular power spectrum. Cell method: The second method is faster. Rather than simulating many redshift-shells, one selects a single redshift $z_{0}$ chosen at the middle of the radial comoving space spanned by the light cone. We generate the corresponding Gaussian field in a comoving volume on a grid at $z=z_{0}$. At this level one needs to include some evolution in the radial direction from the point of view of an observer located at the center of the box. To do so, there are two possibilities.

• •

One can simply think of rescaling the Gaussian field at a comoving radial distance $x(z)$ (from the observer) with the corresponding growth factor $D(z)$ which rules the evolution (Peebles, 1980) of linear matter perturbations. In this way it is clear that as on large scale the power spectrum of the density field will follow the expected evolution in $D^{2}(z)$, however the small scales will be affected in a non trivial way leading to a modification of the shape of the power spectrum,

• •

or one can change the contrast field so that the evolution of the density field will follow the growth factor $D(z)$. For the LN case this would read

[TABLE]

The second option is particularly well suited when generating a density field following a linear evolution. However the first one, although not exact, can allow for the fast computation of spectra evolution for more complex cases as when $D(k,z)$. In the following section we compare, in the linear regime, the shell-method and the two cell-methods in the case of the Log-Normal density field.

4 Application to tomography

In cosmology, several arguments can be put forward to justify a tomographic approach. Unlike the estimation of the power spectrum or the $2$-point correlation function, no fiducial cosmology needs to be assumed in order to estimate the observable (see Bonvin & Durrer, 2011; Montanari & Durrer, 2012; Asorey et al., 2012). Only angular observed positions and measured redshift are required, thus making it a true observable quantity. In addition, the observable is defined on a sphere simplifying its combinations with other cosmological probes such as lensing (Cai & Bernstein, 2012; Gaztañaga et al., 2012) or CMB and H$\alpha$ intensity mapping.

4.1 Angular power spectrum $C_{\ell}$

So far we worked on a Fourier basis but it is usefull to expand the matter perturbations into spherical harmonics (Peebles, 1980) and consider its coefficients

[TABLE]

Assuming that the field is statistically invariant by rotation, i.e. its angular $2$-point correlation function only depends on the angular separation and not on the absolute angular position (the analogue of translational invariance in 3D) on the sky, the two point correlation of the harmonic coefficients depends only on the order $\ell$, thus $C_{\ell}(r,r^{\prime})\equiv\left<\delta_{\ell}^{m}(r){\delta_{\ell}^{m}}^{\star}(r^{\prime})\right>$ is defined as the angular power spectrum between shells $r$ and $r^{\prime}$. We may relate this spectrum to the isotropic 3D one as

[TABLE]

where $j_{\ell}$ is the spherical Bessel function of order $\ell$. However, in this expression there is an explicit dependence on the radial comoving distances $r$ and $r^{\prime}$. One can project the density field of a thick redshift shell with a given weighting function $W(z)$ and define our observable density field as

[TABLE]

The corresponding angular power spectrum 333in the following we only consider auto-correlations between shells can be predicted from

[TABLE]

In practice, the numerical evaluation of equation 31 is not simple and we will use the Angpow software (Campagne et al., 2018) which is fully optimised for this task. This theoretical quantity may then be compared to our simulations by considering the number counts within pixels from samples in the $z_{1}<z<z_{2}$ range (i.e using for $W$ a top-hat window) since in our case the only source of fluctuations is due to the over-density field. We then simply project the objects of the catalogue on the sky, count them, normalize by the mean value within pixels $\bar{N}$, and compute the spherical power-spectrum $\hat{C}_{\ell}$ with the Healpix (Górski et al., 2005) software using the parameter $\text{nside}=2^{11}$. The shot-noise contribution is classically $\tfrac{1}{\bar{N}}$. Note that as for the figure (10), one compute the angular power spectrum at scales smaller than the grid pitch. In a spherical basis, the equivalent of the Nyquist mode is obtained using $\ell_{N}\sim R[z_{\text{mean}}]k_{N}$ where $z_{\text{mean}}$ is the averaged redshift of the particules composing the catalogue. In Figure 11 we compare the estimated angular power spectrum in the $[0.2,0.3]$ redshift range to the predicted one (Eq. 31) with the shell method described in section 3 ($N_{shl}=250)$ for thousand generated light-cones. In the lower panel of the same figure we display the relative difference between the two showing that the agreement is better than the percent level. In order to quantify the impact of the choice of the number of shells in the shell method, we run a test comparison between the cell method and various number of shells $N_{shl}$ in the cell method. Note that since we are using a power spectrum which is evolving linearly across cosmic time we know that the shell method is expected to converge to cell method if we rescale the density field with the linear growing mode $D(z)$ as described in section 3. Figure 12 shows the outcome of this analysis, which shows that in the considered redshift range the shell method is indeed converging to the cell method below the percent level as long as the number of shell is higher than $200$. Finally we make a comparison within the cell method, rescaling the Gaussian field instead of rescaling the density field. This way we can quantify the deviation when assuming that the Gaussian field evolves linearly when instead it is the density field which is evolving linearly. In Figure 13, we show the relative deviation in the two cases with respect to the expected power spectrum. One can see that the deviation, despite being systematic, remains small (around the percent level). Therefore considering these two cell-methods and as stated in the previous section, only the one offering better results will be recommended: the rescaling of the density field (top panel).

4.2 Covariance matrix

In this section we consider the cell method with linear rescaling of the density field. We aim at estimating the covariance matrix of the $\hat{C}_{\ell}$ estimator defined as

[TABLE]

with a high level of precision. Let us first show that the covariance matrix has a similar structure as the one of the power spectrum estimator studied in section 2. By definition the covariance of $\hat{C}_{\ell}$ is

[TABLE]

where we can substitue $\hat{C}_{\ell}$ with its expression (see Eq. 32). One immediately see that the first term of equation 33 will let appear a $4$-point moment which can be expanded (Fry, 1984a, b) in terms of cumulent moments (or connected expectation values). It follows that it takes the general form

[TABLE]

where $\bar{T}_{\ell\ell^{\prime}}$ accounts for non-Gaussian contribution. Instead, in case of $\tilde{\delta}_{\ell}^{m}$ being a Gaussian field one can see that the covariance matrix is diagonal. We generate $N=10000$ realisations and measure the angular power spectrum in each of them, in order to finally estimate the covariance matrix in the same way as described in section 2. In Figure 14 we represent the diagonal of the covariance matrix, the errors on the covariance matrix elements are computed with equation 19. Since in the Gaussian case the relative error expected on the diagonal of the covariance matrix elements is given by $\sqrt{2/(N-1)}$, the interest of using such a large number of realisation is that we expect a $1.4$% precision on the estimation of the diagonal of the covariance matrix and an absolute precision on the correlation coefficients $r_{ij}=C_{ij}/\sqrt{C_{ii}C_{jj}}$ of roughly $0.02$. In the bottom panel of figure 14 we show the relative deviation between the Gaussian prediction and the measured variance of the angular power spectrum, we see that the maximum of deviation is about $15$% at $\ell\sim 600$. It appears that deviations from Gaussianity remain small compared to the deviation obtained for the power spectrum covariance matrix (see section 2).

In addition, in Figure 15 we display some off-diagonal covariance elements with their error bars. Despite some fluctuations it is consistant with zero indicating that the covariance matrix is close to be diagonal as expected in the Gaussian case (at least for the 300 firsts elements of the matrix by counting them following the description in caption). In order to make sure that this is indeed the case, in Figure 16 we show the correlation coefficients $r_{ij}=C_{ij}/\sqrt{C_{ii}C_{jj}}$, we see that the matrix is close to be diagonal only considering $\ell<200$. It therefore confirms that projecting a thick redshift shell onto the sky tends to turn more Gaussian the density field. This is coherent with what one would naively expect from the central limit theorem, since the projection is made by summing over many values of a non-Gaussian field with some weights its appears that the resulting distribution should tend to a Gaussian as the volume of the projection increases. However, for large $\ell$-values we measure a significant amount of correlation typically of order $10$% reaching $30$% at $\ell\sim 600$.

Conclusion

In this paper we refined a known process allowing to generate a non-Gaussian density field with a given p.d.f and power spectrum. We first pointed out the main current mathematical issue arising when one wants to generate a density field with a cut-off scale by filtering its power spectrum. Indeed, we demonstrated that the power spectrum of the Gaussian field that will eventually be transformed into a non-Gaussian is likely to be undefined (i.e. with negative values) on some bandwidth. Despite the fact that, we stressed that there is in principle no way of sorting mathematically this problem we have shown that a simple criterion allows to get around: for a Gaussian filtering, one needs to take a spatial sampling rate $a$ which is larger than twice the cutoff scale $R_{f}$. We demonstrated that taking into account aliasing at the stage of generating the density field in Fourier space is of paramount importance in order to maintain the output power spectrum under control. In addition, we have shown that without imposing an explicit cutoff scale, at the stage of producing a catalogue with a local Poisson process, we introduce an effective filtering of the density field which can be predicted with a subpercent level accuracy. Regarding the Poisson sampling, we proposed a natural extension of the usual Top-Hat method consisting in populating the cubical cells uniformly with objects: we can linearly interpolate the density field between nodes and populate the cells with a probability distribution following the interpolated density field. The interest of this extension is that it allows to get closer to the ideal power spectrum by strongly decreasing the amplitude of aliasing. Regarding the density field, we have shown that one can predict in a perturbative way the expected bi-spectrum and tri-spectrum and provided an analytical approximation allowing to predict the variance of the power spectrum estimated for the non-Gaussian density field. This allowed us to check that the statistical behaviour of our method was going in the expected direction. At the end of section 3 we discussed two different methods to build a light cone out of our catalogues. The shell method can be used no matter if the evolution of the power spectrum is linear or not but involves a large number of redshift shell which is time consuming. The other is much faster, the cell method suits particularly well when the power spectrum evolves linearly but does not allow to keep a perfect control on the power spectrum when it presents a non-linear evolution. Finally, we presented a possible application of this kind of Monte-Carlo catalogues of objects to tomographic analysis. We have shown that the estimated angular power spectrum is in agreement at the percent level with the expected one. This validated the shell and cell methods used to build the light cones. Thanks to the numerical efficiency of this Monte-Carlo we could generate $10000$ realisations allowing to estimate the covariance matrix elements with a percent accuracy. Despite the reduction of non-Gaussianity involved in the projection of the catalogue on the sky, we are still able to detect on small scales clear signature coming from the fact that the catalogues have been generated out of a non-Gaussian density field. Such a Monte-Carlo method might be useful in investigating the dependance of the covariance matrix on the cosmological parameter. As recently done by Lippich et al. (2019) and Blot et al. (2019), in a future work we plan to compare the covariance matrix obtained with this method to the one estimated from cosmological N-body simulations and will include the treatment of redshift space distortions.

Appendix A Some properties of the LN field

Let X follow a Gaussian distribution $X\sim{\cal N}(\mu,\sigma^{2})$ then $Y=e^{X}$ follows a LN distribution. For simplicity we consider in the following that the Gaussian has a null mean $\mu=0$. Its moments can be immediately computed

[TABLE]

In particular

[TABLE]

The idea for cosmology is to ensure a positive energy density (noted $\rho$ in the following) by transforming a Gaussian density contrast ($\delta^{g}$) into

[TABLE]

One recovers the LN contrast using Eq.36

[TABLE]

This is a linear transformation of the pure LN distribution $e^{\delta^{g}}$ so we can compute immediately its first 2 moments:

[TABLE]

The random field is created by considering it a function of spatial coordinates $x_{i}$ and from now on we will use the shorthand $\delta_{i}=\delta(x_{i})$ or $\rho_{i}=\rho(x_{i})$, dropping the ”LN” subscript. Its autocorrelation function, assuming isotropy reads

[TABLE]

Calling $f_{2}(\rho_{1},\rho_{2})$ the 2D density distribution of the LN energy density random field, probability conservation yields

[TABLE]

The covariance matrix of the Gaussian field beeing

[TABLE]

One can then compute its 2-point function in a way similar to moments

[TABLE]

The last line can be obtained from a direct computation or recalling that the generative functional of a multi-dimensional Gaussian is

[TABLE]

$x,t$ representing vectors. Finally, using Eqs.38 and 36, one obtains for the contrast density of the LN field the beautiful result that

[TABLE]

One may check that the variance ($\xi(r=0)$) indeed follows Eq.41.

Appendix B The Mehler formalism

The Mehler transform for bivariate distributions is not a well known tool, while it is particularly convenient to ease computations of 2D integrals involving Gaussian distributions as was demonstrated in Simpson et al. (2013) or more recently in Bel et al. (2016). Let $(X_{1},X_{2})$ follow a central bivariate distribution

[TABLE]

with a covariance matrix

[TABLE]

For convenience we restrict the variance term to 1 , so that the covariance term $\xi_{\nu}$ is the correlation coefficient, and we will show at the end of this appendix how to treat the general case. By denoting, in loose notation, ${\mathcal{N}}(x)$ as the 1D normal distribution, the transform reads

[TABLE]

where $H_{n}$ are the (probabilistic) Hermite polynomials which are orthogonal wrt to the Gaussian measure

[TABLE]

The interest here is that when applying some local transform to a Gaussian field $Y={\mathcal{L}}(X)$ the covariance of the transformed field becomes

[TABLE]

Then if we decompose the local field onto the Hermite polynomials

[TABLE]

and use the orthogonality properties, one obtains the simple expansion

[TABLE]

where

[TABLE]

An important point to notice is that all the coefficients in the expansion are positive. Compare this to the series expansion of the inverse log-transform (Eq. 5). One sees immediately that a field with a $\ln(1+\xi_{X}$) covariance cannot be obtained from a Gaussian one. Let us now reconsider the classical log-normal field (but with $\sigma=1$). The local transform reads

[TABLE]

Then

[TABLE]

where we used the ${\mathcal{N}}(x)$ expression, $H_{0}(x)=1$ and $H_{n}(x+1)=\displaystyle{\sum_{k=0}^{n}}\binom{n}{k}H_{n-k}(x)$. From Eq.56 the autocorrelation of the LN field is

[TABLE]

in agreement with the more classical way to derive it shown in Appendix A. While unnecessary in the LN case, such an approach is very powerful in computing numerically the auto-correlation of any transformed Gaussian field. When the Gaussian field does not have a unit variance $\sigma^{2}=\xi_{X}(0)\neq 1$ which is generally the case, one works with rescaled variables leading to

[TABLE]

Then using the more general local transform discussed in Appendix A,

[TABLE]

one recovers Eq.60.

Appendix C Higher order correlation functions

In this appendix we show how to predict in a perturbative way the bi-spectrum and tri-spectrum a non-gaussian density field generated from the local transformation of a Gaussian field. Let us consider a density field $\epsilon(\@vec{x})$ in configuration space. We can therefore define its Fourier transform as

[TABLE]

As explained is section 2.2 we generate a Gaussian random field in Fourier space (assuming a power spectrum), we inverse Fourier transform it to get its analog in configuration space. We further apply a local transform ${\mathcal{L}}$ to map the Gaussian field into a stochastic field that is characterised by a target p.d.f. Thus, the $N$-point moments can be in principle predicted as soon as the local transform and the target power spectrum have been specified. Be $\nu$ a stochastic field following a centered ($\langle\nu\rangle=0$) reduced ($\sigma_{\nu}^{2}\equiv\langle\nu^{2}\rangle_{c}=1$) Gaussian distribution. From a realisation of this field, one can generate a non-Gaussian density field $\delta$ by applying a local mapping ${\mathcal{L}}$ between the two, hence

[TABLE]

Without lake of generality, one can express the $N$-point moments of the transformed density field with respect to the $2$-point correlation of the Gaussian field as

[TABLE]

where $\mathcal{B}^{(N)}$ is a $N$-variate Gaussian distribution with a $N\times N$ covariance matrix $C_{\nu}$ and sub-indexes are referring to positions $\delta_{1}\equiv\delta(\@vec{x}_{1})$. In practice the computation of equation 66 can be numerically expensive, however as shown in appendix B it can be efficiently computed thanks to the Mehler expansion, at least in the case of the $2$-point moment (see equation 56). Assuming that in the local transform ${\mathcal{L}}$ the amplitude of the coefficients of its Hermite transform (see equation 57) is decreasing with the order $n$, it offers the possibility of ordering the various contributions to the total moment. In order to evaluate equation 66 one can use extensions of the Mehler formula (Carlitz, 1970), for example the third order leads to

[TABLE]

where the correlation functions $\xi_{12}$, $\xi_{13}$ and $\xi_{23}$ are the three off-diagonal elements of the covariance matrix $C_{\nu}$ and the function $G^{N}$ is defined as an $N$-variate Gaussian distribution with a diagonal covariance matrix which values are all set to unity

[TABLE]

At fourth order ($N=4$), the $4$-variate Gaussian can be expressed as

[TABLE]

where again $\xi_{12}$, $\xi_{13}$, $\xi_{14}$, $\xi_{23}$, $\xi_{24}$ and $\xi_{34}$ are the $6$ off-diagonal elements of the covariance matrix $C_{\nu}$. By replacing equations 67 and 69 in expression 66 one can integrate over the $N$ variables $\nu_{1}$ to $\nu_{N}$ and express the $3$- and $4$-points moments as a sum over product of the two point correlation function of the Gaussian field

[TABLE]

and

[TABLE]

the coefficients $c_{i}$ are still the coefficients of the Hermite expansion defined by equation 57. Equations 70 and 71 are particularly useful when one wants to evaluate the $3$- and $4$-point correlation functions of the density field $\delta$ or their Fourier counterparts the bi-spectrum and tri-spectrum. Let us express, first, the power spectrum of the density field $\delta$ with respect to the power spectrum of the Gaussian field $\nu$. By Fourier transforming equation 56 one can obtain

[TABLE]

where the $P^{(n)}(k)$ represent what we will call loop corrections of order $n-1$ and are defined as $P^{(n)}(k)\equiv{\mathcal{F}}[\xi_{\nu}^{n}]$. The leading order or tree-level contribution is given by $c_{1}^{2}P(k)$ which is just a change in amplitude of the power spectrum of the Gaussian field. It represents the change of the power spectrum one would expect if the local transformation ${\mathcal{L}}$ was linear. We now express the $3$-point correlation function $\zeta_{\delta,123}\equiv\langle\delta_{1}\delta_{2}\delta_{3}\rangle_{c}$, dropping terms higher then $1$-loop corrections one would obtain

[TABLE]

Taking the Fourier transform of the above equation 73 one can obtain the expression of the bi-spectrum of the density field as

[TABLE]

where we use the short-hand notations $k_{i}=\@vec{k}_{i}$, $k_{ij}=|\@vec{k}_{i}+\@vec{k}_{j}|$ and

[TABLE]

which can also be expressed in terms of a triple product of the power spectrum at different wave modes

[TABLE]

In the very same way one can also express the $1$-loop prediction of the tri-spectrum, we need to start from the four-point correlation function $\eta_{\delta,1234}\equiv\langle\delta_{1}\delta_{2}\delta_{3}\delta_{4}\rangle_{c}=\langle\delta_{1}\delta_{2}\delta_{3}\delta_{4}\rangle-\xi_{\delta,12}\xi_{\delta,34}-\xi_{\delta,13}\xi_{\delta,24}-\xi_{\delta,14}\xi_{\delta,23}$ (Fry, 1984b), where we need to express the products of $2$-point correlation functions at fourth order, it follows

[TABLE]

Keeping terms of order lower or equal to four in terms of $\xi$ in equation 71 and subtracting permutations of equation 77 one can obtain the $1$-loop expression of the four-point correlation function

[TABLE]

In order to recover the correct permutations, one has to notice that $\xi_{12}$ and $\xi_{34}$, $\xi_{13}$ and $\xi_{24}$, $\xi_{14}$ and $\xi_{23}$ can be interchanged without modification of the coefficients in front, thus in the second line we have four permutations involving the product $\xi_{12}\xi_{34}$ and we can iterate three times by taking the mentioned specific pairs ($\xi_{12}\xi_{34}$, $\xi_{13}\xi_{24}$ and $\xi_{14}\xi_{23}$). The above expression can be transformed into the tri-spectrum by just taking its Fourier transform, it reads

[TABLE]

where we define the fourth order tri-spectrum as

[TABLE]

Of course, the practical evaluation of all the terms in equation 79 is not easy to get, however in order to predict the covariance matrix we only need specific configuration of the tri-spectrum the problem can simplified when trying to predict only the diagonal contribution ($k_{i}=k_{j}$). This has been already considered in literature (Scoccimarro et al., 1999) at the tree-level and we obtain the same expression which is

[TABLE]

where as in Scoccimarro et al. (1999) we approximate the angular averages the power spectrum of the sum of two wave modes as being equal to the power spectrum evaluated at half the modulus of the two wave modes. A simple extension of the previous result can be obtained by neglecting the $B^{(3)}$ and $T^{(4)}$ terms in equation 80, one would get

[TABLE]

which can be used to control the statistical behaviour of our Monte-Carlo density fields.