Density distribution of the cosmological matter field

TL;DR

This paper provides a highly accurate, comprehensive analysis of the matter density PDF in the universe using extensive N-body simulations, revealing its simple shape and dependence on cosmological parameters.

Contribution

It introduces a double-exponential approximation for the matter density PDF that outperforms previous models and combines it with the spherical infall model for precise, wide-range predictions.

Findings

The PDF declines as a power-law with density, P~rho^(-2).

The double-exponential approximation fits simulation data with few percent errors.

The PDF depends on redshift and cosmological density parameter Omega_m.

Abstract

The one-point probability distribution function (PDF) of the matter density field in the universe is a fundamental property that plays an essential role in cosmology for estimates such as gravitational weak lensing, non-linear clustering, massive production of mock galaxy catalogs, and testing predictions of cosmological models. Here we make a comprehensive analysis of the dark matter PDF using a suite of 7000 N-body simulations that covers a wide range of numerical and cosmological parameters. We find that the PDF has a simple shape: it declines with density as a power-law P~rho**(-2), which is exponentially suppressed on both small and large densities. The proposed double-exponential approximation provides an accurate fit to all our N-body results for small filtering scales R< 5Mpc/h with rms density fluctuations sigma>1. In combination with the spherical infall model that works well for small fluctuations sigma<1, the PDF is now approximated with just few percent errors over the range of twelve orders of magnitude -- a remarkable example of precision cosmology. We find that at 5-10% level the PDF explicitly depends on redshift (at fixed sigma) and on cosmological density parameter Omega_m. We test different existing analytical approximations and find that the often used log-normal approximation is always 3-5 times less accurate than either the double-exponential approximation or the spherical infall model.

Tables

Table 1. Table 1: Numerical and cosmological parameters of different simulations. The columns give the simulation identifier, the size of the simulated box in h − 1 ​ Mpc superscript ℎ 1 Mpc h^{-1}\,{\rm Mpc} , the number of particles, the mass per simulation particle m p subscript 𝑚 𝑝 m_{p} in units h − 1 ​ M ⊙ superscript ℎ 1 subscript 𝑀 direct-product h^{-1}\,M_{\odot} , the mesh size N g 3 superscript subscript 𝑁 g 3 \mbox{$N_{\rm g}$}^{3} , the gravitational softening length ϵ italic-ϵ \epsilon in units of h − 1 ​ Mpc superscript ℎ 1 Mpc h^{-1}\,{\rm Mpc} , the number of time-steps N s subscript 𝑁 𝑠 N_{s} , the amplitude of perturbations σ 8 subscript 𝜎 8 \sigma_{8} , the matter density Ω m subscript Ω 𝑚 \Omega_{m} , the number of realisations N r subscript 𝑁 𝑟 N_{r}

Simulations	Box	particles	$m_p$mmm	$N_g^3$	$ϵ$	$N_s$	${\sigma }_8$	${\Omega }_m$	$N_r$
A0.5	5003	12003	$6.16\times {10}^9$	24003	0.208	181	0.822	0.307	680
A1.5	15003	12003	$1.66\times {10}^{11}$	24003	0.625	136	0.822	0.307	4513
A2.5	25003	10003	$1.33\times {10}^{12}$	20003	1.250	136	0.822	0.307	1960
B0.5	5003	16003	$2.66\times {10}^9$	32003	0.156	271	0.828	0.307	5
D0.5	5003	16003	$2.33\times {10}^9$	32003	0.156	271	0.828	0.270	5
MDPL1	10003	38403	$1.5\times {10}^9$	–	0.010	–	0.828	0.307	1
BolshoiP	2503	20483	$1.5\times {10}^8$	–	0.001	–	0.828	0.307	1

Table 2. Table 2: Parameters for the double-exponential model, eq.( 18 ), for redshift z = 0 𝑧 0 z=0 . Different columns give: (1) Name of the simulation, (2) Cell size in h − 1 superscript ℎ 1 h^{-1} Mpc , (3) r ​ m ​ s 𝑟 𝑚 𝑠 rms density fluctuation, (4-7) parameters for eq.( 18 ), (8) r ​ m ​ s 𝑟 𝑚 𝑠 rms of the relative error of the fit, (9) number of cells in 1D, (10) relative force resolution – the number of force resolution elements per each density assignment cell.

Simulation	Cell Size	$\sigma$	$A$	$\alpha +2$	${\rho }_0$	${\rho }_1$	relative error	$N_{\mathrm{cell}}$	$\Delta x/ϵ$
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)
MDPL1	0.167	24.020	1.145E-01	0.0438	0.1001	2080.000	0.066	6000	17
MDPL1	0.250	19.070	1.122E-01	0.0639	0.0925	1225.927	0.058	4000	25
MDPL1	0.313	16.630	1.148E-01	0.0677	0.1104	894.088	0.056	3200	31
MDPL1	0.333	15.960	1.175E-01	0.0670	0.1351	810.545	0.057	3000	33
MDPL1	0.500	12.150	1.306E-01	0.0640	0.1359	451.211	0.056	2000	50
MDPL1	0.625	10.340	1.385E-01	0.0636	0.1238	318.285	0.052	1600	62
A0.5	0.833	7.753	1.651E-01	0.0457	0.1518	180.093	0.052	600	4
MDPL1	1.000	7.140	1.730E-01	0.0441	0.1523	143.084	0.045	1000	100
MDPL1	1.250	5.902	1.983E-01	0.0274	0.1754	95.484	0.040	800	125
B0.5	1.250	5.712	2.059E-01	0.0050	0.1868	103.010	0.033	400	8
MDPL1	1.429	5.243	2.178E-01	0.0123	0.1946	75.355	0.036	700	143
A0.5	1.667	4.422	2.519E-01	-0.0232	0.2281	58.198	0.026	300	8
MDPL1	2.000	3.834	2.843E-01	-0.0185	0.2488	36.196	0.030	500	200
MDPL1	2.500	3.102	3.498E-01	-0.0399	0.2952	21.230	0.027	400	250
B0.5	2.500	3.048	3.587E-01	-0.0455	0.3021	20.465	0.026	200	16
A0.5	3.333	2.299	5.024E-01	-0.0545	0.3756	9.010	0.022	150	16
MDPL1	3.333	2.343	4.852E-01	-0.0520	0.3672	9.604	0.027	300	333
A1.5	5.000	1.555	9.813E-01	-0.0590	0.5324	2.829	0.020	300	8
MDPL1	5.000	1.586	8.902E-01	-0.0500	0.4950	3.100	0.020	200	500
A1.5	10.000	0.836	8.976E+00	-0.0423	1.0600	0.361	0.017	150	16
A2.5	10.000	0.832	9.127E+00	-0.0477	1.0971	0.368	0.013	250	8
A2.5	20.000	0.452	2.198E+03	-0.0340	2.5080	0.052	0.069	125	16

Table 3. Table 3: Parameters for the double-exponential model, eq.( 18 ), for redshifts z ≈ 0.5 − 1 𝑧 0.5 1 z\approx 0.5-1 . Different columns give: (1) Name of the simulation, (2) Cell size in h − 1 superscript ℎ 1 h^{-1} Mpc , (3) r ​ m ​ s 𝑟 𝑚 𝑠 rms density fluctuation, (4-7) parameters for eq.( 18 ), (8) r ​ m ​ s 𝑟 𝑚 𝑠 rms of the relative error of the fit, (9) number of cells in 1D, (10) relative force resolution – the number of force resolution elements per each density assignment cell, (11) redshift.

