# Redshift-Space Distortions in Lagrangian Perturbation Theory

Shi-Fan Chen  Zvonimir Vlah  Emanuele Castorina  Martin White

###### Abstract

We present the one-loop 2-point function of biased tracers in redshift space
computed with Lagrangian perturbation theory, including a full resummation of
both long-wavelength (infrared) displacements and associated velocities. The
resulting model accurately predicts the power spectrum and correlation
function of halos and mock galaxies from two different sets of N-body
simulations at the percent level for quasi-linear scales, including the
damping of the baryon acoustic oscillation signal due to the bulk motions of
galaxies. We compare this full resummation with other, approximate, techniques
including the moment expansion and Gaussian streaming model. We discuss
infrared resummation in detail and compare our Lagrangian formulation with the
Eulerian theory augmented by an infrared resummation based on splitting the
input power spectrum into “wiggle” and “no-wiggle” components. We show that
our model is able to recover unbiased cosmological parameters in mock data
encompassing a volume much larger than what will be available to future galaxy
surveys. We demonstrate how to efficiently compute the resulting expressions
numerically, making available a fast Python code capable of rapidly computing
these statistics in both configuration and Fourier space.

##  1 Introduction

The measured redshifts of galaxies receive a contribution proportional to the
relative velocity between the observer and the emitting source. This term
breaks the isotropy of space, as it singles out the observer location as
special, and causes a distinct anisotropic pattern in the clustering
statistics of biased tracers known as Redshift Space Distortions (RSD) [1, 2,
3, 4]. Being a probe of the velocity field, RSD contain extra cosmological
information compared to the density field only, and they have been shown to be
a useful probe of modified gravity models [5, 6, 7, 8, 9]. Current and
upcoming spectroscopic redshift surveys, like DESI [10] and Euclid [11], will
provide measurements of the power spectrum of galaxies with much better
precision than currently available, making the modeling of RSD of paramount
importance to achieve their science goals. Within the framework of
Perturbation Theory (PT) several different approaches have been put forward to
compute clustering statistics in redshift space. We can divide them into two
main categories: Eulerian PT (EPT) methods, where density and velocity fields
are the relevant degrees of freedom (dof), and Lagrangian PT (LPT) methods,
where the displacements of dark matter particles and galaxies are the
fundamental dof from which observables are computed111In both cases we will
only consider the effective field theory (EFT) approach to PT, see [12, 13,
14, 15, 16, 17] and references therein.. This work focuses on the latter, but
more generally one of our main goals is to clarify the relation between the
two approaches to RSD.

LPT has a long history, especially in the context of RSD [18, 19, 20, 21, 22,
23, 24], and LPT models of the galaxy correlation function, i.e. in position
space, have been successfully applied to data (see e.g. refs. [25, 26, 27, 28,
29, 30] for a sampling of the literature). In Fourier space, while RSD in EPT
can be straightforwardly implemented [31, 32, 33, 34, 35, 36], LPT has posed a
number of technical difficulties that have only recently been overcome [23].
Using an expansion in the moments of the density-weighted pairwise velocities,
ref. [24] presented a derivation of the $1$-loop redshift space galaxy power
spectrum in LPT, finding good agreement in comparison to simulated data.

The main goal of this paper is to compute the $1$-loop redshift space power
spectrum in LPT via direct evaluation of the integrals, without employing the
moment expansion (MOME) [37, 38, 39, 23]. This is not only a practical choice,
but, for example, it will allow us to clarify the effect on the power spectrum
of various resummation schemes for the long-wavelength (IR) displacement modes
[40, 41, 15, 42, 43, 44, 45, 46].

In EPT, IR-resummation is performed a posteriori, after the $1$-loop power
spectrum is computed, and different procedures have been discussed in the
literature. In redshift space in particular, a number of approximations have
been employed to render the calculation more tractable [47, 48, 45, 46]. In
the current implementation of MOME in LPT [24], only the long wavelength
displacements are resummed, leaving the long wavelength velocity effects un-
resummed. If desired, these long wavelength velocity effects can also be
resummed in a posteriori way, as is done in EPT, and as we show in Appendix A.
Alternatively, by truncating the configuration-space velocity cumulants at
second order, the Gaussian streaming model can be used to approximately resum
velocities from linear (Gaussian) modes at the expense of neglecting some one-
loop contributions222See, e.g. Appendix B of ref. [49]. [22].

The direct LPT approach allows us to efficiently resum the long displacement
contributions without relying on any of the above-mentioned approximations.
Some of these different choices for IR-resummation can lead to different
behaviours of the power spectra at small (UV) scales. However, as we shall
discuss further below, these differences can be associated to the different
perturbative expansion parameters they employ. Moreover, even though these
differences arise from long wavelength displacements, which are under
perturbative control, these residual contributions are also approximately
degenerate with, and can thus be absorbed by, the free coefficients of the
effective theory.

The results presented in this paper complement the existing literature on the
one-loop LPT power spectrum in $\Lambda$CDM cosmologies referenced above.
Additionally, our work provides the analytical machinery to better understand
the performance of forward model or density field reconstruction algorithms in
redshift space based on Lagrangian displacements [50, 51, 52].

This paper is structured as follows. Section 2 will introduce the notation,
Section 3 the relevant equations for the computation of the one-loop RSD power
spectrum of biased tracers, whose evaluation is discussed in Section 4.
Numerical fits to N-body simulations and mock catalogs are presented in
Section 5, as well as comparisons between the different PT methods. This
section also demonstrates that the model is able to recover unbiased estimates
of cosmological parameters in a “blind” challenge. Our conclusions are
presented in Section 6. A number of technical points are relegated to
Appendices A, B and C.

##  2 Overview of Lagrangian Perturbation Theory

Our goal in this section is to give a quick overview of Lagrangian
perturbation theory (LPT) as pertains to this paper, both as a review and to
establish our notation and conventions. The reader is referred to the
references in the introduction for further details, and especially to refs.
[18, 53, 19, 21, 15, 22, 49] whose notations we adopt.

Within the Lagrangian picture the gravitational evolution of large-scale
structure is described via the of motion fluid elements starting at initial
(Lagrangian) positions q with trajectories given by
$\textbf{x}(\textbf{q},t)=\textbf{q}+\boldsymbol{\Psi}(\textbf{q},t)$. The
Lagrangian displacements, $\boldsymbol{\Psi}$, obey the equation of motion
$\ddot{\boldsymbol{\Psi}}(\textbf{q})+\mathcal{H}\dot{\boldsymbol{\Psi}}(\textbf{q})=-\nabla_{\textbf{x}}\Phi(\textbf{x})$,
where dots indicate derivatives with respect to the conformal time, and the
gravitational potential $\Phi$ is in turn sourced by the matter overdensity
$\delta_{m}$ given by

$\displaystyle 1+\delta_{m}(\textbf{x})=\int d^{3}\textbf{q}\
\delta_{D}(\textbf{x}-\textbf{q}-\boldsymbol{\Psi}(\textbf{q}))\quad,\quad(2\pi)^{3}\delta_{D}(\textbf{k})+\tilde{\delta}_{m}(\textbf{k})=\int
d^{3}\textbf{q}\ e^{i\textbf{k}\cdot(\textbf{q}+\boldsymbol{\Psi})}$ (2.1)

via Poisson’s equation. In LPT these quantities are solved for order-by-order
in the initial conditions $\delta_{0}(\textbf{q})$, such that the
displacements are given by
$\boldsymbol{\Psi}=\boldsymbol{\Psi}^{(1)}+\boldsymbol{\Psi}^{(2)}+\boldsymbol{\Psi}^{(3)}+...$
Of particular interest is the linear solution
$\boldsymbol{\Psi}^{(1)}=-D(z)\nabla^{-1}_{\textbf{q}}\delta_{0}(\textbf{q})$,
also known as the Zeldovich approximation. The specific forms of the higher-
order solutions are given for example in refs. [54, 55, 56]. These solutions
contain parametrizeable dependences on small-scale physics which are captured
by including additional effective-theory counterterms [14, 15].

In this paper we will be primarily interested in the clustering of biased
tracers of matter like galaxies which are the target of galaxy redshift
surveys. Within the Lagrangian framework biased tracers are modeled as
functionals of the initial conditions $F[\delta_{0}](\textbf{q})$ at their
Lagrangian positions q and advected along with the matter fluid, such that
their observed overdensities are given by number conservation to be

$1+\delta_{g}(\textbf{x})=\int d^{3}\textbf{q}\ F(\textbf{q})\
\delta_{D}(\textbf{x}-\textbf{q}-\boldsymbol{\Psi}).$ (2.2)

The bias functional $F(\textbf{q})$ is a local function of the initial
conditions with effective corrections, and at one-loop order in the power
spectrum includes linear and quadratic density bias, shear and third-order
contributions as well as effective corrections like derivative bias
($\propto\nabla^{2}\delta^{2}$) similar to the dynamical ones described at the
end of the previous paragraph. We will furthermore assume that cold dark
matter and baryons can be treated as a single fluid and therefore discard bias
operators and dynamics proportional to relative density and velocity of the
different fluids [17, 57, 58, 59, 60, 61, 62]. These extra terms are thought
to be small and therefore need only be implemented at leading order in the
(relative) displacements, as discussed in ref. [59]. Our conventions follow
Equation 5.1 in ref. [49].

Finally, galaxy surveys determine the line-of-sight (LOS) position of observed
galaxies via redshifts $z$ whose cosmological and peculiar-velocity
contributions are degenerate. These redshift-space distortions (RSD) can be
accounted for within LPT by boosting displacements along the LOS direction
$\hat{n}$ by their corresponding velocities
$\boldsymbol{\Psi}^{s}=\boldsymbol{\Psi}+(\hat{n}\cdot\textbf{v})\hat{n}/\mathcal{H}$,
where $\mathcal{H}$ is the conformal Hubble parameter. We will use the
superscript $s$ to refer to vectors boosted into redshift space throughout
this work. Within the Einstein-de Sitter approximation (EdS) we have the
further simplification that

$\boldsymbol{\Psi}^{s,(n)}=\boldsymbol{\Psi}^{(n)}+nf\
(\hat{n}\cdot\boldsymbol{\Psi}^{(n)})\ \hat{n}\equiv
R^{(n)}\boldsymbol{\Psi}^{(n)},$ (2.3)

where the matrix $R^{(n)}_{ij}=\delta_{ij}+nf\ \hat{n}_{i}\hat{n}_{j}$ and $f$
is the linear growth rate. We will operate within the EdS approximation for
the remainder of the paper, and further make the distant observer
approximation such that $\hat{n}$ is the same for each galaxy. These
approximations are known to be quite good in the limit of high redshifts and
on scales where higher order perturbation theory is most applicable. A
discussion of violations of these approximations within the LPT context can be
found in refs. [63, 64, 65, 66, 67, 68].

##  3 Redshift-Space Power Spectrum

We now proceed to write down the power spectrum at one loop in Lagrangian
perturbation theory [53, 19, 20, 21, 22, 23, 49]. From Equation 2.1 and its
counterpart for biased tracers we have that the galaxy autospectrum is given
by

$P(\textbf{k})=\int d^{3}\textbf{q}\
e^{i\textbf{k}\cdot\textbf{q}}\left\langle
e^{i\textbf{k}\cdot\Delta}F(\textbf{q}_{1})F(\textbf{q}_{2})\right\rangle_{\textbf{q}=\textbf{q}_{1}-\textbf{q}_{2}}$

where we have defined the pairwise Lagrangian displacement
$\Delta_{i}=\boldsymbol{\Psi}_{i}(\textbf{q}_{1})-\boldsymbol{\Psi}_{i}(\textbf{q}_{2})$.
Setting $F=1$, for matter, the bracketed average can expressed using the
cumulant theorem as [19]

$\ln\left\langle
e^{i\textbf{k}\cdot\Delta}\right\rangle=-\frac{1}{2}k_{i}k_{j}A_{ij}-\frac{i}{6}k_{i}k_{j}k_{k}W_{ijk}+...,$

where we have defined the cumulants of the pairwise displacements as
$A_{ij}=\left\langle\Delta_{i}\Delta_{j}\right\rangle_{c}$ and
$W_{ijk}=\left\langle\Delta_{i}\Delta_{j}\Delta_{k}\right\rangle_{c}$. For
biased tracers one simply needs to compute cumulants with sources like
$J(\textbf{q})\delta(\textbf{q})$ added to the exponent and take functional
derivatives; this produces terms like
$U_{i}=\left\langle\delta_{0}(\textbf{q}_{1})\Delta_{i}\right\rangle.$ The
above calculations can be promoted to redshift space by promoting the
displacements to redshift space ($\Delta^{s}$) order-by-order as in Equation
2.3.

From the above, the one-loop galaxy autospectrum in redshift space is given in
LPT by [53, 19, 20, 21, 22, 23, 49]

$\displaystyle P_{s}(\textbf{k})$ $\displaystyle=\int d^{3}\textbf{q}\
e^{i\textbf{k}\cdot\textbf{q}}\
e^{-\frac{1}{2}k_{i}k_{j}A^{s,<}_{ij}}\Big{\\{}1-\frac{1}{2}k_{i}k_{j}A^{s,>}_{ij}+\frac{1}{8}k_{i}k_{j}k_{k}k_{l}A^{s,>}_{ij}A^{s,>}_{kl}$
$\displaystyle-\frac{1}{2}k_{i}k_{j}A^{s,\rm
loop}_{ij}+\frac{i}{6}k_{i}k_{j}k_{k}W^{s}_{ijk}$
$\displaystyle+2ib_{1}k_{i}(1-\frac{1}{2}k_{i}k_{j}A^{s,>}_{ij})U^{s}_{i}-b_{1}k_{i}k_{j}A^{s,10}_{ij}$
$\displaystyle+b_{1}^{2}(1-\frac{1}{2}k_{i}k_{j}A^{s,>}_{ij})\xi_{\rm
lin}+ib_{1}^{2}k_{i}U^{s,11}_{i}-b_{1}^{2}k_{i}k_{j}U^{s,\rm lin}_{i}U^{s,\rm
lin}_{j}$ $\displaystyle+\frac{1}{2}b_{2}^{2}\xi_{\rm
lin}^{2}+2ib_{1}b_{2}\xi_{\rm lin}k_{i}U^{s,\rm
lin}_{i}-b_{2}k_{i}k_{j}U^{s,\rm lin}_{i}U^{s,\rm
lin}_{j}+ib_{2}k_{i}U^{s,20}_{i}$
$\displaystyle+b_{s}(-k_{i}k_{j}\Upsilon^{s}_{ij}+2ik_{i}V^{s,10}_{i})+2ik_{i}b_{1}b_{s}V^{s,12}_{i}+b_{2}b_{s}\chi+b_{s}^{2}\zeta+2ib_{3}k_{i}U^{s}_{b_{3},i}+2b_{1}b_{3}\theta+...\Big{\\}}$
$\displaystyle+k^{2}(\alpha_{0}+\alpha_{2}\mu^{2}+\alpha_{4}\mu^{4}+\alpha_{6}\mu^{6})P_{s,\rm
Zel}(\textbf{k})+R_{h}^{3}(1+\sigma_{2}k^{2}\mu^{2}+\sigma_{4}k^{4}\mu^{4}).$
(3.1)

Here we have split the second cumulant $A_{ij}$ into long and short linear
components and a loop component, keeping only the long-wavelength piece
$A^{s,<}_{ij}$ exponentiated; we will comment on this further below. In
addition to $A_{ij}$, $W_{ijk}$ and $U_{i}$ defined above, Equation 3.1
contains additional correlators of pairwise Lagrangian displacements and
higher-order bias operators like shear that are defined explicitly in refs.
[19, 21, 22, 49]. The last line of Equation 3.1 includes counterterms
($\alpha_{n}$) and stochastic contributions ($\sigma_{n}$) proportional to the
typical scale of halo/galaxy formation $R_{h}$. These include what are
traditionally referred to as the “shot noise” and “finger of god” (FoG) terms.
The small-scale sensitivities that give rise to these terms are described in
detail333In particular, the advantages of this form for treating fingers of
god (and redshift errors) are discussed in detail in §§4.1.3, 4.2.3, 5.3 and
Appendix C of ref. [49]. in ref. [49].

