# Fast computation of non-linear power spectrum in cosmologies with massive
neutrinos

Hernán E. Noriega  Alejandro Aviles,11footnotetext: Corresponding author.
Sebastien Fromenteau  Mariana Vargas-Magaña

###### Abstract

We compute 1-loop corrections to the redshift space galaxy power spectrum in
cosmologies containing additional scales, and hence kernels different from
Einstein-de Sitter (EdS). Specifically, our method is tailored for cosmologies
in the presence of massive neutrinos and some modified gravity models; in this
article we concentrate on the former case. The perturbative kernels have
contributions that we notice appear either from the logarithmic growth rate
$f(k,t)$, which is scale-dependent because of the neutrino free-streaming, or
from the failure of the commonly used approximation $f^{2}=\Omega_{m}$. The
latter contributions make the computation of loop corrections quite slow,
precluding full-shape analyses for parameter estimation. However, we identify
that the dominant pieces of the kernels come from the growth factor, allowing
us to simplify the kernels but retaining the characteristic free-streaming
scale introduced by the neutrinos’ mass. Moreover, with this simplification
one can exploit FFTLog methods to speed up the computations even more. We
validate our analytical modeling and numerical method with halo catalogs
extracted from the Quijote simulations finding good agreement with the, a
priori, known cosmological parameters. We make public our Python code
FOLPS$\nu$ to compute the redshift space power spectrum in a fraction of
second. Code available at https://github.com/henoriega/FOLPS-nu.

##  1 Introduction

The neutrinos are perhaps the most elusive particles in the Universe, they
couple very weakly to the rest of the standard model particles, and they
always interact or are created in definite flavor states. However, the
observed oscillations of atmospheric and solar neutrinos indicate that their
energy and flavor states are not the same, but they mix, being this only
possible if they have non-zero mass [1]. These observations also put lower
bounds on the sum of their masses, being $0.06\,\text{eV}$ for the normal
hierarchy and $0.11\,\text{eV}$ for the inverted one [2]. On the other hand,
from the energy spectrum of beta electrons emitted in the decay of tritium
into helium-3, the KATRIN experiment has put an upper bound of
$0.8\,\text{eV}$ at 0.90 confidence level (c.l.) in the mass of the lightest
energy state [3, 4], meaning an upper bound of $2.4\,\text{eV}$ in the sum of
the three masses. Hence, neutrinos have excelling small masses. Despite this,
they are so abundant that they contribute considerably to the cosmic energy
budget. They also have specific, well known signatures on the large scale
properties of the Universe which are highly dependent on their masses; see [5,
6] for reviews on neutrino physics in Cosmology. Primeval neutrinos decoupled
early on from the cosmic plasma while they were still relativistic, as they
continued to be until well inside the matter dominated epoch when they have
slowed down sufficiently to become non-relativistic. Since then, primordial
neutrinos have behaved as an additional dark matter component. However, they
are a non-cold component since they still have large velocity dispersion,
preventing them from clustering below a certain scale, named the free-
streaming scale. Since neutrinos do participate in the clustering at large
scale, they leave a well understood signature in the power spectrum,
suppressing it above the free-streaming wave-number $k_{\text{FS}}$ [7, 8].
For realistic masses, this is located at about $0.1\,\,h^{-1}\text{Mpc}$
nowadays, that is, at the onset of the non-linear scales.

Since gravity responds to the total energy momentum tensor, cosmological
probes are particularly sensitive to the sum of the neutrino masses, while the
mass splitting and hierarchy can be safely neglected for current and near
future experiments [9]. Cosmic Microwave Background (CMB) temperature and
polarization anisotropies from Planck 2018 legacy data put stringent
constraints in the neutrino mass of $\sum_{i}m_{\nu,i}\lesssim
0.26\,\text{eV}$ (0.95 c.l.), which is reduced to $0.12\,\text{eV}$ when
combined with DR12 BOSS, MGS and 6dFGS Baryon Acoustic Oscillations (BAO) data
and CMB-lensing [10]; and even more, down to $0.09\,\text{eV}$ when the SDSS
eBOSS DR14 Lyman-alpha forest flux power spectrum data is added [11]. Other
recent analysis as minimal extensions to the base $\Lambda$CDM have yielded
similar bounds, e.g [12, 13, 14, 15, 16].222No minimal extensions as dynamical
dark energy of modified gravity can lead to very different constraints on the
neutrino masses; e.g. [14, 17] That is, the cosmological bounds on the
neutrino mass absolute scale are very impressive and excel those from particle
physics experiments by about one order of magnitude. Clearly, this success
comes with the cost of being indirect measurements and hence very sensitive to
our ability to estimate the other cosmological parameters, particularly the
Hubble rate $H_{0}$ and the amplitude of fluctuations $\sigma_{8}$. Future
galaxy surveys such as DESI333desi.lbl.gov [18, 19],
EUCLID444sci.esa.int/web/euclid [20] and the LSST at the Vera C. Rubin
Observatory555lsst.org/ [21] will be deeper in redshift and wider in subtended
angle, so that it is expected that the neutrino masses will be measured with
unprecedented precision in the near future. This is mainly because of the
reduction of statistical errors in the data. Accordingly, methods to reduce
the systematic errors in the modeling, well below the statistical errors, are
very valuable. This is the subject of this paper.

Recently, the methods of Effective Field Theory (EFT) for large scale
structure formation [22, 23, 24, 25, 26] built on Perturbation Theory (PT)
[27] have reached sufficient maturity to be used to compare the galaxy power
spectrum (as well as other statistics) of real survey data directly to the
analytical modeling of non-linearities [28, 29]; see also [30, 31, 32, 33, 34,
35, 36, 37, 38, 39, 40]. These full-shape analyses capture at the same time
the BAO and the redshift space distortions (RSD) features, reaching similar
constrictive power as Planck experiment for some cosmological
parameters.666Meanwhile, more traditional methods that compress cosmological
information are in continue development, and recent implementations have
similar success than full-shape analyses [41, 42]. Furthermore, with the aid
of FFTLog methods [43], the computation of loop corrections has been speed up
[44, 45, 46, 47, 48] in order to use efficient Markov Chain Monte Carlo (MCMC)
algorithms to draw the posterior distribution of parameters. All the existing
methods have used Einstein-de Sitter (EdS) kernels so far, which is an
excellent approximation as long as the neutrino mass is small ($\sim
0.1\,\text{eV}$), as expected from Planck data. However, the use of
uninformative priors together with EdS kernels can introduce systematic errors
when performing full-shape analyses of data, since one has to explore up to
high mass values. Perhaps the main obstacle, and objection, against the use of
more proper kernels for modeling cosmologies that include massive neutrinos is
the computational time, which for full kernels is prohibitive for parameter
estimation purposes. The objective of this work is to show that under simple
approximations, it is possible to retain the main characteristics of the
proper kernels without sacrificing valuable and expensive computational time.

Before entering into the more technical details in the following sections,
here we present a schematic, but general form of our framework and describe
how we overcome the difficulties.

There are by now, many PT in the presence of massive neutrinos with different
complexity levels and types of approximations when treating non-linearities
[49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61], in this work we adopt
the one developed in [62, 63, 64, 65]. The velocity and density fields of the
combined cold dark matter and baryons fluid $cb$ in the presence of massive
neutrinos are non-local related even at large scales by [63]

$\theta^{(1)}_{cb}(\bm{k})=\frac{f(k)}{f_{0}}\delta^{(1)}_{cb}(\bm{k})$ (1.1)

with $f_{0}(t)=f(k\rightarrow 0,t)$ the cosmic growth rate at very large
scales, or simply the growth rate $f=f_{0}$ in models without massive
neutrinos. Strictly, the above expression holds to linear order in
perturbation theory. However, this property is inherited to higher orders due
to the advection of fields, enforcing the appearance of the free-streaming
scale not only through the linear power spectrum but also via the non-linear
kernels.

The redshift space power spectrum can be written as [64]

$P(k,\mu)=\sum_{m=0}^{\infty}\sum_{n=0}^{m}\mu^{2n}f_{0}^{m}I^{m}_{n}(k)$
(1.2)

where $\mu=\hat{\bm{k}}\cdot\hat{\bm{n}}$ is the cosine angle between the
wave-vector $\bm{k}$ and the line-of-sight direction $\hat{\bm{n}}$. Each
function $I^{m}_{n}$ is obtained via a “loop” integral over internal momenta
$\bm{p}$,

$I(k)=\int d^{3}p\,\mathcal{I}(\bm{k},\bm{p})$ (1.3)

of a function $\mathcal{I}(\bm{p},\bm{k})$. These $\mathcal{I}$ functions are
invariant against rotations, and so they only depend on the length of the
sides of the triangle formed with the momenta $\bm{p}$, $\bm{k}$ and
$\bm{k}-\bm{p}$. In the absence of massive neutrinos, the functions
$\mathcal{I}$ are known analytically, up to the linear power spectra, and the
1-loop integrals in (1.3) can be performed with high precision and very
rapidly using FFTLog methods, since one can write

$\mathcal{I}(\bm{k},\bm{p})=\sum_{a,b,c}\,C_{abc}\,p^{z_{a}}|\bm{k}-\bm{p}|^{z_{b}}k^{z_{c}}$
(1.4)

for a collection of complex powers $z_{a}$, $z_{b}$ and $z_{c}$ and
coefficients $C_{abc}$.

However, in the presence of massive neutrinos the functions
$\mathcal{I}(\bm{k},\bm{p})$ have no analytical form, but instead should be
obtained from a set of differential equations at each point of the
3-dimensional grid in (1.3), that is, at each configuration $k$, $p$ and
$|\bm{k}-\bm{p}|$. This, of course, makes the computation of the loop
integrals quite slow, precluding efficient Markov Chain exploration in order
to draw cosmological parameter constraints via a full-shape analysis.

The perturbative kernels developed in [63, 65] differ from EdS by two kind of
contributions: The first one is due to the growth rate discussed above. The
second is pure non-linear coming from the failure of the approximation
$f^{2}=\Omega_{m}$ that should be obtained by solving differential equations
via Green functions in a similar manner than in $\Lambda$CDM, with the
difference that in our case these corrections are also scale-dependent and so,
the solution to these equations should be obtained at each configuration of
internal and external wave-vectors in loop integrals, and not only one time as
in the massless neutrino case. It turns out that the former contribution,
introducing the free-streaming scale, is much more important than the latter
[65], that plays a role similar to the normalization factors that exist
between $\Lambda$CDM and EdS kernels. If we keep only the growth rates in the
kernels, there is no necessity of solving differential equations and one can
write the functions $\mathcal{I}$ as in eq. (1.4). This approach has two
important advantages over the use of full kernels for efficient evaluation:
there is no need to solve differential equations and one can use FFTLog
methods.

In this paper we present our analytical modeling based on the above mentioned
approximation, and a numerical method with its implementation into the
FOLPS$\nu$ Python code.777Publicly available at
https://github.com/henoriega/FOLPS-nu. We validate the formalism by comparing
against the redshift space power spectrum multipoles of halos obtained from
the Quijote simulations suite
[66].888https://github.com/franciscovillaescusa/Quijote-simulations. We obtain
very satisfactory results, recovering the standard cosmological parameters and
the neutrino mass. We do some unrealistic simplifications in our analyses
since the reach of this paper is to validate our methodology: first we use
halos, instead of galaxies, which together with the use of 100 simulation
boxes each with volume $(1\,h^{-1}\text{Gpc})^{3}$, allow us to have a small
covariance and so, any mismatch in recovering the cosmological parameters of
the simulations is probably more due to errors on our perturbative approach
and not due to statistical errors, or the imperfection on the galaxy-halo
connection modeling that would be present if we use e.g. HODs instead of
halos. Second, by using a large value for the sum of the neutrinos’ three mass
states, $\sum_{i}m_{\nu,i}=0.4\,\text{eV}$, which correspond to the largest
mass available in the Quijote suite, we test the theory in the regimes where
the use of non-EdS kernels are more relevant. Although, we also test against
the mass neutrinos’ case with $\sum_{i}m_{\nu,i}=0.1\,\text{eV}$.

The rest of the paper is organized as follows: In §2, we introduce the
Neutrino free-streaming scale, our key approximations and the linear theory.
§3 reviews the PT theory developed in [63, 64, 65] and introduces the fk-
kernels. In §4 we present the FTTLog method we use which is complemented in
appendix B. The validation of our model is presented in §5 where we test our
modeling to the Quijote suite of simulations; complementary plots are shown in
appendix C. In §6 we give a brief review and some useful information about our
code FOLPS$\nu$. Finally, in §7 we present our conclusions.

##  2 Neutrino free-streaming scale and linear evolution

At very early times, primordial neutrinos are coupled to the primeval plasma,
with a common temperature $T_{\nu}=T$. At this stage, neutrinos are
relativistic following a Fermi-Dirac distribution. As the Universe expands, it
cools down as $T\propto 1/a$ until the weak interaction rate $\Gamma_{\nu}\sim
G^{2}_{\text{F}}T^{5}$ falls below the expansion rate $H\sim
G^{1/2}_{\text{N}}T^{2}$, where $G_{\text{F}}$ and $G_{\text{N}}$ are the
Fermi and Newton constants, and neutrinos decouple from the cosmic plasma when
$\Gamma_{\nu}\sim H$. The decoupling temperature is $T_{\text{dec}}\sim
1\,\text{MeV}$. Since particle physics experiments put a limit on the neutrino
masses $m_{\nu}\lesssim\mathcal{O}(1\,\text{eV})\ll T_{\text{dec}}$, the
neutrinos are still relativistic at decoupling. After that, the neutrinos do
not interact anymore and follow the same Fermi-Dirac distribution, with
$T_{\nu}\propto 1/a$ no longer a “true” temperature, but a parameter of a
distribution that has been freezed out. Neutrinos become non-relativistic when
the mean energy per particle falls below their mass $m_{\nu,\text{i}}$. This
occurs at redshift $1+z_{\text{nr, i}}\approx
1890\,m_{\nu,\text{i}}\,\text{eV}^{-1}$ for each species. Since then, the
neutrinos behave as dark matter, and nowadays their abundance is

$\Omega_{\nu}=\frac{M_{\nu}}{93.14\,h^{2}\,\text{eV}},\qquad\qquad
M_{\nu}\equiv\sum_{\text{i}}m_{\nu,\text{i}}.$ (2.1)

where we defined capital “$M_{\nu}$” as the sum of the masses of the three
species. Despite massive neutrinos become non-relativistic at very early
times, their velocity dispersion is large [67, 5]

$\sigma_{\nu}(z)\simeq\sqrt{\frac{15\zeta(5)}{\zeta(3)}}\frac{T_{\nu}}{m_{\nu}}\simeq
1811\,(1+z)\left(\frac{0.1\,\text{eV}}{m_{\nu}}\right)\,\text{km/s}.$ (2.2)

This velocity dispersion acts as a pressure support counteracting the action
of gravity. As in a Jeans-like mechanism, the neutrinos cannot be confined and
contribute to the formation of structure below a certain scale, named the
free-streaming scale. This corresponds to a wave-number [67]

$k_{\text{FS}}(z)=\sqrt{\frac{4\pi
G\bar{\rho}a^{2}}{\frac{5}{9}\sigma^{2}_{\nu}}}\simeq
0.0908\frac{H(z)/H_{0}}{(1+z)^{2}}\left(\frac{m_{\nu}}{0.1\,\text{eV}}\right)\,h\,\text{Mpc}^{-1}.$
(2.3)

This expression gives us an idea of the involved scales, but is never used in
this work. In fact, there are several definitions for the free-streaming scale
in the literature, all of them having similar numerical values. Notice that
the above equation is valid when neutrinos are non-relativistic, while before
that, their free-streaming was equal to the Hubble horizon. Hence, the free-
streaming has a minimum at the non-relativistic transition redshift, and all
clustering on length scales below it becomes suppressed by the presence of
massive neutrinos [5]. For realistic neutrino masses this scale is about
$1\,h^{-1}\text{Gpc}$, much larger than $1/k_{\text{FS}}(z)$ at recent times.
Further, eq. (2.3) is valid for a single neutrino species, corresponding to
the three different masses. From now on, we will consider equal mass
neutrinos, so $M_{\nu}=3\times m_{\nu,i}$ and there is only one single free-
streaming scale; although, our method does not depend on this choice. Three
equal masses is also commonly assumed in $N$-body simulations. Below the free-
streaming scale $k\ll k_{\text{FS}}$ neutrinos behave as cold dark matter
contributing to the formation of structures, while above it $k\gg
k_{\text{FS}}$ they do not cluster anymore, producing a suppression of the
power spectrum and correlation function for scales below than the free-
streaming.

At late times, inside the matter dominated epoch, the gravitational potential
is related to the matter fields through the Poisson equation

$-\frac{k^{2}}{a^{2}}\Phi(\bm{k},t)=\frac{3}{2}\Omega_{m}(t)H^{2}(t)\big{(}f_{cb}\delta_{cb}(\bm{k},t)+f_{\nu}\delta_{\nu}(\bm{k},t)\big{)},$
(2.4)

where $f_{X}=\Omega_{X}/\Omega_{m}$ is the relative abundance of the matter
component $X$ to the total matter $m$. Since the only clustering components
are assumed to be neutrinos, cold dark matter and baryons, then
$f_{\nu}=1-f_{cb}$. We are primarily interested in the evolution of the $cb$
field, since to a very good approximation galaxies are tracers of it, and not
of the total matter field [68, 69, 70, 71, 72]. Hence, our approach takes the
non-linear neutrino density field as an intermediate quantity that we
approximate as [65]

$\delta_{\nu}(\bm{k},t)=\frac{T_{\nu}(k,t)}{T_{cb}(k,t)}\delta_{cb}(\bm{k},t),$
(2.5)

with $T_{\nu,cb}$ the linear transfer functions of the neutrinos and the $cb$
field. That is, we do not assume that the neutrino field is linear, which
would violate momentum conservation [56], but instead that it is proportional
to the non-linear $cb$ density field. Our approximation was shown to be very
accurate for the total matter power spectrum up to non-linear scales using
simulations (see figure 3 of [65]). The same, or very similar, approximations
were considered in some works [53, 54, 58], but their perturbative theories
differ than our in almost any other aspect.

Using eq. (2.5), the Poisson equation becomes

$-\frac{k^{2}}{a^{2}}\Phi(\bm{k},t)=A(k,t)\delta_{cb}(\bm{k},t)$ (2.6)

with

$A(k,t)\equiv\frac{3}{2}\Omega_{m}(t)H^{2}(t)\left(f_{cb}+f_{\nu}\frac{T_{\nu}(k,t)}{T_{cb}(k,t)}\right).$
(2.7)

At large scales, the ratio of transfer functions is unity and $A\rightarrow
4\pi G\bar{\rho}_{m}$, while the neutrino transfer function falls sharply
about the free-streaming scale and $A\rightarrow 4\pi G\bar{\rho}_{cb}$ at
small scales.

Under our approximation, the $cb$ field evolves according to the fluid
equations,

$\displaystyle\partial_{t}\delta_{cb}(\bm{x},t)+\frac{1}{a}\partial_{i}\big{[}(1+\delta_{cb})v^{i}_{cb}\big{]}$
$\displaystyle=0,$ (2.8)
$\displaystyle\partial_{t}v^{i}_{cb}(\bm{x},t)+\frac{1}{a}v^{j}_{cb}\partial_{j}v^{i}_{cb}+Hv^{i}_{cb}+\frac{1}{a}\partial^{i}\Phi$
$\displaystyle=0,$ (2.9)

where $\bm{v}_{cb}$ is the peculiar velocity of the $cb$ field. Linearizing
the above equations for a longitudinal velocity field, we can solve for the
linear order density field

$\delta_{cb}^{(1)}(k,t)=D_{+}(k,t)\delta_{cb}^{(1)}(k,t_{0}),$ (2.10)

where $D_{+}$ is the scale-dependent growth function given by

$\frac{d^{2}D_{+}(k,t)}{dt^{2}}+2H\frac{dD_{+}(k,t)}{dt}-A(k,t)D_{+}(k,t)=0$
(2.11)

with the initial conditions given in [8] to pick up the growing solution. We
choose the normalization $D_{+}(k\rightarrow 0,t_{0})=1$ with $t_{0}$ the
present time. Notice that although the decomposition of eq. (2.10) looks
arbitrary, with this normalization it becomes unique and coincides with the
massless neutrinos case at large scales.

Figure 1: Growth rate $f(k)$ divided by their large scale limits
$f_{0}=f(k\rightarrow 0)\equiv f^{M_{\nu}=0}$ evaluated at redshift $z=0.5$,
for two cosmologies including massive neutrinos with $M_{\nu}=0.1\,\text{eV}$
(blue lines) and $M_{\nu}=0.4\,\text{eV}$ (red lines). We use eqs. (2.11) and
(2.12) (dotted lines) which we call it “exact”, and the algebraic
approximation of Hu & Eisenstein (1998) given in [7] (solid lines). The
horizontal dot-dashed gray lines show the small scales limit.

