# A Lagrangian theory for galaxy shape statistics

Shi-Fan Chen  Nickolas Kokron

###### Abstract

We formulate the Lagrangian perturbation theory of galaxy intrinsic alignments
and compute the resulting auto and cross power spectra of galaxy shapes,
densities and matter to 1-loop order. Our model represents a consistent
effective-theory description of galaxy shape including the resummation of
long-wavelength displacements which damp baryon acoustic oscillations, and
includes one linear, three quadratic and two cubic dimensionless bias
coefficients at this order, along with counterterms and stochastic
contributions whose structure we derive. We compare this Lagrangian model
against the three-dimensional helicity spectra of halo shapes measured in
N-body simulations by ref. [1] and find excellent agreement on perturbative
scales while testing a number of more restrictive bias parametrizations. The
calculations presented are immediately relevant to analyses of both cosmic
shear surveys and spectroscopic shape measurements, and we make a fast FFTLog-
based code spinosaurus publicly available with this publication.

##  1 Introduction

The last decade of cosmic shear surveys such as the Dark Energy Survey [2],
the Kilo-Degree Survey [3] and the Hyper Suprime Cam Strategic Survey Program
[4] have demonstrated the unique power of weak gravitational lensing from
cosmic shear and its cross correlation with galaxy clustering in probing the
growth of cosmic structure. Cosmic shear surveys rely on measurements of the
distortion of _shapes_ of galaxies, which are subtly sheared and magnified by
the distribution of dark matter along the line of sight. In the weak lensing
paradigm, we are specifically interested in probing how the quadrupole moment
of the projected distribution of galaxy light (also called the galaxy
ellipticity) is altered, and the effect of cosmic shear is teased out through
correlated measurements of these shapes. The effect of shear is described as
an additive contribution to the projected quadrupole moments of galaxies [5]

$g_{ab}(\hat{n})\approx g^{I}_{ab}(\hat{n})+\gamma_{ab}(\hat{n})+\cdots,$

where $g^{I}$ is the assumed intrinsic ellipticity of the galaxy, $\gamma$ is
the shear, $\hat{n}$ is the direction on the sky of this field, and $a,b$ are
the indices of the matrix which describes the galaxy shape ellipsoid. The
shapes of galaxies themselves are not distributed in an uncorrelated fashion,
but rather correlated with large-scale structure. These _intrinsic alignments_
of galaxies are a key contaminant of the cosmic shear signal [6, 7] and their
contributions to the statistics of shapes must be included in order to infer
cosmological parameters in an unbiased fashion [8]. Schematically, the 2-point
function between two maps of measured galaxy ellipticities will be

$\langle
g_{ab}(\hat{n}_{1})g_{cd}(\hat{n}_{2})\rangle=\langle\gamma\gamma\rangle+\langle
g^{I}\gamma\rangle+\langle\gamma g^{I}\rangle+\langle g^{I}g^{I}\rangle,$
(1.1)

such that the cosmic shear spectrum $\left\langle\gamma\gamma\right\rangle$
receives contamination from cross correlations between cosmic shear and galaxy
shape (‘GI’, ‘IG’) as well as galaxy shape autocorrelations (‘II’). Since
shear probes the matter density field, such correlations probe the cross-
spectrum between shapes and matter

$P^{GI}\sim\langle g_{ab}\delta_{m}\rangle,\quad P^{II}\sim\langle
g_{ab}g_{cd}\rangle.$

Studies of cross-correlations between galaxy positions and galaxy shapes
(galaxy–galaxy lensing and ‘3D intrinsic alignment’ measurements) will instead
probe the ‘gI’ terms

$P^{gI}\sim\langle\delta_{g}g_{ab}\rangle.$ (1.2)

While the contributions from cosmic shear are delocalized in space due to the
broad support from the lensing efficiency kernel, the contributions from
intrinsic alignments depend only on the redshift distribution of the source
galaxies whose shapes we are measuring – their projection kernel is equivalent
to that of clustering. This implies that while the challenges of modelling
cosmic shear are significant (as they probe deeply nonlinear small scales),
models of intrinsic alignments may be formulated purely perturbatively and
enjoy comparable success to perturbation theories of the galaxy density used
to analyze large-scale structure surveys.

Despite being considered purely a contaminant in cosmic shear surveys, in
recent years there has been growing interest in intrinsic alignments as
cosmological signals in their own right. In addition to being helpful in
understanding the size of intrinsic alignments in galaxies [9, 10], such
measurements can also be used along with standard density statistics to
improve cosmological constraints. For example, refs. [11, 12] recently
examined the cross correlation between galaxy shapes and densities in the
Sloan Digital Survey to respectively enhance BAO and redshift-space
distortions measurements. Beyond these, since shapes are tensors they can
potentially probe different primordial physics than densities, for example in
the case of anisotropic primordial non-Gaussianity [13, 14, 15, 16, 17].

There is a long history of perturbative models of galaxy intrinsic alignments
(see e.g. [7] for a review). Earlier works typically assumed that elliptical
galaxies aligned with the local tidal field, giving rise to linear alignment
(LA) while spiral galaxies gained their shapes by having their angular
momentum align with halos, which in tidal torquing (TT) theory led to a
dependence on the _square_ of the tidal field [18]. More recent formulations
of galaxy alignments unified these two theories into the TATT (tidal alignment
and tidal torquing) model [19], which was the reference model for intrinsic
alignments used in the Dark Energy Survey Year 3 analysis [2]. In particular,
ref. [20] showed that, accounting for time dependence, the galaxy shape field
can be expressed as a linear combination of all rotationally invariant
combinations of the Hessian of the gravitational potential, and derived this
expansion to second order in perturbation theory. Recently, ref. [21]
introduced an effective-theory description of galaxy shapes based on Eulerian
perturbation theory. This formulation completes the necessary ingredients for
predicting 2-point correlations of intrinsic alignments up to 1-loop order by
including third-order bias contributions and effective-theory corrections
needed to tame small-scale dependences of the theory, and includes the TATT
model as a subset of possible contributions. This full Eulerian EFT for galaxy
shapes has recently been shown to accurately capture the power spectrum of
projected halo shapes on perturbative scales [22].

All analytic approaches to modelling the summary statistics of galaxy shapes
so far in the literature have relied on Eulerian perturbation theory (EPT),
where the shapes of galaxies are expanded perturbatively in terms of their
present environment. An alternative method to model the large scale properties
of galaxies is Lagrangian perturbation theory (LPT), where e.g. galaxy
densities can be written in terms of the initial conditions at their initial
(Lagrangian) positions, which are then “advected” to their final positions by
the flow of dark matter due to gravity. These two approaches can be shown to
be perturbatively equivalent in the case of galaxy densities, and indeed can
be understood as convenient parametrizations of a more general picture where
galaxies respond to their environments along their trajectories across cosmic
time (see e.g. [23]). Nonetheless, LPT provides a complementary perspective to
structure formation, for example by providing a simple physical picture of the
damping of baryon acoustic oscillations by long-wavelength modes and method to
resum them [24, 25, 26, 27], as well as simplifying the relation between
galaxy bias and dynamical contributions from advection [28]. LPT is also a
natural formalism in which to understand redshift-space distortions, since in
it galaxies simply have to be advected by an “extra” displacement proportional
to their velocities [24, 25, 29]. Finally, a method to augment the predictions
of LPT with exact dark matter displacements from simulations—termed Hybrid
Effective Field Theory—has recently been proposed and shows great promise in
future galaxy-galaxy lensing analyses [30, 31, 32, 33].

The main purpose of this work is to extend Lagrangian perturbation theory into
the realm of galaxy shapes and to derive the associated summary statistics,
inspired by the foundational work of ref. [21] in the Eulerian formalism. We
are not the first authors to consider this extension—while this work was in
the early stages of preparation, refs. [34, 35, 36] appeared which detailed
the formal construction of a Lagrangian theory of shapes within the Lagrangian
integrated perturbation theory (iPT)—however to our knowledge we are the first
to explicitly construct the necessary bias expansion of shapes up to 1-loop in
perturbation theory, including effective-theory corrections, and compute
summary statistics within a resummation scheme that respects Galilean
invariance. The paper is structured as follows: in §2 we formulate the bias
expasion of shapes in terms of Lagrangian displacements, describe their
behavior under advection, and relate the resulting contributions to previous
results in the literature. This expansion also includes counter-terms and a
stochasticity tensor whose properties we derive. In §3 we review analytic
forms and computational techniques for the power spectrum of biased tracers in
LPT and extend them to the case of tensor shape fields. Sections 4, 5 and 6
are the main sections of this work and present calculations of the matter-
shape, density-shape and shape-shape correlations in LPT. Finally, we compare
our theoretical predictions to N-body data in §7 before concluding in §8. All
of the expressions derived in this paper are numerically implemented and made
publicly available in the code
spinosaurus111https://github.com/sfschen/spinosaurus, which builds on top of
the framework of velocileptors [37] to efficiently predict power spectra in
LPT.

###  1.1 Notation

Throughout the text we will refer to both Lagrangian $\boldsymbol{q}$ and
Eulerian $\boldsymbol{x}$ coordinates, related via Lagrangian displacements
$\boldsymbol{x}(\boldsymbol{q},\tau)=\boldsymbol{q}+\Psi(\boldsymbol{q},\tau)$.
Our main results will be the observed power spectra of galaxy shapes and
densities in Eulerian space—we will write these cross spectra between fields
$A(\boldsymbol{x})$ and $B(\boldsymbol{x})$, which may be scalars or index-
carrying tensors, as

$P_{AB}(k)=\left\langle
A(\textbf{k})|B(\textbf{k}^{\prime})\right\rangle^{\prime}.$ (1.3)

The ′ denotes the two-point function with a Dirac delta
$(2\pi)^{3}\delta_{D}(\textbf{k}+\textbf{k}^{\prime})$ removed, and the
Fourier transform is with respect to the Eulerian coordinate $\boldsymbol{x}$.

In addition to the observed two-point functions in Eulerian space we will deal
extensively with Lagrangian correlators of the form

$\langle O_{a}(\boldsymbol{q}_{2})\ O_{b}(\boldsymbol{q}_{1})\rangle=\langle
O_{a}(\boldsymbol{q}=\boldsymbol{q}_{2}-\boldsymbol{q}_{1})\
O_{b}(\textbf{0})\rangle=\text{iFT}\\{\langle
O_{a}(\textbf{k})|O_{b}(\textbf{k}^{\prime})\rangle_{\boldsymbol{q}}^{\prime}\\}$
(1.4)

where “iFT” refers to the inverse Fourier transform. The subscript
$\boldsymbol{q}$ refers implies that the Fourier transforms are with respect
to $\boldsymbol{q}$.

We will make extensive use of the Legendre basis $P_{\ell}(\hat{n})$ of
tensors, which Fourier transforms into itself:

$\displaystyle\int\frac{d^{3}\textbf{k}}{(2\pi)^{3}}\
e^{i\textbf{k}\cdot\boldsymbol{q}}A(k)P_{\ell}(\hat{k})=i^{\ell}P_{\ell}(\hat{q})\int\frac{dk\
k^{2}}{2\pi^{2}}A(k)\ j_{\ell}(kq)$ $\displaystyle\int d^{3}\boldsymbol{q}\
e^{-i\textbf{k}\cdot\boldsymbol{q}}A(q)P_{\ell}(\hat{q})=4\pi(-i)^{\ell}P_{\ell}(\hat{k})\
\int dq\ q^{2}\ A(q)j_{\ell}(kq).$ (1.5)

These basis tensors are the unique traceless symmetric tensors constructed
from e.g., $\hat{k}$, such that for any unit vector $\hat{n}$ we have that
$P_{\ell}(\hat{k}\cdot\hat{n})=\hat{n}_{i_{1}}...\hat{n}_{i_{\ell}}P_{\ell,i_{1}\cdots
i_{\ell}}(\hat{k})$, from which the above relations can be proven [38]. We
list them in App. A. The Hankel transforms in Equation 1.5 can be expressed in
terms of generalized correlation functions [39]

$\xi^{\ell}_{n}[F](q)=\int_{0}^{\infty}\frac{dk}{2\pi^{2}}k^{2+n}j_{\ell}(kq)F(k).$
(1.6)

We will often use $\xi^{\ell}_{n}$ without the square brackets to refer
specifically to Hankel transforms of the linear power spectrum $P_{\rm
lin}(k)$. As discussed in [39], the decomposition of perturbation theory
integrals into these particular generalized correlation functions allows for
their efficient evaluation as 1D Fast Fourier Transforms using the FFTLog
algorithm.

The angular structure of tensor correlators can be significantly simplified by
adopting the helicity basis of ref. [21]. This helicity basis is composed of
five $\ell=2$ tensor spherical harmonics
$Y_{\ell=2,ab}^{(m)}(\boldsymbol{\hat{k}})$ with helicities $m=0,\pm 1,\pm 2$.
These are constructed from tensor products of the Fourier wave-vector
$\boldsymbol{\hat{k}}$ and two orthogonal raising/lowering vectors
$\boldsymbol{\hat{e}^{\pm}}$

$Y^{(0)}_{2,ab}=\sqrt{\frac{3}{2}}\left(\hat{k}_{a}\hat{k}_{b}-\frac{1}{3}\delta_{ab}\right),\quad
Y^{(\pm
1)}_{2,ab}=\sqrt{\frac{1}{2}}\left(\hat{k}_{a}\hat{e}^{\pm}_{b}+\hat{k}_{b}\hat{e}^{\pm}_{a}\right),\quad
Y^{(\pm 2)}_{2,ab}=\hat{e}^{\pm}_{a}\hat{e}^{\pm}_{b}.$ (1.7)

We refer the reader to ref. [21] for a detailed derivation of this basis and
its properties.

##  2 The Lagrangian Theory of Shapes

Within the Lagrangian picture, the galaxy density is modeled by first
constructing the initial overdensity of galaxy progenitors in Lagrangian space
at their initial positions $\boldsymbol{q}$, $\delta_{g}(\boldsymbol{q})$,
then advecting these galaxies to their present day positions
$\boldsymbol{x}=\boldsymbol{q}+\Psi(\boldsymbol{q},\tau)$, where the
displacement $\Psi(\boldsymbol{q},\tau)$ defines the trajectory of a fluid
element initially at $\boldsymbol{q}$. The observed galaxy density is then
given by number conservation to be [28]

$1+\delta_{g}(\boldsymbol{x},\tau)=\int d^{3}\boldsymbol{q}\
(1+\delta_{g}(\boldsymbol{q}))\
\delta_{D}(\boldsymbol{x}-\boldsymbol{q}-\Psi(\boldsymbol{q},\tau))=(1+\delta_{m}(\boldsymbol{x},\tau))(1+\delta_{g}(\boldsymbol{q})).$
(2.1)

The final equality holds in the single-stream regime and derives from the fact
that the Jacobian of the transformation between $\boldsymbol{x}$ and
$\boldsymbol{q}$ is given by the matter density. Within this formalism then
there is a separation of galaxy clustering due to biasing, encoded in the
initial Lagrangian density $\delta_{g}(\boldsymbol{q})$, and the clustering
due to the underlying dynamics, quantified by $\Psi$.

We can similarly model the shape density of a population of galaxies or halos
by considering first the shapes of their progenitors. For example, if there
are a collection of objects with number density $\delta_{g}(\boldsymbol{q})$
with moment-of-inertia tensor $I_{ab}(\boldsymbol{q})$ we have that their
shape density is given by

$M_{ab}(\boldsymbol{q})=(1+\delta_{g}(\boldsymbol{q}))I_{ab}(\boldsymbol{q}),\quad
I_{ab}(\boldsymbol{q})=\int d^{3}\textbf{s}\ \rho(\textbf{s},\boldsymbol{q})\
\textbf{s}_{a}\textbf{s}_{b}$ (2.2)

where $\rho(\textbf{s},\boldsymbol{q})$ is the three-dimensional profile of
the object centered at $\boldsymbol{q}$. As structure forms, each volume
element at $\boldsymbol{q}$ is advected to
$\boldsymbol{q}+\Psi(\boldsymbol{q})$, rearranging the distribution of fluid
elements withn the object and this, along with number conservation, gives the
advected shape density

$M_{ab}(\boldsymbol{x})=\int d^{3}\boldsymbol{q}\
\delta_{D}(\boldsymbol{x}-\boldsymbol{q}-\Psi(\boldsymbol{q}))\
R_{ac}(\boldsymbol{q})R_{db}(\boldsymbol{q})M_{cd}(\boldsymbol{q}),\
R_{ab}(\boldsymbol{q})=\delta_{ab}+\nabla_{a}\Psi_{b}.$ (2.3)

We see therefore that the shape density is a symmetric rank two tensor and
transforms as such under the transformation
$\boldsymbol{q}\rightarrow\boldsymbol{q}+\Psi(\boldsymbol{q})$, with again an
additional multiplication by the matter density at that point. This factor,
and the advection integral, implies we always deal with the density-weighted
shape in the Lagrangian formalism, as opposed to the shape on its own—it is
important to note, however, that the former and not the latter is indeed the
physical quantity measured in cosmological surveys, since the spin of an empty
region of space is undefined. We will refer to the integral above with and
without the factors of $R_{ab}(\boldsymbol{q})$ as active and passive
advection; the latter amounts to the parallel transport of galaxy shapes as
they travel along trajectories $\boldsymbol{x}(\boldsymbol{q},\tau)$. As we
will see, active advection induces contributions to the shapes of galaxies
that are entirely degenerate with the Lagrangian bias expansion of
$M_{ab}(\boldsymbol{q})$ itself, and as such we will not explicitly include
the factors of $R_{ab}$ in most of what follows but rather absorb all
contributions to the shape density outside of the $\delta$ function into
$\textbf{M}(\boldsymbol{q})$.

It will be useful in the discussion below to decompose these tensors into a
scalar trace component and a trace-free component

$M_{ab}:\
M=\Tr{M_{ab}},\text{TF}\left\\{M_{ab}\right\\}=M_{ab}-\frac{1}{3}\delta_{ab}M.$
(2.4)

such that the tensor can also be written as $M_{ab}=\left\langle
M\right\rangle((1+\delta_{M})\ \delta_{ab}/3+\gamma_{ab})$, with
$\gamma_{ab}=\text{TF}\left\\{M\right\\}_{ab}/\left\langle M\right\rangle$.
The former, $\delta_{M}$, characterizes the relative fluctuation in galaxy
size density and is a scalar that follows the familiar bias expansion for e.g.
galaxy densities.

###  2.1 Bias Expansion