Simulation	Cell Size	$\sigma$	$A$	$\alpha +2$	${\rho }_0$	${\rho }_1$	relative error	$N_{\mathrm{cell}}$	$\Delta x/ϵ$	redshift
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)
MDPL1	0.167	16.050	1.448E-01	0.0146	0.0714	946.937	0.062	6000	16.7	0.492
MDPL1	0.250	12.780	1.392E-01	0.0346	0.0685	555.436	0.048	4000	25.0	0.492
MDPL1	0.500	8.177	1.683E-01	0.0218	0.1562	213.103	0.045	2000	50.0	0.492
MDPL1	1.000	4.793	2.476E-01	-0.0461	0.2313	72.823	0.031	1000	100.	0.492
MDPL1	1.250	3.961	2.909E-01	-0.0680	0.2668	46.150	0.028	800	125.	0.492
B0.5	1.250	3.837	2.993E-01	-0.0806	0.2764	45.377	0.026	400	8.	0.492
A0.5	1.667	2.923	3.882E-01	-0.0816	0.3300	19.242	0.038	300	8	0.5097
B0.5	2.500	2.068	6.078E-01	-0.0673	0.4228	6.282	0.042	200	16	0.492
MDPL1	2.500	2.112	5.801E-01	-0.0690	0.4130	6.892	0.044	400	250.	0.492
A0.5	3.333	1.575	1.015E+00	-0.0616	0.5398	2.635	0.038	150	16	0.5097
MDPL1	5.000	1.138	2.413E+00	-0.0538	0.7508	0.974	0.033	200	500.	0.492
A2.5	10.00	0.625	7.163E+01	-0.0526	1.6335	0.141	0.017	250	8	0.497
A2.5	20.00	0.347	1.409E+05	-0.0261	3.5731	0.024	0.140	125	16	0.497
B0.5	1.25	2.671	4.516E-01	-0.0597	0.3612	11.71	0.092	400	8	0.99
A0.5	1.667	2.066	6.462E-01	-0.0541	0.4374	5.372	0.090	300	8	0.99
B0.5	2.50	1.475	1.248E+00	-0.0501	0.5876	1.966	0.069	200	16	0.99
A0.5	3.333	1.157	2.484E+00	-0.0460	0.7563	0.941	0.056	150	16	0.99
A1.5	5.00	0.833	1.013E+01	-0.0413	1.1148	0.343	0.043	300	8	1.0
A1.5	10.00	0.487	1.134E+03	-0.0436	2.3475	0.061	0.032	150	16	1.0
A2.5	10.00	0.485	1.297E+03	-0.0403	2.4207	0.0597	0.0564	250	8	1.0
A2.5	20.00	0.274	3.522E+09	-0.0319	6.1154	0.0079	0.0782	125	16	1.0

Equations

$$\rho \equiv 1+{\delta }_{\mathrm{NL}}={\rho }_{\mathrm{DM}}/{\Omega }_m{\rho }_{\mathrm{cr}},$$

$$P(\rho )=\frac{\Delta N_{\mathrm{cell}}}{N_{cell}^3\Delta \rho }.$$

$${\int }_0^{\infty }P(\rho )d\rho =1,{\int }_0^{\infty }\rho P(\rho )d\rho =1.$$

$${\sigma }^2={\int }_0^{\infty }(\rho -1{)}^{2}P(\rho )d\rho =\frac{1}{2{\pi }^2}{\int }_0^{k_{Ny}}{P}_{NL}(k)W^2(k\Delta x)k^{2}dk,$$

$$\rho P_{LN}(\rho )=\frac{1}{\sqrt{2\pi {\sigma }_{LN}^2}}\exp \left(-\frac{[\ln (\rho )+{\sigma }_{LN}^2/2{]}^2}{2{\sigma }_{LN}^2}\right),$$

$${\sigma }_{LN}^2=\ln [1+{\sigma }^2].$$

$$P_{NBN}(\rho )=\frac{1}{\bar{N}}\cdot \frac{\Gamma (N+1/g)}{\Gamma (1/g)\Gamma (N+1)}\cdot \frac{(g\bar{N}{)}^N}{(1+g\bar{N}{)}^{N+1/g}},$$

$$P_{NBN}(\rho )\approx \frac{1}{{\bar{N}}^{1/g}\Gamma (1/g)}{\rho }^{1-1/g}\exp \left(-\frac{\rho }{g}\right).$$

$$\rho P_{GEV}(\rho )=\frac{1}{\ln (10)\beta }\frac{\exp (-z^{-1/k})}{z^{1+1/k}},z\equiv 1+\frac{k}{\beta }{\lg }_{10}\left(\frac{\rho }{{\rho }_0}\right),$$

$${\sigma }_L^2(M)=\frac{1}{2{\pi }^2}\int k^{2}dkP(k)W^2(kR).$$

$$\rho ={\left(1-\frac{{\delta }_L}{{\delta }_c}\right)}^{-{\delta }_c},\rho \equiv \frac{M}{\bar{M}},\bar{M}={\rho }_b\Delta x^3,$$

$${\rho }^{2}P(\rho )=\frac{1}{\sqrt{2\pi }}\exp \left(-\frac{{\delta }_L^2}{2{\sigma }_L^2}\right)\frac{d({\delta }_L/{\sigma }_L)}{d\ln \rho }.$$

$${\rho }^{2}P(\rho )=\frac{1}{\sqrt{2\pi {\sigma }_L^2}}\exp \left(-\frac{{\delta }_L^2}{2{\sigma }_L^2}\right)\left[1-{\delta }_L\left(\frac{1}{{\delta }_c}-\frac{\gamma }{3}\right)\right],$$

$$\gamma =-3\frac{d\ln {\sigma }_L^2}{d\ln M}.$$

$$P(\rho )d\rho =dV/V,$$

$$P(\rho )=-\frac{4\pi r^2(\rho )}{V}\frac{dr(\rho )}{d\rho }.$$

$${\rho }^{2}P(\rho )\propto {\rho }^{1-3/\alpha }.$$

$$P(\rho )=A{\rho }^{-\alpha }\exp \left[-{\left(\frac{{\rho }_0}{\rho }\right)}^{1.1}\right]\exp \left[-{\left(\frac{\rho }{{\rho }_1}\right)}^{0.55}\right],$$

$${\rho }_{\mathrm{one}}={\left(\frac{L}{\Delta x{N}_p}\right)}^3,$$

Full text

1 Introduction

The one-point probability distribution function (PDF) of the matter density field in the universe, and its related statistics the distribution of galaxy counts, have a long and somewhat patchy history in cosmology and extragalactic astronomy. It was Edwin Hubble almost a century ago who found that the counts of about $44,000$ extra-galactic nebulae distributed over a large area of the sky have a probability distribution that is not Gaussian but can be approximated by a log-normal distribution (Hubble, 1934). The statistics of galaxy counts in the Lick survey, in projected cells of size $10^{\prime}\times 10^{\prime}$, was studied by Soneira & Peebles (1978) who also discovered that the distribution of the counts is much broader than the Poisson PDF.