The linear (logarithmic) growth rate is given by

$f(k,t)=\frac{d\ln D_{+}(k,t)}{d\ln a(t)},$ (2.12)

which is again scale-dependent. We define its value at very large scales as
[73]

$f_{0}(t)\equiv f(k\rightarrow 0,t)=f^{M_{\nu}=0}(t),$ (2.13)

which coincides with the scale-independent growth rate $f(t)$ of a cosmology
with the same amount of total matter, but on which the neutrinos are massless.
In figure 1 we show a plot for $f(k,t)$ computed directly using eqs. (2.11)
and (2.12) in dotted lines and using the approximate formulae given in [7] in
solid lines. The discrepancy among both methods is smaller than the $0.1\%$
even for the largest mass of $M_{\nu}=0.4\,\text{eV}$. The gray horizontal
lines show the limit of small scales, $f(k\gg
k_{\text{FS}})\rightarrow(1-\frac{3}{5}f_{\nu})f_{0}$, which corresponds to
the maximum suppression, that is, when the contributions from neutrinos is
completely negligible. The Hu-Eisenstein approximation is used in our code
because otherwise one has to compute the neutrinos and $cb$ transfer functions
for each wave-number and for each cosmology, which is very time-consuming.

As is common, we define the longitudinal piece of the velocity field as

$\theta_{cb}(\bm{x},t)=-\frac{\partial_{i}v^{i}_{cb}}{aHf_{0}}(\bm{x},t).$
(2.14)

With this, we obtain the relation (1.1) between linear velocity and density
fields.

As a recapitulation, our method relies on the approximation for the neutrino
density field given by eq. (2.5). Further, we are assuming that the neutrinos
late times velocity dispersion tensor has no impact on the gravitational
potentials of the metric, and so they become equal $\Psi=\Phi$. It is known
that this is no precise and actually the velocity dispersion leads to the very
definition of free-streaming scale [74] and the generation of vorticity [75].
But, by considering eq. (2.5), we are actually restricting our theory to the
$cb$ field in the presence of massive neutrinos, and we do not deal directly
with the non-linear evolution of the neutrinos themselves. This can be seen as
a disadvantage, since one cannot construct a consistent non-linear total
matter power spectrum, and in this respect other approaches, particularly the
one of [60], are more complete. But as we mentioned, the galaxy field is a
biased tracer of the $cb$ fluid and not of the total matter fluid, and this is
the object of our primary interest in this work. However, for some
applications one deals with the total real space matter power spectrum, as in
weak lensing. In such a case, the approximation of eq. (2.5) implies

$P_{m}(k,t)=P_{cb}(k,t)\left[(1-f_{\nu})^{2}+2(1-f_{\nu})f_{\nu}\frac{T_{\nu}(k,t)}{T_{cb}(k,t)}+f_{\nu}^{2}\left(\frac{T_{\nu}(k,t)}{T_{cb}(k,t)}\right)^{2}\right],$
(2.15)

which still gives very good results [63].

##  3 Perturbative model in massive neutrino cosmologies

The peculiar velocities of galaxies $\bm{v}$ produce a Doppler effect along
the line-of-sight direction $\hat{\bm{n}}$ (we assume the distant observer
approximation where we fix $\hat{\bm{n}}$), breaking the real space isotropy,
such that the power spectrum is given by [76, 77]

$(2\pi)^{3}\delta_{\text{D}}(\bm{k})+P_{s}(\bm{k})=\int
d^{3}x\,e^{-i\bm{k}\cdot\bm{x}}\Big{[}1+\mathcal{M}(\bm{J}=\bm{k},\bm{x})\Big{]},$
(3.1)

in terms of the density-weighted velocity moments generating function

$\mathcal{M}(\bm{J},\bm{x})=\left\langle\big{(}1+\delta(\bm{x}_{1})\big{)}\big{(}1+\delta(\bm{x}_{2})\big{)}e^{-i\bm{J}\cdot\Delta\bm{u}}\right\rangle,$
(3.2)

with $\bm{x}=\bm{x}_{2}-\bm{x}_{1}$ and
$\Delta\bm{u}=\bm{u}(\bm{x}_{2})-\bm{u}(\bm{x}_{1})$. The “velocity” $\bm{u}$
and the peculiar velocities are related by

$\bm{u}(\bm{x})=\hat{\bm{n}}\frac{\bm{v}(\bm{x})\cdot\hat{\bm{n}}}{aH}.$ (3.3)

The $m$-th density weighted velocity field moment of the generation function
is a $m$-rank tensor defined as

$\displaystyle\Xi^{m}_{i_{1}\cdots i_{m}}(\bm{x})$ $\displaystyle\equiv
i^{m}\frac{\partial^{m}}{\partial J_{i_{1}}\cdots\partial
J_{i_{m}}}\big{[}1+\mathcal{M}(\bm{J},\bm{x})\big{]}\Big{|}_{\bm{J}=0}=\langle\big{(}1+\delta_{1}\big{)}\big{(}1+\delta_{2}\big{)}\Delta
u_{i_{1}}\cdots\Delta u_{i_{m}}\rangle,$ (3.4)

with $\delta_{1}=\delta(\bm{x}_{1})$ and $\delta_{2}=\delta(\bm{x}_{2})$.

By expanding in this way eq. (3.1), we get the power spectrum at the moment
expansion approach,

$\displaystyle(2\pi)^{3}\delta_{\text{D}}(\bm{k})+P_{s}(\bm{k})$
$\displaystyle=\sum_{m=0}^{\infty}\frac{(-i)^{m}}{m!}k_{i_{1}}\dots
k_{i_{m}}\tilde{\Xi}_{i_{1}\cdots i_{m}}^{m}(\bm{k}),$ (3.5)

where the $\tilde{\Xi}^{m}_{i_{1}\cdots i_{m}}(\bm{k})$ are the Fourier
moments of the generating function —the Fourier transforms of their
configuration space counterparts, $\Xi^{m}_{i_{1}\cdots i_{n}}(\bm{x})$:

$\tilde{\Xi}_{i_{1}\cdots i_{m}}^{m}(\bm{k})=\int
d^{3}x\,e^{-i\bm{k}\cdot\bm{x}}\,\Xi^{m}_{i_{1}\cdots i_{m}}(\bm{x}).$ (3.6)

For rotational invariant $S(\bm{k},\bm{p})$ function one has [78, 64]

$\displaystyle\int\frac{d^{3}p}{(2\pi)^{3}}(\hat{\bm{p}}\cdot\hat{\bm{n}})^{n}S(\bm{k},\bm{p})$
$\displaystyle=\sum_{m=0}^{n}\mu^{m}\int\frac{d^{3}p}{(2\pi)^{3}}G_{nm}(\hat{\bm{k}}\cdot\hat{\bm{p}})S(\bm{k},\bm{p})$
$\displaystyle=\frac{k^{3}}{4\pi^{2}}\sum_{m=0}^{n}\mu^{m}\int_{0}^{\infty}dr\,r^{2}\int_{-1}^{1}dx\,G_{nm}(x)S(k,r,x),$
(3.7)

with $x=\hat{\bm{k}}\cdot\hat{\bm{p}}$ and $r=p/k$, and

$\displaystyle G_{nm}(x)$
$\displaystyle=\sum_{\ell=0}^{n}\frac{(1+(-1)^{\ell+n})(2\ell+1)}{2(1+\ell+n)}\binom{\ell}{m}\binom{2\ell}{\ell}\binom{\frac{\ell+m-1}{2}}{\ell}$
$\displaystyle\quad\times\,{}_{3}F_{2}\left(\frac{1-\ell}{2},-\frac{\ell}{2},\frac{1}{2}(-1-\ell-n);\frac{1}{2}-\ell,\frac{1}{2}(1-\ell-n);1\right)\mathcal{L}_{\ell}(x),$
(3.8)

with $\mathcal{L}_{\ell}$ the Legendre polynomial of degree $\ell$, and
${}_{3}F_{2}(a_{1},a_{2},a_{3};b_{1},b_{2};z)$ the extended hypergeometric
function of the type $p=3$ and $q=2$.

Using eq.(3), and the properties of $n$-point correlators of homogeneous and
isotropic fields it can be shown that

$\frac{(-i)^{m}}{m!}k_{i_{1}}\dots k_{i_{m}}\tilde{\Xi}_{i_{1}\cdots
i_{m}}^{m}(\bm{k})=\sum_{n=0}^{m}\mu^{2n}f_{0}^{m}I^{m}_{n}(k),$ (3.9)

for only momentum magnitude dependent functions $I^{m}_{n}(k)$, recovering
eq.(1.2). The case that we will use, up to 1-loop corrections, has been shown
in [64].

###  3.1 Perturbation Theory with full kernels

As was mentioned in §2, there is some evidence that galaxies and other cosmic
objects are tracers of the combined dark matter and baryonic fields, and not
of the total matter field that also includes the neutrinos. The underlying
reason is not completely clear but can be understood from peak theory since
primeval density overdensities, from which halos will arise, are formed early
on, before the neutrinos are non-relativistic; and moreover, these proto-halos
do not capture the neutrinos because of the free-streaming. Up to mildly non-
linear scales, the combined CDM and baryonic fields behave essentially as a
single fluid, named $cb$. The scale at which this stops to occur is also not
clear and is highly dependent on how the baryonic physics are modeled [79,
80]; however, the consensus is that at large scales, those reached by
perturbation theory and EFTofLSS, baryonic effects can be neglected to a very
good approximation [81, 79, 80]. Therefore, we are interested in the loop
corrections of the $cb$ power spectrum, and on top of that we model the
tracers with a flexible biasing scheme.

There are several perturbation theory schemes in the literature, we follow
here the one developed recently in [63, 65]. The first step was to compute the
Lagrangian Perturbation Theory (LPT) kernels $L^{(n)}_{i}$, which are easier
to obtain than their Eulerian, Standard Perturbation Theory (SPT) counterparts
when additional scales are involved. Thereafter, these kernels can be mapped
to the Eulerian frame, to obtain the $F_{n}$ and $G_{n}$. In the following, we
summarize the results of [63, 65].

The $cb$ fluid elements follow trajectories given by

$\bm{x}=\bm{q}+\bm{\Psi}(\bm{q},t),$ (3.10)

where $\bm{q}$ is the initial position and $\bm{\Psi}$ is the Lagrangian
displacement, which is obtained through the geodesic equation

$\ddot{\bm{\Psi}}(\bm{q},t)+2H\dot{\bm{\Psi}}(\bm{q},t)=-\nabla\Phi(\bm{x},t)\big{|}_{\bm{x}=\bm{q}+\bm{\Psi}(\bm{q},t)}.$
(3.11)

The Lagrangian displacement is expanded in a Taylor series in Fourier space,
$\bm{\Psi}=\sum_{n}\bm{\Psi}^{(n)}$. At order $n$999For simplicity, we adopt
the notations  $\underset{\bm{k}_{1\cdots
n}=\bm{k}}{\int}=\int\frac{d^{3}\bm{k}_{1}}{(2\pi)^{3}}\cdots\frac{d^{3}\bm{k}_{n}}{(2\pi)^{3}}(2\pi)^{3}\delta_{\text{D}}(\bm{k}_{1\cdots
n}-\bm{k})\quad\mathrm{and}\quad\bm{k}_{1\cdots
n}=\bm{k}_{1}+\cdots+\bm{k}_{n}.$ (3.12)

$\Psi_{i}^{(n)}(\bm{k},t)=\frac{i}{n!}\underset{\bm{k}_{1\cdots
n}=\bm{k}}{\int}L^{(n)}_{i}(\bm{k}_{1},\cdots,\bm{k}_{n};t)\delta_{cb}^{(1)}(\bm{k}_{1},t)\cdots\delta_{cb}^{(1)}(\bm{k}_{n},t).$
(3.13)

In [63] the $L^{(n)}_{i}$ kernels are obtained, and we reproduce them here.
The linear order kernel is

$L_{i}^{(1)}(\bm{k})=\frac{k_{i}}{k^{2}},$ (3.14)

which is identical when neutrinos are massless. To second and third orders,

$\displaystyle L^{(2)}_{i}(\bm{k}_{1},\bm{k}_{2},t)$
$\displaystyle=\frac{3}{7}\frac{k^{i}_{12}}{k^{2}_{12}}\left(\mathcal{A}-\mathcal{B}\,\frac{(\bm{k}_{1}\cdot\bm{k}_{2})^{2}}{k_{1}^{2}k_{2}^{2}}\right),$
(3.15) $\displaystyle L^{(3)}_{i}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3},t)$
$\displaystyle=\frac{k^{i}_{123}}{k^{2}_{123}}\Bigg{\\{}\frac{5}{7}\left(\mathcal{A}^{(3)}-\mathcal{B}^{(3)}\frac{(\bm{k}_{2}\cdot\bm{k}_{3})^{2}}{k^{2}_{2}k^{3}_{2}}\right)\left(1-\frac{(\bm{k}_{1}\cdot\bm{k}_{23})^{2}}{k_{1}^{2}k_{23}^{2}}\right)$
$\displaystyle-\frac{1}{3}\left(\,\mathcal{C}^{(3)}-3\mathcal{D}^{(3)}\frac{(\bm{k}_{2}\cdot\bm{k}_{3})^{2}}{k^{2}_{2}k^{2}_{3}}+2\mathcal{E}^{(3)}\frac{(\bm{k}_{1}\cdot\bm{k}_{2})(\bm{k}_{2}\cdot\bm{k}_{3})(\bm{k}_{3}\cdot\bm{k}_{1})}{k_{1}^{2}k^{2}_{2}k^{2}_{3}}\,\right)\Bigg{\\}}.$
(3.16)

Functions $\mathcal{A}$ and $\mathcal{B}$ in the second order kernels depend
on the wave-vectors $\bm{k}_{1}$ and $\bm{k}_{2}$ and on time $t$, while all
the functions in the third order Lagrangian kernel
$\mathcal{A}^{(3)},\dots,\mathcal{E}^{(3)}$ depend on the same arguments as
$L^{(3)}$, that is on $\bm{k}_{1}$, $\bm{k}_{2}$ and $\bm{k}_{3}$ and on time
$t$. In the case of $\Lambda$CDM with massless neutrinos all these functions
are only time dependent and reduce to unity for EdS evolution. These functions
completely determine the kernels and should be obtained by solving a system of
linear differential equations at each wave-vector configuration and for each
set of cosmological parameters $\Omega_{cb}$, $\Omega_{\nu}$ and $h$. The
expressions for them are not displayed here, not only because they are large
but also because they are not used by the method presented in this work, they
can be found in [63].

We use the functions $C^{(n)}\Gamma_{n}(\bm{k}_{1},\cdots,\bm{k}_{n},t)$ and
$C^{(n)}\Gamma_{n}^{f}(\bm{k}_{1},\cdots,\bm{k}_{n},t)$ as the kernels of the
fields $-i\bm{k}\cdot\bm{\Psi}^{(n)}$ and
$-i(nHf_{0})^{-1}\bm{k}\cdot\dot{\bm{\Psi}}^{(n)}$, respectively [73]. The
numbers $C^{(n)}$ are chosen for algebraical convenience and have the values
$C^{(1)}=C^{(3)}=1$ and $C^{(2)}=3/7$. As usual, we assume the Lagrangian
displacement is longitudinal at early times, and as such the only relevant
kernels for 1-loop, 2-point statistics are the $\Gamma_{n}$ functions. In
terms of the standard Lagrangian kernels $L^{(n)}$, one can write

$\displaystyle C^{(n)}\Gamma_{n}(\bm{k}_{1},\cdots,\bm{k}_{n},t)$
$\displaystyle=k^{i}_{1\cdots n}L_{i}^{(n)}(\bm{k}_{1},\cdots,\bm{k}_{n},t).$
(3.17)

and

$\Gamma^{f}_{n}(\bm{k}_{1},\cdots,\bm{k}_{n},t)=\Gamma_{n}(\bm{k}_{1},\cdots,\bm{k}_{n},t)\frac{f(k_{1})+\cdots+f(k_{n})}{nf_{0}}+\frac{1}{nf_{0}H}\dot{\Gamma}_{n}(\bm{k}_{1},\cdots,\bm{k}_{n},t).$
(3.18)

For the case of EdS evolution, $\Gamma_{n}=\Gamma_{n}^{f}$, and one obtains
the widely used relation $\dot{\bm{\Psi}}^{(n)}=nHf\bm{\Psi}^{(n)}$ [82, 83].

The $cb$ density fluctuation becomes [65]

$\displaystyle\delta_{cb}(\bm{k})$
$\displaystyle=\sum_{m=1}^{\infty}\frac{(-i)^{m}}{m!}k_{i_{1}}\cdots
k_{i_{m}}\underset{\bm{p}_{1\cdots
m}=\bm{k}}{\int}\Psi_{i_{1}}(\bm{p}_{1})\cdots\Psi_{i_{m}}(\bm{p}_{m}).$
(3.19)

which is obtained from the well-known relation $\delta(\bm{k})=\int
d^{3}q\;e^{-i\bm{k}\cdot\bm{q}}\Big{[}e^{-i\bm{k}\cdot\bm{\Psi}(\bm{q})}-1\Big{]}$.
Similarly, considering the peculiar velocity $\bm{v}_{cb}=a\dot{\bm{\Psi}}$,
the $\theta_{cb}$ field is obtained through [65, 84]

$\displaystyle\theta_{cb}(\bm{k})$
$\displaystyle=-\frac{1}{Hf_{0}}\sum_{m=0}^{\infty}\frac{(-i)^{m}}{m!}\int
d^{3}q\;e^{-i\bm{k}\cdot\bm{q}}\big{(}\bm{k}\cdot\bm{\Psi}\big{)}^{m}$
$\displaystyle\quad\times\Big{[}\,\dot{\Psi}_{i,i}+(\delta_{ij}\delta_{ab}-\delta_{ia}\delta_{jb})\Psi_{a,b}\dot{\Psi}_{i,j}+\frac{1}{2}\epsilon_{ikp}\epsilon_{jqr}\Psi_{k,q}\Psi_{p,r}\dot{\Psi}_{i,j}\Big{]},$
(3.20)

with the Lagrangian displacements evaluated at the Lagrangian coordinate
$\bm{q}$, and $\epsilon_{ijk}$ is the Levi-Civita symbol and $\delta_{ij}$ the
Kronocker delta.

On the other hand, the kernels $F_{n}$ and $G_{n}$ are defined via the
following equations for the $n$-th order density fluctuation and velocity

$\displaystyle\delta_{cb}^{(n)}(\bm{k},t)$
$\displaystyle=\underset{\bm{k}_{1\cdots
n}=\bm{k}}{\int}F_{n}(\bm{k}_{1},\cdots,\bm{k}_{n};t)\delta_{cb}^{(1)}(\bm{k}_{1},t)\cdots\delta_{cb}^{(1)}(\bm{k}_{n},t),$
(3.21) $\displaystyle\theta_{cb}^{(n)}(\bm{k},t)$
$\displaystyle=\underset{\bm{k}_{1\cdots
n}=\bm{k}}{\int}G_{n}(\bm{k}_{1},\cdots,\bm{k}_{n};t)\delta_{cb}^{(1)}(\bm{k}_{1},t)\cdots\delta_{cb}^{(1)}(\bm{k}_{n},t).$
(3.22)

Comparing order by order eq. (3.21) with (3.19) and eq. (3.22) with (3.1) one
obtains the Eulerian kernels in terms of the Lagrangian ones. For our
purposes, we are interested in the specific configurations of wave-vectors
that enter in 1-loop integrals of 2-point statistics: for second order
kernels, these are $\bm{k}_{1}=\bm{p}$ and $\bm{k}_{2}=\bm{k}-\bm{p}$, while
for third order $\bm{k}_{1}=\bm{k}$, $\bm{k}_{2}=-\bm{p}$ and
$\bm{k}_{3}=\bm{p}$. We obtain

