# Clustering in Massive Neutrino Cosmologies via Eulerian Perturbation Theory

Alejandro Aviles  Arka Banerjee  Gustavo Niz  Zachary Slepian

###### Abstract

We introduce an Eulerian Perturbation Theory to study the clustering of
tracers for cosmologies in the presence of massive neutrinos. Our approach is
based on mapping recently-obtained Lagrangian Perturbation Theory results to
the Eulerian framework. We add Effective Field Theory counterterms, IR-
resummations and a biasing scheme to compute the one-loop redshift-space power
spectrum. To assess our predictions, we compare the power spectrum multipoles
against synthetic halo catalogues from the Quijote simulations, finding
excellent agreement on scales $k\lesssim 0.25\,h\,\text{Mpc}^{-1}$. One can
obtain the same fitting accuracy using higher wave-numbers, but then the
theory fails to give a correct estimation of the linear bias parameter. We
further discuss the implications for the tree-level bispectrum. Finally,
calculating loop corrections is computationally costly, hence we derive an
accurate approximation wherein we retain only the main features of the
kernels, as produced by changes to the growth rate. As a result, we show how
FFTLog methods can be used to further accelerate the loop computations with
these reduced kernels.

##  1 Introduction

The standard cosmological model, which commonly assumes a massless neutrino
component, has successfully explained most observations, from the exquisite
final Cosmic Microwave Background (CMB) data of the Planck satellite [1] to
the large scale structure in the recent analysis of galaxy survey catalog
release of SDSS eBOSS DR16 [2] or the Y3 DES results (e.g. [3]). However, the
observed oscillations of atmospheric and solar neutrinos show that their
flavour and mass eigenstates are not the same [4, 5, 6], constraining the sum
of neutrino masses to $\sum_{i}m_{\nu_{i}}\gtrsim 0.06\,(0.10)$ eV for a
normal (inverted) hierarchy [7]. Although particle physics experiments have
reached a sub-eV sensitivity on the lightest neutrino mass (e.g. the KATRIN
experiment upper bound $m_{\nu}<0.8\,\text{eV}$ (0.90 c.l.) using tritium
single beta decay [8, 9]), the absolute mass scale, the mass hierarchy and
whether neutrinos are Majorana or Dirac fermions are still open questions.

Indeed, cosmological experiments have indirectly placed the strongest
constraints on the sum of the neutrino masses, due to the non-negligible
impact of these hot dark matter-like components on the overall matter
evolution in the Universe. The Planck temperature-temperature power spectrum
analysis constrains $\sum_{i}m_{\nu,i}\lesssim 0.24$ eV, which is narrowed
down to $\sum_{i}m_{\nu,i}\lesssim 0.12$ eV when combined with Baryonic
Acoustic Oscillations (BAO) data [10]. Moreover, a “full-shape” analysis of
the redshift-space power spectrum using the low-$z$ North Galactic Cap (NGC)
data selection of the SDSS BOSS DR12 galaxy sample, without combining it with
other datasets, yields $\sum_{i}m_{\nu,i}<1.17$ eV (0.95 c.l.) [11] which is
still better than current particle physics experiments. Another promising
avenue using higher-redshift tracers is the 1D power spectrum of the Lyman
alpha forest, which places a bound of $\sum_{i}m_{\nu,i}<0.71$ eV (0.95 c.l.).
This bound shrinks to $\sum_{i}m_{\nu,i}<0.09$ eV (0.95 c.l.) when combined
with BAO data and the CMB temperature, polarization and lensing [12]. It
should be stressed that all the neutrino mass bounds from cosmology have
strong degeneracies with other cosmological parameters.

Primordial neutrinos decouple from baryons at early epochs
($\sim\\!1\,\text{MeV}$), when the weak interaction rate falls below the
Hubble rate and is hence no longer able to maintain them in equilibrium with
the primeval plasma, which at that time is composed mainly of protons,
electrons and photons. Hence, neutrinos decouple while they are still
relativistic, and become non-relativistic at late times, contributing to the
total matter abundance in this late epoch (e.g. [13]). Despite this present-
day non-relativistic nature in the bulk, neutrinos do have a large velocity
dispersion, inherited from their phase-space distribution at freeze-out, which
in turn prevents them from clustering below the free-streaming scale, roughly
defined through a Jeans-like mechanism (e.g. [14]).

Another relevant scale in neutrino cosmology is the maximum comoving distance
that neutrinos can travel over the history of the Universe. For realistic
masses, this distance is about $1\,h^{-1}\,\text{Gpc}$. Above this length-
scale, neutrinos behave essentially as Cold Dark Matter (CDM); while below it,
their density fluctuations are suppressed at some stage of evolution. A
backreaction effect, which also damps the gravitational potentials driving the
growth of structure in the Universe, leads to a suppression of CDM and baryon
fluctuations as well. This implies that both the total matter and the CDM
power spectra become damped below this scale. For a comprehensive review of
the impact of neutrino physics in cosmology see e.g. [15].

Furthermore, the matter velocity fields are also affected in the presence of
massive neutrinos, mainly because the free-streaming scale induces a scale-
dependent growth of overdensities, even in the linear regime. Hence the large-
scale velocity field receives additional contributions compared to that in a
cosmology with massless neutrinos. In particular, the Kaiser boost also
becomes scale-dependent, because of the scale dependence in the logarithmic
growth rate. Since the scales produced by the neutrino mass as discussed above
are within the linear and quasi-linear regimes, Perturbation Theory (PT) is a
valuable tool, complementary to $N$-body simulations, to study the clustering
induced by massive neutrinos. In redshift-space the neutrino effects become
not only more evident but also more degenerate with other PT parameters. In
particular, we see this when looking at the broadband shape of the power
spectrum using a multipole expansion or wedges. In this context, higher-order
statistics may help to break these degeneracies (see e.g. [16, 17, 18]).

As the forthcoming generation of higher-precision galaxy surveys, such as DESI
[19],111desi.lbl.gov/ Euclid [20]222sci.esa.int/web/euclid or the Vera C.
Rubin Observatory [21]333lsst.org/ are expected to reach enough precision to
measure the absolute mass scale of the neutrinos, and potentially also its
hierarchy, exhaustive analytical and semi-analytical methods should be
accordingly developed for cosmological parameter inference. This was recently
done using the BOSS DR12 data in [11, 22, 23] with the use of EPT and the
Effective Field Theory of Large-Scale Structure (EFT)

In the current work, we construct an Eulerian Perturbation Theory (EPT) for
the clustering of tracers in redshift space in the presence of massive
neutrinos, which is self-consistent and well-behaved in both the UV and the
IR. Our modelling is based on the Lagrangian Perturbation Theory (LPT)
constructed recently by two of us in [24]. More precisely, the velocity and
density kernels are obtained by mapping the LPT kernels of [24] to the
Eulerian framework. Furthermore, the mapping method heavily relies on a
previous study [25], where the redshift-space power spectrum with generalized
kernels, beyond Einstein-de Sitter (EdS), was obtained for representative
modified gravity models. The latter work is, in turn, a generalization of the
velocity moments expansion approach, in the form presented in [26]. In
addition, we include the non-linear effects of large-scales in terms of bulk
flows via the IR-resummations techniques of [27, 28, 29], which are not
captured by the EPT and serve to damp the BAO features in the power spectrum.
We also add EFT counterterms to model the small-scales physics, which are out
of reach from the PT, together with the damping of the density and velocity
fields along the line-of-sight due to Redshift Space Distortions (RSD).
Finally, we use the biasing scheme of [30], slightly modified in [25], to
account for additional scales. In summary and to our knowledge, this is the
first work that contains all the above ingredients in redshift-space for
massive neutrino cosmologies.

One important disadvantage of the method we have described is that computing
the loop corrections is computationally very costly; which makes it highly
unlikely that our theory, as it stands, can be used for efficient cosmological
parameter estimation in the near future. The reason for this slowness is that
we have to solve several ordinary differential equations at each volume
element in the loop integrals. However, we show that the main differences of
the perturbative kernels from those in a $\Lambda$CDM cosmology come from the
scale-dependent rate. Hence, we simplify them in a manner where the growth
rate effects are maintained, and obtain simple analytical expressions that
allow the use of FFTLog methods [31, 32, 33], dramatically reducing the
computational time.

At the same time as analytic methods have been improved, significant progress
has been made in capturing the effects of massive neutrinos on structure
formation down to fully nonlinear scales using $N$-body simulations. Various
numerical techniques have been developed over the years to account for
neutrino effects in these simulations [34, 35, 36, 37, 38, 39, 40, 41, 42].
These techniques range from treating neutrinos as an Eulerian fluid on a grid
coupled to the full nonlinear gravitational potential, adding neutrinos as a
separate set of $N$-body particles, and hybrid techniques.

Despite the wide variety of techniques used in these simulations, the results
on the quasi-linear scales of interest in this paper are in agreement with
each other at the percent level. These simulation results serve two purposes
for various analytical and semi-analytical methods—they can be used to
calibrate the analytic methods, e.g. to help fix values of certain free
parameters in the model, and to validate the results from the semi-analytical
methods over a wide range of the parameter space of interest. Therefore, we
compare the multipoles of the power spectrum obtained from our theory against
the dark matter halos obtained from the
Quijote444github.com/franciscovillaescusa/Quijote-simulations suite of
simulations [43] at redshifts $z=0.5$ and $1.0$, and for neutrinos with total
mass $M_{\nu}=\sum_{i}m_{\nu,i}=$ 0, 0.1, 0.2 and 0.4 eV, where the masses are
equally distributed among the three mass eigenstates. We find excellent
agreement up to scales $k\simeq 0.25\,h\,\text{Mpc}^{-1}$ for the monopole and
quadrupole, as is expected generally for EFT theory. The agreement can be
extended to smaller scales at the cost of losing accuracy at small wave-
numbers, and probably such small scales are not in the range where the
broadband shape of the power spectrum provides cosmological information

We would like to end this introductory section by referring the reader to
other relevant works which include the effects of neutrinos using PT,
certainly a topic much less explored than the $\Lambda$CDM case. We start by
mentioning the early work of [44] which describes the linear theory only, but
provides several meaningful physical insights and approximate analytical
formulae for the linear growth and transfer functions. To our knowledge, the
first study of nonlinear PT was done in [45] with the use of Einstein-de
Sitter (EdS) kernels, with the subsequently refined work of [46, 47, 48, 49,
50, 51, 52, 53, 54, 55, 56, 57], which included a more complete description of
non-linearities. We will highlight some of these references in what follows.

The rest of this paper is organized as follows. In §2 we describe our
treatment of linear and non-linear neutrino density fluctuations. §3 is
devoted to obtaining the EPT kernels by mapping the LPT ones obtained in [24].
In §4 we present the theory for the 1-dimensional spectra of density-density
and velocity-velocity fluctuations, as well as the density-velocity cross
spectra; we also present the biasing scheme. In §5, we outline the redshift-
space description of the power spectrum and add EFT counterterms and IR
resummations. We confront our predictions with $N$-body simulations in §6. In
§7 we show how the FFTLog method can be used to accelerate the loop integrals.
We summarize and conclude in §8. In the appendix A we perform a direct
calculation of the second-order Eulerian kernels, complementary to the mapping
performed in §3.

##  2 Neutrino density

Our starting point is the Poisson equation in the presence of massive
neutrinos and at late times when relativistic components can be neglected,
which can be written as

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

where $\Phi$ is the gravitational potential and we have defined the following
quantities:

$A(k,t)\equiv A_{0}\big{[}f_{cb}+f_{\nu}\alpha(k,t)\big{]},\qquad
A_{0}\equiv\frac{3}{2}\Omega_{m}H^{2},\qquad\alpha(k,t)\equiv\frac{\delta_{\nu}(\bm{k},t)}{\delta_{cb}(\bm{k},t)}.$
(2.2)

Therefore, $A(k,t)$ acts as a scale-dependent gravitational strength and
$A_{0}$ is its limit as $k\rightarrow 0$. $\delta_{cb}$ and $\delta_{\nu}$ are
the density fluctuations of the combined CDM + baryon fluid and the massive
neutrinos, respectively, while $f_{cb}=\Omega_{cb}/\Omega_{m}$ and
$f_{\nu}=\Omega_{\nu}/\Omega_{m}=1-f_{cb}$ are the relative abundances with
respect to the total matter fluid $m=cb+\nu$. We are working on sufficiently
large scales that baryonic effects are negligible, and so we treat the baryons
as part of the cold dark matter component; see, e.g., [58]. We note that at
least in linear theory $\delta_{cb}$ and $\delta_{\nu}$ have the same phase,
inherited from the curvature potential at the end of inflation, hence
$\alpha(k,t)$ reduces to a function of $k$ only.

The key approximation in our analysis is to take

$\alpha(k,t)\approx\frac{T_{\nu}(k,t)}{T_{cb}(k,t)},$ (2.3)

where $T_{\nu}$ and $T_{cb}$ are the transfer functions connecting
overdensities imprinted at the end of inflation with linear overdensities at
time $t$. This approximation has been tested using simulations in [24],
finding a very good agreement up to well inside the non-linear regime $k\sim
1\,h\,\text{Mpc}^{-1}$, for neutrino masses up to $M_{\nu}=0.4\,\text{eV}$. In
more detail, the total matter power spectrum can be written as the sum

$P_{mm}(k)=f_{cb}^{2}P_{cb}(k)+2f_{\nu}f_{cb}P_{cb,\nu}(k)+f_{\nu}^{2}P_{\nu}(k),$

and using the approximation of eq. (2.3), we can substitute

$\displaystyle P_{cb,\nu}(k)=\alpha(k)P_{cb}(k)\quad$
$\displaystyle\longrightarrow\quad\frac{T_{\nu}(k)}{T_{cb}(k)}P_{cb}(k),$
$\displaystyle P_{\nu}(k)=\alpha^{2}(k)P_{cb}(k)\quad$
$\displaystyle\longrightarrow\quad\left(\frac{T_{\nu}(k)}{T_{cb}(k)}\right)^{2}P_{cb}(k),$

to obtain

$P_{mm}(k)=\left[f_{cb}^{2}+2f_{\nu}f_{cb}\frac{T_{\nu}(k)}{T_{cb}(k)}+f_{\nu}^{2}\left(\frac{T_{\nu}(k)}{T_{cb}(k)}\right)^{2}\right]P_{cb}(k).$

The analysis of [24] extracts the $P_{mm}^{\text{sims}}$ and
$P_{cb}^{\text{sims}}$ power spectra directly from the Quijote simulations
[43], and constructs an approximation of $P_{mm}^{\text{approx}}$ by first,
substituting in the above equation $P_{cb}\rightarrow P_{cb}^{\text{sims}}$
and second, using the linear transfer functions computed with the code CAMB. A
comparison to the simulated matter power spectrum yields
$|P_{mm}^{\text{approx}}/P_{mm}^{\text{sims}}-1|<0.001$. This result is
consistent with previous work [49, 50], which find a $0.1\%$ error in the
real-space power spectrum when considering the approximation
$\delta_{\nu}=(\delta^{(1)}_{\nu}/\delta^{(1)}_{cb})\delta_{cb}$ for non-
linear neutrino overdensities; in agreement with our approximation [eq. (2.3)]
for adiabatic perturbations, which we assume throughout this work.

These equations show the main aspects of neutrino clustering. Neutrinos become
non-relativistic when their temperature drops below their mass; this happens
at redshift $1+z_{nr,i}\approx 1894\,m_{\nu,i}/\text{eV}$ for each mass
eigenstate $m_{\nu,i}$. Afterwards, they behave as hot dark matter with
abundance

$\Omega_{\nu 0}=\frac{M_{\nu}}{93.14\,h^{2}\,\text{eV}},\qquad
M_{\nu}=\sum_{i=1}^{3}m_{\nu,i},$ (2.4)

that do not cluster below their free-streaming scales, which is given for each
neutrino species by [14]

$k_{\text{FS,i}}(z)\approx
0.0908\frac{H(z)}{(1+z)^{2}}\left(\frac{m_{\nu,i}}{\text{0.1
eV}}\right)\,h\,\text{Mpc}^{-1}.$ (2.5)

For equal mass neutrinos, $m_{\nu,i}=M_{\nu}/3$, the different free-streaming
scales reduce to a single value, that we simply call $k_{\text{FS}}$. $N$-body
simulations with massive neutrinos usually make this assumption, as we also
consider throughout this work, even though our results do not depend on this
choice. At very large scales, $k\ll k_{\text{FS}}$, neutrinos behave as CDM
($T_{\nu}\sim T_{cb}(k)$) and the Poisson equation is that of $\Lambda$CDM. In
the opposite limit ($k\gg k_{\text{FS}}$), one gets that $T_{\nu}\ll
T_{cb}(k)$ and the source of the gravitational potential becomes proportional
to $f_{cb}$, which ends up in the damping of fluctuations at scales smaller
than the free-streaming scale.

To find the evolution of $cb$ fields in PT, we must supplement the Poisson
equation with the continuity and Euler equations and solve them order by order
to obtain the kernels $F_{n}$ and $G_{n}$, leading to the $n$-ordered density
fluctuation and velocity fields

$\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),$
(2.6) $\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),$
(2.7)

where we have adopted the shorthand 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}.$ (2.8)

The $\theta_{cb}$ field is the divergence of the peculiar velocity $\bm{v}$;
which more precisely we define as

$\theta_{cb}(\bm{k},t)=-\frac{i\bm{k}\cdot\bm{v}}{aHf_{0}},$ (2.9)

where $f$ is the logarithmic derivative of the linear growth rate $D_{+}$ with
respect to scale factor $a$,

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

and

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

is the large scale limit of the growth rate, which is not sensitive to the
neutrino mass because on these scales neutrinos behave as cold dark matter.
One should notice that both, the linear growth function and the growth rate,
become scale-dependent due to the additional scale introduced by the free-
streaming. Moreover, the linear growth function $D_{+}$ is the solution to the
equation

