#  Cosmological perturbation theory using generalized Einstein de Sitter
cosmologies

Michael Joyce<EMAIL_ADDRESS>Laboratoire de Physique Nucléaire et de
Hautes Energies, UPMC IN2P3 CNRS UMR 7585, Sorbonne Université, 4, place
Jussieu, 75252 Paris Cedex 05, France  Azrul Pohan<EMAIL_ADDRESS>Laboratoire de Physique Nucléaire et de Hautes Energies, UPMC IN2P3 CNRS UMR
7585, Sorbonne Université, 4, place Jussieu, 75252 Paris Cedex 05, France
Laboratory of Theoretical Physics, Institut Teknologi Sumatera, Lampung 35365,
Indonesia

###### Abstract

The separable analytical solution in standard perturbation theory for an
Einstein de Sitter (EdS) universe can be generalized to the wider class of
such cosmologies (“generalized EdS”, or gEdS) in which a fraction of the
pressure-less fluid does not cluster. We derive the corresponding kernels in
both Eulerian perturbation theory (EPT) and Lagrangian perturbation theory,
generalizing the canonical EdS expressions to a one-parameter family where the
parameter can be taken to be the exponent $\alpha$ of the growing mode linear
amplification $D(a)\propto a^{\alpha}$. For the power spectrum (PS) at one
loop in EPT, the contribution additional to standard EdS is given, for each of
the ‘13’ and ‘22’ terms, as a function of two infra-red safe integrals. In the
second part of the paper we show that the calculation of cosmology-dependent
corrections in perturbation theory in standard (e.g. LCDM-like) models can be
simplified, and their magnitude and parameter dependence better understood, by
relating them to our analytic results for gEdS models. At second order the
time dependent kernels are equivalent to the analytic kernels of the gEdS
model with $\alpha$ replaced by a single redshift dependent effective growth
rate $\alpha_{2}(z)$. At third order the time evolution can be conveniently
parametrized in terms of two additional such effective growth rates. For the
PS calculated at one loop order, the correction to the PS relative to the EdS
limit can be expressed in terms of just $\alpha_{2}(z)$, one additional
effective growth rate function and the four infra-red safe integrals of the
gEdS limit. This is much simplified compared to expressions in the literature
that use six or eight red-shift dependent functions and are not explicitly
infra-red safe. Using the analytic gEdS expression for the PS with
$\alpha=\alpha_{2}(z)$ gives a good approximation (to $\sim 25\%$) for the
exact result.

††preprint: APS/123-QED

##  I Introduction

Understanding the origin of the large-scale structures in the Universe first
revealed several decades ago by early red-shift surveys remains a central and
open problem in cosmology. This is driven by the ever-increasing richness and
accuracy of observational data which will grow in the coming years with
programs such as Euclid, DESI and LSST [1, 2, 3]. In standard cosmological
models the clustering observed at cosmological scales today is the result of
evolution under the gravity of primordial fluctuations, that is highly
constrained by observations of the fluctuations in the cosmic microwave
background, with potential fluctuations $\sim 10^{-5}$ at the time of
decoupling. Calculation of predictions derived from the non-linear equations
governing the evolution of matter fluctuations rests in practice essentially
on numerical solutions using the $N$-body method (for a recent review, see
e.g. [4]). At sufficiently early times and/or large scales, when fluctuations
to uniformity are small, analytical perturbative approaches may be used and
provide an invaluable guideline and control on numerical simulations.
Perturbation theory beyond the leading (linear) order was originally treated
in early works (e.g. [5, 6, 7, 8, 9, 10, 11, 12, 13, 14]) and has been
developed since in an extensive literature (for a review, see [14]). In recent
years perturbation theory has remained highly topical, with development of
various different aspects and cosmological applications (see e.g. [15, 16, 17,
18, 19, 20, 21, 22]), and, more recently, of the so-called Effective Field
Theory (EFT) of cosmological structure formation (see e.g. [15, 23, 24, 25,
26]).

In this article we explore and derive results for the perturbation theory of
cosmological structure formation in a particular set of idealized cosmological
models. We refer to these as “generalized Einstein de Sitter” (gEdS) models as
they are models in which the background cosmological evolution is that of an
EdS cosmology. However, only a fraction of the total energy density driving
the expansion of the background corresponds to clustering matter, the rest
remaining uniform. This is simply equivalent to a smooth matter-like component
(e.g. arising from a homogeneous mode of a scalar field, or from a massive
neutrino-like contribution when perturbations are neglected). We derive
analytical solutions for the key quantities (kernels) in both standard
(Eulerian) perturbation theory (SPT) and also in Lagrangian perturbation
theory (LPT), and use them to derive the expression for the power spectrum at
the leading non-trivial order (one loop). The motivation for studying these
models is twofold. On the one hand, as we will show in the second part of this
paper, these models turn out to be very useful to better understand, and
simplify the calculation of, the small cosmology-dependent corrections to
results obtained using the usual “separable EdS” approximation made in
applying perturbation theory to the standard (i.e. LCDM-like) cosmologies [27,
28, 29, 30, 20, 31, 32, 14]. In particular we will show how expressions for
the non-EdS corrections to the one-loop PS in perturbation theory can be given
(exactly) in terms of just two redshift dependent “effective growth rates”,
and, to a good approximation, in terms of a single such function. On the other
hand, our motivation comes from the intrinsic interest of new analytical
results in perturbation theory in order to test its validity, and in
particular to probe its accuracy in predicting the dependence on expansion
history beyond the usually used approximation. The models we focus on are
particularly suitable for robust numerical tests as in the case of initial
Gaussian fluctuations specified by a simple power-law power spectrum these
gEdS models are scale-free and characterized (at least for some range of
spectra) by “self-similarity” in their evolution [33, 34]. In principle this
property allows one to obtain very accurate numerical results with which to
test perturbation theory (in particular by applying the methods recently
discussed for scale-free EdS models in [35, 36]). An analysis of these
“generalized scale-free” models and numerical tests will be reported in a
separate forthcoming article.

The article is organized as follows. The next section presents first the
family of gEdS models, and then the calculation of perturbation theory kernels
in them for both Eulerian and Lagrangian formulations, giving the density and
velocity kernels up to the third order. In section 3 we study the convergence
properties of the one-loop PS in gEdS cosmologies, considering how the well-
known results for standard EdS are modified. Section 4 considers how the
functions characterising the exact evolution of the PS in perturbation theory
in a class of standard (LCDM-like) cosmological models can be approximated by
interpolation of gEdS models. In the final section we summarize and comment on
further possible exploitation of the results we have derived, notably in
relation to understanding better the ultra-violet divergences in standard
perturbation theory. Some details of our calculations are relegated to
appendices, as well as additional detailed comparison with the previous
results of [27], and also to those of [11, 20, 32]. For readers wishing to
make use of our simplified expression for the cosmological corrections to the
EdS approximation for the one-loop PS in standard cosmologies, given in terms
of just two redshift-dependent functions, without going through the full
derivation given in the body of the article, we give the necessary equations
also in a short appendix.

##  II Perturbation Theory kernels in generalized EdS models

###  II.1 Generalized EdS models

We consider perturbed Friedman-Lemaitre-Robertson-Walker (FLRW) models with a
pressureless perturbed matter component and a smooth component (i.e. without
perturbations). The Friedmann equations can be written as

$\begin{split}\mathcal{H}^{2}&=\frac{8\pi
G}{3}a^{2}\big{[}\rho_{m}+\rho_{s}\big{]},\\\
\frac{\partial{\mathcal{H}}}{\partial{\tau}}&=-\frac{4\pi
G}{3}a^{2}\big{[}\rho_{m}+(1+3w_{s})\rho_{s}\big{]},\end{split}$ (1)

where $\mathcal{H}=d(\log a)/d\tau$ and $\tau$ is the conformal time related
to cosmic time $t$ by $\tau=\int dt/a$, $a$ is the scale factor, $\rho_{m}$ is
the mean density of the perturbed matter component, $\rho_{s}$ is the smooth
“dark energy” component with equation of state $p_{s}=w_{s}\rho_{s}$ (where
$w_{s}$ can be a function of $a$). The case $w_{s}=-1$ corresponds to the
Lambda Cold Dark Matter (LCDM) cosmology. As canonically, we define the
fraction of matter and dark energy as

$\Omega_{m}=\frac{\rho_{m}}{\rho_{tot}},\quad\Omega_{s}=\frac{\rho_{s}}{\rho_{tot}},$
(2)

where $\rho_{tot}=\rho_{m}+\rho_{s}$ is the total energy density (and
$\Omega_{m}+\Omega_{s}=1$). Once $w_{s}$ is specified, the model is fully
characterized by specifying the value at some time of one of these two
parameters. Our analysis of standard cosmologies below applies to this class
of models, with the sole further assumption that the total energy density
asymptotically scales as matter at early times.

The family of cosmologies we will focus on first, those we will refer to as
generalized Einstein de Sitter cosmologies (hereafter gEdS), simply correspond
to the case $w_{s}=0$. In this case [34, 33] the smooth component behaves
exactly like a matter component scaling as $\rho_{s}\propto 1/a^{3}$, so that
both $\Omega_{m}$ and $\Omega_{s}$ are constant. They may thus be
parameterized by this constant value of $\Omega_{m}$ (or $\Omega_{s}$). In
order to avoid confusion below with the general LCDM-type model where
$\Omega_{m}$ is a function of time, we will use, as in [34], the parameter

$\kappa^{2}=\frac{1}{\Omega_{m}}\,.$ (3)

We then have that

$\begin{split}\mathcal{H}^{2}&=\frac{8\pi G}{3}a^{2}\kappa^{2}\rho_{m},\\\
\frac{\partial{\mathcal{H}}}{\partial{\tau}}&=-\frac{4\pi
G}{3}a^{2}\kappa^{2}\rho_{m},\end{split}$ (4)

and the scale factor $a(t)$ evolves as that of an EdS model, with

$\begin{split}a=\Big{(}\frac{\kappa
t}{t_{0}}\Big{)}^{2/3}\quad\text{where}\quad t_{0}=\frac{1}{\sqrt{6\pi
G\rho_{m,0}}}\,,\end{split}$ (5)

where $\rho_{m,0}$ is the matter density at some reference time. At the level
of background cosmology this family of cosmologies are all equivalent to a
standard EdS cosmology driven by a total “matter” density
$\rho_{tot}=\kappa^{2}\rho_{m}$. However, as soon as perturbations to
uniformity are considered, the dependence on the evolution of the scale factor
on $\kappa^{2}$ when expressed in terms of the mean (clustering) mass density
translates into a physical difference of the models as a function of
$\kappa^{2}$. This is seen already, as we will show below, in the evolution of
perturbations at linear order. As in the standard EdS case ($\kappa=1$) there
is a growing mode and a decaying mode, but their temporal evolution is
modified when $\kappa\neq 1$. In particular the growing mode for density
perturbations has the behaviour $D(a)\propto a^{\alpha}$ where

$\alpha=-\frac{1}{4}+\frac{1}{4}\sqrt{1+\frac{24}{\kappa^{2}}}\,.$ (6)

As $\alpha$ is a one-to-one function of $\kappa^{2}$, either can be used to
characterize the one-parameter family of gEdS models. Because of the crucial
role played in perturbation theory by the linear theory growing mode, we will
often find it convenient to give our results as functions of $\alpha$ in our
analysis below.

We note that while it is required that $\rho_{tot}>0$, and $\rho_{m}>0$, it is
not necessary that $\rho_{s}>0$. Thus we consider the family of gEdS
cosmologies to be parameterized either by $\kappa\in[0,\infty]$ or,
alternatively, by $\alpha\in[0,\infty]$. Our motivation for studying these
models here comes, as has been outlined in the introduction, from their
interest as a theoretical tool rather than as relevant physical models. We
note however that, for $\rho_{s}>0$ (i.e. $0<\Omega_{m}<1$), a smooth
component with zero pressure can model approximately in certain regimes a
component in massive neutrinos, or also a contribution from the zero mode of a
scalar field oscillating about the minimum of a quadratic potential or rolling
in a simple exponential potential (see e.g. [37] and references therein). An
alternative physical interpretation, which extends also to models with
$\rho_{s}<0$ (or $\kappa^{2}<1$), is in terms of a model in which the
gravitational coupling is scale-dependent, with the expansion being driven by
a rescaled coupling $\kappa^{2}G$. The limit $\kappa\rightarrow 0$ corresponds
formally to the case of a static universe, in which a smooth component with
negative energy exactly balances the clustering matter (for further detail,
see [34]).

###  II.2 Eulerian Perturbation Theory (EPT)

Our starting point is the usual fluid equations for the evolution, under
gravity only, of perturbations in the matter in an FLRW universe in the
Newtonian (sub-horizon, non-relativistic) limit (see e.g. [14]):

$\displaystyle\frac{\partial\delta(\mathbf{x},\tau)}{\partial\tau}+\nabla\cdot\big{\\{}[1+\delta(\mathbf{x},\tau)]\mathbf{u}(\mathbf{x},\tau)\big{\\}}$
$\displaystyle=$ $\displaystyle 0,$ (7)
$\displaystyle\frac{\partial\mathbf{u}(\mathbf{x},\tau)}{\partial\tau}+\mathcal{H}(\tau)\mathbf{u}(\mathbf{x},\tau)+[\mathbf{u}(\mathbf{x},\tau)\cdot\nabla]\mathbf{u}(\mathbf{x},\tau)$
$\displaystyle=$
$\displaystyle-\nabla\phi(\mathbf{x},\tau)-\frac{1}{\rho}\nabla(\rho\sigma_{ij}),$
(8) $\displaystyle\nabla^{2}\phi(\textbf{x},\tau)$ $\displaystyle=$
$\displaystyle\frac{3}{2}\Omega_{m}\mathcal{H}^{2}\delta,$ (9)

where $\delta(\mathbf{x},\tau)\equiv\rho(x,\tau)/\bar{\rho}(\tau)-1$ is the
matter density fluctuation, with $\rho(x,\tau)$ the matter density and
$\bar{\rho}(\tau)$ its mean value, $\mathbf{u}$ is the (peculiar) velocity,
$\phi(\mathbf{x},\tau)$ is the (Newtonian) gravitational potential, and
$\sigma_{ij}$ is the velocity dispersion tensor. We will henceforth neglect
this last term, treating the fluid in the single flow approximation.

Using the Fourier transform conventions

$\begin{split}f(\textbf{k})&=\int
d^{3}r\exp\big{[}-i\textbf{k}\cdot\textbf{r}\big{]}f(\textbf{r}),\\\
f(\textbf{r})&=\int\frac{d^{3}k}{(2\pi)^{3}}\exp\big{[}i\textbf{k}\cdot\textbf{r}\big{]}f(\textbf{k}),\end{split}$
(10)

Eq. (7) and the divergence of Eq. (8) can be combined to give

$\delta^{\prime}(\textit{{k}},\tau)+\theta(\textit{{k}},\tau)=S_{\tilde{\alpha}}(\textit{{k}},\tau),$
(11)
$\theta^{\prime}(\textit{{k}},\tau)+\mathcal{H}\theta(\textit{{k}},\tau)+\frac{3}{2}\Omega_{m}(a)\mathcal{H}^{2}\delta(\textit{{k}},\tau)=S_{\tilde{\beta}}(\textit{{k}},\tau),$
(12)

where $\theta$ is the divergence of the velocity field $\mathbf{u}$, so that

$\textbf{u}(\textbf{k})=-i\frac{\textbf{k}}{k^{2}}\theta(\textbf{k}).$ (13)

The velocity field can be assumed to have zero vorticity because the latter
can be shown to always decay compared to the irrotational part (see e.g.
[14]). The source terms $S_{\alpha}$ and $S_{\beta}$ in Eqs. (11) and (12) are
given by

$S_{\tilde{\alpha}}(\textit{{k}},\tau)=-\int\frac{\text{d}^{3}q}{(2\pi)^{3}}\tilde{\alpha}(\textit{{q}},\textit{{k}}-\textit{{q}})\theta(\textit{{q}},\tau)\delta(\textit{{k}}-\textit{{q}},\tau),$
(14)
$S_{\tilde{\beta}}(\textit{{k}},\tau)=-\int\frac{\text{d}^{3}q}{(2\pi)^{3}}\tilde{\beta}(\textit{{q}},\textit{{k}}-\textit{{q}})\theta(\textit{{q}},\tau)\theta(\textit{{k}}-\textit{{q}},\tau),$
(15)

where the coupling kernels $\tilde{\alpha}$ and $\tilde{\beta}$ respectively
are defined as

$\tilde{\alpha}(\textit{{q}}_{1},\textit{{q}}_{2})=\frac{\textit{{q}}_{1}.(\textit{{q}}_{1}+\textit{{q}}_{2})}{\textit{q}_{1}^{2}},$
(16)
$\tilde{\beta}(\textit{{q}}_{1},\textit{{q}}_{2})=\frac{1}{2}(\textit{{q}}_{1}+\textit{{q}}_{2})^{2}\frac{\textit{{q}}_{1}.\textit{{q}}_{2}}{\textit{q}_{1}^{2}\textit{q}_{2}^{2}}.$
(17)

For a gEdS cosmology, parametrized by $\kappa^{2}$, combining Eqs. (11) and
(12) with Eq. (9) to obtain a single second-order differential equation for
the density fluctuations $\delta(\textbf{x},\tau)$ and for the divergence of
velocity fluctuations $\theta(\textbf{x},\tau)$, respectively, we have

$\displaystyle\mathcal{H}^{2}\Big{\\{}-a^{2}\partial_{a}^{2}-\frac{3}{2}a\partial_{a}+\frac{3}{2\kappa^{2}}\Big{\\}}\delta(\textit{{k}},a)$
$\displaystyle=$ $\displaystyle
S_{\tilde{\beta}}(\textit{{k}},a)-\mathcal{H}\partial_{a}(aS_{\tilde{\alpha}}(\textit{{k}},a)),$
(18)
$\displaystyle\mathcal{H}\Big{\\{}a^{2}\partial_{a}^{2}+\frac{5}{2}a\partial_{a}+\Big{(}\frac{1}{2}-\frac{3}{2\kappa^{2}}\Big{)}\Big{\\}}\theta(\textit{{k}},a)$
$\displaystyle=$
$\displaystyle\partial_{a}(aS_{\tilde{\beta}}(\textit{{k}},a))-\frac{3}{2\kappa^{2}}\mathcal{H}S_{\tilde{\alpha}}(\textit{{k}},a),$
(19)

where $\partial_{a}$ is the partial derivative respect to scale factor $a$.

We now proceed to solve these equations following the standard method (see
e.g. [14]) in standard Eulerian perturbation theory (EPT). Working to linear
order in $\delta$ and $\theta$, Eq. (18) is just

$-a^{2}\partial_{a}^{2}\delta(\textit{{k}},a)-\frac{3}{2}a\partial_{a}\delta(\textit{{k}},a)+\frac{3}{2\kappa^{2}}\delta(\textit{{k}},a)=0,$
(20)

which has linearly independent decaying and growing mode solutions so that

$\begin{split}\delta(\textit{{k}},a)&=\delta_{+}(\textbf{k},a)+\delta_{-}(\textbf{k},a)\\\
&=D_{+}(a)\delta_{+,0}(\textbf{k})+D_{-}(a)\delta_{-,0}(\textbf{k}),\end{split}$
(21)

where $\delta_{\pm,0}(\textbf{k})$ are constants fixed at some reference time,
$a=1$, and $D_{\pm}=a^{\alpha_{\pm}}$ with

$\alpha_{\pm}=-\frac{1}{4}\pm\frac{1}{4}\sqrt{1+\frac{24}{\kappa^{2}}}.$ (22)

We assume the growing mode initial conditions, corresponding to the solution
at linear order which may be written

