15(\bar{\cal E}_{4}-1)\,,$ (182) $\displaystyle\bar{\cal E}_{8}$
$\displaystyle\leq$ $\displaystyle 105-210\bar{\cal E}_{4}+35\bar{\cal
E}_{4}^{2}+28\bar{\cal E}_{6}\,,$ $\displaystyle\bar{\cal E}_{10}$
$\displaystyle\geq$ $\displaystyle 45\bar{\cal E}_{8}+210\bar{\cal
E}_{4}(15+\bar{\cal E}_{6})-945-1575\bar{\cal E}_{4}^{2}-630\bar{\cal
E}_{6}\,,$ $\displaystyle\bar{\cal E}_{12}$ $\displaystyle\leq$ $\displaystyle
10395+51975\bar{\cal E}_{4}^{2}-5775\bar{\cal E}_{4}^{3}+13860\bar{\cal
E}_{6}+462\bar{\cal E}_{6}^{2}$ $\displaystyle-495\bar{\cal
E}_{4}(105+28\bar{\cal E}_{6}-\bar{\cal E}_{8})-1485\bar{\cal
E}_{8}+66\bar{\cal E}_{10}\,,$

for all values of $c_{\text{max}}$ for which the corresponding expectation
values enter, being $c_{\text{max}}\geq 3,5,7,9,11$ for the five inequalities,
respectively. The first condition is consistent with Eq. (VII.4). These
conditions are however not sufficient.

A _sufficient_ set of stability conditions can be obtained by applying an
algorithm known as Sturm chain. To that end we define a set of polynomials
$P_{n}(\lambda)$, starting from

$P_{0}(\lambda)=q_{c_{\text{max}}}(\lambda),\quad
P_{1}(\lambda)=dq_{c_{\text{max}}}/d\lambda\,.$ (183)

Then, we recursively compute the polynomial quotient $Q_{n}(\lambda)$ from
polynomial division of $P_{n}$ by $P_{n+1}$, and define

$P_{n+2}=Q_{n}P_{n+1}-P_{n}\,,$ (184)

being the rest term up to an overall sign. The order of the polynomials
decreases with increasing $n$, and at some point one obtains a constant, where
the chain is terminated. For all $c_{\text{max}}$ we considered, this is the
case for $n$ being equal to the order of the polynomial
$q_{c_{\text{max}}}(\lambda)$, being

$N_{c_{\text{max}}}\equiv[(c_{\text{max}}+1)/2]\,.$ (185)

Here the square bracket denotes the integer part. The number of roots of
$q_{c_{\text{max}}}(\lambda)$ in the intervall $a<\lambda\leq b$ for some real
values $a<b$ is then given by $\sigma(a)-\sigma(b)$, where $\sigma(\lambda)$
is the number of sign changes in the series
$P_{0}(\lambda),P_{1}(\lambda),\dots,P_{N_{c_{\text{max}}}}(\lambda)$.

Stability requires that all $N_{c_{\text{max}}}$ roots of
$q_{c_{\text{max}}}(\lambda)$ lie in the intervall $-\infty<\lambda\leq 0$. We
therefore consider the choice $a\to-\infty$ and $b=0$, and require
$\sigma(a)-\sigma(b)=N_{c_{\text{max}}}$. The only way how this condition can
be satisfied is if $\sigma(a)=N_{c_{\text{max}}}$ and $\sigma(b)=0$. For
$\lambda\to-\infty$ the coefficient of the highest monomial in each $P_{n}$
determines its sign, and for $\lambda=0$ the constant term. We write the
polynomials in the form

$P_{n}(\lambda)=C_{n}\lambda^{N_{c_{\text{max}}}-n}+\dots+D_{n}\,,$ (186)

where the ellipsis denotes summands with powers $\lambda^{m}$ with
$0<m<N_{c_{\text{max}}}-n$, that are irrelevant here. Due to the alternating
sign of $\lambda^{N_{c_{\text{max}}}-n}$ for $n=0,\dots,N_{c_{\text{max}}}$
and $\lambda<0$, all coefficients $C_{n}$ are required to have the same sign
in order to satisfy $\sigma(-\infty)=N_{c_{\text{max}}}$. Furthermore, by
definition of the characteristic polynomial, we have $C_{0}=+1$. This implies
the conditions that $C_{n}\geq 0$ for all $n=0,\dots,N_{c_{\text{max}}}$.
Similarly, the constraint $\sigma(0)=0$ implies that all $D_{n}$ have to have
the same sign. In principle they could all be positive or all be negative.
However, the latter can be excluded by the following argument: For the
constant polynomial with $n=N_{c_{\text{max}}}$, the term with highest and
lowest power of $\lambda$ are trivially identical, i.e.
$C_{N_{c_{\text{max}}}}=D_{N_{c_{\text{max}}}}$. Since we already obtained the
condition that all $C_{n}$ need to be non-negative, this implies that also all
$D_{n}$ have to be non-negative for stability to hold,

$C_{n}\geq 0,\quad D_{n}\geq 0,\qquad 0\leq n\leq N_{c_{\text{max}}}\,.$ (187)

As mentioned above, the condition is trivially satisfied for $C_{0}=1$ and
degenerate for $C_{n}=D_{n}$ for $n=N_{c_{\text{max}}}$. There are therefore
in general $2N_{c_{\text{max}}}$ distinct conditions. However, some of them
are either trivially satisfied or equivalent. For example, for
$c_{\text{max}}=4$ one has $N_{c_{\text{max}}}=2$ and we find
$C_{1}=2,C_{2}=D_{2}=\frac{5}{9}(2+\bar{\cal E}_{4})$ and
$D_{0}=\frac{5}{9}(3-\bar{\cal E}_{4}),D_{1}=\frac{10}{3}$. The stability
conditions Eq. (187) therefore precisely yield the constraint Eq. (VII.4)
obtained previously for $c_{\text{max}}=4$. The same can be checked for
$c_{\text{max}}=3$.

In addition, the Sturm chain algorithm allows us to obtain stability
conditions for arbitrary truncation order $c_{\text{max}}$. We observe a
number of general patterns: We find that $C_{1}=N_{c_{\text{max}}}$, such that
the condition $C_{1}\geq 0$ is trivially satisfied. From the condition
$C_{2}\geq 0$ we obtain a lower bound on $\bar{\cal E}_{4}$ for all
$c_{\text{max}}\geq 3$,

$\bar{\cal E}_{4}\geq\left\\{\begin{array}[]{l@\qquad
l}-\frac{6}{c_{\text{max}}-2}&c_{\text{max}}\ \text{odd}\,,\\\\[6.45831pt]
-\frac{6}{c_{\text{max}}-1}&c_{\text{max}}\ \text{even}\,.\\\
\end{array}\right.$ (188)

This is a generalization of the lower bound obtained in Eq. (VII.4).222The
same condition can alternatively be derived from the Laguerre-Samuelson rule,
being that the generalized discriminant is positive,
$q_{c_{\text{max}}}^{(2)}q_{c_{\text{max}}}^{(0)}\frac{2N_{c_{\text{max}}}}{N_{c_{\text{max}}}-1}\leq\left(q_{c_{\text{max}}}^{(1)}\right)^{2}$,
where
$q_{c_{\text{max}}}(\lambda)=\sum_{n}q_{c_{\text{max}}}^{(n)}\lambda^{N_{c_{\text{max}}}-n}$.
However, the Laguerre-Samuelson rule is necessary but not sufficient for
stability, while the full set of conditions derived from the Sturm chain are
sufficient. We checked Eq. (188) up to $c_{\text{max}}=12$.

Furthermore, the condition obtained from $D_{0}\geq 0$ is equivalent to the
first line in Eq. (VII.4) for $c_{\text{max}}=3,4$, to the second line for
$c_{\text{max}}=5,6$, the third line for $c_{\text{max}}=7,8$, and so on.
Similarly, $D_{1}\geq 0$ yields the first line in Eq. (VII.4) for
$c_{\text{max}}=5,6$, the second line for $c_{\text{max}}=7,8$, etc.

Figure 7:  Constraints on the non-Gaussianity of the distribution function
from requiring stability of the linear solutions. The top panel shows the
parameter region allowed by stability within the plane spanned by the
dimensionless $4$th and $6$th cumulant expectation values $(\bar{\cal
E}_{4},\bar{\cal E}_{6})$, while the bottom panel shows constraints on the
$4$th and $8$th cumulant expectation values $(\bar{\cal E}_{4},\bar{\cal
E}_{8})$. In each case, these are obtained from the evolution equations for
perturbations modes of cumulants up to $c_{\text{max}}$ as given in the
legend. Note that $\bar{\cal E}_{6}$ is relevant for $c_{\text{max}}\geq 5$,
and $\bar{\cal E}_{8}$ for $c_{\text{max}}\geq 7$. The grey dashed and black
lines show the expectation from the Evans and NFW halo models, see Sec. II.
The fact that the backreaction on linear modes from dispersion and higher
cumulants expected from halos is broadly stable is reassuring, making
linearized a good starting point for a perturbative expansion.

The complete set of stability conditions obtained from the Sturm chain
algorithm for $c_{\text{max}}=5$ reads

$\displaystyle-2$ $\displaystyle\leq$ $\displaystyle\bar{\cal E}_{4}\leq 3\,,$
$\displaystyle 15(\bar{\cal E}_{4}-1)$ $\displaystyle\leq$
$\displaystyle\bar{\cal E}_{6}\leq 10(6-\bar{\cal E}_{4})\,,$ (189)
$\displaystyle 0$ $\displaystyle\leq$ $\displaystyle 100(24+12\bar{\cal
E}_{4}-6\bar{\cal E}_{4}^{2}+5\bar{\cal E}_{4}^{3})$
$\displaystyle-40\bar{\cal E}_{6}(2+3\bar{\cal E}_{4})-\bar{\cal
E}_{6}^{2}\,.$

For $c_{\text{max}}=6$ we obtain

$\displaystyle-\frac{6}{5}$ $\displaystyle\leq$ $\displaystyle\bar{\cal
E}_{4}\leq 3\,,$ $\displaystyle 15(\bar{\cal E}_{4}-1)$ $\displaystyle\leq$
$\displaystyle\bar{\cal E}_{6}\leq 10(2+\bar{\cal E}_{4}/3)\,,$ (190)
$\displaystyle 0$ $\displaystyle\leq$ $\displaystyle 20(216+324\bar{\cal
E}_{4}+90\bar{\cal E}_{4}^{2}+175\bar{\cal E}_{4}^{3})$
$\displaystyle-108\bar{\cal E}_{6}(4+10\bar{\cal E}_{4})-27\bar{\cal
E}_{6}^{2}\,.$

The conditions for $c_{\text{max}}=7,8$ are given in App. G. We checked that
vanishing expectation values $\bar{\cal E}_{2n}=0$ for $2n\geq 4$ do satisfy
the stability conditions in all cases. Therefore, in general, stability sets
an upper limit on the magnitude of the $\bar{\cal E}_{2n}$, i.e. how strongly
non-Gaussian the distribution function can be (on average).

Let us now discuss the impact of stability constraints. For
$c_{\text{max}}=1,2$ the stability conditions are trivially satisfied, and for
$c_{\text{max}}=3,4$ they amount to the constraint on the size of $\bar{\cal
E}_{4}$ given in Eq. (VII.4). For $c_{\text{max}}=5,6$ the perturbation
equations depend on $\bar{\cal E}_{4}$ and $\bar{\cal E}_{6}$, and the
stability constraints given in Eqs. (VII.4, VII.4) are satisfied within a
finite region in the two-dimensional $(\bar{\cal E}_{4},\bar{\cal E}_{6})$
parameter space, shown in Fig. 7 (top panel). We see that the point $\bar{\cal
E}_{4}=\bar{\cal E}_{6}=0$ lies within the stable region, as claimed aboved.
Furthermore, the allowed region is more restricted for $c_{\text{max}}=6$
compared to $c_{\text{max}}=5$.

For $c_{\text{max}}=7,8$, also $\bar{\cal E}_{8}$ enters in the perturbation
equations. Stability yields an allowed region in the three-dimensional
parameter space $(\bar{\cal E}_{4},\bar{\cal E}_{6},\bar{\cal E}_{8})$, that
contains the origin. In Fig. 7 we show the projection of this region on the
$(\bar{\cal E}_{4},\bar{\cal E}_{6})$ as well as $(\bar{\cal E}_{4},\bar{\cal
E}_{8})$ plane in the top and bottom panel, respectively. As expected, the
allowed regions for higher $c_{\text{max}}$ are contained in those for lower
values of $c_{\text{max}}$. Generically, one may expect the normalized
quantities $\bar{\cal E}_{2n}$ to be broadly of order unity, which is allowed
by stability. Nevertheless, the stability conditions impose relevant
constraints, in particular for $\bar{\cal E}_{4}$, as well as correlations
among the relative size of the $\bar{\cal E}_{2n}$.

It is instructive to compare the stability regions of the cumulant expectation
values with those obtained within the Evans and NFW halo models discussed in
Sec. II, that are shown as grey dashed and black lines in Fig. 7,
respectively. Since the stability conditions for the linearized Vlasov
hierarchy and the averaged cumulants obtained from virialized halos are based
on quite distinct physical situations, it is remarkable that the latter are
largely contained within the stable regions for a wide range of halo
concentrations (NFW halos) and different halo shapes (Evans halos). Note the
Evans halos hit the stability boundary when the shape parameter $q$ approaches
the unphysical limit of extreme oblateness $q=1/\sqrt{2}$, where the density
ceases to be positive definite. For NFW halos, the stability boundary is hit
for high mass halos, corresponding to high concentrations which are the most
non-Gaussian (See Fig. 2). But this is unrealistic since in practice we have a
spectrum of halos, and the abundance of high-mass halos is exponentially
suppressed, therefore one must integrate over the mass function to compare
properly. The advantage of using the results of individual halos in Fig. 7 is
that we are insensitive to the initial power spectrum shape, therefore these
results apply broadly to scale-free power spectra as well as CDM spectra. The
dependence on the initial spectrum enters only through the weight given to
different halo masses (or concentrations) by the mass function. That the
constrains on non-Gaussianity of the distribution function are satisfied by
halo estimates implies that the linear solutions within can be considered as a
good starting point for a perturbative analysis for realistic values of the
higher cumulant expectation values.

##  VIII Dispersion in a scaling universe

So far, we treated the background dispersion $\epsilon(\eta)={\cal E}_{2}$ as
well as the expectation value $\omega(\eta)=5{\cal E}_{4}/3$ of the fourth
cumulant, and those of yet higher cumulants (${\cal E}_{2n}$), as external
inputs for solving the equations for perturbation modes up to a certain
cumulant order. In this section we return to the Eqs. (47, V.1, 89, 159) for
the background quantities themselves, that are in turn sourced by the
fluctuations of the perturbation modes, and present self-consistent solutions
of the perturbation and background equations in various approximations. For
illustration, we restrict ourselves to a scaling universe in this work, for
which the differential equations for the background values turn into algebraic
equations, as we shall see. This makes the problem tractable and allows us to
study the dependence on the truncation of the cumulant expansion.

A scaling universe is characterized by a power-law initial spectrum,

$P_{0}(k)=A\,k^{n_{s}}\,,$ (191)

with spectral index $n_{s}$, and an EdS background ($\Omega_{m}=1$). The
linear power spectrum in the SPT approximation is given by

$P^{\text{lin}}_{\text{SPT}}(k,\eta)=e^{2\eta}P_{0}(k)\,,$ (192)

where $e^{2\eta}=D^{2}$ is the square of the conventional linear growth
factor, which in turn equals the scale-factor within EdS. For the
dimensionless power spectrum $\Delta\equiv 4\pi k^{3}P$ this means

$\Delta^{\text{lin}}_{\text{SPT}}(k,\eta)=4\pi
e^{2\eta}A\,k^{n_{s}+3}\equiv\left(\frac{k}{k_{nl}(\eta)}\right)^{n_{s}+3}\,,$
(193)

where we introduced the usual non-linear scale

$k_{nl}(\eta)=k_{nl}\,e^{-2\eta/(n_{s}+3)}\,,$ (194)

with $k_{nl}=(4\pi A)^{-1/(n_{s}+3)}$ being the non-linear scale today
($\eta=0$). The power spectrum obeys a scaling symmetry (for any $r>0$)

$k\to rk,\qquad e^{\eta}\to r^{-(n_{s}+3)/2}e^{\eta}\,,$ (195)

that suggests that the non-linear power spectrum (of any dimensionless
variable) is a function of the ratio

$\Delta(k,\eta)=\Delta(k/k_{nl}(\eta))\,.$ (196)

The background value $\epsilon(\eta)$ of the velocity dispersion defines the
scale

$k_{\sigma}(\eta)\equiv\frac{1}{\sqrt{\epsilon(\eta)}}\,.$ (197)

The scaling symmetry suggests that $k_{\sigma}(\eta)/k_{nl}(\eta)$ is
constant, implying that

$\epsilon=\epsilon_{0}e^{\alpha\eta}\,,$ (198)

follows a power-law with exponent

$\alpha=4/(n_{s}+3)\,.$ (199)

We denote the value today by

$k_{\sigma}\equiv 1/\sqrt{\epsilon_{0}}\,,$ (200)

without time-argument. In addition, when taking higher cumulants into account,
the dimensionless ratio $\bar{\cal E}_{4}=3\bar{\omega}/5$ as well as in
general all $\bar{\cal E}_{2n}$ are constant in time. As anticipated, the
scaling universe therefore provides an example for which the assumptions on
the time-dependence of $\epsilon(\eta)$ and higher cumulant expectation values
taken in Secs. VI and VII are satisfied exactly, and the corresponding linear
kernels for the perturbation variables can be used.

The linear matter power spectrum within is given by

$P_{\delta\delta}^{\text{lin}}(k,\eta)=F_{1,\delta}(k,\eta)^{2}\,e^{2\eta}P_{0}(k)\,,$
(201)

where $F_{1,\delta}$ is the linear kernel obtained from a solution of the
appropriate linearized perturbation equations. The (cross-)power spectra for
any pair $a,b=\delta,\theta,\bar{g},\delta\bar{\epsilon},\dots$ of
dimensionless perturbation modes in is analogously given by

$P_{ab}^{\text{lin}}(k,\eta)=F_{1,a}(k,\eta)F_{1,b}(k,\eta)\,e^{2\eta}P_{0}(k)\,.$
(202)

As shown in Secs. VI and VII, for any truncation order of the cumulant
expansion the linear kernels $F_{1,a}(k,\eta)=F_{1,a}(s)$ depend on time and
scale only via the dimensionless variable

$s\equiv s_{k}(\eta)=\sqrt{3\epsilon(\eta)k^{2}}\,.$ (203)

Therefore we may write the dimensionless power spectrum
$\Delta_{ab}^{\text{lin}}=4\pi k^{3}P_{ab}^{\text{lin}}$ as

$\Delta_{ab}^{\text{lin}}(k,\eta)=x\times\hat{\Delta}_{ab}^{\text{lin}}(s)\,,$
(204)

where

$\displaystyle x$ $\displaystyle\equiv$
$\displaystyle\left(\frac{k_{\sigma}}{\sqrt{3}k_{nl}}\right)^{n_{s}+3}=\left(\frac{1}{3\epsilon_{0}k_{nl}^{2}}\right)^{\frac{n_{s}+3}{2}}\,,$
(205)

is time-independent, and we defined

$\hat{\Delta}_{ab}^{\text{lin}}(s)\equiv F_{1,a}(s)F_{1,b}(s)s^{n_{s}+3}\,.$
(206)

A similar relation holds for the full non-linear power spectrum, computed
within . As in SPT, one can use the perturbative solutions to obtain a loop
expansion,

$\Delta_{ab}(k,\eta)=\sum_{L\geq 0}\Delta_{ab}^{L-\text{loop}}(k,\eta)\,,$
(207)

with $L=0$ being the linear solution. For
$ab=\delta\delta,\delta\theta,\theta\theta$, the loop corrections formally
take a similar form as in SPT, but with non-linear kernels computed based on
the perturbation modes and their non-linear vertices presented here. For a
detailed discussion of loop corrections within we refer to paper II Garny  _et
al._ (2022). Here we restrict ourselves to the general structure of loop
corrections. Since the $L$-loop contribution encompasses $L+1$ factors of
$P_{0}\propto A\propto x$, one has

$\Delta_{ab}^{L-\text{loop}}(k,\eta)=x^{L+1}\times\hat{\Delta}_{ab}^{L-\text{loop}}(s)\,,$
(208)

with the second factor involving loop integrals for $L\geq 1$.

Let us now discuss how to obtain a self-consistent solution for the background
dispersion $\epsilon(\eta)$. Its equation of motion is given in Eq. (47), with
source term from Eq. (V.2). It can be written in terms of dimensionless
quantities as

$\frac{\epsilon^{\prime}(\eta)}{\epsilon(\eta)}+1=\frac{1}{3}\int_{0}^{\infty}\frac{dk}{k}\,\left(\Delta_{\theta\bar{g}}-\Delta_{\theta\delta\bar{\epsilon}}+2\Delta_{w_{i}\bar{\nu}_{i}}+\Delta_{A\bar{\pi}}\right)\,,$
(209)

where $\bar{g}=g/\epsilon$, $\delta\bar{\epsilon}=\delta\epsilon/\epsilon$ and
$\bar{\nu}_{i}=\nu_{i}/\epsilon$ are the dimensionless scalar and vector
perturbation modes of the velocity dispersion tensor, $\bar{\pi}=\pi/\epsilon$
is a scalar mode of the third cumulant, and we have set
$\Omega_{m}/f^{2}\mapsto 1$.

For the time-dependence Eq. (198) expected from scaling symmetry, the left-
hand side of this equation is constant, and equal to
$\alpha+1=(n_{s}+7)/(n_{s}+3)$. Therefore, the ansatz Eq. (198) is consistent
if also the right-hand side is time-independent. Using Eq. (206) in linear
approximation and Eq. (208) in general, we see that this is indeed the case,
since the variable $x$ is time-independent and the power spectrum integrated
over $s\propto k$ as well. Therefore, Eq. (209) turns into an algebraic
equation for $x$ given by

$\frac{n_{s}+7}{n_{s}+3}=\frac{1}{3}\sum_{L\geq
0}x^{L+1}I^{L-\text{loop}}(n_{s})\,,$ (210)

where

$\displaystyle I^{L-\text{loop}}(n_{s})$ $\displaystyle\equiv$
$\displaystyle\int_{0}^{\infty}\frac{ds}{s}\Big{(}\hat{\Delta}_{\theta\bar{g}}(s)-\hat{\Delta}_{\theta\delta\bar{\epsilon}}(s)$
(211)
$\displaystyle+2\hat{\Delta}_{w_{i}\bar{\nu}_{i}}(s)+\hat{\Delta}_{A\bar{\pi}}(s)\Big{)}^{L-\text{loop}}\,.$

In particular, the linear contribution ($L=0$) reads

$\displaystyle I^{\text{lin}}(n_{s})$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}ds\,s^{n_{s}+2}\Big{(}F_{1,\theta}(s)(F_{1,\bar{g}}(s)-F_{1,\delta\bar{\epsilon}}(s))$
(212) $\displaystyle+F_{1,\delta}(s)F_{1,\bar{\pi}}(s)\Big{)}\,,$