$\ddot{D}(k,t)+2H\dot{D}-A(k,t)D=0,$ (2.12)

with the appropriate initial conditions to isolate the growing solution, which
are chosen well inside the matter dominated epoch by the approximated
analytical formula of [44].

In figure 1, we show plots for the scale-dependent linear growth function and
growth rate of the $cb$ fluid, at redshifts $z=0$ and $z=1$, for degenerate
neutrinos with total mass $M_{\nu}=$ 0.1, 0.2 and 0.4 eV, which for
$\Omega_{m}=0.3175$ and $h=0.6711$ correspond to
$f_{\nu}=0.0075,\,0.0150,\,0.0300$, respectively. The linear transfer
functions $T_{\nu}(k,z)$ and $T_{cb}(k,z)$ are computed with the code
CAMB555https://camb.info/. [59]. Alternatively, these function may be obtained
to a very high accuracy with the analytical formulae of [44].

Figure 1: Linear growth function $D_{+}(k,z)$ (left panel) and growth rate
$f(k,z)$ (right panel) for neutrino masses $M_{\nu}=0.1,0.2,0.4$ eV, divided
by the massless neutrino cases. We consider the redshifts $z=1$ (solid lines)
and $z=0$ (dashed lines). In the left panel, the vertical dashed lines show
the free-streaming scales at $z=0$ computed from eq. (2.5). In the right
panel, the dashed horizontal lines show the small scales limit $f(k\gg
k_{\text{FS}})\approx(1-3/5f_{\nu})f^{M_{\nu}=0}$.

Using the linear growth rate $D_{+}$, the linear density fields evolve as

$\delta^{(1)}_{cb}(\bm{k},t)=\frac{D_{+}(k,t)}{D_{+}(k,t_{0})}\delta^{(1)}_{cb}(\bm{k},t_{0}),$
(2.13)

with $t_{0}$ the present time. However, note that this decomposition is
somewhat arbitrary since a full factorization into scale-dependent and time-
dependent pieces is not possible. On the other hand, the velocity field and
the density fluctuations are related by

$\theta_{cb}^{(1)}(\bm{k},t)=\frac{f(k,t)}{f_{0}(t)}\delta_{cb}^{(1)}(\bm{k},t),$
(2.14)

thus coinciding at large scales but resulting in suppressed velocities, by a
factor of $f(k)/f_{0}$, below the free-streaming scale. From above’s equation
and eq. (2.7) we obtain the first order kernels

$F_{1}(\bm{k},t)=1,\qquad\text{and}\qquad
G_{1}(\bm{k},t)=\frac{f(k,t)}{f_{0}(t)}.$ (2.15)

Since the relation between the linear velocity and density fields is
multiplicative in Fourier space, they become non-locally related in
configuration space. However, the $f(k)$ growth rate is close to a step
function, sharply interpolating in $\log(k)$, at a scale $\sim k_{\text{FS}}$,
between $f_{0}$ and $f_{\infty}=f(k\rightarrow\infty)$. Hence, the non-local
effect arising from convoluting the inverse Fourier transform of $f(k)$ with
the configuration space $cb$ overdensity is small. In practice, the density
fields are always smoothed over certain scale; and therefore, if the width of
support of $\tilde{f}$ (or the physical region where it is reasonable
different from zero) turns out to be smaller than the smoothing scale of the
density fields, the relation between the $\theta^{(1)}_{cb}$ and
$\delta^{(1)}_{cb}$ becomes effectively local.

Finding the higher order kernels is quite cumbersome using fluid equations and
eqs. (2.6) and (2.7). However, in this work we take a different approach,
which maps the LPT kernels found in [24] to the required EPT kernels. Before
presenting this mapping (in section 3) let us review the basic ingredients of
the LPT formalism. Moreover, in appendix A, the reader may find the second
order kernels $F_{2}$ and $G_{2}$ directly.

###  2.1 Lagrangian evolution

In a Lagrangian description the observer follows trajectories of the $cb$
fluid elements that have initial coordinates $\bm{q}$ and final, Eulerian
coordinates

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

with $\bm{\Psi}$ the Lagrangian displacement of the $cb$ field, which obeys
the evolution 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}}.$
(2.17)

The non-linear neutrino fluctuations in the Lagrangian space become [24]

$\displaystyle\tilde{\delta}_{\nu}(\bm{k})$
$\displaystyle=\alpha(k)\tilde{\delta}_{cb}(\bm{k})+\frac{1}{A_{0}f_{\nu}}\underset{\bm{k}_{12}=\bm{k}}{\int}\mathcal{K}^{\text{FL$\Psi$}}_{li}(\bm{k}_{1},\bm{k}_{2})\Psi_{l}(\bm{k}_{1})\Psi_{i}(\bm{k}_{2})$
$\displaystyle\quad+\frac{1}{A_{0}f_{\nu}}\underset{\bm{k}_{123}=\bm{k}}{\int}\mathcal{K}^{\text{FL$\Psi$}}_{lij}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})\Psi_{l}(\bm{k}_{1})\Psi_{i}(\bm{k}_{2})\Psi_{j}(\bm{k}_{3})+\cdots,$
(2.18)

where we have omit to write the time dependence. In general, we do not write a
tilde over quantities transformed to Fourier space, since it usually does not
lead to confusions. However, in the equation above we do write a tilde over
the $\nu$ and $cb$ overdensities to point out that they are Fourier transforms
of Eulerian quantities but with respect to the Lagrangian coordinates, that
is, $\tilde{\delta}_{cb,\nu}(\bm{k})=\int
d^{3}q\,e^{-i\bm{k}\cdot\bm{q}}\delta_{cb,\nu}(\bm{x})=\int
d^{3}q\,e^{-i\bm{k}\cdot\bm{q}}\delta_{cb,\nu}(\bm{q}+\bm{\Psi})$. In
contrast, the function $\alpha(k)$ (without a tilde) is indeed the $x$-Fourier
transform of the neutrino to $cb$ fluid overdensities ratio in the Eulerian
framework, and not
$\tilde{\alpha}(k)=\tilde{\delta}_{\nu}(\bm{k})/\tilde{\delta}_{cb}(\bm{k})$.
As a result, eq. (2.1) shows the non-linear nature of the relation between
neutrino and $cb$ overdensities in Lagrangian space, regardless if one takes
$\alpha$ to be equal to the linear transfer functions ratio; for further
details and a detailed derivation we refer the reader to [24]. Actually, this
non-linear relation is encoded in the kernels $\mathcal{K}^{\text{FL}\Psi}$,
namely

$\displaystyle\mathcal{K}^{\text{FL$\Psi$}}_{li}(\bm{k}_{1},\bm{k}_{2})$
$\displaystyle=\big{[}A(k_{1})-A(k_{12})\big{]}k_{1}^{l}k_{1}^{i},$ (2.19)
$\displaystyle\mathcal{K}^{\text{FL$\Psi$}}_{lij}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$
$\displaystyle=-i\big{[}A(k_{123})-A(k_{1})\big{]}k_{1}^{l}k_{1}^{i}k_{1}^{j}$
$\displaystyle\quad+i\big{[}A(k_{123})-A(k_{13})\big{]}k_{13}^{i}\left[k_{1}^{j}k_{1}^{l}+\frac{1}{2}k_{1}^{l}k_{3}^{j}+\frac{1}{2}k_{1}^{j}k_{3}^{l}\right],$
(2.20)

which becomes exactly zero when the function $A(k)$ (or $\alpha(k)$) is scale-
independent. This is expected since it would mean no new scale was introduced
into the theory, as in the limit of zero neutrino masses where the relation of
eq. (2.1) becomes linear. These kernels were introduced and named “frame-
lagging” in [60], in a slightly different mathematical form, to correct the
“lag” between spatial derivatives taken in the Eulerian and Lagrangian frames.

Now, the LPT formalism expands the Lagrangian displacements
$\bm{\Psi}=\bm{\Psi}^{(1)}+\bm{\Psi}^{(2)}+\bm{\Psi}^{(3)}+\cdots$ and solves
eq. (2.17) iteratively, to find

$\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).$
(2.21)

The LPT kernels $\bm{L}^{(n)}$ up to third order for massive neutrino
cosmologies are found in eqs. (3.4), (3.8) and (3.18) of [24]. It is common to
assume that the linear fields are purely longitudinal in these computations,
so that the transverse pieces, of the form
$\bm{\Psi}^{\text{trans}}(\bm{k})=\bm{k}\times\bm{T}$, appear at terms with
orders $n>2$ [61]. However, in one-loop 2-point statistics, these transverse
parts are contracted with external wave-vectors,
$\bm{k}\cdot\bm{\Psi}^{\text{trans}}$, so that they do not contribute.
Therefore, without loss of generality we can use the divergence of the
Lagrangian fields rather than the Lagrangian displacements. Doing this will
ease the computations shown in the next section. To order $n$, we have

$-i\bm{k}\cdot\bm{\Psi}^{(n)}(\bm{k},t)=\frac{1}{n!}\underset{\bm{k}_{1\cdots
n}=\bm{k}}{\int}C^{(n)}\Gamma_{n}(\bm{k}_{1},\cdots,\bm{k}_{n};t)\delta_{cb}^{(1)}(\bm{k}_{1},t)\cdots\delta_{cb}^{(1)}(\bm{k}_{n},t),$
(2.22)

where $C^{(n)}\Gamma_{n}$ is the $n$th order kernel for
$-i\bm{k}\cdot\bm{\Psi}$, and the $C^{(n)}$ are constant coefficients that we
fixed, as in [62], to the values $C^{(1)}=C^{(3)}=1$, and $C^{(2)}=3/7$.
Moreover, eqs. (2.21) and (2.22) imply

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

Equivalently, for the time derivative of $\bm{\Psi}$ we find at order $n$ that

$-i\bm{k}\cdot\dot{\bm{\Psi}}^{(n)}(\bm{k},t)=\frac{1}{n!}nHf_{0}\underset{\bm{k}_{1\cdots
n}=\bm{k}}{\int}C^{(n)}\Gamma^{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),$
(2.24)

with

$\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).$
(2.25)

In the particular case of EdS kernels, one gets that
$\Gamma^{f}_{n}=\Gamma_{n}$, recovering the known relation
$\dot{\bm{\Psi}}^{(n)}=nHf\bm{\Psi}^{(n)}$ [63].

##  3 From Lagrangian to Eulerian kernels

In this section we map the LPT kernels of [24] to the EPT counterparts, to
obtain the $F_{2,3}$ and $G_{2,3}$ kernels that we use in this analysis. We
start with the Eulerian-space overdensity written in terms of Lagrangian
displacements:

$\delta_{cb}(\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{]}.$
(3.1)

To bring the above equation into a complete EPT description, in the form given
in eq. (2.6), we first expand the Lagrangian displacement out of the
exponential, as

$\displaystyle\delta_{cb}(\bm{k})$
$\displaystyle=\sum_{m=1}^{\infty}\frac{(-i)^{m}}{m!}\int
d^{3}q\;e^{-i\bm{k}\cdot\bm{q}}\big{[}\bm{k}\cdot\bm{\Psi}(\bm{q})\big{]}^{m}$
(3.2) $\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.3)

where in the second equality we have Fourier-transformed each Lagrangian
displacement,
$\Psi_{i}(\bm{q})=\int_{\bm{p}}e^{i\bm{p}\cdot\bm{k}}\Psi_{i}(\bm{p})$, and
performed the $\bm{q}$ integral analytically, which yields to a Dirac delta
function ensuring momentum conservation $\bm{p}_{1}+\cdots+\bm{p}_{m}=\bm{k}$.
The next step is to expand the Lagrangian displacements perturbatively as
$\Psi_{i}(\bm{p})=\Psi_{i}^{(1)}(\bm{p})+\cdots=[\bm{p}\cdot\bm{\Psi}^{(1)}(\bm{p})]p_{i}/p^{2}+\cdots$;
then use eq. (2.22), and finally relate the $\Gamma_{n}$ kernels with the
$F_{n}$ kernels using eq. (2.6). This procedure yields to the following
results for $F_{2}$ and $F_{3}$

$\displaystyle F_{2}(\bm{k}_{1},\bm{k}_{2})$
$\displaystyle=\frac{3}{14}\Gamma_{2}(\bm{k}_{1},\bm{k}_{2})+\frac{1}{2}\frac{(\bm{k}_{12}\cdot\bm{k}_{1})(\bm{k}_{12}\cdot\bm{k}_{2})}{k_{1}^{2}k_{2}^{2}},$
(3.4) $\displaystyle F_{3}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$
$\displaystyle=\frac{1}{6}\left[\Gamma_{3}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})+\frac{(\bm{k}_{123}\cdot\bm{k}_{1})(\bm{k}_{123}\cdot\bm{k}_{2})(\bm{k}_{123}\cdot\bm{k}_{3})}{k_{1}^{2}k_{2}^{2}k_{3}^{2}}\right]$
$\displaystyle\quad+\frac{1}{14}\left[\frac{(\bm{k}_{123}\cdot\bm{k}_{12})(\bm{k}_{123}\cdot\bm{k}_{3})}{k_{12}^{2}k_{3}^{2}}\Gamma_{2}(\bm{k}_{1},\bm{k}_{2})+\text{(cyclic)}\right].$
(3.5)

The computation of the $G_{n}$ kernels requires more work. The peculiar
velocity field is given by $\bm{v}=a\dot{\bm{\Psi}}$, or its divergence

$J(\bm{q},t)\frac{1}{a}\frac{\partial\,}{\partial
x^{i}}v^{i}(\bm{x},t)=\dot{\Psi}_{i,i}+\Psi_{j,j}\dot{\Psi}_{i,i}-\Psi_{i,j}\dot{\Psi}_{i,j}+\frac{1}{2}\epsilon_{ikp}\epsilon_{jqr}\Psi_{k,q}\Psi_{p,r}\dot{\Psi}_{i,j},$
(3.6)

where a colon means derivative with respect to Lagrangian coordinates, and
$J=\det(J_{ij})$ with $J_{ij}=\partial x_{i}/\partial q^{j}$ is the Jacobian
matrix of the coordinate transformation of eq. (2.16). In order to arrive to
eq. (3.6) we use $\partial_{\bm{x}^{i}}=(J^{-1})_{ji}\partial_{\bm{q}^{j}}$
and the relation
$J\,(J^{-1})_{ji}=\delta_{ij}+(\delta_{ij}\delta_{ab}-\delta_{ia}\delta_{jb})\Psi_{a,b}+\frac{1}{2}\epsilon_{ikp}\epsilon_{jqr}\Psi_{k,q}\Psi_{p,r}$.
With these identities at hand, we compute

$\displaystyle\theta_{cb}(\bm{k})$ $\displaystyle=-\frac{1}{aHf_{0}}\int
d^{3}x\;e^{-i\bm{k}\cdot\bm{x}}\frac{\partial v_{i}}{\partial
x^{i}}=-\frac{1}{Hf_{0}}\int
d^{3}q\;e^{-i\bm{k}\cdot\bm{q}}e^{-i\bm{k}\cdot\bm{\Psi}(\bm{q},t)}J(\bm{q},t)\frac{1}{a}\frac{\partial
v_{i}}{\partial x^{i}}$
$\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}\Bigg{[}\dot{\Psi}_{i,i}+(\delta_{ij}\delta_{ab}-\delta_{ia}\delta_{jb})\Psi_{a,b}\dot{\Psi}_{i,j}$
$\displaystyle\quad\qquad+\frac{1}{2}\epsilon_{ikp}\epsilon_{jqr}\Psi_{k,q}\Psi_{p,r}\dot{\Psi}_{i,j}\Bigg{]},$
(3.7)

where we use $d^{3}x=J(\bm{q},t)d^{3}q$ in the second equality, and Taylor-
expand the exponential $e^{-i\bm{k}\cdot\bm{\Psi}(\bm{q},t)}$ together with
eq. (3.6) in the third equality. Finally, by using eqs. (2.22) and (2.24) we
can relate the $\Gamma_{n}$ and $\Gamma^{f}_{n}$ LPT kernels to the $G_{n}$
EPT kernels, to obtain the following expression for $G_{2}$

$\displaystyle G_{2}(\bm{k}_{1},\bm{k}_{2})$
$\displaystyle=\frac{3}{7}\Gamma^{f}_{2}(\bm{k}_{1},\bm{k}_{2})+\frac{(\bm{k}_{1}\cdot\bm{k}_{2})^{2}}{k_{1}^{2}k_{2}^{2}}\frac{f_{1}+f_{2}}{2f_{0}}+\frac{1}{2}\frac{\bm{k}_{1}\cdot\bm{k}_{2}}{k_{1}k_{2}}\left(\frac{k_{2}}{k_{1}}\frac{f_{2}}{f_{0}}+\frac{k_{1}}{k_{2}}\frac{f_{1}}{f_{0}}\right),$
(3.8)

where $f_{1}=f(k_{1})$, $f_{2}=f(k_{2})$. The third-order kernel
$G_{3}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$ is much more cumbersome and its
expression is not very illuminating; while it can be found in eq. (C.22) of
[25]. However, it simplifies considerably when it is symmetrized over its
three arguments and evaluated at the wave-vectors $\bm{k}_{1}=\bm{k}$,
$\bm{k}_{3}=-\bm{k}_{2}=\bm{p}$, which correspond to the configuration
appearing in the one-loop computations. By doing this, one arrives at

$\displaystyle
G_{3}(\bm{k},-\bm{p},\bm{p})=\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}{6}\frac{(\bm{k}\cdot\bm{p})^{2}}{p^{4}}\frac{f(k)}{f_{0}}-\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].$
(3.9)

The $F_{3}$ kernel for the same configuration reduces to

$\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.10)

It is easy to show that $F_{3}$ and $G_{3}$ given by eqs. (3.10) and (3)
reduce to the standard kernels for EdS evolution with the use of

