Precision Prediction for the Cosmological Density Distribution

TL;DR

This paper develops an analytic model for the matter density distribution in the universe, focusing on skewness and using a GEV distribution to improve the lognormal approximation, achieving high accuracy across cosmological parameters.

Contribution

It introduces a GEV-based model for the log density field that captures skewness and intrapixel correlations, enhancing the accuracy of cosmological matter distribution predictions.

Findings

GEV distribution models the log density field with 1.7% accuracy

Analytic fits for mean, variance, and skewness of the log density field

Model applicable across various redshifts and smoothing scales

Abstract

The distribution of matter in the universe is, to first order, lognormal. Improving this approximation requires characterization of the third moment (skewness) of the log density field. Thus, using Millennium Simulation phenomenology and building on previous work, we present analytic fits for the mean, variance, and skewness of the log density field $A$. We further show that a Generalized Extreme Value (GEV) distribution accurately models $A$; we submit that this GEV behavior is the result of strong intrapixel correlations, without which the smoothed distribution would tend (by the Central Limit Theorem) toward a Gaussian. Our GEV model yields cumulative distribution functions accurate to within 1.7 per cent for near-concordance cosmologies, over a range of redshifts and smoothing scales.

Equations

$$A=\ln (1+\delta ).$$

$$P_A(k)=b_A^2{P}_{\mathrm{lin}}(k).$$

$${\sigma }_{A,\mathrm{lin}}^2(k_N)\equiv {\int }_0^{k_N}\frac{dk{k}^2}{2{\pi }^2}{P}_{A,\mathrm{lin}}(k),$$

$$b_A^2=\frac{{\sigma }_A^2(k_N)}{{\sigma }_{\mathrm{lin}}^2(k_N)}.$$

$${\sigma }_A^2(k_N)=\mu \ln \left(1+\frac{{\sigma }_{\mathrm{lin}}^2(k_N)}{\mu }\right),$$

$$P_{A\mathrm{meas}}(\mathbf{k})=\sum _{\mathbf{n}\in {\mathbb{Z}}^3}{P}_A(\mathbf{k}+2k_N\mathbf{n})W(\mathbf{k}+2k_N\mathbf{n}{)}^2,$$

$${\sigma }_A^2(\ell )={\int }_{V_k\setminus \{0\}}\frac{d^{3}k}{(2\pi {)}^3}{P}_{A\mathrm{meas}}(\mathbf{k}).$$

$$⟨A⟩=-\lambda \ln \left(1+\frac{{\sigma }_{\mathrm{lin}}^2(k_N)}{2\lambda }\right);$$

$$S_3=\frac{⟨{\delta }^3⟩}{{\sigma }^4}.$$

$$T_3\equiv \frac{⟨{\left(A-⟨A⟩\right)}^3⟩}{{\sigma }_A^4(\ell )}.$$

$${\gamma }_1\equiv \frac{⟨{\left(A-⟨A⟩\right)}^3⟩}{{\sigma }_A^3(\ell )}=T_3\cdot {\sigma }_A(\ell ).$$

$$T_3=T(n_s)\cdot {\left({\sigma }_A^2(\ell )\right)}^{-p(n_s)},$$

$$T(n_s)$$

$$p(n_s)$$

$$\begin{array}[]{ccccccc}a&=&-0.70&&c&=&-0.26\\ b&=&1.25&&d&=&0.06.\end{array}$$

$$\mathcal{P}(A)=\frac{1}{{\sigma }_{\mathrm{GEV}}}t(A{)}^{1+\xi }{e}^{-t(A)},$$

$$t(A)={\left(1+\frac{A-{\mu }_{\mathrm{GEV}}}{{\sigma }_{\mathrm{GEV}}}\xi \right)}^{-1/\xi }.$$

$${\gamma }_1=-\frac{\Gamma (1-3\xi )-3\Gamma (1-\xi )\Gamma (1-2\xi )+2{\Gamma }^3(1-\xi )}{{\left(\Gamma (1-2\xi )-{\Gamma }^2(1-\xi )\right)}^{3/2}};$$

$${\sigma }_{\mathrm{GEV}}=\xi {\sigma }_A(\ell ){\left(\Gamma (1-2\xi )-{\Gamma }^2(1-\xi )\right)}^{-1/2}$$

$${\mu }_{\mathrm{GEV}}=⟨A⟩-{\sigma }_{\mathrm{GEV}}\frac{\Gamma (1-\xi )-1}{\xi }.$$

$${\Psi }_R(\hat{\rho })$$

$$\hat{\rho }(\tau )$$

$${\sigma }^2(R)$$

$${\Psi }_R(\hat{\rho })=\frac{{\nu }^2}{2{\sigma }^2(R_p)}{\left(\frac{R}{R_p}\right)}^{n_s+3}{\left({\hat{\rho }}^{\left(\frac{n_s+3}{6}\right)}-{\hat{\rho }}^{\left(\frac{n_s+3}{6}-\frac{1}{\nu }\right)}\right)}^2.$$