where we used that vorticity and vector modes contribute only starting at one-
loop, and that $A\mapsto\delta$ in linear approximation. When including terms
up to a given loop order $L$, Eq. (210) is a polynomial equation for $x$ of
degree $L+1$, the solution(s) $x_{*}$ of which determine the ratio of scales

$k_{\sigma}/k_{nl}=\sqrt{3}x_{*}^{1/(n_{s}+3)}\,,$ (213)

setting the overall magnitude of the background dispersion
$\epsilon_{0}=1/k_{\sigma}^{2}$. In linear approximation, the solution is
given by

$x_{*}^{\text{lin}}=3(n_{s}+7)/((n_{s}+3)I^{\text{lin}}(n_{s}))\,.$ (214)

While significant changes from non-linear corrections can be expected, we
discuss self-consistent solutions in linear approximation as a proof-of-
principle in this work. Non-linear corrections are presented in paper II Garny
_et al._ (2022). In particular, we are interested in the sensitivity to the
truncation of the cumulant expansion, and start with the case of including
velocity dispersion only.

###  VIII.1 Self-consistent solution in second cumulant approximation

Figure 8:  Velocity dispersion scale $k_{\sigma}=\epsilon^{-1/2}$ relative to
the non-linear scale $k_{nl}$ for power-law initial spectrum $P_{0}\propto
k^{n_{s}}$, obtained when using the linear approximation for the source term
for $\epsilon$, and neglecting third and higher cumulants.

When neglecting third and higher cumulants, we can obtain the self-consistent
solution for the background dispersion by computing the integral Eq. (212)
using the analytical linear kernels for $F_{1,\theta}$ and $F_{1,\bar{g}}$
given in App. F, while $F_{1,\bar{\pi}}=F_{1,\delta\bar{\epsilon}}=0$. The
latter follows since the scalar mode $\delta\epsilon$ has no growing mode in
linear approximation, and at second cumulant order, and the former since the
third cumulant perturbation $\pi$ is neglected presently (see below for the
generalization to higher cumulants).

The self-consistent solution exists provided that the integral
$I^{\text{lin}}(n_{s})$ converges. In the infrared limit, $k\propto s\to 0$,
one has $F_{1,\theta}\to 1$ and $F_{1,\bar{g}}\to 2/(2+\alpha)$. Therefore the
integral is infrared-finite provided that $n_{s}>-3$. In the ultraviolet limit
$k\propto s\gg 1$, the linear kernels have asymptotic form given in Eq. (143).
The integral is absolutely convergent if

$d(n_{s})\equiv
3+n_{s}-\min(d_{\theta}+d_{g},d_{\theta}+e_{g},e_{\theta}+d_{g},e_{\theta}+e_{g})$
(215)

is less than zero. Using results from App. F one finds

$d(n_{s})=\left\\{\begin{array}[]{ccc}-(11n_{s}+53)/24&\mbox{for}&n_{s}<1\,,\\\
-(n_{s}+7)/3&\mbox{for}&n_{s}>1\,,\\\ \end{array}\right.$ (216)

which is less than zero for $n_{s}>-53/11\approx-4.8$. Note that even for a
very blue initial spectrum (large $n_{s}$), the damping due to velocity
dispersion is strong enough to compensate the growth of power at large $k$ and
make the integral converge. On the contrary, the sensitivity to short modes
grows when decreasing $n_{s}$, e.g. $d(-1)\approx-1.8$, $d(-2)\approx-1.3$,
$d(-3)\approx-0.8$. Overall, the whole integral is convergent for all
$n_{s}>-3$. Note that the solution is not valid for $n_{s}<-3$, because then
one would have $\alpha<0$, i.e. a velocity dispersion that was larger in the
past and decays with time, rather than being generated, leading to a
qualitatively different behavior. For the relevant case where $n_{s}\to-3$
from above, one has $\alpha\to\infty$, i.e. velocity dispersion grows very
quickly.

In Fig. 8, the ratio of $k_{\sigma}=\epsilon_{0}^{-1/2}$ and $k_{nl}$, as
determined by the linear approximation, is shown by the black line. The ratio
is always larger than one, and becomes very large for $n_{s}\to-3$. This
implies that non-linear corrections are expected to be more relevant the
smaller $n_{s}$.

###  VIII.2 Self-consistent solution in third and fourth cumulant
approximation

When including perturbation modes of the third and fourth cumulant, the
background value $\omega(\eta)$ of the fourth cumulant has to be taken into
account. Its equation of motion when neglecting fifth cumulant perturbations
is given in Eq. (V.2), and can (for $\Omega_{m}/f^{2}=1$) be rewritten as

$\displaystyle\frac{\omega^{\prime}+2\omega}{\epsilon^{2}}$ $\displaystyle=$
$\displaystyle\frac{1}{3}\int_{0}^{\infty}\frac{dk}{k}\,\Big{\\{}4\Delta_{\theta\bar{\xi}}-\Delta_{\theta\bar{\kappa}}+2\Delta_{\bar{g}\bar{\pi}}$
(217)
$\displaystyle+6\Delta_{\delta\bar{\epsilon}\bar{\pi}}-\frac{8}{5}\Delta_{\bar{g}\bar{\chi}}\Big{\\}}\,,$

where $\bar{\pi},\bar{\chi}$ and $\bar{\xi},\bar{\kappa}$ are the
dimensionless perturbation modes of the third and fourth cumulant,
respectively. Let us show that for a scaling universe, solutions with time-
dependence $\omega(\eta)\propto\epsilon(\eta)^{2}$ are consistent with this
equation of motion. For constant $\bar{\omega}=\omega/\epsilon^{2}$, the left-
hand side of Eq. (217) is equal to
$2(\alpha+1)\bar{\omega}=2(n_{s}+7)\bar{\omega}/(n_{s}+3)$, and is itself
time-independent. The dimensionless power spectra on the right-hand side can
be decomposed in a sum over loop contributions. Each of them satisfies a
relation analogous to Eq. (208), except that the
$\hat{\Delta}_{ab}^{L-\text{loop}}$ can depend on $\bar{\omega}$ in addition
to $s$. Therefore, also the right-hand side of Eq. (217) is constant after
integration over $k\propto s$, if $\bar{\omega}$ is constant. Thus, both the
left- and right-hand side of Eq. (217) are time-independent for constant
$\bar{\omega}$, implying that $\omega(\eta)\propto\epsilon(\eta)^{2}$ is a
consistent ansatz. What remains to be done is to find a solution for the
constant values of $k_{\sigma}/k_{nl}$ and $\bar{\omega}$, or equivalently $x$
and $\bar{\omega}$. For that purpose we rewrite Eq. (210) and Eq. (217) in the
form

$\displaystyle\frac{n_{s}+7}{n_{s}+3}$ $\displaystyle=$
$\displaystyle\frac{1}{3}\sum_{L\geq
0}x^{L+1}I^{L-\text{loop}}(n_{s},\bar{\omega})\,,$ $\displaystyle
2\frac{n_{s}+7}{n_{s}+3}\bar{\omega}$ $\displaystyle=$
$\displaystyle\frac{1}{3}\sum_{L\geq
0}x^{L+1}J^{L-\text{loop}}(n_{s},\bar{\omega})\,,$ (218)

yielding a coupled set of equations for the unknowns $x$ and $\bar{\omega}$,
with a polynomial dependence on $x$ for any given loop order, and an implicit
dependence on $\bar{\omega}$ that can in general only be determined
numerically. Here we defined

$\displaystyle J^{L-\text{loop}}(n_{s},\bar{\omega})$ $\displaystyle\equiv$
$\displaystyle\int_{0}^{\infty}\frac{ds}{s}\Big{(}4\hat{\Delta}_{\theta\bar{\xi}}(s)-\hat{\Delta}_{\theta\bar{\kappa}}(s)+2\hat{\Delta}_{\bar{g}\bar{\pi}}(s)$
(219)
$\displaystyle+6\hat{\Delta}_{\delta\bar{\epsilon}\bar{\pi}}(s)-\frac{8}{5}\hat{\Delta}_{\bar{g}\bar{\chi}}(s)\Big{)}^{L-\text{loop}}\,.$

The integrals $I^{L-\text{loop}}$ are given by the same expression as in Eq.
(211), but with power spectra computed based on kernels and non-linear
vertices including third and fourth cumulants. In turn, this leads to an
implicit dependence on the background value of the fourth cumulant
$\bar{\omega}$, as indicated in the arguments of $I^{L-\text{loop}}$ and
$J^{L-\text{loop}}$.

In linear approximation, we can eliminate $x$ by taking the ratio of both
equations in Eq. (VIII.2), giving an implicit equation for the fourth cumulant
expectation value $\bar{\omega}$,

$\bar{\omega}=\frac{1}{2}\frac{J^{\text{lin}}(n_{s},\bar{\omega})}{I^{\text{lin}}(n_{s},\bar{\omega})}\,.$
(220)

The solution $\bar{\omega}_{*}^{\text{lin}}$ of this equation can be
determined numerically using Eq. (212) for $I^{\text{lin}}$ as well as an
analogous expression for $J^{\text{lin}}\equiv J^{0-\text{loop}}$ and Eq.
(219). The power spectra entering both integrals can be expressed in terms of
the linear kernels $F_{1,a}$ for the dimensionless perturbation variables via
Eq. (206). The $F_{1,a}$ are obtained by numerically solving the linear
evolutions equations Eq. (VI.2) and using Eq. (151). Finally, the background
dispersion scale $k_{\sigma}/k_{nl}$ can be obtained using Eq. (213) with the
linear solution for $x$ given by Eq. (214) with
$I^{\text{lin}}=I^{\text{lin}}(n_{s},\bar{\omega}_{*}^{\text{lin}})$ evaluated
on the solution $\bar{\omega}=\bar{\omega}_{*}^{\text{lin}}$.

When truncating the perturbation modes at third cumulant order, corresponding
to (cum3+) or equivalently $c_{\text{max}}=3$, we find that Eq. (220) indeed
has a solution in linear approximation. The corresponding values of
$\bar{\omega}$ as well as $k_{\sigma}/k_{nl}$ are given in Table 1 for various
spectral indices $n_{s}$. We observe that the dimensionless fourth cumulant
expectation value is of order unity, indicating that higher cumulants are
relevant quantitatively, but of the same order as the background dispersion.
In addition, the shift in the value of $k_{\sigma}/k_{nl}$ compared to the
second cumulant approximation is sizeable, while the overall magnitude is
comparable. This indicates that higher cumulants are important quantitatively,
but do not invalidate the qualitative behavior of the second cumulant
approximation. We find that no self-consistent solutions exist when including
fourth cumulant perturbations, which may be attributed to the shortcomings of
the linear approximation. We investigate the impact of cumulants beyond the
fourth order in the next section.

Table 1: Self-consistent solutions within linear approximation for the velocity dispersion scale $k_{\sigma}=\epsilon^{-1/2}={\cal E}_{2}^{-1/2}$ relative to the non-linear scale, as well as the normalized expectation values $\bar{\cal E}_{2n}={\cal E}_{2n}/\epsilon^{n}$ of higher cumulants. We show results for scaling universes with spectral indices $n_{s}=-1,0,1,2$, and for various truncations of the cumulant expansion, with perturbation modes up to order $c_{\text{max}}$. The cases $c_{\text{max}}=2,3$ are equivalent to (cum2) and (cum3+), respectively, with $\bar{\cal E}_{4}=3\bar{\omega}/5$. | $n_{s}=-1$ | $n_{s}=0$ | $n_{s}=1$ | $n_{s}=2$
---|---|---|---|---
$c_{\text{max}}$ | $k_{\sigma}/k_{nl}$ | $\bar{\cal E}_{4}$ | $\bar{\cal E}_{6}$ | $\bar{\cal E}_{8}$ | $k_{\sigma}/k_{nl}$ | $\bar{\cal E}_{4}$ | $\bar{\cal E}_{6}$ | $\bar{\cal E}_{8}$ | $k_{\sigma}/k_{nl}$ | $\bar{\cal E}_{4}$ | $\bar{\cal E}_{6}$ | $\bar{\cal E}_{8}$ | $k_{\sigma}/k_{nl}$ | $\bar{\cal E}_{4}$ | $\bar{\cal E}_{6}$ | $\bar{\cal E}_{8}$
$2$ | 4.1 | - | - | - | 3.0 | - | - | - | 2.6 | - | - | - | 2.3 | - | - | -
$3$ | 3.4 | 0.45 | - | - | 2.5 | 0.40 | - | - | 2.2 | 0.37 | - | - | 2.0 | 0.35 | - | -
$6$ | 3.8 | 0.37 | 0.86 | - | 2.7 | 0.34 | 0.92 | - | 2.3 | 0.31 | 0.93 | - | 2.1 | 0.29 | 0.92 |
$7$ | 3.8 | 0.36 | 0.78 | 3.5 | 2.7 | 0.36 | 0.94 | 4.5 | 2.3 | 0.35 | 1.03 | 5.1 | 2.1 | 0.34 | 1.08 | 5.3

###  VIII.3 Self-consistent solutions for the full cumulant hierarchy

Figure 9:  Velocity dispersion scale $k_{\sigma}=\epsilon^{-1/2}$ relative to
the non-linear scale $k_{nl}$ for power-law initial spectrum $P_{0}\propto
k^{n_{s}}$, obtained from a self-consistent solution when including cumulant
perturbations up to order $c_{\text{max}}$, and for $n_{s}=-1,0,1,2$,
respectively. The corresponding self-consistent solutions for the fourth,
sixth and eighth cumulant expectation values are shown in Table 1.

The self-consistent solutions of background values for a scaling universe can
be extended to truncations including cumulants beyond the fourth order
following Sec. VII. In linear approximation, the expectation values ${\cal
E}_{2n}$ of the $2n$th cumulant satisfy equations of motion given in Eq.
(284). For $2n=2$ and $2n=4$ they agree with those discussed above for
$\epsilon={\cal E}_{2}$ and $\omega=5{\cal E}_{4}/3$, respectively. For a
scaling universe these equations allow for constant values of the
dimensionless ratios $\bar{\cal E}_{2n}={\cal E}_{2n}/\epsilon^{n}$,
determined by the set of implicit equations (in linear approximation)

$n\frac{n_{s}+7}{n_{s}+3}\bar{\cal E}_{2n}=x\times I_{\bar{\cal
E}_{2n}}^{\text{lin}}(n_{s},\bar{\cal E}_{4},\bar{\cal E}_{6},\dots)\,,$ (221)

with

$\displaystyle I_{\bar{\cal E}_{2n}}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}ds\,s^{n_{s}+2}\sum_{\ell=0}^{n}(\ell+1)\sum_{m_{1},m_{2}=0}^{n-\ell}\delta_{m_{1}+m_{2},n-\ell}$
$\displaystyle\frac{(2n)!(2(m_{2}-m_{1}-\ell)-3)}{(2m_{1})!(2m_{2})!}T_{\ell+1,2m_{1}}(s)T_{\ell,2m_{2}}(s)\,,$

