New fitting formula for cosmic non-linear density distribution

TL;DR

This paper introduces a new fitting formula, Nv2, for the cosmic matter density PDF, which outperforms the traditional log-normal model, especially in underdense regions, and captures non-Gaussian features across different scales and redshifts.

Contribution

The paper proposes the Nv2 distribution as a superior fit for cosmic density PDFs, with detailed analysis of its parameters' evolution across scales and cosmologies.

Findings

Nv2 fits cosmic density PDFs better than log-normal.

Significant improvement in underdense regions.

Parameters evolve with redshift and cosmology.

Abstract

We have measured the probability distribution function (PDF) of cosmic matter density field from a suite of N-body simulations. We propose the generalized normal distribution of version 2 (Nv2) as an alternative fitting formula to the well-known log-normal distribution. We find that Nv2 provides significantly better fit than the log-normal distribution for all smoothing radii (2, 5, 10, 25 [Mpc/h]) that we studied. The improvement is substantial in the underdense regions. The development of non- Gaissianities in the cosmic matter density field is captured by continuous evolution of the skewness and shifts parameters of the Nv2 distribution. We present the redshift evolution of these parameters for aforementioned smoothing radii and various background cosmology models. All the PDFs measured from large and high-resolution N-body simulations that we use in this study can be obtained from a Web site at https://astro.kias.re.kr/jhshin.

Tables

Table 1. Table 1 Numerical coefficients of Fitted functions for α 𝛼 \alpha , κ 𝜅 \kappa , and ξ 𝜉 \xi

Model	Numerical coefficients	$R^2$ goodness
${\Lambda }_{m0}$	$a_1=32.92$	0.992
	$b_1=10.38,b_2=2.118,b_3=3.151$	0.994
	$c_1=-4.082,c_2=44.92$	0.969
${\Lambda }_{m-}$	$a_1=35.48$	0.992
	$b_1=11.85,b_2=3.132,b_3=3.683$	0.993
	$c_1=-4.370,c_2=54.48$	0.985
${\Lambda }_{m+}$	$a_1=30.82$	0.992
	$b_1=9.369,b_2=1.426,b_3=2.787$	0.993
	$c_1=-3.901,c_2=39.85$	0.972
$Q_{w-}$	$a_1=31.54$	0.992
	$b_1=10.22,b_2=2.920,b_3=2.124$	0.996
	$c_1=-3.929,c_2=41.41$	0.980
$Q_{w+}$	$a_1=32.05$	0.994
	$b_1=10.30,b_2=2.890,b_3=2.250$	0.996
	$c_1=-4.024,c_2=45.35$	0.980

Equations

$${\mathcal{N}}_{v2}({\nu }_{\mathrm{sim}})=\frac{\mathcal{N}(y_{\mathrm{sim}})}{\alpha +\kappa ({\nu }_{\mathrm{sim}}-\xi )},$$

$$y_{\mathrm{sim}}=\{\begin{array}{ll}\frac{1}{\kappa }\ln \left[1+\frac{\kappa ({\nu }_{\mathrm{sim}}-\xi )}{\alpha }\right]& \text{if κ=0}\\ \frac{{\nu }_{\mathrm{sim}}-\xi }{\alpha }& \text{otherwise}\end{array},$$

$$P_{\mathrm{LN}}({\delta }_{\mathrm{sim}})=$$

$$\times \exp \left\{\frac{[\ln (1+{\delta }_{\mathrm{sim}})+{\sigma }_1^2/2{]}^2}{2{\sigma }_1^2}\right\},$$

$${\sigma }_1^2=\ln [1+{\sigma }_{\mathrm{sim}}^2].$$

$$P_{\mathrm{SLN}}({\upsilon }_{\mathrm{sim}})=$$

$$+\frac{10}{6!}{T}_3^2{\sigma }_{\Phi }^2{H}_6({\upsilon }_{\mathrm{sim}})]\mathcal{N}({\upsilon }_{\mathrm{sim}}),$$

