# Perturbation theory with dispersion and higher cumulants: non-linear regime

Mathias Garny<EMAIL_ADDRESS>Physik Department T31, Technische
Universität München, James-Franck-Straße 1, D-85748 Garching, Germany  Dominik
Laxhuber<EMAIL_ADDRESS>Physik Department T31, Technische
Universität München, James-Franck-Straße 1, D-85748 Garching, Germany  Román
Scoccimarro<EMAIL_ADDRESS>Center for Cosmology and Particle Physics,
Department of Physics, New York University, NY 10003, New York, USA

(October 14, 2022)

###### Abstract

We present non-linear solutions of Vlasov Perturbation Theory (), describing
gravitational clustering of collisionless dark matter with dispersion and
higher cumulants induced by orbit crossing. We show that can be cast into a
form that is formally analogous to standard perturbation theory (SPT), but
including additional perturbation variables, non-linear interactions, and a
more complex propagation. non-linear kernels have a crucial decoupling
property: for fixed total momentum, the kernels becomes strongly suppressed
when any of the individual momenta cross the dispersion scale into the non-
linear regime. This screening of UV modes allows us to compute non-linear
corrections to power spectra even for cosmologies with very blue power-law
input spectra, for which SPT diverges. We compare predictions for the density
and velocity divergence power spectra as well as the bispectrum at one-loop
order to N-body results in a scaling universe with spectral indices $-1\leq
n_{s}\leq+2$. We find a good agreement up to the non-linear scale for all
cases, with a reach that increases with the spectral index $n_{s}$. We discuss
the generation of vorticity as well as vector and tensor modes of the velocity
dispersion, showing that neglecting vorticity when including dispersion would
lead to a violation of momentum conservation. We verify momentum conservation
when including vorticity, and compute the vorticity power spectrum at two-loop
order, necessary to recover the correct large-scale limit with slope
$n_{w}=2$. Comparing to our N-body measurements confirms the cross-over from
$k^{4}$ to $k^{2}$ scaling on large scales. Our results provide a proof-of-
principle that perturbative techniques for dark matter clustering can be
systematically improved based on the known underlying collisionless dynamics.

††preprint: TUM-HEP 1424/22

###### Contents

  1. I Introduction
  2. II Non-linear kernels with dispersion and higher cumulants
    1. II.1 Setup
    2. II.2 Non-linear kernels
    3. II.3 Analytical results in the limit $\epsilon\to 0$
      1. II.3.1 First order kernels
      2. II.3.2 Second order kernels
      3. II.3.3 Wilson coefficients
    4. II.4 Numerical results
  3. III Symmetry constraints on non-linear kernels
    1. III.1 Galilean invariance
    2. III.2 Mass and momentum conservation
      1. III.2.1 Density contrast
      2. III.2.2 Velocity divergence
  4. IV Generation of vorticity, vector and tensor modes
    1. IV.1 Linear approximation
    2. IV.2 Vorticity generation at second order
      1. IV.2.1 Analytical results for $\epsilon\to 0$
      2. IV.2.2 Numerical results
    3. IV.3 Vorticity generation at third order
    4. IV.4 Generation of tensor modes
  5. V Background dispersion from simulations and perturbation theory
    1. V.1 Self-consistent solution in perturbation theory
    2. V.2 Dispersion from simulated halos
  6. VI  predictions vs Simulations
    1. VI.1 Density power spectrum
      1. VI.1.1 Impact of vorticity, vector and tensor modes
      2. VI.1.2 Impact of higher cumulants and loops
    2. VI.2 Velocity divergence and cross power spectrum
    3. VI.3 Bispectrum
    4. VI.4 Vorticity power spectrum
  7. VII Conclusions
  8. A Tadpoles
  9. B Time integrals
  10. C Recursion relation and Galilean invariance
  11. D Treatment of vorticity, vector and tensor modes
    1. D.1 Simplified algorithm for the one-loop power spectrum
    2. D.2 General case
  12. E Rescaling for power law universe
  13. F N-body simulations
  14. G The halo mass function
    1. G.1 Measuring the Mass Function
    2. G.2 From Isothermal to NFW halos
  15. H Measuring Velocity Divergence and Vorticity

##  I Introduction

Perturbative techniques play a central role in understanding the formation and
evolution of large-scale structures in the universe. The foundation for this
understanding rests on the gravitational clustering of collisionless cold dark
matter (CDM). As such it is fully described by the Vlasov equation, which
follows the evolution of the phase-space distribution function (DF) coupled to
gravity via the Poisson equation at scales below the Hubble radius. The
conventional approach of standard perturbation theory (SPT) is a truncation of
this description, in which only density (zero-th moment of the DF) and
velocity (first moment of the DF) fields are taken into account (see
Bernardeau  _et al._ (2002) for a review).

This truncation is motivated by the fact that CDM has negligible primordial
velocity dispersion (second cumulant of the DF). However, even if absent
initially, second and higher cumulants of the DF are generated by orbit
crossing, and they should play a role in the transition to the non-linear
regime. This is particularly so for initial conditions with significant small-
scale power, i.e. blue spectral indices, for which SPT loop corrections
describing this transition are UV divergent Makino  _et al._ (1992);
Scoccimarro and Frieman (1996a). These UV divergences stand in sharp contrast
to what is seen in N-body simulations Colombi  _et al._ (1996), showing that
for such initial conditions fluctuations in the non-linear regime are most
suppressed compared to linear evolution in SPT.

This motivates to systematically improve the perturbative description by
taking second and higher cumulants into account, which obey a coupled
hierarchy of equations of motion Pueblas and Scoccimarro (2009). Linearized
solutions to this hierarchy were developed in paper I Garny  _et al._ (2022),
showing a structure that is much richer than in SPT as a result of the
backreaction of orbit crossing on linear modes. This leads to a suppression on
small scales that saturates (for a given fixed wavenumber) when truncating the
cumulant expansion at sufficiently high order. This damping of linear UV modes
sets in at a new scale, the dispersion scale, determined by the background
value of the second cumulant, i.e. the velocity dispersion tensor. Physically,
this damping corresponds to the backreaction of halo formation at small scales
on linear modes.

The purpose of the present paper is to present non-linear solutions obtained
by expanding around this new linear theory. We call this new perturbative
expansion Vlasov Perturbation Theory or . We show that can be cast into a form
that is formally analogous to SPT, but including additional perturbation
variables, non-linear interactions, and a more complex propagation. This all
follows from the Vlasov-Poisson equation in a straightforward way, without
having to adhere to any ad hoc assumptions.

The final outcome is to replace the well known non-linear kernels of SPT
describing the perturbative solutions by new kernels that have a crucial
decoupling property: for fixed total momentum ${\bm{k}}$, the kernels become
strongly suppressed when any of the individual momenta cross the dispersion
scale into the non-linear regime. As we shall see, this is precisely what is
needed to bring the UV divergencies of SPT for blue spectra under control.
This screening of UV modes is also what is expected physically to happen in
the non-linear regime, where the predominant structures are fairly stable dark
matter halos. For CDM spectra, this regime is characterized by the “virial
turnover”, where the non-linear power grows with wavenumber less than in the
weakly non-linear regime. The physical reason for this turnover can, at least
approximately, be understood by the stable clustering picture Davis and
Peebles (1977); Peebles (1980), where pairwise velocities cancel the Hubble
flow.

It may be convenient to illustrate here the main ideas behind how accounts for
the screening mechanism of UV modes for the simplest case, the density power
spectrum (analogous results hold for other variables, e.g. velocity
divergence, and other statistics, e.g. the bispectrum). Consider scaling
universes, i.e. scale-free Gaussian initial conditions with initial power
spectra $P_{0}\propto k^{n_{s}}$ and $\Omega_{m}=1$. Then the power spectrum
in follows a perturbative expansion of the form,

${P(k,\eta)\over e^{2\eta}P_{0}}=\sum_{L=0}^{\infty}\Big{(}{k\over k_{\rm
nl}(\eta)}\Big{)}^{(n_{s}+3)L}\
H_{L}(k/k_{\sigma}(\eta),\bar{\omega},\ldots)\,,$ (1)

where $e^{\eta}$ is the scale factor, $L$ is the number of loops, $k_{\rm
nl}(\eta)$ is the usual non-linear scale, $k_{\sigma}(\eta)\equiv
1/\sqrt{\epsilon}$ is the dispersion scale obtained from the background value
of the velocity dispersion tensor $\epsilon_{ij}$ (suitably rescaled, see Eq.
II.1 below), $\bar{\omega}\equiv\omega/\epsilon^{2}$ is a constant that
characterizes the background value of the fourth cumulant $\omega$, and the
dots denote other (even) background values of the cumulants, describing the DF
non-Gaussianity111These are the normalized cumulants $\bar{\cal E}_{2n}$
discussed in paper I Garny  _et al._ (2022).. By self-similarity
$k_{\sigma}/k_{\rm nl}$ is a constant that depends only on $n_{s}$. The
dispersion scale is a new scale induced by orbit crossing, and as such is
absent in SPT.

Figure 1:  Integrand of the one-loop contribution to the density power
spectrum in and SPT, for a scaling universe with power-law input spectrum
$P_{0}\propto k^{n_{s}}$ with spectral indices $n_{s}=2,1,0,-1$ shown in the
four panels. In SPT the integrand grows as
$F_{3}(\bm{k},\bm{q},-\bm{q})P_{0}(q)d^{3}q\sim q^{n_{s}+1}d\ln q$ for large
loop wavenumber $q$, i.e. the well-known suppression $F_{3}\propto
k^{2}/q^{2}$ due to momentum conservation is not sufficient to render the one-
loop integral UV finite for all $n_{s}\geq-1$. In contrast, the physically
expected screening of UV modes within the full collisionless Vlasov-Poisson
dynamics is captured by , compensating the increase for large $q$ for blue
spectra with high $n_{s}$ and leading to a finite integral with asymptotic
decay of the integrand as
$F_{1,\delta}(k,\eta)F_{3,\delta}(\bm{k},\bm{q},-\bm{q},\eta)P_{0}(q)d^{3}q\sim
q^{-2}d\ln q$ independently of the spectral index $n_{s}$. The sum of linear
and one-loop contributions within matches well with N-body results (see Fig.
10, including details on the cumulant hierarchy truncation used for this
figure). Note that the $1/q^{2}$ scaling within is universal and occurs for
all truncations considered in this work, see Fig. 2 and Eq. (II.4) for
details.

The prefactors of $(k/k_{\rm nl})$ in Eq. (1) are the standard loop factors
that appear in SPT, while the functions $H_{L}$ describe the scale-dependence
due to dispersion and higher cumulants that give rise to the screening of UV
modes. For example, $H_{0}$ is the square of the linear kernel which is highly
suppressed when $k\gg k_{\sigma}$. As stressed in paper I Garny  _et al._
(2022), this suppression is physically different from a Jeans scale due to
pressure, or viscosity in a normal fluid. In fact, the linear theory of
collisionless matter is much richer than that of such fluids, with an
effective description in terms of density and velocity divergence that is non-
local in time, which is key to satisfy the cosmic energy equation Garny  _et
al._ (2022). The dispersion scale can in principle be measured independently
from simulations, and as we show in this paper it can be estimated from basic
properties of dark matter halos. Intuitively we know the dispersion scale is
close to the non-linear scale from the “fingers of god” features seen in the
redshift-space matter distribution.

The SPT limit of Eq. (1) is somewhat delicate. Naively, it is reached by
taking the $k_{\sigma}\to\infty$ limit (zero dispersion, including vanishing
higher cumulants). In that case, however, the functions $H_{L}$ for $L\geq 1$
only have a finite limit (giving $k$-independent finite functions of $n_{s}$)
when $n_{s}<-1$222For example, for the one-loop density power spectrum $H_{1}$
becomes the function $\alpha(n_{s})$ given by Eq. (58) and Fig. 9 in
Scoccimarro (1997).. This means that the limit of zero dispersion only exists
for spectra which are sufficiently red-tilted. For such spectra though, there
is still orbit crossing, $k_{\sigma}$ is still finite and corrections to SPT
are predicted.

In this paper we focus on spectra with $n_{s}\geq-1$, in order to
comprehensively assess the applicability of in a situation where the
limitations of SPT are most severe333For an effective field theory description
of these spectra at the one-loop level see Pajer and Zaldarriaga (2013).. We
compare one-loop predictions for the density power spectrum and bispectrum,
the velocity divergence power spectrum and the density-velocity divergence
power spectrum against N-body simulations. In addition, we discuss the
generation of vorticity, presenting solutions up to two-loop in order to
capture the correct scaling of the vorticity power spectrum on large scales.
We also highlight the importance of backreaction of vorticity modes on the
density contrast, both from theoretical considerations as well as numerical
results.

To illustrate the screening of UV modes captured by (as opposed to SPT) we
show in Fig. 1 the _integrand_ of the one-loop density power spectrum for blue
power-law input spectra $P_{0}\propto k^{n_{s}}$ with $n_{s}=2,1,0,-1$. The
normalization is such that $P_{\delta\delta}^{1L}(k,\eta)=\int d\ln
qP_{\delta\delta}^{\text{Integrand}}(k,q,\eta)$. This integrand is closely
related to the one-loop contribution to the response function studied in
Nishimichi  _et al._ (2016,
2017)444$P_{\delta\delta}^{\text{Integrand}}(k,q,\eta)=T^{1L}(k,q)qP_{0}(q,\eta)P_{0}(k,\eta)$
with response $T(k,q)$ in the notation of Nishimichi  _et al._ (2016).. In SPT
one would have

$P_{\delta\delta}^{\text{Integrand}}(k,q,\eta)\big{|}_{\text{SPT}}\propto
q^{n_{s}+1}\,,$ (2)

such that the loop integration is UV divergent for all $n_{s}\geq-1$. In
contrast, the prediction from the Vlasov theory is that the integrand drops as

$P_{\delta\delta}^{\text{Integrand}}(k,q,\eta)\big{|}_{\text{\rm}}\propto
q^{-2}\,,$ (3)

_independently of $n_{s}$_, rendering the integral well-behaved and implying
that perturbation modes with wavenumber far above the dispersion scale (which
is of order of the non-linear scale, as will be discussed in detail below)
give a negligible contribution, in accordance with physical expectations. We
emphasize that also the linear theory in is changed compared to SPT, and we
find that the sum of linear and integrated one-loop contributions matches well
with N-body results.

Previous work in the literature on corrections to SPT implied by the Vlasov
equation spans a number of different fronts. In Pueblas and Scoccimarro (2009)
the equation of motion for the Vlasov cumulant generating function hierarchy
was written down, explaining the conditions under which it matches solution of
the Vlasov equation after shell-crossing. They used measurements of the stress
tensor in simulations to close the hierarchy and estimate the backreaction of
shell-crossing on large-scale density and velocity divergence power spectra,
including a study on the growth of vorticity. In McDonald (2011) a low-$k$
expansion of the linearized hierarchy truncated at the second cumulant is used
to obtain the impact of velocity dispersion on density and velocity
divergence. In addition, an approximate treatment of non-linear evolution is
presented to show that the growth of the velocity dispersion background is
sustained by mode coupling. Still within the second-cumulant truncation,
Aviles (2016) and Erschfeld and Floerchinger (2019) discuss the growth of the
velocity dispersion tensor, with the latter work featuring the evolution of
the background dispersion at the non-linear level while keeping the
fluctuations in the dispersion linear (see also Erschfeld (2021)). Different
strategies for closing the Vlasov hierarchy are presented in Uhlemann (2018).
Complementary approaches from the Lagrangian picture include McDonald and Vlah
(2018) and Cusin  _et al._ (2017) that computes the growth of vorticity driven
by velocity dispersion.

This work is structured as follows. In Section II we discuss the basics of the
Vlasov cumulant hierarchy, splitting cumulants of the DF into their
expectation values and fluctuations and presenting their equations of motion.
We also introduce the non-linear kernels, and study analytically their scaling
in the limit of small wavenumber compared to the dispersion scale. In
addition, we study numerically the case when wavenumbers cross the dispersion
scale highlighting the screening of UV modes. Along the way we discuss the
properties of these kernels in detail for various truncations of the cumulant
hierarchy. Section III discusses the symmetry properties of the non-linear
kernels, including Galilean invariance and mass and momentum conservation. In
Section IV we discuss the generation of vorticity, and vector and tensor modes
of the velocity dispersion tensor. Section V computes the dispersion scale in
two different, complementary ways. First, using perturbation theory, and
second, using knowledge of dark matter profiles together with the halo mass
function measured in our N-body simulations. In Section VI we present the
predictions of against measurements in our simulations, including the density
power spectrum and bispectrum, velocity divergence power spectrum, the cross-
spectrum between density and velocity divergence, and the vorticity power
spectrum. Section VII presents our conclusions.

A number of appendices contain supplementary results that support the
discussion in the main text. Appendix A discusses our approach to tadpoles,
App. B presents analytic results for time integrations over the background
dispersion, and App. C proves the scaling of non-linear kernels in the
squeezed limit due to Galilean invariance. Appendix D discusses our
implementation of vorticity, vector and tensor modes in the numerical scheme,
and App. E shows how to rescale solutions to different values of the
dispersion scale. Appendix F presents our N-body simulations, App. G discusses
the determination of the halo mass function and details about the calculation
of the dispersion scale from halos, and App. H presents our algorithm for
measuring velocity divergence and vorticity in the simulations.

##  II Non-linear kernels with dispersion and higher cumulants

We follow the formalism for perturbation theory introduced in paper I Garny
_et al._ (2022), which presents linear solutions for the density and velocity
perturbations in the presence of a non-zero velocity dispersion tensor and
higher cumulants. We start from the fact that the phase-space distribution
function $f(\tau,\bm{x},\bm{p})$ of non-relativistic matter can be
characterized by its cumulants, obtained by expanding the cumulant generating
function in the auxiliary vector $\bm{l}$ Pueblas and Scoccimarro (2009);
Garny  _et al._ (2022),

$\displaystyle{\cal C}(\tau,\bm{x},\bm{l})$ $\displaystyle=$
$\displaystyle\ln\int d^{3}p\,e^{\bm{l}\cdot\bm{p}/a}\,f(\tau,\bm{x},\bm{p})$
(4) $\displaystyle=$
$\displaystyle\ln(1+\delta)+l_{i}v_{i}+\frac{1}{2}l_{i}l_{j}\sigma_{ij}+\frac{1}{6}l_{i}l_{j}l_{k}{\cal
C}_{ijk}$ $\displaystyle{}+\frac{1}{24}l_{i}l_{j}l_{k}l_{m}{\cal
C}_{ijkm}+\dots\,,$

where $\delta(\tau,\bm{x})$ is the density contrast, $v_{i}(\tau,\bm{x})$ the
peculiar velocity, $\sigma_{ij}(\tau,\bm{x})$ the velocity dispersion tensor,
and ${\cal C}_{ijk}({\cal C}_{ijkm})$ the third (fourth) cumulant.
Furthermore, $a$ is the scale-factor, $\tau$ the conformal time and $\bm{p}$
the comoving momentum per unit particle mass. In this section we first review
the relevant equations and notation, and then describe how to recursively
obtain non-linear kernels. These characterize the non-linear solutions,
generalizing the well known non-linear kernels of standard perturbation theory
(SPT) Fry (1984); Goroff  _et al._ (1986); Bernardeau  _et al._ (2002).