The Lagrangian shape field $\textbf{M}(\boldsymbol{q})$ is a function of the
gravitational environment around which a halo or galaxy forms along its
trajectory $\boldsymbol{x}(t)=\boldsymbol{q}+\Psi(\boldsymbol{q},t)$.
Specifically, it can be written in terms of locally observable quantities such
as the tidal field $\nabla_{i}\nabla_{j}\Phi$ and velocity gradients
$\nabla_{i}v_{j}$—these can in turn all be written in terms of the Lagrangian
shear tensor $L_{ij}=\partial_{i}\Psi_{j}$. Significantly, even for initially
isotropic protohalos these types of shape terms are generated dynamically,
since $R_{ij}=\delta_{ij}+L_{ij}$. This makes sense since gravitational
evolution produces sheets and filaments from close-to homogeneous initial
conditions even in the Zel’dovich approximation. Since galaxy formation is
nonlocal in time—that is, it occurs over the same $\sim H_{0}^{-1}$ time scale
as structure formation itself, it is important to allow contributions at each
order $n$ in perturbation theory to have independent coefficients, since their
relative effect is a function of which epochs each sample of galaxies is
sensitive to. For example, active advection generates new contributions to the
bias expansion by multiplying the Lagrangian shape $M_{ab}$ by $R_{ab}$;
however, the coefficients generated in this way should not be taken as fixed
since their precise values will depend on, e.g. the time at which a galaxy
decouples from large-scale structure following gravitational collapse: the
earlier this happens the smaller the effect of higher-order terms. The
takeaway, rather is that the dynamical advection of galaxies produces
contributions to the shapes of galaxies which are already included in the bias
expansion based on $L_{ij}$ evaluated at each order, much like the case for
galaxy densities where the relevant expansion is based on the scalar functions
of $L_{ij}$ [23].

From the above discussion we can see that a complete bias expansion for galaxy
shapes can be built out of the Lagrangian shear tensor of order $n$, denoted
by $L^{(n)}_{ij}$, along with other rank-2 functions of it obeying the
underlying symmetries. The complete set of such contributions up to third
order in perturbation theory is given by:

$\displaystyle 1^{\rm st}\text{ order: }$ $\displaystyle L^{(1)}_{ij}$
$\displaystyle 2^{\rm nd}\text{ order: }$ $\displaystyle L^{(2)}_{ij},\
(\textbf{L}^{(1)}\textbf{L}^{(1)})_{ij},\ \Tr{\textbf{L}}L^{(1)}_{ij}$
$\displaystyle 3^{\rm rd}\text{ order: }$ $\displaystyle L^{(3)}_{ij},\
(\textbf{L}^{(2)}\textbf{L}^{(1)})_{ij},\ \Tr{\textbf{L}^{(1)}}L^{(2)}_{ij},$
$\displaystyle(\textbf{L}^{(1)}\textbf{L}^{(1)}\textbf{L}^{(1)})_{ij},\
\Tr{\textbf{L}^{(1)}\textbf{L}^{(1)}}L^{(1)}_{ij},\Tr{\textbf{L}^{(1)}}^{2}L^{(1)}_{ij},\
\Tr{\textbf{L}^{(1)}}(\textbf{L}^{(1)}\textbf{L}^{(1)})_{ij}.$

Since the shape tensor is symmetric we implictly refer to the symmetrized
versions of these tensors. In the above we have used that
$\Tr{\textbf{L}^{(n)}}$ can be written in terms of lower-order operators in
order to omit certain linearly dependent contributions—note that this
degeneracy does not exist at the full tensor level [23, 21]. In addition, we
can construct operators proportional to the Kronecker delta, given by the
trace components of the above: $\Tr{O_{ab}}\delta_{ab}$. Up to the Kronecker
delta these terms are equivalent to the scalar bias expansion for galaxy
densities, yielding one, two and four bias parameters at first, second and
third order, respectively [23].

The expansion above is a general bias expansion for galaxy shapes, not just
the traceless part—in order to better distinguish these we can separate an
operator into its trace and trace-free parts:

$O_{ab}(\boldsymbol{q}):\
\Tr{O_{ab}(\boldsymbol{q})},\text{TF}\left\\{O_{ab}(\boldsymbol{q})\right\\}.$
(2.5)

The trace component is described by the Kronecker delta terms in bias
expansion, so performing this split amounts to a linear redefinition of the
bias basis. Since we are primarily interested in galaxy shapes in this paper
we will implicitly refer only to the trace-free component, and drop the “TF”
in most cases, unless explicitly noted otherwise.

In order to make better contact with the literature we note that it is
customary to decompose the linear Lagrangian shear as its trace and trace-free
components

$L^{(1)}_{ij}(\boldsymbol{q})=-\frac{1}{3}\delta(\boldsymbol{q})\delta_{ij}-s_{ij}(\boldsymbol{q})$
(2.6)

as well as to introduce the quadratic operator

$t_{ij}=\Big{(}\frac{\nabla_{i}\nabla_{j}}{\nabla^{2}}-\frac{1}{3}\delta_{ij}\Big{)}(\theta^{(2)}-\delta^{(2)})=\frac{4}{3}\text{TF}\left\\{\textbf{L}^{(2)}\right\\}_{ij}$

where $\theta$ is the velocity potential. In terms of these operators we can
rewrite the above as

$\displaystyle 1^{\rm st}\text{ order: }$ $\displaystyle s_{ij}$
$\displaystyle 2^{\rm nd}\text{ order: }$ $\displaystyle t_{ij},\ s^{2}_{ij},\
\delta s_{ij}$ $\displaystyle 3^{\rm rd}\text{ order: }$ $\displaystyle
L^{(3)}_{ij},\ (st)_{ij},\ \delta t_{ij},\ s^{3}_{ij},\
s^{2}s_{ij},\delta^{2}s_{ij},\ \delta s^{2}_{ij}$

We note that including terms up to second order reproduces the TATT model
[19], and there are in principle seven new operators that contribute at third
order.

###  2.2 Cubic Operators in the 1-loop Power Spectrum

At the one-loop order we consider in this paper, we are interested exclusively
in the 2-point correlations of Lagrangian bias operators $O(\boldsymbol{q})$
and displacements $\Psi(\boldsymbol{q})$ up to second order in the linear
power spectrum. In this regime contributions due to cubic operators can come
only from cross correlation with the linear Lagrangian displacement or
density, with one internal loop

$O_{ab}(\textbf{k})\rightarrow\left(3\int_{\textbf{p}}K^{O}_{ab}(\textbf{k},\textbf{p},-\textbf{p})P_{L}(p)\right)\delta(\textbf{k})=\left(\frac{k_{a}k_{b}}{k^{2}}-\frac{1}{3}\delta_{ab}\right)\
A_{O}(k)\ \delta(\textbf{k})$ (2.8)

where $K^{O}_{ab}$ is the symmetrized kernel for $O$ and $A_{O}(k)$ is a
scalar function of $k$ and form following the last equality follows from
symmetry. Note that the Fourier transform here is in Lagrangian space. For
operators that are built from powers of $\\{\delta,s_{ij}\\}$ these
contributions are entirely degenerate with the linear theory contributions,
i.e. $A_{O}(k)$ is a constant proportional to $\sigma^{2}$ for these
operators, since we can use Wick’s theorem to show e.g.

$\delta^{2}s_{ab}(\boldsymbol{q})\rightarrow\sigma^{2}s_{ab}(\boldsymbol{q})$
(2.9)

where the only nonzero contraction is between the two $\delta$’s due to
$s_{ab}$ being traceless.333The other Zel’dovich cubic operators give
$s^{2}s_{ab}\rightarrow\frac{14}{15}\sigma^{2}s_{ab},\
\delta^{2}s_{ab}\rightarrow\sigma^{2}s_{ab},\ \delta s^{2}_{ab}\rightarrow 0,\
s^{3}_{ab}\rightarrow\frac{7}{15}\sigma^{2}s_{ab}$ (2.10) where the right-hand
side is the equivalent linear operator that gives the (13) correlation. For
the remaining ‘non-Zel’dovich operators’ their contributions must be
equivalent to that of the linear Lagrangian shear $s_{ab}(\boldsymbol{q})$
with the input power spectrum substituted for $A_{O}(k)P_{\rm lin}(k)$.

In order to compute $A_{O}(k)$ for the three cubic non-Zel’dovich operators we
can compute their cross correlation with $\delta$ in Lagrangian space. Of
these we have directly from the correlator of $\Psi^{(3)}$ with $\delta$ [28]

$A_{L^{(3)}}(k)=\frac{5}{21}R_{1}(k),$ (2.11)

where $R_{1}$ is an integral which is given by the mode-coupling structure of
this (13) correlator, defined in ref. [24]. It is part of a broader family of
integrals $R_{n}$ and $Q_{n}$ which will be referenced throughout this text
and whose expressions are found in refs. [24, 28]. The other two cubic
correlators are given by the projections of the tensor (13) integrals

$\displaystyle I^{\delta
t}_{ab}(\textbf{k})=-\frac{4}{7}P(k)\int_{\textbf{p}}\left(\frac{(\textbf{k}-\textbf{p})_{a}(\textbf{k}-\textbf{p})_{b}}{(\textbf{k}-\textbf{p})^{2}}-\frac{1}{3}\delta_{ab}\right)\left(1-\frac{(\textbf{k}\cdot\textbf{p})^{2}}{k^{2}p^{2}}\right)P(p)$
$\displaystyle
I^{st}_{ab}(\textbf{k})=-\frac{4}{7}P(k)\int_{\textbf{p}}\left(\frac{\textbf{p}_{a}\textbf{p}_{c}}{p^{2}}-\frac{\delta_{ac}}{3}\right)\left(\frac{(\textbf{k}-\textbf{p})_{c}(\textbf{k}-\textbf{p})_{b}}{(\textbf{k}-\textbf{p})^{2}}-\frac{1}{3}\delta_{cb}\right)\left(1-\frac{(\textbf{k}\cdot\textbf{p})^{2}}{k^{2}p^{2}}\right)P(p).$

onto the quadrupole tensor $P_{2}^{ab}(\hat{k})$ (Eqn. A.3). Decomposing these
projected integrals into generalized correlation functions as in ref. [39] we
find

$\displaystyle A_{\delta t}(k)=-\frac{8}{735}\int$ $\displaystyle
dqq\,\big{[}42k^{2}\xi^{0}_{0}(q)j_{0}(kq)-56k\xi^{1}_{1}(q)j_{1}(kq)+5(-7k^{2}\xi^{2}_{0}(q)+5\xi^{2}_{2}(q))j_{2}(kq)$
$\displaystyle+42k\xi^{3}_{1}(q)j_{3}(kq)-18\xi^{4}_{2}(q)j_{4}(kq)\big{]}+\frac{8}{105}\sigma^{2}$
$\displaystyle A_{st}(k)=-\frac{8}{2205}\int$ $\displaystyle
dqq\,\big{[}-38k\xi^{1}_{1}(k)j_{1}(kq)+25(k^{2}\xi^{2}_{0}(q)+\xi^{2}_{2})j_{2}(kq)+14k\xi^{3}_{1}(q)j_{3}(kq)$
$\displaystyle-18(k^{2}\xi^{4}_{0}(q)+\xi^{4}_{2})j_{4}(kq)+10k\xi^{5}_{1}(q)j_{5}(kq)\big{]}+\frac{8}{315}\sigma^{2}.$

In the above we have explicitly included the zero-lag piece proportional to
$\sigma^{2}$. In general these contributions can be removed by appropriate
defining a “normal-ordered” operator, equivalent to a linear combination of
the original operator with lower-order ones, with

$:O^{(3)}_{ab}:=O^{(3)}_{ab}-A_{O}(k\rightarrow 0)\ s_{ab},$

which we will implicitly assume throughout the rest of this work.

Finally, as ref. [21] first discussed, various degeneracies between the bias
parameters exist at one-loop level for traceless galaxy shapes; some algebra
shows that

$A_{st}-\frac{1}{6}A_{\delta
t}-\frac{3}{5}A_{L^{(3)}}=\frac{1}{7}\int_{\textbf{p}}\left(\frac{11}{45}-\frac{4}{3}\mu^{2}+\mu^{4}\right)P(p)=0.$

This implies we only need to introduce two bias parameters at this order;
since $A_{L^{(3)}}$ is shared with the density bias case we will adopt it and
$\delta t_{ab}$ as our conventional choice. This is slightly different from
the case of cubic operators in the perturbation theory of scalar fields, where
only one non-degenerate transfer function $\sigma^{2}(k)$ and effective cubic
parameter ‘$\tilde{b}_{3}$’ are introduced [40].

###  2.3 Effective Theory Considerations

In order to complete the list of perturbative contributions to galaxy shapes
we need to include dimensionful coefficients reflecting that our perturbation
theory is an effective one. The first type are the so-called counterterms or
derivative biases, which in the bias expansion reflect the nonlocality of
astrophysical processes, i.e. that the formation and shape of a galaxy depends
not on its environment not at a point but in a neighborhood around it with
size given roughly by its Lagrangian radius $R_{h}$. This gives rises to
contributions proportional to higher derivatives such as

$\delta_{g}(\boldsymbol{q})\ni\alpha_{g}\nabla^{2}\delta(\boldsymbol{q}),\quad
M_{ab}(\boldsymbol{q})\ni\alpha_{M}\nabla^{2}s_{ab}(\boldsymbol{q}),\quad\alpha_{g,M}\sim\mathcal{O}(R_{h}^{2}).$
(2.12)

Similar considerations require additional contributions to the dynamics, i.e.
displacements [27]. At 1-loop order, where we will work, all of these
contributions lead to corrections to the power spectrum $k^{2}P(k)$, with a
tensor structure given either by the scalar $\delta$ (density bias and
dynamics) or trace-free tensor $s_{ab}$ (shape bias).

Beyond the counterterms we also need to account for stochastic contributions,
which quantify the effect of short wavelength physics which we do not
explicitly model. We will be most concerned with contributions uncorrelated
with long wavelengths,

$\delta_{g}\ni\epsilon(\boldsymbol{q}),\quad
M_{ab}\ni\epsilon_{ab}(\boldsymbol{q}).$ (2.13)

Since these encapsulate small-scale physics their correlation will be
localized to $|\boldsymbol{q}_{1}-\boldsymbol{q}_{2}|\lesssim R_{h}$, where
$R_{h}$ is the halo radius. This localization highly restricts the form these
stochastic terms can contribute to the power spectrum—we will return to this
topic in §3.4.444Beyond these contributions the stochastic contribution to
active advection can also help us explain the smallness of the intrinsic
alignment signal—whereas galaxy bias parameters tend to be order one, galaxy
shape biases tend to be of order a few percent when they are detected at all.
Let small-scale astrophysics contribute some additional rotational “kick” to
the active advection so that the shape goes to  $M(1+\gamma^{\rm
PT})\rightarrow\textbf{R}_{\rm stoch}^{T}(\boldsymbol{x})M(1+\gamma^{\rm
PT})\textbf{R}_{\rm stoch}(\boldsymbol{x})$ (2.14) where $M$ and $\gamma$ are
the perturbative size and shape of the galaxy. Since the stochastic terms do
not correlate with the large-scale structure by construction we can compute
their average ahead of time, i.e.  $M(1+\gamma^{\rm PT})_{ij}\rightarrow
M\langle\textbf{R}_{ia,\rm stoch}\textbf{R}_{aj,\rm
stoch}\rangle+M\langle\textbf{R}_{ia,\rm stoch}\textbf{R}_{bj,\rm
stoch}\rangle\gamma^{\rm PT}_{ab}$ (2.15) It is important to note that the
trace and trace-free parts do not transform equivalently under coordinate
transformations, so the change induced by this additional rotation cannot in
general be accounted for by the normalization $M$. Let us take as an example
the case where R imposes a random rotation about a direction $\vec{\omega}$.
In this case the trace part will not be affected, but it is straightforward to
show that the, averaging over orientations, we have
$\gamma\rightarrow\frac{1}{5}\big{(}1+2\cos\omega+2\cos(2\omega)\big{)}\
\gamma$ (2.16) which is to say that in general, the larger the rotation
magnitude, the more washed out the IA signal will tend to be. This is in
contrast to scalar quantities such as the galaxy density—in that case, adding
noise does not change the amplitude of the bias. In general, to lowest order
we can express the mean of the product of $\textbf{R}_{\rm stoch}$ with itself
as a series of Kronecker $\delta$’s , such that under active advection by
short-wavelength dynamics the large-scale intrinsic alignment amplitude gets
multiplicatively reduced.

##  3 Power Spectra in Lagrangian Perturbation Theory

The expressions we have written down so far have been in terms of the full
three-dimensional components of galaxy shape tensors $M_{ab}$. In order to
evaluate the $n$-point correlation functions of these tensors, and their cross
correlations with galaxy density, however, it will be more convenient to
project the tensors onto an irreducible basis under rotations, and to compute
the spectra of these independent components. Such a spherical harmonic basis
was proposed in ref. [21], and we will adopt it in this paper; in particular a
traceless symmmetric tensor like the galaxy shape can be decomposed as a sum
of helicity modes $m$ of $\ell=2$ tensor spherical harmonics

$M_{ab}(\textbf{k})=\sum_{m=-2}^{2}M_{m}(\textbf{k})Y_{2,ab}^{m}(\hat{k}).$

Enforcing rotational symmetry about k requires that the only nonzero
correlations are between components of the same helicity or, as a special
case, the helicity-zero component with a density.

Within this basis it is easy to see that there is only one nonzero component
of the galaxy density-shape power spectrum, which can be written as a
quadrupole operator and a scalar spectrum

$\left\langle\delta_{g}(\textbf{k})|M_{ab}(\textbf{k}^{\prime})\right\rangle^{\prime}=\frac{3}{2}\left(\hat{k}_{a}\hat{k}_{b}-\frac{1}{3}\delta_{ab}\right)P_{gI}(k).$
(3.1)

For the shape-shape autocorrelations, there will tend to be three components
due to the fact that the three helicity modes do not correlate [21]555This
projection is motivated by properties of the spherical tensor basis derived in
ref. [21]. The general power spectrum of two shape fields can be written as
$\left\langle
M_{ab}(\textbf{k})M_{cd}(\textbf{k}^{\prime})\right\rangle=\sum_{m=-2}^{2}Y_{2,ab}^{m}(\hat{k})Y_{2,cd}^{m}(-\hat{k})P_{2m}(k)=\sum_{m=-2}^{2}(-1)^{m}Y_{2,ab}^{m}(\hat{k})Y_{2,cd}^{-m}(\hat{k})P_{2m}(k).$
(3.2) In order to extract the $P_{2m}$ we can dot each spherical tensor with
its complex conjugate, given by
$Y^{m\ast}_{2}(\hat{k})=(-1)^{m}Y_{2}^{-m}(\hat{k})$, where the two factors of
$(-1)^{m}$ from the two conjugated $Y_{2}^{m}$’s cancel. Noting that
$P_{2m}=P_{2(-m)}$ due to parity then yields the symmetried projection and
factor of $(-1)^{m}/2$.

$P_{2m}(k)=(-1)^{m}\frac{1}{2}\left(Y_{2}^{m}(\hat{k})Y_{2}^{-m}(\hat{k})+Y_{2}^{-m}(\hat{k})Y_{2}^{m}(\hat{k})\right)_{abcd}\left\langle
M_{ab}(\textbf{k})|M_{cd}(\textbf{k}^{\prime})\right\rangle^{\prime}.$ (3.3)