The $rms$ of galaxy counts $\sigma$ in cells of size $R$ is an integral over the power spectrum of the galaxy distribution (e.g., Peebles, 1980, Sec.36). As such, in former times, a count-in-cells analysis of the IRAS redshift galaxy survey was performed by Efstathiou et al. (1990) who used the counts as a measure of the two-point clustering statistics on different scales. Once methods to estimate the correlation function and the power spectrum were developed and new large-scale galaxy surveys were available, the count-in-cells as clustering statistics started to play a secondary role. Higher moments of cell counts depend on correlation functions of order larger than two. This means that the whole PDF has information not only on the two-point clustering but also on higher order statistics, which by itself is very valuable information.

At present, a precise description and modeling of the underlying matter density distribution - and biasing prescription that connects the dark matter field with the galaxy distribution - are fundamental to extract cosmological information from current and upcoming large-scale redshift and lensing galaxy surveys (e.g., Taruya et al., 2002; Takahashi et al., 2011; Carron et al., 2015; Clerkin et al., 2017; Manera et al., 2013; Kitaura et al., 2016). For this reason in the last years there has been a rejuvenated interest in the cosmic density distribution both from cosmological $N$-body simulations and galaxy surveys.

Wild et al. (2005) estimated the PDF of galaxies in the 2dF redshift survey using about 200,000 galaxies. Because of a relatively small volume, their analysis was done only for large cells of size 10-30 Mpc. They found that the log-normal distribution fits the data reasonably well, but the noise in the data did not allow them to make accurate measurements of the PDF. The situation was improved by Hurtado-Gil et al. (2017) using the count-in-cells statistics for galaxies in the SDSS main sample. They used $\sim 100,000$ galaxies and estimated the PDF for spheres of radius R=(8-24)\mbox{h^{-1}Mpc}. They found that the log-normal distribution was very inaccurate (a factor of $\sim 2$ errors) for spheres of R=8\mbox{h^{-1}Mpc}. A modification of the log-normal distribution (called log-normal+bias) somewhat improved the fits, but it still had $\sim 50\%$ errors at low number of galaxy counts. The negative binomial distribution was a much better fit for all filtering scales. At hight redshift, Bel et al. (2016) studied the count-in-cells distribution of $\sim 30,000$ galaxies in the VIPERS redshift survey with the typical average number of galaxies per cell of 0.5-5 and spherical cells of radius R=(4-8)\mbox{h^{-1}Mpc}. They found that the skewed log-normal distribution (a modification of the log-normal distribution with 4 more free parameters) was not accurate enough to fit the results of observations. Instead, they found that the negative binomial distribution was much more accurate. Yet, Clerkin et al. (2017) using the DES science verification data confirmed that the log-normal model is a good fit to both the galaxy density contrast and weak lensing convergence PDFs on scales of $(3-10)$ Mpc at median redshift $z=0.3$. In spite of the fact that at present these seems to be the best observational results, the errors and noise in the data are still substantial.

On the theoretical side the situation is also complicated. There are two types of approaches: models that start with some dynamical description of the non-linear evolution of the density field and proceed to make predictions of the matter PDF (e.g., Betancort-Rijo, 1991; Bernardeau, 1994; Kofman et al., 1994; Betancort-Rijo & López-Corredoira, 2002; Ohta et al., 2003; Lam & Sheth, 2008b); and then there are phenomenological approximations that assume a specific analytical form of the PDF (Coles & Jones, 1991; Gaztañaga et al., 2000; Lee et al., 2017; Shin et al., 2017), and find best-fitted parameters of this distribution function to simulation data.

Theoretical models based on the nonlinear dynamics typically use either some variant of the spherical infall model (e.g., Ohta et al., 2003; Lam & Sheth, 2008b; Neyrinck, 2016) or the Zeldovich approximation (e.g., Kofman et al., 1994; Betancort-Rijo & López-Corredoira, 2002). These models have made substantial progress and now can make very accurate predictions for relatively large smoothing scales R\gtrsim 5\mbox{h^{-1}Mpc} and small $\sigma\lesssim 1$ (Lam & Sheth, 2008b) giving errors less than $\sim 10\%$. Not surprisingly, as expected, the models are not very useful and start to fail at larger $rms$ fluctuations $\sigma\gtrsim 1$ (Lam & Sheth, 2008b; Neyrinck, 2016).

One of the disadvantages of the dynamical models is their complexity. They typically require some manipulation of the linear power spectrum, analytical approximations for different terms, and can be quite cumbersome to deal with. This is not a serious impediment to their use, but it is a nuisance. Simple analytical functions can serve as an alternative to more complicated dynamical models.

The log-normal distribution is an example of this approach. It was heavily advocated by Coles & Jones (1991) and is often used for relatively large smoothing scales. There is little justification why the density distribution function should be log-normal. Coles & Jones (1991) argue that under the assumption that the divergence of the peculiar velocity field in Eulerian coordinates grows as the velocities themselves (as given by linear theory) the density field can be expressed as the exponential of a Gaussian field. But while their assumption is acceptable for the Lagrangian divergence, for the Eulerian one there is an additional growth roughly proportional to the cubic root of the normalized density. This leads to a density field that is equal to a Gaussian field to the third power, whose PDF is quite different from a log-normal distribution. Here is the main argument of Coles & Jones (1991): “The lognormal is one of the simplest ways of defining a fully self-consistent random field which always has $\rho>0$ and, most importantly, is one of the few non-Gaussian random fields for which interesting properties are calculable analytically.” This says that the PDF should be log-normal because it can be handled analytically – hardly a serious argument. Another argument is of the same caliber: the log-normal distribution is well-known and frequently used in other fields of science (Ohta et al., 2003). Similar arguments were used for other phenomenological models (Gaztañaga et al., 2000; Lee et al., 2017).

The only real justification for the existing phenomenological approximations (including the log-normal) is that they make a fit to $N$-body results. This is the reason why cosmological $N$-body simulations are important for the field. In this paper we use a very large suite of cosmological simulations to produce accurate estimates of the dark matter distribution functions. Our simulations cover a wide range of numerical and cosmological parameters. We use the estimates to test different dynamical models and approximations for the PDF and to study its dependance on redshift and cosmological parameters.

One of our goals in this paper is to make a comprehensive study of the different effects that can be associated with $N$-body results regarding the matter density distribution, such as mass and force resolution, size of the box, shot noise and cosmic variance. In this regard we find that systematic errors in the PDF can be important. For example, noise related to the discreteness of the density probed by particles is the leading factor of seriously wrong estimates of the PDF in underdense regions.

Generation of mock galaxy catalogs provides a motivation for our study of the density distribution function. One needs to produce thousands of realizations of the dark matter density and velocity fields. This can be done by carefully tuning parameters of simulations and limiting their resolution to a fraction of a megaparsec (see e.g., Tassev et al., 2013; Chuang et al., 2015; Klypin & Prada, 2017). A biasing prescription then connects the dark matter with galaxies. This path requires knowledge of the distribution of dark matter mass on very small scales 100\mbox{h^{-1}kpc}-1\mbox{h^{-1}Mpc}. This is a challenge because the resolution of these simulations is not sufficient to resolve individual halos and subhalos making it difficult to apply existing tools such as Halo Abundance Matching and Halo Occupation Distribution. A path to solve the problem is to map dark matter to galaxies using a biasing scheme (e.g., Kitaura et al., 2016) which requires the understanding details of the density distribution function and finding limitations to its estimates.

Unfortunately, only very few studies in the literature provide PDF results for small smoothing scales \lesssim 1\,\mbox{h^{-1}Mpc} (Platen, 2009; Pandey et al., 2013; Lee et al., 2017). So, we will make an effort to study this regime too. Platen (2009) finds that in the regime of small smoothing and large density the PDF has a power-law shape with a slope of $\approx-1.9$, which is similar to what we find in this work.