$$T_3\equiv \frac{⟨{\Phi }^3⟩}{{\sigma }_{\Phi }^4},T_4\equiv \frac{⟨{\Phi }^4⟩-3{\sigma }_{\Phi }^4}{{\sigma }_{\Phi }^6}.$$

$${\alpha }_{\mathrm{fit}}(R_{\mathrm{th}},D_1)$$

$${\kappa }_{\mathrm{fit}}(R_{\mathrm{th}},D_1)$$

$${\xi }_{\mathrm{fit}}(R_{\mathrm{th}},D_1)$$

Full text

1. Introduction

The inflationary models (Starobinskiǐ, 1979; Starobinsky, 1982; Guth, 1981; Sato, 1981; Linde, 1982; Albrecht & Steinhardt, 1982) of the early Universe predict that the primordial density perturbations generated during inflation (Mukhanov & Chibisov, 1981; Hawking, 1982; Guth & Pi, 1982; Bardeen et al., 1983) must obey nearly Gaussian statistics (Maldacena, 2003; Acquaviva et al., 2003; Creminelli & Zaldarriaga, 2004). This prediction is confirmed by the observations of temperature anisotropies and polarizations of cosmic microwave background radiation (Planck Collaboration et al., 2016), as well as scale-dependent galaxy bias on large-scales measured from galaxies (Giannantonio et al., 2014) and quasars (Leistedt et al., 2014).

The late time nonlinear gravitational evolution, however, induces phase coupling in the cosmic matter density and generates non-Gaussian features in the one-point probability distribution function (PDF) (Peebles, 1980; Juszkiewicz et al., 1993; Bernardeau, 1994). The PDFs measured from cosmological $N$-body simulations show a significant deviation from the Gaussian PDF reflecting the prominent nonlinear structures such as clusters, filaments, and cosmic voids (Hamilton, 1985; Bouchet et al., 1993; Kofman et al., 1994; Taylor & Watts, 2000; Kayo et al., 2001). These late-time non-Gaussian PDFs is directly observable from the cosmic shear measurement of weak lensing surveys (Kruse & Schneider, 2000; Clerkin et al., 2017; Takahashi et al., 2011). Quantifying the cosmic structure with the non-Gaussian PDF in the cosmic density field, therefore, is crucial to understand the nonlinear growth of large-scale structure. Upon exploiting the PDFs, one may tighten cosmological constraints on, for example, dark energy (Tatekawa & Mizuno, 2006; Seo et al., 2012; Codis et al., 2016).

Previous studies on the non-Gaussian PDF have suggested that the distribution of the cosmic density field follows approximately the log-normal PDF (Hubble, 1934; Coles & Jones, 1991; Kofman et al., 1994; Bernadeau & Kofman, 1995; Kayo et al., 2001). Meanwhile, an alternative fitting formula to the log-normal PDF was proposed by Colombi (1994), the so-called skewed log-normal PDF (Ueda & Yokoyama, 1996). Since then, pieces of evidence for the deviations from the log-normal distributions have come into sight relying on improved large-box-size, high-precision cosmological simulations (Szapudi & Pan, 2004; Pandey et al., 2013). More recently, Uhlemann et al. (2016) have analytically calculated the deviation of the logarithmic density PDF from the Gaussian one.

In this paper, we propose a new functional form to fit the non-Gaussian PDF: the generalized normal distribution of version 2 (${\mathcal{N}_{\mathrm{v2}}}$). The ${\mathcal{N}_{\mathrm{v2}}}$ PDF is a three-parameter extension of the Gaussian distribution incorporating the skewness. We show that the PDFs measured from the N-body simulations are well described by this model over a wide range of density, redshift, smoothing kernel radii, and cosmology.

2. Simulation