This decomposition implies the angular structure of shape–shape correlators is
simpler than naively anticipated given the number of tensor components – and a
reduced number of spectra need to be computed in order to describe the two-
point statistics of shapes and their cross-correlations with density. For the
rest of this paper we will compute the quadrupole component density–shape
correlations as in 3.1 and the three components of the shape–shape
correlations as shown above.

###  3.1 The Generating Function

We now describe the computing of power spectra of biased tracers in Lagrangian
Perturbation Theory. From Equation 2.1, the power spectrum between two fields
$\alpha,\beta$ in LPT is given by [28, 25]

$P(\textbf{k})=\int d^{3}\boldsymbol{q}\
e^{i\textbf{k}\cdot\boldsymbol{q}}\langle
e^{i\textbf{k}\cdot\Delta(\boldsymbol{q}_{1},\boldsymbol{q}_{2})}F_{\alpha}(\boldsymbol{q}_{1})F_{\beta}(\boldsymbol{q}_{2})\rangle_{\boldsymbol{q}=\boldsymbol{q}_{2}-\boldsymbol{q}_{1}}$
(3.4)

where $\Delta=\Psi(\boldsymbol{q}_{2})-\Psi(\boldsymbol{q}_{1})$. The
functions $F_{\alpha,\beta}(\boldsymbol{q})$ denote the unadvected values of
these fields, i.e. the proto-halo density $1+\delta_{g}(\boldsymbol{q})$ for
densities and the progenitor shape distribution
$(1+\delta_{g}(\boldsymbol{q}))I_{ab}(\boldsymbol{q})$ for shapes. We have
dropped a delta function at $\textbf{k}=0$ that would appear in the case of
density-density statistics and suppressed tensor indices that would appear in
the case of spin statistics.

The $F_{\alpha,\beta}$ are functions of the local environment that can be
perturbatively expanded. In order to derive expressions for the two-point
function of scalar or tensor biased tracers in LPT we make use of the
generating function666An alternative treatment would be to treat the sources
$\alpha_{i}$ as $\boldsymbol{q}$-dependent functions – that is the field and
coefficients would be integrated as $\int d^{3}\boldsymbol{q}\
\alpha_{i}(\boldsymbol{q})O_{i}(\boldsymbol{q})$ and functional derivatives
taken to isolate particular bias operators. These two approaches are
equivalent, but the functional derivative approach is more natural when
dealing with higher-order statistics where many more $\alpha_{n}$ would be
involved.

$M(\boldsymbol{\alpha})=\left\langle\exp{\big{[}i\textbf{k}\cdot\Delta+\sum_{a}\alpha_{1}^{a}O^{a}(\boldsymbol{q}_{1})+\alpha_{2}^{a}O^{a}(\boldsymbol{q}_{2})\big{]}}\right\rangle.$
(3.5)

which can be evaluated via the cumulant expansion. For example, the
$\boldsymbol{\alpha}=0$ piece is given by

$M(\mathbf{0})=\exp\left[-\frac{1}{2}k_{i}k_{j}A_{ij}(\boldsymbol{q})-i\frac{1}{6}k_{i}k_{j}k_{k}W_{ijk}(\boldsymbol{q})+\cdots\right],\quad
A_{ij}(\boldsymbol{q})=\langle\Delta_{i}\Delta_{j}\rangle,W_{ijk}(\boldsymbol{q})=\langle\Delta_{i}\Delta_{j}\Delta_{k}\rangle.$

Ignoring the self-pairs we have

$\displaystyle\frac{M(\boldsymbol{\alpha})}{M(\mathbf{0})}=$
$\displaystyle\exp{\Big{[}\sum_{a}(\alpha_{1}^{a}+\alpha_{2}^{a})ik_{i}\langle\Delta_{i}O_{1}^{a}\rangle+\sum_{a,b}\alpha_{1}^{a}\alpha_{2}^{b}\langle
O_{1}^{a}O_{2}^{b}\rangle\Big{]}}$
$\displaystyle\exp{\Big{[}-\frac{1}{2}\sum_{a}(\alpha_{1}^{a}+\alpha_{2}^{a})k_{i}k_{j}\langle\Delta_{i}\Delta_{j}O_{1}^{a}\rangle_{c}+\sum_{a,b}i\alpha_{1}^{a}\alpha_{2}^{b}k_{i}\langle\Delta_{i}O_{1}^{a}O_{2}^{b}\rangle_{c}+...\
\Big{]}}$ (3.6)

where we have separated the exponentials by cumulant, truncating at cubic
order. We have also dropped any terms corresponding to more than one $\alpha$
from any given point, or that is more than quadratic in $\alpha$, since these
will not appear in the 2-point function. The correlators in the above equation
are invariant under inversion $\boldsymbol{q}\rightarrow-\boldsymbol{q}$,
which takes $\Psi(\boldsymbol{q})\rightarrow-\Psi(-\boldsymbol{q})$ and thus,
along with a translation, gives $\Delta\rightarrow\Delta$, $O_{1}\rightarrow
O_{2}$.

Given the generating function above, matter–tracer and tracer–tracer
correlation functions are computed by taking derivatives with respect to
$\alpha_{i}$ and then setting these ‘sources’ to zero. Carrying out these
derivatives the power spectra in question are given by

$P_{mt}(\textbf{k})=\sum_{a}b_{a}\int d^{3}\boldsymbol{q}\
e^{i\textbf{k}\cdot\boldsymbol{q}}\left.\frac{\partial
M}{\partial\alpha_{2}^{a}}\right|_{\boldsymbol{\alpha}=0},\quad
P_{tt}=\sum_{a,b}b_{a}b_{b}\int d^{3}\boldsymbol{q}\
e^{i\textbf{k}\cdot\boldsymbol{q}}\left.\frac{\partial^{2}M}{\partial\alpha_{1}^{a}\partial\alpha_{2}^{b}}\right|_{\boldsymbol{\alpha}=0}$
(3.7)

where we have

$\left.\frac{\partial
M}{\partial\alpha_{2}^{a}}\right|_{\boldsymbol{\alpha}=0}=M(\mathbf{0})\Big{(}ik_{i}\langle\Delta_{i}O_{2}^{a}\rangle-\frac{1}{2}k_{i}k_{j}\langle\Delta_{i}\Delta_{j}O_{2}^{a}\rangle+...\Big{)}$
$\mathclap{\left.\frac{\partial^{2}M}{\partial\alpha_{1}^{a}\partial\alpha_{2}^{b}}\right|_{\boldsymbol{\alpha}=0}\quad\quad}=M(\mathbf{0})\Big{(}-k_{i}k_{j}\langle\Delta_{i}O_{1}^{a}\rangle\langle\Delta_{j}O_{2}^{b}\rangle+\langle
O_{1}^{a}O_{2}^{b}\rangle+ik_{i}\langle\Delta_{i}O_{1}^{a}O_{2}^{b}\rangle_{c}+...\Big{)}.$
(3.8)

These expressions are ‘master formulae’ for the power spectra of biased
tracers in LPT. Each of the correlators $\langle\Delta_{i}O\rangle$,
$\langle\Delta_{i}\Delta_{j}O\rangle$ and $\langle O_{1}O_{2}\rangle$ are then
evaluated in pertubation theory up to fourth–order in fields. The complete
order-by-order structure of these correlators which contribute to biased
tracer spectra are enumerated in Appendix B.

###  3.2 IR Resummation and the Eulerian Limit

One of the key features of Lagrangian perturbation theory is that it allows
for a straightforward treatment of long-wavelength linear displacements, which
lead to large nonlinearities such as the damping of baryon acoustic
oscillations (BAO), and hence must be resummed. Within the generating
function, these effects are captured by the linear contribution to
$M(\boldsymbol{\alpha})$, i.e. $A^{\rm
lin}_{ij}(\boldsymbol{q})=\langle\Delta^{(1)}_{i}\Delta^{(1)}_{j}\rangle$. In
particular, it is desirable to split the linear displacements above and below
some exponential cutoff $S(k)=e^{-(k/k_{\rm IR})^{2}/2}$ and resum the long-
wavelength piece while expanding the short-wavelength piece to 1-loop order
such that

$M^{\rm
1loop}(\boldsymbol{0})=\exp\left[-\frac{1}{2}k_{i}k_{j}A^{<}_{ij}(\boldsymbol{q})\right]\left(1-\frac{1}{2}k_{i}k_{j}(A_{ij}^{>}+A_{ij}^{\rm
loop})-\frac{i}{6}k_{i}k_{j}k_{k}W_{ijk}+\frac{1}{8}k_{i}k_{j}k_{k}k_{l}A_{ij}^{>}A_{kl}^{>}\right)$

and similarly

$M^{\rm
1loop}(\boldsymbol{\alpha})=\exp\left[-\frac{1}{2}k_{i}k_{j}A^{<}_{ij}(\boldsymbol{q})\right]\left(1+M^{\rm
linear}(\boldsymbol{\alpha})+M^{\rm loop}(\boldsymbol{\alpha})\right)$ (3.9)

where the linear and loop terms in parentheses have their contributions from
long displacements resummed in the exponential. From the above expressions it
is also possible to extract the un-IR-resummed Eulerian predictions for shape
statistics; these are simply equivalent to setting $k_{\rm IR}=0$, i.e. fully
expanding the exponentiating exponents to 1-loop order [37]. We can also
extract a different resummation based on splitting the linear power spectrum
into its “wiggle” and “no-wiggle” components by noting that any Lagrangian
correlators with the a peak at the BAO scale will pick out the exponentiated
piece in Equation 3.9 in the $\boldsymbol{q}$ integral in Equation 3.4 [41].
Performing this saddle-point approximation yields the IR-resummed form for
shapes derived in Equation 4.24 in ref. [21].

###  3.3 Angular Structure of Lagrangian Correlators

Having laid out the general structure of tensor power spectra in Lagrangian
perturbation theory, much of the remaining work of this paper will be in
calculating the correlators between bias operators and pairwise displacements
that compose the generating function in Equation 3.6. Before we proceed it
will thus be helpful to describe the general structure of these correlators.
We are interested in the Lagrangian space correlators of the form

$\left\langle\Delta_{i_{1}}...\Delta_{i_{n}}O_{ab}(\boldsymbol{q}_{1})O_{cd}(\boldsymbol{q}_{2})\right\rangle,\quad\left\langle\Delta_{i_{1}}...\Delta_{i_{n}}O(\boldsymbol{q}_{1})O_{ab}(\boldsymbol{q}_{2})\right\rangle$
(3.10)

where, in the case of cross correlations with matter we have $O=1$. By
symmetry, any such correlators must be expressible as functions of
$q=|\boldsymbol{q}|$ multiplied by products of Kronecker $\delta$’s and the
unit vector $\hat{q}$. However, since perturbation-theory solutions for these
operators are typically written in Fourier space it will be necessary to first
compute the Fourier transforms of these correlators before converting them
into configuration (Lagrangian) space.

As a first example let us consider the linear density-shear correlator
$E_{ab}(\boldsymbol{q})$. Taking its Fourier transform we have that

$E_{ab}(\boldsymbol{q})=\left\langle\delta(\boldsymbol{q}_{1})s_{ab}(\boldsymbol{q}_{2})\right\rangle=\text{FT}\left\\{\left(\frac{k_{a}k_{b}}{k^{2}}-\frac{1}{3}\delta_{ab}\right)P_{\rm
lin}(k)\right\\}=-\xi^{2}_{0}(q)(\hat{q}_{a}\hat{q}_{b}-\frac{1}{3}\delta_{ab})$
(3.11)

where we have used Equation 1.5 in the final step. This example is
particularly straightforward because it and its Fourier transform can be
written in terms of a single Legendre basis tensor. Typically this will not be
the case and the Fourier-space spectra have to be suitably decomposed.

The most general way to write down these correlators in Fourier space is in
the helicity basis defined in ref. [21]. Since the properties of this basis
are extensively documented in that work we will only summarize the relevant
pieces here. Briefly, the traceless symmetric tensors of rank 2, i.e. shapes,
can be described as the spin 2 component of $1\otimes 1=2\oplus 1\oplus 0$,
which is spanned by $\ket{2,m}$ with $m=0,\pm 1,\pm 2$, represented by the
spherical harmonic tensors $Y_{\ell=2,m}$ defined in Equation 1.7. Scalars,
i.e. densities, are described by the spin 0 part, $Y_{0,0}=1$. For correlators
involving purely shape and density operators, i.e. which do not contain
factors of the displacement $\Psi_{i}$, symmetry under rotations about k
requires that only projections of the same helicity $m$ are nonzero, such that
the correlators can be decomposed into sums of the products
$Y^{m}_{\ell}(\hat{k})Y^{-m}_{\ell^{\prime}}(\hat{k})$. Importantly, these
products themselves also have cancelling azimuthal dependence such that they
can re-projected into the Legendre basis [42] and Fourier-transformed as in
Equation 1.5.

For correlators involving one or more displacements the picture is similar,
except now we must deal with operators of the form

$k_{i_{1}}...k_{i_{m}}\left(\Psi_{i_{1}}...\Psi_{i_{m}}O_{ab}\right)(\boldsymbol{q}).$
(3.12)

where we have kept the powers of the wave-vector k accompanying each
displacement in LPT explicit, in the case of shapes. The case for density
operators is similar. The tensor operator in parentheses is now a rank $m+2$
tensor and as such must now be expressed in the subspace of rank $m+2$ tensors
$1\otimes\cdots\otimes 1$ where the first $m$ indices are symmetric and the
last two are symmetric and traceless. When long-wavelength linear
displacements aren’t resummed, i.e. in the Eulerian theory, we can further
express the dot product with the wavevector in front as dot products with the
angular moment eigenstate $\ket{1,0}$, such that the relevant projected
operators which contribute to the shape statistics with helicity $m$ are
restricted to the $(\bigotimes^{n}\ket{1,0})\otimes\ket{\ell,m}$ components of
the product operator in parentheses. In LPT, however, the distinction between
Lagrangian and Eulerian coordinates and wavevectors due to resummed
displacements makes this observation less useful. Furthermore, as we discuss
in §8 it is of general interest to evaluate the full tensorial structure of
displacement-shape correlators, for example to compute the effect of redshift-
space distortions, rather than restricting to the components selected out by
dotting with the wavevector. This structure however can be slightly
simplified: since any two-point correlators, with or without displacements,
must still have azimuthal symmetry about k, they will nonetheless be forced to
have total helicity zero and be expressible in terms of the Legendre tensor
basis.

Let us now explicitly describe how to evaluate these correlators in terms of
Legendre basis tensors. As an example let us take correlators of the form

$\left\langle
O_{ab}(\textbf{k})|\Psi_{i}(\textbf{k}^{\prime})\right\rangle^{\prime}=F_{3}(k)\mathcal{Q}^{3}_{iab}(\hat{k})+F_{1}(k)\mathcal{Q}^{1}_{iab}(\hat{k}).$
(3.13)

By the above arguments this correlator must be expressible in terms of the
Legendre basis and, due to parity and rank, must be expressible solely in
terms of $P_{1,3}$. The tensor $P_{3}$ can immediately be seen to be a
possible component of the correlator as it is both traceless and symmetric
under all indices; we thus define $\mathcal{Q}_{iab}(\hat{k})=P_{3,iab}$.
Since $P_{1,a}(\hat{k})=\hat{k}_{a}$ is rank 1 it must be multiplied by
factors of the Kronecker $\delta$. The two possible combinations are

$P_{1,i}(\hat{k})\delta_{ab},\quad P_{1,(a}(\hat{k})\delta_{b)i}$

where we have used that the indices $a,b$ must be symmetric. However, since
this basis must additionally be traceless in $a,b$ we have that the only
viable combination is

$\mathcal{Q}^{1}_{iab}(\hat{k})=P_{1,i}\delta_{ab}-\frac{3}{2}P_{1,(a}\delta_{b)i}$
(3.14)

leading to the two components shown above. Since $\mathcal{Q}^{1,3}$ are
written in terms of $P_{\ell}$ this expression can then be straightforwardly
Fourier transformed to yield

$\left\langle
O_{ab}(\boldsymbol{q}_{2})\Psi(\boldsymbol{q}_{1})\right\rangle=-i\xi^{3}_{0}\left[F_{3}\right](q)\mathcal{Q}^{3}_{iab}(\hat{q})+i\xi^{1}_{0}\left[F_{1}\right](q)\mathcal{Q}^{1}_{iab}(\hat{q}).$
(3.15)

In general, we can generate an appropriate Legendre basis for two-point
correlators with $n$ indices as follows:

  1. 1.

For a given Legendre tensor $P_{\ell}(\hat{k})$ of the same parity, form all
tensor products with Kronecker $\delta$’s with combined rank $n$.

  2. 2.

Symmetrize each product under the pairwise and exchange symmetries and keep
the linearly independent pieces.

  3. 3.

Solve for the linear combinations of these pieces that are trace-free in the
shape indices to generate the full basis.

We will make extensive use of this construction throughout the rest of the
paper and denote the components as $\mathcal{Q}^{\ell}_{ijk...abc...}$ where
$ijk$ and $abc$ refer to displacements and shape indices as usual. Note that
the basis elements with different $\ell$ will be orthogonal, that is to say
$\mathcal{Q}^{\ell}_{ijk...abc...}\mathcal{Q}^{\ell^{\prime}}_{ijk...abc...}=0$
for $\ell\neq\ell^{\prime}$ due to $P_{\ell}$ being traceless.

###  3.4 Angular Structure of Stochastic Terms

We can also infer the shape of stochastic contributions to shape-shape and
density-shape power spectra using the mathematical tools developed in this
section. The shape of these contributions follows from the fact that their
correlations are limited in support to small, non-perturbative scales. In the
case of the density this implies
$\left\langle\epsilon(\boldsymbol{q}_{1})\epsilon(\boldsymbol{q}_{2})\right\rangle=Af(q/R_{h})$
for some unit-normalized function $f$. The stochastic power spectrum is then

$\left\langle\epsilon(\textbf{k})|\epsilon(\textbf{k}^{\prime})\right\rangle^{\prime}=A\int
d^{3}\boldsymbol{q}\ e^{i\textbf{k}\cdot\boldsymbol{q}}\
f\left(\frac{q}{R_{h}}\right)=AR_{h}^{3}\left(1-\frac{1}{2}\sigma^{2}R_{h}^{2}k^{2}+\mathcal{O}(R_{h}^{4}k^{4})\right)$
(3.16)

where $\sigma^{2}\sim 1$ is the dimensionelss second moment of $f$.
Equivalently the form above can be derived by approximating the distribution
$f$ as a series of $\delta$ functions. For galaxies, $\epsilon$ is
qualitatively well-approximated by a Poisson-sampled random field [43], in
which case the quantity $AR_{h}^{3}$ is given by the inverse number density
$\bar{n}_{g}^{-1}$. If the galaxies are sufficiently dense then the scale-
dependent corrections tend to be small and can often be neglected.

For galaxy and halo shapes the situation is similar, with the caveats that (1)
the function $f_{abcd}$ is now a tensor and (2) the so-called “shape noise” of
galaxies is typically somewhat larger than the large-scale clustering signal
itself on scales of interest. Accounting for the symmetries of galaxy shapes
the leading-order contribution to the shape noise is given by [21]