While the correlators in Equation 3.1 have been extensively described
elsewhere, it is instructive to elucidate their general perturbative and
angular structure with an example. Let us consider the displacement two-point
function up to one-loop

$A_{ij}\equiv\left\langle\Delta_{i}\Delta_{j}\right\rangle=A^{\rm
lin}_{ij}+A^{\rm loop}_{ij}\quad,\quad A^{\rm
loop}_{ij}=A^{(22)}_{ij}+2A^{(13)}_{ij},$ (3.2)

which is given by a linear piece (lin) from contracting two first-order
displacements and a one-loop piece from contracting two second-order
displacements $(22)$ or one first and one third-order displacement each
$(13)$. To go into redshift space, each of these pieces must be transformed
separately — this is because at each order in perturbation theory the
translation to redshift space depends on $n$, such that for example

$A^{s,(13)}_{ij}=R^{(1)}_{in}R^{(3)}_{im}A^{(13)}_{nm}\quad\mathrm{but}\quad
A^{s,(22)}_{ij}=R^{(2)}_{in}R^{(2)}_{im}A^{(22)}_{nm}\quad.$

The $n$-dependence of these transformations encodes information about beyond-
linear velocities within the RSD spectrum and is the primary complication in
extending the treatment of RSD beyond the Zeldovich approximation, where all
vectors transform via $\textbf{R}=\textbf{R}^{(1)}$ (see refs. [21, 23, 24]).
A similar observation applies to all correlators with vector indices in
Equation 3.1.

Finally, let us comment on our resummation of the linear piece of $A_{ij}$. A
salient feature of cosmologies like $\Lambda$CDM is that large-scale
displacements produce nonlinear damping of spatially localized features in the
power spectrum such as baryon acoustic oscillations (BAO) that cannot be
captured simply with an order-by-order expansion in the linear initial
conditions [69, 70, 18, 53, 71, 72, 73, 19, 74, 75, 76, 40, 77, 42, 43, 78,
44, 79]. This is because, while the dynamics on these large scales are
essentially linear, the size of these displacements on BAO scales can be large
compared to wavenumbers where the BAO wiggles have support. As such, the
effects of these displacements must be manually resummed in order-by-order
expansions such as Eulerian perturbation theory (EPT). On the other hand,
within LPT the exponential in Equation 2.1 and the cumulant theorem for a
Gaussian variable

$\left\langle
e^{i\textbf{k}\cdot\Delta^{(1)}}\right\rangle=e^{-\frac{1}{2}k_{i}k_{j}A^{(11)}_{ij}}$
(3.3)

suggests a natural resummation scheme wherein the linear displacements are
kept exponentiated. These exponentiated displacements should be resummed only
up to an IR scale $k_{\rm IR}$; we will follow the convention in ref. [49] and
perform this split by an exponential cutoff $\exp[-\frac{1}{2}(k/k_{\rm
IR})^{2}]$ in the $A_{ij}$ integral for the “less than” displacements kept
resummed, $A^{s,<}_{ij}$ in Equation 3.1, and the “greater than” displacements
(defined with $1-\exp[-\frac{1}{2}(k/k_{\rm IR})^{2}]$ in the integral),
$A^{s,>}_{ij}$, which are expanded to second order in the curly brackets in
Equation 3.1. Within $\Lambda$CDM-like cosmologies, $A_{ij}$ is close to
saturated on BAO scales and it might be expected that an IR cutoff should make
only small differences in the final theory prediction. This is true for the
density statistics [76, 78, 22, 42, 46]. However, ref. [49] showed that the
higher-order velocity statistics that enter into RSD are especially sensitive
to IR resummation in both broadband and BAO wiggles. This is discussed in
further detail in Appendix A. We conclude this section by noting that,
contrary to EPT where the a posteriori IR-resummation has been implemented
only for equal time correlators, the direct evaluation of the LPT integrals
presented in this work automatically evaluates the power spectrum at unequal
times. In this case, the bulk displacement contributions do not cancel exactly
in the exponent in Eq.(3.3) leading to rapid suppression and decorrelation of
unequal time correlators. This has recently been considered [80] in the
context of weak lensing analyses, that all involve unequal time correlators.
On the same topic, this should also clarify some recent concerns raised in
ref. [81] about the use of perturbation theory for unequal time correlators.

##  4 Numerical Implementations

The primary challenge in evaluating the integral in Equation 3.1 lies in the
angular dependence due to the three vectors, q, k and $\hat{n}$, that enter
the calculation (Fig. 1a). By symmetry, each of the tensor-indexed Lagrangian-
space correlators in Equation 3.1 can be decomposed into components
multiplying products of $\hat{q}$ and the Kronecker delta symbol; for example,
we can write $A_{ij}(\textbf{q})=X(q)\delta_{ij}+Y(q)\hat{q}_{i}\hat{q}_{j}$
[19]. In real space, where the angular dependence is due only to k and q, it
is customary to proceed by defining a coordinate system wherein k points
towards the zenith such that the integrand has azimuthal symmetry. The
resulting dependence on $\mu_{\textbf{q}}=\hat{k}\cdot\hat{q}$, where the
subscript is meant to distinguish it from the familiar LOS angle
$\mu=\hat{k}\cdot\hat{n}$, can then be recast into infinite sums of spherical
Bessel functions using, for example, the identity

$\frac{1}{2}\int d\mu_{\textbf{q}}\
e^{iA\mu_{\textbf{q}}-\frac{1}{2}B\mu_{\textbf{q}}^{2}}=e^{-B/2}\sum_{n=0}\Big{(}\frac{B}{A}\Big{)}^{n}j_{n}(A)$
(4.1)

and its derivatives, with the resulting integrals in $q$ efficiently computed
using the FFTLog algorithm [82, 83]. In the particular case of LPT, we have
$A=kq$ and $B=k^{2}Y^{<}(q)$, such that in the Zeldovich power spectrum
(second line of Equation 3.1 in the $k_{\rm IR}=\infty$ limit) is

$P_{\rm Zel}(k)=4\pi\sum_{n=0}^{\infty}\int dq\,q^{2}\
e^{-\frac{1}{2}k^{2}(X+Y)}\Big{(}\frac{kY}{q}\Big{)}^{n}j_{n}(kq)$ (4.2)

in real space. This form of the integral, suitable for using fast Hankel
transforms, generalizes for higher loop terms, as well as when we go to
redshift space, as we show in the rest of this section. Note that the
computation calls for one transform for each $k$ value. The expansion
converges quickly for the $k$ values of our interest and it is typically
sufficient to keep only the $n<10$ terms in the sum above. Moreover, for
higher $n$ terms the Limber approximation [84, 85] can be used

$j_{\ell}(kq)\approx\sqrt{\frac{\pi}{2\ell+1}}\
\delta^{D}\left(kq-\ell-\frac{1}{2}\right),~{}~{}{\rm
when}~{}~{}\ell\to\infty,$ (4.3)

which provides an accurate approximation for the integral above when used for
$n>3$.

In redshift space, the azimuthal symmetry is broken by the line-of-sight
dependence, as shown in Figure 1a, which selects a preferred plane containing
$\hat{k}$ and $\hat{n}$ (blue). Below, we outline one method to efficiently
perform the integral in redshift space. Following previous work [23] we will
call it Method II. Our development extends the Zeldovich calculations for
matter and biased tracers in refs. [86, 23, 24] to one-loop order. Method II
relies on an active transformation of the wavevector k into a frame more
conducive to evaluating the integral in Equation 3.1. Of course, it is also
possible to directly evaluate the integral within the original frame — this is
the strategy of Method I. This alternative method is described for the
interested reader in Appendix B.

Figure 1: (a) Geometry of the vector and tensor quantities in the integral Eq.
3.1. In the absence of redshift-space distortions ($\hat{n}$ dependence) the
integral over q is azimuthally symmetric; with RSD, a $\phi$ dependence occurs
since $\hat{k}$ and $\hat{n}$ lie on a preferred plane. (b) In Method II,
vectors are boosted into the “Zeldovich” frame where projections along the
line of sight are amplified by the linear growth rate $f$ and the zenith is
redefined to be the thus-boosted $K_{i}=R_{ij}k_{j}$. In both frames,
$\hat{n}$, $\hat{k}$ and $\hat{K}$ are coplanar as shown in blue.

We now outline the rough strategy for Method II. In our expressions
$\boldsymbol{\Psi}$ enters only in the combination
$\textbf{k}\cdot\boldsymbol{\Psi}$. Thus instead of transforming all
displacements into redshift space via
$\boldsymbol{\Psi}^{(n)}\to\textbf{R}^{(n)}\boldsymbol{\Psi}^{(n)}$ we can
instead passively transform the wave vectors, multiplying them by
$\textbf{R}^{T}=\mathbf{R}$. In particular we shall apply the linear theory
transformation $\textbf{R}^{(1)}$ to k to yield the vector
$K_{i}=R_{ij}^{(1)}k_{j}$. We will deal with the non-linear contributions to
$\boldsymbol{\Psi}$ below. In this new “Zeldovich” frame, shown in Figure 1b,
the zenith direction is set to be $\hat{K}$, and we redefine
$\mu_{\textbf{q}}=\hat{K}\cdot\hat{q}$. Thus
$\mathbf{K}\cdot\hat{n}=k\mu(1+f)$ and $K^{2}=k^{2}[1+f(2+f)\mu^{2}]$. Note
that $\hat{k}$, $\hat{K}$ and $\hat{n}$ are coplanar.

Let us begin by reviewing the angular structure of this coordinate choice. We
have

$\displaystyle\hat{n}\cdot\hat{q}=A(\mu)\mu_{\textbf{q}}+B(\mu)\sqrt{1-\mu_{\textbf{q}}^{2}}\,\cos\phi,\quad$
$\displaystyle\textbf{k}\cdot\textbf{q}=kq\Big{(}c(\mu)\mu_{\textbf{q}}-s(\mu)\sqrt{1-\mu_{\textbf{q}}^{2}}\cos\phi\Big{)},$
(4.4)

with the definitions

$\displaystyle A(\mu)$
$\displaystyle=\frac{\mu(1+f)}{\sqrt{1+f(2+f)\mu^{2}}}\quad,\quad
B(\mu)=\sqrt{\frac{1-\mu^{2}}{1+f(2+f)\mu^{2}}}$ $\displaystyle c(\mu)$
$\displaystyle=\frac{1+f\mu^{2}}{\sqrt{1+f(2+f)\mu^{2}}}\quad,\quad
s(\mu)=\frac{f\mu\sqrt{1-\mu^{2}}}{\sqrt{1+f(2+f)\mu^{2}}}\quad.$ (4.5)

Note the square root in the denominators is simply $K/k$. That the azimuthal
dependence always multiplies the sine, $\sqrt{1-\mu_{\textbf{q}}^{2}}$, will
prove a particular convenience in this frame.

In terms of the above, the Zeldovich matter power spectrum can be succintly
expressed as

$\displaystyle P_{s}(\textbf{k})$ $\displaystyle=\int dq\ d\mu_{\textbf{q}}\
q^{2}\ e^{ikqc\mu_{\textbf{q}}-\frac{1}{2}K^{2}(X+Y\mu_{\textbf{q}}^{2})}\
\Bigg{(}\int d\phi\ e^{-ikqs\sqrt{1-\mu_{\textbf{q}}^{2}}\cos\phi}\Bigg{)}\ $

In ref. [23] this integral was shown to be expressable in terms of Bessel
functions via the identity

$I(A,B,C)=\int d\mu_{\textbf{q}}\ d\phi\
e^{-iC\sqrt{1-\mu_{\textbf{q}}^{2}}\cos\phi+iA\mu_{\textbf{q}}+B\mu_{\textbf{q}}^{2}}=4\pi
e^{B}\sum_{\ell=0}^{\infty}\Big{(}\frac{-2}{\rho}\Big{)}^{\ell}\tilde{G}^{(0)}_{0,\ell}(A,B,\rho)j_{\ell}(\rho)$
(4.6)

by substituting $A=kqc$, $B=-\frac{1}{2}K^{2}Y$ and $C=kqs$. The exact form of
the the kernel $\tilde{G}^{(0)}_{0,\ell}$ is given in Appendix C.4. Building
on top of this, any bias contribution involving only the linear (Zeldovich)
displacement (e.g. the linear bias term $ik_{i}U^{\rm lin}_{i}$), simply
acquires powers of $\mu_{\textbf{q}}$ and $K$ (e.g. $iK\mu_{\textbf{q}}$) that
can be evaluated as derivatives of the above with respect to $A$, since
correlators in Lagrangian space are always decomposable into $\delta_{ij}$ and
tensor products of $\hat{q}_{i}$ [24].

The simple Zeldovich angular structure above is, however, broken by the
inclusion of higher-order displacements. This is because these displacements
get boosted along the line of sight by more than linear theory when going to
redshift space. With $\mathbf{R}=\mathbf{R}^{(1)}$

$\dot{\boldsymbol{\Psi}}^{(n)}=\mathbf{R}^{(n)}\boldsymbol{\Psi}^{(n)}=(\mathbf{R}+(n-1)f\
\hat{n}\otimes\hat{n})\boldsymbol{\Psi}^{(n)}.$ (4.7)

In the spirit of the above calculations we can dot the matrix into the
wavevectors and take

$k_{i}\rightarrow K_{i}+f(n-1)k_{\parallel,i}\quad,\quad
k_{\parallel}=(k\mu)\hat{n}.$ (4.8)

Dotting the transformed wavevector with Lagrangian correlators thus simply
requires additional powers of $\hat{n}\cdot\hat{q}$, which conveniently
translates into powers of $\mu_{\textbf{q}}$ and
$\sqrt{1-\mu_{\textbf{q}}^{2}}\cos\phi$, i.e. the coefficients multiplying $A$
and $C$ in Equation 4.6. We thus see that any contribution to the power
spectrum can be evaluated via mixed ($A,C$) derivatives of Equation 4.6. We
refer the reader to Appendix C for further details and an example application
to the one-loop matter power spectrum.

##  5 Results

Having laid out how the one-loop power spectrum can be efficiently computed
within fully-resummed one-loop LPT, our goal in this section is to validate
our model against N-body data. In addition, we compare the performance of our
model with previous models such as the Gaussian streaming model (GSM; [87, 88,
89, 25, 20, 22]) and moment expansion (MOME; [23, 49]) in LPT and resummed
Eulerian perturbation theory (REPT) in both Fourier and configuration space.

###  5.1 Comparison to N-body

Figure 2: Fits to the redshift-space power spectrum wedges (top) and
multipoles (bottom) of the fiducial halo sample with
$10^{12.5}M_{\odot}<M<10^{13.0}M_{\odot}$ at $z=0.8.$ Both statistics were fit
assuming Gaussian covariances using a consistent set of bias parameters and
with linear displacements resummed up to $k_{\rm IR}=0.2\,h\,{\rm Mpc}^{-1}$.
The fiducial LPT model gives an excellent fit to the anisotropic power
spectrum inside the range of fit ($k<0.2\,h\,{\rm Mpc}^{-1}$) well within the
few-percent systematics expected from the N-body data. Shaded regions indicate
wavenumbers beyond the range of fit, with higher multipoles diverging faster
from the data past this point.