where $\delta_{i,j}$ is the Kronecker symbol, and $T_{\ell,2m}(s)$ are the
dimensionless linear kernels for perturbation modes of cumulant order
$\ell+2m$ defined in Eq. (167). The linear kernels are given by numerical
solutions of the equations of motion Eq. (286). They depend on time and scale
via the single variable $s=\sqrt{3\epsilon(\eta)k^{2}}$, which can be seen
using Eq. (172) and that $\alpha=\partial_{\eta}\ln\epsilon$ is constant for a
scaling universe. In addition, the linear kernels depend parametrically on the
spectral index $n_{s}$ due to $\alpha=4/(3+n_{s})$, as well as the backgound
values $\bar{\cal E}_{4},\bar{\cal E}_{6},\dots$. Therefore, Eq. (221) is a
highly implicit and coupled set of equations for the self-consistent values of
$\bar{\cal E}_{2n}$ with $2n\geq 4$ as well as the overall magnitude of the
background dispersion parameterized by the variable $x$ defined in Eq. (205).
By taking the ratio of Eq. (221) for $2n\geq 4$ and Eq. (221) for $2n=2$, and
using that $\bar{\cal E}_{2}=1$ by definition, we obtain a coupled set of
equations for the cumulant expectation values of order $2n\geq 4$ that is
independent of $x$,

$\bar{\cal E}_{2n}=\frac{1}{n}\frac{I_{\bar{\cal
E}_{2n}}^{\text{lin}}(n_{s},\bar{\cal E}_{4},\bar{\cal
E}_{6},\dots)}{I_{\bar{\cal E}_{2}}^{\text{lin}}(n_{s},\bar{\cal
E}_{4},\bar{\cal E}_{6},\dots)}\,.$ (223)

These implicit equations can be viewed as a generalization of Eq. (220) to
beyond the fourth cumulant order. Once a solution to these equations is found,
it can be inserted in Eq. (221) for $n=1$, yielding the solution for $x$ and
thereby for the background dispersion scale

$\frac{k_{\sigma}}{k_{nl}}=\sqrt{3}\left(\frac{n_{s}+7}{(n_{s}+3)I_{\bar{\cal
E}_{2}}^{\text{lin}}}\right)^{1/(n_{s}+3)}\,.$ (224)

For a given truncation of cumulant perturbations at order $c_{\text{max}}$,
all $T_{\ell,2m}(s)$ with $\ell+2m>c_{\text{max}}$ are neglected. In this case
the equations for the linear kernels and therefore also $I_{\bar{\cal
E}_{2n}}^{\text{lin}}$ depend only on the background values $\bar{\cal
E}_{4},\dots,\bar{\cal E}_{c_{\text{max}}+1}$ for odd $c_{\text{max}}$, and
$\bar{\cal E}_{4},\dots,\bar{\cal E}_{c_{\text{max}}}$ for even
$c_{\text{max}}$. For $c_{\text{max}}=2$ only the background dispersion enters
and it is sufficient to solve Eq. (224), with results identical to those from
Sec. VIII.1. For $c_{\text{max}}=3,4$, Eq. (223) reduces to a single equation
and yields values for $\bar{\cal E}_{4}=3\bar{\omega}/5$ consistent with those
obtained in Sec. VIII.2. For $c_{\text{max}}=5,6$, Eq. (223) yields a coupled
set of equations for $\bar{\cal E}_{4},\bar{\cal E}_{6}$. For
$c_{\text{max}}=7,8$, one obtains three coupled equations for $\bar{\cal
E}_{4},\bar{\cal E}_{6},\bar{\cal E}_{8}$. Up to eighth order, we find that
within the linear approximation considered here a joint self-consistent
solution exists only for $c_{\text{max}}=2,3,6,7$, confirming the previous
findings for $c_{\text{max}}=2,3$. The corresponding values are shown in Table
1 for various values of $n_{s}$. We observe that the size of $\bar{\cal
E}_{4}$ is comparable for $c_{\text{max}}=3,6,7$, and the one of $\bar{\cal
E}_{6}$ for $c_{\text{max}}=6,7$, indicating that increasing the truncation
order does not lead to dramatic changes.

The dependence of the dispersion scale $k_{\sigma}$ on the truncation order
$c_{\text{max}}$ is shown in Fig. 9. The largest shift occurs when going from
$c_{\text{max}}=2$ to $c_{\text{max}}=3$, while even higher cumulants have
only a minor impact. This indicates that the relevant contributions to the
source terms Eq. (VIII.3) arise from scales where the dependence of the linear
kernels on $c_{\text{max}}$ is already converged. In addition, we find that
the impact of higher cumulants is smaller for larger $n_{s}$. While spectra
with large values of $n_{s}$ do have a lot of power on small scales initially,
it is more efficiently erased by the suppression due to the buildup of
velocity dispersion, making the linear kernels drop faster, and hence the
source terms less sensitive to the contribution from large wavenumbers.

##  IX Conclusions

In this work we discuss the extension of standard perturbation theory (SPT) to
include higher cumulants of the phase-space distribution function, based on
the underlying Vlasov-Poisson dynamics for collisionless matter, dubbed
_Vlasov Perturbation Theory_ (). This takes into account that even for an
initially perfectly cold dark matter distribution, orbit crossing generates
velocity dispersion and higher cumulants. We provide the explicit form of non-
linear evolution equations when taking up to the fourth cumulant into account,
and derive evolution equations linear in perturbations up to arbitrary order
in the cumulant expansion.

splits cumulants into their average values and fluctuations around them. The
evolution equations for the perturbations depend on the average values, and
vice versa. Since in general the average values are sourced by fluctuations
integrated over all scales, we argue that the background values of even
cumulants should be treated as “${\cal O}(1)$” quantities. This leads to a
consistent, systematic perturbative expansion scheme for the fluctuations, and
allows us to describe the screening of UV modes crucial for improving the
convergence of SPT, as we discuss in paper II Garny  _et al._ (2022). The
resulting equations can be cast into a form that is formally analogous to SPT,
but with an extended set of perturbation variables and non-linear interactions
among them, in presence of a background given by the average values of the
even cumulants.

Our main findings are:

  * i)

Linear is far richer compared to SPT. Even in the simplest approximation,
where linear theory is truncated at the second cumulant, the effective
description in terms of density and velocity divergence is non-local in time,
which is key to satisfy the cosmic energy equation. When a given mode $k$
crosses the dispersion scale $k_{\sigma}$, its growth is suppressed. This back
reaction on modes from small-scale dispersion is completely absent in SPT. The
suppression is only affected by higher cumulants when $k\gg k_{\sigma}$, and
they make the screening mechanism even more efficient. For any given
wavenumber the cumulant expansion converges. This motivates a further study of
perturbation theory with a truncated cumulant expansion beyond the linear
approximation. This is what we carry out in paper II Garny  _et al._ (2022).

  * ii)

The UV screening mechanism is in principle not guaranteed, i.e. the complexity
of the linear theory of collisionless dynamics allows exponential
instabilities. Requiring that these be absent leads to stability conditions
that we derive analytically up to cumulant of order eight. Remarkably, these
stability conditions are independent of the value of the velocity dispersion
and spectral index and only constrain the non-Gaussianity of the distribution
function. The Gaussian case (vanishing average values of fourth and higher
cumulants) is always within the stable domain.

  * iii)

We therefore consider higher cumulants averaged over stationary dark matter
haloes, finding that they are generically of order unity when compared to an
appropriate power of the second cumulant. While the halo analysis is not
required for the development of the framework, it serves as a useful benchmark
and to gain some insight into the distribution function non-Gaussianity.
Interestingly, we find that the cumulant expectation values obtained from the
halo analysis satisfy the stability conditions. Altogether, this implies that
the linear approximation within is a good starting point for a perturbative
analysis with realistic values of the cumulant expectation values.

  * iv)

Finally, we determine self-consistent solutions of the coupled set of
perturbation and background equations. We consider a scaling universe for this
analysis, allowing us to transform the set of equations for the average values
into coupled algebraic equations, that we solve up to cumulant of order eight.
We find that the decoupling of UV modes is more pronounced for larger spectral
index $n_{s}$. This is because very blue initial spectra lead to pronounced
orbit crossing on small scales and therefore quickly generate a large
dispersion. Remarkably, this causes the integral over the power spectra of
cumulant perturbations that source the average values to converge even for
arbitrarily large $n_{s}$. The resulting background values for the cumulants
in this self-consistent approach also satisfy the stability conditions.

In summary, the framework of perturbation theory for dark matter clustering
laid out in this work is directly derived from the underlying fundamental
collisionless Vlasov-Poisson equations. It can, from the conceptual point of
view, be regarded as a straightforward extension of SPT by taking second and
higher cumulants of the distribution function into account. Our results show
that captures physical effects that are neglected in SPT, in particular shell
crossing and the screening of UV modes. It therefore abandons a major
shortcoming of SPT. Furthermore, the framework does by construction neither
contain any ad hoc assumptions nor undetermined free fitting parameters, and
therefore preserves the predictivity of the underlying Vlasov-Poisson
equations. Nevertheless, it allows for systematic and tractable extensions of
SPT. In an accompanying work, we show that the extended framework can be used
beyond the linear approximation as well, and present detailed comparisons with
N-body simulation results (see paper II Garny  _et al._ (2022)).

###### Acknowledgements.

We thank Michael Buehlmann, ChangHoon Hahn, Oliver Hahn and Lam Hui for useful
discussions. MG and RS acknowledge support from the Munich Institute for
Astro- and Particle Physics (MIAPP), and RS thanks the Theory Group at CERN
for hospitality during the Fall of 2014, where this collaboration was started.
We acknowledge support by the Excellence Cluster ORIGINS, which is funded by
the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under
Germany’s Excellence Strategy - EXC-2094 - 390783311.

##  APPENDIX A Further Details on Halo Calculations

In this appendix we provide some extensions of the calculations presented in
Sec. II for the cumulants of the phase-space distribution function and an
alternative approach in the NFW case that gives the shape of the distribution
function itself.

###  A.1 Evans Halos

We first we provide some more details on the halo calculations presented in
Sec. II.1. For the distribution function given by Eq. (4), the expectation
values of the sixth and eighth cumulants are, respectively

$\displaystyle 7\,{\cal E}_{6}$ $\displaystyle=$ $\displaystyle{315\over
8}\mathop{\langle}w_{a}\mathop{\rangle}+{105\over
8}\mathop{\langle}w_{b}\mathop{\rangle}+105\mathop{\langle}w_{c}\mathop{\rangle}+\frac{687}{4}\mathop{\langle}w_{a}^{3}\mathop{\rangle}-\frac{1449}{8}\mathop{\langle}w_{a}^{2}\mathop{\rangle}-\frac{315}{2}\mathop{\langle}w_{a}w_{b}\mathop{\rangle}+\frac{987}{4}\mathop{\langle}w_{a}^{2}w_{b}\mathop{\rangle}-\frac{315}{8}\mathop{\langle}w_{b}^{2}\mathop{\rangle}$
(225)
$\displaystyle+\frac{525}{4}\mathop{\langle}w_{a}w_{b}^{2}\mathop{\rangle}+\frac{105}{4}\mathop{\langle}w_{b}^{3}\mathop{\rangle}-\frac{1785}{4}\mathop{\langle}w_{a}w_{c}\mathop{\rangle}+\frac{987}{2}\mathop{\langle}w_{a}^{2}w_{c}\mathop{\rangle}-\frac{945}{4}\mathop{\langle}w_{b}w_{c}\mathop{\rangle}+525\mathop{\langle}w_{a}w_{b}w_{c}\mathop{\rangle}+\frac{315}{2}\mathop{\langle}w_{b}^{2}w_{c}\mathop{\rangle}$
$\displaystyle-315\mathop{\langle}w_{c}^{2}\mathop{\rangle}+525\mathop{\langle}w_{a}w_{c}^{2}\mathop{\rangle}+315\mathop{\langle}w_{b}w_{c}^{2}\mathop{\rangle}+210\mathop{\langle}w_{c}^{3}\mathop{\rangle}\,,$
$\displaystyle 9\,{\cal E}_{8}$ $\displaystyle=$
$\displaystyle\frac{3465}{16}\mathop{\langle}w_{a}\mathop{\rangle}+\frac{945}{16}\mathop{\langle}w_{b}\mathop{\rangle}+945\mathop{\langle}w_{c}\mathop{\rangle}-\frac{42399}{16}\mathop{\langle}w_{a}^{2}\mathop{\rangle}+\frac{28215}{4}\mathop{\langle}w_{a}^{3}\mathop{\rangle}-\frac{42003}{8}\mathop{\langle}w_{a}^{4}\mathop{\rangle}-\frac{15435}{8}\mathop{\langle}w_{a}w_{b}\mathop{\rangle}+\frac{34965}{4}\mathop{\langle}w_{a}^{2}w_{b}\mathop{\rangle}$
(226)
$\displaystyle-\frac{18549}{2}\mathop{\langle}w_{a}^{3}w_{b}\mathop{\rangle}-\frac{6615}{16}\mathop{\langle}w_{b}^{2}\mathop{\rangle}+\frac{16065}{4}\mathop{\langle}w_{a}w_{b}^{2}\mathop{\rangle}-\frac{26649}{4}\mathop{\langle}w_{a}^{2}w_{b}^{2}\mathop{\rangle}+\frac{2835}{4}\mathop{\langle}w_{b}^{3}\mathop{\rangle}-\frac{4725}{2}\mathop{\langle}w_{a}w_{b}^{3}\mathop{\rangle}-\frac{2835}{8}\mathop{\langle}w_{b}^{4}\mathop{\rangle}$
$\displaystyle-7875\mathop{\langle}w_{a}w_{c}\mathop{\rangle}+21924\mathop{\langle}w_{a}^{2}w_{c}\mathop{\rangle}-18549\mathop{\langle}w_{a}^{3}w_{c}\mathop{\rangle}-3780\mathop{\langle}w_{b}w_{c}\mathop{\rangle}+20790\mathop{\langle}w_{a}w_{b}w_{c}\mathop{\rangle}-26649\mathop{\langle}w_{a}^{2}w_{b}w_{c}\mathop{\rangle}$
$\displaystyle+5670\mathop{\langle}w_{b}^{2}w_{c}\mathop{\rangle}-14175\mathop{\langle}w_{a}w_{b}^{2}w_{c}\mathop{\rangle}-2835\mathop{\langle}w_{b}^{3}w_{c}\mathop{\rangle}-6615\mathop{\langle}w_{c}^{2}\mathop{\rangle}+25515\mathop{\langle}w_{a}w_{c}^{2}\mathop{\rangle}-26649\mathop{\langle}w_{a}^{2}w_{c}^{2}\mathop{\rangle}+14175\mathop{\langle}w_{b}w_{c}^{2}\mathop{\rangle}$
$\displaystyle-28350\mathop{\langle}w_{a}w_{b}w_{c}^{2}\mathop{\rangle}-8505\mathop{\langle}w_{b}^{2}w_{c}^{2}\mathop{\rangle}+11340\mathop{\langle}w_{c}^{3}\mathop{\rangle}-18900\mathop{\langle}w_{a}w_{c}^{3}\mathop{\rangle}-11340\mathop{\langle}w_{b}w_{c}^{3}\mathop{\rangle}-5670\mathop{\langle}w_{c}^{4}\mathop{\rangle}\,.$

As mentioned in Sec. II, for infinite halos $w_{b}$ drops from these
expressions and everything can be written down in terms of $w_{a}$ due to the
$w_{a}+w_{b}+w_{c}=1$ identity. The integrals over the halo of $w_{a}^{n}$ can
be done analytically, resulting in simple analytic expressions for the
expectation value of the cumulants. Defining a real-valued function $f(q)$ for
both prolate ($q>1$) and oblate ($q<1$) halos,

$f(q>1)\equiv\frac{1}{\sqrt{q^{2}-1}}\
\operatorname{arc\,cot}\Big{(}\frac{q}{\sqrt{q^{2}-1}}\Big{)},\ \ \ \ \
f(q<1)\equiv\frac{1}{\sqrt{1-q^{2}}}\
\operatorname{arc\,coth}\Big{(}\frac{q}{\sqrt{1-q^{2}}}\Big{)}\,,$ (227)

we have for the expectation value of the cumulants

${\cal E}_{2}=1-{1\over 3}q^{2}+{q\over 3}(2q^{2}-1)\,f(q),\ \ \ \ \ {\cal
E}_{4}={1\over 5}(1-14q^{2})+{q\over 5}\,(2q^{2}-1)(7q^{2}+6)\,f(q)\,,$ (228)
${\cal E}_{6}=-{3q^{2}\over 28}\,(1-234q^{2}+144q^{4})+{3q\over
28}\,(2q^{2}-1)(108q^{4}-234q^{2}+1)\,f(q)\,,$ (229) ${\cal E}_{8}={q^{2}\over
24}\,(3183-12308q^{2}+5292^{4}-13696q^{6})+{q\over
8}\,(2q^{2}-1)(2568q^{6}+924q^{4}+940q^{2}+1411)\,f(q)\,.$ (230)

These are the expressions used to compute the normalized cumulants $\bar{\cal
E}_{2n}={\cal E}_{2n}/{\cal E}_{2}^{n}$ shown in Fig. 1.

Figure 10:  The distribution function monopole for NFW halos with $\beta=1/2$
(solid) compared to a Gaussian distribution with the same second cumulant
(dashed), at the scale radius of the halo $x\equiv cr/r_{\rm vir}=1$. The left
panel shows $f^{(0)}(p)$ as a function of momentum $p$ in units of
$p_{0}\equiv\sqrt{Gmf(c)c/r_{\rm vir}}$ while the right panel shows $4\pi
p^{2}f^{(0)}(p)$.

###  A.2 NFW Halos

