various approximation schemes (see Table 2). At one-loop, we find that only
one of the roots $x_{*}$ of the quadratic equation (115) is positive, thus
leading only to a single physical solution for $k_{\sigma}/k_{\text{nl}}$.

The results for the self-consistent value of $k_{\sigma}/k_{\text{nl}}$ are
shown in Table 3. Including the one-loop correction leads to an increase in
the background dispersion $\epsilon(\eta)$ relative to the linear
approximation, corresponding to a decrease in $k_{\sigma}$ by about $20-30\%$.
While being significant, the shift is consistent with the one-loop
contribution being a perturbative correction. In addition, we observe that the
self-consistent background dispersion is larger (i.e. $k_{\sigma}$ smaller)
when including third cumulant perturbations (cum3+) compared to the second
cumulant approximation (cum2). For the (cum3+) case, we fix $\bar{\omega}$ to
the value obtained from self-consistently solving the equation of motion of
the fourth cumulant expectation value in linear approximation Garny  _et al._
(2022) as fiducial value (see Table 3), and consider the alternative value
$\bar{\omega}=1$ to assess the dependence on this quantity. A self-consistent
determination of $\bar{\omega}$ at one-loop is beyond the scope of this work,
but would be an interesting extension. Nevertheless, the dependence of
$k_{\sigma}/k_{\text{nl}}$ on $\bar{\omega}$ is relatively mild, such that our
fiducial choice may be considered as a reasonable estimate. Note that
vorticity, vector and tensor modes do not enter in the linear approximation.
At one-loop, taking vorticity and vector modes into account (sv) leads to a
moderate increase in $\epsilon(\eta)$ compared to the case with scalar
perturbations only (s). On the contrary, adding also tensor modes (svt) has a
negligible impact. Finally, we note that we checked the convergence of the
integrals over the power spectra entering $I^{1L}(n_{s})$, with negligible
differences when using $30k_{\sigma}$ versus $40k_{\sigma}$ as cutoff for
$n_{s}=2,1$. For $n_{s}=0$, the difference is at the percent level.

Overall, we find that the self-consistent prediction for the background
dispersion $\epsilon(\eta)$ has the tendency to increase (such that
$k_{\sigma}$ decreases) when enhancing the complexity of the approximation
scheme, i.e. when successively including one-loop, vorticity/vector mode and
higher cumulant effects. The self-consistent determination of $k_{\sigma}$ at
one-loop therefore may be considered as indicative, with the actual value
likely being smaller. In order to bracket the uncertainty on the background
dispersion, we consider in addition an orthogonal approach based on a halo
model in the following.

###  V.2 Dispersion from simulated halos

Let us calculate the velocity dispersion contribution expected from dark
matter halos. If the dispersion is only non-zero inside dark matter halos, we
have

$\epsilon(z)=\mathop{\langle}\epsilon_{ii}(\bm{x},z)/3\mathop{\rangle}=\int\epsilon
P(\epsilon)d\epsilon=\int\epsilon P(\epsilon/h)P(h)dh\,d\epsilon$ (117)

where we used Bayes’ theorem to introduce the conditional probability
$P(\epsilon/h)$, i.e. the probability for $\epsilon(\bm{x},z)$ given that we
are inside a dark matter halo. Since halos are characterized by their mass we
can write, more precisely

$\epsilon(z)=\int_{0}^{\infty}\bar{\epsilon}_{h}(m,z)\,f(m)\,dm\,,$ (118)

where $f(m)$ is the fraction of mass in halos of mass $m$, which plays the
role of $P(h)$, and the conditional average

$\int\epsilon
P(\epsilon/h)\,d\epsilon=\bar{\epsilon}_{h}(m,z)\equiv\Big{(}{3\over 4\pi
r_{\text{vir}}^{3}}\Big{)}\int d^{3}x\ \epsilon_{h,m}(\bm{x},z)\,,$ (119)

where the integral on the RHS is done over the halo. The mass fraction $f(m)$
is related to the mass function $dn/dm$, i.e. the number density of halos of
mass $m$ per unit mass, by $dn/dm\equiv\bar{\rho}f(m)/m$. A convenient
parametrization of the mass function is given by Sheth and Tormen (1999);
Despali  _et al._ (2016)

$\displaystyle m\,f(m)$ $\displaystyle=$ $\displaystyle
A\Big{[}1+{1\over(a\nu)^{p}}\Big{]}\sqrt{a\nu\over 2\pi}{\rm
e}^{-a\nu/2}\,{d\ln\nu\over d\ln m}$ (120) $\displaystyle\equiv$
$\displaystyle\nu f(\nu)\ {d\ln\nu\over d\ln m}\,,$

where $\sqrt{\nu}\equiv\delta_{c}/\sigma(R,z)$, with $\delta_{c}$ the critical
density for collapse at the corresponding time and $\sigma^{2}(R,z)$ the
variance of linear-SPT fluctuations, with $R$ the Lagrangian radius of halos
of mass $m$, i.e. $m=4\pi\bar{\rho}R^{3}/3$ where $\bar{\rho}$ is the total
homogeneous matter density. The constants $a,p,A$ are fit to N-body
simulations. The condition $\nu=1$ defines the non-linear scale $R_{*}$ and
its associated non-linear mass scale $m_{*}$, with $R_{*}$ being the
Lagrangian size of halos of mass $m_{*}$. This mass scale determines when the
abundance of halos is exponentially suppressed from Eq. (120), i.e.
essentially there is not enough time to build significant number of halos more
massive than $m=m_{*}$.

For the application we are interested in this paper, a scale-free universe
with $-1\leq n_{s}\leq 2$ and $\Omega_{m}=1$, $\delta_{c}(z)=1.686$. However,
the usual calculation of $\sigma^{2}(R,z)$ becomes problematic for blue
spectral indices as the standard choice of a top-hat filter gives a divergent
answer for $\sigma^{2}(R,z)$ when $n_{s}\geq 1$. Replacing the linear-SPT
evolution by the linear evolution with dispersion found in paper I Garny  _et
al._ (2022) cures this UV divergence, thus providing an interesting
application of our results. In this case, the variance is given by,

$\sigma^{2}=\Big{(}{k_{\sigma}\over
k_{\text{nl}}}\Big{)}^{n_{s}+3}\int_{0}^{\infty}{dx\over 2}\,W_{\rm
TH}^{2}(\sqrt{x}k_{\sigma}R)\,x^{{n_{s}+1}\over
2}\,F_{1,\delta}^{2}(x,n_{s})\,,$ (121)

where $x\equiv k^{2}\epsilon=(k/k_{\sigma})^{2}$, $W_{\rm TH}$ is the top-hat
window in Fourier space, the linear kernel is given by

$F_{1,\delta}(x,n_{s})\equiv{}_{1}F_{2}\left(\frac{4+\alpha}{3\alpha};1+\frac{2}{\alpha},1+\frac{5}{2\alpha};-\frac{3x}{\alpha^{2}}\right)\,,$
(122)

with ${}_{1}F_{2}$ a hypergeometric function, $\alpha(n_{s})=4/(n_{s}+3)$ as
usual, and $k_{\sigma}/k_{\text{nl}}$ is obtained from the self-consistent
linear solution with dispersion (see Fig. 9 and Table I in paper I Garny  _et
al._ (2022)). For example, setting $\sigma=\delta_{c}$ gives
$k_{\text{nl}}R_{*}=0.83,0.98,1.09,1.15$ for $n_{s}=-1,0,1,2$ respectively.
For reference, using linear-SPT fluctuations gives
$k_{\text{nl}}R_{*}=0.89,1.18$ for $n_{s}=-1,0$.

It is instructive to first consider the contribution to $\bar{\epsilon}_{h}$
in Eq. (118) from an isothermal halo, before we address the more realistic
case of an NFW profile discussed in paper I Garny  _et al._ (2022). In this
case, the contribution to 1D velocity dispersion is simply $Gm/(2r_{\rm vir})$
independent of location, which in the normalization corresponding to
$\bar{\epsilon}_{h}$ means an isothermal halo of mass $m$ contributes to
$\epsilon(z)$ as:

${G\,m\over 2\,r_{\rm vir}\,({\cal H}f)^{2}}={\widehat{m}^{2/3}\over
4f^{2}}\,\frac{\Omega_{m,0}^{2/3}\ \Delta_{\rm
vir}^{1/3}}{\Omega_{m,0}+\Omega_{\Lambda,0}\,a^{3}}\ R_{*}^{2}\,,$ (123)

where we have defined the dimensionless variable $\widehat{m}\equiv m/m_{*}$,
and used that $m=(4\pi/3)\rho_{c,0}\Delta_{\rm vir}r_{\rm vir}^{3}/a^{3}$, the
non-linear mass $m_{*}=(4\pi/3)\rho_{c,0}\Omega_{m,0}R_{*}^{3}$ and assumed a
flat $\Lambda$CDM cosmology that goes to the self-similar case we are mostly
interested in here when $\Omega_{\Lambda,0}=1-\Omega_{m,0}\to 0$. Here
quantities with $0$ subindices correspond to $z=0$, and note that we assume
halos are identified as having a mean physical density of $\Delta_{\rm vir}$
over the critical density $\rho_{c}$ Despali  _et al._ (2016). A simple
estimate is $\Delta_{\rm
vir}(z)=18\pi^{2}+82(\Omega_{m}-1)-39(\Omega_{m}-1)^{2}$ Bryan and Norman
(1998). Since we are interested in predicting $k_{\sigma}/k_{\text{nl}}$ it is
useful to consider

$\Big{(}{G\,m\,k_{\text{nl}}^{2}\over 2\,r_{\rm vir}\,({\cal
H}f)^{2}}\Big{)}^{-1/2}=2f\widehat{m}^{-1/3}\sqrt{\frac{\Omega_{m,0}+\Omega_{\Lambda,0}\,a^{3}}{\Omega_{m,0}^{2/3}\
\Delta_{\rm vir}^{1/3}}}\ {1\over k_{\text{nl}}R_{*}}\,,$ (124)

which is the contribution from an isothermal halo of mass $m$ to
$k_{\sigma}/k_{\text{nl}}$. In the $\Omega_{m}\to 1$ limit, this becomes,

$\Big{(}{G\,m\,k_{\text{nl}}^{2}\over 2\,r_{\rm vir}\,{\cal
H}^{2}}\Big{)}^{-1/2}={2\,\widehat{m}^{-1/3}\over(18\pi^{2})^{1/6}}\ {1\over
k_{\text{nl}}R_{*}}\simeq 0.84\ {\widehat{m}^{-1/3}\over
k_{\text{nl}}R_{*}}\,,$ (125)

which given the values of $k_{\text{nl}}R_{*}$ quoted above says that halos of
order $m\sim m_{*}$ ($\widehat{m}=1$) contribute $k_{\sigma}\sim
k_{\text{nl}}$. More precisely, if we integrate over the mass function we
obtain $(k_{\sigma}/k_{\text{nl}})\simeq 0.71,0.81,0.91,1.04$ for
$n_{s}=2,1,0,-1$ respectively, i.e. a dispersion scale larger than the non-
linear scale. This gives a useful order of magnitude, but isothermal halos
overestimate the velocity dispersion compared to realistic halo profiles.
Therefore, let us now consider NFW halos.

Figure 8:  Mass function $\nu f(\nu)$ fits to simulation measurements (see
App. G) as a function of $\widehat{m}=m/m_{*}$ for $n_{s}=-1,0,1,2$ (top to
bottom from the low-mass end).

From paper I Garny  _et al._ (2022), we have that the 1D velocity dispersion
$m_{2}(r)$ as a function of radius $r$ inside a halo with profile $\rho(r)$
and velocity dispersion tensor anisotropy parameter $\beta$ (assumed to be
constant) is given by

$m_{2}(r)=\frac{(1-2\beta/3)}{\rho\,r^{2\beta}}\int_{r}^{\infty}{d\Phi\over
dr}\rho(r)\,r^{2\beta}\,dr\,,$ (126)

where for the case of NFW halos the potential $\Phi(r)$ can be written as

$\Phi(r)=-{Gm\over r_{\rm vir}}\,c\,g(c)\,\Big{[}{\ln(1+x)\over x}-{1\over
1+c}\Big{]}\,,$ (127)

where $x\equiv c\,r/r_{\rm vir}$, $c$ is the concentration parameter and
$g^{-1}(c)\equiv\ln(1+c)-c/(1+c)$. According to simulations, the concentration
of a halo is a weak function of halo mass, we take
$c(\widehat{m})=9~{}\widehat{m}^{-0.13}$ Bullock  _et al._ (2001). A more
recent study that includes scale-free simulations with red spectral indices
$n_{s}\leq-1$ builds a more accurate model with six free parameters Diemer and
Kravtsov (2015), but their own comparison (see their Fig. 10) shows small
differences at high-mass with Bullock  _et al._ (2001) that become smaller in
the direction of increasing $n_{s}$, which is most relevant to our
application. Given Eqs. (126-127), we can write the contribution to $\epsilon$
in analogy to the results above for an isothermal halo. We have,

$\bar{\epsilon}_{h}(m,z)={\widehat{m}^{2/3}\over
4f^{2}}\,\frac{\Omega_{m,0}^{2/3}\ \Delta_{\rm vir}^{1/3}\
R_{*}^{2}}{\Omega_{m,0}+\Omega_{\Lambda,0}\,a^{3}}\ I(\beta,c)\,,$ (128)

where we have factorized the isothermal halo contribution (see Eq. 123) and
the integral $I$ is given by

$\displaystyle I(\beta,c)$ $\displaystyle\equiv$
$\displaystyle\Big{(}6-4\beta\Big{)}{g(c)\over
c^{2}}\int_{0}^{c}dx\,x^{3-2\beta}(1+x)^{2}$
$\displaystyle\times\int_{x}^{\infty}dy{y^{2\beta-1}\over(1+y)^{2}}\Big{[}{\ln(1+y)\over
y^{2}}-{1\over y(1+y)}\Big{]}\,,$

and is a function of $\widehat{m}$ through the concentration parameter. It can
be computed analytically, see Eqs. (164-165) in App. G where we present
results for the cases of interest, namely $\beta=0$ (isotropic case) and
$\beta=1/2$ (radially biased case, i.e. radial dispersion larger than
tangential dispersion). Simulations show $\beta$ to be scale dependent, with
$\beta\approx 0$ near halo centers and $\beta$ approaching 1/2 at the halo
outskirts. Since the profile is volume weighted we consider $\beta=1/2$ to be
more realistic for our purposes.

Table 4: Velocity dispersion scale (in terms of the non-linear scale) from NFW halos as a function of spectral index. $n_{s}$ | $(k_{\sigma}/k_{\text{nl}})_{\beta=0}$ | $(k_{\sigma}/k_{\text{nl}})_{\beta=1/2}$
---|---|---
$2$ | $0.845$ | $0.863$
$1$ | $0.957$ | $0.977$
$0$ | $1.084$ | $1.105$
$-1$ | $1.219$ | $1.241$

For the case of a scale-free universe we have, integrating over the mass
function and using $d\ln\nu/d\ln m=(n_{s}+3)/3$,

$\epsilon={{n_{s}+3}\over
12}(18\pi^{2})^{1/3}R_{*}^{2}\int_{0}^{\infty}{d\widehat{m}\over\widehat{m}^{1/3}}\,[\nu
f(\nu)]_{\widehat{m}}\ I(\beta,c)\,,$ (130)

which gives us the desired $k_{\text{nl}}/k_{\sigma}$ as a function of
spectral index $n_{s}$

${k_{\sigma}\over
k_{\text{nl}}}=\Big{(}\frac{96}{\pi^{2}}\Big{)}^{1/6}\frac{k_{\text{nl}}R_{*}}{{\sqrt{(}n_{s}+3)\int
d\widehat{m}\,\widehat{m}^{-1/3}[\nu f(\nu)]_{\widehat{m}}\ I(\beta,c)}}\,.$
(131)

The mass function is measured in our simulations for different spectral
indices as discussed in App. G, where we check our mass function measurements
obey the expected self-similarity. As a summary of these measurements, Fig. 8
presents the resulting fits for $[\nu f(\nu)]$ as a function of $\widehat{m}$.
The parameters $A,a,p$ (see Eq. 120) for each spectral index are presented in
Table 6 in App. G. The behavior of the mass function at the high-mass end
agrees with standard expectations, i.e. red spectra with more large-scale
power have enhanced exponential tails. At the small-mass limit, however, blue
spectra with more initial small-scale power show a suppressed halo abundance,
highlighting the effects of increased velocity dispersion which is reflected
already by linear evolution in (see paper I Garny  _et al._ (2022)) and
confirmed by our non-linear analysis.

The results of Eq. (131) after integration over the mass function are given in
Table 4. Overall we see that the dispersion values are smaller (higher
$k_{\sigma}$) than for the corresponding isothermal halos quoted earlier, as
expected. We also see that the dispersion scale changes weakly (by about 2%)
when $\beta$ is varied within reasonable limits. As mentioned earlier (see Eq.
122) our results (as well as the mass function fits) are obtained using the
self-consistent linear solution with dispersion only. Including higher
cumulants in the linear self-consistent solution (for both the computation of
$R_{*}$ and the corresponding $m_{*}$ that appears in the mass function fits)
changes the values in Table 4 by at most $0.1\%$, which is reassuring.

It is worth noting comparing to Table 3 that these halo-based calculations
yield stronger velocity dispersion than the self-consistent one-loop
calculation done in the previous section, but by less than a factor of two.
This is remarkable, given the wide differences in the assumptions:
perturbative dynamics with non-zero dispersion everywhere on the one hand, and
dispersion only inside halos (modeled as spherical NFW profiles) on the other.
Interestingly, we shall see that fitting for $k_{\sigma}$ from the density
power spectrum measured in N-body simulations yields values that are in
between the one-loop perturbative and halo calculations. We discuss such
determinations next.