$\displaystyle\Gamma_{2}^{\text{EdS}}(\bm{k},-\bm{p})=\Gamma_{2}^{f,\,\text{EdS}}(\bm{k},-\bm{p})$
$\displaystyle=1-\frac{(\bm{k}\cdot\bm{p})^{2}}{k^{2}p^{2}},$ (3.11)
$\displaystyle\Gamma_{3}^{\text{EdS}}(\bm{k},-\bm{p},\bm{p})=\Gamma_{2}^{f,\,\text{EdS}}(\bm{k},-\bm{p},\bm{p})$
$\displaystyle=\frac{5}{21}\frac{k^{2}}{|\bm{k}-\bm{p}|^{2}}\left(1-\frac{(\bm{k}\cdot\bm{p})^{2}}{k^{2}p^{2}}\right)^{2}+\,(\bm{p}\rightarrow-\bm{p}).$
(3.12)

###  3.1 Second order kernels and the tree-level bispectrum

We now take a closer look to the second order kernels, which are explicitly
given in eqs. (A.22) and (A.23). It is instructive to write them as

$\displaystyle F_{2}(k_{1},k_{2},x)$
$\displaystyle=\left(\frac{2}{3}+\frac{3\mathcal{A}-\mathcal{B}}{14}\right)\mathcal{L}_{0}(x)+\frac{1}{2}\left(\frac{k_{2}}{k_{1}}+\frac{k_{1}}{k_{2}}\right)\mathcal{L}_{1}(x)+\left(\frac{1}{3}-\frac{\mathcal{B}}{7}\right)\mathcal{L}_{2}(x),$
(3.13) $\displaystyle G_{2}(k_{1},k_{2},x)$
$\displaystyle=\left(\frac{3(3\mathcal{A}-\mathcal{B})(f_{1}+f_{2})+3(3\dot{\mathcal{A}}-\dot{\mathcal{B}})/H+7(f_{1}+f_{2})}{42f_{0}}\right)\mathcal{L}_{0}(x)$
$\displaystyle\quad+\frac{1}{2}\left(\frac{f_{2}}{f_{0}}\frac{k_{2}}{k_{1}}+\frac{f_{1}}{f_{0}}\frac{k_{1}}{k_{2}}\right)\mathcal{L}_{1}(x)+\left(\frac{(3\mathcal{B}-7)(f_{1}+f_{2})+\dot{\mathcal{B}}/H}{21f_{0}}\right)\mathcal{L}_{2}(x),$
(3.14)

where $x\equiv\hat{k}_{1}\cdot\hat{k}_{2}$ and $\mathcal{L}_{L}(x)$ denotes
the Legendre polynomial of degree $L$. The functions $\mathcal{A}$ and
$\mathcal{B}$ depend on the wave-vector magnitudes $k_{1}$, $k_{2}$ and the
cosine of their corresponding opening angle, $x$. These functions show up
because of the scale- and time-dependence of the growth rate $f$; or more
generally, because the relation $f^{2}=\Omega_{m}(t)$ does not hold.
Therefore, one cannot solve analytically the differential equations for the
kernels, as it happens in the simplest case of EdS where both the growth rate
and the matter abundance are unity [64]. As such, these functions are
solutions to second-order linear differential equations, which can be found in
[24] [eqs. (3.9)-(3.11)]; and also in eqs. (A.25) and (A.26) in appendix A. In
$\Lambda$CDM, the functions $\mathcal{A}$ and $\mathcal{B}$ are scale-
independent and evolve very slowly with time, becoming exactly unity for EdS.
However, despite its appearance, notice that the above equations are not
complete Legendre decompositions since $\mathcal{A}$ and $\mathcal{B}$ are
themselves functions of $x$. For example, the dipole coefficients ($L=1$) are

$\displaystyle F_{2}^{L=1}(k_{1},k_{2})$
$\displaystyle=\frac{1}{2}\left(\frac{k_{2}}{k_{1}}+\frac{k_{1}}{k_{2}}\right)+[\Delta
F_{2}(k_{1},k_{2},x)]^{L=1},$ (3.15) $\displaystyle G_{2}^{L=1}(k_{1},k_{2})$
$\displaystyle=\frac{1}{2}\left(\frac{k_{2}}{k_{1}}\frac{f(k_{2})}{f_{0}}+\frac{k_{1}}{k_{2}}\frac{f(k_{1})}{f_{0}}\right)+[\Delta
G_{2}(k_{1},k_{2},x)]^{L=1},$ (3.16)

where the first pieces in the above equations are the so-called gradient
terms, while the pieces on squared brackets are given by

$\displaystyle[\Delta F_{2}(k_{1},k_{2},x)]^{L=1}$
$\displaystyle=\frac{3}{2}\int_{-1}^{1}dx\left[F_{2}(k_{1},k_{2},x)-\frac{x}{2}\left(\frac{k_{2}}{k_{1}}+\frac{k_{1}}{k_{2}}\right)\right]\mathcal{L}_{1}(x),$
(3.17) $\displaystyle[\Delta G_{2}(k_{1},k_{2},x)]^{L=1}$
$\displaystyle=\frac{3}{2}\int_{-1}^{1}dx\left[G_{2}(k_{1},k_{2},x)-\frac{x}{2}\left(\frac{k_{2}}{k_{1}}\frac{f(k_{2})}{f_{0}}+\frac{k_{1}}{k_{2}}\frac{f(k_{1})}{f_{0}}\right)\right]\mathcal{L}_{1}(x).$
(3.18)

Actually, it turns out that $\Delta G_{2}$ and $\Delta F_{2}$ are negligible
contributions, because

$\displaystyle[\Delta F_{2}(k_{1},k_{2},x)]^{L=1}$
$\displaystyle\ll\frac{1}{2}\left(\frac{k_{2}}{k_{1}}+\frac{k_{1}}{k_{2}}\right),$
(3.19) $\displaystyle[\Delta G_{2}(k_{1},k_{2},x)]^{L=1}$
$\displaystyle\ll\frac{1}{2}\left(\frac{k_{2}}{k_{1}}\frac{f(k_{2})}{f_{0}}+\frac{k_{1}}{k_{2}}\frac{f(k_{1})}{f_{0}}\right).$
(3.20)

Figure 2: Ratios of dipole terms $[\Delta F_{2}]^{L=1}$ and $[\Delta
G_{2}]^{L=1}$, defined in eqs. (3.17) and (3.18), over the gradient terms of
$F_{2}$ and $G_{2}$. These plots show numerically that eqs. (3.19) and (3.20)
hold true. We use $M_{\nu}=0.4\,\text{eV}$ for which the ratios are the
largest, reaching at most $2\%$ for the $G_{2}$ kernel within the shown $k$
range. The evaluation redshift is $z=0.5$.

In figure 2, we show density plots of the ratios of the LHS to the RHS sides
in above inequalities, showing that these inequalities indeed hold. In
conclusion, the dipoles for the $F_{2}$ and $G_{2}$ kernels are approximated,
to high accuracy, by

$\displaystyle F_{2}^{L=1}(k_{1},k_{2})$
$\displaystyle=\frac{1}{2}\left(\frac{k_{2}}{k_{1}}+\frac{k_{1}}{k_{2}}\right),$
(3.21) $\displaystyle G_{2}^{L=1}(k_{1},k_{2})$
$\displaystyle=\frac{1}{2}\left(\frac{k_{2}}{k_{1}}\frac{f(k_{2})}{f_{0}}+\frac{k_{1}}{k_{2}}\frac{f(k_{1})}{f_{0}}\right).$
(3.22)

Moreover, when both momenta are larger than the free streaming wavenumber
($k_{1},k_{2}>k_{\text{FS}}\simeq 0.1\,h\,\text{Mpc}^{-1}$), the last
expression reduces to

$G_{2}^{L=1}(k_{1},k_{2})\simeq\frac{1}{2}\left(\frac{k_{2}}{k_{1}}+\frac{k_{1}}{k_{2}}\right)\left(1-\frac{3}{5}f_{\nu}\right),$
(3.23)

where we used the approximation [44, 46] (see also figure 1)

$f(k>k_{\text{FS}})=\left(1-\frac{3}{5}f_{\nu}\right)f_{0}.$ (3.24)

As a result, we obtain a dipole piece in agreement with [18], namely

$\displaystyle[F_{2}^{\text{EdS}}(\bm{k}_{1},\bm{k}_{2})-F_{2}(\bm{k}_{1},\bm{k}_{2})]^{L=1}$
$\displaystyle=0,$ (3.25)
$\displaystyle[G_{2}^{\text{EdS}}(\bm{k}_{1},\bm{k}_{2})-G_{2}(\bm{k}_{1},\bm{k}_{2})]^{L=1}$
$\displaystyle=-\frac{3}{10}\left(\frac{k_{2}}{k_{1}}+\frac{k_{1}}{k_{2}}\right)f_{\nu},$
(3.26)

which as we stated above, is valid for $k_{1},k_{2}>k_{\text{FS}}$. It is
worth noticing that this result was obtained in [18] using the kernels of
[46], which are different from ours since that PT assumes the neutrino density
field is linear; in particular, the gradient terms do not take the form of
eqs. (3.28) and (3.29) below. Despite these differences, the $G_{2}$ and
$F_{2}$ kernels of [46] coincide with ours at two important limits. First, at
$k\ll k_{\text{FS}}$ since in this limit neutrino densities are linear and
behave as CDM; and second, at small scales, $k\gg k_{\text{FS}}$, where the
neutrino overdensities are essentially zero regardless if these are considered
linear or not. Because of the latter, the results in eqs. (3.25) and (3.26)
coincide in both approaches.

On small scales within the free-streaming scale, [18] applied the expression
of [46] for the linear growth rate, and derived leading-order in $f_{\nu}$
corrections to the $F_{2}$ and $G_{2}$ kernels to construct the tree-level
bispectrum. Importantly, they found a change to the dipole of the $G_{2}$
kernel, which dipole is at the pre-cyclic level protected from any changes due
to biasing. This is a signature of how neutrinos truly modify velocities in a
unique manner.

The reader may also notice that the gradient terms on the EPT kernels are
fixed by the advection of density and velocity fields due to the Lagrangian
displacement. In the case of the $F_{2}$ gradient term, it arises from
expanding $\delta(\bm{x}+\bm{\Psi})$ to lowest order, which results in

$\delta^{(1)}(\bm{x}+\bm{\Psi})-\delta^{(1)}(\bm{x})=\Psi^{(1)}_{i}(\bm{x})\partial_{i}\delta^{(1)}(\bm{x})=\partial_{i}\big{[}\nabla^{-2}\delta^{(1)}(\bm{x})\big{]}\partial_{i}\delta^{(1)}(\bm{x}).$
(3.27)

The previous expression is part of the second order perturbation,
$\delta^{(2)}(\bm{x})$, and is the same contribution that one gets in the
massless neutrino $\Lambda$CDM model. At this level the relation with the
gradient terms of $F_{2}(\bm{k}_{1},\bm{k}_{2})$ is not evident, but once we
move to Fourier space and symmetrize over the momenta $\bm{k}_{1}$ and
$\bm{k}_{2}$, the previous expansion becomes

$\frac{x}{2}\left(\frac{k_{2}}{k_{1}}+\frac{k_{1}}{k_{2}}\right),$ (3.28)

which clearly corresponds to the gradients terms in dipole coefficient of eq.
(3.21). (Do not confuse the Eulerian coordinate $\bm{x}$ with the angle’s
cosine $x$.) On the other hand, for the $G_{2}$ one cannot take the $F_{2}$
result unmodified, because of the relation between the linear overdensities
and the velocities (see eq. (2.14)). Therefore, a similar expansion, but for
$\theta^{(1)}(\bm{x}+\bm{\Psi})-\theta^{(1)}(\bm{x})$, in Fourier space now
leads to

$\frac{x}{2}\left(\frac{f(k_{2})}{f_{0}}\frac{k_{2}}{k_{1}}+\frac{f(k_{1})}{f_{0}}\frac{k_{1}}{k_{2}}\right),$
(3.29)

which again is the gradient term in eq. (3.22) times $\mathcal{L}_{1}$. As a
result, the advection of matter fields differs from the massless neutrino case
because of the scale-dependent growth of the linear velocity divergence.

We end this section by discussing other approaches in the literature. In the
seminal article [45], the nonlinear $cb$ fluctuations are obtained neglecting
the neutrino density fluctuation and using the EdS kernels, but with the
linear power spectrum of the $cb$ fluid as input for the loop computations.
Another recurrent approach is to approximate
$\delta_{\nu}=\delta^{(1)}_{\nu}$, obtained from a Boltzmann code, and use the
linear fluctuation as an external source to compute the non-linear $cb$
overdensities [47, 46, 55]. However, the latter approach violates momentum
conservation, as proved in [52]; in particular, one does not obtain the
relation $(P_{\text{NL}}-P_{\text{L}})/P_{\text{L}}\propto k^{2}$ at large
scales. Instead, the prescription of [52] evolves the non-linear neutrino
density fields by truncating the Boltzmann hierarchy at the Euler equation,
and approximating the second moment of the phase-space distribution function
to be proportional to an effective sound speed times the density contrast.
This fluid approach finds justification in the work of [48, 14, 65] (see also
appendix C of [66]). In more recent works [56, 67], the authors use non-linear
perturbations around a Fermi-Dirac massive neutrino distribution and solve the
coupled Boltzmann and CDM density field equations iteratively, by expanding in
powers of $f_{\nu}$ and keeping only the linear order. In addition, it is
worth mentioning a paper published very recently in the subject [57], where
neutrino’s field fluctuations are described by a hybrid Boltzmann-two fluid
model, which then are used to compute the matter real-space power spectrum up
to two loops. Finally, in a similar manner as the renormalization group time
flow scheme used in [49, 50], or the work in [54] which treats $f_{\nu}$ as an
expansion parameter, we improve the description of neutrino evolution by using
the non-linear quantity $\delta_{cb}(\delta_{\nu}^{(1)}/\delta_{cb}^{(1)})$ as
a proxy for $\delta_{\nu}$, instead of consider it a linear quantity or simply
neglect it. However, in our approach, we further find exact kernels, up to a
few functions that are solutions to linear second-order ordinary differential
equations, and which nonetheless are very close to unity. The simple form of
our kernels allows us to draw rapid conclusions of the theory and to develop
efficient algorithms for the computation of loop corrections; as we do in
section 7. To the best of our knowledge, the approach presented here is the
only one that explicitly recovers the gradient terms in eqs. (3.28) and (3.29)
that are fixed by the advection of matter fields.

##  4 1-dimensional spectra

A direct consequence of the results presented so far is the one-loop auto and
cross power spectra for the dentsity and velocity fields. The resulting
expressions for the three possible combinations are

$\displaystyle P^{\text{1-loop}}_{cb,\delta\delta}(k)$
$\displaystyle=P^{L}_{cb,\delta\delta}(k)+2\int_{\bm{p}}\big{[}F_{2}(\bm{p},\bm{k}-\bm{p})\big{]}^{2}P_{L}(p)P_{L}(|\bm{k}-\bm{p}|)$
$\displaystyle\quad+6\int_{\bm{p}}F_{3}(\bm{k},-\bm{p},\bm{p})P_{L}(k)P_{L}(p),$
(4.1) $\displaystyle P^{\text{1-loop}}_{cb,\delta\theta}(k)$
$\displaystyle=P^{L}_{cb,\delta\theta}(k)+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}|)$
$\displaystyle\quad+3\int_{\bm{p}}\left[F_{3}(\bm{k},-\bm{p},\bm{p})\frac{f(k)}{f_{0}}+G_{3}(\bm{k},-\bm{p},\bm{p})\right]P_{L}(k)P_{L}(p),$
(4.2) $\displaystyle P^{\text{1-loop}}_{cb,\theta\theta}(k)$
$\displaystyle=P^{L}_{cb,\theta\theta}(k)+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\quad+6\int_{\bm{p}}G_{3}(\bm{k},-\bm{p},\bm{p})\frac{f(k)}{f_{0}}P_{L}(k)P_{L}(p),$
(4.3)

where, as before, we have adopted the compact integral notation
$\int_{\bm{p}}=\int d^{3}p/(2\pi)^{3}$, and the linear contributions
$P^{L}_{cb,\dots}$ are given by

$\displaystyle P_{L}(k)\equiv
P^{L}_{cb,\delta\delta}(k)=\langle\delta^{(1)}_{cb}(\bm{k})\delta^{(1)}_{cb}(\bm{k}^{\prime})\rangle^{\prime},$
(4.4)

and

$\displaystyle P^{L}_{cb,\delta\theta}(k)=\frac{f(k)}{f_{0}}P_{L}(k),\qquad
P^{L}_{cb,\theta\theta}(k)=\left(\frac{f(k)}{f_{0}}\right)^{2}P_{L}(k).$ (4.5)

Using these results, we depict the one-loop spectra $P_{cb,\delta\delta}(k)$,
$P_{cb,\delta\theta}(k)$ and $P_{cb,\theta\theta}(k)$ in figure 3, at an
intermediate redshift of $z=0.5$ for three degenerate sum of neutrino masses
($M_{\nu}=0.4,0.2,0.1$ eV). It is possible to appreciate how the one-loop
result departs from the linear spectra at small scales; this happens more
noticeable for $\theta$ field than for $\delta$. The departure from linear
theory is more pronounced for more massive neutrinos, and when velocity fields
are involved, because these are affected by the scale dependence of both
linear function and the growth rate, whereas the functions $\mathcal{A}$ and
$\mathcal{B}$ and their third order counterparts play a subdominant role.
Moreover, there is increment in power before one gets the expected suppression
with respect to the linear contribution. Actually, in order to observe the
differences of the neutrino masses on short scales, the cosmological
parameters are the same as the cosmologies used in section 6, but instead of
keeping fixed $\sigma_{8}$ we anchor all the models to have the same
primordial scalar amplitude $A_{s}=2.13\times 10^{-9}$ at a pivotal scale
$k=0.05\,\text{Mpc}^{-1}$. By choosing this normalization, the power spectra
for each massive neutrino case coincide at large scales.