In Sec. II.2 we calculated the cumulants of the phase-space distribution
function for a distribution function with constant anisotropy, Eq. (22), by
integrating the steady Vlasov equation directly, leading to a recursion
relation for the moments in terms of the density profile and potential, Eq.
(LABEL:recrel2). Here we take an alternative approach, which is to compute
directly the shape of the distribution function. The function $f_{E}$ in Eq.
(22) can be computed by a generalization of the so-called Eddington inversion
method Eddington (1916); Binney and Tremaine (2008) originally developed for
the isotropic case ($\beta=0$),

$f_{E}(\varepsilon)={2^{\beta}\over(2\pi)^{3/2}\Gamma(1-\lambda)\Gamma(1-\beta)}\,{d\over
d\varepsilon}\int_{0}^{\varepsilon}\frac{d\psi}{(\varepsilon-\psi)^{\lambda}}\,{d^{n}h\over
d\psi^{n}}\,,$ (231)

where $\varepsilon\equiv-E=\psi-p^{2}/2\geq 0$ is the binding energy per unit
mass, $\psi=-\Phi$, $\lambda=3/2-\beta-n$, $n=[3/2-\beta]$ and $h\equiv
r^{2\beta}\rho$ is the radially weighted density profile expressed as a
function of $\psi$ rather than $r$. When $\beta$ is half-integer
($\beta=1/2,-1/2,\ldots$) this reduces to only derivatives, giving a simple
expression for the distribution function Evans and Williams (2014),

$f(E,L^{2})={1\over 2\pi^{2}}\,{L^{-2\beta}\over(-2\beta)!!}\
{d^{3/2-\beta}h\over d\psi^{3/2-\beta}}\Big{|}_{\psi=\varepsilon}\,.$ (232)

In particular, we are interested in $\beta=1/2$ in which case this expression
reduces to a first derivative. Using these results, one could compute the
cumulants by integration over momentum rather than using the recursion
relation Eq. (LABEL:recrel2), but in practice the latter is actually easier.
However, this approach also gives us explicitly the distribution function. For
$\beta=1/2$ we have,

$f(\bm{p},\bm{x})={1\over 2\pi^{2}}\,{1\over
p\,r\sqrt{1-(\hat{\bm{p}}\cdot\hat{\bm{x}})^{2}}}\ {dh\over
d\psi}\Big{|}_{\psi\to\psi-p^{2}/2},\ \ \ \ \ \ \ \ \ \ f^{(0)}(p,r)={1\over
p\,r}\ {dh\over d\psi}\Big{|}_{\psi\to\psi-p^{2}/2}\,,$ (233)

where $r\equiv|\bm{x}|$. Taking the monopole of this distribution by
integrating over the direction of momentum gives $f^{(0)}$. Figure 10 shows
the result of this calculation for $\beta=1/2$ (solid lines), comparing it to
a Gaussian distribution of the same width (dashed). As expected from the form
of the distribution function, and the normalized cumulants shown in the left
panel of Fig. 2, the distribution is fairly different from a Maxwellian.

##  APPENDIX B Equations of motion up to the second cumulant

In this appendix we give the equations of motion up to the second cumulant
when neglecting third and higher order cumulants in non-linear terms (see App.
D for those). The equations for the density contrast $\delta$, the velocity
divergence $\theta$ and the log-density field $A$ are given in the main text
in Eqs. (69, IV.2, 72), respectively. We use non-bold symbols to denote
wavevectors here and below.

###  B.1 Equation of motion for the vorticity

$\displaystyle
w_{k,i}^{\prime}+\left(\frac{3}{2}\frac{\Omega_{m}}{f^{2}}-1\right)w_{k,i}+k^{2}\nu_{k,i}$
(234) $\displaystyle=$
$\displaystyle\int_{pq}\Bigg{\\{}-\frac{1}{p^{2}}\theta_{p}(k\times(p\times
w_{q}))_{i}+\frac{1}{p^{2}}(k\times((p\times w_{p})\times
w_{q}))_{i}-\frac{p\cdot q}{q^{2}}(p\times q)_{i}A_{p}g_{q}+(p\times
q)_{i}A_{p}\delta\epsilon_{q}$ $\displaystyle{}-\frac{1}{q^{2}}\Big{[}(p\cdot
q)(k\cdot q\delta_{ij}-k_{j}q_{i})-(p\times q)_{i}(p\times
q)_{j}\Big{]}A_{p}\nu_{q,j}-\varepsilon_{ijn}k_{j}p_{m}A_{p}t_{q,nm}\Bigg{\\}}\,.$

###  B.2 Equations of motion for scalar, vector and tensor perturbations of
$\epsilon_{ij}$

$\displaystyle
g_{k}^{\prime}+2\left(\frac{3}{2}\frac{\Omega_{m}}{f^{2}}-1\right)g_{k}-2\epsilon\theta_{k}$
(235) $\displaystyle=$ $\displaystyle\int_{pq}\Bigg{\\{}\left(3\frac{(k\cdot
p)^{2}}{k^{2}p^{2}}-1\right)\theta_{p}\delta\epsilon_{q}{}+\frac{p\cdot
q}{q^{2}}\left(\frac{k^{2}}{p^{2}}+\frac{1}{2}-\frac{3}{2}\frac{(k\cdot
p)^{2}}{k^{2}p^{2}}\right)\theta_{p}g_{q}+3\frac{k\cdot p}{p^{2}k^{2}}(p\times
q)\cdot w_{p}\delta\epsilon_{q}$
$\displaystyle{}+\frac{1}{p^{2}}\left(\frac{3}{2}\frac{(k\cdot
q)^{2}}{k^{2}q^{2}}+3\frac{k\cdot qp\cdot q}{k^{2}q^{2}}-\frac{p\cdot
q}{q^{2}}-\frac{1}{2}\right)(p\times q)\cdot
w_{p}g_{q}{}-\frac{1}{k^{2}p^{2}q^{2}}\left(p\cdot qk^{2}+3k\cdot pk\cdot
q\right)\theta_{p}(p\times q)\cdot\nu_{q}$ $\displaystyle{}-\frac{p\cdot
q}{p^{2}q^{2}}(p\cdot q\delta_{ij}-q_{i}p_{j})w_{p,i}\nu_{q,j}{}-\frac{6k\cdot
q-k^{2}+3p\cdot q}{k^{2}p^{2}q^{2}}(p\times q)\cdot w_{p}(p\times
q)\cdot\nu_{q}$ $\displaystyle{}-\frac{2k^{2}-3p\cdot q-6k\cdot
p}{2k^{2}p^{2}}p_{i}p_{j}\theta_{p}t_{q,ij}+\frac{2k^{2}(p\times
w_{p})_{i}+9(p\times q)\cdot
w_{p}p_{i}}{2k^{2}p^{2}}p_{j}t_{q,ij}\Bigg{\\}}\,,$
$\displaystyle\delta\epsilon_{k}^{\prime}+2\left(\frac{3}{2}\frac{\Omega_{m}}{f^{2}}-1\right)\delta\epsilon_{k}$
(236) $\displaystyle=$ $\displaystyle\int_{pq}\Bigg{\\{}\frac{k\cdot pk\cdot
q}{k^{2}p^{2}}\theta_{p}\delta\epsilon_{q}{}+\frac{p\cdot
q}{2p^{2}q^{2}k^{2}}((k\cdot p)^{2}-p^{2}k^{2})\theta_{p}g_{q}+\frac{k\cdot
q}{p^{2}k^{2}}(p\times q)\cdot w_{p}\delta\epsilon_{q}$
$\displaystyle{}+\frac{k^{2}q^{2}-(k\cdot q)^{2}+2p\cdot qk\cdot
p}{2p^{2}q^{2}k^{2}}(p\times q)\cdot
w_{p}g_{q}{}+\frac{1}{q^{2}}\left(1-\frac{(k\cdot
p)^{2}}{k^{2}p^{2}}\right)\theta_{p}(p\times q)\cdot\nu_{q}$
$\displaystyle{}+\frac{p\cdot q}{p^{2}q^{2}}(p\cdot
q\delta_{ij}-q_{i}p_{j})w_{p,i}\nu_{q,j}{}+\frac{2k\cdot q-k^{2}+p\cdot
q}{k^{2}p^{2}q^{2}}(p\times q)\cdot w_{p}(p\times q)\cdot\nu_{q}$
$\displaystyle{}+\frac{2k\cdot q-p\cdot
q}{2k^{2}p^{2}}p_{i}p_{j}\theta_{p}t_{q,ij}-\frac{2k^{2}(p\times
w_{p})_{i}+3(p\times q)\cdot
w_{p}p_{i}}{2k^{2}p^{2}}p_{j}t_{q,ij}\Bigg{\\}}-Q\delta^{(3)}(k)\,,$
$\displaystyle\nu_{k,i}^{\prime}+2\left(\frac{3}{2}\frac{\Omega_{m}}{f^{2}}-1\right)\nu_{k,i}-\epsilon
w_{k,i}$ (237) $\displaystyle=$
$\displaystyle\int_{pq}\Bigg{\\{}-2\frac{k\cdot p(p\times
q)_{i}}{k^{2}p^{2}}\theta_{p}\delta\epsilon_{q}{}+\frac{p\cdot qk\cdot
p(p\times q)_{i}}{k^{2}p^{2}q^{2}}\theta_{p}g_{q}$
$\displaystyle{}+\frac{k\cdot p(k\cdot p\delta_{ij}-p_{i}k_{j})-(p\times
q)_{i}(p\times q)_{j}}{k^{2}p^{2}}w_{p,j}\delta\epsilon_{q}$
$\displaystyle{}+\frac{k\cdot qp\cdot q(k\cdot
p\delta_{ij}-p_{i}k_{j})+(p\times q)_{i}(p\times
q)_{j}(k^{2}-p^{2})}{k^{2}p^{2}q^{2}}w_{p,j}g_{q}$
$\displaystyle{}+\frac{k^{2}p\cdot q(k\cdot q\delta_{ij}-q_{i}k_{j})-(p\times
q)_{i}(p\times q)_{j}(p^{2}-q^{2})}{k^{2}p^{2}q^{2}}\theta_{p}\nu_{q,j}$
$\displaystyle{}+\frac{k^{2}-p^{2}}{k^{2}p^{2}q^{2}}\left[(k\times(p\times
w_{p}))_{i}(p\times q)\cdot\nu_{q}-(p\times q)\cdot
w_{p}(k\times(q\times\nu_{q}))_{i}\right]$ $\displaystyle{}-2\frac{(p\times
q)_{i}}{k^{2}p^{2}q^{2}}(p\times q)\cdot w_{p}(p\times
q)\cdot\nu_{q}+\frac{(k^{2}-q^{2})\epsilon_{inm}k_{n}-(p\times
q)_{i}p_{m}}{k^{2}p^{2}}p_{j}\theta_{p}t_{q,mj}$
$\displaystyle{}+\frac{2(p\times q)\cdot
w_{p}\epsilon_{inm}k_{n}-(k\times(p\times
w_{p}))_{i}p_{m}}{k^{2}p^{2}}p_{j}t_{q,mj}\Bigg{\\}}\,,$ $\displaystyle
t_{k,ij}^{\prime}+2\left(\frac{3}{2}\frac{\Omega_{m}}{f^{2}}-1\right)t_{k,ij}$
(238) $\displaystyle=$
$\displaystyle\int_{pq}\Bigg{\\{}\Bigg{[}\left(\delta_{ij}-\frac{k_{i}k_{j}}{k^{2}}\right)-2\frac{(p\times
q)_{i}(p\times q)_{j}}{(p\times q)^{2}}\Bigg{]}\frac{(p\times
q)^{2}}{k^{2}p^{2}}\Bigg{[}\theta_{p}\delta\epsilon_{q}-\frac{1}{2}\frac{p\cdot
q}{q^{2}}\theta_{p}g_{q}+\theta_{p}\frac{(p\times
q)\cdot\nu_{q}}{q^{2}}\Bigg{]}$
$\displaystyle{}+\frac{1}{k^{2}p^{2}}\Bigg{[}(q^{2}k_{i}-k^{2}q_{i})(p_{m}\delta_{j\ell}+p_{\ell}\delta_{jm})+(q^{2}k_{j}-k^{2}q_{j})(p_{m}\delta_{i\ell}+p_{\ell}\delta_{im})$
$\displaystyle{}+p\cdot
q\left(\delta_{i\ell}\delta_{jm}+\delta_{im}\delta_{j\ell}-\frac{p_{\ell}p_{m}}{k^{2}}\left(\delta_{ij}-\frac{k_{i}k_{j}}{k^{2}}\right)\right)$
$\displaystyle{}-2\frac{p_{\ell}p_{m}}{k^{2}}\left(k_{i}p_{j}+k_{j}p_{i}-k_{i}k_{j}+q^{2}\delta_{ij}-k\cdot
p\frac{k_{i}k_{j}}{k^{2}}\right)\Bigg{]}\theta_{p}t_{q,\ell m}$
$\displaystyle{}+\frac{1}{k^{2}p^{2}}\Bigg{[}(p\times
w_{p})_{i}\left(k_{j}(k\cdot p)-k^{2}p_{j})\right)+(p\times
w_{p})_{j}\left(k_{i}(k\cdot p)-k^{2}p_{i})\right)$ $\displaystyle{}+(p\times
q)\cdot w_{p}\left(k\cdot
p\left(\delta_{ij}+\frac{k_{i}k_{j}}{k^{2}}\right)-k_{i}p_{j}-k_{j}p_{i}\right)\Bigg{]}\delta\epsilon_{q}$
$\displaystyle{}+\frac{1}{p^{2}q^{2}}\Bigg{[}(p\times w_{p})_{i}(p\cdot
q)\left(k\cdot q\frac{k_{j}}{k^{2}}-q_{j}\right)+(p\times w_{p})_{j}(p\cdot
q)\left(k\cdot q\frac{k_{i}}{k^{2}}-q_{i}\right)$ $\displaystyle{}+(p\times
q)\cdot
w_{p}\Bigg{\\{}\left(\delta_{ij}+\frac{k_{i}k_{j}}{k^{2}}\right)\frac{k\cdot
q(k\cdot q+2p\cdot q)}{2k^{2}}+q_{i}q_{j}$
$\displaystyle{}+\frac{1}{2}(k^{2}-p^{2})\left(\delta_{ij}-\frac{k_{i}k_{j}}{k^{2}}+2\frac{k_{i}q_{j}+k_{j}q_{i}}{k^{2}}\right)\Bigg{\\}}\Bigg{]}g_{q}$
$\displaystyle{}+\frac{p^{2}k_{i}-k^{2}p_{i}}{k^{2}p^{2}q^{2}}\left[(p\times
q)\cdot\nu_{q}(p\times w_{p})_{j}-(p\times q)\cdot
w_{p}(q\times\nu_{q})_{j}\right]$
$\displaystyle{}+\frac{p^{2}k_{j}-k^{2}p_{j}}{k^{2}p^{2}q^{2}}\left[(p\times
q)\cdot\nu_{q}(p\times w_{p})_{i}-(p\times q)\cdot
w_{p}(q\times\nu_{q})_{i}\right]$ $\displaystyle{}+\frac{p\cdot
q}{p^{2}q^{2}}\Bigg{[}(p\times w_{p})_{i}(q\times\nu_{q})_{j}+(p\times
w_{p})_{j}(q\times\nu_{q})_{i}-(p\times
w_{p})\cdot(q\times\nu_{q})\left(\delta_{ij}-\frac{k_{i}k_{j}}{k^{2}}\right)\Bigg{]}$
$\displaystyle{}-\frac{(p\times q)\cdot w_{p}(p\times
q)\cdot\nu_{q}}{p^{2}q^{2}}\Bigg{[}\left(\delta_{ij}+\frac{k_{i}k_{j}}{k^{2}}\right)\frac{3p\cdot
q+2q^{2}}{k^{2}}-2\frac{k_{i}q_{j}+k_{j}q_{i}}{k^{2}}-\left(\delta_{ij}-\frac{k_{i}k_{j}}{k^{2}}\right)\Bigg{]}$
$\displaystyle{}+\frac{1}{2p^{2}}\Bigg{[}2\left(\delta_{ij}-\frac{k_{i}k_{j}}{k^{2}}\right)\left[(p\times
w_{p})_{m}p_{n}+(p\times w_{p})_{n}p_{m}\right]$ $\displaystyle{}-2(p\times
w_{p})_{i}\left(p_{n}\delta_{jm}+p_{m}\delta_{jn}-2\frac{p_{m}p_{n}k_{j}}{k^{2}}\right)-2(p\times
w_{p})_{j}\left(p_{n}\delta_{im}+p_{m}\delta_{in}-2\frac{p_{m}p_{n}k_{i}}{k^{2}}\right)$
$\displaystyle{}+(p\times q)\cdot
w_{p}\Bigg{(}6\frac{p_{m}p_{n}}{k^{2}}\left(\delta_{ij}+\frac{k_{i}k_{j}}{k^{2}}\right)+\delta_{im}\left(\delta_{jn}-4\frac{p_{n}k_{j}}{k^{2}}\right)$
$\displaystyle{}+\delta_{jn}\left(\delta_{ik}-4\frac{p_{m}k_{i}}{k^{2}}\right)+\delta_{in}\left(\delta_{jm}-4\frac{p_{m}k_{j}}{k^{2}}\right)+\delta_{jm}\left(\delta_{in}-4\frac{p_{n}k_{i}}{k^{2}}\right)\Bigg{)}\Bigg{]}t_{q,mn}\Bigg{\\}}\,.$

##  APPENDIX C Vertices up to the second cumulant

In this appendix we collect all vertices $\gamma_{abc}(p,q)$ for perturbation
modes up to the second cumulant. We write $\epsilon$ instead of
$\delta\epsilon$ in the index for simplicity. Furthermore $k\equiv p+q$ in the
vectorial sense.

###  C.1 Vertices involving only scalar perturbations