##  VI  predictions vs Simulations

In this section we present results for density and velocity power spectra as
well as the bispectrum computed within the framework of and compare to N-body
simulations. We focus on a scaling universe with linear input power spectrum
$P_{0}\propto k^{n_{s}}$ and consider the values $n_{s}=2,1,0,-1$, with the
very blue cases included to expose the strong screening of UV modes that
occurs when taking velocity dispersion and higher cumulants into account. Note
that within SPT the one-loop integrals become UV divergent for $n_{s}\geq-1$,
therefore there are no SPT predictions for these spectra beyond tree-level. As
we shall see, this shortcoming is remedied in .

Our perturbative results are based on the framework described above. In
particular, we use (cum3+, svt) as our baseline scheme, and compare the impact
of various approximation schemes (see Table 2) against this fiducial choice.
It includes apart from $\delta$ and $\theta$ the vorticity $w_{i}$, all scalar
($g,\delta\epsilon$) vector ($\nu_{i}$) and tensor ($t_{ij}$) modes of the
second cumulant as well as its background value $\epsilon(\eta)$, and higher
cumulants as detailed in Sec. II (scalar perturbations $\pi,\chi$ of the third
cumulant and background value $\omega(\eta)$ of the fourth cumulant).

###  VI.1 Density power spectrum

Figure 9:  Density power spectrum $P_{\delta\delta}(k,\eta)$ in one-loop
(solid) and linear (dashed) approximation for various values of the background
dispersion $\epsilon_{0}=1/k_{\sigma}^{2}$ for a power law input spectrum
$P_{0}\propto k^{n_{s}}$ with $n_{s}=0$ and
$\epsilon(\eta)=\epsilon_{0}e^{4\eta/(n_{s}+3)}$. The linear and one-loop
result contains the numerical kernels $F_{1,\delta}(k,\eta)$,
$F_{2,\delta}(\bm{k}-\bm{q},\bm{q},\eta)$ and
$F_{3,\delta}(\bm{k},\bm{q},-\bm{q},\eta)$ computed taking up to third
cumulant scalar perturbations as well as vorticity, vector and tensor modes of
the stress tensor into account (cum3+, svt). We set $\bar{\omega}=0.668$ (see
text and Garny  _et al._ (2022)) for the fourth cumulant expectation value
$\omega(\eta)=\bar{\omega}\times\epsilon(\eta)^{2}$.  Figure 10:  Comparison
of $P_{\delta\delta}$ in linear (dashed) and one-loop (solid) approximation in
within the (cum3+, svt) scheme to N-body results (red circles), for
$n_{s}=2,1,0,-1$, as noted. In absence of a precise knowledge of the
background dispersion, we adjusted $k_{\sigma}/k_{\text{nl}}$ by a one-
parameter fit up to $k\leq k_{\text{max}}=0.6k_{\text{nl}}$, with best-fit
values given in the legend for each case. Notably, these values of
$k_{\sigma}/k_{\text{nl}}$ are consistent with our theoretical expectations
given in Sec. V. We use the _same fixed_ values of $k_{\sigma}/k_{\text{nl}}$
as given here when comparing and N-body results for velocity spectra and
bispectra below (for each approximation and $n_{s}$, respectively) such that
for all these statistics there is only a single free parameter that is
adjusted once to a common value.

Figure 11:  Left panel: dependence of the density power spectrum
$P_{\delta\delta}(k,\eta)$ in one-loop (solid) and linear (dashed) order on
the approximation schemes (s), (sw), (sv), (svt) quantifying the impact of
successively including vorticity (sw), vector (sv) and tensor (svt) modes.
Note that (sv) and (svt) are essentially indistinguishable. Right panel:
second versus third cumulant approximations (cum2) vs (cum3+), with three
choices of $\bar{\omega}$ for the latter and a common value
$k_{\sigma}=1.31k_{\text{nl}}$ (see Fig. 10) in all cases.

We start with a discussion of the perturbative prediction of the density power
spectrum within , and then compare to N-body results. We use Eq. (II.2) to
compute the linear and one-loop contribution, with kernels
$F_{1,\delta}(k,\eta)$, $F_{2,\delta}(\bm{k}-\bm{q},\bm{q},\eta)$ and
$F_{3,\delta}(\bm{k},\bm{q},-\bm{q},\eta)$ determined numerically, as detailed
in Sec. II.4, for an input power spectrum $P_{0}\propto k^{n_{s}}$. Scaling
symmetry implies

$\epsilon(\eta)=\epsilon_{0}\,e^{4\eta/(n_{s}+3)},\ \ \ \ \
\bar{\omega}=\omega(\eta)/\epsilon(\eta)^{2}={\rm const.}$ (132)

Therefore, apart from $n_{s}$, the result depends only on the two parameters
$\epsilon_{0}$ and $\bar{\omega}$. As we shall see, the latter has only a
minor impact, and (as already mentioned in Sec. V.1) we use the value obtained
by self-consistently solving the equation of motion for the fourth cumulant in
linear approximation Garny  _et al._ (2022) as fiducial value, given by

$\bar{\omega}=0.579,\,0.616,\,0.668,\,0.752\,,$ (133)

for $n_{s}=2,1,0,-1$, respectively (we assess the dependence on $\bar{\omega}$
below). This leaves $\epsilon_{0}\equiv 1/k_{\sigma}^{2}$ as remaining
dimensionful parameter, which we quantify by the ratio
$k_{\sigma}/k_{\text{nl}}$ of the dispersion to the non-linear scale. We
discussed estimates for this quantity via two different approaches,
perturbation theory and halo modeling, in Sec. V, indicating a value for
$k_{\sigma}/k_{\text{nl}}$ in the range $\sim 1-2$. We shall see that this is
indeed consistent with direct determinations from fitting the density power
spectrum.

In Fig. 9 we show the density power spectrum in linear and one-loop
approximation obtained in based on the (cum3+, svt) scheme (see Table 2), for
three different values of $k_{\sigma}/k_{\text{nl}}$999Note that scaling
symmetry can be used to rescale numerical perturbative results for a given
value of $k_{\sigma}/k_{\text{nl}}$ to any other value, see App. E. and
$n_{s}=0$. The result is normalized to the linear power spectrum in SPT, i.e.
for the perfect pressureless fluid approximation, given by
$P_{0}(k,\eta)=e^{2\eta}P_{0}(k)$. By showing the result for
$P_{\delta\delta}(k,\eta)/P_{0}(k,\eta)$ versus $k/k_{\text{nl}}(\eta)$, the
dependence on time $\eta$ drops out, and we therefore omit the time arguments
for brevity.

From the dashed lines in Fig. 9 one can see that the linear result is
suppressed compared to the linear SPT spectrum, with the suppression being
stronger for larger background dispersion $\epsilon$, i.e. smaller
$k_{\sigma}$, as expected. Adding the one-loop correction (solid lines) leads
to additional suppression for $k\lesssim k_{\text{nl}}$, but enhances the
power spectrum relative to the linear result for $k\gtrsim k_{\text{nl}}$.
Note that the linear and one-loop results are close to each other within the
expected regime of perturbativity, i.e. for $k\lesssim{\cal
O}(k_{\text{nl}})$. This is an important property, indicating perturbative
stability. Note that within SPT the one-loop integrals are UV divergent, and
we therefore do not show any one-loop SPT result. In contrast, we checked that
the one-loop integrals are finite; this is a consequence of the strong
suppression of the non-linear kernels for large loop wavenumber as discussed
in Sec. II.4.

For comparison, we performed N-body simulations for scaling universes with
spectral indices $n_{s}=2,1,0,-1$. Our main suite of N-body simulations
consists of two sets of fixed amplitude initial conditions with opposite
phases to cancel Gaussian cosmic variance, each with $512^{3}$ particles. See
App. F for more details, including tests of self-similarity, initial
conditions, and computation of error bars.

Figure 10 shows the power spectrum $P_{\delta\delta}$ in linear (dashed) and
one-loop (solid) approximation along with the N-body measurements (red
circles). Note that error bars are too small to be noticed in
$P_{\delta\delta}$. In order to assess the level of agreement, and in absence
of a precise independent determination of the background dispersion, we
adjusted the value of $k_{\sigma}/k_{\text{nl}}$ by fitting the perturbative
to the N-body result at low wavenumber, specifically for $k\leq
k_{\text{max}}=0.6k_{\text{nl}}$. Nevertheless, we stress that as a matter of
principle this fitting procedure is not necessary to obtain predictions within
. It hinges on our ability to accurately determine the spatial average of the
stress tensor, being a well-defined quantity with immediate physical meaning.

To further test the viability of this procedure, we use the _same_ value of
$k_{\sigma}/k_{\text{nl}}$ as obtained from the density power spectrum (for a
given approximation scheme and $n_{s}$) when comparing the velocity power
spectrum and bispectrum to N-body results further below. In addition, we
checked that our results are stable under changing $k_{\text{max}}$ as long as
it is below the non-linear scale. Finally, we note that the best-fit values
for $k_{\sigma}/k_{\text{nl}}$ (as quoted in Fig. 10) lie within the expected
range as obtained in Sec. V, giving further support. Notwithstanding, the
fitting procedure likely slightly exaggerates the level of agreement, and we
comment on the impact of missing two-loop contributions below.

With these caveats in mind, we observe that the one-loop density power
spectrum can reproduce the N-body result up to $k\lesssim{\cal
O}(k_{\text{nl}})$. Furthermore, including the one-loop correction
significantly improves the agreement compared to the linear result. To be
conservative, we show the result for the optimal value of
$k_{\sigma}/k_{\text{nl}}$ for the linear approximation as well as for the sum
of the linear and one-loop contributions, respectively. This shows the actual
improvement stemming from one-loop corrections while “marginalizing” over
$k_{\sigma}/k_{\text{nl}}$. Importantly, and in line with the previous
observations, the difference between the two values of
$k_{\sigma}/k_{\text{nl}}$ (see legends in Fig. 10) as well as the difference
between linear and one-loop results for $P_{\delta\delta}$ (dashed and solid
lines in Fig. 10) up to wavenumbers $k\lesssim k_{\text{nl}}$ are reasonably
small, as required for a perturbative expansion.

####  VI.1.1 Impact of vorticity, vector and tensor modes

Figure 12:  One-loop density power spectrum $P_{\delta\delta}(k,\eta)$
compared to N-body results for $n_{s}=2,1,0,-1$ and three approximation
schemes. Taking only scalar perturbation modes into account (cum3+, s) leads
to significant deviations from our fiducial scheme (cum3+, svt) with
vorticity, vector and tensor modes. On the contrary, the second cumulant
approximation (cum2, svt) gives results close to (cum3+, svt). This means the
impact of the third cumulant is mild whereas neglecting the backreaction of
vorticity+vector(+tensor) modes would require a notably different
$k_{\sigma}/k_{\text{nl}}$ which generally makes the agreement worse for
$k\lesssim k_{\text{nl}}$ (grey lines). It is therefore important to take this
backreaction into account, i.e. use the (cum2/3+, svt) (or (cum2/3+, sv), see
Fig. 11) schemes.

Let us now turn to discussing the impact of various approximation schemes
within (see Table 2) on the density power spectrum. We start with backreaction
of vorticity, as well as vector and tensor modes of the velocity dispersion
tensor. They affect the kernel $F_{3,\delta}(\bm{k},\bm{q},-\bm{q},\eta)$
entering the one-loop density power spectrum. To assess the size of this
backreaction, we compare the density power spectrum in (s), (sw), (sv) and
(svt) approximation, all of them within the (cum3+) scheme. In the left panel
of Fig. 11 we show the comparison for a fixed common value of $k_{\sigma}$ for
all cases, and $n_{s}=0$. The linear results are identical.

We observe a significant impact of vorticity backreaction on the third-order
density contrast, i.e. a sizable shift between the one-loop density spectrum
in the (s) and (sw) approximation. For example, at the non-linear scale the
impact of vorticity backreaction is about 10% for $n_{s}=0$ (Fig. 11). We find
the size of this depends on $n_{s}$. For blue indices $n_{s}=1,2$ this grows
to about 15%, while for $n_{s}=-1$ it drops to 3%. Including the vector mode
of the dispersion tensor leads to another small shift, when going from (sw) to
(sv). In contrast to that, the tensor mode has practically no effect, with the
(sv) and (svt) lines in Fig. 11 lying on top of each other.

As discussed in Sec. III.2, including vorticity is necessary to ensure
momentum conservation. Even though the scaling
$F_{3,\delta}(\bm{k},\bm{q},-\bm{q},\eta)\propto k^{2}$ for $k\to 0$ occurs in
all schemes (and even in the momentum violating (s) scheme for the particular
set of arguments entering in the one-loop power spectrum), such that
$P_{\delta\delta}/P_{0}\to 1$ for $k\to 0$, we find that vorticity
backreaction has an important effect for $k\gtrsim 0.3k_{\text{nl}}$.

Note also that the impact of vorticity plus vector and tensor modes on the
density sets in when $q\gtrsim k$ or when $q\gtrsim k_{\sigma}$ (if
$k_{\sigma}>k$). This is the reason why in the left panel of Fig. 11 for
$k\gtrsim 3k_{\text{nl}}$ the relevance of them becomes negligible. The
contribution to the one-loop result from loop wavenumbers for which $q\gtrsim
k$ decreases when $k$ increases due to the UV screening in . Hence vorticity
backreaction eventually peters out for very high $k$ for the one-loop
contribution. Note that this is also in line with vorticity backreaction shown
in Fig. 3 and Fig. 5. Nevertheless, at these large wavenumbers, two- and
higher loops have to become relevant (see below), which should themselves be
affected by vorticity backreaction as well. Altogether, our results show that
within the mildly non-linear regime that is most relevant for perturbative
approaches, vorticity backreaction is significant.

The impact of vorticity backreaction can be seen even more clearly in Fig. 12,
by comparing the one-loop results in the (cum3+, s) and (cum3+, svt) scheme,
respectively. In this figure we adjusted $k_{\sigma}/k_{\text{nl}}$ to
optimally match the N-body result for $k\ll k_{\text{nl}}$ in each scheme. The
best-fit values as well as the power spectra differ significantly, and the
(cum3+, s) scheme shows a markedly worse agreement with N-body results
compared to (cum3+, svt) (except for $n_{s}=-1$, for which the difference
however occurs at wavenumbers where two-loop corrections are expected to
matter, see Sec. VI.1.2).

Also note that in particular for $n_{s}=0,1$ the (cum3+, s) result presents a
poor match to N-body data even at rather small $k$. The reason is that when
varying $k_{\sigma}$ the value of $P_{\delta\delta}$ at a given fixed scale
$k$ features a maximum for some particular $k_{\sigma}$, i.e.
$P_{\delta\delta}$ at one-loop is constrained to be below a certain maximal
possible value for any choice of $k_{\sigma}$. Since for (cum3+, s) and
$n_{s}=0,1$ this maximum is below the N-body result, they cannot be brought
into agreement even at low $k$ when varying $k_{\sigma}$101010The reason for
this maximum of $P_{\delta\delta}$ when varying $k_{\sigma}$ (for $k$ fixed)
is that the linear result decreases when decreasing $k_{\sigma}$
(corresponding to higher background dispersion), while the one-loop increases.
The increase of the one-loop piece can be understood from UV screening within
: the most important role is played be the contribution to the one-loop from
$F_{3,\delta}(\bm{k},\bm{q},-\bm{q},\eta)$ for $k\lesssim q$, which is
negative in SPT. The screening built into makes it less negative or even
positive (see Fig. 2). This explains the opposite behaviour of the linear and
one-loop contributions when changing $k_{\sigma}$. Their interplay then
creates the aforementioned maximum in the sum of linear and one-loop pieces..
This feature emphasizes that the agreement between and N-body results when
taking vorticity backreaction into account (as opposed to when not) is a non-
trivial result.

We conclude that it is vital to take vorticity backreaction into account, both
for physical reasons (momentum conservation, see Sec. III.2), and due to its
quantitative impact on the density contrast. In addition, vector modes of the
stress tensor have a certain effect on $P_{\delta\delta}$, while the
backreaction of tensor modes is negligible at one-loop.

Figure 13:  Velocity divergence and cross power spectra $P_{\delta\theta}$
(first row) and $P_{\theta\theta}$ (second row) for $n_{s}=2$ (left column)
and $n_{s}=0$ (right column). Red circles with error bars show the
corresponding N-body results. The linear (one-loop) approximation within the
(cum3+, svt) scheme of is shown with dashed (solid) lines, with
$k_{\sigma}/k_{\text{nl}}$ fixed to the _same_ value as for $P_{\delta\delta}$
(see Fig. 10) for each case, respectively. Therefore, these can be regarded as
unique predictions of as there is no free parameter being fit to the
measurements.

####  VI.1.2 Impact of higher cumulants and loops

To assess the relevance of cumulants beyond second order we compare the (cum2)
and (cum3+) approximation schemes, adopting (svt) in both cases. For the
latter, we also consider three choices for the fourth cumulant expectation
value, given by $\bar{\omega}=\pm 1$ in addition to our fiducial choice (see
Eq. 133). A comparison for fixed $k_{\sigma}$ and $n_{s}=0$ is shown in the
right panel of Fig. 11. The differences between the linear results for
$P_{\delta\delta}$ are hardly visible (dashed lines). At one-loop, the second
cumulant approximation (cum2) leads to a slightly smaller result compared to
(cum3+) with fiducial choice $\bar{\omega}=0.668$, while (cum3+) with
$\bar{\omega}=1$ is slightly larger. The (cum3+) case with $\bar{\omega}=-1$,
motivated by the values found from halos in paper I Garny  _et al._ (2022),
shows a smaller amplitude of the one-loop correction, close to the second
cumulant approximation. However, the differences are very small for the
relevant range $k\lesssim k_{\text{nl}}$. In addition, and in stark contrast
to vorticity as discussed above, they can largely be compensated by a small
shift in $k_{\sigma}$. This can be seen by comparing the results for (cum2,
svt) and (cum3+, svt) in Fig. 12, with both lines lying on top of each other
for $k\lesssim k_{\text{nl}}$.