We run a suite of cosmological $N$-body simulations using the GOTPM code (Dubinski et al., 2004; Kim et al., 2009, 2011) with $2048^{3}$ particles in a cubic box of $L_{\rm box}=1024~{}h^{-1}$Mpc. The reference cosmology model (hereafter, $\Lambda_{m0}$) adopts the WMAP 5-year cosmology with ($\Omega_{m,0}$, $\Omega_{b,0}$, $\Omega_{DE,0}$, $w$) = ($0.26,~{}0.044,~{}0.74,-1$), $H_{0}=72$ km/s/Mpc, and $\sigma_{8}=0.79$, where $w$ is the equation of state parameter of dark energy. Also, we have run four simulations with spatially flat, but non-standard cosmologies: ($\Omega_{m,0}$, $\Omega_{b,0}$, $\Omega_{DE,0}$, $w$) = ($0.31,~{}0.044,~{}0.69,-1$), ($0.21,~{}0.044,~{}0.79,-1$), ($0.26,~{}0.044,~{}0.74,-1.5$), and ($0.26,~{}0.044,~{}0.74,-0.5$) to highlight the effect of total matter density and the equation of state of dark energy. We name these simulations, $\Lambda_{m-}$, $\Lambda_{m+}$, $Q_{w-}$, and $Q_{w+}$, respectively. Each simulation starts from $z=100$ with the initial conditions generated by the second-order Lagrangian perturbation theory (2LPT; Scoccimarro 1998; McCullagh et al. 2016; L’Huiller et al. 2014) with the linear power spectrum calculated from CAMB (Lewis et al., 2000). In this study, we have used snapshot particle data at six redshifts of $z=0,~{}0.2,~{}0.5,~{}1,~{}2$, and $4$.

From this set of particle data, we have measured the one-point PDF on $2048^{3}$ regular grid points laid over the simulation box with the spherical top-hat kernel with the radius of ${R_{\mathrm{th}}}=2,~{}5,~{}10$, and $25~{}h^{-1}$Mpc. The particle density is directly measured in real space (direct count in the spherical region).

3. One-point density distribution

3.1. Fitting the simulated PDF

After measuring the density, we calculate the probability distribution function of the density contrast $\delta\equiv\Delta\rho/\langle{\rho}\rangle$ as a function of its significance ${\nu_{\mathrm{sim}}}$ (${\nu_{\mathrm{sim}}}\equiv{\delta_{\mathrm{sim}}}/{\sigma_{\mathrm{sim}}}$ where ${\sigma_{\mathrm{sim}}}$ is the standard deviation of the density contrast). Hereafter, the subscript of ‘sim’ refers to the quantity directly measured from simulation data.

Figure 1 shows the PDFs measured at $z=0$ from the simulation with reference cosmology ($\Lambda_{m0}$) for four different ${R_{\mathrm{th}}}=2,~{}5,~{}10~{},25~{}{\rm Mpc}/h$. The PDF of density field smoothed with the narrower kernels deviates more from the Gaussian distribution (gray, solid line); the PDFs are skewed more toward the high-density (right-hand) side and present more kurtosis as ${R_{\mathrm{th}}}$ decreases ($\sigma_{\rm sim}$ gets larger). To facilitate the comparison in the lower density part of the PDF, we also show the same PDFs (but as a function of $1+{\delta_{\mathrm{sim}}}$) in Figure 2.

We fit the measured PDFs to the generalized normal distribution of version 2 (${\mathcal{N}_{\mathrm{v2}}}$),

[TABLE]

in which three parameters $\alpha,~{}\kappa,~{}\xi$ are used to parametrize the deviation from the normal distribution $\mathcal{N}$. Here, the distortion argument (${y_{\mathrm{sim}}}$) is defined as

[TABLE]

and $\alpha$, $\kappa$, and $\xi$ quantify, respectively, the scale, shape, and location of the skewed distribution ${\mathcal{N}_{\mathrm{v2}}}$. A positive (or negative) value of $\kappa$ yields left-skewed (or right-skewed) distributions with a sharp cut-off in the right (or left) distribution wing. Since ${\mathcal{N}_{\mathrm{v2}}}$ approaches $\mathcal{N}$ as $\kappa\rightarrow 0$, ${\mathcal{N}_{\mathrm{v2}}}$ is useful to describe deviations from $\mathcal{N}$ in a continuous manner. In addition, the cumulative distribution of ${\mathcal{N}_{\mathrm{v2}}}$ is the same as that of $\mathcal{N}$ and, consequently, ${\mathcal{N}_{\mathrm{v2}}}$ is a generalized version of the normal distribution. As can be seen in Figure 3, the measured density PDFs depend on redshift and smoothing length; therefore, the ${\mathcal{N}_{\mathrm{v2}}}$ parameters $\alpha$, $\kappa$, and $\xi$ must also be a function of redshift (more specifically, the linear growth factor, $D_{1}$) and ${R_{\mathrm{th}}}$.