###  II.1 Setup

The Vlasov or collisionless Boltzmann equation for the phase-space
distribution in the non-relativistic, sub-horizon limit yields the following
equation of motion for the cumulant generating function Pueblas and
Scoccimarro (2009)

$\partial_{\tau}{\cal C}+{\cal H}(\bm{l}\cdot\nabla_{l}){\cal C}+(\nabla{\cal
C})\cdot(\nabla_{l}{\cal C})+(\nabla\cdot\nabla_{l}){\cal
C}=-\bm{l}\cdot\nabla\Phi\,,$ (5)

where $\Phi$ is the gravitational potential given by
$\nabla^{2}\Phi=\frac{3}{2}{\cal H}^{2}\Omega_{m}\delta$, with conformal
Hubble rate ${\cal H}=\partial_{\tau}\ln a$ and time-dependent matter density
parameter $\Omega_{m}$.

The cumulants satisfy a coupled, non-linear hierarchy of differential
equations, obtained by Taylor expanding Eq. (5) in $\bm{l}$, given by the
continuity and momentum conservation equation at zeroth and first order,

$\displaystyle\partial_{\tau}\delta+\nabla_{i}[(1+\delta)v_{i}]$
$\displaystyle=$ $\displaystyle 0\,,$ (6)
$\displaystyle\partial_{\tau}v_{i}+{\cal
H}v_{i}+v_{j}\nabla_{j}v_{i}+\nabla_{i}\Phi$ $\displaystyle=$
$\displaystyle-\nabla_{j}\sigma_{ij}-\sigma_{ij}\nabla_{j}\ln(1+\delta),$

with the latter containing the velocity dispersion tensor on the right-hand
side. Its equation of motion is obtained when expanding Eq. (5) to second
order in $\bm{l}$, and it is in turn coupled to the third cumulant, whose
equation follows from expanding to third order. They are given by

$\displaystyle\partial_{\tau}\sigma_{ij}$ $\displaystyle+$ $\displaystyle
2{\cal
H}\sigma_{ij}+v_{k}\nabla_{k}\sigma_{ij}+\sigma_{jk}\nabla_{k}v_{i}+\sigma_{ik}\nabla_{k}v_{j}$
$\displaystyle=$ $\displaystyle-\nabla_{k}{\cal C}_{ijk}-{\cal
C}_{ijk}\nabla_{k}\ln(1+\delta)\,,$ $\displaystyle\partial_{\tau}{\cal
C}_{ijk}$ $\displaystyle+$ $\displaystyle 3{\cal H}{\cal
C}_{ijk}+v_{m}\nabla_{m}{\cal C}_{ijk}$ (7) $\displaystyle+$
$\displaystyle\sigma_{km}\nabla_{m}\sigma_{ij}+\sigma_{im}\nabla_{m}\sigma_{kj}+\sigma_{jm}\nabla_{m}\sigma_{ik}$
$\displaystyle+$ $\displaystyle{\cal C}_{jkm}\nabla_{m}v_{i}+{\cal
C}_{ikm}\nabla_{m}v_{j}+{\cal C}_{ijm}\nabla_{m}v_{k}$ $\displaystyle=$
$\displaystyle-\nabla_{m}{\cal C}_{ijkm}-{\cal
C}_{ijkm}\nabla_{m}\ln(1+\delta)\,.$

The framework of SPT is based on neglecting the velocity dispersion in the
momentum conservation equation, which then becomes the standard Euler equation
for a perfect fluid with zero pressure. Here we include the dispersion tensor
and also higher cumulants, systematically taking non-linear terms in their
equations of motion into account within perturbation theory.

In order to be able to deal with the occurrence of the log-density field

$A\equiv\ln(1+\delta)\,,$ (8)

on the right-hand side of the evolution equations, we follow the hybrid scheme
introduced in paper I Garny  _et al._ (2022) and perturbatively solve for both
$\delta$ as well as $A$, complementing the set of evolutions equations by

$\partial_{\tau}A+\nabla_{i}v_{i}+v_{i}\nabla_{i}A=0\,,$ (9)

which follows directly from the continuity equation. It is convenient to
consider the rescaled cumulants,

$\displaystyle u_{i}$ $\displaystyle=$ $\displaystyle\frac{v_{i}}{-{\cal
H}f},\quad\epsilon_{ij}=\frac{\sigma_{ij}}{({\cal H}f)^{2}}\,,$
$\displaystyle\pi_{ijk}$ $\displaystyle=$ $\displaystyle\frac{{\cal
C}_{ijk}}{(-{\cal H}f)^{3}},\quad\Lambda_{ijkm}=\frac{{\cal C}_{ijkm}}{({\cal
H}f)^{4}}\,,$ (10)

where $f=d\ln D/d\ln a$ is the usual growth rate, with $D(a)$ being the
conventional growth factor. In addition, we switch from conformal time $\tau$
to $\eta=\ln(D)$ using $\partial_{\tau}={\cal H}f\partial_{\eta}$ and for any
quantity $X$ we have,

$\displaystyle(\partial_{\tau}+n{\cal H})[(-{\cal H}f)^{n}X]$ (11)
$\displaystyle=$ $\displaystyle-(-{\cal
H}f)^{n+1}\left(\partial_{\eta}+n\left(\frac{3}{2}\frac{\Omega_{m}}{f^{2}}-1\right)\right)X\,.$

Rotational invariance allows all even cumulants to possess a homogeneous
background value, that we parametrize for the second and fourth cumulant by
$\epsilon(\eta)$ and $\omega(\eta)$ via

$\displaystyle\langle\epsilon_{ij}(\eta,\bm{x})\rangle$ $\displaystyle=$
$\displaystyle\delta_{ij}\ \epsilon(\eta)\,,$
$\displaystyle\langle\Lambda_{ijkm}(\eta,\bm{x})\rangle$ $\displaystyle=$
$\displaystyle\left(\delta_{ij}\delta_{km}+2\,{\rm
cyc.}\right)\,\frac{\omega(\eta)}{5}\,.$ (12)

Equations of motion for the background values can be obtained by taking an
ensemble average of the evolution equations of the second and fourth cumulant,
respectively (see paper I Garny  _et al._ (2022)). For the second cumulant, it
is given by

$\partial_{\eta}\epsilon(\eta)+2\left(\frac{3}{2}\frac{\Omega_{m}}{f^{2}}-1\right)\epsilon(\eta)=Q(\eta)\,,$
(13)

with source term

$Q(\eta)=\frac{1}{3}\langle
u_{l}\nabla_{l}\epsilon_{ii}\rangle+\frac{2}{3}\langle\epsilon_{il}\nabla_{l}u_{i}\rangle+\frac{1}{3}\langle\pi_{iil}\nabla_{l}\ln(1+\delta)\rangle\,.$
(14)

Note that while $\langle\delta\rangle=0$ is ensured by mass conservation, the
zeroth cumulant $A=\ln(1+\delta)$ does also have a background value ${\cal
A}=\langle A\rangle$. Nevertheless, since $A$ enters only via
$\nabla_{i}\ln(1+\delta)$ in the equations of motion, the background value is
not needed explicitly Garny  _et al._ (2022).

To obtain equations for the perturbations we expand the cumulants around their
background values, in particular

$\epsilon_{ij}(\eta,\bm{x})=\delta_{ij}\,\epsilon(\eta)+\delta\epsilon_{ij}(\eta,\bm{x})\,.$
(15)

While $\delta,u_{i},{\cal C}_{ijk}$ have no background contribution, a similar
split is necessary for $A$ and analogously for $\Lambda_{ijkm}$. Following the
previous discussion, it is understood in the following that $A$ denotes the
perturbation part only, with ${\cal A}$ subtracted. The equation of motion for
$\delta\epsilon_{ij}$ is obtained by subtracting from Eq. (II.1) its ensemble
average, and similarly for $A$ and higher cumulants.

The velocity can be split into its divergence $\theta=\nabla_{i}u_{i}$ and
vorticity $w_{i}=\varepsilon_{ijk}\nabla_{j}u_{k}$ contributions via

$u_{i}=u_{i}^{S}+u_{i}^{V}=\frac{\nabla_{i}}{\nabla^{2}}\theta-\frac{\varepsilon_{ijk}\nabla_{j}}{\nabla^{2}}w_{k}\,,$
(16)

corresponding to one scalar and two vector modes, while $\delta\epsilon_{ij}$
consists of two scalar modes $g$ and $\delta\epsilon$, two vector modes
contained in the divergence-free vector $\nu_{i}$, and two tensor modes
represented by the symmetric, transverse-traceless matrix $t_{ij}$,

$\displaystyle\delta\epsilon_{ij}$ $\displaystyle=$
$\displaystyle\delta\epsilon_{ij}^{S}+\delta\epsilon_{ij}^{V}+\delta\epsilon_{ij}^{T}\,,$
(17)

where

$\displaystyle\delta\epsilon_{ij}^{S}$ $\displaystyle=$
$\displaystyle\delta_{ij}\,\delta\epsilon+\frac{\nabla_{i}\nabla_{j}}{\nabla^{2}}g\,,$
$\displaystyle\delta\epsilon_{ij}^{V}$ $\displaystyle=$
$\displaystyle-\frac{\varepsilon_{ilk}\nabla_{l}\nabla_{j}}{\nabla^{2}}\nu_{k}-\frac{\varepsilon_{jlk}\nabla_{l}\nabla_{i}}{\nabla^{2}}\nu_{k}\,,$
(18) $\displaystyle\delta\epsilon_{ij}^{T}$ $\displaystyle=$ $\displaystyle
t_{ij}\equiv P_{ij,ls}^{T}\,\delta\epsilon_{ls}\,,$

and the tensor projector $P_{ij,ls}^{T}$ is given in Garny  _et al._ (2022).

The impact of third, fourth and higher cumulants in the linear approximation
has been analyzed in detail in paper I Garny  _et al._ (2022). In this work,
we assess their impact at non-linear level by including the two scalar
perturbations $\pi$ and $\chi$ of the third cumulant,

$\pi_{ijk}^{S}=-\left(\delta_{ij}\frac{\nabla_{k}}{\nabla^{2}}+2\,{\rm
cyc.}\right)\frac{\chi}{5}-\frac{\nabla_{i}\nabla_{j}\nabla_{k}}{\nabla^{4}}(\pi-\chi)\,,$
(19)

as well as the background value $\omega(\eta)$ of the fourth cumulant. Note
that the latter enters the evolution equations of $\pi_{ijk}$ via the last
term in the last line of Eq. (II.1). We neglect non-scalar contributions to
$\pi_{ijk}$ as well as perturbation modes of the fourth cumulant, in order to
obtain a tractable approximation scheme. The impact of higher cumulants will
be quantified by comparing solutions with and without taking $\pi_{ijk}^{S}$
into account.

Apart from the complete set of perturbation modes for the zeroth, first and
second cumulant, it is important to note that our present setup includes all
higher cumulant perturbation modes relevant for the source term of the
background dispersion, given by

$Q(\eta)=\frac{1}{3}\int
d^{3}k\left(P_{\theta\tilde{g}}(k,\eta)+2P_{w_{i}\nu_{i}}(k,\eta)+P_{A\pi}(k,\eta)\right)\,,$
(20)

where $\tilde{g}\equiv g-\delta\epsilon$.

The scalar, vector and tensor perturbation modes can be collected in a
component vector

$\psi=(\psi^{S},\psi^{V},\psi^{T})\,,$ (21)

with

$\displaystyle\psi^{S}$ $\displaystyle=$
$\displaystyle(\delta,\theta,g,\delta\epsilon,A,\pi,\chi)\,,$
$\displaystyle\psi^{V}$ $\displaystyle=$ $\displaystyle(w_{i},\nu_{i})\,,$
$\displaystyle\psi^{T}$ $\displaystyle=$ $\displaystyle(t_{ij})\,.$ (22)

In Fourier space, the evolution equations can be written in the form

$\psi_{k,a}^{\prime}(\eta)+\Omega_{ab}(k,\eta)\psi_{k,b}(\eta)=\int_{pq}\gamma_{abc}(\bm{p},\bm{q})\psi_{p,b}(\eta)\psi_{q,c}(\eta)\,,$
(23)

where $\psi_{k,a}$ stands for the Fourier mode $\bm{k}$ of $\psi$, and $a$
runs over all perturbation modes, including their spatial indices in case of
vectors and tensors, with summation over repeated indices implied.
Furthermore, ${}^{\prime}=d/d\eta$ and $\int_{pq}=\int
d^{3}p\,d^{3}q\,\delta^{(3)}(\bm{k}-\bm{p}-\bm{q})$. The linear evolution is
characterized by the scale- and time-dependent matrix $\Omega_{ab}(k,\eta)$.
Due to rotational invariance it has a block-diagonal structure for scalar,
vector and tensor modes, with the blocks given by

$\displaystyle\Omega^{S}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{ccccccc}&-1\\\
-\frac{3}{2}d&\frac{1}{2}e&k^{2}&k^{2}&k^{2}\epsilon\\\
&-2\epsilon&e&&&1&-\frac{3}{5}\\\ &&&e&&&\frac{1}{5}\\\ &-1\\\
&&-3k^{2}\epsilon&-5k^{2}\epsilon&-k^{2}\omega&\frac{3}{2}e\\\
&&&-5k^{2}\epsilon&-k^{2}\omega&&\frac{3}{2}e\\\ \end{array}\right)\,,$ (31)
$\displaystyle\Omega^{V}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cc}\frac{1}{2}e&k^{2}\\\ -\epsilon&e\\\
\end{array}\right)\,,\quad\Omega^{T}=e\,,$ (34)

with $d(\eta)\equiv\Omega_{m}/f^{2}$ and $e(\eta)\equiv 3d(\eta)-2$, with
$d,e\mapsto 1$ in the EdS approximation. Note that the upper two-by-two part
of $\Omega^{S}$ corresponds to the well-known form within SPT. In addition,
each entry in $\Omega^{V}$ is understood to act on all three vector components
of $w_{i}$ and $\nu_{i}$, respectively, with the off-diagonal entries
describing a mixing of vorticity and vector modes of the dispersion tensor.
Similarly, $\Omega^{T}$ acts equally on all components of $t_{ij}$.

It is worth emphasizing that background values of the second and fourth
cumulant, $\epsilon(\eta)$ and $\omega(\eta)$, are not considered to be a
priori of any particular order in the perturbative expansion as they are
sourced by one-point moments (see Eq. 14) that correspond to integrals of
power spectra over all momenta (see Eq. 20). That is, these quantities are not
proportional to the fluctuation modes $\psi_{k,a}(\eta)$, and we treat them as
any function of time (e.g. $\Omega_{m}$), as we are explicitly agnostic about
the size of velocity dispersion effects in our approach. This is crucial and
allows us, as discussed in paper I Garny  _et al._ (2022), to obtain the
expected decoupling of large-scale modes from halo formation at small scales.
As we shall see below, this is described by the suppression of the non-linear
kernels at high momenta, generalizing the results for the linear kernel
described in detail in paper I Garny  _et al._ (2022). This reduces the
sensitivity of to UV modes when compared to SPT, helping with the convergence
of perturbative calculations. The decoupling effect appears for large
wavenumbers $k\gtrsim k_{\sigma}$ above the scale set by the background
dispersion

$k_{\sigma}(\eta)\equiv 1/\sqrt{\epsilon(\eta)}\,.$ (35)

As shown in paper I Garny  _et al._ (2022), the absence of exponential
instabilities of the linear system for $k\gg k_{\sigma}$ yields a restriction
on the size of the fourth relative to the second cumulant expectation value,
particularly for the dimensionless ratio

$\bar{\omega}\equiv\frac{\omega}{\epsilon^{2}}\,.$ (36)

For our setup considered here, the stability condition is given by

$-10\leq\bar{\omega}\leq 5\,,$ (37)

provided that $\alpha\equiv\partial_{\eta}\ln\epsilon$ and $\bar{\omega}$ are
constant or slowly varying. This is the case in the limit of a scaling
universe, where $\Omega_{m}=1$ and the initial perturbations are scale-free,
$P_{0}(k)\propto k^{n_{s}}$.

All non-linearities are captured by the _vertices_
$\gamma_{abc}(\bm{p},\bm{q})$ in Eq. (23). Their explicit form is given in
paper I Garny  _et al._ (2022), satisfying

$\gamma_{abc}(\bm{p},\bm{q})=\gamma_{acb}(\bm{q},\bm{p})\,,$ (38)

and we take all of them into account in our analysis555Up to the second
cumulant, all vertices obtained from the equations of motion are included. For
vertices involving at least one third cumulant perturbation $\pi$ or $\chi$,
we only include contributions where all three fields are scalar modes. We
checked that this is consistent with the symmetry constraints discussed in
Sec. III.. Note that the equations of motion include at most quadratic terms
when including both $\delta$ and $A$ in the set of perturbation variables.
They include the familiar SPT vertices $\gamma_{\delta\theta\delta}$,
$\gamma_{\delta\delta\theta}$ and $\gamma_{\theta\theta\theta}$, as well as a
large set of extra vertices involving second and third cumulant perturbations
and $A$. It is crucial to take the full set of vertices into account along
with the matrix $\Omega_{ab}$ to ensure that symmetry constraints on the
resulting perturbative solutions from Galilean symmetry as well as mass and
momentum conservation hold, as we explain in detail in Sec. III.

A final comment regarding Eq. (23) is worth pointing out. When deriving the
equation for the perturbation part by subtracting the full from the ensemble
averaged equations of motion, one obtains an equation that contains an
additional contribution compared to Eq. (23). It is proportional to a three-
dimensional Dirac delta supported at $\bm{k}=0$, and contributes only to the
components that possess also a background value, being $\delta\epsilon$ and
$A$ in the setup adopted here. We refer to App. A for a discussion of these
terms, that ensure the cancellation of so-called “tadpole” diagrams in the
perturbative solution. As shown there, these terms can be traded for a
modification of the vertices, that consists in setting
$\gamma_{abc}(\bm{p},\bm{q})$ to zero for $\bm{p}=-\bm{q}$. We adopt this
choice in the following, and use Eq. (23) to generate perturbative solutions.

###  II.2 Non-linear kernels

We consider perturbative solutions of the equation of motion, Eq. (23),
assuming that at some “initial” time $\eta_{\text{ini}}$ long before the onset
of non-linearity, but long after recombination, all perturbation modes are
proportional to the initial density field $\delta_{k0}$. This is the case for
adiabatic initial conditions and when assuming an initially negligibly small
(but non-zero, see below) velocity dispersion, as appropriate for cold dark
matter. The vector $\psi$ of perturbations may then be formally Taylor
expanded in powers of $\delta_{k0}$ as

$\psi_{k,a}(\eta)=\sum_{n\geq
1}\int_{k_{i}}F_{n,a}(\bm{k}_{1},\dots,\bm{k}_{n},\eta)\,e^{n\eta}\,\delta_{k_{1}0}\cdots\delta_{k_{n}0}\,,$
(39)