$\left\langle\epsilon_{ab}(\textbf{k})|\epsilon_{cd}(\textbf{k}^{\prime})\right\rangle^{\prime}=BR_{h}^{3}\left(\delta_{ac}\delta_{bd}+\delta_{ad}\delta_{bc}-\frac{2}{3}\delta_{ab}\delta_{cd}\right)$
(3.17)

i.e. there is only one linearly independent contribution that does not depend
on the vector k. This is equivalent to constructing $\mathcal{Q}^{0}_{abcd}$
with the symmetries of this correlator.

At next to leading order in the expansion around the zero-lag limit of $q\to
0$ we have the scale-dependent contribution

$-\frac{1}{2}k_{i}k_{j}\int d^{3}\boldsymbol{q}\ q_{i}q_{j}\
f_{abcd}\left(\frac{\boldsymbol{q}}{R_{h}}\right)\sim-\frac{R_{h}^{5}}{2}\times
k_{i}k_{j}\left(\delta\delta\delta\right)_{ijabcd}$

where the factor of $(\delta\delta\delta)$ denotes the angular structure of
the second moment of $f_{abcd}$. Again using the symmetries of the problem it
is not hard to see that there are only two independent combinations that
involve at most two powers of the unit wavevector $\hat{k}$, equivalent to
$k^{2}\mathcal{Q}^{0}_{abcd}$ and $k^{2}\mathcal{Q}^{2}_{abcd}$:

$k^{2}\left(\delta_{ac}\delta_{bd}+\delta_{ad}\delta_{bc}-\frac{2}{3}\delta_{ab}\delta_{cd}\right),\quad
k^{2}\left(Y_{2}^{1}(\hat{k})Y_{2}^{-1}(\hat{k})+4Y_{2}^{2}(\hat{k})Y_{2}^{-2}(\hat{k})\right)_{\rm
sym}.$ (3.18)

This implies that the $k^{2}$ stochastic contributions to the three helicity
spectra are not independent but rather described by two free parameters.
Similar logic implies that the cross-stochasticity can be shown to have no
scale-independent piece

$\left\langle\epsilon(\textbf{k})|\epsilon_{ab}(\textbf{k}^{\prime})\right\rangle=CR_{h}^{5}k^{2}\left(\hat{k}_{a}\hat{k}_{b}-\frac{1}{3}\delta_{ab}\right)+...$
(3.19)

such that the leading order contribution is a quadrupole proportional to
$k^{2}$.

###  3.5 Numerical Evaluation

From the above discussion it is clear that the Lagrangian correlators in
Equation 3.6 can be written in terms of Kronecker $\delta$’s and powers of
$\hat{q}_{i}$, multiplied by scalar functions of the separation $q$. When
dotted with the appropriate wave-vectors and spherical harmonics to obtain the
components in Equation 3.3 this therefore yields the form [27, 44, 37]

$\displaystyle\int d^{3}\boldsymbol{q}\ $ $\displaystyle
e^{i\textbf{k}\cdot\boldsymbol{q}-\frac{1}{2}k_{i}k_{j}A^{<}_{ij}(\boldsymbol{q})}\
c_{n}(q)\mu^{2n(+1)}$ $\displaystyle=i^{0(+1)}4\pi\int dq\ q^{2}\
e^{-\frac{1}{2}k^{2}(X^{<}(q)+Y^{<}(q))}\sum_{\ell}f^{2n}_{\ell}(q)\left(\frac{kY^{<}(q)}{q}\right)^{\ell}c_{n}(q)j_{\ell(+1)}(kq).$

Here $c_{n}$ is a scalar function of $q$, $\mu=\hat{q}\cdot\hat{k}$ and we
have defined the kernel using confluent hypergeometric functions of the second
kind

$f^{2n}_{\ell}(q)=\left(\frac{2}{k^{2}Y^{<}}\right)^{n}U\left(-n,\ell-n+1,\frac{k^{2}Y^{<}}{2}\right).$
(3.20)

These integrals then take the form of Hankel transforms and can be rapidly
performed using the FFTLog algorithm [45].

##  4 Matter-Shape Correlators at 1-loop

We now proceed to use the bias expansion and mathematical framework described
in §2 and §3 to compute the matter-shape power spectrum. This correlation
contributes to the cosmic shear auto-spectrum as the correlation between
intrinsic alignments and lensing (‘GI’), and is also a subset of the
contributions in galaxy-galaxy lensing through the galaxy-intrinsic alignment
(‘gI’) cross correlation. We start with matter-shape correlations in the
Zel’dovich approximation before moving on to the full 1-loop calculation. This
section begins with a pedagogical derivation of the Lagrangian linear
alignment contribution to the ‘GI’ correlator in order to establish the
notation and lines of thought used to evaluate the panoply of terms which
contribute to the full one-loop expansion.

###  4.1 Matter-Shape Correlators in the Zel’dovich Approximation

As a warm-up to the general structure of the calculations we will carry out in
the rest of this paper, we begin by examining only the terms which depend on
the Zel’dovich displacement and the linear density and tidal fields. As seen
in equation 3.4 we need to compute correlators of these fields
($\Psi^{(1)},\delta,s_{ab}$) in Lagrangian coordinates $\boldsymbol{q}$. These
correlators between Lagrangian fields and displacements are the building
blocks of predictions in LPT and have been well-studied and computed in
previous works [25]. Their definitions are reviewed in Appendix C.

The linear matter-shape correlator in LPT is given by

$P_{ab}(\textbf{k})=c_{s}\int
d^{3}\boldsymbol{q}e^{i\textbf{k}\cdot\boldsymbol{q}}\langle
e^{i\textbf{k}\cdot\Delta(\boldsymbol{q}_{1},\boldsymbol{q}_{2})}\
s_{ab}(\boldsymbol{q}_{2})\rangle.$ (4.1)

Assuming the Zel’dovich approximation
($\boldsymbol{\Psi}=\boldsymbol{\Psi}^{(1)}$), from the generating function
formalism we see that this is equivalent to

$P_{ab}(\textbf{k})\ni c_{s}\int
d^{3}\boldsymbol{q}e^{i\textbf{k}\cdot\boldsymbol{q}-\frac{1}{2}k_{i}k_{j}A_{ij}(\boldsymbol{q})}\
ik_{i}B_{iab}(\boldsymbol{q}).$ (4.2)

where we have followed ref. [46] and defined
$B_{iab}(\boldsymbol{q})=\langle\Delta^{(1)}_{i}(\boldsymbol{q})s_{ab}(\boldsymbol{q}_{2})\rangle$.
The angular structure of $B_{iab}$ was computed in ref. [47] to be

$B_{iab}(\boldsymbol{q})=\hat{q}_{i}\delta_{ab}\mathcal{J}_{2}(q)+\big{[}\delta_{ia}\hat{q}_{b}+\delta_{ib}\hat{q}_{a}\big{]}\mathcal{J}_{3}(q)+\hat{q}_{i}\hat{q}_{a}\hat{q}_{b}\mathcal{J}_{4}(q)$
(4.3)

where $\mathcal{J}_{n}$ are sums of Hankel transforms of the linear power
spectrum. From the discussion surrounding Equation 3.13 it is clear that
$B_{iab}$ can also be written in terms of the Legendre basis where it is
manifest that only two of $\mathcal{J}_{2,3,4}$ are independendent due to the
trace condition. As discussed in §3 the resulting signal will have a
quadrupole structure by symmetry, and projecting out this component by dotting
it with $ik_{i}P_{2}^{ab}(\hat{k})$ gives us

$\int
d^{3}\boldsymbol{q}e^{i\textbf{k}\cdot\boldsymbol{q}-\frac{1}{2}k_{i}k_{j}A_{ij}(\boldsymbol{q})}\left(-\frac{4i}{15}\mu\xi^{1}_{-1}(q)+\frac{2i}{5}P_{3}(\mu)\xi^{3}_{-1}(q)\right),$
(4.4)

where we have defined $\mu=\hat{k}\cdot\hat{q}$, $P_{\ell}(\mu)$ is the
$\ell$-th scalar Legendre polynomial and $\xi^{\ell}_{n}$ are the generalized
correlation functions defined in §1.1. In the absence of IR-resummation
($A_{ij}\rightarrow 0$) the angular integrals in the integral above can be
done analytically as in Equation 1.5 and we simply obtain the linear theory
prediction of $\frac{2}{3}P_{\rm lin}(k)$.

Let us now look at the contribution of quadratic operators in the Zel’dovich
approximation. These come from the operators $\delta s_{ab}$ and
$s_{ac}s_{cb}$. Noting that both of these operators are products of Gaussian
fields we can directly apply Wick’s theorem in configuration space to compute
the relevant Lagrangian correlators, i.e.

$\displaystyle-\frac{1}{2}k_{i}k_{j}\langle\Delta_{i}\Delta_{j}\delta(\boldsymbol{q}_{2})s_{ab}(\boldsymbol{q}_{2})\rangle_{c}=-k_{i}k_{j}\underbrace{U_{i}}_{\langle\Delta_{i}\delta\rangle}\underbrace{B_{jab}}_{\langle\Delta_{j}s_{ab}\rangle}$
$\displaystyle-\frac{1}{2}k_{i}k_{j}\langle\Delta_{i}\Delta_{j}(s^{2}_{ab})(\boldsymbol{q}_{2})\rangle_{c}=-k_{i}k_{j}\underbrace{B_{iac}}_{\langle\Delta_{i}s_{ac}\rangle}\underbrace{B_{jcb}}_{\langle\Delta_{j}s_{cb}\rangle}.$
(4.5)

For the latter we also need to remove the trace of the right hand side, but
this is accomplished automatically by projecting into the quadrupole
$P_{2}^{ab}(\hat{k})$. Much like the linear contribution these terms can be
expressed as products of $\xi^{\ell}_{n}$ and $\mu$ in the integrand.

Finally, let us examine the contributions of cubic bias operators built from
$\\{\delta,s_{ab}\\}$. As described in §2.2 these operators’ contributions are
degenerate with the linear shape bias—let us see how this occurs in practice
in the case of $\delta^{2}s_{ab}$. In this case the only contribution at one
loop is

$ik_{i}\left\langle\Delta_{i}\delta^{2}s_{ab}(\textbf{0})\right\rangle=ik_{i}\left\langle\Psi^{(1)}_{i}(\boldsymbol{q})(\delta^{2}s_{ab})(\textbf{0})\right\rangle=i\sigma^{2}k_{i}\left\langle\Psi^{(1)}_{i}(\boldsymbol{q})s_{ab}(\textbf{0})\right\rangle$
(4.6)

which is equivalent to the linear contribution with $c_{s}=\sigma^{2}$ or,
equivalently, can be removed by using the normal-ordered operator as we will
do in this paper. The case of the other Zel’dovich cubic operators is entirely
similar.

###  4.2 Terms with One Pairwise Displacement ($\Delta$)

In the previous section we computed contributions to the matter–shape power
spectrum which arose solely due to the linear Zel’dovich displacement. At one
loop in perturbation theory we will also find contributions from the second
and third order displacements to these spectra. The contributions with one
pairwise displacement $\Delta$ will be structured as

$\langle O^{(4-n)}_{ab}(\boldsymbol{q}_{2})\ \Delta^{(n)}_{i}\rangle=-\langle
O^{(4-n)}_{ab}(\boldsymbol{q}_{2})\
\Psi^{(n)}_{i}(\boldsymbol{q}_{1})\rangle.$

This is precisely the correlator whose structure we derived using symmetry
arguments in Equation 3.13. The contribution to the power spectrum is then
computed from contraction of this basis with the hanging factor of $ik_{i}$
and the $P^{2}_{ab}(\hat{k})$ tensor, as previously carried out for the linear
Zel’dovich correlator.

Let us first consider terms due to the second-order displacement. For these
(22) terms we can use that $\Psi^{(2)}_{i}$ is related to the inverse
gradients of $\delta^{2}$ and $s^{2}$ such that

$\langle
O^{(2)}_{ab}(\textbf{k})|\Psi^{(2)}_{i}(\textbf{k}^{\prime})\rangle^{\prime}=-\frac{1}{7}\frac{ik_{i}}{k^{2}}\langle
O^{(2)}_{ab}(\textbf{k})|(\delta^{2}-\frac{3}{2}s^{2})(\textbf{k}^{\prime})\rangle^{\prime}=-\frac{ik_{i}}{k^{2}}\Big{(}\frac{k_{a}k_{b}}{k^{2}}-\frac{1}{3}\delta_{ab}\Big{)}A(k)$
(4.7)

where the last step follows from the tracelessness of $O_{ab}$. From the above
we see that $F_{3}=-\frac{2i}{5k}A$ and $F_{1}=\frac{2i}{15k}A$ are both
described by a single scalar function. Since two of the quadratic spin
operators are simply given by products of $\delta$ and $s_{ij}$ we may compute
their correlators with the left-hand quantity by Wick contraction in
configuration space. The resulting decompositions into generalized correlation
functions are

$A_{s^{2}}(k)=-4\pi\int dq\ q^{2}\ j_{2}(kq)\ \Big{(}\frac{4}{5145}\Big{)}\
\Big{[}49\xi^{0}_{0}\xi^{2}_{0}+130(\xi^{2}_{0})^{2}+36\xi^{2}_{0}\xi^{4}_{0}-45(\xi^{4}_{0})^{2}\Big{]},$
(4.8) $A_{\delta s}=-4\pi\int dq\ q^{2}\ j_{2}(kq)\
\Big{(}-\frac{4}{245}\Big{)}\
\xi^{2}_{0}\Big{[}14\xi^{0}_{0}+5\xi^{2}_{0}-9\xi^{4}_{0}\Big{]}.$ (4.9)

The third, $L^{(2)}_{ab}$, is equal to the left-hand side up to a few factors
of $k$, so that the relevant spectrum is equivalent to the autospectrum of
$\Psi^{(2)}$ which was already computed for the matter power spectrum in ref.
[24]777In fact following the procedure above decomposing $\Psi^{(2)}$ and
$L^{(2)}$ into a derivatives of a scalar we recover the FFTLog expression for
$Q_{1}$ derived in [39].:

$A_{L^{(2)}}(k)=-\frac{9}{98}Q_{1}(k).$ (4.10)

Finally, let us examine the (13) terms. The contributions due to the cubic
bias operators are equivalent to swapping in $A_{O}(k)$ for $P(k)$ when
computing $B_{iab}$ for the Zel’dovich $(1,s_{ab})$ spectrum, as previously
discussed. Similarly, contributions due to the third-order displacement can be
computed by swapping

$A_{s}(k)=\frac{5}{21}R_{1}(k).$ (4.11)

with $P(k)$ when computing e.g.
$\left\langle\Delta_{i}^{(3)}s_{ab}\right\rangle$ rather than $B_{iab}$.

###  4.3 Terms with Two $\Delta$’s

Certain spectra will receive contributions from two powers of the displacement
$\Delta$. These spectra take the form

$\langle\Delta_{i}\Delta_{j}O_{2,ab}\rangle=\langle\Psi_{2,i}\Psi_{2,j}O_{2,ab}\rangle+\langle\Psi_{1,i}\Psi_{1,j}O_{2,ab}\rangle-\langle\Psi_{1,i}\Psi_{2,j}O_{2,ab}\rangle-\langle\Psi_{2,i}\Psi_{1,j}O_{2,ab}\rangle.$
(4.12)

We see that there are three kinds of terms: the first one is a zero lag term,
the second is a correlator of two displacements at one point correlated with a
bias operator at a second, the third is the product of a displacement and bias
operator correlated with another displacement. All three terms become
identical in the zero-separation limit such that the total correlator
vanishes.

As usual it will be easier to compute the Fourier transforms of these objects.
Using the symmetry of the correlator between $i,j$ and $a,b$—note that in this
case we only want the symmetrized component in $i,j$ since we will dot these
with $k_{i}k_{j}$—we can write the Fourier transform in terms of the four-
indexed orthogonal tensor components (§3.3)

$\langle\Delta_{i}\Delta_{j}O_{2,ab}\rangle(\textbf{k})=A_{4}(k)\mathcal{Q}^{4}_{ijab}(\hat{k})+A_{2}(k)\mathcal{Q}^{2A}_{ijab}(\hat{k})+B_{2}(k)\mathcal{Q}^{2B}_{ijab}(\hat{k})+A_{0}(k)\mathcal{Q}_{ijab}^{0}(\hat{k})$
(4.13)

where we have defined

$\displaystyle\mathcal{Q}^{4}_{ijab}(\hat{k})=P^{ijab}_{4}(\hat{k}),$
$\displaystyle\mathcal{Q}^{2A}_{ijab}(\hat{k})=P^{ij}_{2}(\hat{k})\delta_{ab}-\frac{3}{4}\Big{(}P^{ia}_{2}(\hat{k})\delta_{jb}+P^{ib}_{2}(\hat{k})\delta_{aj}+P^{aj}_{2}(\hat{k})\delta_{ib}+P^{bj}_{2}(\hat{k})\delta_{ai}\Big{)},$
$\displaystyle\mathcal{Q}^{2B}_{ijab}(\hat{k})=\delta_{ij}P^{ab}_{2}(\hat{k}),$
$\displaystyle\mathcal{Q}^{0}_{ijab}(\hat{k})=\delta_{ij}\delta_{ab}-\frac{3}{2}\Big{(}\delta_{ia}\delta_{jb}+\delta_{ja}\delta_{ib}\Big{)}.$
(4.14)

The $\mathcal{Q}^{\ell}$ are orthogonal except for

$\mathcal{Q}^{2A}_{ijab}\mathcal{Q}^{2B}_{ijab}=-\frac{9}{2}$

so it is straightforward to obtain the various components by direct
contraction. Note that it is important to subtract the zero-lag term from the
$\mathcal{Q}^{0}$ piece after Fourier transforming to enforce the
$\boldsymbol{q}\rightarrow 0$ limit. We will now derive expressions for the
different two-displacement contributions for linear and quadratic alignment
operators.

####  4.3.1 Contribution to linear alignment – $\langle
s_{1,ab}\Delta^{(1)}_{i}\Delta^{(2)}_{j}\rangle$

For the linear alignment operator we will have contributions from the second-
order displacement vector. The different terms in Equation 4.12 are given by

$\langle
s_{ab}(\textbf{k})\;|\;(\Psi^{(1)}_{i}\Psi^{(2)}_{j})(\textbf{k}^{\prime})\rangle^{\prime}=2S_{ab}(\textbf{k})P(k)\int_{\textbf{p}}L^{(2)}_{i}(\textbf{p},-\textbf{k})L^{(1)}_{j}(-\textbf{p})\
P(p)$
$\langle(\Psi^{(1)}_{i}s_{ab})(\textbf{k})\;|\Psi^{(2)}_{j}(\textbf{k}^{\prime})\rangle^{\prime}=2\int_{\textbf{p}}L^{(1)}_{i}(\textbf{p})S_{ab}(\textbf{k}-\textbf{p})L^{(2)}_{j}(-\textbf{p},-\textbf{k}+\textbf{p})P(p)P(|\textbf{k}-\textbf{p}|).$
$\langle(\Psi^{(2)}_{i}s_{ab})(\textbf{k})\;|\Psi^{(1)}_{j}(\textbf{k}^{\prime})\rangle^{\prime}=2L^{(1)}_{j}(-\textbf{k})P(k)\int_{\textbf{p}}L^{(2)}_{i}(\textbf{p},\textbf{k})\
S_{ab}(-\textbf{p})P(p)$ (4.15)

