# One-point probability distribution function (PDF) from spherical collapse:
Early Dark Energy (EDE) vs. $\Lambda$CDM

Ankush Mandal1,2 and Sharvari Nadkarni-Ghosh1
1Department of Physics, I.I.T. Kanpur, Kanpur, U.P. 208016 India
2Inter-University Centre for Astronomy and Astrophysics, Pune-411007, India
E-mail:ankushm<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We compute the one-point PDF of an initially Gaussian dark matter density
field using spherical collapse (SC). We compare the results to other forms
available in the literature and also compare the PDFs in the $\Lambda$CDM
model with an early dark energy (EDE) model. We find that the skewed log-
normal distribution provides the best fit to the non-linear PDF from SC for
both cosmologies, from $a=0.1$ to 1 and for scales characterized by the
comoving width of the Gaussian: $\sigma_{G}=0.5,1,2$. To elucidate the effect
of cosmology, we examine the linear and non-linear growth rates through test
cases. For overdensities, when the two models have the same initial density
contrast, the differences due to cosmology are amplified in the non-linear
regime, whereas, if the two models have the same linear density contrast
today, then the differences in cosmology are damped in the non-linear regime.
This behaviour is in contrast with voids, where the non-linear growth becomes
‘self-regulatory’ and is less sensitive to cosmology and initial conditions.
To compare the PDFs, we examine the difference of the PDFs and evolution of
the width of the PDF. The trends with scale and redshift are as expected. A
tertiary aim of this paper was to check if the fitting form for the non-linear
density-velocity divergence relation, derived for constant equation of state
($w$) models by Nadkarni-Ghosh holds for the EDE model. We find that it does
with an accuracy of 4%, thus increasing its range of validity.

###### keywords:

cosmology: theory - large-scale structure of Universe - dark energy

††pagerange: One-point probability distribution function (PDF) from spherical
collapse: Early Dark Energy (EDE) vs. $\Lambda$CDM–References

##  1 Introduction

The large scale structure observed today is believed to have grown from tiny
fluctuations laid down during the inflationary era in the history of the
Universe. Observations of the cosmic microwave background suggest that the
one-point probability distribution function (PDF) of these primordial
fluctuations is very close to a Gaussian (Planck
Collaboration111https://www.cosmos.esa.int/web/planck, Planck 2018 results.
IX.). At early epochs, when the evolution is linear, the shape of the PDF
remains Gaussian and hence symmetric about the origin. However, this symmetry
is broken in the non-linear regime. For instance, the density contrast in
overdense virialized regions is of the order of a hundred, but it never drops
below -1 in void regions, since the total density is always positive. Tracking
this (secondary) non-Gaussianity is important because the shape of the PDF
carries information about the evolution history and hence can be used to
constrain cosmological parameters. While the two-point and higher-order
correlation functions give information about the spectral distribution, there
is still observational interest in additionally measuring the matter density
PDF, either directly in terms of the galaxy counts-in-cell distribution or
indirectly through the lensing convergence (Bel et al., 2016; Liu et al.,
2016; Patton et al., 2017; Hurtado-Gil et al., 2017; Clerkin et al., 2017;
Gruen et al., 2018; Friedrich et al., 2018). This has been shown to provide
complementary constraints and help break degeneracies in parameter estimation.

The history of the one-point PDF begins almost a century ago, when Hubble
(1934) examined the frequency distribution of about 44000 galaxies distributed
over three-fourth of the sky and found it to be close to log-normal. Following
this, there were other many other statistical measures proposed to quantify
the departure from Gaussianity, including the extensive work by Neyman (1962)
and Scott (1962). A detailed historical summary can be found in Peebles
(1980). The early 80s and 90s saw the advent of phenomenological models some
of which assumed that gravitational clustering can be described by equilibrium
thermodynamics (Saslaw & Hamilton, 1984; Suto et al., 1990; Lahav et al.,
1993). Another phenomenological model which has remained popular till date is
the log-normal model proposed by Coles & Jones (1991). This model is partly
motivated by the theoretical fact that the full non-linear continuity equation
combined with linear statistics for the velocity field, predicts a log-normal
statistics for the density field. This form has proven to be a good fit in the
quasi-linear regime both theoretically (Kofman, 1993; Coles, Melott, &
Shandarin, 1993; Ohta, Kayo, & Taruya, 2003) and observationally (Colombi et
al., 2000; Szapudi et al., 2000; Wild et al., 2005; Clerkin et al., 2017). It
has also been proposed to fit the weak lensing convergence field Taruya et al.
2002; Hilbert et al. 2011) as well as to predict the void distribution
function (Russell & Pycke, 2017). However, its validity in the non-linear
regime is less established (for e.g., Colombi 1994; Bernardeau & Kofman 1995;
Ueda & Yokoyama 1996; Szapudi & Pan 2004; Repp & Szapudi 2018). Recently,
Xavier et al. (2016) have also pointed out a mathematical inconsistency in
simultaneous modelling the joint distribution of the density and convergence
fields and have proposed modifications to the standard form. Numerical
simulations have also shown deviations from this form and alternate fitting
forms and schemes have been proposed (for e.g., Shin et al. 2017; Klypin et
al. 2018; Tosone et al. 2020).

On the theoretical side, analytical modelling of the one-point PDF has been of
interest for over three decades. In the 90s, Kofman and Bernardeau and their
collaborators used perturbation theory (PT) in the Eulerian and Lagrangian
frames to derive the density PDFs and its higher-order moments (Kofman, 1993;
Kofman et al., 1994; Bernardeau, 1994a, b; Bernardeau & Kofman, 1995; Colombi
et al., 1997). Extensions of this work were done by Fosalba & Gaztanaga
(1998b) who simplified the higher-order perturbative terms using spherical
dynamics so as to apply the formalism to non-Gaussian initial conditions.
Juszkiewicz et al. (1995) and Gaztañaga et al. (2000) combined PT with the
Edgeworth expansion to get better approximations in the non-linear regime.
Taylor & Watts (2000) combined second-order PT with the Chapman-Kolmogorov
equation to understand the evolution of the distribution. Estimates based on
local dynamics such as spherical and ellipsoidal collapse were given by Ohta
et al. (2003, 2004) who formulated the differential equations for the
evolution of the one-point PDFs as well as Lam & Sheth (2008a, b) who used
excursion sets to construct the PDF in both real and redshift space.

One non-perturbative approach was put forth by Valageas who recast the
solution to the collisionless Boltzmann equation as a path integral to be
solved within the saddle-point approximation using the method of steepest
descent. The dynamics near the saddle point depends only on the scale
separation because of the homogeneity and isotropy of the Universe and
spherical symmetry is invoked to provide a unique mapping between the initial
linear state and the final non-linear state of the system (Valageas, 1998,
2001, 2002a, 2002b, 2002c, 2002d; Valageas & Munshi, 2004). This was further
developed and applied to compute the statistic of density profiles by
Bernardeau et al. (2014, 2015). Very recently, Ivanov et al. (2019) have
extended this calculation to compute the effect of fluctuations around the
spherical saddle point, which they argue is necessary for the consistency of
the theory. Bernardeau & Reimberg (2016) have recently proposed another non-
perturbative approach that uses Large Deviations Theory to compute the
cumulant generating functions. This approach also relies on having a unique
mapping between initial and final densities given by local approximations like
spherical or cylindrical collapse (Uhlemann et al., 2016, 2018b). This method
has been applied to a variety of problems: to compute marginal and joint
density and velocity-divergence PDFs (Uhlemann et al., 2017), the weak lensing
convergence PDF (Reimberg & Bernardeau, 2018; Barthelemy et al., 2020) and to
constrain dark energy models (Codis et al., 2016).

Many of the above theoretical methods use spherical collapse dynamics as a
proxy for the non-linear regime. The main advantage is that the equations are
simple second-order ordinary differential equations (ODEs). Analytic solutions
are available for pure matter cosmologies and numerical solutions can be
easily generated not only for $\Lambda$CDM but also for many other non-
standard cosmological models such as those with a constant equation of state
($w\neq-1$) or those with a time-varying $w$. The exceptions are cases where
dark matter and dark energy are coupled or Einstein’s law of gravity is
modified since these systems are higher-order with additional perturbation
variables. In spite of this simplicity, most of the aforementioned
investigations that invoke spherical dynamics have been restricted to pure
matter or $\Lambda$CDM only.

In this paper, we use the spherical collapse model to compute the one-point
PDF and compare the case of an Early Dark Energy (EDE) model with the usual
$\Lambda$CDM cosmology. The physical assumption is that each spatial point
evolves independently (local approximation) as if it were an isolated
spherical perturbation in an expanding background cosmology. The stochasticity
is in the initial conditions: they are drawn from a random Gaussian
distribution. This paper has two primary aims. First, we want to check how
well the non-linear PDF from spherical collapse compares with other fitting
forms in the literature and second, we want to understand the effect of early
dark energy on structure, in particular on the one-point PDF. Nadkarni-Ghosh
(2013) also considered spherical collapse to obtain the density-velocity
divergence relation (DVDR) in the non-linear regime. However, the
investigation was restricted to dark energy models with constant equation of
state ($w=w_{0}$). A third, tertiary aim of this work is to also check if the
fitting form for the non-linear DVDR obtained in N13 also works for the EDE
model with a time-varying $w$.

The paper is organized as follows. §2 sets up the equations for non-linear
evolution of the density field for spherically symmetric perturbations
expanding in a background cosmology consisting of a dark energy component that
neither clusters nor couples to the matter component. §3 discusses the growth
rates in the linear and non-linear regime and compares the behaviour of the
EDE model with the standard $\Lambda$CDM model. Here we accomplish our third
aim of checking how well the non-linear DVDR relation obtained in N13 works
for the EDE model by comparing the theoretical and numerical non-linear growth
rates. §4 gives the details of the numerical runs performed to generate the
PDF and compares the results to fitting forms existing in the literature. We
consider three forms: the theoretical prediction based on perturbation theory
(Bernardeau, 1994a), the skewed log-normal (Colombi, 1994) and the generalized
normal distribution, version 2 (Shin et al., 2017). These estimates are in the
Eulerian frame, whereas, the PDFs obtained through spherical collapse are
Lagrangian PDFs and we discuss the transformation between the two. §5 compares
the PDFs for $\Lambda$CDM and EDE models in terms of (a) the difference of the
PDFs and the (b) evolution of the width of the PDF. We conclude in §6.

##  2 The physical set-up

###  2.1 Non-linear evolution equations in the Spherical Collapse (SC) model

The physical system is a constant density (top-hat) sphere, evolving in a
homogenous and isotropic, flat background universe consisting of pressureless
matter and dark energy with equation of state $w$, which may vary with time.
Matter is allowed to cluster whereas dark energy is assumed to be uniform, and
non-interacting with the matter sector. For such cosmologies, the scale factor
$a$ of the background evolves as

$H^{2}=H_{0}^{2}\left[\frac{\Omega_{m,0}a_{0}^{3}}{a^{3}}+\Omega_{de,0}\exp\left\\{-3\int_{a_{0}}^{a}(1+w(y))\frac{dy}{y}\right\\}\right],$
(1)

where $H={\dot{a}}/a$. The ‘dot’ denotes derivative with respect to coordinate
time, $y$ is the dummy variable used for the integral and $a_{0}$ is the value
of the scale factor today, usually set to 1. The matter and dark energy
parameters are denoted as $\Omega_{m}$ and $\Omega_{de}$ respectively and the
subscript ‘0’ denotes their values today. Using the conservation of energy
density, the density parameters at any epoch are related to their values today
as

$\displaystyle\Omega_{m}(a)$ $\displaystyle=$
$\displaystyle\frac{H_{0}^{2}a_{0}^{3}\Omega_{m,0}}{H^{2}a^{3}}$ (2)
$\displaystyle\Omega_{de}(a)$ $\displaystyle=$
$\displaystyle\frac{H_{0}^{2}}{H^{2}}\Omega_{de,0}\exp\left\\{-3\int_{a_{0}}^{a}(1+w(y))\frac{dy}{y}\right\\}$
(3)

For a flat universe, $\Omega_{m}(a)+\Omega_{de}(a)=1$ for all epochs. Let the
origin be at the centre of the spherical perturbation and let ${\bf r}$ and
${\bf x}$ be the physical and comoving radius of the outer edge. The perturbed
sphere may be surrounded by a compensating region, which is either a mass
shell or a void, depending upon whether the perturbation is underdense or
overdense with respect to the background. Let ${\bar{\rho}}_{m}$ and
$\rho_{m}$ be the density of the background and perturbation respectively. For
the dark energy models considered in this paper, an initial top-hat stays a
top-hat and the dynamics can be solely described by two purely time-dependent
quantities: fractional overdensity ($\delta$) and the scaled divergence of the
peculiar velocity ($\Theta$) defined as