For our main comparisons to N-body data we use halo catalogs from the
simulations in ref. [90] and employ NbodyKit [91] to compute redshift-space
power spectrum wedges, multipoles and correlation function multipoles at
$z=0.8$. These simulations assume a $\Lambda$CDM cosmology with
$\Omega_{m}=0.2648$, $\Omega_{b}h^{2}=0.02258$, $h=0.71$, $n_{s}=0.963$ and
$\sigma_{8}=0.8$. We adopt the mass bin $12.5<\log(M/h^{-1}M_{\odot})<13.0$ as
our fiducial sample but have checked that we get similar results for a higher
mass bin as well as the mock galaxy sample described in ref. [49]. Our
fiducial sample has a number density of $\bar{n}=0.53\times
10^{-3}\,h^{3}\,{\rm Mpc}^{-3}$ and linear (Eulerian) bias of $b\approx 1.7$,
making it slightly sparser but about 40% more biased than the DESI ELG sample
at $z=0.85$ [10]. We have chosen these simulations due to their relatively
large total volume (4 boxes with volume $[4h^{-1}$Gpc$]^{3}$). With such a
large volume the statistical errors on the two-point functions will
necessarily be significantly smaller than galaxy surveys at comparable
redshifts; however, we caution that the use of “derated” time steps in the
running of these simulations may cause systematic errors on the few percent
level, as discussed in further detail in refs. [22, 49]. We have attempted to
mitigate this effect by using only the high redshift catalog at $z=0.8$.

Figure 2 compares our LPT model to the power spectrum wedges and multipoles of
our fiducial halo sample. We fit for the wedges, $P(k,\mu)$, up to $k_{\rm
max}=0.2\,h\,{\rm Mpc}^{-1}$ assuming Gaussian covariances. We use the same
parameters for the multipoles. We find bias parameters of order unity and the
isotropic stochastic contribution $R_{h}^{3}$ comparable to the shot noise,
noting that extending to higher (unperturbative) $k_{\rm max}$ tends to
recover apparently good fits with anomalously large bias and effective
parameters. Our model is in excellent agreement with the power spectrum wedges
at the scales shown, differing from the data at levels comparable to their
statistical uncertainty, with qualitatively similar behavior in the
multipoles, though the anisotropic contributions ($\ell>0$) diverge faster
than the monopoles as expected due to the enhanced nonlinearity of halo
velocities.

In Figure 2, as well as throughout the main body of this work, we have adopted
the fiducial choice of infrared cutoff $k_{\rm IR}=0.2\,h\,{\rm Mpc}^{-1}$. As
discussed in Section 3, compared to density statistics the velocity
statistics’ underlying redshift-space distortions have broadband shapes that
are especially sensitive to the choice of infared cutoff. For example, as
shown in ref. [49] the monopole and quadrupole of the second moment of the
pairwise velocity, responsible for contributions to the power spectrum
proportional to the growth rate ($f$) squared, respectively have broadband
shapes better captured by large and small $k_{\rm IR}$. One might thus hope to
find an intermediate regime wherein both statistics are reasonably captured,
and indeed we find that the choice $k_{\rm IR}=0.2\,h\,{\rm Mpc}^{-1}$
reproduces the hexadecapole better than either the fully-exponentiated limit
($k_{\rm IR}=\infty$) or $k_{\rm IR}=0.$ In principle, the spirit of
perturbation theory should demand that the expanded displacements
$k^{2}\Sigma^{2}_{>}$ be small444Here we define
$\Sigma^{2}_{<}=\frac{2}{3}\int\frac{dk}{2\pi^{2}}P_{\rm lin}(k)\
e^{-(k/k_{\rm
IR})^{2}},\quad\Sigma^{2}_{>}=\frac{2}{3}\int\frac{dk}{2\pi^{2}}P_{\rm
lin}(k)\ (1-e^{-(k/k_{\rm IR})^{2}}),$ (5.1) such that the sum
$\Sigma^{2}=\Sigma^{2}_{<}+\Sigma^{2}_{>}$ is the mean square pairwise
displacement of two distant points in the Zeldovich approximation.  while the
exponentiated ones kept manageable; for $k<0.2\,h\,{\rm Mpc}^{-1}$ this is
satisfied by our choice, though given that the total Zeldovich displacement
for the fiducial cosmology at $z=0.8$ is $\Sigma^{-1}\approx 0.2\,h\,{\rm
Mpc}^{-1}$ this is relatively insensitive to the choice of IR cutoff.
Moreover, while differences exist towards high $k$ and $\mu$, we find that in
general the small-scale differences between the theory’s predictions for
reasonable values of $k_{\rm IR}$ can largely be absorbed by the effective
parameters of the theory; further discussion of the interplay between $k_{\rm
IR}$ choice and our model’s predictions can be found in Appendix A; we intend
to return to this topic in greater depth in future work.

###  5.2 Comparison to Other Models in Fourier and Configuration Space

Figure 3: Redshift-space power spectrum multipoles of the fiducial halo sample
fit using three effective theory models: the fiducial LPT model, the
Lagrangian moment expansion and resummed Eulerian perturbation theory. All
three models are fit as in Figure 2 and are in excellent quantitative
agreement with the N-body data. The three models differ slightly in their
prediction for the hexadecapole broadband; we have explicitly tuned our LPT IR
resummation scheme to provide a good match to the data, though we note the
relatively large statistical uncertainty in the hexadecapole. Figure 4:
Configuration space correlation function multipoles predicted by our LPT
model, the Lagrangian moment expansion, Eulerian perturbation theory and the
Gaussian streaming model compared to N-body data. Each of the models are in
good agreement with the data within the few-percent systematic uncertainties
expected of the simulations, though we note that they all slightly overshoot
the dip around $80\ h^{-1}$ Mpc by around two percent. Note that due to the
high degree of similarity between the theory predictions many of the lines lie
on top of each other even in the fractional residuals in the bottom panel,
particularly when comparing LPT (solid) and MOME (dashed). Black dashed lines
in the lower panels indicate 2 and 5 percent errors.

The main difference between the model presented in this work and previous
effective-theory models of the redshift-space galaxy two-point function lies
in the IR-resummation procedure. Existing LPT formulations typically
incorporate bulk velocities either via streaming model resummations [22, 23,
49] or direct expansions of velocity statistics [23, 49]. While these
approaches have been shown to be sufficiently accurate to model redshift-space
distortions in a variety of contexts, they have nonetheless exhibited a number
of shortcomings. The Gaussian streaming model allows for a partial
resummation. However, the resummation procedure calls for nonlinear mapping of
all the loop contributions, even those that require counterterms in order to
regularize the UV dependence. Even though such mapping could in principle be
restricted to only long wavelength contributions, the model also exhibits a
somewhat cumbersome structure in Fourier space. The moment-expansion approach,
on the other hand, relies on an explicit expansion in the velocity moments.
However, in a similar way to the long wavelength displacement contributions,
long wavelength velocity contributions also affect the BAO feature in a manner
that then prompts the additional resummation of these contributions. In the
Eulerian approach this can be done a posteriori using an ad hoc, wiggle-no-
wiggle splitting of the power spectrum and resumming only the contributions
related to the BAO feature (REPT). Such a procedure could also be performed
for long wavelength velocity contributions in MOME, but this was not done in
ref. [49] where only long wavelength displacements were resummed in the LPT
manner. We outline this resummation of long wavelength velocities in Appendix
A. While the MOME approach has been shown to give excellent predictions for
the Fourier-space power spectrum, configuration space statistics (where the
BAO focus is not merely few-per cent oscillations on top of the broadband) are
expected to be more sensitive to the details of IR resummation and BAO
damping.

More generally, IR resummation has also been extensively studied in the
Eulerian context [42, 43, 92, 46, 93]. Most often, these rely on the wiggle-
no-wiggle splitting procedure, separating the smooth and BAO components of the
linear power spectrum (see, e.g. [70, 43, 93]). This procedure allows for a
simplified treatment of the nonlinear effects of the BAO where typically only
leading effects are captured, neglecting the more intricate structure captured
by LPT. Nonetheless, the controlled approximations that enter into this form
of IR resummation are generally subdominant to higher-order (two-loop)
corrections that have been studied in e.g. refs. [43, 44, 34]. In addition to
these EPT approaches, refs. [47, 48] take an intermediate approach. These rely
on an LPT-like resummation procedure that tries to preserve the unresummed EPT
broadband behaviour, thus effectively retaining an EPT-like perturbative
structure. In this approach, the anisotropic part of the exponent in Eq. (3.3)
is expanded while only the isotropic part is left resummed (see also Appendix
B of ref. [15] for a more detailed connection between the two approaches).

Our goal in this subsection is to investigate how one-loop LPT with long
wavelength velocity contributions fully resummed compares to the approaches
mentioned above, focusing on the anisotropic redshift-space broadband and the
BAO feature in configuration space. Figure 3 compares the Fourier-space
multipoles predicted by our fiducial LPT framework to the LPT moment expansion
(MOME) and one-loop resummed EPT (REPT). All three frameworks are fitted
assuming Gaussian covariances up to $k_{\rm max}=0.2\,h\,{\rm Mpc}^{-1}$ in
$P(k,\mu)$ as in the previous subsection. All three frameworks show excellent
agreement with the data, with any disagreements, including inter-framework
disagreements, well within the few percent systematic errors we expect from
these simulations. The frameworks differ most in the the hexadecapole, with
the pure LPT framework apparently a better fit to the broadband shape over the
scales shown; this should be taken with a grain of salt, however, as the
statistical errors are large and we specifically checked our IR resummation
procedure for the LPT framework using these data. Similarly, in Figure 4 we
fit the correlation function multipoles of the fiducial halo sample using the
three frameworks above as well as the Gaussian streaming model (GSM). All four
are in excellent agreement with regards to both the BAO feature and broadband
shape at quasi-linear scales. Since the correlation function multipoles probe
a slightly different combination of modes than their Fourier-space
counterparts with a hard $k$ cut, we have adjusted the best-fit bias
parameters “by eye” to yield a better fit at $r>30\ h^{-1}\text{Mpc}$, though
we note that directly transforming the previous Fourier-space results still
yield theoretical predictions within the few-percent systematic errors
expected from these simulations. Together, Figures 3 and 4 suggest that,
despite differences in IR resummations schemes, existing effective-theory
frameworks of the redshift-space two-point function offer similar levels of
performance on pertubative scales. Finally, we note that, while in the above
comparisons we have independently fit the bias parameters and EFT corrections
of each model to most favorably evaluate the performance of each, their bias
bases can in principle be mapped onto each other order-by-order; when the bias
parameters are thus fixed, these models will tend to make slightly different
predictions due to differences in resummed IR modes at higher order. We
discuss these differences in Appendix A.

###  5.3 Cosmological Constraints using Blind Challenge Data

Figure 5: Multipoles of the blind challenge power spectrum along with the best
fit one-loop LPT, MOME and REPT models. The top panels shows the (unbinned)
theory curves along with the data. Both the error bars and theory differences
are too small to see except in a few places. The lower panels show the
fractional residuals of each (binned) theory curve, with each $k$ and $\ell$
bin separated by $0.0008\,h\,{\rm Mpc}^{-1}$ for clarity of presentation.

Figure 6: Histograms and two-dimensional contours for the three cosmological
parameters in the blind challenge. (Left): Comparison of the LPT model in this
paper to our previous submissions using the moment expansion (MOME) and
resummed Eulerian perturbation theory (REPT), all at $k_{\rm
max}=0.12\,h\,{\rm Mpc}^{-1}$. The LPT model performs competitively to
existing models and indeed slightly improves upon the LPT-based MOME model’s
constraints on $h$. (Right): The LPT model constraints using three different
scale cuts. All three scale cuts recover the truth on these parameters to
within $2\,\sigma$.

As a final test of our LPT model, we use it to model the “blind challenge”
data described in ref. [94]555The data and more information about the blind
challenge can be found at
https://www2.yukawa.kyoto-u.ac.jp/~takahiro.nishimichi/data/PTchallenge/..
These are redshift-space power spectra for a BOSS-like HOD sample at $z\simeq
0.6$ constructed from ten N-body boxes each with sidelength
$L=3.84\,h^{-1}$Gpc sampled with $3072^{3}$ equal mass particles. Since the
total volume amounts to about 100 times the volume of the BOSS DR12 sample
[95], the statistical error associated with these data are expected to be far
below any realizable galaxy survey at this redshift. These data were designed
for a blind challenge wherein three cosmological parameters
($\Omega_{M},h,\ln[10^{10}A_{s}]$) need to be fit assuming these tiny
statistical errors while the baryon fraction, $f_{b}$, and spectral tilt,
$n_{s}$, are fixed to the values used in the simulations. The challenge was
designed to evaluate the performance of different PT approaches. Any group
wishing to enter the challenge submits their best fit cosmological parameters,
without knowing the true ones, to Takahiro Nishimichi and collaborators. After
submission one discovers if the model provides an unbiased estimate of the
parameters.

We had previously submitted best fit parameters for MOME and REPT models that
included scales up to $k_{\rm max}=0.12\,\,h\,{\rm Mpc}^{-1}$, and in both
cases the models “passed”: our REPT submission yielded means for the
cosmological parameters well within $1\,\sigma$ of the truth, while MOME
yielded $1.1\,\sigma$ and $1.8\,\sigma$ deviations for $\Omega_{m}$ and $h$,
both well within errors expected for realistic galaxy surveys666Indeed, MOME
also yields errors below $1\,\sigma$ for $k_{\rm max}=0.14\,h\,{\rm
Mpc}^{-1}$, though we did not know this prior to submission and unblinding.
and possibly consistent with fluctuations in the challenge data themselves. We
have repeated the same exercise with the direct LPT model discussed in this
work, using the same set of parameters777These are
$b_{1},b_{2},b_{s},\alpha_{0},\alpha_{2},R_{h}^{3},\sigma_{2}$ in the
Lagrangian basis and the equivalent set mapped onto the Eulerian basis. We
have dropped the Lagrangian third-order bias because it is expected to be
small and somewhat degenerate with other terms, $\alpha_{4,6}$ because we fit
only up to the quadrupole and $\sigma_{4}$ because it was not necessary to fit
the data at the scales we fit.. As we have already participated in the
challenge, we now know the true cosmological parameters. However this should
not affect the evaluation of our new LPT approach, since the analysis pipeline
is the same one we adopted for MOME and REPT and we did not change the model
in any way from that described in previous sections in order to participate in
the challenge except to use the unblinded values as a seed in the MCMC to more
quickly reach the maximum likelihood region. As in our previous submissions,
uninformative priors were placed on all of the model parameters.

Figure 5 shows the measurements of the multipoles of the power spectrum along
with the best fit one-loop LPT, MOME and EPT models. It’s worth noticing that
since we are also fitting for cosmological parameters, compared to the
previous section where the linear power spectrum was held fixed and we varied
only the bias parameters, the model has to include the Alcock-Paczynski (AP)
effect888Specifically, we use Equations 42-46 of ref. [57], though without the
factors of $r_{s}$. [96, 57]. Since this plot is just for visual comparison,
we only show the best fit model with $k_{\rm max}=0.12\,h\,{\rm Mpc}^{-1}$.
These data, produced using a different N-body code, halo finder and HOD
prescription at a different redshift than the simulations in §5.1, act as an
additional test of the three PT models, and indeed the agreement between the
models and with the data is remarkable up to the smallest scales included in
the fit. All three models have $\chi^{2}/\text{dof}\approx 13/(24-10)$,
demonstrating good fits compared to their degrees of freedom.