Figure 3: $cb$ power spectra in the presence of neutrinos with mass
$M_{\nu}=0.1,0.2,0.4\,\text{eV}$ over massless neutrinos $cb$ spectra at
redshift $z=0.5$. Left panel shows $P_{\delta\delta}$, middle panel
$P_{\delta\theta}$, and right panel $P_{\theta\theta}$. We plot the linear
(dashed lines) and one-loop (solid lines) perturbative results computed from
eqs. (4.1)-(4.3). The spectra are normalized to have the same primordial
amplitude $A_{s}$.

###  4.1 Biasing the cb field

Large scale bias is primarily sensitive to the clustering of CDM and baryons,
with neutrinos playing a subdominant role because their large velocity
dispersions prevent them to cluster at the relevant scales and redshifts for
halos and galaxies formation. The consequence is that tracers are biased
objects of the $cb$ field, and not of the total matter field [68, 69, 70, 71].
Secondary effects, that we do not consider here, arise due to the critical
overdensity for collapse dependence on the neutrino mass, affecting the halo
abundance and its response to changes in the background density (and hence the
biases) [72, 73]. However, these effects are subdominant and in the following
we assume that the biasing of the $cb$ field is universal and constant at
large scales.666Actually [24] found a preference for a non-zero small
curvature bias for $M_{\nu}=0.4\,\text{eV}$, but consistency with zero for
$M_{\nu}=0.1$ and $0.2\,\text{eV}$.

We employ the biasing scheme of [30], recently generalized to cosmologies that
pose additional scales in [25], which assumes that tracer velocities follow
exactly the $cb$ velocities, while tracer density fluctuations $\delta$ are
biased by

$\displaystyle\delta(\bm{x})$
$\displaystyle=c_{\delta}\delta_{cb}+\frac{1}{2}c_{\delta^{2}}\delta_{cb}^{2}+c_{s^{2}}s^{2}+\frac{1}{6}c_{\delta^{3}}\delta_{cb}^{3}+\frac{1}{2}c_{\delta
s^{2}}\delta_{cb}s^{2}+c_{\psi}\psi+c_{st}st+\frac{1}{6}c_{s^{3}}s^{3}$
$\displaystyle\quad+\text{stochastic terms},$ (4.6)

where $s^{2}(\bm{x})=s_{ij}s_{ij}$, $s^{3}(\bm{x})=s_{ij}s_{jk}s_{ki}$ and

$s_{ij}(\bm{x})=\left(\nabla_{i}\nabla_{j}\nabla^{-2}-\frac{1}{3}\delta_{ij}\right)\delta_{cb}(\bm{x}).$
(4.7)

A useful operator called $\eta(\bm{x})$ [30], defined here through its Fourier
representation

$\eta(\bm{k})=\theta_{cb}(\bm{k})-\frac{f(k)}{f_{0}}\delta_{cb}(\bm{k}),$
(4.8)

manifestly vanishes at first order by means of eq. (2.14). With $\eta$ in
hand, we can construct the following bias operators

$\displaystyle t_{ij}(\bm{k})$
$\displaystyle=\left(\frac{k_{i}k_{j}}{k^{2}}-\frac{1}{3}\delta_{ij}\right)\eta,$
(4.9) $\displaystyle\psi(\bm{k})$
$\displaystyle=\eta-\frac{f(k)}{f_{0}}\left[\frac{2}{7}s^{2}-\frac{4}{21}\delta_{cb}^{2}\right],$
(4.10)

from which we also obtained the third order operator $st$ as $st\equiv
s_{ij}t_{ij}$. It is important to stress that the $\psi$ bias is third order
in EdS, but has a second order contribution in the $\Lambda$CDM model with
massless neutrinos. Since this second order piece is degenerate with the
$s^{2}$ and $\delta^{2}$ operators, $\psi$ can be thought to be a pure third
order quantity in $\Lambda$CDM [74]. However, in the presence of massive
neutrinos, the second order contribution to $\psi$ cannot be absorbed into
other operators, but because it is quite small, in the following, we will
consider it as a third order operator.

In the massless neutrino case, the third order bias contributions to the one-
loop power spectrum (corresponding to the third order bias operators $st$,
$\delta s^{2}$ and $\psi$) differ by a constant that can be absorbed by the
stochastic noise. Hence, they can be collected into a single component given
by $b_{3nl}\sigma^{2}_{3}(k)P_{L}(k)$, where $b_{3nl}$ is a non-local bias
parameter of third order and $\sigma_{3}(k)$ takes the form [30, 75]

$\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],$
(4.11)

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}.$
(4.12)

In contrast, this is no longer possible when additional scales enter into the
theory, as in our case the neutrinos’ mass scale. Indeed all the contributions
have different shapes and one would have to deal with three of them. However,
in [25] it is shown that adding two bias operators $\nabla^{2}\delta$ and
$\nabla^{2}\theta$ suffices to absorb a significant part of the differences in
the third order bias contributions which are inherent to a non-EdS evolution.
This renormalization-like procedure is not strictly necessary because one can
keep the different biasing shape contributions in the analysis. However, for
fitting purposes it is convenient to deal with just one third-order bias
parameter, $b_{3nl}$, instead of having three of them that are closely,
although not completely, degenerate. Actually, this is the only place where we
use the curvature biasing, since we do not use it to construct the scale
dependent linear bias or to fit the data. As a matter of fact, curvature
biases are degenerate with the EFT counterterms that will be introduced later,
and hence terms of the form $k^{2}P_{L}(k)$ do enter into the theory.

Collecting all the bias contributions and after renormalization, one obtains
the tracer spectra to one-loop [75]

$\displaystyle P_{\delta\delta}(k)$
$\displaystyle=b_{1}^{2}P^{\text{1-loop}}_{cb,\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}_{cb,\delta\delta}(k),$
(4.13) $\displaystyle P_{\delta\theta}(k)$
$\displaystyle=b_{1}P^{\text{1-loop}}_{cb,\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),$
(4.14) $\displaystyle P_{\theta\theta}(k)$
$\displaystyle=P^{\text{1-loop}}_{cb,\theta\theta}(k),$ (4.15)

where the one-loop auto and cross spectra
$P^{\text{1-loop}}_{cb,\delta\delta}(k)$,
$P^{\text{1-loop}}_{cb,\delta\theta}(k)$ and
$P^{\text{1-loop}}_{cb,\theta\theta}(k)$ are given by eqs. (4.1), (4.2) and
(4.3), while the remaining pieces are defined by

$\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}|),$
(4.16) $\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}|),$
(4.17) $\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{]},$
(4.18) $\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],$
(4.19) $\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]$
(4.20) $\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}|),$
(4.21) $\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}|).$
(4.22)

Notice that these expressions take the same form as in $\Lambda$CDM, and the
only differences rely in the relations between the bare and the renormalized
biases [25].

One may think of reducing the number of bias parameters by using co-evolution
theory [76, 77, 75], which mathematically implies that

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

These relations assume an EdS evolution, so in the present context of massive
neutrinos they are not entirely satisfied. Yet, they only exhibit trends,
since its derivation assume that the initial tracer fluctuations are expanded
in densities, while tidal and non-local contributions are subsequently
generated by the non-linear gravitational evolution. Actually the last
statement is not completely true, and indeed a non-vanishing Lagrangian tidal
bias has been measured in simulations [78]. However, we believe that for the
statistical level present in our halo catalogues, in principle, the use of the
co-evolution relations (4.23) should be sufficient for parameter estimation.
We will discuss further on this when comparing the theory to the synthetic
data, but before doing so let us study the effect introduced by RSD.

##  5 Redshift-space power spectrum

Peculiar velocities distort the observed positions of objects in the sky
through the Doppler effect. For non-relativistic particles, the effect is on
the direction of observation, such that in the distant-observer approximation,
the map between the true $\bm{x}$ and the observed $\bm{s}$ positions is

$\bm{s}=\bm{x}+\bm{u},$ (5.1)

where the line-of-sight velocity

$\bm{u}=\hat{\bm{n}}\frac{\hat{\bm{n}}\cdot\bm{v}}{aH},$ (5.2)

is constructed from the peculiar velocity, $\bm{v}=a\dot{\bm{x}}$, and the
angular direction to the sample, $\hat{\bm{n}}$. Conservation of number of
objects against the coordinate transformation (5.1) gives the redshift-space
power spectrum [79]

$(2\pi)^{3}\delta_{\text{D}}(\bm{k})+P_{s}(\bm{k})=\int
d^{3}xe^{-i\bm{k}\cdot\bm{x}}\Big{\langle}\big{(}1+\delta(\bm{x}_{1})\big{)}\big{(}1+\delta(\bm{x}_{2})\big{)}e^{-i\bm{k}\cdot\Delta\bm{u}}\Big{\rangle},$
(5.3)

with $\bm{x}=\bm{x}_{2}-\bm{x}_{1}$ and
$\Delta\bm{u}=\bm{u}(\bm{x}_{2})-\bm{u}(\bm{x}_{1})$. By Taylor expanding the
exponential inside the correlation of above’s equation, the power spectrum
becomes [26]

$\displaystyle(2\pi)^{3}\delta_{\text{D}}(\bm{k})+P_{s}(\bm{k})$
$\displaystyle=\sum_{m=0}^{\infty}\frac{(-i)^{m}}{m!}(k\mu)^{m}\tilde{\Xi}^{(m)}_{\hat{\bm{n}}}(\bm{k}),$
(5.4)

where $\mu=\hat{\bm{k}}\cdot\hat{\bm{n}}$, and
$\tilde{\Xi}^{(m)}_{\hat{\bm{n}}}$ are the Fourier transforms of the pairwise
velocity moments along the line-of-sight,

$\displaystyle\Xi^{(m)}_{\hat{\bm{n}}}(\bm{x})$
$\displaystyle\equiv\langle\big{(}1+\delta(\bm{x}_{1})\big{)}\big{(}1+\delta(\bm{x}_{2})\big{)}\big{(}u(\bm{x}_{2})-u(\bm{x}_{1})\big{)}^{m}\rangle,$
(5.5)

with $u=\hat{\bm{n}}\cdot\bm{u}=|\bm{u}|$. Notice that the difference of
peculiar velocities of tracers located at $\bm{x}_{2}$ and $\bm{x}_{1}$
projected along the line-of-sight is given by
$(aH)\hat{\bm{n}}\cdot\Delta\bm{u}=\hat{\bm{n}}\cdot(\bm{v}(\bm{x}_{2})-\bm{v}(\bm{x}_{1}))$.
To linear order in perturbation, one truncates the sum in eq. (5.4) at $m=2$,
to obtain the (linear) Kaiser power spectrum [80]

$P^{K}_{s}(k,\mu)=\left(b_{1}+f(k)\mu^{2}\right)^{2}P_{L}(k),$ (5.6)

which has the same structure as in the massless neutrinos case, but now the
Kaiser boost becomes scale-dependent due to the linear growth rate.

###  5.1 One-loop power spectrum

We generalize the Kaiser spectrum by cutting the sum in eq. (5.4) again at
$m=2$, but keeping terms up to one-loop, obtaining the non-linear Kaiser power
spectrum

$P^{K,\text{NL}}_{s}(k,\mu)=P_{\delta\delta}(k)+f_{0}\mu^{2}P_{\delta\theta}(k)+f_{0}^{2}\mu^{4}P_{\theta\theta}(k),$
(5.7)

with $P_{\delta\delta}$, $P_{\delta\theta}$ and $P_{\theta\theta}$ given by
eqs. (4.1), (4.2) and (4.3). However, this expression does not contain all
one-loop corrections. In fact, the moment $\Xi^{(m)}_{\hat{\bm{n}}}$ involves
correlations of at least $m$ velocity fields and density fluctuations;
therefore, the complete expression for the one-loop EPT power spectrum is
obtained by truncating the sum at $m=4$, yielding [81, 25]

$P^{\text{EPT}}_{s}(k,\mu)=\sum_{m=0}^{4}\frac{(-i)^{m}}{m!}(k\mu)^{m}\tilde{\Xi}^{(m)}_{\hat{\bm{n}}}(\bm{k})=P^{K,\text{NL}}_{s}(k,\mu)+A^{\text{TNS}}(k,\mu)+D(k,\mu),$
(5.8)

while higher order moments are pure two or higher loops contributions. The
functions $A^{\text{TNS}}(k,\mu)$ and $D(k,\mu)$ are [25, 81]

$\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})\,,$
(5.9) $\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{\\}},$
(5.10)

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}}$, and

$\displaystyle B_{\sigma}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$
$\displaystyle=\left\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]\right\rangle^{\prime},$
(5.11) $\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{]}.$
(5.12)

The decomposition in the second equality of eq. (5.8) makes sense since
$P^{K,\text{NL}}$ is constructed out of correlations of only two fields,
either overdensities or velocities, $A$ from three fields, and $D$ from four
fields.

This model is similar to TNS [81], though the latter does not contain the term
$P^{K}_{s}(|\bm{k}-\bm{p}|,\mu_{\bm{k}-\bm{p}})$ in the squared brackets of
eq. (5.10), which is necessary to have a well behaved theory, free of both IR
and UV divergences once assuming typical linear power spectra [25].
Furthermore, the TNS model have a phenomenological Gaussian suppression factor
$\exp[k^{2}\mu^{2}f^{2}\sigma^{2}_{v}]$, providing the damping of the power
spectrum along the line-of-sight direction. This term can be interpreted also
as a partial factorization of the Fingers-of-God (FoG) and Kaiser boost
effects.

Another key aspect to consider is that several of the functions involved in
the line-of-sight moments $\Xi^{(m)}_{\hat{\bm{n}}}$ contain non-vanishing
zero-lag correlators sensible to small-scales physics, so there is a need to
add counterterms to these expressions. Specifically, these come from

$\displaystyle-\sigma^{2}_{v}f_{0}^{2}(\mu^{2}+2f\mu^{4}+f^{2}\mu^{6})k^{2}P_{L}(k)\in
D(k,\mu),$ (5.13)
$\displaystyle-\sigma^{2}_{\Psi}(1+2f\mu^{2}+f^{2}\mu^{4})k^{2}P_{L}(k)\in
P_{s}^{K,\text{NL}}(k,\mu),$ (5.14)

with $\sigma^{2}_{\Psi}=\int_{0}^{\infty}dp\;P_{L}(p)/(6\pi^{2})$ the
1-dimensional variance of Lagrangian displacements, and
$\sigma^{2}_{v}=\int_{0}^{\infty}dp\;P^{L}_{\theta\theta}(p)/(6\pi^{2})$ the
velocity variance. Notice that these two correlators are identical for the
massless neutrino case, and are simply named $\sigma^{2}_{v}$ in that
instance. Therefore, the EFT contribution to the power spectrum becomes [82,
83]

$P^{\text{EFT}}_{s}(k,\mu)=P^{\text{EPT}}_{s}(k,\mu)+(\alpha_{0}+\alpha_{2}\mu^{2}+\alpha_{4}\mu^{4}+\alpha_{6}\mu^{6})k^{2}P_{L}(k)+P_{\epsilon\epsilon}(k,\mu),$
(5.15)

where we have neglected the growth rate scale-dependence in eqs. (5.13) and
(5.14), since at the scales the counterterms become important, it has reached
its small scale limit $f\approx(1-3f_{\nu}/5)f_{0}$. However, adding the
$f(k)$ factors to the above equation can be done without complications. Below,
we will compute multipoles of the power spectrum up to $\ell=4$, thus the
counterterm $\alpha_{6}$ would be redundant and not considered in this work.
One should notice that in real space, EFT contributions are included to model
the small scale physics out of the reach of PT, and to tame UV-divergencies
(at one-loop present in $P_{13}$). In redshift space the counterterms further
model the non-linear relation between overdensities in real and redshift
spaces. Being this similar to the approach of some earlier works in the
subject, which add to the velocity dispersion variance, $\sigma^{2}_{v}$, a
phenomenological free parameter, $\sigma^{2}_{\text{FoG}}$, to better
modelling the FoG.

With respect to a possible stochastic terms, uncorrelated with long wave-
length fluctuations, we use white noise, plus a tilt proportional to
$(k\mu)^{2}$, resulting in the additional power spectra contribution [82, 83,
84, 85, 86]

$P_{\epsilon}(k,\mu)=P_{shot}\big{[}\alpha^{shot}_{0}+\alpha^{shot}_{2}(k\mu)^{2}\big{]},$
(5.16)

with the standard Poisson process shot-noise as normalization factor
$P_{shot}=1/\bar{n}_{X}$, where $\bar{n}_{X}$ is the mean number density of
tracers. The departure from a white noise arises because stochasticity is not
localized at a single point [30, 87] and because the stochastic nature of
peculiar velocities at small scales [83]. Moreover, stochastic contributions
of the form $k^{2n}$ are also necessary for the renormalization of contact
terms that typically diverge in Fourier space [88, 89, 90].

Finally, and to complete the full bias framework, we need to introduce biasing
in functions $A^{\text{TNS}}$ and $D$. These are obtained by weighting them
with the correct power 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}),$
(5.17) $\displaystyle D(k,\mu;f_{0})$
$\displaystyle\,\longrightarrow\,b_{1}^{4}D(k,\mu;f_{0}/b_{1}).$ (5.18)

It turns out that the biasing for the $D$ function is exact, because $D$ is
constructed out of only linear fields. Instead, the function $A^{\text{TNS}}$
is also biased by the non-linear $\delta^{2}_{cb}$ and $s^{2}$ operators.
However, their contributions are very small so that we neglect them; see also
[91] for using this approximation in the reduction in bias parameters for the
TNS model. As a reference, the complete biased expression of
$A^{\text{TNS}}(k,\mu)$, including $b_{2}$ and $b_{s^{2}}$ parameters, can be
found in appendix A.1 of [25].