This paper is organized as follows. In Sec. 2 we define quantities related with the PDF and provide details of some analytical approximations in Sec. 3.1. The spherical infall and the double-exponential models are introduced in Sec. 3.2. Numerical simulations used in this paper are discussed in Sec. 4. Sec. 5 also presents main features of the PDF. Accuracy of different approximations are discussed in Sec. 6. Summary of results is given in Sec. 7. Finally, numerical effects are discussed in Appendix A, and Appendix B presents tables of parameters for the double-exponential approximation.

2 Definitions

In order to estimate the density distribution function $P(\rho)$ from $N$-body simulations we split the computational volume $L^{3}$, where $L$ is the box size, with a $3D$ mesh of size $N^{3}_{\rm cell}$ and use the Cloud-In-Cell (CIC) density assignment scheme to estimate the density $\rho$ at each grid point of the mesh. The cell size of the grid $\Delta x=L/N_{\rm cell}$ defines the smoothing length. The density is normalized to the average matter density $\rho_{av}\equiv\Omega_{m}\rho_{\rm cr}$, i.e.

[TABLE]

where $\delta_{\rm NL}$ is the matter density contrast or overdensity. The index $NL$ highlights the fact that $\rho$ is a non-linear quantity and it can be distinguished from the density contrast $\delta$ as estimated by the linear theory. Throughout the paper we use the quantity $\rho$ as “density” in spite of the fact that it is really a normalized density – a dimensionless quantity as shown in eq. (1). This is done for convenience to avoid repeating $1+\delta_{\rm NL}$ in most plots and equations.

The values of density $\rho$ are binned using logarithmically spaced bins with width $\Delta\log_{10}(\rho)=0.025-0.050$. The density distribution function - PDF of the cosmic density field - is then defined as a normalized number of cells with density in the range $[\rho,\rho+\Delta\rho]$:

[TABLE]

The PDF can have a surprisingly large range of values. For example, density can reach values larger than $10^{5}$ for hundreds of cells when we use a large mesh of $\sim 3000^{3}$ cells in high-resolution simulations. That gives $P(\rho)\sim 10^{-12}$. At the same time the number of cells at low densities can be millions for a small density bin leading to a large PDF value $P(\rho)\gtrsim 1$. In order to avoid a large dynamical range of quantities, we typically plot $\rho^{2}P(\rho)$.

By design, the density distribution function is normalized to have the total volume and total mass equal to unity:

[TABLE]

The second moment of $P(\rho)$ gives the $rms$ density fluctuation of the density field $\sigma$, and as such it is related to the non-linear power spectrum $P_{NL}(k)$ of density perturbations:

[TABLE]

where $W^{2}(k\Delta x)$ is the power spectrum of the CIC filter with the width $\Delta x$, and the integral is truncated at the Nyquist frequency of the mesh $k_{Ny}=\pi/\Delta x$.

Figure 1 shows our first results on the structure of the PDF. More detailed discussion is given later in Sec. 5. The plot demonstrates the main trend of the shape of the PDF. In the linear regime of the growth of fluctuations, when $\sigma\ll 1$, the PDF is a Gaussian distribution that quickly acquires a skewed shape as $\sigma$ increases. With the further increase of $\sigma$ the PDF becomes progressively wider and becomes a power-law with the slope close to -2 that is smoothly suppressed on small and large densities.

3 Methodology

Different analytical approximations and theoretical models, as mentioned in Sec. 1, are used to fit and make predictions for the PDF estimates obtained from numerical simulations. Here we describe both approaches.

3.1 PDF analytical approximations: log-normal,

Negative Binomial, and Generalized Extreme Value

We introduce the log-normal, negative binomial and the generalized extreme value distributions which have been traditionally adopted as analytical approximations for the PDF in many works.

The log-normal (LN) distribution function $P_{LN}$ is defined as:

[TABLE]

where

[TABLE]

is the only free parameter. $\sigma^{2}_{LN}$ can be obtained from results of simulations, and thus should be considered as fixed.

Because the log-normal distribution does not provide accurate fits to numerical simulations, a number of modifications have been proposed (Hamilton, 1985; Colombi, 1994; Shin et al., 2017). None of those modifications extend the approximation to large densities, hence we do not discuss them in this paper.

The negative binomial (NBN) distribution (Betancort-Rijo, 2000; Gaztañaga et al., 2000; Bel et al., 2016) is defined as a discrete distribution. It is the probability $P_{N}(V)$ to find $N$ particles in a cell of volume $V$ with the average number of particles $\bar{N}$. It can be re-written as a distribution function of density contrast $\rho=N/\bar{N}$:

[TABLE]

where $g$ is a parameter that is defined by $rms$ fluctuations of counts: $g=(\sigma_{N}^{2}-\bar{N})/\bar{N}^{2}$. The average number of particles per cell $\bar{N}$ in general is not an integer number and is defined by the average density and cell volume: $\bar{N}=\bar{\rho}V$. So, the two parameters $g$ and $\bar{N}$ that define the NBN distribution are not free and can be fixed from simulations. However, for PDFs with not too small $rms$ fluctuations ($\sigma\gtrsim 1$) the fits produced by this distribution are not very accurate. As the result, we decided to treat both $\bar{N}$ and $g$ as free parameters.

While the negative bimodal distribution formally has two parameters, there is little change in $P_{NBN}$ when $\bar{N}>10$. For cases where NBN makes some reasonable fits, our simulations have typical values of $\bar{N}$ of many hundreds. So, in practice the NBN PDF depends only on one parameter $g$, which defines the width of the distribution function: the larger is $g$, the wider is $P(\rho)$.

At large average number of objects in a cell $\bar{N}\gg 1$ and for large $N>\bar{N}$ the NBN approximation predicts that the density distribution changes with density $\rho=N/\bar{N}$ as follows

[TABLE]

This expression is very different compared to the behavior of the PDF observed in $N$-body simulations for large $\rho$ and $\sigma$: in that regime $P(\rho)\propto\rho^{-2}\exp\left(-C\rho^{0.5}\right)$. Thus, the NBN approximation predicts too steep a decline with density and lacks the power-law regime of the $N$-body PDF.

A generalized extreme value (GEV) distribution as an approximation for the density distribution function was used by Lee et al. (2017). It can be written as:

[TABLE]

where $k,\rho_{0},$ and $\beta$ are free parameters.

3.2 Spherical Infall and double-exponential PDF models

Approximations discussed so far were not based on any dynamical models. They simply make a guess regarding the functional form of $P(\rho)$ and then proceed to finding parameters that produce the best fit. The guess is not based on any insights from the dynamics of clustering either. Spherical infall models are different because they are theoretical predictions for the density distribution function that are based on simplified approximations of the non-linear evolution of the density field. Here we closely follow the theoretical framework developed by Betancort-Rijo & López-Corredoira (2002) and Lam & Sheth (2008b). We assume that the linear density field was smoothed with a top-hat filter with radius $R_{f}$ and corresponding mass $M$. The variance of the smoothed field is equal to

[TABLE]

In the spherical infall model the mapping from linear density contrast $\delta_{L}$ to the nonlinear overdensity $\rho$ is approximated by the following relation (Betancort-Rijo, 1991; Bernardeau, 1994):

[TABLE]