where we have defined
$S_{ab}(\textbf{k})=(\hat{k}_{a}\hat{k}_{b}-\frac{1}{3}\delta_{ab}$). In order
to maintain the symmetry properties of the basis above we need to evaluate the
symmetrized version $\langle
s_{1,ab}\Delta^{(1)}_{(i}\Delta^{(2)}_{j)}\rangle$, which is in fact what
enters into the power spectrum since each $\Delta$ is always contracted with
an identical $k$. Each term is computed by a separate technique that we have
already introduced – reduction to a previous result, reduction to
configuration-space Wick contractions or direct evaluation into sums of
products of generalized correlation functions.

The first term – $\langle
s_{ab}(\textbf{k})\;|\;(\Psi^{(1)}_{i}\Psi^{(2)}_{j})(\textbf{k}^{\prime})\rangle^{\prime}$–
can be inferred from a related correlator in the galaxy power spectrum [28]

$\displaystyle\langle
s_{ab}(\textbf{k})\;|\;(\Psi^{(1)}_{i}\Psi^{(2)}_{j})(\textbf{k}^{\prime})\rangle^{\prime}$
$\displaystyle=\Big{(}\frac{k_{a}k_{b}}{k^{2}}-\frac{1}{3}\delta_{ab}\Big{)}\langle\delta(\textbf{k})|(\Psi^{(1)}_{i}\Psi^{(2)}_{j})(\textbf{k}^{\prime})\rangle^{\prime}$
$\displaystyle=-\Big{(}\frac{k_{a}k_{b}}{k^{2}}-\frac{1}{3}\delta_{ab}\Big{)}\Big{(}-\frac{3}{14}\frac{\delta_{ij}}{k^{2}}R_{1}(k)+\frac{3}{14}\frac{k_{i}k_{j}}{k^{4}}(R_{1}+2R_{2})\Big{)}.$
(4.16)

The second two terms both involve correlators of products of displacements and
bias operators with a single displacement; since displacements in LPT are
inverse gradients we can pull out a power of $ik^{-1}\hat{k}$ from both of
them such that

$\displaystyle\langle(\Psi^{(1)}_{i}s_{ab})(\textbf{k})\;|\Psi^{(2)}_{j}(\textbf{k}^{\prime})\rangle^{\prime}$
$\displaystyle=\frac{k_{j}}{k^{2}}\Big{(}Q^{s}_{3}(k)\mathcal{Q}^{3}_{iab}(\hat{k})-Q^{s}_{1}(k)\mathcal{Q}^{1}_{iab}(\hat{k})\Big{)}$
$\displaystyle\langle(\Psi^{(2)}_{i}s_{ab})(\textbf{k})\;|\Psi^{(1)}_{j}(\textbf{k}^{\prime})\rangle^{\prime}$
$\displaystyle=\frac{\,k_{j}}{k^{2}}\left[R_{3}^{s}(k)\mathcal{Q}^{3}_{iab}+R_{1}^{s}(k)\mathcal{Q}^{1}_{iab}\right].$
(4.17)

The first of these can be evaluated using the inverse-gradient form of
$\boldsymbol{\Psi}^{(2)}_{i}$ as in Equation 4.7. Explicit expressions for
these kernels can be found in Appendix D. It is worth noting that these two
integrals are not independent, in that
$R_{3}^{s}=\frac{1}{4}A_{st}-\frac{1}{2}R_{1}^{s}$. Recasting into the four-
indexed $\mathcal{Q}_{ijab}$ basis we find the final expression for the
components of the two-displacement linear alignment correlator.

$\displaystyle
A_{0}^{s}(k)=\frac{2}{105k^{2}}(R_{1}+2R_{2})+\frac{2}{3k}(Q^{s}_{1}-R^{s}_{1})$
$\displaystyle
A_{2}^{s}(k)=\frac{8}{147k^{2}}(R_{1}+2R_{2})+\frac{4}{21k}(7Q^{s}_{1}+Q^{s}_{3}-7R^{s}_{1}+R^{s}_{3})$
$\displaystyle
B_{2}^{s}(k)=\frac{4}{49k^{2}}(3R_{1}-R_{2})-\frac{10}{21k}(Q^{s}_{3}+R^{s}_{3})$
$\displaystyle
A_{4}^{s}(k)=-\frac{24}{245k^{2}}(R_{1}+2R_{2})-\frac{8}{7k}(Q^{s}_{3}+R^{s}_{3}).$

Figure 1: Contributions to the matter–shape spectrum at one-loop order.
_Left:_ The individual contributions to the matter-tide correlator compared to
the linear Eulerian prediction (in blue), the Zel’dovich approximation
(orange) and the 1-loop LPT spectrum including higher-order displacements. The
corrections are $\sim 20\%$ at $k=0.3\,h{\rm Mpc}^{-1}$. _Right:_ The
remaining matter–shape spectra at one-loop order, with all shape bias
coefficients set to one.

####  4.3.2 Quadratic Bias

We now proceed to compute the correlations between quadratic alignment
operators and two $\Delta$ operators. There are three operators at quadratic
order

$O^{(2)}_{ab}\ni(s_{ak}s_{kb},\delta s_{ab},L^{(2)}_{ab}).$

The first two of these are straightforward products of linear operators and
their contributions were already obtained in Equation 4.1. The third one,
$L^{(2)}_{ab}$, however is nonlocal and we have to compute the (13) and (22)
integrals

$\displaystyle\langle(\Psi^{(1)}_{i}L^{2}_{ab})(\textbf{k})\;|\Psi^{(1)}_{j}(\textbf{k}^{\prime})\rangle^{\prime}$
$\displaystyle=2L^{(1)}_{i}(-\textbf{k})P(k)\int_{\textbf{p}}L^{(1)}_{j}(\textbf{p})\
L^{(2)}_{ab}(-\textbf{p},\textbf{k})P(p)$ $\displaystyle\langle
L^{(2)}_{ab}(\textbf{k})\;|\;(\Psi^{(1)}_{i}\Psi^{(1)}_{j})(\textbf{k}^{\prime})\rangle^{\prime}$
$\displaystyle=2\int_{\textbf{p}}L^{(1)}_{i}(\textbf{p})L^{(1)}_{j}(\textbf{k}-\textbf{p})\
L^{(2)}_{ab}(\textbf{p},\textbf{k}-\textbf{p})P(p)P(|\textbf{k}-\textbf{p}|)$
(4.18)

where we have dropped a term using $L^{(2)}_{ab}(\textbf{p},-\textbf{p})=0$.
We can decompose the former into the form of Equation 4.17 where in this case
we substitute in the corresponding kernels $R_{1,3}^{L^{(2)}}$ whose explicit
definitions are given in Appendix D. Again, these integrals are not
independent, but are rather linear combinations of previous cubic
contributions, $kR_{1}^{L^{(2)}}=-1/10A_{\delta t}+\frac{3}{5}A_{st}$ and
$kR_{1}^{L^{(2)}}=-1/5A_{\delta t}-\frac{3}{10}A_{st}$.

For the latter, we note these are equivalent to results obtained when
computing cubic correlators of displacements [28]

$\langle
L^{(2)}_{ab}(\textbf{k})\;|\;(\Psi^{(1)}_{i}\Psi^{(1)}_{j})(\textbf{k}^{\prime})\rangle^{\prime}=\frac{3}{14}\Big{(}\frac{-k_{a}k_{b}}{k^{4}}\Big{)}\Big{(}Q_{1}(k)\delta_{ij}-(Q_{1}(k)+2Q_{2}(k))\hat{k}_{i}\hat{k}_{j}\Big{)}.$
(4.19)

Both of these integrals can then be Fourier-transformed back into real space
by first projecting into the $\mathcal{Q}$ basis in Equation 4.14 as done for
the linear alignment contribution. Combining again all the terms we have

$\displaystyle
A_{0}^{L}(k)=-\frac{1}{105k^{2}}(Q_{1}+2Q_{2})-\frac{2}{3k}R_{1}^{L^{(2)}}$
$\displaystyle
A_{2}^{L}(k)=-\frac{4}{147k^{2}}(Q_{1}+2Q_{2})-\frac{4}{21k}(7R^{L^{(2)}}_{1}-R^{L^{(2)}}_{3})$
$\displaystyle
B_{2}^{L}(k)=-\frac{2}{49k^{2}}(3Q_{1}-Q_{2})-\frac{10}{21k}R^{L^{(2)}}_{3}$
$\displaystyle
A_{4}^{L}(k)=+\frac{12}{245k^{2}}(Q_{1}+2Q_{2})-\frac{8}{7k}R^{L^{(2)}}_{3}$

This concludes the derivation of the matter–shape spectrum to one-loop order
in perturbation theory. In Figure 1 we show predictions for all matter–shape
spectra derived above at $z=0.51$, chosen to match the N-body measurements in
§7, as computed by spinosaurus, the code we present as part of this
publication.

###  4.4 IR Resummation and Shape Correlation Functions

Figure 2: The impact of the resummation scale $k_{\rm IR}$ on density–shear
observables from 1-loop (non-resummed) EPT to fully-resummed LPT. _Left:_ the
matter–shape power spectrum as $k_{\rm IR}$ is changed, assuming linear
Lagrangian shape bias. _Right:_ The matter–shape correlation function
quadrupole as a function of resummation scale, compared to the linear
alignment prediction. Note the failure of the un-resummed theory in capturing
the BAO feature, as well as the insensitivity of this feature to $k_{\rm IR}$
in LPT.

As previously mentioned, a key advantage of formulating a theory of intrinsic
alignments within the framework of LPT is the ability to naturally capture the
smoothing of the BAO feature due to long-wavelength displacements. To
illustrate this effect on intrinsic alignment observables, we calculate the
matter–linear shape power spectrum and associated correlation function (given
by $\xi^{2}_{0}$) in comparison to linear theory and the un-resummed 1-loop
prediction, shown in Figure 2. The BAO features manifests as a “dip” in the
correlation function since the signal has a quadrupolar angular dependence,
much like the BAO feature in the redshift-space quadrupole of the galaxy
density. We show these spectra as the IR resummation scale, $k_{\rm IR}$, is
varied from no-resummation to fully resummed LPT. While the broadband shape of
the power spectrum is mildly affected by $k_{\rm IR}$ at small scales, a fit
to data in practice would adjust for EFT corrections (stochastic and
counterterms) that would make up for some of the broadband difference. The
impact of this resummation on the BAO feature and correlation function is on
the other hand stark, with the BAO damped in the resummed versions but
persisting out to high $k$ and leading to wiggly features in configuration
space when unresummed. The damping is realtively robust to the choice of
resummation scale, which plays little effect on the configuration-space BAO
feature shape, reflecting that most of the damping is done by long-wavelength
modes near the BAO scale.

##  5 Galaxy Density-Shape Cross Spectrum

The final scalar–shape set of correlators we must consider are the galaxy
density–shape cross-spectra. These spectra are the building blocks for the
‘gI’ contribution to intrinsic alignments in studies of galaxy–galaxy lensing,
as well as studies cross-correlating 3D catalogs of galaxy shapes with their
positions. These cross-correlations will be of the form $(\delta_{g},M_{ab}).$
In the ‘linear alignment’ paradigm the only operators that contribute to this
cross-correlation are linear in density – $\delta$ and $s_{ab}$. However, at
1-loop this cross-spectrum relies also on correlations of higher order
operators, whose cross correlations with density bias operators must be
included in a complete model of shape statistics. We proceed to enumerate
these below. Since each bias operator includes at least one additional power
of the density field, the calculations are simpler and the vast majority can
be done via configuration-space Wick contractions.

###  5.1 Terms with no displacements

Let us list contributing spectra in terms of ordered pairs
$O_{a}(\boldsymbol{q}_{2}),O_{b}(\boldsymbol{q}_{1})$, giving expressions of
the cross-spectra for each pair. The linear term is simply

  * •

$(\delta,s_{ab})$: $E_{ab}=-\frac{2}{3}P^{2}_{ab}(\hat{q})\xi^{2}_{0}.$

These bias operators also produce a term due to the square of the first
cumulant (Eqn. 3.5)
$\langle\delta(\boldsymbol{q}_{1})s_{ab}(\boldsymbol{q}_{2})\Delta_{i}\Delta_{j}\rangle=2U_{i}B_{jab}$.
This is the only contribution due to two pairwise displacements in the non-
matter density-shape correlators.

For pure contributions due to two quadratic bias operators we have the
following

  * •

$(\delta s_{ab},\delta^{2}):\ 2\xi E_{ab}$

  * •

$(s^{2}_{ab},\delta^{2}):\ 2E_{ac}E_{cb}$

  * •

$(\delta s_{ab},s^{2}):\ 2E_{cd}C_{abcd}$

  * •

$(s^{2}_{ab},s^{2}):\ 2C_{acde}C_{cbde}$.

There are in addition to the above two contributions from $L^{(2)}_{ab}$ that
can be computed as Hankel transforms of LPT kernels using that
$L_{ab}=\partial_{i}\Psi_{j}$ [44]:

  * •

$(L^{(2)}_{ab},\delta^{2}):\
\frac{2}{7}P_{2}(\hat{q})\xi^{2}_{0}\left[Q_{8}\right]$

  * •

$(L^{(2)}_{ab},s^{2}):\
\frac{2}{21}P_{2}(\hat{q})\xi^{2}_{0}\left[Q_{s^{2}}\right]$

where we have dropped the trace components. Finally, the cross correlation of
cubic bias operators with the linear shear field is given by computing
$E_{ab}$ using $A_{O}(k)$.

Figure 3: The full set of density–shape correlators at one-loop order. Note
there is significant variation in y-axes between panels. _Left:_ Spectra
between the linear density operator and galaxy shape operators. _Center:_
Spectra between the quadratic density operator and shape operators. _Right:_
Spectra between the tidal field strength and shape operators.

###  5.2 Terms with One $\Delta$

In the case of one pairwise linear displacement we have the following
contributions

  * •

$(s_{ab},\delta^{2}):2ik_{i}U_{i}E_{ab}$

  * •

$(s_{ab},s^{2}):2ik_{i}B_{icd}C_{abcd}$

  * •

$(\delta s_{ab},\delta):ik_{i}(\xi B_{iab}+U_{i}E_{ab})$

  * •

$(s^{2}_{ab},\delta):2ik_{i}B_{iac}E_{cb}$

  * •

$(L^{(2)}_{ab},\delta):\langle
L^{(2)}_{ab,2}\Psi^{(1)}_{i,2}\delta_{1}\rangle-\langle
L^{(2)}_{ab,2}\Psi^{(1)}_{i,1}\delta_{1}\rangle$.

  * •

$(s_{ab},\delta):\langle
s_{ab}(\boldsymbol{q}_{2})\Psi_{i}^{(2)}(\boldsymbol{q}_{2})\delta(\boldsymbol{q}_{1})\rangle-\langle
s_{ab}(\boldsymbol{q}_{2})\Psi_{i}^{(2)}(\boldsymbol{q}_{1})\delta(\boldsymbol{q}_{1})\rangle$

Both terms in the penultimate expression involve previously-computed
integrals. In particular using Equation 4.18 we have

$\langle
L^{(2)}_{ab,2}\Psi^{(1)}_{i,2}\delta_{1}\rangle=\xi^{3}_{0}\left[R_{3}^{L^{(2)}}(k)\right]\mathcal{Q}^{3}_{iab}(\hat{q})-\xi^{1}_{0}\left[R_{1}^{L^{(2)}}\right]\mathcal{Q}^{1}_{iab}(\hat{q})$

and from correlators previously computed for galaxy densities [28] we have

$\langle
L^{(2)}_{ab,2}\Psi^{(1)}_{i,1}\delta_{1}\rangle=\frac{3}{7}\Big{(}\frac{2}{5}\xi^{3}_{-1}\left[Q_{5}\right]\mathcal{Q}^{3}_{iab}(\hat{q})+\frac{2}{15}\xi^{1}_{-1}\left[Q_{5}\right]\mathcal{Q}^{1}_{iab}(\hat{q})\Big{)}.$
(5.1)

The final item is a contribution due to the second-order displacement given by

$\displaystyle\langle s_{ab,2}\Psi_{i,2}^{(2)}\delta_{1}\rangle$
$\displaystyle=\xi^{3}_{0}\left[R_{3}^{s}\right]\mathcal{Q}^{3}_{iab}(\hat{q})-\xi^{1}_{0}\left[R_{1}^{s}\right]\mathcal{Q}^{1}_{iab}(\hat{q})$
$\displaystyle\langle s_{ab,2}\Psi_{i,1}^{(2)}\delta_{1}\rangle$
$\displaystyle=\frac{3}{7}\Big{(}-\frac{2}{5}\xi^{3}_{-1}\left[R_{1}+R_{2}\right]\mathcal{Q}^{3}_{iab}(\hat{q})-\frac{2}{15}\xi^{1}_{-1}\left[R_{1}+R_{2}\right]\mathcal{Q}^{1}_{iab}(\hat{q})\Big{)}$
(5.2)

where in the first line we have used our result from Equation 4.17. This
concludes the derivation of all density–shape correlators which contribute at
one-loop order. We show these spectra, computed at $z=0.51$, in Figure 3.
Except for the pure linear-Lagrangian bias term $(\delta,s_{ab})$ spectrum,
all the other spectra are 1-loop at leading order and amount to perturbative
corrections to the overall clustering signal over the range shown.

##  6 Shape–shape correlations

Finally let us compute the shape-shape two-point function of galaxies. These
correlations are particularly relevant for analyses of cosmic shear since they
directly contribute to shear-shear correlations through the measured
ellipticities of galaxies. Unlike in the density-shape correlations calculated
in the previous two subsections these shape-shape correlations involve three
distinct components given by the independent spectra of the $m=0,1,2$ helicity
components. As was the case for the galaxy density-shape correlators all of
the Lagrangian correlators needed for these calculations can be obtained by
re-arranging previous results.

###  6.1 Advection Mixes Modes

Before we enumerate all the contributions to shape-shape spectra at 1-loop let
us first return to the illustrative example of linear Lagrangian shape bias in
the Zel’dovich approximation, where all large-scale correlations are due to
the Lagrangian shear operator $s_{ab}(\boldsymbol{q})$. This operator has only
$m=0$ modes in Lagrangian space by definition.

Let us examine the generation of nonlinearities due to advection. For
simplicity we will operate in the Zel’dovich approximation and resum all
linear displacements (i.e. $k_{\rm IR}=\infty$). In this regime the leading
contribution is

$P_{s_{ab}s_{cd}}=\int d^{3}\boldsymbol{q}\
e^{i\textbf{k}\cdot\boldsymbol{q}-\frac{1}{2}k_{i}k_{j}A_{ij}}C_{abcd}(\boldsymbol{q}).$
(6.1)