We find the best-fitting parameters ($\alpha_{\rm fit}$, $\kappa_{\rm fit}$ and $\xi_{\rm fit}$) by applying the $\chi^{2}$-minimization method with a thousand density bins. Hereafter, the subscript of ‘fit’ refers to the best-fitting quantities. The resulting best-fitting PDFs to the reference simulation at $z=0$ are shown in Figure 1. As shown there, the overall shape of the simulated density PDF ($P_{\mathrm{sim}}$) is well fitted with ${\mathcal{N}_{\mathrm{v2}}}$ for a wide range of ${\nu_{\mathrm{sim}}}$ and ${R_{\mathrm{th}}}$.

We also compare $P_{\mathrm{sim}}$ with the log-normal and the skewed log-normal PDFs in Figure 2. The log-normal PDF ($P_{\mathrm{LN}}$) is defined as

[TABLE]

where the variance $\sigma_{1}$ can be derived by

[TABLE]

The skewed log-normal PDF ($P_{\mathrm{SLN}}$) combines the log-normal distribution and the Edgeworth expansion (e.g. Juszkiewicz et al. 1995). At third-order approximation (Colombi, 1994), $P_{\mathrm{SLN}}$ reads

[TABLE]

where $\upsilon\equiv\Phi/\sigma_{\Phi}$, $\Phi\equiv\ln(1+\delta)-\langle\ln(1+\delta)\rangle$, and $\sigma_{\Phi}$ is the variance of the log-density field $\Phi$. $H_{m}(\upsilon)$ is the Hermite polynomial of degree $m$, and $T_{3}$ and $T_{4}$ are the renormalized skewness and kurtosis of the field $\Phi$, respectively:

[TABLE]

While $P_{\mathrm{sim}}$ is well reproduced by all the ${\mathcal{N}_{\mathrm{v2}}}$, $P_{\mathrm{LN}}$, and $P_{\mathrm{SLN}}$ in the high-density regions, the low-density cliffs are better fitted by ${\mathcal{N}_{\mathrm{v2}}}$ and $P_{\mathrm{SLN}}$ than by $P_{\mathrm{LN}}$. For the smaller ${R_{\mathrm{th}}}$, the deviation between $P_{\mathrm{sim}}$ and $P_{\mathrm{LN}}$ in the underdense regions becomes more prominent. It is consistent with the analysis by Ueda & Yokoyama (1996) and the perturbation theory presented by Bernadeau & Kofman (1995). Although the fits by ${\mathcal{N}_{\mathrm{v2}}}$ also differ from $P_{\mathrm{sim}}$ in the underdense regions, in particular for the smaller ${R_{\mathrm{th}}}$, the deviation is much milder than that of $P_{\mathrm{LN}}$.

3.2. Fitting to $\alpha$, $\kappa$, and $\xi$

All best-fitting parameters ($\alpha_{\rm fit}$, $\kappa_{\rm fit}$, and $\xi_{\rm fit}$) vary with redshifts (or $D_{1}$) and ${R_{\mathrm{th}}}$. We therefore compile these fitting values at six redshifts ($z=$ 0, 0.2, 0.5, 1, 2, and 4), four smoothing radii (${R_{\mathrm{th}}}=$ 2, 5, 10, and 25 $h^{-1}$Mpc) for five different simulated models ($\Lambda_{m0}$, $\Lambda_{m-}$, $\Lambda_{m+}$, $Q_{w-}$, and $Q_{w+}$). We then find functions that incorporate the redshift and smoothing-scale dependence of the best-fitting parameters by using the Eureqa software 111http://www.nutonian.com/products/eureqa.