$F_{1}(k)=1,\qquad G_{1}(k)=\frac{f(k)}{f_{0}},$ (3.23) $\displaystyle
F_{2}(\bm{p},\bm{k}-\bm{p})$
$\displaystyle=\frac{3}{14}\Gamma_{2}(\bm{p},\bm{k}-\bm{p})+\frac{1}{2}\frac{(\bm{k}\cdot\bm{p})(\bm{k}\cdot(\bm{k}-\bm{p}))}{p^{2}|\bm{k}-\bm{p}|^{2}},$
(3.24) $\displaystyle G_{2}(\bm{p},\bm{k}-\bm{p})$
$\displaystyle=\frac{3}{7}\Gamma^{f}_{2}(\bm{p},\bm{k}-\bm{p})+\frac{(\bm{p}\cdot(\bm{k}-\bm{p}))^{2}}{p^{2}|\bm{k}-\bm{p}|^{2}}\frac{f(|\bm{k}-\bm{p}|)+f(p)}{2f_{0}}$
$\displaystyle\quad+\frac{1}{2}\frac{\bm{p}\cdot(\bm{k}-\bm{p})}{p\,|\bm{k}-\bm{p}|}\left(\frac{|\bm{k}-\bm{p}|}{p}\frac{f(|\bm{k}-\bm{p}|)}{f_{0}}+\frac{p}{|\bm{k}-\bm{p}|}\frac{f(p)}{f_{0}}\right),$
(3.25)

and,

$\displaystyle F_{3}(\bm{k},-\bm{p},\bm{p})$
$\displaystyle=\frac{1}{6}\Gamma_{3}(\bm{k},-\bm{p},\bm{p})+\frac{1}{7}\frac{\bm{k}\cdot(\bm{k}-\bm{p})\,\bm{k}\cdot\bm{p}}{p^{2}|\bm{k}-\bm{p}|^{2}}\Gamma_{2}(\bm{k},-\bm{p})-\frac{1}{6}\frac{(\bm{k}\cdot\bm{p})^{2}}{p^{4}},$
(3.26) $\displaystyle G_{3}(\bm{k},-\bm{p},\bm{p})$
$\displaystyle=\frac{1}{2}\Gamma^{f}_{3}(\bm{k},-\bm{p},\bm{p})+\frac{2}{7}\frac{\bm{k}\cdot\bm{p}}{p^{2}}\Gamma^{f}_{2}(\bm{k},-\bm{p})+\frac{1}{7}\frac{f(p)}{f_{0}}\Gamma_{2}(\bm{k},-\bm{p})\frac{\bm{k}\cdot(\bm{k}-\bm{p})}{|\bm{k}-\bm{p}|^{2}}$
$\displaystyle\quad-\frac{1}{7}\left[2\Gamma_{2}^{f}(\bm{k},-\bm{p})+\Gamma_{2}(\bm{k},-\bm{p})\frac{f(p)}{f_{0}}\right]\left[1-\frac{(\bm{p}\cdot(\bm{k}-\bm{p}))^{2}}{p^{2}|\bm{k}-\bm{p}|^{2}}\right]$
$\displaystyle\quad-\frac{1}{6}\frac{(\bm{k}\cdot\bm{p})^{2}}{p^{4}}\frac{f(k)}{f_{0}}.$
(3.27)

where the $\Gamma_{3}$ kernels used in the above equations are constructed by
first symmetrizing $L^{(3)}_{i}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$ in eq.
(3.16) over its three arguments and thereafter evaluating at
$\bm{k}_{1}=\bm{k}$, $\bm{k}_{2}=-\bm{p}=-\bm{k}_{3}$. For $\Omega_{cb}(a)=1$,
one has $f(k)=f_{0}$ and all the scale- and time-dependent functions
$\mathcal{A},...$ are unity, reducing the above expressions to the standard
EdS kernels.

The integrals for the loop corrections have exactly the same form as in the
standard massless neutrino case, but with the use of the corrected kernels.
The $\delta\delta$, $\delta\theta$ and $\theta\theta$ power spectra are

$P^{\text{1-loop}}_{ab}(k)=P^{L}_{ab}(k)+P_{ab}^{22}(k)+P_{ab}^{13}(k),$
(3.28)

where $a$ and $b$ refer to $\delta$ or $\theta$ fields, and linear power
spectra $P^{L}_{ab}(k)$,

$P^{L}_{\delta\delta}(k)\equiv P_{L}(k),\quad
P^{L}_{\delta\theta}(k)=\frac{f(k)}{f_{0}}P_{L}(k),\quad
P^{L}_{\theta\theta}(k)=\left(\frac{f(k)}{f_{0}}\right)^{2}P_{L}(k),$ (3.29)

with $P_{L}$ the linear $cb$ real space power spectrum, and leading non-linear
contributions given by101010Hereafter we use the notation
$\int_{\bm{p}}=\int\frac{d^{3}p}{(2\pi)^{3}}.$ (3.30)

$\displaystyle P_{\delta\delta}^{22}(k)$
$\displaystyle=2\int_{\bm{p}}\big{[}F_{2}(\bm{p},\bm{k}-\bm{p})\big{]}^{2}P_{L}(p)P_{L}(|\bm{k}-\bm{p}|),$
(3.31) $\displaystyle P_{\delta\theta}^{22}(k)$
$\displaystyle=2\int_{\bm{p}}F_{2}(\bm{p},\bm{k}-\bm{p})G_{2}(\bm{p},\bm{k}-\bm{p})P_{L}(p)P_{L}(|\bm{k}-\bm{p}|),$
(3.32) $\displaystyle P_{\theta\theta}^{22}(k)$
$\displaystyle=2\int_{\bm{p}}\big{[}G_{2}(\bm{p},\bm{k}-\bm{p})\big{]}^{2}P_{L}(p)P_{L}(|\bm{k}-\bm{p}|),$
(3.33) $\displaystyle P_{\delta\delta}^{13}(k)$
$\displaystyle=6P_{L}(k)\int_{\bm{p}}F_{3}(\bm{k},-\bm{p},\bm{p})P_{L}(p),$
(3.34) $\displaystyle P_{\delta\theta}^{13}(k)$
$\displaystyle=3P_{L}(k)\int_{\bm{p}}\big{[}F_{3}(\bm{k},-\bm{p},\bm{p})\frac{f(k)}{f_{0}}+G_{3}(\bm{k},-\bm{p},\bm{p})\big{]}P_{L}(p),$
(3.35) $\displaystyle P_{\theta\theta}^{13}(k)$
$\displaystyle=6\frac{f(k)}{f_{0}}P_{L}(k)\int_{\bm{p}}G_{3}(\bm{k},-\bm{p},\bm{p})P_{L}(p).$
(3.36)

These integrals can be performed straightforwardly with standard methods in
perturbation theory. However, the computational time becomes quite large
because one has to solve a set of differential equations for the functions
$\mathcal{A},\mathcal{B},\mathcal{A}^{(3)},\mathcal{B}^{(3)}$ at each volume
element in the 2-dimensional quadratures, since these do not have analytical
expressions for arbitrary momenta, making the computations considerably slower
than in the massless neutrinos case. This makes the use of these kernels not
viable for parameters estimation, where the integrals must be evaluated
thousands of times.

We will come back later in §3.3 to the rest of the necessary functions
entering eq. (1.2) and the biasing scheme. But before that we are showing the
approximation we use to overcome the difficulties explained in the above
paragraph.

###  3.2 fk-kernels

The fk-kernels method consists on fixing the functions $\mathcal{A},\cdots$
equal to their $\Lambda$CDM or EdS values, but keeping the linear growth rates
$f(k)$. More precisely, one constructs the fk-kernels by taking the
$\Gamma_{n}$ kernels [eq. (3.17)] as their EdS or $\Lambda$CDM counterparts,
which are known analytically (we use EdS in this paper) or computed by solving
a system of differential equations only once (for $\Lambda$CDM),111111For
$\Lambda$CDM, $\mathcal{A}=\mathcal{B}$, $\mathcal{A}^{(3)}=\mathcal{B}^{(3)}$
and $\mathcal{C}^{(3)}=\mathcal{D}^{(3)}=\mathcal{E}^{(3)}$, and are all
functions of time only. Thus, they can be solved just once and thereafter one
performs the loop integrals. and in the $\Gamma_{n}^{f}$ kernels of eq. (3.18)
one keeps the $f(k)$ rates. Thereafter, one uses these $\Gamma_{n}$ into the
SPT kernels $F_{n}$ and $G_{n}$ expressions of eqs. (3.24)-(3.27).

In this way, the perturbative kernels inherit the free-streaming scale
imprinted by the large neutrinos’ velocity dispersion, such that also the
advection of fields is well modeled at non-linear orders. For example, the
difference of linear velocity fields
$\theta^{(1)}(\bm{x}+\bm{\Psi})-\theta^{(1)}(\bm{x})$ is a second order effect
that yields the contribution

$\underset{\bm{k}_{12}=\bm{k}}{\int}\frac{\hat{\bm{k}}_{1}\cdot\hat{\bm{k}}_{2}}{2}\left(\frac{k_{1}}{k_{2}}\frac{f(k_{1})}{f_{0}}+\frac{k_{2}}{k_{1}}\frac{f(k_{2})}{f_{0}}\right)\delta_{cb}^{(1)}(\bm{k}_{1})\delta_{cb}^{(1)}(\bm{k}_{2})\in\theta^{(2)}(\bm{k}),$
(3.37)

to the velocity field, where the rates $f(k)/f_{0}$ appear because of the non-
trivial linear relation between density and velocity fields in eq. (1.1).

In ref. [65] it is shown that this approximation is quite accurate and the
growth rates constitute the most important additional pieces for kernels of
the $cb$ fields; see figure 9 of that paper.

A second advantage of the use of fk-kernels is that it is possible to write
them as the sum of power laws

$p^{2n_{1}}|\bm{k}-\bm{p}|^{2n_{2}}k^{-2(n_{1}+n_{2})}$ (3.38)

multiplied by growth rates evaluated either at wave-numbers $p$,
$|\bm{k}-\bm{p}|$ or $k$. This allows us to use FFTLog methods [43] to speed
up even more the calculations [44, 45, 46, 47, 48].

For example, the second order kernels are
$F_{2}^{\text{{fk}}}=F_{2}^{\text{EdS}}$ and

$\displaystyle G_{2}^{\text{{fk}}}(\bm{p},\bm{k}-\bm{p})$
$\displaystyle=\left(\frac{k^{4}}{14|\bm{k}-\bm{p}|^{2}p^{2}}+\frac{3k^{2}}{28|\bm{k}-\bm{p}|^{2}}-\frac{k^{2}}{7p^{2}}-\frac{5p^{2}}{28|\bm{k}-\bm{p}|^{2}}+\frac{|\bm{k}-\bm{p}|^{2}}{14p^{2}}+\frac{3}{28}\right)\frac{f(p)}{f_{0}}$
$\displaystyle\quad+\big{(}\,p\leftrightarrow|\bm{k}-\bm{p}|\,\big{)},$ (3.39)

and clearly for $f(p)=f(|\bm{k}-\bm{p}|)=f_{0}$, one has
$G_{2}^{\text{{fk}}}(\bm{p},\bm{k}-\bm{p})=G_{2}^{\text{EdS}}$. We will come
back to the numerical method to speed up the loop integrals in §4

###  3.3 EFT galaxy power spectrum multipoles

We use the EFT galaxy power spectrum as given by [85, 76, 86, 64, 65, 87, 37,
77]

$\displaystyle P^{\text{EFT}}_{s}(k,\mu)$
$\displaystyle=P_{\delta\delta}(k)+2f_{0}\mu^{2}P_{\delta\theta}(k)+f_{0}^{2}\mu^{4}P_{\theta\theta}(k)+A^{\text{TNS}}(k,\mu)+D(k,\mu)$
$\displaystyle\quad+(\alpha_{0}+\alpha_{2}\mu^{2}+\alpha_{4}\mu^{4}+\alpha_{6}\mu^{6})k^{2}P_{L}(k)+\tilde{c}(f_{0}\sigma_{v}k\mu)^{4}P_{s}^{K}(k,\mu)$
$\displaystyle\quad+P_{\text{shot}}\big{[}\alpha^{shot}_{0}+\alpha^{shot}_{2}(k\mu)^{2}\big{]}.$
(3.40)

The first line of the above equation is obtained by using the density-weighted
velocity momentum expansion (3.5) up to 1-loop in perturbation theory:
$P_{\delta\delta}$, $P_{\delta\theta}$ and $P_{\theta\theta}$ are the tracers
1-loop real space power spectra for the velocity and density fields, as given
below in §3.3.1. The functions $A^{\text{TNS}}(k,\mu)$ and $D(k,\mu)$ are [64,
86]

$\displaystyle A^{\text{TNS}}(k,\mu)$ $\displaystyle=2k\mu
f_{0}\int_{\bm{p}}\frac{\bm{p}\cdot\hat{\bm{n}}}{p^{2}}B_{\sigma}(\bm{p},-\bm{k},\bm{k}-\bm{p})\,,$
(3.41) $\displaystyle D(k,\mu)$ $\displaystyle=(k\mu
f_{0})^{2}\int_{\bm{p}}\Big{\\{}F(\bm{p})F(\bm{k}-\bm{p})$
$\displaystyle\quad+\frac{(\bm{p}\cdot\hat{\bm{n}})^{2}}{p^{4}}P_{\theta\theta}^{L}(p)\big{[}P^{K}_{s}(|\bm{k}-\bm{p}|,\mu_{\bm{k}-\bm{p}})-P^{K}_{s}(k,\mu)\big{]}\Big{\\}},$
(3.42)

with $\mu_{\bm{k}-\bm{p}}$ the cosine angle between the wave-vector
$\bm{k}-\bm{p}$ and the line-of-sight direction $\hat{\bm{n}}$. Here, the
bispectrum $B_{\sigma}$ and function $F$ are

$\displaystyle(2\pi)^{3}\delta_{\text{D}}(\bm{k}_{1}+\bm{k}_{2}+\bm{k}_{3})B_{\sigma}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})=$
$\displaystyle\qquad\qquad\Big{\langle}\theta(\bm{k}_{1})\left[\delta(\bm{k}_{2})+f_{0}\frac{(\bm{k}_{2}\cdot\hat{\bm{n}})^{2}}{k_{2}^{2}}\theta(\bm{k}_{2})\right]\left[\delta(\bm{k}_{3})+f_{0}\frac{(\bm{k}_{3}\cdot\hat{\bm{n}})^{2}}{k_{3}^{2}}\theta(\bm{k}_{3})\right]\Big{\rangle},$
(3.43)

and

$\displaystyle F(\bm{p})$
$\displaystyle=\frac{\bm{p}\cdot\hat{\bm{n}}}{p^{2}}\Big{[}P_{\delta\theta}(p)+f_{0}\frac{(\bm{p}\cdot\hat{\bm{n}})^{2}}{p^{2}}P_{\theta\theta}(p)\Big{]}.$
(3.44)

While

$P^{K}_{s}(k,\mu)=\big{(}1+\mu^{2}f(k)\big{)}^{2}P_{L}(k)$ (3.45)

is the linear Kaiser power spectrum [85]. In Appendix A we will write the
functions $A^{\text{TNS}}(k,\mu)$ and $D(k,\mu)$ in the form of eq. (1.2)
which is more suitable for numerical analyses.

In the second line of eq. (3.3), we add EFT parameters $\alpha_{0}$,
$\alpha_{2}$, $\alpha_{4}$ and $\alpha_{6}$ that model the backreaction of
small scales over large scales and the non-linear map between real and
redshift spaces. Further, along the line-of-sight direction, 2-point
statistics are dominated by Fingers of God as a non-linear coupling between
the velocity and density fields, with a characteristic scale given by the
velocity dispersion $\sigma_{v}$, yielding the next-to-leading order
counterterm $\tilde{c}$ [29].

The third line in eq. (3.3) includes the noise parameters, uncorrelated with
long wave-length fluctuations. $P_{\text{shot}}$ is considered a constant that
can be set equal to the Poisson process shot noise
$P_{\text{Poisson}}=1/\bar{n}_{X}$ with $\bar{n}_{X}$ the number density of
tracers, or to any other constant whose value is not relevant since it is
completely degenerate with $\alpha_{0}^{shot}$. We have added a tilt
proportional to $(k\mu)^{2}$ [88, 89, 87, 90, 36, 91, 37]. The shot noise
departs from a pure white noise (when $\alpha^{shot}_{2}\neq 0$) because
stochasticity is not localized at a single point [88, 89] and because of the
stochastic nature of peculiar velocities at small scales [90].

Regarding the set of EFT and noise parameters introduced above, not all are
necessary at the same time. As we already said, the value of $P_{\text{shot}}$
is not relevant and only alters the numerical values of $\alpha^{shot}_{0}$
and $\alpha^{shot}_{2}$, so it can take any value, e.g. the Poissonian shot
noise. On the other hand, one usually leads with only the first two or three
non-vanishing multipoles of the power spectrum, hence not all the EFT
parameters are necessary: for example, below we will compute only the monopole
and the quadrupole, thus the counterterms $\alpha_{4}$ and $\alpha_{6}$ would
be redundant and not considered. Finally, the introduced functional dependence
by $\tilde{c}$ is $k^{4}P_{L}(k)$, which is approximately $\propto k^{2}$ at
high-$k$, and hence degenerate with $\alpha_{2}^{shot}$. And so, the common is
to choose between the use of either $\tilde{c}$ or $\alpha_{2}^{shot}$. In
this work, we use the latter as a nuisance parameter.

####  3.3.1 Biased tracers

We adopt the bias prescription of McDonald [92, 88, 93], which has been
tweaked to work for theories with additional scales [63]. Notice the
introduction of the neutrino mass impedes to have a complete bias expansion at
any given PT order; however, curvature and higher order biases can be used to
partially tame this inconsistency [89]. Ultimately, we obtain good numerical
results with the same number of bias parameters as in $\Lambda$CDM. That is,
we consider the bias parameters $b_{1}$, $b_{2}$, $b_{s^{2}}$ and $b_{3nl}$.
We do not add the curvature bias $b_{\nabla^{2}\delta}$, because it is
absorbed by the Effective Field Theory parameters (introduced below) since
they are degenerate. The isotropic power spectra for velocity and density
fields become

$\displaystyle P_{\delta\delta}(k)$
$\displaystyle=b_{1}^{2}P^{\text{1-loop}}_{\delta\delta}(k)+2b_{1}b_{2}P_{b_{1}b_{2}}(k)+2b_{1}b_{s^{2}}P_{b_{1}b_{s^{2}}}(k)+b_{2}^{2}P_{b_{2}^{2}}(k)$
$\displaystyle\quad+2b_{2}b_{s^{2}}P_{b_{2}b_{s^{2}}}(k)+b_{s^{2}}^{2}P_{b_{s^{2}}^{2}}(k)+2b_{1}b_{3nl}\sigma^{2}_{3}(k)P_{L}(k),$
(3.46) $\displaystyle P_{\delta\theta}(k)$
$\displaystyle=b_{1}P^{\text{1-loop}}_{\delta\theta}(k)+b_{2}P_{b_{2},\theta}(k)+b_{s^{2}}P_{b_{s^{2}},\theta}(k)+b_{3nl}\sigma^{2}_{3}(k)P^{L}_{cb,\delta\theta}(k),$
(3.47) $\displaystyle P_{\theta\theta}(k)$
$\displaystyle=P^{\text{1-loop}}_{\theta\theta}(k),$ (3.48)

where the quantities $P^{\text{1-loop}}_{\delta\delta}(k)$,
$P^{\text{1-loop}}_{\delta\theta}(k)$ and
$P^{\text{1-loop}}_{\theta\theta}(k)$ represent the 1-loop power spectra
obtained from eqs. (3.28) and (3.31)–(3.36). The other contributions are
defined as [93]

$\displaystyle P_{b_{1}b_{2}}(k)$
$\displaystyle=\int_{\bm{p}}F_{2}(\bm{p},\bm{k}-\bm{p})P_{L}(p)P_{L}(|\bm{k}-\bm{p}|),$
(3.49) $\displaystyle P_{b_{1}b_{s^{2}}}(k)$
$\displaystyle=\int_{\bm{p}}F_{2}(\bm{p},\bm{k}-\bm{p})S_{2}(\bm{p},\bm{k}-\bm{p})P_{L}(p)P_{L}(|\bm{k}-\bm{p}|),$
(3.50) $\displaystyle P_{b_{2}^{2}}(k)$
$\displaystyle=\frac{1}{2}\int_{\bm{p}}P_{L}(p)\big{[}P_{L}(|\bm{k}-\bm{p}|)-P_{L}(p)\big{]},$
(3.51) $\displaystyle P_{b_{2}b_{s^{2}}}(k)$
$\displaystyle=\frac{1}{2}\int_{\bm{p}}P_{L}(p)\left[P_{L}(|\bm{k}-\bm{p}|)S_{2}(\bm{p},\bm{k}-\bm{p})-\frac{2}{3}P_{L}(p)\right],$
(3.52) $\displaystyle P_{b_{s^{2}}^{2}}(k)$
$\displaystyle=\frac{1}{2}\int_{\bm{p}}P_{L}(p)\left[P_{L}(|\bm{k}-\bm{p}|)[S_{2}(\bm{p},\bm{k}-\bm{p})]^{2}-\frac{4}{9}P_{L}(p)\right],$
(3.53) $\displaystyle P_{b_{2},\theta}(k)$
$\displaystyle=\int_{\bm{p}}G_{2}(\bm{p},\bm{k}-\bm{p})P_{L}(p)P_{L}(|\bm{k}-\bm{p}|),$
(3.54) $\displaystyle P_{b_{s^{2}},\theta}(k)$
$\displaystyle=\int_{\bm{p}}G_{2}(\bm{p},\bm{k}-\bm{p})S_{2}(\bm{p},\bm{k}-\bm{p})P_{L}(p)P_{L}(|\bm{k}-\bm{p}|),$
(3.55)