$$\mathcal{P}(A)=\hat{\rho }{\left[\frac{{\Psi }_R^{\prime \prime }(\hat{\rho })+{\Psi }_R^{\prime }(\hat{\rho })/\hat{\rho }}{2\pi }\right]}^{1/2}\exp \left(-{\Psi }_R(\hat{\rho })\right),$$

$$\mathcal{P}(A)=\frac{\nu }{\sigma (R_p)\sqrt{2\pi }}{\left(\frac{R}{R_p}\right)}^{\frac{n_s+3}{2}}{\hat{\rho }}^{\left(\frac{n_s+3}{6}\right)}$$

Full text

††pagerange: Precision Prediction for the Cosmological Density Distribution–Appendix

1 Introduction

The fluctuations in the cosmic microwave background (CMB) are remarkably Gaussian, implying a similarly Gaussian matter distribution at decoupling. It is the subsequent non-linear action of gravity that has produced the highly non-Gaussian distribution observable today (Fry & Peebles 1978; Szapudi et al. 1992; Gaztañaga 1994).

This distribution is approximately lognormal, a fact which emerges naturally under the assumption that peculiar velocities grow in linear fashion (Coles & Jones 1991). Lognormality is however only a rough first-order approximation (see Fig. 4), and in the framework of high-precision cosmology it is important to better characterize this distribution.

One reason for doing so is that forecasts of future surveys’ effectiveness require assumptions about the distribution of the underlying matter field – and a standard assumption is that this field is Gaussian (e.g. Wang et al. 2010). This assumption can result in significant overestimates obtainable from a survey – up to an order of magnitude for amplitude-like parameters (Repp et al. 2015).

Accurate characterization of the matter distribution is essential to recovery of this information. It is the non-Gaussianity of the overdensity field that causes a sizable fraction of the information to escape its power spectrum (Rimes & Hamilton 2005; Neyrinck & Szapudi 2007; Carron 2011; Carron & Neyrinck 2012; Wolk et al. 2013). By considering instead an “optimal observable” (Carron & Szapudi 2013), one can recapture this otherwise-lost information. However, to predict the power spectrum of the optimal observable for galaxy surveys (denoted $A^{*}$ – see Carron & Szapudi 2014), one must know the underlying dark matter distribution (Repp & Szapudi, submitted; see Section 6). Furthermore, knowledge of the matter distribution allows determination of the galaxy bias function itself (Szapudi & Pan 2004). Hence, to access all of the Fisher information in galaxy surveys, it is necessary to precisely describe the cosmological density distribution.

Various fits to the matter distribution have appeared in recent literature. Lee et al. (2017) use multiple Generalized Extreme Value (GEV) distributions to fit the density at various scales and redshifts. Shin et al. (2017) use a generalized normal distribution to describe the density for five cosmological models. Klypin et al. (2017) fit the matter distribution using a power law exponentially suppressed at both ends, considering two values of $\sigma_{8}$ and of $\Omega_{m}$, and they extend the fits to extremely small scales (80$h^{-1}$ kpc). These fits, however, do not explicitly express the distribution parameters in terms of cosmology.

Uhlemann et al. (2016) employ a first-principles approach to this problem; their results (from applying the Large Deviation Principle to spherical collapse) provide a good description of the density field in the low-variance limit. They find it helpful to consider the log density – which we denote $A=\ln(1+\delta)$ – rather than the overdensity $\delta$ itself. Repp & Szapudi (2017) also consider the log field, but they employ a phenomenological approach (using the Millennium Simulation) to predict the power spectrum $P_{A}(k)$; they find that $P_{A}(k)$ essentially equals the linear power spectrum $P_{\mathrm{lin}}(k)$, biased by the ratio of the linear variance $\sigma^{2}_{\mathrm{lin}}$ to the $A$-variance $\sigma_{A}^{2}$.

Since $\delta$ is approximately lognormal, it is reasonable to suppose in extracting the cosmological information) that the first three moments of $A$ characterize the log density. Repp & Szapudi (2017) have already characterized the variance of $A$; in this work we provide fits for its mean $\langle A\rangle$ and skewness $T_{3}$. These fits allow prediction of these moments for any (near-concordance) set of cosmological parameter values, including combinations of parameters that were not in the original fit.

We also show that the Generalized Extreme Value (GEV) distribution specified by these moments describes well the actual distribution of $A$. In particular, the GEV distribution determined by the predicted moment values is accurate to within 1.7 per cent (in near-concordance cosmologies) for redshifts up to $z\sim 2$ and for smoothing scales (pixel side lengths) from 2 to $30h^{-1}\mathrm{Mpc}$.

We structure this paper as follows: in Section 2 we reiterate our prescription for the variance of the log field and clarify the passage from $\sigma_{A}^{2}(k)$ to $\sigma_{A}^{2}(\ell)$. In Section 3 we present a similar prescription for the mean $\langle A\rangle$; and in Section 4 we present our prescription for the skewness $T_{3}$ of the log density field. In Section 5 we present the GEV model for the $A$-probability distribution and quantify its accuracy. Discussion, including comparison to Uhlemann et al. (2016), follows in Section 6; we summarize and conclude in Section 7.