$\delta(\textit{{k}},a)=D(a)\delta^{(1)}(\textbf{k}),$ (23)

where $D(a)=a^{\alpha}$, and $\delta^{(1)}(\textbf{k})$ is the fluctuation at
$a=1$. At the same linear order we have the relation between density
fluctuations and velocity

$\theta(\textit{{k}},\tau)=-\delta^{\prime}(\textit{{k}},\tau)\,,$ (24)

from which the solution at linear order for the velocity divergence follows,

$\theta(\textit{{k}},a)=-\mathcal{H}(a)\alpha D(a)\delta^{(1)}(\textbf{k}).$
(25)

For the case $\kappa=1$ (and $\alpha=1$) we recover the standard expression
for the EdS model. These expressions are just special cases of the standard
ones for FLRW cosmologies, which correspond to Eq. (23) and

$\theta(\textit{{k}},a)=-\mathcal{H}(a)\alpha_{1}(a)D_{1}(a)\delta^{(1)}(\textbf{k}),$
(26)

where $D_{1}(a)$ is the appropriate growth factor for the linear theory
growing mode in the cosmology, and

$\alpha_{1}(a)=\frac{d\ln D_{1}}{d\ln a}.$ (27)

Taking leading order solutions as Eqs. (23) and (25), we now proceed to solve
Eqs. (18) and (19) perturbatively to higher orders taking the ansatz

$\begin{split}\delta(\textbf{k},a)&=\sum_{i=1}^{\infty}D^{i}(a)\,\delta^{(i)}(\textbf{k}),\\\
\theta(\textbf{k},a)&=-\mathcal{H}(a)\alpha\sum_{i=1}^{\infty}D^{i}(a)\,\theta^{(i)}(\textbf{k}),\end{split}$
(28)

where $D^{i}(a)\equiv(a^{\alpha})^{i}$, and $\delta^{(i)}$,
$\theta^{(i)}(\textbf{k})$ are the contributions at order $i$ to each of these
quantities. Such a “separable” ansatz means that all dependence on the
cosmological model disappears when $\delta(\textbf{k})$ and
$\theta(\textbf{k})$ are expressed in terms of their values at linear order.
Such an ansatz solves the equations exactly in perturbation for the EdS model
(see e.g. [14]). For the family of perturbed FLRW cosmologies with a perturbed
matter component, one can consider the same ansatz with $\alpha$ replaced by
$\alpha_{1}(a)$, and it has been shown [38, 39] that the condition for it be a
solution is that the ratio

$\frac{\Omega_{m}}{\alpha_{1}^{2}},$ (29)

is constant in time. While this condition is not satisfied in LCDM-like
models, it is in gEdS models (in which both $\Omega_{m}$ and $\alpha_{1}$ are
individually constant).

Using this ansatz (28) the n-th order solutions for the density and velocity
fields can be expressed in exactly the same form as in standard EdS as

$\displaystyle\delta^{(n)}(\textbf{k})$ $\displaystyle=$
$\displaystyle\prod_{m=1}^{n}\Big{\\{}\int\frac{d^{3}q_{m}}{(2\pi)^{3}}\delta^{(1)}(\textbf{q}_{m})\Big{\\}}$
(30)
$\displaystyle\times(2\pi)^{3}\delta_{D}(\textbf{k}-\textbf{q}|_{1}^{n})F_{n}(\textbf{q}_{1},\dots,\textbf{q}_{n}),$
$\displaystyle\theta^{(n)}(\textbf{k})$ $\displaystyle=$
$\displaystyle\prod_{m=1}^{n}\Big{\\{}\int\frac{d^{3}q_{m}}{(2\pi)^{3}}\delta^{(1)}(\textbf{q}_{m})\Big{\\}}$
(31)
$\displaystyle\times(2\pi)^{3}\delta_{D}(\textbf{k}-\textbf{q}|_{1}^{n})G_{n}(\textbf{q}_{1},\dots,\textbf{q}_{n}).$

Working to second-order in the expansion, we have

$\displaystyle\mathcal{H}^{2}\Big{\\{}-a^{2}\partial_{a}^{2}-\frac{3}{2}a\partial_{a}+\frac{3}{2\kappa^{2}}\Big{\\}}a^{2\alpha}\delta^{(2)}(\textit{{k}})$
(32) $\displaystyle=$
$\displaystyle-\int\frac{\text{d}^{3}q}{(2\pi)^{3}}\tilde{\beta}(\textit{{q}},\textit{{k}}-\textit{{q}})\delta^{(1)}(\textit{{q}})\delta^{(1)}(\textit{{k}}-\textit{{q}})\times\big{[}\mathcal{H}^{2}\alpha^{2}a^{2\alpha}\big{]}$
$\displaystyle-\int\frac{\text{d}^{3}q}{(2\pi)^{3}}\tilde{\alpha}(\textit{{q}},\textit{{k}}-\textit{{q}})\delta^{(1)}(\textit{{q}})\times\delta^{(1)}(\textit{{k}}-\textit{{q}})$
$\displaystyle\times\mathcal{H}^{2}\big{[}\frac{\alpha}{2}+2\alpha^{2}\big{]}a^{2\alpha},$

from which we obtain the second-order density kernel
$\textit{F}_{2}(\textit{{q}}_{1},\textit{{q}}_{2})$ as

$\Big{(}3\alpha^{2}+\frac{\alpha}{2}\Big{)}\textit{F}_{2}(\textit{{q}}_{1},\textit{{q}}_{2})=\alpha^{2}\tilde{\beta}(\textit{{q}}_{1},\textit{{q}}_{2})+\Big{(}\frac{\alpha}{2}+2\alpha^{2}\Big{)}\tilde{\alpha}(\textit{{q}}_{1},\textit{{q}}_{2}),$
(33)

and thus finally

$\textit{F}_{2}(\textit{{q}}_{1},\textit{{q}}_{2})=d_{2}\tilde{\alpha}(\textit{{q}}_{1},\textit{{q}}_{2})+(1-d_{2})\tilde{\beta}(\textit{{q}}_{1},\textit{{q}}_{2}),$
(34)

where

$d_{2}=\Big{(}\frac{1+4\alpha}{1+6\alpha}\Big{)}.$ (35)

For the velocity field divergence at the same order we have

$\displaystyle\mathcal{H}\Big{\\{}a^{2}\partial_{a}^{2}+\frac{5}{2}a\partial_{a}+(\frac{1}{2}-\frac{3}{2\kappa^{2}})\Big{\\}}\big{[}-\mathcal{H}\alpha
a^{2\alpha}\theta^{(2)}(\textit{{k}})\big{]}$ (36) $\displaystyle=$
$\displaystyle-\int\frac{\text{d}^{3}q}{(2\pi)^{3}}\tilde{\beta}(\textit{{q}},\textit{{k}}-\textit{{q}})\delta^{(1)}(\textit{{q}})\delta^{(1)}(\textit{{k}}-\textit{{q}})$
$\displaystyle\times\big{[}\mathcal{H}^{2}2\alpha^{3}a^{2\alpha}\big{]}-\int\frac{\text{d}^{3}q}{(2\pi)^{3}}\tilde{\alpha}(\textit{{q}},\textit{{k}}-\textit{{q}})$
$\displaystyle\times\delta^{(1)}(\textit{{q}})\delta^{(1)}(\textit{{k}}-\textit{{q}})\mathcal{H}^{2}\big{[}\frac{\alpha^{2}}{2}+\alpha^{3}\big{]}a^{2\alpha},$

which leads to the solution

$\textit{G}_{2}(\textit{{q}}_{1},\textit{{q}}_{2})=\tilde{d}_{2}\tilde{\alpha}(\textit{{q}}_{1},\textit{{q}}_{2})+(1-\tilde{d}_{2})\tilde{\beta}(\textit{{q}}_{1},\textit{{q}}_{2})$
(37)

where

$\tilde{d}_{2}=\Big{(}\frac{1+2\alpha}{1+6\alpha}\Big{)}.$ (38)

Continuing to third-order in the same manner we obtain

$\displaystyle\textit{F}_{3}(\textit{{q}}_{1},\textit{{q}}_{2},\textit{{q}}_{3})$
$\displaystyle=$
$\displaystyle\frac{1}{2}\bigg{[}d_{3}\tilde{\alpha}(\textit{{q}}_{1},\textit{{q}}_{2}+\textit{{q}}_{3})\textit{F}_{2}(\textit{{q}}_{2},\textit{{q}}_{3})$
(39)
$\displaystyle+(1-d_{3})\tilde{\beta}(\textit{{q}}_{1},\textit{{q}}_{2}+\textit{{q}}_{3})\textit{G}_{2}(\textit{{q}}_{2},\textit{{q}}_{3})$
$\displaystyle+\big{[}d_{3}\tilde{\alpha}(\textit{{q}}_{1}+\textit{{q}}_{2},\textit{{q}}_{3})$
$\displaystyle+(1-d_{3})\tilde{\beta}(\textit{{q}}_{1}+\textit{{q}}_{2},\textit{{q}}_{3})\big{]}$
$\displaystyle\times\textit{G}_{2}(\textit{{q}}_{1},\textit{{q}}_{2})\bigg{]},$

and

$\displaystyle\textit{G}_{3}(\textit{{q}}_{1},\textit{{q}}_{2},\textit{{q}}_{3})$
$\displaystyle=$
$\displaystyle\frac{1}{2}\bigg{[}\tilde{d}_{3}\tilde{\alpha}(\textit{{q}}_{1},\textit{{q}}_{2}+\textit{{q}}_{3})\textit{F}_{2}(\textit{{q}}_{2},\textit{{q}}_{3})$
(40)
$\displaystyle+(1-\tilde{d}_{3})\tilde{\beta}(\textit{{q}}_{1},\textit{{q}}_{2}+\textit{{q}}_{3})\textit{G}_{2}(\textit{{q}}_{2},\textit{{q}}_{3})$
$\displaystyle+\big{[}\tilde{d}_{3}\tilde{\alpha}(\textit{{q}}_{1}+\textit{{q}}_{2},\textit{{q}}_{3})$
$\displaystyle+(1-\tilde{d}_{3})\tilde{\beta}(\textit{{q}}_{1}+\textit{{q}}_{2},\textit{{q}}_{3})\big{]}$
$\displaystyle\times\textit{G}_{2}(\textit{{q}}_{1},\textit{{q}}_{2})\bigg{]},$

where

$d_{3}=\frac{1+6\alpha}{1+8\alpha},\quad\tilde{d}_{3}=\frac{1+2\alpha}{1+8\alpha}.$
(41)

Higher-order kernels can be inferred by continuing this recursion.111Further
details of these calculations are given in Appendix A.

Figure 1: Variation of the functions $d_{2}$ and $d_{3}$ determining the EPT
kernels at second and third order relative to their values in EdS, as a
function of the logarithmic growth rate $\alpha$ of gEdS. The dashed
horizontal lines bracket approximately the region in which the effective
growth rate $\alpha_{1}(a)$ varies in LCDM-like models.

Fig. 1 shows how the kernels $F_{2}$ and $G_{2}$ vary relative to those in the
standard EdS case as a function of $\alpha$ (and $F_{3}$ and $G_{3}$ have very
similar behaviours). We see that the dependence on $\alpha$ is very weak,
other than at values of $\alpha$ very much less than unity. Indeed expanding
about $\alpha=1$ we have

$\displaystyle
d_{2}=d_{2}(1)\left[1+\frac{2}{35}(1-\alpha)+O((1-\alpha)^{2})\right],$
$\displaystyle
d_{3}=d_{3}(1)\left[1+\frac{2}{63}(1-\alpha)+O((1-\alpha)^{2})\right],$ (42)

and, over the whole range of $\alpha$, varying from $0$ to $\infty$, we have

$\displaystyle\frac{7}{5}>\frac{d_{2}}{d_{2}(1)}>\frac{14}{15},$ (43)
$\displaystyle\frac{9}{7}>\frac{d_{3}}{d_{3}(1)}>\frac{27}{28}.$ (44)

A naive guess would be that, in a generic FLRW cosmology, the kernel at any
time can be approximated by those in a gEdS model, with the growth rate
corresponding to the instantaneous value of its effective logarithmic growth
rate $\alpha_{1}(a)$. For standard-type cosmological models, as we will
discuss in detail in the next section, such an effective growth rate decreases
relative to EdS to at most $\alpha\approx 0.5$ at $z=0$ as dark energy
contributions accelerate the expansion. From Fig. 1 we see that this would
correspond to a change of about $5\%$ ($3\%$) in $d_{2}$($d_{3}$) relative to
EdS. In practice, as we will see in the second part of the paper, such a
mapping to an approximate mapping may indeed to used, but the relevant
effective growth exponent is given by this naive guess only for a sufficiently
slowly varying smooth component. For LCDM-like models it is much closer to the
initial value. The very small sub-percent corrections calculated for LCDM-like
models can thus be understood as a combination of the very weak sensitivity of
the gEdS kernels to the growth rate of fluctuations, and the rapid evolution
of dark energy at low redshift.

###  II.3 Lagrangian Perturbation Theory (LPT)

In LPT one characterizes the evolution of the fluid using its displacement
field $\Psi(\textbf{q},t)$ as a function of its initial Lagrangian position q
of the particle, in terms of which the Eulerian position is given by

$\textbf{x}(\textbf{q},\tau)=\textbf{q}+\mathbf{\Psi}(\textbf{q},\tau).$ (45)

It follows in particular that

$\big{[}1+\delta(\textbf{x})\big{]}\text{d}^{3}\text{x}=\text{d}^{3}\text{q},$
(46)

and that the determinant $J$ of the Jacobian matrix is given by

$J=\Big{|}\frac{d^{3}x}{d^{3}q}\Big{|}=\text{det}\Big{[}\delta_{ij}^{(K)}+\Psi_{i,j}\Big{]},$
(47)

where

$\Psi_{i,j}=\frac{\partial\Psi_{i}}{\partial q_{j}},$ (48)

so that

$\delta(\textbf{x})=\frac{1}{J}-1.$ (49)

The equation of motion for the fluid particle is

$\frac{\text{d}^{2}\mathbf{\Psi}}{\text{d}\tau^{2}}+\mathcal{H}\frac{\text{d}\mathbf{\Psi}}{\text{d}\tau}=-\nabla_{x}\phi.$
(50)

Taking the divergence of this equation, and using the Poisson Eq. (9), one
obtains

$\begin{split}J\bm{\nabla}_{\textbf{x}}\cdot\Big{[}\frac{\text{d}^{2}\mathbf{\Psi}}{\text{d}\tau^{2}}+\mathcal{H}\frac{\text{d}\mathbf{\Psi}}{\text{d}\tau}\Big{]}&=\frac{3}{2}\Omega_{m}\mathcal{H}^{2}[J-1].\end{split}$
(51)

This equation order is then solved order by order using the perturbative
expansion for

$\mathbf{\Psi}=\mathbf{\Psi}^{(1)}+\mathbf{\Psi}^{(2)}+\mathbf{\Psi}^{(3)}+\dots$
(52)

and

$J=1+J^{(1)}+J^{(2)}+J^{(3)}+\dots$ (53)

where the latter are related to the former up to third order by the following
expressions (for a detailed derivation, see e.g. [40]):

$\begin{split}J^{(1)}&=\sum_{i}\Psi_{i,i}^{(1)},\\\
J^{(2)}&=\sum_{i}\Psi_{i,i}^{(2)}+\frac{1}{2}\sum_{i\neq
j}\Big{\\{}\Psi_{i,i}^{(1)}\Psi_{j,j}^{(1)}-\Psi_{i,j}^{(1)}\Psi_{j,i}^{(1)}\Big{\\}},\\\
J^{(3)}&=\sum_{i}\Psi_{i,i}^{(3)}+\sum_{i\neq
j}\Big{\\{}\Psi_{i,i}^{(2)}\Psi_{j,j}^{(1)}-\Psi_{i,j}^{(2)}\Psi_{j,i}^{(1)}\Big{\\}}+\text{det}\Psi_{i,j}^{(1)}.\end{split}$
(54)

Using the chain rule

$\nabla_{x_{i}}=(\delta_{ij}^{(K)}+\Psi_{i,j})^{-1}\nabla_{q_{j}},$ (55)

one then solves Eq. (51) order by order. At linear order we have simply

$\frac{\text{d}^{2}\Psi_{i,i}^{(1)}}{\text{d}\tau^{2}}+\mathcal{H}\frac{\text{d}\Psi_{i,i}^{(1)}}{\text{d}\tau}-\frac{3}{2}\Omega_{m}\mathcal{H}^{2}\Psi_{i,i}^{(1)}=0,$
(56)

and thus, choosing the growing mode solution,

$\Psi_{i,i}^{(1)}(\textbf{x},\tau)=-\delta^{(1)}(\textbf{x},\tau),$ (57)

which gives in Fourier space

$\Psi^{(1)}(\textbf{k},\tau)=-\frac{i\textbf{k}}{k^{2}}\delta^{(1)}(\textbf{k})D_{1}(\tau).$
(58)

This is the solution for any perturbed FLRW cosmology and thus for gEdS in
particular with $D_{1}(a)=D(a)=a^{\alpha}$.

Considering now gEdS cosmologies, we proceed to solve Eq. (51) to non-linear
order making the separable ansatz analogous to that in EPT:

$\Psi^{(n)}(\textbf{k},\tau)=D^{n}(\tau)\Psi^{(n)}(\textbf{k})\,.$ (59)

We have then

$\begin{split}&(1+J^{(1)}+J^{(2)}+J^{(3)}+\dots)(\delta_{ij}^{(K)}-\Psi_{i,j}+\Psi_{i,l}\Psi_{l,j}+\dots)\\\
&\times
n\alpha\Big{[}n\alpha+\frac{1}{2}\Big{]}\Psi_{i,j}=\frac{3}{2\kappa^{2}}[J^{(1)}+J^{(2)}+J^{(3)}+\dots],\end{split}$
(60)

where we have used
$\frac{\text{d}\mathcal{H}}{\text{d}a}=-\frac{\mathcal{H}}{2a}$. At second
order we have therefore

$\Psi_{i,i}^{(2)}=-\frac{2\alpha+1}{2(6\alpha+1)}\sum_{i,j}\big{[}\Psi_{i,i}^{(1)}\Psi_{j,j}^{(1)}-\Psi_{i,j}^{(1)}\Psi_{i,j}^{(1)}\Big{]}.$
(61)

The result at third-order, of which further details are given in Appendix A.2,
is

$\displaystyle\Psi_{i,i}^{(3)}$ $\displaystyle=$
$\displaystyle\frac{1}{(8\alpha+1)}\Big{[}(4\alpha+1)\Big{\\{}-\Psi_{i,i}^{(2)}\Psi_{j,j}^{(1)}+\Psi_{i,j}^{(2)}\Psi_{j,i}^{(1)}\Big{\\}}$
(62)
$\displaystyle+(2\alpha+1)\Big{\\{}\frac{1}{2}\Psi_{i,i}^{(1)}\Psi_{i,j}^{(1)}\Psi_{j,i}^{(1)}$
$\displaystyle-\frac{1}{6}\Psi_{i,i}^{(1)}\Psi_{j,j}^{(1)}\Psi_{k,k}^{(1)}-\frac{1}{3}\Psi_{i,k}^{(1)}\Psi_{k,j}^{(1)}\Psi_{j,i}^{(1)}\Big{\\}}\Big{]}.$

At n-th order the displacement field in Fourier space can be written in the
form

$\displaystyle\Psi^{(n)}(\textbf{k})$ $\displaystyle=$
$\displaystyle-\frac{i}{n}\prod_{i=1}^{n}\Big{\\{}\int\frac{\text{d}^{3}q_{i}}{(2\pi)^{3}}\delta^{(1)}(\textbf{q}_{i})\Big{\\}}\textbf{L}_{n}(\textbf{q}_{1},\dots,\textbf{q}_{n})$
(63)
$\displaystyle\times(2\pi)^{3}\delta^{(D)}(\textbf{k}-\textbf{q}_{1}\dots\textbf{q}_{n}),$