where $\delta_{c}$ is the linear theory prediction for the critical overdensity of collapse. Here we will use $\delta_{c}=5/3$ as suggested by Betancort-Rijo & López-Corredoira (2002). For the initial Gaussian fluctuations this model gives the density distribution function:

[TABLE]

This is the same as eq.(6) in Lam & Sheth (2008b). We also use a modification of the spherical infall model, which is based on the excursion set model (Sheth, 1998; Lam & Sheth, 2008b):

[TABLE]

where $\gamma(\rho)$ is:

[TABLE]

Note that in these equations $\sigma_{L}$ is a function of filtering mass $M$, which in turn depends on density $\rho=M/\bar{M}$. So, $\sigma_{L}=\sigma_{L}(\rho)$.

In order to apply the models, we need to adjust the top-hat filtering scale $R_{f}$ used in the spherical infall model so that it matches the Cloud-In-Cell filtering scale used in our simulations. This is done by matching the power spectra of both filters at wavenumbers $k<0.7k_{Ny}$ by applying $R_{f}=\Delta x/\sqrt{1.3}$, where $k_{Ny}=\pi/\Delta x$ is the Nyquist frequency of the density grid used for density assignment, and $\Delta x$ is the grid size. Specifically, for every bin with density $\rho$ we find the mass $M=\rho\bar{M}$ and then the top-hat filtering radius is found as $R_{f}=\sqrt{2}\left[3M/4\pi\rho_{cr}\Omega_{m}\right]^{1/3}$. When applying the relations given in eqs.(10–14), we integrate numerically eq.(4) with the top-hat filter $W(kR_{f})$ and use finite differences to estimate the derivatives in eq.(12) and eq.(14).

As was found previously (e.g., Betancort-Rijo & López-Corredoira, 2002; Lam & Sheth, 2008b, a; Neyrinck, 2016), the spherical infall model provides good approximations for the PDF in those cases with large filtering scales where the $rms$ fluctuation in the simulation box $\sigma(M)$ is relatively small $\sigma\lesssim 1$. This is consistent with our results, which will be presented later. At larger $\sigma$ (large densities) the model provides results that are much less accurate. Even so, the spherical infall model predicts trends that track our $N$-body results. This is somewhat unexpected because at large overdensities $\rho\gtrsim 10^{2}$ the model is clearly outside of the limits of its dynamical applicability. After all, a simplistic treatment of non-linear evolution used by the model cannot be valid for densities that are appropriate for collapsed and virialized halos.

It is interesting and instructive to find what makes the spherical infall model much better than expected. The last factor in eq. (12) is just a correction to the leading terms that are the power-law $P\propto\rho^{-2}$ (see also Figure 1) and the exponential term on the right-hand-side of eq. (12). The exponential term originates from the assumption that the density distribution function of primordial fluctuations is Gaussian. The $\rho^{-2}$ term comes from integrating the Gaussian distribution function over mass and then by writing it in differential form. This is basically the same logic as in the Press-Schechter derivation of the mass function of dark matter halos.

More in detail, we see that the exponential term provides the truncation of the $P\propto\rho^{-2}$ behavior both on the low-density $\rho<1$ and on the high-density $\rho\gg 1$ regimes. At low densities (small masses, large $\sigma_{L}$) the truncation is mostly due to large negative values of $\delta_{L}$. At large $\rho$ (large masses and small $\sigma_{L}$) the decline is related to the combination of decreasing $\sigma_{L}$ and increasing $\delta_{L}$. However, the increase of $\delta_{L}$ is limited: it can not exceed $\delta_{L,{\rm max}}=\delta_{c}=5/3$, and the decline in $\sigma_{L}$ is not strong enough by itself to produce a substantial suppression of the PDF. So, the $\rho^{2}P(\rho)$ shape becomes nearly flat at very large densities $\rho>100$ and small smoothing scales \Delta x\lesssim 1\mbox{h^{-1}Mpc}. This is clearly seen in Figure 2.

At very small smoothing scales \Delta x\ll 1\mbox{h^{-1}Mpc} and large densities there is another regime, that sheds light on how the density distribution function should behave at very large densities. When density is larger than $\sim 100$ we are likely dealing with interiors of collapsed dark matter halos. In this regime the distribution function $P(\rho)$ is a sum over the distribution functions of individual halos. Assuming that the density profile of a halo can be approximated by the Navarro-Frank-White (NFW) profile, and for a very small filtering scale, we can can derive the PDF provided by a single halo. If $\rho(r)$ is the halo density profile, then the density distribution function given in eq. (2) can be written as

[TABLE]

where $dV/V$ is a fraction of volume with density in the range $(\rho,\rho+d\rho)$. The density dependance on radius can be inverted to give us the radius at a given density $r=r(\rho)$. Then the PDF in eq.(15) takes the form

[TABLE]

This can be applied, for example, to the NFW profile. It is easy to see the trend, if the density profile is a power-law $\rho\propto r^{-\alpha}$ with the slope $\alpha$. In this case

[TABLE]

In the outer regions of a dark matter halo the density declines as a power-law with the slope $\alpha\approx 3$, which implies that $\rho^{2}P(\rho)\approx{\it constant}$. It is easy to invert the NFW profile numerically. Results show that for the halo masses in the range M=(10^{12}-10^{15})\mbox{h^{-1}M_{\odot} } the product $\rho^{2}P(\rho)$ is nearly constant for overdensities $\rho=10^{2}-10^{4}$ and declines at much larger densities. This decline is consistent with the fact that the slope $\alpha$ becomes smaller at radii comparable with the characteristic scale radius $r_{s}$ of the NFW profile. Adding results for many halos with different masses and concentrations will change the behavior of the $\rho^{2}P(\rho)$ trend at $\rho>10^{4}$, but not for the range $\rho=10^{2}-10^{4}$, where it should remain flat because it is also flat for each halo.

In summary, both the spherical infall model and the dark matter halo profiles indicate that the leading term in the density distribution function should be $P(\rho)\propto\rho^{-2}$. The trend should be modified by adding suppression on large and small scales. However insightful, the spherical infall model or the NFW results at large densities cannot be used to produce accurate results for the density distribution function. We use those hints to construct our own approximation for the PDF $P(\rho)$.

Motivated by these results and by our simulations we design our own model. It is nearly a power-law $P(\rho)\propto\rho^{-\alpha}$ with the slope $\alpha\approx 2$ that is truncated with exponents on both the small and large densities. We call this model double-exponential. We tested different shapes for the exponential terms and find that the following expression provides errors less than few percent at all redshifts, smoothing scales, and cosmologies that we study:

[TABLE]

where $A,\,\alpha$,  $\rho_{0}$, and $\rho_{1}$ are free parameters. As noticed above the slope $\alpha\approx 2$. The slopes in the exponential terms 0.55 and 1.1 are results of the fitting of numerical PDFs at different smoothing scales and redshifts. One may expect that adding two more free parameters to the approximation (i.e., the slopes in the exponential terms) may further improve the quality of the fits. We find that this is not the case: the data prefer the same slopes regardless of the value of $\sigma$.

The double-exponential model has four formal free parameters. One may use three constraints to limit the parameters: the total mass and volume should be equal to unity (see eq. (3)), and the second moment of the PDF should be equal to $\sigma$ measured in simulations (see eq. (4)). Note that there must be a degree of freedom left after fixing the constraints otherwise the model would not be able to reproduce the numerical results that show that the PDF is not defined solely by $\sigma$, and depends on both the redshift and $\Omega_{m}$. The double-exponential model has this additional degree of freedom.

