<figure><img src="image_1.png" /></figure>

# New Fitting Formula for Cosmic Nonlinear Density Distribution

Jihye Shin¹, Juhan Kim², Christophe Pichon¹,³, Donghui Jeong⁴, and Changbom Park¹

¹School of Physics, Korea Institute for Advanced Study, 85 Hoegiro, Dongdaemun-gu, Seoul 130-722, Korea

²Center for Advanced Computation, Korea Institute for Advanced Study, 85 Hoegiro-ro, Dongdaemun-gu, Seoul 130-722, Korea

³CNRS and UPMC Université Paris 06, UMR 7095, Institut d'Astrophysique de Paris, 98 bis Boulevard Arago, Paris F-75014, France

⁴Department of Astronomy and Astrophysics, and Institute for Gravitation and the Cosmos, The Pennsylvania State University, University Park, PA 16802, USA

*Received 2017 February 4; revised 2017 May 1; accepted 2017 May 18; published 2017 July 5*

## Abstract

We have measured the probability distribution function (PDF) of a cosmic matter density field from a suite of N-body simulations. We propose the generalized normal distribution of version 2 ($N_{v2}$) as an alternative fitting formula to the well-known log-normal distribution. We find that $N_{v2}$ provides a significantly better fit than that of the log-normal distribution for all smoothing radii (2, 5, 10, 25 [Mpc $h^{-1}$]) that we studied. The improvement is substantial in the underdense regions. The development of non-Gaussianities in the cosmic matter density field is captured by continuous evolution of the skewness and shift parameters of the $N_{v2}$ 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 the web site https://astro.kias.re.kr/jhshin.

**Key words:** cosmology: theory – large-scale structure of universe – methods: numerical

## 1. Introduction