with

$S_{2}(\bm{k}_{1},\bm{k}_{2})=\frac{(\bm{k}_{1}\cdot\bm{k}_{2})^{2}}{k_{1}^{2}k_{2}^{2}}-\frac{1}{3}.$
(3.56)

Additionally, the function $\sigma^{2}_{3}(k)$ is given by

$\sigma^{2}_{3}(k)=\frac{105}{16}\int_{\bm{p}}P_{L}(p)\left[S_{2}(\bm{p},\bm{k}-\bm{p})\left(\frac{2}{7}S_{2}(-\bm{p},\bm{k})-\frac{4}{21}\right)+\frac{8}{63}\right].$
(3.57)

Alternatively, one can use the bias expansion with $b_{\mathcal{G}_{2}}$ and
$b_{\Gamma_{3}}$ [94], instead of $b_{s^{2}}$ and $b_{3nl}$ [88].

####  3.3.2 IR-resummations

Coherent flows of the matter field stream over a scale settled by the variance
of Lagrangian displacements $\sigma^{2}_{\Psi}$. Because the scale is
comparable in size to the BAO peak width, overdense regions in the Universe
become partially depleted, while underdense regions partially populated,
broadening the acoustic peak [95, 96, 97]. This effect is well described by
LPT, since bulk flows are captured by the displacement fields even at linear
order [82, 83, 25, 90, 37]. In contrast, in SPT the convergence is very slow
[98], and therefore non-perturbative methods, called IR-resummations, must be
used [99]. We follow the prescription of [100, 101, 32], that splits the
linear power spectrum in a piece that does not contain the BAO (the non-wiggle
power spectrum, $P_{nw}$) and a wiggle piece $P_{w}$, such that the real space
linear power spectrum can be written as $P_{L}=P_{nw}+P_{w}$. As a result of
this splitting, the 1-loop IR-resummed EFT redshift space power spectrum
becomes [101]

$\displaystyle P_{s}^{\text{IR}}(k,\mu)$
$\displaystyle=e^{-k^{2}\Sigma^{2}_{\text{tot}}(k,\mu)}P_{s}^{\text{EFT}}(k,\mu)+\big{(}1-e^{-k^{2}\Sigma^{2}_{\text{tot}}(k,\mu)}\big{)}P_{s,nw}^{\text{EFT}}(k,\mu)$
$\displaystyle\quad+e^{-k^{2}\Sigma^{2}_{\text{tot}}(k,\mu)}P_{w}(k)k^{2}\Sigma^{2}_{\text{tot}}(k,\mu),$
(3.58)

where the wiggle piece $P_{s}^{\text{EFT}}(k,\mu)$ is the 1-loop power
spectrum computed using eq. (3.3). The non-wiggle part,
$P_{s,nw}^{\text{EFT}}(k,\mu)$ is also computed with eq. (3.3) but using as
input the non-wiggle linear power spectrum $P_{nw}$. The function
$\Sigma^{2}_{\text{tot}}$ is given by

$\Sigma^{2}_{\text{tot}}(k,\mu)=\big{[}1+f(k)\mu^{2}\big{(}2+f(k)\big{)}\big{]}\Sigma^{2}+f^{2}(k)\mu^{2}(\mu^{2}-1)\delta\Sigma^{2},$
(3.59)

with

$\displaystyle\Sigma^{2}$
$\displaystyle=\frac{1}{6\pi^{2}}\int_{0}^{k_{s}}dp\,P_{nw}(p)\left[1-j_{0}\left(p\,\ell_{\text{BAO}}\right)+2j_{2}\left(p\,\ell_{\text{BAO}}\right)\right],$
(3.60) $\displaystyle\delta\Sigma^{2}$
$\displaystyle=\frac{1}{2\pi^{2}}\int_{0}^{k_{s}}dp\,P_{nw}(p)j_{2}\left(p\,\ell_{\text{BAO}}\right),$
(3.61)

where $\ell_{\text{BAO}}\simeq 105\,h^{-1}\text{Mpc}$ is the BAO peak scale
and $j_{0}$ and $j_{2}$ are the spherical Bessel functions of degree 0 and 2.
The scale $k_{s}$ splits the long and short modes, whose choice is somewhat
arbitrary. We use the value $k_{s}=0.4\,h\,\text{Mpc}^{-1}$.

####  3.3.3 Multipoles

The IR-resummed EFT power spectrum of eq. (3.3.2) is the one we compare the
simulations to. More precisely, we take its monopole, quadrupole and
hexadecapole multipoles from

$P_{\ell}(k)=\frac{2\ell+1}{2}\int_{-1}^{1}d\mu\;P_{s}^{\text{IR}}(k,\mu)\mathcal{L}_{\ell}(\mu),$
(3.62)

where $\mathcal{L}_{\ell}$ are the Legendre polynomial of degree $\ell$.

##  4 FFTLog formalism

FFTLog algorithms give us the 1-dimensional Fast Fourier Transform or Fast
Hankel Transform of a function evaluated over a set of log-spaced intervals
[102, 43]. Its usefulness in cosmology arise since one has to manipulate
functions with relevant features spanning several orders of magnitude
[43].121212https://jila.colorado.edu/$\sim$ajsh/FFTLog/ For example, it serves
to obtain the real space matter 2-point correlation function from the matter
power spectrum faster, with higher accuracy, and notably reduced ringing and
aliasing than when using direct integration. More recently, these methods have
been reintroduced to help to bring the problem of integrating loop expressions
onto the problem of matrix multiplication, which can be handled very fast with
the use of proper numerical methods.

The starting point is to decompose the linear real space matter power spectrum
in to a series of power laws

$\bar{P}_{L}(k)=\sum_{m=-N/2}^{N/2}c_{m}\,k^{\nu+i\eta_{m}},$ (4.1)

where $N$ is the number of sampling points logarithmic spaced over the
interval $[k_{\text{min}},\dots,k_{\text{max}}]$. The Fourier coefficients
$c_{m}$ and exponents $\eta_{m}$ are given by

$\displaystyle c_{m}$
$\displaystyle=\frac{W_{m}}{N}k_{\text{min}}^{-(\nu+i\eta_{m})}\sum_{l=0}^{N-1}P_{L}(k_{l})\left(\frac{k_{l}}{k_{\text{min}}}\right)^{-\nu}e^{-2\pi
iml/N},$ (4.2) $\displaystyle\eta_{m}$ $\displaystyle=\frac{N-1}{N}\frac{2\pi
m}{\log(k_{\text{max}}/k_{\text{min}})},$ (4.3)

with $W_{m}=1$ except $W_{\pm N/2}=1/2$. The parameter $\nu$ is an arbitrary
real number known as the FFTLog bias. It serves as a leverage to obtain better
convergence of the loop integrals. That is, the $c_{m}$ are not the FFT of the
set $P_{L}(k_{m})$, but of $P_{L}(k_{m})k_{m}^{-\nu}$. Notice $\bar{P}_{L}$ in
eq. (4.1) denotes the approximation of the linear power spectrum, while the
$P_{L}$ appearing in eq. (4.2) is the real one as obtained from an Einstein-
Boltzmann code. For brevity, we will not continue to use this “bar” notation.

In the presence of massive neutrinos one must distinguish between the linear
power spectra $\delta\delta$, $\delta\theta$ and $\theta\theta$, since these
are not equal due to eq. (1.1). Hence, similar expansions to eq. (4.1) are
constructed,

$\displaystyle P_{\delta\theta}^{L}(k)$
$\displaystyle=\frac{f(k)}{f_{0}}P_{L}(k)=\sum_{m=-N/2}^{N/2}c_{m}^{f}\,k^{\nu+i\eta_{m}},$
(4.4) $\displaystyle P_{\theta\theta}^{L}(k)$
$\displaystyle=\left(\frac{f(k)}{f_{0}}\right)^{2}P_{L}(k)=\sum_{m=-N/2}^{N/2}c_{m}^{ff}\,k^{\nu+i\eta_{m}},$
(4.5)

where the coefficients $c_{m}^{f}$ and $c_{m}^{ff}$ are obtained analogously
to $c_{m}$ in eq. (4.2) but with the substitution of $P_{L}$ by
$P^{L}_{\delta\theta}$ and $P^{L}_{\theta\theta}$, respectively.

###  4.1 Velocity power spectrum

All the $I^{m}_{n}(k)$ functions are either of the “$P_{22}$-type” or
“$P_{13}$-type”. The former type means that one has to integrate over a kernel
that can be written as a power series of $p^{2}$, $k^{2}$ and
$|\bm{k}-\bm{p}|^{2}$, multiplied by growth rates and by one linear power
spectrum evaluated at $p$ and other evaluated at $|\bm{k}-\bm{p}|$. For the
$P_{13}$-type case is similar but with the power spectra evaluated one at the
internal momentum $p$ and the other at external momentum $k$.

In this subsection, we consider the velocity power spectrum
$P^{\text{loop}}_{\theta\theta}(k)$ as an example. The rest of the functions
$I^{m}_{n}$ are computed in the same form, but the algebra is lengthy, so we
show the results in appendix B.

We start with the “$P_{22}$-type” piece given by eq. (3.33), where

$\displaystyle P^{22}_{\theta\theta}(k)$
$\displaystyle=2\int_{\bm{p}}\big{[}G_{2}(\bm{p},\bm{k}-\bm{p})\big{]}^{2}P_{L}(p)P_{L}(|\bm{k}-\bm{p}|)$
$\displaystyle=2\int_{\bm{p}}K^{f_{\bm{k}-\bm{p}}f_{\bm{k}-\bm{p}}}_{\theta\theta}(\bm{p},\bm{k}-\bm{p})\frac{f^{2}(|\bm{k}-\bm{p}|)}{f^{2}_{0}}P_{L}(p)P_{L}(|\bm{k}-\bm{p}|)$
$\displaystyle\quad+2\int_{\bm{p}}K^{f_{\bm{p}}f_{\bm{p}}}_{\theta\theta}(\bm{p},\bm{k}-\bm{p})\frac{f^{2}(p)}{f^{2}_{0}}P_{L}(p)P_{L}(|\bm{k}-\bm{p}|)$
$\displaystyle\quad+2\int_{\bm{p}}K^{f_{\bm{k}-\bm{p}}f_{\bm{p}}}_{\theta\theta}(\bm{p},\bm{k}-\bm{p})\frac{f(|\bm{k}-\bm{p}|)f(p)}{f^{2}_{0}}P_{L}(p)P_{L}(|\bm{k}-\bm{p}|),$
(4.6)

with kernels

$\displaystyle K^{f_{\bm{p}}f_{\bm{p}}}_{\theta\theta}(\bm{p},\bm{k}-\bm{p})$
$\displaystyle=\frac{k^{8}}{196|\bm{k}-\bm{p}|^{4}p^{4}}+\frac{3k^{6}}{196|\bm{k}-\bm{p}|^{4}p^{2}}-\frac{k^{6}}{49|\bm{k}-\bm{p}|^{2}p^{4}}-\frac{3k^{4}}{196|\bm{k}-\bm{p}|^{2}p^{2}}$
$\displaystyle\quad-\frac{11k^{4}}{784|\bm{k}-\bm{p}|^{4}}+\frac{3k^{4}}{98p^{4}}-\frac{15k^{2}p^{2}}{392|\bm{k}-\bm{p}|^{4}}-\frac{k^{2}|\bm{k}-\bm{p}|^{2}}{49p^{4}}+\frac{29k^{2}}{392|\bm{k}-\bm{p}|^{2}}$
$\displaystyle\quad-\frac{3k^{2}}{196p^{2}}-\frac{15p^{2}}{392|\bm{k}-\bm{p}|^{2}}+\frac{25p^{4}}{784|\bm{k}-\bm{p}|^{4}}+\frac{3|\bm{k}-\bm{p}|^{2}}{196p^{2}}+\frac{|\bm{k}-\bm{p}|^{4}}{196p^{4}}$
$\displaystyle\quad-\frac{11}{784},$ (4.7) $\displaystyle
K^{f_{\bm{k}-\bm{p}}f_{\bm{p}}}_{\theta\theta}(\bm{p},\bm{k}-\bm{p})$
$\displaystyle=\frac{k^{8}}{196|\bm{k}-\bm{p}|^{4}p^{4}}-\frac{k^{6}}{196|\bm{k}-\bm{p}|^{4}p^{2}}+\frac{37k^{4}}{784|\bm{k}-\bm{p}|^{2}p^{2}}-\frac{9k^{4}}{196|\bm{k}-\bm{p}|^{4}}$
$\displaystyle\quad+\frac{13k^{2}p^{2}}{196|\bm{k}-\bm{p}|^{4}}-\frac{13k^{2}}{196|\bm{k}-\bm{p}|^{2}}-\frac{9p^{2}}{392|\bm{k}-\bm{p}|^{2}}-\frac{5p^{4}}{196|\bm{k}-\bm{p}|^{4}}+\frac{19}{392}$
$\displaystyle\quad+\big{(}\,p\leftrightarrow|\bm{k}-\bm{p}|\,\big{)},$ (4.8)
$\displaystyle
K^{f_{\bm{k}-\bm{p}}f_{\bm{k}-\bm{p}}}_{\theta\theta}(\bm{p},\bm{k}-\bm{p})$
$\displaystyle=K^{f_{\bm{p}}f_{\bm{p}}}_{\theta\theta}(\bm{k}-\bm{p},\bm{p}).$
(4.9)

In the second equality of eq. (4.1) we have split
$\big{[}G_{2}(\bm{p},\bm{k}-\bm{p})\big{]}^{2}$ in three pieces, each of them
proportional to the growth rate evaluated at different wave-numbers.
Interchanging the integration variable $\bm{p}\rightarrow\bm{k}-\bm{p}$ of the
first of these three terms, we can rewrite eq. (4.1) as

$\displaystyle P^{22}_{\theta\theta}(k)$
$\displaystyle=4\int_{\bm{p}}K^{f_{\bm{p}}f_{\bm{p}}}_{\theta\theta}(\bm{p},\bm{k}-\bm{p})P^{L}_{\theta\theta}(p)P_{L}(|\bm{k}-\bm{p}|)$
$\displaystyle\quad+2\int_{\bm{p}}K^{f_{\bm{k}-\bm{p}}f_{\bm{p}}}_{\theta\theta}(\bm{p},\bm{k}-\bm{p})P^{L}_{\delta\theta}(p)P^{L}_{\delta\theta}(|\bm{k}-\bm{p}|),$
(4.10)

with

$K^{f_{\bm{p}}f_{\bm{p}}}_{\theta\theta}(\bm{p},\bm{k}-\bm{p})\equiv\displaystyle\sum\limits_{n_{1},n_{2}=-2}^{2}f^{f_{\bm{p}}f_{\bm{p}}}_{22,\theta\theta}(n_{1},n_{2})k^{-2(n_{1}+n_{2})}p^{2n_{1}}|\bm{k}-\bm{p}|^{2n_{2}},$
(4.11)
$f^{f_{\bm{p}}f_{\bm{p}}}_{22,\theta\theta}(n_{1},n_{2})=\blockarray{cccccc}\text{\scriptsize$n_{2}=-2$}&\text{\scriptsize$-1$}\text{\scriptsize$0$}\text{\scriptsize$1$}\text{\scriptsize$2$}\\\
\block{(ccccc)c}\frac{1}{196}-\frac{1}{49}\frac{3}{98}-\frac{1}{49}\frac{1}{196}\text{\scriptsize$n_{1}=-2$}\\\
\frac{3}{196}-\frac{3}{196}-\frac{3}{196}\frac{3}{196}0\text{\scriptsize$-1$}\\\
-\frac{11}{784}\frac{29}{392}-\frac{11}{784}00\text{\scriptsize$0$}\\\
-\frac{15}{392}-\frac{15}{392}000\text{\scriptsize$1$}\\\
\frac{25}{784}0000\text{\scriptsize$2$}\\\ ,$ (4.12)

and kernel

$K^{f_{\bm{k}-\bm{p}}f_{\bm{p}}}_{\theta\theta}(\bm{p},\bm{k}-\bm{p})\equiv\displaystyle\sum\limits_{n_{1},n_{2}=-2}^{2}f^{f_{\bm{k}-\bm{p}}f_{\bm{p}}}_{22,\theta\theta}(n_{1},n_{2})k^{-2(n_{1}+n_{2})}p^{2n_{1}}|\bm{k}-\bm{p}|^{2n_{2}},$
(4.13)

with matrix

$f^{f_{\bm{k}-\bm{p}}f_{\bm{p}}}_{22,\theta\theta}(n_{1},n_{2})=\blockarray{cccccc}\text{\scriptsize$n_{2}=-2$}&\text{\scriptsize$-1$}\text{\scriptsize$0$}\text{\scriptsize$1$}\text{\scriptsize$2$}\\\
\block{(ccccc)c}\frac{1}{98}-\frac{1}{196}-\frac{9}{196}\frac{13}{196}-\frac{5}{196}\text{\scriptsize$n_{1}=-2$}\\\
-\frac{1}{196}\frac{37}{392}-\frac{13}{196}-\frac{9}{392}0\text{\scriptsize$-1$}\\\
-\frac{9}{196}-\frac{13}{196}\frac{19}{196}00\text{\scriptsize$0$}\\\
\frac{13}{196}-\frac{9}{392}000\text{\scriptsize$1$}\\\
-\frac{5}{196}0000\text{\scriptsize$2$}\\\ ,$ (4.14)

where we have written the expansions (4.7) and (4.8) compactly since all their
summands have the form $k^{-2(n_{1}+n_{2})}p^{2n_{1}}|\bm{k}-\bm{p}|^{2n_{2}}$
with $n_{1},n_{2}\in\\{-2,-1,0,1,2\\}$.

Using the FFTLog decompositions of eq. (4.1), (4.4) and (4.5), and the above
defined matrices we approximate

$\displaystyle P^{22}_{\delta\theta}(k)$
$\displaystyle=4\sum_{m_{1},m_{2}}c^{ff}_{m_{1}}c_{m_{2}}\sum_{n_{1},n_{2}=-2}^{2}f^{f_{\bm{p}}f_{\bm{p}}}_{22,\theta\theta}(n_{1},n_{2})k^{-2(n_{1}+n_{2})}\int_{\bm{p}}\frac{1}{p^{2\nu_{1}-2n_{1}}|\bm{k}-\bm{p}|^{2\nu_{2}-2n_{2}}}$
$\displaystyle+2\sum_{m_{1},m_{2}}c_{m_{1}}^{f}c_{m_{2}}^{f}\sum_{n_{1},n_{2}=-2}^{2}f^{f_{\bm{k}-\bm{p}}f_{\bm{p}}}_{22,\theta\theta}(n_{1},n_{2})k^{-2(n_{1}+n_{2})}\int_{\bm{p}}\frac{1}{p^{2\nu_{1}-2n_{1}}|\bm{k}-\bm{p}|^{2\nu_{2}-2n_{2}}},$
(4.15)

where the sums over $m_{1}$ and $m_{2}$ run from $-N/2$ to $N/2$, and for
simplicity we defined

$\nu_{1}=-\frac{1}{2}(\nu+i\eta_{m_{1}}),\qquad\nu_{2}=-\frac{1}{2}(\nu+i\eta_{m_{2}}).$
(4.16)

To further simplify the previous expression, we define

$\displaystyle M_{22,\theta\theta}^{f_{\bm{p}}f_{\bm{p}}}(\nu_{1},\nu_{2})$
$\displaystyle=4\sum_{n_{1},n_{2}=-2}^{2}f^{f_{\bm{p}}f_{\bm{p}}}_{22,\theta\theta}(n_{1},n_{2})I(\nu_{1}-n_{1},\nu_{2}-n_{2}),$
(4.17) $\displaystyle
M_{22,\theta\theta}^{f_{\bm{k}-\bm{p}}f_{\bm{p}}}(\nu_{1},\nu_{2})$
$\displaystyle=2\sum_{n_{1},n_{2}=-2}^{2}f^{f_{\bm{k}-\bm{p}}f_{\bm{p}}}_{22,\theta\theta}(n_{1},n_{2})I(\nu_{1}-n_{1},\nu_{2}-n_{2}),$
(4.18)