Turning to the cosmological parameters, the left panel of Figure 6 shows the
$1\,\sigma$ and $2\,\sigma$ constraints obtained by fitting the monopole and
quadrupole of $P(k)$ up to $k_{\rm max}=0.12\,h\,{\rm Mpc}^{-1}$ using our
LPT, MOME and REPT models. This was the scale cut we chose in submitting
results using MOME and REPT to the blind challenge, conservatively selected
given the unusually low statistical uncertainty of the sample, as well as the
main case analyzed in ref. [94]. The LPT model performs slightly better than
the other two, providing unbiased constraints on the three cosmological
parameters. In particular the bias in the Hubble constant, $h$, is reduced in
LPT compared to both MOME and REPT. While the difference is less than
$2\,\sigma$ and therefore well within the realm of possible statistical
fluctuations in the N-body data, the improvement in our $\Omega_{M}$ and $h$
constraints, particularly relative to MOME in which bulk velocities are not
fully resummed, suggests that our IR resummation scheme is correctly capturing
the effects of large scale modes on both the BAO feature and broadband shape.

The right panel in Figure 6 shows the two-dimensional confidence intervals for
different choices of $k_{\rm max}=0.12$, 0.16, $0.20\,h\,{\rm Mpc}^{-1}$ for
LPT. The constraints are within $2\,\sigma$ of the truth for each scale cut,
but at $k_{\rm max}=0.20\,h\,{\rm Mpc}^{-1}$ the Hubble parameter $h$ has a
mean very close to $2\,\sigma$ away from the truth while all parameters are
well within $1\,\sigma$ at $0.12\,h\,{\rm Mpc}^{-1}$, suggesting growing
systematic bias at higher scale cuts where higher-order corrections are
expected to play a more significant role. However, it should be noted that
since these error bars are derived from the covariance of the sample itself —
within the Gaussian approximation no less — the exhibited errors are well
within the realm of statistical possibility and we cannot conclusively
determine that any model is biased. These results are summarized in Figure 7,
which shows the shift between the inferred and true bias parameters as a
function of $k_{\rm max}$. The shaded regions indicate $10\times$ the derived
standard deviations, i.e. approximately the expected errors for a survey like
BOSS, or a redshift slice of $\Delta z=0.25$ at $z=0.6$ for DESI. The LPT
model presented in this work correctly recovers the underlying cosmology for
all the scale cuts shown to well within the expected errors of surveys like
BOSS and DESI. We anticipate that the model would perform even better at
higher redshift where the degree of non-linearity is smaller. We therefore
conclude that, at least for the cosmological parameters probed in the
challenge, our LPT model should provide an accurate tool for modeling RSD in
upcoming surveys.

Figure 7: Parameter constraints using the LPT model for the blind challenge as
a function of scale cut $k_{\rm max}$. The shaded blue region show errors
scaled to a survey of $10\times$ less volume, i.e. BOSS or a $\Delta z=0.1$
slice of DESI at $z=1.2$. All constraints shown are within $2\,\sigma$ of the
truth, and within $1\,\sigma$ when $k_{\rm max}\leq 0.14\,h\,{\rm Mpc}^{-1}$,
well within error bars expected for future surveys as well as the realm of
possible statistical fluctuations for this simulated sample.

##  6 Conclusions

The anisotropic galaxy clustering observed by spectroscopic surveys probes
density and velocity fields on large scales, enabling us to test the growth of
structure in the quasilinear regime as predicted by General Relativity. In
addition, the baryon acoustic oscillations in the galaxy clustering signal
provide geometric information that constrain the cosmological expansion
history. Perturbation theory is an ab initio approach with clear physical
assumptions, and it is therefore the preferred tool for a precise and rigorous
mapping between cosmological parameters and the observed clustering signal, in
the quest for a better understanding of the cosmological model and in the
search for new physics.

The main purpose of this work was to further develop the modeling of the
redshift-space two-point function within Lagrangian perturbation theory. Two
critical aspects of any perturbation-theory model of the power spectrum or
correlation function are the treatment of the nonlinear damping of the BAO
signal due to large scale motions of galaxies and the inclusion of redshift-
space distortions due to a degeneracy between the observed line-of-sight
positions and peculiar velocities of galaxies. By expanding directly in
displacements, Lagrangian perturbation theory naturally treats both phenomena
within the same framework without relying on the wiggle-no-wiggle splitting
procedure to damp the BAO oscillations. In particular, LPT exponentiates
linear, or Zeldovich, displacements via the cumulant theorem, thereby
resumming the long-wavelength (IR) modes primarily responsible for nonlinear
BAO damping.

In this work, we extend the numerical techniques developed in refs. [23, 49]
for the Zeldovich approximation to calculate the one-loop LPT power spectrum
of biased tracers with both bulk displacements and velocities resummed. This
is in contrast to previous efforts to model redshift-space distortions within
effective LPT that focused on predicting velocity statistics for which only
the displacements were exponentiated999Earlier work, e.g. ref. [19], performed
a similar calculation as the one here for the correlation function only, but
exponentiated the modes coming from the one-loop contributions as well.. In
Sections 2 and 3, we outlined the effective formalism developed in those works
and used them to write down the redshift-space power spectrum including third-
order biasing, counterterms and stochastic contributions, of which the latter
two in part play the role of “finger-of-god” terms in traditional RSD models.
Then, in Section 4, we developed the required numerical framework for the
implementation of the model. We have publically released our implementation of
Method II as part of the velocileptors
code.101010https://github.com/sfschen/velocileptors Additional details of
these calculations can be found in Appendices B and C.

Finally, we compare the LPT model developed in this work to N-body data and
previous models in Section 5. First, we fit the power spectrum of a halo
sample drawn from a set of N-body simulations at $z=0.8$ with masses
$12.5<\log(M/M_{\odot})<13.0$ assuming Gaussian covariances using our model in
Figure 2, finding excellent agreement in both the wedges and first three
multipoles for a consistent set of bias and effective parameters. Our model
has an extra degree of freedom in the IR cutoff $k_{\rm IR}$, which dictates
the wavelength beyond which displacement modes are resummed. In ref. [49] it
was shown that the higher-order velocity statistics that enter into redshift-
space distortions have broadband shapes that are especially sensitive to this
choice, and indeed we find that a choice of $k_{\rm IR}=0.2\,h\,{\rm
Mpc}^{-1}$, which lies between the broadband predictions of full-expanded LPT
and EPT, gives the best match to the hexadecapole; this choice is explored
further in Appendix A, where we also show (Fig. 10) that the configuration-
space anisotropic BAO feature is remarkably robust to this hyperparameter.

We then compare the performance of our LPT model with other existing
effective-theory models. In Figure 3, we compare the aforementioned power
spectrum predictions to the Lagrangian moment expansion (MOME) and resummed
Eulerian perturbation theory, finding that all three can fit the data to
within the expected systematic error of the simulations. Similar results for
the configuration-space multipoles are shown in Figure 4, where we also
compare to the Lagrangian Gaussian streaming model. As a last numerical test
we checked whether the full one-loop LPT model presented in this work can
recover unbiased cosmological parameters in a data analysis challenge using a
different N-body code, halo finder, halo occupation distribution and redshift.
We find that LPT performs better than MOME and REPT at the reference scale cut
chosen in ref. [94], and is able to recover the true cosmology with negligible
biases up to $k_{\rm max}=0.2\,h\,{\rm Mpc}^{-1}$ in a volume approximately
one hundred times that of DESI or Euclid at the same redshift. Results of this
test are given in Figure 6.

Let us conclude by noting some possible extensions of our model and numerical
implementation. There has been considerable recent interest in extending the
bias expansion of galaxies to include anisotropic selection effects and,
indeed, the redshift-space galaxy density can be decomposed into a generalized
expansion of operators with LOS symmetry [97, 33, 98, 99]. Since the angular
dependencies of these operators will in general involve only tensor products
of $\hat{k}$, $\hat{n}$ and $\hat{q}$ their 2-point functions (with IR
displacements resummed) should follow straightforwardly from our calculations.
Our calculations should also be straightforwardly extendable to modeling the
reconstructed galaxy power spectrum at one-loop order, especially the “Rec-
Sym” scheme [49] which features an identical structure, with the only
difference being a larger set of terms involved. Concerning extensions of the
$\Lambda$CDM model, perhaps the simplest one to implement is massive neutrino
cosmologies. It is well known that halos and galaxies are biased tracers of
the dark matter and baryon fluids only [100, 101, 102, 103, 104, 105, 106],
which implies the bias expansion presented in Section 2 will still be valid
with the trivial replacement of the total matter field with the dark matter
plus baryon one. Extra care should be taken with RSD since the growth rate is
now scale dependent and Equation 2.3 is no longer valid [107] (and similarly
in modified gravity theories [108, 109]). However given the smallness of
neutrino masses, this complication is usually neglected when evaluating loop
integrals. We intend to return to these, admittedly more involved,
calculations in future work.

## Acknowledgments

We thank Jahmour Givans, Marko Simonović and the anonymous referee for helping
us clarify various points in the draft. S.C. is supported by the National
Science Foundation Graduate Research Fellowship (Grant No. DGE 1106400) and by
the UC Berkeley Theoretical Astrophysics Center Astronomy and Astrophysics
Graduate Fellowship. M.W. is supported by the U.S. Department of Energy and
the NSF. Z.V. is supported by the Kavli Foundation. This research has made use
of NASA’s Astrophysics Data System and the arXiv preprint server. This
research used resources of the National Energy Research Scientific Computing
Center (NERSC), a U.S. Department of Energy Office of Science User Facility
operated under Contract No. DE-AC02-05CH11231.

##  Appendix A Infrared Resummation and the Broadband

Figure 8: Power spectrum predictions of the three models (LPT, MOME, REPT)
given the same set of bias parameters, but with all counterterms and
stochastic parameters adjusted in the MOME and REPT curves to match the LPT
prediction with no counterterms or stochastic terms. While much of the
differences between the three formalisms can be soaked up by the counterterms,
the fractional residuals (lower panel) at high $k$ and $\mu$ begin to show
less trivial behavior, suggesting non-negligible theory error.

Our goal in the main body of this paper was to develop the fully-resummed LPT
model and compare its performance to existing models such MOME and REPT. To be
as fair as possible to each of these models, we have independently fit for the
bias parameters in counterterms in each. However, it should be noted that in
principle the bias bases for LPT and MOME are identical, and equivalent up to
a mapping to the Eulerian basis employed in REPT111111See, for example,
Equation 4.8 in ref. [49]. However, while all three models should be
equivalent order-by-order under these mappings, they tend to make somewhat
different predictions, especially towards small scales, due to the different
IR resummation schemes involved. Namely, each of the three schemes organizes
the perturbative expansion in slightly different expansion parameters. While
LPT resums all the two-point long wavelength displacement and velocity
contributions, in EPT these are considered perturbative and are accordingly
expanded. Nevertheless, in EPT what is resummed are the contributions to the
BAO feature from the relative motions of the long modes. In the MOME
expansion, RSD contributions are organized following the moment expansion [37,
38, 39, 23], while each of the contributions is then evaluated in LPT [49].
This seemingly leaves the long velocity contributions expanded, contrary to
the full LPT approach. In ref. [49] these were left un-resummed, although a
straight forward approach to add these would follow the EPT procedure, just
excluding the displacement contributions which have already been resummed. For
MOME we can thus write

$\displaystyle P^{s,IR}_{\rm 1-loop}(\textbf{k})$ $\displaystyle\approx
P^{s,{\rm nw}}_{\rm
1-loop}(\textbf{k})+e^{-\frac{1}{2}\Sigma_{s}^{2}(\mu)k^{2}}\left(1+\tfrac{1}{2}\Sigma_{s}^{2}(\mu)k^{2}\right)P^{s,{\rm
w}}_{\rm
lin}(\textbf{k})+e^{-\frac{1}{2}\Sigma_{s}^{2}(\mu)k^{2}}\left(P^{s}_{\rm
loop}(\textbf{k})-P^{s,{\rm nw}}_{\rm loop}(\textbf{k})\right),$

where the wiggle and no-wiggle $P^{s}_{\rm 1-loop}$ (and similarly the
$P^{s}_{\rm loop}$ by dropping the linear Kaiser part) are given by MOME
predictions computed in ref. [49]. The difference with the EPT resummation
scheme is in the definition of $\Sigma_{\rm s}$ which now contains only the
velocity contributions

$\Sigma_{s}^{2}(\mu)=f(f+2)\mu^{2}\Sigma^{2},$ (A.1)

where $\Sigma^{2}$ is the velocity dispersion due to the long wavelength
modes.

Figure 8 shows the predictions for $P(k,\mu)$ of LPT, MOME and REPT when the
bias parameters ($b_{1}$, $b_{2}$, $b_{s}$, $b_{3}$) are fixed to the best-fit
values for the fiducial halo sample in LPT, with counterterms and stochastic
contributions in the MOME and REPT cases adjusted to fit the LPT result. The
three schemes differ systematically towards higher $k$ and $\mu$. Compared to
its Eulerian counterpart, LPT shows significant supression of power towards
high $k$, particularly along the line of sight; this suppression is absorbed
by adjusting counterterms and stochastic contributions in Figure 8, though the
the theories nonetheless begin to diverge at the percent level around
$k=0.15\,h\,{\rm Mpc}^{-1}$, especially towards higher $\mu$. The LPT and MOME
schemes are quite similar at low $\mu$ since they differ only in the inclusion
of higher-order velocities along the line of sight in the former121212In fact,
MOME also shows slightly more suppression at low $\mu$ since we have followed
the main text of ref. [49] and the velocileptors code in not including any IR
cutoffs in the expanded velocity moments.—and both are damped relative to
REPT—but closer to the line of sight their oscillatory components begin to
differ significantly more than between LPT and REPT since long velocity modes
are not resummed in MOME.

The above differences in the three schemes have their origin in the fact that
the expansion parameters in these schemes do not match exactly. As we
mentioned before, in LPT the long wavelength displacement and velocity two-
point contributions are directly resummed, while in REPT only the relative
effects of these are resummed and thus affect only the BAO. The MOME scheme,
on the other hand, takes a hybrid approach between the previous two. Since the
predictions of the three schemes, as plotted in Figure 8, are equal up to
second (one-loop) order in the linear power spectrum when expanded order-by-
order; the apparent differences therefore reflect differences at two-loop
order or beyond, even though these residuals are due to the long mode
contributions and are nominally under the perturbative control. At low $k$
they manifest as contributions proportional to the wavenumber squared and can
be largely absorbed by existing counterterms $\propto k^{2}\mu^{2n}P(k)$ and
stochastic contributions. The residual deviations at higher $k$ and $\mu$, can
also provide rough estimates of the theory error of these common perturbative
schemes, indicating the range of validity of current perturbative models. We
intend to return to a more in-depth comparison of these schemes in a future
work.

Figure 9: Power spectrum multipole predictions for a range of IR cutoffs
$k_{\rm IR}$. Including more IR modes by upping $k_{\rm IR}$ damps the total
power at high $k$, especially in the higher multipoles. The limit of $k_{\rm
IR}\rightarrow 0$ corresponds to unresummed EPT (thin solid lines), with
noticeably larger BAO oscillations at both linear and one-loop order, though
even $k_{\rm IR}=0.05\,h\,{\rm Mpc}^{-1}$ dramatically reduces these.

In addition, as noted in the main body of the text, the choice of infrared
cutoff $k_{\rm IR}$ has a significant effect on the broadband power of the
anisotropic power spectrum within the LPT model itself. Figure 9 shows this
effect on the power spectrum multipoles, with bias parameters fixed to those
obtained by fitting the fiducial sample to LPT with $k_{\rm IR}=0.2\,h\,{\rm
Mpc}^{-1}$. In all three multipoles shown, decreasing the IR cutoff results in
increasing power at high $k$, with higher $\ell$ increasingly sensitive to
these nonlinearities. Indeed, by $k=0.1\,h\,{\rm Mpc}^{-1}$ the hexadecapoles
with $k_{\rm IR}$ equal to zero is close to twice as large as that with no IR
cutoff, i.e. $k_{\rm IR}=\infty$. In the $k_{\rm IR}\rightarrow 0$ limit the
LPT prediction is equal to that in EPT without any additional IR resummation
and with the bias parameters appropriately mapped. Of course, we caution that
while the differences shown in Figure 9 are intended to demonstrate the full
spectrum of resummations possible in our scheme, in reality much of these
differences can be absorbed by effective corrections as in Figure 8.