$\displaystyle\gamma_{\delta\theta\delta}(p,q)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\alpha_{pq}=\frac{(p+q)\cdot p}{2p^{2}}\,,$
$\displaystyle\gamma_{\theta\theta\theta}(p,q)$ $\displaystyle=$
$\displaystyle\beta_{pq}=\frac{(p+q)^{2}p\cdot q}{2p^{2}q^{2}}\,,$
$\displaystyle\gamma_{A\theta A}(p,q)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\frac{q\cdot p}{p^{2}}\,,$
$\displaystyle\gamma_{\theta Ag}(p,q)$ $\displaystyle=$
$\displaystyle-\frac{1}{2}(p+q)\cdot q\ \frac{q\cdot p}{q^{2}}\,,$
$\displaystyle\gamma_{\theta A\epsilon}(p,q)$ $\displaystyle=$
$\displaystyle-\frac{1}{2}(p+q)\cdot p\,,$ (239) $\displaystyle\gamma_{g\theta
g}(p,q)$ $\displaystyle=$ $\displaystyle\frac{1}{2}\frac{p\cdot
q}{q^{2}}\left(\frac{k^{2}}{p^{2}}+\frac{1}{2}-\frac{3}{2}\frac{(k\cdot
p)^{2}}{k^{2}p^{2}}\right)\,,$ $\displaystyle\gamma_{g\theta\epsilon}(p,q)$
$\displaystyle=$ $\displaystyle\frac{1}{2}\left(3\frac{(k\cdot
p)^{2}}{k^{2}p^{2}}-1\right)\,,$ $\displaystyle\gamma_{\epsilon\theta g}(p,q)$
$\displaystyle=$ $\displaystyle\frac{1}{2}\frac{p\cdot
q}{2p^{2}q^{2}k^{2}}((k\cdot p)^{2}-p^{2}k^{2})\,,$
$\displaystyle\gamma_{\epsilon\theta\epsilon}(p,q)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\frac{k\cdot pk\cdot q}{k^{2}p^{2}}\,.$ (240)

###  C.2 Vertices involving at most one vorticity, vector or tensor
perturbation

$\displaystyle\gamma_{\delta w_{i}\delta}(p,q)$ $\displaystyle=$
$\displaystyle\gamma_{Aw_{i}A}(p,q)=\frac{1}{2}\frac{(p\times
q)_{i}}{p^{2}}\,,$ $\displaystyle\gamma_{\theta w_{i}\theta}(p,q)$
$\displaystyle=$ $\displaystyle\frac{1}{2}\left(1+\frac{2p\cdot
q}{q^{2}}\right)\frac{(p\times q)_{i}}{p^{2}}\,,$ $\displaystyle\gamma_{\theta
A\nu_{i}}(p,q)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(1+\frac{2p\cdot q}{q^{2}}\right)(p\times
q)_{i}\,,$ $\displaystyle\gamma_{\theta At_{ij}}(p,q)$ $\displaystyle=$
$\displaystyle-\frac{1}{2}p_{i}p_{j}\,,$ (241)
$\displaystyle\gamma_{gw_{i}g}(p,q)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\frac{1}{p^{2}}\left(\frac{3}{2}\frac{(k\cdot
q)^{2}}{k^{2}q^{2}}+3\frac{k\cdot qp\cdot q}{k^{2}q^{2}}-\frac{p\cdot
q}{q^{2}}-\frac{1}{2}\right)(p\times q)_{i}\,,$
$\displaystyle\gamma_{gw_{i}\epsilon}(p,q)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\,3\frac{k\cdot p}{p^{2}k^{2}}(p\times q)_{i}\,,$
$\displaystyle\gamma_{g\theta\nu_{i}}(p,q)$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\frac{1}{k^{2}p^{2}q^{2}}\left(p\cdot qk^{2}+3k\cdot
pk\cdot q\right)(p\times q)_{i}\,,$ $\displaystyle\gamma_{g\theta
t_{ij}}(p,q)$ $\displaystyle=$ $\displaystyle-\frac{1}{2}\frac{2k^{2}-3p\cdot
q-6k\cdot p}{2k^{2}p^{2}}p_{i}p_{j}\,,$ (242) $\displaystyle\gamma_{\epsilon
w_{i}g}(p,q)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\frac{k^{2}q^{2}-(k\cdot q)^{2}+2p\cdot qk\cdot
p}{2p^{2}q^{2}k^{2}}(p\times q)_{i}\,,$ $\displaystyle\gamma_{\epsilon
w_{i}\epsilon}(p,q)$ $\displaystyle=$ $\displaystyle\frac{1}{2}\frac{k\cdot
q}{p^{2}k^{2}}(p\times q)_{i}\,,$
$\displaystyle\gamma_{\epsilon\theta\nu_{i}}(p,q)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\frac{1}{q^{2}}\left(1-\frac{(k\cdot
p)^{2}}{k^{2}p^{2}}\right)(p\times q)_{i}\,,$
$\displaystyle\gamma_{\epsilon\theta t_{ij}}(p,q)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\frac{2k\cdot q-p\cdot q}{2k^{2}p^{2}}p_{i}p_{j}\,,$
(243) $\displaystyle\gamma_{w_{i}Ag}(p,q)$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\frac{p\cdot q}{q^{2}}(p\times q)_{i}\,,$
$\displaystyle\gamma_{w_{i}A\epsilon}(p,q)$ $\displaystyle=$
$\displaystyle\frac{1}{2}(p\times q)_{i}\,,$
$\displaystyle\gamma_{\nu_{i}\theta g}(p,q)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\frac{p\cdot qk\cdot p(p\times
q)_{i}}{k^{2}p^{2}q^{2}}\,,$
$\displaystyle\gamma_{\nu_{i}\theta\epsilon}(p,q)$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\,2\frac{k\cdot p(p\times q)_{i}}{k^{2}p^{2}}\,,$
(244) $\displaystyle\gamma_{t_{ij}\theta g}(p,q)$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\frac{(p\times q)^{2}\,p\cdot
q}{2k^{2}p^{2}q^{2}}\Bigg{\\{}\left(\delta_{ij}-\frac{k_{i}k_{j}}{k^{2}}\right)-2\frac{(p\times
q)_{i}(p\times q)_{j}}{(p\times q)^{2}}\Bigg{\\}}\,,$
$\displaystyle\gamma_{t_{ij}\theta\epsilon}(p,q)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\frac{(p\times
q)^{2}}{k^{2}p^{2}}\Bigg{\\{}\left(\delta_{ij}-\frac{k_{i}k_{j}}{k^{2}}\right)-2\frac{(p\times
q)_{i}(p\times q)_{j}}{(p\times q)^{2}}\Bigg{\\}}\,.$ (245)

###  C.3 Vertices involving two vorticity, vector or tensor perturbations

$\displaystyle\gamma_{\theta w_{i}w_{j}}(p,q)$ $\displaystyle=$
$\displaystyle-\frac{(p\times q)_{i}(p\times q)_{j}}{p^{2}q^{2}}\,,$
$\displaystyle\gamma_{\epsilon w_{i}\nu_{j}}(p,q)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\,\frac{(p\cdot q-p^{2}+q^{2})(p\times q)_{i}(p\times
q)_{j}}{p^{2}q^{2}k^{2}}\,,$ $\displaystyle{}+\frac{1}{2}\frac{k^{2}(p\cdot
q)\left(\delta_{ij}\,p\cdot q-p_{j}q_{i}\right)}{p^{2}q^{2}k^{2}}\,,$
$\displaystyle\gamma_{gw_{i}\nu_{j}}(p,q)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\,\frac{(p^{2}-5q^{2}-7p\cdot q)(p\times
q)_{i}(p\times q)_{j}}{p^{2}q^{2}k^{2}}\,,$
$\displaystyle{}-\frac{1}{2}\frac{k^{2}(p\cdot q)\left(\delta_{ij}\,p\cdot
q-p_{j}q_{i}\right)}{p^{2}q^{2}k^{2}}\,,$
$\displaystyle\gamma_{gw_{i}t_{jm}}(p,q)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\,\frac{9(p\times
q)_{i}p_{j}+2k^{2}\varepsilon_{ijl}p_{l}}{2k^{2}p^{2}}p_{m}\,,$
$\displaystyle\gamma_{\epsilon w_{i}t_{jm}}(p,q)$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\,\frac{3(p\times
q)_{i}p_{j}+2k^{2}\varepsilon_{ijl}p_{l}}{2k^{2}p^{2}}p_{m}\,,$ (246)
$\displaystyle\gamma_{w_{i}\theta w_{j}}(p,q)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\,\frac{\delta_{ij}p\cdot k-p_{i}p_{j}}{p^{2}}\,,$
$\displaystyle\gamma_{w_{i}A\nu_{j}}(p,q)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\,\frac{(p\times q)_{i}(p\times q)_{j}+(p\cdot
q)\left(p_{j}q_{i}-\delta_{ij}(k\cdot q)\right)}{q^{2}}\,,$
$\displaystyle\gamma_{w_{i}At_{nm}}(p,q)$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\varepsilon_{ijn}k_{j}p_{m}\,,$
$\displaystyle\gamma_{\nu_{i}w_{j}g}(p,q)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\,\frac{\left(\delta_{ij}p\cdot
k-p_{i}q_{j}\right)(q\cdot k)(p\cdot q)}{p^{2}q^{2}k^{2}}$
$\displaystyle{}+\frac{1}{2}\frac{(p\times q)_{i}(p\times
q)_{j}(k^{2}-p^{2})}{p^{2}q^{2}k^{2}}\,,$
$\displaystyle\gamma_{\nu_{i}w_{j}\epsilon}(p,q)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\,\frac{\left(\delta_{ij}p\cdot
k-p_{i}q_{j}\right)(p\cdot k)-(p\times q)_{i}(p\times q)_{j}}{p^{2}k^{2}}\,,$
$\displaystyle\gamma_{\nu_{i}\theta\nu_{j}}(p,q)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\,\frac{\left(\delta_{ij}q\cdot
k-p_{j}q_{i}\right)(p\cdot q)k^{2}+(p\times q)_{i}(p\times
q)_{j}(q^{2}-p^{2})}{p^{2}q^{2}k^{2}}\,,$ $\displaystyle\gamma_{\nu_{i}\theta
t_{mj}}(p,q)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\,\frac{(k^{2}-q^{2})\varepsilon_{inm}k_{n}-(p\times
q)_{i}p_{m}}{k^{2}p^{2}}p_{j}\,,$ $\displaystyle\gamma_{t_{ij}w_{\ell}g}(p,q)$
$\displaystyle=$
$\displaystyle\frac{1}{2}\,\frac{1}{p^{2}q^{2}}\Bigg{[}\varepsilon_{in\ell}\,p_{n}(p\cdot
q)\left(k\cdot
q\frac{k_{j}}{k^{2}}-q_{j}\right)+\varepsilon_{jn\ell}\,p_{n}(p\cdot
q)\left(k\cdot q\frac{k_{i}}{k^{2}}-q_{i}\right)$
$\displaystyle{}+\varepsilon_{\ell
nm}\,p_{n}q_{m}\Bigg{\\{}\left(\delta_{ij}+\frac{k_{i}k_{j}}{k^{2}}\right)\frac{k\cdot
q(k\cdot q+2p\cdot q)}{2k^{2}}+q_{i}q_{j}$
$\displaystyle{}+\frac{1}{2}(k^{2}-p^{2})\left(\delta_{ij}-\frac{k_{i}k_{j}}{k^{2}}+2\frac{k_{i}q_{j}+k_{j}q_{i}}{k^{2}}\right)\Bigg{\\}}\Bigg{]}\,,$
$\displaystyle\gamma_{t_{ij}w_{\ell}\epsilon}(p,q)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\,\frac{p_{n}}{k^{2}p^{2}}\Bigg{[}\varepsilon_{in\ell}\left(k_{j}(k\cdot
p)-k^{2}p_{j})\right)+\varepsilon_{jn\ell}\left(k_{i}(k\cdot
p)-k^{2}p_{i})\right)$ $\displaystyle{}+\varepsilon_{\ell
nm}\,q_{m}\left(k\cdot
p\left(\delta_{ij}+\frac{k_{i}k_{j}}{k^{2}}\right)-k_{i}p_{j}-k_{j}p_{i}\right)\Bigg{]}\,,$
$\displaystyle\gamma_{t_{ij}\theta\nu_{\ell}}(p,q)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\,\frac{\varepsilon_{mn\ell}\,p_{m}q_{n}(p\times
q)^{2}}{k^{2}p^{2}q^{2}}\Bigg{[}\left(\delta_{ij}-\frac{k^{i}k^{j}}{k^{2}}\right)-2\frac{(p\times
q)_{i}(p\times q)_{j}}{(p\times q)^{2}}\Bigg{]}\,,$
$\displaystyle\gamma_{t_{ij}\theta t_{\ell m}}(p,q)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\,\frac{1}{k^{2}p^{2}}\Bigg{[}(q^{2}k_{i}-k^{2}q_{i})(p_{m}\delta_{j\ell}+p_{\ell}\delta_{jm})+(q^{2}k_{j}-k^{2}q_{j})(p_{m}\delta_{i\ell}+p_{\ell}\delta_{im})$
(248) $\displaystyle{}+p\cdot
q\left(\delta_{i\ell}\delta_{jm}+\delta_{im}\delta_{j\ell}-\frac{p_{\ell}p_{m}}{k^{2}}\left(\delta_{ij}-\frac{k_{i}k_{j}}{k^{2}}\right)\right)$
$\displaystyle{}-2\frac{p_{\ell}p_{m}}{k^{2}}\left(k_{i}p_{j}+k_{j}p_{i}-k_{i}k_{j}+q^{2}\delta_{ij}-k\cdot
p\frac{k_{i}k_{j}}{k^{2}}\right)\Bigg{]}\,.$

###  C.4 Vertices involving three vorticity, vector or tensor perturbations

$\displaystyle\gamma_{w_{i}w_{j}w_{\ell}}(p,q)$ $\displaystyle=$
$\displaystyle\frac{\varepsilon_{imj}\,p_{\ell}p_{m}+\delta_{i\ell}(p\times
q)_{j}}{2p^{2}}+\frac{\varepsilon_{im\ell}\,q_{j}q_{m}-\delta_{ij}(p\times
q)_{\ell}}{2q^{2}}\,,$ $\displaystyle\gamma_{\nu_{i}w_{j}\nu_{\ell}}(p,q)$
$\displaystyle=$
$\displaystyle\frac{1}{2}\,\frac{1}{p^{2}q^{2}k^{2}}\Bigg{\\{}\left[\delta_{i\ell}q\cdot
k-p_{\ell}q_{i}\right](p\times q)_{j}(k^{2}-p^{2})$
$\displaystyle{}-\left[\delta_{ij}p\cdot k-p_{i}q_{j}\right](p\times
q)_{\ell}(k^{2}-p^{2})$ $\displaystyle{}-2(p\times q)_{i}(p\times
q)_{j}(p\times q)_{\ell}\Bigg{\\}}\,,$
$\displaystyle\gamma_{\nu_{i}w_{\ell}t_{mj}}(p,q)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\,\frac{2(p\times
q)_{\ell}\,\varepsilon_{inm}k_{n}-p_{m}\left(p_{i}q_{\ell}-k\cdot
p\delta_{i\ell}\right)}{k^{2}p^{2}}p_{j}\,,$ (249)
$\displaystyle\gamma_{t_{ij}w_{\ell}\nu_{m}}(p,q)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\,\frac{p^{2}k_{i}-k^{2}p_{i}}{k^{2}p^{2}q^{2}}\left[\varepsilon_{mns}\,\varepsilon_{jr\ell}\,p_{r}-\varepsilon_{\ell
ns}\,\varepsilon_{jrm}\,q_{r}\right]p_{n}q_{s}$
$\displaystyle{}+\frac{1}{2}\frac{p^{2}k_{j}-k^{2}p_{j}}{k^{2}p^{2}q^{2}}\left[\varepsilon_{mns}\,\varepsilon_{ir\ell}\,p_{r}-\varepsilon_{\ell
ns}\,\varepsilon_{irm}\,q_{r}\right]p_{n}q_{s}$
$\displaystyle{}+\frac{1}{2}\frac{p\cdot
q}{p^{2}q^{2}}p_{n}q_{r}\Bigg{[}\varepsilon_{in\ell}\,\varepsilon_{jrm}+\varepsilon_{jn\ell}\,\varepsilon_{irm}-\varepsilon_{sn\ell}\,\varepsilon_{srm}\left(\delta_{ij}-\frac{k_{i}k_{j}}{k^{2}}\right)\Bigg{]}$
$\displaystyle{}-\frac{1}{2}\frac{\varepsilon_{\ell
ns}\,\varepsilon_{mkr}\,p_{n}p_{k}q_{s}q_{r}}{p^{2}q^{2}}\Bigg{[}\left(\delta_{ij}+\frac{k_{i}k_{j}}{k^{2}}\right)\frac{3p\cdot
q+2q^{2}}{k^{2}}-2\frac{k_{i}q_{j}+k_{j}q_{i}}{k^{2}}-\left(\delta_{ij}-\frac{k_{i}k_{j}}{k^{2}}\right)\Bigg{]}\,,$
$\displaystyle\gamma_{t_{ij}w_{\ell}t_{mn}}(p,q)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\,\frac{1}{2p^{2}}\Bigg{[}2\left(\delta_{ij}-\frac{k_{i}k_{j}}{k^{2}}\right)\left[\varepsilon_{mr\ell}\,p_{n}+\varepsilon_{nr\ell}\,p_{m}\right]p_{r}$
(250)
$\displaystyle{}-2\varepsilon_{ir\ell}\,p_{r}\left(p_{n}\delta_{jm}+p_{m}\delta_{jn}-2\frac{p_{m}p_{n}k_{j}}{k^{2}}\right)-2\varepsilon_{jr\ell}\,p_{r}\left(p_{n}\delta_{im}+p_{m}\delta_{in}-2\frac{p_{m}p_{n}k_{i}}{k^{2}}\right)$
$\displaystyle{}+\varepsilon_{\ell
rs}\,p_{r}q_{s}\Bigg{(}6\frac{p_{m}p_{n}}{k^{2}}\left(\delta_{ij}+\frac{k_{i}k_{j}}{k^{2}}\right)+\delta_{im}\left(\delta_{jn}-4\frac{p_{n}k_{j}}{k^{2}}\right)$
$\displaystyle{}+\delta_{jn}\left(\delta_{ik}-4\frac{p_{m}k_{i}}{k^{2}}\right)+\delta_{in}\left(\delta_{jm}-4\frac{p_{m}k_{j}}{k^{2}}\right)+\delta_{jm}\left(\delta_{in}-4\frac{p_{n}k_{i}}{k^{2}}\right)\Bigg{)}\Bigg{]}\,.$

##  APPENDIX D Third and fourth cumulant

###  D.1 Evolution matrix for scalar perturbations

The scalar part of the block-diagonal evolution matrix $\Omega_{ab}(k,\eta)$
for the perturbation modes up to the fourth cumulant

$\psi^{S}=(\delta,\theta,g,\delta\epsilon,A,\pi,\chi,\kappa,\xi,\psi)\,,$
(251)

in the approximation $\Omega_{m}/f^{2}\to 1$ is given by