We thus find no indication that higher cumulants beyond the velocity
dispersion tensor invalidate the approach for $k\lesssim{\cal
O}(k_{\text{nl}})$. On the contrary, the good agreement between (cum2) and
(cum3+) as well as small sensitivity on the fourth cumulant expectation value
$\bar{\omega}$ can be taken as an indication that the truncation of higher
cumulants leads to an acceptable uncertainty for the one-loop density power
spectrum within the mildly non-linear regime and for the considered values of
$n_{s}$. This finding is remarkable since e.g. for a single shell crossing,
higher cumulants of all orders are generated at once. Nevertheless, the
gradual build-up of an average velocity dispersion due to the superposition of
many shell crossings, as well as the question of what impact do higher
cumulants have on the large-scale density contrast complicate this simple
argument, rendering the relevance of third and higher cumulants unclear a
priori. It’s also worth pointing out that the picture that emerges from this
one-loop calculation is consistent with that of the linear approximation
detailed in paper I Garny  _et al._ (2022), i.e. while higher cumulant
perturbations are important to capture the suppression of the linear kernel
for $k\gg k_{\sigma}$, they only mildly affect the transition region between
the low-$k$ regime and the onset of suppression describing the screening of UV
modes (see Fig. 6 in Garny  _et al._ (2022)). For a more detailed discussion
of the convergence of the cumulant expansion in the linear approximation we
refer the reader to paper I Garny  _et al._ (2022).

However, while the (cum2) scheme is apparently a reasonable approximation for
$P_{\delta\delta}$ at one-loop, higher cumulants likely have a larger impact
at higher order in perturbation theory, where backreaction effects can be more
pronounced, and deserve further study in the future. Still, in order to
explore the perturbative stability we checked in an explorative study that
adding the two-loop correction for (cum2, sv) leads to only minor changes in
$k_{\sigma}/k_{\text{nl}}$ and an agreement with N-body results that is
comparable or slightly improved compared to the one-loop (cum3+, svt) case.
Interestingly, by fitting the two-loop (cum2, sv) approximation to the N-body
measurements we obtain the values:

$k_{\sigma}/k_{\text{nl}}=1.32,\,1.36,\,1.39,\,1.40\,,$ (134)

for $n_{s}=2,1,0,-1$, respectively. Note that here $k_{\sigma}/k_{\text{nl}}$
_increases_ when decreasing $n_{s}$. This is the same trend observed in Sec. V
both for the self-consistent and the halo determinations of
$k_{\sigma}/k_{\text{nl}}$. We therefore conclude that the lack of a similar
trend for the one-loop fits is most likely a limitation of the one-loop
approximation. In addition, note that the values in Eq. (134) lie in between
the independent estimates from self-consistent perturbation theory and halo
calculations (compare to Tables 3 and 4). Overall, this is compatible with the
fact that through our fitting procedure we are determining the spatial average
of the dispersion tensor, a well-defined quantity with clear physical meaning.

Figure 14:  Equilateral bispectrum at tree-level (dashed) and one-loop (solid)
compared to N-body results for $n_{s}=2,1,0,-1$, respectively. The bispectra
are normalized to the SPT tree-level bispectrum. In each case
$k_{\sigma}/k_{\text{nl}}$ is fixed to the _same_ value as for
$P_{\delta\delta}$ (see Fig. 10). Therefore, these can be regarded as unique
predictions of as there is no free parameter being fit to the measurements.

###  VI.2 Velocity divergence and cross power spectrum

Measuring the power spectrum of the velocity divergence and its cross
correlation with the density contrast from N-body simulations requires a
reconstruction of the velocity field. We describe our procedure to accomplish
this in App. H.

Within , $P_{\theta\theta}$ and $P_{\delta\theta}$ can be computed analogously
to $P_{\delta\delta}$, using Eq. (II.2) at one-loop with $a=b=\theta$ or
$a=\delta,b=\theta$, respectively. We use the _same_ values of
$k_{\sigma}/k_{\text{nl}}$ as determined in Sec. VI.1 by matching and N-body
results for $P_{\delta\delta}$, for each approximation scheme and $n_{s}$,
respectively. Therefore, there are no free parameters entering
$P_{\theta\theta}$ and $P_{\delta\theta}$, and they are predicted uniquely.
The corresponding linear and one-loop results are shown in Fig. 13 for
$n_{s}=2$ as well as $n_{s}=0$. We find that the one-loop and N-body results
are very close to each other for $k\lesssim k_{\text{nl}}$, with the agreement
being somewhat better for $n_{s}=2$. Taking the uncertainty in the velocity
reconstruction into account, the level of agreement shows that yields relevant
results also for the velocity divergence. This is an important consistency
check of the formalism. We stress again that no free parameters were adjusted
to obtain the velocity and cross power spectra, and the same is also true for
the next consistency check that we discuss, the bispectrum.

Figure 15:  Vorticity power spectrum $P_{w_{i}w_{i}}$ within in one-loop
(dashed) and two-loop (solid) approximation, for $n_{s}=2,1,0,-1$. The scaling
$P_{w_{i}w_{i}}\propto k^{2}$ occurs only starting at two-loop order, while
the one-loop result scales as $k^{4}$ for small $k$. The figure shows the
dependence of the one- and two-loop pieces on $k_{\sigma}/k_{\text{nl}}$, see
Eq. (137). The error bars correspond to the Monte Carlo integration error
(visible only at low $k$).

###  VI.3 Bispectrum

To further test the viability of , we compare perturbative and N-body results
for the bispectrum. The tree-level and one-loop bispectrum can be computed
using the formal expressions as in SPT, but replacing the SPT kernels $F_{n}$
by the numerically computed kernels $F_{n,\delta}$ within . A subtlety is that
also the first order kernel $F_{1,\delta}$ is non-trivial and therefore needs
to be taken into account, giving the expressions

$\displaystyle
B^{\text{tree}}(k_{1},k_{2},k_{3},\eta)=2e^{4\eta}F_{2,\delta}(\bm{k}_{1},\bm{k}_{2},\eta)F_{1,\delta}(k_{1},\eta)$
$\displaystyle\times$ $\displaystyle
F_{1,\delta}(k_{2},\eta)P_{0}(k_{1})P_{0}(k_{2})+2\ \text{perm.}\,,$
$\displaystyle B^{1L}(k_{1},k_{2},k_{3},\eta)=e^{6\eta}\int
d^{3}qP_{0}(q)\Big{\\{}8F_{2,\delta}(\bm{q}+\bm{k}_{1},-\bm{q},\eta)$ (135)
$\displaystyle\times$ $\displaystyle
F_{2,\delta}(\bm{q}-\bm{k}_{2},-\bm{q}-\bm{k}_{1},\eta)F_{2,\delta}(\bm{q},\bm{k}_{2}-\bm{q},\eta)P_{0}(|\bm{q}+\bm{k}_{1}|)$
$\displaystyle\times$ $\displaystyle
P_{0}(|\bm{q}-\bm{k}_{2}|)+6\big{[}F_{1,\delta}(k_{3},\eta)F_{3,\delta}(\bm{q}-\bm{k}_{2},-\bm{q},-\bm{k}_{3},\eta)$
$\displaystyle\times$ $\displaystyle
F_{2,\delta}(\bm{q},\bm{k}_{2}-\bm{q},\eta)P_{0}(|\bm{q}-\bm{k}_{2}|)P_{0}(k_{3})+F_{2,\delta}(\bm{k}_{2},\bm{k}_{3},\eta)$
$\displaystyle\times$ $\displaystyle
F_{1,\delta}(k_{2},\eta)F_{3,\delta}(\bm{k}_{3},\bm{q},-\bm{q},\eta)P_{0}(k_{2})P_{0}(k_{3})+5\
\text{perm.}\big{]}$ $\displaystyle+$ $\displaystyle
12\big{[}F_{1,\delta}(k_{2},\eta)F_{1,\delta}(k_{3},\eta)F_{4,\delta}(\bm{k}_{2},\bm{k}_{3},\bm{q},-\bm{q},\eta)$
$\displaystyle\times$ $\displaystyle P_{0}(k_{2})P_{0}(k_{3})+2\
\text{perm.}\big{]}\Big{\\}}\,.$

For the non-linear kernels entering the one-loop bispectrum, the inclusion of
vorticity and the vector mode of the dispersion tensor requires to use the
advanced algorithm described in App. D.2. We do not include the tensor mode of
the dispersion for the bispectrum for simplicity, since its impact on the
power spectrum was found to be tiny. This corresponds to the (cum3+, sv)
scheme. Furthermore, we use the algorithm described in App. B of Floerchinger
_et al._ (2019) to rewrite the integrand such that contributions that are
enhanced for $q\to 0$ cancel in the sum of all terms at the integrand level.

We normalize our results to the tree-level SPT bispectrum
$B^{\text{tree}}_{\text{SPT}}(k_{1},k_{2},k_{3},\eta)=2e^{4\eta}F_{2}(\bm{k}_{1},\bm{k}_{2})P_{0}(k_{1})P_{0}(k_{2})+2\
$perm., with EdS-SPT kernel $F_{2}(\bm{k}_{1},\bm{k}_{2})$. In the equilateral
configuration $k_{1}=k_{2}=k_{3}$, the relevant SPT kernels are given by
$F_{2}(\bm{k}_{1},\bm{k}_{2})=F_{2}(\bm{k}_{1},\bm{k}_{3})=F_{2}(\bm{k}_{2},\bm{k}_{3})=2/7$.
Note that for the considered initial power spectra with $n_{s}\geq-1$ the SPT
one-loop integrals would be UV divergent, while the result is finite, as for
the power spectrum.

In Fig. 14 we compare the tree-level and one-loop equilateral bispectrum
within to N-body results. As above, we use the _same_ value for
$k_{\sigma}/k_{\text{nl}}$ as in Sec. VI.1 for each case, such that the
bispectrum is predicted without any additional free parameters. We find that
even the tree-level bispectrum within (dashed lines) yields a rather
reasonable result, accounting for the suppression relative to the SPT tree-
level prediction (the latter corresponding to the line at $1.0$ in Fig. 14).
Adding the one-loop contribution (solid lines) further improves the agreement
with the N-body results.

Importantly, we find that despite of the large values for $n_{s}$ the one-loop
correction remains small in the perturbative regime, indicating perturbative
stability of also for the bispectrum. Indeed, the one-loop correction is
smaller for larger $n_{s}$. This can be explained by the stronger damping of
the non-linear kernels for larger $n_{s}$, which overcompensates the larger
weight of large wavenumbers due to $P_{0}(q)\propto q^{n_{s}}$. The same trend
can also be observed regarding the agreement between one-loop and N-body
results, which reaches out to rather large wavenumbers for $n_{s}=2$. We
checked that the results for the bispectrum are very similar for the (cum2)
scheme, with differences to (cum3+) comparable in size as for the power
spectrum.

Thus, we find that the formalism developed in this work yields an accurate
description not only of the density and velocity divergence power spectra, but
also of the bispectrum. This adds another piece of evidence that is able to
capture the relevant dynamics. It further demonstrates its predictivity,
since, as for the velocity divergence spectra, _no_ free parameters were
adjusted for obtaining the bispectrum.

Figure 16:  Vorticity power spectrum $P_{w_{i}w_{i}}$ within in one-loop
(dashed) and two-loop (solid) approximation, for $n_{s}=0$. The generic
scaling $P_{w_{i}w_{i}}\propto k^{2}$ for low $k$ is obtained starting at two-
loop only, since the one-loop piece is suppressed as $k^{4}$. The figure shows
the difference between the (cum2, sw) and (cum2, sv) scheme (see Table 2),
i.e. the backreaction of the vector modes of the velocity dispersion on
vorticity.

Figure 17:  Comparison of the vorticity power spectrum $P_{w_{i}w_{i}}$ in at
two-loop order (solid line) to N-body results (red circles) for $n_{s}=0$
(top) and $n_{s}=-1$ (bottom). The cross-over from $k^{4}$ to $k^{2}$ scaling
for low $k$ is clearly observed in both the N-body and two-loop results. Note
that the absolute normalization of the N-body vorticity power spectrum is
sensitive to mass resolution (see App. H). For $n_{s}=-1$ we rescaled the
N-body result by a factor $4$ in order to ease comparison of the shape. The
dashed and dot-dashed lines show the one- and two-loop contributions,
respectively. The value of $k_{\sigma}/k_{\text{nl}}$ is fixed to the one
obtained from $P_{\delta\delta}$ in the same approximation, two-loop (cum2,
sv), see Eq. (134).

###  VI.4 Vorticity power spectrum

The vorticity power spectrum $P_{w_{i}w_{i}}(k,\eta)$ is a probe that is
particularly sensitive to velocity dispersion, since no vorticity would be
generated when neglecting second and higher cumulants, as is the case in SPT.
Within cumPT, the leading contribution arises at one-loop, since vorticity
vanishes at the linear level. Apart from the overall size, the scaling with
$k$ for small wavenumber is also an interesting issue. As discussed in Sec.
IV.3, we find a scaling

$P_{w_{i}w_{i}}(k,\eta)\propto k^{2}\ \text{for}\ k\to 0\,.$ (136)

However, this scaling occurs only _starting at two-loop order_ , while the
one-loop contribution scales as $k^{4}$. Therefore, the two-loop is expected
to eventually dominate in the limit $k\to 0$, since the one-loop is
“accidentally” suppressed there. Note that, in contrast to the opposite
regime, this behaviour on large scales is definitely not a sign that
perturbation theory breaks down. It merely arises due to the strong
suppression of the one-loop piece. Starting at two loops all higher loop
corrections generically scale as $k^{2}$ for $k\to 0$.

Nevertheless, we expect that the $k^{2}$ scaling can only be observed at small
enough $k$, since also the one-loop piece contributes on intermediate scales.
To test this expectation, we compute the vorticity power spectrum up to two
loops. Given the mild differences between (cum2) and (cum3+) as well as the
tiny impact of tensor modes, we restrict ourselves to the second cumulant
approximation including vorticity as well as scalar and vector modes of the
dispersion tensor for this analysis, i.e. the (cum2, sv) scheme (see Table 2).
As for the one-loop bispectrum, computing the two-loop vorticity power
spectrum requires the algorithm described in App. D.2, where both modes of the
vorticity in the plane perpendicular to the corresponding Fourier wavevector
are included explicitly, and similarly for the vector mode.

In Fig. 15 we show the sum of the one- and two-loop contributions to the
vorticity power spectrum (solid lines), compared to the one-loop result only
(dashed lines), for $n_{s}=2,1,0,-1$. The different scaling of the one-
($\propto k^{4}$) and two-loop ($\propto k^{2}$) results for $k\ll
k_{\text{nl}}$ is clearly visible. As expected, the power spectrum is
dominated by the two-loop piece at low wavenumber, specifically for $k\lesssim
0.5k_{\text{nl}}$. The $k^{2}$ scaling is approached only gradually, occurring
for $k\lesssim 0.1-0.2k_{\text{nl}}$. The panels of Fig. 15 show the
dependence on $k_{\sigma}/k_{\text{nl}}$ for each $n_{s}$. The different
scaling with $k_{\sigma}$ of the one- and two-loop results can also be
understood analytically. The dimensionless power spectrum
$k^{3}P_{w_{i}w_{i}}$ (as well as $P_{ww}/P_{0}$ shown in Fig. 15) can only
depend on dimensionless ratios, that we can take as $k/k_{\sigma}$ and
$k_{\sigma}/k_{\text{nl}}$. Furthermore, the power of $k_{\text{nl}}$ is fixed
at a given order in perturbation theory, since the $L$-loop contribution
contains $L+1$ factors of $P_{0}=Ak^{n_{s}}$, and $A\propto
1/k_{\text{nl}}^{n_{s}+3}$. Using furthermore that $P_{w_{i}w_{i}}^{1L}\propto
k^{4}$ and $P_{w_{i}w_{i}}^{2L}\propto k^{2}$ for small $k$ implies

$\displaystyle k^{3}P_{w_{i}w_{i}}^{1L}$ $\displaystyle\propto$
$\displaystyle\left(\frac{k}{k_{\sigma}}\right)^{4+3}\left(\frac{k_{\sigma}}{k_{\text{nl}}}\right)^{2(n_{s}+3)}\propto
k_{\sigma}^{2n_{s}-1}\,,$ $\displaystyle k^{3}P_{w_{i}w_{i}}^{2L}$
$\displaystyle\propto$
$\displaystyle\left(\frac{k}{k_{\sigma}}\right)^{2+3}\left(\frac{k_{\sigma}}{k_{\text{nl}}}\right)^{3(n_{s}+3)}\propto
k_{\sigma}^{3n_{s}+4}\,,$ (137)
$\displaystyle\qquad\qquad\qquad\qquad\text{for}\ k\ll
k_{\text{nl}},k_{\sigma}\,.$

For $n_{s}=2,1$ ($n_{s}=0,-1$) this explains the increase (decrease) of
$P_{w_{i}w_{i}}^{1L}$ when increasing $k_{\sigma}$ for some fixed
$k/k_{\text{nl}}\ll 1$ (dashed lines in each panel of Fig. 15). In addition,
the second line in Eq. (137) explains the increase of $P_{w_{i}w_{i}}^{2L}$
for all $n_{s}=2,1,0,-1$, which is least pronounced for $n_{s}=-1$ (solid
lines, note that the sum of one- and two-loop is dominated by the latter for
small $k$ as discussed above).

The impact of vector modes of the stress tensor is shown in Fig. 16, where we
compare the (sw) and (sv) schemes. While both include vorticity, the vector
mode is taken into account only in the latter case. We find a very mild
difference within the perturbative regime $k\lesssim{\cal O}(k_{\text{nl}})$.