where the $I(\nu_{1},\nu_{2})$ function is defined as [103, 104]

$\displaystyle I(z_{1},z_{2})$ $\displaystyle\equiv
k^{-3+2z_{12}}\int_{\bm{p}}\frac{1}{p^{2z_{1}}|\bm{k}-\bm{p}|^{2z_{2}}}=\frac{1}{8\pi^{3/2}}\frac{\Gamma\left(\frac{3}{2}-z_{1}\right)\Gamma\left(\frac{3}{2}-z_{2}\right)\Gamma\left(z_{12}-\frac{3}{2}\right)}{\Gamma{(z_{1})}\Gamma{(z_{2})}\Gamma{(3-z_{12})}},$
(4.19)

for complex numbers $z_{1}$, $z_{2}$ and $z_{12}\equiv z_{1}+z_{2}$, and
$\Gamma(z)$ is the Gamma function.

Then, using the “$M$” matrices we can write the approximation for
$P^{22}_{\theta\theta}(k)$ as

$\displaystyle P^{22}_{\theta\theta}(k)$
$\displaystyle=k^{3}\displaystyle\sum\limits_{m_{1},m_{2}}c^{ff}_{m_{1}}k^{-2\nu_{1}}M^{f_{\bm{p}}f_{\bm{p}}}_{22,\theta\theta}(\nu_{1},\nu_{2})\,c_{m_{2}}k^{-2\nu_{2}}$
$\displaystyle\quad+k^{3}\displaystyle\sum\limits_{m_{1},m_{2}}c^{f}_{m_{1}}k^{-2\nu_{1}}M^{f_{\bm{k}-\bm{p}}f_{\bm{p}}}_{22,\theta\theta}(\nu_{1},\nu_{2})\,c^{f}_{m_{2}}k^{-2\nu_{2}},$
(4.20)

with

$\displaystyle M^{f_{\bm{p}}f_{\bm{p}}}_{22,\theta\theta}(\nu_{1},\nu_{2})$
$\displaystyle=\Big{[}98\nu_{1}^{3}\nu_{2}+7\nu_{1}^{2}\big{(}2\nu_{2}(7\nu_{2}-8)+1\big{)}-\nu_{1}\big{(}2\nu_{2}(7\nu_{2}+17)+53\big{)}-12(1-2\nu_{2})^{2}\Big{]}$
$\displaystyle\quad\times\frac{2\nu_{12}-3}{98\nu_{1}(\nu_{1}+1)\nu_{2}(\nu_{2}+1)(2\nu_{2}-1)}I(\nu_{1},\nu_{2}),$
(4.21) $\displaystyle
M^{f_{\bm{k}-\bm{p}}f_{\bm{p}}}_{22,\theta\theta}(\nu_{1},\nu_{2})$
$\displaystyle=\Big{[}7\nu_{1}^{2}(7\nu_{2}+3)+\nu_{1}\big{(}7\nu_{2}(7\nu_{2}-1)-10\big{)}+\nu_{2}(21\nu_{2}-10)-37\Big{]}$
$\displaystyle\quad\times\frac{2\nu_{12}-3}{98\nu_{1}(\nu_{1}+1)\nu_{2}(\nu_{2}+1)}I(\nu_{1},\nu_{2}).$
(4.22)

Before continuing with the rest of the terms in
$P^{\text{loop}}_{\theta\theta}$ let us discuss a little bit on what we have
done. The formalism developed above allows us to transform the evaluation of
the loop integrals into the contraction of cosmology-independent matrices $M$
and cosmology-dependent terms $c_{m}$ (and $c_{m}^{f}$ and $c_{m}^{ff}$ when
massive neutrinos are present), speeding up the loop calculations compared to
the usual, direct integration, because one is able to perform matrix
multiplications with fast algorithms which are very common in standard
software library for numerical linear algebra. Moreover, since the matrices do
not depend on cosmology, they can be pre-computed and stored, while the
cosmology-dependent terms can be found quickly through Fast Fourier Transform
(FFT). Naively, one may think of that for EdS the number of computations were
reduced drastically, about half, since in this case
$c_{m}=c_{m}^{f}=c_{m}^{ff}$. However, this is not correct since some of the
same computations are shared for the different $I^{m}_{n}$ functions and we
take advantage of this fact as much as possible. Hence, we have to perform a
total of 26 $P_{22}$-type matrix multiplications with our method, instead of
24 if EdS is used, and all of the multiplications take about the same amount
of time. For $P_{13}$-type the numbers of multiplications is reduced to 7 from
11. Since most of the time the code is performing $P_{22}$-type matrix
multiplications, the time spent with the use of EdS kernels is only marginally
smaller than with our method.

Finally, we develop the approximation to $P^{13}_{\theta\theta}(k)$,

$\displaystyle P^{13}_{\theta\theta}(k)$
$\displaystyle=6\frac{f(k)}{f_{0}}P_{L}(k)\int_{\bm{p}}G_{3}(\bm{k},-\bm{p},\bm{p})P_{L}(p)$
$\displaystyle=6\left(\frac{f(k)}{f_{0}}\right)^{2}P_{L}(k)\int_{\bm{p}}K^{f_{\bm{k}}}_{\theta\theta}(\bm{k},-\bm{p},\bm{p})P_{L}(p)+6\frac{f(k)}{f_{0}}P_{L}(k)\int_{\bm{p}}K^{f_{\bm{p}}}_{\theta\theta}(\bm{k},-\bm{p},\bm{p})P^{L}_{\delta\theta}(p).$
(4.23)

Observe that we have decomposed $G_{3}(\bm{k},-\bm{p},\bm{p})$ in two pieces,
one of them proportional to $f(k)/f_{0}$ and the other proportional to
$f(p)/f_{0}$. The former involves the kernel

$\displaystyle K^{f_{\bm{k}}}_{\theta\theta}(\bm{k},-\bm{p},\bm{p})$
$\displaystyle=-\frac{k^{6}}{252|\bm{k}-\bm{p}|^{2}p^{4}}+\frac{k^{4}}{63|\bm{k}-\bm{p}|^{2}p^{2}}-\frac{11k^{4}}{252p^{4}}+\frac{19k^{2}|\bm{k}-\bm{p}|^{2}}{168p^{4}}-\frac{k^{2}}{42|\bm{k}-\bm{p}|^{2}}$
$\displaystyle\quad-\frac{41k^{2}}{504p^{2}}-\frac{p^{2}}{504k^{2}}+\frac{p^{2}}{63|\bm{k}-\bm{p}|^{2}}-\frac{p^{4}}{252k^{2}|\bm{k}-\bm{p}|^{2}}+\frac{13|\bm{k}-\bm{p}|^{2}}{126p^{2}}+\frac{5|\bm{k}-\bm{p}|^{2}}{168k^{2}}$
$\displaystyle\quad-\frac{5|\bm{k}-\bm{p}|^{4}}{63p^{4}}-\frac{19|\bm{k}-\bm{p}|^{4}}{504k^{2}p^{2}}+\frac{|\bm{k}-\bm{p}|^{6}}{72k^{2}p^{4}}-\frac{5}{126}$
$\displaystyle\equiv\displaystyle\sum\limits_{n_{1}=-2}^{2}\displaystyle\sum\limits_{n_{2}=-1}^{3}f^{f_{\bm{k}}}_{13,\theta\theta}(n_{1},n_{2})k^{-2(n_{1}+n_{2})}p^{2n_{1}}|\bm{k}-\bm{p}|^{2n_{2}},$
(4.24)

with matrix

$f^{f_{\bm{k}}}_{13,\theta\theta}(n_{1},n_{2})=\blockarray{cccccc}\text{\scriptsize$n_{2}=-1$}&\text{\scriptsize$0$}\text{\scriptsize$1$}\text{\scriptsize$2$}\text{\scriptsize$3$}\\\
\block{(ccccc)c}-\frac{1}{252}-\frac{11}{252}\frac{19}{168}-\frac{5}{63}\frac{1}{72}\text{\scriptsize$n_{1}=-2$}\\\
\frac{1}{63}-\frac{41}{504}\frac{13}{126}-\frac{19}{504}0\text{\scriptsize$-1$}\\\
-\frac{1}{42}-\frac{5}{126}\frac{5}{168}00\text{\scriptsize$0$}\\\
\frac{1}{63}-\frac{1}{504}000\text{\scriptsize$1$}\\\
-\frac{1}{252}0000\text{\scriptsize$2$}\\\ ,$ (4.25)

while the latter involves the kernel

$\displaystyle K^{f_{\bm{p}}}_{\theta\theta}(\bm{k},-\bm{p},\bm{p})$
$\displaystyle=-\frac{k^{6}}{126|\bm{k}-\bm{p}|^{2}p^{4}}+\frac{k^{4}}{72|\bm{k}-\bm{p}|^{2}p^{2}}+\frac{k^{4}}{72p^{4}}+\frac{k^{2}|\bm{k}-\bm{p}|^{2}}{168p^{4}}+\frac{k^{2}}{168|\bm{k}-\bm{p}|^{2}}$
$\displaystyle\quad+\frac{k^{2}}{252p^{2}}-\frac{5p^{2}}{126k^{2}}-\frac{11p^{2}}{504|\bm{k}-\bm{p}|^{2}}+\frac{5p^{4}}{504k^{2}|\bm{k}-\bm{p}|^{2}}+\frac{11|\bm{k}-\bm{p}|^{2}}{504p^{2}}+\frac{5|\bm{k}-\bm{p}|^{2}}{84k^{2}}$
$\displaystyle\quad-\frac{11|\bm{k}-\bm{p}|^{4}}{504p^{4}}-\frac{5|\bm{k}-\bm{p}|^{4}}{126k^{2}p^{2}}+\frac{5|\bm{k}-\bm{p}|^{6}}{504k^{2}p^{4}}+\frac{11}{504}$
$\displaystyle\equiv\displaystyle\sum\limits_{n_{1}=-2}^{2}\displaystyle\sum\limits_{n_{2}=-1}^{3}f^{f_{\bm{p}}}_{13,\theta\theta}(n_{1},n_{2})k^{-2(n_{1}+n_{2})}p^{2n_{1}}|\bm{k}-\bm{p}|^{2n_{2}},$
(4.26)

with matrix

$\
f^{f_{\bm{p}}}_{13,\theta\theta}(n_{1},n_{2})=\blockarray{cccccc}\text{\scriptsize$n_{2}=-1$}&\text{\scriptsize$0$}\text{\scriptsize$1$}\text{\scriptsize$2$}\text{\scriptsize$3$}\\\
\block{(ccccc)c}-\frac{1}{126}\frac{1}{72}\frac{1}{168}-\frac{11}{504}\frac{5}{504}\text{\scriptsize$n_{1}=-2$}\\\
\frac{1}{72}\frac{1}{252}\frac{11}{504}-\frac{5}{126}0\text{\scriptsize$-1$}\\\
\frac{1}{168}\frac{11}{504}\frac{5}{84}00\text{\scriptsize$0$}\\\
-\frac{11}{504}-\frac{5}{126}000\text{\scriptsize$1$}\\\
\frac{5}{504}0000\text{\scriptsize$2$}\\\ .$ (4.27)

Then, the approximation for $P^{13}_{cb,\theta\theta}(k)$ gives

$\displaystyle
P^{13}_{\theta\theta}(k)=k^{3}\frac{f(k)}{f_{0}}P_{L}(k)\left(\frac{f(k)}{f_{0}}\displaystyle\sum\limits_{m_{1}}c_{m_{1}}k^{-2\nu_{1}}M^{f_{\bm{k}}}_{13,\theta\theta}(\nu_{1})+\displaystyle\sum\limits_{m_{1}}c^{f}_{m_{1}}k^{-2\nu_{1}}M^{f_{\bm{p}}}_{13,\theta\theta}(\nu_{1})\right),$
(4.28)

where

$\displaystyle M^{f_{\bm{k}}}_{13,\theta\theta}(\nu_{1})$
$\displaystyle=-\frac{\tan(\nu_{1}\pi)}{14\pi(\nu_{1}+1)\nu_{1}(\nu_{1}-1)(\nu_{1}-2)(\nu_{1}-3)},$
(4.29) $\displaystyle M^{f_{\bm{p}}}_{13,\theta\theta}(\nu_{1})$
$\displaystyle=\frac{9\nu_{1}-7}{4}\frac{\tan(\nu_{1}\pi)}{28\pi(\nu_{1}+1)\nu_{1}(\nu_{1}-1)(\nu_{1}-2)(\nu_{1}-3)}.$
(4.30)

In Appendix B, we present the FFTLog approximations to obtain all the
functions $I^{m}_{n}(k)$ in eq. (1.2).

###  4.2 UV and IR corrections

Notice the function $I(z_{1},z_{2})$ vanishes if one of the arguments is zero
or a negative integer. For example, consider $z_{2}=0$. In this case, eq.
(4.19) becomes $\int_{\bm{p}}p^{-2z_{1}}=0$. On the other hand, it is well-
known that $\int_{\bm{p}}P_{L}(p)$ diverges. It seems contradictory because
the linear power spectrum can be decomposed into power laws, and as a
consequence of applying eq. (4.19) we would find zero as the final result. The
latter suggests that $I(z_{1},z_{2})$ gives inappropriate results when the
integral is divergent. This is because $I(z_{1},z_{2})$ calculates only the
finite part of the loop integral. Therefore, if the integral we are interested
in has an ultraviolet (UV) or an infrared (IR) divergence, to get the correct
answer, one simply has to add the UV/IR contribution by hand [48, 98].

Now, let us discuss the convergence properties of the loop integrals (3.31) –
(3.36). The behavior of fk-kernels under the UV $(k\ll p)$ and IR $(p\ll k)$
regimes are the same as $\Lambda$CDM, with the difference of some
multiplications by $f_{0}$ and $f(k\rightarrow\infty)$ factors. Then, when
considering the FFTLog decomposition, eq. (4.1), the convergence of the loop
integrals is determined by the bias parameter $\nu$, as much as in the
massless neutrino case. For example, if we take the UV/IR limit and consider
$P_{L}(k)\sim k^{\nu}$, one can find that the $P^{22}_{ab}$ contributions are
UV convergent for $\nu<1/2$ and IR convergent for $-1<\nu$. Then for
$-1<\nu<1/2$ the $P^{22}_{ab}$ integrals are convergent, and consequently, the
use of eq. (4.19) return the same results as with the traditional direct
computation. Nevertheless, when we choose values for the bias $\nu$ outside
the convergence range, the $P^{22}_{ab}$ integrals become UV or IR divergent.
Then, eq. (4.19) does not guarantee the correct answer because the divergent
pieces are not captured by dimensional regularization. Therefore, to obtain
the correct answer, the corresponding UV or IR piece must be added to the
result obtained through the FFTLog formalism [48, 98].

UV and IR corrections
---
| UV | IR
$P^{22}$ |  |
$\tfrac{1}{2}<\nu<\tfrac{3}{2}$ | $-3<\nu<-1$
$P^{22,\text{UV}}_{\delta\delta}(k)=\frac{9k^{4}}{196\pi^{2}}\int_{0}^{\infty}dp\,\frac{P^{2}_{L}(p)}{p^{2}}$ | $P^{22,\text{IR}}_{\delta\delta}(k)=P_{L}(k)k^{2}\sigma_{\Psi}^{2}$
$P^{22,\text{UV}}_{\delta\theta}(k)=-\frac{3k^{4}}{196\pi^{2}}\int_{0}^{\infty}dp\,\frac{f(p)}{f_{0}}\frac{P^{2}_{L}(p)}{p^{2}}$ | $P^{22,\text{IR}}_{\delta\theta}(k)=P^{L}_{\delta\theta}(k)k^{2}\sigma_{\Psi}^{2}$
$P^{22,\text{UV}}_{\theta\theta}(k)=\frac{k^{4}}{196\pi^{2}}\int_{0}^{\infty}dp\,\left(\frac{f(p)}{f_{0}}\right)^{2}\frac{P^{2}_{L}(p)}{p^{2}}$ | $P^{22,\text{IR}}_{\theta\theta}(k)=P^{L}_{\theta\theta}(k)k^{2}\sigma_{\Psi}^{2}$
$P^{13}$ |  |
$-1<\nu<1$ | $-3<\nu<-1$
$P^{13,\text{UV}}_{\delta\delta}(k)=-\frac{61}{105}P_{L}(k)k^{2}\sigma^{2}_{\Psi}$ | $P^{13,\text{IR}}_{\delta\delta}(k)=-P_{L}(k)k^{2}\sigma_{\Psi}^{2}$
$P^{13,\text{UV}}_{\delta\theta}(k)=-\left(\frac{23}{21}\frac{f(k)}{f_{0}}\sigma^{2}_{\Psi}+\frac{2}{21}\sigma^{2}_{v}\right)k^{2}P_{L}(k)$ | $P^{13,\text{IR}}_{\delta\theta}(k)=-P^{L}_{\delta\theta}(k)k^{2}\sigma_{\Psi}^{2}$
$P^{13,\text{UV}}_{\theta\theta}(k)=-\left(\frac{169}{105}\frac{f(k)}{f_{0}}\sigma^{2}_{\Psi}+\frac{4}{21}\sigma^{2}_{v}\right)k^{2}P^{L}_{\delta\theta}(k)$ | $P^{13,\text{IR}}_{\theta\theta}(k)=-P^{L}_{\theta\theta}(k)k^{2}\sigma_{\Psi}^{2}$
Table 1:  Leading UV/IR contributions for different values of bias $\nu$

A similar analysis for contributions $P^{13}_{ab}$ finds these loop integrals
are divergent for $\nu>-1$ and $\nu<-1$, i.e. the integrals never converge.
Thus, based on the $\nu$ value, we have to add the UV or IR divergence, as
appropriate. For example, for $\nu>-1$ the $P^{13}_{ab}$ loop integrals are UV
divergent. In this limit the leading contribution to $P^{13}_{\theta\theta}$
is

$P^{13,\text{UV}}_{\theta\theta}(k)=-\left(\frac{169}{105}\frac{f(k)}{f_{0}}\sigma^{2}_{\Psi}+\frac{4}{21}\sigma^{2}_{v}\right)k^{2}P^{L}_{\delta\theta}(k),$
(4.31)

with

$\sigma^{2}_{\Psi}\equiv\frac{1}{6\pi^{2}}\int_{0}^{\infty}dp\,P_{L}(p),\qquad\sigma^{2}_{v}\equiv\frac{1}{6\pi^{2}}\int_{0}^{\infty}dp\,P^{L}_{\delta\theta}(p).$
(4.32)

Notice that we have kept only the leading UV contribution, which is valid in
the range $-1<\nu<1$. For biases $\nu>1$, one has to consider subleading UV
corrections.

Similarly, if we consider a bias $\nu<-1$, the IR corrections are required.
For these biases, the leading IR contribution of $P^{13}_{\theta\theta}$ is

$P^{13,\text{IR}}_{\theta\theta}(k)=-P^{L}_{\theta\theta}(k)k^{2}\sigma_{\Psi}^{2},$
(4.33)

which holds for $-3<\nu<-1$. For smaller values of the bias, the subleading IR
contributions are necessary.

Returning to $P^{22}_{ab}$, we find that for the velocity spectra, the leading
UV and IR contributions are

$\displaystyle P^{22,\text{UV}}_{\theta\theta}(k)$
$\displaystyle=\frac{k^{4}}{196\pi^{2}}\int_{0}^{\infty}dp\,\left(\frac{f(p)}{f_{0}}\right)^{2}\frac{P^{2}_{L}(p)}{p^{2}},\qquad(\tfrac{1}{2}<\nu<\tfrac{3}{2})$
(4.34) $\displaystyle P^{22,\text{IR}}_{\theta\theta}(k)$
$\displaystyle=P^{L}_{\theta\theta}(k)k^{2}\sigma_{\Psi}^{2},\qquad\qquad\qquad\qquad\quad(-3<\nu<-1)$
(4.35)

In table 1 we summarize all the leading UV and IR contributions with the
corresponding range where they become necessary. Notice that $P^{22}_{ab}$ and
$P^{13}_{ab}$ are individually IR divergent for $\nu<-1$. Furthermore, notice
in the range $-3<\nu<-1$ the IR-contributions from $P^{22}_{ab}$ and
$P^{13}_{ab}$ cancel out, so they do not must to be added, as much as when
neutrinos are massless [105].

##  5 Model validation