where the first, second, and third-order solutions correspond to the following
kernels:

$\displaystyle\textbf{L}_{1}$ $\displaystyle=$
$\displaystyle\frac{\textbf{k}}{k^{2}},$ (64) $\displaystyle\textbf{L}_{2}$
$\displaystyle=$
$\displaystyle\Big{(}\frac{2\alpha+1}{6\alpha+1}\Big{)}\frac{\textbf{k}}{k^{2}}\Big{[}1-\frac{(\textbf{q}_{1}\cdot\textbf{q}_{2})^{2}}{q_{1}^{2}q_{2}^{2}}\Big{]},$
(65) $\displaystyle\textbf{L}_{3}$ $\displaystyle=$ $\displaystyle
3\Big{(}\frac{4\alpha+1}{8\alpha+1}\Big{)}\Big{(}\frac{2\alpha+1}{6\alpha+1}\Big{)}\frac{\textbf{k}}{k^{2}}\Big{[}1-\Big{(}\frac{\textbf{q}_{1}\cdot\textbf{q}_{2}}{q_{1}q_{2}}\Big{)}^{2}\Big{]}$
(66)
$\displaystyle\times\Big{[}1-\Big{(}\frac{(\textbf{q}_{1}+\textbf{q}_{2})\cdot\textbf{q}_{3}}{|(\textbf{q}_{1}+\textbf{q}_{2}|q_{3}}\Big{)}^{2}\Big{]}-\Big{(}\frac{2\alpha+1}{8\alpha+1}\Big{)}\frac{\textbf{k}}{k^{2}}$
$\displaystyle\times\Big{[}1-3\Big{(}\frac{\textbf{q}_{1}\cdot\textbf{q}_{2}}{q_{1}q_{2}}\Big{)}^{2}$
$\displaystyle+2\frac{(\textbf{q}_{1}\cdot\textbf{q}_{2})(\textbf{q}_{2}\cdot\textbf{q}_{3})(\textbf{q}_{3}\cdot\textbf{q}_{1})}{q_{1}^{2}q_{2}^{2}q_{3}^{2}}\Big{]}.$

As a check on our calculation we can use the known relations between the EPT
and LPT kernels (see [29, 41, 42]). At first-order we have simply

$F_{1}(\textbf{k})=\textbf{k}\cdot\textbf{L}_{1}(\textbf{k})=\textbf{k}\cdot\frac{\textbf{k}}{k^{2}}=1.$
(67)

The second-order density kernels are related by

$\displaystyle F_{2}(\textbf{q}_{1},\textbf{q}_{2})$ $\displaystyle=$
$\displaystyle\frac{1}{2}\Big{[}\textbf{k}\cdot\textbf{L}_{2}(\textbf{q}_{1},\textbf{q}_{2})+[\textbf{k}\cdot\textbf{L}_{1}(\textbf{q}_{1}][\textbf{k}\cdot\textbf{L}_{1}(\textbf{q}_{2})]\Big{]}$
(68) $\displaystyle=$
$\displaystyle\frac{4\alpha+1}{6\alpha+1}+\frac{1}{2}\frac{\textit{{q}}_{1}\cdot\textit{{q}}_{2}}{q_{1}q_{2}}\Big{[}\frac{q_{2}}{q_{1}}+\frac{q_{1}}{q_{2}}\Big{]}$
$\displaystyle+\Big{(}\frac{2\alpha}{6\alpha+1}\Big{)}\frac{(\textit{{q}}_{1}\cdot\textit{{q}}_{2})^{2}}{q_{1}^{2}q_{2}^{2}},$

while at third order we have

$\displaystyle F_{3}(\textbf{q}_{1},\textbf{q}_{2},\textbf{q}_{3})$
$\displaystyle=$
$\displaystyle\frac{1}{3!}\Big{[}\textbf{k}\cdot\textbf{L}_{3}^{(s)}(\textbf{q}_{1},\textbf{q}_{2},\textbf{q}_{3})$
(69)
$\displaystyle+\big{\\{}[\textbf{k}\cdot\textbf{L}_{1}(\textbf{q}_{1})][\textbf{k}\cdot\textbf{L}_{2}(\textbf{q}_{2},\textbf{q}_{3})]$
$\displaystyle+\text{cyc}\big{\\}}$
$\displaystyle+[\textbf{k}\cdot\textbf{L}_{1}(\textbf{q}_{1})][\textbf{k}\cdot\textbf{L}_{1}(\textbf{q}_{2})]$
$\displaystyle\times[\textbf{k}\cdot\textbf{L}_{1}(\textbf{q}_{3})]\Big{]}.$

It is straightforward to check that we indeed recover the results in EPT of
the previous subsection.

##  III Power spectrum in gEdS models

We define the power spectrum $P(\vec{k})\equiv P(k)$ ($k=|\textbf{k}|$) of the
(assumed) statistically homogeneous and isotropic stochastic density field by

$\langle\delta(\textbf{k},a)\delta(\textbf{k}^{\prime},a)\rangle=(2\pi)^{3}\delta^{(D)}(\textbf{k}+\textbf{k}^{\prime})P(|\textbf{k}|,a),$
(70)

where $\langle\cdots\rangle$ denotes the ensemble average. Assuming that the
fluctuations are Gaussian at linear order, one obtains, retaining terms up to
fourth order in the linear order solution,

$\displaystyle\langle\delta(\textbf{k},a)\delta(\textbf{k}^{\prime},a)\rangle$
$\displaystyle=$
$\displaystyle\langle\delta^{(1)}(\textbf{k},a)\delta^{(1)}(\textbf{k}^{\prime},a)\rangle$
(71)
$\displaystyle+2\langle\delta^{(1)}(\textbf{k},a)\delta^{(3)}(\textbf{k}^{\prime},a)\rangle$
$\displaystyle+\langle\delta^{(2)}(\textbf{k},a)\delta^{(2)}(\textbf{k}^{\prime},a)\rangle,$

and thus the expression for the “one-loop” power spectrum as

$\textit{P}_{1-\text{loop}}(k,a)=P_{\text{\rm
lin}}(k,a)+P_{22}(k,a)+2P_{13}(k,a),$ (72)

where $P_{\text{lin}}(k,a)$ is the linear power spectrum (i.e. defined by Eq.
(70) with $\delta(\vec{k},a)=D(a)\delta^{(1)}(\vec{k}$), and the two other
terms correspond to

$\displaystyle P_{22}(k,a)$ $\displaystyle=$ $\displaystyle
2\int\frac{d^{3}q}{(2\pi)^{3}}P_{\text{lin}}(q,a)P_{\text{lin}}(|\textbf{k}-\textbf{q}|,a)$
(73) $\displaystyle\times|F_{2}^{(s)}(\textbf{k}-\textbf{q},\textbf{q})|^{2}.$

and

$\displaystyle P_{13}(k,a)$ $\displaystyle=$ $\displaystyle
3P_{\text{lin}}(k,a)\int\frac{d^{3}q}{(2\pi)^{3}}P_{\text{lin}}(q,a)$ (74)
$\displaystyle\times F_{3}^{(s)}(\textbf{k},\textbf{q},-\textbf{q}),$

where $F_{2}^{(s)}$ and $F_{3}^{(s)}$ are the symmetrized form of $F_{2}$ and
$F_{3}$ (with respect to permutation of their arguments). These expressions
are formally identical to those in standard EdS, and share the property of any
separable solution that the perturbative corrections to the PS at any time are
functionals only of the linear power spectrum at that time. The modifications
of gEdS relative to the standard EdS model are expressed solely through the
$\alpha$-dependence of the kernels $F_{2}$ and $F_{3}$ which we have derived
above. Using these to make the $\alpha$-dependence explicit we obtain

$\displaystyle P_{22}$ $\displaystyle=$ $\displaystyle
M_{0}+\frac{1+4\alpha}{1+6\alpha}M_{1}+\left(\frac{1+4\alpha}{1+6\alpha}\right)^{2}M_{2},$
(75) $\displaystyle 2P_{13}$ $\displaystyle=$ $\displaystyle
N_{0}+\frac{1+2\alpha}{1+8\alpha}N_{1}+\frac{2\alpha(1+2\alpha)}{(1+6\alpha)(1+8\alpha)}N_{2},$
(76)

where the $M_{i}$ are the integrals

$\displaystyle M_{i}(k,a)$ $\displaystyle=$
$\displaystyle\frac{1}{8\pi^{2}}k^{3}\int_{0}^{\infty}drP_{\text{lin}}(kr,a)$
(77)
$\displaystyle\times\int_{-1}^{1}d\mu\frac{P_{\text{lin}}(k\sqrt{1+r^{2}-2\mu
r},a)}{(1+r^{2}-2\mu r)^{2}}\,$ $\displaystyle\times m_{i}(r,\mu),$

with

$\displaystyle m_{0}(r,\mu)$ $\displaystyle=$ $\displaystyle(\mu-r)^{2},$ (78)
$\displaystyle m_{1}(r,\mu)$ $\displaystyle=$ $\displaystyle
4r(\mu-r)(1-\mu^{2}),$ (79) $\displaystyle m_{2}(r,\mu)$ $\displaystyle=$
$\displaystyle 4r^{2}(1-\mu^{2})^{2},$ (80)

and $N_{i}$ are the integrals

$N_{i}(k,a)=\frac{1}{8\pi^{2}}k^{3}P_{\rm lin}(k,a)\int_{0}^{\infty}drP_{\rm
lin}(kr,a)n_{i}(r)$ (81)

with

$\begin{split}n_{0}(r)&=-\frac{4}{3},\\\
n_{1}(r)&=1+\frac{8}{3}r^{2}-r^{4}+\frac{(r^{2}-1)^{3}}{2r}\ln\frac{|1+r|}{|1-r|},\\\
n_{2}(r)&=\frac{1}{r^{2}}(1-\frac{8}{3}r^{2}-r^{4})+\frac{(r^{2}-1)^{3}}{2r^{3}}\ln\frac{|1+r|}{|1-r|},\end{split}$
(82)

and we have defined $\mu=\textbf{k}\cdot\textbf{q}/kq$ and $r=q/k$.

###  III.1 Convergence properties of the one-loop power spectrum

| IR-divergent | UV-divergent
---|---|---
$M_{0}$ | $n\leq-1$ | $n\geq 1/2$
$M_{1}$ | $n\leq-3$ | $n\geq 1/2$
$M_{2}$ | $n\leq-3$ | $n\geq 1/2$
$N_{0}$ | $n\leq-1$ | $n\geq-1$
$N_{1}$ | $n\leq-3$ | $n\geq-1$
$N_{2}$ | $n\leq-3$ | $n\geq-1$
$N_{0}+M_{0}$ | $n\leq-3$ | $n\geq-1$
Table 1: Convergence properties of the integrals $M_{i}$ and $N_{i}$
($i=0,1,2$). In each case the bound on $n$ is that obtained assuming that
$P_{\rm lin}\sim k^{n}$ in the relevant ($k\rightarrow 0$ or
$k\rightarrow\infty$) limit. Integrals are “infra-red safe” if they converge
for $n>-3$, i.e. when $P_{\rm lin}(k)$ itself is integrable (in three
dimensions) as $k\rightarrow 0$. The cancellation of the divergences for
$-1\geq n>-3$ in the sum $N_{0}+M_{0}$ corresponds to the well-known
cancellation between the full $P_{22}$ and $P_{13}$ contributions in an EdS
cosmology. As expected this cancellation generalizes to the gEdS case, but
without the need for any cancellation in the $\alpha$-dependent terms. We will
show below that the corrections in standard (e.g. LCDM-like) cosmologies to
the one-loop PS are expressed in terms of the four infra-red safe integrals
$M_{1},M_{2},N_{1},N_{2}$.

We now consider the convergence properties of these one-loop contributions to
the PS in the gEdS model i.e. for what asymptotic behaviours of the linear PS,
$P_{\text{lin}}(k,a)$, the integrals converge. The results of this analysis
are summarized in Table 1. These are obtained straightforwardly by considering
the $r\rightarrow 0$, and $r\rightarrow\infty$, behaviours of the integrands
for the infra-red and ultra-violet limits, respectively. In the $M_{i}$
integrals there is also a divergence in the angular integral at $\mu=1$ for
$r=1$, corresponding to $|\bf{q}-\bf{k}|\rightarrow 0$, but because of the
symmetry of exchange of the variables $\bf{k}$ and $\bf{q}-\bf{k}$ in the
integral (as noted e.g. in [13]), its contribution is identical to that from
$r\rightarrow 0$ and is thus taken into account by multiplying it by a factor
of two.

The leading contribution to $M_{0}$ as $r\rightarrow 0$ corresponds to
$m_{0}=\mu^{2}$, which when integrated over angle, and multiplied by two,
gives a factor of $4/3$ which can be seen to exactly cancel the corresponding
leading contribution from $N_{0}$. These contributions to each integral are
proportional to $\int P_{\rm lin}(k)dk$ i.e. to the variance of the
displacement field. As understood in early works on perturbation theory (e.g.
[6, 43], see also [44] for a useful discussion), these are unphysical
divergences if the system (as is the case here, just as in standard EdS) is
Galilean invariant, and their cancellation is a consequence of this
invariance. The next term in the expansion as $r\rightarrow 0$ does not cancel
but is proportional to $\int P_{\rm lin}(k)k^{2}dk$ i.e. to the variance of
the fluctuation field (which is Galilean invariant). Their sum $M_{0}+N_{0}$
is then said to be “infra-red safe” in this case: the integral converges for
any $P_{\text{lin}}(k,a)$ which is integrable at $k\rightarrow 0$.

We see, on the other hand, that the integrals $M_{1},M_{2},N_{1}$ and $N_{2}$
are all infra-red safe: for $r\rightarrow 0$, the functions
$m_{1},m_{2},n_{1},n_{2}$ all have, after integration over angles, a leading
behaviour $\sim r^{2}$. This means that infra-red safety of all the
$\alpha$-dependent contributions is obtained term by term, without any
cancellation between different terms. This is evidently a necessary condition
to obtain infra-red safety given that the $\alpha$-dependence in the $F_{2}$
and $F_{3}$ kernels are given through independent non-linear functions of
$\alpha$. As we will see below it turns out that we can write the cosmological
corrections (due to departures from EdS) of the one-loop PS in LCDM-like
models in terms of the same integrals. Infra-red safety of these cosmological
corrections is thus explicit and their numerical calculation simplified,
without the additional manipulations needed for the leading EdS approximated
contribution (see e.g. [45]).

Although we will not discuss them further in this paper, it is interesting to
consider also the ultraviolet divergences. These are physical divergences
marking a real breakdown of standard perturbation theory. It is
straightforward to show (see Appendix B for details) that combining the
leading contributions in the six integrals above we obtain

$P_{22}(k,a)=\frac{7+36\alpha+92\alpha^{2}}{30(1+6\alpha)^{2}}k^{4}\int\frac{d^{3}q}{(2\pi)^{3}}\frac{P_{\text{lin}}^{2}(q)}{q^{4}}+\cdots,$
(83)

and

$\displaystyle 2P_{13}(k)$ $\displaystyle=$
$\displaystyle\frac{7-14\alpha-176\alpha^{2}}{15(1+6\alpha)(1+8\alpha)}k^{2}P_{\text{lin}}(k)\int\frac{d^{3}q}{(2\pi)^{3}}\frac{P_{\text{lin}}(q)}{q^{2}}$
(84)
$\displaystyle+\frac{12(8\alpha-1)(1+2\alpha)}{105(1+6\alpha)(1+8\alpha)}k^{4}P_{\text{lin}}(k)$
$\displaystyle\times\int\frac{d^{3}q}{(2\pi)^{3}}\frac{P_{\text{lin}}(q)}{q^{4}}+\cdots$

where the dots now indicate terms proportional to integrals which converge
more rapidly in the ultra-violet than these leading terms. The first
expression corresponds to the sum of the leading contribution in the integrals
$M_{i}$, which we see (as indicated in Table 1) diverge for $n\geq 1/2$, while
the second shows the sum of the leading contributions from the $N_{i}$, the
first diverging for $n\geq-1$ (as indicated in Table 1) and the next one
diverging for $n\geq 1$.

The integrals in these expressions are the same as those in the standard EdS
models, the only difference being in the coefficients which are manifestly all
explicitly $\alpha$ dependent. The leading divergence overall is that in
$P_{13}$, diverging for $n\geq-1$. We note that this remains true in the
family of gEdS models, except that there is a specific value,
$\alpha=0.1635\cdots$, at which the coefficient of this leading term,
proportional to $7-14\alpha-176\alpha^{2}$, vanishes. For this specific value
of $\alpha$ therefore the one-loop PS is ultra-violet convergent for $n<1/2$.
We will not explore further in this article the physical significance of this
result, nor more generally the $\alpha$ dependence of the ultra-violet
divergences. To address these questions it is necessary to consider the
regularisation of these divergences, using for example, the RPT or EFT
approaches (see references cited in the introduction above). We will explore
these issues in future work.

##  IV LCDM approximated as an interpolation of gEdS

In the rest of this article we explore the relation of these models to
standard (e.g. LCDM-like) cosmological models. The class of models considered
are the FLRW models described at the beginning of Section II.1, i.e., with a
clustering matter component and a smooth component with a generic equation of
state. Further to this we will assume only that cosmology is matter dominated
at early times, i.e. that it is EdS, or gEdS if there remains a smooth matter-
like component (e.g. massive neutrinos or matter-like dark energy).
Specifically we explore whether we can make use of the analytical results for
perturbation theory in gEdS to calculate, exactly or approximately, results in
standard models (see references further below). More precisely we wish to
consider whether the functions, $Q(z)$ say, characterising the evolution of
fluctuations in perturbation theory in an LCDM model as a function of redshift
$z$ may be interpolated on the family of gEds models, so that we can write

$Q(z)\approx Q_{gEdS}(\alpha=\alpha_{eff}(z)),$ (85)

where $\alpha_{eff}(z)$ is an effective growth exponent which can be
calculated given the parameters of the LCDM model. Further we will determine
whether the approximation

$Q(z)\approx Q_{gEdS}(\kappa^{2}=1/\Omega_{m}(z)),$ (86)

is valid, i.e., whether we can take

$\alpha_{eff}(z)=-\frac{1}{4}+\frac{1}{4}\sqrt{1+24\Omega_{m}(z)},$ (87)

which means that we match both the instantaneous expansion rate and growth
rate of the FLRW cosmology to that of a gEdS cosmology. We will say in this
case that the interpolation is adiabatic, because we expect such an
approximation to be valid when we can neglect the temporal variation of the
smooth component relative to that of matter. Indeed we will see in all our
equations below for effective growth exponents that the approximation Eq. (87)
becomes exact as $w_{s}\rightarrow 0$, which corresponds to this limit.

We will apply here this analysis up to third order in perturbation theory,
which is what is needed to calculate the one-loop power spectrum. We have
chosen for our analysis to follow closely that of [27], who have calculated
explicitly the corrections in LCDM and provided phenomenological fits to the
relevant redshift dependent functions for standard type models. As we will
show, our method of analysis in seeking to relate these results to those in
gEdS, leads to a simplification of the formulation and calculations of the
results of [27], and other treatments given in the literature (e.g. [11, 28,
20, 32]) Further in so doing it provides greater insight into the reason why
these cosmology-dependent corrections are typically so small, and how they
depend on cosmological parameters. We compare our numerical results in detail
to the fits of [27] in Appendix D. For completeness we detail also a
comparison of our analysis, and numerical results when possible, with that of
[32] in Appendix E, and with those of [11] and [20] in Appendix F.

###  IV.1 Linear growth rate interpolated on gEdS

By combining Eqs. (11) and (12) at linear order (i.e. neglecting the source
terms $S_{\tilde{\alpha}}$ and $S_{\tilde{\beta}}$) we obtain easily the
equation obeyed by the linear growth factor $D_{1}(a)$ which can be cast as

$\displaystyle\frac{d^{2}}{d^{2}\ln a^{2}}\frac{D_{1}}{a}$ $\displaystyle+$
$\displaystyle\Big{(}4+\frac{d\ln H}{d\ln a}\Big{)}\frac{d}{d\ln
a}\frac{D_{1}}{a}$ (88) $\displaystyle+$ $\displaystyle\Big{(}3+\frac{d\ln
H}{d\ln a}-\frac{3}{2}\Omega_{m}(a)\Big{)}\frac{D_{1}}{a}=0.$

Writing

$D_{1}=\exp\left(\int_{\ln a_{0}}^{\ln a}\alpha_{1}(a)\,d(\ln a_{0})\right),$
(89)

where $a_{0}$ to be the reference time at which $D_{1}=1$, and

$\alpha_{1}=\frac{d\ln D_{1}}{d\ln a},$ (90)

we have, trivially, an interpolation in the sense of Eq. (85) taking
$\alpha_{eff}(z)=\alpha_{1}$.

Let us now consider the accuracy of an adiabatic gEdS interpolation, in the
sense defined by Eq. (86) above, i.e., the approximation

$\alpha_{1}(z)\approx\alpha_{10}(z)=\frac{1}{4}\big{[}-1+\sqrt{1+24\Omega_{m}(z)}\big{]}.$
(91)

Fig. 2 shows a comparison of the exact growth rate $\alpha_{1}$ (obtained by
solving Eq. (88) numerically) and $\alpha_{10}$, for models with constant
equation of state $w_{s}=w_{0}$, for $w_{0}=-0.5,-1,-1.5$. The left panel
shows the two quantities, and the right panel their ratio. We see that, even
for the LCDM case ($w_{0}=-1$), the approximation is quite good down to $z=0$
(when $\Omega_{m}\approx 0.3$).

Figure 2: Left panel: Exact solution for the growth rate $\alpha_{1}$ compared
with its “adiabatic gEdS approximation” $\alpha_{10}$, as a function of
$\Omega_{m}$ and for different indicated values of $w_{0}$, parametrizing the
constant equation of state of the smooth (dark energy) component. Right panel:
the relative difference of the same quantities.

This result is somewhat surprising as one would anticipate that it would hold
only for $|w_{0}|$ small compared to unity. To understand what is observed
better, we re-express Eq. (88) as an equation for $\alpha_{1}$:

$\left[\frac{d\alpha_{1}}{d\ln
a}-\frac{3}{2}w\alpha_{1}\right]+\alpha_{1}^{2}+\frac{1}{2}\alpha_{1}-\frac{3}{2}\Omega_{m}=0,$
(92)

where

$w=-1-\frac{2}{3}\frac{d\ln H}{d\ln a}=\frac{1}{3}\frac{d\ln\Omega_{m}}{d\ln
a},$ (93)

and for the specific case of a smooth (dark energy) component with the
constant equation of state, $w_{s}=w_{0}$, we have

$w=w_{0}(1-\Omega_{m}).$ (94)

The approximate solution Eq. (91) corresponds to neglecting the first two
terms in square brackets i.e., setting $w=0$ and neglecting the time
derivative, and taking the growing mode solution of the remaining quadratic
equation. It is thus formally the solution at leading order in a gradient
expansion in $\Omega_{m}(z)$ for which $w$ (or $w_{0}$ for a constant equation
of state) is the control parameter. Given that e.g. in LCDM $w=-0.7$ when
$\Omega_{m}=0.3$, it is indeed, as we have noted, surprising that the leading
$w=0$ approximation is so good. To see why this is the case more explicitly we
define

$\alpha_{1}=\alpha_{10}(1+\epsilon_{1}(z))\,.$ (95)

and rewrite Eq. (92), neglecting the single quadratic term in $\epsilon_{1}$,
as a linearized equation in $\epsilon_{1}$:

$\frac{d\epsilon_{1}}{d\ln
a}+\Big{[}-S(\Omega_{m})w+(\frac{1}{2}+2\alpha_{10})\Big{]}\epsilon_{1}=S(\Omega_{m})w,$
(96)

where

$S(\Omega_{m})=\frac{3}{2}-\frac{1}{w}\frac{d\ln\alpha_{10}}{d\ln
a}=\frac{3}{2}-\frac{9\Omega_{m}}{\alpha_{10}(1+4\alpha_{10})},$ (97)

is the coefficient multiplying $w$ which gives the source term for
$\epsilon_{1}$, the correction to the leading adiabatic gEdS approximation.

Why $\epsilon_{1}$ remains small for relatively large $w$ is just a
consequence of the smallness of the source term $S$ relative to the
coefficient of the $\epsilon_{1}$ term on the left-hand size of the equation.
This is shown explicitly in Fig. 3 which shows the numerical value of
$S(\Omega_{m})$ plotted against that of $(\frac{1}{2}+2\alpha_{10})$.

Figure 3: The coefficient $S(\Omega_{m})$ of the source term in Eq. (96) as a
function of $(1/2+2\alpha_{10})$, and as a function of $\Omega_{m}$ (upper
$x$-axis). We see that the ratio of $|S(\Omega_{m})|/(1/2+2\alpha_{10})$
remains quite small as dark energy starts to dominate, so that even for
$|w|\sim 1$ the adiabatic gEdS approximation for the linear growth factor can
remain good.

The relative smallness of $S$ arises from an approximate cancellation of the
two terms in Eq. (97), which corresponds to that between the two terms in
square brackets in Eq. (92). Indeed at high redshift we have
$S\approx\frac{3}{2}-\frac{9}{5}=-0.3$.

A corollary of this observed accuracy of the leading gEdS interpolation is
that the approximation for $\alpha_{1}$ obtained by solving the linear Eq.
(96) is an accurate one. Fig. 4 shows a comparison of the value of the growth
rate obtained in this way with the exact $\alpha_{1}$. We see that it is
indeed accurate well below the sub-percent level even for a LCDM model at
$z=0$.

Figure 4: Ratio of the growth rate $\alpha_{11}$ obtained using the linearized
approximation Eq. (96) for $\epsilon_{1}=(\alpha_{1}/\alpha_{10})-1$, compared
to the exact solution from Eq. (92).

###  IV.2 Second-order growth rates interpolated on gEdS models

In treating LCDM models in perturbation theory at second and higher order, a
considerable additional complexity arises because the separability of the EdS
solution is not valid. In practice it has been shown (see e.g. [14, 28, 20]
and references therein) that the associated corrections are very small for
standard type cosmological models, and for most purposes can be done using the
EdS kernels and the assumption of separability, replacing simply the EdS
linear theory growth rate by that of the model. Given the ever-improving
precision of cosmological observations, however, even small theoretical errors
arising from cosmology dependence are of interest and have been discussed in a
number of recent works (see e.g. [46, 47, 48, 49, 20, 26, 30]).

We discuss here the calculations of such corrections in light of our results
in gEdS models. More specifically we consider whether the time-dependent
functions, additional to $D_{1}(a)$, which appear in describing the cosmology
dependence of the evolved density field can be interpolated on gEdS models in
the sense we have defined at the beginning of this section. In this subsection
we consider the density fluctuation at second order, and in the following one
at third order. For conciseness, we will not present here the full derivation
of the results we use from EPT including cosmological corrections as it has
been given in several other works, both for LCDM-type models ([27, 20, 32]) or
in a broader class of models (e.g. [28, 17]). As anticipated above we will
follow in the main text the treatment of [27], but we also provide in
Appendices E-F detailed comparison with several other works [32, 11, 20].

The solution for the density perturbation in a general FLRW cosmology at
second order is given by [27] as

$\delta_{2}(\textbf{k},a)=D_{2A}A(\textbf{k})+D_{2B}B(\textbf{k}),$ (98)

where

$\displaystyle A(\textbf{k})=\frac{5}{7}\hat{A}(\textbf{k})=\frac{5}{7}\int
d^{3}\textbf{q}\tilde{\alpha}(\textbf{q},\textbf{k}-\textbf{q})\delta_{1}(\textbf{q})\delta_{1}(\textbf{k}-\textbf{q}),$
$\displaystyle B(\textbf{k})=\frac{2}{7}\hat{B}(\textbf{k})=\frac{2}{7}\int
d^{3}\textbf{q}\tilde{\beta}(\textbf{q},\textbf{k}-\textbf{q})\delta_{1}(\textbf{q})\delta_{1}(\textbf{k}-\textbf{q}),$

and $D_{2A}$ and $D_{2B}$ are time-dependent functions which are solutions to
the equations

$\displaystyle\frac{d^{2}}{d\ln a^{2}}\frac{D_{2}}{a^{2}}+\Big{(}6+\frac{d\ln
H}{d\ln a}\Big{)}\frac{d}{d\ln a}\frac{D_{2}}{a^{2}}$ (100)
$\displaystyle+\Big{[}8+2\frac{d\ln H}{d\ln
a}-\frac{3}{2}\Omega_{m}(a)\Big{]}\frac{D_{2}}{a^{2}}$ $\displaystyle=$
$\displaystyle\begin{cases}\frac{7}{5}\Bigg{[}\Big{(}\frac{dD_{1}}{da}\Big{)}^{2}+\frac{3}{2}\Omega_{m}(a)\Big{(}\frac{D_{1}}{a}\Big{)}^{2}\Bigg{]}\text{
for $D_{2A}$},\\\ \frac{7}{2}\Big{(}\frac{dD_{1}}{da}\Big{)}^{2}\text{ for
$D_{2B}$}.\end{cases}$

These equations have been written in a form adapted to their comparison with
the standard EdS model, for which the right-hand side vanishes giving the
solutions $D_{2A}=D_{2B}=a^{2}$. It is straightforward to verify that in the
gEdS cosmology they admit the solutions

$D_{2A}=\frac{7}{5}\frac{1+4\alpha}{1+6\alpha}a^{2\alpha},\quad
D_{2B}=\frac{7}{2}\frac{2\alpha}{1+6\alpha}a^{2\alpha},$ (101)

corresponding exactly to the separable solutions we have derived in the
previous section.

In light of these solutions it is natural to define the functions

$d_{2A}=\frac{5}{7}\frac{D_{2A}}{D_{1}^{2}},\quad
d_{2B}=\frac{2}{7}\frac{D_{2B}}{D_{1}^{2}},$ (102)

in terms of which

$\delta_{2}(\textbf{k},a)=d_{2A}D_{1}^{2}\hat{A}(\textbf{k})+d_{2B}D_{1}^{2}\hat{B}(\textbf{k})\,.$
(103)

When the functions $d_{2A}$ and $d_{2B}$ are constant in time,
$\delta_{2}(\textbf{k},a)$ is a time-independent functional of
$\delta_{1}(\textbf{k},a)$. Non-separability, on the contrary, corresponds to
a time-dependence of the functions $d_{2A}$ and $d_{2B}$. We first note that
Eqs. (100) can be written, making use of the definitions of $\alpha_{1}$ and
$w$ in the previous section, and using Eq. (92) to eliminate terms in
$d\ln\alpha_{1}/d\ln a$, as

${\cal{D}}^{(2)}d_{2A}=\alpha_{1}^{2}+\frac{3}{2}\Omega_{m},\quad{\cal{D}}^{(2)}d_{2B}=\alpha_{1}^{2},$
(104)

where

${\cal{D}}^{(2)}=\frac{d^{2}}{d\ln
a^{2}}+\Big{[}\frac{1}{2}(1-3w)+4\alpha_{1}\Big{]}\frac{d}{d\ln
a}+(2\alpha_{1}^{2}+\frac{3}{2}\Omega_{m}).$ (105)

We see immediately that ${\cal{D}}^{(2)}(d_{2A}+d_{2B}-1)=0$, from which it
follows that the evident property

$d_{2A}+d_{2B}=1,$ (106)

of the usual EdS solution (with $d_{2A}=5/7$ and $d_{2B}=2/7$), shared by the
gEdS cosmologies (with $d_{2A}=\frac{1+4\alpha}{1+6\alpha}$ and
$d_{2B}=\frac{2\alpha}{1+6\alpha}$), will hold in a general LCDM-like
cosmology, if at asymptotically early times the cosmology approaches EdS (or
gEdS). There is therefore in this case only one real function (in addition to
$D_{1}$) is needed to describe the evolution of the density perturbation at
second order, just as in the original treatment of [11]. Given our gEdS
solutions (101), it is natural to define, without loss of generality, the time
dependence using the function $\alpha_{2}(a)$ defined by

$d_{2A}(a)=\frac{1+4\alpha_{2}(a)}{1+6\alpha_{2}(a)},\quad
d_{2B}(a)=1-d_{2A}(a)\,.$ (107)

$\alpha_{2}(a)$ is just then the value of $\alpha$ in the gEdS model with the
same instantaneous values of $d_{2A}$ and $d_{2B}$, i.e., it is the effective
growth exponent required to interpolate (exactly) these functions on gEdS
models. Anticipating the fact that we will consider in particular here the
small corrections due to the non-separability of the LCDM model relative to
the EdS approximation (with $d_{2A}=5/7$ and $d_{2B}=2/7$), we also define

$d_{2A}=\frac{5}{7}\Big{[}1+\frac{2}{35}\gamma_{2}\Big{]},\quad
d_{2B}=\frac{2}{7}\Big{[}1-\frac{1}{7}\gamma_{2}\Big{]},$ (108)

where

$\gamma_{2}=\frac{1-\alpha_{\rm 2}}{1+\frac{6}{7}(1-\alpha_{\rm 2})}\,,$ (109)

i.e. the coefficients in the definition of the perturbation $\gamma_{2}$ have
been chosen so that, in the limit of small $\gamma_{2}$, we obtain
$\gamma_{2}=1-\alpha_{2}$.

From Eq. (104) we have that the equation obeyed by $\gamma_{2}$ is

${\cal{D}}^{(2)}\gamma_{2}=\frac{21}{2}\Big{(}\Omega_{m}-\alpha_{1}^{2}\Big{)}.$
(110)

It is instructive to rewrite this equation as

$\displaystyle\frac{d^{2}\gamma_{2}}{d\ln a^{2}}$ $\displaystyle+$
$\displaystyle\Big{[}\frac{1}{2}(1-3w)+4\alpha_{1}\Big{]}\frac{d\gamma_{2}}{d\ln
a}$ (111) $\displaystyle+$
$\displaystyle\Big{(}2\alpha_{1}^{2}+\frac{3}{2}\Omega_{m}\Big{)}[\gamma_{2}-\gamma_{2}^{(0)}]=0,$

where

$\gamma_{2}^{(0)}=\frac{21(\Omega_{m}-\alpha_{1}^{2})}{4\alpha_{1}^{2}+3\Omega_{m}}\,.$
(112)

We recover explicitly from this equation the result [38, 39] that the solution
is separable if and only if the ratio $\Omega_{m}/\alpha_{1}^{2}$ is constant.
Indeed, as we have noted, separability corresponds to a time-independent
solution for $\gamma_{2}$, which here is admitted if and only if
$\gamma_{2}^{(0)}$ is constant in time, and therefore if and only if the ratio
$\Omega_{m}/\alpha_{1}^{2}$ is so. More generally we note that solving in the
adiabatic gEdS approximation, i.e. neglecting all terms proportional to $w$
and all time derivatives, and taking, in the same approximation,
$\alpha_{1}=\alpha_{10}$, we then obtain at leading order

$\gamma_{2}=7\frac{1-\alpha_{10}}{1+6\alpha_{10}},$ (113)

which corresponds to

$\alpha_{\rm 2}=\alpha_{10},$ (114)

i.e. we recover as expected, as $w\rightarrow 0$, the adiabatic gEdS
approximation for $d_{2A}$ (and $d_{2B}$).

Figure 5: The effective growth rate $\alpha_{2}$ of a gEdS model with the same
instantaneous value of the functions $d_{2A}$ (and $d_{2B}$), for dark energy
models with the constant equation of state parametrized by the indicated
values of $w_{0}$. We see that, for very small absolute values of $w_{0}$,
$\alpha_{2}$ is well approximated by the adiabatic linear growth rate
$\alpha_{10}$ so that, as expected, $d_{2A}$ and $d_{2B}$ are well
approximated by an adiabatic interpolation of gEdS models. For larger absolute
values of $w_{0}$, $\alpha_{2}$ stays much closer to its initial value because
of rapid evolution which causes effective damping in its evolution equation.

Fig. 5 shows the numerical solutions obtained for $\alpha_{2}$ as a function
of $\Omega_{m}$, for models with the different indicated values of $w_{0}$.
These are obtained by solving Eq. (110) for $\gamma_{2}$, with initial
conditions corresponding to the asymptote to EdS as $z\rightarrow\infty$:

$\gamma_{2}(0)=0,\quad\frac{d\gamma_{2}}{d\ln a}(0)=0\,.$ (115)

In practice, rather than extrapolating to an initial time $a_{i}\ll 1$ at
which $\Omega_{m}\approx 1$, we can avoid the very early time transients by
taking directly as the initial condition

$\gamma_{2}(a_{i})=1-\alpha_{10}(a_{i}),\quad\frac{d\gamma_{2}}{d\ln
a}(a_{i})=-\frac{d\alpha_{10}}{d\ln a}(a_{i})\,,$ (116)

since $\gamma_{2}(a)=1-\alpha_{10}(a)$, as we have noted above, becomes a good
approximation to the exact evolution in the limit $w\ll 1$. Indeed this can be
seen in Fig. 5, in which this adiabatic approximation $\alpha_{2}=\alpha_{10}$
is shown (solid line). We see that, for very small absolute values of $w_{0}$,
$\alpha_{2}$ is very well approximated by $\alpha_{10}$. In other words as
expected, $d_{2A}$ and $d_{2B}$ are then well approximated by an adiabatic
interpolation of gEdS models. However, unlike what we observed for
$\alpha_{1}(z)$, this adiabatic approximation fails for increasing values of
$|w_{0}|$. In this case we see that $\alpha_{2}$ stays much closer to its
initial value. This can be understood simply as a result of the increasing
contribution of the dark energy, which increases (through the dependence in
$w$) an effective damping term in Eq. (110) which slows greatly the
“relaxation” to $\gamma_{2}=\gamma_{2}^{(0)}$ occurring for small $|w_{0}|$.

###  IV.3 Third order growth rates

Following again the treatment of [27], we write the density perturbation to
third-order as

$\displaystyle\delta_{3}(\textbf{k},a)$ $\displaystyle=$ $\displaystyle
D_{3AA}(a)C_{AA}(\textbf{k})+D^{\prime}_{3AA}(a)C^{\prime}_{AA}(\textbf{k})$
$\displaystyle+D_{3AB}(a)C_{AB}(\textbf{k})+D^{\prime}_{3AB}(a)C^{\prime}_{AB}(\textbf{k})$
$\displaystyle+D_{3BA}(a)C_{BA}(\textbf{k})+D_{3BB}(a)C_{BB}(\textbf{k}),$

where

$\displaystyle C_{AA}(\textbf{k})$ $\displaystyle=$
$\displaystyle\frac{7}{18}\int
d^{3}\textbf{q}\tilde{\alpha}(\textbf{q},\textbf{k}-\textbf{q})\delta_{1}(\textbf{q})A(\textbf{k}-\textbf{q}),$
$\displaystyle C_{AA}^{\prime}(\textbf{k})$ $\displaystyle=$
$\displaystyle\frac{7}{30}\int
d^{3}\textbf{q}\tilde{\alpha}(\textbf{q},\textbf{k}-\textbf{q})\delta_{1}(\textbf{k}-\textbf{q})A(\textbf{q}),$
$\displaystyle C_{AB}(\textbf{k})$ $\displaystyle=$
$\displaystyle\frac{7}{18}\int
d^{3}\textbf{q}\tilde{\alpha}(\textbf{q},\textbf{k}-\textbf{q})\delta_{1}(\textbf{q})B(\textbf{k}-\textbf{q}),$
$\displaystyle C_{AB}^{\prime}(\textbf{k})$ $\displaystyle=$
$\displaystyle\frac{7}{9}\int
d^{3}\textbf{q}\tilde{\alpha}(\textbf{q},\textbf{k}-\textbf{q})\delta_{1}(\textbf{k}-\textbf{q})B(\textbf{q}),$
$\displaystyle C_{BA}(\textbf{k})$ $\displaystyle=$
$\displaystyle\frac{2}{15}\int
d^{3}\textbf{q}\tilde{\beta}(\textbf{q},\textbf{k}-\textbf{q})\delta_{1}(\textbf{q})A(\textbf{k}-\textbf{q}),$
$\displaystyle C_{BB}(\textbf{k})$ $\displaystyle=$
$\displaystyle\frac{4}{9}\int
d^{3}\textbf{q}\tilde{\beta}(\textbf{q},\textbf{k}-\textbf{q})\delta_{1}(\textbf{q})B(\textbf{k}-\textbf{q}).$

and the six functions $D_{3XX}$ are the solutions of the four differential
equations

$\displaystyle\frac{d^{2}}{d\ln a^{2}}\frac{D_{3}}{a^{3}}+\Big{(}8+\frac{d\ln
H}{d\ln a}\Big{)}\frac{d}{d\ln a}\frac{D_{3}}{a^{3}}$ (124)
$\displaystyle+\Big{[}15+3\frac{d\ln H}{d\ln
a}-\frac{3}{2}\Omega_{m}(a)\Big{]}\frac{D_{3}}{a^{3}}$ $\displaystyle=$
$\displaystyle\begin{cases}\frac{18}{7}\Big{[}2\frac{dD_{1}}{da}+\frac{3}{2}\Omega_{m}(a)\frac{D_{1}}{a}\Big{]}\frac{D_{2A,B}}{a^{2}}\\\
+\frac{18}{7}a\frac{dD_{1}}{da}\frac{d}{da}\frac{D_{2A,B}}{a^{2}}\text{ for
$D_{3AA,3AB}$}\\\
15\frac{dD_{1}}{da}\Big{[}a\frac{d}{da}\frac{D_{2A}}{a^{2}}+2\frac{D_{2A}}{a^{2}}-\frac{7}{5}\frac{D_{1}}{a}\frac{dD_{1}}{da}\Big{]}\text{
for $D_{3BA}$}\\\
\frac{9}{2}\frac{dD_{1}}{da}\Big{[}a\frac{d}{da}\frac{D_{2B}}{a^{2}}+2\frac{D_{2B}}{a^{2}}\Big{]}\text{
for $D_{3BB}$},\end{cases}$

and the two remaining functions are determined by the relations

$\displaystyle\frac{5}{18}D_{3AA}+\frac{2}{9}D_{3AB}^{\prime}=\frac{1}{2}D_{1}^{3},$
$\displaystyle\frac{1}{6}D_{3AA}^{\prime}+\frac{1}{9}D_{3AB}+\frac{2}{21}D_{3BA}+\frac{8}{63}D_{3BB}=\frac{1}{2}D_{1}^{3}.$

Analogously to the previous section we now define

$\displaystyle d_{3AA}$ $\displaystyle=$
$\displaystyle\frac{5}{9}\frac{D_{3AA}}{D_{1}^{3}},\quad
d_{3AA}^{\prime}=\frac{1}{3}\frac{D_{3AA}^{\prime}}{D_{1}^{3}},$
$\displaystyle d_{3AB}$ $\displaystyle=$
$\displaystyle\frac{2}{9}\frac{D_{3AB}}{D_{1}^{3}},\quad
d_{3AB}^{\prime}=\frac{4}{9}\frac{D_{3AB}^{\prime}}{D_{1}^{3}},$
$\displaystyle\quad d_{3BA}$ $\displaystyle=$
$\displaystyle\frac{2}{21}\frac{D_{3BA}}{D_{1}^{3}},\quad
d_{3BB}=\frac{8}{63}\frac{D_{3BB}}{D_{1}^{3}}.$ (126)

For the case of the gEdS cosmologies, it is simple to verify using the kernels
$F_{2}$, $G_{2}$ and $F_{3}$ derived in Section II that we obtain

$\displaystyle d_{3AA}$ $\displaystyle=$ $\displaystyle
d_{3}d_{2}=\frac{1+4\alpha}{1+8\alpha},$ $\displaystyle d^{\prime}_{3AA}$
$\displaystyle=$ $\displaystyle
d_{3}\tilde{d}_{2}=\frac{1+2\alpha}{1+8\alpha},$ $\displaystyle d_{3AB}$
$\displaystyle=$ $\displaystyle d_{3}(1-d_{2})=\frac{2\alpha}{1+8\alpha},$
$\displaystyle d^{\prime}_{3AB}$ $\displaystyle=$ $\displaystyle
d_{3}(1-\tilde{d}_{2})=\frac{4\alpha}{1+8\alpha},$ $\displaystyle d_{3BA}$
$\displaystyle=$
$\displaystyle(1-d_{3})\tilde{d}_{2}=\frac{2\alpha(1+2\alpha)}{(1+6\alpha)(1+8\alpha)},$
$\displaystyle d_{3BB}$ $\displaystyle=$
$\displaystyle(1-d_{3})(1-\tilde{d}_{2})=\frac{8\alpha^{2}}{(1+6\alpha)(1+8\alpha)}.$

It is straightforward to show that using Eqs. (124), and $d_{2B}=1-d_{2A}$, we
obtain

$\displaystyle{\cal{D}}^{(3)}d_{3AA}$ $\displaystyle=$
$\displaystyle(4\alpha_{1}^{2}+3\Omega_{m})d_{2A}+2\alpha_{1}\frac{d\,d_{2A}}{d\ln
a},$ $\displaystyle{\cal{D}}^{(3)}d_{3AB}$ $\displaystyle=$
$\displaystyle(4\alpha_{1}^{2}+3\Omega_{m})-(4\alpha_{1}^{2}+3\Omega_{m})d_{2A}$
$\displaystyle-2\alpha_{1}\frac{d\,d_{2A}}{d\ln a},$
$\displaystyle{\cal{D}}^{(3)}d_{3BA}$ $\displaystyle=$
$\displaystyle-2\alpha_{1}^{2}+4\alpha_{1}^{2}d_{2A}+2\alpha_{1}\frac{d\,d_{2A}}{d\ln
a},$ $\displaystyle{\cal{D}}^{(3)}d_{3BB}$ $\displaystyle=$ $\displaystyle
4\alpha_{1}^{2}-4\alpha_{1}^{2}d_{2A}-2\alpha_{1}\frac{d\,d_{2A}}{d\ln a},$
(128)

where

${\cal{D}}^{(3)}=\frac{d^{2}}{d\ln
a^{2}}+\Big{[}\frac{1}{2}(1-3w)+6\alpha_{1}\Big{]}\frac{d}{d\ln
a}+3(2\alpha_{1}^{2}+\Omega_{m})\,.$ (129)

The two constraints Eqs. (IV.3) may be written as

$\displaystyle d_{3AA}+d_{3AB}^{\prime}=1,$ (130) $\displaystyle
d_{3AB}+d_{3AA}^{\prime}+2(d_{3BA}+d_{3BB})=1.$ (131)

Setting all derivatives to zero in Eq. (IV.3), we recover, as required, the
expression for the gEdS cosmologies Eq. (IV.3) above.

Taking the sum of Eqs. (IV.3) we note that we obtain also

${\cal{D}}^{(3)}(d_{3AA}+d_{3AB}+d_{3BA}+d_{3BB}-1)=0,$ (132)

and therefore, given that this quantity is zero in an EdS (and gEdS)
cosmology, we obtain the additional constraint

$\displaystyle d_{3AA}+d_{3AB}+d_{3BA}+d_{3BB}=1,$ (133)

when we assume that the cosmology is asymptotically EdS (or gEdS). Further
writing Eq. (104) for $d_{2A}$ as

${\cal{D}}^{(3)}d_{2A}=\alpha_{1}^{2}+\frac{3}{2}\Omega_{m}+(4\alpha_{1}^{2}+\frac{3}{2}\Omega_{m})d_{2A}+2\alpha_{1}\frac{d\,d_{2A}}{d\ln
a},$ (134)

we see also that

${\cal{D}}^{(3)}(2d_{2A}-d_{3AA}+d_{3BB}-1)=0,$ (135)

so that we can infer (again assuming the cosmology to be asymptotically EdS or
gEdS) the additional relation

$2d_{2A}-d_{3AA}+d_{3BB}=1.$ (136)

The full set of constraints now reads

$\displaystyle d_{3AA}+d_{3AB}^{\prime}=1,$ (137) $\displaystyle
d_{3AB}^{\prime}+d_{3AA}^{\prime}+d_{3BA}+d_{3BB}=1,$ (138) $\displaystyle
d_{3AB}+d_{3AA}+d_{3BA}+d_{3BB}=1,$ (139) $\displaystyle
2d_{2A}-d_{3AA}+d_{3BB}=1.$ (140)

We can thus conclude that the perturbation at the third order can therefore be
described fully by only just two independent functions in addition to
$d_{2A}$. This result agrees again, as we saw at second order, with the
original analysis of [11] for the asymptotic EdS case. The exact mapping
between our functions and those of [11] is given in Appendix F. Our analysis
above shows that this same result (and either our parametrization or that of
[11]) generalizes to the broader class of cosmologies which are asymptotically
gEdS.

We choose now (arbitrarily) to use, in addition to $d_{2A}$, the two functions
$d_{3AA}$ and $d_{3AB}$ to describe the cosmology-dependent corrections to the
third-order perturbation To describe them in terms of interpolation on gEdS
models we thus introduce, following the treatment of the second order
perturbation, effective growth exponents as follows:

$\displaystyle d_{3AA}(a)$ $\displaystyle=$
$\displaystyle\frac{1+4\alpha_{3AA}(a)}{1+8\alpha_{3AA}(a)}=\frac{5}{9}\Big{[}1+\frac{4}{45}\gamma_{3AA}(a)\Big{]},$
$\displaystyle d_{3AB}(a)$ $\displaystyle=$
$\displaystyle\frac{2\alpha_{3AB}(a)}{1+8\alpha_{3AB}(a)}=\frac{2}{9}\Big{[}1-\frac{1}{9}\gamma_{3AB}(a)\Big{]},$

where, as for the second order case, the numerical coefficients have been
chosen so that we will have, in the limit of small $\gamma_{3XX}$, that
$\alpha_{3XX}\approx 1-\gamma_{3XX}$.

The evolution of the $\gamma_{3AA}$ and $\gamma_{3AB}$ are now given by the
equations

$\displaystyle{\cal{D}}^{(3)}\gamma_{3AA}$ $\displaystyle=$
$\displaystyle\frac{135}{14}\Big{(}\Omega_{m}-\alpha_{1}^{2}\Big{)}$ (142)
$\displaystyle+\frac{81}{98}\Big{[}(4\alpha_{1}^{2}+3\Omega_{m})\gamma_{2}+2\alpha_{1}\frac{d\,\gamma_{2}}{d\ln
a}\Big{]},\ $ $\displaystyle{\cal{D}}^{(3)}\gamma_{3AB}$ $\displaystyle=$
$\displaystyle-\frac{54}{7}\Big{(}\Omega_{m}-\alpha_{1}^{2}\Big{)}$
$\displaystyle+\frac{81}{49}\Big{[}(4\alpha_{1}^{2}+3\Omega_{m})\gamma_{2}+2\alpha_{1}\frac{d\,\gamma_{2}}{d\ln
a}\Big{]},\ $

which can also be written as

$\begin{split}\frac{d^{2}\gamma_{3XX}}{d\ln
a^{2}}+&\Big{[}\frac{1}{2}(1-3w)+6\alpha_{1}\Big{]}\frac{d\gamma_{3XX}}{d\ln
a}+3\Big{[}2\alpha_{1}^{2}+\Omega_{m}\Big{]}\\\
&\times(\gamma_{3XX}-\gamma_{3XX}^{(0)})=c_{3XX}\Big{[}(4\alpha_{1}^{2}+3\Omega_{m})\\\
&\times(\gamma_{2}-\gamma_{2}^{(0)})+2\alpha_{1}\frac{d\,\gamma_{2}}{d\ln
a}\Big{]},\end{split}$ (144)

where $c_{3AA}=\frac{81}{98},c_{3AB}=\frac{81}{49}$ and

$\gamma_{3AA}^{(0)}=\gamma_{3AB}^{(0)}=\frac{9(\Omega_{m}-\alpha_{1}^{2})}{2\alpha_{1}^{2}+\Omega_{m}}.$
(145)

For purposes of comparison it is useful to rewrite Eq. (111) for $\gamma_{2}$
as

$\begin{split}{\cal{D}}^{(3)}\gamma_{2}&=\frac{21}{2}\Big{(}\Omega_{m}-\alpha_{1}^{2}\Big{)}+\Big{[}(4\alpha_{1}^{2}+\frac{3}{2}\Omega_{m})\gamma_{2}+2\alpha_{1}\frac{d\,\gamma_{2}}{d\ln
a}\Big{]},\end{split}$ (146)

and

$\begin{split}\frac{d^{2}\gamma_{2}}{d\ln
a^{2}}+&\Big{[}\frac{1}{2}(1-3w)+6\alpha_{1}\Big{]}\frac{d\gamma_{2}}{d\ln
a}+3\Big{[}2\alpha_{1}^{2}+\Omega_{m}\Big{]}\\\
&\times(\gamma_{2}-\gamma_{2}^{(0)})=\Big{[}(4\alpha_{1}^{2}+\frac{3}{2}\Omega_{m})(\gamma_{2}-\gamma_{2}^{(0)})\\\
&+2\alpha_{1}\frac{d\,\gamma_{2}}{d\ln a}\Big{]}.\end{split}$ (147)

Solving Eqs. (146) and (142) in the adiabatic gEdS approximation, i.e.
neglecting all time derivatives in these equations, we can see that we obtain
the solution Eq. (113) for $\gamma_{2}$ and now also

$\gamma_{3AA}=\gamma_{3AB}=9\frac{1-\alpha_{10}}{1+8\alpha_{1,0}},$ (148)

which corresponds, as anticipated, to the leading adiabatic solutions

$\alpha_{2}(a)=\alpha_{3AA}(a)=\alpha_{3AB}(a)=\alpha_{10}(a).$ (149)

Figure 6: Evolution of the effective growth rate $\alpha_{2}$, $\alpha_{3AA}$
and $\alpha_{3AB}$ which describes in general the cosmology-dependent
corrections to the separable EdS approximation at third order in perturbation
theory. The different sets of curves correspond to the different indicated
values of $w_{0}$. Also shown is the adiabatic linear growth rate
$\alpha_{10}$. We see that for the smallest values the adiabatic approximation
is good, but that it degrades increasingly as $w$ increases. Figure 7: Same as
in the previous figure but for models with larger absolute values of $w_{0}$.
We see that for these values of $w_{0}$ becomes larger, one of the three
exponents has quite a different value from the $\alpha_{2}$, so an
interpolation of the functions describing deviations from separability on a
single gEdS the model with $\alpha=\alpha_{2}$ is not necessarily a valid
approximation in general, although it turns out to be for the calculation of
the one-loop PS.

Fig. 6 and Fig. 7 show the results for the three effective growth exponents
$\alpha_{2}$ (as in Fig. 5), $\alpha_{3AA}$ and $\alpha_{3AB}$ as a function
of $\Omega_{m}$ obtained by numerically integrating these equations for
different dark energy models with the fixed equation of state with the same
initial conditions as for $\gamma_{2}$ (corresponding to EdS at high
redshift). For the smallest value of $w_{0}$ we see, as expected, all three
exponents well approximated by the adiabatic solution Eq. (149). Then as
$|w_{0}|$ increases we see that not only do the exponents differ from this
solution but they start differing, albeit to a lesser extent, from one
another. More specifically we note that $\alpha_{2}\approx\alpha_{3AA}$, while
the third exponent differs. The origin of these behaviours can be seen easily
in the coefficients of the equations above: those in the equations for
$\gamma_{2}$ and $\alpha_{3AA}$ are almost identical, while those in the
equation $\alpha_{3AB}$ differ not only in magnitude (by about a factor of
$2$) but also in sign. In Appendix D we report a detailed comparison of our
numerical results with those of [27], and adapt the fitting functions provided
by it to infer corresponding fits for the two functions $\gamma_{3AA}$ and
$\gamma_{3AB}$.

###  IV.4 The PS at one loop

We now derive our simplified expressions for the PS at one loop in the class
of FLRW cosmologies specified in Section II.1 (clustering matter plus a
generic dark energy fluid component), making use of the reduction of the
number of independent functions and expressing these in terms of the effective
growth exponents we have defined above. We also consider finally to what
accuracy this exact result can be approximated by a direct interpolation of
the PS in gEdS models (i.e. by just replacing $\alpha$ in the analytic gEdS
kernels by a single effective growth rate). In our analysis we start from the
expressions we need (for each of the ‘22’ and ‘13’ contributions) in a form
easily comparable (up to notation) to the expressions given in [27]. This
facilitates the explanation of the simplifications we obtain.

####  IV.4.1 ‘22’ term

Starting from Eq. (98), we obtain the full time-dependent $P_{22}$ directly as

$\displaystyle P_{22}$ $\displaystyle=$ $\displaystyle
d_{2A}^{2}[M_{0}+M_{1}+M_{2}]+2d_{2A}d_{2B}[M_{0}+\frac{1}{2}M_{1}]$ (150)
$\displaystyle+d_{2BB}^{2}M_{0}$ $\displaystyle=$
$\displaystyle(d_{2A}+d_{2B})^{2}M_{0}+d_{2A}(d_{2A}+d_{2B})M_{1}$
$\displaystyle+d_{2A}^{2}M_{2}$

where the integrals $M_{i}$ are those defined above in Eq. (77). The first
expression is, up to simple notational differences, identical to that given by
[27].222An equivalent expression is given in [28], also in terms of two
redshift dependent functions and three integrals which are linear combinations
of $M_{0}$, $M_{1}$ and $M_{2}$. Making use of the additional relation Eq.
(106) we have derived above, this simplifies to

$P_{22}=M_{0}+d_{2A}M_{1}+d_{2A}^{2}M_{2},$ (151)

i.e. to a form evidently equivalent to that in the gEdS model Eq. (75) when we
parametrize $d_{2A}$ by $\alpha_{2}$ or $\gamma_{2}$ as defined in the
previous section. Defining now

$\Delta P_{22}=P_{22}-P_{22}^{EdS},$ (152)

where $P_{22}^{EdS}$ is the result in the separable EdS approximation

$P_{22}^{EdS}=M_{0}+\frac{5}{7}M_{1}+\frac{25}{49}M_{2},$ (153)

the cosmology dependent correction $\Delta P_{22}$ is expressed explicitly as
a function of $\gamma_{2}(a)$ by

$\Delta
P_{22}=\frac{2}{49}\gamma_{2}(M_{1}+\frac{10}{7}M_{2})+(\frac{2}{49}\gamma_{2})^{2}M_{2}.$
(154)

Recalling our discussion below we see that this correction is thus given in
terms just of integrals that are explicitly infra-red safe, which is not the
case if we do not make use of the relation Eq. (106).

Table 2: A summary of the calculation of the three “effective growth rates” parametrizing the second and third-order density kernels using the different approximations we have discussed and also the exact calculation. Approximation | $\alpha_{2}$ | $\alpha_{3AA}$ | $\alpha_{3AB}$ | Precision
---|---|---|---|---
EdS | $1$ | $1$ | $1$ |  $\lesssim 0.5\%$ for PS at $z=0$ (LCDM)
adiabatic gEdS |  $\alpha_{1}=\frac{d\ln D_{1}}{d\ln a}$, Eq. (92) | $\alpha_{1}$ | $\alpha_{1}$ | valid for small $|w_{0}|$ (cf. Figs. 5, 6, and 7)
gEdS |  $\alpha_{2}$, Eq. (111) | $\alpha_{2}$ | $\alpha_{2}$ |  $\lesssim 0.1\%$ for PS at $z=0$ (LCDM)
exact |  $\alpha_{2}$, Eq. (111) |  $\alpha_{3AA}$ , Eq.(144) |  $\alpha_{3AB}$, Eq. (144) | —

####  IV.4.2 ‘13’ term

Starting from Eq. (IV.3), the expression for $P_{13}$ can be obtained directly
as

$\displaystyle 2P_{13}$ $\displaystyle=$ $\displaystyle
d_{3AA}(N_{0}+N_{3})+d_{3AA}^{\prime}(N_{1}-N_{3})$ (155)
$\displaystyle+d_{3AB}(N_{0}+N_{3})-d_{3AB}^{\prime}N_{3}+d_{3BA}(N_{0}+N_{2})$
$\displaystyle+d_{3BB}N_{0}$ $\displaystyle=$
$\displaystyle(d_{3AA}+d_{3AB}+d_{3BA}+d_{3BB})N_{0}$
$\displaystyle+d_{3AA}^{\prime}N_{1}+d_{3BA}N_{2}+(d_{3AA}-d_{3AA}^{\prime}$
$\displaystyle+d_{3AB}-d_{3AB}^{\prime})N_{3},$

where the integrals $N_{0},N_{1},N_{2}$ are those defined by Eq. (81) and Eq.
(82). $N_{3}$ is an additional integral defined also by Eq. (81) but with
$n_{3}=-\frac{4}{3}r^{2}$, and which, it is simple to verify, is infra-red
safe. The first expression has been written to be easily compared to that
given by [27] in terms of the six functions $d_{3XX}$ and six integrals.333The
equivalent expression in [28] is given in terms of six redshift dependent
functions and five integrals. Making use now of the relations Eqs. (137)-(140)
obtained from the EdS (or gEdS) boundary condition, we find that the
coefficient of $N_{3}$ vanishes, and further that the expression can be
simplified to be a function of only of $d_{2A}$ and only one of the initial
six functions $d_{3XX}$:

$2P_{13}=N_{0}-N_{1}+2d_{2A}N_{1}+d_{3BA}(N_{2}-N_{1}),$ (156)

Thus, as for $P_{22}$, we obtain an expression involving only the same three
integrals as in the gEdS one-loop PS, and as in that case, the cosmology-
dependent terms involve only the two infra-red safe integrals $N_{1}$ and
$N_{2}$. We note that this latter property is obtained using the relation Eq.
(139). Conversely, without employing this relation, which we have shown here
follows from the condition that the cosmology asymptotes to EdS (or gEdS) at
early times, it is not possible to recover explicitly the infra-red safety of
the cosmological correction.

The cosmological correction in $P_{13}$ can thus be taken to depend in general
only on $\gamma_{2}$, as defined in Eq. (108), and on just one additional
function. It is convenient then to introduce the function
$\gamma_{3BA}=\gamma_{3}$, defined, following the same treatment as in Section
IV, as

$d_{3BA}=\frac{2\alpha_{3}(1+2\alpha_{3})}{(1+6\alpha_{3})(1+8\alpha_{3})}=\frac{2}{21}\Big{[}1+\frac{5}{63}\gamma_{3}\Big{]},$
(157)

so that $\alpha_{3}\approx 1-\gamma_{3}$ for $\gamma_{3}\ll 1$. $\gamma_{3}$
is given in terms of $\gamma_{2}$ and the two functions $\gamma_{3AA}$ and
$\gamma_{3BA}$ analysed in the Section IV, by

$\displaystyle 15\gamma_{3}=162\gamma_{2}-196\gamma_{3AA}+49\gamma_{3AB}.$
(158)

It is thus a solution to the equation

${\cal{D}}^{(3)}\gamma_{3}=-\frac{189}{5}\Big{(}\Omega_{m}-\alpha_{1}^{2}\Big{)}+\frac{189}{35}\Big{[}4\alpha_{1}^{2}\gamma_{2}+2\alpha_{1}\frac{d\,\gamma_{2}}{d\ln
a}\Big{]},$ (159)

with appropriate boundary conditions. Defining

$2\Delta P_{13}=2P_{13}-2P_{13}^{EdS},$ (160)

where

$2P_{13}^{EdS}=N_{0}+\frac{1}{3}N_{1}+\frac{2}{21}N_{2},$ (161)

is the usual result in the separable EdS approximation, the cosmological
correction is given as a function of $\gamma_{2}$ and $\gamma_{3}$ as

$2\Delta
P_{13}=\frac{4}{49}\big{[}\gamma_{2}N_{1}+\frac{5}{54}\gamma_{3}(N_{2}-N_{1})\big{]}.$
(162)

####  IV.4.3 The cosmological correction to the PS at one loop

Combining Eqs. (154) and (162), we get the final expression for the full
cosmological (i.e. non-EdS) correction to the one-loop PS as

$\displaystyle\Delta P_{1-loop}$ $\displaystyle=$
$\displaystyle\frac{2}{49}\gamma_{2}(M_{1}+\frac{10}{7}M_{2})+(\frac{2}{49}\gamma_{2})^{2}M_{2}$
(163)
$\displaystyle+\frac{4}{49}\big{[}\gamma_{2}N_{1}+\frac{5}{54}\gamma_{3}(N_{2}-N_{1})\big{]}.$

Figure 8: Left panel: The one loop cosmological (non-EdS) correction $\Delta
P_{1-loop}$ as a fraction of $P_{1-loop}$, the full one-loop PS in EdS, at
z=0, for $w_{0}=-1$ and $w_{0}=-0.5$ (and a typical standard like cosmological
PS, see text) calculated using the exact expression Eq. (163), and also using
the exact gEdS kernels with $\alpha$ replaced by the single effective growth
exponent $\alpha_{2}(z)$. Right panel: ratio of the gEdS approximation to the
exact result.

The left panel in Fig. 8 shows $\Delta P_{1-loop}$ relative to the full one
loop PS in EdS, at $z=0$, in a standard cosmological model (specifically, with
the cosmological parameters of [50] and calculated using the approximation to
the transfer function of [51]). The model labeled by $w_{0}=-0.5$ is identical
but for this parameter. As previously documented in the literature (e.g. [27])
the corrections are at the sub percent level in the former case but increase
as $|w_{0}|$ decreases. The reason for this behaviour is very evident in light
of our analysis: a smaller $|w|$ is closer to the adiabatic limit, which
corresponds to an effective $\alpha$ closer to the instantaneous linear growth
exponent $\alpha_{1}=d\log D_{1}/d\log a$ (i.e. to a larger value of
$|\gamma_{2}|$ and $|\gamma_{3}|$.) Also shown are the approximation to
$\Delta P_{1-loop}$ obtained by directly interpolating with a gEdS model i.e.,
by replacing $\alpha$ by $\alpha_{2}(z)$ in Eqs. (75) and (76) (and
subtracting the EdS result from their sum). The right panel of Fig. 8 shows
the ratio of this approximation to the exact result. As might be anticipated
from the expression Eq. (163) — which depends on the function $\gamma_{3}$
only in one term with a small pre-factor — we find that interpolation provides
a good approximation (within about $25\%$). In practice it is not
significantly more difficult to calculate the exact result, but the result
shows that for understanding the origin of the smallness of the calculated
correction, the approximation to gEdS is very instructive. We underline,
however, that the adiabatic approximation $\alpha_{2}\approx\alpha_{10}$ is
not valid, so this is not an adiabatic interpolation in the sense we have
defined. The effective growth exponent $\alpha_{2}$ can be roughly thought of
a time-averaged behaviour of the exact linear growth exponent $\alpha_{1}$,
which is itself a good approximation to $\alpha_{2}$ for small $|w|$.

To summarize we present in Table 2 a synthesis of the different approximations
for the calculation of density perturbations up to third order, in
cosmological models with a self-gravitating matter component and a smooth
time-dependent component (given through an effective equation of state
$p=w(a)\rho$). To parametrize the temporal evolution of the kernels we have
defined the three functions $\alpha_{2}$, $\alpha_{3AA}$ and $\alpha_{3AB}$
(and also equivalent functions $\gamma_{2}$, $\gamma_{3AA}$ and
$\gamma_{3AB}$). The table indicates how each of these three functions is
calculated in each of the three approximations we have discussed for the exact
result.

##  V Discussion and conclusions

We have derived in detail the kernels for both standard EPT and LPT in a
family of EdS models, characterized by their constant linear growth rate
$\alpha$. While such results for this family of models are implicit in
previous treatments of perturbation theory in the FLRW cosmologies (e.g.
[27]), or in a broader class of cosmologies (e.g. [28, 32]), the simple
analytical results in this case which generalize those of the standard EdS
model have not (to our knowledge) been discussed in the literature previously.
For the one-loop power spectrum in EFT we have used these kernels to show that
the expected infra-red convergence properties of the corrections relative to a
standard EdS model are obtained, without the need for cancellation between the
two contributing integrals. Further we have noted that the coefficients of
ultra-violet divergent terms have a strong dependence on $\alpha$, and the
leading divergence can even change sign (and is zero at a certain value).

In the second part of the paper we have shown these analytical results for the
kernels in these models can be exploited to obtain a simplified formulation of
the calculation of the cosmology-dependent corrections to the usual separable
EdS approximation in any standard type non-EdS cosmology, and in particular in
LCDM-like cosmologies. We do so by parametrizing the relevant time-dependent
functions characterizing this cosmology dependence as effective growth rates
in gEdS models. Assuming only that the cosmology at asymptotically early times
is gEdS, second-order corrections are parametrized in terms of one such
exponent, and third-order corrections in terms of it and two more such
exponents. Thus, for example, the bispectrum calculated at leading non-trivial
order in perturbation theory, depending on the second-order perturbation, is,
at any time, exactly equal to that in a gEdS model. For the power spectrum we
show the results on the infra-red convergence behaviour of the gEdS model
generalize, and we derive explicit expressions for the corrections to the EdS
model expressed in terms of infra-red safe integrals. This expression turns
out to depend on only one of the two effective exponents needed at the third
order. Nevertheless the PS is in fact given to a relatively good approximation
(up to $\sim 25\%$) by that obtained with the analytic gEdS kernels on
replacing the fixed growth rate $\alpha$ of gEdS by the single effective
exponent $\alpha_{2}(z)$ defining the exact evolution at second order. Our
exact expression is much simplified compared to previous (equivalent)
expressions derived in [27, 32] which involve six or eight redshift-dependent
functions and do not explicitly recover the infra-red convergent properties as
we have done.

Our formulation of the cosmological corrections in this way makes it simple to
understand why they are so very small in LCDM-like models (typically $<1\%$
even at $z=0$). The effective growth exponents we introduce to map the time-
dependent kernels to those of the gEdS models obey equations in which one can
clearly identify an adiabatic limit, of a dark energy component scaling like
matter (i.e. $w=0$) in which these exponents converge to the instantaneous
logarithmic linear growth rate. The more rapid evolution of dark energy at low
redshift acts like a damping of this evolution, leading to the effective
exponent much closer to its initial EdS value. Combined with the very weak
$\alpha$ dependence shown by the analytical expression for the gEdS kernels,
the corrections turn out to be as small as they are. Only for models
accelerating much faster than $LCDM$ at $z=0$ can such corrections be
comparable to the EdS approximated loop corrections, and even in the
asymptotic limit that $\Omega_{m}\rightarrow 0$ they remain finite and at most
of the order unity.

In a forthcoming paper, we will exploit the results for the kernels in gEdS to
obtain exact analytical solutions for the PS in what we will call “generalized
scale-free models” i.e. for a PS which is a simple power law and the cosmology
is gEdS [34, 33]. These models share the property of self-similarity of usual
scale-free models (with a standard EdS expansion law). This makes them ideal
for numerical testing for the accuracy of the analytical predictions, as the
property of self-similarity allows the extraction of highly accurate results
from simulations (see [35, 36]). In further work we plan also to explore the
ultra-violet divergences in perturbation theory in these models, and to
examine their regularization by methods such as effective field theory. This
may, we anticipate, provide insights into the application of such techniques
in standard cosmologies. It would also be interesting to extend the analysis
we have presented here to determine analogous formulations (in terms of
effective growth rate functions) of analyses previously reported in the
literature for cosmological corrections, both for other statistics (two-point
velocity statistics, pairwise velocities, trispectrum etc.) at one-loop and
two-loop, and also for a broader range of cosmological models. In particular
it would be interesting to explore such a formulation of perturbation theory
with a massive neutrino component including its perturbations (see e.g. [48,
49, 20, 52]). In this case the kernels have also a scale dependence that might
be expected to be conveniently formulated in terms of (scale-dependent)
effective growth exponents analogous to those we have employed here.

###### Acknowledgements.

A.P. is supported by the Indonesia Endowment Fund for Education (LPDP). The
numerical calculations in this research used the python packages Matplotlib
[53], NumPy [54], and SciPy [55].

##  Appendix A PT kernels at third-order in gEdS

###  A.1 EPT

Eq. (18) at third order in $\delta^{(1)}$ gives

$\displaystyle\mathcal{H}^{2}\Big{\\{}-a^{2}\partial_{a}^{2}-\frac{3}{2}a\partial_{a}$
$\displaystyle+$
$\displaystyle\frac{3}{2\kappa^{2}}\Big{\\}}\delta^{(3)}(\textit{{k}},a)=S^{(3)}_{\tilde{\beta}}(\textit{{k}},a)-\mathcal{H}\partial_{a}(aS^{(3)}_{\tilde{\alpha}}(\textit{{k}},a)),$
(164)

where

$\displaystyle
S_{\tilde{\alpha}}^{(3)}(\textit{{k}},a)=-\int\frac{\text{d}^{3}q}{(2\pi)^{3}}\tilde{\alpha}(\textit{{q}},\textit{{k}}-\textit{{q}},a)\theta^{(1)}(\textit{{q}},a)\delta^{(2)}(\textit{{k}}-\textit{{q}},a)-\int\frac{\text{d}^{3}q}{(2\pi)^{3}}\tilde{\alpha}(\textit{{q}},\textit{{k}}-\textit{{q}},a)\theta^{(2)}(\textit{{q}},a)\delta^{(1)}(\textit{{k}}-\textit{{q}},a),$
(165) $\displaystyle
S_{\tilde{\beta}}^{(3)}(\textit{{k}},a)=-\int\frac{\text{d}^{3}q}{(2\pi)^{3}}\tilde{\beta}(\textit{{q}},\textit{{k}}-\textit{{q}},a)\theta^{(1)}(\textit{{q}},a)\theta^{(2)}(\textit{{k}}-\textit{{q}},a)-\int\frac{\text{d}^{3}q}{(2\pi)^{3}}\tilde{\beta}(\textit{{q}},\textit{{k}}-\textit{{q}},a)\theta^{(2)}(\textit{{q}},a)\theta^{(1)}(\textit{{k}}-\textit{{q}},a).$
(166)

We can write Eq. (165) as

$\begin{split}S^{(3)}_{\tilde{\alpha}}(\textit{{k}},a)=\mathcal{H}fD^{3}(a)(A_{1}+A_{2}),\end{split}$
(167)

where

$\displaystyle A_{1}$ $\displaystyle=$
$\displaystyle\int\frac{\text{d}^{3}q}{(2\pi)^{3}}\frac{\text{d}^{3}q_{1}}{(2\pi)^{3}}\frac{\text{d}^{3}q_{2}}{(2\pi)^{3}}\delta^{(1)}(\textit{{q}})\delta^{(1)}(\textit{{q}}_{1})\delta^{(1)}(\textit{{q}}_{2})\tilde{\alpha}(\textit{{q}},\textit{{k}}-\textit{{q}})\textit{F}_{2}(\textit{{q}}_{1},\textit{{q}}_{2})(2\pi)^{3}\delta^{(D)}(\textit{{k}}-\textit{{q}}-\textit{{q}}_{1}-\textit{{q}}_{2}),$
(168) $\displaystyle A_{2}$ $\displaystyle=$
$\displaystyle\int\frac{\text{d}^{3}q}{(2\pi)^{3}}\frac{\text{d}^{3}q_{1}}{(2\pi)^{3}}\frac{\text{d}^{3}q_{2}}{(2\pi)^{3}}\delta^{(1)}(\textit{{q}})\delta^{(1)}(\textit{{q}}_{1})\delta^{(1)}(\textit{{q}}_{2})\tilde{\alpha}(\textit{{q}},\textit{{k}}-\textit{{q}})\textit{G}_{2}(\textit{{q}}_{1},\textit{{q}}_{2})(2\pi)^{3}\delta^{(D)}(\textit{{q}}-\textit{{q}}_{1}-\textit{{q}}_{2}).$
(169)

and Eq. (166) as

$\begin{split}S^{(3)}_{\tilde{\beta}}(\textit{{k}},a)=-\mathcal{H}^{2}f^{2}D^{3}(a)(B_{1}+B_{2}),\end{split}$
(170)

where

$\displaystyle B_{1}$ $\displaystyle=$
$\displaystyle\int\frac{\text{d}^{3}q}{(2\pi)^{3}}\frac{\text{d}^{3}q_{1}}{(2\pi)^{3}}\frac{\text{d}^{3}q_{2}}{(2\pi)^{3}}\delta^{(1)}(\textit{{q}})\delta^{(1)}(\textit{{q}}_{1})\delta^{(1)}(\textit{{q}}_{2})\tilde{\beta}(\textit{{q}},\textit{{k}}-\textit{{q}})\textit{G}_{2}(\textit{{q}}_{1},\textit{{q}}_{2})(2\pi)^{3}\delta^{(D)}(\textit{{k}}-\textit{{q}}-\textit{{q}}_{1}-\textit{{q}}_{2}),$
(171) $\displaystyle B_{2}$ $\displaystyle=$
$\displaystyle\int\frac{\text{d}^{3}q}{(2\pi)^{3}}\frac{\text{d}^{3}q_{1}}{(2\pi)^{3}}\frac{\text{d}^{3}q_{2}}{(2\pi)^{3}}\delta^{(1)}(\textit{{q}})\delta^{(1)}(\textit{{q}}_{1})\delta^{(1)}(\textit{{q}}_{2})\tilde{\beta}(\textit{{q}},\textit{{k}}-\textit{{q}})\textit{G}_{2}(\textit{{q}}_{1},\textit{{q}}_{2})(2\pi)^{3}\delta^{(D)}(\textit{{q}}-\textit{{q}}_{1}-\textit{{q}}_{2}).$
(172)

The right-hand side of Eq. (164) can then be written as

$\displaystyle
S^{(3)}_{\tilde{\beta}}-\mathcal{H}\partial_{a}\big{(}aS^{(3)}_{\tilde{\alpha}})=-\mathcal{H}^{2}\alpha^{2}a^{3\alpha}(B_{1}+B_{2})-\mathcal{H}^{2}a^{3\alpha}\Big{[}\frac{1}{2}\alpha+3\alpha^{2}\Big{]}(A_{1}+A_{2}),$

from which we obtain

$\displaystyle\delta^{(3)}(\textbf{k})=\frac{1}{2(8\alpha+1)}\bigg{[}2\alpha(B_{1}+B_{2})+(1+6\alpha)(A_{1}+A_{2})\bigg{]},$

which corresponds to

$\displaystyle\textit{F}_{3}(\textit{{q}}_{1},\textit{{q}}_{2},\textit{{q}}_{3})$
$\displaystyle=$
$\displaystyle\frac{1}{2(8\alpha+1)}\bigg{[}(6\alpha+1)\tilde{\alpha}(\textit{{q}}_{1},\textit{{q}}_{2}+\textit{{q}}_{3})\times\textit{F}_{2}(\textit{{q}}_{2},\textit{{q}}_{3})+(2\alpha)\tilde{\beta}(\textit{{q}}_{1},\textit{{q}}_{2}+\textit{{q}}_{3})\textit{G}_{2}(\textit{{q}}_{2},\textit{{q}}_{3})$
$\displaystyle+\big{[}(6\alpha+1)\tilde{\alpha}(\textit{{q}}_{1}+\textit{{q}}_{2},\textit{{q}}_{3})+(2\alpha)\tilde{\beta}(\textit{{q}}_{1}+\textit{{q}}_{2},\textit{{q}}_{3})\big{]}\textit{G}_{2}(\textit{{q}}_{1},\textit{{q}}_{2})\bigg{]}.$

###  A.2 LPT

The left-hand side of Eq. (60) at third order can be written

$\displaystyle\frac{(18\alpha^{2}+3\alpha)}{2}\Psi_{i,i}^{(3)}-(4\alpha^{2}+\alpha)\Psi_{i,j}^{(1)}\Psi_{i,j}^{(2)}-\frac{(2\alpha^{2}+\alpha)}{2}\Psi_{i,j}^{(2)}\Psi_{i,j}^{(1)}+(4\alpha^{2}+\alpha)\Psi_{i,i}^{(1)}\Psi_{j,j}^{(2)}$
$\displaystyle-\frac{(2\alpha^{2}+\alpha)}{2}\Psi_{i,i}^{(1)}\Psi_{i,j}^{(1)}\Psi_{i,j}^{(1)}+\frac{(2\alpha^{2}+\alpha)}{2}\Psi_{i,i}^{(2)}\Psi_{i,j}^{(1)}+\frac{(2\alpha^{2}+\alpha)}{4}\Big{\\{}\Psi_{i,i}^{(1)}\Psi_{j,j}^{(1)}\Psi_{k,k}^{(1)}-\Psi_{i,i}^{(1)}\Psi_{j,k}^{(1)}\Psi_{k,j}^{(1)}\Big{\\}}$
$\displaystyle+\frac{(2\alpha^{2}+\alpha)}{2}\Psi_{i,l}^{(1)}\Psi_{l,j}^{(1)}\Psi_{i,j}^{(1)},$
(176)

and the right-hand side as

$\displaystyle\frac{(2\alpha^{2}+\alpha)}{2}\Psi_{i,i}^{(3)}+\frac{(2\alpha^{2}+\alpha)}{2}\big{\\{}\Psi_{i,i}^{(2)}\Psi_{j,j}^{(1)}-\Psi_{i,j}^{(2)}\Psi_{j,i}^{(1)}\big{\\}}+\frac{(2\alpha^{2}+\alpha)}{2}\Big{\\{}\frac{1}{6}\Psi_{i,i}^{(1)}\Psi_{j,j}^{(1)}\Psi_{k,k}^{(1)}$
$\displaystyle-\frac{1}{2}\Psi_{i,i}^{(1)}\Psi_{j,k}^{(1)}\Psi_{k,j}^{(1)}+\frac{1}{3}\Psi_{i,l}^{(1)}\Psi_{l,j}^{(1)}\Psi_{i,j}^{(1)}\Big{\\}}.$
(177)

Collecting up the terms we have

$\displaystyle(8\alpha+1)\Psi_{i,i}^{(3)}$ $\displaystyle=$
$\displaystyle-(4\alpha+1)\Psi_{i,i}^{(2)}\Psi_{j,j}^{(1)}+(4\alpha+1)\Psi_{i,j}^{(2)}\Psi_{j,i}^{(1)}+\frac{(2\alpha+1)}{2}\Psi_{i,i}^{(1)}\Psi_{i,j}^{(1)}\Psi_{j,i}^{(1)}-\frac{(2\alpha+1)}{6}\Psi_{i,i}^{(1)}\Psi_{j,j}^{(1)}\Psi_{k,k}^{(1)}$
(178)
$\displaystyle-\frac{(2\alpha+1)}{3}\Psi_{i,k}^{(1)}\Psi_{k,j}^{(1)}\Psi_{j,i}^{(1)},$

which gives Eq. (62).

##  Appendix B Ultra-violet divergences in one-loop PS for gEdS

We detail here a little more the analysis leading to the expressions Eq. (83)
and Eq. (84) which allows us to infer the ultraviolet convergence properties
of these integrals i.e. the asymptotic large $k$ behaviour of $P_{\rm lin}(k)$
for which they converge.

For the $P_{22}$ the leading divergence is obtained simply by replacing
$P_{\rm lin}(|\textbf{k}-\textbf{q}|)$ by $P_{\rm lin}(q)$ and taking the
leading term in the series expansion in $k/q$ around $(k/q)=0$:

$\displaystyle|F_{2}^{(s)}(\textbf{k}-\textbf{q},\textbf{q})|^{2}$
$\displaystyle=$
$\displaystyle\frac{1}{4}\frac{k^{4}}{q^{4}}+\Big{(}\frac{4\alpha+1}{6\alpha+1}\Big{)}\Big{[}(\mu^{2}-1)\frac{k^{4}}{q^{4}}\Big{]}+\Big{(}\frac{4\alpha+1}{6\alpha+1}\Big{)}^{2}\Big{[}(\mu^{2}-1)^{2}\frac{k^{4}}{q^{4}}\Big{]}+\dots$
$\displaystyle=$
$\displaystyle\Big{(}\frac{2\alpha-(8\alpha+2)\mu^{2}+1}{2(6\alpha+1)}\Big{)}^{2}\frac{k^{4}}{q^{4}}+\cdots$

where $\mu=\frac{\textbf{k}\cdot\textbf{q}}{kq}$ is the angular variable.
Using Eqs. (73) and (B) we then obtain, after integrating $\mu$, Eq. (83).

Likewise expanding $\int_{-1}^{1}d\mu
F_{3}^{(s)}(\textbf{k},\textbf{q},-\textbf{q})$ in a power series in $k/q$
around $(k/q)=0$ we obtain

$\displaystyle\lim_{k/q\to 0}\int_{-1}^{1}d\mu
F_{3}^{(s)}(\textbf{k},\textbf{q},-\textbf{q})$ $\displaystyle=$
$\displaystyle-\frac{1}{9}\frac{k^{2}}{q^{2}}+\Big{(}\frac{2\alpha+1}{8\alpha+1}\Big{)}\Big{[}-\frac{4}{105}\frac{k^{4}}{q^{4}}+\times\frac{4}{15}\frac{k^{2}}{q^{2}}\Big{]}+\Big{(}\frac{2\alpha+1}{8\alpha+1}\Big{)}\Big{(}\frac{2\alpha}{6\alpha+1}\Big{)}\Big{[}\frac{4}{15}\frac{k^{4}}{q^{4}}-\frac{4}{9}\frac{k^{2}}{q^{2}}\Big{]}+\dots$
(180) $\displaystyle=$
$\displaystyle\frac{7-14\alpha-176\alpha^{2}}{45(1+6\alpha)(1+8\alpha)}\frac{k^{2}}{q^{2}}+\frac{4(8\alpha-1)(1+2\alpha)}{105(1+6\alpha)(1+8\alpha)}\frac{k^{4}}{q^{4}}+\dots$

which leads directly to Eq. (84).

##  Appendix C Summary of calculation of cosmology dependent corrections to
the EdS approximation for the one loop PS in standard perturbation theory

We have shown that this correction is given by

$\displaystyle\Delta
P_{1-loop}=\frac{2}{49}\gamma_{2}(M_{1}+\frac{10}{7}M_{2})+(\frac{2}{49}\gamma_{2})^{2}M_{2}+\frac{4}{49}\big{[}\gamma_{2}N_{1}+\frac{5}{54}\gamma_{3}(N_{2}-N_{1})\big{]},$
(181)

where the four integrals are given in Eqs. (77) and (81). These expressions
are valid for a generic FLRW cosmology with a clustering matter component and
a smooth component that scales at asymptotically early times as matter i.e.
with $w\rightarrow 0$ as $a\rightarrow 0$ (where
$w=\frac{1}{3}\frac{d\Omega_{m}(a)}{d\ln a}$). The two functions
$\gamma_{2}(a)$ and $\gamma_{3}(a)$ are obtained by solving the coupled
equations:

$\displaystyle\frac{d\alpha_{1}}{d\ln a}$ $\displaystyle=$
$\displaystyle\frac{3}{2}w\alpha_{1}-\alpha_{1}^{2}-\frac{1}{2}\alpha_{1}+\frac{3}{2}\Omega_{m},$
(182) $\displaystyle\frac{d^{2}\gamma_{2}}{d\ln a^{2}}$ $\displaystyle=$
$\displaystyle-\Big{[}\frac{1}{2}(1-3w)+4\alpha_{1}\Big{]}\frac{d\gamma_{2}}{d\ln
a}-\Big{(}2\alpha_{1}^{2}+\frac{3}{2}\Omega_{m}\Big{)}\gamma_{2}+\frac{21}{2}(\Omega_{m}-\alpha_{1}^{2}),$
(183) $\displaystyle\frac{d^{2}\gamma_{3}}{d\ln a^{2}}$ $\displaystyle=$
$\displaystyle-\Big{[}\frac{1}{2}(1-3w)+6\alpha_{1}\Big{]}\frac{d\gamma_{3}}{d\ln
a}-3(2\alpha_{1}^{2}+\Omega_{m})\gamma_{3}-\frac{189}{5}\Big{(}\Omega_{m}-\alpha_{1}^{2}\Big{)}+\frac{189}{35}\Big{[}4\alpha_{1}^{2}\gamma_{2}+2\alpha_{1}\frac{d\gamma_{2}}{d\ln
a}\Big{]},$ (184)

subject to the boundary conditions, given at $a=0$, by

$\alpha_{1}(0)=\alpha_{2}(0)=\alpha_{3}(0)=\alpha_{10}(0),\,\frac{d\alpha_{2}}{da}(0)=\frac{d\alpha_{3}}{da}(0)=0,$
(185)

or, more explicitly,

$\displaystyle\gamma_{2}(0)$ $\displaystyle=$ $\displaystyle
7\frac{1-\alpha_{10}(0)}{1+6\alpha_{10}(0)},$ $\displaystyle\gamma_{3}(0)$
$\displaystyle=$
$\displaystyle\frac{63}{5}\frac{(1-\alpha_{10}(0))(6\alpha_{10}(0)-1)}{(1+6\alpha_{10}(0))1+8\alpha_{10}(0))},$
(186)

and

$\frac{d\gamma_{2}}{da}(0)=0\,\,,\frac{d\gamma_{3}}{da}(0)=0,$ (187)

and where $\alpha_{10}(0)$ is given by Eq. (91) with $\Omega_{m}$ given by the
asymptotic (constant) fraction of clustering matter at high redshift. Given
that $\alpha_{1}=\alpha_{2}=\alpha_{3}(0)=\alpha_{10}(a)$ is the solution at
leading order in $w$, we can integrate in practice from $0<a_{i}\ll 1$
assuming exactly the initial conditions given by this solution. Note that in a
standard EdS cosmology we have of course $\alpha_{10}(0)=1$, and the case
$\alpha_{10}(0)\neq 1$ corresponds to the more general gEdS case where part of
the matter component at high redshift is smooth (e.g. a massive neutrino
component or matter-like dark energy component).

##  Appendix D Detailed comparison with results of [27]

Expressions for the cosmological corrections, at second and third order in
EPT, and for the one-loop PS (as in the previous appendix), have been given in
[27], in terms of six time-dependent functions. Further this paper provides
phenomenological fits for the parameter dependence in the case of a smooth
component with a constant equation of state. We make explicit here how these
six functions are related to one another by the three additional constraints
we have derived (which reduces the number of independent functions to three,
of which only two are required to calculate the PS at one loop). We also check
that there is indeed good numerical agreement between our results and those of
[27].

###  D.1 Second-order growth factor

[27] gives results for the second-order growth factor in terms of two
numerically fitted functions as

$\displaystyle T_{2A}$ $\displaystyle=$
$\displaystyle\frac{D_{2A}}{D_{1}^{2}}-1=|\ln{\Omega_{m}}|\Big{(}\frac{5.54\times
10^{-3}}{|w_{0}|}-\frac{3.40\times 10^{-3}}{\sqrt{|w_{0}|}}\Big{)},$ (188)
$\displaystyle T_{2B}$ $\displaystyle=$
$\displaystyle\frac{D_{2B}}{D_{1}^{2}}-1=|\ln{\Omega_{m}}|\Big{(}-\frac{1.384\times
10^{-2}}{|w_{0}|}+\frac{8.50\times 10^{-3}}{\sqrt{|w_{0}|}}\Big{)}.$ (189)

These are related (cf. Eqs. (108)) to our effective exponent $\gamma_{2}$ by

$T_{2A}=\frac{2}{35}\gamma_{2},\qquad T_{2B}=-\frac{1}{7}\gamma_{2},$ (190)

and thus $T_{2A}=-(2/5)T_{2B}$. Comparing the coefficients we see that this
relation is indeed verified within the level of numerical precision quoted by
[27].444What we have denoted by $w_{0}$ here corresponds to $w$ in [27].

Figure 9: The ratio of our numerical solutions for $\gamma_{2}$ to the
parametric fit of [27] $\gamma_{2,Tak}$ for the same quantity, for dark energy
with constant equation of state $w_{0}=-0.5,-1.0$, and $-1.5$.

Figure 10: Left panel: Ratio of the function $\gamma_{3AA}$ given by our
direct numerical solution of Eq. (142) to its value inferred from the fits of
[27], for a smooth component with constant equation of state ($w_{0}=-0.5,-1$
and $-1.5$). Right panel: Same comparison, but only for the model with
$w_{0}=-1$, and for $\gamma_{3}=\gamma_{3BA}$ (from numerical solution of Eq.
(184)). The two different curves correspond to the two different possible
numerical fits which can be inferred from those of [27]: the full line
corresponds to that inferred from $T_{3BA}$ directly, and the dashed line to
that obtained using the linear combination on the right-hand side of Eq.
(200)). The level of accuracy is less good than for $\gamma_{2}$ and
$\gamma_{3AA}$, but quite adequate for calculation notably of the PS at any
practically relevant level of precision.

Fig. 9 shows the ratio between the results we obtain for $\gamma_{2}$ by
numerical integration of Eqs.(182) and (183) with the fit for the same
quantity inferred from [27]:

$\displaystyle\gamma_{2,Tak}=(35/2)T_{2A}\approx|\ln{\Omega_{m}}|\Big{(}\frac{9.7\times
10^{-2}}{|w_{0}|}-\frac{5.95\times 10^{-2}}{\sqrt{|w_{0}|}}\Big{)},$

We see that our numerical solution for $\gamma_{2}$ well agree very accurately
with Eq. (D.1) for $w_{0}=-0.5,-1$ and well, albeit slightly less so, for
$w_{0}=-1.5$.

###  D.2 Third-order growth factor

At this order [27] gives results in terms of the four functions

$\displaystyle T_{3AA}=\frac{D_{3AA}}{D_{1}^{3}}-1$ $\displaystyle\simeq$
$\displaystyle|\ln{\Omega_{m}}|\Big{(}\frac{8.21\times
10^{-3}}{|w_{0}|}-\frac{5.14\times 10^{-3}}{\sqrt{|w_{0}|}}\Big{)},$ (192)
$\displaystyle T_{3AB}=\frac{D_{3AB}}{D_{1}^{3}}-1$ $\displaystyle\simeq$
$\displaystyle|\ln{\Omega_{m}}|^{1.5+0.4\ln{|w_{0}|}}|\Omega_{m}|^{0.7|w_{0}|}\Big{(}-\frac{9.16\times
10^{-3}}{|w_{0}|}+\frac{8.95\times 10^{-3}}{\sqrt{|w_{0}|}}\Big{)},$ (193)
$\displaystyle T_{3BA}=\frac{D_{3BA}}{D_{1}^{3}}-1$ $\displaystyle\simeq$
$\displaystyle|\ln{\Omega_{m}}|^{1.06-0.5\ln{|w_{0}|}}\Big{(}7.68\times
10^{-3}|w_{0}|-1.130\times 10^{-2}\sqrt{|w_{0}|}\Big{)},$ (194) $\displaystyle
T_{3BB}=\frac{D_{3BB}}{D_{1}^{3}}-1$ $\displaystyle\simeq$
$\displaystyle|\ln{\Omega_{m}}|\Big{(}-\frac{2.641\times
10^{-2}}{|w_{0}|}+\frac{1.582\times 10^{-2}}{\sqrt{|w_{0}|}}\Big{)},$ (195)

and two additional functions $D_{3AA}^{\prime}$ and $D_{3AB}^{\prime}$
determined by the relations in Eqs. (IV.3).

These are related to the functions we have introduced by

$\displaystyle T_{3AA}$ $\displaystyle=$
$\displaystyle\frac{4}{45}\gamma_{3AA},\quad
T_{3AB}=-\frac{1}{9}\gamma_{3AB},$ $\displaystyle T_{3BA}$ $\displaystyle=$
$\displaystyle\frac{5}{63}\gamma_{3BA},\quad
T_{3BB}=-\frac{16}{63}\gamma_{3BB}.$ (196)

From Eqs. (130)-(131) it is straightforward to show that our additional
constraints relating to these functions are

$\displaystyle 98\gamma_{3AA}-49\gamma_{3AB}+15\gamma_{3BA}-64\gamma_{3BB}=0,$
(197) $\displaystyle 49\gamma_{3AA}+32\gamma_{3BB}=81\gamma_{2}.$ (198)

Using these we can infer the relations

$\displaystyle T_{3BB}$ $\displaystyle=$
$\displaystyle\frac{9}{2}T_{2B}+\frac{35}{8}T_{3AA},$ (199) $\displaystyle
T_{3BA}$ $\displaystyle=$
$\displaystyle-\big{(}\frac{4}{3}T_{3BB}+\frac{35}{6}T_{3AA}+\frac{7}{3}T_{3AB}\big{)}.$
(200)

For the first relation it is easy to check the accuracy with which it is
satisfied by the expressions above from [27], because the same functional form
has been used to fit the functions on both sides: comparing the two fitted
coefficients on the left and right-hand side of the relation Eq. (199), we
find they agree at the sub-percent level. Further to check the numerical
accuracy of these fits against ours we plot in the left panel of Fig. 10 the
ratio between our own numerical solution for $\gamma_{3AA}$ and
$\gamma_{3AA,Tak}=(45/4)T_{3AA}$. The agreement is at a comparable level to
that for $\gamma_{2}$.

To assess both the agreement of Eq. (200) with the fits of [27] given above,
and the accuracy of these fits in describing the exact solutions, we plot in
the right panel of Fig. 10 the ratio of our numerical solutions for
$\gamma_{3}=\gamma_{3BA}$ to $\gamma_{3BA,Tak1}=(63/5)T_{3BA}$ and
$\gamma_{3BA,Tak2}=(63/5)T_{3BA}^{\prime}$ respectively, where
$T_{3BA}^{\prime}$ is the linear combination on the right-hand side of Eq.
(200). These curves are for a model with $w_{0}=-1$. We can infer that Eq.
(200) is obeyed by the fitting formulae at the $10\%$ level, and that the
agreement with our numerical results is slightly better if we use the fit for
$T_{3BA}$. We have checked carefully the origin of these small discrepancies
and find that they arise from transients from initial conditions which
disappear. In practice however, notably for the PS corrections which depend
only very weakly on $\gamma_{3}$, the fits provided by [27] are accurate
enough for any current practical application.

##  Appendix E Comparison with analysis of [32]

Figure 11: Our numerical solutions for $\gamma_{2}$, and
$\gamma_{3}=\gamma_{3BA}$ (solid black lines) in LCDM ($w_{0}=-1$) compared
with the power series expressions in Eqs. (203)-(204) obtained from [32]. The
different lines correspond to including terms in the power series up to the
indicated order.

In [32] the time dependence of the second order kernel is parametrized by two
functions $\lambda_{2}^{(1)}$ and $\lambda_{2}^{(2)}$, which are identical to
the two functions we have used:

$\displaystyle\lambda_{2}^{(1)}=d_{2A},\quad\lambda_{2}^{(2)}=d_{2B}.$ (201)

The six functions used to parametrize the third-order density kernels just
differ by a factor of $2$:

$\displaystyle\lambda_{3}^{(1)}=\frac{1}{2}d_{3AA},\quad\lambda_{3}^{(2)}=\frac{1}{2}d_{3AB},$
$\displaystyle\lambda_{3}^{(3)}=\frac{1}{2}d_{3AA}^{\prime},\quad\lambda_{3}^{(4)}=\frac{1}{2}d_{3AB}^{\prime},$
$\displaystyle\lambda_{3}^{(5)}=\frac{1}{2}d_{3BA},\quad\lambda_{3}^{(6)}=\frac{1}{2}d_{3BB}.$
(202)

For the specific case of a $LCDM$ cosmology, [32] provide analytical solutions
for these eight functions (and further ones parametrizing the fourth and fifth
order kernels) as power series in the parameter $\zeta=\frac{\Omega_{\Lambda
0}}{\Omega_{m0}}e^{3\eta}$, where $\eta\equiv\ln D_{+}$ ($D_{+}$ the linear
growth factor normalized to unity at $z=0$, and $\Omega_{\Lambda 0}$ and
$\Omega_{m0}$ are the dark energy and matter fractions at $z=0$,
respectively). It is straightforward to verify that the expressions given in
Appendix A of [32], up to third order in $\eta$, satisfy the five relations
Eq. (106) and Eqs. (137)-(140).

To check the numerical agreement of our results with those of [32] we
calculate the functions we have used to parametrize the PS and find

$\displaystyle\gamma_{2}$ $\displaystyle=$
$\displaystyle\frac{7}{2}(7\lambda_{2}^{(1)}-5)=-\frac{7c_{1}}{26}\zeta-\frac{2c_{2}}{19}\zeta^{2}-\frac{7c_{3}}{125}\zeta^{3},$
(203) $\displaystyle\gamma_{3BA}$ $\displaystyle=$
$\displaystyle\frac{63}{5}(21\lambda_{3}^{(5)}-1)=\frac{462c_{1}}{1625}\zeta+\frac{51c_{2}}{266}\zeta^{2}+\frac{2576c_{3}}{20625}\zeta^{3},$
(204)

where

$c_{1}=-\frac{3}{32},\quad c_{2}=-\frac{141}{4114},\quad
c_{3}=-\frac{9993}{1040842}$ (205)

Figure 11 shows the comparison between our numerical solutions to the exact
equations for the functions $\gamma_{2}$ and $\gamma_{3}=\gamma_{3BA}$ and the
power series solutions of [32] for the same quantities, in LCDM models, up to
the indicated orders. We observe excellent agreement where the power series
solutions appear to converge well.

##  Appendix F Comparison with analysis of [11, 20]

[11] parametrized the time dependence of second order and third order density
kernels in terms of three functions. Comparing their definitions directly to
ours (equivalent to those of [27, 32]) we obtain

$\displaystyle d_{2A}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}+\frac{3}{4}\nu_{2},$ (206) $\displaystyle d_{2B}$
$\displaystyle=$ $\displaystyle\frac{3}{2}-\frac{3}{4}\nu_{2},$ (207)
$\displaystyle d_{3AA}$ $\displaystyle=$
$\displaystyle-\frac{\lambda_{3}}{2}-\frac{3\nu_{2}}{4}+\frac{3\nu_{3}}{8}+\frac{1}{2},$
(208) $\displaystyle d_{3AA}^{\prime}$ $\displaystyle=$
$\displaystyle\frac{\lambda_{3}}{6}-\frac{3\nu_{2}}{4}+\frac{3\nu_{3}}{8}+\frac{1}{6},$
(209) $\displaystyle d_{3AB}$ $\displaystyle=$
$\displaystyle\frac{7\lambda_{3}}{6}+\frac{3\nu_{2}}{4}-\frac{3\nu_{3}}{8}+\frac{1}{6},$
(210) $\displaystyle d_{3AB}^{\prime}$ $\displaystyle=$
$\displaystyle\frac{\lambda_{3}}{2}+\frac{3\nu_{2}}{4}-\frac{3\nu_{3}}{8}+\frac{1}{2},$
(211) $\displaystyle d_{3BA}$ $\displaystyle=$
$\displaystyle-\frac{\lambda_{3}}{6}+\frac{9\nu_{2}}{4}-\frac{3\nu_{3}}{8}-\frac{13}{6},$
(212) $\displaystyle d_{3BB}$ $\displaystyle=$
$\displaystyle-\frac{\lambda_{3}}{2}-\frac{9\nu_{2}}{4}+\frac{3\nu_{3}}{8}+\frac{5}{2}.$
(213)