Finally, a comparison to N-body results for the vorticity power spectrum is
presented in Fig. 17. The extraction of vorticity from N-body data requires a
careful treatment, see discussion in App. H for details. In particular, as is
well-known Pueblas and Scoccimarro (2009); Hahn  _et al._ (2015); Jelic-Cizmek
_et al._ (2018), the overall amplitude of the vorticity power spectrum is
rather sensitive to the mass resolution in the simulations, but it’s _shape_
is rather robust. We verify this as well for our scale-free simulations in
App. H, showing in particular that our $n_{s}=1,2$ simulations do not have
enough resolution to capture the amplitude reliably. We therefore show results
here only for $n_{s}=0,-1$, cautioning that in particular the former case may
not be fully converged yet (see discussion in App. H). Similarly, our two-loop
calculation of the vorticity power spectrum has some uncertainty as well, as
it only includes the second cumulant. A full two-loop calculation would
include higher cumulants, including their vector modes as well. This is well
beyond the scope of the present paper. We therefore primarily base our
comparison on the $k$-dependence of the vorticity power spectrum. To aid this
comparison note that we have multiplied the N-body measurements for $n_{s}=-1$
by a factor of four, while the $n_{s}=0$ measurement is unchanged.

Given these limitations, we find a good agreement of the shape with our
perturbative result. The solid line in Fig. 17 shows the sum of one- and two-
loop contributions, while the dashed and dot-dashed lines show both of them
individually. The dip of the two-loop correction at $k\sim 1.5k_{\text{nl}}$
for $n_{s}=0$ is due to an (almost) cancellation among the contributions
$P^{2L}_{w_{1}w_{1}}$ and $P^{2L}_{w_{2}w_{2}}$ of the two vorticity degrees
of freedom (see App. D.2) at this point, while the sum is positive below and
above, as is the one-loop result. Note that in Fig. 17 we fixed
$k_{\sigma}/k_{\text{nl}}$ to the value obtained from $P_{\delta\delta}$
within the same approximation scheme, and going up to two-loop order (see Eq.
134), but our conclusion is independent of the precise value. We observe in
particular that adding the two-loop contribution significantly improves the
agreement of the shape of $P_{w_{i}w_{i}}$ with N-body results at small
wavenumbers, where the transition from $k^{4}$ to $k^{2}$ scaling occurs, as
captured by the two-loop part. It is remarkable that this transition is seen
both in two-loop and N-body results.

Thus, even though the overall normalization of the vorticity power spectrum is
subject of further study, the _shape_ of the perturbative and N-body results
agrees very well in the weakly non-linear regime. This constitutes yet another
non-trivial consistency check of our approach, and demonstrates the
versatility and predictivity of .

The low-$k$ behavior of the vorticity power spectrum has been the subject of
some study in the literature. In the context of $\Lambda$CDM simulations at
$z=0$, Jelic-Cizmek  _et al._ (2018) measure a low-$k$ asymptotic behavior
$P_{ww}\propto k^{n_{w}}$ with $n_{w}=2.55\pm 0.02$ (for
$k\mathrel{\mathchoice{\lower
3.6pt\vbox{\halign{$\mathsurround=0pt\displaystyle\hfill#\hfil$\cr<\crcr\sim\crcr}}}{\lower
3.6pt\vbox{\halign{$\mathsurround=0pt\textstyle\hfill#\hfil$\cr<\crcr\sim\crcr}}}{\lower
3.6pt\vbox{\halign{$\mathsurround=0pt\scriptstyle\hfill#\hfil$\cr<\crcr\sim\crcr}}}{\lower
3.6pt\vbox{\halign{$\mathsurround=0pt\scriptscriptstyle\hfill#\hfil$\cr<\crcr\sim\crcr}}}}0.4h/\mathrm{Mpc}$)
in agreement with Hahn  _et al._ (2015) which quote $n_{w}\simeq 5/2$ for
$k\simeq 0.1h/\mathrm{Mpc}$. Note that these works use small simulation boxes,
$256\mathrm{Mpc}/h$ and $100-300\mathrm{Mpc}/h$, respectively, thus most
likely they have not reached the regime where predicts a transition to $k^{2}$
scaling, i.e. $k\simeq(0.1-0.2)\,k_{\rm nl}$. On the other hand, from the
perturbative side, there have been two predictions based on EFT and one based
on Lagrangian perturbation theory (LPT). From EFT, Carrasco  _et al._ (2014)
predict a “leading order” contribution at the one-loop level with
$n_{w}=7+3n_{s}$, and a stochastic component with $n_{w}=4$. On the other
hand, Mercolli and Pajer (2014) predict $n_{w}=7+2n_{s}$ and $n_{w}=4$. From
LPT, Cusin  _et al._ (2017) also predict $n_{w}=4$. For our simulations with
$n_{s}\geq-1$, all these three predictions would lead to a low-$k$ slope of
$n_{w}=4$. As discussed above, the collisionless Vlasov dynamics captured by
yields $n_{w}=2$, in agreement with our N-body results. This result is
independent of the input power spectrum, and thus also holds for $\Lambda$CDM
power spectra.

##  VII Conclusions

In this paper we presented non-linear perturbative solutions to the Vlasov
hierarchy, a cumulant expansion of the Vlasov-Poisson equations of motion
describing collisionless matter. This Vlasov perturbation theory () represents
a systematic and predictive extension of the well-known framework of standard
perturbation theory (SPT).

Formally, this approach is similar to SPT, in that the calculation of density
and velocity divergence fields has the same expressions but with different
non-linear kernels $F_{n}/G_{n}\to F_{n,a}$ and an enlarged set of
perturbation modes $a=\delta,\theta,w_{i},...$. One key modification is that
the linear kernel ($n=1$) is not unity, except in the $k\to 0$ limit, with a
scale-dependence describing how linear growth is suppressed at small scales
once $k$ crosses the (time-dependent) dispersion scale $k_{\sigma}$ set by the
expectation value of the velocity dispersion tensor. In contrast with SPT,
there are additional degrees of freedom corresponding to vorticity $w_{i}$,
the velocity dispersion tensor $\epsilon_{ij}$ and higher cumulants ${\cal
C}_{ijk\dots}$ of the DF, which provide an important new window to testing the
theory against simulations.

For the non-linear kernels appearing in loop corrections, the suppression
compared to the equivalent SPT kernels is present even when their total
momentum $k<k_{\sigma}$ as long as the loop momenta $q_{i}$ cross the
dispersion scale. Indeed, when such crossing takes place the linear modes
cease to grow (as described by the linear kernel) and therefore they cease to
source the non-linear kernels. This reduces the UV sensitivity of loops
compared to SPT, capturing the expected decoupling of UV modes. This allows us
to compute non-linear corrections to power spectra even for cosmologies with
very blue power-law input spectra, for which SPT does not converge.
Furthermore, it realizes the screening of UV modes for any input spectrum,
considerably reducing the sensitivity of non-linear corrections to highly non-
linear scales.

Along the way, we obtained a number of noteworthy results, namely:

  * i)

We provide analytical results when expanding in powers of the average
background value of the velocity dispersion tensor (Sec. II.3). These
solutions illustrate the type of corrections to SPT, but are not sufficient to
capture the UV screening.

  * ii)

We characterize the impact of higher cumulants of the DF on the non-linear
kernels, and in addition we study the impact of scalar, vector and tensor
modes of the second cumulant, i.e. the velocity dispersion tensor (Sec. II.4).

  * iii)

We demonstrate screening of UV modes, finding a universality that can be
explained analytically. We derive an asymptotic scaling relation Eq. (II.4)
showing how non-linear density and velocity divergence kernels are suppressed
relative to their SPT counterparts in the limit $k=|\sum_{i}\bm{k}_{i}|\ll
k_{\sigma}\ll|\bm{k}_{i}|$,

$\displaystyle\ \ \ \ \
F_{n,\delta/\theta}(\bm{k}_{1},\dots,\bm{k}_{n},\eta)\sim\frac{k^{2}}{q^{2}}\,\left(\frac{k_{\sigma}}{k_{1}}\right)^{2/\alpha}\\!\\!\\!\cdots\left(\frac{k_{\sigma}}{k_{n}}\right)^{2/\alpha}\\!\\!\\!\,.$

This makes the UV screening manifest for large internal wavenumbers $k_{i}$,
with an extra suppression on top of the $k^{2}/q^{2}$ scaling imposed by
momentum conservation well-known from SPT kernels, $F_{n}/G_{n}\sim
k^{2}/q^{2}$ (where $q\sim\text{max}_{i}|\bm{k}_{i}|$). The exponent is
related to the growth rate of the background dispersion
$\alpha\equiv\partial_{\eta}\ln\epsilon$, and given by $2/\alpha=(n_{s}+3)/2$
for a scaling universe with spectral index $n_{s}$. This result is universal
(e.g. valid in all truncation schemes we considered) and only hinges on the
stalling of usual linear growth once a perturbation mode enters the
“dispersion horizon”, i.e. for
$k_{i}>k_{\sigma}(\eta)=1/\epsilon(\eta)^{1/2}$. This scaling implies that the
one-loop integral Eq. (II.2) is UV finite for any spectral index $n_{s}$ and
explains Eq. (3).

  * iv)

We derive constraints on the non-linear kernels from underlying symmetries,
such as Galilean invariance and momentum conservation (Sec. III). These
generalize well-known relations to the case when including dispersion and
higher cumulants, and provide non-trivial checks and implications for viable
approximation schemes that truncate the hierarchy. In particular, we show that
neglecting vorticity when including dispersion leads to a violation of
momentum conservation. However, when including vorticity the truncations we
explore do respect momentum conservation and Galilean invariance.

  * v)

We discuss the generation of vorticity as well as vector and tensor modes of
the velocity dispersion at second and higher order in perturbation theory
(Sec. IV). We find that the vorticity power spectrum generated non-linearly
scales as $k^{2}$ for low wavenumber in general. However, the leading (one-
loop) contribution is accidentally suppressed as $k^{4}$, such that the
$k^{2}$ scaling appears only starting at two-loop order.

  * vi)

We provide estimates of the dispersion scale within two complementary
frameworks. First, we use a purely perturbative approach based on self-
consistent solutions of the coupled equations for the perturbations and the
average dispersion (Sec. V.1). Second, we use the halo mass function measured
in simulations together with the calculations of the velocity dispersion
profile in NFW halos to compute the dispersion scale from virialized regions
(Sec. V.2). Remarkably, both of these approaches yield estimates that are
close to each other, while not allowing for a precise determination.

  * vii)

The linear theory gives for the first time a way to calculate the non-linear
halo mass for blue spectra (Sec. V.2), avoiding UV divergences that appear for
spectral indices $n_{s}\geq 1$. This is relevant as well for CDM spectra when
considering biased tracers, to avoid UV divergences that appear when dealing
with derivatives of the smoothed density field in some bias schemes, see e.g.
Paranjape  _et al._ (2013).

  * viii)

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$ (Sec. VI). We use the
density power spectrum to fix the precise value of the background dispersion,
which then completely determines all other power- and bispectra. 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}$. Even though, naively (e.g. within
SPT), the loop integration would be UV dominated for large $n_{s}$, we find
that this is in fact not the case when including higher cumulants due to the
screening of small-scale modes that becomes more efficient with larger
$n_{s}$. As a result of this, reconciles perturbation theory behavior with the
observation that blue spectra in simulations show most suppression of non-
linear power, explaining the dependence of non-linear power on linear power
with spectral index.

  * ix)

Interestingly, we find the values for the dispersion scale fixed from the two-
loop density power spectrum (Eq. 134) is in between the perturbative estimate
and the halo estimate (compare to Tables 3 and 4, respectively). Furthermore,
the trend of these values with spectral index is also consistent. This is
compatible with the fact that by matching the density power spectrum we are
determining the spatial average of the dispersion tensor, a well-defined
quantity with clear physical meaning.

  * x)

We observe a significant impact of vorticity backreaction on the (third- and
higher order) density contrast. At the non-linear scale the impact of
vorticity backreaction on the density power spectrum is $\sim 10\%$ for
$n_{s}=0$ (Fig. 11). For blue indices $n_{s}=1,2$ this grows to about 15%,
while for $n_{s}=-1$ it drops to 3%. The backreaction of the vector mode of
the dispersion tensor also has a certain impact, while its tensor modes yield
a negligible backreaction up to one-loop.

  * xi)

A nontrivial prediction of is the vorticity power spectrum, which we compute
up to two-loop order in order to recover the correct large-scale limit.
Comparing to our N-body measurements confirms the cross-over from $k^{4}$ to
$k^{2}$ scaling on large scales (Sec. VI.2).

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. In particular, this deterministic approach features a
screening of UV modes and thereby abandons one of the main shortcomings of the
ubiquitous SPT approximation. This motivates to develop the framework further
and study also its application in the context of $\Lambda$CDM cosmologies.
Another obvious extension is to carry out a systematic study of two-loop
corrections, and check the convergence of , which given the UV screening may
be expected to hold best for blue spectra, the opposite of what happens for
SPT. Finally, along the same lines, we plan to study power spectra response
functions, i.e. the dependence of the non-linear power spectrum on the initial
spectrum. The suppression of this reponse seen in N-body simulations
Nishimichi  _et al._ (2016, 2017) compared to SPT is another manifestation of
the screening of UV modes (see Fig. 1), and a quantitative comparison of this
provides yet another test of .

###### Acknowledgements.

We thank Michael Buehlmann, ChangHoon Hahn, Oliver Hahn and Lam Hui for useful
discussions, and E. Sefusatti for useful exchange on fixed-amplitude initial
conditions and sharing his unpublished results with A. Oddo Oddo and Sefusatti
(2022). We acknowledge support by the Excellence Cluster ORIGINS, which is
funded by the Deutsche Forschungsgemeinschaft (DFG, German Research
Foundation) under Germany’s Excellence Strategy - EXC-2094 - 390783311.

##  APPENDIX A Tadpoles

In this appendix we discuss subtleties related to the splitting between
background value and perturbation. For example, for the velocity dispersion
the equation of motion for the perturbations
$\delta\epsilon_{ij}=\epsilon_{ij}-\langle\epsilon_{ij}\rangle$ are obtained
by subtracting the Eq. (13) for the average value
$\langle\epsilon_{ij}\rangle=\epsilon(\eta)\delta_{ij}$ from the evolution
equation for $\epsilon_{ij}=\sigma_{ij}/(f{\cal H})^{2}$ obtained from the
second moment of the Vlasov equation Eq. (II.1). The same holds for the
perturbation of $A=\ln(1+\delta)$ around its average ${\cal A}=\langle
A\rangle$, that we denote by $\delta A\equiv A-\langle A\rangle$ in this
section for clarity (in the rest of the work we use the symbol $A$ also for
the perturbation part in a slight abuse of notation). The equation of motion
for its background value reads ${\cal A}^{\prime}=Q_{A}(\eta)$ with source
term $Q_{A}(\eta)\equiv-\int d^{3}kP_{\theta A}(k,\eta)$ Garny  _et al._
(2022).

The evolution equation for the vector Eq. (21) of perturbation modes obtained
in this way reads

$\psi_{k,a}^{\prime}(\eta)+\Omega_{ab}(k,\eta)\psi_{k,b}(\eta)=\int_{pq}\hat{\gamma}_{abc}(\bm{p},\bm{q})\psi_{p,b}(\eta)\psi_{q,c}(\eta)-Q_{a}(\eta)\delta^{(3)}(\bm{k})\,,$
(139)

which differs from Eq. (23) given in the main text by the last term, and
contains vertices $\hat{\gamma}_{abc}(\bm{p},\bm{q})$, that, within this
section, denote the vertices as derived from the Vlasov equation (we use the
hatted notation to discriminate them from modified vertices to be defined
below). The last term arises from subtracting the equations for the average
values from the original evolution equations for $\delta\epsilon$ and $\delta
A$, respectively. We use a vector notation with $Q_{a}\equiv Q$ (see Eq. (20))
for $a=\delta\epsilon$, $Q_{a}\equiv Q_{A}$ for $a=\delta A$, and $Q_{a}=0$
for all other perturbation modes considered here, that do not possess average
values. The last term in Eq. (139) ensures that
$\langle\psi_{k,a}(\eta)\rangle=0$, i.e. that $\delta\epsilon$ and $\delta A$
indeed correspond to the perturbation around the average.

In this appendix we show that the term $-Q_{a}(\eta)\delta^{(3)}(\bm{k})$ can
be skipped in the evolution equations provided one replaces the vertices
$\hat{\gamma}_{abc}(\bm{p},\bm{q})$ as derived from the Vlasov equation by
vertices $\gamma_{abc}(\bm{p},\bm{q})$ given by