In order to extract the helicity components we need to project the tensor
expression above into the helicity basis as in Equation 3.3. This projection
leads to angular integrals that can be written in terms of Legendre
polynomials of $\mu=\hat{k}\cdot\hat{q}$

  * •

$m=0$: $\int d^{3}\boldsymbol{q}\
e^{i\textbf{k}\cdot\boldsymbol{q}-\frac{1}{2}k_{i}k_{j}A_{ij}}\big{(}\frac{2}{15}\xi^{0}_{0}(q)-\frac{4}{21}\xi^{2}_{0}(q)P_{2}(\mu)+\frac{12}{35}\xi^{4}_{0}(q)P_{4}(\mu)\big{)}$

  * •

$m=1$: $-\frac{1}{2}\int d^{3}\boldsymbol{q}\
e^{i\textbf{k}\cdot\boldsymbol{q}-\frac{1}{2}k_{i}k_{j}A_{ij}}\big{(}-\frac{4}{15}\xi^{0}_{0}(q)+\frac{4}{21}\xi^{2}_{0}(q)P_{2}(\mu)+\frac{16}{35}\xi^{4}_{0}(q)P_{4}(\mu)\big{)}$

  * •

$m=2$: $\frac{1}{2}\int d^{3}\boldsymbol{q}\
e^{i\textbf{k}\cdot\boldsymbol{q}-\frac{1}{2}k_{i}k_{j}A_{ij}}\big{(}\frac{4}{15}\xi^{0}_{0}(q)+\frac{8}{21}\xi^{2}_{0}(q)P_{2}(\mu)+\frac{4}{35}\xi^{4}_{0}(q)P_{4}(\mu)\big{)}$

In the absence of the displacements in $A_{ij}$, Equation 1.5 gives that the
$m=1,2$ components would vanish, as they do in linear theory. With advection,
however, they will tend to be nonzero, as shown in Figure 4. Indeed, the
generation of nonzero helicity modes by advection represents a sizeable
contribution at the tens-of-percents level compared to the linear theory
prediction. This is analogous to the generation of B modes by CMB lensing
[48]. The contributions of cubic spin operators is again accomplished by
simply exchanging
$\xi^{\ell}_{n}\rightarrow\xi^{\ell}_{n}\left[A_{O}(k)\right]$ in the above.

Figure 4: The linear shear–shear correlator pre and post-advection. The blue
line shows the $m=0$ (and only) contribution to the pre-advection correlator
and the subsequent lines show the generation of helicity 1 and 2 spectra from
advection and non-linear evolution. The correlators are computed with $k_{\rm
IR}=0.2\,h{\rm Mpc}^{-1}$.

The angular structure above is fully general to any advected correlator
$f_{abcd}$ between two shape operators. Fundamentally this is because any such
operators can be decomposed into the three Legendre basis tensors

$\displaystyle\mathcal{Q}^{4}_{abcd}=P_{4}^{abcd}(\hat{q})$
$\displaystyle\mathcal{Q}^{2}_{abcd}=\delta_{ab}P_{2}^{cd}(\hat{q})+\delta_{cd}P_{2}^{ab}(\hat{q})-\frac{3}{4}\left(\delta_{ac}P_{2}^{bd}(\hat{q})+\delta_{ad}P_{2}^{bc}(\hat{q})+\delta_{bd}P_{2}^{ac}(\hat{q})+\delta_{bc}P_{2}^{ad}(\hat{q})\right)$
$\displaystyle\mathcal{Q}^{0}_{abcd}(\hat{q})=\delta_{ab}\delta_{cd}-\frac{3}{2}\Big{(}\delta_{ac}\delta_{bd}+\delta_{ad}\delta_{bc}\Big{)}$
(6.2)

as in §3.3, which in this case is equivalent to noting that there are three
independent helicity components describing any such correlator. These basis
tensors are orthogonal to their Fourier-space counterparts
$\mathcal{Q}^{\ell}_{abcd}(\hat{k})$ when $\ell\neq\ell^{\prime}$ and dot into
Legendre polynomials $P_{\ell}(\mu)$ otherwise, implying that the helicity
spectra are each composed of the same three functions multiplying Legendre
polynomials with coefficients depending only on $m$.

###  6.2 $(s_{ab},s_{cd})$

Beyond the lowest order contribution we have two further contributions

$P_{s_{ab}s_{cd}}=\int d^{3}\boldsymbol{q}\
e^{i\textbf{k}\cdot\boldsymbol{q}-\frac{1}{2}k_{i}k_{j}A_{ij}}\Big{(}ik_{i}\langle\Delta_{i}^{(2)}s_{2,ab}s_{1,cd}\rangle-
k_{i}k_{j}B_{iab}B_{jcd}\Big{)}$ (6.3)

The former term is given by

$\displaystyle\langle\Delta_{i}^{(2)}s_{2,ab}s_{1,cd}\rangle$
$\displaystyle=\langle\Psi_{i}^{(2)}(\boldsymbol{q}_{2})s_{ab}(\boldsymbol{q}_{2})s_{cd}(\boldsymbol{q}_{1})\rangle-\langle\Psi_{i}^{(2)}(\boldsymbol{q}_{1})s_{ab}(\boldsymbol{q}_{2})s_{cd}(\boldsymbol{q}_{1})\rangle$
$\displaystyle=\langle\Psi_{i}^{(2)}(\boldsymbol{q}_{2})s_{ab}(\boldsymbol{q}_{2})s_{cd}(\boldsymbol{q}_{1})\rangle+\langle\Psi_{i}^{(2)}(\boldsymbol{q}_{2})s_{cd}(\boldsymbol{q}_{2})s_{ab}(\boldsymbol{q}_{1})\rangle$

and can be expressed in terms of the integrals in Equation 4.17

$\displaystyle\langle\Psi_{i}^{(2)}s_{ab}|s_{cd}\rangle^{\prime}=i\left[R^{s}_{3}(k)\mathcal{Q}^{3}_{iab}(\hat{k})+R^{s}_{1}(k)\mathcal{Q}^{1}_{iab}(\hat{k})\right]\left(\frac{k_{c}k_{d}}{k^{2}}-\frac{1}{3}\delta_{cd}\right).$
(6.4)

In order to Fourier-transform this correlator, symmetrized between the two
pairs of indices $(a,b)$, $(c,d)$, we can construct again a Legendre basis for
rank-5 tensors satisfying that the last four indices are composed of two
traceless and symmetric pairs using the method described in § 3.3. This
produces a set of five basis tensors
$\mathcal{Q}^{5},\mathcal{Q}^{3A},\mathcal{Q}^{3B},\mathcal{Q}^{1A},\mathcal{Q}^{1B}$
whose construction follows the procedure described in § 3.3888As this basis is
composed of a large number of terms, we present its construction and
manipulation in the basis_iabcd.nb Mathematica notebook in the spinosaurus git
repository.. Note that this term has no helicity one or two contributions in
EPT due to the factors of $\hat{k}_{a}\hat{k}_{b}-\delta_{ab}/3$, so that the
only Eulerian contribution comes from the square of $B_{iab}$.

Figure 5: All components of shape–shape spectra at one-loop order. Each row
shows a separate helicity $m=\\{0,1,2\\}$ of the shape–shape correlators, in
descending order.

###  6.3 Linear x Quadratic

Let us now look at correlators of the form

$ik_{i}\left\langle\Delta_{i}s_{2,ab}O^{(2)}_{1,cd}\right\rangle=ik_{i}\left(\left\langle\Psi_{2,i}s_{2,ab}O^{(2)}_{1,cd}\right\rangle-\left\langle
s_{2,ab}\Psi_{1,i}O^{(2)}_{1,cd}\right\rangle\right).$ (6.5)

Two of these are straightforward to write down in terms of Zel’dovich
correlators

  * •

$(s_{ab},\delta s_{cd})$: $ik_{i}(U_{i}C_{abcd}+E_{ab}B_{icd})$

  * •

$(s_{ab},s^{2}_{cd})$: $ik_{i}(C_{abce}B_{ied}+C_{abed}B_{iec})$

where we need also to remove the $cd$ trace of the latter. Removing the $c,d$
trace ensures that the resulting correlator is symmetric under exchange of
$ab$, $cd$, though this is also accomplished directly when projecting into the
helicity basis. Additionally, these operators do not have $m=1$ helicity modes
prior to advection, as was first pointed out by ref. [1]. Finally we have

$(s_{ab},L^{(2)}_{cd}):\ ik_{i}\left\langle
s_{ab}(\boldsymbol{q}_{1})L^{(2)}_{cd}(\boldsymbol{q}_{2})\Delta_{i}\right\rangle=ik_{i}\left(\left\langle\Psi_{i}(\boldsymbol{q}_{2})L^{(2)}_{cd}(\boldsymbol{q}_{2})s_{ab}(\boldsymbol{q}_{1})\right\rangle+\left\langle\Psi_{i}(\boldsymbol{q}_{2})s_{ab}(\boldsymbol{q}_{2})L^{(2)}_{cd}(\boldsymbol{q}_{1})\right\rangle\right)$
(6.6)

Both of the correlators above can be related to ones we have already
previously computed (Eqn. 4.18)

$\displaystyle\left\langle(\Psi_{i}s_{ab})(\textbf{k})|\text{TF}\left\\{L^{(2)}\right\\}_{cd}(\textbf{k}^{\prime})\right\rangle^{\prime}=-i\left(Q^{s}_{3}(k)\mathcal{Q}^{3}_{iab}(\hat{k})-Q^{s}_{1}(k)\mathcal{Q}^{1}_{iab}(\hat{k})\right)\left(\frac{k_{c}k_{d}}{k^{2}}-\frac{1}{3}\delta_{cd}\right)$
$\displaystyle\left\langle(\Psi_{i}\text{TF}\left\\{L^{(2)}\right\\}_{cd})(\textbf{k})|s_{ab}(\textbf{k}^{\prime})\right\rangle^{\prime}=i\left(R^{L^{(2)}}_{3}(k)\mathcal{Q}^{3}_{icd}(\hat{k})+R^{L^{(2)}}_{1}(k)\mathcal{Q}^{1}_{icd}(\hat{k})\right)\left(\frac{k_{a}k_{b}}{k^{2}}-\frac{1}{3}\delta_{ab}\right)$

where we have explicitly written the trace-free component in both cases. In
order to project these onto the $Q^{\ell}_{iabcd}$ basis, it is necessary to
further symmetrize the trace-removed correlators above under pair swap—while
these correlators are not themselves symmetric under pair swap, their Fourier
transform after dotting with $k_{i}$ and IR resummation must be due to
parity999It is straightforward to check that the Fourier-space correlators
above are also symmetric in this way after contracting with $k_{i}$., so this
step can be safely performed, though it is also sufficient to project with the
$Q^{\ell}_{iabcd}$ basis directly to remove offending components. Note also
that these spectra again have no helicity one or two contributions in EPT due
to the factors of $\hat{k}_{a}\hat{k}_{b}-\delta_{ab}/3$.

###  6.4 Quadratic x Quadratic

Finally let us consider the case of terms involving two quadratic bias
operators. In this case at 1-loop order no further displacements can be
involved so we need only consider Lagrangian correlators of the form
$\left\langle O_{1,ab}O_{2,cd}\right\rangle$ whose angular structure is
described by the basis in Equation 6.2.

Several of the terms with two quadratic operators can be written purely in
terms of Zel’dovich correlators

  * •

$(\delta s_{ab},\delta s_{cd})$: $\xi C_{abcd}+E_{ab}E_{cd}$

  * •

$(\delta s_{ab},s^{2}_{cd})$: $E_{ce}C_{abed}+E_{ed}C_{abce}$

  * •

$(s^{2}_{ab},s^{2}_{cd})$: $C_{aecf}C_{ebfd}+C_{aefd}C_{ebcf}$

while the remaining ones can be reduced to what we’ve already computed
previously, i.e.

$\left\langle
O^{(2)}_{ab}(\textbf{k})|L^{(2)}_{cd}(\textbf{k}^{\prime})\right\rangle^{\prime}=-\left(\frac{k_{a}k_{b}}{k^{2}}-\frac{1}{3}\delta_{ab}\right)\left(\frac{k_{c}k_{d}}{k^{2}}-\frac{1}{3}\delta_{cd}\right)A_{O^{(2)}}(k),$
(6.7)

where $A_{O^{(2)}}$ have been previously derived in the correlators of eqns.
4.8, 4.9 and 4.10. Note that these latter correlators have only $m=0$
components prior to advection.

The entire set of helicity spectra for shape–shape correlators at one loop is
shown in Figure 5. Each row shows the $m=\\{0,1,2\\}$ helicity spectra
separated by columns of shape operator. The first column shows shape–shape
correlators with dependence on the linear tidal field $s_{ab}$, the second
those with dependence on the quadratic term $\delta s_{ab}$ and finally those
with dependence on $s_{ac}s_{cb}$. The $(\delta s_{ab},\delta s_{ab})$ seem to
have the highest amplitude for higher helicity spectra, and indeed we will see
in §7 that this is also the most important contribution to the higher helicity
spectra in the so-called local Lagrangian bias ansatz [1]. We will return to
this point shortly when we compare our results to $N$-body simulations.

##  7 Comparison to N-body

Figure 6: Perturbation theory fits to halo shape spectra of helicities
$m=0,1,2$ with shape noise (SN) subtracted for clarity. The higher helicity
spectra are suppressed relative to $m=0$ and are multiplied by constant
factors to enhance their visibility. Solid and dotted lines show the results
when using the full Lagrangian effective model with (LPT) and without (CSN)
scale-dependent shape noise, while the dashed line shows the fit to the
(non-)linear alignment model, also known as NLA. For lower mass halos both LPT
and CSN represent excellent fits with reduced $\chi^{2}$ around one, while the
scale-dependent shape noise is required to obtain a good fit for more massive
halos. The NLA model fails to fit the $m=0$ spectrum past $k=0.1\,h{\rm
Mpc}^{-1}$ and predicts zero for the higher-helicity spectra.

In order to validate our perturbation theory we compare its predictions to the
clustering of the three-dimensional shapes of halos measured in N-body
simulations. In particular we will use the helicity auto spectra measured in
ref. [1], which were taken from 8 boxes each with volume
$(1.5h^{-1}\text{Gpc})^{3}$. Unlike the projected $E$ and $B$ mode statistics
used in previous works validating perturbation theory predictions for galaxy
shapes, these spectra fully describe the three-dimensional clustering without
mixing independent components, and moreover do not contain projection effects
which break the isotropy of galaxy shapes in real space [22]. Due to the
orthogonality of the helicity components the $m=0,1,2$ spectra can be treated
as diagonal, and we will assume their covariance can be approximated as
Gaussian.

The left and right panels of Figure 6 show the best fits of the theory to
halos in the mass bins $12.0<\log_{10}(Mh\,/M_{\odot})<12.5$ and
$13.0<\log_{10}(Mh\,/M_{\odot})<13.5$ at redshift $z=0.51$, along with the
shape-noise subtracted helicity spectra from N-body simulations. We also show
a comparison with the standard (non)-linear alignment model adopted by stage-
III cosmic shear surveys such as KiDS-1000 [3]. The NLA model is constructed
by taking the linear Eulerian prediction for $(s_{ab},s_{cd})$ and mapping the
power spectrum to the non-linear power spectrum using a fitting function such
as Halofit [49].

Lagrangian perturbation theory at 1-loop order, with all the independent bias
and effective-theory parameters free, provides an excellent fit to all three
helicity spectra for both the high and low mass N-body halos up to $k=0.3h$
Mpc-1, as shown by the solid lines, with a reduced $\chi^{2}$ in both cases of
approximately one. We note that these fits are performed assuming the full
statistical power of the N-body volume ($27h^{-3}$ Gpc3). The NLA model, shown
in dashed lines, on the other hand, cannot fit the $m=0$ spectrum past scales
of $\sim 0.1\,h{\rm Mpc}^{-1}$ despite using the fully nonlinear matter power
spectrum. Additionally, since the fiducial NLA modelling of shear surveys uses
a fixed value for the shape noise (assumed to be the Poisson value) that is
subtracted from the spectrum, with no additional freedom, the NLA model wholly
fails to capture the higher helicity spectra – which have been detected to be
non-zero at high significance even after shot noise subtraction.

The dotted lines show the prediction of perturbation theory without the scale-
dependent ($k^{2}$) stochastic term. This term is rather non-negligible for
the high mass bin, particularly in the $m=1,2$ spectra which are 1-loop at
leading order, and setting it to zero in our fits incurs a $\chi^{2}$ penalty
mostly from these two curves. At first sight this may be surprising, given
that halo densities in the same mass bins are well-described assuming this
term is negligible [37, 50]. However, the importance of this term can be
understood by noting that the deterministic bias contributions for shapes are
suppressed by more than an order of magnitude compared to those of densities,
such that the leading (constant) stochastic contribution ($\propto R_{h}^{3}$)
is a much larger fraction of the overall signal, as shown by the black dot-
dashed line. Once this is taken into account, the scale dependence of the
stochastic contribution ($\propto R_{h}^{5}k^{2}$) is then roughly of the size
preferred by the data. It is worth noting again that while there are three
helicity spectra there are only two independent $k^{2}$ terms, which are
enough to fit the data. Nonetheless, the prominence of the scale-dependent
contribution for the $m=1,2$ spectra suggests that shape spectra observed in
galaxy surveys, particularly for the $B$ modes, cannot be perturbatively
modeled well past what is shown here even at higher orders in perturbation
theory. On the other hand, these B mode contributions are all extremely
subdominant to the shape noise, making their signal relatively suppressed and
hard to detect in realistic surveys.

Figure 7: Fits to halo helicity spectra, as in Fig. 6. The dashed lines now
show the ‘local Lagrangian bias’ model fit to these catalogs. The $\chi^{2}$
for LLB is approximately 5 (17)% higher for the low (high) mass bins compared
to the full LPT fit.

The above discussion centered on fits to galaxy shape spectra involving all
free parameters at 1-loop in perturbation theory. Recently, ref. [1] showed
that halos are, at least approximately, well described as Lagrangian tracers
of the initial shear field, up to multiplication by the (scalar) density. This
suggests that an approximate biasing scheme for halos can be written as

$M_{ab}(\boldsymbol{q})=\left(1+b_{1}\delta(\boldsymbol{q})+\frac{1}{2}b_{2}(\delta^{2}(\boldsymbol{q})-\left\langle\delta^{2}\right\rangle)+b_{s}(s^{2}(\boldsymbol{q})-\left\langle
s^{2}\right\rangle)+...\right)\ c_{s}s_{ab}(\boldsymbol{q}).$ (7.1)