It is straightforward to verify that, using these expressions, the five
relations Eq. (106) and Eqs. (137)-(140) are satisfied. As we have shown in
Section IV, this corresponds to assuming that the eight functions obey
conditions valid for EdS boundary conditions, as assumed by [11]. As we have
seen in Section IV, these same conditions are also appropriate for the more
general case of gEdS boundary conditions.

The three effective growth rate functions ($\alpha_{2}$, $\alpha_{3AA}$, and
$\alpha_{3AB}$) that we have chosen to characterize the second order and third
order kernels can be inferred from

$\displaystyle\nu_{2}$ $\displaystyle=$
$\displaystyle\frac{2}{3}\left(2d_{2A}+1\right),$ (214) $\displaystyle\nu_{3}$
$\displaystyle=$ $\displaystyle\frac{2}{3}(4d_{2A}+7d_{3AA}+3d_{3AB}-2),$
(215) $\displaystyle\lambda_{3}$ $\displaystyle=$
$\displaystyle\frac{1}{2}(3d_{3AA}+3d_{3AB}-2).$ (216)

by substituting the appropriate definitions, Eq. (107) and Eqs. (LABEL:def-
gammas-3). In terms of the functions $\gamma_{2}$, $\gamma_{3AA}$, and
$\gamma_{3AB}$ we have