The inflationary models (Starobinskiǐ 1979, 1982; Guth 1981; Sato 1981; Albrecht & Steinhardt 1982; Linde 1982) of the early universe predict that the primordial density perturbations generated during inflation (Mukhanov & Chibisov 1981; Guth & Pi 1982; Hawking 1982; Bardeen et al. 1983) must obey nearly Gaussian statistics (Acquaviva et al. 2003; Maldacena 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 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 are directly observable from the cosmic shear measurement of weak-lensing surveys (Kruse & Schneider 2000; Takahashi et al. 2011; Clerkin et al. 2017). Quantifying the cosmic structure with the non-Gaussian PDF in the cosmic density field, therefore, is crucial for understanding 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; Bernardeau & 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 ($N_{v2}$). The $N_{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 described well 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³ particles in a cubic box of $L_{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$^{-1}$ Mpc$^{-1}$, 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; L'Huillier et al. 2014; McCullagh et al. 2016), and 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.
<figure><img src="image_2.png" /><figcaption>**Figure 1.** First panel: PDFs with various smoothing radii at $z=0$ in the reference model of $\Lambda_{m0}$. The filled circles are obtained from the simulation ($P_{sim}$), while the solid lines are the corresponding fitted curves ($P_{fit}$), a generalized normal distribution (Ver. 2). The log-normal distribution is overplotted as the dashed lines. In the legend, we show the radius of the smoothing top-hat kernel. For a reference, we show the Gaussian distribution in the thick gray curve. Second panel: error bars comparing $P_{sim}$ to the corresponding smoothing radii in the first panel. For a clear appearance, offsets are used in the x-axis. Third and fourth panels: relative residuals of $P_{fit}$ (third) and log-normal distributions (fourth) comparing $P_{sim}$ to the corresponding smoothing radii in the first panel.</figcaption></figure>

From this set of particle data, we have measured the one-point PDF on 2048³ regular grid points laid over the simulation box with the spherical top-hat kernel with a radius of $R_{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 PDF of the density contrast $\delta \equiv \Delta\rho/\langle\rho\rangle$ as a function of its significance $\nu_{sim}$ ($\nu_{sim} \equiv \delta_{sim}/\sigma_{sim}$, where $\sigma_{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_{th} = 2, 5, 10, 25 h^{-1}$ Mpc smoothed with the narrower kernels deviates more from the Gaussian distribution (gray, solid line); the PDFs are skewed more toward the high-density (right) side and present more kurtosis as $R_{th}$ decreases ($\sigma_{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_{sim}$) in Figure 2.

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Same as the top panel of Figure 1 but as a function of $1 + \delta_{sim}$. The x-axis is scaled in log 1 + $\delta_{sim}$ to show the PDFs of the underdense regions in detail. It clearly shows how the fit fails in the underdense nonlinear region (blue curve).</figcaption></figure>

We fit the measured PDFs to the generalized normal distribution of version 2 ($N_{v2}$),

$$N_{v2}(\nu_{sim}) = \frac{N(y_{sim})}{\alpha + \kappa(\nu_{sim} - \xi)}, \qquad (1)$$

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

$$y_{sim} = \begin{cases} \frac{1}{\kappa} \ln \left[ 1 + \frac{\kappa(\nu_{sim} - \xi)}{\alpha} \right] & \text{if } \kappa \neq 0 \\ \frac{\nu_{sim} - \xi}{\alpha} & \text{otherwise,} \end{cases} \qquad (2)$$

and $\alpha$, $\kappa$, and $\xi$ quantify, respectively, the scale, shape, and location of the skewed distribution $N_{v2}$. A positive (or negative) value of $\kappa$ yields left-skewed (or right-skewed) distributions with a sharp cutoff in the right (or left) distribution wing. Since $N_{v2}$ approaches $N$ as $\kappa \to 0$, $N_{v2}$ is useful for describing deviations from $N$ in a continuous manner. In addition, the cumulative distribution of $N_{v2}$ is the same as that of $N$, and consequently, $N_{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 $N_{v2}$ parameters $\alpha$, $\kappa$, and $\xi$ must also be a function of redshift (more specifically, the linear growth factor, $D_1$) and $R_{th}$.

We find the best-fitting parameters ($\alpha_{fit}$, $\kappa_{fit}$, and $\xi_{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_{sim}$) is well fitted, with $N_{v2}$ for a wide range of $\nu_{sim}$ and $R_{th}$.
<figure><img src="image_4.png" /><figcaption>**Figure 3.** PDFs with various smoothing radii at $z = 0, 1$, and 4 for the five simulated models. $P_{sim}$ (filled circles) are compared with the corresponding predicted PDFs, $N_{v2}(\alpha_{fit}, \kappa_{fit}, \xi_{fit})$ (solid lines), and the log-normal distributions (dashed lines). In the first row, cumulative distribution functions (CDFs) of $P_{sim}$ (filled circles) are compared with the corresponding CDFs of $N_{v2}(\alpha_{fit}, \kappa_{fit}, \xi_{fit})$, which are the same as the CDFs of $N(y_{sim})$.</figcaption></figure>
<table><caption>**Table 1**
Numerical Coefficients of Fitted Functions for $\alpha$, $\kappa$, and $\xi$</caption><thead><tr><th>Model</th><th>Numerical Coefficients</th><th>R<sup>2</sup> Goodness</th></tr></thead><tbody><tr><td rowspan="3">Λ<sub>m0</sub></td><td>a<sub>1</sub> = 32.92</td><td>0.992</td></tr><tr><td>b<sub>1</sub> = 10.38, b<sub>2</sub> = 2.118, b<sub>3</sub> = 3.151</td><td>0.994</td></tr><tr><td>c<sub>1</sub> = -4.082, c<sub>2</sub> = 44.92</td><td>0.969</td></tr><tr><td rowspan="3">Λ<sub>m-</sub></td><td>a<sub>1</sub> = 35.48</td><td>0.992</td></tr><tr><td>b<sub>1</sub> = 11.85, b<sub>2</sub> = 3.132, b<sub>3</sub> = 3.683</td><td>0.993</td></tr><tr><td>c<sub>1</sub> = -4.370, c<sub>2</sub> = 54.48</td><td>0.985</td></tr><tr><td rowspan="3">Λ<sub>m+</sub></td><td>a<sub>1</sub> = 30.82</td><td>0.992</td></tr><tr><td>b<sub>1</sub> = 9.369, b<sub>2</sub> = 1.426, b<sub>3</sub> = 2.787</td><td>0.993</td></tr><tr><td>c<sub>1</sub> = -3.901, c<sub>2</sub> = 39.85</td><td>0.972</td></tr><tr><td rowspan="3">Q<sub>w-</sub></td><td>a<sub>1</sub> = 31.54</td><td>0.992</td></tr><tr><td>b<sub>1</sub> = 10.22, b<sub>2</sub> = 2.920, b<sub>3</sub> = 2.124</td><td>0.996</td></tr><tr><td>c<sub>1</sub> = -3.929, c<sub>2</sub> = 41.41</td><td>0.980</td></tr><tr><td rowspan="3">Q<sub>w+</sub></td><td>a<sub>1</sub> = 32.05</td><td>0.994</td></tr><tr><td>b<sub>1</sub> = 10.30, b<sub>2</sub> = 2.890, b<sub>3</sub> = 2.250</td><td>0.996</td></tr><tr><td>c<sub>1</sub> = -4.024, c<sub>2</sub> = 45.35</td><td>0.980</td></tr></tbody></table>

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

$$P_{LN}(\delta_{sim}) = \frac{1}{(2\pi\sigma_1^2)^{1/2}} \frac{1}{1 + \delta_{sim}} \times \exp\left\{\frac{[\ln(1 + \delta_{sim}) + \sigma_1^2/2]^2}{2\sigma_1^2}\right\}, \quad (3)$$

where the variance $\sigma_1^2$ can be derived by

$$\sigma_1^2 = \ln[1 + \sigma_{sim}^2]. \quad (4)$$

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

$$P_{SLN}(v_{sim}) = \left[ 1 + \frac{1}{3!} T_3 \sigma_{\Phi} H_3(v_{sim}) + \frac{1}{4!} T_4 \sigma_{\Phi}^2 H_4(v_{sim}) + \frac{10}{6!} T_3^2 \sigma_{\Phi}^2 H_6(v_{sim}) \right] N(v_{sim}), \quad (5)$$

where $v \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(v)$ is the Hermite polynomial of degree $m$, and $T_3$ and $T_4$ are the renormalized skewness and kurtosis of the field $\Phi$, respectively:

$$T_3 \equiv \frac{\langle \Phi^3 \rangle}{\sigma_{\Phi}^4}, \quad T_4 \equiv \frac{\langle \Phi^4 \rangle - 3\sigma_{\Phi}^4}{\sigma_{\Phi}^6}. \quad (6)$$

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

<figure><img src="image_5.png" /><figcaption>**Figure 4.** $S_3$ values of $N_{v2}(\alpha_{fit}, \kappa_{fit}, \xi_{fit})$ (red line) and $P_{LN}$ (black line) as a function of $1+z$. The horizontal lines indicate $S_3$ values in an EdS universe with $n = -3$ (dotted line) and in the CDM universe (dashed line), predicted by the EPT. The black filled circles represent $S_3$ values measured from our simulation data, while the blue filled square indicates that from the 2LPT density field at $z=19$. Here, the same smoothing radius as $R_{th} = 10 \text{ Mpc h}^{-1}$ is adopted for all the $S_3$ values.</figcaption></figure>

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

All best-fitting parameters ($\alpha_{fit}$, $\kappa_{fit}$, and $\xi_{fit}$) vary with redshifts (or $D_1$) and $R_{th}$. We therefore compile these fitting values at six redshifts ($z=0, 0.2, 0.5, 1, 2$, and $4$) and four smoothing radii ($R_{th} = 2, 5, 10$, and $25 \text{ h}^{-1} \text{ 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 using the Eureqa software.<fnref n="5" />

Among all possible fitting forms to $\alpha_{fit}(R_{th}, D_1)$, $\kappa_{fit}(R_{th}, D_1)$, and $\xi_{fit}(R_{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_{fit} \to 1$, $\kappa_{fit} \to 0$, and $\xi_{fit} \to 0$ as $R_{th} \to \infty$ and/or $D_1 \to 0$, (3) the $R^2$ 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

$$\alpha_{fit}(R_{th}, D_1) = \frac{R_{th}^2}{R_{th}^2 + a_1 D_1^2 + R_{th} D_1^2}, \quad (7)$$

$$\kappa_{fit}(R_{th}, D_1) = \frac{b_1 D_1}{b_2 + R_{th} + b_3 D_1}, \quad (8)$$

$$\xi_{fit}(R_{th}, D_1) = \frac{c_1 R_{th} D_1}{R_{th} + R_{th}^2 + c_2 D_1^2}, \quad (9)$$

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 $N_{v2}$ by $\alpha_{fit}(R_{th}, D_1)$, $\kappa_{fit}(R_{th}, D_1)$, and $\xi_{fit}(R_{th}, D_1)$, hereafter $N_{v2}(\alpha_{fit}, \kappa_{fit}, \xi_{fit})$, are compared with the corresponding $P_{sim}$ and the log-normal distribution in Figure 3. $N_{v2}(\alpha_{fit}, \kappa_{fit}, \xi_{fit})$ reproduces the

<footnote n="5">⁵ http://www.nutonian.com/products/eureqa</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Differences between PDFs ($P - P_{\Lambda_{m0}}$) for various simulated models ($\Lambda_{m-}$, $\Lambda_{m+}$, $Q_{w-}$, $Q_{w+}$) relative to the reference model ($\Lambda_{m0}$). The thick solid lines are obtained from the simulations ($P_{sim}$), while the dashed lines are from the corresponding fitted curves by $N_{v2}$ ($P_{fit}$). For smaller radii and lower redshift, the model fails to capture the difference between the different cosmologies and the reference ΛCDM model.</figcaption></figure>

overall shape of $P_{sim}$ for a wide range of $\nu_{sim}$, $R_{th}$, redshift, and cosmology. The PDFs over the entire density scale are better reproduced by $N_{v2}(\alpha_{fit}, \kappa_{fit}, \xi_{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 $N_{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 \langle\delta^3\rangle/\sigma^4$ in an Einstein de Sitter (EdS) universe approaches $\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 $N_{v2}(\alpha_{fit}, \kappa_{fit}, \xi_{fit})$ and $P_{LN}$ for the ΛCDM model ($\Lambda_{m0}$) at $R_{th} = 10$ Mpc h$^{-1}$. Although the $N_{v2}(\alpha_{fit}, \kappa_{fit}, \xi_{fit})$ is chosen to converge to the normal distribution at early times (second condition of Section 3.2), the $S_3$ value approaches a non-zero value of $\sim 3.6$. The $S_3$ values directly measured from our simulation data (black circles) closely follow that of $N_{v2}(\alpha_{fit}, \kappa_{fit}, \xi_{fit})$ rather than that of $P_{LN}$. To calculate the $S_3$ value at 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 $N_{v2}(\alpha_{fit}, \kappa_{fit}, \xi_{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 ΛCDM model ($\Lambda_{m0}$) are shown in Figure 5. The differences ($P - P_{\Lambda_{m0}}$) compiled by both $P_{sim}$ and $P_{fit}$ are compared to each other in order to check how $P_{fit}$ captures the different models. $P - P_{\Lambda_{m0}}$ are well-reproduced at high redshift and/or large smoothing. However,
$P - P_{\Lambda m0}$ for smaller redshift or small smoothing show significant deviations from that of $P_{sim}$. Thus, the $N_{v2}$ fits are not accurate enough to make a distinction between the models for the strongly nonlinear regime. The failure is due to poor fits of the PDFs in the underdense region (see Figure 2).

## 4. Summary and Discussion

In this paper, we presented 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 ($N_{v2}$) provides a significantly better fit compared to the log-normal distribution. In particular, $N_{v2}$ reproduces the overall PDFs for a wide range of density, smoothing kernel, redshift, and cosmology, except in strongly nonlinear regimes. The improvement by $N_{v2}$ is substantial in the underdense regions, and is also achieved by the skewed log-normal distribution, which is the third-order Edgeworth expansion of the log-normal distribution (Colombi 1994; Ueda & Yokoyama 1996).

As the $N_{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 $N_{v2}$ for various smoothing kernels ($R_{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 https://astro.kias.re.kr/jhshin.

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. C.P.'s research is partially funded by Spin(e) ANR-13-BS05-0005. D.J. acknowledges support from National Science Foundation grant AST-1517363.

## References

* Acquaviva, V., Bartolo, N., Matarrese, S., & Riotto, A. 2003, NuPhB, 667, 119

* Albrecht, A., & Steinhardt, P. J. 1982, PhRvL, 48, 1220

* Bardeen, J. M., Steinhardt, P. J., & Turner, M. S. 1983, PhRvD, 28, 679

* Bernardeau, F. 1994, ApJ, 433, 1

* Bernardeau, F., Colombi, S., Gaztañaga, E., & Scoccimarro, R. 2002, PhR, 367, 1

* Bernardeau, F., & Kofman, L. 1995, ApJ, 443, 479

* Bouchet, F., Strauss, M. A., Davis, M., et al. 1993, ApJ, 417, 36

* Clerkin, L., Kirk, D., Manera, M., et al. 2017, MNRAS, 466, 1444

* Codis, S., Pichon, C., Bernardeau, F., Uhlemann, C., & Prunet, S. 2016, MNRAS, 460, 1549

* Coles, P., & Jones, B. 1991, MNRAS, 248, 1

* Colombi, S. 1994, ApJ, 435, 536

* Creminelli, P., & Zaldarriaga, M. 2004, JCAP, 10, 006

* Dubinski, J., Kim, J., Park, C., & Humble, R. 2004, NewA, 9, 111

* Fry, J. N. 1984, ApJ, 279, 499

* Fry, J. N., & Scherrer, R. J. 1994, ApJ, 429, 36

* Giannantonio, T., Ross, A. J., Percival, W. J., et al. 2014, PhRvD, 89, 023511

* Guth, A. H. 1981, PhRvD, 23, 347

* Guth, A. H., & Pi, S.-Y. 1982, PhRvL, 49, 1110

* Hawking, S. W. 1982, PhLB, 115, 295

* Hamilton, A. J. S. 1985, ApJL, 292, L35

* Hubble, E. R. 1934, ApJ, 79, 8

* Juszkiewicz, R., Bouchet, F., & Colombi, S. 1993, ApJL, 419, L9

* Juszkiewicz, R., Weinberg, D. H., Amsterdamski, P., Chodorowski, M., & Bouchet, F. 1995, ApJ, 442, 39

* Kayo, I., Taruya, A., & Suto, Y. 2001, ApJ, 561, 22

* Kofman, L., Bertschinger, E., Gelb, J. M., Nusser, A., & Dekel, A. 1994, ApJ, 420, 44

* Kim, J., Park, C., Gott, J. R., III, & Dubinski, J. 2009, ApJ, 701, 1547

* Kim, J., Park, C., Lee, S. M., & Gott, J. R., III 2011, JKAS, 44, 217

* Kruse, G., & Schneider, P. 2000, MNRAS, 318, 321

* Leistedt, B., Peiris, H. V., & Roth, N. 2014, PhRvL, 113, 221301

* Lewis, A., Challinor, A., & Lasenby, A. 2000, ApJ, 538, 473

* L'Huillier, B., Park, C., & Kim, J. 2014, NewA, 30, 79

* Linde, A. D. 1982, PhLB, 108, 389

* Maldacena, J. 2003, JHEP, 5, 013

* McCullagh, N., Jeong, D., & Szalay, A. S. 2016, MNRAS, 455, 2945

* Mukhanov, V. F., & Chibisov, G. V. 1981, JETPL, 33, 532

* Neyrinck, M. C. 2013, MNRAS, 428, 141

* Pandey, B., White, S. D. M., Springel, V., & Angulo, R. E. 2013, MNRAS, 435, 2968

* Peebles, P. J. E. 1980, The Large-Scale Structure of the Universe (Princeton, NJ: Princeton Univ. Press)

* Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 594 A17

* Sato, K. 1981, MNRAS, 195, 467

* Scoccimarro, R. 1998, MNRAS, 299, 1097

* Seo, H.-J., Sato, M., Takada, M., & Dodelson, S. 2012, ApJ, 748, 57

* Starobinski, A. A. 1979, JETPL, 30, 682

* Starobinski, A. A. 1982, PhLB, 117, 175

* Szapudi, I., & Pan, J. 2004, ApJ, 602, 26

* Takahashi, R., Oguri, M., Sato, M., & Hamana, T. 2011, ApJ, 742, 15

* Taylor, A. N., & Watts, P. I. R. 2000, MNRAS, 314, 92

* Tatekawa, T., & Mizuno, S. 2006, JCAP, 2, 006

* Ueda, H., & Yokoyama, J. 1996, MNRAS, 280, 754

* Uhlemann, C., Codis, S., Pichon, C., Bernardeau, F., & Reimberg, P. 2016, MNRAS, 460, 1529