###  5.2 IR-resummations

Long wave-length, bulk displacements of matter tend to degrade the BAO in the
power spectrum, because coherent flows stream over a scale settled by (2
times) the 1-dimensional variance of Lagrangian displacements
$\sigma^{2}_{\Psi}$, which is comparable in size to the BAO peak width. As a
result, overdense regions are partially depleted, while underdense regions are
partially populated, broadening the acoustic peak [92, 93, 94]. This effect is
reasonably well described by LPT, even at its first order, the Zeldovich
Approximation (ZA), since bulk flows are captured by the linear displacement
fields [63, 95, 96, 83, 86]. In contrast, in the EPT the convergence is very
slow (see e.g. [97]), thus non-perturbative methods are commonly used to
account for this effect. Here, we employ IR-resummation methods [27] to model
the degradation of the BAO features. We follow the prescription of [28, 29,
98], 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}$;777We perform this decomposition using the fast sine transform recipe
of [99]. 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 one-loop IR-resummed
EFT redshift-space power spectrum becomes [29]

$\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),$
(5.19)

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

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

where $f$ is the scale-dependent growth rate in the presence of massive
neutrinos, and

$\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],$
(5.21) $\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),$
(5.22)

where $\ell_{\text{BAO}}\simeq 105\,h^{-1}\text{Mpc}$ is the BAO peak scale
and $j_{n}$ are the spherical Bessel functions of degree $n$. The scale
$k_{s}$ separates the long and short modes, whose choice is somewhat
arbitrary. However, final results depend weakly on it if $k_{s}\gtrsim
0.05\,h\,\text{Mpc}^{-1}$. Here we use the value
$k_{s}=0.2\,h\,\text{Mpc}^{-1}$, following [98]. We use the same standard IR-
resummed method as in the prescription of [29], with the only difference that
the damping (5.20) is scale-dependent because of the $f(k)$ presence. Notice
that we have no formal derivation of this IR-resummation in the presence of
massive neutrinos, however, we expect to be at the same level of approximation
than in [29] given that the resummed long-wavelengths behave essentially as in
the massless neutrino case.

The IR-resummed EFT power spectrum of eq. (5.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),$
(5.23)

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

##  6 Numerical results for tracers

Figure 4: Comparison of the halo power spectrum multipoles obtained from the
Quijote simulations and the predictions from the perturbation theory: monopole
(black lines), quadrupole (red) and hexadecapole (green). The different panels
show the cases $M_{\nu}=0$ and $0.4$ eV at redshifts $z=0.5$ and $z=1$ (the
cases $M_{\nu}=0.1$ and $0.2$ eV give similar results). The fittings to data
are performed up to $k_{\text{max}}=0.25\,h\,\text{Mpc}^{-1}$ (for comparison
of different choices of $k_{\text{max}}$ see fig. 7). The analytical monopole
results lie within the 1% from the data points and the quadrupole within the
$3\%$. The Poisson noise $P_{shot}$ has been subtracted from the monopole.

Figure 5: Posterior contours for combinations of bias parameters in the
massive neutrino cosmology with $M_{\nu}=0.1$ eV at redshift $z=0.5$. In the
left panel we show the contour plot at 0.68 and 0.95 C.L. for the subspace
$b_{s^{2}}$-$b_{3nl}$. The gray and vertical bands show the expected values
from co-evolution in EdS. In the right panel we show the contour plot for the
subspace $b_{1}$-$b_{2}$.

Figure 6: Power spectrum for tracers: monopole (top panels) and quadrupole
(bottom panels) at $z=0.5$ for $M_{\nu}=0.0\,\text{eV}$ (solid purple line),
$M_{\nu}=0.1\,\text{eV}$ (dotted blue), $M_{\nu}=0.2\,\text{eV}$ (dot-dashed
brown) and $M_{\nu}=0.4\,\text{eV}$ (dashed red) compared to $N$-body
simulations data. For visualization purposes, we have multiplied the results
by constant factors
$A_{s}^{M_{\nu}=0}/A_{s}^{M_{\nu}=0,0.1,0.2,0.4\,\text{eV}}$, such that the
corresponding unbiased multipoles in all models and their linear real space
power spectra have the same primordial amplitude, instead of the same
$\sigma_{8}$ normalization. We do this because for fixed $\sigma_{8}$, the
differences among the different massive neutrino cases are located at large
scales where the sample variance is large. The dotted red lines show the
linear Kaiser result for the $M_{\nu}=0.4\,\text{eV}$ case. The right panels
show the ratios over the best fit massless neutrino case, together with a
dispersion around the central values by sampling 500 times the parameter
posterior distributions.  Figure 7: Comparison when the fitting to data of the
halo power spectra is performed up to $k_{\text{max}}=0.25$, $0.35$ and
$0.45\,h\,\text{Mpc}^{-1}$, for sum of neutrino masses
$M_{\nu}=0.1\,\text{eV}$ at redshift $z=0.5$. In the left panels, we show the
relative differences of the monopole and quadrupole of the theory and the data
from simulations. We observe that the fitting is similar in all cases. In the
right panel, we show the Zeldovich Approximation correlation function using
the linear biases obtained from the different power spectrum fittings. We
observe that as larger is $k_{\rm max}$, the poorer is the correlation
function fitting to the data.  Figure 8: Posterior contours and 1-dimensional
marginalized posterior densities for the linear bias parameter $b_{1}$ and the
shot noise parameter. This plot is for the massive neutrino cosmology with
$M_{\nu}=0.1$ eV at $z=0.5$. The units of $\alpha^{shot}_{2}$ are
$h^{-2}\,\text{Mpc}^{2}$.

We are now in position to compare our analytical predictions to the halo
redshift space power spectrum multipoles obtained from $N$-body simulations.
To this end we use the Quijote suite of 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. We
use $N_{r}=100$ realizations for each cosmology, each one with a volume
$(1\,\text{Gpc}\,h^{-1})^{3}$, containing $512^{3}$ CDM particles and
$512^{3}$ neutrino particles. All simulations have the same $\sigma_{8}$
value, corresponding to scalar primordial amplitudes
$A_{s}^{M_{\nu}=0}=2.13\times 10^{-9}$,
$A_{s}^{M_{\nu}=0.1\,\text{eV}}=2.25\times 10^{-9}$,
$A_{s}^{M_{\nu}=0.2\,\text{eV}}=2.40\times 10^{-9}$ and
$A_{s}^{M_{\nu}=0.4\,\text{eV}}=2.74\times 10^{-9}$. 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 [100] over the CDM particles only, with a linking
length parameter $b=0.2$. We choose to show results for the redshifts $z=0.5$
and $z=1$, though the case $z=0$ has also been checked, finding similar
behaviours.

To compare our analytical model of eq. (5.23) against the simulated data we
let free nine parameters: four bias parameters
$\\{b_{1},b_{2},b_{s^{2}},b_{3nl}\\}$, three counterterms
$\\{\alpha_{0},\alpha_{2},\alpha_{4}\\}$, and two stochastic parameters
$\\{\alpha^{shot}_{0}-1,\alpha^{shot}_{2}\\}$. Given that we subtracted the
Poissonian shot noise $P_{shot}$ from the monopole data, hence the parameter
to fit is $\alpha^{shot}_{0}-1$ instead of $\alpha^{shot}_{0}$. Moreover, for
our fittings we only use the diagonal elements of the covariance matrix,
because the non-diagonal ones are dominated by statistical noise, even for
different multipoles at the same $k$-bin. In order to sample the parameter
space, we run Monte Carlo Markov Chains (MCMC) using the emcee sampler [101]
and the GetDist Python package to present the posteriors [102]. We choose
uniform priors over sufficiently wide intervals because apriori all sensible
values are equally likely. (Notice that some recent work, which also estimate
cosmological parameters, choose Gaussian priors around expected values
obtained through $N$-body simulations and co-evolution [22, 84, 98].)

The first three non-vanishing multipoles ($\ell=0,2,4$) with the best fit
parameters for redshifts $z=0.5$ and $z=1$, together with those obtained from
the simulations, are depicted in figure 4. The error bars denote the RMS
errors for the $N_{r}$ realizations on each model. We choose to show the
fiducial cosmology $M_{\nu}=0$ and the massive neutrino
$M_{\nu}=0.4\,\text{eV}$, the latter presenting the largest deviation from the
massless neutrino cosmology. With the exception of the hexadecapole, which is
very noisy, the analytical results lie within the RMS error bars at all wave-
number bins. We add that we have found very similar agreement for the rest of
massive neutrino cosmology cases: $M_{\nu}=0.1$ and $0.2$ eV. In the left
panel of figure 5 we plot the posterior distribution contour for the subspace
$b_{s^{2}}$-$b_{3nl}$ at $1\sigma$ and $2\sigma$, together with their values
given by co-evolution in EdS, given by eqs. (4.23) (vertical gray bands, with
width given by the dispersion of $b_{1}$). Although the dispersion in
$b_{s^{2}}$-$b_{3nl}$ is large, the final results are not affected by this:
indeed, we also fixed the bias parameters $b_{s^{2}}$ and $b_{3nl}$ to their
co-evolution values before performing the MCMC, and find no differences in the
precision of the final result; and importantly, the linear bias estimations
remains almost equal. In the right panel of figure 5 we show the contour plots
for the subspace $b_{1}$-$b_{2}$. The observed positive correlation is
expected from Peak-Background-Split with a Sheth-Tormen (or Press-Schechter)
mass function for halos in this mass range [103, 104, 87]. Also, the
differences in the large-scale bias is small, consistent with the almost
universal halo mass function for cosmologies that have the same amplitude
$\sigma_{8}$ [68]. We find the best fit values for the linear local bias
$b_{1}=1.72,1.71,1.69,1.67$, for $M_{\nu}=0,0.1,0.2,0.4\,\text{eV}$
respectively, when fitting up to $k_{max}=0.25\,h\,\text{Mpc}^{-1}$ for
redshift $z=0.5$. Similar values are found for the same halo catalogue in [24]
using the 2-point correlation function; see Table 1 of that paper. Notice that
in that work Lagrangian linear biases are reported, so they are related to the
Eulerian linear bias used here as $b_{1}=1+b_{1}^{L}$.

In figure 6 we show the suppression due to the free-streaming for the
different neutrino masses and for multipoles $\ell=0$ and $2$. In the left
panels we show the power spectrum multipoles multiplied by $k$, and we further
multiply by factors $A^{M_{\nu}=0}/A^{M_{\nu}}$, such that at large scales the
different massive cases differ only because of the different linear biases
(which nonetheless are very similar for all the models). Remember that the
simulations are normalized to have the same $\sigma_{8}$, hence by multiplying
by the primordial amplitudes we can observe the main differences at scales
around and above the free-streaming, which is better to visualize the power
spectrum suppression. In the right panels we display the dispersion in the
model by sampling the free parameters posterior distribution. We also show the
monopole and quadrupole in the Kaiser power spectrum plus the Poissonian shot
noise, $P_{K}=(b_{1}+\mu^{2}f)^{2}P_{L}+1/\bar{n}_{X}$, in order to observe
the scales at which the linear theory departs from the data. Notice here that
we are not using any damping factor, so the differences between Kaiser and EFT
theories start at quite small $k$.

We want to discern up to what extent the theory is reliable. To this end, we
perform the parameters sampling up to different maximum wave-numbers $k_{max}$
($=0.25$, $0.35$ and $0.45\,h\,\text{Mpc}^{-1}$). In the left panels of figure
7, we plot the ratios of the analytical predictions to the simulated data for
multipoles $\ell=0$ and $2$, where we use the massive neutrino cosmology with
$M_{\nu}=0.1\,\text{eV}$ at redsfhift $z=0.5$. We notice that the accuracy of
the results are similar for the three cases, all presenting a reasonable
fitting to the data up to the different $k_{max}$. However, the estimated
values of the free parameters change, which can be very harmful when one
estimates the cosmological parameters. Since we are only sampling the
nuisances, we put attention to the linear local bias $b_{1}$, for which we
obtain the best fits $b_{1}=1.71$, $1.74$ and $1.78$ for $k_{max}=0.25$,
$0.35$ and $0.45$, respectively. These are considerable differences in the
estimation of a physical parameter that is intrinsic to the halo catalogue;
indeed, we can use different statistics to measure it. In particular, in [24]
using the correlation function it is found the Lagrangian bias value
$b_{1}^{L}=0.70$ (corresponding to $b_{1}=1.70$) for the same catalogue. In
the right panel of figure 7 we show the halo ZA real space correlation
function for the three estimations of $b_{1}$, together with the data from the
simulations; we omit to show the error bars, but they are sufficiently large
to encompass the three analytical curves; see figure 7 of [24]. We observe
that as larger is $k_{max}$, the fitting to the correlation function becomes
poorer. We interpret it as a biased estimation of the linear bias; which is
particularly important in this context since $b_{1}$ is very degenerate with
the total mass of the neutrinos. Two considerations are worthy to mention:
First, the correlation function is insensitive to the shot noise because
stochastic parameters are uncorrelated to density fields at large scales in
configuration space, hence measuring $b_{1}$ using the correlation function is
more direct than with the power spectrum. Second, at large scales the ZA
describes correctly the correlation function, including the BAO damping due to
the bulk flows of CDM particles. Hence we use linear order in LPT to show the
right-panel on figure 7; moreover, moving to higher orders in PT, beyond the
ZA, introduces additional Lagrangian bias parameters that are non-locally
related to the Eulerian biases, with no general prescription to compare
between them. We believe the main reason for this biased estimation of $b_{1}$
is the shot noise influence: at moderate non-linear scales the halo power
spectrum follows essentially a power law, and it is largely dominated by the
$\alpha^{shot}_{2}\mu^{2}k^{2}$ term, which becomes tightly constrained by the
data for large $k_{max}$. In turn, this constrains even more the constant shot
noise, that is degenerate with $b_{1}$ at large scales. Notice also that at
high-$k$ the EFT parameters are quite degenerate with $\alpha^{shot}_{0}$,
because the former scales as $k^{2}P_{L}(k)$ which is close to a constant at
small scales; however this degeneracy becomes less important at large scales
where the degeneracy between $b_{1}$ and the constant piece of the shot noise
is more relevant. This analysis suggests that the one-loop EFT is not valid
for wave-numbers $\gtrsim 0.25\,h\,\text{Mpc}^{-1}$, as expected. Moreover, at
the scales where the $(\mu k)^{2}$ piece of the shot noise dominates, the RSD
are driven mainly by stochastic velocities and probably no cosmological
information can be extracted [83]. In figure 8 we show triangular posterior
contour plots for the shot parameters and $b_{1}$, for the three analyzed
$k_{max}$ showing the degeneracies and tendencies above mentioned.

Finally, we notice that we may add the next-to-leading order EFT correction
introduced in [11], with counterterm $\tilde{c}$, that serves for a better
modeling of the power spectrum damping along the line-of-sight direction due
to the FoG. However, this parameter does not help for a better estimation of
$b_{1}$, while at the same time reaching good accuracy at smaller scales. This
is because 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}$.

In the end, our main interest is to estimate the neutrino mass scale from
data. However, as it stands, our theory is quite slow and not suitable for an
MCMC cosmological parameter exploration. We have devoted this section to only
fit the nuisance parameters, which implies that one would only need to compute
the loop integrals once for each model. Indeed, in our preliminary numerical
computations, we have also used the EdS kernels to fit the data and find
results as good as with our full kernels. However, we have seen in figure 7
that a good fit to the data does not imply a correct estimation of parameters.
We showed this for the linear bias, but it would not be a complete surprise if
a similar situation occurs with the cosmological parameters when using a
perturbative theory that is not complete, especially when the data standard
deviations are large, and one had to explore, e.g., large neutrino masses,
where the EdS kernels and full kernels differ even more. However, to perform a
complete MCMC analysis, we need a method that accelerates the computation of
loop integrals; such endeavour is the subject of the following section.

##  7 FFTLog method

Obtaining loop corrections with the full massive neutrino kernels requires the
computation of several integrals of the form $I(k)=\int_{\bm{p}}K_{I}(k,p,x)$.
Performing these loop calculations is computationally slow, since at each
volume element of the integration we need to solve a system of differential
equations for second and third order growth functions to find the values of
$\mathcal{A},\mathcal{B}(k,p,x)$ and their third order counterparts for
constructing $K_{I}(k,p,x)$. The slowness of this procedure precludes the use
of efficient parameter sampling algorithms for estimation of cosmological
parameters.

However, the dominant contributions to the loop corrections come from the
growth rates instead of these computationally costly functions. The latter
calculations mainly play a normalization role, and their effect is similar to
how they affect the $\Lambda$CDM kernels, as oppose to the use of EdS kernels.
Therefore, one may use this hierarchy of calculations in order to speed up the
one-loop computations. For this purpose, we first define the fk kernels as

$\displaystyle F_{2}^{\texttt{fk}}(\bm{k}_{1},\bm{k}_{2})$
$\displaystyle=F_{2}(\bm{k}_{1},\bm{k}_{2})\Big{|}_{\mathcal{A}=\mathcal{B}=\mathcal{A}^{f_{\nu}=0}},$
(7.1) $\displaystyle G_{2}^{\texttt{fk}}(\bm{k}_{1},\bm{k}_{2})$
$\displaystyle=G_{2}(\bm{k}_{1},\bm{k}_{2})\Big{|}_{\mathcal{A}=\mathcal{B}=\mathcal{A}^{f_{\nu}=0}},$
(7.2)

and similar for the third order kernels. In other words, we fix the functions
$\mathcal{A}$ and $\mathcal{B}$ to have their large scale ($f_{\nu}=0$) value
$\mathcal{A}^{f_{\nu}=0}$, which is the same for both functions and coincides
with the $\Lambda$CDM value, namely

$\mathcal{A}^{\text{$\Lambda$CDM}}(t)=\mathcal{A}^{f_{\nu}=0}(t)=\mathcal{B}^{f_{\nu}=0}(t)=\frac{D_{\mathcal{A}}(t)}{\frac{3}{7}D^{2}_{+}(t)},$
(7.3)