Since the cubic operators in this product are all degenerate with the linear
bias $c_{s}$ this leads to a parameterization of the shape field with only two
free dimensionless bias operators: $s_{ab}$ and $\delta s_{ab}$. The dashed
lines in Fig. 7 show the fits using this reduced parameter set, along with
effective-theory contributions, when compared to the full LPT model described
above. Evidently the fit is not perfect and we are unable to fully recover the
shape power spectra on large scales; on the other hand, the resulting
difference with N-body data is only at the $2\%$ or smaller level in terms of
the linear bias. The shift in $\chi^{2}$ from this local Lagrangian bias (LLB)
model is $\sim$ 5 (17) % for the lower (higher) mass bins, respectively. The
shifts in $\chi^{2}$ are higher than the change in degrees of freedom from
restricting the parameter space of the alignment model. In addition, we note
that the our fits to the full model are consistent with a Lagrangian bias
model wherein only the Zeldovich terms ($s_{ab},\delta s_{ab},s^{2}_{ab}$), in
addition to the EFT corrections, are nonzero. This result is in good agreement
with ref. [1], though it should be treated with caution since we have not
performed a proper MCMC analysis of the simulated data. We intend to return to
this topic in the near term.

Nevertheless, since the intrinsic alignment signal is a small contaminant in
cosmic shear data for source/lens redshift distributions that do not possess
significant overlap, it seems likely that a reduced bias parameterization like
LLB or Zeldovich-only can be used in lieu of the full model, assuming that the
bias relations of halos continue to hold for galaxies. The net result would be
a reduced parametrization for the galaxy shapes that is nonetheless more
accurate than the TATT or NLA models alone. Analyses with such a reduced model
could have implications for even current cosmic shear analyses [51], where
degeneracies between intrinsic alignment model parameters (both bias
coefficient amplitudes and their redshift-dependence) currently limit
constraining power.

##  8 Conclusions

The cosmological distribution of galaxy shapes, known as intrinsic alignments,
play an important role in current and upcoming galaxy surveys, both as one of
the most significant contaminants in cosmic shear surveys and as a probe of
large scale structure in their own right. In recent years, significant
progress has been made towards the modeling these large-scale shape
correlations, particularly within the language of perturbation theory, leading
to consecutive development of the (non)linear alignment (NLA), quadratic tidal
alignment and tidal torquing (TATT) and finally effective field theory (EFT)
models [6, 20, 21]. The EFT model, in particular, represents a consistent
treatment of perturbative contributions to intrinsic alignments up to a given
order in the perturbation theory, and was shown in ref. [22] to be in
excellent agreement with the shapes of halos in N-body simulations.

In this work we have continued this program, particularly as set out in ref.
[21], by formulating a description of intrinsic alignments within Lagrangian
perturbation theory (LPT). Within this description the dynamics of galaxies
are computed in terms of their displacements $\Psi$ from initial Lagrangian
positions $\boldsymbol{q}$, and their shapes described by a bias expansion in
terms of their local initial conditions written as spatial derivatives of
these displacements. This advective description of large-scale structure
allows for a natural way to understand the large-scale displacements that damp
the baryon acoustic oscillation (BAO) peak and other oscillatory features in
the power spectrum, whose effect we resum in our calculations. We enumerate
the contributions to this bias expansion up to cubic order and derive the
1-loop power spectrum of both galaxy shape–shape and galaxy density–shape
correlations—the latter includes the special case of matter when all the
Lagrangian bias coefficients are set to zero—representing a complete list of
the relevant contributions to cosmic shear (‘GI’ and ‘II’) and galaxy-galaxy
lensing (‘gI’) measurements as well as the real-space contributions to
spectroscopic shape analyses. At this order we require one linear, three
quadratic and two cubic bias parameters. Our formalism also includes
effective-theory counterterms and stochastic contributions necessary to
control the impact of small-scale physics, including two next-to-leading order
scale-dependent contributions to the latter. We make spinosaurus, a Python
code that performs these calculations using Hankel transforms and FFTLog,
publically available.

As a test of the theory model developed in this work, we fit our theory
predictions to the three-dimensional helicity spectra of simulated halos
measured in ref. [1]. These fits capture the full angular structure of galaxy
spin statistics at the 2-point function level and are distinct from previous
comparisons in the literature which used only the statistics of projected E
and B modes. Perturbation theory at 1-loop order gives an excellent fit to the
shape-shape autospectra of these halos, with $\chi^{2}$ per degree of freedom
close to unity, out to $k=0.3h$ Mpc-1 at $z=0.51$, and we find a strong
preference for including a $k^{2}$ scale dependence to the shot noise,
particularly in order to fit the $m=2$ spectrum. We also test the “local
Lagrangian bias” ansatz, previously seen to be a good approximation at the
field [1] and projected 2-point function [22, 52] levels. This ansatz involves
only two free bias parameters in addition to the counter and stochastic
terms—we find that it is able to reasonably describe the data in qualitative
terms but results in a mis-estimation of the linear bias at the few percent
level when fit over the same range of scales, leading to a non-negligible
$\chi^{2}$ difference. Nonetheless, we note that our fits are performed for
full three-dimensional measurements with total volume $27h^{-3}$ Gpc3,
representing greater constraining power than near-term cosmic shear or
spectroscopic shape surveys, such that any of the schemes we have introduced
in this work would likely be acceptable in realistic applications despite the
full model being the only one that adequately fits the simulations at their
true statistical power.

Let us close by mentioning a few possible extensions to this work. Most
directly, our model can be applied to analyses of existing weak lensing data
sets for both cosmic shear alone and galaxy–galaxy lensing, though scale cuts
would necessarily need to be applied in order to limit the contributions of
intrinsic alignments to these data to perturbative scales. Recently, ref. [52]
proposed a hybrid model [30] of intrinsic alignments wherein Lagrangian bias
operators are advected by the non-perturbative displacements of dark matter
solved for in N-body simulations, with the aim of extending the reach of the
perturbative bias expansion beyond the dynamical nonlinear scale. The success
of such an approach is related to the ability of the limited Lagrangian model
described above to qualitatively describe the data. Our calculations in this
paper provide direct predictions for the components of the hybrid model and
will be essential in producing a viable emulator for it to interpolate across
different cosmologies in data analyses. Our calculations similarly serve as
building blocks for extensions of the technique of ‘Zel’dovich control
variates’ [53, 54, 55] to the realm of intrinsic alignment statistics –
allowing for more accurate and precise measurements of IA power spectra from
$N$-body simulations.

Recently, there has been a renewed interest in measuring intrinsic alignments
not just as a contaminant in lensing surveys but as a cosmological signal in
their own right, particularly in the context of spectroscopic surveys. Perhaps
the most exciting extension of this work is to the modeling of galaxy shape
statistcs in redshift space, where these surveys operate. In such surveys
anisotropies known as redshift-space distortions arise due to the degeneracy
between cosmological redshifts and peculiar velocities along the line of
sight, and manifest as a boost $\textbf{u}=\textbf{v}_{\parallel}/\mathcal{H}$
to the line-of-sight positions of galaxies. Within Lagrangian perturbation
this amounts to boosting the order $n$ displacement by [24, 25]

$\Psi^{(n),s}_{i}(\boldsymbol{q})=\Psi^{(n)}_{i}(\boldsymbol{q})+\mathbf{u}_{i}^{(n)}=R^{(n)}_{ij}\Psi^{(n)}_{j}(\boldsymbol{q})$
(8.1)

where we have defined the matrix
$R^{(n)}_{ij}=\delta_{ij}+nf(a)\hat{n}_{i}\hat{n}_{j}$, such that the redshift
space power spectrum becomes

$P_{s,abcd}(\textbf{k})=\int d^{3}\boldsymbol{q}\
e^{i\textbf{k}\cdot\boldsymbol{q}}\left\langle e^{i\textbf{k}\cdot\Delta_{s}}\
F_{ab}(\boldsymbol{q}_{1})F_{cd}(\boldsymbol{q}_{2})\right\rangle.$ (8.2)

where $\Delta_{s,i}^{(n)}=R^{(n)}_{ij}\Delta^{(n)}_{j}$ is the pairwise
distortion in redshift space. Calculating the power spectrum in redshift space
is then equivalent to boosting the pairwise-displacement and bias operator
correlators we have already computed at each order $n$ with the matrix
$R^{(n)}_{ij}$, though we caution that the resulting angular integrals will be
slightly more complicated than in this work [56, 29] – the resulting power
spectra will now only be symmetric about rotations about the line-of-sight and
thus require a different tensor basis [34, 35, 36]. Alternatively the
displacements can be perturbatively expanded to yield a moment expansion in
the velocities, also known as the distribution function approach to redshift-
space distortions [57, 56, 37]. In this case we can write

$P_{s,abcd}(\textbf{k})=\sum_{n=0}^{\infty}\frac{i^{n}}{n!}k_{i_{1}}...k_{i_{n}}\int
d^{3}\boldsymbol{q}\ e^{i\textbf{k}\cdot\boldsymbol{q}}\left\langle
e^{i\textbf{k}\cdot\Delta}\Delta\textbf{u}_{i_{1}}...\Delta\textbf{u}_{i_{n}}F_{ab}(\boldsymbol{q}_{1})F_{cd}(\boldsymbol{q}_{2})\right\rangle$
(8.3)

where again $\Delta\textbf{u}$ are simply scalar multiples of $\Delta$ given
by their time derivatives, precisely what we have painstakingly computed in
this paper to 1-loop order. We intend to return to this topic in future work.

## Acknowledgments

We thank Joe DeRose, Zvonimir Vlah, Martin White for useful discussions and
Elisabeth Krause for helpful comments on a draft of this work. We especially
thank Kazuyuki Akitsu for frequent conversations throughout the preparation of
this work that greatly aided our understanding of the subject, as well as for
providing his measurements of halo shape spectra. S.C. is supported by the
Bezos Membership at the Institute for Advanced Study. N.K. would like to thank
the International Center for Theoretical Physics – South American Institute
for Fundamental Research (ICTP-SAIFR) for their hospitality during part of the
completion of this work. ICTP-SAIFR is funded by FAPESP grant 2021/14335-0.
Figures and code to produce them in this work have been made using SciPy Stack
[58, 59, 60]. This research has made use of NASA’s Astrophysics Data System
and the arXiv preprint server.

##  Appendix A Tensor analogs of Legendre polynomials

We collect definitions of the tensor Legendre polynomials used throughout the
text, up to fourth order. The fifth order polynomial contains many terms and
we refer to the included Mathematica notebooks for its presentation and use.