where $F_{n,a}$ is the $n$th order non-linear kernel for component $a$, and
$\int_{k_{i}}\equiv\int d^{3}k_{1}\cdots
d^{3}k_{n}\,\delta^{(3)}(\bm{k}-\sum_{i}\bm{k}_{i})$, with the constraint on
the sum over all wavevectors arising from spatial translation symmetry of the
equations of motion. Even if initially all higher cumulant _perturbations_ are
taken to be zero, they are generated by time-evolution in presence of a
background dispersion (and background values of higher cumulants in general).
As discussed in paper I Garny  _et al._ (2022), the perturbations in turn
source the background values, and we come back to the possibility of
determining self-consistent solutions of this system in Sec. V.1. For the
moment, we consider the background values $\epsilon(\eta)$ and $\omega(\eta)$
entering via Eq. (31) as given and discuss the implications for the
perturbations $\psi_{k,a}(\eta)$.

As usual, the $F_{n,a}$ can be taken to be symmetric under arbitrary exchanges
of wavenumbers in their arguments, which we assume in the following. The
normalization is chosen such that, when ignoring second and higher cumulants
and using the EdS approximation, $F_{n,\delta}$ and $F_{n,\theta}$ would
become time-independent and coincide with the conventional EdS-SPT kernels
$F_{n}$ and $G_{n}$, with $F_{1}=G_{1}=1$. When taking second and higher
cumulants into account, the kernels do depend on time, and even the linear
solutions described by $F_{1,a}(k,\eta)$ become non-trivial, featuring a
power-like suppression for $k\gg k_{\sigma}$, see paper I Garny  _et al._
(2022).

The expansion in Eq. (39) can be used to perturbatively compute the
(cross-)power spectra

$\langle\psi_{k,a}(\eta)\psi_{k^{\prime},b}(\eta)\rangle=\delta^{(3)}(\bm{k}+\bm{k}^{\prime})P_{ab}(k,\eta)\,,$
(40)

where the subscript $ab$ labels the various modes, with for example
$a=b=\delta$ corresponding to the density power spectrum, and
$a=\delta,b=\theta$ the density/velocity divergence cross spectrum. Following
the standard scenario, the initial density field $\delta_{k0}$ can be taken to
be drawn from a Gaussian random field fully characterized by its initial power
spectrum
$\langle\delta_{k0}\,\delta_{k^{\prime}0}\rangle=\delta^{(3)}(\bm{k}+\bm{k}^{\prime})P_{0}(k)$.
By inserting Eq. (39) into Eq. (40) and using the Wick theorem, we obtain a
loop expansion in close analogy to SPT,

$P_{ab}(k,\eta)=P_{ab}^{\text{lin}}(k,\eta)+P_{ab}^{1L}(k,\eta)+P_{ab}^{2L}(k,\eta)+\dots\,,$
(41)

but with modified non-linear kernels. For example, the linear and one-loop
contributions are given by

$\displaystyle P_{ab}^{\text{lin}}(k,\eta)$ $\displaystyle=$ $\displaystyle
e^{2\eta}F_{1,a}(k,\eta)F_{1,b}(k,\eta)P_{0}(k)\,,$ $\displaystyle
P_{ab}^{1L}(k,\eta)$ $\displaystyle=$ $\displaystyle e^{4\eta}\int
d^{3}q\Big{\\{}2F_{2,a}(\bm{k}-\bm{q},\bm{q},\eta)F_{2,b}(\bm{k}-\bm{q},\bm{q},\eta)$
(42) $\displaystyle{}\times P_{0}(|\bm{k}-\bm{q}|)P_{0}(q)$ $\displaystyle+$
$\displaystyle
3F_{1,a}(k,\eta)F_{3,b}(\bm{k},\bm{q},-\bm{q},\eta)P_{0}(k)P_{0}(q)$
$\displaystyle+$ $\displaystyle
3F_{3,a}(\bm{k},\bm{q},-\bm{q},\eta)F_{1,b}(k,\eta)P_{0}(k)P_{0}(q)\Big{\\}}\,.$

The expressions can be used to numerically compute loop corrections, with
$P_{0}(k)$ being the usual linear matter power spectrum obtained from
Boltzmann solvers such as CLASS Blas  _et al._ (2011) or CAMB Lewis  _et al._
(2000). As in SPT, the integration region with $|\bm{q}|\ll k\equiv|\bm{k}|$
yields large contributions to the individual summands, that cancel in their
sum (see Sec. III.1). It is therefore advantageous to ensure their
cancellation at the integrand level by multiplying the first summand in the
integrand of $P_{ab}^{1L}$ with $2\Theta(|\bm{k}-\bm{q}|-|\bm{q}|)$ and then
symmetrizing it with respect to the substitution $\bm{q}\to-\bm{q}$ Blas  _et
al._ (2013). Here $\Theta(x)$ is the Heaviside function. This operation leaves
the integral unchanged, but is preferable for numerical Monte Carlo
integration. The two-loop contribution to the power spectrum can be computed
analogously (see Sec. VI.4 and Blas  _et al._ (2013, 2013); Garny and Taule
(2021)). Similarly, the bispectrum can be computed using the modified non-
linear kernels and employing the same algorithm as in SPT to ensure the
cancellation of large contributions for low loop wavenumber at the integrand
level (see Sec. VI.3 and Floerchinger  _et al._ (2019)).

The difference in compared to SPT enters via the non-linear kernels $F_{n,a}$.
Equations of motion for the kernels can be obtained by inserting Eq. (39) into
Eq. (23) and collecting all terms with a given power of initial density
fields. This yields a set of coupled, ordinary differential equations given by

$\displaystyle(\partial_{\eta}\delta_{ab}+n\delta_{ab}+\Omega_{ab}(k,\eta))F_{n,b}(\bm{k}_{1},\dots,\bm{k}_{n},\eta)$
$\displaystyle=$
$\displaystyle\sum_{m=1}^{n-1}\Big{\\{}\gamma_{abc}(\bm{q}_{1}+\cdots+\bm{q}_{m},\bm{q}_{m+1}+\cdots+\bm{q}_{n})$
$\displaystyle\times
F_{m,b}(\bm{q}_{1},\dots,\bm{q}_{m},\eta)F_{n-m,c}(\bm{q}_{m+1},\dots,\bm{q}_{n},\eta)\Big{\\}}^{s}\,,$

where $\\{\cdots\\}^{s}=\sum_{\text{perm}}\\{\cdots\\}/|\text{perm}|$ denotes
an average over all $|\text{perm}|=n!/m!/(n-m)!$ possibilities to choose the
subset of wavevectors $\\{\bm{q}_{1},\dots,\bm{q}_{m}\\}$ from
$\\{\bm{k}_{1},\dots,\bm{k}_{n}\\}$, and $k\equiv|\sum_{i}\bm{k}_{i}|$. When
neglecting second and higher cumulants (i.e. restricting to
$a,b=\delta,\theta$) and in the EdS approximation these equations have time-
independent solutions for $F_{n,\delta}$ and $F_{n,\theta}$. Indeed, by
combining the equations for these two kernels in that limit and using time-
independence one recovers the well-known algebraic recursion relations Goroff
_et al._ (1986).

When taking second and higher cumulants into account, one obtains a coupled
system for all components of $F_{n,a}$, with $\Omega_{ab}(k,\eta)$ being
scale- and time-dependent. In this case, solutions can in general only be
found numerically by integrating the set of coupled equations, and for this
purpose we follow the strategy developed in Blas  _et al._ (2015); Garny and
Taule (2021). Nevertheless, it is instructive to consider analytical results
in the limit when the impact of velocity dispersion and higher cumulants is
small (but non-zero), which we discuss before presenting our numerical
results. In the following we set $\Omega_{m}/f^{2}\mapsto 1$ for simplicity.

###  II.3 Analytical results in the limit $\epsilon\to 0$

For wavenumbers that satisfy $\epsilon k_{i}^{2}\ll 1$ it is possible to find
approximate analytical solutions by solving for the non-linear kernels
perturbatively in powers of the background dispersion $\epsilon$. Using the
counting $\omega\propto\epsilon^{2}$ the kernels $F_{n,a}$ for
$a=\delta,\theta,A$ start at order $\epsilon^{0}$, the kernels for
$a=w_{i},g,\delta\epsilon,\nu_{i},t_{ij}$ at order $\epsilon^{1}$, and for
$a=\pi,\chi$ at order $\epsilon^{2}$, in accordance with the scaling derived
in paper I Garny  _et al._ (2022). The kernels
$F_{n,\delta}|_{\epsilon^{0}}=F_{n}$ and $F_{n,\theta}|_{\epsilon^{0}}=G_{n}$
coincide with the EdS-SPT kernels at lowest order in $\epsilon$, and for all
$n$.

####  II.3.1 First order kernels

Let us first derive the expansion of the linear kernels $F_{1,a}$ in
$\epsilon$, for which the right-hand side of Eq. (II.2) is zero. At linear
order one has $F_{1,A}=F_{1,\delta}$ and
$F_{1,w_{i}}=F_{1,\nu_{i}}=F_{1,t_{ij}}=0$, see Sec. IV for more details on
vorticity, vector and tensor modes. At lowest order in the expansion in
$\epsilon$ one recovers the EdS-SPT result
$F_{1,\delta}|_{\epsilon^{0}}=F_{1,\theta}|_{\epsilon^{0}}=1$. The evolution
equation for $g$ related to the second line of Eq. (31) yields

$\displaystyle(\partial_{\eta}+2)F_{1,g}(k,\eta)$ $\displaystyle=$
$\displaystyle 2\epsilon(\eta)F_{1,\theta}(k,\eta)-F_{1,\pi}(k,\eta)$ (44)
$\displaystyle{}+\frac{3}{5}F_{1,\chi}(k,\eta)\,.$

Using that the third cumulant modes $\pi,\chi$ contribute only starting at
order $\epsilon^{2}$ and inserting the lowest order expression for
$F_{1,\theta}$ yields

$F_{1,g}(k,\eta)=2E_{2}(\eta)+{\cal O}(\epsilon^{2})\,,$ (45)

where we define the weighted time integral of $\epsilon(\eta)$,

$E_{m}(\eta)\equiv\int^{\eta}d\eta^{\prime}\,e^{m(\eta^{\prime}-\eta)}\epsilon(\eta^{\prime})\,,$
(46)

which e.g. for $\epsilon=\epsilon_{0}\,e^{\alpha\eta}$ gives
$E_{m}(\eta)=\epsilon(\eta)/(m+\alpha)$.

Similarly one finds $F_{1,\delta\epsilon}={\cal O}(\epsilon^{2})$, i.e. the
linear kernel for $\delta\epsilon$ starts only at order $\epsilon^{2}$. Note
that the $\epsilon^{2}$ contribution is generated only by the third cumulant.
Indeed, when neglecting third and higher cumulants $F_{1,\delta\epsilon}$
would be zero at all orders in $\epsilon$, because the evolution equation for
$\delta\epsilon$ becomes trivial in that limit and corresponds to a decaying
mode. However, below we shall see that starting from second order, the kernels
for $\delta\epsilon$ have contributions at order $\epsilon^{1}$, as expected
from the general scaling discussed above.

The evolution equations for $\delta$ and $\theta$ can be written as

$\left(\partial_{\eta}+{\bf
1}+\Omega_{\text{SPT}}\right)\left(F_{1,\delta}\atop
F_{1,\theta}\right)=-k^{2}\left(0\atop F_{1,g}+F_{1,\delta\epsilon}+\epsilon
F_{1,A}\right)\,,$ (47)

where

$\Omega_{\text{SPT}}=\left(\begin{array}[]{cc}0&-1\\\
-\frac{3}{2}&\frac{1}{2}\end{array}\right)\,,$ (48)

is the standard EdS-SPT evolution matrix, ${\bf 1}$ denotes the unit matrix,
and the right-hand side contains the impact of velocity dispersion. We can
integrate out this equation to write the formal solution as

$\left(F_{1,\delta}\atop F_{1,\theta}\right)=\left(1\atop
1\right)+\int^{\eta}d\eta^{\prime}\,e^{\eta^{\prime}-\eta}g_{\text{SPT}}(\eta-\eta^{\prime})\left(0\atop{\cal
S}(k,\eta^{\prime})\right)\,,$ (49)

where we used SPT initial conditions in the far past (which specifies the
lower limit of integration), and $g_{\text{SPT}}(\eta-\eta^{\prime})$ is the
standard EdS-SPT linear Green function Scoccimarro (1998a) with ${\cal
S}(k,\eta)\equiv-k^{2}(F_{1,g}+F_{1,\delta\epsilon}+\epsilon F_{1,A})$. We
obtain the ${\cal O}(\epsilon)$ corrections to $F_{1,\delta}$ and
$F_{1,\theta}$ by inserting the lowest order results for ${\cal S}$, giving

$\displaystyle F_{1,\delta}(k,\eta)$ $\displaystyle=$ $\displaystyle
1-k^{2}I_{\delta}(\eta)+{\cal O}(\epsilon^{2})\,,$ $\displaystyle
F_{1,\theta}(k,\eta)$ $\displaystyle=$ $\displaystyle
1-k^{2}I_{\theta}(\eta)+{\cal O}(\epsilon^{2})\,,$ (50)

where

$\displaystyle I_{\delta}(\eta)$ $\displaystyle\equiv$
$\displaystyle\frac{2}{5}\int^{\eta}d\eta^{\prime}\,\left(1-e^{5(\eta^{\prime}-\eta)/2}\right)$
$\displaystyle{}\left(\epsilon(\eta^{\prime})+2\int^{\eta^{\prime}}d\eta^{\prime\prime}\,e^{2(\eta^{\prime\prime}-\eta^{\prime})}\epsilon(\eta^{\prime\prime})\right)\,,$
$\displaystyle I_{\theta}(\eta)$ $\displaystyle\equiv$
$\displaystyle\frac{2}{5}\int^{\eta}d\eta^{\prime}\,\left(1+\frac{3}{2}e^{5(\eta^{\prime}-\eta)/2}\right)$
(51)
$\displaystyle{}\left(\epsilon(\eta^{\prime})+2\int^{\eta^{\prime}}d\eta^{\prime\prime}\,e^{2(\eta^{\prime\prime}-\eta^{\prime})}\epsilon(\eta^{\prime\prime})\right)\,.$

Since $\epsilon\geq 0$ and $\eta^{\prime}\leq\eta$ we observe that the
correction term has a negative sign irrespective of the precise time-
dependence of $\epsilon(\eta)$, such that velocity dispersion necessarily
leads to a _suppression_ relative to SPT at first order in $\epsilon$. The
integrals can be simplified using the relation

$\displaystyle\int^{\eta}d\eta^{\prime}\,e^{a(\eta^{\prime}-\eta)}\int^{\eta^{\prime}}d\eta^{\prime\prime}\,e^{b(\eta^{\prime\prime}-\eta^{\prime})}\,f(\eta^{\prime\prime})$
$\displaystyle=-\frac{1}{a-b}\int^{\eta}d\eta^{\prime}\,(e^{a(\eta^{\prime}-\eta)}-e^{b(\eta^{\prime}-\eta)})\,f(\eta^{\prime})\,,$
(52)

for $a\not=b$. Using the definition of Eq. (46) of the time-weighted
background dispersion one finds,

$\displaystyle I_{\delta}(\eta)$ $\displaystyle=$
$\displaystyle\frac{4}{5}E_{0}(\eta)-2E_{2}(\eta)+\frac{6}{5}E_{5/2}(\eta)\,,$
$\displaystyle I_{\theta}(\eta)$ $\displaystyle=$
$\displaystyle\frac{4}{5}E_{0}(\eta)+2E_{2}(\eta)-\frac{9}{5}E_{5/2}(\eta)\,.$
(53)

Proceeding similarly for the third cumulant perturbations gives

$\displaystyle F_{1,\pi}$ $\displaystyle=$ $\displaystyle
F_{1,\chi}+3k^{2}\int^{\eta}d\eta^{\prime}\,e^{5(\eta^{\prime}-\eta)/2}\epsilon(\eta^{\prime})F_{1,g}(k,\eta^{\prime})\,,$
$\displaystyle F_{1,\chi}$ $\displaystyle=$ $\displaystyle
k^{2}\int^{\eta}d\eta^{\prime}\,e^{5(\eta^{\prime}-\eta)/2}(5\epsilon(\eta^{\prime})F_{1,\delta\epsilon}(k,\eta^{\prime})$
(54) $\displaystyle{}+\omega(\eta^{\prime})F_{1,A}(k,\eta^{\prime}))\,,$

which readily yields the leading ${\cal O}(\epsilon^{2})$ contribution when
inserting the lowest order expressions $F_{1,A}|_{\epsilon^{0}}=1$,
$F_{1,\delta\epsilon}|_{\epsilon^{1}}=0$ and $F_{1,g}|_{\epsilon^{1}}$ from
Eq. (45) on the right-hand side and recalling the power counting
$\omega\propto\epsilon^{2}$.

####  II.3.2 Second order kernels

At second order in perturbation theory we can proceed analogously, taking the
vertices on the right-hand side of Eq. (II.2) into account. For example, the
differential equation for the $g$ mode reads

$\displaystyle(\partial_{\eta}+3)F_{2,g}(\bm{p},\bm{q},\eta)=2\epsilon(\eta)F_{2,\theta}(\bm{p},\bm{q},\eta)-F_{2,\pi}(\bm{p},\bm{q},\eta)$
$\displaystyle+\frac{3}{5}F_{2,\chi}(\bm{p},\bm{q},\eta)+\gamma_{gbc}(\bm{p},\bm{q})F_{1,b}(p,\eta)F_{1,c}(q,\eta)\,,$

where summation over all $b,c$ is implied, with contributions coming in
general from $\gamma_{g\theta g}$, $\gamma_{g\theta\epsilon}$,
$\gamma_{gA\pi}$, $\gamma_{gA\chi}$, $\gamma_{gw_{i}g}$,
$\gamma_{gw_{i}\epsilon}$, $\gamma_{g\theta\nu_{i}}$, $\gamma_{g\theta
t_{ij}}$, $\gamma_{gw_{i}\nu_{j}}$, $\gamma_{gw_{i}t_{jk}}$, as well as
corresponding contributions with the second two entries in the subscript
flipped. All of these vertices are elementary functions of their arguments and
are given explicitly in paper I Garny  _et al._ (2022). Note that $\epsilon$
in the subscript stands for the mode $\delta\epsilon$. At second order, none
of the vertices involving vorticity, vector or tensor modes contribute since
$F_{1,a}=0$ for them. In addition, for the leading (linear) contribution in
the expansion in $\epsilon$ we can omit third cumulant contributions and also
those involving $\delta\epsilon$, since $F_{1,\delta\epsilon}$ starts at
$\epsilon^{2}$. This means only

$\gamma_{g\theta
g}(\bm{p},\bm{q})=\frac{1}{2}\frac{\bm{p}\cdot\bm{q}}{q^{2}}\left(\frac{(\bm{p}+\bm{q})^{2}}{p^{2}}+\frac{1}{2}-\frac{3}{2}\frac{((\bm{p}+\bm{q})\cdot\bm{p})^{2}}{(\bm{p}+\bm{q})^{2}p^{2}}\right)\,,$
(56)