We are now in position to validate our analytical model and numerical recipes
by comparing against the Quijote simulations. Quijote is a suite of 44,100
full N-body simulations whose fiducial cosmology is
$\\{\Omega_{m}=0.3175,\Omega_{b}=0.049,h=0.6711,n_{s}=0.9624,\sigma_{8}=0.834,M_{\nu}=0\\}$,
and contains additional neutrino cosmologies with total mass $M_{\nu}=0.1$,
$0.2$ and $0.4$ eV equally distributed among the three neutrino species. In
this work, we use mainly $M_{\nu}=0.4\,\text{eV}$ since this is the most
massive case, and hence the effects of having kernels beyond EdS are more
important. But also, we compare against the more likely value of
$M_{\nu}=0.1\,\text{eV}$. We use $N_{T}=100$ realizations for each cosmology,
containing $512^{3}$ cold dark matter ($cdm$) particles and $512^{3}$ neutrino
particles. We consider halos with masses in the range
$13.1<\log(M/h^{-1}M_{\odot})<13.5$ identified with a Friends-of-Friends
algorithm [106] over the $cdm$ particles, with a linking length parameter
$b=0.2$. These are the same halos utilized in [63, 65], which facilitates the
comparisons to previous works. Each realization has a volume of
$1\,(\,h^{-1}\text{Gpc})^{3}$, hence the total volume is
$100\,(\,h^{-1}\text{Gpc})^{3}$, allowing us to have very small statistical
errors and hence any departure from recovering the correct parameters from the
simulations are expected to come from effects that are not present or not well
modeled in our theory, e.g. 2-loop corrections, halo identifications, biasing
and effective field theory counterterms, and more important for us, to test
the fk-kernels method, since even small inaccuracies in the modeling can lead
to very large biased constraints on the parameters.

Figure 2: Contour plots for the posterior distributions at 68 and 95% confidence level computed with different covariance matrices, corresponding to $V_{1}$, $V_{5}$, $V_{25}$ and $V_{100}$ effective volume data sets. This is the case of mass $M_{\nu}=0.4\,\text{eV}$ at redshift $z=0.5$. The fittings are performed using the monopole and quadrupole of the power spectrum up to $k_{\text{max}}=0.2\,\,h\,\text{Mpc}^{-1}$. Vertical and horizontal dashed lines show the values of the simulations. The triangular plot with all the parameters, including nuisances, is shown in figure 10.  | Parameter | $\bm{V}_{1}$ | $\bm{V}_{5}$ | $\bm{V}_{25}$ | $\bm{V}_{100}$
---|---|---|---|---
$h$ | $0.700\pm 0.035$ | $0.686\pm 0.015$ | $0.6763\pm 0.0096$ | $0.6736\pm 0.0046$
$\Omega_{cdm}$ | $0.272\pm 0.023$ | $0.271\pm 0.013$ | $0.2643\pm 0.0073$ | $0.2626^{+0.0042}_{-0.0038}$
$M_{\nu}$ [eV] | — | $0.71^{+0.33}_{-0.47}$ | $0.39^{+0.16}_{-0.30}$ | $0.31\pm 0.11$
$\log(10^{10}A_{s})$ | $3.10^{+0.19}_{-0.28}$ | $3.25\pm 0.12$ | $3.288\pm 0.057$ | $3.305\pm 0.035$
$b_{1}$ | $1.83^{+0.28}_{-0.25}$ | $1.70\pm 0.13$ | $1.645\pm 0.066$ | $1.626\pm 0.038$
$b_{2}$ | $-0.3\pm 1.8$ | $-0.4\pm 1.1$ | $-0.46^{+0.56}_{-0.68}$ | $-0.53\pm 0.34$
$\alpha_{0}$ [$h^{-2}\,\text{Mpc}^{2}$] | $12\pm 81$ | $15\pm 35$ | $3^{+17}_{-20}$ | $1.8^{+8.9}_{-10}$
$\alpha_{2}$ [$h^{-2}\,\text{Mpc}^{2}$] | $0^{+31}_{-62}$ | $-26^{+19}_{-22}$ | $-28.9\pm 7.2$ | $-29.9\pm 3.9$
$\alpha^{shot}_{0}-1$ | $-0.91^{+1.0}_{-0.80}$ | $-1.01\pm 0.51$ | $-0.92\pm 0.25$ | $-0.90^{+0.16}_{-0.13}$
$\alpha^{shot}_{2}$ [$h^{-2}\,\text{Mpc}^{2}$] | $-8.5^{+5.9}_{-5.2}$ | $-7.7\pm 2.7$ | $-8.1\pm 1.1$ | $-8.07\pm 0.54$
Table 2:  1-dimensional constraints for the different covariance matrices,
corresponding to $V_{1}$, $V_{5}$, $V_{25}$ and $V_{100}$ effective volume
data sets. We show the mean of the posterior distributions with the 68%
confidence level intervals. We consider the case of massive neutrinos with
$M_{\nu}=0.4\,\text{eV}$ at redshift $z=0.5$. The fittings are performed using
the monopole and quadrupole of the power spectrum up to
$k_{\text{max}}=0.2\,\,h\,\text{Mpc}^{-1}$. The em dash “—” in the neutrino
mass for $V_{1}$ means no conclusive results because the posterior
distribution saturates the prior of $[0,2]\,\text{eV}$. The corresponding
triangular plot is shown in figure 2.

We construct four data sets to compare to. All of them having the average of
the 100 realizations as central points, but with different covariance matrices

$C_{N_{r}}=\frac{1}{N_{r}}C_{1},\qquad\text{with}\qquad
N_{r}=1,\,5,\,25,\,\text{and}\,100,$ (5.1)

with $C_{1}$ the covariance of just one realization, but using the average of
the $N_{T}=100$ realizations:

$C_{1}(k_{i},k_{j})=\left\langle\Big{(}\bar{P}(k_{i})-\hat{P}(k_{i})\Big{)}\Big{(}\bar{P}(k_{j})-\hat{P}(k_{j})\Big{)}\right\rangle$
(5.2)

where $\hat{P}(k_{i})$ is the value of the power spectrum of a single
realization at bin $k_{i}$ and $\bar{P}(k_{i})=\langle\hat{P}(k_{i})\rangle$
is the mean over the ensemble of realizations at bin $k_{i}$. It may be more
precise the use of $N_{r}\leq N_{T}$ random subsamples of realizations and fit
to their average [107], instead of fitting to the average of the $N_{T}$
realizations as we do. However, we have made tests down to $N_{r}=5$ and found
no significant difference between the two approaches. Ultimately, our interest
is to have a consistent way to compare between different effective volumes to
validate our modeling. Another way to understand our approach is under the
assumption of a Gaussian covariance matrix, scaling with the total volume.

We choose four different volumes, corresponding to _1_) the volume of one
realization $V_{1}=1\,h^{-1}\text{Gpc}$ ($N_{r}=1$), _2_) the volume of the
$N_{T}=100$ realizations $V_{100}=100\,h^{-1}\text{Gpc}$ ($N_{r}=N_{T}=100$),
_3_) we also consider a “BOSS-like” volume of $V_{5}=5\,h^{-1}\text{Gpc}$
($N_{r}=5$), and _4_) a “DESI-like” volume $V_{25}=25\,h^{-1}\text{Gpc}$
($N_{r}=25$).

Figure 3: Comparison on the constraints of the sum of the neutrino masses with
the use of fk\- and EdS-kernels for the volumes $V_{25}$ (left panel) and
$V_{100}$ (right panel).

Our baseline fitting consist on letting free the cosmological parameters
$\Omega_{cdm}$, $h$, $\log(10^{10}A_{s})$ and $M_{\nu}$, where
$\Omega_{cb}=\Omega_{cdm}+\Omega_{b}$ with the baryon abundance fixed to the
Quijote fiducial value $\omega_{b}=0.02207$. We consider also the bias
parameters $b_{1}$, $b_{2}$, the EFT counterterms $\alpha_{0}$, $\alpha_{2}$,
and the shot noise $\alpha_{0}^{shot}$ and $\alpha_{2}^{shot}$. The tidal and
non-local bias are fixed by co-evolution theory [108, 109, 93], implying that

$b_{s^{2}}=-\frac{4}{7}(b_{1}-1),\qquad b_{3nl}=\frac{32}{315}(b_{1}-1).$
(5.3)

These expressions neglect an early time, Lagrangian tidal bias and are
obtained using EdS evolution. However, we find that they yield good results
also in the presence of massive neutrinos when the biased field is $cb$, so we
use them to reduce the number of free parameters of the theory. Later we will
relax this assumption and let free the biases $b_{s^{2}}$ and $b_{3nl}$, as
well as the baryon abundance $\omega_{b}$. We consider uninformative priors on
the cosmological and nuisances, with the exception of the neutrino mass prior,
which is uniform over the interval $[0,2]\,\text{eV}$. The shot noise constant
is fixed to $P_{shot}=1/\bar{n}_{x}=4719.7\,h^{-3}\,\text{Mpc}^{3}$, hence any
departure from the Poissonian process is encoded in the values of
$\alpha_{0}^{shot}$ and the tilt $\alpha_{2}^{shot}$. The linear $cb$ power
spectrum is obtained from the Einstein-Boltzmann code
CLASS131313https://lesgourg.github.io/class_public/class.html [110] which
serves as an input of FOLPS$\nu$ giving us the galaxy power spectrum
multipoles that we compare against the simulations using a Gaussian
likelihood, $L\propto\exp(-\chi^{2}/2)$, where $\chi^{2}$ is defined as
$\chi^{2}\equiv D^{T}C_{N_{r}}^{-1}D$ and $D$ is the residual of the data
vector and the model and $C_{N_{r}}$ is the covariance matrix in eq. (5.1).

Figure 4: Effects on the posterior distributions of the cosmological
parameters of considering different values of the maximum wave-number
$k_{\text{max}}=0.15,\,0.2,\,0.25\,h\,\text{Mpc}^{-1}$. We use the DESI-like
volume $V_{25}$. Vertical and horizontal dashed lines show the values of the
simulations. The triangular plot with all the parameters, including nuisances,
is shown in figure 13.  Figure 5: Effects on the posterior distributions of
the cosmological parameters of letting free the biases $b_{s^{2}}$ and
$b_{3nl}$ instead of fixing them to co-evolution theory [93]. Also, we vary
the baryon abundance with a Gaussian prior $\omega_{b}=0.02207\pm 0.00015$. We
use the DESI-like volume $V_{25}$. Vertical and horizontal dashed lines show
the values of the simulations.

To sample the parameter space, we run the code
emcee141414https://emcee.readthedocs.io/ [111] which utilizes the affine-
invariant ensemble sampler method for MCMC [112]. Finally, we obtain the plots
and confidence intervals using the GetDist Python package [113].

We compare against the monopole and quadrupole of the power spectrum up to a
wave-number $k_{\text{max}}=0.20\,h\,\text{Mpc}^{-1}$ for redshift $z=0.5$. In
table 2 we show the mean values and 0.68 c.l. intervals of the 1-dimensional
marginalized posterior distributions, where we considered the four volumes
$V_{1}$, $V_{5}$, $V_{25}$, $V_{100}$. A triangular plot for these fits is
shown in figure 2 where we only show the cosmological parameters (the plot for
all the parameters, including nuisances, is given in Appendix C, figure 10).
We notice that for the smallest volume, $V_{1}$, we are not able to recover
the mass of the neutrinos, but the posterior saturates the prior of
$[0,2]\,\text{eV}$, as it is represented with an em dash in table 2. However,
for the other effective volume cases we indeed recover all the cosmological
parameters inside the 0.68 confidence intervals, it is easily noticed that
even for the largest volume, where we expect that the modeling systematics
dominate the error, we are able to obtain back all the parameters.

Figure 6: Contour plots for the posterior distributions at 68 and 95%
confidence level computed with different covariance matrices, corresponding to
$V_{25}$ and $V_{100}$ effective volume data sets. This is the case of mass
$M_{\nu}=0.1\,\text{eV}$ at redshift $z=0.5$. The fittings are performed using
the monopole and quadrupole of the power spectrum up to
$k_{\text{max}}=0.2\,\,h\,\text{Mpc}^{-1}$. Vertical and horizontal dashed
lines show the values of the simulations.

We perform the same exercise using EdS kernels, by means of simply setting
$f(k)\rightarrow f_{0}$ in our FOLPS$\nu$ code, and although the fittings are
good as well, when our method of fk-kernels is used the results for $M_{\nu}$
are modestly better. This is shown in figure 3 where we plot the marginalized
1-dimensional posterior for $M_{\nu}$ for volumes $V_{25}$ and $V_{100}$. The
estimation for the sum of the neutrino masses is within the 68% limits for
both cases when using the fk-kernels, but it fails when using EdS with volume
$V_{100}$, although is recovered within the 95% limits. The corresponding
triangular plot is shown in Appendix C, figures 11 and 12, where all the rest
of the cosmological parameters are recovered within the 0.68 c.l. intervals
both with EdS and with fk kernels.

| Parameter | $\bm{V}_{25}$ | $\bm{V}_{100}$
---|---|---
$h$ | $0.6800\pm 0.0099$ | $0.6744^{+0.0048}_{-0.0056}$
$\Omega_{cdm}$ | $0.2763\pm 0.0077$ | $0.2709^{+0.0042}_{-0.0052}$
$M_{\nu}$ [eV] | $<0.493$ | $<0.219$
$\log(10^{10}A_{s})$ | $3.166^{+0.057}_{-0.065}$ | $3.168^{+0.027}_{-0.032}$
$b_{1}$ | $1.696\pm 0.073$ | $1.667\pm 0.041$
$b_{2}$ | $-0.48^{+0.63}_{-0.52}$ | $-0.28\pm 0.34$
$\alpha_{0}$ [$h^{-2}\,\text{Mpc}^{2}$] | $12^{+18}_{-23}$ | $-0.9^{+8.2}_{-9.7}$
$\alpha_{2}$ [$h^{-2}\,\text{Mpc}^{2}$] | $-25.5\pm 7.7$ | $-28.0\pm 3.4$
$\alpha^{shot}_{0}-1$ | $-1.00\pm 0.22$ | $-0.95\pm 0.15$
$\alpha^{shot}_{2}$ [$h^{-2}\,\text{Mpc}^{2}$] | $-8.0\pm 1.1$ | $-7.98\pm 0.53$
Table 3:  1-dimensional constraints for the different covariance matrices, corresponding to $V_{25}$ and $V_{100}$ effective volume data sets. We show the mean of the posterior distributions with the 68% confidence level intervals. We consider the case of massive neutrinos with $M_{\nu}=0.1\,\text{eV}$ at redshift $z=0.5$. The fittings are performed using the monopole and quadrupole of the power spectrum up to $k_{\text{max}}=0.2\,\,h\,\text{Mpc}^{-1}$. The corresponding triangular plot is shown in figure 6. Figure 7: Contour plots for the posterior distributions at 0.68 and 0.95 c.l. computed with different covariance matrices, corresponding to $V_{25}$ and $V_{100}$ effective volume data sets. This is the case of mass $M_{\nu}=0.4\,\text{eV}$ at redshift $z=1.0$. The fittings are performed using the monopole and quadrupole of the power spectrum up to $k_{\text{max}}=0.2\,\,h\,\text{Mpc}^{-1}$. Vertical and horizontal dashed lines show the values of the simulations.  | Parameter | $\bm{V}_{25}$ | $\bm{V}_{100}$
---|---|---
$h$ | $0.6726\pm 0.0089$ | $0.6694^{+0.0056}_{-0.0049}$
$\Omega_{cdm}$ | $0.2596\pm 0.0062$ | $0.2578^{+0.0043}_{-0.0038}$
$M_{\nu}$ [eV] | $<0.460$ | $0.26^{+0.12}_{-0.19}$
$\omega_{b}$ | $0.02208^{+0.00014}_{-0.00020}$ | $0.02209^{+0.00018}_{-0.00016}$
$\log(10^{10}A_{s})$ | $3.347^{+0.090}_{-0.081}$ | $3.340^{+0.061}_{-0.051}$
$b_{1}$ | $2.27^{+0.12}_{-0.14}$ | $2.270^{+0.070}_{-0.084}$
$b_{2}$ | $0.4^{+1.7}_{-1.5}$ | $-0.8^{+1.6}_{-1.4}$
$b_{s^{2}}$ | $2.4^{+3.3}_{-2.8}$ | $-1.0^{+2.2}_{-3.5}$
$b_{3nl}$ | $0.23^{+0.38}_{-0.46}$ | $0.14^{+0.20}_{-0.25}$
$\alpha_{0}$ [$h^{-2}\,\text{Mpc}^{2}$] | $-10.5^{+1.6}_{-1.8}$ | $-7.6^{+4.5}_{-5.1}$
$\alpha_{2}$ [$h^{-2}\,\text{Mpc}^{2}$] | $-14.3\pm 2.7$ | $-15.2^{+8.9}_{-11}$
$\alpha^{shot}_{0}-1$ | $-0.71^{+0.29}_{-0.34}$ | $-0.78^{+0.15}_{-0.20}$
$\alpha^{shot}_{2}$ [$h^{-2}\,\text{Mpc}^{2}$] | $-5.25^{+0.84}_{-1.1}$ | $-5.10\pm 0.83$
Table 4:  1-dimensional constraints for the different covariance matrices,
corresponding to $V_{25}$ and $V_{100}$ effective volume data sets. We show
the mean of the posterior distributions with the 68% confidence level
intervals. We consider the case of massive neutrinos with
$M_{\nu}=0.4\,\text{eV}$ at redshift $z=1.0$. The fittings are performed using
the monopole and quadrupole of the power spectrum up to
$k_{\text{max}}=0.2\,\,h\,\text{Mpc}^{-1}$. The corresponding triangular plot
is shown in figure 7.

So far, we have fixed the maximum wave-number value to
$k_{\text{max}}=0.2\,h\,\text{Mpc}^{-1}$ since we found it optimal under our
testings. In figure 4 we show contour and 1-dimensional posterior
distributions for three different values, $k_{\text{max}}=0.15,\,0.2$ and
$0.25\,h\,\text{Mpc}^{-1}$, using the DESI-like volume $V_{25}$. The case of
$k_{\text{max}}=0.25\,h\,\text{Mpc}^{-1}$ underestimates the neutrino mass
$M_{\nu}$ considerably and overestimates the amplitude of fluctuations
$A_{s}$, or alternatively underestimate the linear local bias which was
measured for this halo sample to be $b_{1}=1.665$ using the power spectrum
[65] and the correlation function [63].151515In those works, the cosmological
parameters are fixed to the Quijote fiducial values, while the nuisances are
fitted to the simulations.  On the other hand, using either
$k_{\text{max}}=0.15$ or $0.20\,h\,\text{Mpc}^{-1}$ we recover the
cosmological parameters within the 68% limits, but
$k_{\text{max}}=0.20\,h\,\text{Mpc}^{-1}$ does it with smaller standard
deviations. Hence, the justification of our choice of $k_{\text{max}}$. The
failure to recover the correct parameters when using large values for
$k_{\text{max}}$ has been observed in previous works, for example in [63]
where only the nuisances were left free. In that work it was clear that as
larger was $k_{\text{max}}$, the more inaccurate was the inferred value for
$b_{1}$. Since this parameter is largely degenerate with the primordial
fluctuation amplitude $A_{s}$, it is natural to expect a wrong estimation of
the cosmological parameters as well. We attribute this to the fact that at
small scales the form of the analytical power spectrum is largely dominated by
the shot-noise and the EFT counterterms, which on the other hand become
tightly constrained by the data large $k$ modes (see the red contours in the
complementary figure 13), leaving little room to the cosmological parameters
to find their “true” values.

Now, we present the results when letting free the biases $b_{s^{2}}$ and
$b_{3nl}$ instead of fixing them to the co-evolution relations (5.3). We show
the triangular plots of posterior distribution using the DESI-like volume
$V_{25}$ in figure 5 with (blue coloured) and without (green coloured) co-
evolution, noticing no significant change in the estimation of parameters.
Also, in this figure we show the results when we let free the baryon abundance
with a Gaussian prior centered at the fiducial value of the simulations and
standard deviation from Planck data 2018 $\omega_{b}=0.02207\pm 0.00015$. In
this case, we also obtain very similar results, but with slight indications of
improvement in the estimation for the $M_{\nu}$ best fit, as can be inferred
by looking at the position of the maxima of the 1-dimensional posterior
distributions of $M_{\nu}$ in figure 5 with and without $\omega_{b}$ fixed. We
think this is because $\Omega_{\nu}\ll\Omega_{b}$, and hence even a small
allowed change in the baryon abundance permits significant variations of the
neutrino abundance.