$\gamma_{abc}(\bm{p},\bm{q})\equiv\left\\{\begin{array}[]{ll}\hat{\gamma}_{abc}(\bm{p},\bm{q})&\bm{p}+\bm{q}\not=0\,,\\\
0&\bm{p}+\bm{q}=0\,.\end{array}\right.$ (140)

Note that, apart from this section, we do not discriminate between
$\gamma_{abc}$ and $\hat{\gamma}_{abc}$, but it is implicitly understood that
vertices are set to zero for $\bm{p}+\bm{q}=0$.

The statement from above can be most easily understood using a diagrammatic
representation of the perturbative expansion, see e.g. Bernardeau  _et al._
(2002). We note that it can be extended in a straightforward way to the
generalized case of a vector of perturbations considered here. Starting point
is the integral representation of the equation of motion

$\psi_{k,a}(\eta)=\psi^{\text{lin}}_{k,a}(\eta)+\int^{\eta}d\eta^{\prime}g_{aa^{\prime}}(\eta,\eta^{\prime},k)\left[\int_{pq}\hat{\gamma}_{a^{\prime}bc}(\bm{p},\bm{q})\psi_{p,b}(\eta^{\prime})\psi_{q,c}(\eta^{\prime})-Q_{a^{\prime}}(\eta)\delta^{(3)}(\bm{k})\right]\,,$
(141)

where $g_{aa^{\prime}}(\eta,\eta^{\prime},k)$ is the retarded propagator,
defined by the solution of the linear equation
$\partial_{\eta}g_{aa^{\prime}}(\eta,\eta^{\prime},k)+\Omega_{ab}(k,\eta)g_{ba^{\prime}}(\eta,\eta^{\prime},k)=0$
with boundary condition
$g_{aa^{\prime}}(\eta^{\prime},\eta^{\prime},k)=\delta_{aa^{\prime}}$ for
$\eta=\eta^{\prime}$, and $g_{aa^{\prime}}(\eta,\eta^{\prime},k)=0$ for
$\eta<\eta^{\prime}$. The perturbative expansion Eq. (39) can be generated by
iteratively inserting the linear solution
$\psi^{\text{lin}}_{k,a}(\eta)=e^{\eta}F_{1,a}(k,\eta)\delta_{k0}$ inside the
non-linear terms in Eq. (141), und using a lower integration boundary formally
shifted to $-\infty$.

The propagator and vertices lead to a diagrammatic expansion that is formally
analogous to the SPT case Bernardeau  _et al._ (2002), except that the
propagator depends on both time arguments as well as on wavenumber, has more
components, and more vertices contribute. In addition, the source terms
$Q_{a}$ can formally be represented by insertions given by
$-Q_{a^{\prime}}(\eta)\delta^{(3)}(\bm{k})$ from which a single propagator
line with wavenumber $\bm{k}=0$ emerges. An example will be shown below. As we
will see, their role is to precisely cancel “tadpoles diagrams”. Similar
cancellations are well-known in the context of effective actions considered in
quantum field theory.

Tadpole diagrams contain some loops that are connected only via a _single_
line to the rest of the diagram. By momentum conservation, this line is
restricted to carry vanishing wavenumber. Such diagrams can arise for
perturbations that also possess a homogeneous background value, i.e. $\delta
A$ and $\delta\epsilon$ in the present case.

An example for a tadpole contribution to the one-loop power spectrum is shown
in Fig. 18 (left). The tadpole loop can be non-zero if the vertex
$\gamma_{abc}(p,q)$ belonging to the loop has a non-vanishing limit for
$\bm{p}+\bm{q}\to 0$. In the following we restrict the discussion to scalar
perturbations up to the second cumulant for simplicity, but note that the
arguments below can be extended to the general case. Within this restriction,
the vertices with a non-zero limit for $\bm{p}+\bm{q}\to 0$ are

$\gamma_{A\theta A}(p,q)\to-\frac{1}{2}\,,\gamma_{g\theta
g}(p,q)\to\frac{1}{4}(3c_{p}^{2}-1)\,,\gamma_{g\theta\epsilon}(p,q)\to\frac{1}{2}(3c_{p}^{2}-1)\,,\gamma_{\epsilon\theta
g}(p,q)\to\frac{1}{4}s_{p}^{2}\,,\gamma_{\epsilon\theta\epsilon}(p,q)\to-\frac{1}{2}c_{p}^{2}\,.$
(142)

Here $c_{p}(s_{p})$ denotes the cosine (sine) of the angle between between
$\bm{p}$ and $\bm{k}\equiv\bm{p}+\bm{q}$. The tadpole loop is given by

$T_{a}(\bm{k},\eta)\equiv\delta^{(3)}(\bm{k})\int
d^{3}p\,\gamma_{abc}(\bm{p},\bm{q})\Big{|}_{\bm{p}+\bm{q}\to
0}\,P_{bc}(p,\eta)\,,$ (143)

where $P_{bc}(p,\eta)$ denotes the linear power spectrum at time $\eta$ (see
below for a generalization to the non-linear case). Since $P_{bc}(p,\eta)$
does not depend on the direction of $\bm{p}$, we can consider the angle-
averaged vertices

$\int\frac{d\Omega}{4\pi}\gamma_{A\theta
A}(\bm{p},\bm{q})\to-\frac{1}{2}\,,\int\frac{d\Omega}{4\pi}\gamma_{g\theta
g}(\bm{p},\bm{q})\to
0\,,\int\frac{d\Omega}{4\pi}\gamma_{g\theta\epsilon}(\bm{p},\bm{q})\to
0\,,\int\frac{d\Omega}{4\pi}\gamma_{\epsilon\theta
g}(\bm{p},\bm{q})\to\frac{1}{6}\,,\int\frac{d\Omega}{4\pi}\gamma_{\epsilon\theta\epsilon}(\bm{p},\bm{q})\to-\frac{1}{6}\,.$
(144)

The contribution from tadpoles containing $\gamma_{g\theta g}$ and
$\gamma_{g\theta\epsilon}$ averages to zero, consistent with the fact that $g$
does not have a homogeneous background value.

Figure 18:  Tadpole diagram contributing to the one-loop power spectrum
(left), and a corresponding diagram with an insertion of the source term
$Q(\eta)\delta^{(3)}(\bm{k})$ contributing to the equation of motion for the
perturbations (right). Both diagrams cancel each other, to ensure that the
perturbation around the background maintains a vanishing mean value. It is
equivalent to skip the right diagram, and use the prescription to set the
vertices to zero when the ingoing lines _exactly_ add up to zero, thereby
eliminating also the left diagram.

Let us consider as an example the tadpole $T_{a}$ for $a=\delta\epsilon$.
Using the angle-averaged vertices, the tadpole Eq. (143) is given by

$T_{\delta\epsilon}(\bm{k},\eta)=\delta^{(3)}(\bm{k})\,\frac{1}{3}\int
d^{3}p\,(P_{\theta g}(p,\eta)-P_{\theta\,\delta\epsilon}(p,\eta))\,.$ (145)

The tadpole is _exactly_ cancelled by the contribution
$-Q(\eta)\delta^{(3)}(\bm{k})$ from the source term to the perturbation
equation for $\delta\epsilon$ Eq. (20) (considering the restriction to scalar
modes up to the second cumulant for the purpose of this discussion), as
expected. The corresponding contribution to the power spectrum is shown in the
right part of Fig. 18.

Note that this cancellation also works at higher loop order. We can decompose
any loop diagram in tadpole contributions, that are attached to some internal
or external lines of the full diagram. The tadpole itself can also contain
loops, that precisely yield the full non-linear power spectrum
$P_{bc}(p,\eta)$ in the integrand in Eq. (143). On the other hand, the source
term $Q$ also contains the integral over the corresponding full non-linear
power spectrum, such that the cancellation occurs at any order in perturbation
theory. The same argument applies to the tadpole Eq. (143) for $a=\delta A$.
It can also be generalized to perturbations of the fourth, sixth, eighth, or
any even cumulant, that may possess an average value (see paper I Garny  _et
al._ (2022)), but are not considered in this work.

One may wonder whether the vertex connecting the line coming from the tadpole
to the rest of the diagram can potentially be singular for vanishing ingoing
wavenumber. One can check that either this singularity is removed by a factor
of $k^{2}$ coming from the internal propagating line connecting with the
tadpole, or that the vertex is not singular.

Therefore, we conclude that with the prescription of setting vertices to zero
when the ingoing lines _exactly_ add up to zero effectively eliminates all
tadpoles, such that we can also ignore diagrams with insertions of the sources
$Q$ or $Q_{A}$. This is equivalent to using the equation of motion Eq. (23)
instead of Eq. (139), with vertices Eq. (140), as claimed above.

##  APPENDIX B Time integrals

The integrals appearing in the analytic results Eq. (II.3.2) for the non-
linear kernels $F_{2,\delta}$ and $F_{2,\theta}$ when expanding to first order
in the power of the background dispersion $\epsilon(\eta)$ are given by

$\displaystyle J^{\delta}_{1}$ $\displaystyle\equiv$
$\displaystyle\frac{4}{5}\int^{\eta}d\eta^{\prime}\,\left(e^{\eta^{\prime}-\eta}-e^{7(\eta^{\prime}-\eta)/2}\right)E_{3}(\eta^{\prime})=\frac{2}{5}(E_{1}(\eta)-5E_{3}(\eta)+4E_{7/2}(\eta))\,,$
$\displaystyle J^{\delta}_{2}$ $\displaystyle\equiv$
$\displaystyle\frac{4}{5}\int^{\eta}d\eta^{\prime}\,\left(e^{\eta^{\prime}-\eta}-e^{7(\eta^{\prime}-\eta)/2}\right)(E_{2}(\eta^{\prime})-E_{3}(\eta^{\prime}))=\frac{2}{15}(3E_{1}(\eta)-10E_{2}(\eta)+15E_{3}(\eta)-8E_{7/2}(\eta))\,,$
$\displaystyle J^{\delta}_{3}$ $\displaystyle\equiv$
$\displaystyle\frac{2}{5}\int^{\eta}d\eta^{\prime}\,\left(e^{\eta^{\prime}-\eta}-e^{7(\eta^{\prime}-\eta)/2}\right)\epsilon(\eta^{\prime})=\frac{2}{5}(E_{1}(\eta)-E_{7/2}(\eta))\,,$
$\displaystyle J^{\delta}_{4}$ $\displaystyle\equiv$
$\displaystyle\frac{2}{5}\int^{\eta}d\eta^{\prime}\,\left(e^{\eta^{\prime}-\eta}-e^{7(\eta^{\prime}-\eta)/2}\right)I_{\theta}(\eta^{\prime})=\frac{2}{105}(12E_{0}(\eta)-70E_{2}(\eta)+63E_{5/2}(\eta)-5E_{7/2}(\eta))\,,$
$\displaystyle J^{\delta}_{5}$ $\displaystyle\equiv$ $\displaystyle
J^{\delta}_{1}+J^{\delta}_{2}=\frac{4}{15}(3E_{1}(\eta)-5E_{2}(\eta)+2E_{7/2}(\eta))\,,$
$\displaystyle J^{\delta}_{6}$ $\displaystyle\equiv$
$\displaystyle\frac{2}{5}\int^{\eta}d\eta^{\prime}\,\left(\frac{3}{2}e^{\eta^{\prime}-\eta}+e^{7(\eta^{\prime}-\eta)/2}\right)I_{\theta}(\eta^{\prime})=\frac{2}{21}(6E_{0}(\eta)-7E_{2}(\eta)+E_{7/2}(\eta))\,,$
$\displaystyle J^{\delta}_{7}$ $\displaystyle\equiv$
$\displaystyle\frac{2}{5}\int^{\eta}d\eta^{\prime}\,\left(\frac{3}{2}e^{\eta^{\prime}-\eta}+e^{7(\eta^{\prime}-\eta)/2}\right)I_{\delta}(\eta^{\prime})=\frac{2}{105}(30E_{0}(\eta)-63E_{1}(\eta)+35E_{2}(\eta)-2E_{7/2}(\eta))\,,$
(146)

and

$\displaystyle J^{\theta}_{1}$ $\displaystyle\equiv$
$\displaystyle\frac{4}{5}\int^{\eta}d\eta^{\prime}\,\left(e^{\eta^{\prime}-\eta}+\frac{3}{2}e^{7(\eta^{\prime}-\eta)/2}\right)E_{3}(\eta^{\prime})=\frac{2}{5}(E_{1}(\eta)+5E_{3}(\eta)-6E_{7/2}(\eta))\,,$
$\displaystyle J^{\theta}_{2}$ $\displaystyle\equiv$
$\displaystyle\frac{4}{5}\int^{\eta}d\eta^{\prime}\,\left(e^{\eta^{\prime}-\eta}+\frac{3}{2}e^{7(\eta^{\prime}-\eta)/2}\right)(E_{2}(\eta^{\prime})-E_{3}(\eta^{\prime}))=\frac{2}{5}(E_{1}(\eta)-5E_{3}(\eta)+4E_{7/2}(\eta))\,,$
$\displaystyle J^{\theta}_{3}$ $\displaystyle\equiv$
$\displaystyle\frac{2}{5}\int^{\eta}d\eta^{\prime}\,\left(e^{\eta^{\prime}-\eta}+\frac{3}{2}e^{7(\eta^{\prime}-\eta)/2}\right)\epsilon(\eta^{\prime})=\frac{1}{5}(2E_{1}(\eta)+3E_{7/2}(\eta))\,,$
$\displaystyle J^{\theta}_{4}$ $\displaystyle\equiv$
$\displaystyle\frac{2}{5}\int^{\eta}d\eta^{\prime}\,\left(e^{\eta^{\prime}-\eta}+\frac{3}{2}e^{7(\eta^{\prime}-\eta)/2}\right)I_{\theta}(\eta^{\prime})=\frac{1}{35}(16E_{0}(\eta)-21E_{5/2}(\eta)+5E_{7/2}(\eta))\,,$
$\displaystyle J^{\theta}_{5}$ $\displaystyle\equiv$ $\displaystyle
J^{\theta}_{1}+J^{\theta}_{2}=\frac{4}{5}(E_{1}(\eta)-E_{7/2}(\eta))\,,$
$\displaystyle J^{\theta}_{6}$ $\displaystyle\equiv$
$\displaystyle\frac{2}{5}\int^{\eta}d\eta^{\prime}\,\left(\frac{3}{2}e^{\eta^{\prime}-\eta}-\frac{3}{2}e^{7(\eta^{\prime}-\eta)/2}\right)I_{\theta}(\eta^{\prime})=\frac{1}{35}(12E_{0}(\eta)-70E_{2}(\eta)+63E_{5/2}(\eta)-5E_{7/2}(\eta))\,,$
$\displaystyle J^{\theta}_{7}$ $\displaystyle\equiv$
$\displaystyle\frac{2}{5}\int^{\eta}d\eta^{\prime}\,\left(\frac{3}{2}e^{\eta^{\prime}-\eta}-\frac{3}{2}e^{7(\eta^{\prime}-\eta)/2}\right)I_{\delta}(\eta^{\prime})=\frac{2}{35}(6E_{0}(\eta)-21E_{1}(\eta)+35E_{2}(\eta)-21E_{5/2}(\eta)+E_{7/2}(\eta))\,.$
(147)

##  APPENDIX C Recursion relation and Galilean invariance

In this appendix we prove the relation Eq. (79), by generalizing the proof
within EdS-SPT presented in Sugiyama and Futamase (2013). The crucial
observation is that the non-linear vertices $\gamma_{abc}(\bm{p},\bm{q})$
satisfy the property Eq. (78) in the limit $\bm{p}\to 0$.

We proceed by induction and assume that Eq. (79) is satisfied for all $m\leq
n$. For the moment we also assume that it is satisfied for the starting point
$n=1$, and show that this is indeed the case further below. Now we want to
proove Eq. (79) for $F_{n+1,a}$. The differential equation Eq. (II.2) for the
kernel can be written as

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

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}\\}$. The symmetrization with respect to the
vector $\bm{p}$ is _not_ included in $\\{\cdots\\}^{s}$ in the equation above,
but explicitly accounted for by the two sums over $m$ in the second and third
line and the two contributions in the last line, respectively, as well as
their prefactors.

In this form, the limit $\bm{p}\to 0$ can be performed, using Eq. (79) for the
kernels $F_{m+1,b}(\bm{p},\bm{q}_{1},\dots,\bm{q}_{m},\eta)$ and
$F_{n-m+1,c}(\bm{p},\bm{q}_{m+1},\dots,\bm{q}_{n},\eta)$, respectively, as
well as

$\displaystyle\gamma_{abc}(\bm{p}+\bm{q}_{1}+\cdots+\bm{q}_{m},\bm{q}_{m+1}+\cdots+\bm{q}_{n})$
$\displaystyle\to$
$\displaystyle\gamma_{abc}(\bm{q}_{1}+\cdots+\bm{q}_{m},\bm{q}_{m+1}+\cdots+\bm{q}_{n})\,,$
$\displaystyle\gamma_{abc}(\bm{q}_{1}+\cdots+\bm{q}_{m},\bm{p}+\bm{q}_{m+1}+\cdots+\bm{q}_{n})$
$\displaystyle\to$
$\displaystyle\gamma_{abc}(\bm{q}_{1}+\cdots+\bm{q}_{m},\bm{q}_{m+1}+\cdots+\bm{q}_{n})\,.$
(149)

For the first term in the last line of Eq. (C) we use the property Eq. (78) of
the vertices, as well as $F_{1,\theta}(\bm{p},\eta)\to 1$ for $\bm{p}\to 0$
and $F_{1,w_{i}}(\bm{p},\eta)=0$. At this point we have to assume that
$|\bm{p}|\ll k_{\sigma}$ in addition to $|\bm{p}|\ll|\bm{k}_{i}|$. Using that
the vertices are symmetrized, see Eq. (38), the second term in the last line
of Eq. (C) is identical to the first one. Finally, we use that
$\Omega_{ab}(|\bm{k}+\bm{p}|,\eta))\to\Omega_{ab}(\bm{k},\eta))$. Altogether,
we obtain

$\displaystyle(\partial_{\eta}\delta_{ab}+(n+1)\delta_{ab}+\Omega_{ab}(k,\eta))F_{n+1,b}(\bm{p},\bm{k}_{1},\dots,\bm{k}_{n},\eta)$
(150) $\displaystyle=$
$\displaystyle\frac{1}{n+1}\frac{\bm{k}\cdot\bm{p}}{p^{2}}\sum_{m=1}^{n-1}\Big{\\{}\gamma_{abc}(\bm{q}_{1}+\cdots+\bm{q}_{m},\bm{q}_{m+1}+\cdots+\bm{q}_{n})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}$
$\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)\
+{\cal O}(p^{0})\,.$

The sum in the middle line is precisely the right-hand side of the
differential equation Eq. (79) for the kernel
$F_{n,a}(\bm{k}_{1},\dots,\bm{k}_{n},\eta)$. Therefore, taking also the last
line into account, we can also write this equation as