with

$D_{\mathcal{A}}^{f_{\nu}=0}(t)=\frac{3}{7}D_{+}^{2}(t)+\frac{4}{7}\left(\frac{d^{2}\,}{dt^{2}}+2H\frac{d\,}{dt}-\frac{3}{2}\Omega_{m}H^{2}\right)^{-1}\left[\frac{3}{2}\Omega_{m}H^{2}\left(1-\frac{f^{2}}{\Omega_{m}}\right)\right].$
(7.4)

In the simpler case of EdS, $\Omega_{m}=1=f$, and the second term in the rhs
of the eq. (7.4) vanishes, reducing the second order kernels to the well-known
EdS ones with $\mathcal{A}=\mathcal{B}=1$. For $\Lambda$CDM, or equivalently
$f_{\nu}=0$, and in the case of observationally allowed cosmologies, one
obtains numerically $\mathcal{A}^{f_{\nu}=0}(z=0)\approx 1.01$.

In order to asses the accuracy of simplified kernels, we measure the relative
difference of the matter power spectra multipoles for $\Lambda$CDM or fk
kernels with respect to the computation using the full kernels. That is, in
figure 9 we plot the variable

$\frac{\Delta P_{\ell}}{P_{\ell}}=\frac{P_{\ell}^{\text{$\Lambda$CDM kernels,
fk kernels}}}{P_{\ell}^{\text{Full kernels}}}-1,$ (7.5)

for multipoles $\ell=0,2,4$. In the previous expression the multipoles
$P_{\ell}^{\text{fk kernels}}$ use the fk kernels, which we remind the reader
are obtained by considering the pieces of the full kernels with the scale-
dependent growth rates only, while for the pieces that are solution to second
order differential equations we assign their massless neutrino values. We find
that the use of the fk kernels is always very close to the full kernels case
(within $0.3\%$) for $k<0.5\,h\,$Mpc-1. In summary, this accuracy of the fk
kernels allow us to use the FFTLog methods [31, 105, 32, 106, 33, 98] to
further accelerate the computation of loop integrals.

Figure 9: Relative differences of the multipoles of the redshift-space matter
power spectrum when using $\Lambda$CDM kernels (solid lines) and fk kernels
(dot-dashed lines) instead of the full kernels, as given by eq. (7.5).

It is not our intention to work out the FFTLog procedure in detail for all the
expressions in this work (which will be presented elsewhere), but to show
explicitly how this can be done for a particular example. We choose the case
of $P^{22}_{\delta\theta}$, given by the second term in eq. (4.2). In fact, we
also show some numerical results for the one-dimensional $P_{\delta\theta}$
and $P_{\theta\theta}$ spectra. For simplicity we assume
$\mathcal{A}^{f_{\nu}=0}=\mathcal{A}^{\text{EdS}}=1$; though the computation
without this restriction is straightforward and it is outlined at the end of
this section.

First, we expand the linear power spectrum as a sum of scale invariant spectra
with complex powers as [107, 31]

$\bar{P}_{L}(k)=\sum_{m=-N/2}^{N/2}c_{m}k^{\nu+i\eta_{m}},\qquad\eta_{m}=\frac{N-1}{N}\frac{2\pi
m}{\ln(k_{\text{max}}/k_{\text{min}})},$ (7.6)

where we have split an interval $[k_{\text{min}},k_{\text{max}}]$ in $N$
logarithmic spaced wave-numbers. The coefficients $c_{m}$ comes from the
discrete log-Fourier transform

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

with the weights $W_{m}=1$, except for the end points, for which
$W_{-N/2}=W_{N/2}=1/2$. The so-called bias $\nu$ is in principle any real
number, but its value is chosen to have a better convergence for loop
integrals [33].

We need a similar expansion to eq. (7.6) for the linear density-velocity
cross-spectra $P^{L}_{\delta\theta}=(f(k)/f_{0})P_{L}(k)$, namely

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

where the coefficients $c_{m}^{f}$ are computed with eq. (7.7) by substituting
$P_{L}(k)$ by $(f(k)/f_{0})P_{L}(k)$.

In the standard method, the kernels are taken to be EdS and these are written
in powers of $k$, $p$ and $|\bm{k}-\bm{p}|$. The difficulty in our case is the
presence of growth rate functions in $G_{2}$. Hence, we split the integrand of
the loop piece of $P^{22}_{\delta\theta}$ into a contribution with a linear
growth rate evaluated at the wavenumber $p$ and a further piece evaluated at
$|\bm{k}-\bm{p}|$. The splitting results in

$\displaystyle P^{22}_{\delta\theta}$
$\displaystyle=2\int_{\bm{p}}F_{2}^{\text{fk}}(\bm{p},\bm{k}-\bm{p})G_{2}^{\text{fk}}(\bm{p},\bm{k}-\bm{p})P_{L}(p)P_{L}(|\bm{k}-\bm{p}|)$
$\displaystyle=2\int_{\bm{p}}K^{f_{\bm{p}}}_{\delta\theta}(\bm{p},\bm{k}-\bm{p})\frac{f(p)}{f_{0}}P_{L}(p)P_{L}(|\bm{k}-\bm{p}|)$
$\displaystyle\quad+2\int_{\bm{p}}K^{f_{\bm{k}-\bm{p}}}_{\delta\theta}(\bm{p},\bm{k}-\bm{p})P_{L}(p)\frac{f(|\bm{k}-\bm{p}|)}{f_{0}}P_{L}(|\bm{k}-\bm{p}|),$
(7.9)

such that the two functions

$\displaystyle K^{f_{\bm{k}-\bm{p}}}_{\delta\theta}(\bm{p},\bm{k}-\bm{p})$
$\displaystyle=\frac{k^{8}}{196p^{4}|\bm{k}-\bm{p}|^{4}}-\frac{k^{6}}{392p^{2}|\bm{k}-\bm{p}|^{4}}+\frac{3k^{6}}{196p^{4}|\bm{k}-\bm{p}|^{2}}-\frac{9k^{4}}{392|\bm{k}-\bm{p}|^{4}}$
$\displaystyle\quad+\frac{23k^{4}}{784p^{2}|\bm{k}-\bm{p}|^{2}}-\frac{11k^{4}}{784p^{4}}+\frac{13k^{2}}{392p^{-2}|\bm{k}-\bm{p}|^{4}}-\frac{15k^{2}}{392p^{4}|\bm{k}-\bm{p}|^{-2}}$
$\displaystyle\quad-\frac{5k^{2}}{98|\bm{k}-\bm{p}|^{2}}+\frac{11k^{2}}{196p^{2}}-\frac{5}{392p^{-4}|\bm{k}-\bm{p}|^{4}}+\frac{25}{784p^{4}|\bm{k}-\bm{p}|^{-4}}$
$\displaystyle\quad+\frac{5}{784p^{-2}|\bm{k}-\bm{p}|^{2}}-\frac{65}{784p^{2}|\bm{k}-\bm{p}|^{-2}}+\frac{45}{784},$
(7.10)

and

$\displaystyle K^{f_{\bm{p}}}_{\delta\theta}(\bm{p},\bm{k}-\bm{p})$
$\displaystyle=\frac{k^{8}}{196p^{4}|\bm{k}-\bm{p}|^{4}}+\frac{3k^{6}}{196p^{2}|\bm{k}-\bm{p}|^{4}}-\frac{k^{6}}{392p^{4}|\bm{k}-\bm{p}|^{2}}-\frac{11k^{4}}{784|\bm{k}-\bm{p}|^{4}}$
$\displaystyle\quad+\frac{23k^{4}}{784p^{2}|\bm{k}-\bm{p}|^{2}}-\frac{9k^{4}}{392p^{4}}-\frac{15k^{2}}{392p^{-2}|\bm{k}-\bm{p}|^{4}}+\frac{13k^{2}}{392p^{4}|\bm{k}-\bm{p}|^{-2}}$
$\displaystyle\quad+\frac{11k^{2}}{196|\bm{k}-\bm{p}|^{2}}-\frac{5k^{2}}{98p^{2}}+\frac{25}{784p^{-4}|\bm{k}-\bm{p}|^{4}}-\frac{5}{392p^{4}|\bm{k}-\bm{p}|^{-4}}$
$\displaystyle\quad-\frac{65}{784p^{-2}|\bm{k}-\bm{p}|^{2}}+\frac{5}{784p^{2}|\bm{k}-\bm{p}|^{-2}}+\frac{45}{784},$
(7.11)

are independent of the input cosmology. As a result of the previous
expressions, one gets that

$K^{f_{\bm{p}}}_{\delta\theta}(\bm{p},\bm{k}-\bm{p})+K^{f_{\bm{k}-\bm{p}}}_{\delta\theta}(\bm{p},\bm{k}-\bm{p})=F_{2}^{\text{EdS}}(\bm{p},\bm{k}-\bm{p})G_{2}^{\text{EdS}}(\bm{p},\bm{k}-\bm{p}).$
(7.12)

Since all the summands in the kernels above are of 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\\}$, it is convenient to define the matrix
$f^{f_{\bm{p}}}_{22,\delta\theta}(n_{1},n_{2})$, with components given by the
coefficient of each term in eq. (7). Similarly, one constructs the matrix
$f^{f_{\bm{k}-\bm{p}}}_{22,\delta\theta}(n_{1},n_{2})$ using eq. (7).
Therefore, using these matrices and with the aid of the discrete Fourier
transforms [eqs. (7.6) and (7.8)], the $P_{22}$ of eq. (7) can be approximated
by

$\displaystyle P^{22}_{\delta\theta}(k)$
$\displaystyle=2\sum_{m_{1},m_{2}}c^{f}_{m_{1}}c_{m_{2}}\sum_{n_{1},n_{2}=-2}^{2}f^{f_{\bm{p}}}_{22,\delta\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}}c^{f}_{m_{2}}\sum_{n_{1},n_{2}=-2}^{2}f^{f_{\bm{k}-\bm{p}}}_{22,\delta\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}}},$
(7.13)

with

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

To further simplify the previous expression, we define

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

and [108, 109]

$\displaystyle I(a,b)\equiv
k^{-3+2(a+b)}\int_{\bm{p}}\frac{1}{p^{2a}|\bm{k}-\bm{p}|^{2b}}=\frac{1}{8\pi^{3/2}}\frac{\Gamma(\frac{3}{2}-a)\Gamma(\frac{3}{2}-b)\Gamma(a+b-\frac{3}{2})}{\Gamma(a)\Gamma(b)\Gamma(3-a-b)}.$
(7.17)

Using the above equations, one can write the approximation for
$P^{22}_{\delta\theta}$ as

$\displaystyle\bar{P}^{22}_{\delta\theta}(k)$
$\displaystyle=k^{3}\sum_{m_{1},m_{2}}c_{m_{1}}^{f}k^{-2\nu_{1}}\,M_{22,\delta\theta}^{f_{\bm{p}}}(\nu_{1},\nu_{2})\,c_{m_{2}}k^{-2\nu_{2}}$
$\displaystyle\quad+k^{3}\sum_{m_{1},m_{2}}c_{m_{1}}k^{-2\nu_{1}}\,M_{22,\delta\theta}^{f_{\bm{k}-\bm{p}}}(\nu_{1},\nu_{2})\,c^{f}_{m_{2}}k^{-2\nu_{2}}.$
(7.18)

Figure 10: Comparison of the FFTLog and Direct computations for the one-loop
contributions to $P_{\delta\delta}$, $P_{\delta\theta}$ and
$P_{\theta\theta}$. The parameters are
$k_{\text{min}}=10^{-6}\,\,h\,\text{Mpc}^{-1}$,
$k_{\text{max}}=10\,\,h\,\text{Mpc}^{-1}$, $N=256$ and bias $\nu=-0.1$.

We notice that the $M_{22}$ matrices defined in eqs. (7.15) and (7.16) are not
symmetric. However, the interchange $\bm{p}\rightarrow\bm{k}-\bm{p}$ symmetry,
which explicitly holds in the first line of eq. (7), implies
$M_{22,\delta\theta}^{f_{\bm{p}}}(\nu_{1},\nu_{2})=M_{22,\delta\theta}^{f_{\bm{k}-\bm{p}}}(\nu_{2},\nu_{1})$.
This equality permits us to reduce the previous expression to

$\displaystyle\bar{P}^{22}_{\delta\theta}(k)$
$\displaystyle=2k^{3}\sum_{m_{1},m_{2}}c_{m_{1}}^{f}k^{-2\nu_{1}}\,M_{22,\delta\theta}^{f_{\bm{p}}}(\nu_{1},\nu_{2})\,c_{m_{2}}k^{-2\nu_{2}}.$
(7.19)

Note that the EdS matrix is recovered as
$M_{22,\delta\theta}^{\text{EdS}}=M_{22,\delta\theta}^{f_{\bm{p}}}+M_{22,\delta\theta}^{f_{\bm{k}-\bm{p}}}$.

The computation of the rest of loop contributions to construct the one-
dimensional spectra of eqs. (4.1), (4.2) and (4.3) is similar, but lengthy so
we do not show them here. In figure 10 we show the results for
$P_{\delta\delta}$, $P_{\delta\theta}$ ($\ni P^{22}_{\delta\theta}$) and
$P_{\theta\theta}$ computed using FFTLog and the Direct integration methods.
To enhance the effects of massive neutrinos we work the case of
$M_{\nu}=0.4\,\text{eV}$, for which the introduction of the scale dependent
growth rate is more important. In conclusion, we show that the FFTLog method,
using a bias of $\nu=-0.1$,888After a few tests we found this number yield
better convergence. For bias $\nu>-1$ one has to add the UV limits of the
$P_{13}$ functions. works well to a very high accuracy, smaller than the
$0.5\%$ over the interval $\,0.001<k<1\,h\,\text{Mpc}^{-1}$.

Finally, let us discuss how to introduce the correct time dependence of the
$\mathcal{A}$ function that we have so far neglected. Notice that the
coefficients in the kernels of eqs. (7) and (7) would become multiplied by
factors $\epsilon^{m}d^{n}$ (with $m+n=0,1,2$), where

$\epsilon=\mathcal{A}^{f_{\nu}=0}(t)-1\lesssim 0.01,\qquad\text{and}\qquad
d=\frac{\dot{\mathcal{A}}^{f_{\nu}=0}(t)}{f_{0}H}\lesssim 0.01.$ (7.20)

As a result, the $M_{22}$ matrices can be split as

$M_{22}=M_{22}\big{|}_{\epsilon=d=0}+\epsilon\Big{[}M_{22}-M_{22}\big{|}_{\epsilon=d=0}\Big{]}_{\epsilon=1,d=0}+d\Big{[}M_{22}-M_{22}\big{|}_{\epsilon=d=0}\Big{]}_{\epsilon=0,d=1},$
(7.21)

where we neglect terms of order $\epsilon^{m}d^{n}$ with $m+n=2$.

This method can also be applied to find the $\Lambda$CDM model with the exact
time dependence (or general models with non-clustering dark energy), at the
cost of introducing two additional matrices. In such a case, one would have
$M_{22}^{\text{$\Lambda$CDM}}\big{|}_{\epsilon=d=0}=M^{\text{EdS}}_{22}$.

##  8 Summary and Conclusions

Despite neutrinos being non-relativistic particles at late times and behaving
analogous to other matter components, rather than radiation, they are not cold
components due to their large velocity dispersions. Depending on the scales of
interest, neutrinos present two distinct behaviors. In regions below the free-
streaming scale, a Jeans-like mechanism with a pressure supported by the
neutrino thermal velocity, prevents their confinement. In contrast, for scales
well above this scale, they cluster in the same way as the CDM component. As a
consequence, matter statistics become suppressed below the free-streaming when
compared to the equivalent massless neutrino scenario.

In the present work, we built upon the work of refs. [24] and [25] to
construct an Eulerian Perturbation Theory (EPT) for the Large-Scale clustering
in massive neutrino cosmologies. The most important assumption we made is
given by eq. (2.3), where the non-linear neutrino density is approximated by
the non-linear $cb$ density fluctuation times the ratio of the neutrinos to
$cb$ transfer functions. As a result, one obtains a theory that exhibits a
scale-dependent gravitational strength, which has been partially motivated by
the recent developments in PT for modified gravity models (e.g. [110, 60]). Of
particular relevance for this work are the findings of [25], where a
generalization of the kernels beyond the EdS approximation were developed for
the density-weighted velocity moments (see [26]).

We do not try to obtain the EPT kernels directly from the continuity, Euler
and Poisson equations, but instead, we map the Lagrangian Perturbation Theory
(LPT) kernels of [24] into the Eulerian frame. This approach is considerably
more straightforward, and for theories with additional scales, the computation
of LPT kernels is generally simpler than for the EPT counterparts. However,
and for the sake of completeness, we show how to derive the second order
kernels directly in appendix A, matching the other approach. Moreover, our EPT
kernels are well behaved in the UV and IR limits, and reduce to the
$\Lambda$CDM expressions when using the large scale limit in the function
$A(k,t)$, introduced in eq. (2.2).

In order to asses the accuracy of our predictions we contrast our modelling
with synthetic data. In particular, we compare the first three non-vanishing
multipoles of the power spectrum to a halo catalogue from the Quijote
simulations. We find good agreement up to the commonly expected validity scale
in the EFT ($k\sim 0.25\,h\,\text{Mpc}^{-1}$). We can extend to higher wave-
numbers with increasing precision, but by doing so we are incorrectly
estimating the linear bias. One need to remember that this bias is physical,
intrinsic to the halo catalogue and degenerate with the neutrino masses, hence
this is potentially harmful for the sum of neutrino mass estimations.