We show our results when we compare against the massive neutrinos case
$M_{\nu}=0.1\,\text{eV}$ at redshift $z=0.5$ and using a maximum wave-number
$k_{\text{max}}=0.2\,h\,\text{Mpc}^{-1}$. These are shown in the triangular
plots of figure 6 and table 3. We found upper bounds for the neutrino mass,
being $M_{\nu}<0.493\,\text{eV}$ for the effective volume $V_{25}$ and
$M_{\nu}<0.219\,\text{eV}$ for $V_{100}$, both at 0.68 c.l. For the rest of
the cosmological parameters we notice that we are not able to recover all of
them for $V_{100}$ within the 68% limits, although they are inside the 95%
limits, while for $V_{25}$ the discrepancies are smaller, and the parameters
lie inside the 0.68 c.l. intervals, except the case of the primordial
fluctuations amplitude $A_{s}$, which is overestimated. We think the reason
for this mismatch is the well-known degeneracy that exists between $A_{s}$ and
$M_{\nu}$. This leads us to believe of that neutrino mass will not be measured
with ongoing and near future surveys using 2-point functions alone, but we
should use additional information coming from higher order statistics as the
bispectrum [114, 115, 116] or even other summary statistics as the marked
power spectrum and correlation function [117, 118, 119, 78].

Finally, we test our modeling and code against the massive neutrinos case
$M_{\nu}=0.4\,\text{eV}$, but at a higher redshift of $z=1.0$, since these
redshifts will be covered by current stage IV experiments as with the DESI-
ELGs. We fit up to a maximum wave-number
$k_{\text{max}}=0.2\,h\,\text{Mpc}^{-1}$ and let free the biases $b_{3nl}$ and
$b_{s^{2}}$, and the baryonic matter abundance $\omega_{b}$. As with the case
of $M_{\nu}=0.1\,\text{eV}$ at redshift $z=0.5$, we test the cases of $V_{25}$
and $V_{100}$. Our results are shown in the triangular plots of figure 7 and
Table 4. We notice we are able to recover all the parameters for the two
effective volumes at 0.68 c.l., except the case of $M_{\nu}$ for the effective
volume $V_{100}$, which shows a small mismatch, but still recovered at 0.95
c.l.

##  6 Code FOLPS$\nu$

Fast One Loop Power Spectrum in the presence of massive neutrinos (FOLPS$\nu$)
is a Python code that computes the redshift space power spectrum multipoles in
a fraction of second. It receives as input the linear power spectrum of matter
or $cb$ field and the set of cosmological parameters
$\\{\Omega_{m},\,h,\,M_{\nu}\\}$ and nuisance parameters
$\\{b_{1},b_{2},b_{s^{2}},b_{3nl},\alpha_{0},\alpha_{2},\alpha_{4},\tilde{c},\alpha_{0}^{shot},\alpha_{2}^{shot}\\}$.
(Alternatively, one can use the bias parameters $b_{\mathcal{G}_{2}}$ and
$b_{\Gamma_{3}}$, instead of $b_{s^{2}}$ and $b_{3nl}$.) The code computes the
EFT galaxy redshift space power spectrum given in eq. (3.3) for both the
wiggle and non-wiggle (also computed by the code) linear power spectra.
Thereafter, it uses the IR-resummation expression in eq. (3.3.2). Finally, the
code computes the Legendre multipoles using the integral of eq. (3.62). All
the loop integrals are computed using the formulae of §4 and Appendix B. The
code obtains $f_{0}$ by solving directly the differential equation for the
linear growth $D_{+}$ at very large scales, using as input the matter
abundance $\Omega_{m}$ and neglecting relativistic species. The Hubble
dimensionless rate $h$ is used to obtain the non-wiggle linear power spectrum
using the fast sine transform recipe of [120]. The SPT kernel
$G_{1}(k)=f(k)/f_{0}$ is obtained with the Hu-Eisenstein approximation of ref.
[7].

The code obtains $f_{0}$ by solving directly the differential equation for
$D_{+}$ at very large scales, as in the massless neutrino case, using as input
the matter abundance $\Omega_{m}$ and neglecting relativistic species.

FOLPS$\nu$ employs the standard libraries NumPy161616We recommend to use NumPy
versions $\geq$ 1.20.0. For older versions, one needs to rescale $c_{m}$,
$c_{m}^{f}$ and $c_{m}^{ff}$ by a factor $1/N$. and SciPy. Other inputs are
the FFTLog parameters $k_{\text{min}}$, $k_{\text{max}}$, $N$ and the bias
$\nu$ introduced in 4. In addition, the code has the capability to switch
between fk\- and EdS-kernels, and it can also account for the Alcock-Paczynski
effect if required [121].

To compute the Fast Fourier Transforms and construct the $M$ matrices,
FOLPS$\nu$ takes $N$ points logarithmically spaced over the interval
$(k_{\text{min}}=10^{-7}\,h\,\text{Mpc}^{-1},k_{\text{max}}=100\,h\,\text{Mpc}^{-1})$.
In the following we test the use of different values $N=64$, $128$, $256$ and
$1000$. We compute the redshift space power spectrum multipoles for a set of
cosmological and nuisance parameters, first using direct integration (with a
very high precision) and then through the FOLPS$\nu$ code. Figure 8 shows the
relative difference between these two approaches for the different numbers of
sampling points. As expected, by sampling with $N=64$ we obtain the highest
relative difference with respect to the direct integration, while the cases of
$N=$ 128, 256 and 1000 provide similar accuracy. In particular, $N=128$
produces a maximal error of $\sim 0.1\%$ for the monopole and quadrupole, and
$\sim 0.2\%$ for hexadecapole up to $k=0.35\,h\,\text{Mpc}^{-1}$, which is
even smaller when sampling up to smaller wave-numbers. In all the
computations, we use the default bias parameter $\nu=-0.1$ for all the
functions except those required only for biased tracers, that is for eqs.
(3.49) – (3.57), for which we use $\nu=-1.51$. We have tested that these
choices yield good convergence.

Figure 8: Accuracy of the code FOLPS$\nu$ when using the sampling points
$N=64$, $128$, $256$ and $1000$. The comparisons are against the direct
integration, using as much precision as possible.

Let us now discuss the performance of our code. The computational times
presented here were obtained with a standard personal computer with the
following specifications: ASUS ZenBook with Cpu: 11th Gen Intel(R) Core(TM)
i7-1165G7 @ 2.80GHz, RAM: 16 GB, OS: Ubuntu 20.04.4 LTS.

Figure 9 shows the estimated end-to-end time taken by the FOLPS$\nu$ code as a
function of the number of sampling points $N$, here is not included the time
taken by the classy module of the code CLASS to compute the linear $cb$ power
spectrum. To these computations we have used 191 $k$ points over the interval
$[0.01,0.5]\,h\,\text{Mpc}^{-1}$. The black dots represent the $N$ values
where the code was evaluated, while the black solid line refers to the average
execution time, and the gray regions are the standard deviations of 100 runs.
In table 5 we decompose some of these times into the different pieces executed
during the code evaluation, the symbol “—” in the cells means that the time is
unchanged since those pieces do not depend on the FFTLog parameter $N$. From
the table we observe that the Hu-Eisenstein approximation is the fastest piece
of the code, the next on this list are the inclusion of the IR-resummations
and the computation of multipoles from eq. (3.62). Also note that the
cosmology-dependent terms $c_{m}$ ($c^{f}_{m}$ and $c^{ff}_{m}$ for massive
neutrinos cosmologies) are computed quickly through the use of the FFT
algorithms, while most of the time is spent calculating matrix
multiplications, especially the $P_{22}$-type.

We notice the computational time of our Python code is competitive to others,
as for example, the C language code class-pt, see table 1 in [32]. (Other
codes, as PyBird [122], FAST-PT [44, 45] and velocileptors [90] run in a
similar amount of time as class-pt.) We think that this is because our
expressions entering eq. (1.2) have been analytically reduced to perform a
smaller number of computations. Further, despite using fk-kernels, the number
of $M_{22}$ matrix multiplications is equal to $2\times 26$ against $2\times
24$ when EdS are used (the factor of 2 comes from the fact that we must
compute all the matrix multiplications using both wiggle and non-wiggle linear
power spectra). Hence our method is expected to run in a similar amount of
time than EdS. If one considers that this code is thought to be used inside a
pipeline containing a Boltzmann code and an MCMC sampler, the saving time when
using EdS kernels is negligible.

Finally, since the matrices $M$ do not depend on the cosmology, they can be
pre-computed and stored. However, it may be convenient to allow the user to
select other different values of $N$ from those presented here. For that
reason, our code computes the matrices $M$, this is the slowest part of the
code, and we do not include it in figure 9 and table 5 because this
calculation is performed only one time. That is, the first time the code is
called, it computes the $M$ matrices and stores them for the rest of the runs,
which can be of the order of thousands in parameter estimations. This is shown
in a Jupyter notebook released with the code. For the cases of $N=64$, $128$
and $256$ the computational time for all the $M$ matrices are 0.118, 0.400,
and 1.251 seconds, respectively.

Figure 9: Performance: Estimated end-to-end time taken by the FOLPS$\nu$ code as a function of the number of sampling points $N$.  $N$ | $f(k)/f_{0}$ | $P_{nw}(k)$ |  $c_{m}$, $c^{f}_{m}$, $c^{ff}_{m}$ |  $P_{22}$-type |  $P_{13}$-type |  $P_{s}^{\text{IR}}(k,\mu)$, $P_{\ell}(k)$ | Total
---|---|---|---|---|---|---|---
64 | 0.00136 | 0.0535 | 0.0104 | 0.0477 | 0.0154 | 0.00787 | 0.136
128 | — | — | 0.0179 | 0.137 | 0.0183 | — | 0.236
256 | — | — | 0.0317 | 0.435 | 0.0254 | — | 0.555
Table 5:  Computing time in [s] for the different terms for $N=64$, $128$, and
$256$. The “—” means that the computing time is almost unchanged compared to
the case $N=64$. The total number of computed of wave-numbers $k$ is 191.

##  7 Conclusions

During the upcoming years, the increasing constraining power of galaxy surveys
will probably yield the tightest bounds on the sum of neutrino masses. The
most important cosmological effect is the suppression of the power spectrum
and correlation function below the free-streaming scale, which arises because
neutrinos cannot be confined at those scales because their large velocity
dispersion counteract the gravitational collapse and prevent neutrinos from
forming structure. The amplitude and location of these effects depend on the
neutrino masses in a well understood physical process. However, the free-
streaming lies at the onset of non-linearities, and as such the modeling of it
requires the use of perturbation theory and non-linear effects as Effective
Field Theory and IR-resummations. On the other hand, as surveys become wider
and deeper they cover larger volumes and the statistical errors become smaller
and systematic errors dominate. As such, improvements on the theoretical
description would be very beneficial for data analyses. In this work we have
continued the development of a comprehensive perturbation theory for
cosmologies in the presence of massive neutrinos building on the works [63,
64, 65] which compute kernels beyond EdS for the combined cold dark matter and
baryonic fluid which is biased to obtain galaxy statistics. The use of more
adequate kernels than EdS has been studied in several works, and its utility
was shown, particularly the Lagrangian approach, when used in COmoving
Lagrangian Acceleration (COLA) implementations [123, 59] as well as other
approximate methods to generate mocks [124], and to initialize full N-body
simulations [125]. However, their use for comparing against data has been used
very little, if any. There are two important reasons for this: _1_) the so
called full-shape methods that compare directly the consistent Effective field
theory for Large Scales Structure to real data have reached sufficient
masterization only recently [29, 28] in order to be used to constrain
cosmological parameters; and, _2_) the use of complete kernels for cosmologies
that contains additional scales, the neutrino mass in our case, requires
lengthy computations of 1-loop integrals with kernels that are not known
analytically, but should be found through differential equations at each bin
of the 2-dimensional integrals.171717We want to mention that this situation is
encountered in some modified gravity models, such as $f(R)$ theories,
symmetron fields, or more generally in a large sector of the Hordenski sector;
see e.g., [126, 127, 62, 128, 64] and Appendix B of [129]. In those cases
there is an additional scale given by the associated scalar field (effective)
mass $m_{eff}$. Clearly, there is no free-streaming effect, but a universal
fifth force of range $1/m_{eff}$. However, the presence of the additional
scale leads to a function $A(k,t)$ different than in $\Lambda$CDM (where
$A(t)=3/2\,\Omega_{m}H^{2}$) and ultimately to a $k$-dependence on the linear
growth function and growth rate. For this reason the method developed in this
work, particularly the fk-kernels formalism, is potentially useful also for
modified gravity theories. We are currently exploring this possibility and
hope to present the results elsewhere in the near future.

We want to finish by mentioning that the results found in this work can be
applied also to some modified gravity models, such as $f(R)$ theories,
symmetron fields, and more generally to a large sector of the Hordenski theory
(see e.g., [126, 127, 62, 128, 64] and Appendix B of [129]). In those cases,
there is an additional scale given by the associated scalar field (effective)
mass $m_{eff}$. Clearly, there is no free-streaming effect, but a universal
fifth force of range $1/m_{eff}$. However, the presence of the additional
scale leads to a function $A(k,t)$ different than $=3/2\,\Omega_{m}H^{2}$, and
ultimately to a $k$-dependence on the growth rate. For this reason the method
developed in this work, particularly the fk-kernels formalism, is potentially
useful also for modified gravity theories. We are currently exploring this
possibility and hope to present the results elsewhere in a near future.

However, the presence of the additional scale leads to a function $A(k,t)$
different than $=3/2\,\Omega_{m}H^{2}$, and ultimately to a $k$-dependence on
the growth rate. For this reason the method developed in this work,
particularly the fk-kernels formalism, is potentially useful also for modified
gravity theories. We are currently exploring this possibility and hope to
present the results elsewhere in a near future.

In this work, following the results of [65], we have developed a method that
keeps the scale-dependent growth rates $f(k,t)$ as inherited from the linear
theory non-local relation between velocity and density fields given in eq.
(1.1) [see e.g. eq. (3.37)], while the rest of the pieces are maintaining
equal to their EdS counterparts. In this way, we are able to construct kernels
with precise analytical expressions (up to the growth rate), but yet the free-
streaming scale is present also at the non-linear level and not only through
the input linear power spectrum. The effects of having these more precise
kernels are about, or smaller than, the 1 per-cent level for neutrino masses
around $0.1\,\text{eV}$, but when spanning over uninformative priors one can
go up to large masses of $\gtrsim 1\,\text{eV}$, where the use of EdS kernels
is at least questionable. Hence, we argue that for a correct sampling of the
parameter space, one needs to use a more complete theory. The main objection
against it, is not that such a theory is missing, since several PTs have
existed from long time [50, 51], but that the computational time of loop
integrals is prohibitive to be used in MCMC algorithms. However, with the
method presented in this work, named by us as fk, we are able to compute the
loop integrals in a time comparable to that in EdS. Further, since our kernels
can now be written as the sum of terms
$p^{2n_{1}}|\bm{k}-\bm{p}|^{2n_{2}}k^{-2(n_{1}+n_{2})}$ times growth rates, we
are able to use FFTLog methods that improve the accuracy and reduce even more
the time spent when computing the loop integrals.

We release the code FOLPS$\nu$ that is fed with the linear $cb$ power spectrum
and a set of cosmological and nuisance parameters, and it computes the
galaxy/tracers redshift space power spectrum, including 1-loop corrections,
effective field theory counterterms, biasing and IR-resummations. Despite our
code is written in Python and goes beyond EdS, it is competitive in time with
other existing codes. We attribute this to the fact that our final expansion
is different, since we are able to write the loop power spectrum in the form
of eq. (1.2), analytically reducing the number of computations performed by
our code. Further, our code does not require any considerably more amount of
time than the use of EdS, since actually the number of matrix multiplications
is only 10% larger when using fk-kernels.

In §5 we validated our model and code by comparing against the Quijote suite
of simulations. We use as our baseline model the most massive case available
in these simulations, $M_{\nu}=0.4\,\text{eV}$, because in this case the use
of a more proper theory becomes more relevant. Further, we use halos and a
very large volume of $(100\,\,h^{-1}\text{Gpc})^{3}$, corresponding to the 100
realizations in the simulations, in order to have very small statistical
errors. As such, any failure in recovering the cosmological parameters of the
simulations is likely to be a cause of a wrong modeling in our theory. The
most important plot of the model validation is shown in figure 2, where we fit
to four different data sets, all with the same central points but the
corresponding covariance scaling with the inverse of the volume, for effective
volumes of $1,\,5,\,25\,\,\text{and}\,\,100\,h^{-3}\,\text{Gpc}^{3}$, named
$V_{1}$, $V_{5}$, $V_{25}$ and $V_{100}$, respectively. We find that we are
able to recover all the cosmological parameters within the 68% confidence
limits, even for the largest volume where the errors are very small. We
further compare our constraints with those obtained with the use of EdS,
finding modest but noticeable improvements when using our method. Curiously,
this happens only for the $M_{\nu}$ parameter, with EdS failing to recover it
at 0.68 c.l. for $V_{100}$.

In our baseline fitting, we fixed the baryon abundance to the simulation’s
fiducial value. However, we also test the effect of varying it with a Gaussian
prior with width given by Planck 2018\. By doing this, being quite small, the
neutrino abundance has enough room to move a little bit without affecting the
baryon abundance considerably, and ultimately allowing a slight improvement
for the best fit.

We also compare against the case of $M_{\nu}=0.1\,\text{eV}$ finding not as
good agreement than the other case. Particularly, we found the mass of the
neutrinos is overestimated, and at the same time, the fluctuations primordial
amplitude $A_{s}$ is overestimated as well, but in an even larger amount.
Since these parameters are highly degenerate, we argue that the galaxy 2-point
functions alone will not be sufficient to measure the absolute mass scale of
the neutrinos with current and upcoming surveys, which is a somewhat believed
result; e.g. [114].

We believe that calibrations in our method as the choice of more optimal
parameters should be done by using mocks and ultimately comparisons with real
data. This is a work in progress that we hope to present soon elsewhere.
Particularly with the view of upcoming results for the first year DESI data.

## Acknowledgments

We would like to thank Willem Elbers, Carlos Frenk, Baojiu Li, Martin White
and Arka Banerjee for useful discussions and suggestions. AA and HEN are
supported by CONACyT Ciencia de Frontera grant No.319359 and CONACyT grant
283151. AA also acknowledges CONACyT grant 102958. MV, HEN and SF acknowledges
PAPIIT IN108321 and PAPIIT IA103421. MV acknowledges CONACyT grant A1-S-1351.
This research was partially supported through computational and human
resources provided by the LAMOD UNAM project through the clusters Atocatl and
Tochtli. LAMOD is a collaborative effort between the IA, ICN and IQ institutes
at UNAM.

##  Appendix A Functions A and D

Function $A^{\text{TNS}}(k,\mu)$ was originally introduced in [86] for the TNS
(Taruya-Nishimichi-Saito) model. It is given by eq. (3.41) and can be written
as

$\displaystyle A^{\text{TNS}}(k,\mu)$
$\displaystyle=\mu^{2}\big{[}f_{0}I^{1,udd}_{1}(k)+f_{0}^{2}I^{2,uud}_{1}(k)\big{]}+\mu^{4}\big{[}f_{0}^{2}I^{2,uud}_{2}(k)$
$\displaystyle\quad+f_{0}^{3}I^{3,uuu}_{2}(k)\big{]}+\mu^{6}f_{0}^{3}I^{3,uuu}_{3}(k),$
(A.1)

where

$\displaystyle I^{m,ijk}_{n}(k)$ $\displaystyle=2\int_{\bm{p}}{\text{\bf
a}}_{nm}(\bm{k},\bm{p})P_{L}(k)P_{L}(p)+\int_{\bm{p}}\tilde{\mathcal{A}}_{nm}(\bm{k},\bm{p})P_{L}(p)P_{L}(|\bm{k}-\bm{p}|),$
(A.2) $\displaystyle\equiv I^{m,ijk}_{n,a}(k)+I^{m,ijk}_{n,b}(k).$ (A.3)

for the following index combination :

Index combination
---
$m=1$ | $n=1$ | $ijk=udd$
$m=2$ | $n=1,2$ | $ijk=uud$
$m=3$ | $n=2,3$ | $ijk=uuu$
Table 6:  Index combination for $I^{m,ijk}_{n}$ functions.

That is, the labels $uuu$, $uud$, $\dots$, serve to remind us where each term
comes from. For example, the function $I^{2,uud}_{2}(k)$ comes from a
bispectrum composed of two velocity fields ($u$) and one density field ($d$).
It should not be confused with the function $I^{2,uudd}_{2}(k)$, introduced
below, which is obtained from a 4-point correlation $\langle uudd\rangle$.
Actually, the total function $I^{m}_{n}$ with $m=2,n=2$, is
$I^{2}_{2}(k)=P_{\theta\theta}(k)+I^{2,uud}_{2}(k)+I^{2,uudd}_{2}(k)$.