$(\partial_{\eta}\delta_{ab}+(n+1)\delta_{ab}+\Omega_{ab}(k,\eta))\left[F_{n+1,b}(\bm{p},\bm{k}_{1},\dots,\bm{k}_{n},\eta)-\frac{1}{n+1}\frac{\bm{k}\cdot\bm{p}}{p^{2}}F_{n,b}(\bm{k}_{1},\dots,\bm{k}_{n},\eta)\right]={\cal
O}(p^{0})\,.$ (151)

For $a=\delta,\theta$, the kernels $F_{n,a}$ approach the EdS-SPT kernels at
early times, which are known to satisfy Eq. (79) Sugiyama and Futamase (2013).
For all other perturbation modes, the kernels can be assumed to vanish at
early times, and are generated only by time evolution. Therefore, the square
bracket can be taken to be zero in the limit $\eta\to-\infty$ for all kernels,
up to corrections of order $p^{0}$. Consequently, the unique solution of the
differential equation is that the square bracket is zero (up to $p^{0}$ terms)
at _all_ times. This completes the proof for all $n\geq 2$. What is still
missing is to show Eq. (79) for the starting point of the induction, $n=1$.
This can be achieved by realizing that when proceeding analogously as above
for the kernels $F_{2,a}(\bm{p},\bm{k},\eta)$, the second and third line in
Eq. (C) is absent, and only the fourth line contributes. Therefore, the
desired relation Eq. (79) follows directly from the vertex property Eq. (78)
in that case, when proceeding analogously as in the general case.

##  APPENDIX D Treatment of vorticity, vector and tensor modes

In Fourier space, the vorticity $w_{i}$ satisfies the transversality condition
$\bm{k}\cdot\bm{w}=0$, such that only two of its three components are
independent. The same holds for the vector mode $\nu_{i}$ of the second
cumulant, while its tensor mode $t_{ij}$ satisfies
$k_{i}t_{ij}=0,t_{ij}=t_{ji},t_{ii}=0$, corresponding to $3\times 3-(3+3+1)=2$
independent degrees of freedom. In our implementation we use a scheme taking
advantage of these constraints. For example, for the vorticity perturbation
$w_{i}$ at some wavenumber $\bm{k}$, we can choose the two independent degrees
of freedom in the plane perpendicular to $\bm{k}$. There is a remaining
freedom to select a basis in this plane. It turns out that for the one-loop
power spectrum a suitable choice of this basis allows us to track only a
single perturbation mode for the vorticity. The same applies to the vector and
tensor mode, respectively. For the two-loop power spectrum or one-loop
bispectrum, we need to include all independent degrees of freedom. We describe
both algorithms below, and checked their exact equivalence whenever the
simpler one is applicable.

Before that, we discuss some general properties. In this appendix we refer to
scalar modes with the letter $s$ (with
$s\in\\{\delta,\theta,g,\delta\epsilon,A,\pi,\chi\\}$ up to (cum3+)), vector
modes with $v_{i}$ (with $v_{i}\in\\{w_{i},\nu_{i}\\}$ encompassing _both_
vorticity and the vector mode of the velocity dispersion tensor), and tensor
modes with $t_{ij}$. Vector and tensor modes can appear in two different
roles:

  * ($i$)

when computing power spectra of vector or tensor modes themselves (involving
non-linear kernels $F_{n,v_{i}}$ or $F_{n,t_{ij}}$),

  * ($ii$)

as _backreaction_ contributions to the power- (and bi-)spectra of scalar modes
via non-linear kernels $F_{n,s}$.

In a diagrammatic representation of the perturbative solution of the equations
of motion (see e.g. Bernardeau  _et al._ (2002) and App. A), vector and tensor
modes appear as _external, outgoing lines_ in case ($i$) and as _internal_
lines in both ($i$) and ($ii$), while “ingoing” lines (that emerge from
insertions of the initial power spectrum $P_{0}$ in loop diagrams) are always
(growing) scalar modes.

Furthermore, statistical isotropy implies rotational invariance of power
spectra involving vector and tensor modes, which implies that all cross
spectra $P_{v_{i}s}=P_{t_{ij}s}=P_{v_{i}t_{jk}}=0$ vanish. In addition,
rotational invariance requires

$P_{v_{i}v_{j}}(k,\eta)=\frac{1}{2}\left(\delta_{ij}-\frac{k_{i}k_{j}}{k^{2}}\right)\sum_{n}P_{v_{n}v_{n}}(k,\eta)\,,\qquad
P_{t_{ij}t_{ls}}(k,\eta)=\frac{1}{2}P^{T}_{ij,ls}(k)\sum_{m,n}P_{t_{nm}t_{nm}}(k,\eta)\,,$
(152)

where we explicitly indicated the summation over $n,m$ here for clarity (while
being implicitly understood in the rest of this work), and $i,j,l,s$ are not
summed over. The tensor projector $P^{T}_{ij,ls}(k)$ can be found in Garny
_et al._ (2022). These results can be derived by making the ansatz that e.g.
$P_{v_{i}v_{j}}(k,\eta)$ is a general linear combination of all rotationally
covariant objects with two indices built from the vector $\bm{k}$, being
$\delta_{ij}$ and $k_{i}k_{j}$, and using the transversality condition
$\bm{k}\cdot\bm{v}=0$. An analogous derivation holds for the tensor mode and
the cross spectra. Note that the structure of $P_{v_{i}v_{j}}$ applies
correspondingly to the three cases
$P_{w_{i}w_{j}},P_{w_{i}\nu_{j}},P_{\nu_{i}\nu_{j}}$. Altogether, this implies
that it is sufficient to compute the vector and tensor power spectra on the
right-hand sides of Eq. (152), where a summation over all indices is
performed. Our numerical implementations take advantage of this observation.

Figure 19:  Generation of vorticity/vector perturbations at one-loop (left
diagram) and backreaction on scalar power spectra (right diagram). Here
$v_{i}\in\\{w_{i},\nu_{i}\\}$ corresponds to dashed lines, and solid lines
represent any scalar perturbations
$s\in\\{\delta,\theta,g,\delta\epsilon,A,\pi,\chi\\}$ (possibly of different
type for each line). The open square denotes the initial power spectrum, and
the filled circles are non-linear vertices.

###  D.1 Simplified algorithm for the one-loop power spectrum

The algorithm described here can be used for all scalar, vector and tensor
power spectra up to one-loop order, including backreaction of vector and
tensor modes for the former. It exploits that, at one-loop, it is sufficient
to track only one of the two independent degrees of freedom for each vector
perturbation $v_{i}$ and tensor perturbation $t_{ij}$, respectively. The
general case is discussed in Sec. D.2.

When being interested in the one-loop power spectrum, vector and tensor modes
can appear only in two places: ($i$) in the kernel
$F_{2,v_{i}}(\bm{q},\bm{k}-\bm{q},\eta)$ or
$F_{2,t_{ij}}(\bm{q},\bm{k}-\bm{q},\eta)$ for the vector- or tensor power
spectra $P_{v_{i}v_{i}}$ or $P_{t_{ij}t_{ij}}$, respectively, and ($ii$) in
the kernels $F_{3,s}(\bm{k},\bm{q},-\bm{q},\eta)$ contributing to power
spectra of scalar modes. These two possibilities are specific examples for the
generic cases ($i$) and ($ii$) from above.

In the following we discuss an algorithm that is sufficient to compute these
kernels. Let us start with the vector modes $v_{i}\in\\{w_{i},\nu_{i}\\}$. We
generically refer to non-linear vertices $\gamma_{abc}(\bm{k}_{1},\bm{k}_{2})$
with e.g. $a$ representing a vector mode and $b,c$ any two scalar modes
$s,s^{\prime}$ as $\gamma_{v_{i}ss^{\prime}}(\bm{k}_{1},\bm{k}_{2})$, and
analogously for other cases. The only relevant vertex combinations containing
vector modes when computing one-loop power spectra are

$\gamma_{v_{i}ss^{\prime}}(\bm{q},\bm{k}-\bm{q})\gamma_{v_{i}^{\prime}s^{\prime\prime}s^{\prime\prime\prime}}(\bm{q},\bm{k}-\bm{q})\,,\qquad\gamma_{ss^{\prime}v_{i}}(\bm{q},\bm{k}-\bm{q})\gamma_{v_{i}^{\prime}s^{\prime\prime}s^{\prime\prime\prime}}(-\bm{q},\bm{k})\,,$
(153)

which contribute to case ($i$) and ($ii$), respectively (see Fig. 19), and
with summation over the spatial index $i=x,y,z$ in both cases. Here $\bm{k}$
is the external wavenumber, and $\bm{q}$ the loop integration variable. Note
that there are only those two independent directions. In addition, we observe
that all vertices of the form
$\gamma_{v_{i}ss^{\prime}}(\bm{k}_{1},\bm{k}_{2})$ as well as
$\gamma_{sv_{i}s^{\prime}}(\bm{k}_{1},\bm{k}_{2})=\gamma_{ss^{\prime}v_{i}}(\bm{k}_{2},\bm{k}_{1})$
are proportional to $\bm{k}_{1}\times\bm{k}_{2}$. Therefore, the vorticity has
to lie along the axis $\bm{k}\times\bm{q}$. Thus, it is sufficient to track a
single degree of freedom for each vector mode. In practice, we can replace the
three kartesian by a single component, $w_{i}\mapsto w_{\text{eff}}$ and
$\nu_{i}\mapsto\nu_{\text{eff}}$, and replace the factor
$(\bm{k}_{1}\times\bm{k}_{2})_{i}$ contained in the relevant vertices by
$\sigma_{\bm{k}_{1}\bm{k}_{2}}|\bm{k}_{1}\times\bm{k}_{2}|$. Here
$\sigma_{\bm{k}_{1}\bm{k}_{2}}\in\\{\pm 1,0\\}$ keeps track of the correct
sign, and is chosen such that
$(\bm{k}_{1}\times\bm{k}_{2})_{i}=\sigma_{\bm{k}_{1}\bm{k}_{2}}|\bm{k}_{1}\times\bm{k}_{2}|(\bm{k}\times\bm{q})_{i}/|\bm{k}\times\bm{q}|$
for all relevant vectors
$\bm{k}_{1,2}\in\\{\pm\bm{q},\pm\bm{k},\pm(\bm{k}-\bm{q}),\pm(\bm{k}+\bm{q})\\}$
that occur when evaluating the integrand of the one-loop power spectrum,
including the symmetrization $\bm{q}\leftrightarrow-\bm{q}$ to ensure
cancellation of infrared singularities at the integrand level. We use
$\sigma_{\bm{k}_{1}\bm{k}_{2}}=\text{sgn}((\bm{k}_{1}\times\bm{k}_{2})\cdot(\bm{k}\times\bm{q}))$.

Altogether, the effective single-component vorticity and vector modes
correctly take into account the backreaction on the scalar kernels
$F_{3,s}(\bm{k},\bm{q},-\bm{q},\eta)$, and allow for the computation of the
vector power spectra according to
$P_{w_{i}w_{i}}=P_{w_{\text{eff}}w_{\text{eff}}},P_{w_{i}\nu_{i}}=P_{w_{\text{eff}}\nu_{\text{eff}}},P_{\nu_{i}\nu_{i}}=P_{\nu_{\text{eff}}\nu_{\text{eff}}}$.

For the tensor mode, the relevant combinations of vertices are analogous to
Eq. (153) with $v_{i},v_{i}^{\prime}$ replaced by $t_{ij},t_{ij}^{\prime}$. We
observe that for vertices of the form
$\gamma_{t_{ij}ss^{\prime}}(\bm{k}_{1},\bm{k}_{2})$ the dependence on the
indices $i,j$ is universally given by the factor

$f_{ij}(\bm{k}_{1},\bm{k}_{2})\equiv\delta_{ij}-(\bm{k}_{1}+\bm{k}_{2})_{i}(\bm{k}_{1}+\bm{k}_{2})_{j}/(\bm{k}_{1}+\bm{k}_{2})^{2}-2(\bm{k}_{1}\times\bm{k}_{2})_{i}(\bm{k}_{1}\times\bm{k}_{2})_{j}/(\bm{k}_{1}\times\bm{k}_{2})^{2}\,,$
(154)

which satisfies the transversality and trace conditions
$0=(\bm{k}_{1}+\bm{k}_{2})_{i}f_{ij}(\bm{k}_{1},\bm{k}_{2}),0=f_{ii}$.
Furthermore, for all vertices of the form
$\gamma_{st_{ij}s^{\prime}}(\bm{k}_{1},\bm{k}_{2})=\gamma_{ss^{\prime}t_{ij}}(\bm{k}_{2},\bm{k}_{1})$
the dependence on $i,j$ is captured by
$g_{ij}(\bm{k}_{1},\bm{k}_{2})\equiv(\bm{k}_{2})_{i}(\bm{k}_{2})_{j}$. We find
that for the one-loop power spectrum we can again represent the tensor mode by
a single degree of freedom, $t_{ij}\mapsto t_{\text{eff}}$, with vertices
related to the full ones via
$\gamma_{t_{ij}ss^{\prime}}(\bm{k}_{1},\bm{k}_{2})=f_{ij}(\bm{k}_{1},\bm{k}_{2})\gamma_{t_{\text{eff}}ss^{\prime}}(\bm{k}_{1},\bm{k}_{2})$
and
$k_{2}^{2}\sin^{2}(\bm{k}_{1},\bm{k}_{2})\gamma_{st_{ij}s^{\prime}}(\bm{k}_{1},\bm{k}_{2})=g_{ij}(\bm{k}_{1},\bm{k}_{2})\gamma_{st_{\text{eff}}s^{\prime}}(\bm{k}_{1},\bm{k}_{2})$,
where
$\sin^{2}(\bm{k}_{1},\bm{k}_{2})=(\bm{k}_{1}\times\bm{k}_{2})^{2}/(k_{1}^{2}k_{2}^{2})$.
With these definitions one has

$\displaystyle\gamma_{t_{ij}ss^{\prime}}(\bm{q},\bm{k}-\bm{q})\gamma_{t_{ij}^{\prime}s^{\prime\prime}s^{\prime\prime\prime}}(\bm{q},\bm{k}-\bm{q})$
$\displaystyle=$ $\displaystyle
2\gamma_{t_{\text{eff}}ss^{\prime}}(\bm{q},\bm{k}-\bm{q})\gamma_{t_{\text{eff}}^{\prime}s^{\prime\prime}s^{\prime\prime\prime}}(\bm{q},\bm{k}-\bm{q})\,,$
$\displaystyle\gamma_{ss^{\prime}t_{ij}}(\bm{q},\bm{k}-\bm{q})\gamma_{t_{ij}^{\prime}s^{\prime\prime}s^{\prime\prime\prime}}(-\bm{q},\bm{k})$
$\displaystyle=$
$\displaystyle\gamma_{ss^{\prime}t_{\text{eff}}}(\bm{q},\bm{k}-\bm{q})\gamma_{t_{\text{eff}}^{\prime}s^{\prime\prime}s^{\prime\prime\prime}}(-\bm{q},\bm{k})\,.$
(155)

The second line implies that the effective tensor mode correctly takes into
account the backreaction of tensor perturbations on the scalar kernels
$F_{3,s}(\bm{k},\bm{q},-\bm{q},\eta)$. Using the first line we find that the
tensor power spectrum at one-loop order can be computed according to
$P_{t_{ij}t_{ij}}(k,\eta)=2P_{t_{\text{eff}}t_{\text{eff}}}(k,\eta)$.

###  D.2 General case

For situations where the simplified treatment described above is not
applicable, we also implemented a full treatment of vorticity and vector mode
perturbations taking both of the two independent degrees of freedom for each
of them into account. This is relevant in particular for the two-loop
vorticity power spectrum (see Sec. VI.4) and the one-loop matter density
bispectrum (see Sec. VI.3). In practice, to make use of the constraint
$\bm{p}\cdot\bm{v}=0$ for a vector mode $v_{i}\in\\{w_{i},\nu_{i}\\}$ of
wavenumber $\bm{p}$, we project $v_{i}$ on a basis that depends on $\bm{p}$.
In particular, we use the $\bm{p}$-dependent orthogonal basis vectors

$\bm{b}_{\bm{p}1}\equiv{\cal N}_{\bm{p}1}(\bm{p}\times\bm{Q}_{\bm{p}}),\
\bm{b}_{\bm{p}2}\equiv{\cal
N}_{\bm{p}2}(\bm{p}\times(\bm{p}\times\bm{Q}_{\bm{p}}))={\cal
N}_{\bm{p}2}\left(\bm{p}\,(\bm{p}\cdot\bm{Q}_{\bm{p}})-\bm{Q}_{\bm{p}}\,p^{2}\right),\
\bm{b}_{\bm{p}3}\equiv{\cal N}_{\bm{p}3}\,\bm{p}\,,$ (156)

where $\bm{Q}_{\bm{p}}$ is an a priori arbitrary reference vector that can be
chosen by convenience for each $\bm{p}$ with the only condition that
$\bm{p}\times\bm{Q}_{\bm{p}}\not=0$. Furthermore, ${\cal N}_{\bm{p}1,2,3}$ are
normalization factors. A natural choice is such that the basis vectors have
norm equal to unity, but this is not required in our implementation. We refer
to this basis as the _transverse basis_ , and decompose the vorticity for
wavevector $\bm{p}$ as

$\bm{w}=w_{\bm{p}1}\bm{b}_{\bm{p}1}+w_{\bm{p}2}\bm{b}_{\bm{p}2}\,,$ (157)

where $w_{\bm{p}1,2}$ are the two tranverse vorticity components, and the
projection along $\bm{b}_{\bm{p}3}$ vanishes by construction. The
decomposition of the vector mode is analogous. In the list of perturbation
variables, we correspondingly include $w_{\bm{p}\alpha}$ and
$\nu_{\bm{p}\alpha}$ for $\alpha=1,2$. The linear part of their evolution is
not affected by the basis choice. For the non-linear part, the vertices in the
transverse basis can be obtained from the generic kartesian ones via