Among all possible fitting forms to $\alpha_{\mathrm{fit}}({R_{\mathrm{th}}},D_{1})$, $\kappa_{\mathrm{fit}}({R_{\mathrm{th}}},D_{1})$, and $\xi_{\mathrm{fit}}({R_{\mathrm{th}}},D_{1})$, we select those which satisfy the following criteria: (1) the functional form must be the same for all the cosmological models, (2) the resulting PDF must asymptote to the normal distribution at early times and for large smoothing radius; that is, $\alpha_{\mathrm{fit}}\to 1$, $\kappa_{\mathrm{fit}}\to 0$ and $\xi_{\mathrm{fit}}\to 0$ as ${R_{\mathrm{th}}}\rightarrow\infty$ and/or $D_{1}\rightarrow 0$, (3) the R2 goodness of fit should be larger than $0.9$, and (4) the fitting equation has the least number of coefficients. Note that the empirical relations that we find here do not necessarily reflect the physical origin of the underlying function. The final functions that we find are

[TABLE]

respectively, where $a_{i}$, $b_{i}$, and $c_{i}$ are numerical coefficients.

Table 1 lists the coefficients and their $R^{2}$ goodness of fit value for the simulated models. The predicted ${\mathcal{N}_{\mathrm{v2}}}$ by $\alpha_{\mathrm{fit}}({R_{\mathrm{th}}},D_{1})$, $\kappa_{\mathrm{fit}}({R_{\mathrm{th}}},D_{1})$, and $\xi_{\mathrm{fit}}({R_{\mathrm{th}}},D_{1})$, hereafter ${\mathcal{N}_{\mathrm{v2}}}(\alpha_{\mathrm{fit}},\kappa_{\mathrm{fit}},\xi_{\mathrm{fit}})$, are compared with the corresponding $P_{\mathrm{sim}}$ and the log-normal distribution in Figure 3. ${\mathcal{N}_{\mathrm{v2}}}(\alpha_{\mathrm{fit}},\kappa_{\mathrm{fit}},\xi_{\mathrm{fit}})$ reproduce well the overall shape of $P_{\mathrm{sim}}$ for a wide range of ${\nu_{\mathrm{sim}}}$, ${R_{\mathrm{th}}}$, redshift, and cosmology. The PDFs over the entire density scale are better reproduced by ${\mathcal{N}_{\mathrm{v2}}}(\alpha_{\mathrm{fit}},\kappa_{\mathrm{fit}},\xi_{\mathrm{fit}})$ than with the log-normal distribution.

3.3. Skewness of fitted PDFs

The density fluctuations in the very early universe are known to be indistinguishable from Gaussian to within measurement error. However, gravity is expected to skew the density distribution, making a lognormal, skewed lognormal, or ${\mathcal{N}_{\mathrm{v2}}}$ a better fit than a Gaussian even at early times (Peebles, 1980; Fry , 1984; Juszkiewicz et al., 1993; Bernardeau, 1994). It should be noted, though, that for small variance at early times, the skewness has a negligible effect on the actual density PDF. Eulerian perturbation theory (EPT) predicts that a reduced skewness parameter $S_{3}\equiv\left<\delta^{3}\right>/{\sigma}^{4}$ in an Einstein de sitter (EdS) universe approaches to $\sim 34/7-(n+3)$ at early times, where $n$ is an index of power-law spectrum (Peebles, 1980; Juszkiewicz et al., 1993; Bernardeau, 1994; Fry & Scherrer , 1994). Since the log-normal PDF has a non-zero skewness as $S_{3}=3$ at early times, the log-normal distribution has been proposed as a better fit than a Gaussian to the initial PDF (Coles & Jones, 1991; Colombi, 1994; Neyrinck, 2013). Figure 4 compares the $S_{3}$ values of ${\mathcal{N}_{\mathrm{v2}}}(\alpha_{\mathrm{fit}},\kappa_{\mathrm{fit}},\xi_{\mathrm{fit}})$ and $P_{\mathrm{LN}}$ for the $\Lambda$CDM model ($\Lambda_{m0}$) at $R_{\mathrm{th}}=10$ Mpc/h. Although the ${\mathcal{N}_{\mathrm{v2}}}(\alpha_{\mathrm{fit}},\kappa_{\mathrm{fit}},\xi_{\mathrm{fit}})$ is chosen to converge to the normal distribution at early times (second condition of §3.2), the $S_{3}$ value approach to a non-zero value of $\sim 3.6$. The $S_{3}$ values directly measured from our simulation data (black circles) closely follow that of ${\mathcal{N}_{\mathrm{v2}}}(\alpha_{\mathrm{fit}},\kappa_{\mathrm{fit}},\xi_{\mathrm{fit}}$) rather than of $P_{\mathrm{LN}}$. To calculate the $S_{3}$ value at the higher redshift, we generate an initial Gaussian density field at $z=19$ evolved by the 2LPT. Following the $S_{3}$ trend of the simulation data, the $S_{3}$ value at $z=19$ (blue filled square) results in $S_{3}\sim 3.4$, which is slightly smaller than that of ${\mathcal{N}_{\mathrm{v2}}}(\alpha_{\mathrm{fit}},\kappa_{\mathrm{fit}},\xi_{\mathrm{fit}}$).