and $\gamma_{gg\theta}(\bm{p},\bm{q})=\gamma_{g\theta g}(\bm{q},\bm{p})$ need
to be kept at ${\cal O}(\epsilon)$. Using in addition that
$F_{2,\theta}(\bm{p},\bm{q},\eta)|_{\epsilon^{0}}=G_{2}(\bm{p},\bm{q})$ is
equal to the usual EdS-SPT second order velocity kernel and
$F_{1,g}|_{\epsilon^{1}}$ from (45) yields

$\displaystyle F_{2,g}(\bm{p},\bm{q},\eta)$ $\displaystyle=$ $\displaystyle
2G_{2}(\bm{p},\bm{q})E_{3}(\eta)+2(\gamma_{g\theta g}(\bm{p},\bm{q})$ (57)
$\displaystyle{}+\gamma_{g\theta g}(\bm{q},\bm{p}))(E_{2}(\eta)-E_{3}(\eta))$
$\displaystyle{}+{\cal O}(\epsilon^{2})\,.$

Analogously one finds

$\displaystyle F_{2,\delta\epsilon}(\bm{p},\bm{q},\eta)$ $\displaystyle=$
$\displaystyle 2(\gamma_{\epsilon\theta
g}(\bm{p},\bm{q})+\gamma_{\epsilon\theta
g}(\bm{q},\bm{p}))(E_{2}(\eta)-E_{3}(\eta))$ (58) $\displaystyle{}+{\cal
O}(\epsilon^{2})\,,$

where Garny  _et al._ (2022)

$\gamma_{\epsilon\theta
g}(\bm{p},\bm{q})=\frac{1}{2}\frac{\bm{p}\cdot\bm{q}}{2p^{2}q^{2}(\bm{p}+\bm{q})^{2}}((\bm{p}\cdot\bm{q})^{2}-p^{2}q^{2})\,.$
(59)

Note that, while $F_{1,\delta\epsilon}$ starts only at order $\epsilon^{2}$,
the second-order kernel $F_{2,\delta\epsilon}$ has a non-zero contribution
already at (leading) order $\epsilon^{1}$. The same is true for higher-order
kernels of $\delta\epsilon$.

For the density and velocity divergence kernels we use

$\left(\partial_{\eta}+2\cdot{\bf
1}+\Omega_{\text{SPT}}\right)\left(F_{2,\delta}\atop
F_{2,\theta}\right)=\left(\gamma_{\delta bc}F_{1,b}F_{1,c}\atop\gamma_{\theta
bc}F_{1,b}F_{1,c}+{\cal S}_{2}\right)\,,$ (60)

where we omit the arguments for brevity and set ${\cal
S}_{2}\equiv-(\bm{p}+\bm{q})^{2}(F_{2,g}+F_{2,\delta\epsilon}+\epsilon
F_{2,A})$. Possible vertices are $\gamma_{\delta\theta\delta},\gamma_{\delta
w_{i}\delta}$ and $\gamma_{\theta\theta\theta}$, $\gamma_{\theta Ag}$,
$\gamma_{\theta A\epsilon}$, $\gamma_{\theta w_{i}\theta}$, $\gamma_{\theta
w_{i}w_{i}}$, $\gamma_{\theta A\nu_{i}}$, $\gamma_{\theta At_{ij}}$ given in
paper I Garny  _et al._ (2022) along with the flipped contributions. As
before, only scalar vertices contribute to $F_{2,\delta}$ and $F_{2,\theta}$,
while contributions involving $F_{1,\delta\epsilon}$ matter only at order
$\epsilon^{2}$. Up to first order in $\epsilon$ the relevant vertices
therefore amount to the two SPT vertices
$\gamma_{\delta\theta\delta}(\bm{p},\bm{q})=\alpha(\bm{p},\bm{q})/2=(\bm{p}+\bm{q})\cdot\bm{p}/(2p^{2})$
and
$\gamma_{\theta\theta\theta}(\bm{p},\bm{q})=\beta(\bm{p},\bm{q})=(\bm{p}+\bm{q})^{2}\bm{p}\cdot\bm{q}/(2p^{2}q^{2})$
as well as

$\gamma_{\theta
Ag}=-((\bm{p}+\bm{q})\cdot\bm{q})(\bm{p}\cdot\bm{q})/(2q^{2})\,.$ (61)

Writing the equation in integral form

$\displaystyle\left(F_{2,\delta}\atop F_{2,\theta}\right)$ $\displaystyle=$
$\displaystyle\int^{\eta}d\eta^{\prime}\,e^{2(\eta^{\prime}-\eta)}g_{\text{SPT}}(\eta-\eta^{\prime})$
(62) $\displaystyle{}\times\left(\gamma_{\delta
bc}F_{1,b}F_{1,c}\atop\gamma_{\theta bc}F_{1,b}F_{1,c}+{\cal
S}_{2}\right)\Bigg{|}_{\eta^{\prime}}\,,$

and expanding the right-hand side up to linear order in $\epsilon$ yields

$\displaystyle F_{2,\delta}(\bm{p},\bm{q},\eta)$ $\displaystyle=$
$\displaystyle
F_{2}(\bm{p},\bm{q})-(\bm{p}+\bm{q})^{2}\sum_{j=1}^{7}\Gamma_{j}(\bm{p},\bm{q})J^{\delta}_{j}(\eta)$
$\displaystyle{}+{\cal O}(\epsilon^{2})\,,$ $\displaystyle
F_{2,\theta}(\bm{p},\bm{q},\eta)$ $\displaystyle=$ $\displaystyle
G_{2}(\bm{p},\bm{q})-(\bm{p}+\bm{q})^{2}\sum_{j=1}^{7}\Gamma_{j}(\bm{p},\bm{q})J^{\theta}_{j}(\eta)$
(63) $\displaystyle{}+{\cal O}(\epsilon^{2})\,,$

where the time-dependence at linear order in $\epsilon$ is given by the
integrals $J_{j}^{\delta/\theta}(\eta)$, with explicit form given in App. B,
and the dependence on wavenumber arising from the various non-linear vertices
can be written as

$\displaystyle\Gamma_{1}(\bm{p},\bm{q})$ $\displaystyle\equiv$ $\displaystyle
G_{2}(\bm{p},\bm{q})\,,$ $\displaystyle\Gamma_{2}(\bm{p},\bm{q})$
$\displaystyle\equiv$ $\displaystyle\gamma_{g\theta
g}(\bm{p},\bm{q})+\gamma_{g\theta g}(\bm{q},\bm{p})+\gamma_{\epsilon\theta
g}(\bm{p},\bm{q})$ $\displaystyle{}+\gamma_{\epsilon\theta
g}(\bm{q},\bm{p})\,,$ $\displaystyle\Gamma_{3}(\bm{p},\bm{q})$
$\displaystyle\equiv$ $\displaystyle
F_{2,A}(\bm{p},\bm{q},\eta)|_{\epsilon^{0}}$ $\displaystyle=$
$\displaystyle\left[G_{2}(\bm{p},\bm{q})+\gamma_{A\theta
A}(\bm{p},\bm{q})+\gamma_{A\theta A}(\bm{q},\bm{p})\right]/2\,,$
$\displaystyle\Gamma_{4}(\bm{p},\bm{q})$ $\displaystyle\equiv$
$\displaystyle(p^{2}+q^{2})\gamma_{\theta\theta\theta}(\bm{p},\bm{q})/(\bm{p}+\bm{q})^{2}\,,$
$\displaystyle\Gamma_{5}(\bm{p},\bm{q})$ $\displaystyle\equiv$
$\displaystyle-(\gamma_{\theta Ag}(\bm{p},\bm{q})+\gamma_{\theta
Ag}(\bm{q},\bm{p}))/(\bm{p}+\bm{q})^{2}\,,$
$\displaystyle\Gamma_{6}(\bm{p},\bm{q})$ $\displaystyle\equiv$
$\displaystyle(p^{2}\gamma_{\delta\theta\delta}(\bm{p},\bm{q})+q^{2}\gamma_{\delta\theta\delta}(\bm{q},\bm{p}))/(\bm{p}+\bm{q})^{2}$
$\displaystyle=$ $\displaystyle 1/2\,,$
$\displaystyle\Gamma_{7}(\bm{p},\bm{q})$ $\displaystyle\equiv$
$\displaystyle(q^{2}\gamma_{\delta\theta\delta}(\bm{p},\bm{q})+p^{2}\gamma_{\delta\theta\delta}(\bm{q},\bm{p}))/(\bm{p}+\bm{q})^{2}\,,$
(64)

with $\gamma_{A\theta A}(\bm{p},\bm{q})=\bm{p}\cdot\bm{q}/(2q^{2})$. Each of
the $\Gamma_{j}$ approaches a constant when the sum of wavenumbers goes to
zero, leading to a scaling of $F_{2,\delta}$ in accordance with the
requirement from mass and momentum conservation (see Sec. III.2).

Table 1: Expansion coefficients in Eq. (65) of the functions $\Gamma_{j}(\bm{p},\bm{q})$ in terms of the four basis functions given by Eq. (II.3.2). $j$ | $1$ | $2$ | $3$ | $4$ | $5$ | $6$ | $7$
---|---|---|---|---|---|---|---
$\beta_{j,a}$ | $0$ | $-1$ | $\frac{1}{2}$ | $1$ | $-\frac{1}{2}$ | $-\frac{1}{2}$ | $\frac{3}{2}$
$\beta_{j,b}$ | $1$ | $2$ | $\frac{1}{2}$ | $0$ | $\frac{1}{2}$ | $\frac{1}{2}$ | $-\frac{1}{2}$
$\beta_{j,c}$ | $0$ | $-\frac{3}{14}$ | $\frac{5}{14}$ | $\frac{5}{7}$ | $\frac{1}{7}$ | $-\frac{5}{14}$ | $\frac{1}{14}$
$\beta_{j,d}$ | $\frac{2}{7}$ | $\frac{13}{14}$ | $\frac{1}{7}$ | $0$ | $-\frac{1}{7}$ | $-\frac{1}{7}$ | $\frac{1}{7}$

Using the explicit expressions for the vertices, one can check that each

$\Gamma_{j}=\beta_{j,a}\Delta_{a}+\beta_{j,b}\Delta_{b}+\beta_{j,c}\Delta_{c}+\beta_{j,d}\Delta_{d}\,,$
(65)

can be expressed in terms the four basis functions

$\displaystyle\Delta_{a}(\bm{p},\bm{q})$ $\displaystyle\equiv$
$\displaystyle(p^{2}+q^{2})F_{2}(\bm{p},\bm{q})/(\bm{p}+\bm{q})^{2}\,,$
$\displaystyle\Delta_{b}(\bm{p},\bm{q})$ $\displaystyle\equiv$ $\displaystyle
F_{2}(\bm{p},\bm{q})\,,$ $\displaystyle\Delta_{c}(\bm{p},\bm{q})$
$\displaystyle\equiv$
$\displaystyle(p^{2}+q^{2})K(\bm{p},\bm{q})/(\bm{p}+\bm{q})^{2}\,,$
$\displaystyle\Delta_{d}(\bm{p},\bm{q})$ $\displaystyle\equiv$ $\displaystyle
K(\bm{p},\bm{q})\equiv(\bm{p}\cdot\bm{q})^{2}/(p^{2}q^{2})-1\,,\ $ (66)

with coefficients $\beta_{j,k}$ given in Table 1. This basis has been
introduced in Eggemeier  _et al._ (2019), and the present result proves
equation (109) therein. Here $K$ is the kernel corresponding to the second
order Galileon operator ${\cal G}_{2}$ in Eggemeier  _et al._ (2019), which
physically represents the tidal field.

The decomposition in Eq. (65) implies that one can write the first-order
correction in $\epsilon$ to the second-order density kernel $F_{2,\delta}$
also as

$\displaystyle F_{2,\delta}(\bm{p},\bm{q},\eta)$ $\displaystyle=$
$\displaystyle
F_{2}(\bm{p},\bm{q})-(\bm{p}+\bm{q})^{2}\sum_{k=a,b,c,d}\Delta_{k}(\bm{p},\bm{q})\beta^{\delta}_{k}(\eta)$
$\displaystyle{}+{\cal O}(\epsilon^{2})\,,$ $\displaystyle
F_{2,\theta}(\bm{p},\bm{q},\eta)$ $\displaystyle=$ $\displaystyle
G_{2}(\bm{p},\bm{q})-(\bm{p}+\bm{q})^{2}\sum_{k=a,b,c,d}\Delta_{k}(\bm{p},\bm{q})\beta^{\theta}_{k}(\eta)$
(67) $\displaystyle{}+{\cal O}(\epsilon^{2})\,,$

with coefficients $\beta^{\delta/\theta}_{k}(\eta)$ for $k=a,b,c,d$ given by
linear combinations of the $J^{\delta/\theta}_{j}(\eta)$ with coefficients
$\beta_{j,k}$ from Table 1. Using the results from App. B we find

$\displaystyle\beta^{\delta}_{a}(\eta)$ $\displaystyle=$
$\displaystyle\frac{2}{5}(2E_{0}(\eta)-6E_{1}(\eta)+5E_{2}(\eta)+3E_{5/2}(\eta)$
$\displaystyle{}-5E_{3}(\eta)+E_{7/2}(\eta))\,,$
$\displaystyle\beta^{\delta}_{b}(\eta)$ $\displaystyle=$
$\displaystyle\frac{2}{5}(6E_{1}(\eta)-10E_{2}(\eta)+5E_{3}(\eta)-E_{7/2}(\eta))\,,$
$\displaystyle\beta^{\delta}_{c}(\eta)$ $\displaystyle=$
$\displaystyle\frac{1}{35}(3E_{1}(\eta)-20E_{2}(\eta)+30E_{5/2}(\eta)-15E_{3}(\eta)$
$\displaystyle{}+2E_{7/2}(\eta))\,,$ $\displaystyle\beta^{\delta}_{d}(\eta)$
$\displaystyle=$
$\displaystyle\frac{3}{35}(3E_{1}(\eta)-10E_{2}(\eta)+15E_{3}(\eta)-8E_{7/2}(\eta))\,,$

and

$\displaystyle\beta^{\theta}_{a}(\eta)$ $\displaystyle=$
$\displaystyle\frac{1}{10}(8E_{0}(\eta)-24E_{1}(\eta)+40E_{2}(\eta)-33E_{5/2}(\eta)$
$\displaystyle{}+20E_{3}(\eta)-6E_{7/2}(\eta))\,,$
$\displaystyle\beta^{\theta}_{b}(\eta)$ $\displaystyle=$
$\displaystyle\frac{1}{10}(24E_{1}(\eta)-20E_{2}(\eta)+15E_{5/2}(\eta)-20E_{3}(\eta)$
$\displaystyle{}+6E_{7/2}(\eta))\,,$ $\displaystyle\beta^{\theta}_{c}(\eta)$
$\displaystyle=$
$\displaystyle\frac{1}{70}(16E_{0}(\eta)+6E_{1}(\eta)+60E_{2}(\eta)-81E_{5/2}(\eta)$
$\displaystyle{}+30E_{3}(\eta)-6E_{7/2}(\eta))\,,$
$\displaystyle\beta^{\theta}_{d}(\eta)$ $\displaystyle=$
$\displaystyle\frac{1}{35}(9E_{1}(\eta)+20E_{2}(\eta)-15E_{5/2}(\eta)-45E_{3}(\eta)$
(69) $\displaystyle{}+36E_{7/2}(\eta))\,.$

The kernels for the third cumulant perturbations $\pi$ and $\chi$ start at
order $\epsilon^{2}$, and can be obtained by inserting the results up to first
order from above in the right-hand side of their equations of motion

$\displaystyle\left(\partial_{\eta}+\frac{7}{2}\right)F_{2,\pi}$
$\displaystyle=$ $\displaystyle 3k^{2}\epsilon F_{2,g}+5k^{2}\epsilon
F_{2,\delta\epsilon}+k^{2}\omega F_{2,A}$ $\displaystyle{}+\gamma_{\pi
bc}F_{1,b}F_{1,c}\,,$
$\displaystyle\left(\partial_{\eta}+\frac{7}{2}\right)F_{2,\chi}$
$\displaystyle=$ $\displaystyle 5k^{2}\epsilon
F_{2,\delta\epsilon}+k^{2}\omega F_{2,A}+\gamma_{\chi bc}F_{1,b}F_{1,c}\,,$

and using the respective vertices.

In addition, at second order, also vorticity, vector and tensor modes are
generated. The corresponding kernels are discussed in Sec. IV and Sec. IV.4,
respectively.

In principle, with the same strategy one can obtain the correction terms to
the third-order kernels $F_{3,a}$ at linear order in $\epsilon$. Similarly, it
is possible to systematically expand to higher powers in $\epsilon$.
Nevertheless, the expressions become too lengthy to be shown, and we revert to
a numerical treatment, that does _not_ require any expansion in $\epsilon$ and
therefore remains valid when $\epsilon k_{i}^{2}$ is sizeable.

####  II.3.3 Wilson coefficients

Note that the four basis functions in Eq. (II.3.2) appearing in the order
$\epsilon$ correction to the density kernel
$F_{2,\delta}(\bm{p},\bm{q},\eta)$, Eq. (II.3.2), are equivalent to the shape
functions introduced in the context of the effective field theory (EFT)
treatment of the bispectrum Baldauf  _et al._ (2015); Angulo  _et al._
(2015)666Using the basis functions $E_{1,2,3}$ and $\Gamma$ in the convention
of Baldauf  _et al._ (2021) one has
$E_{1}=k^{2}(-\Delta_{a}+\Delta_{b}-\frac{5}{7}\Delta_{c}-\frac{2}{7}\Delta_{d})$,
$E_{2}=k^{2}(-\frac{2}{3}\Delta_{a}+\frac{2}{3}\Delta_{b}-\frac{10}{21}\Delta_{c}+\frac{17}{21}\Delta_{d})$,
$E_{3}=k^{2}(-\frac{1}{3}\Delta_{a}+\frac{1}{3}\Delta_{b}+\frac{11}{42}\Delta_{c}-\frac{2}{21}\Delta_{d})$,
$\Gamma=k^{2}(\frac{2}{11}\Delta_{a}+\frac{9}{11}\Delta_{b}-\frac{1}{7}\Delta_{c}-\frac{6}{77}\Delta_{d})$,
where $k^{2}\equiv(\bm{p}+\bm{q})^{2}$.. Indeed, within the effective field
theory language the coefficients
$\beta^{\delta}_{a},\beta^{\delta}_{b},\beta^{\delta}_{c},\beta^{\delta}_{d}$
are the set of four counter-terms (or rather Wilson coefficients) at second
order in perturbation theory and at leading order in the derivative expansion.
The result obtained here can be viewed as a first-principle _determination_ of
those Wilson coefficients by matching the EFT to the exact UV theory, being
the Vlasov equation (at least when considering cold dark matter only).
Nevertheless, this does not mean that their values are identical to the
corresponding counter-terms obtained when complementing an SPT loop
computation with correction terms of the form of Eq. (II.3.2) and fitting them
to simulation results, since these counter-terms have to absorb also the
unphysical contributions from the UV region of the SPT loop integration (or
alternatively account for the missing modes when imposing a cutoff). Instead,
the Wilson coefficients
$\beta^{\delta}_{a},\beta^{\delta}_{b},\beta^{\delta}_{c},\beta^{\delta}_{d}$
obtained here capture the actual impact of non-linearly induced dark matter
velocity dispersion on the second-order kernel. Similarly, $I^{\delta}$
appearing in the $\epsilon$ correction to $F_{1,\delta}$, see Eq. (II.3.1),
may be viewed as the result of matching the EFT correction at first order in
perturbations and leading order in gradients to the Vlasov theory.