$\displaystyle\nu_{2}$ $\displaystyle=$
$\displaystyle\frac{2}{147}(119+4\gamma_{2}),$ (217) $\displaystyle\nu_{3}$
$\displaystyle=$
$\displaystyle\frac{2}{11907}(1372\gamma_{3AA}-294\gamma_{3AB}+648\gamma_{2}+21483),$
(218) $\displaystyle\lambda_{3}$ $\displaystyle=$
$\displaystyle\frac{1}{54}(4\gamma_{3AA}-2\gamma_{3AB}+9).$ (219)

[20] give (Appendix A) numerical values for $\nu_{2}$, $\nu_{3}$ and
$\lambda_{3}$ in an LCDM model with $\Omega_{m}=0.3$ at $z=0$. For this case
our numerical integration gives $\gamma_{2}=0.0452\dots$,
$\gamma_{3AA}=0.0417\dots$ and $\gamma_{3AB}=0.00142\dots$. Using Eqs.
(217)-(219) we find numerical agreement up to at least four significant
decimal places.

We note also that our result, that the one loop PS can be written in terms of
$d_{2A}$ and $d_{3BA}$ only, implies that it can equivalently be written only
in terms of $\nu_{2}$ and the linear combination $(4\lambda_{3}+9\nu_{3})$.
Finally we note that in the gEdS limit, corresponding to
$\alpha_{2}=\alpha_{3AA}=\alpha_{3AB}=\alpha$ we obtain