In practice, we use all four parameters to fit the numerical data. We typically find that the best-fit parameters provide PDF approximations that within 1-2% conserve the mass and match well the numerical value of $\sigma$ measured in the simulations. The volume is conserved within 1-5% accuracy.

4 Simulations

Numerical parameters of our simulations are presented in Table 1, which gives box size, number of particles, mass of a particle $m_{p}$, number of mesh points $N_{g}^{3}$ (if relevant), cell size of the density/force mesh $\epsilon$, the number of time-steps $N_{s}$, cosmological parameters $\sigma_{8}$ and $\Omega_{m}$, and number of realizations $N_{r}$.

In order to estimate the density distribution function, we split each simulation box with a $3-D$ mesh with the cell size $\Delta x$. We then use the Cloud-In-Cell (CIC) density assignment to generate the density field. Many filtering sizes $\Delta x$ were used for each simulation and snapshot.

Different codes were used to make those simulations. The MultiDark $1~{}Gpc/h$ simulation (MDPL1) (Klypin et al., 2016) was done with the GADGET-2 code (Springel, 2005). The ART code (Kravtsov et al., 1997) was used to produce the BolshoiP simulation (Klypin et al., 2011). These two simulations have the largest resolution and the largest number-density of particles. However, there are only two of these simulations because they are very expensive computationally. Other simulations were carried out with the parallel Particle-Mesh code GLAM (Klypin & Prada, 2017). Because the GLAM code is much faster, we have many realizations of the same cosmological and numerical parameters. Simulations B0.5 and D0.5 were designed for testing possible dependance of the PDF on the matter density $\Omega_{m}$. These simulations have the same random seeds to make comparisons of results easier.

All the simulations were started at initial redshift $z_{\rm init}=100$ using the Zeldovich approximation. The simulations span three orders of magnitude in mass resolution, a factor of hundred in force resolution and differ by a factor of 500 in effective volume (see Table 1). Altogether, we use about 7000 simulations to study the density distribution function. To our knowledge this is the largest set of simulations available today for this type of analysis.

Analysis of different numerical effects in the estimates of the PDF is presented in Appendix A. In summary, our results are mostly dominated by systematics, not by the finite-volume simulation variance. When dealing with individual simulations such as MDPL1 or BolshoiP we use only bins with more than $N>100$ cells per bin. For the large sets of simulations A0.5, A1.5, A2.5 we accept bins with more than 10 cells. The discreteness of density assignment may become an issue at small densities while the force resolution may affect the high-density tail of the PDF. We find limits on numerical parameters that should be satisfied to produce the $P(\rho)$ with errors less than few percent: (1) the filtering scale $\Delta x$ must be resolved with not less than 8 force resolution elements: $\Delta x>8\epsilon$, and (2) the number of particles per filtering cell should be not less than 10-20.

5 PDF: main trends

The overall dependance of $P(\rho)$ on filtering scale and $rms$ density fluctuation $\sigma$ is illustrated in Figures 1 and 2. Full curves in the plots show the results from our simulations. At small $\sigma$ and large $\Delta x$ the PDF has a peak at $\rho\approx 1$ that shifts to smaller densities as $\sigma$ increases (and smoothing scale $\Delta x$ decreases). At the same time the distribution function becomes extremely wide and develops a distinct $\rho^{-2}$ power-law trend that expands both to large and small densities. By scaling out the $\rho^{-2}$ dependance (Figures 2) we reduce the dynamical range and can see better details of the PDF. In particular, we find a very steep decline of $\rho^{2}P(\rho)$ on both high and low density tails. The double-exponential model eq.(18) was tuned to find the shape of both declines. Our results for different filtering scales and different redshifts show that the decline at the large-density limit is $\propto\exp[-(\rho/\rho_{1})^{\nu}]$ and at the small-density limit it is $\propto\exp[-(\rho/\rho_{0})^{-2\nu}]$ with $\nu\approx 0.55$.

Note that at very large $\sigma\gtrsim 10$ the peak of $\rho^{2}P(\rho)$ has a nearly constant amplitude $\rho^{2}P(\rho)\approx 0.11$ but the position of the peak shifts to larger values of $\rho$. Because the total mass must be preserved ($\int\rho P(\rho)d\rho=1$), this implies that at intermediate scales $\rho\approx 1-100$ the PDF $P(\rho)$ should decline when $\sigma$ increases, and that the slope $\alpha$ in eq.(18) must become slightly shallower with increasing $\sigma$.

It is often taken for granted that the PDF depends only on the amplitude of the density perturbations on a given filtering scale $\sigma$. Indeed, this is the dominant behavior of $P(\rho)$. However, this is not exactly correct. Our simulations have such a good accuracy that now we can test the dependance of the PDF on redshift at fixed $\sigma$ and on cosmological parameters.

We first select redshifts and smoothing scales in such a way that $\sigma$ for two different redshifts are nearly identical. Figure 3 presents two examples of such cases – one for relatively low $\sigma\approx 1$ and another for larger $\sigma\approx 4$. The differences between PDFs at the same $\sigma$ and different $z$ are not large: 5-10% depending on the density where the differences are measured. Nevertheless, the differences clearly exist. We also estimate the differences using the spherical infall model and present the results in the right panel of Figure 3.

The PDF also slightly depends on parameters of the cosmological model. In the left panel of Figure 4 we compare $z=0$ results for the B0.5Gpc and D0.5Gpc simulations that differ only by the matter density parameter $\Omega_{m}$. Again, differences are small but clearly present at $\sim 5-10\%$ level. The right panel shows predictions for the spherical infall model with the same basic conclusion: PDF does depend on $\Omega_{m}$.

One can understand why the PDF depends on $z$ and $\Omega_{m}$, if one realizes that at any given density $\rho$ the value of $P(\rho)$ is formally a functional on the non-linear power spectrum. This means that it depends not only on $\sigma(R)$ but also on the whole shape of the power spectrum. We can analyze the situation by assuming that $P(\rho)$ depends just on $\sigma_{L}(R)$ (see eq. (10)) and on its local logarithmic derivative $\gamma$ at scale $R$ as defined by eq. (14). Then the dependence of the PDF on redshift for a fixed value of $\sigma_{L}(R)$ can be explained, because at larger redshits the given value of $\sigma_{L}$ is attained at a smaller scale $R$, where $\sigma_{L}(R)$ is less steep for a CDM-type power spectrum.

If we look at different terms for $\rho^{2}P(\rho)$ in eq. (12), then we note that the $\sigma_{L}$ in the argument of the exponential term is also a function of $\rho$ (falling faster with $\rho$ as $\gamma$ grows). So, for larger values of $\gamma$ the PDF will be smaller for large $\rho$ values, where the exponential behavior dominates, while for small values of $\rho$ it will be larger. The derivative in the right-hand-side of eq. (12) also depends on $\gamma$: it will be larger for larger $\gamma$ values. However, it is only important for the intermediate values of density, where the behavior can qualitatively be explained by the behavior in the extremes and the conservation of probability. In short, this all implies a smaller PDF in the low $\rho$ limit and a larger one in the large $\rho$ limit as compared to the PDF at smaller redshift.

The same qualitative behavior would be observed for smaller values of $\Omega_{m}$ when studying the dependence of the PDF on $\Omega_{m}$ at a given redshift: as $\Omega_{m}$ decreases the power spectrum flattens (up to scales of the order of the horizon at the start of matter domination), leading to a less steep $\sigma_{L}.$