Apart from that, we stress that we do not expand in powers of $\epsilon$ in
our numerical analysis (see below), and therefore our results go beyond the
four correction terms in Eq. (II.3.2). For example, formally, the
contributions to $F_{2,\delta}$ at order $\epsilon^{2}$ would correspond to
additional higher-derivative operators in the EFT with two additional spatial
gradients compared to those captured by Eq. (II.3.2). At order $\epsilon^{3}$,
the EFT matching would require to include operators with four additional
gradients, and so on. Furthermore, the correction to the third-order kernel
$F_{3,\delta}$ would correspond to EFT operators at third order in
perturbation theory, etc. All of these contributions are implicitly contained
in the kernels $F_{n,a}$ considered here.

###  II.4 Numerical results

Table 2: Background values and perturbation modes taken into account in various approximation schemes of . | s | sw | sv | svt
---|---|---|---|---
cum2 | $\epsilon(\eta)$ | $\delta,\theta,g,\delta\epsilon,A$ | $+w_{i}$ | $+\nu_{i}$ | $+t_{ij}$
cum3+ | $\epsilon(\eta),\omega(\eta)$ | $\delta,\theta,g,\delta\epsilon,A,\pi,\chi$ | $+w_{i}$ | $+\nu_{i}$ | $+t_{ij}$

Figure 2:  Non-linear kernels $F_{2,a}(\bm{k}-\bm{q},\bm{q},\eta)$ (left
column) and $F_{3,a}(\bm{k},\bm{q},-\bm{q},\eta)$ (right column) versus
$q=|\bm{q}|$ in (blue lines) compared to SPT (black dashed) for the density
contrast $a=\delta$ (top row) and velocity divergence $a=\theta$ (bottom row).
For the background dispersion we chose
$\epsilon(\eta)=\epsilon_{0}\,e^{\alpha\eta}$ as well as
$\omega(\eta)=\bar{\omega}\times\epsilon(\eta)^{2}$, with parameters as
indicated in the figures, and $\eta=0$. We show the dependence on $q$ relative
to the dispersion scale $k_{\sigma}$, such that the result applies to any
value of $\epsilon_{0}=1/k_{\sigma}^{2}$. The blue lines show second and third
cumulant approximations (see Table 2), which are both suppressed compared to
the EdS-SPT kernels for $q\gtrsim k_{\sigma}$ (shown in black dashed lines).
The thin black solid line shows the analytical result of Eq. (II.3.2) at first
order in $\epsilon$.  Figure 3:  Non-linear kernels
$F_{3,a}(\bm{k},\bm{q},-\bm{q},\eta)$ versus $q=|\bm{q}|$ for
$a=\delta,\theta$ as in Fig. 2, but comparing the impact of various
approximations related to vorticity, as well as vector and tensor modes of the
velocity dispersion tensor (see Table 2). Including vorticity backreaction has
a significant impact on $F_{3,\delta}$ within the most relevant regime
$k_{\sigma}\lesssim q\lesssim 2k_{\sigma}$ where the suppression of relative
to SPT sets in ((sw) vs (s) in upper panel).

In our analysis we solve the system of coupled differential equations, Eq.
(II.2), for the non-linear kernels numerically using a recursive algorithm
Blas  _et al._ (2015); Garny and Taule (2021).

We consider the approximation schemes given in Table 2. The most inclusive
scheme is (cum3+/svt), for which apart from $\delta$ and $\theta$ also all
scalar, vector and tensor perturbations of the dispersion tensor
$\epsilon_{ij}$, its background expectation value $\epsilon(\eta)$, as well as
the scalar modes $\pi$ and $\chi$ of the third cumulant are taken into
account. Furthermore, their evolution equations contain the expectation value
$\omega(\eta)$ of the fourth cumulant, that we also include. For comparison,
we consider approximations (cum2) that omit the third cumulant perturbations.
In that case, $\omega(\eta)$ does not enter in the perturbation equations. In
addition, we consider various approximations related to scalar, vector and
tensor modes. For approximation scheme (s), only scalar modes are included.
For (sw) we include in addition vorticity. Within the schemes (sv) and (svt)
also the vector and tensor modes of $\epsilon_{ij}$ are added successively.

Since the second and higher cumulants are expected to play a role only at low
redshift, we initialize the kernels at some finite time $\eta_{\text{ini}}$,
using EdS-SPT kernels for the density and velocity divergence, and setting all
other kernels to zero (except $A$, which is initialized as $\delta$),

$\displaystyle F_{n,\delta}(\bm{k}_{1},\dots,\bm{k}_{n},\eta_{\text{ini}})$
$\displaystyle\equiv$ $\displaystyle F_{n}(\bm{k}_{1},\dots,\bm{k}_{n})\,,$
$\displaystyle F_{n,\theta}(\bm{k}_{1},\dots,\bm{k}_{n},\eta_{\text{ini}})$
$\displaystyle\equiv$ $\displaystyle G_{n}(\bm{k}_{1},\dots,\bm{k}_{n})\,,$
$\displaystyle F_{n,A}(\bm{k}_{1},\dots,\bm{k}_{n},\eta_{\text{ini}})$
$\displaystyle\equiv$ $\displaystyle F_{n}(\bm{k}_{1},\dots,\bm{k}_{n})\,,$
$\displaystyle F_{n,a}(\bm{k}_{1},\dots,\bm{k}_{n},\eta_{\text{ini}})$
$\displaystyle\equiv$ $\displaystyle 0,\qquad
a=w_{i},g,\delta\epsilon,\nu_{i},t_{ij},\pi,\chi\,.$

The most important part of the initial condition are the first two lines, and
in particular for the linear kernels $n=1$. These conditions must be
supplemented by an initial condition on the background dispersion $\epsilon$
of infinitesimal value, which makes the Vlasov cumulant hierarchy evolution
consistent with the solution of the Vlasov equation777See Sec. III.B in
Pueblas and Scoccimarro (2009) for more discussion on this., otherwise no
dispersion and higher cumulants would be generated by time evolution out of
Eqs. (II.4). As long a $\eta_{\text{ini}}$ is chosen early enough, the results
at low redshift are insensitive to the choice of the initial condition for the
higher kernels with $n>1$, as well as for the perturbations of $A$ and the
second and higher cumulants, since the dominant contribution is generated by
the time-evolution at late times. In practice we use $\eta_{\text{ini}}=-20$,
and checked that our results are insensitive to the precise choice of the
initial time, as expected since transients will be down by at least
$\exp(-20)$. Our implementation for the vector and tensor modes is described
in App. D.

Figure 4:  Non-linear kernels $F_{3,a}(\bm{k},\bm{q},-\bm{q},\eta)$ versus
$k=|\bm{k}|$ for $a=\delta,\theta$ (top) and the scalar perturbations
$a=g,\delta\epsilon$ of the stress tensor (bottom). $F_{3,\delta}\propto
k^{2}$ for $k\to 0$ is ensured by mass and momentum conservation, while for
$F_{3,\theta}$ it is a consequence of the symmetry $\bm{q}\to-\bm{q}$, see
Sec. III.2. In contrast, the kernels of the velocity dispersion modes approach
a constant at low $k$.

The numerical result for the kernels $F_{2,a}(\bm{k}-\bm{q},\bm{q},\eta)$ and
$F_{3,a}(\bm{k},\bm{q},-\bm{q},\eta)$, that enter in the one-loop density
($a=\delta$) and velocity divergence ($a=\theta$) power spectra, Eq. (II.2),
are shown in Fig. 2 (see also Fig. 3). We show the dependence on the magnitude
of the loop wavenumber $q=|\bm{q}|$ for fixed $k=|\bm{k}|$ and
$\cos_{kq}=\mu=\bm{k}\cdot\bm{q}/(kq)$, as well as assuming a power-law time-
dependence $\epsilon(\eta)=\epsilon_{0}\,e^{\alpha\eta}$ and
$\omega(\eta)=\bar{\omega}\times\epsilon(\eta)^{2}$, with parameters as stated
in the figure. Note that the kernels are dimensionless, such that the result
is independent of $\epsilon_{0}=1/k_{\sigma}^{2}$ when plotted against the
ratio $q/k_{\sigma}$. The numerical results for the density kernels (in blue)
approach the corresponding EdS-SPT kernels (black dashed) in the limit when
both $q\ll k_{\sigma}$ and $k\ll k_{\sigma}$. Due to our choice
$k/k_{\sigma}=0.2$, the blue and black dashed curves therefore lie (almost) on
top of each other for $q/k_{\sigma}\ll 1$, as expected. Furthermore, we show
the analytical result at order $\epsilon$ (Eq. II.3.2) for $F_{2,\delta}$ and
$F_{2,\theta}$. It agrees with the numerical result for $q/k_{\sigma}\lesssim
1$, and captures the onset of the deviation of from SPT.

In the limit of large loop wavenumber $q$ we observe a substantial
_suppression_ of the kernels relative to EdS-SPT. This implies that, when
integrating over the loop wavenumber $q$, the one-loop correction to the
density power spectrum is less sensitive to the UV regime as compared to SPT.
This result is in line with theoretical expectations Peebles (1980) and in
qualitative agreement with numerical findings on the response of the non-
linear power spectrum to small-scale initial fluctuations Nishimichi  _et al._
(2016, 2017), which reflects the decoupling of small-scale modes. The
inclusion of velocity dispersion and higher cumulants in  therefore addresses
one of the major shortcomings of SPT, namely the sensitivity to UV modes, by
taking the physically expected screening of those perturbation modes into
account.

The asymptotic slope of the kernels
$F_{n,\delta/\theta}(\bm{k}_{1},\dots,\bm{k}_{n},\eta)$ in the limit
$k\equiv|\sum\bm{k}_{i}|\ll k_{\sigma}\ll|\bm{k}_{i}|$ can be understood
analytically. Within SPT momentum conservation implies $F_{n}/G_{n}\propto
k^{2}/q^{2}$ in the corresponding limit $k\ll|\bm{k}_{i}|$, where
$q\equiv\text{max}_{i}|\bm{k}_{i}|$, and the same factor arises in . In
addition, within the dispersion scale $k_{\sigma}$ is relevant. For all
wavenumbers with $|\bm{k}_{i}|\gg k_{\sigma}$, the corresponding modes
contributing to the kernel have already entered the “dispersion horizon”. The
moment of entry corresponds to time $\eta=\eta_{k_{i}}$ satisfying the
condition $\epsilon(\eta_{k_{i}})k_{i}^{2}=1$. For
$\epsilon(\eta)=\epsilon_{0}e^{\alpha\eta}$ this corresponds to
$e^{\eta_{k_{i}}}=(k_{\sigma}^{2}/k_{i}^{2})^{1/\alpha}$. Now, the most
important effect is that the linear modes effectively stop growing in the
usual way (i.e. $\propto e^{\eta}$) once entering the dispersion horizon. This
implies that

$\displaystyle
F_{n,\delta/\theta}(\bm{k}_{1},\dots,\bm{k}_{n},\eta)\sim\frac{k^{2}}{q^{2}}\,e^{\eta_{k_{1}}}\,e^{\eta_{k_{2}}}\cdots
e^{\eta_{k_{n}}}$ $\displaystyle\text{for}\quad k\ll k_{\sigma}\ll q\,,$ (72)

within . Using $e^{\eta_{k_{i}}}=(k_{\sigma}^{2}/k_{i}^{2})^{1/\alpha}$ this
yields a power-law screening of UV modes. Indeed, applying this argument to
$F_{2,a}(\bm{k}-\bm{q},\bm{q},\eta)\propto q^{-2-4/\alpha}$ and
$F_{3,a}(\bm{k},\bm{q},-\bm{q},\eta)\propto q^{-2-4/\alpha}$ yields the slope
by which these kernels decay at very large $q$. Using that for a power-law
universe $\alpha=4/(n_{s}+3)$ this furthermore explains the universal
$1/q^{2}$ decay stated in Eq. (3) and shown in Fig. 1. Moreover, the argument
leading to this scaling is independent of the precise way on how the linear
growth stalls when entering the dispersion horizon, and therefore shows that
the asymptotic scaling of the kernels is universal, independent of the precise
approximation scheme. We shall also confirm this with our numerical solutions
in the following.

In Fig. 2 we compare the (cum2) and (cum3+) approximation schemes (see Table
2). For $F_{2,\delta(\theta)}$ the difference is very small (with both lines
being almost on top of each other), while for $F_{3,\delta(\theta)}$ a more
noticeable difference arises. But overall, the impact of third cumulant
perturbations is small within the regime $q\sim{\cal O}(k_{\sigma})$ where the
suppression relative to SPT sets in.

In Fig. 3, we show the dependence of the kernel
$F_{3,\delta(\theta)}(\bm{k},\bm{q},-\bm{q},\eta)$ on the approximation scheme
for vorticity, vector and tensor modes (see Table 2; note that only scalar
modes contribute to $F_{2,\delta(\theta)}$). Within the regime $q\lesssim
2k_{\sigma}$, where the suppression relative to SPT becomes important, the
(sw), (sv) and (svt) schemes are all close to each other for $F_{3,\delta}$,
while the (s) approximation differs significantly. This implies that the
backreaction of vorticity on the third-order density kernel is sizeable, and
leads to a substantial suppression of the $F_{3,\delta}$ kernel. The
difference between (s) and (sw) captures the impact of vorticity, while
including the vector mode of the stress tensor in addition within the (sv)
scheme leads to further differences at large $q\gtrsim 2k_{\sigma}$. On the
other hand, we find that the impact of the tensor perturbation of the stress
tensor on both $F_{3,\delta}$ and $F_{3,\theta}$ is very minor, as quantified
by the difference between (sv) and (svt). We find a qualitatively similar
behaviour, including in particular the suppression at large $q$ relative to
SPT, independently of the precise choice of the parameters $k$, $\cos_{kq}$,
$\alpha$ and $\bar{\omega}$.

The kernels $F_{3,a}(\bm{k},\bm{q},-\bm{q},\eta)$ of the two scalar
perturbation modes $a=g,\delta\epsilon$ of the velocity dispersion tensor are
shown in Fig. 4, as compared to those for $a=\delta,\theta$. Note that we show
the dependence on $k$ for fixed $q$ here. For $a=\delta,\theta$ the kernels
scale as $k^{2}$ for small $k$, as in SPT, and guaranteed by mass and momentum
conservation for the density, as well as by symmetry $\bm{k}\to-\bm{k}$ for
$\theta$ (see Sec. III). For large $k\gtrsim k_{\sigma}$ the kernels are
suppressed compared to their SPT values, similarly as for the dependence on
$q$.

The behaviour for the stress tensor perturbations $a=g,\delta\epsilon$ is
qualitatively different for small $k$, where their kernels approach a constant
value. This implies that their linear theory values are modified also in the
limit $k\to 0$, as quantified e.g. by a “bias” $b_{a}\equiv
P_{a\delta}/P_{\delta\delta}|_{k\to 0}$. In the limit $k\to 0$, the bias is
given by

$\displaystyle b_{a}$ $\displaystyle=$ $\displaystyle F_{1,a}(k=0,\eta)+3\int
d^{3}q\,F_{3,a}(\bm{k},\bm{q},-\bm{q},\eta)\,P_{0}(q)|_{k\to 0}$ (73)
$\displaystyle{}+\dots\,,$

where the ellipsis stand for two and higher loop contributions. For $a=g$ the
linear contribution is given by
$b_{g}^{\text{lin}}=F_{1,g}(k=0,\eta)=2E_{2}(\eta)$ (see Eq. 45), while for
the case $a=\delta\epsilon$ we have $b_{\delta\epsilon}^{\text{lin}}=0$888Note
that the ${\cal O}(\epsilon^{2})$ contributions to $F_{1,a}(k,\eta)$ have to
vanish for $a=g,\delta\epsilon$ in the limit $k\to 0$ for dimensional reasons,
being proportional to some time integral of either
$\epsilon(\eta^{\prime})\epsilon(\eta^{\prime\prime})k^{2}$ or
$\omega(\eta^{\prime})k^{2}$..

Overall, we conclude that vorticity as well as scalar and vector modes of the
velocity dispersion tensor should be taken into account when computing the
one-loop correction to the density power spectrum, while third cumulant
perturbations and especially the tensor mode are less relevant. Taking
vorticity into account is also important to guarantee momentum conservation,
as we shall see in the next section.

##  III Symmetry constraints on non-linear kernels

The symmetries underlying the Vlasov-Poisson system lead to constraints on the
behaviour of the non-linear kernels
$F_{n,a}(\bm{k}_{1},\dots,\bm{k}_{n},\eta)$ in particular limiting cases. In
the following we discuss the implications of

  * •

Galilean invariance in the limit where one of the arguments $\bm{k}_{i}\to 0$,
corresponding to the impact of a very large scale mode on the propagation on
smaller scales, and

  * •

mass and momentum conservation in the limit when the total wavenumber
$\bm{k}\equiv\sum_{i}\bm{k}_{i}\to 0$, while the individual $\bm{k}_{i}$
remain finite, corresponding to the impact of small-scale modes onto
perturbation modes on larger scales.

The corresponding constraints are well-known within EdS-SPT Zeldovich (1965);
Peebles (1980); Goroff  _et al._ (1986), and here we discuss and show how they
are generalized when including velocity dispersion and higher cumulants. Since
these symmetry constraints guarantee delicate cancellations among various
terms at a given order in perturbation theory, it is important to ensure their
validity for a given approximation scheme. For example, as we will see, and
unlike what happens in SPT, momentum conservation requires to take vorticity
into account. Finally, the symmetry constraints also provide a rather non-
trivial check of the formalism and implementation.

###  III.1 Galilean invariance

The Vlasov-Poisson system within the non-relativistic limit features a shift
symmetry

$\bm{x}\to\bm{x}^{\prime}=\bm{x}+\bm{n}(\tau),\quad\tau\to\tau^{\prime}=\tau\,,$
(74)

with an arbitrary time-dependent shift function $\bm{n}(\tau)$, as pointed out
in Kehagias and Riotto (2013); Peloso and Pietroni (2013) generalizing the
earlier work in Scoccimarro and Frieman (1996b). More broadly, this symmetry
is a remnant of those induced by general relativistic coordinate invariance,
being the only one that has a non-trivial Newtonian limit Creminelli  _et al._
(2013); Horn  _et al._ (2014), corresponding to the equivalence principle. It
can be viewed as a generalized Galilean invariance (in the absence of gravity
and expansion of the universe, ${\bm{n}}$ becomes linear in time), with
gravitational potential transforming as
$\Phi(\tau,\bm{x})\to\Phi(\tau^{\prime},\bm{x}^{\prime})=\Phi^{\prime}(\tau,\bm{x})-(d^{2}{\bm{n}}/d\tau^{2}+{\cal
H}d{\bm{n}}/d\tau)\cdot\bm{x}$. While the density contrast transforms like a
scalar under this symmetry,
$\delta(\tau,\bm{x})\to\delta(\tau^{\prime},\bm{x}^{\prime})=\delta^{\prime}(\tau,\bm{x})$,
the velocity field is shifted by a time-dependent, but spatially constant
vector,
$\bm{v}(\tau,\bm{x})\to\bm{v}(\tau^{\prime},\bm{x}^{\prime})=\bm{v}^{\prime}(\tau,\bm{x})+d\bm{n}/d\tau$.