The increasing effect of the IR cutoff on redshift space distortions at high
$\mu$ can be understood intuitively within the language of the moment
expansion. In general, higher-order velocities will tend to be more sensitive
to the IR cutoff. This can be seen as follows: the velocity statistics of
interest for RSD can be schematically written as $M_{n}=\left\langle
X^{n}e^{iX}\right\rangle$ where $X=\textbf{k}\cdot\Delta$. Approximating $X$
to be Gaussian with variance $\sigma^{2}$ we can write the even moments as

$M_{2n}=\langle X^{2n}\rangle(1+a_{2}\sigma^{2}+a_{4}\sigma^{4}+...\ )\
\exp\left[-\sigma^{2}/2\right].$ (A.2)

For $n>0$ we always have $a_{2}<0$, i.e. $M_{2n}$ damps faster than the
exponential damping in $M_{0}$. This is easily understood: higher $M_{n}$
receive more contributions from larger values of $X$, where the complex
exponential oscillates rapidly, and are suppressed by $X^{n}$ at small $X$
where the exponential varies slowly. Indeed, this effect was observed in ref.
[49], where it was noted that the broadband of the second moment of the
pairwise velocity, $\sigma_{\rm
12,ij}=\left\langle(1+\delta_{1})(1+\delta+2)\dot{\Delta}_{i}\dot{\Delta}_{j}\right\rangle$
is very sensitive to cutoff choice, with its monopole and quadrupole
respectively being better predicted by higher and lower values of $k_{\rm
IR}$. Since the second moment’s quadrupole is the leading $\mu^{4}$
contribution to the anisotropic power spectrum, one might expect that $P_{4}$
should in turn be very sensitive to this choice. In light of the effects of
$k_{\rm IR}$ on the second moment, in this paper we have chosen the
“intermediate” value of $0.2\,h\,{\rm Mpc}^{-1}$ as our fiducial IR cutoff (at
$z\approx 0.8$), though we caution that further investigation is warranted
when operating at significantly higher or lower redshifts or with highly
biased tracers.

Finally, let us note that, in contrast to the anisotropic broadband, the
corresponding BAO features in the correlation function monopole and
quadrupole, shown in Figure 10, are rather insensitive to the specific choice
of $k_{\rm IR}$. Indeed, even $k_{\rm IR}=0.05\,h\,{\rm Mpc}^{-1}$, which is
almost identical to EPT in its broadband, demonstrates significant damping of
the BAO feature. This figure also shows the unresummed EPT ($k_{\rm IR}=0$)
prediction, which clearly illustrates the non-convergence of the
configuration-space BAO feature in one-loop EPT that necessitates the a
posteriori IR resummation implemented in these models.

Figure 10: Correlation function multipole predictions with a range of IR
cutoffs $k_{\rm IR}$ for bias parameters fixed to those obtained from the
fiducial halo sample. The BAO features in both the monopole and quadrupole are
rather robust, showing little change for $k_{\rm IR}>0.05\,h\,{\rm Mpc}^{-1}$,
despite significant differences in the corresponding power spectrum
broadbands. For reference, the unresummed EPT ($k_{\rm IR}=0$) prediction,
which clearly illustrates the non-convergence of the nonlinear configuration-
space BAO feature absent IR resummation, is also shown (thin solid lines).

##  Appendix B Method I

In this section we extend Method I, first presented in [23], to include terms
up to one-loop order. In contrast to Method II, described in the main body of
the text, this method does not rely on boosting the wavevector k into a more
convenient frame but rather evaluates the integral in Equation 3.1 directly in
the observed frame.

From the form of the RSD operator
$R^{(n)}_{ij}=\delta_{ij}+nf\hat{n}_{i}\hat{n}_{j},$ where $n$ index counts a
given order in PT, we can see that all integrals in Eq.(3.1) can be written in
terms of scalar functions and dot products between three unit vectors
($\hat{q}$, $\hat{k}$ and $\hat{n}$). The angular structure is given in terms
of products

$\hat{n}\cdot\hat{k}=\mu,\quad\hat{q}\cdot\hat{k}=\mu_{\textbf{q}},\quad\hat{q}\cdot\hat{n}=\mu_{\textbf{q}}\mu+\sqrt{1-\mu_{\textbf{q}}^{2}}\,\sqrt{1-\mu^{2}}\,\cos\phi,$
(B.1)

where $\phi$ is the azimuthal angle in a polar coordinate system where the
zenith is given by $\hat{k}$ and the plane $\phi=0$ is spanned by $\hat{k}$
and $\hat{q}$. The effect of RSD operators $\mathbf{R}^{(n)}$ can then be
captured by looking how it acts on the tensor basis comprised of
$\hat{q}_{i},\delta_{ij},$. We have

$\displaystyle k_{i}k_{j}\delta_{ij}$ $\displaystyle\rightarrow
k_{i}k_{j}R^{(m_{1})}_{in}R^{(m_{2})}_{jn}=k^{2}\left[1+f(m_{1}+m_{2}+m_{1}m_{2}f)\mu^{2}\right],$
(B.2) $\displaystyle k_{i}\hat{q}_{i}$ $\displaystyle\rightarrow
k_{i}R^{(m)}_{ij}\hat{q}_{j}=k\mu_{\textbf{q}}\left[1+mf\mu^{2}+mf\mu^{2}\gamma(\mu_{\textbf{q}},\mu)\cos\phi\right],$
(B.3)

where
$\gamma(\mu_{\textbf{q}},\mu)=\sqrt{1-\mu_{\textbf{q}}^{2}}\,\sqrt{1-\mu^{2}}/\mu_{\textbf{q}}\mu.$
The azimuthal dependence of the exponentiated linear displacements
$A^{s}_{ij}$ requires us to calculate polar-coordinate integrals of the form
[23]

$I_{n}(f,\mu_{\textbf{q}},\mu)=\int^{2\pi}_{0}\frac{d\phi}{2\pi}\
e^{-\frac{1}{2}k^{2}Y(\alpha_{2}\gamma\cos\phi+\alpha_{3}\gamma^{2}\cos^{2}\phi)\mu_{\textbf{q}}^{2}},$
(B.4)

where $\alpha_{2}=f\mu^{2}(m_{1}+m_{2}+fm_{1}m_{2}\mu^{2})$ and
$\alpha_{3}=m_{1}m_{2}f^{2}\mu^{4}$. This expression can be used to compute
all the other loop contributions (except the $W_{ijk}$ term). These can be
calculated by taking derivatives in either $\alpha$ or $\beta$ of the identity

$I_{\phi}\left(\alpha,\beta,\mu_{\textbf{q}}\right)=\int_{0}^{2\pi}\frac{d\phi}{2\pi}~{}e^{\alpha\mu_{\textbf{q}}\sqrt{1-\mu_{\textbf{q}}^{2}}\cos\phi+\beta(1-\mu_{\textbf{q}}^{2})\cos^{2}\phi}=\sum_{\ell=0}^{\infty}F_{\ell}(\alpha,\beta)\left(\alpha^{2}\mu_{\textbf{q}}^{2}/\beta\right)^{\ell},~{}$
(B.5)

where

$\displaystyle F_{\ell}(\alpha,\beta)$
$\displaystyle=\sum_{m=0}^{\ell}\frac{\Gamma(m+\frac{1}{2})}{\pi^{1/2}\Gamma(m+1)\Gamma(1+2m-\ell)\Gamma(2\ell-2m+1)}\left(-\frac{\beta^{2}}{\alpha^{2}}\right)^{m}$
$\displaystyle\hskip 68.28644pt\times
M\left(\ell-2m;\ell-m+\tfrac{1}{2};\frac{\alpha^{2}}{4\beta}\right)M\left(m+\frac{1}{2};m+1;\beta\right)$
(B.6)

and $M(a,b,z)$ are hypergeometric functions of the first kind.

In order to capture the contribution of $W_{ijk}$ a slight generalisation of
the integrals above is required. The integral we need is of the form

$\int_{0}^{2\pi}\frac{d\phi}{2\pi}~{}\left(\sqrt{1-\mu_{\textbf{q}}^{2}}\cos\phi\right)e^{\alpha\mu_{\textbf{q}}\sqrt{1-\mu_{\textbf{q}}^{2}}\cos\phi+\beta(1-\mu_{\textbf{q}}^{2})\cos^{2}\phi}=\frac{1}{\mu_{\textbf{q}}}\frac{d}{d\alpha}I_{\phi}\left(\alpha,\beta,\mu_{\textbf{q}}\right).$
(B.7)

However, since $F_{0}(\alpha,\beta)$ does not depend on $\alpha$ we have

$\frac{1}{\mu_{\textbf{q}}}\frac{d}{d\alpha}I_{\phi}\left(\alpha,\beta,\mu_{\textbf{q}}\right)=\mu_{\textbf{q}}\frac{1}{\beta}\frac{d}{\
d\alpha}\sum_{\ell=0}^{\infty}\alpha^{2}F_{\ell+1}(\alpha,\beta)\left(\alpha^{2}\mu_{\textbf{q}}^{2}/\beta\right)^{\ell}.$
(B.8)

##  Appendix C Method II

###  C.1 General Mathematical Structure

The workhorse integral of Method II is [23]

$I(A,B,C)=\int d\mu_{\textbf{q}}\ d\phi\
e^{-iC\sqrt{1-\mu_{\textbf{q}}^{2}}\cos\phi+iA\mu_{\textbf{q}}+B\mu_{\textbf{q}}^{2}}=4\pi
e^{B}\sum_{\ell=0}^{\infty}\Big{(}\frac{-2}{\rho}\Big{)}^{\ell}\tilde{G}^{(0)}_{0,\ell}(A,B,\rho)j_{\ell}(\rho),$

where $\rho^{2}=A^{2}+C^{2}.$ The exact form of the kernel
$\tilde{G}^{(0)}_{0,\ell}$ is given in Section C.4.

To compute the redshift-space power spectrum for biased tracers we will need
the integrals derived from $A,C$ derivatives of $I$, i.e.

$\displaystyle I_{n,m}(A,B,C)$ $\displaystyle=i^{m-n}\int d\mu_{\textbf{q}}\
d\phi\
e^{-iC\sqrt{1-\mu_{\textbf{q}}^{2}}\cos\phi+iA\mu_{\textbf{q}}+B\mu_{\textbf{q}}^{2}}\big{(}\sqrt{1-\mu_{\textbf{q}}^{2}}\cos\phi\big{)}^{n}\mu_{\textbf{q}}^{m}$
(C.1) $\displaystyle=4\pi
e^{B}\sum_{\ell=0}^{\infty}\Big{(}\frac{-2}{\rho}\Big{)}^{\ell}\tilde{G}^{(m)}_{n,\ell}(A,B,\rho)j_{\ell}(\rho).$
(C.2)

These satisfy the recursion relations

$\tilde{G}^{(m)}_{n,\ell}=\frac{d\tilde{G}^{(m)}_{n,\ell}}{dA}+\frac{A}{2}\
\tilde{G}^{(m-1)}_{n,\ell-1},\quad\tilde{G}^{(m)}_{n,l}=\frac{\partial\tilde{G}^{(m)}_{n,\ell}}{\partial
C}+\frac{C}{2}\tilde{G}^{(m-1)}_{n,\ell-1}.$ (C.3)

For convenience we list the first few $\tilde{G}^{(m)}_{0,\ell}$:

$\displaystyle\tilde{G}^{(1)}_{0,\ell}$
$\displaystyle=\frac{\partial\tilde{G}^{(0)}_{0,\ell}}{\partial
A}+\frac{A}{2}\tilde{G}^{(0)}_{0,\ell-1}$
$\displaystyle\tilde{G}^{(2)}_{0,\ell}$
$\displaystyle=\frac{\partial^{2}\tilde{G}^{(0)}_{0,\ell}}{\partial
A^{2}}+A\frac{\partial\tilde{G}^{(0)}_{0,\ell-1}}{\partial
A}+\frac{1}{2}\tilde{G}^{(0)}_{0,\ell-1}+\frac{A^{2}}{4}\tilde{G}^{(0)}_{0,\ell-2}$
$\displaystyle\tilde{G}^{(3)}_{0,\ell}$
$\displaystyle=\frac{\partial^{3}\tilde{G}^{(0)}_{0,\ell}}{\partial
A^{3}}+\frac{3A}{2}\frac{\partial^{2}\tilde{G}^{(0)}_{0,\ell-1}}{\partial
A^{2}}+\frac{3}{2}\frac{\partial\tilde{G}^{(0)}_{0,\ell-1}}{\partial
A}+\frac{3A^{2}}{4}\frac{\partial\tilde{G}^{(0)}_{0,\ell-2}}{\partial
A}+\frac{3A}{4}\tilde{G}^{(0)}_{0,\ell-2}+\frac{A^{3}}{8}\tilde{G}^{(0)}_{0,\ell-3}.$
(C.4)

The derivatives with respect to $C$ are entirely analogous, swapping $C$ for
$A$ and $m$ for $n$. In addition, we will need

$\tilde{G}^{(1)}_{1,\ell}=\frac{\partial^{2}\tilde{G}^{(0)}_{0,\ell}}{\partial
A\partial C}+\frac{C}{2}\frac{\partial\tilde{G}^{(0)}_{0,\ell-1}}{\partial
A}+\frac{A}{2}\frac{\partial\tilde{G}^{(0)}_{0,\ell-1}}{\partial
C}+\frac{AC}{4}\tilde{G}^{(0)}_{0,\ell-2}$
$\displaystyle\tilde{G}^{(2)}_{1,\ell}$
$\displaystyle=\frac{\partial^{3}\tilde{G}^{(0)}_{0,\ell}}{\partial
A^{2}\partial
C}+\frac{C}{2}\frac{\partial^{2}\tilde{G}^{(0)}_{0,\ell-1}}{\partial
A^{2}}+A\frac{\partial^{2}\tilde{G}^{(0)}_{0,\ell-1}}{\partial A\partial
C}+\frac{1}{2}\frac{\partial\tilde{G}^{(0)}_{0,\ell-1}}{\partial C}$
$\displaystyle+\frac{AC}{2}\frac{\partial\tilde{G}^{(0)}_{0,\ell-2}}{\partial
A}+\frac{A^{2}}{4}\frac{\partial\tilde{G}^{(0)}_{0,\ell-2}}{\partial
C}+\frac{C}{4}\tilde{G}^{(0)}_{0,\ell-2}+\frac{A^{2}C}{8}\tilde{G}^{(0)}_{0,\ell-3}$

###  C.2 Example: One-Loop Matter Power Spectrum in Redshift Space

Let us consider the one-loop matter power spectrum as an example for the kinds
of angular terms that can arise. For simplicity, we focus on what happens to
the un-exponentiated one-loop contribution $A^{(22)}_{ij}$. In this case we
have that the relevant quantity in redshift space is

$\displaystyle k_{i}k_{j}(\mathbf{R}^{(1)}+f\
\hat{n}\otimes\hat{n})_{in}(\mathbf{R}^{(1)}+f\
\hat{n}\otimes\hat{n})_{jm}A^{(22)}_{nm}=\big{(}K_{i}K_{j}+2fk_{\parallel,i}K_{j}+f^{2}k_{\parallel,i}k_{\parallel,j}\big{)}A^{(22)}_{ij}.$
(C.5)

The first piece ($K_{i}K_{j}$) is identical in angular structure to those in
the Zeldovich case, so we restrict our attention to the other two.

Let’s begin with the term proportional to $f$ in Equation C.5. We have
(dropping the $(22)$ for brevity)