$\displaystyle P_{0}(\hat{k})=1$ (A.1) $\displaystyle
P_{1}^{i}(\hat{k})=\hat{k}_{i}$ (A.2) $\displaystyle
P_{2}^{ij}(\hat{k})=\frac{3}{2}\left(\hat{k}_{i}\hat{k}_{j}-\frac{1}{3}\delta_{ij}\right)$
(A.3) $\displaystyle
P_{3}^{iab}(\hat{k})=\frac{1}{2}\left(5\hat{k}_{i}\hat{k}_{a}\hat{k}_{b}-(\hat{k}_{i}\delta_{ab}+\hat{k}_{(a}\delta_{b)i}\right)$
(A.4) $\displaystyle
P_{4}^{ijab}(\hat{k})=\frac{1}{8}\left(35\hat{k}_{i}\hat{k}_{j}\hat{k}_{a}\hat{k}_{b}-5(\hat{k}_{i}\hat{k}_{j}\delta_{ab}+\hat{k}_{a}\hat{k}_{b}\delta_{ij}+\hat{k}_{i}\hat{k}_{(a}\delta_{b)j}+\hat{k}_{j}\hat{k}_{(a}\delta_{b)i})+(\delta_{ij}\delta_{ab}+\delta_{i(a}\delta_{b)j})\right)$
(A.5)

##  Appendix B Order-by-order contributions to derivatives of the generating
function at 1-loop

In this appendix, we enumerate contributions to the generating function, or
its derivatives, to 1-loop order in perturbation theory.

###  B.1 Matter Tracer Cross Spectrum

In the case of the matter-tracer cross spectrum we need to look for terms with
only one power of $\alpha$, i.e.

$\left.\frac{\partial
M}{\partial\alpha_{2}^{a}}\right|_{\boldsymbol{\alpha}=0}=M(\mathbf{0})\Big{(}ik_{i}\langle\Delta_{i}O_{2}^{a}\rangle-\frac{1}{2}k_{i}k_{j}\langle\Delta_{i}\Delta_{j}O_{2}^{a}\rangle+...\Big{)}$

The contributions are

$\displaystyle\langle\Delta_{i}O_{2}^{a}\rangle\ni\langle\Delta_{i}^{(1)}O_{2}^{a,(1)}\rangle,\langle\Delta_{i}^{(2)}O_{2}^{a,(2)}\rangle,\langle\Delta_{i}^{(1)}O_{2}^{a,(3)}\rangle,\langle\Delta_{i}^{(3)}O_{2}^{a,(1)}\rangle$
$\displaystyle\langle\Delta_{i}\Delta_{j}O_{2}^{a}\rangle\ni\langle\Delta_{i}^{(2)}\Delta_{j}^{(1)}O_{2}^{a,(1)}\rangle,\langle\Delta_{i}^{(1)}\Delta_{j}^{(1)}O_{2}^{a,(2)}\rangle$
(B.1)

###  B.2 Tracer-Tracer Spectrum

For the tracer–tracer spectrum two derivatives are needed

$\left.\frac{\partial^{2}M}{\partial\alpha_{1}^{a}\partial\alpha_{2}^{b}}\right|_{\boldsymbol{\alpha}=0}=M(\mathbf{0})\Big{(}-k_{i}k_{j}\langle\Delta_{i}O_{1}^{a}\rangle\langle\Delta_{j}O_{2}^{b}\rangle+\langle
O_{1}^{a}O_{2}^{b}\rangle+ik_{i}\langle\Delta_{i}O_{1}^{a}O_{2}^{b}\rangle_{c}+...\Big{)}$

Again we can enumerate the relevant contributions by order:

$\displaystyle\langle\Delta_{i}O_{1}^{a}\rangle\langle\Delta_{j}O_{2}^{b}\rangle\ni\langle\Delta_{i}^{(1)}O_{1}^{a,(1)}\rangle\langle\Delta_{j}^{(1)}O_{2}^{b,(1)}\rangle$
$\displaystyle\langle O_{1}^{a}O_{2}^{b}\rangle\ni\langle
O_{1}^{a,(1)}O_{2}^{b,(1)}\rangle,\langle
O_{1}^{a,(1)}O_{2}^{b,(3)}\rangle,\langle
O_{1}^{a,(3)}O_{2}^{b,(1)}\rangle,\langle O_{1}^{a,(2)}O_{2}^{b,(2)}\rangle$
$\displaystyle\langle\Delta_{i}O_{1}^{a}O_{2}^{b}\rangle_{c}\ni\langle\Delta_{i}^{(1)}O_{1}^{a,(1)}O_{2}^{b,(2)}\rangle_{c},\langle\Delta_{i}^{(1)}O_{1}^{a,(2)}O_{2}^{b,(1)}\rangle_{c},\langle\Delta_{i}^{(2)}O_{1}^{a,(1)}O_{2}^{b,(1)}\rangle_{c}$
(B.2)

##  Appendix C Definitions of LPT correlators

In this section we review the Lagrangian correlators defined in refs. [47,
46]. Specifically we will list the correlators including at least one power of
the Lagrangian tidal field $s_{ab}(\boldsymbol{q})$. Readers interested in the
remaining correlators (e.g.$\ A_{ij},U_{i},\xi_{L}$ etc.) are directed to ref.
[25]. The first correlator as discussed in §3.3 is

$E_{ab}(\boldsymbol{q})=\left\langle\delta(\boldsymbol{q}_{1})s_{ab}(\boldsymbol{q}_{2})\right\rangle=\mathcal{J}_{1}(q)(\delta_{ab}/3-\hat{q}_{a}\hat{q}_{b})=-\frac{2}{3}\xi^{2}_{0}(q)P_{2,ab}(\hat{q}).$
(C.1)

where $\mathcal{J}_{1}=\xi^{2}_{0}$. Ref. [47] defined a series of functions
$\mathcal{J}_{n}$ via sums of generalized correlation functions, and we will
follows that convention here. The next is the shear-displacement correlator

$B_{iab}(\boldsymbol{q})=\left\langle\Delta
s_{ab}(\boldsymbol{q}_{2})\right\rangle=\delta_{ab}\hat{q}_{i}\mathcal{J}_{2}+(\delta_{ia}\hat{q}_{b}+\delta_{ib}\hat{q}_{a})\mathcal{J}_{3}+\hat{q}_{i}\hat{q}_{a}\hat{q}_{b}\mathcal{J}_{4}$

where

$\mathcal{J}_{2}=\frac{2}{15}\xi^{1}_{-1}-\frac{1}{5}\xi^{3}_{-1},\
\mathcal{J}_{3}=-\frac{1}{5}\xi^{1}_{-1}-\frac{1}{5}\xi^{3}_{-1},\
\mathcal{J}_{4}=\xi^{3}_{-1}.$ (C.2)

As discussed in §3.3 this correlator can equivalently be written in the
Legendre basis adopted in this work as

$B_{iab}(\boldsymbol{q})=\frac{2}{5}\xi^{3}_{-1}(q)\mathcal{Q}^{3}_{iab}(\hat{q})+\frac{2}{15}\xi^{1}_{-1}(q)\mathcal{Q}^{1}_{iab}(\hat{q}).$
(C.3)

Finally the shear-shear correlator is

$\displaystyle C_{abcd}(\boldsymbol{q})=\left\langle
s_{ab}(\boldsymbol{q}_{1})s_{cd}(\boldsymbol{q}_{2})\right\rangle=\ $
$\displaystyle\delta_{ab}\delta_{cd}\mathcal{J}_{5}+(\delta_{ac}\delta_{bd}+\delta_{ad}\delta_{bc})\mathcal{J}_{6}$
$\displaystyle+(\delta_{ab}\hat{q}_{c}\hat{q}_{d}+\delta_{cd}\hat{q}_{a}\hat{q}_{b})\mathcal{J}_{7}$
$\displaystyle+(\delta_{ac}\hat{q}_{b}\hat{q}_{d}+\delta_{bd}\hat{q}_{a}\hat{q}_{c}+\delta_{ad}\hat{q}_{b}\hat{q}_{c}+\delta_{bc}\hat{q}_{a}\hat{q}_{d})\mathcal{J}_{8}$
$\displaystyle+\hat{q}_{a}\hat{q}_{b}\hat{q}_{c}\hat{q}_{d}\mathcal{J}_{9}$

where

$\mathcal{J}_{5}=-\frac{1}{315}(14\xi^{0}_{0}+40\xi^{2}_{0}-9\xi^{4}_{0}),\
\mathcal{J}_{6}=\frac{1}{105}(7\xi^{0}_{0}+10\xi^{2}_{0}+3\xi^{4}_{0})$
$\mathcal{J}_{7}=\frac{1}{21}(4\xi^{2}_{0}-3\xi^{4}_{0}),\
\mathcal{J}_{8}=-\frac{1}{7}(\xi^{2}_{0}+\xi^{4}_{0}),\
\mathcal{J}_{9}=\xi^{4}_{0}.$ (C.4)

Again we can equivalently express this correlator in terms of the Legendre
basis in which

$C_{abcd}(\boldsymbol{q})=\frac{8}{35}\xi^{4}_{0}(q)\mathcal{Q}^{4}_{abcd}(\hat{q})+\frac{8}{63}\xi^{2}_{0}(q)\mathcal{Q}^{2}_{abcd}(\hat{q})-\frac{2}{45}\xi^{0}_{0}(q)\mathcal{Q}^{0}_{abcd}(\hat{q}).$
(C.5)

##  Appendix D Explicit Expressions for Loop Integrals

###  D.1  $\left\langle s_{1,ab}\Delta^{(1)}_{i}\Delta^{(2)}_{j}\right\rangle$

For the first 22-type term we can write

$\displaystyle\langle(\Psi^{(1)}_{i}s_{ab})(\textbf{k})\;|\Psi^{(2)}_{j}(\textbf{k}^{\prime})\rangle^{\prime}$
$\displaystyle=-\frac{1}{7}\frac{ik_{j}}{k^{2}}\langle(\Psi^{(1)}_{i}s_{ab})(\textbf{k})\;|\;\big{(}\delta^{2}-\frac{3}{2}s^{2}\big{)}(\textbf{k}^{\prime})\rangle^{\prime}$
$\displaystyle=-\frac{1}{7}\frac{ik_{j}}{k^{2}}\text{FT}\Big{\\{}2U_{i}(\boldsymbol{q})E_{ab}(\boldsymbol{q})-3B_{icd}(\boldsymbol{q})C_{abcd}(\boldsymbol{q})\Big{\\}}$
$\displaystyle=-\frac{ik_{j}}{k^{2}}\text{FT}\Big{\\{}\tilde{F}^{s}_{3}(q)\mathcal{Q}^{3}_{iab}(\hat{q})+\tilde{F}^{s}_{1}(q)\mathcal{Q}^{1}_{iab}(\hat{q})\Big{\\}}$
$\displaystyle\equiv\frac{k_{j}}{k^{2}}\Big{(}Q_{3}^{s}(k)P_{3}(\hat{k})-Q_{1}^{s}(k)(\hat{k}_{i}\delta_{ab}-\frac{3}{2}\hat{k}_{(a}\delta_{b)i})\Big{)}$
(D.1)

where $\tilde{F}^{s}_{1,3}$ can be obtained by projecting the correlators in
curly brackets in the second line into the $\mathcal{Q}(\hat{q})$ basis, and
$Q^{s}_{1,3}=\xi^{1,3}_{0}[\tilde{F}^{s}_{1,3}]$.

The integrals for the 13-type term can similalry be directly projected onto
the $\mathcal{Q}(\hat{k})$ basis such that the integrals can be rewritten in
terms of generalized correlation functions:

$\displaystyle R_{3}^{s}(k)=\frac{2}{5}P(k)\int_{0}^{\infty}dq\ q\ $
$\displaystyle\Big{[}-\frac{2}{35}k\xi^{0}_{0}(q)j_{0}(kq)+\frac{18}{245}\xi^{1}_{1}(q)j_{1}(kq)+\frac{10}{49}k\xi^{2}_{0}(q)j_{2}(kq)$
$\displaystyle-\frac{22}{105}\xi^{3}_{1}(q)j_{3}(kq)-\frac{36}{245}k\xi^{4}_{0}(q)j_{4}(kq)+\frac{20}{147}\xi^{5}_{1}(q)j_{5}(kq)\Big{]}$
(D.2) $\displaystyle R_{1}^{s}(k)=\frac{1}{15}P(k)\int_{0}^{\infty}dq\ q\ $
$\displaystyle\Big{[}\frac{4}{35}k\xi^{0}_{0}(q)j_{0}(kq)+\frac{12}{35}\xi^{1}_{1}(q)j_{1}(kq)-\frac{20}{49}k\xi^{2}_{0}(q)j_{2}(kq)-\frac{12}{35}\xi^{3}_{1}(q)j_{3}(kq)$
$\displaystyle+\frac{72}{245}k\xi^{4}_{0}(q)j_{4}(kq)\Big{]}.$ (D.3)

###  D.2 $\left\langle
O^{(2)}_{ab}\Delta^{(1)}_{i}\Delta^{(1)}_{j}\right\rangle$

Again pulling out a factor of $ik^{-1}\hat{k}$ we can write

$\left\langle\Psi^{(1)}_{i}L^{(2)}_{ab}(\textbf{k})|\Psi^{(1)}_{j}(\textbf{k}^{\prime})\right\rangle=\frac{k_{j}}{k^{2}}\left[R_{3}^{L^{(2)}}(k)\mathcal{Q}^{3}_{iab}(\hat{k})+R_{1}^{L^{(2)}}(k)\mathcal{Q}^{1}_{iab}(\hat{k})\right]$
(D.4)

where we have defined the kernels

$\displaystyle R_{3}^{L^{(2)}}(k)$ $\displaystyle=-P(k)\int dq\ q\
\Big{[}\frac{12}{175}k^{2}\xi^{1}_{-1}(q)j_{1}(kq)-\frac{12}{175}k^{2}\xi^{3}_{-1}(q)j_{3}(kq)+\frac{8}{175}k\xi^{0}_{0}(q)j_{0}(kq)j_{2}(kq)$
$\displaystyle-\frac{8}{49}k\xi^{2}_{0}(q)+\frac{144k\xi^{4}_{0}(q)j_{4}(kq)}{1225}-\frac{36\xi^{1}_{1}(q)j_{1}(kq)}{1225}+\frac{44}{525}\xi^{3}_{1}(q)j_{3}(kq)-\frac{8}{147}\xi^{5}_{1}(q)j_{5}(kq)\Big{]}$
$\displaystyle R_{1}^{L^{(2)}}(k)=-P(k)\int$ $\displaystyle dq\ q\
\Big{[}-\frac{4}{175}k^{2}\xi^{1}_{-1}(q)j_{1}(kq)+\frac{4}{175}k^{2}\xi^{3}_{-1}(q)j_{3}(kq)+\frac{32}{525}k\xi^{0}_{0}(q)j_{0}(kq)$
$\displaystyle-\frac{8}{147}k\xi^{2}_{0}(q)j_{2}(kq)-\frac{8k\xi^{4}_{0}(q)j_{4}(kq)}{1225}-\frac{4}{175}\xi^{1}_{1}(q)j_{1}(kq)+\frac{4}{175}\xi^{3}_{1}(q)j_{3}(kq)\Big{]}.$

## References

  * [1] K. Akitsu, Y. Li and T. Okumura, _Quadratic shape biases in three-dimensional halo intrinsic alignments_ , _arXiv e-prints_ (2023) arXiv:2306.00969 [2306.00969].
  * [2] DES Collaboration, T.M.C. Abbott, M. Aguena, A. Alarcon, S. Allam, O. Alves et al., _Dark Energy Survey Year 3 Results: Cosmological Constraints from Galaxy Clustering and Weak Lensing_ , _arXiv e-prints_ (2021) arXiv:2105.13549 [2105.13549].
  * [3] C. Heymans, T. Tröster, M. Asgari, C. Blake, H. Hildebrandt, B. Joachimi et al., _KiDS-1000 Cosmology: Multi-probe weak gravitational lensing and spectroscopic galaxy clustering constraints_ , _A &A_ 646 (2021) A140 [2007.15632].
  * [4] H. Miyatake, S. Sugiyama, M. Takada, T. Nishimichi, X. Li, M. Shirasaki et al., _Hyper Suprime-Cam Year 3 Results: Cosmology from Galaxy Clustering and Weak Lensing with HSC and SDSS using the Emulator Based Halo Model_ , _arXiv e-prints_ (2023) arXiv:2304.00704 [2304.00704].
  * [5] M. Kilbinger, _Cosmology with cosmic shear observations: a review_ , _Reports on Progress in Physics_ 78 (2015) 086901.
  * [6] C.M. Hirata and U. Seljak, _Intrinsic alignment-lensing interference as a contaminant of cosmic shear_ , _Phys. Rev. D_ 70 (2004) 063526 [astro-ph/0406275].
  * [7] M.A. Troxel and M. Ishak, _The intrinsic alignment of galaxies and its impact on weak gravitational lensing in an era of precision cosmology_ , _Phys. Rep._ 558 (2015) 1 [1407.6990].
  * [8] E. Krause, T. Eifler and J. Blazek, _The impact of intrinsic alignment on current and future cosmic shear surveys_ , _Monthly Notices of the Royal Astronomical Society_ 456 (2015) 207.
  * [9] R. Mandelbaum, C.M. Hirata, M. Ishak, U. Seljak and J. Brinkmann, _Detection of large-scale intrinsic ellipticity-density correlation from the Sloan Digital Sky Survey and implications for weak lensing surveys_ , _MNRAS_ 367 (2006) 611 [astro-ph/0509026].
  * [10] S. Singh, R. Mandelbaum and S. More, _Intrinsic alignments of SDSS-III BOSS LOWZ sample galaxies_ , _MNRAS_ 450 (2015) 2195 [1411.1755].
  * [11] K. Xu, Y.P. Jing, G.-B. Zhao and A.J. Cuesta, _Evidence for baryon acoustic oscillations from galaxy-ellipticity correlations_ , 2023.
  * [12] T. Okumura and A. Taruya, _First Constraints on Growth Rate from Redshift-space Ellipticity Correlations of SDSS Galaxies at 0.16 < z < 0.70_, _ApJ_ 945 (2023) L30 [2301.06273].
  * [13] F. Schmidt, N.E. Chisari and C. Dvorkin, _Imprint of inflation on galaxy shape correlations_ , _Journal of Cosmology and Astroparticle Physics_ 2015 (2015) 032.
  * [14] K. Akitsu, T. Kurita, T. Nishimichi, M. Takada and S. Tanaka, _Imprint of anisotropic primordial non-Gaussianity on halo intrinsic alignments in simulations_ , _Phys. Rev. D_ 103 (2021) 083508 [2007.03670].
  * [15] K. Kogai, K. Akitsu, F. Schmidt and Y. Urakawa, _Galaxy imaging surveys as spin-sensitive detector for cosmological colliders_ , _J. Cosmology Astropart. Phys_ 2021 (2021) 060 [2009.05517].
  * [16] K. Akitsu, Y. Li and T. Okumura, _Gravitational wave fossils in nonlinear regime: Halo tidal bias and intrinsic alignments from gravitational wave separate universe simulations_ , _Phys. Rev. D_ 107 (2023) 063531 [2209.06226].
  * [17] T. Kurita and M. Takada, _Constraints on anisotropic primordial non-Gaussianity from intrinsic alignments of SDSS-III BOSS galaxies_ , _arXiv e-prints_ (2023) arXiv:2302.02925 [2302.02925].
  * [18] P. Catelan, M. Kamionkowski and R.D. Blandford, _Intrinsic and extrinsic galaxy alignment_ , _Monthly Notices of the Royal Astronomical Society_ 320 (2001) L7.
  * [19] J.A. Blazek, N. MacCrann, M. Troxel and X. Fang, _Beyond linear galaxy alignments_ , _Physical Review D_ 100 (2019) .
  * [20] D. Schmitz, C. Hirata, J. Blazek and E. Krause, _Time evolution of intrinsic alignments of galaxies_ , _Journal of Cosmology and Astroparticle Physics_ 2018 (2018) 030.
  * [21] Z. Vlah, N.E. Chisari and F. Schmidt, _An EFT description of galaxy intrinsic alignments_ , _JCAP_ 01 (2020) 025 [1910.08085].
  * [22] T. Bakx, T. Kurita, N.E. Chisari, Z. Vlah and F. Schmidt, _Effective Field Theory of Intrinsic Alignments at One Loop Order: a Comparison to Dark Matter Simulations_ , _arXiv e-prints_ (2023) arXiv:2303.15565 [2303.15565].
  * [23] V. Desjacques, D. Jeong and F. Schmidt, _Large-scale galaxy bias_ , _Phys. Rep._ 733 (2018) 1 [1611.09787].
  * [24] T. Matsubara, _Resumming cosmological perturbations via the lagrangian picture: One-loop results in real space and in redshift space_ , _Physical Review D_ 77 (2008) .
  * [25] J. Carlson, B. Reid and M. White, _Convolution Lagrangian perturbation theory for biased tracers_ , _MNRAS_ 429 (2013) 1674 [1209.0780].
  * [26] L. Senatore and M. Zaldarriaga, _The IR-resummed Effective Field Theory of Large Scale Structures_ , _J. Cosmology Astropart. Phys_ 2015 (2015) 013 [1404.5954].
  * [27] Z. Vlah, M. White and A. Aviles, _A Lagrangian effective field theory_ , _J. Cosmology Astropart. Phys_ 2015 (2015) 014 [1506.05264].
  * [28] T. Matsubara, _Nonlinear perturbation theory with halo bias and redshift-space distortions via the Lagrangian picture_ , _Phys. Rev. D_ 78 (2008) 083519 [0807.1733].
  * [29] S.-F. Chen, Z. Vlah, E. Castorina and M. White, _Redshift-space distortions in lagrangian perturbation theory_ , _Journal of Cosmology and Astroparticle Physics_ 2021 (2021) 100.
  * [30] C. Modi, S.-F. Chen and M. White, _Simulations and symmetries_ , _Monthly Notices of the Royal Astronomical Society_ 492 (2020) 5754.
  * [31] N. Kokron, J. DeRose, S.-F. Chen, M. White and R.H. Wechsler, _The cosmology dependence of galaxy clustering and lensing from a hybrid n-body–perturbation theory model_ , _Monthly Notices of the Royal Astronomical Society_ 505 (2021) 1422.
  * [32] B. Hadzhiyska, C. Garcí a-García, D. Alonso, A. Nicola and A. Slosar, _Hefty enhancement of cosmological constraints from the DES y1 data using a hybrid effective field theory approach to galaxy bias_ , _Journal of Cosmology and Astroparticle Physics_ 2021 (2021) 020.
  * [33] J. DeRose, N. Kokron, A. Banerjee, S.-F. Chen, M. White, R. Wechsler et al., _Aemulus $\nu$: precise predictions for matter and biased tracer power spectra in the presence of neutrinos_, _Journal of Cosmology and Astroparticle Physics_ 2023 (2023) 054.
  * [34] T. Matsubara, _The integrated perturbation theory for cosmological tensor fields I: Basic formulation_ , _arXiv e-prints_ (2022) arXiv:2210.10435 [2210.10435].
  * [35] T. Matsubara, _The integrated perturbation theory for cosmological tensor fields II: Loop corrections_ , _arXiv e-prints_ (2022) arXiv:2210.11085 [2210.11085].
  * [36] T. Matsubara, _The integrated perturbation theory for cosmological tensor fields III: Projection effects_ , _arXiv e-prints_ (2023) arXiv:2304.13304 [2304.13304].
  * [37] S.-F. Chen, Z. Vlah and M. White, _Consistent modeling of velocity statistics and redshift-space distortions in one-loop perturbation theory_ , _J. Cosmology Astropart. Phys_ 2020 (2020) 062 [2005.00523].
  * [38] A. Guth, _Traceless symmetric tensor approach to legendre polynomials and spherical harmonics_ , 2012.
  * [39] M. Schmittfull, Z. Vlah and P. McDonald, _Fast large scale structure perturbation theory using one-dimensional fast fourier transforms_ , _Physical Review D_ 93 (2016) .
  * [40] P. McDonald and A. Roy, _Clustering of dark matter tracers: generalizing bias for the coming era of precision LSS_ , _J. Cosmology Astropart. Phys_ 2009 (2009) 020 [0902.0991].
  * [41] Z. Vlah, U. Seljak, M. Yat Chu and Y. Feng, _Perturbation theory, effective field theory, and oscillations in the power spectrum_ , _J. Cosmology Astropart. Phys_ 2016 (2016) 057 [1509.02120].
  * [42] J.J. Sakurai, _Modern quantum mechanics; rev. ed._ , Addison-Wesley, Reading, MA (1994).
  * [43] T. Baldauf, U. Seljak, R.E. Smith, N. Hamaus and V. Desjacques, _Halo stochasticity from exclusion and nonlinear clustering_ , _Physical Review D_ 88 (2013) .
  * [44] Z. Vlah, E. Castorina and M. White, _The gaussian streaming model and convolution lagrangian effective field theory_ , _Journal of Cosmology and Astroparticle Physics_ 2016 (2016) 007.
  * [45] A.J.S. Hamilton, _Uncorrelated modes of the non-linear power spectrum_ , _MNRAS_ 312 (2000) 257 [astro-ph/9905191].
  * [46] N. Kokron, S.-F. Chen, M. White, J. DeRose and M. Maus, _Accurate predictions from small boxes: variance suppression via the Zel’dovich approximation_ , _J. Cosmology Astropart. Phys_ 2022 (2022) 059 [2205.15327].
  * [47] M. White, _The zel'dovich approximation_ , _Monthly Notices of the Royal Astronomical Society_ 439 (2014) 3630.
  * [48] M. Zaldarriaga and U. Seljak, _Gravitational lensing effect on cosmic microwave background polarization_ , _Phys. Rev. D_ 58 (1998) 023003 [astro-ph/9803150].
  * [49] R. Takahashi, M. Sato, T. Nishimichi, A. Taruya and M. Oguri, _REVISING THE HALOFIT MODEL FOR THE NONLINEAR MATTER POWER SPECTRUM_ , _The Astrophysical Journal_ 761 (2012) 152.
  * [50] M. Schmittfull, M. Simonović, M.M. Ivanov, O.H.E. Philcox and M. Zaldarriaga, _Modeling galaxies in redshift space at the field level_ , _J. Cosmology Astropart. Phys_ 2021 (2021) 059 [2012.03334].
  * [51] L. Secco et al., _Dark energy survey year 3 results: Cosmology from cosmic shear and robustness to modeling uncertainty_ , _Physical Review D_ 105 (2022) .
  * [52] F. Maion, R.E. Angulo, T. Bakx, N.E. Chisari, T. Kurita and M. Pellejero-Ibáñez, _HYMALAIA: A Hybrid Lagrangian Model for Intrinsic Alignments_ , _arXiv e-prints_ (2023) arXiv:2307.13754 [2307.13754].
  * [53] N. Kokron, S.-F. Chen, M. White, J. DeRose and M. Maus, _Accurate predictions from small boxes: variance suppression via the zel'dovich approximation_ , _Journal of Cosmology and Astroparticle Physics_ 2022 (2022) 059.
  * [54] J. DeRose, S.-F. Chen, N. Kokron and M. White, _Precision redshift-space galaxy power spectra using zel'dovich control variates_ , _Journal of Cosmology and Astroparticle Physics_ 2023 (2023) 008.
  * [55] B. Hadzhiyska, M.J. White, X. Chen, L.H. Garrison, J. DeRose, N. Padmanabhan et al., _Mitigating the noise of desi mocks using analytic control variates_ , 2023.
  * [56] Z. Vlah and M. White, _Exploring redshift-space distortions in large-scale structure_ , _J. Cosmology Astropart. Phys_ 2019 (2019) 007 [1812.02775].
  * [57] U. Seljak and P. McDonald, _Distribution function approach to redshift space distortions_ , _J. Cosmology Astropart. Phys_ 2011 (2011) 039 [1109.1888].