In order to determine the transformation properties of higher cumulants, we
use that the phase-space distribution function itself is a scalar under
Galilean transformations, $f(\tau,\bm{x},\bm{p})\to
f(\tau^{\prime},\bm{x}^{\prime},\bm{p}^{\prime})=f^{\prime}(\tau,\bm{x},\bm{p})$,
where $\bm{p}^{\prime}=\bm{p}+a(\tau)d\bm{n}/d\tau$. From this, one finds that
the cumulant generating function, Eq. (4), transforms as

${\cal C}(\tau,\bm{x},\bm{l})\to{\cal
C}(\tau^{\prime},\bm{x}^{\prime},\bm{l})={\cal
C}^{\prime}(\tau,\bm{x},\bm{l})+\bm{l}\cdot d\bm{n}/d\tau\,,$ (75)

which implies that all but the first cumulant, i.e. the velocity, are scalars
under Galilean transformations. In particular,

$\displaystyle\sigma_{ij}(\tau,\bm{x})$ $\displaystyle\to$
$\displaystyle\sigma_{ij}(\tau^{\prime},\bm{x}^{\prime})=\sigma_{ij}^{\prime}(\tau,\bm{x})\,,$
$\displaystyle{\cal C}_{ijk\cdots}(\tau,\bm{x})$ $\displaystyle\to$
$\displaystyle{\cal C}_{ijk\cdots}(\tau^{\prime},\bm{x}^{\prime})={\cal
C}_{ijk\cdots}^{\prime}(\tau,\bm{x})\,,$ (76)

for the velocity dispersion tensor, the third and all higher cumulants. Note
that also the divergence of the velocity $\theta$ as well as the vorticity
$w_{i}$ transform as scalars, or in general any gradient $\nabla_{i}v_{j}$ of
the velocity field, as well as second gradients $\nabla_{i}\nabla_{j}\Phi$ of
the gravitational potential.

Galilean invariance of the equations of motion of the density contrast, the
velocity divergence and vorticity, as well as the dispersion tensor and higher
cumulants therefore implies that only covariant quantities may enter these
equations. This requires in particular that only gradients of the velocity
field, and second gradients of the gravitational potential may appear, as can
readily be checked to be the case for the continuity and Euler equations (6),
and also for the equations of motion (II.1) of the dispersion tensor and the
higher cumulants (note that the equations for $\theta$ and $w_{i}$ contain an
extra gradient compared to the Euler equation for $\bm{v}$). The only
exception is the time-derivative. Galilean invariance requires that it appears
in the covariant combination $\partial_{\tau}+{\bf v}\cdot\nabla$. This is
indeed the case in Eqs. (6) and (II.1). In terms of the rescaled time variable
$\eta$ and velocity $\bm{u}$, the covariant derivative reads

$D_{\eta}\equiv\partial_{\eta}-\bm{u}\cdot\nabla=\partial_{\eta}-\left(\frac{\nabla\theta}{\nabla^{2}}-\frac{\nabla\times\bm{w}}{\nabla^{2}}\right)\cdot\nabla\,.$
(77)

Galilean invariance therefore implies that the covariant time derivative is
the only place where the velocity field may appear _without_ a gradient acting
on it, when considering equations of motion for quantities that are scalars
under Galilean transformations. As shown above, this is the case for all
perturbation variables considered here. When transforming the equations of
motion into Fourier space as in Eq. (23), this implies that the non-linear
vertices satisfy the property

$\gamma_{abc}(\bm{p},\bm{q})=\delta_{ac}\left(\delta_{b\theta}\frac{\bm{q}\cdot\bm{p}}{2p^{2}}-\delta_{bw_{i}}\frac{(\bm{q}\times\bm{p})_{i}}{2p^{2}}\right)+{\cal
O}(p^{0})\,,$ (78)

in the limit $\bm{p}\to 0$, generalizing the result in Crocce and Scoccimarro
(2006) to include vorticity. As shown in App. C, this implies the relation

$\displaystyle F_{n+1,a}(\bm{p},\bm{k}_{1},\dots,\bm{k}_{n},\eta)$
$\displaystyle=$
$\displaystyle\frac{1}{n+1}\frac{\bm{k}\cdot\bm{p}}{p^{2}}F_{n,a}(\bm{k}_{1},\dots,\bm{k}_{n},\eta)$
(79) $\displaystyle{}+{\cal O}(p^{0})\ \ \ \text{for}\ \bm{p}\to 0\,,$

for the squeezed limit of the kernels, where $\bm{k}=\sum_{i}\bm{k}_{i}$. This
property is well-known for the EdS-SPT kernels, see e.g. Sugiyama and Futamase
(2013). Galilean invariance ensures its validity also for the non-linear
kernels, including those of higher cumulants. The limit is approached when
$|\bm{p}|\ll|\bm{k}_{i}|$ _and_ $|\bm{p}|\ll k_{\sigma}$, i.e. the size $\sim
1/\bm{p}$ of the large-scale mode has to be larger than the dispersion length-
scale $\sqrt{\epsilon}$. We checked that the relation above is satisfied for
the analytical results of the second-order kernels at order $\epsilon$. In
addition we checked its validity analytically for the third-order kernels
$F_{3,a}$ for $a=\delta,\theta,g,\delta\epsilon,A$ expanded to order
$\epsilon$.

Equation (79) ensures that, as in SPT Scoccimarro and Frieman (1996b); Jain
and Bertschinger (1996), all contributions to the equal-time power spectra and
bispectra for which the kernels become singular in the limit of small loop
wavenumber cancel. This cancellation occurs only after adding all
contributions at a given order in perturbation theory, i.e. for example among
the various lines in the one-loop expression, Eq. (II.2). It is easy to see
that Eq. (79) also implies that the cancellation can be made manifest on the
level of the loop integrand by the same strategy as in SPT, for the power
spectrum Blas  _et al._ (2013, 2013) and bispectrum Floerchinger  _et al._
(2019), respectively. We already described its implementation at one-loop
above, and use it throughout this work.

For unequal times, one can check that Eq. (79) gives rise to the standard
$(1/p)$ pole present in the so-called consistency relations between the
squeezed $(N+1)$-point functions of density fluctuations and its regular
$N$-point function Kehagias and Riotto (2013); Peloso and Pietroni (2013).

###  III.2 Mass and momentum conservation

Within SPT, it is well-known that $F_{n}(\bm{k}_{1},\dots,\bm{k}_{n})\propto
k^{2}$ and $G_{n}(\bm{k}_{1},\dots,\bm{k}_{n})\propto k^{2}$ in the large-
scale limit when the sum $\bm{k}\equiv\sum_{i}\bm{k}_{i}$ of wavevectors goes
to zero (with the magnitude of at least two of the individual wavevectors
remaining constant). Here we show that, as expected from mass and momentum
conservation, this property is preserved for the kernels $F_{n,\delta}$ of the
density contrast also in provided vorticity is taken into account, as hinted
at already above. However, as we show below, the scaling of the velocity
divergence kernels $F_{n,\theta}$ in is in general different compared to SPT,
but still conforms with all symmetry requirements.

####  III.2.1 Density contrast

Figure 5:  Non-linear kernels
$F_{3,a}(\bm{k}-\bm{p}-\bm{q},\bm{p},\bm{q},\eta)$ for $a=\delta$ (upper
panel) and $a=\theta$ (lower panel) versus $k=|\bm{k}|$, for fixed $\bm{p}$
and $\bm{q}$ at $\eta=0$, and assuming $\epsilon=\epsilon_{0}\,e^{\alpha\eta}$
with $\alpha=4/3$. The blue lines show the (s), (sw) and (sv) approximation
schemes, respectively, for (cum3+) with $\bar{\omega}=1$. The black dashed
lines show the EdS-SPT kernels $F_{3}$ and $G_{3}$ for comparison. For
$F_{3,\delta}$, a spurious linear scaling with $k$ occurs in the limit $k\to
0$ when vorticity is neglected (s), while the $k^{2}$ scaling ensured by mass
and momentum conservation is realized for all other approximation schemes that
include vorticity. For $F_{3,\theta}$, the linear scaling with $k$ is a
physical effect, and occurs independently of the approximation scheme as long
as $\bm{p}+\bm{q}\not=0$. For the figure we chose cosines
$c_{kp}=0.3,c_{kq}=0.5,c_{pq}=0.875,p=0.5k_{\sigma},q=0.5k_{\sigma}$ where
$k_{\sigma}=1/\sqrt{\epsilon_{0}}$.

Let us start with the kernels $F_{n,\delta}$ of the density contrast. To study
their asymptotic behaviour it is convenient to write the continuity equation
in the form

$\partial_{\tau}\delta+\nabla_{i}P_{i}=0\,,$ (80)

where

$P_{i}(\tau,\bm{x})\equiv(1+\delta)v_{i}\,,$ (81)

is the momentum field. Using the Euler equation, its equation of motion can be
written as

$\partial_{\tau}P_{i}+{\cal H}P_{i}+\nabla_{j}T_{ij}=0\,,$ (82)

which takes the form of a conservation equation for the comoving momentum
$aP_{i}$, with

$\displaystyle T_{ij}$ $\displaystyle\equiv$
$\displaystyle(1+\delta)(\sigma_{ij}+v_{i}v_{j})+\Phi\delta_{ij}$
$\displaystyle+$ $\displaystyle\frac{1}{3{\cal
H}^{2}\Omega_{m}}(\nabla_{i}\Phi\nabla_{j}\Phi-\Phi\nabla_{i}\nabla_{j}\Phi+\Phi\delta_{ij}\nabla^{2}\Phi)\,.$

This is of the expected form coming from momentum conservation for free
particles described by a distribution function in general relativity, where
the stress-energy-momentum tensor ${\cal T}^{\mu\nu}=\int
dV_{p}f({\bm{p}})p^{\mu}p^{\nu}$, and $dV_{p}\equiv d^{3}p/\sqrt{-g}\,p^{0}$
is the invariant 3-volume in phase-space and $g$ the metric’s determinant Ma
and Bertschinger (1995). Indeed, in the weak-field limit for non-relativistic
velocities and scalar-only metric perturbations described by $\Phi$, the
covariant derivative in momentum conservation induces the $\Phi$-dependent
terms inside the Newtonian limit $T_{ij}$ that satisfy
$\nabla_{j}T^{\Phi}_{ij}=(1+\delta)\nabla_{i}\Phi$, as required. Combining
Eqs. (80) and (82) gives

$\partial_{\tau}^{2}\delta+{\cal
H}\partial_{\tau}\delta=\nabla_{i}\nabla_{j}T_{ij}\,,$ (84)

which is valid non-linearly and including the stress tensor contribution
$\Sigma_{ij}\equiv(1+\delta)\sigma_{ij}$. Its form also does not change when
allowing for third or higher cumulants. Inserting the perturbative expansion
of the density contrast readily yields the claimed scaling due to the two
overall derivatives in front of $T_{ij}$,

$F_{n,\delta}(\bm{k}_{1},\dots,\bm{k}_{n},\eta)\propto
k^{2},\qquad\bm{k}=\sum_{i}\bm{k}_{i}\to 0\,.$ (85)

This proves that the quadratic scaling of the kernels $F_{n,\delta}$ due to
mass and momentum conservation is a general property, as expected.

Let us now discuss in how far this constraint is compatible with the various
approximation schemes considered in this work. One can check that Eq. (82)
holds independently of whether or not third and higher cumulants are taken
into account, i.e. in particular for both the (cum2) and (cum3+)
approximations. Furthermore, it also holds when including only a subset of the
velocity dispersion tensor perturbations, as is the case in the (sw) or (sv)
approximations. Therefore also the scaling property is preserved for all of
these approximation schemes.

However, we find that _neglecting vorticity breaks momentum conservation_. The
reason for this is that the velocity dispersion tensor contribution to the
Euler equation sources vorticity. Neglecting it discards part of the velocity
field, and therefore also some part of the momentum, leading to artificial
violation of momentum conservation. Therefore, the (s) approximation scheme
should not be used in general, and we do indeed always include vorticity in
this work.

Nevertheless, it’s worth pointing out that when neglecting vorticity, the
spurious linear terms do not affect the one-loop density power spectrum.
Technically, this can be seen by applying the scalar projector
$P^{s}_{ij}=\nabla_{i}\nabla_{j}/\nabla^{2}$ to the Euler equation, and using
this projected equation to derive an equation for the momentum. This would
lead to a contribution to the momentum equation (82) of the form

$(1+\delta)\left[\frac{1}{1+\delta}\partial_{j}[(1+\delta)\epsilon_{ij}]\right]^{s}\,,$
(86)

where the superscript $s$ stands for applying the scalar projector,
$[X_{i}]^{s}\equiv P^{s}_{ik}X_{k}$. This implies that the two factors of
$1+\delta$ in front do not cancel in that case, and therefore do not yield an
expression with an overall derivative in front. This, in turn, would lead to a
violation of the $k^{2}$ scaling when neglecting vorticity. When taking
vorticity into account, one effectively adds the corresponding contribution
with a vector instead of scalar projection, $[X_{i}]^{v}\equiv
P^{v}_{ik}X_{k}$. Due to $P^{v}_{ij}+P^{s}_{ij}=\delta_{ij}$, the projectors
drop out and one obtains again Eq. (82), i.e. a $k^{2}$ scaling of the density
kernels. However, when neglecting vorticity, one obtains an extra term on the
right-hand side of Eq. (84) given by

$\nabla_{i}\left[\delta\left([\nabla_{j}\sigma_{ij}]^{s}-\nabla_{j}\sigma_{ij}+\left[\sigma_{ij}\nabla_{j}A\right]^{s}-\sigma_{ij}\nabla_{j}A\right)\right]\,,$
(87)

which only contains a single overall derivative and would therefore
potentially lead to terms linear in $\bm{k}$ in the large scale limit of
$F_{n,\delta}$ when neglecting vorticity. The first two summands cancel for
the scalar contributions to the stress tensor, and are at least of second
order for the vector and tensor contributions. The last two summands are at
least of second order as well. Therefore, taking the additional factor
$\delta$ in front of the parenthesis into account, the spurious linear scaling
in $\bm{k}$ can only appear starting at third order in perturbation theory,
i.e. for $F_{3,\delta}$. Furthermore, due to the property
$F_{n,\delta}(\bm{k}_{1},\dots,\bm{k}_{n},\eta)=F_{n,\delta}(-\bm{k}_{1},\dots,-\bm{k}_{n},\eta)$,
the kernel $F_{3,\delta}(\bm{k},\bm{q},-\bm{q},\eta)$ that appears in the one-
loop power spectrum is odd under $\bm{k}\to-\bm{k}$, and therefore any linear
term in $\bm{k}$ has to cancel.

Thus, the spurious linear scaling that occurs when neglecting vorticity can
only affect the power spectrum starting at two-loop order, e.g. via the kernel
$F_{3,\delta}(\bm{k}-\bm{p}-\bm{q},\bm{p},\bm{q},\eta)$. For this reason we
include the (s) approximation scheme in our analysis below at one-loop order
for illustration and in order to compare to the approximation schemes (sw),
(sv) and (svt), that are all compatible with momentum conservation. However,
we do not use the (s) scheme when going beyond one-loop.

For illustration, we show the kernel
$F_{3,\delta}(\bm{k}-\bm{p}-\bm{q},\bm{p},\bm{q},\eta)$ in the upper panel of
Fig. 5. For the (sw) and (sv) approximations, that take vorticity into
account, the kernel scales as $k^{2}$ for small $k$. However, in the (s)
approximation scheme, a spurious linear scaling with $\propto k$ appears. Note
that, similarly as before, the kernels are strongly suppressed for large $k$
as compared to SPT, and also differ for small $k$ due to the finite, fixed
size of $\bm{p}$ and $\bm{q}$. The behaviour is qualitatively similar for the
(cum2) and (cum3+) approximations, respectively, and we therefore only show
the latter.

Finally, we comment on why neglecting vorticity within the SPT approximation
does not lead to spurious linear terms. Within SPT, there is no vorticity
source term in the Euler equation. The remaining non-linear term in the Euler
equation has the property
$[v^{s}_{j}\partial_{j}v^{s}_{i}]^{s}=v^{s}_{j}\partial_{j}v^{s}_{i}$, which
implies the well-known property that no vorticity is generated in the perfect
fluid approximation. This has the consequence that, when neglecting the
velocity dispersion tensor altogether, the omission of vorticity does not
spoil momentum conservation. Therefore, the $k^{2}$ scaling of the density
kernels holds in SPT even in absence of vorticity. However, this is not true
any longer as soon as taking the velocity dispersion tensor into account, in
which case it is mandatory to include also vorticity in order to guarantee
that momentum conservation holds.

####  III.2.2 Velocity divergence

In contrast to the density, the scaling of the velocity divergence kernels is
indeed modified when including the velocity dispersion tensor and higher
cumulants along with vorticity,

$F_{n,\theta}(\bm{k}_{1},\dots,\bm{k}_{n},\eta)\propto
k,\qquad\bm{k}=\sum_{i}\bm{k}_{i}\to 0\,.$ (88)

This scaling is in accordance with mass and momentum conservation, i.e.
allowed by the symmetries of the system. However, we observe that for $n=2$ it
does not occur, i.e. $F_{2,\theta}\propto k^{2}$, independently of the
approximation scheme. Furthermore, due to
$F_{n,\theta}(\bm{k}_{1},\dots,\bm{k}_{n},\eta)=F_{n,\theta}(-\bm{k}_{1},\dots,-\bm{k}_{n},\eta)$
the linear term cancels in the one-loop kernel
$F_{3,\theta}(\bm{k},\bm{q},-\bm{q},\eta)\propto k^{2}$. Therefore, the linear
scaling does not affect the one-loop power spectrum of the velocity
divergence, or its cross spectrum with the density contrast. The simplest
kernel that displays linear scaling is
$F_{3,\theta}(\bm{k}-\bm{p}-\bm{q},\bm{p},\bm{q},\eta)$, for
$\bm{p}+\bm{q}\not=0$. This kernel is shown in the lower panel of Fig. 5. In
contrast to the density kernel, we find linear scaling $\propto k$ in the
limit $k\to 0$ also when taking vorticity into account, i.e. for the (sw) and
(sv) approximation schemes. Furthermore, the linear scaling occurs both in the
(cum2) and (cum3+) approximation schemes (only the latter is shown in Fig. 5).

Proceeding analogously as in Sec. II.3 we find the following analytical result
for the ${\cal O}(\epsilon)$ contribution to the third-order velocity
divergence kernel, when keeping the leading term in a Taylor expansion in $k$,