$\displaystyle K_{i}k_{\parallel,j}A_{ij}$ $\displaystyle=(K\cdot
k_{\parallel})X+(\hat{q}\cdot K)(\hat{q}\cdot k_{\parallel})Y$
$\displaystyle=k^{2}\mu^{2}(1+f)\ X+(K\mu_{\textbf{q}})\
k\mu(A(\mu)\mu_{\textbf{q}}+B(\mu)\sqrt{1-\mu_{\textbf{q}}^{2}}\cos\phi)\ Y.$

The piece proportional to $X$ poses no problem since it has no angular
dependence. The term proportional to $Y$ has a piece proportional to
$\mu_{\textbf{q}}^{2}$ , which can be computed via two derivatives of Equation
4.6 w.r.t. $A$, and another with $\phi$ dependence calculable via a $C$
derivative; both are of the form C.2. The term proportional to $f^{2}$ is
similar and involves

$k_{\parallel,i}k_{\parallel,j}A_{ij}=k^{2}\mu^{2}\left[X+(\hat{q}\cdot\hat{n})^{2}Y\right].$

This piece proportional to $Y$ then involves up to two $C$ derivatives.

###  C.3 General Angular Structure of Bias Contributions

Let us now list all the possible angular dependencies at one-loop order,
organized in powers $\mu_{\textbf{q}}^{a}(\hat{n}\cdot\hat{q})^{b}$. This
format is convenient because each such power can be readily integrated in
$\phi$ and $\mu_{\textbf{q}}$ to give

$\displaystyle\mu_{\textbf{q}}^{a}(\hat{n}\cdot\hat{q})^{b}=\sum_{n=0}^{b}\binom{b}{n}$
$\displaystyle
A^{n}(\mu)B^{b-n}(\mu)\mu_{\textbf{q}}^{a+n}\big{(}\sqrt{1-\mu_{\textbf{q}}^{2}}\cos\phi\big{)}^{b-n}\rightarrow$
$\displaystyle
4\pi\sum_{n=0}^{b}\binom{b}{n}A^{n}(\mu)B^{b-n}(\mu)\Big{(}\frac{-2}{kq}\Big{)}^{\ell}\tilde{G}_{b-n,\ell}^{a+n}\
j_{\ell}(kq)$

The simplest case involves correlators with one order $n$ displacement, of the
form $U^{(n)}_{i}=U(q)\hat{q}_{i}$:

$\displaystyle
k_{i}U^{s,(n)}_{i}=[K\mu_{\textbf{q}}+f(n-1)k\mu(\hat{q}\cdot\hat{n})]\ U(q)$
(C.6)

Then we have terms involving two displacements with order $n,m$, which we can
write as $A^{(n,m)}_{ij}=X\delta_{ij}+Y\hat{q}_{i}\hat{q}_{j}$:

$\displaystyle k_{i}k_{j}A^{s,(n,m)}_{ij}=$ $\displaystyle
K^{2}[X(q)+Y(q)\mu_{\textbf{q}}^{2}]+(n+m-2)fk\mu[X(q)(\hat{K}\cdot\hat{n})+Y(q)\mu_{\textbf{q}}(\hat{q}\cdot\hat{n})]$
$\displaystyle+(n-1)(m-1)f^{2}k^{2}\mu^{2}[X(q)+Y(q)(\hat{q}\cdot\hat{n})^{2}].$
(C.7)

Finally, at one-loop order there is one term involving three displacements
involving their (112) bispectrum
$W^{(112)}_{ijk}=V_{1}(\hat{q}_{i}\delta_{jk}+\hat{q}_{j}\delta_{ik})+V_{3}\hat{q}_{k}\delta_{ij}+T\hat{q}_{i}\hat{q}_{j}\hat{q}_{k}$:

$\displaystyle k_{i}k_{j}k_{k}W^{(112)}_{ijk}$
$\displaystyle=2\big{(}K^{3}+fK(k\mu)^{2}(1+f)\big{)}\mu_{\textbf{q}}\
V_{1}(q)$
$\displaystyle+K^{2}\big{(}K\mu_{\textbf{q}}+f(k\mu)(\hat{n}\cdot\hat{q})\big{)}\
V_{3}(q)$
$\displaystyle+K^{2}\big{(}K\mu_{\textbf{q}}^{3}+f(k\mu)\mu_{\textbf{q}}^{2}(\hat{n}\cdot\hat{q})\big{)}\
T(q).$ (C.8)

The full list of non-Zeldovich angular dependences required for the one-loop
power spectrum is given in Table 1.

bias |  $(n,m,l)$: correlator
---|---
$1$ |  $(1,3)$: $A^{(13)}_{ij}$, $(2,2)$: $A^{(22)}_{ij}$, $(1,1,2)$: $W^{(112)}_{ijk}$
$b_{1}$ |  $(3)$: $U^{(3)}_{i}$, $(1,2)$: $A^{10}_{ij}$
$b_{1}^{2}$ |  $(2)$: $U^{11}_{i}$
$b_{2}$ |  $(2)$: $U^{20}_{i}$
$b_{s}$ |  $(2)$: $V^{10}_{i}$
Table 1: Contributions to the one-loop power spectrum and the perturbative
order of the displacements they contain.

###  C.4 $\tilde{G}^{(0)}_{0,m}$ and Its Derivatives

The basic kernel for Method II is the function

$\tilde{G}^{(0)}_{m}(A,B,\rho)=\sum_{n=m}^{\infty}f_{nm}\Big{(}\frac{BA^{2}}{\rho^{2}}\Big{)}^{n}{}_{2}F_{1}\Big{(}\frac{1}{2}-n,-n;\frac{1}{2}-m-n;\frac{\rho^{2}}{A^{2}}\Big{)},$
(C.9)

where $\rho=\sqrt{A^{2}+C^{2}}$, ${}_{2}F_{1}$ is the ordinary hypergeometric
function and $f_{nm}$ is

$f_{nm}=\frac{\Gamma(m+n+\frac{1}{2})}{\Gamma(m+1)\Gamma(n+\frac{1}{2})\Gamma(1-m+n)}\quad.$
(C.10)

The angular dependences in Method II require us to take $A$ and $C$
derivatives of the above. The first three derivatives of $\tilde{G}^{(0)}_{0}$
with respect to $A$ are given by

$\displaystyle\frac{d\tilde{G}^{(0)}_{0,m}}{dA}=\sum_{n=m}^{\infty}\Big{(}\frac{BA^{2}}{\rho^{2}}\Big{)}^{n}$
$\displaystyle
f_{nm}\Bigg{[}\Big{(}\frac{2n}{A}-\frac{2nA}{\rho^{2}}\Big{)}{}_{2}F_{1}\Big{(}\frac{1}{2}-n,-n;\frac{1}{2}-m-n;\frac{\rho^{2}}{A^{2}}\Big{)}$
$\displaystyle+$
$\displaystyle\Big{(}-\frac{2\rho^{2}}{A^{3}}+\frac{2}{A}\Big{)}\frac{(\frac{1}{2}-n)(-n)}{(\frac{1}{2}-m-n)}{}_{2}F_{1}\Big{(}\frac{3}{2}-n,1-n;\frac{3}{2}-m-n;\frac{\rho^{2}}{A^{2}}\Big{)}\Bigg{]}$
(C.11)
$\displaystyle\frac{d^{2}\tilde{G}^{(0)}_{0,m}}{dA^{2}}=\sum_{n=m}^{\infty}\Big{(}\frac{BA^{2}}{\rho^{2}}\Big{)}^{n}$
$\displaystyle
f_{nm}\Big{(}\frac{\rho^{2}-A^{2}}{\rho^{4}}\Big{)}\Bigg{[}(2m-1-4n(m+1))\
{}_{2}F_{1}\Big{(}\frac{1}{2}-n,-n;\frac{1}{2}-m-n;\frac{\rho^{2}}{A^{2}}\Big{)}$
$\displaystyle+(1-4n^{2}+m(4n-2))\
{}_{2}F_{1}\Big{(}\frac{3}{2}-n,-n;\frac{1}{2}-m-n;\frac{\rho^{2}}{A^{2}}\Big{)}\Bigg{]}.$
(C.12)
$\displaystyle\frac{d^{3}\tilde{G}^{(0)}_{0,m}}{dA^{3}}=\frac{C^{2}}{A\rho^{6}}\sum_{n=m}^{\infty}\Big{(}\frac{BA^{2}}{\rho^{2}}\Big{)}^{n}$
$\displaystyle
f_{nm}\Bigg{[}\Big{(}(2(1-m)(1-2m)+8n(2-m)(1+m)+8n^{2}(1+m))A^{2}$
$\displaystyle-(1-2m+4n(1+m))C^{2}\Big{)}\
{}_{2}F_{1}\Big{(}\frac{1}{2}-n,-n;\frac{1}{2}-m-n;\frac{\rho^{2}}{A^{2}}\Big{)}$
$\displaystyle-(1-2n)\Big{(}2(1-2m+2n)(1-m+n)A^{2}$
$\displaystyle-(1-2m+4n(1+m))C^{2}\Big{)}\
{}_{2}F_{1}\Big{(}\frac{3}{2}-n,1-n;\frac{1}{2}-m-n;\frac{\rho^{2}}{A^{2}}\Big{)}\Bigg{]}$

The derivatives with respect to $C$ are

$\displaystyle\frac{d\tilde{G}^{(0)}_{0,m}}{dC}=-\frac{C}{\rho^{2}}\sum_{n=m}^{\infty}\Big{(}\frac{BA^{2}}{\rho^{2}}\Big{)}^{n}f_{nm}\Big{[}$
$\displaystyle{}_{2}F_{1}\Big{(}\frac{1}{2}-n,-n;\frac{1}{2}-m-n;\frac{\rho^{2}}{A^{2}}\Big{)}$
$\displaystyle-(1-2n)\
{}_{2}F_{1}\Big{(}\frac{3}{2}-n,1-n;\frac{1}{2}-m-n;\frac{\rho^{2}}{A^{2}}\Big{)}\Big{]}$
(C.13)
$\displaystyle\frac{d^{2}\tilde{G}^{(0)}_{0,m}}{dC^{2}}=\rho^{-4}\sum_{n=m}^{\infty}\Big{(}$
$\displaystyle\frac{BA^{2}}{\rho^{2}}\Big{)}^{n}f_{nm}\Big{[}\big{(}(1+2m-4n(1+m))A^{2}+2C^{2}\big{)}{}_{2}F_{1}\Big{(}\frac{1}{2}-n,-n;\frac{1}{2}-m-n;\frac{\rho^{2}}{A^{2}}\Big{)}$
$\displaystyle-(1-2n)\big{(}(1+2m-2n)A^{2}+2C^{2}\big{)}{}_{2}F_{1}\Big{(}\frac{3}{2}-n,-n;\frac{1}{2}-m-n;\frac{\rho^{2}}{A^{2}}\Big{)}\Big{]}.$
(C.14)

Note that we can use $dG/dA=-(C/A)\ dG/dC$ to recast the first derivative
w.r.t. $A$ in a convenient form as well.

In addition, we need two mixed derivatives
$\partial_{C}\partial_{A}^{(1,2)}G$. These are

$\displaystyle\frac{\partial^{2}\tilde{G}^{(0)}_{0,m}}{\partial C\partial
A}=-\frac{C}{A\rho^{4}}\sum_{n=m}^{\infty}\Big{(}$
$\displaystyle\frac{BA^{2}}{\rho^{2}}\Big{)}^{n}f_{nm}\Big{(}(2(m-2n(1+m))A^{2}+C^{2}){}_{2}F_{1}(\frac{1}{2}-n,-n;\frac{1}{2}-m-n;\frac{\rho^{2}}{A^{2}})$
$\displaystyle-(1-2n)(2(m-n)A^{2}+C^{2}){}_{2}F_{1}(\frac{3}{2}-n,-n;\frac{1}{2}-m-n;\frac{\rho^{2}}{A^{2}})\Big{)}$
(C.15) $\displaystyle\frac{\partial^{3}\tilde{G}^{(0)}_{0,m}}{\partial
C\partial A^{2}}=\frac{C}{\rho^{6}}\sum_{n=m}^{\infty}\Big{(}$
$\displaystyle\frac{BA^{2}}{\rho^{2}}\Big{)}^{n}f_{nm}\Big{(}\big{(}2(m-2m^{2}-4n(1-m^{2})-4n^{2}(1+m))A^{2}$
$\displaystyle+3(1-2m+4n(1+m))C^{2}\big{)}\
{}_{2}F_{1}(\frac{1}{2}-n,-n;\frac{1}{2}-m-n;\frac{\rho^{2}}{A^{2}})$
$\displaystyle-(1-2n)\big{(}2(1-2m+2n)(m-n)A^{2}$
$\displaystyle+(3-6m+8n+4mn)C^{2}\big{)}\
{}_{2}F_{1}(\frac{3}{2}-n,-n;\frac{1}{2}-m-n;\frac{\rho^{2}}{A^{2}})\Big{)}$
(C.16)

###  C.5 Implementation in Python

Our implementation of Method II, lpt_rsd_fftw.py, is available as part of
velocileptors131313https://github.com/sfschen/velocileptors, a Python package
for the one-loop redshift-space power spectrum that also includes modules
implementing the moment expansion, Gaussian streaming model, and resummed
Eulerian perturbation theory. The LPT module includes auxiliary functions to
compute multipoles, add Alcock-Paczynski parameters not equal to one and
combine the various bias contributions. We also include a sample jupyter
notebook containing example usage.

Our IR-resummation procedure is inherently rather numerically involved because
of the angular dependence of the resummed displacements. To speed up this
calculation (in Python) we can take advantage of the fact that the derivatives
$\partial^{b}_{C}\partial^{a}_{A}\tilde{G}^{(0)}_{0,\ell}$ can be written as

$\sum_{n=\ell}^{\infty}\Big{[}-\frac{1}{2}K^{2}Y^{\rm
lin}(q)\Big{]}^{n}(kq)^{-(a+b)}c^{2n}F_{n}(f,\mu),$

of which the only vector operations involve multiplying by $q$ and $Y(q)$. The
remaining factors are independent of $k$ and $q$ and can be tabulated for each
value of $\mu$. Further, the only special functions we need are
${}_{2}F_{1}(\frac{1}{2}-n,-n;\frac{1}{2}-m-n;x)$,
${}_{2}F_{1}(\frac{3}{2}-n,-n;\frac{1}{2}-m-n;x)$, and the $\Gamma$-functions
in $f_{nm}$, which can all be tabulated in advance and do not have to be
calculated at each wavenumber. With these simplifications, it takes about two
and a half seconds to compute the mutipoles at 50 $k$ points between $0.01$
and $0.25\,h\,{\rm Mpc}^{-1}$; applying a cubic spline to interpolate between
these points is sufficient to achieve sub-percent accuracy for any $k$ in this
range. Similarly, it takes less than a second to compute $P(k,\mu)$ over the
same number of points for a fixed $\mu$.

###  C.6 An alternative formulation of Method II

In this section we provide an alternative numerical solution to the direct
evaluation of the RSD integrals in LPT. For simplicity we will present only
the leading order term, i.e. the Zeldovich approximation, but the main result
is trivially extended to one-loop, for which we will provide the necessary
ingredients. We begin by writing the Zeldovich RSD power spectrum as

$\displaystyle P_{s,\rm
Zel}(k,\mu)=2\pi\int\differential{q}q^{2}e^{-\frac{1}{2}K^{2}(X(q)+Y(q))}\int_{-1}^{1}d\mu_{\mathbf{q}}~{}e^{i\mu_{\mathbf{q}}A+(\mu_{\mathbf{q}}^{2}-1)B}J_{0}(C\sqrt{1-\mu_{\mathbf{q}}})$
(C.17)

with $A\equiv kqc$, $B\equiv-K^{2}Y(q)/2$ and $C\equiv kqs$. We then Taylor
series expand in $B$ and integrate the $A$ piece by parts $n$ times, when $n$
goes to infinity eventually, to rewrite the integral over $\mu_{\mathbf{q}}$
as