$\displaystyle\nu_{2}$ $\displaystyle=$
$\displaystyle\frac{2}{3}\Big{(}\frac{3+14\alpha}{1+6\alpha}\Big{)},$ (220)
$\displaystyle\nu_{3}$ $\displaystyle=$
$\displaystyle\frac{2}{3}\frac{\left(236\alpha^{2}+96\alpha+9\right)}{(1+6\alpha)(1+8\alpha)},$
(221) $\displaystyle\lambda_{3}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\Big{(}\frac{1+2\alpha}{1+8\alpha}\Big{)}.$ (222)

## References

  * Laureijs  _et al._ [2011] R. Laureijs, J. Amiaux, S. Arduini, J. L. Auguères, J. Brinchmann, R. Cole, M. Cropper, C. Dabin, L. Duvet, A. Ealet, _et al._ , arXiv e-prints , arXiv:1110.3193 (2011), arXiv:1110.3193 [astro-ph.CO] .
  * Dey  _et al._ [2019] A. Dey, D. J. Schlegel, D. Lang, R. Blum, K. Burleigh, X. Fan, J. R. Findlay, D. Finkbeiner, D. Herrera, S. Juneau, _et al._ , Astron. J. 157, 168 (2019), arXiv:1804.08657 [astro-ph.IM] .
  * Ivezić  _et al._ [2019] Ž. Ivezić, S. M. Kahn, J. A. Tyson, B. Abel, E. Acosta, R. Allsman, D. Alonso, Y. AlSayyad, S. F. Anderson, J. Andrew, _et al._ , Astrophys. J. 873, 111 (2019), arXiv:0805.2366 [astro-ph] .
  * Angulo and Hahn [2022] R. E. Angulo and O. Hahn, Living rev. comput. astrophys. 8, 1 (2022), arXiv:2112.05165 [astro-ph.CO] .
  * Juszkiewicz [1981] R. Juszkiewicz, Mon. Not. Roy. Astron. Soc. 197, 931 (1981).
  * Vishniac [1983] E. T. Vishniac, Mon. Not. Roy. Astron. Soc. 203, 345 (1983).
  * Goroff  _et al._ [1986] M. H. Goroff, B. Grinstein, S. J. Rey, and M. B. Wise, Astrophys. J. 311, 6 (1986).
  * Suto and Sasaki [1991] Y. Suto and M. Sasaki, Phys. Rev. Lett. 66, 264 (1991).
  * Bertschinger and Jain [1994] E. Bertschinger and B. Jain, Astrophys. J. 431, 486 (1994), arXiv:astro-ph/9307033 [astro-ph] .
  * Bernardeau [1994a] F. Bernardeau, Astrophys. J. 427, 51 (1994a), arXiv:astro-ph/9311066 [astro-ph] .
  * Bernardeau [1994b] F. Bernardeau, Astrophys. J. 433, 1 (1994b), arXiv:astro-ph/9312026 [astro-ph] .
  * Catelan  _et al._ [1995] P. Catelan, F. Lucchin, S. Matarrese, and L. Moscardini, Mon. Not. Roy. Astron. Soc. 276, 39 (1995), arXiv:astro-ph/9411066 [astro-ph] .
  * Makino  _et al._ [1992] N. Makino, M. Sasaki, and Y. Suto, Phys. Rev. D 46, 585 (1992).
  * Bernardeau  _et al._ [2002] F. Bernardeau, S. Colombi, E. Gaztañaga, and R. Scoccimarro, Phys. Rept. 367, 1 (2002), arXiv:astro-ph/0112551 [astro-ph] .
  * Crocce and Scoccimarro [2006] M. Crocce and R. Scoccimarro, Phys. Rev. D 73, 063519 (2006), arXiv:astro-ph/0509418 [astro-ph] .
  * Carlson  _et al._ [2009] J. Carlson, M. White, and N. Padmanabhan, Phys. Rev. D 80, 043531 (2009), arXiv:0905.0479 [astro-ph.CO] .
  * Taruya [2016] A. Taruya, Phys. Rev. D 94, 023504 (2016), arXiv:1606.02168 [astro-ph.CO] .
  * Modi  _et al._ [2017] C. Modi, M. White, and Z. Vlah, JCAP 2017 (8), 009, arXiv:1706.03173 [astro-ph.CO] .
  * Chen  _et al._ [2020] S.-F. Chen, Z. Vlah, and M. White, JCAP 2020 (7), 062, arXiv:2005.00523 [astro-ph.CO] .
  * Garny and Taule [2021] M. Garny and P. Taule, JCAP 2021 (1), 020, arXiv:2008.00013 [astro-ph.CO] .
  * Chen  _et al._ [2021] S.-F. Chen, Z. Vlah, and M. White, JCAP 2021 (5), 053, arXiv:2103.13498 [astro-ph.CO] .
  * Chen  _et al._ [2022] S.-F. Chen, M. White, J. DeRose, and N. Kokron, arXiv e-prints , arXiv:2204.10392 (2022), arXiv:2204.10392 [astro-ph.CO] .
  * McDonald [2007] P. McDonald, Phys. Rev. D 75, 043514 (2007), arXiv:astro-ph/0606028 [astro-ph] .
  * Carrasco  _et al._ [2012] J. J. M. Carrasco, M. P. Hertzberg, and L. Senatore, JHEP 2012, 82, arXiv:1206.2926 [astro-ph.CO] .
  * Pajer and Zaldarriaga [2013] E. Pajer and M. Zaldarriaga, JCAP 2013 (8), 037, arXiv:1301.7182 [astro-ph.CO] .
  * Steele and Baldauf [2021] T. Steele and T. Baldauf, Phys. Rev. D 103, 023520 (2021), arXiv:2009.01200 [astro-ph.CO] .
  * Takahashi [2008] R. Takahashi, Prog. Theor. Phys. 120, 549 (2008), arXiv:0806.1437 [astro-ph] .
  * Fasiello and Vlah [2016] M. Fasiello and Z. Vlah, Phys. Rev. D 94, 063516 (2016), arXiv:1604.04612 [astro-ph.CO] .
  * Aviles  _et al._ [2018] A. Aviles, M. A. Rodriguez-Meza, J. De-Santiago, and J. L. Cervantes-Cota, JCAP 2018 (11), 013, arXiv:1809.07713 [astro-ph.CO] .
  * Baldauf  _et al._ [2021] T. Baldauf, M. Garny, P. Taule, and T. Steele, Phys. Rev. D 104, 123551 (2021), arXiv:2110.13930 [astro-ph.CO] .
  * Alkhanishvili  _et al._ [2022] D. Alkhanishvili, C. Porciani, E. Sefusatti, M. Biagetti, A. Lazanu, A. Oddo, and V. Yankelevich, Mon. Not. R. Astron. Soc. 512, 4961 (2022), arXiv:2107.08054 [astro-ph.CO] .
  * Fasiello  _et al._ [2022] M. Fasiello, T. Fujita, and Z. Vlah, Phys. Rev. D 106, 123504 (2022), arXiv:2205.10026 [astro-ph.CO] .
  * Benhaiem [2013] D. Benhaiem,  _Non linear gravitational clustering in scale free cosmological models_ , Ph.D. thesis, Université Pierre et Marie Curie (2013).
  * Benhaiem  _et al._ [2014] D. Benhaiem, M. Joyce, and B. Marcos, Mon. Not. R. Astron. Soc. 443, 2126 (2014), arXiv:1309.2753 [astro-ph.CO] .
  * Joyce  _et al._ [2021] M. Joyce, L. Garrison, and D. Eisenstein, Mon. Not. R. Astron. Soc. 501, 5051 (2021), arXiv:2004.07256 [astro-ph.CO] .
  * Maleubre  _et al._ [2022] S. Maleubre, D. Eisenstein, L. H. Garrison, and M. Joyce, Mon. Not. R. Astron. Soc. 512, 1829 (2022), arXiv:2109.04397 [astro-ph.CO] .
  * Ferreira and Joyce [1998] P. G. Ferreira and M. Joyce, Phys. Rev. D 58, 023503 (1998), arXiv:astro-ph/9711102 [astro-ph] .
  * Martel and Freudling [1991] H. Martel and W. Freudling, Astrophys. J. 371, 1 (1991).
  * Scoccimarro  _et al._ [1998] R. Scoccimarro, S. Colombi, J. N. Fry, J. A. Frieman, E. Hivon, and A. Melott, Astrophys. J. 496, 586 (1998), arXiv:astro-ph/9704075 [astro-ph] .
  * Bouchet  _et al._ [1995] F. R. Bouchet, S. Colombi, E. Hivon, and R. Juszkiewicz, Astron. Astrophys. 296, 575 (1995), arXiv:astro-ph/9406013 [astro-ph] .
  * Matsubara [2011] T. Matsubara, Phys. Rev. D 83, 083518 (2011), arXiv:1102.4619 [astro-ph.CO] .
  * Rampf and Buchert [2012] C. Rampf and T. Buchert, JCAP 2012 (6), 021, arXiv:1203.4260 [astro-ph.CO] .
  * Scoccimarro and Frieman [1996] R. Scoccimarro and J. Frieman, Astrophys. J. 105, 37 (1996), arXiv:astro-ph/9509047 [astro-ph] .
  * Peloso and Pietroni [2013] M. Peloso and M. Pietroni, JCAP 2013 (5), 031, arXiv:1302.0223 [astro-ph.CO] .
  * Carrasco  _et al._ [2014] J. J. M. Carrasco, S. Foreman, D. Green, and L. Senatore, JCAP 2014 (7), 056, arXiv:1304.4946 [astro-ph.CO] .
  * Baldauf [2013] T. Baldauf,  _Modelling large scale structure statistics for precision cosmology_ , Ph.D. thesis, University of Zurich (2013).
  * Baldauf  _et al._ [2015] T. Baldauf, L. Mercolli, M. Mirbabayi, and E. Pajer, JCAP 2015 (5), 007, arXiv:1406.4135 [astro-ph.CO] .
  * Aviles and Banerjee [2020] A. Aviles and A. Banerjee, JCAP 2020 (10), 034, arXiv:2007.06508 [astro-ph.CO] .
  * Aviles  _et al._ [2021] A. Aviles, A. Banerjee, G. Niz, and Z. Slepian, JCAP 2021 (11), 028, arXiv:2106.13771 [astro-ph.CO] .
  * Planck Collaboration  _et al._ [2014] Planck Collaboration, P. A. R. Ade, N. Aghanim, C. Armitage-Caplan, M. Arnaud, M. Ashdown, F. Atrio-Barandela, J. Aumont, C. Baccigalupi, A. J. Banday, R. B. Barreiro, _et al._ , Astron. Astrophys. 571, A16 (2014), arXiv:1303.5076 [astro-ph.CO] .
  * Eisenstein and Hu [1998] D. J. Eisenstein and W. Hu, Astrophys. J. 496, 605 (1998), arXiv:astro-ph/9709112 [astro-ph] .
  * Garny and Taule [2022] M. Garny and P. Taule, JCAP 2022 (9), 054, arXiv:2205.11533 [astro-ph.CO] .
  * Hunter [2007] J. D. Hunter, Comput. Sci. Eng. 9, 90 (2007).
  * Harris  _et al._ [2020] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, _et al._ , Nature 585, 357 (2020), arXiv:2006.10256 [cs.MS] .
  * Virtanen  _et al._ [2020] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, _et al._ , Nat. Methods 17, 261 (2020), arXiv:1907.10121 [cs.MS] .