The usual functions in the TNS model are $A_{nm}$, $\tilde{A}_{nm}$ and
$a_{nm}$ [86], instead of our bold faced ${\text{\bf a}}_{nm}$ and
$\mathcal{\tilde{A}}_{nm}$. But they are actually the same, as is manifestly
under the identifications ${\text{\bf
a}}_{nm}(\bm{k},\bm{p})=r^{-2}a_{nm}(\bm{k},\bm{p})$ and
$\mathcal{\tilde{A}}_{nm}(\bm{k},\bm{p})=r^{-2}\tilde{A}_{nm}(\bm{k},\bm{p})/(1+r^{2}-2rx)^{2}$
(with $r=p/k$ and $x=\hat{\bm{p}}\cdot\hat{\bm{k}}$), and only differ in the
used kernels $G_{2}$ and rates $f(k)/f_{0}$ [64]. Also, the functions
$\mathcal{A}_{nm}$ are unnecessary because
$\mathcal{A}_{nm}(\bm{k},\bm{k}-\bm{p})={\text{\bf a}}_{nm}(\bm{k},\bm{p})$,
and hence the factor “2” in eq. (A.2).

The functions ${\text{\bf a}}_{nm}(\bm{k},\bm{p})$,
$\tilde{\mathcal{A}}_{nm}(\bm{k},\bm{p})$ are given by

$\displaystyle{\text{\bf a}}_{11}(\bm{k},\bm{p})$
$\displaystyle=2\frac{\bm{k}\cdot\bm{p}}{p^{2}}G_{1}(p)F_{2}(-\bm{k},\bm{p})+2\frac{\bm{k}\cdot(\bm{k}-\bm{p})}{|\bm{k}-\bm{p}|^{2}}G_{2}(-\bm{k},\bm{p}),$
(A.4) $\displaystyle\tilde{\mathcal{A}}_{11}(\bm{k},\bm{p})$
$\displaystyle=4\frac{\bm{k}\cdot\bm{p}}{p^{2}}G_{1}(p)F_{2}(\bm{p},\bm{k}-\bm{p}),$
(A.5)

$\displaystyle{\text{\bf a}}_{12}(\bm{k},\bm{p})$
$\displaystyle=\frac{1}{|\bm{k}-\bm{p}|^{2}}\left[\frac{(\bm{k}\cdot\bm{p})^{2}}{p^{2}}-k^{2}\right]G_{1}(p)G_{2}(-\bm{k},\bm{p}),$
(A.6) $\displaystyle\tilde{\mathcal{A}}_{12}(\bm{k},\bm{p})$
$\displaystyle=\frac{1}{|\bm{k}-\bm{p}|^{2}}\left[\frac{(\bm{k}\cdot\bm{p})^{2}}{p^{2}}-k^{2}\right]G_{1}(p)G_{1}(|\bm{k}-\bm{p}|)F_{2}(\bm{p},\bm{k}-\bm{p}),$
(A.7)

$\displaystyle{\text{\bf a}}_{33}(\bm{k},\bm{p})$
$\displaystyle=\frac{1}{|\bm{k}-\bm{p}|^{2}}\left[k^{2}+2\frac{k^{2}}{p^{2}}(\bm{k}\cdot\bm{p})-3\frac{(\bm{k}\cdot\bm{p})^{2}}{p^{2}}\right]G_{1}(k)G_{1}(p)G_{2}(-\bm{k},\bm{p}),$
(A.8) $\displaystyle\tilde{\mathcal{A}}_{33}(\bm{k},\bm{p})$
$\displaystyle=\frac{1}{|\bm{k}-\bm{p}|^{2}}\left[k^{2}+2\frac{k^{2}}{p^{2}}(\bm{k}\cdot\bm{p})-3\frac{(\bm{k}\cdot\bm{p})^{2}}{p^{2}}\right]G_{1}(p)G_{1}(|\bm{k}-\bm{p}|)G_{2}(\bm{p},\bm{k}-\bm{p}),$
(A.9)

$\displaystyle{\text{\bf a}}_{22}(\bm{k},\bm{p})$
$\displaystyle=\frac{{\text{\bf a}}_{33}(\bm{k},\bm{p})}{G_{1}(k)}+{\text{\bf
a}}_{11}(\bm{k},\bm{p})G_{1}(k),$ (A.10)
$\displaystyle\tilde{\mathcal{A}}_{22}(\bm{k},\bm{p})$
$\displaystyle=\frac{1}{|\bm{k}-\bm{p}|^{2}}\left[k^{2}+2\frac{k^{2}}{p^{2}}(\bm{k}\cdot\bm{p})-3\frac{(\bm{k}\cdot\bm{p})^{2}}{p^{2}}\right]G_{1}(p)G_{1}(|\bm{k}-\bm{p}|)F_{2}(\bm{p},\bm{k}-\bm{p})$
$\displaystyle\quad+4\frac{\bm{k}\cdot\bm{p}}{p^{2}}G_{1}(p)G_{2}(\bm{p},\bm{k}-\bm{p}),$
(A.11)

$\displaystyle{\text{\bf a}}_{23}(\bm{k},\bm{p})$ $\displaystyle={\text{\bf
a}}_{12}(\bm{k},\bm{p})G_{1}(k),$ (A.12)
$\displaystyle\tilde{\mathcal{A}}_{23}(\bm{k},\bm{p})$
$\displaystyle=\frac{1}{|\bm{k}-\bm{p}|^{2}}\left[\frac{(\bm{k}\cdot\bm{p})^{2}}{p^{2}}-k^{2}\right]G_{1}(p)G_{1}(|\bm{k}-\bm{p}|)G_{2}(\bm{p},\bm{k}-\bm{p}).$
(A.13)

For EdS kernels, all these functions coincide with those reported in TNS paper
[86], but here we have reduced the number of different functions $A_{mn}$,
$a_{mn}$ and $\tilde{A}_{mn}$ from being 15 to 8, plus two more (${\text{\bf
a}}_{22}$ and ${\text{\bf a}}_{23}$) that are trivially obtained from the
rest: that is, they differ only by factors $G_{1}(k)$ that can be pulled out
of the internal momentum $\bm{p}$ integrals as shown in eqs. (B.38) and
(B.39).

Now, the $D(k,\mu)$ function in eq. (3.42) can be expressed as [65]

$\displaystyle D(k,\mu)$
$\displaystyle=\mu^{2}f^{2}_{0}I^{2,uudd}_{1}(k)+\mu^{4}\big{[}f^{2}_{0}I^{2,uudd}_{2}(k)+f^{3}_{0}I^{3,uuud}_{2}(k)+f^{4}_{0}I^{4,uuuu}_{2}(k)\big{]}$
$\displaystyle\quad+\mu^{6}\big{[}f^{3}_{0}I^{3,uuud}_{3}(k)+f^{4}_{0}I^{4,uuuu}_{3}(k)\big{]}+\mu^{8}f^{4}_{0}I^{4,uuuu}_{4}(k),$
(A.14)

with

$\displaystyle I^{2,uudd}_{n}(k)$
$\displaystyle=I^{2,uudd}_{n,D}(k)-\delta_{n1}k^{2}\sigma^{2}_{vv}P_{\delta\delta}(k),\hskip
28.45274pt(n=1,2)$ (A.15) $\displaystyle I^{3,uuud}_{n}(k)$
$\displaystyle=I^{3,uuud}_{n,D}(k)-\delta_{n2}2k^{2}\sigma^{2}_{vv}P_{\delta\theta}(k),\qquad(n=2,3)$
(A.16) $\displaystyle I^{4,uuuu}_{n}(k)$
$\displaystyle=I^{4,uuuu}_{n,D}(k)-\delta_{n3}k^{2}\sigma^{2}_{vv}P_{\theta\theta}(k),\hskip
25.6073pt(n=2,3,4)$ (A.17)

where the label $D$ means the sum of contributions $B$ and $C$. The velocity
dispersion $\sigma^{2}_{vv}$ is given by

$\sigma^{2}_{vv}\equiv\frac{1}{6\pi^{2}}\int^{\infty}_{0}dp\,P^{L}_{\theta\theta}(p).$
(A.18)

Explicitly, the above $I^{m}_{n}$ functions are

$\displaystyle
I^{2,uudd}_{n,B}(k)=\int_{\bm{p}}\mathcal{B}^{n}_{11}(\bm{k},\bm{p})P^{L}_{\delta\theta}(p)P^{L}_{\delta\theta}(|\bm{k}-\bm{p}|),$
(A.19)

with

$\displaystyle\mathcal{B}^{1}_{11}(\bm{k},\bm{p})$
$\displaystyle=\frac{1}{2|\bm{k}-\bm{p}|^{2}}\left[\frac{(\bm{k}\cdot\bm{p})^{2}}{p^{2}}-k^{2}\right],$
(A.20) $\displaystyle\mathcal{B}^{2}_{11}(\bm{k},\bm{p})$
$\displaystyle=\frac{1}{2|\bm{k}-\bm{p}|^{2}}\left[k^{2}+2\frac{k^{2}}{p^{2}}(\bm{k}\cdot\bm{p})-3\frac{(\bm{k}\cdot\bm{p})^{2}}{p^{2}}\right],$
(A.21)

and

$\displaystyle
I^{2,uudd}_{n,C}(k)=2\int_{\bm{p}}\mathcal{C}^{n}_{11}(\bm{k},\bm{p})P_{L}(p)P^{L}_{\theta\theta}(|\bm{k}-\bm{p}|),$
(A.22)

with

$\displaystyle\mathcal{C}^{1}_{11}(\bm{k},\bm{p})$
$\displaystyle=\frac{1}{4|\bm{k}-\bm{p}|^{4}}\left[k^{2}p^{2}-(\bm{k}\cdot\bm{p})^{2}\right],$
(A.23) $\displaystyle\mathcal{C}^{2}_{11}(\bm{k},\bm{p})$
$\displaystyle=\frac{1}{4|\bm{k}-\bm{p}|^{4}}\left[2k^{4}-k^{2}p^{2}-4k^{2}(\bm{k}\cdot\bm{p})+3(\bm{k}\cdot\bm{p})^{2}\right].$
(A.24)

Notice these expressions have different power spectrum product combinations,
so they can not be written into a single term, as we do in the following for
the rest of the $I^{m}_{n}$ functions for $D(k,\mu)$.

Functions with $m=3$ are

$\displaystyle
I^{3,uuud}_{n,D}(k)=-2\int_{\bm{p}}\mathcal{D}^{n}_{21}(\bm{k},\bm{p})P^{L}_{\delta\theta}(p)P^{L}_{\theta\theta}(|\bm{k}-\bm{p}|),$
(A.25)

with

$\displaystyle\mathcal{D}^{2}_{21}(\bm{k},\bm{p})$
$\displaystyle=\frac{k^{4}}{2|\bm{k}-\bm{p}|^{4}}\left[2-3\frac{(\bm{k}\cdot\bm{p})}{k^{2}}\right]\left[1-\frac{(\bm{k}\cdot\bm{p})^{2}}{k^{2}p^{2}}\right],$
(A.26) $\displaystyle\mathcal{D}^{3}_{21}(\bm{k},\bm{p})$
$\displaystyle=-\frac{k^{4}}{2|\bm{k}-\bm{p}|^{4}}\left[2+\left(\frac{2}{p^{2}}-\frac{3}{k^{2}}\right)(\bm{k}\cdot\bm{p})-6\frac{(\bm{k}\cdot\bm{p})^{2}}{k^{2}p^{2}}+5\frac{(\bm{k}\cdot\bm{p})^{3}}{k^{4}p^{2}}\right].$
(A.27)

Finally, functions with $m=4$ are given by

$\displaystyle
I^{4,uuuu}_{n,D}(k)=\int_{\bm{p}}\mathcal{D}^{n}_{22}(\bm{k},\bm{p})P^{L}_{\theta\theta}(p)P^{L}_{\theta\theta}(|\bm{k}-\bm{p}|),$
(A.28)

with

$\displaystyle\mathcal{D}^{2}_{22}(\bm{k},\bm{p})$
$\displaystyle=\frac{3}{16|\bm{k}-\bm{p}|^{4}}\left[k^{2}-\frac{(\bm{k}\cdot\bm{p})^{2}}{p^{2}}\right]^{2},$
(A.29) $\displaystyle\mathcal{D}^{3}_{22}(\bm{k},\bm{p})$
$\displaystyle=\frac{k^{2}}{8p^{2}|\bm{k}-\bm{p}|^{4}}\left[k^{2}-\frac{(\bm{k}\cdot\bm{p})^{2}}{p^{2}}\right]\left[2k^{2}-3p^{2}-12(\bm{k}\cdot\bm{p})+15\frac{(\bm{k}\cdot\bm{p})^{2}}{k^{2}}\right],$
(A.30) $\displaystyle\mathcal{D}^{4}_{22}(\bm{k},\bm{p})$
$\displaystyle=\frac{k^{4}}{16p^{2}|\bm{k}-\bm{p}|^{4}}\bigg{[}3p^{2}-4k^{2}+24(\bm{k}\cdot\bm{p})+\left(\frac{12}{p^{2}}-\frac{30}{k^{2}}\right)(\bm{k}\cdot\bm{p})^{2}-40\frac{(\bm{k}\cdot\bm{p})^{3}}{k^{2}p^{2}}$
$\displaystyle\quad+35\frac{(\bm{k}\cdot\bm{p})^{4}}{k^{4}p^{2}}\bigg{]}.$
(A.31)

We have reduced the $m=3,4$ contributions to a single term expression. Then,
for these cases, the $B$ and $C$ contributions boil down to only one matrix
multiplication for each $n$ value when applying the FFTLog formalism.

To complete the full framework, we need to introduce the bias into functions
$A^{\text{TNS}}$ and $D$. These are obtained by substituting $f_{0}\rightarrow
f_{0}/b_{1}$ and weighting the functions with powers of $b_{1}$, namely

$\displaystyle A^{\text{TNS}}(k,\mu;f_{0})$
$\displaystyle\,\longrightarrow\,b_{1}^{3}A^{\text{TNS}}(k,\mu;f_{0}/b_{1}),$
(A.32) $\displaystyle D(k,\mu;f_{0})$
$\displaystyle\,\longrightarrow\,b_{1}^{4}D(k,\mu;f_{0}/b_{1}).$ (A.33)

The biasing for the $D$ function is exact, because $D$ is constructed from
linear fields only. Instead, the function $A^{\text{TNS}}$ is also biased by
the non-linear $\delta^{2}_{cb}$ and $s^{2}$ operators. The biased expression
of $A^{\text{TNS}}(k,\mu)$, includes also $b_{2}$ and $b_{s^{2}}$ but the
resulting corrections are subdominant, even smaller in amplitude than the
other contributions carrying these bias parameters, so we neglect them. The
complete expression including second order linear and tidal biases can be
found in Appendix A.1 of [64].

##  Appendix B FFTLog contributions

In this appendix we present explicit expressions for FFTLog contributions,
including their matrices or vectors and UV or IR corrections if necessary. We
employ two different values for the bias parameter $\nu$. For the biasing
power spectra $P_{b_{1}b_{2}},\dots$ we use $\nu=-1.51$, while for the rest of
the functions we use $\nu=-0.1$.

In the following, remind that $\nu_{1}$ and $\nu_{2}$ should not be confused
with the bias $\nu$, but these are defined through eqs. (4.16).

### Overdensity and velocity power spectra

We start with the leading non-linear contributions (3.31) – (3.36), whose
approximations are

$\displaystyle
P^{22}_{\delta\delta}(k)=k^{3}\displaystyle\sum\limits_{m_{1},m_{2}}c_{m_{1}}k^{-2\nu_{1}}M_{22,\delta\delta}(\nu_{1},\nu_{2})\,c_{m_{2}}k^{-2\nu_{2}},$
(B.1) $\displaystyle
M_{22,\delta\delta}(\nu_{1},\nu_{2})=\Big{[}\nu_{1}\nu_{2}\left(98\nu^{2}_{12}-14\nu_{12}+36\right)-91\nu^{2}_{12}+3\nu_{12}+58\Big{]}$
$\displaystyle\quad\qquad\qquad\qquad\times\frac{\left(\frac{3}{2}-\nu_{12}\right)\left(\frac{1}{2}-\nu_{12}\right)}{196\nu_{1}(1+\nu_{1})\left(\frac{1}{2}-\nu_{1}\right)\nu_{2}(1+\nu_{2})\left(\frac{1}{2}-\nu_{2}\right)}I(\nu_{1},\nu_{2}).$
(B.2)

$\displaystyle
P^{22}_{\delta\theta}(k)=2k^{3}\displaystyle\sum\limits_{m_{1},m_{2}}c^{f}_{m_{1}}k^{-2\nu_{1}}M^{f_{\bm{p}}}_{22,\delta\theta}(\nu_{1},\nu_{2})\,c_{m_{2}}k^{-2\nu_{2}},$
(B.3) $\displaystyle
M^{f_{\bm{p}}}_{22,\delta\theta}(\nu_{1},\nu_{2})=\Big{[}7(7\nu_{1}+3)\nu_{2}^{2}+(7\nu_{1}(7\nu_{1}-1)-38)\nu_{2}-21\nu_{1}-23\Big{]}$
$\displaystyle\quad\qquad\qquad\qquad\times\frac{(2\nu_{12}-3)(2\nu_{12}-1)}{196\nu_{1}(\nu_{1}+1)\nu_{2}(\nu_{2}+1)(2\nu_{2}-1)}I(\nu_{1},\nu_{2}).$
(B.4)

$\displaystyle
P^{22}_{\theta\theta}(k)=k^{3}\displaystyle\sum\limits_{m_{1},m_{2}}c^{ff}_{m_{1}}k^{-2\nu_{1}}M^{f_{\bm{p}}f_{\bm{p}}}_{22,\theta\theta}(\nu_{1},\nu_{2})\,c_{m_{2}}k^{-2\nu_{2}}$
$\displaystyle\quad\qquad+k^{3}\displaystyle\sum\limits_{m_{1},m_{2}}c^{f}_{m_{1}}k^{-2\nu_{1}}M^{f_{\bm{k}-\bm{p}}f_{\bm{p}}}_{22,\theta\theta}(\nu_{1},\nu_{2})\,c^{f}_{m_{2}}k^{-2\nu_{2}},$
(B.5) $\displaystyle
M^{f_{\bm{p}}f_{\bm{p}}}_{22,\theta\theta}(\nu_{1},\nu_{2})=\Big{[}98\nu_{1}^{3}\nu_{2}+7\nu_{1}^{2}\big{(}2\nu_{2}(7\nu_{2}-8)+1\big{)}-\nu_{1}\big{(}2\nu_{2}(7\nu_{2}+17)+53\big{)}$
$\displaystyle\quad\qquad\qquad\qquad-12(1-2\nu_{2})^{2}\Big{]}\frac{2\nu_{12}-3}{98\nu_{1}(\nu_{1}+1)\nu_{2}(\nu_{2}+1)(2\nu_{2}-1)}I(\nu_{1},\nu_{2}),$
(B.6) $\displaystyle
M^{f_{\bm{k}-\bm{p}}f_{\bm{p}}}_{22,\theta\theta}(\nu_{1},\nu_{2})=\Big{[}7\nu_{1}^{2}(7\nu_{2}+3)+\nu_{1}\big{(}7\nu_{2}(7\nu_{2}-1)-10\big{)}+\nu_{2}(21\nu_{2}-10)-37\Big{]}$
$\displaystyle\quad\qquad\qquad\qquad\times\frac{2\nu_{12}-3}{98\nu_{1}(\nu_{1}+1)\nu_{2}(\nu_{2}+1)}I(\nu_{1},\nu_{2}).$
(B.7)

$\displaystyle
P^{13}_{\delta\delta}(k)=k^{3}P_{L}(k)\displaystyle\sum\limits_{m_{1}}c_{m_{1}}k^{-2\nu_{1}}M_{13,\delta\delta}(\nu_{1}),$
(B.8) $\displaystyle
M_{13,\delta\delta}(\nu_{1})=\frac{1+9\nu_{1}}{4}\frac{\tan(\nu_{1}\pi)}{28\pi(\nu_{1}+1)\nu_{1}(\nu_{1}-1)(\nu_{1}-2)(\nu_{1}-3)}.$
(B.9)

$\displaystyle
P^{13}_{\delta\theta}(k)=\frac{1}{2}k^{3}P_{L}(k)\displaystyle\sum\limits_{m_{1}}\left(\frac{f(k)}{f_{0}}c_{m_{1}}+c^{f}_{m_{1}}\right)k^{-2\nu_{1}}M^{f_{\bm{k}}}_{13,\delta\theta}(\nu_{1}),$
(B.10) $\displaystyle
M^{f_{\bm{k}}}_{13,\delta\theta}(\nu_{1})=\frac{9\nu_{1}-7}{4}\frac{\tan(\nu_{1}\pi)}{28\pi(\nu_{1}+1)\nu_{1}(\nu_{1}-1)(\nu_{1}-2)(\nu_{1}-3)}.$
(B.11)

$\displaystyle