$\gamma_{aw_{\bm{p}\alpha}c}(\bm{p},\bm{q})=\gamma_{aw_{i}c}(\bm{p},\bm{q})\bm{b}_{\bm{p}\alpha,i},\
\gamma_{abw_{\bm{q}\alpha}}(\bm{p},\bm{q})=\gamma_{abw_{i}}(\bm{p},\bm{q})\bm{b}_{\bm{q}\alpha,i},\
\gamma_{w_{\bm{k}\alpha}bc}(\bm{p},\bm{q})=\bm{b}_{\bm{p}+\bm{q}\alpha,i}/|\bm{b}_{\bm{p}+\bm{q}\alpha}|^{2}\gamma_{w_{i}bc}(\bm{p},\bm{q})\,,$
(158)

where summation over $i=x,y,z$ is implied. Note that the appropriate
transverse basis for each of the wavevectors $\bm{p}$, $\bm{q}$ or
$\bm{p}+\bm{q}$ is used in the three cases. For vertices involving the vector
mode $\nu_{i}$ the projection is analogous, and for vertices where more than
one of the indices $a,b,c$ corresponds to a vorticity or vector mode, the
projection is applied for each of them correspondingly. For example,

$\gamma_{w_{\bm{p}+\bm{q}\alpha}w_{\bm{p}\beta}w_{\bm{q}\gamma}}(p,q)=b_{\bm{p}+\bm{q}\alpha,i}/|\bm{b}_{\bm{p}+\bm{q}\alpha}|^{2}\gamma_{w_{i}w_{j}w_{k}}(\bm{p},\bm{q})\bm{b}_{\bm{p}\beta,j}\bm{b}_{\bm{q}\gamma,k}\,,$
(159)

with summation over $i,j,k=x,y,z$ and for $\alpha,\beta,\gamma=1,2$. Therefore
the $www$ vertex corresponds in the transverse basis to $2^{3}=8$ distinct
vertices, while in the kartesian basis it would have in general $3^{3}=27$
components. This illustrates the advantage of using the transverse basis.

We note that in the transverse basis, and using well-known relations from
vector algebra, the projected vertices can be reduced to expressions that
depend solely on scalar products $\bm{p}\cdot\bm{q}$,
$\bm{p}\cdot\bm{Q}_{\bm{r}}$ and $\bm{q}\cdot\bm{Q}_{\bm{r}}$, as well as on

$D_{\bm{r}\bm{s}}\equiv\text{det}(\bm{s},\bm{r},\bm{Q}_{r})=(\bm{r}\times\bm{Q}_{r})\cdot\bm{s}\,,$
(160)

with $\bm{r},\bm{s}\in\\{\bm{p},\bm{q},\bm{p}+\bm{q}\\}$. This is advantageous
when following the algorithm outlined in Blas  _et al._ (2013), and we
precompute these values for a given configuration of wavevectors corresponding
to a single evaluation of the integrand for the loop evaluation. In
particular, the relevant elementary wavevectors that may appear are the
external wavevector $\bm{k}$ for the power spectrum (or
$\bm{k}_{1},\bm{k}_{2}$ for the bispectrum) as well as loop wavenumbers
$\bm{Q}_{1},\dots,\bm{Q}_{L}$ at $L$-loop order. The most general set of
vectors that may appear as arguments of non-linear kernels as well as of
vertices has then the form $\sum_{n=1}^{L}c_{n}Q_{n}+c_{L+1}\bm{k}$ (or
$\sum_{n=1}^{L}c_{n}Q_{n}+c_{L+1}\bm{k}_{1}+c_{L+2}\bm{k}_{2}$ for the
bispectrum), with some coefficients $c_{n}=0,\pm 1$. For each of those vectors
$\bm{p}$, we define the associated vector $\bm{Q}_{\bm{p}}$ entering the
definition of the corresponding transverse basis as the vector within the set
${\cal B}=\\{\bm{k},\bm{Q}_{1},\dots,\bm{Q}_{L}\\}$ (or ${\cal
B}=\\{\bm{k}_{1},\bm{k}_{2},\bm{Q}_{1},\dots,\bm{Q}_{L}\\}$ for the
bispectrum) that has the largest projection within the plane perpendicular to
$\bm{p}$, i.e. for which
$s_{\bm{p}\bm{Q}}^{2}=1-\frac{(\bm{p}\cdot\bm{Q})^{2}}{\bm{p}^{2}\bm{Q}^{2}}$
is maximal out of all $\bm{Q}\in{\cal B}$. This allows for an efficient
implementation of vertices involving vector modes within the algorithm used in
this work, based on Blas  _et al._ (2013, 2015); Floerchinger  _et al._
(2019); Garny and Taule (2021).

##  APPENDIX E Rescaling for power law universe

For a scaling universe with linear input power spectrum $P_{0}(k)=Ak^{n_{s}}$
and EdS background, the average value of the velocity dispersion is given by
$\epsilon(\eta)=\epsilon_{0}e^{\alpha\eta}$ with $\alpha=4/(n_{s}+3)$, and the
dimensionless ratio $\bar{\omega}=\omega(\eta)/\epsilon^{2}(\eta)$ involving
the fourth cumulant expectation value is constant (as are dimensionless ratios
of higher order cumulants Garny  _et al._ (2022)). The only dimensionful
scales are therefore $k_{\sigma}\equiv 1/\sqrt{\epsilon_{0}}$ and
$k_{\text{nl}}\equiv 1/(4\pi A)^{1/(n_{s}+3)}$. We indicate the dependence of
power spectra on $\epsilon_{0}$ in this appendix using the notation
$P_{ab}(k,\eta;\epsilon_{0})$, for any perturbation modes $a$ and $b$ (e.g.
$a=b=\delta$ for the density spectrum), and taking $A$ as given.

Scaling symmetry can be used to rescale power pectra computed for a given
reference value $\epsilon_{0}^{\text{ref}}$ to any other value $\epsilon_{0}$,
while keeping all dimensionless ratios of higher cumulants fixed, specifically
$\bar{\omega}$ (the dependence on which we suppress for brevity). To see this,
we consider for a moment power spectra
$P_{\bar{a}\bar{b}}(k,\eta;\epsilon_{0})=P_{ab}(k,\eta;\epsilon_{0})/\epsilon^{d_{a}+d_{b}}(\eta)$
of dimensionless perturbation variables $\bar{a}=a/\epsilon^{d_{a}}(\eta)$ and
$\bar{b}=b/\epsilon^{d_{b}}(\eta)$, with appropriate powers $d_{a},d_{b}$. For
example, $d_{\delta}=d_{\theta}=d_{A}=d_{w_{i}}=0$,
$d_{g}=d_{\delta\epsilon}=d_{\nu_{i}}=d_{t_{ij}}=d_{\pi}=d_{\chi}=1$. Then the
dimensionless power spectrum $4\pi k^{3}P_{\bar{a}\bar{b}}$ can depend only on
dimensionless ratios of $k$, $k_{\sigma}$ and $k_{\text{nl}}$, which can be
taken to be $k/k_{\sigma}$ and $k_{\sigma}/k_{\text{nl}}$. In addition, within
perturbation theory, the $L$-loop contribution involves exactly $L+1$ powers
of $P_{0}$, i.e. is proportional to $A^{L+1}$. This implies that $4\pi
k^{3}P^{L-\text{loop}}_{\bar{a}\bar{b}}(k,\eta;\epsilon_{0})=(k_{\sigma}/k_{\text{nl}})^{(n_{s}+3)(L+1)}\bar{\Delta}_{\bar{a}\bar{b}}^{L}(k/k_{\sigma},\eta)$,
with some dimensionless function
$\bar{\Delta}_{\bar{a}\bar{b}}^{L}(k/k_{\sigma},\eta)$ that depends on $k$
only via the ratio $k/k_{\sigma}$. From this we obtain the rescaling relation
for the original power spectrum,

$P_{ab}^{L-\text{loop}}(k,\eta;\epsilon_{0})=\left(\frac{\epsilon_{0}^{\text{ref}}}{\epsilon_{0}}\right)^{\frac{(n_{s}+3)(L+1)}{2}-\frac{3}{2}-d_{a}-d_{b}}P_{ab}^{L-\text{loop}}(k\times(\epsilon_{0}/\epsilon_{0}^{\text{ref}})^{1/2},\eta;\epsilon_{0}^{\text{ref}})\,.$
(161)

Similarly, for the bispectrum we find

$B_{abc}^{L-\text{loop}}(k_{1},k_{2},k_{3},\eta;\epsilon_{0})=\left(\frac{\epsilon_{0}^{\text{ref}}}{\epsilon_{0}}\right)^{\frac{(n_{s}+3)(L+2)}{2}-\frac{6}{2}-d_{a}-d_{b}-d_{c}}B_{abc}^{L-\text{loop}}(k_{1}^{\prime},k_{2}^{\prime},k_{3}^{\prime},\eta;\epsilon_{0}^{\text{ref}})\,,$
(162)

where
$k_{i}^{\prime}=k_{i}\times(\epsilon_{0}/\epsilon_{0}^{\text{ref}})^{1/2}$ for
$i=1,2,3$.

In addition, the dependence on time (or equivalently $\eta$) is also fixed by
scaling symmetry. To see this it is useful to introduce generalized time-
dependent scales $k_{\sigma}(\eta)\equiv 1/\sqrt{\epsilon(\eta)}$,
$k_{\text{nl}}(\eta)\equiv 1/(4\pi Ae^{2\eta})^{1/(n_{s}+3)}$. The
dimensionless power spectrum $4\pi
k^{3}P_{\bar{a}\bar{b}}(k,\eta;\epsilon_{0})$ can depend on time only via the
ratios $k/k_{\sigma}(\eta)$ and $k_{\sigma}(\eta)/k_{\text{nl}}(\eta)$.
Observing that the latter ratio is constant in time implies that
$\bar{\Delta}_{\bar{a}\bar{b}}^{L}(k/k_{\sigma},\eta)=\bar{\Delta}_{\bar{a}\bar{b}}^{L}(k/k_{\sigma}(\eta))$
is a function of a single variable. This yields
$P_{ab}(k,\eta;\epsilon_{0})=(\epsilon(\eta^{\prime})/\epsilon(\eta))^{-3/2-d_{a}-d_{b}}P_{ab}(k^{\prime},\eta^{\prime};\epsilon_{0})$
and
$B_{abc}(k_{1},k_{2},k_{3},\eta;\epsilon_{0})=(\epsilon(\eta^{\prime})/\epsilon(\eta))^{-6/2-d_{a}-d_{b}-d_{c}}B_{abc}(k_{1}^{\prime},k_{2}^{\prime},k_{3}^{\prime},\eta^{\prime};\epsilon_{0})$
where now
$k^{\prime}=k\times(\epsilon(\eta)/\epsilon(\eta^{\prime}))^{1/2}=ke^{2(\eta-\eta^{\prime})/(n_{s}+3)}$
and analogously for $k_{i}^{\prime}$.

##  APPENDIX F N-body simulations

Our main suite of N-body simulations (see Tab. 5) consists of 4 realizations
of $512^{3}$ particles for each power-law spectral index $n_{s}=-1,0,1,2$, run
with $\Omega_{m}=1$ using the Gadget code Springel (2005). The four
realizations in each case correspond to two sets of fixed amplitude initial
conditions with opposite phases Angulo and Pontzen (2016), to cancel Gaussian
cosmic variance. The initial conditions are set at scale factor $a_{\rm
ic}=0.001$ and evolved until $a=1$ with outputs at $a=10^{s}$ with $s=-2$ to
$s=0$ in steps of $\Delta s=0.2$. The setup of initial conditions necessarily
induces transients in subsequent evolution, due to imperfect dynamics
Scoccimarro (1998b); Valageas (2002); Crocce  _et al._ (2006); Jenkins (2010);
Reed  _et al._ (2013); Michaux  _et al._ (2021) (the use of perturbation
theory to compute initial displacements) and special configurations of
discrete particles Marcos  _et al._ (2006); Joyce and Marcos (2007a, b);
Garrison  _et al._ (2016) (cartesian grid out of which the initial
displacements are imposed) that are not exact solutions of the equations of
motion.