$\Omega^{S}=\left(\begin{array}[]{cccccccccc}&-1\\\
-3/2&1/2&k^{2}&k^{2}&k^{2}\epsilon\\\ &-2\epsilon&1&&&1&-3/5\\\ &&&1&&&1/5\\\
&-1\\\ &&-3k^{2}\epsilon&-5k^{2}\epsilon&-k^{2}\omega&3/2&&&-k^{2}\\\
&&&-5k^{2}\epsilon&-k^{2}\omega&&3/2&-5/2k^{2}&5/2k^{2}&4k^{2}\\\
&-4\omega&&&&4\epsilon&&2\\\ &-16/5\omega&&&&4\epsilon&-4/5\epsilon&&2\\\
&&&&&&&&&2\end{array}\right)\,.$ (252)

###  D.2 Vertices involving scalar modes $\pi$ and $\chi$ of the third
cumulant

$\displaystyle\gamma_{gA\pi}(p,q)$ $\displaystyle=$
$\displaystyle\frac{1}{4}\frac{p\cdot q}{q^{2}}-\frac{3}{4}\frac{(k\cdot
q)^{2}p\cdot q}{k^{2}q^{4}}\,,$ $\displaystyle\gamma_{gA\chi}(p,q)$
$\displaystyle=$ $\displaystyle-3\gamma_{\epsilon A\chi}(p,q)\,,$
$\displaystyle\gamma_{\epsilon A\pi}(p,q)$ $\displaystyle=$
$\displaystyle-\frac{1}{4}\frac{p\cdot q}{q^{2}}+\frac{1}{4}\frac{(k\cdot
q)^{2}p\cdot q}{k^{2}q^{4}}\,,$ $\displaystyle\gamma_{\epsilon A\chi}(p,q)$
$\displaystyle=$ $\displaystyle\frac{1}{4}\left(\frac{p\cdot
q}{5q^{2}}+\frac{2}{5}\frac{k\cdot pk\cdot q}{k^{2}q^{2}}-\frac{(k\cdot
q)^{2}p\cdot q}{k^{2}q^{4}}\right)\,,$ (253)

$\displaystyle\gamma_{\pi gg}(p,q)$ $\displaystyle=$
$\displaystyle\frac{k^{2}(p\cdot q)^{2}}{p^{2}q^{2}}+\frac{1}{2}p\cdot
q\left(\frac{k\cdot q}{q^{2}}+\frac{k\cdot p}{p^{2}}\right)\,,$
$\displaystyle\gamma_{\pi\epsilon g}(p,q)$ $\displaystyle=$
$\displaystyle\frac{1}{2}k\cdot q\left(3+\frac{5p\cdot q}{q^{2}}\right)\,,$
$\displaystyle\gamma_{\pi\epsilon\epsilon}(p,q)$ $\displaystyle=$
$\displaystyle\frac{5}{2}k^{2}\,,$ $\displaystyle\gamma_{\pi\theta\pi}(p,q)$
$\displaystyle=$ $\displaystyle\frac{1}{2}\frac{k^{2}p\cdot
q}{p^{2}q^{2}}+\frac{k\cdot q(p\cdot q)^{2}}{p^{2}q^{4}}\,,$
$\displaystyle\gamma_{\pi\theta\chi}(p,q)$ $\displaystyle=$
$\displaystyle\frac{1}{5}\frac{k\cdot q}{q^{2}}+\frac{2}{5}\frac{k\cdot
pp\cdot q}{p^{2}q^{2}}-\frac{k\cdot q(p\cdot q)^{2}}{p^{2}q^{4}}\,,$ (254)

$\displaystyle\gamma_{\chi gg}(p,q)$ $\displaystyle=$
$\displaystyle-\frac{15}{4}\frac{k\cdot qk\cdot pp\cdot
q}{p^{2}q^{2}}+\frac{5}{2}\gamma_{\pi gg}(p,q)\,,$
$\displaystyle\gamma_{\chi\epsilon g}(p,q)$ $\displaystyle=$
$\displaystyle-\frac{15}{4}\left(\frac{k\cdot qp\cdot q}{q^{2}}+\frac{(k\cdot
q)^{3}}{k^{2}q^{2}}\right)+\frac{5}{2}\gamma_{\pi\epsilon g}(p,q)\,,$
$\displaystyle\gamma_{\chi\epsilon\epsilon}(p,q)$ $\displaystyle=$
$\displaystyle-\frac{15}{4}k^{2}+\frac{5}{2}\gamma_{\pi\epsilon\epsilon}(p,q)\,,$
$\displaystyle\gamma_{\chi\theta\pi}(p,q)$ $\displaystyle=$
$\displaystyle-\frac{5}{4}\frac{(k\cdot q)^{2}p\cdot q(k^{2}+2k\cdot
p)}{k^{2}p^{2}q^{4}}+\frac{5}{2}\gamma_{\pi\theta\pi}(p,q)\,,$
$\displaystyle\gamma_{\chi\theta\chi}(p,q)$ $\displaystyle=$
$\displaystyle-\frac{3}{4}\frac{k^{2}p\cdot
q}{p^{2}q^{2}}-\frac{3}{2}\frac{(k\cdot p)^{2}k\cdot
q}{k^{2}p^{2}q^{2}}+\frac{5}{4}\frac{(k\cdot q)^{2}p\cdot q(k^{2}+2k\cdot
p)}{k^{2}p^{2}q^{4}}+\frac{5}{2}\gamma_{\pi\theta\chi}(p,q)\,.$ (255)

##  APPENDIX E The local cosmic energy equation

In this section we review the evolution equation for the total energy (see
e.g. Binney and Tremaine (2008)), given by the sum of kinetic and potential
energy,

$\displaystyle E_{\text{kin}}$ $\displaystyle=$ $\displaystyle\frac{1}{2}\int
d^{3}p\frac{p^{2}}{a^{2}}f(\tau,\bm{x},\bm{p})=\frac{1}{2}(1+\delta)(\sigma_{ii}+v_{i}v_{i})\,,$
$\displaystyle E_{\text{pot}}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\Phi\delta\,.$ (256)

Using the evolution equation for $\sigma_{ij}$, the continuity, Euler and
Poisson equations one obtains

$\displaystyle\partial_{\tau}E_{\text{kin}}+2{\cal
H}E_{\text{kin}}+\nabla_{i}J_{i}^{\text{kin}}$ $\displaystyle=$
$\displaystyle-\Phi\partial_{\tau}\delta\,,$
$\displaystyle\partial_{\tau}E_{\text{pot}}+q{\cal
H}E_{\text{pot}}+\nabla_{i}J_{i}^{\text{pot}}$ $\displaystyle=$
$\displaystyle+\Phi\partial_{\tau}\delta\,,$ (257)

where

$\displaystyle J_{i}^{\text{kin}}$ $\displaystyle=$ $\displaystyle
E_{\text{kin}}v_{i}+(1+\delta)\left(\sigma_{ij}v_{j}+v_{i}\Phi+\frac{1}{2}{\cal
C}_{ijj}\right)\,,$ $\displaystyle J_{i}^{\text{pot}}$ $\displaystyle=$
$\displaystyle\frac{1}{8\pi
G\bar{\rho}_{0}a^{2}}\left(\Phi\partial_{\tau}\nabla_{i}\Phi-(\nabla_{i}\Phi)\partial_{\tau}\Phi\right)\,,$
(258)

with third cumulant ${\cal C}_{ijk}$ and

$q\equiv-\frac{d\ln(a^{2}\bar{\rho}_{0})}{d\ln a}\to 1\,,$ (259)

for the usual scaling $\bar{\rho}_{0}\propto a^{-3}$ of the matter rest energy
density, that we assume throughout. The evolution of the total energy is
therefore given by

$\partial_{\tau}E+{\cal
H}(2E_{\text{kin}}+E_{\text{pot}})+\nabla_{i}J_{i}=0\,,$ (260)

with $J_{i}=J_{i}^{\text{kin}}+J_{i}^{\text{pot}}$, and recovering the energy
conservation law on a static background. Equation (260) is the local version
of the cosmic energy equation Peebles (1980). Rescaling to $\eta=\ln(D)$ and
$\hat{X}=X/(f{\cal H})^{2}$ for $X=\Phi,E,E_{\text{kin}},E_{\text{pot}}$ and
$\hat{X}_{i}=X_{i}/(-f{\cal H})^{3}$ for $X=J,J^{\text{kin}},J^{\text{pot}}$
gives

$\displaystyle\partial_{\eta}\hat{E}_{\text{kin}}+2\left(\frac{3}{2}\frac{\Omega_{m}}{f^{2}}-1\right)\hat{E}_{\text{kin}}-\nabla_{i}\hat{J}_{i}^{\text{kin}}$
$\displaystyle=$ $\displaystyle-\hat{\Phi}\partial_{\eta}\delta\,,$
$\displaystyle\partial_{\eta}\hat{E}_{\text{pot}}+2\left(\frac{3}{2}\frac{\Omega_{m}}{f^{2}}-1-\frac{1}{2f}\right)\hat{E}_{\text{pot}}-\nabla_{i}\hat{J}_{i}^{\text{pot}}$
$\displaystyle=$ $\displaystyle+\hat{\Phi}\partial_{\eta}\delta\,,$ (261)

where

$\displaystyle\hat{J}_{i}^{\text{kin}}$ $\displaystyle=$
$\displaystyle\hat{E}_{\text{kin}}u_{i}+(1+\delta)\left(\epsilon_{ij}u_{j}+u_{i}\hat{\Phi}+\frac{1}{2}\pi_{ijj}\right)\,,$
$\displaystyle\hat{J}_{i}^{\text{pot}}$ $\displaystyle=$
$\displaystyle\frac{f^{2}}{3\Omega_{m}}\left(\hat{\Phi}\partial_{\eta}\nabla_{i}\hat{\Phi}-(\nabla_{i}\hat{\Phi})\partial_{\eta}\hat{\Phi}\right)\,,$
(262)

with $\pi_{ijk}={\cal C}_{ijk}/(-f{\cal H})^{3}$ and

$\displaystyle\hat{E}_{\text{kin}}$ $\displaystyle=$
$\displaystyle\frac{1}{2}(1+\delta)(\epsilon_{ii}+u_{i}u_{i})\,,$
$\displaystyle\hat{E}_{\text{pot}}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\hat{\Phi}\delta\,,$

so that

$\partial_{\eta}\hat{E}+2\left(\frac{3}{2}\frac{\Omega_{m}}{f^{2}}-1\right)\hat{E}-\nabla_{i}\hat{J}_{i}=\frac{1}{f}\hat{E}_{\text{pot}}\,,$
(264)

where $\hat{J}_{i}=\hat{J}_{i}^{\text{kin}}+\hat{J}_{i}^{\text{pot}}$.

Let us check that the local cosmic energy equation is indeed satisfied at
lowest order in perturbation theory, when including second and higher cumulant
contributions as done in this work. We insert the expansion
$\epsilon_{ij}=\epsilon(\eta)\delta_{ij}+\delta\epsilon_{ij}$ as well as the
decomposition of $\delta\epsilon_{ij}$ and $\pi_{ijk}$ from Eqs. (57, V.2).
Keeping terms up to linear order in perturbations (except for contributions
involving $\delta\epsilon$, see below) yields

$\displaystyle\hat{E}_{\text{kin}}$ $\displaystyle=$
$\displaystyle\frac{3}{2}(1+\delta)(\epsilon(\eta)+\delta\epsilon)+\frac{1}{2}g\,,$
$\displaystyle\nabla_{i}\hat{J}_{i}$ $\displaystyle=$
$\displaystyle\frac{5}{2}\epsilon\nabla_{i}u_{i}+\frac{1}{2}\nabla_{i}\pi_{ijj}=\frac{5}{2}\epsilon(\eta)\theta-\frac{1}{2}\pi\,,$
(265)

while $\hat{E}_{\text{pot}}=0$ at this order. In the following we use the EdS
approximation $\Omega_{m}=f=1$. One can check that, as expected, Eq. (264) is
indeed satisfied up to terms of higher order in perturbation theory when using
the linear equations of motion (see Sec. V)

$\partial_{\eta}\delta=\theta\,,\quad(\partial_{\eta}+1)\delta\epsilon=-\chi/5-Q\,,\quad(\partial_{\eta}+1)g=2\theta\epsilon(\eta)-\pi-3\chi/5\,,$
(266)

as well as the background dispersion evolution equation
$(\partial_{\eta}+1)\epsilon(\eta)=Q$. Note that here it is important to
retain the contribution from $Q$ in the equation for $\delta\epsilon$, which
has to cancel with corresponding contributions from the equation for
$\epsilon(\eta)$. Alternatively, this cancellation can be seen directly by
realizing that they appear in the combination $\epsilon+\delta\epsilon$, and
that $Q$ drops out when adding the corresponding equations of motion.

##  APPENDIX F Kernels in linear approximation

The linear solution for the velocity divergence $\theta_{k}$ and scalar
perturbation mode $\bar{g}_{k}=g_{k}/\epsilon$ of the velocity dispersion
tensor in the second-cumulant approximation and for a power-law dependence
$\epsilon=\epsilon_{0}e^{\alpha\eta}$ is given by Eq. (VI.1), with linear
kernels

$\displaystyle F_{1,\theta}(k,\eta)$ $\displaystyle=$ $\displaystyle
F_{1,\delta}(k,\eta)-\frac{2(4+\alpha)k^{2}\epsilon}{(2+\alpha)(5+2\alpha)}\,\;{}_{1}F_{2}\left(\frac{4+4\alpha}{3\alpha};2+\frac{2}{\alpha},2+\frac{5}{2\alpha};\frac{-3k^{2}\epsilon}{\alpha^{2}}\right)\,,$
$\displaystyle F_{1,\bar{g}}(k,\eta)$ $\displaystyle=$
$\displaystyle\frac{2}{2+\alpha}\Bigg{[}\;{}_{1}F_{2}\left(\frac{4+\alpha}{3\alpha};2+\frac{2}{\alpha},1+\frac{5}{2\alpha};\frac{-3k^{2}\epsilon}{\alpha^{2}}\right)$
(267)
$\displaystyle-\frac{(4+\alpha)k^{2}\epsilon}{(1+\alpha)(5+2\alpha)}\,\;{}_{1}F_{2}\left(\frac{4+4\alpha}{3\alpha};3+\frac{2}{\alpha},2+\frac{5}{2\alpha};\frac{-3k^{2}\epsilon}{\alpha^{2}}\right)\Bigg{]}\,.$

The coefficients for the expansion Eq. (143) for large $k^{2}\epsilon$ are
given in Table 2.

Table 2: Coefficients in the expansion Eq. (143) of the linear kernels $F_{1,a}(k,\eta)$ for large $k^{2}\epsilon$. $a$ | $d_{a}$ | $e_{a}$ | $D_{a}$ | $E_{a}$ | $\varphi_{a}$
---|---|---|---|---|---
$\delta$ | $\frac{2}{3\alpha}(4+\alpha)$ | $\frac{19+7\alpha}{6\alpha}$ | $\frac{5\alpha^{4\frac{2-\alpha}{3\alpha}}\Gamma(\frac{2}{\alpha})\Gamma(\frac{5}{2\alpha})}{\Gamma(\frac{4\alpha+7}{6\alpha})\Gamma(\frac{2(1+\alpha)}{3\alpha})}$ | $-\frac{5\alpha^{\frac{19-5\alpha}{6\alpha}}\Gamma(\frac{2}{\alpha})\Gamma(\frac{5}{2\alpha})}{\sqrt{\pi}\Gamma(\frac{4+\alpha}{3\alpha})}$ | $\frac{5\alpha-19}{12\alpha}\pi$
$\theta$ | $\frac{2}{3\alpha}(4+\alpha)$ | $\frac{19+\alpha}{6\alpha}$ | $-\frac{5\alpha^{4\frac{2-\alpha}{3\alpha}}(\alpha+1)\Gamma(\frac{2}{\alpha})\Gamma(\frac{5}{2\alpha})}{3\Gamma(\frac{4\alpha+7}{6\alpha})\Gamma(\frac{2(1+\alpha)}{3\alpha})}$ | $-\frac{5\alpha^{\frac{19-5\alpha}{6\alpha}}\Gamma(\frac{2}{\alpha})\Gamma(\frac{5}{2\alpha})}{\sqrt{\pi}\Gamma(\frac{4+\alpha}{3\alpha})}$ | $\frac{11\alpha-19}{12\alpha}\pi$
$\bar{g}$ | $\frac{2}{3\alpha}(4+\alpha)$ | $\frac{19+7\alpha}{6\alpha}$ | $-\frac{5\alpha^{4\frac{2-\alpha}{3\alpha}}\Gamma(\frac{2}{\alpha})\Gamma(\frac{5}{2\alpha})}{\Gamma(\frac{4\alpha+7}{6\alpha})\Gamma(\frac{2(1+\alpha)}{3\alpha})}$ | $-\frac{10\alpha^{\frac{19-5\alpha}{6\alpha}}\Gamma(\frac{2}{\alpha})\Gamma(\frac{5}{2\alpha})}{\sqrt{\pi}\Gamma(\frac{4+\alpha}{3\alpha})}$ | $\frac{5\alpha-19}{12\alpha}\pi$

The full linear solution of the coupled system Eq. (VI.1) for
$\delta_{k},\theta_{k},g_{k}$, that includes all eigenmodes, can be found by
making the ansatz

$\delta_{k}=e^{\lambda\eta}\sum_{n\geq
0}c_{n}(\epsilon(\eta)k^{2})^{n}\,,\qquad\theta_{k}=e^{\lambda\eta}\sum_{n\geq
0}d_{n}(\epsilon(\eta)k^{2})^{n}\,,\qquad
k^{2}g_{k}=e^{\lambda\eta}\sum_{n\geq 0}e_{n}(\epsilon(\eta)k^{2})^{n}\,,$
(268)

with some exponent $\lambda$ and coefficients $c_{n},d_{n},e_{n}$. Inserting
this ansatz into Eq. (VI.1) and assuming $\epsilon\propto e^{\alpha\eta}$
yields the recursions

$c_{n}(\lambda+\alpha n)=d_{n}\,,\qquad d_{n}(\lambda+\alpha
n+\frac{1}{2})=\frac{3}{2}c_{n}-c_{n-1}-e_{n}\,,\qquad e_{n}(\lambda+\alpha
n+1)=2d_{n-1}\,,$ (269)

for $n\geq 1$. They can be combined into

$c_{n}=-c_{n-1}\frac{3(\lambda+\alpha n)+1-2\alpha}{(\lambda+\alpha
n-1)(\lambda+\alpha n+1)(\lambda+\alpha n+\frac{3}{2})}\,,$ (270)

with explicit expression given by

$c_{n}=\left(\frac{-3}{\alpha^{2}}\right)^{n}c_{0}\frac{\Gamma(n+p_{1})}{\Gamma(p_{1})}\frac{\Gamma(q_{1})}{\Gamma(n+q_{1})}\frac{\Gamma(q_{2})}{\Gamma(n+q_{2})}\frac{\Gamma(q_{3})}{\Gamma(n+q_{3})}\,,$
(271)