6 PDF: Testing different approximations and models.

We start our analysis of different approximations by testing the spherical infall model and the log-normal distribution. Both models are expected to work and typically used for relatively low rms fluctuations $\sigma\lesssim 1$.

Figure 5 shows results for two modifications of the spherical infall model. We select five configurations with different filtering scales, $\sigma$ and redshifts. There are some differences between the pure spherical infall model in eq. (12) and the excursion set model in eq. (13). For example, the excursion set model produces smaller errors at $\rho<1$ and $\sigma\approx 1$. At the same time, it makes visibly larger errors at $\rho>1$. For this reason we prefer the standard spherical infall model. It provides smaller than 10% errors for points that are larger than 0.1 of the maximum of the PDF. The error increases substantially in the peripheral regions. If better accuracy is required for the small values of the PDF, one would need to use results of $N$-body simulations, not the approximations.

Results for the log-normal distribution are shown in Figure 6, where in the left panel we present $z=0$ results, and results for different redshifts are shown in the right panel. It is clear that the log-normal distribution produces much worse fits as compared with the spherical infall model. For example, errors in the central region are less than 10% only for $\sigma<0.5$. They become dramatically worse for even slightly larger $\sigma$. For $\sigma>1$ the lognormal distribution makes typically $\sim 50\%$ errors and, which is even worse, predicts a wrong shape of the PDF. It predicts a wrong position of the maximum; slopes of the declining PDF are not correct.

The log-normal distribution has two advantages as compared with the spherical infall model: (a) It is simple; does not require the machinery of handling the power spectrum and numerical derivatives. (b) it makes mediocre predictions that do not totally fail. We definitely recommend it for quick-and-dirty applications, but not for accurate estimates.

Left panels in Figure 7 show results for the Negative Binomial distribution given by eq. (7). For the large smoothing scale \Delta x=20\mbox{h^{-1}Mpc} it provides a reasonable fit with errors $\sim 10\%$ for densities $\rho=0.4-2.5$. However, it is clear that it has a wrong shape: too steep at large densities and too shallow at small densities. This becomes a serious issue for large values of $\sigma$. For example, we could not find any good fit for \Delta x=5\mbox{h^{-1}Mpc} shown in the bottom panel. It may be not fare to use the NBN for the dark matter PDF because we are in a regime that was not favorable for the NBN: study configurations with very large number of particles per cell while the NBN was designed to handle very small number $\bar{N}$.

The GEV approximation eq. (9) scored much better, as illustrated in the right panels of Figure 7. Indeed, it provided excellent fits for $\sigma\lesssim 2$ with errors less than 10% and even for larger $\sigma$ it gives very good accuracy, but it starts to fail catastrophically at very large densities. Lee et al. (2017) tested this approximation for PDF for the Bolshoi simulation. Their results are compatible with ours. However, Lee et al. (2017) seems did not pay attention to the situation at large densities, where the GEV becomes unacceptable. It still can be useful for low densities, but the problem is that there is no a priory estimate at what density and $\sigma$ the GEV fails.

Comparison of the double-exponential model with the $N$-body results has been already presented in Figure 1. More detailed analysis is shown in Figure 8. Parameters of the fits for different smoothing scales and redshifts are given in Appendix B. Of all approximations studied in this paper the double exponential model is by far the best. It provides very accurate (few percent) fits for $\sigma\gtrsim 1$ with densities $\rho\approx 10^{-1}-10^{5}$. With somewhat larger errors it still works down to $\sigma\approx 0.5$. To some degree the success of the approximation is not surprising because it was designed to reproduce the main features of the spherical infall model and the PDF of individual NFW halo profiles that predict the $P(\rho)\propto\rho^{-2}$ trend for large densities. And, of course, the parameters of the approximation – the two power-law slopes – were tuned to produce best fits.

However, it is unexpected that the exponential terms in eq. (18) should have the same power-law slopes 0.55 and 1.11 for all densities and smoothing scales. Indeed, we tried different combinations of the slopes – even with changing values of the slopes for different $\sigma$ – and did not find them to improve the fits. There is one problem with the constant slopes, though: the approximation cannot work for a very small $\sigma$ where the PDF must become a Gaussian. The eq. (18) does not allow a transition to a Gaussian distribution. This is not a serious issue because at $\sigma\lesssim 1$ the spherical infall model provides an adequate approximation for the dark matter density distribution function.

7 Conclusions

Using a large suite of cosmological $N$-body simulations we study the shape and the evolution of the dark matter one-point probability distribution function. Unlike most of other studies, we cover a very large range of smoothing scales R=100\,\mbox{h^{-1}kpc}-20\,\mbox{h^{-1}Mpc} and $rms$ density fluctuations $\sigma$. We find that as $\sigma$ increases, the PDF becomes a power-law $P\propto\rho^{-2}$ that is truncated with exponential terms on both small and large densities. This trend is consistent with the extrapolation of both the spherical infall model and the PDF expected at high densities for the NFW density profile of dark matter halos.

The PDF weakly depends on redshift (at fixed $\sigma$) and on matter density $\Omega_{m}$. The effect is relatively small ($\sim 5-10\%$) but is clearly observed in simulations. The spherical infall model also has the same trend. This behavior contradicts analytical approximations such as the log-normal or the GEV that assume that the PDF should depend only on $\sigma$.

The basic trend $P\propto\rho^{-2}$ gives us a motivation to construct a new model given in eq. (18), which we call double-exponential distribution. It formally has four free parameters of which two can be fixed by requiring that the total volume and mass must be equal to unity. The model works only for $\sigma\gtrsim 1$ and does not allow the transition to a Gaussian distribution as expected for $\sigma\ll 1$. Nevertheless, for $\sigma\gtrsim 1$ the model gives the best performance of all approximations that we tested in this work, with errors of just few percent for $P(\rho)$ when the PDF changes by 12 orders of magnitude. Parameters of the double-exponential model are provided in Appendix B. The model potentially may be modified by adding few extra parameters to allow for accurate treatment at small $\sigma$. We did not try to do it for two reasons: (1) the spherical infall model gives accurate enough treatment for this regime, and (2) it is cheap to make adequately accurate $N$-body simulations for $\sigma<1$ if needed.

The spherical infall model provides accurate predictions for the $N$-body PDF results for low values of $\sigma<1$, but it becomes less reliable for larger $\sigma$, which is expected for this model. The combination of the double-exponential model at large $rms$ fluctuation with the spherical infall model in the small $rms$ regime yields a remarkably accurate density distribution for all regimes of clustering of the cosmological matter field.

We also tested different analytical approximations. The often used log-normal distribution, as was shown before, does not make good fits and needs significant modifications before it can provide accurate results. It is pretty much useless for large densities because it does not provide a path to explain the main trend at large densities, i.e. the power-law trend $P\propto\rho^{-2}$. It clearly has an advantage of being simple, and it does not fail catastrophically as shown in Figure 6. Our results for $\sigma<1$ demonstrate that the log-normal approximation always made significantly worse fits as compared with the spherical infall model predictions. In addition, the log-normal approximation cannot accomodate the dependance of the PDF on redshift and $\Omega_{m}$, while the spherical infall model nicely reproduces the effect.

The GEV approximation scores much better than the log-normal distribution. Even for a large $\sigma\approx 10$ it gives remarkably accurate results for densities up to $\rho\approx 10^{3}$. However, the approximation fails catastrophically at larger densities, and there is no obvious way to predict at what density it still works or fails.