For our application of interest, namely testing 3D gravitational clustering in
the regime where orbit-crossing is most significant, we would like to focus on
spectra with blue spectral indices, hence our choice of pushing $n_{s}$ to be
as high as $n_{s}=+2$, which as far as we know has never before been
considered in the literature of scale-free 3D cosmological simulations Colombi
_et al._ (1996); Peacock and Dodds (1996); Smith  _et al._ (2003); Bagla  _et
al._ (2009); Widrow  _et al._ (2009); Joyce  _et al._ (2021); Leroy  _et al._
(2021). This gives rise to some challenges from the numerical simulation point
of view. In particular, since we use perturbation theory to set up initial
conditions, we would like to have the maximum amplitude of fluctuations to be
initially less than unity, which for such blue spectra makes the fluctuation
amplitude at the scale of the box very small. This in turn makes following the
evolution of fluctuations more difficult by the tree algorithm due to
accumulation of numerical errors (e.g. Reed  _et al._ (2013); Angulo and Hahn
(2022). To strike a balance, the amplitude of initial fluctuations at the
Nyquist frequency of the particle grid is an increasing function of $n_{s}$,
with $\Delta_{\rm Nyq}=(k_{\rm Nyq}/k_{\rm nl})^{n_{s}+3}=0.01,0.01,0.05,0.25$
for $n_{s}=-1,0,1,2$. This gives an initial amplitude at the scale of the box
of $\Delta_{\rm box}=1.5\times 10^{-7},6\times 10^{-10},1.2\times
10^{-11},2.3\times 10^{-13}$, respectively. Even then, we see some numerical
noise in the evolution of the large-scale modes down to late times,
particularly for the bispectrum when $n_{s}=2$ (see Fig. 14).

This limitation imposes constraints on the order of perturbation theory that
can be safely used in setting up initial conditions. Indeed, at the
$\Delta_{\rm Nyq}$ values mentioned above, our linear solutions to the
equations of motion with dispersion and higher cumulants tell us that already
the corrections over the unperturbed initial spectrum $P_{0}\propto k^{n_{s}}$
are about $0,3,15,30\%$ at the Nyquist frequency for $n_{s}=-1,0,1,2$
respectively at the initial conditions (see Fig. 10). This means that for blue
spectral indices $n_{s}=1,2$, Lagrangian perturbation theory convergence is
doubtful, as one would have to take into account dispersion and higher
cumulants. For this reason in such cases we set initial conditions using the
Zel’dovich approximation (ZA), rather than second-order Lagrangian
perturbation theory (2LPT).

The use of ZA initial conditions introduces violations of self-similarity at
early times due to transients, which we correct using tree-level standard
perturbation theory (SPT), tested against $n_{s}=-1,0$ 2LPT initial conditions
simulations for which we have also run their ZA counterparts. For the
bispectrum for equilateral triangles, tree-level SPT gives Scoccimarro (1998b)

$B_{\rm eq}(k,a)=\Big{\\{}\frac{12}{7}-\frac{27}{20}\Big{(}{a_{\rm ic}\over
a}\Big{)}+\frac{27}{70}\Big{(}{a_{\rm ic}\over
a}\Big{)}^{7/2}\Big{\\}}\Big{[}P_{\rm lin}^{\rm SPT}(k,a)\Big{]}^{2}\ ,$ (163)

which at the initial conditions ($a=a_{\rm ic}$) gives the ZA tree-level
bispectrum, and for $a\gg a_{\rm ic}$ gives the tree-level equilateral
bispectrum of the exact dynamics in the absence of dispersion (which is
correctly predicted by 2LPT). In between, transients (spurious $a/a_{\rm ic}$
dependence in Eq. 163) induce violations to self-similar evolution. At early
times, when our simulations are weakly non-linear, we find that Eq. (163)
matches the differences between the ZA and 2LPT initial conditions simulations
very well, as expected, and in agreement with previous tests of higher-order
statistics in the literature Scoccimarro (1998b); Crocce  _et al._ (2006). In
fact, we have checked this for all triangle shapes, not just equilateral. The
use of Eq. (163) at late times to compute the relative corrections due to
transients is also justified in practice because, although tree-level SPT is
not valid, the size of the transient corrections is rather small, so a larger
error in estimating it is acceptable compared to the statistical errors. More
specifically, we find that the maximum deviation between ZA initial conditions
simulations corrected for transients by SPT compared to the 2LPT initial
conditions simulations for all triangle shapes and simulation outputs is about
$1.5\%$ for $n_{s}=-1$ and $0.5\%$ for $n_{s}=0$, suggesting our application
to $n_{s}=1,2$ should yield sub-percent systematic errors.

Table 5: N-Body simulations as a function of spectral index. $n_{s}$ | $N_{\rm particles}$ | $N_{\rm realizations}$ | Initial Conditions
---|---|---|---
$2$ | $256^{3},512^{3}$ | 6 (3 paired fixed-amplitude) | ZA
$1$ | $256^{3},512^{3}$ | 4 (2 paired fixed-amplitude) | ZA
$0$ | $256^{3},512^{3}$ | 4 (2 paired fixed-amplitude) | ZA
$0$ | $384^{3},768^{3}$ | 1 fixed amplitude | ZA
$0$ | $512^{3}$ | 4 (2 paired fixed-amplitude) | 2LPT
$-1$ | $512^{3}$ | 4 (2 paired fixed-amplitude) | ZA
$-1$ | $512^{3}$ | 4 (2 paired fixed-amplitude) | 2LPT

The situation for power spectrum measurements is less challenging, since ZA
initial conditions reproduce the correct linear spectrum in the absence of
dispersion, with transients entering only trough loop corrections in the
context of SPT Crocce  _et al._ (2006). By comparing our ZA initial conditions
simulations to those run with 2LPT for the cases $n_{s}=-1,0$, we see that
correcting power spectrum measurements for transients is in practice not
necessary, as deviations are at most $0.5\%$ for $n_{s}=-1$ and $0.3\%$ for
$n_{s}=0$ once the power spectra from different outputs are binned into
functions of $k/k_{\rm nl}$.

Apart from transients, there are other spurious effects that break self-
similarity Colombi  _et al._ (1996). N-body simulations cannot follow modes
with wavelength larger than the size of the box, and at small scales the
gravitational force is softened. The former means in particular that at late
times, when the non-linear scale becomes comparable to the box size (e.g.
$k_{\rm nl}L\mathrel{\mathchoice{\lower
3.6pt\vbox{\halign{$\mathsurround=0pt\displaystyle\hfill#\hfil$\cr<\crcr\sim\crcr}}}{\lower
3.6pt\vbox{\halign{$\mathsurround=0pt\textstyle\hfill#\hfil$\cr<\crcr\sim\crcr}}}{\lower
3.6pt\vbox{\halign{$\mathsurround=0pt\scriptstyle\hfill#\hfil$\cr<\crcr\sim\crcr}}}{\lower
3.6pt\vbox{\halign{$\mathsurround=0pt\scriptscriptstyle\hfill#\hfil$\cr<\crcr\sim\crcr}}}}0.1$),
the simulation cannot follow the true dynamics as mode-coupling to relevant
modes is missing. In addition, particles in simulations can be thought of as
having an effective size equal to the softening length, and to be in the fluid
limit (avoid spurious discreteness) requires that the non-linear scale be
larger than the mean interparticle separation, or equivalently that the
amplitude of fluctuations at the Nyquist frequency of the particles be larger
than unity. This is of course the opposite of the requirement that
perturbation theory be valid (as required for initial conditions). In practice
this means that between the initial conditions and the first output that can
be analyzed one must leave sufficient time for both conditions to be met. This
also minimizes transients induced by simplified dynamics in initial conditions
generation and associated initial grid effects. Therefore, to obtain
measurements as a function of $k/k_{\rm nl}$ with a broad range of values, we
use measurements starting at a scale factor 25 times larger than at the
initial conditions.

Finally, another issue that can potentially impact our simulations is the use
of fixed-amplitude initial conditions to reduce cosmic variance Angulo and
Pontzen (2016). Because the amplitude of the Fourier modes are fixed, these
initial conditions are actually non-Gaussian. That’s precisely the point when
it comes to the initial four-point function, as its modification leads to
Gaussian cosmic variance cancellation. However this also affects the
bispectrum. At tree-level in SPT, the bispectrum gets modified by order unity
only for triangles of the form $B(\bm{k},\bm{k},-2\bm{k})$, which involves the
initial four-point function of pairs of opposite modes (same that enters into
the power spectrum cosmic variance). Since we present measurements of the
bispectrum for equilateral triangles in this paper, this is of no concern,
except for the fact that we use a wide bin in $k$ with $\Delta k=6k_{F}$, with
$k_{F}$ the fundamental model of the simulation box, and the first bin at
$k=6\pm\Delta k/2$ does have contributions from such fundamental triangles.
Therefore, to avoid this issue we do not use the first two equilateral
triangle bins in any of our bispectrum results.

At one-loop order, however, it’s easy to see that these fixed-amplitude
initial conditions induce bispectrum modifications to triangles of all shapes
(as well as to the power spectrum). This is so because given any choice of
external momenta (triangle shape) one can always match the configuration of
initial four-point functions affected in the expectation values over linear
fields for a suitable choice of the loop momentum that is being integrated
over. However, since only a special configuration is involved, the loop
integral collapses to a particular case of discrete modes. By dimensional
analysis the loop then gives a contribution proportional to $k_{F}^{3}P_{\rm
lin}(k)\propto(k/k_{\rm nl})^{n_{s}}(k_{F}/k_{\rm nl})^{3}$, which means that
fixed-amplitude initial conditions induce violations of self-similarity
through loop corrections. Following Oddo and Sefusatti (2022), we computed
these for the power spectrum and bispectrum, but found that these violations
are too small (sub-percent even for the latest outputs used) to affect our
measurements.

To calculate error bars for our power spectrum and bispectrum measurements, we
first build an independent realization from the two fixed-amplitude initial
conditions simulations with opposite phases, and do the same for the second
random seed for each spectral index. These two independent realizations
contain several outputs (different scale factors), whose measurements (power
spectrum and bispectrum) are then binned into functions of $k/k_{\rm nl}$ by
self-similarity. The error bars are then computed from the scatter of these
measurements at fixed $k/k_{\rm nl}$. Thus our error bars contain scatter from
two independent random seeds, as well as scatter from several different
outputs that land at the same $k/k_{\rm nl}$. To illustrate the extent to
which self-similarity is satisfied for the density power spectrum, we show in
Fig. 20 the density power spectrum in units of the initial spectrum for the
four spectral indices ($n_{s}=-1$ to $n_{s}=2$) in terms of the self-similar
variable $k/k_{\rm nl}$. Different outputs (each corresponding to a different
color) lie on top of each other, as expected from self-similar scaling. In
this plot, each output corresponds to the average of one fixed-amplitude
realization over the two opposite phases, and we have removed the high-$k$
portion of the spectrum when violation of self-similarity occurs (e.g. due to
the softening length of the simulations) for clarity.

##  APPENDIX G The halo mass function

###  G.1 Measuring the Mass Function

To measure the mass function in our simulations, we run the Rockstar halo
finder Behroozi  _et al._ (2013). We use the virial mass definition, and
consider only halos with more than 50 particles. The halos are then binned in
200 logarithmically spaced bins, and the fixed-amplitude realizations of
opposite phases are averaged together. This gives us 2 paired fixed-amplitude
realizations for each output which are then combined using self-similarity.

The scatter among these is used to compute the error bars, and a fit of the
form given by Eq. (120) is used to find the best fit parameters for $A$, $a$
and $p$ which are given in Tab. 6. Figure 8 shows the resulting fits as a
function of the ratio of mass to the non-linear mass. We have verified
comparing the halo mass function with 2LPT initial conditions to those with ZA
initial conditions for $n_{s}=-1,0$ (see Tab. 5) that transients from ZA
initial conditions do not affect our mass function fits. Figure 21 shows how
well our simulations obey self-similarity for the mass function. For each
spectral index, different colors denote different outputs. For this figure, we
have binned halos into 50 logarithmically spaced bins and suppressed error
bars for clarity. Each line shows the average over two paired fixed-amplitude
realizations.

###  G.2 From Isothermal to NFW halos

The contribution to the expectation value of the stress tensor $\epsilon$ of
an NFW halo can be written in terms of its isothermal halo counterpart as Eq.
(128), where $I(\beta,c)$ is a 2D integral given by Eq. (V.2) that represents
the correction over the isothermal value. Here we present the analytic result
for this integral in the two cases of interest, $\beta=0$ (isotropic) and
$\beta=1/2$ (radially biased dispersion). We have, as a function of halo
concentration $c$,

$\displaystyle I(\beta=0,c)$ $\displaystyle=$ $\displaystyle{g(c)\over
20c^{2}}\Big{[}c\\{6-c[6+c(25+2c(59+35c))]+c^{3}[15+2c(12+5c)]\pi^{2}\\}-c^{4}[15+2c(12+5c)]\ln(1+c)$
(164)
$\displaystyle+(1+c)\\{-6+c(1+c)[12+c(3+2c(-28+5c))]\\}\ln(1+c)+3\\{-1+c^{4}[15+2c(12+5c)]\\}[\ln(1+c)]^{2}$
$\displaystyle+6c^{4}[15+2c(12+5c)]{\rm Li}_{2}(-c)\Big{]}$

and

$\displaystyle I(\beta=1/2,c)$ $\displaystyle=$ $\displaystyle{2g(c)\over
45c^{2}}\Big{[}c^{2}\\{-3+9c[7+4c]-c[10+3c(5+2c)]\pi^{2}\\}+6c(1+c)[1+6c(1+c)]\ln(1+c)$
(165)
$\displaystyle-3(1+c)^{3}(1-3c+6c^{2})[\ln(1+c)]^{2}-6c^{3}[10+3c(5+2c)]{\rm
Li}_{2}(-c)\Big{]}$

respectively. In these expressions ${\rm Li}_{2}$ represents the dilogarithm
function.

##  APPENDIX H Measuring Velocity Divergence and Vorticity

Measuring the velocity field divergence and vorticity is somewhat challenging
in N-body simulations, as the velocity field is sampled by particles (which
are themselves a sampling of the phase-space distribution function $f$). This
means that in low-density regions the sampling is sparser than in high-density
regions. In particular, in voids there is little velocity information at
scales smaller than the void size. One would like to construct a grid where
the velocity field is defined, then Fourier transformed, and
divergence/vorticity and their spectra computed. But this can only be done at
all grid points if there is always a contribution from nearby particles. In
practical terms, the issue manifests itself when defining the peculiar
velocity field,

$\bm{v}(\bm{x})=\frac{1}{1+\delta}\int
d^{3}p\,\frac{\bm{p}}{a}\,f(\bm{x},\bm{p})\,.$ (166)

because at points $\bm{x}$ where $f$ vanishes (and $\delta=-1$), i.e. voids,
this leads to a 0/0 indeterminacy. A number of strategies have been developed
in the literature to overcome this limitation, among them Voronoi and Delaunay
tessellations Bernardeau and van de Weygaert (1996); Romano-Díaz and van de
Weygaert (2007); Pueblas and Scoccimarro (2009), phase-space constructions
Hahn  _et al._ (2015); Buehlmann and Hahn (2019) and rescaling earlier time
velocities Jelic-Cizmek  _et al._ (2018). In this work, we have developed an
alternative approach based on a multigrid method.

To avoid the indeterminacy problem, we start by interpolating the particle
velocities to a coarse grid (e.g. $64^{3}$), where all grid points receive
velocity information. To do this, we use the Piecewise Cubic Spline (PCS)
algorithm discussed in Sefusatti  _et al._ (2015) which employs fourth-order
interpolation in interlaced grids to reduce aliasing. Fourth-order
interpolation has a broad kernel and thus helps in giving velocity information
to all grid points. Then, we recursively make the grid finer by factors of two
in each dimension ($128^{3},256^{3}$, etc). As the grid becomes finer,
eventually there will be grid points (particularly at late times) for which
the indeterminacy happens somewhere inside a void, even with the PCS
algorithm. In that case, we look up the information from the coarser grid, and
perform an interpolation from the nearby points in the coarser grid one level
above (which has no indeterminacy problem) to the finer grid. In this way we
recursively make the grid finer and finer to the desired level, i.e. enough to
resolve the vorticity power spectrum peak at small scales (see below), ending
up at a $1024^{3}$ grid which is then Fourier transformed. The divergence and
vorticity are then calculated in Fourier space, and their power spectrum
computed as usual. Note that as we make the grid finer, if no grid points have
indeterminacy the coarser level grid information is never used, and gets
discarded. At each stage in the algorithm there are thus only two concurrent
grids, and the finer grid becomes the coarser grid in the next iteration of
the recursive algorithm.

A second issue is that the velocity field, as seen from Eq. (166), is the
ratio of two interpolated quantities (momentum over density). Normally, it is
easy to correct for the interpolation scheme used by dividing by the
corresponding Fourier transform of the interpolation window in Fourier space.
This corrects frequencies close to the Nyquist frequency of the grid, which
get damped the most. However, since the velocity is a ratio of interpolated
quantities it is difficult to reliably correct for the interpolation kernel
(i.e. correction of the interpolated density is sensitive to noise), although
being a ratio means that the damping caused by the interpolation window is not
as severe. In practice we minimize this issue by using frequencies up to one
tenth of the Nyquist frequency, comparison between results for different grid
sizes show that this leads to no noticeable window interpolation effects.

The most challenging of our measurements correspond to the vorticity power
spectrum, which shows significant deviations of self-similarity except at late
times, as we detail below. For the velocity divergence power spectrum, Fig. 22
shows a test of self-similarity by plotting each output for a given spectral
index in a different color. We see comparing to Fig. 20 that deviations from
self-similarity are more significant than for the density power spectrum case,
as expected given the challenges discussed above. However, in comparison to
the vorticity power spectrum (see Fig. 23), self-similarity of the velocity
divergence power spectrum is fairly well satisfied.

As has been known from CDM simulations, the vorticity power spectrum is
sensitive to mass resolution, and poor mass resolution results in an
artificial amplification of the vorticity power Pueblas and Scoccimarro
(2009); Hahn  _et al._ (2015); Jelic-Cizmek  _et al._ (2018). For CDM initial
conditions, the vorticity power spectrum typically has a peak at $k_{\rm
peak}\sim 1h/\mathrm{Mpc}$ at $z=0$, and at higher redshifts this peak shifts
to higher $k_{\rm peak}$ as expected. A reasonable criterion for convergence
is to characterize the mass resolution (in terms of the mean interparticle
separation $\lambda\equiv L_{\rm box}/N_{\rm par}$) needed to “resolve” this
peak, i.e. the characteristic scale where most vorticity is produced,
extracted from the convergence studies mentioned above. We obtain

$\lambda\,k_{\rm peak}\mathrel{\mathchoice{\lower
3.6pt\vbox{\halign{$\mathsurround=0pt\displaystyle\hfill#\hfil$\cr<\crcr\sim\crcr}}}{\lower
3.6pt\vbox{\halign{$\mathsurround=0pt\textstyle\hfill#\hfil$\cr<\crcr\sim\crcr}}}{\lower
3.6pt\vbox{\halign{$\mathsurround=0pt\scriptstyle\hfill#\hfil$\cr<\crcr\sim\crcr}}}{\lower
3.6pt\vbox{\halign{$\mathsurround=0pt\scriptscriptstyle\hfill#\hfil$\cr<\crcr\sim\crcr}}}}0.25,0.27,0.19,$
(167)

respectively from the $z=0$ vorticity power spectra convergence studies from
Fig. 3 in Pueblas and Scoccimarro (2009), Fig. 13 in Hahn  _et al._ (2015) and
Fig. 6 in Jelic-Cizmek  _et al._ (2018) and the mass resolutions of their
simulations.

For scale-free simulations, we have a similar situation to CDM, in that the
vorticity power spectrum has a peak at small scales. It peaks at $k_{\rm
peak}/k_{\text{nl}}\simeq 5,4,3.5,2.75$ for $n_{s}=-1,0,1,2$ respectively.
Table 7 shows the resulting values of $\lambda\,k_{\rm peak}$ for the last
four outputs of our $512^{3}$ simulations. We see comparing these values to
the simple convergence criterion in Eq. (167) that we would expect convergence
to be achieved only in the case of $n_{s}=-1$ simulations for the last two or
perhaps three outputs. The next best situation is for $n_{s}=0$ last output.

Figure 23 shows that these estimates are quite reasonable, perhaps on the
conservative side. The left panel shows the vorticity power spectrum for
$n_{s}=-1$ for different outputs. Each line shows a single output (averaged
over two pair-fixed realizations), starting from early outputs (highest
amplitude at low-$k$) down to the last output. We see that the last three to
four outputs show a reasonably converged answer to a self-similar result. As
expected from CDM simulations, and results in Tab. 7, earlier outputs
increasingly overestimate the vorticity power spectrum amplitude, although its
shape is fairly robust. The right panel shows the same convergence study for
$n=0$. As expected, the convergence is less solid, although the last two
outputs are fairly close. The situation is less convincing for $n_{s}=1,2$
(not shown), as expected from the values quoted in Tab. 7, thus we won’t
present any vorticity results for these spectra in this work.

As a result of this analysis, we proceed as follows to build the vorticity
power spectrum that are compared to our predictions in the main text. For
$n_{s}=-1$, we compute the mean and scatter over the last four outputs (8
pair-fixed measurements). For $n_{s}=0$, we do the same for the last two
outputs, with the caveat that the result may be not be fully converged. Still,
since the shape of the power spectra are robust to resolution, it provides
useful information in comparison to our theoretical predictions.

For the velocity divergence power spectra the situation is less delicate. The
large-scale part of the spectrum, which is what tests our perturbative
predictions is rather insensitive to mass resolution; there is more
sensitivity to mass resolution at small scales ($k\mathrel{\mathchoice{\lower
3.6pt\vbox{\halign{$\mathsurround=0pt\displaystyle\hfill#\hfil$\cr>\crcr\sim\crcr}}}{\lower
3.6pt\vbox{\halign{$\mathsurround=0pt\textstyle\hfill#\hfil$\cr>\crcr\sim\crcr}}}{\lower
3.6pt\vbox{\halign{$\mathsurround=0pt\scriptstyle\hfill#\hfil$\cr>\crcr\sim\crcr}}}{\lower
3.6pt\vbox{\halign{$\mathsurround=0pt\scriptscriptstyle\hfill#\hfil$\cr>\crcr\sim\crcr}}}}k_{\text{nl}}$)
in agreement with previous results in the literature Pueblas and Scoccimarro
(2009); Hahn  _et al._ (2015); Jelic-Cizmek  _et al._ (2018). Therefore, in
this case we combine all outputs into a self-similar spectrum, with error bars
describing the scatter between two pair-fixed realizations, and deviations
from self-similarity.

Figure 20:  Self-similarity of the density power spectrum for different spectral indices ($n_{s}=-1$ to $n_{s}=2$ from top to bottom at high-$k$). For each spectral index, different colors denote different outputs.  Figure 21:  Self-similarity of the halo mass function for different spectral indices ($n_{s}=-1$ to $n_{s}=2$ from top to bottom at low $\hat{m}=m/m_{*}$). For each spectral index, different colors denote different outputs.  Table 6: Mass function best fit parameters (see Eq. 120) as a function of spectral index. $n_{s}$ | $A$ | $a$ | $p$
---|---|---|---