Given the small discrepancies of the power spectrum predictions between EdS
and full kernels for realistic neutrino masses, the functional liberty allowed
by counterterms, biasing and stochasticity, as well as the expected precision
of current and upcoming galaxy surveys, one may be skeptical about the
necessity of using a more complicated theory, which is also much more
demanding computationally. However, recent analyses of large scale structure
BOSS data [11, 12] have shown that the constraining power on the neutrino
masses of surveys without adding further data-sets is limited, with the need
of sampling parameters up to large masses (of a few eV), where the use of EdS
kernels is questionable. Another reason to use neutrino kernels is that higher
order statistics may break degeneracies between the galaxy bias and the
neutrino masses. In particular, recent results for the redshift space
bispectrum show that the use of a proper description of $G_{2}$, away from the
EdS approximation, may help breaking degeneracies between $b_{1}$ and
$f_{\nu}$ [18]. In spite of these motivations, the computation of loop
corrections with full kernels is quite expensive, which makes the use of this
theory unviable for routinely cosmological parameter estimation. For this
reason, we reduce the information content of the kernels, but keep its main
feature, the inherited —from linear theory— effect of having a scale-dependent
growth rate. The advantage of these reduced kernels is twofold: first, there
is no need of solving differential equations at each step of the loop
correction quadratures; and second, these kernels are suitable for an FFTLog
method. We make some developments in the latter and find very good agreement
with the direct integration method for the three first density-weighted
velocity moments.

On the other hand, we have discussed at the end of section 6 how the PT–EFT
theory stops being that accurate at around $k=0.25\,h\,\text{Mpc}^{-1}$.
Actually, one is able to fit the data with high precision at higher
wavelengths, but it results in an inaccurate estimation of the linear bias.
This is better illustrated in figure 7 where both, the correlation function
and the power spectrum are fit to the data. It would not be a surprise if this
same biased estimation when using large wavelengths happens with the
cosmological parameters. In our early numerical results we found that by
fitting the nuisance parameters with the use of EdS kernels, instead of the
exact ones, one obtains good results, with slightly different best-fitting
nuisance parameters. But in the same line of thought as above, we can
speculate further and conceive that the use of EdS kernels in fitting the
cosmological parameters may lead to the wrong estimation of parameters,
particularly when the standard deviations of measurements are large and one
has to explore large neutrino masses. Therefore, exploring the whole space of
parameters with the use of MCMC algorithms may turn out to be very useful. In
this context, the full FFTLog analysis is a natural avenue to carry on this
exploration and we leave it for future work.

## Acknowledgments

We thank Hernán E. Noriega and Mario A. Rodriguez-Meza for useful discussions,
together with the Instituto Avanzado de Cosmología A. C. and the DCI-UG
DataLab for academic and computational resources. AA and GN acknowledge the
support of CONACyT, specially through the projects 283151, 286897 and Ciencia
de Frontera No. 102958. GN also appreciates the grants of UG-DAIP. AB was
supported by the Fermi Research Alliance, LLC under Contract No. DE-
AC02-07CH11359 with the U.S. Department of Energy, and the U.S. Department of
Energy (DOE) Office of Science Distinguished Scientist Fellow Program

##  Appendix A Direct computation of second order Eulerian kernels

In this appendix we find the second order EPT kernels directly from the fluid
equations and show that these coincide with those obtained with the mapping
between Lagrangian and Eulerian frames used in section 3. The continuity and
Euler equations subjected to the scale dependent gravitational strength $A(k)$
are

$\displaystyle\frac{1}{H}\frac{\partial\delta(\bm{k})}{\partial
t}-f_{0}\theta(\bm{k})$
$\displaystyle=f_{0}\underset{\bm{k}_{12}=\bm{k}}{\int}\alpha(\bm{k}_{1},\bm{k}_{2})\theta(\bm{k}_{1})\delta(\bm{k}_{2}),$
(A.1) $\displaystyle\frac{1}{H}\frac{\partial f_{0}\theta(\bm{k})}{\partial
t}+\left(2+\frac{\dot{H}}{H^{2}}\right)f_{0}\theta(\bm{k})$
$\displaystyle-\frac{A(k)}{H^{2}}\delta(\bm{k})=f_{0}^{2}\underset{\bm{k}_{12}=\bm{k}}{\int}\beta(\bm{k}_{1},\bm{k}_{2})\theta(\bm{k}_{1})\theta(\bm{k}_{2}),$
(A.2)

with

$\alpha(\bm{k}_{1},\bm{k}_{2})=1+\frac{\bm{k}_{1}\cdot\bm{k}_{2}}{k_{1}^{2}},\qquad\beta(\bm{k}_{1},\bm{k}_{2})=\frac{k_{12}^{2}(\bm{k}_{1}\cdot\bm{k}_{2})}{2k_{1}^{2}k_{2}^{2}},$
(A.3)

where we used the Poisson equation $(k/a)^{2}\Phi=-A(k,t)\delta$ and the
definition of the rescaled velocity divergence $\theta=-ik_{i}v^{i}/(aHf_{0})$
in terms of the peculiar velocity $v^{i}$.

We introduce the EPT kernels through