where

$p_{1}\equiv\frac{3\lambda+1+\alpha}{3\alpha}\,,\qquad q_{1}\equiv
1+\frac{\lambda-1}{\alpha}\,,\qquad q_{2}\equiv
1+\frac{\lambda+1}{\alpha}\,,\qquad q_{3}\equiv
1+\frac{\lambda+\frac{3}{2}}{\alpha}\,.$ (272)

The sum over all $n$ yields a solution in terms of generalized hypergeometric
functions. To find the allowed values for $\lambda$ it is sufficient to
consider Eq. (VI.1) in the limit $\epsilon\to 0$, which yields a linear set of
three algebraic equations for $(c_{0},d_{0},e_{0})$,

$\left(\begin{array}[]{ccc}\lambda&-1&0\\\
-\frac{3}{2}&\lambda+\frac{1}{2}&1\\\
0&0&\lambda+1\end{array}\right)\left(\begin{array}[]{c}c_{0}\\\ d_{0}\\\
e_{0}\end{array}\right)=0\,.$ (273)

Non-trivial solutions exist for $\lambda=+1$ (with $c_{0}=d_{0},e_{0}=0$),
$\lambda=-3/2$ (with $c_{0}=-\frac{2}{3}d_{0},e_{0}=0$) and $\lambda=-1$ (with
$c_{0}=-d_{0}=e_{0}$). Each possibility yields one of the three linearly
independent solutions. Inserting the values for $\lambda$ in the recursion
relation, and building a generic linear combination of the three solutions,
yields the general solution Eq. (VI.1). The solution for all three
perturbation modes can be written as

$\left(\begin{array}[]{c}\delta_{k}\\\ \theta_{k}\\\
g_{k}\end{array}\right)={\cal M}_{k}(\eta)\left(\begin{array}[]{c}A\\\ B\\\
C\end{array}\right)\,,$ (274)

with free coefficients $A,B,C$ and $3\times 3$ matrix given by the tensor
product

${\cal M}_{k}(\eta)\equiv\left(\begin{array}[]{c}1\\\ \partial_{\eta}\\\
(-\partial_{\eta}^{2}-\frac{1}{2}\partial_{\eta}+\frac{3}{2}-\epsilon
k^{2})/k^{2}\end{array}\right)\otimes\left(e^{\eta}\,{}_{1}F_{2}\left(p_{1};q_{2},q_{3};x\right),\,e^{-\frac{3}{2}\eta}\,{}_{1}F_{2}\left(p_{1};q_{1},q_{2};x\right),\,e^{-\eta}\,{}_{1}F_{2}\left(p_{1};q_{1},q_{3};x\right)\right)\,,$
(275)

with $x\equiv\frac{-3k^{2}\epsilon(\eta)}{\alpha^{2}}$, and $p_{i},q_{i}$
evaluated for $\lambda=1,-3/2,-1$ for the three terms, respectively. This
result can be used to give an analytic expression for the linear propagator

$G_{k}(\eta,\eta^{\prime})={\cal M}_{k}(\eta)[{\cal
M}_{k}(\eta^{\prime})]^{-1}\,.$ (276)

A full treatment of the scalar modes in second cumulant approximation requires
to include also the mode $\delta\epsilon_{k}$ as well as $A_{k}$. Including
$\delta\epsilon_{k}$ leads to an additional decaying mode solution
$g_{k}=-\delta\epsilon_{k}=De^{-\eta}$ and $\delta_{k}=\theta_{k}=A_{k}=0$,
with free coefficient $D$. This solution remains valid also for $\epsilon
k^{2}\gg 1$, and can easily be included in the linear propagator by extending
Eq. (274) by a fourth row. The solutions given previously remain valid when
including $\delta\epsilon_{k}$.

Finally, taking the $A_{k}$ mode into account formally yields an additional
eigenvalue $\lambda=0$, related to the freedom to choose different initial
conditions for $\delta_{k}$ and $A_{k}$. This is irrelevant for the linear
evolution, but the additional linearly independent solution enters in the
linear propagator. The corresponding additional solution is given by

$\delta_{k}=E\,\left[{}_{2}F_{3}\left(1,p_{1};q_{1},q_{2},q_{3};x\right)-1\right]\,,$
(277)

with free coefficient $E$ and $p_{i},q_{i}$ evaluated for $\lambda=0$, and
$A_{k}=E+\delta_{k},\theta_{k}=\partial_{\eta}\delta_{k}$,
$k^{2}g_{k}=(-\partial_{\eta}^{2}-\frac{1}{2}\partial_{\eta}+\frac{3}{2}-\epsilon
k^{2})\delta_{k}$, and $\delta\epsilon_{k}=0$. The most general solution for
all five scalar modes is therefore given by

$\left(\begin{array}[]{c}\delta_{k}\\\ \theta_{k}\\\ g_{k}\\\
\delta\epsilon_{k}\\\ A_{k}\end{array}\right)={\cal M}_{k}^{(5\times
5)}(\eta)\left(\begin{array}[]{c}A\\\ B\\\ C\\\ D\\\ E\end{array}\right)\,,$
(278)

with ${\cal M}_{k}^{(5\times 5)}(\eta)$ given by the $5\times 5$ matrix

$\left(\begin{array}[]{ccc|cc}&&&0&{}_{2}F_{3}\left(1,p_{1};q_{1},q_{2},q_{3};x\right)-1\\\
&{\cal
M}_{k}(\eta)&&0&\partial_{\eta}\,[{}_{2}F_{3}\left(1,p_{1};q_{1},q_{2},q_{3};x\right)-1]\\\
&&&-e^{-\eta}&D_{\eta}[{}_{2}F_{3}\left(1,p_{1};q_{1},q_{2},q_{3};x\right)-1]\\\
\hline\cr 0&0&0&e^{-\eta}&0\\\
e^{\eta}\,{}_{1}F_{2}\left(p_{1};q_{2},q_{3};x\right)&e^{-\frac{3}{2}\eta}\,{}_{1}F_{2}\left(p_{1};q_{1},q_{2};x\right)&e^{-\eta}\,{}_{1}F_{2}\left(p_{1};q_{1},q_{3};x\right)&0&{}_{2}F_{3}\left(1,p_{1};q_{1},q_{2},q_{3};x\right)\end{array}\right)\,,$
(279)

where
$D_{\eta}\equiv(-\partial_{\eta}^{2}-\frac{1}{2}\partial_{\eta}+\frac{3}{2}-\epsilon
k^{2})/k^{2}$. The corresponding linear propagator reads

$G_{k}^{(5\times 5)}(\eta,\eta^{\prime})={\cal M}_{k}^{(5\times
5)}(\eta)[{\cal M}_{k}^{(5\times 5)}(\eta^{\prime})]^{-1}\,.$ (280)

##  APPENDIX G Evolution equations for the full hierarchy of cumulants

The linear evolution equation obtained from expanding Eq. (157) in powers of
$L$ is given by

$\displaystyle\left[\partial_{\eta}+1+\left(\frac{3}{2}\frac{\Omega_{m}}{f^{2}}-1\right)\left(\ell+2n\right)\right]{\cal
C}_{\ell,2n}$ $\displaystyle=$
$\displaystyle\frac{\\{1,k^{2}\\}}{2\ell+1}{\cal
R}_{\ell,2n}+\frac{1}{2}\delta_{\ell
1}\delta_{n0}\frac{\Omega_{m}}{f^{2}}{\cal C}_{0,0}\,,$ (281)

where $\\{A,B\\}=A$ for even $\ell$, and $B$ for odd $\ell$, and

$\displaystyle{\cal R}_{\ell,2n}$ $\displaystyle\equiv$
$\displaystyle\sum_{m=0}^{n}\frac{(2n)!}{(2m+1)!(2n-2m)!}{\cal
E}_{2m+2}\left((\ell+1)(2n-2m)(2n-2m-1){\cal C}_{\ell+1,2(n-m-1)}-\ell{\cal
C}_{\ell-1,2(n-m)}\right)$ (282)
$\displaystyle+(\ell+1)\left(2\ell+3+2n\right){\cal
C}_{\ell+1,2n}-\frac{\ell}{2n+1}{\cal C}_{\ell-1,2(n+1)}\,.$

We solve this system of equations with growing mode initial conditions

$\displaystyle{\cal C}_{0,0}$ $\displaystyle\to$ $\displaystyle 1\,,$
$\displaystyle{\cal C}_{1,0}$ $\displaystyle\to$ $\displaystyle\frac{1}{3}\,,$
$\displaystyle{\cal C}_{\ell,2n}$ $\displaystyle\to$ $\displaystyle
0,\qquad\ell+2n\geq 2\,,$ (283)

for $\eta\to-\infty$, as appropriate for cold dark matter and adiabatic
initial conditions.

The evolution equation for the expectation values of even cumulants obtained
from Eq. (89) and Taylor expanding Eq. (159) in $L$ is given by

$\left[\partial_{\eta}+\left(\frac{3}{2}\frac{\Omega_{m}}{f^{2}}-1\right)2n\right]{\cal
E}_{2n}=Q_{{\cal E}_{2n}}\,,$ (284)

with

$\displaystyle Q_{{\cal E}_{2n}}$ $\displaystyle=$ $\displaystyle
4\pi\int_{0}^{\infty}dk\,k^{2}\,e^{2\eta}\,P_{0}(k)\sum_{\ell=0}^{n}(\ell+1)\sum_{m_{1},m_{2}=0}^{n-\ell}\delta_{m_{1}+m_{2},n-\ell}\frac{(2n)!}{(2m_{1})!(2m_{2})!}$
(285) $\displaystyle(2m_{2}-2m_{1}-2\ell-3){\cal C}_{\ell+1,2m_{1}}{\cal
C}_{\ell,2m_{2}}\,.$

The evolution equations for the rescaled variables $T_{\ell,2n}={\cal
C}_{\ell,2n}/\epsilon^{n+[\ell/2]}$ and $\bar{\cal E}_{2n}={\cal
E}_{2n}/\epsilon^{n}$ read

$\displaystyle\left[\partial_{\eta}+1+\left(\frac{3}{2}\frac{\Omega_{m}}{f^{2}}-1\right)\left(\ell+2n\right)+(n+\\{\ell/2,(\ell-1)/2\\})(\partial_{\eta}\ln\epsilon)\right]T_{\ell,2n}$
$\displaystyle=$ $\displaystyle\frac{\\{1,\epsilon k^{2}\\}}{2\ell+1}\bar{\cal
R}_{\ell,2n}+\frac{1}{2}\delta_{\ell
1}\delta_{n0}\frac{\Omega_{m}}{f^{2}}T_{0,0}\,,$ (286)
$\left[\partial_{\eta}+\left(\frac{3}{2}\frac{\Omega_{m}}{f^{2}}-1\right)2n+n(\partial_{\eta}\ln\epsilon)\right]\bar{\cal
E}_{2n}=\bar{Q}_{{\cal E}_{2n}}\,,$ (287)

where $\bar{\cal R}_{\ell,2n}={\cal R}_{\ell,2n}|_{{\cal C}\to T,{\cal
E}\to\bar{\cal E}}$ and $\bar{Q}_{{\cal E}_{2n}}=Q_{{\cal E}_{2n}}|_{{\cal
C}\to T}$.

Table 3: Factor $q_{c_{\text{max}}}(\lambda)$ contributing to the characteristic polynomial related to the asymptotic behavior of the linear kernels for large $\epsilon k^{2}$ when taking scalar perturbation modes of cumulants up to order $c_{\text{max}}$ into account. Stability requires that all roots of $q_{c_{\text{max}}}(\lambda)$ lie on the negative real axis or are zero. We set $X_{10}\equiv 945+1575\bar{\cal E}_{4}^{2}+630\bar{\cal E}_{6}-210\bar{\cal E}_{4}(15+\bar{\cal E}_{6})-45\bar{\cal E}_{8}+\bar{\cal E}_{10}$ and $X_{12}\equiv 10395+51975\bar{\cal E}_{4}^{2}-5775\bar{\cal E}_{4}^{3}+13860\bar{\cal E}_{6}+462\bar{\cal E}_{6}^{2}-495\bar{\cal E}_{4}(105+28\bar{\cal E}_{6}-\bar{\cal E}_{8})-1485\bar{\cal E}_{8}+66\bar{\cal E}_{10}-\bar{\cal E}_{12}$. $c_{\text{max}}$ | $q_{c_{\text{max}}}(\lambda)$
---|---
$3$ | $\lambda^{2}+2\lambda+\frac{1}{9}(3-\bar{\cal E}_{4})$
$4$ | $\lambda^{2}+\frac{10}{3}\lambda+\frac{5}{9}(3-\bar{\cal E}_{4})$
$5$ | $\lambda^{3}+5\lambda^{2}+\frac{5}{3}\lambda(3-\bar{\cal E}_{4})+\frac{1}{27}(15-15\bar{\cal E}_{4}+\bar{\cal E}_{6})$
$6$ | $\lambda^{3}+7\lambda^{2}+\frac{35}{9}\lambda(3-\bar{\cal E}_{4})+\frac{7}{27}(15-15\bar{\cal E}_{4}+\bar{\cal E}_{6})$
$7$ | $\lambda^{4}+\frac{28}{3}\lambda^{3}+\frac{70}{9}\lambda^{2}(3-\bar{\cal E}_{4})+\frac{28}{27}\lambda(15-15\bar{\cal E}_{4}+\bar{\cal E}_{6})+\frac{1}{81}(105-210\bar{\cal E}_{4}+35\bar{\cal E}_{4}^{2}+28\bar{\cal E}_{6}-\bar{\cal E}_{8})$
$8$ | $\lambda^{4}+12\lambda^{3}+14\lambda^{2}(3-\bar{\cal E}_{4})+\frac{28}{9}\lambda(15-15\bar{\cal E}_{4}+\bar{\cal E}_{6})+\frac{1}{9}(105-210\bar{\cal E}_{4}+35\bar{\cal E}_{4}^{2}+28\bar{\cal E}_{6}-\bar{\cal E}_{8})$
$9$ | $\lambda^{5}+15\lambda^{4}+\frac{70}{3}\lambda^{3}(3-\bar{\cal E}_{4})+\frac{70}{9}\lambda^{2}(15-15\bar{\cal E}_{4}+\bar{\cal E}_{6})+\frac{5}{9}\lambda(105-210\bar{\cal E}_{4}+35\bar{\cal E}_{4}^{2}+28\bar{\cal E}_{6}-\bar{\cal E}_{8})+\frac{1}{243}X_{10}$
$10$ | $\lambda^{5}+\frac{55}{3}\lambda^{4}+\frac{110}{3}\lambda^{3}(3-\bar{\cal E}_{4})+\frac{154}{9}\lambda^{2}(15-15\bar{\cal E}_{4}+\bar{\cal E}_{6})+\frac{55}{27}\lambda(105-210\bar{\cal E}_{4}+35\bar{\cal E}_{4}^{2}+28\bar{\cal E}_{6}-\bar{\cal E}_{8})+\frac{11}{243}X_{10}$
$11$ | $\lambda^{6}+22\lambda^{5}+55\lambda^{4}(3-\bar{\cal E}_{4})+\frac{308}{9}\lambda^{3}(15-15\bar{\cal E}_{4}+\bar{\cal E}_{6})+\frac{55}{9}\lambda^{2}(105-210\bar{\cal E}_{4}+35\bar{\cal E}_{4}^{2}+28\bar{\cal E}_{6}-\bar{\cal E}_{8})+\frac{22}{81}\lambda X_{10}+\frac{1}{729}X_{12}$
$12$ | $\lambda^{6}+26\lambda^{5}+\frac{715}{9}\lambda^{4}(3-\bar{\cal E}_{4})+\frac{572}{9}\lambda^{3}(15-15\bar{\cal E}_{4}+\bar{\cal E}_{6})+\frac{143}{9}\lambda^{2}(105-210\bar{\cal E}_{4}+35\bar{\cal E}_{4}^{2}+28\bar{\cal E}_{6}-\bar{\cal E}_{8})+\frac{286}{243}\lambda X_{10}+\frac{13}{729}X_{12}$

The stability conditions for $c_{\text{max}}=7$ are given by those for
$c_{\text{max}}=5$, and in addition

$\displaystyle 0$ $\displaystyle\leq$ $\displaystyle 6+5\bar{\cal E}_{4}\,,$
$\displaystyle 0$ $\displaystyle\leq$ $\displaystyle 105-210\bar{\cal
E}_{4}+35\bar{\cal E}_{4}^{2}+28\bar{\cal E}_{6}-\bar{\cal E}_{8}\,,$
$\displaystyle 0$ $\displaystyle\leq$ $\displaystyle 630-525\bar{\cal
E}_{4}-35\bar{\cal E}_{4}^{2}+21\bar{\cal E}_{6}+\bar{\cal E}_{8}\,,$
$\displaystyle 0$ $\displaystyle\leq$ $\displaystyle 70(216+108\bar{\cal
E}_{4}-378\bar{\cal E}_{4}^{2}+215\bar{\cal E}_{4}^{3})-14(-108+180\bar{\cal
E}_{4}+95\bar{\cal E}_{4}^{2})\bar{\cal E}_{6}+63\bar{\cal
E}_{6}^{2}+(66+95\bar{\cal E}_{4}+3\bar{\cal E}_{6})\bar{\cal E}_{8}\,,$
$\displaystyle 0$ $\displaystyle\leq$ $\displaystyle 70(72+132\bar{\cal
E}_{4}+2\bar{\cal E}_{4}^{2}+85\bar{\cal E}_{4}^{3})-1680\bar{\cal
E}_{4}\bar{\cal E}_{6}-63\bar{\cal E}_{6}^{2}+(6+5\bar{\cal E}_{4})\bar{\cal
E}_{8}\,,$ $\displaystyle 0$ $\displaystyle\leq$ $\displaystyle
34300(5184+5184\bar{\cal E}_{4}-8208\bar{\cal E}_{4}^{2}+7008\bar{\cal
E}_{4}^{3}-9204\bar{\cal E}_{4}^{4}+16340\bar{\cal E}_{4}^{5}+1445\bar{\cal
E}_{4}^{6})$ (288) $\displaystyle-54880(432+1296\bar{\cal E}_{4}-2016\bar{\cal
E}_{4}^{2}+4100\bar{\cal E}_{4}^{3}+1805\bar{\cal E}_{4}^{4})\bar{\cal
E}_{6}+1372(-4104+13500\bar{\cal E}_{4}+15390\bar{\cal
E}_{4}^{2}+2275\bar{\cal E}_{4}^{3})\bar{\cal E}_{6}^{2}$
$\displaystyle-592704\bar{\cal E}_{6}^{3}-64827\bar{\cal
E}_{6}^{4}+\big{[}-7840(-108-324\bar{\cal E}_{4}-441\bar{\cal
E}_{4}^{2}-290\bar{\cal E}_{4}^{3}+170\bar{\cal
E}_{4}^{4})+4704(-72-150\bar{\cal E}_{4}+95\bar{\cal E}_{4}^{2})\bar{\cal
E}_{6}$ $\displaystyle+882(-54+35\bar{\cal E}_{4})\bar{\cal
E}_{6}^{2}\big{]}\bar{\cal E}_{8}-7(396+1140\bar{\cal E}_{4}+335\bar{\cal
E}_{4}^{2}+72\bar{\cal E}_{6})\bar{\cal E}_{8}^{2}-\bar{\cal E}_{8}^{3}\,.$