While it is useful and insightful to have analytical models for the PDF, one does not really need them if $\sim 1$% accuracy is a requirement and $\sigma\gtrsim 0.5$. Then there is no alternative to $N$-body simulations: one can get very accurate, fast and cheap results. The cost of one GLAM A0.5 simulation is just 3.5 hrs of wall clock on a modest data server (Klypin & Prada, 2017). For the B0.5 run it is 10 hrs. In this paper we analyse thousands of realizations, but this was to prove that one does not need those to make accurate PDF: just a few realizations is enough. However, one needs to be carefull about numerical effects when using these fast Particle-Mesh simulations. Klypin & Prada (2017) provide detailed description of constraints for the simulations, and Appendix A in this paper gives prescriptions on how to use the $N$-body simulations to reach $\sim 1\%$ accuracy in the PDF.

Acknowledgements

A.K. acknowledges support of the Fulbright Foundation and support of the Instituto de Fisica Teorica, CSIC, Madrid, Spain. F.P. acknowledges support from the Spanish MINECO grant AYA2014-60641-C2-1-P. FDA acknowledges financial support from “la Caixa”-Severo Ochoa doctoral fellowship. We thank Johan Comparat (IFT, Madrid), Claudia Scoccola (IAC, Tenerife), and Sergio Rodriguez-Torres (UAM, Madrid) for comments and fruitful discussions. The PPM-GLAM simulations have been performed on the FinisTarrae II supercomputer at CESGA in Galicia, supported by the Xunta de Galicia, CISIC, MINECO and EU-ERDF.

Appendix A Numerical effects: finite-volume variance,

mass and force resolution

Numerical effects may play an important role for the estimates of the density distribution function. We start with the analysis of the effects of the variance due to the fine-volume of the simulations. We use the GLAM simulations A0.5 with small cell \Delta x=0.83\mbox{h^{-1}Mpc} and A1.5 with larger cell \Delta x=0.5\mbox{h^{-1}Mpc} to estimate the level of fluctuations in different realizations. Bottom panels in Figure 9 show the average values of $\rho^{2}P(\rho)$ and the statistical fluctuations of a single realization. As expected, the fluctuations due to the fine-volume simulation variance are larger for very large densities and become very small for $\rho\approx 0.5-10$.

The middle panels present relative fluctuations $\Delta P/\langle P\rangle$, where $\langle P\rangle$ is the average PDF over an ensamble of realizations and $\Delta P$ is the $rms$ deviation. The vertical dotted lines at large $\rho$ show the density bin with 100 cells in a single realization. At this density the level of statistical fluctuations $\Delta P/P$ is about 0.1, which is consistant with the expected shot noise. We clarify this situation by plotting in the top panels the relative fluctuations scaled with $N^{1/2}$, where $N$ is the number of cells of given density in a bin. Indeed, the fluctuations are defined by the number of cells in a bin for large densities.

The situation is different at small densities $\rho\lesssim 1$ where the fluctuations become substantially stronger than the Gaussian $\Delta P/P=N^{-1/2}$. This is likely related with the increasing noise in the density field due to too few particles per cell $n$. The vertical dotted lines on low densities show the bin at which $n=10$.

In spite of being strongly non-Poissonian at small densities, the errors are still very small. For example, at $\rho=0.1$ for the A0.5 simulations the errors are just $\sim 2$% for a single realization (right panels). Note that the fluctuations plotted in Figure 9 provide average deviations of a single realization from the ensamble average. Errors of the average $P(\rho)$ are significantly smaller. For example, for the A1.5 simulations (left panels) the error of the mean at $\rho=60$, $N=100$ is just $\sim 0.1$%.

In summary, our results are mostly dominated by systematics, not by the variance. When dealing with individual simulations such as MDPL1 or BolshoiP we use only bins with more than $N>100$ per bin. For large sets of simulations A0.5, 1.5, 2.5 we accept bins with more than 10 cells.

In order to evaluate other possible numerical effects, we select two filtering scales \Delta x=1.25\mbox{h^{-1}Mpc} and \Delta x=5\mbox{h^{-1}Mpc} and analyze $P(\rho)$ at $z=0$ obtained from different simulations. The two filtering scales probe different dynamical regimes. For \Delta x=1.25\mbox{h^{-1}Mpc} the $rms$ density fluctuation is large $\sigma\approx 5$. So, we are testing very nonlinear regime with densities up to 1000. The larger filtering scale \Delta x=5\mbox{h^{-1}Mpc} probes more modest fluctuations with $\sigma\approx 1.5$ and densities $\rho=0.1-100$.

By comparing these results we test the effects of the finite box size (ranging from 500\mbox{h^{-1}Mpc} to 2500\mbox{h^{-1}Mpc}), the effects of force resolution $\Delta x/\epsilon$ (ranging from 1 to 500) and the discreteness effects associated with the finite number of dark matter particles. The later can be characterized by the density $\rho_{\rm one}$ produced by a single particle placed at the node of a cell:

[TABLE]

where $N_{p}^{3}$ is the number of particles and $L$ is the size of the simulation box. Vertical lines in the left panel of Figure 10 show $\rho_{\rm one}$ for different simulations. Values of different parameters are also given in Figure 10.

We use the left panel in Figure 10 to demonstrate two numerical effects. At large densities $\rho\gtrsim 50$ the discreteness effects are small even for the A2.5 simulations and the results are dominated by the force resolution. Here we find a trend that is expected for simulations with low force resolution: the PDF increases with increasing force resolution. However, there is little difference between simulations when the resolution becomes sufficiently high: the PDF for B0.5 simulations with 8 resolution elements across one cell( $\Delta x/\epsilon=8$) is nearly the same as for the much high resolution simulation MDPL1 with $\Delta x/\epsilon=125$. This is a signature that the results have converged.

The discreteness of density assignment becomes an issue at small densities as is clearly seen in Figure 10 for $\rho\lesssim 10$. The A2.5 simulations provide a good example on how the particle noise affects the PDF. There is a large bump in $\rho^{2}P(\rho)$ at densities slightly below $\rho_{\rm one}$. At much smaller densities $\rho\approx 0.1\rho_{\rm one}$ the PDF falls much below the real one. The effects of the particle noise extend to densities above $\rho_{\rm one}$ but quickly die out beyond $\rho\approx 10\rho_{\rm one}$.

The right panel in Figure 10 shows better convergence of $P(\rho)$ because we select simulations that have better force resolution and smaller level of the particle noise. In order to see the differences more clearly we make a fit to the data, and on the top panel plot we show only the deviations from the fit. The A2.5 data still fall below the more accurate MDPL1 results both at small densities (effects of the particle noise) and large densities (effects of the force resolution). The A1.5 data show much smaller errors.

Using the results presented in Figure 10 and similar results of the comparison between different simulations (e.g., comparison of BolshoiP and MDPL1) we find limits on numerical parameters that should be satisfied to produce a PDF $P(\rho)$ with errors less than a few percent:

• •

the filtering scale $\Delta x$ must be resolved with not less than 8 force resolution elements: $\Delta x>8\epsilon$

• •

the number of particles per filtering cell should be not less than 10-20: $\rho>(10-20)\rho_{\rm one}$.

Appendix B Parameters of the double-exponential model for the Dark Matter Density Distribution Function

Tables 2 and 3 give parameters for the double-exponential model given in eq.(18) for different simulations, filtering scales and redshifts.