$\displaystyle\delta$ $\displaystyle=$
$\displaystyle\frac{\rho_{m}}{{\bar{\rho}}_{m}}-1,$ (4) $\displaystyle\Theta$
$\displaystyle=$ $\displaystyle\frac{1}{H}\nabla_{r}\cdot{\bf v},$ (5)

where ${\bar{\rho}}_{m}$ and $\rho_{m}$ are the densities of the background
and perturbation respectively and ${\bf v}={\dot{\bf r}}-H{\bf r}=a{\dot{\bf
x}}$ is the peculiar velocity. For a pure matter cosmology, it is possible to
invoke Birkoff’s theorem and treat the perturbed sphere as a separate universe
whose evolution is dictated by an equation analogous to eq. (1), with a
different scale factor and density parameter. However, extensions to dark
energy cosmologies are not obvious. Here we adopt the approach of starting
with the hydrodynamic equations: the generalized continuity and Euler
equations along with the Poisson equation for the gravitational potential
(Lima et al., 1997; Abramo et al., 2009; Pace et al., 2010, 2017). These are
valid for general non-interacting dark energy cosmologies, in the sub-horizon
limit and in the absence of shell-crossing. The resulting equations for
$\delta$ and $\Theta$ are

$\displaystyle\delta^{\prime}$ $\displaystyle=$
$\displaystyle-(1+\delta)\Theta,$ (6) $\displaystyle\Theta^{\prime}$
$\displaystyle=$
$\displaystyle-\frac{3}{2}\Omega_{m}(a)\delta-\Theta\left(2+\frac{H^{\prime}}{H}\right)-\frac{\Theta^{2}}{3}.$
(7)

The prime denotes derivative with respect to $\ln a$. The details of the
derivation are in appendix A. Combining eq. (6) and eq. (7) gives one second
order equation for the non-linear evolution of $\delta$.

$\delta^{\prime\prime}+\left(2+\frac{H^{\prime}}{H}\right)\delta^{\prime}-\frac{4}{3}\frac{\delta^{\prime
2}}{(1+\delta)}-\frac{3}{2}\Omega_{m}(a)\delta(1+\delta)=0.$ (8)

For the flat dark energy models, with varying equation of state, eq. (8)
reduces to

$\delta^{\prime\prime}+\frac{1}{2}\left[1-3w(a)\Omega_{de}(a)\right]\delta^{\prime}-\frac{4}{3}\frac{\delta^{\prime
2}}{(1+\delta)}-\frac{3}{2}\Omega_{m}(a)\delta(1+\delta)=0.$ (9)

This system is also valid if the density has a radial variation: in that case
$\delta$ has to be replaced by its spherically averaged value (N13).

###  2.2 Initial conditions

Equation (9) can be completely solved by specifying $\delta$ and
$\delta^{\prime}$ at some initial time $a_{i}$. From eq. (6), we have
$\delta^{\prime}_{i}=-(1+\delta_{i})\Theta_{i}$. In principle, $\delta_{i}$
and $\Theta_{i}$ are independent however, usually cosmological initial
conditions assume the Zeldovich approximation (Zeldovich, 1970). This states
that the velocity field is proportional to the acceleration field; the
proportionality constant is chosen so as to ensure that there are no
perturbations at the big bang time. In linear, Eulerian perturbation theory,
for the EdS cosmology, this is equivalent to setting the coefficient of the
decaying mode to zero and the growing mode solution is $\delta(a)\propto a$,
which gives $\delta^{\prime}=\delta$ or $\Theta=-\delta$. In this paper, we
use the generalization given by Linder (2005) for dark energy cosmologies:
$\Theta_{i}=\Omega_{m,i}^{0.55}\delta_{i}$. For sufficiently early epochs and
sufficiently small $\delta_{i}$, this reduces to
$\delta^{\prime}_{i}=\delta_{i}$. We will discuss extensions of this linear
relation in §3.

###  2.3 Cosmological models

In this paper, we primarily consider two models: (a) $\Lambda$CDM with
$\Omega_{m,0}=0.29$, $\Omega_{de,0}=0.71$ and (b) Early Dark Energy (EDE
hereafter) model with varying equation of state

$w(a)=\frac{w_{0}}{(1-b\ln a)^{2}}$ (10)

with

$b=-\frac{3w_{0}}{\text{ln}\Big{(}\frac{1-\Omega_{de,e}}{\Omega_{de,e}}\Big{)}+\text{ln}\Big{(}\frac{1-\Omega_{m,0}}{\Omega_{m,0}}\Big{)}},$
(11)

where $\Omega_{de,e}$ is the value of $\Omega_{de}$ at very early time
($a\rightarrow 0$), $w_{0}$ is the value of equation of state parameter at
present day. This three-parameter model was first proposed by Wetterich (2004)
and later also considered by Doran & Robbers (2006). The EDE model contains a
low but non-vanishing amount of dark energy at early epoch. It can be
characterized by the quantity

${\bar{\Omega}_{de,sf}}=(\ln a_{tr}-\ln a_{eq})^{-1}\int_{\ln a_{eq}}^{\ln
a_{tr}}\Omega_{de}(a)d\ln a,$ (12)

where $a_{eq}$ is the epoch of matter radiation equality and $a_{tr}$ is the
epoch where dark energy starts to dominate ($a_{tr}\sim 1/3$).
${\bar{\Omega}_{de,sf}}$ denotes the average fraction of dark energy during
the structure formation/matter domination era and using a combination of data
sets it is constrained to be ${\bar{\Omega}_{de,sf}}\leq 0.04$ (Doran,
Schwindt, & Wetterich, 2001, 2007). This constraint was imposed by Bartelmann,
Doran, & Wetterich (2006) as well as Grossi & Springel (2008) who investigated
the impact of EDE on non-linear structure. They considered two EDE models with
${\bar{\Omega}_{de,sf}}=0.04$, but other parameters randomly picked from a
Monte Carlo chain. The two EDE models did not give qualitatively different
results, so, we choose one of them: $\Omega_{de,e}=8\times 10^{-4}$ and
$w_{0}=-0.99$. The value for $\Omega_{m,0}$ was chosen to be the same as the
$\Lambda$CDM model and the requirement of spatial flatness sets
$\Omega_{de,0}$. An alternative approach is to consider the best fit values
for an EDE cosmology, which in general will be different than those in
$\Lambda$CDM (Jennings et al., 2010), but we do not adopt this approach here.
We do not consider a larger range of EDE models because the main aim is to
qualitatively contrast the growth of the PDF in $\Lambda$CDM and EDE
scenarios.

Another commonly used parameterization for dark energy is the two-parameter
CPL parameterization (Chevallier & Polarski, 2001; Linder, 2003)

$w_{CPL}(a)=w_{0}+w_{a}(1-a),$ (13)

where $w_{0}$ and $w_{a}$ are two free parameters to be set by data. This
parameterization corresponds to a first order Taylor expansion of the equation
of state around $a=1$, with $w_{0}$ being the value at $a=1$. Performing a
similar expansion on eq. (10) gives

$w_{EDE}(a)=w_{0}-2bw_{0}(1-a).$ (14)

Matching eqs. (13) and (14) gives $w_{a}=-2bw_{0}$. Most recent observational
constraint from Planck Collaboration et al. (2018) is $w_{0}=-0.957$ and
$w_{a}=-0.29$ combining the Planck, SNe and BAO data. A more detailed
comparison of the two models can be found in appendix B.

Figure 1 shows the evolution of the matter density for the models considered
in this paper. For pedagogical reasons, we also include here and in the next
section the EdS (flat, pure matter) cosmology with $\Omega_{m}=1$. The left
panel shows the density parameters $\Omega_{m}$ and $\Omega_{de}$ in all three
cosmologies. The EDE cosmology has a non-negligible amount of dark energy at
early times. Assuming flatness, this amounts to a lower value of $\Omega_{m}$
at early epochs. The EdS cosmology has $\Omega_{m}=1$ at all epochs. The right
panel shows the ratio of the matter density parameters for the dark energy
cosmologies: $\Omega_{m,{\rm EDE}}/\Omega_{m,\Lambda{\rm CDM}}$. The ratio is
minimum at around $a\sim 0.4$. The difference in growth rates between the two
models is sensitive to this ratio; see figure 4.

Figure 1: Evolution of density parameters. Left panel: the matter density
parameter is unity for the EdS cosmology. The EDE model has non-negligible
dark energy at early epochs as compared to $\Lambda$CDM, but the $\Omega_{m}$
is still close to unity for both models. Right panel: the ratio of the matter
density parameters in EDE relative to the $\Lambda$CDM cosmology. The ratio is
minimum around $a\sim 0.4$. The vertical blue and red lines in both panels
indicate the epoch of matter-dark energy equality for EDE and $\Lambda$CDM
respectively. They are very close and almost indistinguishable on a log-axis.

##  3 Linear and non-linear growth rates

The non-linear evolution of $\delta$ in the SC model is given by eq. (9). The
growth rate $f$ is defined as

$f=\frac{d\ln\delta}{d\ln a}.$ (15)

This definition is valid in the linear as well as non-linear regime. In the
linear regime, $\delta\propto D(a)$, where $D(a)$ is the growth factor (see
appendix A) and $f=D^{\prime}(a)/D(a)$. Using the continuity equation (eq.
(6)), $f$ can also be expressed as

$f=-\frac{1+\delta}{\delta}\Theta.$ (16)

In non-linear regime, $f$ is a function of the instantaneous $\delta$. The
higher the $\delta$, the higher the growth rate. Knowing the relation between
the $\delta$ and $\Theta$, it is possible to compute the growth rate as a
function of density. This relation has been studied extensively by various
authors in the past, using perturbation theory (e.g., Bernardeau 1992) or
numerical simulations. Bilicki & Chodorowski (2008) gave a fit based on
analytical solutions of SC. N13, investigated the $w$ dependence of this
relation, also using SC, but without relying on the analytic solutions. The
main idea was to impose the condition that there be ‘no perturbations at the
big bang’. This is a generalization of the idea of ‘no decaying mode’, which
is the basis of the ‘Zeldovich approximation’ (Zeldovich, 1970; Peebles,
1980). This is also sometimes referred to as the ‘slaving condition’. The
assumption that the initial peculiar velocity field is ‘slaved’ to the initial
gravitational field has found various applications in the literature. For
example, Brenier et al. (2003) combined this condition with the MAK
reconstruction formalism to evolve perturbations back in time to obtain
information about the primordial initial conditions. Buchert (1992) showed
that this condition was a special case of a more general first-order
Lagrangian Perturbation Theory (LPT). Nadkarni-Ghosh & Chernoff (2011), while
investigating the convergence properties of LPT, found that this condition,
imposed in the SC model for an EdS cosmology, traces out a unique curve in the
two dimensional $\delta-\Theta$ phase space. This curve, called the ‘Zeldovich
curve’ is an invariant set of the dynamics and this property was exploited to
extend the range of validity of LPT. Nadkarni-Ghosh (2013)222In N13,
$\delta_{v}$ was used instead of $\Theta$; $\Theta=3\delta_{v}$. extended the
phase space analysis to include constant equation of state cosmologies and
showed that the invariant property of the ‘Zeldovich curve’ still holds -
perturbations that start with the slaving condition, maintain the condition
and those that do not start out slaved, eventually get slaved. Thus, this
condition acts like an attractor of the dynamical system defined by eq. (6)
and (7) and it corresponds to the late time (non-linear) density velocity-
divergence relation (DVDR). The fitting form obtained in N13 based constant
$w$ cosmologies is:

$\Theta_{fit}(\Omega_{m},w,\delta)=\left\\{\begin{array}[]{cc}3A(\Omega_{m},w)\left(1-(1+\delta)^{B(\Omega_{m},w)}\right),&-1\leq\delta\leq
1\\\ \\\
3\Omega_{m}^{\gamma_{1}(w)+\gamma_{2}(w)}\left[(1+\delta)^{1/6}-(1+\delta)^{1/2}\right],&1\leq\delta\leq
10.\end{array}\right.$ (17)

This fit is an extension of the forms given by Bernardeau (1992) and Bilicki &
Chodorowski (2008). Substituting from eq. (17) gives

$f_{fit}(\Omega_{m},w,\delta)=\left\\{\begin{array}[]{cc}-3A(\Omega_{m},w)\delta^{-1}\left((1+\delta)-(1+\delta)^{B(\Omega_{m},w)+1}\right),&-1\leq\delta\leq
1\\\ \\\
-3\Omega_{m}^{\gamma_{1}(w)+\gamma_{2}(w)}\delta^{-1}\left[(1+\delta)^{7/6}-(1+\delta)^{3/2}\right],&1\leq\delta\leq
10\end{array}\right.$ (18)

where

$A(\Omega_{m},w)=\frac{1}{2}\Omega_{m}^{\gamma_{1}(w)},\;B(\Omega_{m},w)=\frac{2}{3}\Omega_{m}^{\gamma_{2}(w)},\;\gamma_{1}(w)=0.56(-w)^{-0.08}\;\;{\rm
and}\;\;\gamma_{2}(w)=-0.01(-w)^{-1.18}.$ (19)

and $\Omega_{m}$ is evaluated at the epoch of interest. In the linear regime,
the fitting function reduces to

$f_{lin,fit}(\Omega_{m},w)=\Omega_{m}^{\gamma_{1}(w)+\gamma_{2}(w)},$ (20)

which for the $\Lambda$CDM case, agrees with the commonly used
parameterization of Linder (2005) $f_{Linder}\approx\Omega_{m}^{0.55}$.

One of the aims of this paper is to check the validity of eq. (18) for the EDE
model. To do this, we compare this form given by eq. (18) to the numerically
computed growth factor obtained by solving eq. (9) for Zeldovich initial
conditions. We consider the following cases:

  * •

Case (1a): Overdensity with same initial amplitude for all three cosmologies
[$\delta(a=0.001)=1.52\times 10^{-3}$].

  * •

Case (1b): Void with same initial amplitude for all three cosmologies
[$\delta(a=0.001)=-1.52\times 10^{-3}$].

  * •

Case (2a): Overdensity with same linear amplitude today for all three
cosmologies [$\delta(a=1)=1$].

  * •

Case (2b): Void with same "linear" amplitude today for all three cosmologies
[$\delta(a=1)=-1$].

For cases (2a) and (2b), the initial values at $a_{i}=0.001$ are obtained by
multiplying by the linear growth factor normalized to unity at $a=1$. This
factor is denoted as ${\tilde{D}}(a)$ and expression for a general cosmology
is in appendix A. This gives $\delta_{i,EDE}=\pm 1.52\times
10^{-3},\delta_{i,\Lambda CDM}=\pm 1.29\times 10^{-3},\delta_{i,EdS}=\pm
0.001$. We note that the limit $\delta_{L}=-1$ does not correspond to the
linear regime, but the generated initial amplitudes are well defined. The
rationale behind considering two classes of initial conditions is the
following. It is clear that the presence of early dark energy will reduce the
growth rate of structure. The simplest way to illustrate this is to consider
the same initial conditions for all three cosmologies. This corresponds to
cases (1a) and (1b). However, observations constrain the linear power spectrum
today i.e., at $a=1$. Thus, it is important to understand the differences in
non-linear growth, subject to this constraint. This corresponds to cases (2a)
and (2b). By construction, the EDE solution is the same for classes (1) and
(2).

Figure 2: Solutions for $\delta$ as a function of scale factor for the four
cases defined in the text. (1a) and (1b) have same initial conditions for all
three cosmologies whereas, (2a) and (2b) have same conditions at $a=1$. The
dependence on cosmology is clear in cases (1a) and (1b) whereas, in cases (2a)
and (2b), a lower $\Omega_{m}$ for the EDE model gets compensated by a higher
initial amplitude and the deviations due to cosmology are less obvious. The
dark energy models in both (2a) and (2b) are more non-linear than the pure
matter cosmology. Note also that for voids, the non-linear solutions grows
slower than in the linear regime. As a check, it is useful to note that in the
pure matter (EdS) cosmology, $\delta\propto a$ in the linear regime.

Figure 2 shows the evolution of $\delta$ as a function of $a$ for the four
cases. The following features can be noted. First, consider the overdensities.
The differences between cosmologies are very clear in case (1a), whereas, they
are much lower in (2a); in fact, the solutions for the $\Lambda$CDM and EDE
cosmologies, almost coincide. This is because in case (2a), the lower
$\Omega_{m}$ in the EDE cosmology is compensated with a higher initial
amplitude. The effect is even more pronounced when comparing EdS and dark
energy cases. The initial amplitude is low enough that the value of $\delta$
is always less than the dark energy cases; however, the growth rate is higher
so the solution evolves to the same $\delta$ value at $a=1$. For voids, we
first note that the linear evolution is ‘faster’ than the non-linear
evolution. This is because a void is underdense as compared to the background
and hence grows slower. Mathematically, in the linear regime,
$|\delta^{\prime}|\propto|\Theta|$ whereas, in the non-linear regime
$|\delta^{\prime}|\propto|(1+\delta)\Theta|$. Since $(1+\delta)<1$, the non-
linear growth is slower than the linear growth. Comparing cases (1a) and (1b),
we note a similar effect as seen in overdensities. In case (1b), the lower
value of $\Omega_{m}$ gives a slower evolution in EDE cosmology as compared to
$\Lambda$CDM. This difference is reduced in case (2b) because the lower value
of $\Omega_{m}$ is compensated by a higher initial amplitude (absolute value)
in the EDE.

Figure 3:  The linear and non-linear growth rates for the four cases discussed
in the text. The solid lines indicate the fit given by eq. (18). The EdS model
is plotted for reference. The linear growth rate varies with cosmology, but is
independent of the initial conditions. The growth rates for EDE and
$\Lambda$CDM match at $a=1$ because the two models have the same $\Omega_{m}$
today. In overdensities, in case (1a), the difference in growth rates is
purely due to the difference in $\Omega_{m}$ in the three cosmologies. The
cosmology dependence between $\Lambda$CDM and EDE is damped in case (2a)
because the EDE cosmology has a higher initial amplitude which compensates the
effect of a lower $\Omega_{m}$. Voids, on the other hand, show a different
behaviour. A higher $\Omega_{m}$ means that the linear growth rate is higher
and a void goes non-linear earlier. However, unlike overdensities, the growth
rate slows down in the non-linear regime
($f\propto\Omega_{m}^{\gamma_{1}}(1+\delta)$). Thus, the growth rate is ‘self-
regulatory’ in voids. This explains the small differences between the three
cosmologies in case (1b). In case (2b), the initial amplitude for the
$\Lambda$CDM (and EdS) is less negative than the corresponding value in case
(1b). Thus, at initial epochs, the differences in cosmology are more
pronounced, but at late epochs, the differences reduce due to the ‘self-
regulation’.

Figure 3 shows the growth rates for the above four cases. We note the
following features. The linear growth rate $d\ln\delta/d\ln a$ does not depend
on the initial amplitude but depends only on the cosmological parameters.
Thus, it remains unchanged in all four panels. On the other hand, the non-
linear growth rate depends on the instantaneous $\delta$, which in turn
depends on the initial amplitude as well as the cosmological parameters. For
overdensities, the difference between EDE and $\Lambda$CDM is more pronounced
when the constraint of ‘same linear $\delta$ today’ is not imposed. This is
similar to the effect seen in figure 2 for the same reason.

However, the case for voids is somewhat counter-intuitive when contrasted with
overdensities and figure 2. The growth rate for voids does not show stark
differences in behaviour when comparing cases (1b) and (2b). This can be
understood as follows. As noted earlier, the more underdense the void, the
slower it grows. From eq. (18), the void limit ($\delta\rightarrow-1$) gives,
$f_{void}\propto\Omega_{m}^{\gamma_{1}}(1+\delta)$. Thus, for the same initial
conditions, a higher $\Omega_{m}$ implies a faster growth which implies a
higher amplitude underdensity which in turn reduces the growth rate i.e., the
two factors act in opposite directions. This ‘self-regulation’ does not occur
in overdensities:
$f_{overdense}\propto\Omega_{m}^{\gamma_{1}+\gamma_{2}}\sqrt{\delta}$, and the
two factors act in tandem. Thus, for the cases of ‘same initial conditions’
(cases (1a) and (1b)), the differences in cosmology are apparent for
overdensities but not for voids. In case (2b), the initial amplitude for
$\Lambda$CDM and EdS is smaller (less negative) than their values in case
(1b). The growth rate is therefore initially larger making the effect of
cosmology slightly more pronounced at early epochs, but the difference is
reduced at later epochs when the rate slows down as the underdensity grows.
Overall, we conclude that cosmology dependence of the growth factor in the
overdense case is more sensitive to how initial conditions are set as compared
to the void case.

The solid lines in figure 3 indicate the fits given by eq. (18), where
$\delta$ is evaluated using the numerical solution. We found that over the
full range of epochs, the maximum relative error between the numerical value
and the value given by the fit was less than 4 %, indicating that the form
given by eq. (18) also works well for the EDE model considered here. From the
fitting functions, it is clear that the dependence on $w$ is very weak, even
in the non-linear regime. In the linear regime, it is usually neglected to
give $f_{lin}\sim\Omega_{m}^{0.55}$ for most dark energy cosmologies with
standard gravity. We found that, for the EDE cosmology, ignoring the $w$
dependence (set $w=-1$) gave only an error of 0.8% (0.4%) in the non-linear
(linear) growth rate. This indicates that the difference in growth rates
between $\Lambda$CDM and EDE cosmology is primarily from the $\Omega_{m}$
dependence. Figure 4 shows the ratios of the linear growth rates in the EDE to
the $\Lambda$CDM cosmology. The black vertical line indicates the position
where the ratio $\Omega_{m,EDE}/\Omega_{m,\Lambda CDM}$ is minimum i.e.,
$a\sim 0.4$. The ratio of the linear growth rates also reaches its minimum
close to $a\sim 0.4$.

Figure 4: Ratio of the linear growth rates. This ratio is independent of the
initial conditions since it depends only on the cosmological parameters. The
ratio is minimum near $a_{min}\sim 0.4$, which corresponds roughly to the
epoch where the ratio of matter density parameters is minimum; the slight
offset indicates a weak dependence on the equation of state $w$.

##  4 One-point PDF of the density field

###  4.1 Numerical Runs

The numerical data was generated as follows. An ensemble of points was drawn
from a Gaussian random field with mean $0$ and standard deviation
$\sigma_{G}$. This corresponds to the linear distribution at $a=1$. The
distribution at any other epoch is obtained by scaling each data point by the
normalized linear growth factor, ${\tilde{D}}(a)$, given in appendix A. The
resulting distribution is Gaussian with the linearly scaled width

$\sigma_{L}(a)=\sigma_{G}{\tilde{D}}(a).$ (21)

The realizations at $a=1$ are the same for all three cosmologies, but differ
at the initial epoch due to differences in the growth factors. This set-up is
similar to cases (2a,b) in the previous section. We chose the initial epoch at
$a=0.001$. We considered three values for $\sigma_{G}:0.5,\,1$ and $2$. For a
BBKS power spectrum with $n_{s}=1$ and $\sigma_{8}=0.9$, this corresponds
roughly to 15, 7 and 3 $h^{-1}$ Mpc. Each initial condition of the ensemble
was evolved using eq. (9) and the final density was read out at $a=0.1,\,0.4$
and $1$ corresponding to redshift $z=9,1.5$ and $0$ respectively. This
procedure gives an ensemble of final densities, whose PDF is constructed by
first imposing a cut-off $\delta\leq 5$ and then binning with 40 bins
equispaced on the log-scale. The PDF was defined as the fractional number of
points in each bin divided by the length of the bin. Thus, the resulting PDF
is normalized by construction. For each value of $\sigma_{G}$, five
realizations each containing $5\times 10^{4}$ points were evolved. The final
PDF plotted below is the average over these five realizations and error bars
correspond to the standard deviation.

The PDF obtained through this procedure assumes that the mass of each sphere
is unchanged throughout the evolution and hence gives a Lagrangian estimate of
the PDF. The numerical forms given below have been obtained using Eulerian
perturbation theory or are based on simulations and provide the Eulerian PDF.
The transformation between the Lagrangian to Eulerian PDF is outlined in
appendix C. However, it is important to note that we are not comparing our
results to data from N-body codes. In view of this, it is possible to treat
the fitting forms given below as purely analytic expressions whose parameters
are estimated from our data generated using spherical dynamics. One can then
also compare both the Lagrangian and Eulerian PDFs to these forms
independently. The exception is B94, which is explicitly Eulerian and does not
have any free parameters. It needs to be converted appropriately to facilitate
comparison.

###  4.2 Fitting forms for comparison

In this paper, we have considered three fitting forms (1) the distribution
based on Eulerian perturbation theory given by Bernardeau (1994a); hereafter
referred to as the B94 distribution, (2) the skewed log-normal distribution
(SLN) which combines the Edgeworth series and the log-normal form (Colombi,
1994) and (3) the generalized normal distribution ($\mathcal{N}_{v2}$) based
on numerical simulations given by (Shin et al., 2017). We compare the analytic
expressions to the PDF generated from our numerical runs.

#### Estimates based on Eulerian Perturbation Theory

The B94 is a piecewise distribution for density field and given by,

$P(\delta)=\left\\{\begin{array}[]{cc}P_{B94}^{I}(\delta)&-1\leq\delta<0.3\\\
\\\ P_{B94}^{II}(\delta)&\delta>0.3,\end{array}\right.$ (22)

where

$P_{B94}^{I}(\delta)d\delta=\Big{(}\frac{7-5(1+\delta)^{\frac{2}{3}}}{4\pi\sigma_{\delta}^{2}}\Big{)}^{\frac{1}{2}}(1+\delta)^{-\frac{5}{3}}\,\text{exp}\Big{[}-\frac{9}{8\sigma_{\delta}^{2}}\Big{(}-1+\frac{1}{(1+\delta)^{\frac{2}{3}}}\Big{)}^{2}\Big{]}d\delta$
(23)
$P_{B94}^{II}(\delta)d\delta=f_{c}\frac{3a_{s\delta}\sigma_{\delta}}{4\sqrt{\pi}}(1+\delta)^{-\frac{5}{2}}\,\text{exp}\Big{[}\frac{-|y_{s\delta}|\delta+|\phi_{s\delta}|}{\sigma_{\delta}^{2}}\Big{]}d\delta,$
(24)

where $\sigma_{\delta}=\sigma_{L}$, where $\sigma_{L}$ is given by eq. (21),
$a_{s\delta}=1.84$, $y_{s\delta}=-0.184$, $\phi_{s\delta}=-0.03$. We have
chosen the $n=-3$ values for the parameters $a_{s}$, $y_{s}$, $\phi_{s}$
(B94). The correction factor $f_{c}$ is

$f_{c}=[1+2(0.8-\sigma_{\delta})\sigma_{\delta}^{-1.3}(1+\delta)^{-0.5}]$ (25)

accounts for the fact that the PDF based on perturbation theory does not
perform well at high $\delta$.

The B94 distribution is derived in the Eulerian frame. Converting to its
Lagrangian form involves inverting the Lagrangian to Eulerian mapping given by
eq. (61) in the appendix. This is not straightforward and as a simplification
we use the approximation $P_{L}(\delta)=P_{E}(\delta)(1+\delta)$ (Nadkarni-
Ghosh & Singhal, 2016).

#### Skewed log-normal distribution (SLN)

The skewed log-normal PDF ($P_{SLN}$) proposed by Colombi (1994) combines the
log-normal PDF and the Edgeworth series. Define $\nu\equiv\Phi/\sigma_{\Phi}$,
with $\Phi\equiv\ln(1+\delta)-\langle\ln(1+\delta)\rangle$,
$\sigma_{\Phi}\equiv\sqrt{\langle\Phi^{2}\rangle}$, where the brackets denote
the average over the $\phi$ values. $P_{SLN}$, which uses the third order
Edgeworth expansion, is defined as

$P_{SLN}(\nu)d\nu=\Big{[}1+\frac{1}{3!}T_{3}\sigma_{\Phi}H_{3}(\nu)+\frac{1}{4!}T_{4}\sigma_{\Phi}^{2}H_{4}(\nu)+\frac{10}{6!}T_{3}^{2}\sigma_{\Phi}^{2}H_{6}(\nu)\Big{]}\mathcal{N}(\nu),$
(26)

where

$T_{3}\equiv\frac{\langle\Phi^{3}\rangle}{\sigma_{\Phi}^{4}},$ (27)
$T_{4}\equiv\frac{\langle\Phi^{4}\rangle-3\sigma_{\Phi}^{4}}{\sigma_{\Phi}^{6}},$
(28)
$\mathcal{N}(\nu)=\frac{1}{\sqrt{2\pi}}\exp\Big{[}-\frac{\nu^{2}}{2}\Big{]}$
(29)

and $H_{m}(\nu)$ is the Hermite polynomial of degree $m$. This distribution is
not always positive definite, i.e., it can be negative valued hence is not a
real PDF (Colombi, 1994). But, N-body simulations show that this approximation
successfully describes the late-time evolution of the dark matter field (Ueda
& Yokoyama, 1996; Szapudi & Pan, 2004).

#### Generalised normal distribution of version 2 ($\mathcal{N}_{v2}$)

This distribution has been recently proposed by Shin et al. (2017) based on
numerical simulations of five different cosmological models. The
$\mathcal{N}_{v2}$ distribution consists of three parameters namely, $\alpha$
(scale), $\kappa$ (shape) and $\xi$ (location). $\alpha$ is related to the
width of the distribution, $\kappa$ quantifies skewness and $\xi$ gives the
location of the peak. In terms of these parameters, the $\mathcal{N}_{v2}$ is
defined as,

$\mathcal{N}_{v2}(\delta)=\frac{\mathcal{N}(y)}{\alpha-\kappa(\delta-\xi)},$
(30)

where $N(y)$ is the usual Gaussian distribution with mean 0 and variance 1.
Here $y$ is defined as,

$y=\left\\{\begin{array}[]{cc}-\frac{1}{\kappa}\ln\Big{[}1-\frac{\kappa(\delta-\xi)}{\alpha}\Big{]},&\text{
if $\kappa\neq 0$}\\\ \\\ \frac{\delta-\xi}{\alpha},&\text{ if
$\kappa=0$},\end{array}\right.$ (31)

A positive (or negative) value of $\kappa$ yields a left-skewed (or right-
skewed) distribution bounded to the right (or left). The mean ($\mu$), median
(${\tilde{\delta}}$) and skewness ($\gamma_{1}$) of the distribution are
related to $\alpha$, $\kappa$ and $\xi$ as

$\displaystyle\mu$ $\displaystyle=$
$\displaystyle\langle\delta\rangle=\xi-\frac{\alpha}{\kappa}\left(e^{\kappa^{2}/2}-1\right)$
(32) $\displaystyle{\tilde{\delta}}$ $\displaystyle=$ $\displaystyle\xi$ (33)
$\displaystyle\gamma_{1}$ $\displaystyle=$
$\displaystyle\frac{\langle(\delta-\mu)^{3}\rangle}{(\langle(\delta-\mu)^{2}\rangle)^{\frac{3}{2}}}=\frac{3e^{\kappa^{2}}-e^{3\kappa^{2}}+2}{\left(e^{\kappa^{2}}-1\right)^{\frac{3}{2}}}\rm{sign}(\kappa).$
(34)

Here the ‘$\langle\cdot\rangle$’ denotes average over the domain of $\delta$.
$\mathcal{N}_{v2}$ approaches the usual Gaussian normal distribution in the
limit that $\kappa$ tends to zero. In this case, $\mu=\xi$ and $\gamma_{1}=0$.
There are two ways to find the parameters $\alpha$, $\kappa$ and $\xi$. One
way is to compute $\mu,{\tilde{x}}$ and $\gamma_{1}$ from the data and solve
eqs. (32), (33) and (34) for the three parameters. The other way is to find
the best fit values by optimizing the fit to the data. We implemented the
latter method using Python’s SciPy inbuilt function. This method requires an
initial guess for $\alpha$, $\kappa$ and $\xi$. We found that the errors were
minimum when the initial guess values were computed analytically using the
first method. Table 1 in appendix D shows the values of parameters for the
best fit.

###  4.3 Results and discussion

Figure 5: Variation of the one-point Lagrangian space PDF with epoch.
$P(\delta)$ vs. $(1+\delta)$ for $\sigma_{G}=0.5$ for the EDE and $\Lambda$CDM
models at $a=0.1$, $0.4$, $1$ ($z=9$, $1.5$, $0$). Simulated PDFs (filled red
circles) are compared with the predicted PDFs based on theory and N-body
simulations. At $a=0.1$, the Gaussian shape is retained, but non-Gaussian
features appear at later epochs.

Figure 6: Dependence of one-point Lagrangian space PDF on $\sigma_{G}$.
$P(\delta)$ vs. $(1+\delta)$ at $a=0.4$ ($z=1.5$) for the EDE and $\Lambda$CDM
models for three values of $\sigma_{G}=0.5,1,2$. Simulated PDFs (filled red
circles) are compared with the predicted PDFs based on theory and N-body
simulations. Higher $\sigma_{G}$ distribution is more non-Gaussian at the same
epoch.

We compute the numerical Lagrangian ($P_{L}$) and Eulerian ($P_{E}$) one-point
PDFs as described in §4.1 and compare with the fitting forms described above.
The shape of the PDF evolves with time and depends upon $\sigma_{G}$. Figure 5
shows the variation of the PDFs with epoch for $\sigma_{G}=0.5$ and figure 6
shows the dependence on $\sigma_{G}$ for a fixed epoch ($a=0.4$). The top two
panels of each figure show $P_{E}$ and the bottom two show $P_{L}$. From
figure 5, it is clear that at $a=0.1$, the evolution is still approximately
linear and the width is given by $\sigma_{L}$ defined in eq. (21). Non-
Gaussianity sets in at later epochs. The width of the PDF increases with
increasing epoch for a given $\sigma_{G}$. At $a=1$, in the non-linear regime,
the SLN distribution provides the best fit to $P_{E}$ in both the overdense
and void regions. SLN also fits $P_{L}$ in the overdense region, but does not
match $P_{L}$ in voids. From figure 6, it is clear that the SLN form continues
to be a good fit to $P_{E}$ for higher values of $\sigma_{G}$ in both
overdense and void regions. As was the case in figure 5, SLN fits $P_{L}$ in
overdensities also for higher values of $\sigma_{G}$ but fails in voids. The
B94 form fits $P_{L}$ for small values of $\sigma_{G}$, over the entire
domain, but fails for larger values of $\sigma_{G}$. This breakdown of B94 is
expected since Eulerian PT and possibly the approximation for $P_{E}$ to
$P_{L}$ conversion both break down in the non-linear regime. Overall, we
conclude that the SLN form fits the Eulerian PDF best over the entire range of
mass scales and redshifts. The Lagrangian PDF matches B94 for small values of
$\sigma_{G}$ over the entire redshift range, but for larger values of
$\sigma_{G}$ no particular form fits the Lagrangian PDF well.

The SC is a local model, i.e., each perturbation is assumed to evolve
independently. It does not account for tidal forces which cause shear,
rotation and accretion. It also ignores non-linear mode coupling, which is a
non-local effect. Yet, it recovers the shape of the PDF given by perturbative
methods and simulations reasonably well, over a wide range of epochs, scales
and cosmologies. The SC model has been successful in other contexts as well.
The non-linear relation between the density and velocity field obtained by
simulations is a scatter plot, but the ‘mean’ relation is well described by
the SC model (Bernardeau et al., 1999; Kudlicki et al., 2000; Bilicki &
Chodorowski, 2008). It has also been successfully used to treat the shell
crossing regime in Lagrangian perturbation theory Kitaura & Hess (2013). Some
of this success can be attributed to a partial cancellation between various
effects that have been neglected. One illustration of this feature is to
consider the effect of rotation and shear. Both terms appear in the
Raychaudhuri equation, which dictates the evolution of $\Theta$ and
consequently the full non-linear evolution of $\delta$. But they appear with
opposite signs (Peebles 1980). Spherical symmetry ignores both terms, giving
rise to a partial cancellation. For example, in modelling the non-linear DVDR,
the result of simulations agree better with spherical collapse rather than
ellipsoidal collapse (Nadkarni-Ghosh & Singhal, 2016). This is because
ellipsoidal dynamics accounts for only for the shear and does not benefit from
the partial cancellation. These arguments are akin to reasons ‘why the Press-
Schecter formalism works so well’ (Monaco 1999).

##  5 The PDF: $\Lambda$CDM vs EDE

Many past investigations have looked at the effect of an early dark energy
component on the non-linear growth of structure, particularly for the type of
EDE model considered here. Bartelmann et al. (2006) considered the SC model
and computed the critical density for collapse $\delta_{c}$, an important
ingredient in excursion set based estimates of the mass function. They found
that $\delta_{c}$ is lower for EDE models and hence for the same observed
structure today, the structure must grow earlier in EDE models. They also
noticed that, while there were more clusters at an earlier epoch in EDE than
$\Lambda$CDM, the growth rate of structure from redshift 1 to 0 was slower in
EDE than $\Lambda$CDM. Grossi & Springel (2008) performed numerical
simulations to test the conclusions based on semi-analytic estimates. They
found that no modification was required to the standard mass function
formalism and the sensitivity of the halo mass function to EDE was much less
than that predicted based on analytic estimates. Similar conclusions were
reached by Francis et al. (2008). Fontanot et al. (2012) coupled dark matter
simulations with semi-analytic models of gas physics to model galaxy formation
in EDE scenarios and again found that galaxy properties could only weakly
distinguish EDE from $\Lambda$CDM. More recently, Shi & Baugh (2016) performed
more consistent simulations, which accounted for the fact that the EDE
parameters that best fit the Planck data are different from those for
$\Lambda$CDM. They found that the cleanest signature of dark energy was on the
shape of the power spectrum and this difference was elucidated only by
adopting consistent linear theory initial conditions. In this paper, we
examine the possibility of exploiting the PDF to discriminate between the EDE
and $\Lambda$CDM cosmologies. We consider two measures based on the PDF: the
difference of the PDFs, particularly in the bulk and the width of the PDF.

###  5.1 Difference of the PDFs

Figure 7: Differences between the Eulerian PDF of $\Lambda$CDM and EDE. The
grey area denotes the ‘bulk’ region, which is most sensitive to the
differences as compared to the tails. The differences near $\delta\sim 0$ are
positive at early epochs and reverse sign as evolution progresses. The peak
magnitude is higher as non-linear evolution proceeds, but non-monotonic due to
the constraint that both cosmologies have the same linear realization at
$a=1$.

Figure 7 shows the difference of the Eulerian PDF between the $\Lambda$CDM and
EDE cosmology for $\sigma_{G}=0.5,1,2$. The grey area shows the bulk which is
defined as the region around maximum probability which includes 87% of the
area under the curve at $a=1$333The bulk was defined with respect to the
$\Lambda$CDM PDF. We found that data points in the void region which were 10%
less than the peak had huge error bars. Therefore, in defining the bulk, we
cut-off the data in the void region which lies below 10% of the peak PDF. . We
use this definition for all values of $\sigma_{G}$. It is seen that, at an
early redshift, for small values of $\delta$, the difference is positive i.e.,
the PDF for $\Lambda$CDM is greater than that of EDE and the sign changes with
increasing epoch. This can be understood as follows. The initial conditions
were chosen such that the realization was the same for the two cosmologies at
$a=1$; the values at the starting epoch ($a=0.001$), were chosen by
multiplying by the appropriate linear growth factor. Hence, for any point in
the realization, the initial value for $\Lambda$CDM is smaller than the
corresponding value for EDE. Thus, given the same PDF at $a=1$, the
distribution for $\Lambda$CDM at $a=0.001$ will have more points in the
$\delta$ bin near zero and hence the difference is positive at the initial
time and early epochs. The growth rate is higher for $\Lambda$CDM than EDE and
as evolution proceeds, the difference in PDFs for small values of $\delta$
reduces and eventually becomes negative (fewer points in a small $\delta$ bin
for $\Lambda$CDM).

At early epochs, the differences peak (and dip) near $\delta\sim 0$, but the
peak moves left as non-linear evolution proceeds. This means the differences
between cosmologies are larger in the voids than in the overdense regions.
This observation is somewhat counter-intuitive when compared to the
discussions in §3, where we conclude that in the non-linear regime, the growth
factor in voids is only mildly sensitive to cosmology. Furthermore, for the
same fixed linear density contrast today, there was no stark difference
between the void and overdense cases. This paradox can be resolved by
understanding that, although the ‘self-regulatory’ behaviour holds for an
individual data point, the cut-off $\delta\geq-1$ means that when an initially
symmetric distribution evolves to the non-linear regime, the distribution is
no longer symmetric around the origin. The bins in the void region have more
data points than those in overdense regions (the overdense points are more
spread out) and hence small differences get amplified more in voids. As the
evolution approaches $a\sim 1$, the amplitude of the peak gets higher. At a
given red-shift, the peak magnitude is higher for a higher $\sigma_{G}$
(smaller scale) since there are more points which are non-linear where the
differences between EDE and $\Lambda$CDM are more pronounced.

Recently, Shin et al. (2017) also examined the differences of the matter
density PDF for four different cosmologies (see their figure 5) and we observe
similar magnitudes for the differences at $a=1$. The differences are observed
to be larger in voids and the trends with decreasing scale (increasing
$\sigma_{G}$) are also similar. However, in their case, the magnitudes
increase with redshift monotonically, which differs from the behaviour
observed here. This is because their method to set up initial conditions
differs from ours. The use of the bulk of the PDF to constrain cosmological
parameters has also been recently advocated by Friedrich et al. (2019) in the
context of constraining primordial non-Gaussianity. Our definition of the bulk
is motivated by this paper and matches it for $R=15h^{-1}$ Mpc. Since we are
not using our method to constrain real parameters, we do not rigorously
examine the effects of this choice.

###  5.2 Width of the PDF

Figure 8: Width of the PDF (top panel) and scaled width (bottom panel) as a
function of time. At very early epochs, the growth is linear and the scaled
width is close to 1. For $\sigma_{G}=0.2$, the evolution stays close to linear
with the $\sigma_{NL}$ deviating from the linear value by only 10%. For
$\sigma_{G}=0.5$, the growth is faster at early epochs but slows down the dark
energy terms starts to dominate.

The width of the PDF is the standard deviation defined as

$\sigma_{NL}=\sqrt{\langle(\delta-\langle\delta\rangle)^{2}\rangle}.$ (35)

The cut-off in the data was imposed at $\delta\leq 20$. In order to highlight
the non-linear evolution, we also compute the scaled standard deviation:
$\sigma_{NL}$ divided by $\sigma_{L}$, where $\sigma_{L}$ is defined in eq.
(21). Figure 8 shows $\sigma_{NL}$ (top) and its scaled value (bottom) for two
values of $\sigma_{G}=0.2$ and $\sigma_{G}=0.5$ (the results for
$\sigma_{G}=1,2$ are subject to the $\delta$ cut-off and are shown in appendix
E). We evaluate the width using the Lagrangian PDF at ten epochs between
$a=0.1$ and $a=1$ (the Eulerian PDF showed similar trends). The top panel
shows that at early epochs, $\sigma_{NL}$ grows linearly with $a$, for both
scales, as expected, since the universe is matter dominated. For
$\sigma_{G}=0.2$, the fluctuations stay mostly linear, since $\sigma_{NL}$
deviates only about 10% from its linear value even at late epochs. For
$\sigma_{G}=0.5$, the fluctuations become non-linear at late epochs. The
growth rate (slope) increases from its linear value from $a\sim 0.4$ to $a\sim
0.8$, as is expected for non-linear evolution and it slows down after $a\sim
0.8$ due to the presence of the dark energy terms. For both scales considered
here, the EDE and $\Lambda$CDM cases have the same qualitative trends. The
fractional change between the two cases ($\sigma_{NL,{\rm
EDE}}/\sigma_{NL,\Lambda CDM}-1$) at $a=1$ is about 6% and 10% for the
$\sigma_{G}=0.2$ and $\sigma_{G}=0.5$ respectively.

Constraining cosmological parameters is not always straightforward. In
general, more than one cosmological parameters can have the same effect on a
given observable. These degeneracies imply that parameters cannot always be
constrained uniquely. For example, at a single red-shift the effect of EDE can
also be mimicked by a lower $\Omega_{m,0}$. Other parameters that can be
degenerate with the EDE equation of state parameters are the tilt of the power
spectrum $n_{s}$ or its normalization $\sigma_{8}$. One interesting way to
break this degeneracy is to use phase space analysis. It was shown in N13, the
simulated density and velocity profiles of galaxies were degenerate with
respect to $\sigma_{8}$ and $w$. However, plotting the same data points on a
$\delta-\Theta$ phase space plot, could break this degeneracy since the late
time non-linear relation between $\delta-\Theta$ depends only on the
parameters that dictate evolution: i.e., the $\Omega$s and $w$. It is
independent of the primordial parameters such as $n_{s}$ and $\sigma_{8}$
because the relation is an attracting set of the dynamics. Other degeneracies
may arise due to new physics - the presence of massive neutrinos or
modifications to gravity and this may complicate parameter estimation.
However, investigations of the PDF that include neutrinos suggest that the
signature of clustering neutrinos on the PDF is not generate with changes in
$\sigma_{8}$ (Uhlemann et al., 2020). In the current framework, the initial
conditions are characterized only by the initial scale ($\sigma_{G}$) and it
has no free parameters other than the fractional densities and equation of
state. Thus, an in-depth analysis of these degeneracies in the context of EDE
cosmology is beyond the scope of this work.

##  6 Conclusion

In this paper, we computed the one-point PDF of the dark matter density field
using spherical collapse as a proxy for the non-linear dynamics. We focussed
primarily on comparing the standard $\Lambda$CDM model with a three-parameter
Early Dark Energy (EDE) model (Wetterich, 2004; Doran & Robbers, 2006). The
initial distribution was assumed to be Gaussian characterized by its width
$\sigma_{G}$. The non-linear dynamics of any point in the initial realization
is governed by eq. (9) and the PDF at later epochs is numerically computed
from knowing the individual densities at the desired final time. Since the
mass of the sphere is assumed to be a constant of motion, the estimated
distribution corresponds to the Lagrangian PDF and it has to be appropriately
transformed to the Eulerian PDF (appendix C). This procedure is purely local;
it does not account for any non-local interactions between spheres. There were
three main aims: (1) to test the efficacy of this method by comparing with
existing closed forms available in the literature, (2) to compare the PDFs for
EDE and $\Lambda$CDM models and (3) to examine if the density-velocity
divergence relation (DVDR) derived by Nadkarni-Ghosh in N13 for constant $w$
cosmologies can also be applied to the EDE model with time-varying $w$. We
considered three values of $\sigma_{G}:0.5,1$ and $2$, which roughly
correspond to a smoothing scale of $R_{f}=15,7$ and 3 $h^{-1}$ Mpc for the
BBKS power spectrum and may vary slightly for other choices444The algorithm
does not require a specification of the power spectrum.. To compare the PDFs,
we considered three fitting functions: the theoretical estimate given by
Bernardeau, referred to as B94 (Bernardeau, 1994a), the skewed lognormal (SLN)
form given by Colombi (1994) and the generalized normal distribution version 2
($\mathcal{N}_{v2}$) recently proposed by Shin et al. 2017, based on numerical
simulations. There are other forms based on large deviation statistics
(Uhlemann et al., 2016) or excursion sets (Lam & Sheth, 2008b) but we do not
compare to these forms for want of a closed form fitting function. To compare
the EDE and $\Lambda$CDM models, we compute the PDF differences in the bulk
and the PDF width as two measures. The main conclusions are:

  * •

The SLN form matches the Eulerian PDF for the entire range of redshifts (from
early epochs to today), scales and cosmologies, in both the overdense and void
regions. It also provides a good fit to the Lagrangian PDF in the overdense
regions but fails in voids. The B94 form (transformed to the Lagrangian frame)
matches the Lagrangian PDF for small values of $\sigma_{G}$ in both
overdensities and voids, but for larger values of $\sigma_{G}$, no form
provides a good fit to the entire Lagrangian PDF.

  * •

The differences between the PDF, were examined in the bulk of the PDF, which
contains about 87% of the probability. The peak differences did not increase
monotonically with redshift because the initial conditions were chosen such
that the two cosmologies have the same linear distribution at $a=1$. Imposing
this constraint damps the effect of the differences in non-linear growth rates
in the two models. The peak height also increases with increasing
$\sigma_{G}$. At early epochs, the peak of the difference is at $\delta\sim
0$, but shifts to the void as non-linear evolution progresses. This is because
of the inherent asymmetry in the dynamics: $\delta\geq-1$ in voids, but it has
no bound in the overdense limit. The temporal evolution of this peak will
potentially different for other initial conditions.

The width of the PDF, defined to be the standard deviation, was computed for
$\sigma_{G}=0.2$ and 0.5. The observed differences between the two models were
about 10%. Both measures, differences and the width can be used as potential
probes to discriminate between models, however, the differences have a slight
advantage. Within this simple model, defining the width for higher values of
$\sigma$ is subject to imposed cut-offs. More realistically, theoretical
modelling in the tails is complicated by other effects (for example baryonic
physics) and observationally, the tails are harder to measure and are ridden
by many uncertainties such as galaxy bias and sample variance. On the other
hand, the differences, particularly in the bulk, have information encoded in
the peak as well as the amplitude and this additional information can be used
to break parameter degeneracies (see for example Friedrich et al. 2019).

  * •

We computed the growth rate of perturbations in the linear and non-linear
regime; the growth as expected is slower in EDE than dark energy. One would
naively expect that the difference in the growth rates would be magnified in
the non-linear regime as compared to the linear regime. However, the linear
growth rate depends only on cosmological parameters whereas the non-linear
growth rate depends not only on the cosmological parameters, but also on the
instantaneous density contrast and hence in turn is sensitive to initial
conditions. Thus, we found that in overdensities, if initial conditions are
the same in both models then the growth rates in the non-linear regime are a
very sensitive probe of cosmology. However, if the initial conditions are
chosen such that the linear contrast is the same for both models at the
current epoch, then the effect of cosmology is damped by this constraint. In
voids, the non-linear growth rate in a void is proportional to
$\Omega_{m}^{0.56}(1+\delta)$. These two factors act in opposite giving rise
to a ‘self-regulatory’ behaviour. A higher value of the matter density
parameter implies a faster transition to the non-linear regime, however, a
higher non-linearity implies a lower value of $(1+\delta)$ bringing down the
net growth rate. Thus, the non-linear growth rate in voids is neither
significantly altered by the choice of initial conditions, nor is it too
sensitive to cosmology. Nevertheless, the small differences between models get
amplified in the PDF because the non-linear PDF is skewed with a cut-off at
$\delta\geq-1$. For initially symmetric initial conditions (as given by a
Gaussian), the differences in the non-linear PDF are higher in voids. The case
may potentially be different in the presence of a primordial non-Gaussianity.

  * •

The form of the non-linear DVDR given by N13 was tested by examining the
observed and predicted non-linear growth rates. The predicted growth rate was
computed by combining the definition given by eq. (16) with the fitting form
of $\Theta$ from N13, given by eq. (17). The error between observed and
predicted rate is, less than 4% confirming the validity of the N13 formula for
EDE models. This suggests that the formula is robust even though it was
derived based on constant $w$ dark energy models.

Observationally, the one-point PDF is measured either using galaxy catalogs or
using weak lensing convergence maps. The survey either measures positions of
galaxies or observed ellipticities to compute a shear map. The density map is
constructed using the Cloud-In-Cell (CIC) method which involves defining 2D or
3D cells in space, adding up the matter content in each cell and allotting a
density to the cell depending upon the cell volume (for e.g., Bernardeau et
al. 2002). This method introduces shot noise because of the finite cell size.
This noise is often modelled by a Poisson distribution in the case of galaxies
and a Gaussian in the case of the weak lensing convergence field (for e.g.,
Clerkin et al. 2017). In order to compare to data, the theoretically predicted
continuous distribution has to be convolved with the noise distribution. In
the case of galaxies, one needs to further account for the fact that they are
biased tracers of the underlying dark matter distribution. Both these effects
may introduce effects that are degenerate with the effect of changing
cosmology (Szapudi & Pan, 2004; Manera & Gaztañaga, 2011; Jee et al., 2012;
Uhlemann et al., 2018a).

The use of spherical dynamics to model the density one-point PDF is not new.
Many investigations have coupled SC with some non-local techniques like
perturbation theory, Edgeworth expansion, excursion sets, large deviation
theory, path integral formalism, all of which account for the interaction
between modes. However, in this analysis, we do not couple SC with any non-
local scheme, which makes this a rather simple technique. Our study is
motivated by the future goal of finding the one-point PDF for non-standard
dark energy models (for e.g., Bagla et al. 2003; Pratap Rajvanshi & Bagla
2018; Sangwan et al. 2018; Gupta et al. 2015; Devi & Sen 2011; Chakraborty et
al. 2019) or modified gravity. Complications arise in these situations because
the dynamical variables are no longer just density and velocity. In the case
of modified gravity, the ratio of the two gravitational potentials also
becomes a dynamical variable (Nadkarni-Ghosh et al, in preparation), whereas
in the case of dark energy coupled to dark matter, the dark energy
perturbation introduces an additional degree of freedom. Top-hats do not stay
as top-hats (e.g., Dai et al. 2008; Borisov et al. 2012; Kopp et al. 2013) and
defining the density PDF will involve taking a spatial average of the density
for each initial perturbation. These issues can be better understood, with a
simple scheme. Our work is perhaps most similar in spirit to that of Ohta et
al. (2003) who obtained an analytical expression for the evolution equation of
the PDF whereas, we follow the evolution of the PDF numerically. The new
features of our analysis are the one-point PDF and the expression for the non-
linear growth rate, which is derived from the density velocity divergence
relation (eq. 18).

The SC model has numerous limitations. It is based on local dynamics. Thus,
the two-point correlation functions and higher-order spatial statistics cannot
be computed solely within this framework. Effect of the environment which
induces tidal forces, rotation and accretion is ignored and so is mode-
coupling and galaxy bias (see Fosalba & Gaztanaga (1998a) for a detailed
analysis of ignoring the tidal fields in SC while modelling cumulants).
Estimates based on perturbative schemes such as Lagrangian perturbation theory
(Nadkarni-Ghosh & Chernoff, 2011, 2013) or other approximate methods Monaco
(2016) can account for some of the neglected effects. Nevertheless, the method
presented in this paper allows for a quick comparative study of the evolution
of the one-point PDF in various cosmologies.

##  7 Acknowledgements

We would like to thank the anonymous referee for very detailed and valuable
comments which improved this work significantly. SN would like to acknowledge
the Department of Science and Technology (DST), Govt. of India for the grant
no. (SR/WOS-A/PM-21/2018).

##  Appendix A The dynamical equations

The physical system comprises of a perturbed sphere evolving in a background
cosmology consisting of matter and dark energy with an equation of state $w$.
The equation of state can be varying or time-dependent. The background
evolution is completely described by a scale factor $a(t)$ given by eq. (1).

Let the origin be at the center of the inner sphere (perturbation) and ${\bf
r}$ and ${\bf x}$ denote the physical and comoving coordinate of the outer
edge of the sphere. Let ${\bf u}={\dot{\bf x}}$ and $\delta$ denote the
fractional overdensity. For pure matter cosmologies, it is possible to write
down an analogous Friedmann equation by replacing the scale factor of the
background by that of the perturbation with a different matter density
parameter. However, the extension to other dark energy models is not obvious.
Therefore, it is preferable to start with the equations of hydrodynamics,
generalized to include dark energy. This set corresponds to the spherically
symmetric, sub-horizon limit of the perturbed Einsteins’ equations along with
the conservation of the stress-energy tensor. In the weak-field, sub-horizon
limit, for pure matter perturbations, the equations are:

$\displaystyle{\dot{\delta}}$ $\displaystyle=$
$\displaystyle-(1+\delta)\nabla_{x}\cdot{\bf u}$ (36) $\displaystyle{\dot{\bf
u}}+2H{\bf u}$ $\displaystyle=$ $\displaystyle=-\frac{1}{a^{2}}\nabla_{x}\psi$
(37) $\displaystyle\nabla_{x}^{2}\psi$ $\displaystyle=$
$\displaystyle\frac{3}{2}\ a^{2}\Omega_{m}(a)H^{2}(a)\delta,$ (38)

where the ‘dot’ refers to the total derivative w.r.t. cosmic time $t$
($d/dt=\partial/\partial t+{\bf u}\cdot\nabla_{x}$) and the partial spatial
derivatives are with respect to the comoving coordinate ${\bf x}$. Taking
divergence of eq. (37), gives

$\nabla_{x}\cdot{\dot{\bf u}}+2H\nabla_{x}\cdot{\bf
u}=-\frac{3}{2}H^{2}\Omega_{m}(a)\delta.$ (39)

For the spherically symmetric system considered here, we have

$\frac{d(\nabla_{x}\cdot{\bf u})}{dt}=\nabla_{x}\cdot{\dot{\bf
u}}-\frac{(\nabla_{x}\cdot{\bf u})^{2}}{3}$ (40)

Define the scaled velocity divergence as

$\Theta=\frac{\nabla_{x}\cdot{\bf u}}{H}.$ (41)

Equations (36) and (39) give the coupled system

$\displaystyle{\dot{\delta}}$ $\displaystyle=$
$\displaystyle-H(1+\delta)\Theta$ (42) $\displaystyle{\dot{\Theta}}$
$\displaystyle=$
$\displaystyle-H\left[\frac{3}{2}\Omega_{m}(a)\delta+\left(2+\frac{\dot{H}}{H^{2}}\right)\Theta+\frac{\Theta^{2}}{3}\right].$
(43)

Combining the two first order equations gives the second order system

${\ddot{\delta}}+2H{\dot{\delta}}-\frac{4}{3}\frac{{\dot{\delta}}^{2}}{(1+\delta)}-\frac{3}{2}H^{2}\Omega_{m}(a)\delta(1+\delta)=0.$
(44)

Using $\ln a$ as the time variable instead of $t$ ($d/dt=Hd/d\ln a$), gives
the coupled system

$\displaystyle\delta^{\prime}$ $\displaystyle=$
$\displaystyle-(1+\delta)\Theta$ (45) $\displaystyle\Theta^{\prime}$
$\displaystyle=$
$\displaystyle-\frac{3}{2}\Omega_{m}(a)\delta-\Theta\left(2+\frac{H^{\prime}}{H}\right)-\frac{\Theta^{2}}{3}$
(46)

and the second order version

$\delta^{\prime\prime}+\left(2+\frac{H^{\prime}}{H}\right)\delta^{\prime}-\frac{4}{3}\frac{\delta^{\prime
2}}{(1+\delta)}-\frac{3}{2}\Omega_{m}(a)\delta(1+\delta)=0.$ (47)

For the dark energy cosmologies considered in this paper:

$\frac{1}{H}\frac{dH}{d\ln
a}=-\frac{3}{2}\Omega_{m}(a)-\frac{3}{2}(1+w(a))\Omega_{de}(a).$ (48)

For a flat universe: $\Omega_{m}(a)+\Omega_{de}(a)=1$ at all epochs. So the
factor

$\left(2+\frac{H^{\prime}}{H}\right)=\frac{1}{2}\left[1-3w(a)\Omega_{de}(a)\right].$
(49)

This gives

$\delta^{\prime\prime}-\frac{4}{3}\frac{\delta^{\prime
2}}{(1+\delta)}+\frac{1}{2}\left[1-3w(a)\Omega_{de}(a)\right]\delta^{\prime}-\frac{3}{2}\Omega_{m}(a)\delta(1+\delta)=0.$
(50)

This second order non-linear equation for the evolution of $\delta$ has been
considered in the past by Fosalba & Gaztanaga (1998b) as well as Ohta et al.
(2003), but applied to pure matter cosmologies.

In the linear regime, eq. (47) reduces to

$\delta^{\prime\prime}+\left(2+\frac{H^{\prime}}{H}\right)\delta^{\prime}-\frac{3}{2}\Omega_{m}(a)\delta=0$
(51)

The general solution to the linear system is (see equations 7.73, chapter 7
Dodelson 2003)

$\delta_{lin}(a)=C_{1}H(a)+C_{2}D(a),$ (52)

where

$D(a)\propto H(a)\int_{0}^{a}\frac{da^{\prime}}{a^{\prime
3}H(a^{\prime})^{3}}.$ (53)

Requiring that there be no decaying modes, sets $C_{1}=0$. Thus, the linear
solution is

$\delta_{lin}(a)=\delta(a_{i})\frac{D(a)}{D(a_{i})},$ (54)

where $\delta(a_{i})$ and $D(a_{i})$ are the initial values of the density and
the growth factor. Define the scaled growth factor

${\tilde{D}}=\frac{D(a)}{D(a_{i})}.$ (55)

Figure 9 shows the growth factor and various ratios for the cosmologies
discussed in this paper. Growth factor is different from the growth rate which
is defined as

$f=\frac{d\ln\delta}{d\ln a}=\frac{d\ln D(a)}{d\ln a},$ (56)

and is discussed in detail in the main text.

Figure 9: Left panel: the growth factor ${\tilde{D}}$, normalized to unity at
the initial time. For the EdS cosmology it is proportional to $a$.Until the
dark energy terms become prominent, the other models too behave like the $EdS$
model. Right Panel: Ratios of the growth factors.The EDE and $\Lambda$CDM have
a lower growth factor as compared to the EdS case; the EDE cosmology has a
slower growth as compared to $\Lambda$CDM.

##  Appendix B Comparison between the EDE and CPL models

The equation of state for the EDE model considered in this paper is given by
eq. (10). The EDE model contains a low but non-vanishing amount of dark energy
at early epoch. It can be characterized by the quantity
${\bar{\Omega}_{de,sf}}$ defined in eq. (12) (Doran et al., 2001, 2007) which
is the average fraction of dark energy during the structure formation/matter
domination era. Using a combination of data sets, including WMAP-3 year data,
Doran, Robbers, & Wetterich (2007) found ${\bar{\Omega}_{de,sf}}\leq 0.04$. A
randomly chosen model that satisfied these constraints is
$\Omega_{de,e}=8\times 10^{-4}$ and $w_{0}=-0.99$ Bartelmann, Doran, &
Wetterich (2006) . Another popular parameterization of dynamical dark energy
models is the CPL parameterization (Chevallier & Polarski, 2001; Linder, 2003)

$w(a)=w_{0}+w_{a}(1-a)$ (57)

where $w_{0}$ and $w_{a}$ are two free parameters and $w_{0}$ is to the value
of $w$ at present time. Most recent observational constraint from Planck
Collaboration et al. (2018), which combines the Planck, SNe and BAO data is
$w_{0}=-0.957$ and $w_{a}=-0.29$.

Figure 10: The left panel shows the evolution of the equation of state
parameter $w$ and as a function of scale factor $a$. The right panel shows the
evolution of density parameters $\Omega_{m}(a)$ (dashed line) and
$\Omega_{de}(a)$ (solid line) as a function of $a$. The blue and the black
lines correspond to the EDE model (with $w_{0}=-0.99$, $\Omega_{de,e}=8\times
10^{-4}$) and CPL model (with $w_{0}=-0.954$, $w_{a}=-0.29$) respectively. The
red line in both panels shows the behaviour in standard $\Lambda$CDM
cosmology.

Figure 10shows the equation of state as well as the matter and dark energy
fractional densities for the EDE and CPL models. It is clear that the two
parameterizations have very different behaviour in the early epoch limit. In
CPL models, $w$ is not constant but the dark energy density tends to zero when
$a$ tends to zero. In contrast, the EDE models have $w(a)\rightarrow 0$ and
the dark energy density is finite when $a\rightarrow 0$. The CPL
parameterization corresponds to a Taylor expansion around $a=1$. Performing a
similar expansion for the EDE equation of state gives

$w(a)=\frac{w_{0}}{[1-b\ln[1-(1-a)]]^{2}}\sim\frac{w_{0}}{[1+b(1-a)]^{2}}\sim
w_{0}-2bw_{0}(1-a)+\mathcal{O}((1-a)^{2})$ (58)

Comparing with eq. (57), we can relate $w_{a}=-2bw_{0}$, since $w_{0}$ denotes
the same quantity in both expressions. For our choice of EDE parameters, we
have $-2bw_{0}=0.73$. Thus, although, eq. (57) and eq. (58) and have the same
form, the value of the parameters are different if consider the recent
constraints from the observational data. In particular, the signs are
different and the equation $w$ is less (EDE) and more (CPL) negative as
compared to the $\Lambda$CDM value of -1.

##  Appendix C Converting to the Eulerian PDF

Given a discrete distribution of matter with mean density ${\bar{\rho}}$,
there are two ways to define the density field. Fix the volume $V_{0}$ around
each (fixed grid) point and estimate the matter enclosed inside it to define
the overdensity $\delta$. Alternately, fix the mass $M_{0}$ associated with
each point and estimate the volume it occupies. The former corresponds to the
Eulerian prescription whereas the latter to the Lagrangian prescription. The
fixed mass $M_{0}$ is related to the scale $R_{f}$ through

$M_{0}=\frac{4\pi}{3}R_{f}^{3}{\bar{\rho}},$ (59)

and the $\sigma_{G}$ is related to the scale $R_{f}$ through

$\sigma^{2}(R_{f})=\int P(k)W^{2}(kR_{f})d^{3}k,$ (60)

where $P(k)$ is the power spectrum. The Eulerian and Lagrangian PDFs are
related as

$(1+\delta)P_{E}(\delta,\sigma_{G})=\left\\{P_{L}(\delta,\sigma_{G})-\frac{d\sigma_{G}}{dR_{f}}\frac{dR_{f}}{d\delta}\int_{\delta}^{\infty}\frac{dP_{L}(\delta^{\prime},\sigma_{G})}{d\sigma_{G}}d\delta^{\prime}\right\\}$
(61)

This expression is equivalent to one derived by Bernardeau (1994a), equation
33, and has been re-derived in the context of local dynamics by Nadkarni-Ghosh
& Singhal (2016), appendix E. The PDF obtained by the method based on
spherical collapse is Lagrangian since it assumes that the mass of each sphere
in the ensemble is the same (set by the scale $\sigma_{G}$) and it remains
constant throughout the evolution. However, the fitting forms based on
Eulerian perturbation theory or numerical simulations correspond to the
Eulerian PDF. The Lagrangian PDF was converted to the Eulerian one using the
above formula. However, we do not know $P_{L}(\delta^{\prime},\sigma_{G})$
analytically. Thus, to find the derivative of $P_{L}$ in the integrand, we
performed numerical runs at $\sigma_{G}=0.6,1.1$ and $2.1$. The data for all
$\sigma_{G}$ was binned in the same range and discrete derivatives were
computed in each $\delta$ bin. The definition
$\delta=M_{0}/{\bar{\rho}}V_{0}-1$ gives the pre-factor
$\frac{dR_{f}}{d\delta}=R_{f}/(1+\delta)$ (for fixed $V_{0}$). The initial
realization is simply an ensemble of points drawn from a Gaussian and does not
correspond to any power spectrum (the spherical dynamics assumes only local
collapse; there is no reference to spatial correlations). Ideally, to compute
the pre-factor $d\sigma_{G}/dR_{F}$, one needs to know the initial power
spectrum that is chosen by a particular simulation, whose results one wants to
compare. However, in this paper, we do not aim to compare our results to a
particular output of a simulation. It is important to keep in mind that the
forms given below only provide the analytic dependence on $\delta$. The
parameters of the fitting functions (for e.g., $\sigma_{\phi}$ in SLN or
$\alpha$, $\kappa$ and $\xi$ of $N_{v2}$) are computed from the data generated
by spherical dynamics and not the N-body simulation and we only aim to check
which form fits best. In view of this, we approximate the pre-factor
$d\sigma_{G}/dR_{F}$ using a BBKS power spectrum under the assumption that the
order of magnitude of the pre-factor will not be significantly different with
any other more realistic power spectrum. The transformed PDF does not satisfy
$\int P_{E}(\delta)d\delta=1$; it is explicitly normalized by dividing by the
norm. A more accurate transformation for general PDFs based on 1d systems is
given in Pajer & van der Woude (2018).

##  Appendix D Best fit parameters for the $\mathcal{N}_{v2}$ distribution for
the three cosmologies.

Model | $\sigma_{G}$ | $a$ | $\alpha$ | $\kappa$ | $\xi$
---|---|---|---|---|---
Lagrangian | Eulerian | Lagrangian | Eulerian | Lagrangian | Eulerian
$\Lambda$CDM | 0.5 | 0.1 | 0.064 | 0.064 | -0.108 | -0.113 | 0.999 | 0.997
0.4 | 0.254 | 0.234 | -0.389 | -0.379 | 1.000 | 0.958
1.0 | 0.504 | 0.351 | -0.670 | -0.599 | 0.988 | 0.856
| 1.0 | 0.1 | 0.128 | 0.125 | -0.201 | -0.200 | 1.001 | 0.987
| 0.4 | 0.476 | 0.342 | -0.645 | -0.586 | 0.993 | 0.824
| 1.0 | 0.832 | 0.323 | -0.948 | -0.749 | 0.990 | 0.567
| 2.0 | 0.1 | 0.256 | 0.230 | -0.389 | -0.365 | 0.9998 | 0.945
| 0.4 | 0.766 | 0.305 | -0.987 | -0.716 | 0.904 | 0.543
| 1.0 | 0.848 | 0.272 | -1.223 | -0.936 | 0.743 | 0.346
EDE | 0.5 | 0.1 | 0.067 | 0.067 | -0.113 | -0.119 | 0.999 | 0.997
0.4 | 0.245 | 0.228 | -0.377 | -0.372 | 1.000 | 0.962
1.0 | 0.539 | 0.345 | -0.671 | -0.583 | 1.078 | 0.870
1.0 | 0.1 | 0.135 | 0.132 | -0.213 | -0.210 | 1.001 | 0.986
0.4 | 0.495 | 0.342 | -0.665 | -0.579 | 1.008 | 0.838
1.0 | 0.728 | 0.331 | -0.867 | -0.742 | 0.977 | 0.592
2.0 | 0.1 | 0.269 | 0.240 | -0.406 | -0.384 | 0.999 | 0.939
0.4 | 0.612 | 0.307 | -0.840 | -0.707 | 0.863 | 0.549
1.0 | 0.851 | 0.233 | -1.134 | -0.874 | 0.812 | 0.327
Table 1: Values of the parameters $\alpha$, $\kappa$ and $\xi$ for best fitted
$\mathcal{N}_{v2}$ distribution:

##  Appendix E Width of the one-point PDF for higher values of $\sigma_{G}$.

Figure 11: Width of the PDF (top panel) and scaled width (bottom panel) as a
function of time. The flattening for $\sigma_{G}=2$ is an effect of a finite
cutoff imposed on the $\delta$ values in computing the PDF.

The upper panel of figure 11 shows the variation of $\sigma_{NL}$ as a
function of $a$ ($a=0.1,0.4,1$). It is clear that for all values of
$\sigma_{G}$ the width of the PDF increases with time. Between $a=0.1$ and
$a=0.4$, the rate of increase is higher for higher $\sigma_{G}$, which is
expected because the non-linear growth rate is higher for higher $\delta$
values. From $a=0.4$ to $a=1$, there is an apparent slow-down in the growth
rate for $\sigma_{G}=1,2$. While slow down is expected at higher epochs due to
the presence of dark energy, the one observed slow down is an artefact of the
cut-off in $\delta$ imposed in the construction of the PDF. For
$\sigma_{G}=2$, that the change (increase) in $\sigma_{NL}$ from $a=0.4$ to
$a=1$ was higher when the cut-off $\delta$ was increased from 5 (plot not
shown) to 20 (this plot).

## References

  * Abramo et al. (2009) Abramo L. R., Batista R. C., Liberato L., Rosenfeld R., 2009, Physical Review D, 79, 023516
  * Bagla et al. (2003) Bagla J. S., Jassal H. K., Padmanabhan T., 2003, Physical Review D, 67, 063504, arXiv: astro-ph/0212198
  * Bartelmann et al. (2006) Bartelmann M., Doran M., Wetterich C., 2006, Astronomy and Astrophysics, 454, 27
  * Barthelemy et al. (2020) Barthelemy A., Codis S., Uhlemann C., Bernardeau F., Gavazzi R., 2020, Monthly Notices of the Royal Astronomical Society, 492, 3420
  * Bel et al. (2016) Bel J. et al., 2016, Astronomy and Astrophysics, 588, A51
  * Bernardeau (1992) Bernardeau F., 1992, The Astrophysical Journal, 390, L61
  * Bernardeau (1994a) Bernardeau F., 1994a, Astronomy and Astrophysics, 291, 697
  * Bernardeau (1994b) Bernardeau F., 1994b, The Astrophysical Journal, 433, 1
  * Bernardeau et al. (1999) Bernardeau F., Chodorowski M. J., Lokas E. L., Stompor R., Kudlicki A., 1999, Monthly Notices of the Royal Astronomical Society, 309, 543
  * Bernardeau et al. (2015) Bernardeau F., Codis S., Pichon C., 2015, Monthly Notices of the Royal Astronomical Society, 449, L105
  * Bernardeau et al. (2002) Bernardeau F., Colombi S., Gaztaï¿œaga E., Scoccimarro R., 2002, Physics Reports, 367, 1
  * Bernardeau & Kofman (1995) Bernardeau F., Kofman L., 1995, The Astrophysical Journal, 443, 479
  * Bernardeau et al. (2014) Bernardeau F., Pichon C., Codis S., 2014, Physical Review D, 90, 103519
  * Bernardeau & Reimberg (2016) Bernardeau F., Reimberg P., 2016, Physical Review D, 94, 063520
  * Bilicki & Chodorowski (2008) Bilicki M., Chodorowski M. J., 2008, Monthly Notices of the Royal Astronomical Society, 391, 1796
  * Borisov et al. (2012) Borisov A., Jain B., Zhang P., 2012, Physical Review D, 85, 63518
  * Brenier et al. (2003) Brenier Y., Frisch U., HÃ©non M., Loeper G., Matarrese S., Mohayaee R., SobolevskiÄ A., 2003, Monthly Notices of the Royal Astronomical Society, 346, 501
  * Buchert (1992) Buchert T., 1992, Monthly Notices of the Royal Astronomical Society, 254, 729
  * Chakraborty et al. (2019) Chakraborty A., Banerjee N., Ghosh A., 2019, General Relativity and Gravitation, 51, 5
  * Chevallier & Polarski (2001) Chevallier M., Polarski D., 2001, International Journal of Modern Physics D, 10, 213
  * Clerkin et al. (2017) Clerkin L. et al., 2017, Monthly Notices of the Royal Astronomical Society, 466, 1444
  * Codis et al. (2016) Codis S., Pichon C., Bernardeau F., Uhlemann C., Prunet S., 2016, Monthly Notices of the Royal Astronomical Society, 460, 1549
  * Coles & Jones (1991) Coles P., Jones B., 1991, Monthly Notices of the Royal Astronomical Society, 248, 1
  * Coles et al. (1993) Coles P., Melott A. L., Shandarin S. F., 1993, Monthly Notices of the Royal Astronomical Society, 260, 765
  * Colombi (1994) Colombi S., 1994, The Astrophysical Journal, 435, 536
  * Colombi et al. (1997) Colombi S., Bernardeau F., Bouchet F. R., Hernquist L., 1997, Monthly Notices of the Royal Astronomical Society, 287, 241
  * Colombi et al. (2000) Colombi S., Szapudi I., Jenkins A., Colberg J., 2000, Monthly Notices of the Royal Astronomical Society, 313, 711
  * Dai et al. (2008) Dai D.-C., Maor I., Starkman G., 2008, Physical Review D, 77, 64016
  * Devi & Sen (2011) Devi N. C., Sen A. A., 2011, Monthly Notices of the Royal Astronomical Society, 413, 2371
  * Dodelson (2003) Dodelson S., 2003, Modern Cosmology, 2nd edn. Academic Press
  * Doran & Robbers (2006) Doran M., Robbers G., 2006, Journal of Cosmology and Astro-Particle Physics, 2006, 026
  * Doran et al. (2007) Doran M., Robbers G., Wetterich C., 2007, Physical Review D, 75, 023003
  * Doran et al. (2001) Doran M., Schwindt J., Wetterich C., 2001, Phys. Rev. D, 64, 123520
  * Fontanot et al. (2012) Fontanot F., Springel V., Angulo R. E., Henriques B., 2012, Monthly Notices of the Royal Astronomical Society, 426, 2335
  * Fosalba & Gaztanaga (1998a) Fosalba P., Gaztanaga E., 1998a, Monthly Notices of the Royal Astronomical Society, 301, 503
  * Fosalba & Gaztanaga (1998b) Fosalba P., Gaztanaga E., 1998b, Monthly Notices of the Royal Astronomical Society, 301, 535
  * Francis et al. (2008) Francis M. J., Lewis G. F., Linder E. V., 2008, 0808.2840, mon.Not.Roy.Astron.Soc.394:605-614,2008
  * Friedrich et al. (2018) Friedrich O., Gruen D., DeRose J., Collaboration D., 2018, Phys. Rev. D, 98, 023508
  * Friedrich et al. (2019) Friedrich O., Uhlemann C., Villaescusa-Navarro F., Baldauf T., Manera M., Nishimichi T., 2019, arXiv e-prints, 1912, arXiv:1912.06621
  * Gaztañaga et al. (2000) Gaztañaga E., Fosalba P., Elizalde E., 2000, The Astrophysical Journal, 539, 522
  * Grossi & Springel (2008) Grossi M., Springel V., 2008, 0809.3404
  * Gruen et al. (2018) Gruen D. et al., 2018, Physical Review D, 98, 023507
  * Gupta et al. (2015) Gupta G., Rangarajan R., Sen A. A., 2015, Physical Review D, 92, 123003
  * Hilbert et al. (2011) Hilbert S., Hartlap J., Schneider P., 2011, Astronomy and Astrophysics, 536, A85
  * Hubble (1934) Hubble E., 1934, The Astrophysical Journal, 79, 8
  * Hurtado-Gil et al. (2017) Hurtado-Gil L., MartÃnez V. J., Arnalte-Mur P., Pons-BorderÃa M.-J., Pareja-Flores C., Paredes S., 2017, Astronomy and Astrophysics, 601, A40
  * Ivanov et al. (2019) Ivanov M. M., Kaurov A. A., Sibiryakov S., 2019, Journal of Cosmology and Astro-Particle Physics, 2019, 009
  * Jee et al. (2012) Jee I., Park C., Kim J., Choi Y.-Y., Kim S. S., 2012, The Astrophysical Journal, 753, 11
  * Jennings et al. (2010) Jennings E., Baugh C. M., Angulo R. E., Pascoli S., 2010, Monthly Notices of the Royal Astronomical Society, 401, 2181
  * Juszkiewicz et al. (1995) Juszkiewicz R., Weinberg D. H., Amsterdamski P., Chodorowski M., Bouchet F., 1995, Astrophysical Journal, 442, 39
  * Kitaura & Hess (2013) Kitaura F. S., Hess S., 2013, Monthly Notices of the Royal Astronomical Society, 435, L78
  * Klypin et al. (2018) Klypin A., Prada F., Betancort-Rijo J., Albareti F. D., 2018, Monthly Notices of the Royal Astronomical Society, 481, 4588
  * Kofman (1993) Kofman L., 1993, arXiv e-prints, astro
  * Kofman et al. (1994) Kofman L., Bertschinger E., Gelb J. M., Nusser A., Dekel A., 1994, The Astrophysical Journal, 420, 44
  * Kopp et al. (2013) Kopp M., Appleby S. A., Achitouv I., Weller J., 2013, Physical Review D, 88, 084015
  * Kudlicki et al. (2000) Kudlicki A., Chodorowski M., Plewa T., Różyczka M., 2000, Monthly Notices of the Royal Astronomical Society, 316, 464
  * Lahav et al. (1993) Lahav O., Itoh M., Inagaki S., Suto Y., 1993, The Astrophysical Journal, 402, 387
  * Lam & Sheth (2008a) Lam T. Y., Sheth R. K., 2008a, Monthly Notices of the Royal Astronomical Society, 389, 1249
  * Lam & Sheth (2008b) Lam T. Y., Sheth R. K., 2008b, Monthly Notices of the Royal Astronomical Society, 386, 407
  * Lima et al. (1997) Lima J. A. S., Zanchin V., Brandenberger R., 1997, Monthly Notices of the Royal Astronomical Society, 291, L1
  * Linder (2003) Linder E. V., 2003, Physical Review Letters, 90, 091301
  * Linder (2005) Linder E. V., 2005, Physical Review D, 72, 43529
  * Liu et al. (2016) Liu J., Hill J. C., Sherwin B. D., Petri A., BÃ¶hm V., Haiman Z., 2016, Physical Review D, 94, 103501
  * Manera & Gaztañaga (2011) Manera M., Gaztañaga E., 2011, MNRAS, 415, 383
  * Monaco (1999) Monaco P., 1999, Observational Cosmology: The Development of Galaxy Systems, 176, 186
  * Monaco (2016) Monaco P., 2016, Galaxies, 4, 53
  * Nadkarni-Ghosh (2013) Nadkarni-Ghosh S., 2013, Monthly Notices of the Royal Astronomical Society, 428, 1166
  * Nadkarni-Ghosh & Chernoff (2011) Nadkarni-Ghosh S., Chernoff D. F., 2011, Monthly Notices of the Royal Astronomical Society, 410, 1454
  * Nadkarni-Ghosh & Chernoff (2013) Nadkarni-Ghosh S., Chernoff D. F., 2013, Monthly Notices of the Royal Astronomical Society, 431, 799
  * Nadkarni-Ghosh & Singhal (2016) Nadkarni-Ghosh S., Singhal A., 2016, Monthly Notices of the Royal Astronomical Society, 457, 2773
  * Neyman (1962) Neyman J., 1962, Problems of Extra-Galactic Research, 15, 294
  * Ohta et al. (2003) Ohta Y., Kayo I., Taruya A., 2003, The Astrophysical Journal, 589, 1
  * Ohta et al. (2004) Ohta Y., Kayo I., Taruya A., 2004, The Astrophysical Journal, 608, 647
  * Pace et al. (2017) Pace F., Meyer S., Bartelmann M., 2017, Journal of Cosmology and Astro-Particle Physics, 10, 040
  * Pace et al. (2010) Pace F., Waizmann J.-C., Bartelmann M., 2010, Monthly Notices of the Royal Astronomical Society, 406, 1865
  * Pajer & van der Woude (2018) Pajer E., van der Woude D., 2018, Journal of Cosmology and Astroparticle Physics, 05, 039
  * Patton et al. (2017) Patton K., Blazek J., Honscheid K., Huff E., Melchior P., Ross A. J., Suchyta E., 2017, Monthly Notices of the Royal Astronomical Society, 472, 439
  * Peebles (1980) Peebles P. J. E., 1980, The large-scale structure of the universe. Princeton University Press
  * Planck Collaboration et al. (2018) Planck Collaboration et al., 2018, arXiv e-prints, arXiv:1807.06209
  * Pratap Rajvanshi & Bagla (2018) Pratap Rajvanshi M., Bagla J. S., 2018, Journal of Cosmology and Astroparticle Physics, 06, 018
  * Reimberg & Bernardeau (2018) Reimberg P., Bernardeau F., 2018, Physical Review D, 97, 023524
  * Repp & Szapudi (2018) Repp A., Szapudi I., 2018, Monthly Notices of the Royal Astronomical Society, 473, 3598
  * Russell & Pycke (2017) Russell E., Pycke J.-R., 2017, The Astrophysical Journal, 835, 69
  * Sangwan et al. (2018) Sangwan A., Mukherjee A., Jassal H. K., 2018, Journal of Cosmology and Astroparticle Physics, 01, 018
  * Saslaw & Hamilton (1984) Saslaw W. C., Hamilton A. J. S., 1984, The Astrophysical Journal, 276, 13
  * Scott (1962) Scott E. L., 1962, Problems of Extra-Galactic Research, 15, 269
  * Shi & Baugh (2016) Shi D., Baugh C. M., 2016, Monthly Notices of the Royal Astronomical Society, 459, 3540
  * Shin et al. (2017) Shin J., Kim J., Pichon C., Jeong D., Park C., 2017, ArXiv e-prints, 1705, arXiv:1705.06863
  * Suto et al. (1990) Suto Y., Itoh M., Inagaki S., 1990, The Astrophysical Journal, 350, 492
  * Szapudi et al. (2000) Szapudi I., Colombi S., Jenkins A., Colberg J., 2000, Monthly Notices of the Royal Astronomical Society, 313, 725
  * Szapudi & Pan (2004) Szapudi I., Pan J., 2004, The Astrophysical Journal, 602, 26
  * Taruya et al. (2002) Taruya A., Takada M., Hamana T., Kayo I., Futamase T., 2002, The Astrophysical Journal, 571, 638
  * Taylor & Watts (2000) Taylor A. N., Watts P. I. R., 2000, Monthly Notices of the Royal Astronomical Society, 314, 92
  * Tosone et al. (2020) Tosone F., Neyrinck M. C., Granett B. R., Guzzo L., Vittorio N., 2020, arXiv e-prints, 2006, arXiv:2006.04653
  * Ueda & Yokoyama (1996) Ueda H., Yokoyama J., 1996, Monthly Notices of the Royal Astronomical Society, 280, 754
  * Uhlemann et al. (2017) Uhlemann C., Codis S., Hahn O., Pichon C., Bernardeau F., 2017, Monthly Notices of the Royal Astronomical Society, 469, 2481
  * Uhlemann et al. (2016) Uhlemann C., Codis S., Pichon C., Bernardeau F., Reimberg P., 2016, Monthly Notices of the Royal Astronomical Society, 460, 1529
  * Uhlemann et al. (2018a) Uhlemann C. et al., 2018a, Monthly Notices of the Royal Astronomical Society, 473, 5098
  * Uhlemann et al. (2020) Uhlemann C., Friedrich O., Villaescusa-Navarro F., Banerjee A., Codis S., 2020, Monthly Notices of the Royal Astronomical Society
  * Uhlemann et al. (2018b) Uhlemann C., Pichon C., Codis S., L’Huillier B., Kim J., Bernardeau F., Park C., Prunet S., 2018b, Monthly Notices of the Royal Astronomical Society, 477, 2772
  * Valageas (1998) Valageas P., 1998, Astronomy and Astrophysics, 337, 655
  * Valageas (2001) Valageas P., 2001, Astronomy and Astrophysics, 379, 8
  * Valageas (2002a) Valageas P., 2002a, Astronomy and Astrophysics, 382, 412
  * Valageas (2002b) Valageas P., 2002b, Astronomy and Astrophysics, 382, 431
  * Valageas (2002c) Valageas P., 2002c, Astronomy and Astrophysics, 382, 450
  * Valageas (2002d) Valageas P., 2002d, Astronomy and Astrophysics, 382, 477
  * Valageas & Munshi (2004) Valageas P., Munshi D., 2004, Monthly Notices of the Royal Astronomical Society, 354, 1146
  * Wetterich (2004) Wetterich C., 2004, Physics Letters B, 594, 17
  * Wild et al. (2005) Wild V. et al., 2005, Monthly Notices of the Royal Astronomical Society, 356, 247
  * Xavier et al. (2016) Xavier H. S., Abdalla F. B., Joachimi B., 2016, Monthly Notices of the Royal Astronomical Society, 459, 3693
  * Zeldovich (1970) Zeldovich Y. B., 1970, Astronomy and Astrophysics, 5, 84