2 Variance

The over/underdensity of a location in the Universe is $\delta=\rho/\overline{\rho}-1$, where $\rho$ denotes density. We here consider the log density field

[TABLE]

Since $\delta$ exhibits an approximately lognormal distribution, $A$ is Gaussian to first approximation.

Furthermore, the log transformation effectively erases non-linear evolution from the power spectrum (Neyrinck et al. 2009) – despite $A$ being only approximately Gaussian. Based on this surprising empirical fact, one can write the following relationship between the power spectrum of $A$ and the linear power spectrum predicted, for instance, by CAMB (Code for Anisotropies in the Microwave Background:111http://camb.info/ Lewis & Challinor 2002):

[TABLE]

Repp & Szapudi (2017) show that this relation holds to better than 10 per cent; the addition of a small slope modulation parameter reduces the error to a few per cent, comparable to that of the standard Smith et al. (2003) fit.

It follows immediately that the bias $b_{A}^{2}$ equals the ratio between the variance of $A$ and the linear variance. Furthermore, since these variances occur in the context of Equation 2, they are functions of the wavenumber $k$ and are obtained by integrating the power spectra against a top-hat filter in $k$-space. We thus write these quantities in terms of the Nyquist wavenumber $k_{N}=\pi/\ell$, where $\ell$ is the side length of the cubical pixels; and we denote these quantities as $\sigma^{2}_{A}(k_{N})$ and $\sigma^{2}_{\mathrm{lin}}(k_{N})$, defined by

[TABLE]

so that

[TABLE]

Repp & Szapudi (2017) use snapshots222http://gavo.mpa-garching.mpg.de/Millennium/ from the Millennium Simulation (Springel et al. 2005) to show that this $A$-variance is a simple function of linear variance:

[TABLE]

where the best-fitting value of $\mu$ is 0.73.

We stress that $\sigma_{A}^{2}(k_{N})$ is not the variance one would compute from counts-in-cells in a cosmological simulation. To elucidate this point, one can consider the $A$-field of the Millennium Simulation at $z=0$. Applying Equation 3 to the measured power spectrum $P_{A}(k)$ of this field yields a variance of 0.966; the actual counts-in-cells variance is 1.49.

Two factors are responsible for this discrepancy. First, the counts-in-cells procedure introduces convolution with a cubical top-hat filter (mass-assignment function) in real space, which would require a factor of $W(k)^{2}$ in Equation 3. Second, and more subtly, this convolution occurs at only a finite number of points and is thus liable to alias effects from wavenumbers that are even-integer multiples of the Nyquist frequency. Both effects – the convolution with the mass-assignment function and the finite sampling of the convolved density field – are significant near the Nyquist wavenumber.

To account for these effects, we follow Jing (2005) in writing

[TABLE]

where the sum runs over all three-dimensional integer vectors $\mathbf{n}$, and the Nyquist wavenumber $k_{N}=\pi/\ell$, $\ell$ being the distance between neighbouring grid points. (Note that we find it sufficient to consider only $|\mathbf{n}|<3$. Note also that this equation is valid only for $|\mathbf{k}|\leq k_{N}$.)

We now use $\sigma^{2}_{A}(\ell)$ to denote the measured counts-in-cells variance – the result of applying a cubical top-hat filter in real space. Then

[TABLE]

Note that the region of integration (denoted $V_{k}\setminus\!\{0\}$) is the set of $k$-vectors corresponding to the real-space volume under consideration; since the mean of $A$ does not vanish, we must exclude the zero-wavenumber power $P_{A}(0)$. (Technically one should also exclude modes larger than the real-space length scale; in practice the impact of doing so is negligible.)

In particular, if the spatial volume is cubical (as in the Millennium Simulation), then $V_{k}$ is a cube in $k$-space extending from $-k_{N}$ to $k_{N}$ along each axis. In such a case, integration over a sphere (i.e., replacing $d^{3}k$ with $4\pi\,k^{2}\,dk$ and integrating up to $k_{N}$) yields a value for $\sigma^{2}_{A}(\ell)$ that is too small; obtaining the correct value requires integration over the entire cube.

One additional subtlety enters here, in that $V_{k}$ includes some vectors with magnitude greater than $k_{N}$, for which Equation 6 does not hold. Comparison with direct measurements from the simulation shows that a power law continuation of Equation 6 yields sufficiently accurate values of $\sigma_{A}^{2}(\ell)$.

Jing (2005) also provides an iterative algorithm to de-alias and deconvolve the measured power, thus reconstructing the original spectrum. Fig. 1 shows that this procedure is necessary to match the measured power spectrum $P(k)$ to that generated by CAMB. It is somewhat surprising that this procedure works for the $A$-power spectrum as well, since measuring $P_{A}(k)$ involves first smoothing and then taking the log, whereas Equation 2 implicitly assumes that we first take the log – to get the ‘true’ $P_{A}(k)$ – and then smooth it (to the measured pixel scale). These operations do not in general commute; however, fig. 2 of Repp & Szapudi (2017) shows that Jing’s prescription succeeds in recovering the theoretical curve, whether one applies it to the measured spectrum of $\delta$ or of $A$.

Thus, returning to the variance measurements quoted following Equation 5, if we calculate the variance by integrating the measured $A$-spectrum (i.e., by using $P_{A\,\mathrm{meas}}$ instead of $P_{A}(k)$ in Equation 3), we obtain 0.966. If we use the deconvolved and antialiased $P_{A}(k)$, Equation 3 gives us $\sigma_{A}^{2}(k_{N})=1.18$. However, the measured counts-in-cells variance $\sigma_{A}^{2}(\ell)=1.49$, and this is the variance which Equation 7 recovers by integrating over the cube in $k$-space. The value of $\ell=\pi/k_{N}$ is the side length of the cubical pixels, which in this case is $1.95h^{-1}$ Mpc.

It is this variance $\sigma^{2}_{A}(\ell)$ which partially determines the probability distribution in Section 5. To predict it, one first obtains $\sigma^{2}_{A}(k_{N})$ from Equation 5, noting that $k_{N}$ depends on the side length of a cubical pixel in real space. One then obtains $P_{A}(k)$ from Equations 2–4, and then $\sigma^{2}_{A}(\ell)$ from Equations 6 and 7.

3 Mean

Since a simple relation connects $\sigma_{A}^{2}(k_{N})$ and $\sigma_{\mathrm{lin}}(k_{N})$, it is reasonable to hope that a similar relation might connect $\langle A\rangle$ and $\sigma_{\mathrm{lin}}^{2}(k_{N})$.

We begin by assuming that $\delta$ is lognormal in the low-variance limit, so that $\langle A\rangle=(-1/2)\sigma_{A}^{2}\sim(-1/2)\sigma^{2}_{\mathrm{lin}}.$ Hence we attempt to fit the mean value of $A$ using the function

[TABLE]

once again, $k_{N}$ is the Nyquist wavenumber associated with the side length of the cubical pixels.

We find this form to be indeed a good fit for the mean of $A$ (see Fig. 2); least squares optimization on points with $\sigma^{2}_{\mathrm{lin}}>2$ (to insure accuracy of the fit at high variances) yields a best value of $\lambda=0.65$.

By combining Equations 5 and 8, we can express $\langle A\rangle$ directly as a function of $\sigma^{2}_{A}(k_{N})$. When we do so (Fig. 3), we see that in the low-variance regime ($\sigma_{A}^{2}(k_{N})\la 0.1$) the mean is (by design) what one would expect from the lognormal approximation. At higher variances we see distinct departure from lognormal behavior.

4 Skewness

For a truly lognormal distribution, the mean and variance of $A$ would be its only nonzero moments. However, it is well-known (e.g., Colombi 1994) that the actual distribution of $A$ is noticeably skewed (see Fig. 4). We thus turn now to characterizing this skewness.

The standard measure for the skewness of the density field is $S_{3}$, the third moment of $\delta$ scaled by the square of the variance:

[TABLE]

We (following Colombi 1994) denote the analogous measure for the log field as $T_{3}$:

[TABLE]

A related measure of skewness is the Pearson’s moment coefficient $\gamma_{1}$:

[TABLE]

When we plot values of $T_{3}$ from the Millennium Simulation (see Fig. 5), we note that it is roughly constant at each scale; however, both our results and the simulations of Uhlemann et al. (2016) show a mild dependence of $T_{3}$ on the variance. We also consider that, according to a standard result from perturbation theory, $S_{3}=34/7-(n_{s}+3)$, where $n_{s}$ is the slope of the linear power spectrum at the scale being considered. While one cannot carry over the derivation of this result into log space, it is reasonable to expect a similar dependence for $T_{3}$.

We thus propose an expression for $T_{3}$ that depends both on $\sigma_{A}^{2}(\ell)$ (which varies with redshift) and on the linear power spectrum slope $n_{s}$ at the scale $\ell$. By rebinning the Millennium Simulation into pixels with side lengths 2, 4, 8, and 16 times the original pixel size (and then calculating $A$ for the rebinned volumes), we obtain $A$-fields at a variety of redshifts (from 0 to 2.1) and multiple smoothing scales. These results are available in three “near-concordance cosmologies,” (viz., the original Millennium Simulation cosmology, the WMAP7 cosmology, and the Planck 2013 cosmology), using the publicly-available rescalings from the Angulo & White (2010) algorithm.

We now plot both $T_{3}$ and $\gamma_{1}$ as functions of $\sigma^{2}_{A}$ for the various redshifts, scales, and cosmologies described above; the results appear in Fig. 5.

As noted before, $T_{3}$ is roughly constant at each smoothing scale; however, both its magnitude and its dependence on $\sigma_{A}^{2}$ are sensitive to this scale. It is at first puzzling that (for small scales) $T_{3}$ appears to rise at low variances, but this fact simply indicates that the third moment of $A$ approaches zero slightly more slowly that $\sigma_{A}^{4}(\ell)$; if we reduce by one the power of $\sigma_{A}(\ell)$ in the denominator, we obtain $\gamma_{1}$, which approaches zero as a power-law.

We thus write $T_{3}$ as a power law in $\sigma_{A}^{2}(\ell)$ with a coefficient $T(n_{s})$; since $T_{3}$ is roughly constant at each scale, we expect the power $-p(n_{s})$ to be relatively small. Thus,

[TABLE]

where $n_{s}$ is the slope of the linear power spectrum. Note that we find it necessary to use the no-wiggle linear power spectrum of Eisenstein & Hu (1998), since the inclusion of baryonic oscillations significantly affects the scale-dependence of the slope. Roughly speaking, $T(n_{s})$ is the constant part of the skewness, and the power expresses evolution with redshift.

It remains to write expressions for $T(n_{s})$ and $p(n_{s})$. We find that $T(n_{s})$ depends linearly upon $n_{s}$, and that $p(n_{s})$ is well-fit by a logarithmic function:

[TABLE]

Note that $b$ is the value in the limit of $n_{s}=-3$, corresponding to $34/7$ in the standard result for $S_{3}$.

We then perform MCMC optimization (over the measured values of $T_{3}$ for our six redshifts, five smoothing scales, and three near-concordance cosmologies) to obtain the following best fit values of $a$, $b$, $c$, and $d$:

[TABLE]

Note that relative to the expression for $S_{3}$, the coefficient of $n_{s}+3$ has changed from $-1$ to $-0.7$, and the constant term has decreased from $34/7$ to 1.25 – indicating the reduction in skewness accomplished by the log transformation.

We plot this fit (for each of our three cosmologies) in Fig. 5. While noting that it captures the trend of the data, there is scope for improvement in two areas. First, the algorithm of Angulo & White (2010) – by which the WMAP and Planck results were derived from the original Millennium Simulation – rescales the simulation volume and reassigns the redshifts of each snapshot to match the variance in a range of redshifts and scales. This procedure reproduces the power spectrum of the target cosmology quite accurately. It is also responsible for the fact that at each rebinning factor, the data points from all three cosmologies line up in Fig. 5: for instance, $z=0$ corresponds to a different snapshot in each of the three cosmologies, but these three snapshots come from the same evolutionary sequence in the same simulation. However, the closely parallel (but not coincident) curves in Fig. 5 show that our use of $n_{s}$ to parametrize the scale-dependence introduces a small cosmology-dependence that differs slightly from Angulo & White’s rescaling. Understanding (and, if necessary, correcting) this difference requires further testing in other simulations (see Section 6).

Second, the low-variance data in Fig. 5 suggest a possible downturn in $T_{3}$ as $\sigma_{A}^{2}$ approaches zero. We considered including a factor in Equation 12 to model this downturn; however, doing so would introduce additional free parameters only in order to model a regime dominated by large error bars. Thus we leave the existence of and modeling of any such downturn for future work, given that Equation 12 is sufficiently accurate for predicting the distribution of $A$.

However, it is worth noting once again how dramatically the log transform reduces the effects of non-linear evolution. The perturbative expression for $S_{3}$, viz. $34/7-(n_{s}+3)$, fails catastrophically at the scale of the full-resolution Millennium Simulation: at this scale, the measured value of $S_{3}$ (at $z=0$ in the original cosmology) is $8.4\pm 0.3$, twice as large as the value of 4.2 predicted by the perturbative expression. At the largest scale we consider (cubical pixels of side $31.25h^{-1}$ Mpc), the perturbative prediction for $S_{3}$ is better, yielding 3.4 as opposed to the measured value of $2.8\pm 0.2$. The analogous expression for $T_{3}$, with no explicit scale-dependence, is Equation 13. On our smallest scale it predicts a value of 0.81, compared to a measured value of $0.75$ (uncertainty $<.01$); on the largest scale, it predicts a value of 0.21, consistent with the measured value of $0.16\pm.09$. In this regard the log-transformed distribution exhibits much more nearly linear behavior even down to scales of $2h^{-1}$ Mpc.

In addition, Neyrinck (2013) notes that in a variety of distributions, including the lognormal, the log transform reduces the skewness by three (so that $T_{3}=S_{3}-3$). In our fits we find behavior consistent with this relationship at the largest scale ($S_{3}=2.9\pm 0.2$, $T_{3}=.16\pm.09$), though not at the smallest ($S_{3}=8.3\pm 0.3$, $T_{3}=0.75$). Furthermore, Uhlemann et al. (2016) show that the disparity between $T_{3}$ and $S_{3}-3$ is proportional to $\sigma_{A}^{2}$ (plus terms of higher order in $\sigma_{A}^{2}$). This dependence on the $A$-variance explains the fact that the smallest scale ($\sigma_{A}^{2}=1.5$) exhibits significant departure from the relationship, whereas the largest scale ($\sigma_{A}^{2}=0.18$) does not.

5 Probability Distribution

With expressions in hand for the mean, variance, and skewness of the log density field, it remains to model the actual probability distribution of $A$. We begin by noting that one does not directly sample the actual, continuous $A$ (or $\delta$) distribution; instead one first smooths the densities across each pixel before sampling. This smoothing is tantamount to summing the underlying density across the points within that pixel.

If these intrapixel densities were uncorrelated, then we would expect a Gaussian result by the Central Limit Theorem. In constrast, we know that densities are well-correlated on scales much lower than our pixel size, and therefore one would expect a different limiting distribution. Bertin & Clusel (2006) demonstrate that extreme value statistics can describe sums of correlated variables, even if the underlying process is not extremal.

There are three classes of such extreme-value distributions, namely, the Gumbel, Fréchet, and reversed Wiebull distributions (see, e.g., Gumbel 1958; Coles 2001; Leadbetter et al. 1983). Of these three distributions, the Gumbel exhibits the simplest form and possesses the advantage of support at all real numbers. Antal et al. (2009) note that the galaxy counts from SDSS seem to follow a Gumbel distribution, and this distribution is also a fairly good fit to the Millennium Simulation distribution of $A$ at low redshifts (without any rebinning).

However, the skewness ($\gamma_{1}$) of the Gumbel distribution is fixed at $\sim 1.1$, and thus it is not a good fit for $A$ at higher redshifts or at different smoothing scales. We obtain a better fit with the reversed Wiebull class of Generalized Extreme Value (GEV) distributions, which allows for variable skewness. The GEV distributions depend on three parameters – a location parameter $\mu_{\mathrm{GEV}}$, a positive scale parameter $\sigma_{\mathrm{GEV}}$, and a shape parameter $\xi$. The probability distribution is then

[TABLE]

where

[TABLE]

The Gumbel distribution is the limit of this expression as $\xi$ approaches zero, and it is clear that a double exponential is the most salient feature of the distribution in this limit.

The reversed Wiebull subclass of the GEV denotes distributions for which the shape parameter $\xi<0$. In this case, the distribution parameters depend on the moments of $A$ as follows:

[TABLE]

from this equation one can obtain $\xi$ numerically, and then the following two equations complete the specification of $\mathcal{P}(A)$:

[TABLE]

[TABLE]

This distribution has support only for $A\leq\mu_{\mathrm{GEV}}-\sigma_{\mathrm{GEV}}/\xi$; it is defined to be zero for any larger values of $A$. Thus this model implies an upper bound on a region’s density at a given epoch and scale.

This GEV distribution, together with our prescriptions for the mean, variance, and skewness of $A$, form a complete model for the log density distribution. We show the GEV fits (using moment values from Equations 7, 8, and 12) in Fig. 4. Inspection of this figure shows that the GEV distribution is a significantly better model than the Gaussian. And of course, given an accurate prescription for $A$, it is trivial to write the probability density for $\delta=e^{A}-1$.

To investigate the accuracy of this model, we compare its cumulative distribution function (CDF) to that of the Millennium Simulation realizations of $A$; we do so for our set of 6 redshifts, 5 smoothing scales, and 3 near-concordance cosmologies. Fig. 6 shows one such comparison; we note again the superiority of the GEV model to the Gaussian.

We explore two means of quantifying the differences between the true and predicted CDFs. We first consider the maximum of the absolute difference between the CDFs; this maximum is in fact the Kolmogorov-Smirnov statistic. (See discussion of the K-S test in Section 6.) Fig. 7 shows these maximum residuals, and we note that the most extreme deviation is 1.7 per cent. This deviation compares quite favorably with the Gaussian distribution, whose worst deviation exceeds 7 per cent.

We next consider the root mean square (rms) values of the CDF differences (for the same set of cosmologies, redshifts, and smoothing scales). To prevent the vanishing residuals in the tails from dominating the rms, we consider only the values of $A$ for which the simulation CDFs fall between 0.05 and 0.95. The results appear in Fig. 8; the worst GEV rms is 0.8 per cent, whereas the worst Gaussian rms is around 4 per cent. In this sense we can say that our prescription describes the log density distribution to sub-per cent levels of accuracy.

6 Discussion

In several ways our approach is complementary to that of Uhlemann et al. (2016): they employ an a priori approach, starting from the assumption of spherical collapse and employing the Large Deviation Principle; our approach on the other hand is phenomenological.

In addition, their prescription applies to the low-variance limit, whereas ours is derived from $N$-body simulations on scales where the variance is of order unity. We have already noted (Section 4) one indication that our prescription for $T_{3}$ might require correction in the low-variance regime. One might also deduce this fact from Figs. 7 and 8, which show that at large scales ($\sim 30h^{-1}$ Mpc) the GEV prescription in general does no better than the Gaussian. This scale ($k\sim 0.1$) represents the transition to the linear regime. It is this low-variance limit which Uhlemann et al. (2016) have well-characterized, and our prescription thus complements theirs by describing the mid- to high-variance regimes.

However, the two results are similar in that both prescriptions predict double-exponential behavior for $\mathcal{P}(A)$ (see Appendix). And although the functional forms differ significantly, the results yield similar predictions for the distribution of $\delta$ (compare Fig. 9 with their fig. 8).

In addition, our prescription approaches that of Uhlemann et al. in the low-variance limit (see left-hand panel of Fig. 10, in particular the $z=2.1$ distribution). Thus, for situations in which their approach is valid, our GEV prescription reproduces the phenomenology predicted by their approach.

We next note three areas in which our fit is amenable to improvement. First, we have already observed (in Section 4) the possibility of modifying our $T_{3}$-prediction at low variances.

Second, Fig. 6 shows that our model overpredicts the probability for small values of $A$ and underpredicts it for intermediate and large values. This effect is not peculiar to the redshift shown in Fig. 6; one can observe it across the entire range of redshifts and smoothing scales. It is worth noting that when we have an extremely large number of simulation data points (as in the non-rebinned snapshots), the Kolmogorov-Smirnov test discriminates decisively between the GEV and the actual $A$ distribution. This result confirms that while the GEV is a good fit (with error less than 2 per cent), it is not a perfect fit.

Third, this prescription allows us to predict the fourth moment $T_{4}$ of $A$. Doing so, we find that the predicted values differ from the measured values by at most $\sim 0.2$. In many cases, this discrepancy is comparable to the size of the measurment error bars, and in the cases for which it is not, it represents an overprediction of about 20 per cent. We conclude that our prescription yields values for $T_{4}$ that are of the correct order of magnitude but cannot form the basis for a precise prediction.

Turning now to the underlying rationale for the GEV behavior of $A$, at least two possible scenarios could produce such a distribution. Both rely on the fact that our density measurements are averages over the points in a given pixel.

The first possibility (indicated by the name of the distribution) is a process that selects extreme values from the multiple realizations (at the various points in each pixel) of the underlying distribution. If the extremes of density dominate the log smoothed density in each cell, then one could expect GEV behavior. However, in that case one would also expect this behavior to be more pronounced at larger smoothing scales, which integrate over a larger volume and thus include more-pronounced extremal values. In practice, we observe the opposite behavior: at large scales, $A$ exhibits a more nearly Gaussian distribution than at small scales. Furthermore, some simple investigations (into the effect of clipping extreme values before smoothing) indicate that the extreme values are not by themselves responsible for the GEV behavior of the log density distribution.

We have already mentioned (in Section 5) a second potential source of GEV behavior, namely, correlation of densities within a pixel. According to the Central Limit Theorem, the average of multiple realizations of independent identically distributed (and well-behaved) random variables approaches a Gaussian; however, correlations between the variables (violating independence) can yield extreme value statistics instead (Bertin & Clusel 2006). In this case, one would expect more-nearly Gaussian behavior on large scales as opposed to small scales, given that the correlations are weaker on larger scales. The Millennium Simulation pixel side lengths ($\sim 2h^{-1}$ Mpc) have the same order of magnitude as the typical correlation scale, and thus it is likely that the strong intrapixel correlations are the source of the GEV behavior in $A$. We leave further investigation of this possibility (using finer-grained simulations) to future work.

We next note that the parameters $\mu$ and $\lambda$ (in Equations 5 and 8) characterize non-perturbatively the deviations from lognormality. The variance is larger than that of a lognormal distribution, and the difference between $\mu$ and $\lambda$ corresponds to a distortion of the lognormal relationship between the mean and the variance. The Millennium Simulation (with its rescalings to two other near-concordance cosmologies) suggests that $\mu$ and $\lambda$ are universal, cosmology independent quantities characterizing the non-linear evolution. While this statement must be checked in wider range of cosmologies, it certainly appears to be true near the concordance model. These investigations we also leave to future research.

Finally, two distinct biases intervene between knowledge of the matter distribution and direct comparison to galaxy counts: the first reflects the impact of discreteness, and the second reflects the processes of galaxy formation and evolution.

Regarding the former, it is $A^{*}$, the discrete analog of $A$, which is the essentially optimal observable for discrete fields (Carron & Szapudi 2014). Since cosmological surveys count galaxies rather than directly measuring matter, it is the $A^{*}$ statistic that will ultimately reveal the information which escapes the standard power spectrum at high wavenumbers. However, the power spectrum of $A^{*}$ is biased with respect to that of $A$. Fig. 11 shows the distribution and power spectrum of $A$ as measured in the Millennium Simulation at $z=0$ in cubical pixels of side length $1.95h^{-1}$ Mpc. To generate a corresponding discrete realization, we first choose an average pixel number density of $\overline{N}=1$ count per pixel; then for each pixel we randomly assign its number of counts $N$ from a Poisson distribution with mean $\overline{N}e^{A}=\overline{N}(1+\delta)$. After calculating $A^{*}(N)$ for each cell, we obtain the one-point distribution and power spectrum of $A^{*}$, both shown in the figure. In determining the power spectra, we use Jing (2005)’s prescription to remove pixel-window and alias effects.

One notices first that $\sigma_{A^{*}}^{2}(\ell)$ is about 30 per cent lower than $\sigma_{A}^{2}(\ell)$, as expected given that much of the negative tail of the $A$-distribution collapses into the discrete $N=0$ spike. This reduced variance manifests itself in a power spectrum for $A^{*}$ that is biased with respect to that of $A$. One notices also that $P_{A^{*}}(k)$ begins to level off at high wavenumbers (an effect similar but not identical to the standard $1/\overline{n}$ discreteness plateau). Upon integrating the power spectra to obtain the variances, this upward bend partially ameliorates the effect of the bias; and thus the bias is even more pronounced than one might expect from the ratio of the variances. Approximate expressions for this bias exist (Wolk et al. 2015), but they require further refinement. We have shown (Repp & Szapudi, submitted) that the description of the $A$ field appearing in this work allows quantitative prediction of this bias $b_{A^{*}}^{2}$; and this bias is necessary to characterize the power spectrum of $A^{*}$. In addition, by predicting the probability distribution $\mathcal{P}(\delta)$, our prescription allows us to predict the distribution $\mathcal{P}(N)$ in simulations with lower particle density than the Millennium, and thus to further compare our results to those of other simulations.

This further comparison is important: thus far we have based all our results on the Millennium Simulation because its high particle density allows us to ignore discreteness effects. Unfortunately, the relatively small size of this simulation makes it susceptible to cosmic variance. Reference to other simulations will both alleviate this problem and expand the range of cosmologies which we can test.

A second bias is the galaxy bias function defined by $\delta_{g}=f(\delta)$, where $\delta_{g}$ is the galaxy overdensity. Szapudi & Pan (2004) have outlined and tested a method for determining this bias function from galaxy catalogs given the underlying dark matter distribution; one issue with which they contended was how to reconstruct (or what to assume for) that underlying distribution. Our prescription obviates this problem, allowing accurate characterization of the galaxy bias function; in particular, it permits such characterization apart from the assumption of a particular form (e.g., linear) for this bias.

7 Conclusion

The lognormal distribution is a reasonable first approximation to the matter distribution of the Universe. However, it is not an extremely accurate approximation, nor has there been until now a robust means of predicting the distribution parameters. This work remedies both problems, measuring and fitting the first three moments of the log distribution and showing that the GEV accurately describes $A=\ln(1+\delta)$. Together with the GEV prescription, these fits allow the prediction of these moments – and the one-point matter density distribution – for any (near-concordance) set of cosmological parameter values.

Both the mean $\langle A\rangle$ and the variance $\sigma_{A}^{2}(k)$ depend in a simple fashion on the linear variance $\sigma^{2}_{\mathrm{lin}}(k)$ (Equations 8 and 3–5), and the dependence reduces to the expected lognormal behavior at low variances. However, the variance $\sigma_{A}^{2}(k)$ in Equation 5 assumes a top-hat filter in $k$ space; to convert it to a measured counts-in-cells variance $\sigma_{A}^{2}(\ell)$, one must account for both the mass-assignment function and alias effects. Equations 6 and 7 quote a prescription for doing so from Jing (2005).

The skewness $T_{3}$ of the $A$-distribution depends in turn on both the variance $\sigma_{A}^{2}(\ell)$ and the slope $n_{s}$ of the no-wiggle linear power spectrum at the smoothing scale $\ell$. Describing this dependence requires four free parameters; the description appears in Equations 12 through 15.

Having expressions for the first three moments of $A$, we show that a generalized extreme value (GEV) distribution of the reversed Wiebull type matches the actual $A$-distribution well (Equations 16 through 20). When we compare the cumulative distribution functions of the actual $A$-distribution to those generated by our GEV prescription, we find the maximum difference to be less than 2 per cent, with rms differences at most 0.8 per cent.

Having thus characterized the distribution of the log density field $A$, it becomes trivial to write the distribution for the actual density perturbation field $\delta$. However, it is the power spectrum of $A$, not of $\delta$, which captures the majority of the cosmological information at small scales (Carron & Szapudi 2013) and which hence will allow us to make full use of the data from upcoming galaxy surveys.

Acknowledgements

The Millennium Simulation data bases used in this work and the web application providing online access to them were constructed as part of the activities of the German Astrophysical Virtual Observatory (GAVO). IS acknowledges support from National Science Foundation (NSF) award 1616974. We also thank Mark Neyrinck for his comments and suggestions, which significantly strengthened this paper.

Appendix

We here briefly demonstrate that both Uhlemann et al. (2016) and our GEV model predict double-exponential behavior.

We begin with equations 5, 12, and 13 of Uhlemann et al., which give us

[TABLE]

Note that $\hat{\rho}$ is the normalized density; $\hat{\rho}=\rho/\overline{\rho}=1+\delta$ in our notation. From these expressions we obtain

[TABLE]

From their equation 11 we can write

[TABLE]

so that after some algebra,

[TABLE]

Since by definition $\hat{\rho}=e^{A}$, the final line of the above equation (our Equation Appendix) will be a double exponential.

For the GEV distribution, it is clear that the function $t(A)$ in Equation 17 becomes exponential for low values of $\xi$. For the cases we consider, the values of $\xi$ are typically quite small, causing double-exponential behavior in our Equation 16.