$\displaystyle\sum_{n=0}^{\infty}(-1)^{n}\int_{-1}^{1}d\mu_{\mathbf{q}}~{}\frac{e^{i\mu_{\mathbf{q}}A}}{(iA/B)^{n}}2^{n}\frac{d^{n}}{d\mu_{\mathbf{q}}^{n}}\left(\frac{(\mu_{\mathbf{q}}^{2}-1)^{n}}{2^{n}n!}J_{0}(C\sqrt{1-\mu_{\mathbf{q}}^{2}})\right)\,.$
(C.18)

Then use 10.1.48 of ref. [110] to expand the $J_{0}$,

$\displaystyle
J_{0}(C\sqrt{1-\mu_{\mathbf{q}}^{2}})=\sum_{\alpha=0}^{\infty}(4\alpha+1)\frac{(2\alpha)!}{2^{2\alpha}(\alpha!)^{2}}j_{2\alpha}(C)\mathcal{P}_{2\alpha}(\mu_{\mathbf{q}})$
(C.19)

to arrive, using the plane wave expansion of the exponential, at

$\displaystyle=\sum_{n,\alpha,\ell}(-1)^{n}\frac{2^{n}(2\alpha)!}{2^{2\alpha}(\alpha!)^{2}}(4\alpha+1)(2\ell+1)\frac{(i)^{\ell}j_{\ell}(A)j_{2\alpha}(C)}{(iA/B)^{n}}$
(C.20)
$\displaystyle\times\int_{-1}^{1}d\mu_{\mathbf{q}}~{}\mathcal{P}_{\ell}(\mu_{\mathbf{q}})\frac{d^{n}}{d\mu_{\mathbf{q}}^{n}}\left(\frac{(\mu_{\mathbf{q}}^{2}-1)^{n}}{2^{n}n!}\mathcal{P}_{2\alpha}(\mu_{\mathbf{q}})\right)\,.$
(C.21)

Now it turns out that

$\displaystyle=\int_{-1}^{1}d\mu_{\mathbf{q}}~{}\mathcal{P}_{\ell}(\mu_{\mathbf{q}})\sum_{k=0}^{n}\binom{n}{k}\mathcal{P}_{n}^{-k}(\mu_{\mathbf{q}})\mathcal{P}_{2\alpha}^{k}(\mu_{\mathbf{q}})$
(C.22)
$\displaystyle=2\sum_{k=0}^{n}\binom{n}{k}\sqrt{\frac{(2\alpha+k)!(n-k)!}{(n+k)!(2\alpha-k)!}}\begin{pmatrix}\ell&n&2\alpha\\\
0&-k&k\end{pmatrix}\begin{pmatrix}\ell&n&2\alpha\\\ 0&0&0\end{pmatrix}$

in terms of 3-$j$ symbols. Putting the above equations together, the angular
part of the Zeldovich RSD integral can be computed analytically

$\displaystyle=\sum_{n,\alpha,\ell}(-1)^{n}\frac{2^{n}(2\alpha)!}{2^{2\alpha}(\alpha!)^{2}}(4\alpha+1)(2\ell+1)\frac{(i)^{\ell}j_{\ell}(A)j_{2\alpha}(C)}{(iA/B)^{n}}$
(C.23) $\displaystyle\times
2\sum_{k}\binom{n}{k}\sqrt{\frac{(2\alpha+k)!(n-k)!}{(n+k)!(2\alpha-k)!}}\begin{pmatrix}\ell&n&2\alpha\\\
0&-k&k\end{pmatrix}\begin{pmatrix}\ell&n&2\alpha\\\ 0&0&0\end{pmatrix}$ (C.24)
$\displaystyle\equiv\sum_{\ell,\alpha}c_{\ell,\alpha}(k,\mu,q)j_{\ell}(A)j_{2\alpha}(C)\,,$
(C.25)

and we can rewrite ZA power spectrum in redshift space as

$\displaystyle P_{s,\rm Zel}(k,\mu)=$ $\displaystyle
2\pi\sum_{\ell,\alpha=0}^{\infty}\int\frac{\differential{q}}{q}q^{3}e^{-1/2K^{2}(X(q)+Y(q))}c_{\ell,\alpha}(k,\mu,q)j_{\ell}(A)j_{2\alpha}(C)\,.$
(C.26)

The two remaining sums over $\alpha$ and $\ell$ run from zero to infinity, but
in practice only the first 5 terms are relevant for sub-% precision. Upon
expanding the non oscillatory part of the integrand above in complex power
laws (FFTlog) [82, 111],

$\displaystyle
q^{3}e^{-1/2K^{2}(X(q)+Y(q))}c_{\ell,\alpha}(k,\mu,q)\equiv\sum_{n}d_{\ell,\alpha,n}(k,\mu)q^{\nu_{n}}\,$
(C.27)

and then using the following analytic integral

$\displaystyle\int\frac{\differential{q}}{q}q^{\nu_{n}}j_{\ell}(kqc)j_{2\alpha}(kqs)$
(C.28) $\displaystyle=$ $\displaystyle(kc)^{-\nu_{n}}\frac{\pi
2^{\nu_{n}-3}t^{2\alpha}\Gamma\left(\frac{1}{2}(2\alpha+\ell+\nu_{n})\right)\,_{2}\tilde{F}_{1}\left(\frac{1}{2}(2\alpha-\ell+\nu_{n}-1),\frac{1}{2}(2\alpha+\ell+\nu_{n});2\alpha+\frac{3}{2};t^{2}\right)}{\Gamma\left(\frac{1}{2}(-2\alpha+\ell-\nu_{n}+3)\right)}$
(C.29) $\displaystyle\equiv(kc)^{-\nu_{n}}I(\ell,\alpha,\nu_{n},t)$ (C.30)

with $t\equiv s/c$, and ${}_{2}\tilde{F}_{1}$ a regularized Hypergeometric
function, we can perform the remaining integral in Eq. C.26 141414See for
example p. 401 (13.4) of ref. [112]. The integral exists for any $\nu_{n}>-2$,
which is the case by an appropriate choice of the FFTlog parameters.. The
Zeldovich RSD power spectrum can then be computed as

$\displaystyle P_{\rm
ZA}(k,\mu)=2\pi\sum_{\ell,\alpha,n}d_{\ell,\alpha,n}(k,\mu)(kc)^{-\nu_{n}}I(\ell,\alpha,\nu_{n},t)\,.$
(C.31)

For a given set of FFTlog parameters, the function $I(\ell,\alpha,\nu_{n},t)$
can be tabulated in advance and the sum above performed quickly.

Loop integrals and bias terms introduce two main complications. First, new
terms proportional to $\mu_{\mathbf{q}}^{\beta}$ appear in the angular
integral in Equation C.17. They can be easily included by rewriting them as
derivatives with respect to $A$ or $B$, analogously to what is done in the
real space calculation [15]. Second, the loop structure will replace the
$J_{0}$ in Equation C.17 with more complicated functions. For the one-loop
calculation at hand, we have to perform the angular integral in Equation C.18
with the following two terms

$\displaystyle-i\sqrt{1-\mu_{\mathbf{q}}^{2}}\
J_{1}(C\sqrt{1-\mu_{\mathbf{q}}^{2}})\quad\;\text{and}\;\quad\frac{1}{2}(1-\mu_{\mathbf{q}}^{2})\left[J_{0}(C\sqrt{1-\mu_{\mathbf{q}}^{2}})-J_{2}(C\sqrt{1-\mu_{\mathbf{q}}^{2}})\right]$
(C.32)

instead of the $J_{0}$. Those are easy to deal with by noticing that

$\displaystyle-i\sqrt{1-\mu_{\mathbf{q}}^{2}}\
J_{1}(C\sqrt{1-\mu_{\mathbf{q}}^{2}})=i\partial_{C}J_{0}(C\sqrt{1-\mu_{\mathbf{q}}^{2}})$
$\displaystyle\longrightarrow i\partial_{C}j_{2\alpha}(C)$ (C.33)
$\displaystyle=i\frac{2aj_{2a}(C)}{C}-ij_{2a+1}(C)$ (C.34)

and

$\displaystyle\frac{1}{2}(1-\mu_{\mathbf{q}}^{2})\left[J_{0}(C\sqrt{1-\mu_{\mathbf{q}}^{2}})-J_{2}(C\sqrt{1-\mu_{\mathbf{q}}^{2}})\right]$
$\displaystyle=-\partial^{2}_{C}J_{0}(C\sqrt{1-\mu_{\mathbf{q}}^{2}})\longrightarrow-\partial^{2}_{C}j_{2\alpha}(C)$
(C.35)
$\displaystyle=-\frac{\left(4a^{2}-2a-C^{2}\right)j_{2a}(C)+2Cj_{2a+1}(C)}{C^{2}}\,,$
(C.36)

which boil down to a reshuffling of the coefficients in the sums of Equation
C.25.