$\displaystyle\delta^{(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_{L}(\bm{k}_{1},t)\cdots\delta_{L}(\bm{k}_{n},t),$
$\displaystyle\theta^{(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_{L}(\bm{k}_{1},t)\cdots\delta_{L}(\bm{k}_{n},t).$
(A.4)

Hence, at first order,

$\displaystyle F_{1}(\bm{k})$ $\displaystyle=1,\qquad\text{and}\qquad
G_{1}(\bm{k})=\frac{f(k)}{f_{0}},$ (A.5)

with the growth rate $f(k,t)=d\ln D_{+}/d\ln a$, $f_{0}=f(k\rightarrow 0)$,
and the linear growth function solving the differential equation
$\ddot{D}_{+}+2H\dot{D}_{+}=A(k)D_{+}$, with appropriate initial conditions
that pick out the fastest growing solution.

To second order, the fluid equations are

$\displaystyle H^{-1}\frac{\partial\delta^{(2)}(\bm{k})}{\partial
t}-f_{0}\theta^{(2)}(\bm{k})$
$\displaystyle=f_{0}\underset{\bm{k}_{12}=\bm{k}}{\int}\alpha(\bm{k}_{1},\bm{k}_{2})\theta^{(1)}(\bm{k}_{1})\delta^{(1)}(\bm{k}_{2})$
$\displaystyle=\frac{1}{2}\underset{\bm{k}_{12}=\bm{k}}{\int}\big{[}\alpha(\bm{k}_{1},\bm{k}_{2})f(k_{1})+\alpha(\bm{k}_{2},\bm{k}_{1})f(k_{2})\big{]}\delta^{(1)}(\bm{k}_{1})\delta^{(1)}(\bm{k}_{2}),$
(A.6) $\displaystyle H^{-1}\frac{\partial f_{0}\theta^{(2)}(\bm{k})}{\partial
t}+f_{0}\left(2+\frac{\dot{H}}{H^{2}}\right)\theta^{(2)}(\bm{k})$
$\displaystyle-\frac{A(k)}{H^{2}}\delta^{(2)}(\bm{k})=f_{0}^{2}\underset{\bm{k}_{12}=\bm{k}}{\int}\beta(\bm{k}_{1},\bm{k}_{2})\theta^{(1)}(\bm{k}_{1})\theta^{(1)}(\bm{k}_{2}),$
$\displaystyle=\underset{\bm{k}_{12}=\bm{k}}{\int}\beta(\bm{k}_{1},\bm{k}_{2})f(k_{1})f(k_{2})\delta^{(1)}(\bm{k}_{1})\delta^{(1)}(\bm{k}_{2}),$
(A.7)

where we have used $\theta^{(1)}(\bm{k})=(f(k)/f_{0})\delta^{(1)}(\bm{k})$,
and inside the integral of the rhs of eq. (A.6) we have symmetrized over.

The second order density fluctuation and velocity fields are

$\displaystyle\delta^{(2)}(\bm{k})$
$\displaystyle=\underset{\bm{k}_{12}=\bm{k}}{\int}F_{2}(\bm{k}_{1},\bm{k}_{2})D_{+}(\bm{k}_{1},t)D_{+}(\bm{k}_{2},t)\delta_{0}(\bm{k}_{1})\delta_{0}(\bm{k}_{2}),$
$\displaystyle\theta^{(2)}(\bm{k})$
$\displaystyle=\underset{\bm{k}_{12}=\bm{k}}{\int}G_{2}(\bm{k}_{1},\bm{k}_{2})D_{+}(\bm{k}_{1},t)D_{+}(\bm{k}_{2},t)\delta_{0}(\bm{k}_{1})\delta_{0}(\bm{k}_{2})$
(A.8)

with $\delta_{0}(\bm{k})=\delta^{(1)}(\bm{k},t_{0})$. Inserting these
expressions into eqs. (A.6) and (A.7),

$\displaystyle\frac{1}{HD_{1}D_{2}}\frac{d\,\,}{dt}(F_{2}D_{1}D_{2})-f_{0}G_{2}=\frac{1}{2}(\alpha_{12}f_{1}+\alpha_{21}f_{2})$
(A.9)
$\displaystyle\frac{1}{HD_{1}D_{2}}\frac{d\,\,}{dt}(f_{0}G_{2}D_{1}D_{2})+\left(2+\frac{\dot{H}}{H^{2}}\right)f_{0}G_{2}-\frac{A(k)}{H^{2}}F_{2}=\beta_{12}f_{1}f_{2},$
(A.10)

with $f_{1,2}=f(k_{1,2})$, $D_{1,2}=D_{+}(\bm{k}_{1,2},t)$,
$\alpha_{12}=\alpha(\bm{k}_{1},\bm{k}_{2})$,
$\alpha_{21}=\alpha(\bm{k}_{2},\bm{k}_{1})$ and
$\beta_{12}=\beta(\bm{k}_{1},\bm{k}_{2})$. We rewrite the above equations as

$\displaystyle\frac{1}{H}\frac{dF_{2}}{dt}+F_{2}(f_{1}+f_{2})-f_{0}G_{2}=\frac{1}{2}(\alpha_{12}f_{1}+\alpha_{21}f_{2}),$
(A.11)
$\displaystyle\frac{1}{H}\frac{df_{0}G_{2}}{dt}+f_{0}G_{2}(f_{1}+f_{2})+\left(2+\frac{\dot{H}}{H^{2}}\right)f_{0}G_{2}-\frac{A(k)}{H^{2}}F_{2}=\beta_{12}f_{1}f_{2}.$
(A.12)

Taking the time derivative of eq. (A.11) and using eq. (A.12) we obtain a
second order equation for $F_{2}$,

$\displaystyle\frac{1}{H^{2}}\ddot{F}_{2}+\frac{2}{H}(1+f_{1}+f_{2})\dot{F}_{2}+\left[\frac{1}{H}(\dot{f}_{1}+\dot{f}_{2})+(f_{1}+f_{2})\left(f_{1}+f_{2}+2+\frac{\dot{H}}{H^{2}}\right)-\frac{A(k)}{H^{2}}\right]F_{2}=$
$\displaystyle\qquad\frac{1}{2H}(\alpha_{12}\dot{f}_{1}+\alpha_{21}\dot{f}_{2})+\frac{1}{2}(\alpha_{12}f_{1}+\alpha_{21}f_{2})\left(f_{1}+f_{2}+2+\frac{\dot{H}}{H^{2}}\right)+\beta_{12}f_{1}f_{2},$
(A.13)

Now, the growth rate $f(k_{1})$ evolves as

$\dot{f}_{1}=\frac{A(k_{1})}{H}-H\left(2+\frac{\dot{H}}{H^{2}}\right)f_{1}-Hf_{1}^{2},$
(A.14)

and an equivalent expression for $f_{2}$. Substituting it in eq. (A.13),

$\displaystyle\frac{1}{H^{2}}\ddot{F}_{2}+\frac{2}{H}(1+f_{1}+f_{2})\dot{F}_{2}+\left[2f_{1}f_{2}+\frac{A(k_{1})+A(k_{2})-A(k)}{H^{2}}\right]F_{2}=$
$\displaystyle\qquad\frac{1}{2}\alpha_{12}\frac{A(k_{1})}{H^{2}}+\frac{1}{2}\alpha_{21}\frac{A(k_{2})}{H^{2}}+\frac{1}{2}f_{1}f_{2}(\alpha_{12}+\alpha_{21})+\beta_{12}f_{1}f_{2}.$
(A.15)

Now, let us define a second order growth function as

$D^{(2)}(\bm{k}_{1},\bm{k}_{2},t)\equiv D_{12}\equiv
2D_{1}D_{2}F_{2}-\chi_{12},$ (A.16)

and

$\displaystyle\chi_{12}$
$\displaystyle\equiv\alpha_{12}+\alpha_{21}-\gamma_{12},\qquad\text{with}\qquad\gamma_{12}\equiv
1-\frac{(\bm{k}_{1}\cdot\bm{k}_{2})^{2}}{k_{1}^{2}k_{2}^{2}}.$ (A.17)

We will now find a differential equation for $D_{12}$. First, the second order
$F_{2}$ kernel is

$F_{2}=\frac{D_{12}}{2D_{1}D_{2}}+\frac{1}{2}\chi_{12}.$ (A.18)

Now, taking time derivatives of the above equation,

$\frac{1}{H}\dot{F}_{2}=\frac{1}{2D_{1}D_{2}}\left(\frac{1}{H}\dot{D}_{12}-D_{12}(f_{1}+f_{2})\right),$
(A.19)

$\displaystyle\frac{1}{H^{2}}\ddot{F}_{2}$
$\displaystyle=\frac{1}{2D_{1}D_{2}}\Bigg{[}\frac{1}{H^{2}}\ddot{D}_{12}-\frac{2}{H}(f_{1}+f_{2})\dot{D}_{12}$
$\displaystyle\quad\qquad+D_{12}\left(2(f_{1}^{2}+f_{2}^{2}+f_{1}f_{2}+f_{1}+f_{2})-\frac{1}{H^{2}}(A(k_{1})+A(k_{2}))\right)\Bigg{]},$
(A.20)

where we used eq. (A.14) and $\ddot{D}+2H\dot{D}=A(k)D$. Substituting the
above equations into eq. (A),

$\displaystyle\ddot{D}_{12}+2H\dot{D}_{12}-A(k)D_{12}$
$\displaystyle=\Bigg{[}A(k)+(A(k)-A(k_{2}))\frac{\bm{k}_{1}\cdot\bm{k}_{2}}{k^{2}_{1}}+(A(k)-A(k_{1}))\frac{\bm{k}_{1}\cdot\bm{k}_{2}}{k^{2}_{2}}$
$\displaystyle\quad\qquad-(A(k_{1})+A(k_{2})-A(k))\frac{(\bm{k}_{1}\cdot\bm{k}_{2})^{2}}{k_{1}^{2}k_{2}^{2}}\Bigg{]}D_{1}D_{2}.$
(A.21)

Hence, using eq. (A.18) we obtain

$\displaystyle F_{2}(\bm{k}_{1},\bm{k}_{2})$
$\displaystyle=\frac{1}{2}+\frac{3}{14}\mathcal{A}+\left(\frac{1}{2}-\frac{3}{14}\mathcal{B}\right)\frac{(\bm{k}_{1}\cdot\bm{k}_{2})^{2}}{k_{1}^{2}k_{2}^{2}}+\frac{\bm{k}_{1}\cdot\bm{k}_{2}}{2k_{1}k_{2}}\left(\frac{k_{2}}{k_{1}}+\frac{k_{1}}{k_{2}}\right),$
(A.22)

and from eq. (A.11),

$\displaystyle G_{2}(\bm{k}_{1},\bm{k}_{2})$
$\displaystyle=\frac{3\mathcal{A}(f_{1}+f_{2})+3\dot{\mathcal{A}}/H}{14f_{0}}+\left(\frac{f_{1}+f_{2}}{2f_{0}}-\frac{3\mathcal{B}(f_{1}+f_{2})+3\dot{\mathcal{B}}/H}{14f_{0}}\right)\frac{(\bm{k}_{1}\cdot\bm{k}_{2})^{2}}{k_{1}^{2}k_{2}^{2}}$
$\displaystyle\quad+\frac{\bm{k}_{1}\cdot\bm{k}_{2}}{2k_{1}k_{2}}\left(\frac{f_{2}}{f_{0}}\frac{k_{2}}{k_{1}}+\frac{f_{1}}{f_{0}}\frac{k_{1}}{k_{2}}\right),$
(A.23)

with $\mathcal{A}$ and $\mathcal{B}$ given by

$\mathcal{A}(\bm{k}_{1},\bm{k}_{2},t)=\frac{7D^{(2)}_{\mathcal{A}}(\bm{k}_{1},\bm{k}_{2},t)}{3D_{+}(k_{1},t)D_{+}(k_{2},t)},\qquad\mathcal{B}(\bm{k}_{1},\bm{k}_{2},t)=\frac{7D^{(2)}_{\mathcal{B}}(\bm{k}_{1},\bm{k}_{2},t)}{3D_{+}(k_{1},t)D_{+}(k_{2},t)},$
(A.24)

with second order growth functions $D^{(2)}_{\mathcal{A},\mathcal{B}}$ the
solutions to the Green problem:

$\displaystyle
D^{(2)}_{\mathcal{A}}=\big{(}\mathbf{\hat{\mathcal{T}}}-A(k)\big{)}^{-1}\Big{[}A(k)$
$\displaystyle+(A(k)-A(k_{1}))\frac{\bm{k}_{1}\cdot\bm{k}_{2}}{k_{2}^{2}}$
$\displaystyle+(A(k)-A(k_{2}))\frac{\bm{k}_{1}\cdot\bm{k}_{2}}{k_{1}^{2}}\Big{]}D_{+}(k_{1})D_{+}(k_{2}),$
(A.25) $\displaystyle
D^{(2)}_{\mathcal{B}}=\big{(}\mathbf{\hat{\mathcal{T}}}-A(k)\big{)}^{-1}\Big{[}A(k_{1})$
$\displaystyle+A(k_{2})-A(k)\Big{]}D_{+}(k_{1})D_{+}(k_{2}),$ (A.26)

with $k=|\bm{k}_{1}+\bm{k}_{2}|$, and [61]

$\mathbf{\hat{\mathcal{T}}}=\frac{d^{2}\,}{dt^{2}}+2H\frac{d\,}{dt}.$ (A.27)

From eqs. (A.22) and (A.23) one recovers eqs. (3.13) and (3.14), obtained from
LPT. And by inverting eqs. (3.4) and (3.8), one obtains the LPT kernels

$\displaystyle\Gamma_{2}(\bm{k}_{1},\bm{k}_{2})$
$\displaystyle=\mathcal{A}-\mathcal{B}\frac{(\bm{k}_{1}\cdot\bm{k}_{2})^{2}}{k_{1}^{2}k_{2}^{2}},$
(A.28) $\displaystyle\Gamma_{2}^{f}(\bm{k}_{1},\bm{k}_{2})$
$\displaystyle=\left(\mathcal{A}-\mathcal{B}\frac{(\bm{k}_{1}\cdot\bm{k}_{2})^{2}}{k_{1}^{2}k_{2}^{2}}\right)\frac{f(k_{1})+f(k_{2})}{2f_{0}}+\frac{1}{2Hf_{0}}\left(\dot{\mathcal{A}}-\dot{\mathcal{B}}\frac{(\bm{k}_{1}\cdot\bm{k}_{2})^{2}}{k_{1}^{2}k_{2}^{2}}\right),$
(A.29)

which are those found in [24] by a direct LPT approach.

## References

  * [1] Planck collaboration, N. Aghanim et al., _Planck 2018 results. I. Overview and the cosmological legacy of Planck_ , _Astron. Astrophys._ 641 (2020) A1, [1807.06205].
  * [2] eBOSS collaboration, S. Alam et al., _Completed SDSS-IV extended Baryon Oscillation Spectroscopic Survey: Cosmological implications from two decades of spectroscopic surveys at the Apache Point Observatory_ , _Phys. Rev. D_ 103 (2021) 083533, [2007.08991].
  * [3] DES collaboration, S. Pandey et al., _Dark Energy Survey Year 3 Results: Constraints on cosmological parameters and galaxy bias models from galaxy clustering and galaxy-galaxy lensing using the redMaGiC sample_ , 2105.13545.
  * [4] Super-Kamiokande collaboration, Y. Fukuda et al., _Evidence for oscillation of atmospheric neutrinos_ , _Phys. Rev. Lett._ 81 (1998) 1562–1567, [hep-ex/9807003].
  * [5] SNO collaboration, Q. R. Ahmad et al., _Measurement of the rate of $\nu_{e}+d\to p+p+e^{-}$ interactions produced by 8B solar neutrinos at the Sudbury Neutrino Observatory_, _Phys. Rev. Lett._ 87 (2001) 071301, [nucl-ex/0106015].
  * [6] SNO collaboration, Q. R. Ahmad et al., _Direct evidence for neutrino flavor transformation from neutral current interactions in the Sudbury Neutrino Observatory_ , _Phys. Rev. Lett._ 89 (2002) 011301, [nucl-ex/0204008].
  * [7] I. Esteban, M. Gonzalez-Garcia, A. Hernandez-Cabezudo, M. Maltoni and T. Schwetz, _Global analysis of three-flavour neutrino oscillations: synergies and tensions in the determination of $\theta_{23}$, $\delta_{CP}$, and the mass ordering_, _JHEP_ 01 (2019) 106, [1811.05487].
  * [8] KATRIN collaboration, M. Aker et al., _Improved Upper Limit on the Neutrino Mass from a Direct Kinematic Method by KATRIN_ , _Phys. Rev. Lett._ 123 (2019) 221802, [1909.06048].
  * [9] M. Aker et al., _First direct neutrino-mass measurement with sub-eV sensitivity_ , 2105.08533.
  * [10] Planck collaboration, N. Aghanim et al., _Planck 2018 results. VI. Cosmological parameters_ , _Astron. Astrophys._ 641 (2020) A6, [1807.06209].
  * [11] M. M. Ivanov, M. Simonović and M. Zaldarriaga, _Cosmological Parameters from the BOSS Galaxy Power Spectrum_ , _JCAP_ 05 (2020) 042, [1909.05277].
  * [12] N. Palanque-Delabrouille, C. Yèche, N. Schöneberg, J. Lesgourgues, M. Walther, S. Chabanier et al., _Hints, neutrino bounds and WDM constraints from SDSS DR14 Lyman- $\alpha$ and Planck full-survey data_, _JCAP_ 04 (2020) 038, [1911.09073].
  * [13] Z. Slepian and S. K. N. Portillo, _Too hot to handle? Analytic solutions for massive neutrino or warm dark matter cosmologies_ , _Mon. Not. Roy. Astron. Soc._ 478 (2018) 516–529, [1710.01785].
  * [14] M. Shoji and E. Komatsu, _Massive Neutrinos in Cosmology: Analytic Solutions and Fluid Approximation_ , _Phys. Rev. D_ 81 (2010) 123516, [1003.0942].
  * [15] J. Lesgourgues and S. Pastor, _Massive neutrinos and cosmology_ , _Phys. Rept._ 429 (2006) 307–379, [astro-ph/0603494].
  * [16] C. Hahn, F. Villaescusa-Navarro, E. Castorina and R. Scoccimarro, _Constraining $M_{\nu}$ with the bispectrum. Part I. Breaking parameter degeneracies_, _JCAP_ 03 (2020) 040, [1909.11107].
  * [17] C. Hahn and F. Villaescusa-Navarro, _Constraining $M_{\nu}$ with the bispectrum. Part II. The information content of the galaxy bispectrum monopole_, _JCAP_ 04 (2021) 029, [2012.02200].
  * [18] F. Kamalinejad and Z. Slepian, _A Non-Degenerate Neutrino Mass Signature in the Galaxy Bispectrum_ , 2011.00899.
  * [19] DESI collaboration, A. Aghamousa et al., _The DESI Experiment Part I: Science,Targeting, and Survey Design_ , 1611.00036.
  * [20] R. Laureijs et al., _Euclid definition study report_ , 1110.3193.
  * [21] LSST collaboration, v. Ivezić et al., _LSST: from Science Drivers to Reference Design and Anticipated Data Products_ , _Astrophys. J._ 873 (2019) 111, [0805.2366].
  * [22] G. D’Amico, J. Gleyzes, N. Kokron, K. Markovic, L. Senatore, P. Zhang et al., _The Cosmological Analysis of the SDSS/BOSS data from the Effective Field Theory of Large-Scale Structure_ , _JCAP_ 05 (2020) 005, [1909.05271].
  * [23] T. Colas, G. D’amico, L. Senatore, P. Zhang and F. Beutler, _Efficient Cosmological Analysis of the SDSS/BOSS data from the Effective Field Theory of Large-Scale Structure_ , _JCAP_ 06 (2020) 001, [1909.07951].
  * [24] A. Aviles and A. Banerjee, _A Lagrangian Perturbation Theory in the presence of massive neutrinos_ , _JCAP_ 10 (2020) 034, [2007.06508].
  * [25] A. Aviles, G. Valogiannis, M. A. Rodriguez-Meza, J. L. Cervantes-Cota, B. Li and R. Bean, _Redshift space power spectrum beyond Einstein-de Sitter kernels_ , _JCAP_ 04 (2021) 039, [2012.05077].
  * [26] Z. Vlah and M. White, _Exploring redshift-space distortions in large-scale structure_ , _JCAP_ 1903 (2019) 007, [1812.02775].
  * [27] L. Senatore and M. Zaldarriaga, _The IR-resummed Effective Field Theory of Large Scale Structures_ , _JCAP_ 02 (2015) 013, [1404.5954].
  * [28] T. Baldauf, M. Mirbabayi, M. Simonović and M. Zaldarriaga, _Equivalence Principle and the Baryon Acoustic Peak_ , _Phys. Rev. D_ 92 (2015) 043514, [1504.04366].
  * [29] M. M. Ivanov and S. Sibiryakov, _Infrared Resummation for Biased Tracers in Redshift Space_ , _JCAP_ 07 (2018) 053, [1804.05080].
  * [30] P. McDonald and A. Roy, _Clustering of dark matter tracers: generalizing bias for the coming era of precision LSS_ , _JCAP_ 0908 (2009) 020, [0902.0991].
  * [31] J. E. McEwen, X. Fang, C. M. Hirata and J. A. Blazek, _FAST-PT: a novel algorithm to calculate convolution integrals in cosmological perturbation theory_ , _JCAP_ 09 (2016) 015, [1603.04826].
  * [32] M. Schmittfull, Z. Vlah and P. McDonald, _Fast large scale structure perturbation theory using one-dimensional fast Fourier transforms_ , _Phys. Rev. D_ 93 (2016) 103528, [1603.04405].
  * [33] M. Simonović, T. Baldauf, M. Zaldarriaga, J. J. Carrasco and J. A. Kollmeier, _Cosmological perturbation theory using the FFTLog: formalism and connection to QFT loop integrals_ , _JCAP_ 04 (2018) 030, [1708.08130].
  * [34] J. Brandbyge and S. Hannestad, _Grid based linear neutrino perturbations in cosmologicaln-body simulations_ , _Journal of Cosmology and Astroparticle Physics_ 2009 (May, 2009) 002–002.
  * [35] M. Viel, M. G. Haehnelt and V. Springel, _The effect of neutrinos on the matter distribution as probed by the intergalactic medium_ , _Journal of Cosmology and Astroparticle Physics_ 2010 (Jun, 2010) 015–015.
  * [36] Y. Ali-Haïmoud and S. Bird, _An efficient implementation of massive neutrinos in non-linear structure formation simulations_ , _Monthly Notices of the Royal Astronomical Society_ 428 (Nov, 2012) 3375–3389.
  * [37] A. Banerjee and N. Dalal, _Simulating nonlinear cosmological structure formation with massive neutrinos_ , _JCAP_ 11 (2016) 015, [1606.06167].
  * [38] D. Inman and U.-L. Pen, _Cosmic neutrinos: A dispersive and nonlinear fluid_ , _Physical Review D_ 95 (Mar, 2017) .
  * [39] A. Banerjee, D. Powell, T. Abel and F. Villaescusa-Navarro, _Reducing noise in cosmological n-body simulations with neutrinos_ , _Journal of Cosmology and Astroparticle Physics_ 2018 (Sep, 2018) 028–028.
  * [40] S. Bird, Y. Ali-Haïmoud, Y. Feng and J. Liu, _An efficient and accurate hybrid method for simulating non-linear neutrino structure_ , _Monthly Notices of the Royal Astronomical Society_ 481 (08, 2018) 1486–1500, [https://academic.oup.com/mnras/article-pdf/481/2/1486/25716087/sty2376.pdf].
  * [41] J. Dakin, J. Brandbyge, S. Hannestad, T. HaugbØlle and T. Tram, _$\nu$ concept: cosmological neutrino simulations from the non-linear boltzmann hierarchy_, _Journal of Cosmology and Astroparticle Physics_ 2019 (Feb, 2019) 052–052.
  * [42] A. E. Bayer, A. Banerjee and Y. Feng, _A fast particle-mesh simulation of non-linear cosmological structure formation with massive neutrinos_ , _JCAP_ 01 (2021) 016, [2007.13394].
  * [43] F. Villaescusa-Navarro et al., _The Quijote simulations_ , 1909.05273.
  * [44] W. Hu and D. J. Eisenstein, _Small scale perturbations in a general MDM cosmology_ , _Astrophys. J._ 498 (1998) 497, [astro-ph/9710216].
  * [45] S. Saito, M. Takada and A. Taruya, _Impact of massive neutrinos on nonlinear matter power spectrum_ , _Phys. Rev. Lett._ 100 (2008) 191301, [0801.0607].
  * [46] Y. Y. Y. Wong, _Higher order corrections to the large scale matter power spectrum in the presence of massive neutrinos_ , _JCAP_ 0810 (2008) 035, [0809.0693].
  * [47] S. Saito, M. Takada and A. Taruya, _Nonlinear power spectrum in the presence of massive neutrinos: perturbation theory approach, galaxy bias and parameter forecasts_ , _Phys. Rev._ D80 (2009) 083528, [0907.2922].
  * [48] M. Shoji and E. Komatsu, _Third-order Perturbation Theory With Non-linear Pressure_ , _Astrophys. J._ 700 (2009) 705–719, [0903.2669].
  * [49] J. Lesgourgues, S. Matarrese, M. Pietroni and A. Riotto, _Non-linear Power Spectrum including Massive Neutrinos: the Time-RG Flow Approach_ , _JCAP_ 0906 (2009) 017, [0901.4550].
  * [50] A. Upadhye, R. Biswas, A. Pope, K. Heitmann, S. Habib, H. Finkel et al., _Large-Scale Structure Formation with Massive Neutrinos and Dynamical Dark Energy_ , _Phys. Rev._ D89 (2014) 103515, [1309.5872].
  * [51] H. Dupuy and F. Bernardeau, _Describing massive neutrinos in cosmology as a collection of independent flows_ , _JCAP_ 01 (2014) 030, [1311.5487].
  * [52] D. Blas, M. Garny, T. Konstandin and J. Lesgourgues, _Structure formation with massive neutrinos: going beyond linear theory_ , _JCAP_ 1411 (2014) 039, [1408.2995].
  * [53] F. Führer and Y. Y. Y. Wong, _Higher-order massive neutrino perturbations in large-scale structure_ , _JCAP_ 03 (2015) 046, [1412.2764].
  * [54] M. Levi and Z. Vlah, _Massive neutrinos in nonlinear large scale structure: A consistent perturbation theory_ , 1605.09417.
  * [55] B. S. Wright, H. A. Winther and K. Koyama, _COLA with massive neutrinos_ , _JCAP_ 10 (2017) 054, [1705.08165].
  * [56] L. Senatore and M. Zaldarriaga, _The Effective Field Theory of Large-Scale Structure in the presence of Massive Neutrinos_ , 1707.04698.
  * [57] M. Garny and P. Taule, _Loop corrections to the power spectrum for massive neutrino cosmologies with full time- and scale-dependence_ , _JCAP_ 01 (2021) 020, [2008.00013].
  * [58] N. E. Chisari et al., _Modelling baryonic feedback for survey cosmology_ , _Open J. Astrophys._ 2 (2019) 4, [1905.06082].
  * [59] A. Lewis, A. Challinor and A. Lasenby, _Efficient computation of CMB anisotropies in closed FRW models_ , _Astrophys. J._ 538 (2000) 473–476, [astro-ph/9911177].
  * [60] A. Aviles and J. L. Cervantes-Cota, _Lagrangian perturbation theory for modified gravity_ , _Phys. Rev._ D96 (2017) 123526, [1705.10719].
  * [61] T. Matsubara, _Recursive Solutions of Lagrangian Perturbation Theory_ , _Phys. Rev._ D92 (2015) 023534, [1505.01481].
  * [62] G. Valogiannis, R. Bean and A. Aviles, _An accurate perturbative approach to redshift space clustering of biased tracers in modified gravity_ , _JCAP_ 2001 (2020) 055, [1909.05261].
  * [63] T. Matsubara, _Resumming Cosmological Perturbations via the Lagrangian Picture: One-loop Results in Real Space and in Redshift Space_ , _Phys. Rev. D_ 77 (2008) 063530, [0711.2521].
  * [64] F. Bernardeau, S. Colombi, E. Gaztanaga and R. Scoccimarro, _Large scale structure of the universe and cosmological perturbation theory_ , _Phys. Rept._ 367 (2002) 1–248, [astro-ph/0112551].
  * [65] E. Castorina, C. Carbone, J. Bel, E. Sefusatti and K. Dolag, _DEMNUni: The clustering of large-scale structures in the presence of massive neutrinos_ , _JCAP_ 07 (2015) 043, [1505.07148].
  * [66] A. Aviles, _Dark matter dispersion tensor in perturbation theory_ , _Phys. Rev. D_ 93 (2016) 063517, [1512.07198].
  * [67] R. de Belsunce and L. Senatore, _Tree-Level Bispectrum in the Effective Field Theory of Large-Scale Structure extended to Massive Neutrinos_ , _JCAP_ 02 (2019) 038, [1804.06849].
  * [68] F. Villaescusa-Navarro, F. Marulli, M. Viel, E. Branchini, E. Castorina, E. Sefusatti et al., _Cosmology with massive neutrinos I: towards a realistic modeling of the relation between matter, haloes and galaxies_ , _JCAP_ 03 (2014) 011, [1311.0866].
  * [69] E. Castorina, E. Sefusatti, R. K. Sheth, F. Villaescusa-Navarro and M. Viel, _Cosmology with massive neutrinos II: on the universality of the halo mass function and bias_ , _JCAP_ 02 (2014) 049, [1311.1212].
  * [70] S. Vagnozzi, T. Brinckmann, M. Archidiacono, K. Freese, M. Gerbino, J. Lesgourgues et al., _Bias due to neutrinos must not uncorrect’d go_ , _JCAP_ 09 (2018) 001, [1807.04672].
  * [71] A. Banerjee, E. Castorina, F. Villaescusa-Navarro, T. Court and M. Viel, _Weighing neutrinos with the halo environment_ , _JCAP_ 06 (2020) 032, [1907.06598].
  * [72] M. LoVerde, _Halo bias in mixed dark matter cosmologies_ , _Phys. Rev. D_ 90 (2014) 083530, [1405.4855].
  * [73] M. LoVerde, _Spherical collapse in $\nu\Lambda$CDM_, _Phys. Rev. D_ 90 (2014) 083518, [1405.4858].
  * [74] Y. Donath and L. Senatore, _Biased Tracers in Redshift Space in the EFTofLSS with exact time dependence_ , _JCAP_ 10 (2020) 039, [2005.04805].