3.4. Sensitivity of fitted PDFs to cosmology

Relative differences of PDFs for four non-standard models ($\Lambda_{m-}$, $\Lambda_{m+}$, $Q_{w-}$, $Q_{w+}$) relative to the $\Lambda$CDM model ($\Lambda_{m0}$) are shown in Figure 5. The differences ($P-P_{\Lambda_{m0}}$) compiled by both $P_{\mathrm{sim}}$ and $P_{\mathrm{fit}}$ are compared to each other in order to check how $P_{\mathrm{fit}}$ capture the different models. $P-P_{\Lambda_{m0}}$ are well reproduced at high redshift and/or the large smoothing. However, $P-P_{\Lambda_{m0}}$ for smaller redshift or the small smoothing show significant deviations from that of $P_{\mathrm{sim}}$. Thus, the ${\mathcal{N}_{\mathrm{v2}}}$ fits are not accurate enough to make the distinction between the models for the strongly nonlinear regime. The failure is due to poor fits of the PDFs in the underdense region (see Fig. 2).

4. Summary & Discussion

In this paper, we presented the one-point PDFs measured from cosmological $N$-body simulations and showed that the new fitting formula based on the generalized normal distribution of version 2 (${\mathcal{N}_{\mathrm{v2}}}$) provides significantly better fit compared to the log-normal distribution. In particular, ${\mathcal{N}_{\mathrm{v2}}}$ reproduces well the overall PDFs for a wide range of density, smoothing kernel, redshift, and cosmology, except in strongly nonlinear regimes. The improvement by the ${\mathcal{N}_{\mathrm{v2}}}$ is substantial in the under-dense regions, which is also achieved by the skewed log-normal distribution, the third-order Edgeworth expansion of the log-normal distribution (Colombi, 1994; Ueda & Yokoyama, 1996).

As the ${\mathcal{N}_{\mathrm{v2}}}$ distribution can accommodate a continuous transition from the initial Gaussian distribution function, the result we present here should pave the way to modeling the density PDF in the quasi-linear regimes where perturbation theory (Bernardeau et al., 2002) captures the nonlinear evolution of cosmic density fields.

The simulated PDFs and their fitted curves by ${\mathcal{N}_{\mathrm{v2}}}$ for various smoothing kernels (${R_{\mathrm{th}}}=2$, 5, 10, and 25 $h^{-1}$Mpc), redshifts ($z=0$, 0.2, 0.5, 1, 2 and 4), and cosmologies ($\Lambda_{m0}$, $\Lambda_{m-}$, $\Lambda_{m+}$, $Q_{w-}$, and $Q_{w+}$) can be obtained from the Web site of the first author at https://astro.kias.re.kr/jhshin.

We appreciate the anonymous referee for his/her helpful comments that improved our manuscript. We thank the Korea Institute for Advanced Study for providing computing resources (KIAS Center for Advanced Computation) for this work. We thank C. Uhlemann for comments. CP’s research is partially funded by Spin(e) ANR-13-BS05-0005. DJ acknowledges support from National Science Foundation grant AST-1517363.