## References

  * [1] N. Kaiser, _Clustering in real space and in redshift space_ , _MNRAS_ 227 (1987) 1.
  * [2] A. J. S. Hamilton, _Measuring Omega and the real correlation function from the redshift correlation function_ , _ApJL_ 385 (1992) L5.
  * [3] J. A. Peacock, _Cosmological Physics_. Jan., 1999.
  * [4] S. Dodelson, _Modern cosmology_. 2003.
  * [5] L. Guzzo et al., _A test of the nature of cosmic acceleration using galaxy redshift distortions_ , _Nature_ 451 (2008) 541 [0802.1944].
  * [6] D. H. Weinberg, M. J. Mortonson, D. J. Eisenstein, C. Hirata, A. G. Riess and E. Rozo, _Observational probes of cosmic acceleration_ , _PhysRep_ 530 (2013) 87 [1201.2434].
  * [7] A. Joyce, B. Jain, J. Khoury and M. Trodden, _Beyond the cosmological standard model_ , _PhysRep_ 568 (2015) 1 [1407.0059].
  * [8] L. Amendola, S. Appleby, A. Avgoustidis, D. Bacon, T. Baker, M. Baldi et al., _Cosmology and fundamental physics with the Euclid satellite_ , _Living Reviews in Relativity_ 21 (2018) 2 [1606.00180].
  * [9] S. Alam, A. Aviles, R. Bean, Y.-C. Cai, M. Cautun, J. L. Cervantes-Cota et al., _Testing the theory of gravity with DESI: estimators, predictions and simulation requirements_ , _arXiv e-prints_ (2020) arXiv:2011.05771 [2011.05771].
  * [10] DESI Collaboration, A. Aghamousa, J. Aguilar, S. Ahlen, S. Alam, L. E. Allen et al., _The DESI Experiment Part I: Science,Targeting, and Survey Design_ , _ArXiv e-prints_ (2016) [1611.00036].
  * [11] R. Laureijs, J. Amiaux, S. Arduini, J. . Auguères, J. Brinchmann, R. Cole et al., _Euclid Definition Study Report_ , _ArXiv e-prints_ (2011) [1110.3193].
  * [12] D. Baumann, A. Nicolis, L. Senatore and M. Zaldarriaga, _Cosmological non-linearities as an effective fluid_ , _JCAP_ 7 (2012) 51 [1004.2488].
  * [13] J. J. M. Carrasco, M. P. Hertzberg and L. Senatore, _The effective field theory of cosmological large scale structures_ , _Journal of High Energy Physics_ 9 (2012) 82 [1206.2926].
  * [14] R. A. Porto, L. Senatore and M. Zaldarriaga, _The Lagrangian-space Effective Field Theory of large scale structures_ , _JCAP_ 5 (2014) 022 [1311.2168].
  * [15] Z. Vlah, M. White and A. Aviles, _A Lagrangian effective field theory_ , _JCAP_ 9 (2015) 014 [1506.05264].
  * [16] L. Senatore, _Bias in the effective field theory of large scale structures_ , _JCAP_ 11 (2015) 007 [1406.7843].
  * [17] R. Angulo, M. Fasiello, L. Senatore and Z. Vlah, _On the statistics of biased tracers in the Effective Field Theory of Large Scale Structures_ , _JCAP_ 9 (2015) 029 [1503.08826].
  * [18] T. Matsubara, _Resumming cosmological perturbations via the Lagrangian picture: One-loop results in real space and in redshift space_ , _PRD_ 77 (2008) 063530 [0711.2521].
  * [19] J. Carlson, B. Reid and M. White, _Convolution Lagrangian perturbation theory for biased tracers_ , _MNRAS_ 429 (2013) 1674 [1209.0780].
  * [20] L. Wang, B. Reid and M. White, _An analytic model for redshift-space distortions_ , _MNRAS_ 437 (2014) 588 [1306.1804].
  * [21] M. White, _The Zel’dovich approximation_ , _MNRAS_ 439 (2014) 3630 [1401.5466].
  * [22] Z. Vlah, E. Castorina and M. White, _The Gaussian streaming model and convolution Lagrangian effective field theory_ , _JCAP_ 12 (2016) 007 [1609.02908].
  * [23] Z. Vlah and M. White, _Exploring redshift-space distortions in large-scale structure_ , _JCAP_ 2019 (2019) 007 [1812.02775].
  * [24] S.-F. Chen, Z. Vlah and M. White, _The reconstructed power spectrum in the Zeldovich approximation_ , _JCAP_ 2019 (2019) 017 [1907.00043].
  * [25] B. A. Reid, L. Samushia, M. White, W. J. Percival, M. Manera, N. Padmanabhan et al., _The clustering of galaxies in the SDSS-III Baryon Oscillation Spectroscopic Survey: measurements of the growth of structure and expansion rate at z = 0.57 from anisotropic clustering_ , _MNRAS_ 426 (2012) 2719 [1203.6641].
  * [26] L. Samushia, B. A. Reid, M. White, W. J. Percival, A. J. Cuesta, G.-B. Zhao et al., _The clustering of galaxies in the SDSS-III Baryon Oscillation Spectroscopic Survey: measuring growth rate and geometry with anisotropic clustering_ , _MNRAS_ 439 (2014) 3504 [1312.4899].
  * [27] S. Satpathy, S. Alam, S. Ho, M. White, N. A. Bahcall, F. Beutler et al., _The clustering of galaxies in the completed SDSS-III Baryon Oscillation Spectroscopic Survey: on the measurement of growth rate using galaxy correlation functions_ , _MNRAS_ 469 (2017) 1369 [1607.03148].
  * [28] A. Tamone, A. Raichoor, C. Zhao, A. de Mattia, C. Gorgoni, E. Burtin et al., _The completed SDSS-IV extended baryon oscillation spectroscopic survey: growth rate of structure measurement from anisotropic clustering analysis in configuration space between redshift 0.6 and 1.1 for the emission-line galaxy sample_ , _MNRAS_ 499 (2020) 5527 [2007.09009].
  * [29] Y. Zhang, A. R. Pullen, S. Alam, S. Singh, E. Burtin, C.-H. Chuang et al., _Testing General Relativity on cosmological scales at redshift z $\sim$ 1.5 with quasar and CMB lensing_, _MNRAS_ (2020)  [2007.12607].
  * [30] J. E. Bautista, R. Paviot, M. Vargas Magaña, S. de la Torre, S. Fromenteau, H. Gil-Marín et al., _The completed SDSS-IV extended Baryon Oscillation Spectroscopic Survey: measurement of the BAO and growth rate of structure of the luminous red galaxy sample from the anisotropic correlation function between redshifts 0.6 and 1_ , _MNRAS_ 500 (2021) 736 [2007.08993].
  * [31] L. F. de la Bella, D. Regan, D. Seery and S. Hotchkiss, _The matter power spectrum in redshift space using effective field theory_ , _JCAP_ 11 (2017) 039 [1704.05309].
  * [32] L. Fonseca de la Bella, D. Regan, D. Seery and D. Parkinson, _Impact of bias and redshift-space modelling for the halo power spectrum: Testing the effective field theory of large-scale structure_ , _JCAP_ 07 (2020) 011 [1805.12394].
  * [33] V. Desjacques, D. Jeong and F. Schmidt, _The galaxy power spectrum and bispectrum in redshift space_ , _JCAP_ 2018 (2018) 035 [1806.04015].
  * [34] A. Chudaykin, M. M. Ivanov, O. H. E. Philcox and M. Simonović, _Nonlinear perturbation theory extension of the Boltzmann code CLASS_ , _PRD_ 102 (2020) 063533 [2004.10607].
  * [35] 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_ 2020 (2020) 005 [1909.05271].
  * [36] J. Tomlinson, H. S. G. Gebhardt and D. Jeong, _Fast calculation of the nonlinear redshift-space galaxy power spectrum including selection bias_ , _Phys. Rev. D_ 101 (2020) 103528 [2004.03629].
  * [37] U. Seljak and P. McDonald, _Distribution function approach to redshift space distortions_ , _JCAP_ 11 (2011) 039 [1109.1888].
  * [38] Z. Vlah, U. Seljak, P. McDonald, T. Okumura and T. Baldauf, _Distribution function approach to redshift space distortions. Part IV: perturbation theory applied to dark matter_ , _JCAP_ 11 (2012) 009 [1207.0839].
  * [39] Z. Vlah, U. Seljak, T. Okumura and V. Desjacques, _Distribution function approach to redshift space distortions. Part V: perturbation theory applied to dark matter halos_ , _JCAP_ 10 (2013) 053 [1308.6294].
  * [40] L. Senatore and M. Zaldarriaga, _The IR-resummed Effective Field Theory of Large Scale Structures_ , _JCAP_ 2 (2015) 13 [1404.5954].
  * [41] L. Senatore and G. Trevisan, _On the IR-Resummation in the EFTofLSS_ , _JCAP_ 05 (2018) 019 [1710.02178].
  * [42] T. Baldauf, M. Mirbabayi, M. Simonović and M. Zaldarriaga, _Equivalence principle and the baryon acoustic peak_ , _PRD_ 92 (2015) 043514 [1504.04366].
  * [43] Z. Vlah, U. Seljak, M. Yat Chu and Y. Feng, _Perturbation theory, effective field theory, and oscillations in the power spectrum_ , _Journal of Cosmology and Astro-Particle Physics_ 2016 (2016) 057 [1509.02120].
  * [44] D. Blas, M. Garny, M. M. Ivanov and S. Sibiryakov, _Time-sliced perturbation theory II: baryon acoustic oscillations and infrared resummation_ , _JCAP_ 2016 (2016) 028 [1605.02149].
  * [45] Z. Ding, H.-J. Seo, Z. Vlah, Y. Feng, M. Schmittfull and F. Beutler, _Theoretical systematics of Future Baryon Acoustic Oscillation Surveys_ , _MNRAS_ 479 (2018) 1021 [1708.01297].
  * [46] M. M. Ivanov and S. Sibiryakov, _Infrared resummation for biased tracers in redshift space_ , _JCAP_ 2018 (2018) 053 [1804.05080].
  * [47] M. Lewandowski, L. Senatore, F. Prada, C. Zhao and C.-H. Chuang, _EFT of large scale structures in redshift space_ , _Phys. Rev. D_ 97 (2018) 063526 [1512.06831].
  * [48] A. Perko, L. Senatore, E. Jennings and R. H. Wechsler, _Biased Tracers in Redshift Space in the EFT of Large-Scale Structure_ , 1610.09321.
  * [49] S.-F. Chen, Z. Vlah and M. White, _Consistent modeling of velocity statistics and redshift-space distortions in one-loop perturbation theory_ , _JCAP_ 2020 (2020) 062 [2005.00523].
  * [50] H.-M. Zhu, Y. Yu and U.-L. Pen, _Nonlinear reconstruction of redshift space distortions_ , _Phys. Rev. D_ 97 (2018) 043502 [1711.03218].
  * [51] C. Modi, M. White, A. Slosar and E. Castorina, _Reconstructing large-scale structure with neutral hydrogen surveys_ , _JCAP_ 11 (2019) 023 [1907.02330].
  * [52] M. Schmittfull, M. Simonović, M. M. Ivanov, O. H. E. Philcox and M. Zaldarriaga, _Modeling Galaxies in Redshift Space at the Field Level_ , _arXiv e-prints_ (2020) arXiv:2012.03334 [2012.03334].
  * [53] T. Matsubara, _Nonlinear perturbation theory with halo bias and redshift-space distortions via the Lagrangian picture_ , _PRD_ 78 (2008) 083519 [0807.1733].
  * [54] T. Matsubara, _Recursive solutions of Lagrangian perturbation theory_ , _PRD_ 92 (2015) 023534 [1505.01481].
  * [55] V. Zheligovsky and U. Frisch, _Time-analyticity of Lagrangian particle trajectories in ideal fluid flow_ , _Journal of Fluid Mechanics_ 749 (2014) 404 [1312.6320].
  * [56] C. Rampf, _The recursion relation in Lagrangian perturbation theory_ , _JCAP_ 2012 (2012) 004 [1205.5274].
  * [57] F. Beutler, U. Seljak and Z. Vlah, _Constraining the relative velocity effect using the Baryon Oscillation Spectroscopic Survey_ , _ArXiv e-prints_ (2016) [1612.04720].
  * [58] F. Schmidt, _Effect of relative velocity and density perturbations between baryons and dark matter on the clustering of galaxies_ , _PRD_ 94 (2016) 063508 [1602.09059].
  * [59] S.-F. Chen, E. Castorina and M. White, _Biased tracers of two fluids in the Lagrangian picture_ , _Journal of Cosmology and Astro-Particle Physics_ 2019 (2019) 006 [1903.00437].
  * [60] A. Barreira, G. Cabass, D. Nelson and F. Schmidt, _Baryon-CDM isocurvature galaxy bias with IllustrisTNG_ , _JCAP_ 02 (2020) 005 [1907.04317].
  * [61] H. Khoraminezhad, T. Lazeyras, R. E. Angulo, O. Hahn and M. Viela, _Quantifying the impact of baryon-CDM perturbations on halo clustering and baryon fraction_ , 2011.01037.
  * [62] C. Rampf, C. Uhlemann and O. Hahn, _Cosmological perturbations for two cold fluids in $\Lambda$CDM_, _arXiv e-prints_ (2020) arXiv:2008.09123 [2008.09123].
  * [63] C. Rampf, B. Villone and U. Frisch, _How smooth are particle trajectories in a $\Lambda$CDM Universe?_, _MNRAS_ 452 (2015) 1421 [1504.00032].
  * [64] M. Fasiello and Z. Vlah, _Nonlinear fields in generalized cosmologies_ , _PRD_ 94 (2016) 063516 [1604.04612].
  * [65] T. Fujita and Z. Vlah, _Perturbative description of biased tracers using consistency relations of LSS_ , _JCAP_ 10 (2020) 059 [2003.10114].
  * [66] E. Castorina and M. White, _Beyond the plane-parallel approximation for redshift surveys_ , _MNRAS_ 476 (2018) 4403 [1709.09730].
  * [67] E. Castorina and M. White, _The Zeldovich approximation and wide-angle redshift-space distortions_ , _MNRAS_ (2018)  [1803.08185].
  * [68] A. Taruya, S. Saga, M.-A. Breton, Y. Rasera and T. Fujita, _Wide-angle redshift-space distortions at quasi-linear scales: cross-correlation functions from Zel’dovich approximation_ , _Mon. Not. Roy. Astron. Soc._ 491 (2020) 4162 [1908.03854].
  * [69] S. Bharadwaj, _The Evolution of Correlation Functions in the Zeldovich Approximation and Its Implications for the Validity of Perturbation Theory_ , _ApJ_ 472 (1996) 1 [arXiv:astro-ph/9606121].
  * [70] D. J. Eisenstein, H.-J. Seo and M. White, _On the Robustness of the Acoustic Scale in the Low-Redshift Clustering of Matter_ , _ApJ_ 664 (2007) 660 [astro-ph/0604361].
  * [71] M. Crocce and R. Scoccimarro, _Nonlinear evolution of baryon acoustic oscillations_ , _PRD_ 77 (2008) 023533 [0704.2783].
  * [72] N. Padmanabhan, M. White and J. D. Cohn, _Reconstructing baryon oscillations: A Lagrangian theory perspective_ , _PRD_ 79 (2009) 063523 [0812.2905].
  * [73] Y. Noh, M. White and N. Padmanabhan, _Reconstructing baryon oscillations_ , _PRD_ 80 (2009) 123501 [0909.1802].
  * [74] S. Tassev and M. Zaldarriaga, _Towards an optimal reconstruction of baryon oscillations_ , _JCAP_ 10 (2012) 006 [1203.6066].
  * [75] N. McCullagh and A. S. Szalay, _Nonlinear Behavior of Baryon Acoustic Oscillations from the Zel’dovich Approximation Using a Non-Fourier Perturbation Approach_ , _ApJ_ 752 (2012) 21 [1202.1306].
  * [76] M. White, _The Zel’dovich approximation_ , _MNRAS_ 439 (2014) 3630 [1401.5466].
  * [77] M. Schmittfull, Y. Feng, F. Beutler, B. Sherwin and M. Y. Chu, _Eulerian BAO reconstructions and N -point statistics_ , _PRD_ 92 (2015) 123522 [1508.06972].
  * [78] M. McQuinn and M. White, _Cosmological perturbation theory in 1+1 dimensions_ , _JCAP_ 1 (2016) 043 [1502.07389].
  * [79] H.-J. Seo, F. Beutler, A. J. Ross and S. Saito, _Modeling the reconstructed BAO in Fourier space_ , _MNRAS_ 460 (2016) 2453 [1511.00663].
  * [80] N. E. Chisari and A. Pontzen, _Unequal time correlators and the Zeldovich approximation_ , _arXiv e-prints_ (2019) arXiv:1905.02078 [1905.02078].
  * [81] L. F. de la Bella, N. Tessore and S. Bridle, _The unequal-time matter power spectrum: impact on weak lensing observables_ , 2011.06185.
  * [82] A. J. S. Hamilton, _Uncorrelated modes of the non-linear power spectrum_ , _MNRAS_ 312 (2000) 257 [astro-ph/9905191].
  * [83] Z. Vlah, U. Seljak and T. Baldauf, _Lagrangian perturbation theory at one loop order: Successes, failu res, and improvements_ , _PRD_ 91 (2015) 023508 [1410.1617].
  * [84] D. N. Limber, _The Analysis of Counts of the Extragalactic Nebulae in Terms of a Fluctuating Density Field._ , _ApJ_ 117 (1953) 134.
  * [85] M. LoVerde and N. Afshordi, _Extended Limber Approximation_ , _Phys. Rev. D_ 78 (2008) 123506 [0809.5112].
  * [86] A. N. Taylor and A. J. S. Hamilton, _Non-linear cosmological power spectra in real and redshift space_ , _MNRAS_ 282 (1996) 767 [astro-ph/9604020].
  * [87] P. J. E. Peebles, _The large-scale structure of the universe_. 1980.
  * [88] K. B. Fisher, _On the Validity of the Streaming Model for the Redshift-Space Correlation Function in the Linear Regime_ , _ApJ_ 448 (1995) 494 [astro-ph/9412081].
  * [89] B. A. Reid and M. White, _Towards an accurate model of the redshift-space clustering of haloes in the quasi-linear regime_ , _MNRAS_ 417 (2011) 1913 [1105.4165].
  * [90] T. Sunayama, N. Padmanabhan, K. Heitmann, S. Habib and E. Rangel, _Efficient construction of mock catalogs for baryon acoustic oscillation surveys_ , _JCAP_ 2016 (2016) 051 [1510.06665].
  * [91] N. Hand, Y. Feng, F. Beutler, Y. Li, C. Modi, U. Seljak et al., _nbodykit: An Open-source, Massively Parallel Toolkit for Large-scale Structure_ , _AJ_ 156 (2018) 160 [1712.05834].
  * [92] J. A. Blazek, J. E. McEwen and C. M. Hirata, _Streaming Velocities and the Baryon Acoustic Oscillation Scale_ , _Physical Review Letters_ 116 (2016) 121303 [1510.03554].
  * [93] S.-F. Chen, Z. Vlah and M. White, _Modeling features in the redshift-space halo power spectrum with perturbation theory_ , _JCAP_ 2020 (2020) 035 [2007.00704].
  * [94] T. Nishimichi, G. D’Amico, M. M. Ivanov, L. Senatore, M. Simonović, M. Takada et al., _Blinded challenge for precision cosmology with large-scale structure: results from effective field theory for the redshift-space galaxy power spectrum_ , _arXiv e-prints_ (2020) arXiv:2003.08277 [2003.08277].
  * [95] S. Alam, M. Ata, S. Bailey, F. Beutler, D. Bizyaev, J. A. Blazek et al., _The clustering of galaxies in the completed SDSS-III Baryon Oscillation Spectroscopic Survey: cosmological analysis of the DR12 galaxy sample_ , _MNRAS_ 470 (2017) 2617 [1607.03155].
  * [96] C. Alcock and B. Paczynski, _An evolution free test for non-zero cosmological constant_ , _Nature_ 281 (1979) 358.
  * [97] C. M. Hirata, _Tidal alignments as a contaminant of redshift space distortions_ , _MNRAS_ 399 (2009) 1074 [0903.4929].
  * [98] J. J. Givans and C. M. Hirata, _Redshift-space streaming velocity effects on the Lyman- $\alpha$ forest baryon acoustic oscillation scale_, _PRD_ 102 (2020) 023515 [2002.12296].
  * [99] A. Obuljen, W. J. Percival and N. Dalal, _Detection of anisotropic galaxy assembly bias in BOSS DR12_ , _JCAP_ 2020 (2020) 058 [2004.07240].
  * [100] 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].
  * [101] 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].
  * [102] 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].
  * [103] F. Villaescusa-Navarro, A. Banerjee, N. Dalal, E. Castorina, R. Scoccimarro, R. Angulo et al., _The imprint of neutrinos on clustering in redshift-space_ , _Astrophys. J._ 861 (2018) 53 [1708.01154].
  * [104] M. LoVerde, _Halo bias in mixed dark matter cosmologies_ , _Phys. Rev. D_ 90 (2014) 083530 [1405.4855].
  * [105] J. B. Muñoz and C. Dvorkin, _Efficient Computation of Galaxy Bias with Neutrinos and Other Relics_ , _Phys. Rev. D_ 98 (2018) 043503 [1805.11623].
  * [106] C. Fidler, N. Sujata and M. Archidiacono, _Relativistic bias in neutrino cosmologies_ , _JCAP_ 06 (2019) 035 [1812.09266].
  * [107] A. Aviles and A. Banerjee, _A Lagrangian Perturbation Theory in the presence of massive neutrinos_ , _JCAP_ 10 (2020) 034 [2007.06508].
  * [108] A. Aviles and J. L. Cervantes-Cota, _A Lagrangian perturbation theory for modify gravity_ , _ArXiv e-prints_ (2017) [1705.10719].
  * [109] G. Valogiannis, R. Bean and A. Aviles, _An accurate perturbative approach to redshift space clustering of biased tracers in modified gravity_ , _JCAP_ 01 (2020) 055 [1909.05261].
  * [110] M. Abramowitz and I. A. Stegun, _Handbook of Mathematical Functions_. 1972.
  * [111] V. Assassi, M. Simonović and M. Zaldarriaga, _Efficient evaluation of angular power spectra and bispectra_ , _JCAP_ 11 (2017) 054 [1705.05022].
  * [112] G. N. Watson, _A Treatise on the Theory of Bessel Functions, 2nd ed_. 1966\.