$\displaystyle
F_{3,\theta}(\bm{k}-\bm{p}-\bm{q},\bm{p},\bm{q},\eta)=G_{3}(\bm{k}-\bm{p}-\bm{q},\bm{p},\bm{q})$
(89) $\displaystyle-$
$\displaystyle\left[\frac{(\bm{k}\cdot\bm{p})(\bm{q}\cdot(\bm{q}+2\bm{p}))s_{pq}^{2}}{(\bm{p}+\bm{q})^{2}}+(\bm{p}\leftrightarrow\bm{q})\right]$
$\displaystyle\times$
$\displaystyle\left(\frac{10}{7}E_{2}(\eta)+\frac{2}{21}E_{3}(\eta)-2E_{5/2}(\eta)+\frac{10}{21}E_{7/2}(\eta)\right)$
$\displaystyle+{\cal}{O}(\epsilon^{2},k^{2})\,,$

where $s_{pq}^{2}\equiv 1-(\bm{p}\cdot\bm{q})^{2}/(p^{2}q^{2})$. The second
line scales linearly with $\bm{k}$, and we recall that $E_{n}={\cal
O}(\epsilon)$. Furthermore, one can check that the second line vanishes in the
limit $\bm{p}+\bm{q}\to 0$, giving a $k^{2}$ scaling for the third-order
kernel entering the one-loop power spectrum, as stated above.

Note that, at first order in $\epsilon$, third and higher cumulants do not
contribute, and therefore (89) is independent of the truncation, such that it
is identical for (cum2) and (cum3+), and would not change when going to higher
cumulant order. Furthermore, we find that tensor modes of the stress tensor do
not contribute to the linear term in $k$, such that it holds for both the (sv)
approximation as well as the full (svt) case. (In (sw) approximation only the
prefactor would change). From this we conclude that the linear scaling with
$\bm{k}$ for the velocity divergence in the limit $\bm{k}\to 0$ cannot be an
artifact of the approximation scheme, but is present in the full Vlasov
theory.

As a further check, we demonstrate that the linear term cancels when
considering the momentum $P_{i}$. The divergence of the momentum field can be
written as

$\nabla_{i}P_{i}=-{\cal
H}f\left(\theta+\nabla_{i}\left(\delta\frac{\nabla_{i}\theta}{\nabla^{2}}\right)-\nabla_{i}\left(\delta\frac{\varepsilon_{ijk}\nabla_{j}w_{k}}{\nabla^{2}}\right)\right)\,.$
(90)

Due to the momentum conservation equation (82), we know that the non-linear
kernels obtained when perturbatively expanding $\nabla_{i}P_{i}$ have to scale
as $k^{2}$ in the limit where the sum of all wavenumbers goes to zero. In
particular, at third order this requires that in the combination of kernels

$\displaystyle
F_{3,\theta}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3},\eta)+\frac{1}{3}\Bigg{[}\Bigg{(}\frac{\bm{k}\cdot\bm{k}_{3}}{k_{3}^{2}}F_{2,\delta}(\bm{k}_{1},\bm{k}_{2},\eta)F_{1,\theta}(\bm{k}_{3},\eta)$
$\displaystyle-$
$\displaystyle\frac{(\bm{k}\times(\bm{k}_{2}+\bm{k}_{3}))_{i}}{(\bm{k}_{1}+\bm{k}_{3})^{2}}F_{1,\delta}(\bm{k}_{1},\eta)F_{2,w_{i}}(\bm{k}_{2},\bm{k}_{3},\eta)$
$\displaystyle+$
$\displaystyle\frac{\bm{k}\cdot(\bm{k}_{2}+\bm{k}_{3})}{(\bm{k}_{2}+\bm{k}_{3})^{2}}F_{1,\delta}(\bm{k}_{1},\eta)F_{2,\theta}(\bm{k}_{2},\bm{k}_{3},\eta)\Bigg{)}+2\,\text{perm.}\Bigg{]}\,,$

the linear term in $\bm{k}=\sum_{i}\bm{k}_{i}$ for $\bm{k}\to 0$ cancels.
Using the analytical results for the density, velocity divergence and
vorticity kernels at ${\cal O}(\epsilon)$ obtained in Sec. II.3 and Sec. IV.2
as well as Eq. (89) this can indeed be shown to be the case. This provides a
non-trivial check of the linear scaling for the velocity divergence kernels.

Since the leading contribution to the power spectrum affected by the linear
scaling in Eq. (89) is the two-loop velocity divergence power spectrum, it is
difficult to notice its impact in practice. In particular, note that the
propagator-like terms (obtained from correlating one linear and one non-linear
field) contain kernels
$F_{n,\theta}(\bm{k},\bm{q}_{1},-\bm{q}_{1},\bm{q}_{2},-\bm{q}_{2},\dots,\eta)$
where the linear terms in $\bm{k}$ drop out for symmetry reasons, at all loop
orders. These propagator-like contributions therefore yield a contribution to
the power spectrum $P_{\theta\theta}(k,\eta)$ that scales as $k^{2}P_{0}(k)$
for $k\to 0$, similarly as in SPT. However, the linear scaling of
$F_{3,\theta}(\bm{k}-\bm{p}-\bm{q},\bm{p},\bm{q},\eta)$ produces an additional
contribution arising from mode-coupling terms starting at two-loop order. It
yields a contribution to $P_{\theta\theta}(k,\eta)$ that scales as $k^{2}$ and
should therefore dominate over $k^{2}P_{0}(k)$ on large enough scales.
Nevertheless, within this regime, $P_{\theta\theta}(k,\eta)$ is dominated by
the linear contribution to the power spectrum, and thus the anomalous behavior
is hard to verify in practice. However, as we will see below, a similar
effects occurs for the vorticity, where it is potentially visible given that
the vorticity power spectrum has no linear contribution.

##  IV Generation of vorticity, vector and tensor modes

It is well known that within the perfect, pressureless fluid approximation
vorticity decays as $(\nabla\times{\bf v})\propto 1/a$ within the linear
regime Peebles (1980), and is not generated by non-linear evolution due to the
neglect of orbit crossing, in contradiction to the physical expectation. While
the former property essentially remains true when including velocity
dispersion and higher cumulants, the latter does not Pueblas and Scoccimarro
(2009). In this section we discuss the generation of vorticity as well as
vector and tensor modes of the stress tensor via non-linear effects within the
framework underlying the present work. The generation of vorticity and vector
modes is coupled to each other, and we discuss them jointly in the following.
The tensor mode is discussed in Sec. IV.4.

Let us first recall the evolution equation (23) for the vorticity ${\bf w}$,
that is coupled to the vector mode ${\bf\nu}$ of the velocity dispersion
tensor,

$\displaystyle\partial_{\eta}w_{k,i}+\left(\frac{3}{2}\frac{\Omega_{m}}{f^{2}}-1\right)w_{k,i}+k^{2}\nu_{k,i}$
$\displaystyle=$
$\displaystyle\int_{pq}\gamma_{w_{i}bc}\psi_{p,b}\psi_{q,c}\,,$
$\displaystyle\partial_{\eta}\nu_{k,i}+2\left(\frac{3}{2}\frac{\Omega_{m}}{f^{2}}-1\right)\nu_{k,i}-\epsilon
w_{k,i}$ $\displaystyle=$
$\displaystyle\int_{pq}\gamma_{\nu_{i}bc}\psi_{p,b}\psi_{q,c}\,,$

with vertices listed below.

###  IV.1 Linear approximation

Within the linear approximation, one obtains a coupled set of ordinary
differential equations for each component of the vorticity and ${\bf\nu}$
vector, respectively. In the limit $\epsilon\to 0$, and setting
$\Omega_{m}/f^{2}\mapsto 1$ for simplicity, the second equation has the
approximate solution
$\nu_{k,i}(\eta)=\nu_{k,i}(\eta_{0})e^{-(\eta-\eta_{0})}$, being a purely
decaying mode. Inserting it in the linearized vorticity equation yields

$\displaystyle w_{k,i}(\eta)$ $\displaystyle=$
$\displaystyle(w_{k,i}(\eta_{0})-2k^{2}\nu_{k,i}(\eta_{0}))e^{-\frac{1}{2}(\eta-\eta_{0})}$
(93) $\displaystyle{}+2k^{2}\nu_{k,i}(\eta_{0})e^{-(\eta-\eta_{0})}\,,$

being a superposition of two decaying modes, that corresponds to a decay
$\propto 1/a$ and $\propto a^{-3/2}$ for the unrescaled vorticity $-f{\cal
H}{\bf w}$. Therefore, the vector mode of the dispersion tensor leads to an
additional decaying mode in the vorticity at linear level. By inserting the
solution for ${\bf w}$ from above into the linearized equation for the vector
mode $\nu$, one can derive corrections at first order in $\epsilon$, that are
relatively suppressed by factors of $\epsilon k^{2}$ as long as this quantity
is small, i.e. at early times. In the opposite limit the coupled set of
linearized equations needs to be solved, and the background dispersion leads
to suppression as well. For example, for a power-law time dependence
$\epsilon=\epsilon_{0}e^{\alpha\eta}$, a closed-form solution of the linear
equations can be found in terms of hypergeometric functions (assuming
$\alpha>1/2$), given by

$\displaystyle w_{k,i}(\eta)$ $\displaystyle=$ $\displaystyle
A\,e^{-\frac{1}{2}\eta}\
{}_{0}F_{1}\left(1+\frac{1}{2\alpha};-\frac{\epsilon(\eta)k^{2}}{\alpha^{2}}\right)$
(94) $\displaystyle{}+B\,e^{-\eta}\
{}_{0}F_{1}\left(1-\frac{1}{2\alpha};-\frac{\epsilon(\eta)k^{2}}{\alpha^{2}}\right)\,,$

with coefficients $A,B$ related to the initial vorticity and vector mode. For
$\epsilon k^{2}\ll 1$ one recovers Eq. (93) since ${}_{0}F_{1}\to 1$ in that
limit. For $\epsilon k^{2}\gg 1$ one finds oscillating behavior with a damped
envelope,

$\displaystyle w_{k,i}(\eta)$ $\displaystyle\to$ $\displaystyle
e^{-\frac{3+\alpha}{4}\eta}\,\Bigg{[}A^{\prime}\sin\left(2\sqrt{\epsilon}k+\frac{\alpha-1}{4\alpha}\pi\right)$
(95)
$\displaystyle{}+B^{\prime}\sin\left(2\sqrt{\epsilon}k+\frac{\alpha+1}{4\alpha}\pi\right)\Bigg{]}\,,$

where

$\displaystyle A^{\prime}$ $\displaystyle\equiv$ $\displaystyle
A\left(\frac{\epsilon_{0}k^{2}}{\alpha^{2}}\right)^{-\frac{1}{4}-\frac{1}{4\alpha}}\frac{\Gamma\left(1+\frac{1}{2\alpha}\right)}{\sqrt{\pi}}\,,$
$\displaystyle B^{\prime}$ $\displaystyle\equiv$ $\displaystyle
B\left(\frac{\epsilon_{0}k^{2}}{\alpha^{2}}\right)^{-\frac{1}{4}+\frac{1}{4\alpha}}\frac{\Gamma\left(1-\frac{1}{2\alpha}\right)}{\sqrt{\pi}}\,.$
(96)

Since the linear solutions correspond to decaying modes, initial vorticity or
$\nu$ modes can be neglected. Indeed, as stated already above, this means that
the first-order kernels vanish,

$F_{1,w_{i}}=F_{1,\nu_{i}}=0\,.$ (97)

This implies in particular that the linear contribution to the vorticity and
vector power spectra is zero.

###  IV.2 Vorticity generation at second order

When including velocity dispersion and higher cumulants, vorticity is
generated at the non-linear level Pueblas and Scoccimarro (2009). The leading
contribution arises at second order in perturbation theory. Two scalar modes
can generate vorticity as well as vector perturbations of the stress tensor.

At second order, the evolution equations for the ${\bf w}$ and $\nu$ kernels
read

$\displaystyle\left(\partial_{\eta}+\frac{5}{2}\right)$ $\displaystyle
F_{2,w_{i}}(\bm{p},\bm{q},\eta)=$
$\displaystyle-(\bm{p}+\bm{q})^{2}F_{2,\nu_{i}}(\bm{p},\bm{q},\eta)$
$\displaystyle+\gamma_{w_{i}bc}(\bm{p},\bm{q})F_{1,b}(p,\eta)F_{1,c}(q,\eta)\,,$
$\displaystyle\left(\partial_{\eta}+3\right)$ $\displaystyle
F_{2,\nu_{i}}(\bm{p},\bm{q},\eta)=$
$\displaystyle\epsilon(\eta)F_{2,w_{i}}(\bm{p},\bm{q},\eta)$
$\displaystyle+\gamma_{\nu_{i}bc}(\bm{p},\bm{q})F_{1,b}(p,\eta)F_{1,c}(q,\eta)\,.$

The list of potential vertices for vorticity is $\gamma_{w_{i}Ag}$,
$\gamma_{w_{i}A\epsilon}$, $\gamma_{w_{i}A\nu_{j}}$, $\gamma_{w_{i}At_{jk}}$,
$\gamma_{w_{i}\theta w_{j}}$, $\gamma_{w_{i}w_{j}w_{k}}$, and for the vector
mode $\gamma_{\nu_{i}\theta g}$, $\gamma_{\nu_{i}\theta\epsilon}$,
$\gamma_{\nu_{i}w_{j}g}$, $\gamma_{\nu_{i}w_{j}\epsilon}$,
$\gamma_{\nu_{i}\theta\nu_{j}}$, $\gamma_{\nu_{i}\theta t_{jk}}$,
$\gamma_{\nu_{i}w_{j}\nu_{k}}$, $\gamma_{\nu_{i}w_{j}t_{k\ell}}$. Since only
scalar modes have non-zero $F_{1,a}$, only the vertices $\gamma_{w_{i}Ag}$,
$\gamma_{w_{i}A\epsilon}$ and $\gamma_{\nu_{i}\theta g}$,
$\gamma_{\nu_{i}\theta\epsilon}$ contribute at second order, and describe the
generation of vorticity and vector modes from a pair of scalar perturbations
at non-linear level. The leading contribution to the vorticity power spectrum
is therefore given by

$\displaystyle P_{w_{i}w_{i}}^{1L}(k,\eta)$ $\displaystyle=$ $\displaystyle
2e^{4\eta}\int d^{3}q\,F_{2,w_{i}}(\bm{k}-\bm{q},\bm{q},\eta)$
$\displaystyle{}\times
F_{2,w_{i}}(\bm{k}-\bm{q},\bm{q},\eta)P_{0}(|\bm{k}-\bm{q}|)P_{0}(q)\,.$

The power spectra of the vector mode and their cross spectrum are given by
analogous expressions. Below we will see that it is actually necessary to
include also the two-loop contribution to capture the correct scaling for
$k\to 0$, see Sec. IV.3 and Sec. VI.4, but first discuss the one-loop part.

####  IV.2.1 Analytical results for $\epsilon\to 0$

Figure 6:  Non-linear kernels describing the generation of vorticity
$F_{n,w_{i}}$ (left column) and the vector mode of the stress tensor
$k^{2}F_{n,\nu_{i}}$ (right column), from non-linear coupling of two ($n=2$,
upper row) or three ($n=3$, lower row) perturbation modes of the initial
linear density field. The kernels are shown for various approximation schemes
as indicated in the legend (see Table 2). For $n=2$ we also show the
analytical results at first order in $\epsilon$ (black solid line), as well as
a line $\propto k^{2}$ (black dashed) for comparison. Note that $F_{2,w_{i}}$
points in the direction perpendicular to the plane spanned by $\bm{k}$ and
$\bm{q}$, and we show its projection on
$\bm{k}\times\bm{q}/|\bm{k}\times\bm{q}|$. For $n=3$ we show two components
$F_{3,w_{i}}$, $i=1,2$, perpendicular to $\bm{k}$ (see App. D for details),
finding $F_{3,w_{i}}\propto k$ for $k\to 0$, similarly to $\theta$ (see Sec.
III.2). The impact of the vector mode on vorticity ((sw) vs (svt) in the left
panels) is noticeable at $k\gtrsim k_{\sigma}$, while the impact of third
cumulant perturbations is relatively mild for vorticity as well as vector mode
kernels ((cum2) vs (cum3+) in the upper row). The analytical results Eqs.
(100, 102) at first order in $\epsilon$ are shown as thin black lines, and
agree well for $k\ll k_{\sigma}$. In SPT all these kernels are zero. For the
figure we chose $\epsilon=\epsilon_{0}\,e^{\alpha\eta}$ with $\alpha=4/3$, and
$c_{kp}=0.3,c_{kq}=0.5,c_{pq}=0.875,p=k_{\sigma},q=0.5k_{\sigma}$, where
$k_{\sigma}=1/\sqrt{\epsilon_{0}}$.

Let us first discuss an approximate analytical solution for the second-order
kernels in the limit $\epsilon\to 0$, similar as above for the scalar kernels.
Both vorticity and the vector mode are generated starting at first order in
$\epsilon$. Therefore, the $\epsilon F_{2,w_{i}}$ term in the equation for the
vector mode can be neglected at lowest order in $\epsilon$, and the equation
for the vector kernel can be solved independently at this order. Using
furthermore that $F_{1,\delta\epsilon}$ starts at order $\epsilon^{2}$, only
the vertex $\gamma_{\nu_{i}\theta g}$ is relevant, and we find

$\displaystyle F_{2,\nu_{i}}(\bm{p},\bm{q},\eta)$ $\displaystyle=$
$\displaystyle 2(\gamma_{\nu_{i}\theta g}(\bm{p},\bm{q})+\gamma_{\nu_{i}\theta
g}(\bm{q},\bm{p}))$ (100)
$\displaystyle{}\times\int^{\eta}d\eta^{\prime}e^{3(\eta^{\prime}-\eta)}\int^{\eta^{\prime}}d\eta^{\prime\prime}\,e^{2(\eta^{\prime\prime}-\eta^{\prime})}\epsilon(\eta^{\prime\prime})$
$\displaystyle{}+{\cal O}(\epsilon^{2})\,.$

Inserting this result in the equation for vorticity, and using that only
$\gamma_{w_{i}Ag}$ contributes at first order in $\epsilon$, one finds

$\displaystyle F_{2,w_{i}}(\bm{p},\bm{q},\eta)$ $\displaystyle=$
$\displaystyle-(\bm{p}+\bm{q})^{2}\int^{\eta}d\eta^{\prime}e^{\frac{5}{2}(\eta^{\prime}-\eta)}\,F_{2,\nu_{i}}(\bm{p},\bm{q},\eta^{\prime})$
(101)
$\displaystyle{}+2(\gamma_{w_{i}Ag}(\bm{p},\bm{q})+\gamma_{w_{i}Ag}(\bm{q},\bm{p}))$
$\displaystyle{}\times\int^{\eta}d\eta^{\prime}e^{\frac{5}{2}(\eta^{\prime}-\eta)}\int^{\eta^{\prime}}d\eta^{\prime\prime}\,e^{2(\eta^{\prime\prime}-\eta^{\prime})}\epsilon(\eta^{\prime\prime})$
$\displaystyle{}+{\cal O}(\epsilon^{2})\,.$