The stability conditions for $c_{\text{max}}=8$ are given by those for
$c_{\text{max}}=6$, and in addition

$\displaystyle 0$ $\displaystyle\leq$ $\displaystyle 6+7\bar{\cal E}_{4}\,,$
$\displaystyle 0$ $\displaystyle\leq$ $\displaystyle 105-210\bar{\cal
E}_{4}+35\bar{\cal E}_{4}^{2}+28\bar{\cal E}_{6}-\bar{\cal E}_{8}\,,$
$\displaystyle 0$ $\displaystyle\leq$ $\displaystyle 210-105\bar{\cal
E}_{4}-35\bar{\cal E}_{4}^{2}-7\bar{\cal E}_{6}+\bar{\cal E}_{8}\,,$
$\displaystyle 0$ $\displaystyle\leq$ $\displaystyle 210(24+36\bar{\cal
E}_{4}-34\bar{\cal E}_{4}^{2}+49\bar{\cal E}_{4}^{3})-14(-12+252\bar{\cal
E}_{4}+91\bar{\cal E}_{4}^{2})\bar{\cal E}_{6}-49\bar{\cal
E}_{6}^{2}+(78+147\bar{\cal E}_{4}+7\bar{\cal E}_{6})\bar{\cal E}_{8}\,,$
$\displaystyle 0$ $\displaystyle\leq$ $\displaystyle 14(72+204\bar{\cal
E}_{4}+90\bar{\cal E}_{4}^{2}+203\bar{\cal E}_{4}^{3})-784\bar{\cal
E}_{4}\bar{\cal E}_{6}-49\bar{\cal E}_{6}^{2}+(6+7\bar{\cal E}_{4})\bar{\cal
E}_{8}\,,$ $\displaystyle 0$ $\displaystyle\leq$ $\displaystyle
980(1728+5184\bar{\cal E}_{4}+3024\bar{\cal E}_{4}^{2}+7776\bar{\cal
E}_{4}^{3}+2628\bar{\cal E}_{4}^{4}+22932\bar{\cal E}_{4}^{5}+5887\bar{\cal
E}_{4}^{6})$ (289) $\displaystyle-4704(144+720\bar{\cal E}_{4}+384\bar{\cal
E}_{4}^{2}+2436\bar{\cal E}_{4}^{3}+1183\bar{\cal E}_{4}^{4})\bar{\cal
E}_{6}+588(-120+1932\bar{\cal E}_{4}+1274\bar{\cal E}_{4}^{2}+147\bar{\cal
E}_{4}^{3})\bar{\cal E}_{6}^{2}$ $\displaystyle+21952\bar{\cal
E}_{6}^{3}-7203\bar{\cal E}_{6}^{4}+\big{[}-672(-36-180\bar{\cal
E}_{4}-411\bar{\cal E}_{4}^{2}-294\bar{\cal E}_{4}^{3}+203\bar{\cal
E}_{4}^{4})+672(-24-42\bar{\cal E}_{4}+91\bar{\cal E}_{4}^{2})\bar{\cal
E}_{6}$ $\displaystyle+2058(-2+3\bar{\cal E}_{4})\bar{\cal
E}_{6}^{2}\big{]}\bar{\cal E}_{8}-3(156+588\bar{\cal E}_{4}+259\bar{\cal
E}_{4}^{2}+56\bar{\cal E}_{6})\bar{\cal E}_{8}^{2}-\bar{\cal E}_{8}^{3}\,.$

## References

  * Bernardeau  _et al._ (2002a) F. Bernardeau, S. Colombi, E. Gaztanaga, and R. Scoccimarro, Phys. Rep. 367, 1 (2002a).
  * Colombi  _et al._ (1996) S. Colombi, F. R. Bouchet, and L. Hernquist, ApJ 465, 14 (1996).
  * Makino  _et al._ (1992) N. Makino, M. Sasaki, and Y. Suto, Phys. Rev. D 46, 585 (1992).
  * Scoccimarro and Frieman (1996) R. Scoccimarro and J. A. Frieman, ApJ 473, 620 (1996), arXiv:astro-ph/9602070 .
  * Cooray and Sheth (2002) A. Cooray and R. Sheth, Phys. Rep. 372, 1 (2002).
  * Nishimichi  _et al._ (2016) T. Nishimichi, F. Bernardeau, and A. Taruya, Physics Letters B 762, 247 (2016), arXiv:1411.2970 .
  * Nishimichi  _et al._ (2017) T. Nishimichi, F. Bernardeau, and A. Taruya, Phys. Rev. D 96, 123515 (2017), arXiv:1708.08946 [astro-ph.CO] .
  * Carrasco  _et al._ (2012) J. J. M. Carrasco, M. P. Hertzberg, and L. Senatore, ArXiv e-prints (2012), arXiv:1206.2926 [astro-ph.CO] .
  * Baumann  _et al._ (2012) D. Baumann, A. Nicolis, L. Senatore, and M. Zaldarriaga, J. Cosmology Astropart. Phys 7, 051 (2012), arXiv:1004.2488 [astro-ph.CO] .
  * Pueblas and Scoccimarro (2009) S. Pueblas and R. Scoccimarro, Phys. Rev. D 80, 043504 (2009), arXiv:0809.4606 .
  * Garny  _et al._ (2022) M. Garny, D. Laxhuber, and R. Scoccimarro, to appear (2022).
  * McDonald (2011) P. McDonald, J. Cosmology Astropart. Phys 2011, 032 (2011), arXiv:0910.1002 [astro-ph.CO] .
  * Erschfeld and Floerchinger (2019) A. Erschfeld and S. Floerchinger, J. Cosmology Astropart. Phys 2019, 039 (2019), arXiv:1812.06891 [astro-ph.CO] .
  * Erschfeld (2021) A. A. Erschfeld,  _Functional methods for cosmic structure formation._ , Ph.D. thesis, U. Heidelberg (main) (2021).
  * Buchert and Domínguez (1998) T. Buchert and A. Domínguez, A&A 335, 395 (1998), arXiv:astro-ph/9702139 [astro-ph] .
  * Adler and Buchert (1999) S. Adler and T. Buchert, A&A 343, 317 (1999), arXiv:astro-ph/9806320 [astro-ph] .
  * Tatekawa  _et al._ (2002) T. Tatekawa, M. Suda, K.-I. Maeda, M. Morita, and H. Anzai, Phys. Rev. D 66, 064014 (2002), arXiv:astro-ph/0205017 [astro-ph] .
  * Morita and Tatekawa (2001) M. Morita and T. Tatekawa, MNRAS 328, 815 (2001), arXiv:astro-ph/0108289 [astro-ph] .
  * Colombi (2015) S. Colombi, MNRAS 446, 2902 (2015), arXiv:1411.4165 .
  * Aviles (2016) A. Aviles, Phys. Rev. D 93, 063517 (2016), arXiv:1512.07198 [astro-ph.CO] .
  * Cusin  _et al._ (2017) G. Cusin, V. Tansella, and R. Durrer, Phys. Rev. D 95, 063527 (2017), arXiv:1612.00783 [astro-ph.CO] .
  * Taruya and Colombi (2017) A. Taruya and S. Colombi, MNRAS 470, 4858 (2017), arXiv:1701.09088 [astro-ph.CO] .
  * Rampf and Frisch (2017) C. Rampf and U. Frisch, MNRAS 471, 671 (2017), arXiv:1705.08456 [astro-ph.CO] .
  * McDonald and Vlah (2018) P. McDonald and Z. Vlah, Phys. Rev. D 97, 023508 (2018), arXiv:1709.02834 [astro-ph.CO] .
  * Saga  _et al._ (2018) S. Saga, A. Taruya, and S. Colombi, Phys. Rev. Lett. 121, 241302 (2018), arXiv:1805.08787 [astro-ph.CO] .
  * Halle  _et al._ (2020) A. Halle, T. Nishimichi, A. Taruya, S. Colombi, and F. Bernardeau, MNRAS 499, 1769 (2020), arXiv:2001.10417 [astro-ph.CO] .
  * Rampf and Hahn (2021) C. Rampf and O. Hahn, MNRAS 501, L71 (2021), arXiv:2010.12584 [astro-ph.CO] .
  * Valageas (2010) P. Valageas, ArXiv e-prints (2010), arXiv:1009.0106 [astro-ph.CO] .
  * Valageas and Nishimichi (2011a) P. Valageas and T. Nishimichi, A&A 527, A87 (2011a), arXiv:1009.0597 [astro-ph.CO] .
  * Valageas and Nishimichi (2011b) P. Valageas and T. Nishimichi, A&A 532, A4 (2011b), arXiv:1102.0641 [astro-ph.CO] .
  * Valageas  _et al._ (2013) P. Valageas, T. Nishimichi, and A. Taruya, Phys. Rev. D 87, 083522 (2013), arXiv:1302.4533 [astro-ph.CO] .
  * Seljak and Vlah (2015) U. Seljak and Z. Vlah, Phys. Rev. D 91, 123516 (2015), arXiv:1501.07512 [astro-ph.CO] .
  * Abel  _et al._ (2012) T. Abel, O. Hahn, and R. Kaehler, MNRAS 427, 61 (2012), arXiv:1111.3944 [astro-ph.CO] .
  * Hahn  _et al._ (2013) O. Hahn, T. Abel, and R. Kaehler, MNRAS 434, 1171 (2013), arXiv:1210.6652 [astro-ph.CO] .
  * Colombi and Touma (2014) S. Colombi and J. Touma, MNRAS 441, 2414 (2014).
  * Hahn and Angulo (2016) O. Hahn and R. E. Angulo, MNRAS 455, 1115 (2016), arXiv:1501.01959 [astro-ph.CO] .
  * Sousbie and Colombi (2016) T. Sousbie and S. Colombi, Journal of Computational Physics 321, 644 (2016), arXiv:1509.07720 [physics.comp-ph] .
  * Stucker  _et al._ (2020) J. Stucker, O. Hahn, R. E. Angulo, and S. D. M. White, Mon. Not. Roy. Astron. Soc. 495, 4943 (2020), arXiv:1909.00008 [astro-ph.CO] .
  * Colombi (2021) S. Colombi, A&A 647, A66 (2021), arXiv:2012.04409 [astro-ph.CO] .
  * Saga  _et al._ (2022) S. Saga, A. Taruya, and S. Colombi, A&A 664, A3 (2022), arXiv:2111.08836 [astro-ph.CO] .
  * Angulo and Hahn (2022) R. E. Angulo and O. Hahn, Living Reviews in Computational Astrophysics 8, 1 (2022), arXiv:2112.05165 [astro-ph.CO] .
  * Widrow and Kaiser (1993) L. M. Widrow and N. Kaiser, ApJ 416, L71+ (1993).
  * Schive  _et al._ (2014) H.-Y. Schive, T. Chiueh, and T. Broadhurst, Nature Phys. 10, 496 (2014), arXiv:1406.6586 [astro-ph.GA] .
  * Kopp  _et al._ (2017) M. Kopp, K. Vattis, and C. Skordis, Phys. Rev. D 96, 123532 (2017), arXiv:1711.00140 [astro-ph.CO] .
  * Garny and Konstandin (2018) M. Garny and T. Konstandin, J. Cosmology Astropart. Phys 2018, 009 (2018), arXiv:1710.04846 [astro-ph.CO] .
  * Uhlemann (2018) C. Uhlemann, J. Cosmology Astropart. Phys 2018, 030 (2018), arXiv:1807.07274 [astro-ph.CO] .
  * Mocz  _et al._ (2018) P. Mocz, L. Lancaster, A. Fialkov, F. Becerra, and P.-H. Chavanis, Phys. Rev. D 97, 083519 (2018), arXiv:1801.03507 [astro-ph.CO] .
  * Li  _et al._ (2019) X. Li, L. Hui, and G. L. Bryan, Phys. Rev. D 99, 063509 (2019), arXiv:1810.01915 [astro-ph.CO] .
  * Uhlemann  _et al._ (2019) C. Uhlemann, C. Rampf, M. Gosenca, and O. Hahn, Phys. Rev. D 99, 083524 (2019), arXiv:1812.05633 [astro-ph.CO] .
  * Garny  _et al._ (2020) M. Garny, T. Konstandin, and H. Rubira, JCAP 04, 003 (2020), arXiv:1911.04505 [astro-ph.CO] .
  * May and Springel (2021) S. May and V. Springel, Mon. Not. Roy. Astron. Soc. 506, 2603 (2021), arXiv:2101.01828 [astro-ph.CO] .
  * Evans (1993) N. W. Evans, MNRAS 260, 191 (1993).
  * Evans (1994) N. W. Evans, MNRAS 267, 333 (1994).
  * Frieman and Scoccimarro (1994) J. Frieman and R. Scoccimarro, ApJ 431, L23 (1994), arXiv:astro-ph/9312043 .
  * Alcock  _et al._ (1995) C. Alcock, R. A. Allsman, T. S. Axelrod, D. P. Bennett, K. H. Cook, N. W. Evans, K. C. Freeman, K. Griest, J. Jijina, M. Lehner, S. L. Marshall, S. Perlmutter, B. A. Peterson, M. R. Pratt, P. J. Quinn, A. W. Rodgers, C. W. Stubbs, W. Sutherland, and MACHO Collaboration, ApJ 449, 28 (1995), arXiv:astro-ph/9411019 [astro-ph] .
  * Calcino  _et al._ (2018) J. Calcino, J. García-Bellido, and T. M. Davis, MNRAS 479, 2889 (2018), arXiv:1803.09205 [astro-ph.CO] .
  * Navarro  _et al._ (1997) J. Navarro, C. Frenk, and S. White, ApJ 490, 493 (1997).
  * Chua  _et al._ (2019) K. T. E. Chua, A. Pillepich, M. Vogelsberger, and L. Hernquist, MNRAS 484, 476 (2019), arXiv:1809.07255 [astro-ph.GA] .
  * Prada  _et al._ (2019) J. Prada, J. E. Forero-Romero, R. J. J. Grand, R. Pakmor, and V. Springel, MNRAS 490, 4877 (2019), arXiv:1910.04045 [astro-ph.GA] .
  * Binney and Tremaine (2008) J. Binney and S. Tremaine,  _Galactic Dynamics: Second Edition_ (2008).
  * Hansen and Moore (2006) S. H. Hansen and B. Moore, New A 11, 333 (2006), arXiv:astro-ph/0411473 [astro-ph] .
  * Wojtak  _et al._ (2008) R. Wojtak, E. L. Łokas, G. A. Mamon, S. Gottlöber, A. Klypin, and Y. Hoffman, MNRAS 388, 815 (2008), arXiv:0802.0429 [astro-ph] .
  * Sparre and Hansen (2012) M. Sparre and S. H. Hansen, J. Cosmology Astropart. Phys 2012, 049 (2012), arXiv:1210.2392 [astro-ph.CO] .
  * Osipkov (1979) L. P. Osipkov, Soviet Astronomy Letters 5, 42 (1979).
  * Merritt (1985) D. Merritt, AJ 90, 1027 (1985).
  * Cole and Lacey (1996) S. Cole and C. Lacey, MNRAS 281, 716 (1996), arXiv:astro-ph/9510147 [astro-ph] .
  * Sheth  _et al._ (2001) R. Sheth, L. Hui, A. Diaferio, and R. Scoccimarro, MNRAS 325, 1288 (2001).
  * Pichon and Bernardeau (1999) C. Pichon and F. Bernardeau, A&A 343, 663 (1999), arXiv:astro-ph/9902142 [astro-ph] .
  * Szapudi and Kaiser (2003) I. Szapudi and N. Kaiser, ApJ 583, L1 (2003), arXiv:astro-ph/0211065 .
  * Neyrinck  _et al._ (2009) M. C. Neyrinck, I. Szapudi, and A. S. Szalay, ApJ 698, L90 (2009), arXiv:0903.4693 [astro-ph.CO] .
  * Wang  _et al._ (2011) X. Wang, M. Neyrinck, I. Szapudi, A. Szalay, X. Chen, J. Lesgourgues, A. Riotto, and M. Sloth, Astrophys. J. 735, 32 (2011), arXiv:1103.2166 [astro-ph.CO] .
  * Carron and Szapudi (2013) J. Carron and I. Szapudi, MNRAS 434, 2961 (2013), arXiv:1306.1230 [astro-ph.CO] .
  * Rubira and Voivodic (2021) H. Rubira and R. Voivodic, JCAP 03, 070 (2021), arXiv:2011.12280 [astro-ph.CO] .
  * Scoccimarro (2001) R. Scoccimarro, ArXiv:astro-ph/0008277, Annals New York Academy Sciences 927, 13 (2001).
  * Goroff  _et al._ (1986) M. Goroff, B. Grinstein, S.-J. Rey, and M. Wise, ApJ 311, 6 (1986).
  * Bernardeau  _et al._ (2002b) F. Bernardeau, S. Colombi, E. Gaztañaga, and R. Scoccimarro, Phys. Rep. 367, 1 (2002b), astro-ph/0112551 .
  * Blas  _et al._ (2015) D. Blas, S. Floerchinger, M. Garny, N. Tetradis, and U. A. Wiedemann, JCAP 11, 049 (2015), arXiv:1507.06665 [astro-ph.CO] .
  * Garny and Taule (2021) M. Garny and P. Taule, JCAP 01, 020 (2021), arXiv:2008.00013 [astro-ph.CO] .
  * Abolhasani  _et al._ (2016) A. A. Abolhasani, M. Mirbabayi, and E. Pajer, JCAP 05, 063 (2016), arXiv:1509.07886 [hep-th] .
  * Peebles (1980) P. Peebles,  _The large-scale structure of the universe_ (Princeton University Press, 1980).
  * Eddington (1916) A. S. Eddington, MNRAS 76, 572 (1916).
  * Evans and Williams (2014) N. W. Evans and A. A. Williams, MNRAS 443, 791 (2014), arXiv:1406.3730 [astro-ph.GA] .