Inserting the explicit form of the vertices Garny  _et al._ (2022), one finds

$F_{2,w_{i}}(\bm{p},\bm{q},\eta)=(\bm{p}\times\bm{q})_{i}(\bm{p}\cdot\bm{q})\left(\frac{1}{p^{2}}-\frac{1}{q^{2}}\right)J^{w}(\eta)+{\cal
O}(\epsilon^{2})\,,$ (102)

with

$\displaystyle J^{w}(\eta)$ $\displaystyle=$
$\displaystyle\int^{\eta}d\eta^{\prime}e^{\frac{5}{2}(\eta^{\prime}-\eta)}\Bigg{[}\int^{\eta^{\prime}}d\eta^{\prime\prime}\,e^{2(\eta^{\prime\prime}-\eta^{\prime})}\epsilon(\eta^{\prime\prime})$
(103)
$\displaystyle+\int^{\eta^{\prime}}d\eta^{\prime\prime}e^{3(\eta^{\prime\prime}-\eta^{\prime})}\int^{\eta^{\prime\prime}}d\eta^{\prime\prime\prime}\,e^{2(\eta^{\prime\prime\prime}-\eta^{\prime\prime})}\epsilon(\eta^{\prime\prime\prime})\Bigg{]}$
$\displaystyle=$ $\displaystyle 4E_{2}(\eta)-6E_{5/2}(\eta)+2E_{3}(\eta)\,.$

Note that the contribution in the second line would be missed when neglecting
the impact of the vector mode of the dispersion tensor. As expected, the
vorticity and vector mode are perpendicular to the plane spanned by the
wavevectors ${\bf p}$ and ${\bf q}$. Furthermore,
$F_{2,w_{i}}(\bm{k}-\bm{q},\bm{q},\eta)$ scales as $k^{2}$ for $\bm{k}\to 0$,
which can be readily checked for the analytical result at order $\epsilon$
from above. Inserting this first order result in $\epsilon$ into Eq. (IV.2)
yields a dependence on wavevectors identical to the integrand of Eq. (70) in
Cusin  _et al._ (2017); in our numerical treatment (see below) we do not
expand in $\epsilon$ and also fully include the impact of the vector mode on
vorticity, while Eq. (102) captures only the lowest order in $\epsilon$.

####  IV.2.2 Numerical results

As for the scalar kernels, we solve the differential equation (IV.2) for the
vorticity and vector kernels numerically in our analysis, thereby resumming
terms of all powers in the background dispersion $\epsilon$. The numerical
results for a particular configuration of $F_{2,w_{i}}$ and $F_{2,\nu_{i}}$
are shown in the upper left and right panels of Fig. 6, respectively. The
numerical and analytical results agree well in the limit $k\ll k_{\sigma}$, as
expected. Furthermore, the scaling
$F_{2,w_{i}}(\bm{k}-\bm{q},\bm{q},\eta)\propto k^{2}$ for $\bm{k}\to 0$ holds
also for the numerical result. For the vector mode, we show the dimensionless
quantity $k^{2}F_{2,\nu}$, which also scales as $k^{2}$, meaning that
$F_{2,\nu}$ approaches a constant, analogously as for the scalar modes of the
dispersion tensor. The simplest approximation capturing vorticity generation
is (sw), for which the vector mode $\nu$ is neglected. As can be seen in the
left upper panel of Fig. 6, this approximation leads to a small shift of
$F_{2,w_{i}}$ for $k\lesssim k_{\sigma}$ as compared to the case (svt) with
vector mode, and a more shallow decline for $k\gg k_{\sigma}$. Note that at
second order the tensor mode is irrelevant for $F_{2,w_{i}}$ and
$F_{2,\nu_{i}}$, such that there is no difference between (sv) and (svt), and
only the latter is shown. This is also the reason why only the (svt) result is
shown for $F_{2,\nu_{i}}$ in the right upper panel of Fig. 6, noting in
addition that $\nu$ is not taken into account within the (sw) scheme. Finally,
the upper panels show the difference between (cum2) and (cum3+), both
evaluated for (svt). The impact of higher cumulants is negligibly small for
$k\lesssim k_{\sigma}$, while it somewhat impacts the asymptotic behaviour for
$k\gg k_{\sigma}$.

###  IV.3 Vorticity generation at third order

Similarly as for the velocity divergence, we find that the non-linear kernels
for the vorticity in general scale only with the first power of the total
wavenumber,

$F_{n,w_{i}}(\bm{k}_{1},\dots,\bm{k}_{n},\eta)\propto
k,\qquad\bm{k}=\sum_{i}\bm{k}_{i}\to 0\,.$ (104)

As for $\theta$, this linear scaling appears first at third order, i.e. for
$n\geq 3$, while $F_{2,w_{i}}\propto k^{2}$ as discussed above.

In the lower left panel of Fig. 6 we show the kernel
$F_{3,w_{i}}(\bm{k}-\bm{p}-\bm{q},\bm{p},\bm{q},\eta)$ versus $k$ in the (sw)
and (sv) approximations, and for two components of $w_{i}$ perpendicular to
${\bf k}$ (see App. D for details on the basis choice and projection). The
linear scaling in $k$ can clearly be seen for all cases. As for $\theta$, the
scaling would be $\propto k^{2}$ for $\bm{p}+\bm{q}=0$. This implies that the
vorticity power spectrum scales as

$P_{w_{i}w_{i}}(k,\eta)\propto k^{2}\,,$ (105)

for $k\ll k_{\sigma}$. However, since $F_{2,w_{i}}\propto k^{2}$, _this
scaling can only be observed starting at two-loop order_ , while
$P_{w_{i}w_{i}}^{1L}(k,\eta)\propto k^{4}$ at one-loop. We will confirm this
behavior against measurements of the vorticity power spectrum in numerical
simulations in Sec. VI.

In contrast to the vorticity, the kernels $F_{n,\nu_{i}}$ approach constant
values for $k\to 0$, such that the dimensionless combination
$k^{2}F_{n,\nu_{i}}$ scales as $k^{2}$ also for $n=3$. This is confirmed by
our numerical results shown in the lower right panel of Fig. 6.

###  IV.4 Generation of tensor modes

Figure 7:  Second order kernel of the tensor perturbation $t_{ij}$ of the
velocity dispersion that describes the generation of tensor modes from
coupling two scalar modes, captured by . We show the dimensionless kernel
$k^{2}F_{2,t_{\text{eff}}}(\bm{k}-\bm{q},\bm{q})$ versus $k$ in two
approximations (see legend), as well as the analytical result at first order
in $\epsilon$ (black solid line). The kernel scales as $k^{2}$ for $k\to 0$
(black dashed line).

As discussed above, tensor perturbations of the stress tensor have only a very
small impact on the non-linear kernels for the density contrast and velocity
divergence (see Fig. 3), as opposed to impact of vorticity and vector modes on
the density. Therefore, tensor modes could be neglected in practice when being
interested in density and velocity power spectra (see Sec. VI). Nevertheless,
as a consequence of non-linear mode coupling, a power spectrum of tensor modes
is generated, analogously to the case of vorticity and vector modes discussed
above. The leading contribution arises from the second-order kernel,

$\displaystyle\left(\partial_{\eta}+3\right)$ $\displaystyle
F_{2,t_{ij}}(\bm{p},\bm{q},\eta)=\gamma_{t_{ij}bc}(\bm{p},\bm{q})F_{1,b}(p,\eta)F_{1,c}(q,\eta)\,,$

giving rise to a tensor power spectrum at one-loop,

$\displaystyle P_{t_{ij}t_{ij}}^{1L}(k,\eta)$ $\displaystyle=$ $\displaystyle
2e^{4\eta}\int d^{3}q\,F_{2,t_{ij}}(\bm{k}-\bm{q},\bm{q},\eta)$
$\displaystyle{}\times
F_{2,t_{ij}}(\bm{k}-\bm{q},\bm{q},\eta)P_{0}(|\bm{k}-\bm{q}|)P_{0}(q)\,.$

Vertices that could potentially contribute (disregarding third and higher
cumulant perturbations) are Garny  _et al._ (2022) $\gamma_{t_{ij}\theta g}$,
$\gamma_{t_{ij}\theta\epsilon}$, $\gamma_{t_{ij}w_{k}g}$,
$\gamma_{t_{ij}w_{k}\epsilon}$, $\gamma_{t_{ij}\theta\nu_{k}}$,
$\gamma_{t_{ij}\theta t_{k\ell}}$, $\gamma_{t_{ij}w_{k}\nu_{\ell}}$,
$\gamma_{t_{ij}w_{k}t_{\ell m}}$. Since $F_{1,b}$ and $F_{1,c}$ are non-zero
only for scalar modes, only the vertices contributing to $F_{2,t_{ij}}$ are

$\displaystyle\gamma_{t_{ij}\theta g}(\bm{p},\bm{q})$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\frac{(\bm{p}\times\bm{q})^{2}\,\bm{p}\cdot\bm{q}}{2(\bm{p}+\bm{q})^{2}p^{2}q^{2}}{\cal
T}_{ij}\,,$ $\displaystyle\gamma_{t_{ij}\theta\epsilon}(\bm{p},\bm{q})$
$\displaystyle=$
$\displaystyle\frac{1}{2}\frac{(\bm{p}\times\bm{q})^{2}}{(\bm{p}+\bm{q})^{2}p^{2}}{\cal
T}_{ij}\,,$ (108)

where we defined

${\cal
T}_{ij}\equiv\delta_{ij}-\frac{(\bm{p}+\bm{q})_{i}(\bm{p}+\bm{q})_{j}}{(\bm{p}+\bm{q})^{2}}-2\frac{(\bm{p}\times\bm{q})_{i}(\bm{p}\times\bm{q})_{j}}{(\bm{p}\times\bm{q})^{2}}\,.$
(109)

Note that $0={\cal T}_{ii}=(\bm{p}+\bm{q})_{i}{\cal T}_{ij}={\cal
T}_{ij}(\bm{p}+\bm{q})_{j}$, as required for a tensor mode.

Let us first discuss the analytical result when expanding at first order in
$\epsilon$. In this case, only the vertex $\gamma_{t_{ij}\theta g}$ is
relevant since $F_{1,\delta\epsilon}\propto\epsilon^{2}$ (we remind the reader
that the index $\epsilon$ in $\gamma_{abc}$ stands for the mode
$\delta\epsilon$). Furthermore, at ${\cal O}(\epsilon)$ we use Eq. (45) for
$F_{1,g}$ and can set $F_{1,\theta}\to 1$, giving

$\displaystyle F_{2,t_{ij}}(\bm{p},\bm{q},\eta)$ $\displaystyle=$
$\displaystyle 2(\gamma_{t_{ij}\theta g}(\bm{p},\bm{q})+\gamma_{t_{ij}\theta
g}(\bm{q},\bm{p}))$ (110)
$\displaystyle{}\times\int^{\eta}d\eta^{\prime}e^{3(\eta^{\prime}-\eta)}\int^{\eta^{\prime}}d\eta^{\prime\prime}\,e^{2(\eta^{\prime\prime}-\eta^{\prime})}\epsilon(\eta^{\prime\prime})$
$\displaystyle{}+{\cal O}(\epsilon^{2})\,.$

This analytical result is compared to the full numerical kernel in Fig. 7.
Here we defined

$F_{2,t_{ij}}(\bm{p},\bm{q},\eta)=F_{2,t_{\text{eff}}}(\bm{p},\bm{q},\eta)\,{\cal
T}_{ij}\,,$ (111)

factoring out the dependence on the indices captured by ${\cal T}_{ij}$ (see
App. D for details). The figure shows the dimensionless combination
$k^{2}F_{2,t_{\text{eff}}}(\bm{k}-\bm{q},\bm{q})$, that scales as $k^{2}$ for
$k\to 0$, analogously to the vector perturbation of the stress tensor. As for
the latter, the same scaling is expected at higher orders, and therefore no
qualitative change is expected when including contributions beyond one-loop
order. Thus, $P_{t_{ij}t_{ij}}(k,\eta)$ approaches a constant for $k\to 0$ at
all loop orders. For the dimensionless tensor power spectrum, this means that
$k^{4}P_{t_{ij}t_{ij}}(k,\eta)\propto k^{4}$ for $k\to 0$.

The tensor modes of the stress tensor generated by non-linearities constitute
a source for a stochastic gravitational wave background. Since their frequency
is related to the wavenumber $k$, its maximal amplitude is expected to occur
for ultra-low frequencies related to the dispersion scale $f_{\text{gw}}\sim
ck_{\sigma}/(2\pi)\simeq 1.5\cdot 10^{-15}\,h$ Hz
$\times(k_{\sigma}/(1h/\text{Mpc}))$. In addition, due to the non-relativistic
velocities of order ${\cal H}\times\sqrt{\epsilon}\simeq
100\text{km}/\text{s}\times\sqrt{\epsilon}/(\text{Mpc}/h)$ the amplitude of
the gravitational wave background is expected to be tiny, although the small
velocities can be somewhat compensated by the large masses that are involved.
For a more detailed assessment, we refer to future work.

##  V Background dispersion from simulations and perturbation theory

In this section we discuss various estimates of the background dispersion
$\epsilon(\eta)$, that plays a central role in . We discuss two complementary
approaches, one based on a self-consistent solution of the equation of motion
for $\epsilon(\eta)$ within perturbation theory (Eq. 13), and one based on
accounting for the dispersion generated in halos following paper I Garny  _et
al._ (2022), and taking the halo profile as well as the halo mass function
measured in numerical N-body simulations into account.

For simplicity, we focus on a scaling universe with linear input power
spectrum

$P_{0}(k)=Ak^{n_{s}}\equiv\frac{1}{4\pi
k^{3}}\left(\frac{k}{k_{\text{nl}}}\right)^{n_{s}+3}\,,$ (112)

with constant power law index $n_{s}$ and amplitude $A$ related to the non-
linear scale $k_{\text{nl}}$, for an EdS background cosmology. Scaling
symmetry ensures that the dimensionless ratio of any dimensionfull quantity to
an appropriate power of the time-dependent non-linear scale
$k_{\text{nl}}(\eta)=k_{\text{nl}}e^{-2\eta/(n_{s}+3)}$ is constant in time.
This implies Garny  _et al._ (2022)

$\epsilon(\eta)=\epsilon_{0}\,e^{\alpha\eta},\quad\alpha=\frac{4}{n_{s}+3}\,.$
(113)

For given $n_{s}$, this leaves a single free parameter $\epsilon_{0}\equiv
1/k_{\sigma}^{2}$ for the background dispersion, that we can conveniently
parameterize by the ratio

$k_{\sigma}/k_{\text{nl}}\equiv(\epsilon_{0}k_{\text{nl}}^{2})^{-1/2}\,.$
(114)

We note that when defining a time-dependent dispersion scale by
$\epsilon(\eta)\equiv 1/k_{\sigma}^{2}(\eta)$, the time dependence given above
indeed implies that
$k_{\sigma}(\eta)/k_{\text{nl}}(\eta)=k_{\sigma}/k_{\text{nl}}$ is constant.
Similarly, scaling symmetry implies a constant dimensionless ratio
$\bar{\omega}\equiv\omega(\eta)/\epsilon(\eta)^{2}$ of the fourth cumulant
expectation value to the dispersion squared.

###  V.1 Self-consistent solution in perturbation theory

Table 3: Self-consistent solutions in linear and one-loop approximation for the velocity dispersion scale $k_{\sigma}=\epsilon_{0}^{-1/2}$ relative to the non-linear scale, for scaling universes with spectral indices $n_{s}=2,1,0$ and in various approximation schemes (see Table 2). Note that vorticity+vector (sv) and tensor (svt) modes enter starting at one-loop only. For (cum3+), we show results for two choices $\bar{\omega}=\bar{\omega}_{\text{fid}}=0.579,0.616,0.668$ for $n_{s}=2,1,0$ and $\bar{\omega}=1$, respectively (see text for details). $n_{s}$ |  $k_{\sigma}/k_{\text{nl}}$ (linear) |  $k_{\sigma}/k_{\text{nl}}$ (one-loop)
---|---|---
| (cum2) | (cum3+) | (cum2) | (cum3+)
| - | - | - | (s) | (sv) | (svt) | (svt) | (svt)
| - | $\bar{\omega}_{\text{fid}}$ | $\bar{\omega}=1$ | - | - | - | $\bar{\omega}_{\text{fid}}$ | $\bar{\omega}=1$
$2$ | $2.34$ | $2.00$ | $1.85$ | $1.80$ | $1.74$ | $1.74$ | $1.62$ | $1.51$
$1$ | $2.55$ | $2.17$ | $2.00$ | $1.96$ | $1.86$ | $1.87$ | $1.69$ | $1.64$
$0$ | $2.97$ | $2.50$ | $2.29$ | $2.28$ | $2.11$ | $2.11$ | $2.13$ | $1.97$

The equation of motion (13) for the background dispersion can for a scaling
universe and EdS background be turned into an algebraic equation for
$k_{\sigma}/k_{\text{nl}}$, that can be written as (see paper I Garny  _et
al._ (2022))

$\frac{\epsilon^{\prime}(\eta)}{\epsilon(\eta)}+1=\frac{1}{3}\left(xI^{\text{lin}}(n_{s})+x^{2}I^{1L}(n_{s})+\dots\right)\,,$
(115)

where $x\equiv(k_{\sigma}/(\sqrt{3}k_{\text{nl}}))^{n_{s}+3}$ and the ellipsis
denote two and higher loop contributions. Using Eq.g (113), the left-hand side
evaluates to a constant, $\epsilon^{\prime}/\epsilon+1=(n_{s}+7)/(n_{s}+3)$.
Furthermore, the coefficients $I^{\text{lin}}(n_{s})$ and $I^{1L}(n_{s})$ are
independent of $x$ by virtue of the scaling symmetry, and are given by
integrals over the sum of the (suitably rescaled) power spectra
$P_{\theta\tilde{g}}=P_{\theta g}-P_{\theta\delta\epsilon}$,
$2P_{w_{i}\nu_{i}}$ and $P_{A\pi}$ evaluated in linear and one-loop
approximation, respectively Garny  _et al._ (2022). These power spectra arise
from the source term, Eq. (20), entering the equation of motion (13) for
$\epsilon(\eta)$. When including one-loop and linear terms on the right-hand
side of Eq. (115) one therefore obtains a quadratic equation for $x$. Its
solution(s) $x_{*}$ determine the self-consistent value for

$k_{\sigma}/k_{\text{nl}}\big{|}_{\text{self-
consistent}}=\sqrt{3}x_{*}^{1/(n_{s}+3)}\,,$ (116)

within perturbation theory up to one-loop order. Correspondingly, dropping
also the term $x^{2}I^{1L}(n_{s})$ yields self-consistent values in linear
approximation. Note that since $I^{\text{lin}}(n_{s})$ and $I^{1L}(n_{s})$ are
independent of $x$, and therefore of $k_{\sigma}$, the power spectra entering
these expressions can be computed at an arbitrary reference value
$k_{\sigma}^{\text{ref}}$ (see App. E for details). We evaluate them in