# The integrated perturbation theory for cosmological tensor fields I: Basic
formulation

Takahiko Matsubara<EMAIL_ADDRESS>Institute of Particle and Nuclear
Studies, High Energy Accelerator Research Organization (KEK), Oho 1-1, Tsukuba
305-0801, Japan  The Graduate University for Advanced Studies (SOKENDAI),
Tsukuba 305-0801, Japan

###### Abstract

In order to extract maximal information about cosmology from the large-scale
structure of the Universe, one needs to use every bit of signals that can be
observed. Beyond the spatial distributions of astronomical objects, the
spatial correlations of tensor fields, such as galaxy spins and shapes, are
ones of promising sources that one can access in the era of large surveys in
near future. The perturbation theory is a powerful tool to analytically
describe the behaviors and evolutions of correlation statistics on large
scales for a given cosmology. In this paper, we formulate a nonlinear
perturbation theory of tensor fields in general, based on the formulation of
integrated perturbation theory for the scalar-valued bias, generalizing it to
include the tensor-valued bias. To take advantage of rotational symmetry, the
formalism is constructed on the basis of irreducible decomposition of tensors,
identifying physical variables which are invariant under the rotation of the
coordinates system.

##  I Introduction

The large-scale structure of the Universe is an inexhaustible fountain of the
cosmological information and has been playing important roles in cosmology.
The various patterns in the spatial distribution of galaxies have emerged from
initial density fluctuations of the Universe, which origin is believed to be
generated by quantum fluctuations in the inflationary phase of the very early
Universe [1]. Considerable information on the evolution of the Universe is
also contained in the large-scale structure. For example, the spatial
clustering of galaxies is a promising probe to reveal the nature of dark
energy [2] and the initial condition in the early Universe such as the
primordial spectrum and non-Gaussianity etc. [3]. While the large-scale
structure is used as a probe of the density field in the relatively late
Universe, the distribution of observable objects, such as galaxies, quasars,
and 21cm lines, etc., are biased tracers of the underlying matter distribution
of the Universe. The difference between distributions of galaxies and matter
is referred to as “galaxy bias” [4]. The phenomena of biasing are not
restricted to the galaxies, but are unavoidable when we use astronomical
objects as tracers of underlying matter distribution.

The number density of galaxies is the primary target to probe the properties
of large-scale structure. Theoretical predictions for the spatial distribution
of galaxies is investigated by various methods, including numerical
simulations, analytic modeling, and the perturbation theories. On one hand,
the dynamical evolutions of the mass density fields on relatively large scales
are analytically described by the perturbation theory [5, 6, 7, 8], which is
valid as long as the density fluctuations are small and in the quasi-linear
regime. On the other hand, the dynamical structure formation on relatively
small scales, including the formation of galaxies and other astronomical
objects, is not analytically tractable because of the full nonlinearity of the
problem. One should resort to numerical simulations and nonlinear modeling to
describe and understand the small-scale structure formation. For example, the
halo model [9, 10, 11] statistically predicts spatial distributions of halos
depending on the initial conditions, using simple assumptions based on the
Press-Schechter theory [12] and its extensions [13, 14] of the nonlinear
structure formation.

The biasing of tracers affects not only spatial distributions on small scales,
but also those on large scales. However, the biasing effects are not so
complicated on large scales as those on small scales. For example, the density
contrast $\delta_{\mathrm{g}}$ of biased tracers on sufficiently large scales,
where the dynamics are described by linear theory, is proportional to the
linear density contrast $\delta_{\mathrm{L}}$, and we have a simple relation
$\delta_{\mathrm{g}}=b\,\delta_{\mathrm{L}}$, where a proportional constant
$b$ is called the linear bias parameter [15, 16]. While nonlinear dynamics
complicate the situation, the nonlinear perturbation theory to describe the
spatial clustering of biased tracers has been also developed [8, 4]. Since the
perturbation theory cannot describe the nonlinear structure formation which
takes place on small scales, it is inevitable to introduce some unknown
parameters or functions whose values strongly depend on the small-scale
nonlinear physics which cannot be predicted from the first principle. These
nuisance parameters are called bias parameters or bias functions in the
nonlinear perturbation theory. When a model of bias is defined in some way,
one can calculate the evolution of the biased tracers by applying the
nonlinear perturbation theory.

There are several different ways of characterizing the bias model in the
nonlinear perturbation theory [4]. Among others, the integrated perturbation
theory (iPT) [17, 19, 18, 20], which is based on the Lagrangian perturbation
theory [21, 22] and orthogonal decomposition of the bias [23], provides a
method to generically include any model of bias formulated in Lagrangian
space, irrespective of whether the bias is local or nonlocal. In this
formalism, the concept of the renormalized bias functions [17] is introduced.
Dynamical evolutions on large scales which are described by perturbation
theory are separated from complicated nonlinear processes of structure
formation on small scales. The latter complications are all encoded in the
renormalized bias functions in the formalism of iPT. The renormalized bias
functions are given for any models of bias, whether they are local or nonlocal
in general. For example, the halo bias is one of the typical models of
Lagrangian bias, and the renormalized bias functions are uniquely determined
and calculated from the model [19].

The number density of biased objects is not the only quantity we can observe.
For example, the position and shape of galaxies are simultaneously observed in
imaging surveys of galaxies, which are essential in observations of weak
lensing fields to probe the nature of dark matter and dark energy [24]. Before
the lensing effects, the galaxy shapes are more or less correlated to the mass
density field through, e.g., gravitational tidal fields and other
environmental structures in the underlying density field. The shapes of
galaxies, i.e., sizes and intrinsic alignments [25, 26, 27, 28, 29, 30, 31,
32, 33] are expected to be promising probes of cosmological information [34]
in near future when unprecedentedly large imaging surveys are going to take
place. For example, the shape statistics of galaxies offer a new probe of
particular features in primordial non-Gaussianity generated during the
inflation in the presence of higher-spin fields [35], which features are
dubbed “cosmological collider physics.”

Motivated by recent progress in observational techniques of imaging surveys,
analytical modelings of galaxy shape statistics by the nonlinear perturbation
theory have recently been considered [36, 37, 38, 39, 40, 41]. Among others,
Ref. [39] developed a formalism based on spherical tensors to describe the
galaxy shapes, with which rotation and parity symmetries in the shape
statistics are respected. They use an approach of what is called Effective
Field Theory of the Large-Scale Structure (EFTofLSS) with the Eulerian
perturbation theory, and galaxy biasing are modeled by introducing a set of
semi-local terms that are generally allowed by conceivable symmetries of the
problem and the coefficients of these terms are free parameters in the theory.

The galaxy shape statistics are usually characterized by the second-order
moments of galaxy images in the two-dimensional sky, which correspond to the
projection of the three-dimensional second-order moments of galaxy shapes to
the two-dimensional sky. In principle, the galaxy images contain information
about other moments higher than the second order. The higher-order moments in
galaxy images turn out to be possible probes for the higher-spin fields in the
context of cosmological collider physics through angle-dependent primordial
non-Gaussianity. Ref. [42] showed that the shape statistics of the order-$s$
moment turn out to be a probe for spin-$s$ fields.

In this paper, inspired by the above developments, we consider a
generalization of the iPT formalism to predict correlations of tensor fields
of rank $l$ in general with the perturbation theory. In order to have
rotational symmetry apparent, we decompose the tensor fields into irreducible
tensors on the basis of spherical tensors, just similarly in the pioneering
work of Refs. [39, 40]. However, unlike the previous work, we do not fix the
coordinates system to the directions of wavevectors of Fourier modes, so as to
make the theory fully covariant with respect to the rotational symmetry in
three-dimensional space. It turns out that the generalizing iPT formalism on
the basis of irreducible tensors results in an elegant way of describing the
perturbation theory only with rotationally invariant quantities in the
presence of tensor-valued bias. Various methods developed for the theory of
angular momentum in quantum mechanics, such as the $3nj$-symbols and their
orthogonality relations, sum rules, the Wigner-Eckart theorem, and so on, are
effectively utilized in the course of calculations for predicting statistical
quantities of tensor fields.

This paper defines a basic formalism of the tensor generalization of the
methods of iPT. Several applications of the fully general formalism to
relatively simple examples of calculating statistics are also illustrated,
i.e., the lowest-order predictions of the power spectrum of tensor fields in
real space and in redshift space, lowest-order effects of primordial non-
Gaussianity in the power spectrum, and tree-level predictions of the
bispectrum of tensor fields in real space. The main purpose of this paper is
to show the fundamental formulation of tensor-valued iPT, and to give useful
methods and formulas for future applications to calculate statistics of
concrete tensor fields, such as galaxy spins, shapes, etc. Further
developments for calculating loop corrections in the higher-order perturbation
theory are described in a separate paper in the series, Paper II [43], which
follows this paper. Since the three-dimensional tensor fields themselves are
not directly observable, instead, the projected two-dimensional tensors on the
sky are observable. The projection effects for predicting observable tensor
components are studied in another separate paper in the series, Paper III
[44].

This paper is organized as follows. In Sec. II, various notations regarding
the spherical basis are introduced and defined, and their fundamental
properties are derived and summarized. In Sec. III, the iPT formalism
originally for scalar-valued fields is generalized to those for tensor-valued
fields, and rotationally invariant components of renormalized bias functions
and higher-order propagators are identified. Additional symmetries in them are
also described. Explicit forms of several lower-order propagators are
calculated and presented. In Sec. IV, the calculation of the power spectrum,
correlation function and bispectrum of tensor fields by our formalism is
illustrated and presented in several cases which are not relatively
complicated. In Sec. V, we define a class of bias models, which we call semi-
local models, in which the Lagrangian bias of the tensor field is determined
only by derivatives of the gravitational potential at the same position in
Lagrangian space. Conclusions are given in Sec. VI. In Appendix A, explicit
expressions of higher-order spherical tensor basis are derived and given.
Formulas regarding spherical harmonics which are repeatedly used throughout
this paper and subsequent papers in the series are summarized in Appendix B
for readers’ convenience. In Appendix C is summarized some properties of
Wigner’s $3j$-, $6j$\- and $9j$-symbols which are useful throughout the series
of papers.

##  II Decomposition of tensors by spherical basis

In this paper, we extensively use the irreducible representations of the
tensor field in spherical basis. We employ a set of unusual conventions for
representations of rotation group, which are sometimes different from those in
standard literature. The reason is that the standard conventions, which are
optimized for applications to the quantum theory of angular momentum, turn out
to be inconvenient for our applications to the perturbation theory and
statistics of large-scale structure. While the contents of this section are
not new, the purpose is to introduce and define our notations and conventions
of various quantities which are repeatedly used in the series of papers. They
are sometimes different from conventional notations found in standard
literature, as we find the conventional ones are quite inconvenient for our
purpose and make many equations unnecessarily tedious and complicated.

For a given set of orthonormal basis vectors $\hat{\mathbf{e}}_{i}$
($i=1,2,3$) in a flat, three-dimensional space, the spherical basis
($\mathbf{e}_{0}$, $\mathbf{e}_{\pm}$) is defined by [45, 46]

$\mathbf{e}_{0}=\hat{\mathbf{e}}_{3},\quad\mathbf{e}_{\pm}=\mp\frac{\hat{\mathbf{e}}_{1}\pm
i\hat{\mathbf{e}}_{2}}{\sqrt{2}}.$ (1)

It is convenient to introduce the dual basis with an upper index by taking
complex conjugates of the spherical basis:

$\mathbf{e}^{0}\equiv\mathbf{e}_{0}^{*}=\hat{\mathbf{e}}_{3},\quad\mathbf{e}^{\pm}\equiv\mathbf{e}_{\pm}^{*}=\mp\frac{\hat{\mathbf{e}}_{1}\mp
i\hat{\mathbf{e}}_{2}}{\sqrt{2}}.$ (2)

The orthogonality relations among these bases are given by

$\mathbf{e}_{m}\cdot\mathbf{e}_{m^{\prime}}=g^{(1)}_{mm^{\prime}},\quad\mathbf{e}_{m}\cdot\mathbf{e}^{m^{\prime}}=\delta_{m}^{m^{\prime}},\quad\mathbf{e}^{m}\cdot\mathbf{e}^{m^{\prime}}=g_{(1)}^{mm^{\prime}},$
(3)

where $m,m^{\prime}$ are azimuthal indices of the spherical basis which
represent one of $(-,0,+)$. The spherical metric tensors are defined by

$g^{(1)}_{mm^{\prime}}=g_{(1)}^{mm^{\prime}}=(-1)^{m}\delta_{m,-m^{\prime}},$
(4)

where $\delta_{m}^{m^{\prime}}=\delta_{mm^{\prime}}$ is the Kronecker’s delta
symbol which is unity for $m=m^{\prime}$ and zero for $m\neq m^{\prime}$. The
matrices $g_{(1)}^{mm^{\prime}}$ and $g^{(1)}_{mm^{\prime}}$ are inverse
matrices of each other,

$g_{(1)}^{mm^{\prime\prime}}g^{(1)}_{m^{\prime\prime}m^{\prime}}=\delta^{m}_{m^{\prime}},$
(5)

where the repeated index $m^{\prime\prime}$ is assumed to be summed over,
omitting summation symbols as commonly adopted in the convention of tensor
analysis (Einstein summation convention). Throughout this paper, the Einstein
convention of summation is always assumed, unless otherwise stated. They act
as metric tensors of the spherical basis,

$\mathbf{e}^{m}=g_{(1)}^{mm^{\prime}}\mathbf{e}_{m^{\prime}},\quad\mathbf{e}_{m}=g^{(1)}_{mm^{\prime}}\mathbf{e}^{m^{\prime}}.$
(6)

Any vector can be represented by $\bm{V}\propto
V^{m}\mathbf{e}_{m}=V_{m}\mathbf{e}^{m}$ up to normalization constant, and
thus $\mathbf{e}_{m}$ is considered as the basis of contravariant spherical
vectors, and $\mathbf{e}^{m}$ is considered as the basis of covariant
spherical vectors.

One of the advantages of the spherical basis over the Cartesian basis is that
the Cartesian tensors are reduced into irreducible tensors on the spherical
basis under a rotation of coordinates system [47]. We define the spherical
basis of traceless, symmetric tensors [48, 39] as111 The normalization of the
spherical tensors in this paper are somehow different from those in Refs. [48,
39]. Our normalization is chosen so that Eq. (11) below should hold. In
addition, in the previous literature, the axis $\mathbf{e}^{0}$ of the
coordinates system is specifically chosen to be parallel to the separation
vector $\bm{r}$ of the correlation function [48], or, to the wavevector
$\bm{k}$ of Fourier modes [39], while our axis can be arbitrarily chosen, so
that fully rotational symmetry is explicitly kept in our formalism, not
choosing a special coordinates system.

$\mathsf{Y}^{(0)}=1,\quad\mathsf{Y}^{(m)}_{i}={\mathrm{e}^{m}}_{i}$ (7)

for tensors of rank 0 and 1, and

$\displaystyle\mathsf{Y}^{(0)}_{ij}=\sqrt{\frac{3}{2}}\left({\mathrm{e}^{0}}_{i}\,{\mathrm{e}^{0}}_{j}-\frac{1}{3}\delta_{ij}\right),$
(8) $\displaystyle\mathsf{Y}^{(\pm
1)}_{ij}=\sqrt{2}\,{\mathrm{e}^{0}}_{(i}\,{\mathrm{e}^{\pm}}_{j)},$ (9)
$\displaystyle\mathsf{Y}^{(\pm
2)}_{ij}={\mathrm{e}^{\pm}}_{i}\,{\mathrm{e}^{\pm}}_{j}$ (10)

for tensors of rank 2, where ${\mathrm{e}^{m}}_{i}=[\mathbf{e}^{m}]_{i}$ is
the Cartesian components of the spherical basis, and round brackets in the
indices of the right-hand side (rhs) of Eq. (9) indicate symmetrization with
respect to the indices inside the brackets, e.g.,
$x_{(i}y_{j)}=(x_{i}y_{j}+x_{j}y_{i})/2$. These bases satisfy orthogonality
relations,

$\mathsf{Y}^{(0)}\mathsf{Y}^{(0)}=1,\quad\mathsf{Y}^{(m)}_{i}\mathsf{Y}^{(m^{\prime})}_{i}=g_{(1)}^{mm^{\prime}},\quad\mathsf{Y}^{(m)}_{ij}\mathsf{Y}^{(m^{\prime})}_{ij}=g_{(2)}^{mm^{\prime}},$
(11)

where, in the last equation,
$g_{(2)}^{mm^{\prime}}=(-1)^{m}\delta_{m,-m^{\prime}}$ is the $5\times 5$
spherical metric, and the indices $m,m^{\prime}$ run over the integers
$-2,-1,0,+1,+2$ in this case. The Einstein summation convention is also
applied for the Cartesian indices in the above, so that $i$ and/or $j$ are
summed over on the left-hand side (lhs) of the above equations. Taking complex
conjugates of the spherical basis virtually lowers the azimuthal indices,

$\mathsf{Y}^{(0)*}=\mathsf{Y}^{(0)},\quad\mathsf{Y}^{(m)*}_{i}=g^{(1)}_{mm^{\prime}}\mathsf{Y}^{(m^{\prime})}_{i},\quad\mathsf{Y}^{(m)*}_{ij}=g^{(2)}_{mm^{\prime}}\mathsf{Y}^{(m^{\prime})}_{ij}.$
(12)

For a unit vector $\bm{n}=(\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta)$
whose direction is given by spherical coordinates $(\theta,\phi)$, the
spherical tensors are related to the spherical harmonics $Y_{lm}(\theta,\phi)$
as

$\displaystyle Y_{00}(\theta,\phi)$
$\displaystyle=\frac{1}{\sqrt{4\pi}}\mathsf{Y}^{(0)*},\quad
Y_{1m}(\theta,\phi)=\sqrt{\frac{3}{4\pi}}\,\mathsf{Y}^{(m)*}_{i}n_{i},$
$\displaystyle Y_{2m}(\theta,\phi)$
$\displaystyle=\sqrt{\frac{15}{8\pi}}\,\mathsf{Y}^{(m)*}_{ij}n_{i}n_{j}.$ (13)

where the normalization convention of the spherical harmonics are given by a
standard one,

$Y_{lm}(\theta,\phi)=\sqrt{\frac{2l+1}{4\pi}}\sqrt{\frac{(l-m)!}{(l+m)!}}P_{l}^{m}(\cos\theta)\,e^{im\phi},$
(14)

and $P_{l}^{m}(x)$ are the associated Legendre polynomials. The Condon-
Shortley phase is included in the latter, namely,

$P_{l}^{m}(x)=\frac{(-1)^{m}}{2^{l}\,l!}\left(1-x^{2}\right)^{m/2}\frac{d^{l+m}}{dx^{l+m}}\left(1-x^{2}\right)^{l}.$
(15)

Similarly, spherical bases for higher-rank symmetric tensors are uniquely
defined by imposing the similar properties of Eq. (13),

$Y_{lm}(\theta,\phi)=\sqrt{\frac{(2l+1)!!}{4\pi\,l!}}\,\mathsf{Y}^{(m)*}_{i_{1}i_{2}\cdots
i_{l}}n_{i_{1}}n_{i_{2}}\cdots n_{i_{l}}.$ (16)

General procedures to construct spherical bases of arbitrary rank are given
and explicit forms of the spherical bases of the rank 0 to 4 are derived in
Appendix A.

The standard normalization of spherical harmonics turns out to be not so
convenient to represent various equations in this work. Instead, the spherical
harmonics with Racah’s normalization, defined by

$C_{lm}(\theta,\phi)\equiv\sqrt{\frac{4\pi}{2l+1}}Y_{lm}(\theta,\phi)=\sqrt{\frac{(l-m)!}{(l+m)!}}P_{l}^{m}(\cos\theta)\,e^{im\phi},$
(17)

turn out to simplify various equations. One of the reason is that the last
normalization is closely related to the standard normalization of Legendre
polynomials, as we have $C_{l0}(\theta,\phi)=P_{l}(\cos\theta)$, where
$P_{l}(x)=P_{l}^{0}(x)$ are the ordinary Legendre polynomials. With the
spherical harmonics with Racah’s normalization, a general relation of Eq. (16)
is given by

$\mathsf{Y}^{(m)}_{i_{1}i_{2}\cdots i_{l}}n_{i_{1}}n_{i_{2}}\cdots
n_{i_{l}}=A_{l}\,C_{lm}^{*}(\theta,\phi),$ (18)

where we define a factor,

$A_{l}\equiv\sqrt{\frac{l!}{(2l-1)!!}},$ (19)

which repeatedly appear in later applications.

The spherical bases are irreducible representations of the rotation group
SO(3) in the same way as the spherical harmonics are. We consider a passive
rotation of the Cartesian basis,

$\hat{\mathbf{e}}_{i}\xrightarrow{\mathbb{R}}\hat{\mathbf{e}}_{i}^{\prime}=R_{ij}\hat{\mathbf{e}}_{j},$
(20)

where $R_{ij}=\hat{\mathbf{e}}_{i}^{\prime}\cdot\hat{\mathbf{e}}_{j}$ are
components of a real orthogonal matrix which satisfies $R^{\mathrm{T}}R=I$.
Cartesian components of a unit vector $\bm{n}$, where $|\bm{n}|=1$, transform
as

$n_{i}\xrightarrow{\mathbb{R}}n_{i}^{\prime}=(R^{-1})_{ij}n_{j}=n_{j}R_{ji},$
(21)

because we have
$n_{i}\hat{\mathbf{e}}_{i}=n_{i}^{\prime}\hat{\mathbf{e}}_{i}^{\prime}$.
Denoting the spherical coordinates of $n_{i}$ and $n_{i}^{\prime}$ as
$(\theta,\phi)$ and $(\theta^{\prime},\phi^{\prime})$, respectively, the
spherical harmonics transform as

$C_{lm}(\theta^{\prime},\phi^{\prime})=C_{lm^{\prime}}(\theta,\phi)D_{(l)m}^{m^{\prime}}(R),$
(22)

where $D_{(l)m}^{m^{\prime}}(R)$ is the Wigner’s rotation matrix of the
passive rotation $R$. Therefore, the construction of the tensor bases, Eq.
(16), indicates that they transform as

$\displaystyle\mathsf{Y}^{(m)}_{i_{1}i_{2}\cdots
i_{l}}\xrightarrow{\mathbb{R}}\mathsf{Y}^{(m)\prime}_{i_{1}i_{2}\cdots i_{l}}$
$\displaystyle=R_{i_{1}j_{1}}R_{i_{2}j_{2}}\cdots
R_{i_{l}j_{l}}\mathsf{Y}^{(m)}_{j_{1}j_{2}\cdots j_{l}}$
$\displaystyle=D^{m}_{(l)m^{\prime}}(R^{-1})\mathsf{Y}^{(m^{\prime})}_{i_{1}i_{2}\cdots
i_{l}},$ (23)

where Wigner’s matrix of the inverse rotation is given by

$D^{m}_{(l)m^{\prime}}(R^{-1})=D^{m^{\prime}*}_{(l)m}(R).$ (24)

The orthogonality relations of Eq. (11) are generalized for higher-rank
tensors, and are given by

$\mathsf{Y}^{(m)}_{i_{1}i_{2}\cdots
i_{l}}\mathsf{Y}^{(m^{\prime})}_{i_{1}i_{2}\cdots
i_{l}}=g_{(l)}^{mm^{\prime}},$ (25)

where $(2l+1)\times(2l+1)$ matrix $g_{(l)}^{mm^{\prime}}$ is defined just as
Eq. (4) but the indices $m,m^{\prime}$ run over integers from $-l$ to $+l$.
The inverse matrix is similarly denoted by $g^{(l)}_{mm^{\prime}}$, which
elements are the same as $g_{(l)}^{mm^{\prime}}$. Thus, similar relations as
in Eqs. (4) and (5) hold for $g_{(l)}^{mm^{\prime}}$ and
$g^{(l)}_{mm^{\prime}}$:

$g_{(l)}^{mm^{\prime}}=g^{(l)}_{mm^{\prime}}=(-1)^{m}\delta_{m,-m^{\prime}},\quad
g_{(l)}^{mm^{\prime\prime}}g^{(l)}_{m^{\prime\prime}m^{\prime}}=\delta^{m}_{m^{\prime}}.$
(26)

The relations to the $1jm$-symbol introduced by Wigner [49] are given by

$(-1)^{l}g^{(l)}_{mm^{\prime}}=\begin{pmatrix}l\\\ m\ \
m^{\prime}\end{pmatrix}=(-1)^{l}g_{(l)}^{mm^{\prime}}.$ (27)

Taking the complex conjugate of the basis virtually lowers the azimuthal index
as in Eq. (12), i.e., we have

$\mathsf{Y}^{(m)*}_{i_{1}i_{2}\cdots
i_{l}}=g^{(l)}_{mm^{\prime}}\mathsf{Y}^{(m^{\prime})}_{i_{1}i_{2}\cdots
i_{l}},\quad\mathsf{Y}^{(m)}_{i_{1}i_{2}\cdots
i_{l}}=g_{(l)}^{mm^{\prime}}\mathsf{Y}^{(m^{\prime})*}_{i_{1}i_{2}\cdots
i_{l}}.$ (28)

The same properties also apply to spherical harmonics:

$C_{lm}^{*}(\theta,\phi)=g_{(l)}^{mm^{\prime}}C_{lm^{\prime}}(\theta,\phi),\quad
C_{lm}(\theta,\phi)=g^{(l)}_{mm^{\prime}}C_{lm^{\prime}}^{*}(\theta,\phi).$
(29)

Consequently, the orthogonality relations of Eq. (25) is also equivalent to

$\mathsf{Y}^{(m)*}_{i_{1}i_{2}\cdots
i_{l}}\mathsf{Y}^{(m^{\prime})}_{i_{1}i_{2}\cdots
i_{l}}=\delta_{m}^{m^{\prime}},\quad\mathsf{Y}^{(m)*}_{i_{1}i_{2}\cdots
i_{l}}\mathsf{Y}^{(m^{\prime})*}_{i_{1}i_{2}\cdots
i_{l}}=g^{(l)}_{mm^{\prime}}.$ (30)

We decompose a scalar function $S$ and a vector function $V_{i}$ into
spherical components as

$S=S_{00}\mathsf{Y}^{(0)},\quad V_{i}=iV_{1m}\mathsf{Y}^{(m)}_{i}.$ (31)

The imaginary unit $i$ in front of the vector decomposition is a normalization
convention that turn out to be convenient for our purpose below. Due to the
orthogonality relations of Eq. (11), the inverted relations of the above are
given by

$S_{00}=S\mathsf{Y}^{(0)*},\quad V_{1m}=-iV_{i}\mathsf{Y}^{(m)*}_{i}.$ (32)

Since the spherical basis is a traceless tensor, a symmetric rank-2 tensor
function $T_{ij}$, which is not necessarily traceless, is decomposed into
spherical components of the trace part and traceless part as

$T_{ij}=\frac{1}{3}\delta_{ij}T_{00}\mathsf{Y}^{(0)}-\sqrt{\frac{2}{3}}\,T_{2m}\mathsf{Y}^{(m)}_{ij},$
(33)

where

$T_{00}=T_{ii}\mathsf{Y}^{(0)*},\quad
T_{2m}=-\sqrt{\frac{3}{2}}\,T_{ij}\mathsf{Y}^{(m)*}_{ij}.$ (34)

The prefactor $-\sqrt{2/3}$ in front of the second term on the rhs of Eq. (33)
is again our normalization convention. The general rules of our normalization
convention is shortly explained at the end of this section.

Higher-rank tensors can be similarly decomposed into spherical tensors. We
first note that any symmetric tensor of rank $l$ can be decomposed into [50]

$T_{i_{1}i_{2}\cdots i_{l}}=T^{(l)}_{i_{1}i_{2}\cdots
i_{l}}+\frac{l(l-1)}{2(2l-1)}\delta_{(i_{1}i_{2}}T^{(l-2)}_{i_{3}\cdots
i_{l})}+\cdots,$ (35)

where $T^{(l)}_{i_{1}i_{2}\cdots i_{l}}$ is a symmetric traceless tensor of
rank $l$, $T^{(l-2)}_{i_{3}\cdots i_{l}}$ is traceless part of
$T_{jji_{3}i_{4}\cdots i_{l}}$, and so forth (one should note that the
Einstein summation convention is applied also to Cartesian indices, and thus
$j=1,2,3$ is summed over in the last tensor). The traceless part of
$T_{i_{1}\cdots i_{l}}$ is generally given by a formula [50, 51]

$T^{(l)}_{i_{1}\cdots
i_{l}}=\frac{l!}{(2l-1)!!}\sum_{k=0}^{[l/2]}\frac{(-1)^{k}(2l-2k-1)!!}{2^{k}k!(l-2k)!}\\\
\times\delta_{(i_{1}i_{2}}\cdots\delta_{i_{2k-1}i_{2k}}T^{(l:k)}_{i_{2k+1}\cdots
i_{l})},$ (36)

where $[l/2]$ is the Gauss symbol, i.e., $[l/2]=l/2$ if $l$ is even and
$[l/2]=(l-1)/2$ if $l$ is odd, and

$\displaystyle T^{(l:k)}_{i_{2k+1}\cdots i_{l}}$
$\displaystyle=\delta_{i_{1}i_{2}}\delta_{i_{3}i_{4}}\cdots\delta_{i_{2k-1}i_{2k}}T_{i_{1}\cdots
i_{2k}i_{2k+1}\cdots i_{l}}$ $\displaystyle=T_{j_{1}j_{1}j_{2}j_{2}\cdots
j_{k}j_{k}i_{2k+1}\cdots i_{l}}$ (37)

is the traces of the original tensor taken $k$ times. Tensors of rank 0 and 1
correspond to a scalar and a vector, respectively, and are already traceless
as they do not have any trace component from the first place. The traceless
part of the tensor components forms an irreducible representation of the
rotation group SO(3). Decomposition of higher-rank tensors can be similarly
obtained according to the procedure described in Ref. [50]. Recursively using
Eq. (36), one obtains explicit expressions of Eq. (35). Up to the fourth rank,
they are given by

$\displaystyle T_{ij}$
$\displaystyle=T^{(2)}_{ij}+\frac{1}{3}\delta_{ij}T^{(0)},$ (38)
$\displaystyle T_{ijk}$
$\displaystyle=T^{(3)}_{ijk}+\frac{3}{5}\delta_{(ij}T^{(1)}_{k)},$ (39)
$\displaystyle T_{ijkl}$
$\displaystyle=T^{(4)}_{ijkl}+\frac{6}{7}\delta_{(ij}T^{(2)}_{kl)}+\frac{1}{5}\delta_{(ij}\delta_{kl)}T^{(0)}.$
(40)

The traceless part of a tensor is decomposed by spherical bases as

$\displaystyle T^{(l)}_{i_{1}i_{2}\cdots i_{l}}$
$\displaystyle=i^{l}A_{l}T_{lm}\mathsf{Y}^{(m)}_{i_{1}i_{2}\cdots i_{l}},$
(41) $\displaystyle T_{lm}$
$\displaystyle=\frac{(-i)^{l}}{A_{l}}T^{(l)}_{i_{1}i_{2}\cdots
i_{l}}\mathsf{Y}^{(m)*}_{i_{1}i_{2}\cdots i_{l}}.$ (42)

The prefactor $i^{l}A_{l}$ on the rhs of Eq. (41) is our normalization
convention, where $A_{l}$ is given by Eq. (19). The merits of this particular
choice of normalization for the irreducible tensor $T_{lm}$ in general appear
in Sec. V below. However, the equations presented in Secs. III–IV of this
paper are not affected by the choice of normalization $A_{l}$, because this
normalization factor simultaneously appears on both sides and cancels to each
other in most of the equations.

The traceless part of a traced tensor is similarly decomposed as

$\displaystyle T^{(l-2k)}_{i_{2k+1}\cdots i_{l}}$
$\displaystyle=i^{l-2k}A_{l-2k}T_{l-2k,m}\mathsf{Y}^{(m)}_{i_{2k+1}\cdots
i_{l}},$ (43) $\displaystyle T_{l-2k,m}$
$\displaystyle=\frac{(-i)^{l-2k}}{A_{l-2k}}T^{(l)}_{i_{2k+1}\cdots
i_{l}}\mathsf{Y}^{(m)*}_{i_{2k+1}\cdots i_{l}}.$ (44)

Under the passive rotation of the basis, Eq. (20), Cartesian tensors of rank
$l$ transform as

$T^{(l)}_{i_{1}i_{2}\cdots
i_{l}}\xrightarrow{\mathbb{R}}{T^{\prime}}^{(l)}_{i_{1}i_{2}\cdots
i_{l}}=T^{(l)}_{j_{1}j_{2}\cdots j_{l}}R_{j_{1}i_{1}}\cdots R_{j_{l}i_{l}},$
(45)

and the corresponding transformation of spherical tensors is given by

$T_{lm}\xrightarrow{\mathbb{R}}T^{\prime}_{lm}=T_{lm^{\prime}}D^{m^{\prime}}_{(l)m}(R),$
(46)

as derived from the transformation of the spherical base, Eq. (23). The
contravariant components of the spherical tensors are defined by

$T_{l}^{\,m}=g_{(l)}^{mm^{\prime}}T_{lm^{\prime}},$ (47)

and transform under the rotation as

$T_{l}^{\,m}\xrightarrow{\mathbb{R}}{T^{\prime}}_{l}^{\,m}=D^{m}_{(l)m^{\prime}}(R^{-1})T_{l}^{\,m^{\prime}}.$
(48)

##  III  The integrated perturbation theory of tensor fields

###  III.1  Formulating the iPT of tensor fields

The integrated perturbation theory (iPT) is a systematic method of
perturbation theory to describe the quasi-nonlinear evolution of the large-
scale structure. This method is based on the Lagrangian scheme of the biasing
and nonlinear structure formation. The formalism of this method is explicitly
developed to describe the scalar fields, such as the number density fields of
galaxies or other biased objects in the Universe [17, 18]. In this section,
the formalism is generalized to be able to describe quasi-nonlinear evolutions
of objects which generally have tensor values, such as the angular momentum of
galaxies (vector), the second moment of intrinsic alignment of galaxy shapes
(rank-2 tensor), or higher moments of intrinsic alignment (rank-3 or higher
tensors), etc.

The fundamental formalism of iPT to calculate spatial correlations of biased
objects is described in Ref. [17]. Basically we can just follow the formalism
with a generalization of assigning the tensor values to the objects. We denote
the observable objects $X$, where $X$ is the class of objects selected in a
given cosmological surveys, such as a certain type of galaxies in galaxy
surveys. We consider each object of $a\in X$ has a tensor value
$F^{a}_{i_{1}i_{2}\cdots}$ in general, and define the tensor field by

$F_{X\,i_{1}i_{2}\cdots}(\bm{x})=\frac{1}{\bar{n}_{X}}\sum_{a\in
X}F^{a}_{i_{1}i_{2}\cdots}\delta_{\mathrm{D}}^{3}\left(\bm{x}-\bm{x}_{a}\right),$
(49)

where $\bm{x}_{a}$ is the location of a particular object $a$, $\bar{n}_{X}$
is the mean number density of objects $X$, and $\delta_{\mathrm{D}}^{3}$ is
the three-dimensional Dirac’s delta function. It is crucial to note that the
above definition of the tensor field is weighted by the number density of the
objects $X$. In the above, although the field really depends on the time
$F_{X\,i_{1}i_{2}\cdots}(\bm{x},t)$, the argument of time is omitted in the
notation for simplicity. Similarly, every fields in this paper depend on time,
although they are not explicitly written in the arguments.

In general, any type of tensor can be decomposed into symmetric tensors [50],
and therefore we assume that the tensor field $F^{a}_{i_{1}i_{2}\cdots}$ is a
symmetric tensor without loss of generality. We decompose the symmetric tensor
of each object $F^{a}_{i_{1}i_{2}\cdots}$ into irreducible tensors
$F^{a}_{lm}$ according to the scheme described in the previous section.
Accordingly, corresponding components of Eq. (49) are given by

$F_{Xlm}(\bm{x})=\frac{1}{\bar{n}_{X}}\sum_{a\in
X}F^{a}_{lm}\delta_{\mathrm{D}}^{3}\left(\bm{x}-\bm{x}_{a}\right),$ (50)

As this tensor field is a number density-weighted quantity by construction, we
also consider the unweighted field $G_{Xlm}(\bm{x})$ defined by

$F_{Xlm}(\bm{x})=\left[1+\delta_{X}(\bm{x})\right]G_{Xlm}(\bm{x}),$ (51)

where $\delta_{X}(\bm{x})$ is the number density contrast of the objects $X$,
defined by

$1+\delta_{X}(\bm{x})=\frac{1}{\bar{n}_{X}}\sum_{a\in
X}\delta_{\mathrm{D}}^{3}\left(\bm{x}-\bm{x}_{a}\right).$ (52)

Substituting Eq. (52) into Eq. (51) and comparing with Eq. (50), the
unweighted field satisfies

$G_{Xlm}(\bm{x}_{a})=F^{a}_{Xlm}$ (53)

at each position of object $a$, as naturally expected. The function
$G_{Xlm}(\bm{x})$ can arbitrarily be interpolated where $\bm{x}\neq\bm{x}_{a}$
for all $a\in X$, because the field $F_{Xlm}(\bm{x})$ is not supported at
those points.

The iPT utilizes the Lagrangian perturbation theory (LPT) for the dynamical
nonlinear evolutions (see Ref. [52] for notations of LPT employed in this
paper). The Eulerian coordinates $\bm{x}$ and the Lagrangian coordinates
$\bm{q}$ of a mass element are related by

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

where $\bm{\Psi}(\bm{q})$ is the displacement field. The number density of the
objects in Lagrangian space $n^{\mathrm{L}}_{X}(\bm{q})$ is given by the
Eulerian counterpart by a continuity relation,

$n_{X}(\bm{x})d^{3}\\!x=n^{\mathrm{L}}_{X}(\bm{q})d^{3}\\!q,$ (55)

i.e., the Lagrangian number density is defined so that the Eulerian positions
of objects are displaced back into the Lagrangian positions, and the number
density in Lagrangian space also depends on the time of evaluating
$n_{X}(\bm{x})$. The mean number density of objects $\bar{n}_{X}$ is the same
in Eulerian and Lagrangian spaces by construction, thus we have
$1+\delta_{X}(\bm{x})=[1+\delta^{\mathrm{L}}_{X}(\bm{q})]/J(\bm{q})$, where
$J(\bm{q})=|\partial(\bm{x})/\partial(\bm{q})|$ is the Jacobian of the mapping
from Lagrangian space to Eulerian space. The density contrasts between these
spaces are therefore related by

$1+\delta_{X}(\bm{x})=\int
d^{3}\\!q\left[1+\delta^{\mathrm{L}}_{X}(\bm{q})\right]\delta_{\mathrm{D}}^{3}\left[\bm{x}-\bm{q}-\bm{\Psi}(\bm{q})\right].$
(56)

As in the Eulerian space, the tensor field in Lagrangian space is also
naturally defined with the weight of the number density of objects. The
unweighted tensor field $G^{\mathrm{L}}_{Xlm}(\bm{q})$ in Lagrangian space is
simply defined by taking the same value displaced back from the Eulerian
position of Eq. (51),

$G^{\mathrm{L}}_{Xlm}(\bm{q})\equiv G_{Xlm}(\bm{x}),$ (57)

where $\bm{x}$ on the right-hand side (rhs) is just given by Eq. (54) and the
number density-weighted tensor field in Lagrangian space,
$F^{\mathrm{L}}_{lm}(\bm{q})$, is defined by

$F^{\mathrm{L}}_{Xlm}(\bm{q})=\left[1+\delta^{\mathrm{L}}_{X}(\bm{q})\right]G^{\mathrm{L}}_{Xlm}(\bm{q}).$
(58)

As obviously derived from Eqs. (56) and (57), the number density-weighted
tensor fields in Eulerian space and Lagrangian space is related by

$F_{Xlm}(\bm{x})=\int
d^{3}\\!qF^{\mathrm{L}}_{Xlm}(\bm{q})\delta_{\mathrm{D}}^{3}\left[\bm{x}-\bm{q}-\bm{\Psi}(\bm{q})\right].$
(59)

In cosmology, the properties of observable quantities are determined by the
initial condition of density fluctuations in the Universe, and the growing
mode of linear density contrast, $\delta_{\mathrm{L}}$, is a representative of
initial conditions. Therefore, the tensor field $F_{Xlm}$ is considered as a
functional of the linear density field $\delta_{\mathrm{L}}$. The ensemble
average of the $n$th functional derivatives of the evolved field is called
multi-point propagators [23, 53, 54], which are defined in equivalently
various manners. Here we explain a comprehensive definition according to Ref.
[23], starting from a configuration space:

$\Gamma^{(n)}_{Xlm}\left(\bm{x}-\bm{x}_{1},\ldots,\bm{x}-\bm{x}_{n}\right)\equiv\left\langle\frac{\delta^{n}{F}_{Xlm}(\bm{x})}{\delta\delta_{\mathrm{L}}(\bm{x}_{1})\cdots\delta\delta_{\mathrm{L}}(\bm{x}_{n})}\right\rangle.$
(60)

where $\delta_{\mathrm{L}}(\bm{x})$ is the linear density contrast,
$\delta/\delta\delta_{\mathrm{L}}(\bm{x})$ is the functional derivative with
respect to the linear density contrast, and $\langle\cdots\rangle$ represents
the ensemble average of the field. Because of the translational invariance of
the Universe in a statistical sense, the rhs should be invariant under the
homogeneous translation, $\bm{x}\rightarrow\bm{x}+\bm{x}_{0}$ and
$\bm{x}_{i}\rightarrow\bm{x}_{i}+\bm{x}_{0}$ for $i=1,\ldots,n$ with any fixed
displacement $\bm{x}_{0}$, and thus the arguments on the lhs of Eq. (60)
should be given as those in the expression.

We apply the Fourier transform of the linear density contrast and the tensor
field,

$\displaystyle\tilde{\delta}_{\mathrm{L}}(\bm{k})$ $\displaystyle=\int
d^{3}\\!x\,e^{-i\bm{k}\cdot\bm{x}}\delta_{\mathrm{L}}(\bm{x}),$ (61)
$\displaystyle\tilde{F}_{Xlm}(\bm{k})$ $\displaystyle=\int
d^{3}\\!x\,e^{-i\bm{k}\cdot\bm{x}}F_{Xlm}(\bm{x}).$ (62)

and the functional derivative in Fourier space is given by

$\frac{\delta}{\delta\tilde{\delta}_{\mathrm{L}}(\bm{k})}=\frac{1}{(2\pi)^{3}}\int
d^{3}\\!x\,e^{i\bm{k}\cdot\bm{x}}\frac{\delta}{\delta\delta_{\mathrm{L}}(\bm{x})}.$
(63)

Using the above equations, the Fourier counterpart of the rhs of Eq. (60) is
calculated, and we have

$\left\langle\frac{\delta^{n}\tilde{F}_{Xlm}(\bm{k})}{\delta\tilde{\delta}_{\mathrm{L}}(\bm{k}_{1})\cdots\delta\tilde{\delta}_{\mathrm{L}}(\bm{k}_{n})}\right\rangle\\\
=(2\pi)^{3-3n}\delta_{\mathrm{D}}^{3}(\bm{k}_{1}+\cdots+\bm{k}_{n}-\bm{k})\tilde{\Gamma}^{(n)}_{Xlm}(\bm{k}_{1},\ldots\bm{k}_{n}),$
(64)

where

$\tilde{\Gamma}^{(n)}_{Xlm}(\bm{k}_{1},\ldots\bm{k}_{n})=\int
d^{3}\\!x_{1}\cdots
d^{3}\\!x_{n}\,e^{-i(\bm{k}_{1}\cdot\bm{x}_{1}+\cdots+\bm{k}_{n}\cdot\bm{x}_{n})}\\\
\times\Gamma^{(n)}_{Xlm}(\bm{x}_{1},\ldots\bm{x}_{n})$ (65)

is the Fourier transform of the propagator. The appearance of the delta
function on the rhs of Eq. (64) is the consequence of the statistical
homogeneity of the Universe. Integrating Eq. (64), we have

$\tilde{\Gamma}^{(n)}_{Xlm}(\bm{k}_{1},\ldots\bm{k}_{n})=(2\pi)^{3n}\int\frac{d^{3}k}{(2\pi)^{3}}\left\langle\frac{\delta^{n}\tilde{F}_{Xlm}(\bm{k})}{\delta\tilde{\delta}_{\mathrm{L}}(\bm{k}_{1})\cdots\delta\tilde{\delta}_{\mathrm{L}}(\bm{k}_{n})}\right\rangle.$
(66)

The Fourier transform of Eq. (59) is given by

$\tilde{F}_{Xlm}(\bm{k})=\int
d^{3}\\!q\,e^{-i\bm{k}\cdot\bm{q}-i\bm{k}\cdot\bm{\Psi}(\bm{q})}F^{\mathrm{L}}_{Xlm}(\bm{q}).$
(67)

In the iPT, the displacement field $\bm{\Psi}(\bm{q})$ is expanded as a
perturbative series of Lagrangian perturbation theory (LPT),

$\bm{\Psi}(\bm{q})=\sum_{n=1}^{\infty}\frac{i}{n!}\int\frac{d^{3}k_{1}}{(2\pi)^{3}}\cdots\frac{d^{3}k_{n}}{(2\pi)^{3}}e^{i(\bm{k}_{1}+\cdots+\bm{k}_{n})\cdot\bm{q}}\\\
\times\bm{L}_{n}(\bm{k}_{1},\ldots,\bm{k}_{n})\tilde{\delta}_{\mathrm{L}}(\bm{k}_{1})\cdots\tilde{\delta}_{\mathrm{L}}(\bm{k}_{n}).$
(68)

The perturbation kernels $\bm{L}_{n}$ are recursively obtained both for the
model of Einstein-de Sitter Universe [55] and for general models [52]. The
propagators of Eq. (66) are systematically evaluated by applying diagrammatic
rules in the iPT. The detailed derivation of the diagrammatic rules for the
number density field $\delta_{X}$ is given in Ref. [17]. Exactly the same
derivation applies to the tensor field, simply replacing
$1+\delta_{X}\rightarrow F_{Xlm}$ and $1+\delta^{\mathrm{L}}_{X}\rightarrow
F^{\mathrm{L}}_{Xlm}$ in that reference. In this method, the renormalized bias
functions play important roles. In the present case of the tensor field, the
renormalized bias functions $c^{(n)}_{Xlm}(\bm{k}_{1},\ldots,\bm{k}_{n})$ in
Fourier space are defined by

$\left\langle\frac{\delta^{n}\tilde{F}^{\mathrm{L}}_{Xlm}(\bm{k})}{\delta\tilde{\delta}_{\mathrm{L}}(\bm{k}_{1})\cdots\delta\tilde{\delta}_{\mathrm{L}}(\bm{k}_{n})}\right\rangle\\\
=(2\pi)^{3-3n}\delta_{\mathrm{D}}^{3}(\bm{k}_{1}+\cdots+\bm{k}_{n}-\bm{k})c^{(n)}_{Xlm}(\bm{k}_{1},\ldots\bm{k}_{n}),$
(69)

or equivalently,

$c^{(n)}_{Xlm}(\bm{k}_{1},\ldots\bm{k}_{n})=(2\pi)^{3n}\int\frac{d^{3}k}{(2\pi)^{3}}\left\langle\frac{\delta^{n}\tilde{F}^{\mathrm{L}}_{Xlm}(\bm{k})}{\delta\tilde{\delta}_{\mathrm{L}}(\bm{k}_{1})\cdots\delta\tilde{\delta}_{\mathrm{L}}(\bm{k}_{n})}\right\rangle.$
(70)

As can be seen by comparing the above definition with Eqs. (64) and (66), the
renormalized bias functions are the counterparts of the propagator for the
biasing in Lagrangian space.

The diagrammatic rules of iPT are summarized in the Appendix of Ref. [18]. The
same rules apply, simply replacing
$\Gamma_{X}^{(n)}\rightarrow\Gamma_{Xlm}^{(n)}$, $c_{X}^{(n)}\rightarrow
c_{Xlm}^{(n)}$, etc. For the terms without $c^{(n)}_{X}$ in that reference, an
additional factor $c^{(0)}_{Xlm}$ should be multiplied, which is unity in the
original formulation for scalar fields. The propagators are given in a form,

$\tilde{\Gamma}^{(n)}_{Xlm}(\bm{k}_{1},\ldots,\bm{k}_{n})=\Pi(\bm{k})\,\hat{\Gamma}^{(n)}_{Xlm}(\bm{k}_{1},\ldots,\bm{k}_{n}),$
(71)

where

$\Pi(\bm{k})=\left\langle
e^{-i\bm{k}\cdot\bm{\Psi}}\right\rangle=\exp\left[\sum_{n=2}^{\infty}\frac{(-i)^{n}}{n!}\left\langle(\bm{k}\cdot\bm{\Psi})^{n}\right\rangle_{\mathrm{c}}\right]$
(72)

is the vertex resummation factor in terms of the displacement field, and
$\langle\cdots\rangle_{\mathrm{c}}$ denotes the $n$-point connected part of
the random variables. The propagator $\hat{\Gamma}^{(n)}_{Xlm}$ without the
resummation factor $\Pi(\bm{k})$ is called reduced propagators in the
following.

Below, we give explicit forms of lower-order propagators, which can be used to
estimate power spectra, bispectra, etc. The vertex resummation factor in the
one-loop approximation of iPT is given by

$\Pi(\bm{k})=\exp\left\\{-\frac{1}{2}\int\frac{d^{3}p}{(2\pi)^{3}}\left[\bm{k}\cdot\bm{L}_{1}(\bm{p})\right]^{2}P_{\mathrm{L}}(p)\right\\},$
(73)

where $P_{\mathrm{L}}(k)$ is the linear power spectrum defined by

$\left\langle\tilde{\delta}_{\mathrm{L}}(\bm{k})\tilde{\delta}_{\mathrm{L}}(\bm{k})\right\rangle=(2\pi)^{3}\delta_{\mathrm{D}}^{3}(\bm{k}+\bm{k}^{\prime})P_{\mathrm{L}}(k),$
(74)

where the delta function on the rhs appears due to the statistical homogeneity
of the space. The reduced propagator of the first order, up to the one-loop
approximation, is given by

$\hat{\Gamma}_{Xlm}^{(1)}(\bm{k})=c_{Xlm}^{(1)}(\bm{k})+\left[\bm{k}\cdot\bm{L}_{1}(\bm{k})\right]c_{Xlm}^{(0)}\\\
+\int\frac{d^{3}p}{(2\pi)^{3}}P_{\mathrm{L}}(p)\biggl{\\{}\left[\bm{k}\cdot\bm{L}_{1}(-\bm{p})\right]c_{Xlm}^{(2)}(\bm{k},\bm{p})\\\
+\left[\bm{k}\cdot\bm{L}_{1}(-\bm{p})\right]\left[\bm{k}\cdot\bm{L}_{1}(\bm{k})\right]c_{Xlm}^{(1)}(\bm{p})\\\
+\left[\bm{k}\cdot\bm{L}_{2}(\bm{k},-\bm{p})\right]c_{Xlm}^{(1)}(\bm{p})\\\
+\frac{1}{2}\left[\bm{k}\cdot\bm{L}_{3}(\bm{k},\bm{p},-\bm{p})\right]c_{Xlm}^{(0)}\\\
+\left[\bm{k}\cdot\bm{L}_{1}(\bm{p})\right]\left[\bm{k}\cdot\bm{L}_{2}(\bm{k},-\bm{p})\right]c_{Xlm}^{(0)}\biggr{\\}}.$
(75)

The first line of the above equation corresponds to the lowest, or tree-level
approximation. The reduced propagator of second order in the tree-level
approximation is given by

$\hat{\Gamma}^{(2)}_{Xlm}(\bm{k}_{1},\bm{k}_{2})=c^{(2)}_{Xlm}(\bm{k}_{1},\bm{k}_{2})\\\
+\left[\bm{k}_{12}\cdot\bm{L}_{1}(\bm{k}_{1})\right]c^{(1)}_{Xlm}(\bm{k}_{2})+\left[\bm{k}_{12}\cdot\bm{L}_{1}(\bm{k}_{2})\right]c^{(1)}_{Xlm}(\bm{k}_{1})\\\
+\left\\{\left[\bm{k}_{12}\cdot\bm{L}_{1}(\bm{k}_{1})\right]\left[\bm{k}_{12}\cdot\bm{L}_{1}(\bm{k}_{2})\right]+\bm{k}_{12}\cdot\bm{L}_{2}(\bm{k}_{1},\bm{k}_{2})\right\\}c^{(0)}_{Xlm},$
(76)

where $\bm{k}_{12}=\bm{k}_{1}+\bm{k}_{2}$, and the reduced propagator of the
third order, again in the tree-level approximation, is given by

$\hat{\Gamma}^{(3)}_{Xlm}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})=c^{(3)}_{Xlm}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})\\\
+\left[\bm{k}_{123}\cdot\bm{L}_{1}(\bm{k}_{1})\right]c^{(2)}_{Xlm}(\bm{k}_{2},\bm{k}_{3})+\,\mathrm{cyc.}\\\
+\Bigl{\\{}\left[\bm{k}_{123}\cdot\bm{L}_{1}(\bm{k}_{1})\right]\left[\bm{k}_{123}\cdot\bm{L}_{1}(\bm{k}_{2})\right]+\left[\bm{k}_{123}\cdot\bm{L}_{2}(\bm{k}_{1},\bm{k}_{2})\right]\Bigr{\\}}\\\
\hskip 120.0pt\times c^{(1)}_{Xlm}(\bm{k}_{3})+\,\mathrm{cyc.}\\\
+\Bigl{\\{}\left[\bm{k}_{123}\cdot\bm{L}_{1}(\bm{k}_{1})\right]\left[\bm{k}_{123}\cdot\bm{L}_{1}(\bm{k}_{2})\right]\left[\bm{k}_{123}\cdot\bm{L}_{1}(\bm{k}_{3})\right]\\\
+\left[\bm{k}_{123}\cdot\bm{L}_{1}(\bm{k}_{1})\right]\left[\bm{k}_{123}\cdot\bm{L}_{2}(\bm{k}_{2},\bm{k}_{3})\right]+\,\mathrm{cyc.}\\\
+\bm{k}_{123}\cdot\bm{L}_{3}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})\Bigr{\\}}c^{(0)}_{Xlm},$
(77)

where $\bm{k}_{123}=\bm{k}_{1}+\bm{k}_{2}+\bm{k}_{3}$, and $+\,\mathrm{cyc.}$
represents two terms which are added with cyclic permutations of each previous
term.

In real space, the kernels of LPT in the standard theory of gravity (in the
Newtonian limit) are given by [56, 57]

$\displaystyle\bm{L}_{1}(\bm{k})=\frac{\bm{k}}{k^{2}},$ (78)
$\displaystyle\bm{L}_{2}(\bm{k}_{1},\bm{k}_{2})=\frac{3}{7}\frac{\bm{k}_{12}}{{k_{12}}^{2}}\left[1-\left(\frac{\bm{k}_{1}\cdot\bm{k}_{2}}{k_{1}k_{2}}\right)^{2}\right],$
(79)
$\displaystyle\bm{L}_{3}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})=\frac{1}{3}\left[\tilde{\bm{L}}_{3}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})+\mathrm{cyc.}\right],$
(80) $\displaystyle\tilde{\bm{L}}_{3}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$
$\displaystyle\quad=\frac{\bm{k}_{123}}{{k_{123}}^{2}}\left\\{\frac{5}{7}\left[1-\left(\frac{\bm{k}_{1}\cdot\bm{k}_{2}}{k_{1}k_{2}}\right)^{2}\right]\left[1-\left(\frac{\bm{k}_{12}\cdot\bm{k}_{3}}{{k}_{12}k_{3}}\right)^{2}\right]\right.$
$\displaystyle\qquad\quad\left.-\frac{1}{3}\left[1-3\left(\frac{\bm{k}_{1}\cdot\bm{k}_{2}}{k_{1}k_{2}}\right)^{2}+\,2\frac{(\bm{k}_{1}\cdot\bm{k}_{2})(\bm{k}_{2}\cdot\bm{k}_{3})(\bm{k}_{3}\cdot\bm{k}_{1})}{{k_{1}}^{2}{k_{2}}^{2}{k_{3}}^{2}}\right]\right\\}$
$\displaystyle\qquad+\frac{3}{7}\frac{\bm{k}_{123}}{{k_{123}}^{2}}\times\frac{(\bm{k}_{1}\times\bm{k}_{23})(\bm{k}_{1}\cdot\bm{k}_{23})}{{k_{1}}^{2}{k_{23}}^{2}}\left[1-\left(\frac{\bm{k}_{2}\cdot\bm{k}_{3}}{k_{2}k_{3}}\right)^{2}\right].$
(81)

In the above, weak dependencies on the time $t$ in the kernels are neglected
[8, 52]. In Ref. [52], complete expressions of the displacement kernels of LPT
up to the seventh order including the transverse parts are explicitly given,
together with a general way of recursively deriving the kernels including weak
dependencies on the time in general cosmology and subleading growing modes.
This method is also generalized to obtain kernels of LPT in modified theories
of gravity [58].

One of the benefits of the Lagrangian picture is that redshift-space
distortions are relatively easy to incorporate into the theory. The
displacement vector in redshift space $\bm{\Psi}^{\mathrm{s}}$ is obtained
from a linear mapping of that in real space by
$\bm{\Psi}^{\mathrm{s}}=\bm{\Psi}+H^{-1}(\hat{\bm{z}}\cdot\dot{\bm{\Psi}})\hat{\bm{z}}$,
where $\dot{\bm{\Psi}}=\partial\bm{\Psi}/\partial t$ is the time derivative of
the displacement field, $\hat{\bm{z}}$ is a unit vector directed to the line
of sight, $H=\dot{a}/a$ is the time-dependent Hubble parameter and $a(t)$ is
the scale factor. A displacement kernel in redshift space
$\bm{L}^{\mathrm{s}}_{n}$ is simply related to the kernel in real space in the
same order by a linear mapping [59]

$\bm{L}_{n}\rightarrow\bm{L}^{\mathrm{s}}_{n}=\bm{L}_{n}+nf\left(\hat{\bm{z}}\cdot\bm{L}_{n}\right)\hat{\bm{z}},$
(82)

where $f=d\ln D/d\ln a=\dot{D}/HD$ is the linear growth rate, $D(t)$ is the
linear growth factor. Throughout this paper, the distant-observer
approximation is assumed in redshift space and the unit vector $\hat{\bm{z}}$
denotes the line-of-sight direction. The expressions of Eqs. (73)–(77) apply
as well in redshift space when the displacement kernels $\bm{L}_{n}$ are
replaced by those in redshift space $\bm{L}^{\mathrm{s}}_{n}$.

###  III.2  Symmetries

####  III.2.1  Complex conjugate

We assume the Cartesian tensor field of Eq. (49) is a real tensor:

$F^{*}_{X\,i_{1}i_{2}\cdots i_{l}}(\bm{x})=F_{X\,i_{1}i_{2}\cdots
i_{l}}(\bm{x}),$ (83)

as it should be for physically observable quantities. For the spherical
decomposition of the traceless part, $F_{Xlm}(\bm{x})$, Eqs. (28), (42), (47)
and the above equation indicates

$F^{*}_{Xlm}(\bm{x})=(-1)^{l}g_{(l)}^{mm^{\prime}}F_{Xlm^{\prime}}(\bm{x})=(-1)^{l}F_{Xl}^{\phantom{X}m}(\bm{x}),$
(84)

and thus the Fourier transform of Eq. (62) satisfies

$\tilde{F}^{*}_{Xlm}(\bm{k})=(-1)^{l}g_{(l)}^{mm^{\prime}}\tilde{F}_{Xlm^{\prime}}(-\bm{k})=(-1)^{l}\tilde{F}_{Xl}^{\phantom{X}m}(-\bm{k}).$
(85)

That is, the complex conjugate of the irreducible tensor field in Fourier
space raises the azimuthal index, inverts the sign of wavevector in the
argument, and put a phase factor $(-1)^{l}$. The necessity of the last phase
factor comes from our convention of the $i^{l}$ factor in Eq. (41). The linear
density contrast $\delta_{\mathrm{L}}(\bm{x})$ is also a real field, and its
Fourier transform satisfies
${\tilde{\delta}_{\mathrm{L}}}^{\,*}(\bm{k})=\tilde{\delta}_{\mathrm{L}}(-\bm{k})$.
Therefore, the renormalized bias functions satisfy

$\displaystyle c^{(n)*}_{Xlm}(\bm{k}_{1},\ldots,\bm{k}_{n})$
$\displaystyle=(-1)^{l}g_{(l)}^{mm^{\prime}}c^{(n)}_{Xlm^{\prime}}(-\bm{k}_{1},\ldots,-\bm{k}_{n})$
$\displaystyle=(-1)^{l}c^{(n)m}_{Xl}(-\bm{k}_{1},\ldots,-\bm{k}_{n}).$ (86)

Similarly, the tensor propagator of Eq. (66) satisfies

$\displaystyle\tilde{\Gamma}^{(n)*}_{Xlm}(\bm{k}_{1},\ldots,\bm{k}_{n})$
$\displaystyle=(-1)^{l}g_{(l)}^{mm^{\prime}}\tilde{\Gamma}^{(n)}_{Xlm^{\prime}}(-\bm{k}_{1},\ldots,-\bm{k}_{n})$
$\displaystyle=(-1)^{l}\tilde{\Gamma}^{(n)m}_{Xl}(-\bm{k}_{1},\ldots,-\bm{k}_{n}).$
(87)

In redshift space, the propagators also depend on the direction of lines of
sight, $\hat{\bm{z}}$. The above property of the complex conjugate equally
applies to those in redshift space, and we have

$\displaystyle\tilde{\Gamma}^{(n)*}_{Xlm}(\bm{k}_{1},\ldots,\bm{k}_{n};\hat{\bm{z}})$
$\displaystyle=(-1)^{l}g_{(l)}^{mm^{\prime}}\tilde{\Gamma}^{(n)}_{Xlm^{\prime}}(-\bm{k}_{1},\ldots,-\bm{k}_{n};\hat{\bm{z}})$
$\displaystyle=(-1)^{l}\tilde{\Gamma}^{(n)m}_{Xl}(-\bm{k}_{1},\ldots,-\bm{k}_{n};\hat{\bm{z}}).$
(88)

####  III.2.2  Rotation

Under the passive rotation of Eq. (20), vector components $V_{i}$ generally
transform as

$V_{i}\rightarrow V_{i}^{\prime}=(R^{-1})_{ij}V_{j}=V_{j}R_{ji}.$ (89)

We denote the rotation of the position and wavevector components as, e.g.,
$\bm{x}\rightarrow\bm{x}^{\prime}=R^{-1}\bm{x}$ and
$\bm{k}\rightarrow\bm{k}^{\prime}=R^{-1}\bm{k}$, respectively. This notation
does not mean the physical rotations of vectors and does mean the rotation of
components, $x_{i}\rightarrow x_{i}^{\prime}=x_{j}R_{ji}$ and
$k_{i}\rightarrow k_{i}^{\prime}=k_{j}R_{ji}$, respectively, as a result of
the passive rotation of Cartesian basis, Eq. (20). The Cartesian tensor field
of Eq. (49) transforms as

$F_{X\,i_{1}i_{2}\cdots}(\bm{x})\xrightarrow{\mathbb{R}}F_{X\,i_{1}i_{2}\cdots}^{\prime}(\bm{x}^{\prime})=F_{X\,j_{1}j_{2}\cdots}(\bm{x})R_{j_{1}i_{1}}R_{j_{2}i_{2}}\cdots,$
(90)

and correspondingly the irreducible tensor of rank $l$ transforms as

$F_{Xlm}(\bm{x})\xrightarrow{\mathbb{R}}F_{Xlm}^{\prime}(\bm{x}^{\prime})=F_{Xlm^{\prime}}(\bm{x})D_{(l)m}^{m^{\prime}}(R).$
(91)

The corresponding transformation in Fourier space, Eq. (62), is given by

$\tilde{F}_{Xlm}(\bm{k})\xrightarrow{\mathbb{R}}\tilde{F}_{Xlm}^{\prime}(\bm{k}^{\prime})=\tilde{F}_{Xlm^{\prime}}(\bm{k})D_{(l)m}^{m^{\prime}}(R),$
(92)

because the volume element $d^{3}\\!x$ and the inner product
$\bm{k}\cdot\bm{x}$ are invariant under the rotation.

The same transformations apply to the tensor fields in Lagrangian space,
$F^{\mathrm{L}}_{Xlm}(\bm{q})$. Therefore, the renormalized bias functions,
Eq. (70), transform as

$c^{(n)}_{Xlm}(\bm{k}_{1},\ldots,\bm{k}_{n})\xrightarrow{\mathbb{R}}c^{(n)\prime}_{Xlm}(\bm{k}_{1}^{\prime},\ldots\bm{k}_{n}^{\prime})\\\
=c^{(n)}_{Xlm^{\prime}}(\bm{k}_{1},\ldots,\bm{k}_{n})D_{(l)m}^{m^{\prime}}(R).$
(93)

Similarly, the propagators of Eq. (66) transform as

$\tilde{\Gamma}^{(n)}_{Xlm}(\bm{k}_{1},\ldots\bm{k}_{n})\xrightarrow{\mathbb{R}}\tilde{\Gamma}^{(n)\prime}_{Xlm}(\bm{k}_{1}^{\prime},\ldots\bm{k}_{n}^{\prime})\\\
=\tilde{\Gamma}^{(n)}_{Xlm^{\prime}}(\bm{k}_{1},\ldots\bm{k}_{n})D_{(l)m}^{m^{\prime}}(R).$
(94)

In redshift space, the lines of sight also rotate,
$\hat{\bm{z}}\rightarrow\hat{\bm{z}}^{\prime}=R^{-1}\hat{\bm{z}}$, although
simultaneous rotations of wavevectors and lines of sight do not change the
scalar products in the propagators, as in Eqs. (73)–(82). Thus we have

$\tilde{\Gamma}^{(n)}_{Xlm}(\bm{k}_{1},\ldots\bm{k}_{n};\hat{\bm{z}})\xrightarrow{\mathbb{R}}\tilde{\Gamma}^{(n)\prime}_{Xlm}(\bm{k}_{1}^{\prime},\ldots\bm{k}_{n}^{\prime};\hat{\bm{z}}^{\prime})\\\
=\tilde{\Gamma}^{(n)}_{Xlm^{\prime}}(\bm{k}_{1},\ldots\bm{k}_{n};\hat{\bm{z}})D_{(l)m}^{m^{\prime}}(R).$
(95)

####  III.2.3  Parity

Next, we consider the property of parity symmetry. Keeping the left-handed
coordinates system, we consider an active parity transformation of the
physical system, instead of flipping the axes of coordinates system. With the
active parity transformation, the field values at a position $\bm{x}$ are
mapped into those at $-\bm{x}$, and the functional form of the tensor field is
transformed as

$F_{X\,i_{1}\cdots i_{l}}(\bm{x})\xrightarrow{\mathbb{P}}F_{X\,i_{1}\cdots
i_{l}}^{\prime}(\bm{x})=(-1)^{p_{X}+l}F_{X\,i_{1}\cdots i_{l}}(-\bm{x}),$ (96)

where $p_{X}=0$ for ordinary tensors and $p_{X}=1$ for pseudotensors. The
angular momentum is a typical example of a pseudotensor of rank 1. The
corresponding transformation for the irreducible tensor is given by

$F_{Xlm}(\bm{x})\xrightarrow{\mathbb{P}}F_{Xlm}^{\prime}(\bm{x})=(-1)^{p_{X}+l}F_{Xlm}(-\bm{x})$
(97)

in configuration space, and

$\tilde{F}_{Xlm}(\bm{k})\xrightarrow{\mathbb{P}}\tilde{F}_{Xlm}^{\prime}(\bm{k})=(-1)^{p_{X}+l}\tilde{F}_{Xlm}(-\bm{k})$
(98)

in Fourier space. The same applies to tensor fields in Lagrangian space.
Accordingly, the renormalized bias functions transform as

$c^{(n)}_{Xlm}(\bm{k}_{1},\ldots,\bm{k}_{n})\xrightarrow{\mathbb{P}}c^{(n)\prime}_{Xlm}(\bm{k}_{1},\ldots,\bm{k}_{n})\\\
=(-1)^{p_{X}+l}\,c^{(n)}_{Xlm}(-\bm{k}_{1},\ldots,-\bm{k}_{n}),$ (99)

for the linear density field transforms as
$\tilde{\delta}_{\mathrm{L}}(\bm{k})\rightarrow\tilde{\delta}_{\mathrm{L}}(-\bm{k})$.

Similarly, the propagators transform as

$\tilde{\Gamma}^{(n)}_{Xlm}(\bm{k}_{1},\ldots,\bm{k}_{n})\xrightarrow{\mathbb{P}}\tilde{\Gamma}^{(n)\prime}_{Xlm}(\bm{k}_{1},\ldots,\bm{k}_{n})\\\
=(-1)^{p_{X}+l}\,\tilde{\Gamma}^{(n)}_{Xlm}(-\bm{k}_{1},\ldots,-\bm{k}_{n}).$
(100)

In redshift space, we have

$\tilde{\Gamma}^{(n)}_{Xlm}(\bm{k}_{1},\ldots,\bm{k}_{n};\hat{\bm{z}})\xrightarrow{\mathbb{P}}\tilde{\Gamma}^{(n)\prime}_{Xlm}(\bm{k}_{1},\ldots,\bm{k}_{n};\hat{\bm{z}})\\\
=(-1)^{p_{X}+l}\,\tilde{\Gamma}^{(n)}_{Xlm}(-\bm{k}_{1},\ldots,-\bm{k}_{n};-\hat{\bm{z}}).$
(101)

####  III.2.4  Interchange of arguments

An obvious symmetry of the renormalized bias functions and propagators is that
they are invariant under permutation of the wavevectors in the argument. For
the renormalized bias function, we have

$c^{(n)}_{Xlm}(\bm{k}_{\sigma(1)},\ldots,\bm{k}_{\sigma(n)})=c^{(n)}_{Xlm}(\bm{k}_{1},\ldots\bm{k}_{n}),$
(102)

where $\sigma\in\mathcal{S}_{n}$ is a permutation of the symmetric group
$\mathcal{S}_{n}$ of order $n$. Any permutation can be realized by a series of
interchange operations of adjacent arguments:

$c^{(n)}_{Xlm}(\bm{k}_{1},\ldots,\bm{k}_{i},\bm{k}_{i+1},\ldots\bm{k}_{n})=c^{(n)}_{Xlm}(\bm{k}_{1},\ldots,\bm{k}_{i+1},\bm{k}_{i},\ldots\bm{k}_{n}),$
(103)

with arbitrary $i=1,\ldots,n-1$. Similarly, we have

$\tilde{\Gamma}^{(n)}_{Xlm}(\bm{k}_{1},\ldots,\bm{k}_{i},\bm{k}_{i+1},\ldots\bm{k}_{n})=\tilde{\Gamma}^{(n)}_{Xlm}(\bm{k}_{1},\ldots,\bm{k}_{i+1},\bm{k}_{i},\ldots\bm{k}_{n}),$
(104)

for the propagators in real space, and

$\tilde{\Gamma}^{(n)}_{Xlm}(\bm{k}_{1},\ldots,\bm{k}_{i},\bm{k}_{i+1},\ldots\bm{k}_{n};\hat{\bm{z}})\\\
=\tilde{\Gamma}^{(n)}_{Xlm}(\bm{k}_{1},\ldots,\bm{k}_{i+1},\bm{k}_{i},\ldots\bm{k}_{n};\hat{\bm{z}}),$
(105)

for the propagators in redshift space.

###  III.3  Renormalized bias functions of tensor fields

In the following of this paper, we mostly work in Fourier space. For
notational simplicity, we omit the tildes above variables in Fourier space,
and denote $\delta_{\mathrm{L}}(\bm{k})$, $F_{Xlm}(\bm{k})$,
$\Gamma^{(n)}_{Xlm}(\bm{k}_{1},\ldots)$ etc. instead of
$\tilde{\delta}_{\mathrm{L}}(\bm{k})$, $\tilde{F}_{Xlm}(\bm{k})$,
$\tilde{\Gamma}^{(n)}_{Xlm}(\bm{k}_{1},\ldots)$, respectively, as long as when
they are not confusing.

One of the essential parts of iPT is the introduction of the renormalized bias
function defined by Eq. (70). They characterize complicated fully nonlinear
processes in the formation of astronomical objects, such as various types of
galaxies and clusters, etc. If these complicated processes are described by
some kind of models and the function $F^{\mathrm{L}}_{Xlm}$ in Lagrangian
space is analytically given in terms of the linear density field
$\delta_{\mathrm{L}}$, then we can analytically calculate the renormalized
bias function $c^{(n)}_{Xlm}$ according to the model. In the case of the
number density of halos, the Press-Schechter model and its extensions are
applied to the concrete calculations along this line [19, 18, 20]. When this
kind of models are not known, the renormalized bias functions have infinite
degrees of freedom which is difficult to calculate from the fundamental level.
However, the rotational symmetry considered in the previous subsection places
some constraints on the functional form of the renormalized bias functions as
we show below.

We consider the renormalized bias functions given in the form of Eq. (70). In
the simplified notation without tildes,

$c^{(n)}_{Xlm}(\bm{k}_{1},\ldots\bm{k}_{n})=(2\pi)^{3n}\int\frac{d^{3}k}{(2\pi)^{3}}\left\langle\frac{\delta^{n}F^{\mathrm{L}}_{Xlm}(\bm{k})}{\delta\delta_{\mathrm{L}}(\bm{k}_{1})\cdots\delta\delta_{\mathrm{L}}(\bm{k}_{n})}\right\rangle.$
(106)

The angular dependencies of the wavevectors in the arguments can be decomposed
in a series expansion with spherical harmonics with Racah’s normalization, Eq.
(17), as

$c^{(n)}_{Xlm}(\bm{k}_{1},\ldots\bm{k}_{n})=\sum_{l_{1},\ldots,l_{n}}c^{(n)\,l_{1}\cdots
l_{n}}_{Xlm;m_{1}\cdots m_{n}}(k_{1},\ldots,k_{n})\\\ \times
C_{l_{1}m_{1}}^{*}(\hat{\bm{k}}_{1})\cdots
C_{l_{n}m_{n}}^{*}(\hat{\bm{k}}_{n}),$ (107)

where $k_{i}=|\bm{k}_{i}|$ are absolute values of the wavevectors, and
$\hat{\bm{k}}_{i}=(\theta_{ki},\phi_{ki})$ represents the spherical
coordinates for the direction of wavevectors $\bm{k}_{i}$. The repeated
indices $m$, $m_{1},\ldots,m_{n}$ in Eq. (107) are summed over according to
the Einstein summation convention as in the previous section. Due to the
orthonormality relation for the spherical harmonics, Eq. (380), the expansion
of Eq. (107) can be inverted and we have

$c^{(n)\,l_{1}\cdots l_{n}}_{Xlm;m_{1}\cdots
m_{n}}(k_{1},\ldots,k_{n})=(2l_{1}+1)\cdots(2l_{n}+1)\\\
\times\int\frac{d^{2}\hat{k}_{1}}{4\pi}\cdots\frac{d^{2}\hat{k}_{n}}{4\pi}c^{(n)}_{Xlm}(\bm{k}_{1},\ldots\bm{k}_{n})\\\
\times C_{l_{1}m_{1}}(\hat{\bm{k}}_{1})\cdots
C_{l_{n}m_{n}}(\hat{\bm{k}}_{n}),$ (108)

where $d^{2}\hat{k}_{i}=\sin\theta_{i}\,d\theta_{i}\,d\phi_{i}$ represents the
angular integration of the wavevector $\bm{k}_{i}$. Under the passive rotation
of Eq. (20), the function of Eq. (108) covariantly transforms as

$c^{(n)\,l_{1}\cdots l_{n}}_{Xlm;m_{1}\cdots
m_{n}}(k_{1},\ldots,k_{n})\xrightarrow{\mathbb{R}}c^{(n)\,l_{1}\cdots
l_{n}}_{Xlm^{\prime};m_{1}^{\prime}\cdots
m_{n}^{\prime}}(k_{1},\ldots,k_{n})\\\ \times
D_{(l)m}^{m^{\prime}}(R)D_{(l_{1})m_{1}}^{m_{1}^{\prime}}(R)\cdots
D_{(l_{n})m_{n}}^{m_{n}^{\prime}}(R),$ (109)

as obviously indicated by the lower positions of the azimuthal indices
$m,m_{1},\ldots m_{n}$.

For the statistically isotropic Universe, the functions of Eq. (108) are
rotationally invariant as these functions are given by ensemble averages of
the field variables as seen by Eqs. (106) and (108), and they should not
depend on the choice of coordinates system $\hat{\mathbf{e}}_{i}$ to describe
wavevectors $\bm{k}_{1},\ldots,\bm{k}_{n}$ of the renormalized bias functions
$c^{(n)}_{Xlm}(\bm{k}_{1},\ldots,\bm{k}_{n})$. Therefore, they are invariant
under the transformation of Eq. (109), and we have

$c^{(n)\,l_{1}\cdots l_{n}}_{Xlm;m_{1}\cdots
m_{n}}(k_{1},\ldots,k_{n})=c^{(n)\,l_{1}\cdots
l_{n}}_{Xlm^{\prime};m_{1}^{\prime}\cdots
m_{n}^{\prime}}(k_{1},\ldots,k_{n})\\\
\times\frac{1}{8\pi^{2}}\int[dR]D_{(l)m}^{m^{\prime}}(R)D_{(l_{1})m_{1}}^{m_{1}^{\prime}}(R)\cdots
D_{(l_{n})m_{n}}^{m_{n}^{\prime}}(R).$ (110)

where

$\int[dR]\cdots=\int_{0}^{2\pi}d\alpha\int_{0}^{\pi}d\beta\,\sin\beta\int_{0}^{2\pi}d\gamma\cdots$
(111)

represents the integrals over Euler angles $(\alpha\beta\gamma)$ of the
rotation $R$. A product of two Wigner’s rotation matrices of the same rotation
reduces to a single matrix through vector-coupling coefficients as [45]

$D_{(l_{1})m_{1}}^{m_{1}^{\prime}}(R)D_{(l_{2})m_{2}}^{m_{2}^{\prime}}(R)\\\
=\sum_{l,m,m^{\prime}}(2l+1)\begin{pmatrix}l_{1}&l_{2}&l\\\
m_{1}^{\prime}&m_{2}^{\prime}&m^{\prime}\end{pmatrix}\begin{pmatrix}l_{1}&l_{2}&l\\\
m_{1}&m_{2}&m\end{pmatrix}D_{(l)m}^{m^{\prime}\,*}(R),$ (112)

using the Wigner’s $3j$-symbol. The complex conjugate of the rotation matrix
is given by
$D_{(l)m_{1}}^{m_{2}*}(R)=g_{(l)}^{m_{1}m_{1}^{\prime}}g^{(l)}_{m_{2}m_{2}^{\prime}}D_{(l)m_{1}^{\prime}}^{m_{2}^{\prime}}(R)$
using our metric tensor for spherical basis.

In the following, we use a simplified notation for the $3j$-symbols,

$\left(l_{1}\,l_{2}\,l_{3}\right)_{m_{1}m_{2}m_{3}}\equiv\begin{pmatrix}l_{1}&l_{2}&l_{3}\\\
m_{1}&m_{2}&m_{3}\end{pmatrix},$ (113)

which is nonzero only when $m_{1}+m_{2}+m_{3}=0$. We consider that the
azimuthal indices can be raised or lowered by spherical metric, e.g.,

$\displaystyle\left(l_{1}\,l_{2}\,l_{3}\right)_{m_{1}}^{\phantom{m_{1}}m_{2}m_{3}}$
$\displaystyle=g_{(l_{2})}^{m_{2}m_{2}^{\prime}}g_{(l_{3})}^{m_{3}m_{3}^{\prime}}\left(l_{1}\,l_{2}\,l_{3}\right)_{m_{1}m_{2}^{\prime}m_{3}^{\prime}}$
$\displaystyle=(-1)^{m_{2}+m_{3}}\begin{pmatrix}l_{1}&l_{2}&l_{3}\\\
m_{1}&-m_{2}&-m_{3}\end{pmatrix},$ (114)

and so forth. The properties and formulas for Wigner’s $3j$-symbol with our
notation, which are repeatedly used in this paper below and in subsequent
papers in the series, are summarized in Appendix C.

With our notation of the $3j$-symbol, Eq. (112) is equivalent to an
expression,

$D_{(l_{1})m_{1}}^{m_{1}^{\prime}}(R)D_{(l_{2})m_{2}}^{m_{2}^{\prime}}(R)=(-1)^{l_{1}+l_{2}}\sum_{l=0}^{\infty}(-1)^{l}(2l+1)\\\
\times\left(l_{1}\,l_{2}\,l\right)_{m_{1}m_{2}}^{\phantom{m_{1}m_{2}}m}\left(l_{1}\,l_{2}\,l\right)^{m_{1}^{\prime}m_{2}^{\prime}}_{\phantom{m_{1}^{\prime}m_{2}^{\prime}}m^{\prime}}D_{(l)m}^{m^{\prime}}(R),$
(115)

in which the rotational covariance is explicit. We use a property of Eq. (404)
for the $3j$-symbol to derive the above equation.

Consecutive applications of Eq. (115) to Eq. (110), the dependence on the
rotation $R$ in the integrand on the rhs can be represented by a linear
combination of products of less than or equal to three rotation matrices and
the coefficients of the expression are given by products of $3j$-symbols. The
averages of products of rotation matrices of three or fewer are given by

$\displaystyle\frac{1}{8\pi^{2}}\int[dR]D_{(l)m}^{m^{\prime}}(R)=\delta_{l0}\delta_{m0}\delta_{m^{\prime}0},$
(116)
$\displaystyle\frac{1}{8\pi^{2}}\int[dR]D_{(l_{1})m_{1}}^{m_{1}^{\prime}}(R)D_{(l_{2})m_{2}}^{m_{2}^{\prime}}(R)=\frac{\delta_{l_{1}l_{2}}}{2l_{1}+1}g^{(l_{1})}_{m_{1}m_{2}}g_{(l_{1})}^{m_{1}^{\prime}m_{2}^{\prime}},$
(117)
$\displaystyle\frac{1}{8\pi^{2}}\int[dR]D_{(l_{1})m_{1}}^{m_{1}^{\prime}}(R)D_{(l_{2})m_{2}}^{m_{2}^{\prime}}(R)D_{(l_{3})m_{3}}^{m_{3}^{\prime}}(R)$
$\displaystyle\hskip
60.0pt=(-1)^{l_{1}+l_{2}+l_{3}}\left(l_{1}\,l_{2}\,l_{3}\right)_{m_{1}m_{2}m_{3}}\left(l_{1}\,l_{2}\,l_{3}\right)^{m_{1}^{\prime}m_{2}^{\prime}m_{3}^{\prime}}.$
(118)

Eqs. (116) and (117) are also derived from Eq. (118) noting
$D_{(0)0}^{0}(R)=1$ and Eqs. (410). Also, Eqs. (117) and (118) are derived
from further applications of Eq. (115) to reduce the number of products and
finally using Eq. (116).

As a result of the above procedure, averaging over the rotation $R$, the
dependence of the indices of $m$’s in Eq. (110) are all represented by
spherical metrics and $3j$-symbols. First, we specifically derive the results
for lower orders of $n$. For $n=0$, using Eqs. (110) and (116), we derive

$c^{(0)}_{Xlm}=c^{(0)}_{X00}\delta_{l0}\delta_{m0}.$ (119)

This result is trivial because the corresponding function of Eq. (106) for
$n=0$ does not depend on the direction of the wavevector. For notational
simplicity, we simply denote $c^{(0)}_{X}\equiv c^{(0)}_{X00}$, and we have

$c^{(0)}_{Xlm}=\delta_{l0}\delta_{m0}c^{(0)}_{X}.$ (120)

For a scalar function, $F_{X}$, the coefficient $c^{(0)}_{X}$ just corresponds
to the mean value, $c^{(0)}_{X}=c^{(0)}_{X00}=\langle F_{X00}\rangle=\langle
F_{X00}\rangle\mathsf{Y}^{(0)}=\langle F_{X}\rangle$, as seen from Eqs. (7)
and (31). Due to the assumption that the field $F_{X}$ is a real field, the
coefficient $c^{(0)}_{X}$ has a real value:

$c^{(0)*}_{X}=c^{(0)}_{X}.$ (121)

The parity transformation is given by

$c^{(0)}_{X}\xrightarrow{\mathbb{P}}(-1)^{p_{X}}c^{(0)}_{X}.$ (122)

In the Universe with parity symmetry, we should have $p_{X}=0$, because the
mean value of a scalar field $\langle F_{X}\rangle$ vanish if the field is a
pseudo-scalar in the Universe with parity symmetry.

For $n=1$, using Eqs. (110) and (117), we derive

$c^{(1)\,l_{1}}_{Xlm;m_{1}}(k)=\frac{\delta_{l_{1}l}}{2l+1}g^{(l)}_{mm_{1}}g_{(l)}^{m^{\prime}m_{1}^{\prime}}c^{(1)\,l}_{Xlm^{\prime};m_{1}^{\prime}}(k),$
(123)

and therefore both indices $m,m_{1}$ appear only in the metric
$g^{(l)}_{mm_{1}}$. Defining a rotationally invariant function,

$c^{(1)}_{Xl}(k)\equiv\frac{(-1)^{l}}{\sqrt{2l+1}}g_{(l)}^{mm_{1}}c^{(1)\,l}_{Xlm;m_{1}}(k),$
(124)

one finds that Eq. (123) reduces to a simple form,

$c^{(1)\,l_{1}}_{Xlm;m_{1}}(k)=\frac{(-1)^{l}}{\sqrt{2l+1}}\delta_{ll_{1}}g^{(l)}_{mm_{1}}c^{(1)}_{Xl}(k).$
(125)

Substituting this expression into Eq. (107) with $n=1$, we derive

$c^{(1)}_{Xlm}(\bm{k})=\frac{(-1)^{l}}{\sqrt{2l+1}}\,c^{(1)}_{Xl}(k)C_{lm}(\hat{\bm{k}}).$
(126)

The invariant coefficient $c^{(1)}_{Xl}(k)$ can also be read off from an
expression of the renormalized bias function $c^{(1)}_{Xlm}(\bm{k})$ in which
the angular dependence on $\bm{k}$ is expanded by the spherical harmonics
$C_{lm}(\hat{\bm{k}})$. One can also explicitly invert the above equation by
using the orthonormality relation of spherical harmonics.

Combining the properties of complex conjugate for the spherical harmonics, Eq.
(377), and for the renormalized bias function, Eq. (86) with $n=1$, we have

$c^{(1)*}_{Xl}(k)=c^{(1)}_{Xl}(k),$ (127)

and thus the reduced function is a real function. Combining the properties of
parity for the spherical harmonics, Eq. (376), and for the renormalized bias
function, Eq. (99) with $n=1$, we have

$c^{(1)}_{Xl}(k)\xrightarrow{\mathbb{P}}(-1)^{p_{X}}c^{(1)}_{Xl}(k).$ (128)

Thus the Universe with parity symmetry suggests $p_{X}=0$ to have a non-zero
value for the function. This means that pseudotensors do not have the first-
order renormalized bias function in the presence of parity symmetry. This
property is partly because we consider the case that only scalar perturbations
$\delta_{\mathrm{L}}$ are responsible for the properties of tensor fields. For
a typical example that is a manifestation of the last conclusion, angular
momenta of galaxies are not generated by the first-order effect in linear
perturbation theory [60].

For $n=2$, using Eqs. (110) and (118), we derive

$c^{(2)\,l_{1}l_{2}}_{Xlm;m_{1}m_{2}}(k_{1},k_{2})=\left(l\,l_{1}\,l_{2}\right)_{mm_{1}m_{2}}c^{(2)\,l}_{Xl_{1}l_{2}}(k_{1},k_{2}),$
(129)

where we define a rotationally invariant function,

$c^{(2)\,l}_{Xl_{1}l_{2}}(k_{1},k_{2})\equiv(-1)^{l+l_{1}+l_{2}}\left(l\,l_{1}\,l_{2}\right)^{mm_{1}m_{2}}c^{(2)\,l_{1}l_{2}}_{Xlm;m_{1}m_{2}}(k_{1},k_{2}).$
(130)

The above Eq. (129) essentially corresponds to the Wigner-Eckart theorem of
representation theory and quantum mechanics [47], and the dependencies on
azimuthal indices $m,m_{1},m_{2}$ appear only through $3j$-symbols.
Substituting Eq. (129) into Eq. (107) with $n=2$, we derive

$c^{(2)}_{Xlm}(\bm{k}_{1},\bm{k}_{2})=\sum_{l_{1},l_{2}}c^{(2)\,l}_{Xl_{1}l_{2}}(k_{1},k_{2})\left(l\,l_{1}\,l_{2}\right)_{m}^{\phantom{m}m_{1}m_{2}}\\\
\times C_{l_{1}m_{1}}(\hat{\bm{k}}_{1})C_{l_{2}m_{2}}(\hat{\bm{k}}_{2}).$
(131)

The last expression contains the bipolar spherical harmonics [46] which is
defined by Eq. (383) in Appendix B with our notation. Rather than using the
conventional notation for the bipolar spherical harmonics,
$\\{Y_{l_{1}}(\hat{\bm{k}}_{1})\otimes Y_{l_{2}}(\hat{\bm{k}}_{2})\\}_{lm}$,
we introduce a notation with convenient normalization of the bipolar spherical
harmonics,

$X^{l_{1}l_{2}}_{lm}(\hat{\bm{k}}_{1},\hat{\bm{k}}_{2})\equiv\frac{(-1)^{l}\,4\pi}{\sqrt{\\{l\\}\\{l_{1}\\}\\{l_{2}\\}}}\left\\{Y_{l_{1}}(\hat{\bm{k}}_{1})\otimes
Y_{l_{2}}(\hat{\bm{k}}_{2})\right\\}_{lm},$ (132)

where $\\{l\\}=2l+1$, $\\{l_{1}\\}=2l_{1}+1$ $\\{l_{2}\\}=2l_{2}+1$. These
factors frequently appear in the equations below, and we use the same notation
for arbitrary non-negative integer $L$,

$\\{L\\}\equiv 2L+1.$ (133)

throughout the paper henceforth. The bipolar spherical harmonics satisfy the
orthonormality relation, Eq. (388), and product formula, Eq. (390), as shown
in Appendix B. With the above notation, Eq. (131) is concisely given by

$c^{(2)}_{Xlm}(\bm{k}_{1},\bm{k}_{2})=\sum_{l_{1},l_{2}}c^{(2)\,l}_{Xl_{1}l_{2}}(k_{1},k_{2})X^{l_{1}l_{2}}_{lm}(\hat{\bm{k}}_{1},\hat{\bm{k}}_{2}).$
(134)

Combining the properties of complex conjugate for the bipolar spherical
harmonics, Eq. (387), and for the renormalized bias function, Eq. (86) with
$n=2$, we have

$c^{(2)\,l*}_{Xl_{1}l_{2}}(k_{1},k_{2})=c^{(2)\,l}_{Xl_{1}l_{2}}(k_{1},k_{2})$
(135)

and thus the reduced function is a real function. Combining the properties of
parity for the bipolar spherical harmonics, Eq. (386), and for the
renormalized bias function, Eq. (99) with $n=2$, we have

$c^{(2)\,l}_{Xl_{1}l_{2}}(k_{1},k_{2})\xrightarrow{\mathbb{P}}(-1)^{p_{X}+l+l_{1}+l_{2}}c^{(2)\,l}_{Xl_{1}l_{2}}(k_{1},k_{2}).$
(136)

Thus the Universe with parity symmetry suggests
$p_{X}+l+l_{1}+l_{2}=\mathrm{even}$ to have a non-zero value for the function.
The interchange symmetry of Eq. (103) in this case is given by
$c^{(2)}_{Xlm}(\bm{k}_{1},\bm{k}_{2})=c^{(2)}_{Xlm}(\bm{k}_{2},\bm{k}_{1})$.
Due to the symmetry of $3j$-symbols of Eq. (406), the bipolar spherical
harmonics satisfy

$X^{l_{1}l_{2}}_{lm}(\hat{\bm{k}}_{1},\hat{\bm{k}}_{2})=(-1)^{l+l_{1}+l_{2}}X^{l_{2}l_{1}}_{lm}(\hat{\bm{k}}_{2},\hat{\bm{k}}_{1}),$
(137)

and therefore the invariant function satisfies

$c^{(2)\,l}_{Xl_{2}l_{1}}(k_{2},k_{1})=(-1)^{l+l_{1}+l_{2}}c^{(2)\,l}_{Xl_{1}l_{2}}(k_{1},k_{2}).$
(138)

For $n=3$, using Eqs. (110), (115) and (118), we derive

$c^{(3)\,l_{1}l_{2}l_{3}}_{Xlm;m_{1}m_{2}m_{3}}(k_{1},k_{2},k_{3})=\sum_{L}(-1)^{L}\sqrt{\\{L\\}}\left(l\,l_{1}\,L\right)_{mm_{1}M}\\\
\times\left(L\,l_{2}\,l_{3}\right)^{M}_{\phantom{M}m_{2}m_{3}}c^{(3)\,l;L}_{Xl_{1}l_{2}l_{3}}(k_{1},k_{2},k_{3}),$
(139)

where we define

$c^{(3)\,l;L}_{Xl_{1}l_{2}l_{3}}(k_{1},k_{2},k_{3})\equiv(-1)^{l+l_{1}+l_{2}+l_{3}+L}\sqrt{\\{L\\}}\left(l\,l_{1}\,L\right)^{mm_{1}}_{\phantom{mm_{1}}M}\\\
\times\left(L\,l_{2}\,l_{3}\right)^{Mm_{2}m_{3}}c^{(3)\,l_{1}l_{2}l_{3}}_{Xlm;m_{1}m_{2}m_{3}}(k_{1},k_{2},k_{3}).$
(140)

Substituting Eq. (139) into Eq. (107) with $n=3$, we derive

$c^{(3)}_{Xlm}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})=\sum_{l_{1},l_{2},l_{3},L}(-1)^{L}\sqrt{\\{L\\}}\\\
\times
c^{(3)\,l;L}_{Xl_{1}l_{2}l_{3}}(k_{1},k_{2},k_{3})\left(l\,l_{1}\,L\right)_{m}^{\phantom{m}m_{1}M}\left(L\,l_{2}\,l_{3}\right)_{M}^{\phantom{M}m_{2}m_{3}}\\\
\times
C_{l_{1}m_{1}}(\hat{\bm{k}}_{1})C_{l_{2}m_{2}}(\hat{\bm{k}}_{2})C_{l_{3}m_{3}}(\hat{\bm{k}}_{3}).$
(141)

The last expression contains the tripolar spherical harmonics [46]. Similar to
Eq. (132) of bipolar spherical harmonics, we introduce a notation for the
tripolar spherical harmonics,

$X^{l_{1}l_{2}l_{3}}_{L;lm}(\hat{\bm{k}}_{1},\hat{\bm{k}}_{2},\hat{\bm{k}}_{3})\equiv\frac{(-1)^{l}\,(4\pi)^{3/2}}{\sqrt{\\{l\\}\\{l_{1}\\}\\{l_{2}\\}\\{l_{3}\\}}}\\\
\times\left\\{Y_{l_{1}}(\hat{\bm{k}}_{1})\otimes\left\\{Y_{l_{2}}(\hat{\bm{k}}_{2})\otimes
Y_{l_{3}}(\hat{\bm{k}}_{3})\right\\}_{L}\right\\}_{lm}.$ (142)

Explicit form of the above harmonics is defined by Eq. (391), and their
orthonormality relation, Eq. (395) and product formula, Eq. (396) are given in
Appendix B with our notations for convenience. Thereby Eq. (141) is concisely
represented as

$c^{(3)}_{Xlm}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})=\sum_{l_{1},l_{2},l_{3},L}c^{(3)\,l;L}_{Xl_{1}l_{2}l_{3}}(k_{1},k_{2},k_{3})X^{l_{1}l_{2}l_{3}}_{L;lm}(\hat{\bm{k}}_{1},\hat{\bm{k}}_{2},\hat{\bm{k}}_{3}).$
(143)

Combining the properties of complex conjugate for the tripolar spherical
harmonics, Eq. (394), and for the renormalized bias function, Eq. (86) with
$n=3$, we have

$c^{(3)\,l;L\,*}_{Xl_{1}l_{2}l_{3}}(k_{1},k_{2},k_{3})=c^{(3)\,l;L}_{Xl_{1}l_{2}l_{3}}(k_{1},k_{2},k_{3}),$
(144)

and thus the reduced function is a real function. Combining the properties of
parity for the tripolar spherical harmonics, Eq. (393), and for the
renormalized bias function, Eq. (99) with $n=3$, we have

$c^{(3)\,l;L}_{Xl_{1}l_{2}l_{3}}(k_{1},k_{2},k_{3})\xrightarrow{\mathbb{P}}(-1)^{p_{X}+l+l_{1}+l_{2}+l_{3}}c^{(3)\,l;L}_{Xl_{1}l_{2}l_{3}}(k_{1},k_{2},k_{3}).$
(145)

Thus the Universe with parity symmetry suggests
$p_{X}+l+l_{1}+l_{2}+l_{3}=\mathrm{even}$ to have a non-zero value for the
function. The interchange symmetries of Eq. (103) in terms of invariant
functions in this case are derived from properties of interchanging arguments
in tripolar spherical harmonics. For an interchange of the last two arguments,
$2\leftrightarrow 3$, we have

$X^{l_{1}l_{2}l_{3}}_{L;lm}(\hat{\bm{k}}_{1},\hat{\bm{k}}_{2},\hat{\bm{k}}_{3})=(-1)^{l_{2}+l_{3}+L}X^{l_{1}l_{3}l_{2}}_{L;lm}(\hat{\bm{k}}_{1},\hat{\bm{k}}_{3},\hat{\bm{k}}_{2}),$
(146)

just similarly in the case of Eq. (137). For an interchange of the first two
arguments, $1\leftrightarrow 2$, however, a recoupling of the $3j$-symbols
appears, and we have

$X^{l_{1}l_{2}l_{3}}_{L;lm}(\hat{\bm{k}}_{1},\hat{\bm{k}}_{2},\hat{\bm{k}}_{3})=(-1)^{l_{1}+l_{2}}\sum_{L^{\prime}}\\{L^{\prime}\\}\begin{Bmatrix}l_{1}&l&L\\\
l_{2}&l_{3}&L^{\prime}\end{Bmatrix}\\\ \times
X^{l_{2}l_{1}l_{3}}_{L^{\prime};lm}(\hat{\bm{k}}_{2},\hat{\bm{k}}_{1},\hat{\bm{k}}_{3}),$
(147)

where the factor in front of the tripolar spherical harmonics on the rhs is a
$6j$-symbol defined by Eq. (417) or equivalently by Eq. (418). Useful
properties of $6j$-symbols are summarized in Eqs. (419), (420) and (426) of
Appendix C. The above relation of Eq. (147) is derived by applying a sum rule
of the $3j$-symbols, Eq. (426), to the definition of the tripolar spherical
harmonics, Eq. (391). Correspondingly, the interchange symmetries of invariant
coefficients of Eq. (143) are given by

$c^{(3)l;L}_{Xl_{l}l_{3}l_{2}}(k_{1},k_{3},k_{2})=(-1)^{l_{2}+l_{3}+L}c^{(3)l;L}_{Xl_{l}l_{2}l_{3}}(k_{1},k_{2},k_{3}),$
(148)

and

$c^{(3)l;L}_{Xl_{2}l_{1}l_{3}}(k_{2},k_{1},k_{3})=(-1)^{l_{1}+l_{2}}\sum_{L^{\prime}}\\{L^{\prime}\\}\begin{Bmatrix}l_{1}&l&L\\\
l_{2}&l_{3}&L^{\prime}\end{Bmatrix}\\\ \times
c^{(3)l;L^{\prime}}_{Xl_{l}l_{2}l_{3}}(k_{1},k_{2},k_{3}).$ (149)

Combining Eqs. (148) and (149), all the other symmetries concerning the
permutation of $(1,2,3)$ in the subscripts of the arguments are
straightforwardly obtained. A similar kind of considerations of the
interchange symmetry in a formulation of angular trispectrum of the cosmic
microwave background are found in Ref. [61] to enforce permutation symmetry in
the angular trispectrum.

In the same way, one derives the results for general orders $n>3$. The
expansion coefficient of order $n$ is given by

$c^{(n)\,l_{1}\cdots l_{n}}_{Xlm;m_{1}\cdots m_{n}}(k_{1},\ldots,k_{n})\\\
=\sum_{L_{2},\ldots,L_{n-1}}(-1)^{L_{2}+\cdots
L_{n-1}}\sqrt{\\{L_{2}\\}\cdots\\{L_{n-1}\\}}\left(l\,l_{1}\,L_{2}\right)_{mm_{1}M_{2}}\\\
\times\left(L_{2}\,l_{2}\,L_{3}\right)^{M_{2}}_{\phantom{M_{2}}m_{2}M_{3}}\cdots\left(L_{n-2}\,l_{n-2}\,L_{n-1}\right)^{M_{n-2}}_{\phantom{M_{n-2}}m_{2}M_{n-1}}\\\
\times\left(L_{n-1}\,l_{n-1}\,l_{n}\right)^{M_{n-1}}_{\phantom{M_{n-1}}m_{n-1}m_{n}}c^{(n)\,l;L_{2}\cdots
L_{n-1}}_{Xl_{1}\cdots l_{n}}(k_{1},\ldots,k_{n}),$ (150)

where we define

$c^{(n)\,l;L_{2}\cdots L_{n-1}}_{Xl_{1}\cdots
l_{n}}(k_{1},\ldots,k_{n})\equiv(-1)^{l+l_{1}+\cdots+l_{n}}\\\
\times(-1)^{L_{2}+\cdots+L_{n-1}}\sqrt{\\{L_{2}\\}\cdots\\{L_{n-1}\\}}\left(l\,l_{1}\,L_{2}\right)^{mm_{1}}_{\phantom{mm_{1}}M_{2}}\\\
\times\left(L_{2}\,l_{2}\,L_{3}\right)^{M_{2}m_{2}}_{\phantom{M_{2}m_{2}}M_{3}}\cdots\left(L_{n-2}\,l_{n-2}\,L_{n-1}\right)^{M_{n-2}m_{n-2}}_{\phantom{M_{n-2}m_{n-2}}M_{n-1}}\\\
\times\left(L_{n-1}\,l_{n-1}\,l_{n}\right)^{M_{n-1}m_{n-1}m_{n}}c^{(n)\,l_{1}\cdots
l_{n}}_{Xlm;m_{1}\cdots m_{n}}(k_{1},\ldots,,k_{n}).$ (151)

Thus the azimuthal indices of the decomposed renormalized bias functions, Eq.
(108) are all represented by combinations of $3j$-symbols, and the physical
contents of the renormalized bias functions are given by the invariant
functions, $c^{(n)\,l;L_{2}\cdots L_{n-1}}_{Xl_{1}\cdots
l_{n}}(k_{1},\ldots,k_{n})$. The fully rotational symmetry is the reason why
the azimuthal indices of decomposed renormalized bias function,
$c^{(n)\,l;l_{1}\cdots l_{n}}_{Xm;m_{1}\cdots m_{n}}(\ldots)$ are represented
by the combinations of $3j$-symbols. Due to the symmetry of $3j$-symbols, we
have $m+m_{1}+\cdots+m_{n}=0$, which is a manifestation of the axial symmetry
of the rotation around the third axis $\hat{\mathbf{e}}_{3}$.

The renormalized bias function of Eq. (107) is given by

$c^{(n)}_{Xlm}(\bm{k}_{1},\ldots,\bm{k}_{n})=\sum_{\begin{subarray}{c}l_{1},\ldots,l_{n}\\\
L_{2},\ldots,L_{n-1}\end{subarray}}(-1)^{L_{2}+\cdots+L_{n-1}}\sqrt{\\{L_{2}\\}\cdots\\{L_{n-1}\\}}\\\
\times c^{(n)\,l;L_{2}\cdots L_{n-1}}_{Xl_{1}\cdots
l_{n}}(k_{1},\ldots,k_{n})\left(l\,l_{1}\,L_{2}\right)_{m}^{\phantom{m}m_{1}M_{2}}\\\
\times\left(L_{2}\,l_{2}\,L_{3}\right)_{M_{2}}^{\phantom{M_{2}}m_{2}M_{3}}\cdots\left(L_{n-2}\,l_{n-2}\,L_{n-1}\right)_{M_{n-2}}^{\phantom{M_{n-2}}m_{n-2}M_{n-1}}\\\
\times\left(L_{n-1}\,l_{n-1}\,l_{n}\right)_{M_{n-1}}^{\phantom{M_{n-1}}m_{n-1}m_{n}}C_{l_{1}m_{1}}(\hat{\bm{k}}_{1})\cdots
C_{l_{n}m_{n}}(\hat{\bm{k}}_{n}).$ (152)

The last expression contains the polypolar spherical harmonics defined by Eq.
(397), and their orthonormality relation, Eq. (401) and product formula, Eq.
(402) are given in Appendix B with our notations for convenience. Similar to
Eqs. (132) and (142) of bipolar and tripolar spherical harmonics, we introduce
a notation for the polypolar spherical harmonics,

$X^{l_{1}\cdots l_{n}}_{L_{2}\cdots
L_{n-1};lm}(\hat{\bm{k}}_{1},\ldots,\hat{\bm{k}}_{n})\equiv\frac{(-1)^{l}\,(4\pi)^{n/2}}{\sqrt{\\{l\\}\\{l_{1}\\}\cdots\\{l_{n}\\}}}\\\
\times\left\\{Y_{l_{1}}(\hat{\bm{k}}_{1})\otimes\left\\{Y_{l_{2}}(\hat{\bm{k}}_{2})\otimes\left\\{\cdots\otimes
Y_{l_{n}}(\hat{\bm{k}}_{n})\right\\}_{L_{n-1}}\cdots\right\\}_{L_{2}}\right\\}_{lm},$
(153)

and Eq. (152) is concisely represented as

$c^{(n)}_{Xlm}(\bm{k}_{1},\cdots,\bm{k}_{n})=\sum_{\begin{subarray}{c}l_{1},\ldots,l_{n}\\\
L_{2},\ldots,L_{n-1}\end{subarray}}c^{(n)\,l;L_{2}\cdots
L_{n-1}}_{Xl_{1}\cdots l_{n}}(k_{1},\ldots,k_{n})\\\ \times X^{l_{1}\cdots
l_{n}}_{L_{2}\cdots L_{n-1};lm}(\hat{\bm{k}}_{1},\ldots,\hat{\bm{k}}_{n}).$
(154)

Combining the properties of complex conjugate for the polypolar spherical
harmonics, Eq. (400), and for the renormalized bias function, Eq. (86), we
have

$c^{(n)\,l;L_{2}\cdots L_{n-1}*}_{Xl_{1}\cdots
l_{n}}(k_{1},\ldots,k_{n})=c^{(n)\,l;L_{2}\cdots L_{n-1}}_{Xl_{1}\cdots
l_{n}}(k_{1},\ldots,k_{n}),$ (155)

and thus the reduced function is a real function. Combining the properties of
parity for the polypolar spherical harmonics, Eq. (399), and for the
renormalized bias function, Eq. (99), we have

$c^{(n)\,l;L_{2}\cdots L_{n-1}}_{Xl_{1}\cdots l_{n}}(k_{1},\ldots,k_{3})\\\
\xrightarrow{\mathbb{P}}(-1)^{p_{X}+l+l_{1}+\cdots+l_{n}}c^{(n)\,l;L_{2}\cdots
L_{n-1}}_{Xl_{1}\cdots l_{n}}(k_{1},\ldots,k_{3}).$ (156)

Thus the Universe with parity symmetry suggests
$p_{X}+l+l_{1}+\cdots+l_{n}=\mathrm{even}$ to have a non-zero value for the
function.

The interchange symmetries of Eq. (103) in terms of invariant function in this
case are also derived following a similar way in the case of $n=3$. For the
interchange of the last two arguments, we have

$c^{(n)l;L_{2}\cdots L_{n-2}L_{n-1}}_{Xl_{l}\cdots
l_{n-2}l_{n}l_{n-1}}(k_{1},\ldots,k_{n-2},k_{n},k_{n-1})\\\
=(-1)^{l_{n-1}+l_{n}+L_{n-1}}c^{(n)l;L_{2}\cdots L_{n-1}}_{Xl_{l}\cdots
l_{n}}(k_{1},\ldots,k_{n}),$ (157)

and for the interchange of the other adjacent arguments, $i\leftrightarrow
i+1$, we have

$c^{(n)\,l;L_{2}\cdots L_{i-1}L_{i+1}L_{i}L_{i+2}\cdots L_{n-1}}_{Xl_{1}\cdots
l_{i-1}l_{i+1}l_{i}l_{i+2}\cdots
l_{n}}(k_{1},\ldots,k_{i-1},k_{i+1},k_{i},k_{i+2},\ldots,k_{n})\\\
=(-1)^{l_{i}+l_{i+1}}\sum_{L^{\prime}}\\{L^{\prime}\\}\begin{Bmatrix}l_{i}&L_{i}&L_{i+1}\\\
l_{i+1}&L_{i+2}&L^{\prime}\end{Bmatrix}\\\ \times c^{(n)\,l;L_{2}\cdots
L_{i}L^{\prime}L_{i+2}\cdots L_{n-1}}_{Xl_{1}\cdots
l_{n}}(k_{1},\ldots,k_{n}),$ (158)

where $i=1,\ldots,n-2$. Combining Eqs. (157) and (158), all the other
symmetries concerning the permutation of $(1,\ldots,n)$ in the subscripts of
the arguments are straightforwardly obtained.

###  III.4  Propagators of tensor fields

The renormalized bias functions describe the bias mechanisms in the physical
space, and thus there is not any corresponding concept of redshift space.
However, the apparent clustering in redshift space introduces anisotropy in
the propagators in Eulerian space. In the following, we separately consider
the propagators in real space and redshift space.

####  III.4.1  Real space

In real space, where the fully rotational symmetry is satisfied even in
Eulerian space, the reduced propagators $\hat{\Gamma}^{(n)}_{Xlm}$ are
similarly decomposed into rotationally invariant coefficients as in the case
of renormalized bias functions explained in the previous subsection. The
derivation of the invariant propagators is mostly parallel to the latter case.
In real space, formally replacing the renormalized bias functions
$c^{(n)}_{Xlm}$ with normalized propagators $\hat{\Gamma}^{(n)}_{Xlm}$ in Sec.
III.3 is enough to derive necessary equations of invariant propagators. We
summarize some important equations explicitly below for later use. The
derivations are exactly the same as those of renormalized bias functions.

The propagator of first order, $n=1$, is given by

$\hat{\Gamma}^{(1)}_{Xlm}(\bm{k})=\frac{(-1)^{l}}{\sqrt{2l+1}}\,\hat{\Gamma}^{(1)}_{Xl}(k)C_{lm}(\hat{\bm{k}}),$
(159)

with the invariant propagator $\hat{\Gamma}^{(1)}_{Xl}(k)$ which is a real
function. If the Universe is statistically invariant under the parity, the
first-order propagator is nonzero only when $p_{X}=0$ and vanishes for
pseudotensors. The parity transformation of the invariant function is given by

$\hat{\Gamma}^{(1)}_{Xl}(k)\xrightarrow{\mathbb{P}}(-1)^{p_{X}}\hat{\Gamma}^{(1)}_{Xl}(k).$
(160)

The propagator of second order, $n=2$, is given by

$\hat{\Gamma}^{(2)}_{Xlm}(\bm{k}_{1},\bm{k}_{2})=\sum_{l_{1},l_{2}}\hat{\Gamma}^{(2)\,l}_{Xl_{1}l_{2}}(k_{1},k_{2})X^{l_{1}l_{2}}_{lm}(\hat{\bm{k}}_{1},\hat{\bm{k}}_{2}).$
(161)

with the invariant propagator $\hat{\Gamma}^{(2)}_{Xl_{1}l_{2}}(k_{1},k_{2})$
which is a real function. The parity transformation of the invariant function
is given by

$\hat{\Gamma}^{(2)\,l}_{Xl_{1}l_{2}}(k_{1},k_{2})\xrightarrow{\mathbb{P}}(-1)^{p_{X}+l+l_{1}+l_{2}}\hat{\Gamma}^{(2)\,l}_{Xl_{1}l_{2}}(k_{1},k_{2}).$
(162)

If the Universe is statistically invariant under the parity, the function is
nonzero only when $p_{X}+l+l_{1}+l_{2}=\mathrm{even}$. The interchange
symmetry for the second-order propagator is given by

$\hat{\Gamma}^{(2)\,l}_{Xl_{2}l_{1}}(k_{2},k_{1})=(-1)^{l+l_{1}+l_{2}}\hat{\Gamma}^{(2)\,l}_{Xl_{1}l_{2}}(k_{1},k_{2}).$
(163)

The propagator of third order, $n=3$, is given by

$\hat{\Gamma}^{(3)}_{Xlm}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})=\sum_{l_{1},l_{2},l_{3},L}\hat{\Gamma}^{(3)\,l;L}_{Xl_{1}l_{2}l_{3}}(k_{1},k_{2},k_{3})\\\
\times
X^{l_{1}l_{2}l_{3}}_{L;lm}(\hat{\bm{k}}_{1},\hat{\bm{k}}_{2},\hat{\bm{k}}_{3}),$
(164)

with the invariant propagator
$\hat{\Gamma}^{(3)\,l;L}_{Xl_{1}l_{2}l_{3}}(k_{1},k_{2},k_{3})$ which is a
real function. The parity transformation of the invariant function is given by

$\hat{\Gamma}^{(3)\,l;L}_{Xl_{1}l_{2}l_{3}}(k_{1},k_{2},k_{3})\xrightarrow{\mathbb{P}}(-1)^{p_{X}+l+l_{1}+l_{2}+l_{3}}\hat{\Gamma}^{(3)\,l;L}_{Xl_{1}l_{2}l_{3}}(k_{1},k_{2},k_{3}).$
(165)

If the Universe is statistically invariant under the parity, the function is
nonzero only when $p_{X}+l+l_{1}+l_{2}+l_{3}=\mathrm{even}$. The interchange
symmetries for the third-order propagator are given by

$\hat{\Gamma}^{(3)l;L}_{Xl_{l}l_{3}l_{2}}(k_{1},k_{3},k_{2})=(-1)^{l_{2}+l_{3}+L}\hat{\Gamma}^{(3)l;L}_{Xl_{l}l_{2}l_{3}}(k_{1},k_{2},k_{3}),$
(166)

and

$\hat{\Gamma}^{(3)l;L}_{Xl_{2}l_{1}l_{3}}(k_{2},k_{1},k_{3})=(-1)^{l_{1}+l_{2}}\sum_{L^{\prime}}\\{L^{\prime}\\}\begin{Bmatrix}l_{1}&l&L\\\
l_{2}&l_{3}&L^{\prime}\end{Bmatrix}\\\
\times\hat{\Gamma}^{(3)l;L^{\prime}}_{Xl_{l}l_{2}l_{3}}(k_{1},k_{2},k_{3}).$
(167)

The propagators of order $n>3$ are given by

$\hat{\Gamma}^{(n)}_{Xlm}(\bm{k}_{1},\cdots,\bm{k}_{n})=\sum_{\begin{subarray}{c}l_{1},\ldots,l_{n}\\\
L_{2},\ldots,L_{n-1}\end{subarray}}\hat{\Gamma}^{(n)\,l;L_{2}\cdots
L_{n-1}}_{Xl_{1}\cdots l_{n}}(k_{1},\ldots,k_{n})\\\ \times X^{l_{1}\cdots
l_{n}}_{L_{2}\cdots L_{n-1};lm}(\hat{\bm{k}}_{1},\ldots,\hat{\bm{k}}_{n}),$
(168)

with the invariant propagators $\hat{\Gamma}^{(n)\,l;L_{2}\cdots
L_{n-1}}_{Xl_{1}\cdots l_{n}}(k_{1},\ldots,k_{n})$ which are real functions.
The parity transformation is given by

$\hat{\Gamma}^{(n)\,l;L_{2}\cdots L_{n-1}}_{Xl_{1}\cdots
l_{n}}(k_{1},\ldots,k_{3})\\\
\xrightarrow{\mathbb{P}}(-1)^{p_{X}+l+l_{1}+\cdots+l_{n}}\hat{\Gamma}^{(n)\,l;L_{2}\cdots
L_{n-1}}_{Xl_{1}\cdots l_{n}}(k_{1},\ldots,k_{3}).$ (169)

The interchange symmetries for the higher-order propagators are given by

$\hat{\Gamma}^{(n)l;L_{2}\cdots L_{n-2}L_{n-1}}_{Xl_{l}\cdots
l_{n-2}l_{n}l_{n-1}}(k_{1},\ldots,k_{n-2},k_{n},k_{n-1})\\\
=(-1)^{l_{n-1}+l_{n}+L_{n-1}}\hat{\Gamma}^{(n)l;L_{2}\cdots
L_{n-1}}_{Xl_{l}\cdots l_{n}}(k_{1},\ldots,k_{n}),$ (170)

and

$\hat{\Gamma}^{(n)\,l;L_{2}\cdots L_{i-1}L_{i+1}L_{i}L_{i+2}\cdots
L_{n-1}}_{Xl_{1}\cdots l_{i-1}l_{i+1}l_{i}l_{i+2}\cdots
l_{n}}(k_{1},\ldots,k_{i-1},k_{i+1},k_{i},k_{i+2},\ldots,k_{n})\\\
=(-1)^{l_{i}+l_{i+1}}\sum_{L^{\prime}}\\{L^{\prime}\\}\begin{Bmatrix}l_{i}&L_{i}&L_{i+1}\\\
l_{i+1}&L_{i+2}&L^{\prime}\end{Bmatrix}\\\
\times\hat{\Gamma}^{(n)\,l;L_{2}\cdots L_{i}L^{\prime}L_{i+2}\cdots
L_{n-1}}_{Xl_{1}\cdots l_{n}}(k_{1},\ldots,k_{n}).$ (171)

Symmetric property of all the other permutations are obtained by combining the
above two cases.

####  III.4.2  Redshift space

In redshift space, the propagators also depend on the direction of the line of
sight, $\hat{\bm{z}}$. In this paper, we approximately consider the direction
of the line of sight is fixed in the observed space. This approximation is
usually called the plane-parallel, or distant-observer approximation, and is
known as a good approximation for most of the practical purposes. In order to
keep the rotational covariance in the theory, the line of sight,
$\hat{\bm{z}}$, is considered to be arbitrarily directed in the coordinates
system. In literature, it is a common practice to fix the direction of line of
sight in the third axis $\hat{\bf{e}}_{3}$ in the distant-observer
approximation. In this paper, we do not fix the direction, and we have
$\hat{\bm{z}}\neq\hat{\mathbf{e}}_{3}$ in general.

In redshift space, the propagators depend on the direction of the line of
sight $\hat{\bm{z}}$ as well as the wavevectors. We decompose the directional
dependencies on directions of both the line of sight $\hat{\bm{z}}$ and
wavevectors $\hat{\bm{k}}_{1},\ldots,\hat{\bm{k}}_{n}$ by spherical harmonics
in a similar manner of Eq. (107) in the case of renormalized bias functions.
However, it is often useful not to expand all the angular dependencies, and
keeping the angular dependencies in the sum of the wavevectors
$\bm{k}=\bm{k}_{1}+\cdots+\bm{k}_{n}$ unexpanded. The sum of wavevectors
corresponds to the external wavevector [17], and is not integrated over when
the predictions of nonlinear spectra are calculated. Therefore, the amplitude
$k=|\bm{k}|$ and directional cosine with respect to the line of sight
$\mu=\hat{\bm{z}}\cdot\hat{\bm{k}}$ can be fixed in the calculation of
propagators in redshift space. Some part of the directional dependencies of
the propagators on the line of sight $\hat{\bm{z}}$ and the wavevectors
$\hat{\bm{k}}_{1},\hat{\bm{k}}_{2},\ldots$, which is included in the
parameters $k$ and $\mu$ are optionally not need to be expanded in the
spherical harmonics. While it is possible to expand every angular dependencies
into spherical harmonics, the calculations are more complicated than keeping
the dependencies in $k$ and $\mu$ unexpanded as much as possible. In the
following, we consider the general case that the propagators conveniently have
the explicit dependence on $k$ and $\mu$.

The expansion is given by

$\hat{\Gamma}^{(n)}_{Xlm}(\bm{k}_{1},\ldots\bm{k}_{n};\hat{\bm{z}};k,\mu)\\\
=\sum_{l_{z},l_{1},\ldots,l_{n}}\hat{\Gamma}^{(n)\,l_{z};l_{1}\cdots
l_{n}}_{Xlm;m_{z};m_{1}\cdots m_{n}}(k_{1},\ldots,k_{n};k,\mu)\\\ \times
C_{l_{z}m_{z}}^{*}(\hat{\bm{z}})C_{l_{1}m_{1}}^{*}(\hat{\bm{k}}_{1})\cdots
C_{l_{n}m_{n}}^{*}(\hat{\bm{k}}_{n}).$ (172)

As the parameters $k$ and $\mu$ are functions of $\hat{\bm{z}}$ and
$\bm{k}_{1},\ldots,\bm{k}_{n}$, the above expansion is not unique in general.
Nevertheless, when which part of the dependencies on the parameters $k$ and
$\mu$ is consistently fixed in the expression on the lhs of Eq. (172)
throughout, the expansion is unique, at least artificially.

Using the expansion of Eq. (172), deriving the invariant propagators in
redshift space is quite similar to the previous cases of renormalized bias
functions and propagators in real space. Comparing Eq. (172) with Eq. (108),
the similarity is apparent, while there are additional dependencies on the
direction of line of sight, $\hat{\bm{z}}$, and also possible dependencies on
$k$ and $\mu$. Bearing these differences in mind, we apply the similar
derivation in Sec. III.3 to the present case below.

The transformation under the passive rotation of Eq. (20) is given by

$\hat{\Gamma}^{(n)\,l_{z};l_{1}\cdots l_{n}}_{Xlm;m_{z};m_{1}\cdots
m_{n}}(k_{1},\ldots,k_{n};k,\mu)\\\
\xrightarrow{\mathbb{R}}\hat{\Gamma}^{(n)\,l_{z};l_{1}\cdots
l_{n}}_{Xlm^{\prime};m_{z}^{\prime};m_{1}^{\prime}\cdots
m_{n}^{\prime}}(k_{1},\ldots,k_{n};k,\mu)D_{(l)m}^{m^{\prime}}(R)\\\ \times
D_{(l_{z})m_{z}}^{m_{z}^{\prime}}(R)D_{(l_{1})m_{1}}^{m_{1}^{\prime}}(R)\cdots
D_{(l_{n})m_{n}}^{m_{n}^{\prime}}(R).$ (173)

When we rotate the direction of the line of sight together with the coordinate
system, the propagators should be invariant under the rotation of the
coordinate system in the Universe which satisfies the cosmological principle.
Thus we have

$\hat{\Gamma}^{(n)\,l_{z};l_{1}\cdots l_{n}}_{Xlm;m_{z};m_{1}\cdots
m_{n}}(k_{1},\ldots,k_{n};k,\mu)\\\ =\hat{\Gamma}^{(n)\,l_{z};l_{1}\cdots
l_{n}}_{Xlm^{\prime};m_{z}^{\prime};m_{1}^{\prime}\cdots
m_{n}^{\prime}}(k_{1},\ldots,k_{n};k,\mu)\frac{1}{8\pi^{2}}\int[dR]D_{(l)m}^{m^{\prime}}(R)\\\
\times
D_{(l_{z})m_{z}}^{m_{z}^{\prime}}(R)D_{(l_{1})m_{1}}^{m_{1}^{\prime}}(R)\cdots
D_{(l_{n})m_{n}}^{m_{n}^{\prime}}(R).$ (174)

Following the same procedure below Eq. (110) in the case of renormalized bias
functions, one can represent the dependence of propagators in redshift space
by spherical metric, $3j$-symbols, and invariant variables. The main
difference is that an angular dependence on the line of sight and dependencies
on $k$ and $\mu$ additionally appear in the expansion.

For the propagator of first order, $n=1$, we have $k=k_{1}$ and one can omit
the dependence on $k_{1}$. We thus have

$\hat{\Gamma}^{(1)\,l_{z};l_{1}}_{Xlm;m_{z};m_{1}}(k,\mu)=\left(l\,l_{z}\,l_{1}\right)_{mm_{z}m_{1}}\hat{\Gamma}^{(1)\,l\,l_{z}}_{Xl_{1}}(k,\mu),$
(175)

with

$\hat{\Gamma}^{(1)l\,l_{z}}_{Xl_{1}}(k,\mu)\equiv(-1)^{l+l_{z}+l_{1}}\left(l\,l_{z}\,l_{1}\right)^{mm_{z}m_{1}}\hat{\Gamma}^{(1)\,l_{z};l_{1}}_{Xlm;m_{z};m_{1}}(k,\mu).$
(176)

The propagator is given by

$\hat{\Gamma}^{(1)}_{Xlm}(\bm{k};\hat{\bm{z}};k,\mu)=\sum_{l_{z},l_{1}}\hat{\Gamma}^{(1)\,l\,l_{z}}_{Xl_{1}}(k,\mu)X^{l_{z}l_{1}}_{lm}(\hat{\bm{z}},\hat{\bm{k}}).$
(177)

Combining the properties of complex conjugate for the bipolar spherical
harmonics, Eq. (387), and for the propagator in redshift space, Eq. (88) with
$n=1$, we have

$\hat{\Gamma}^{(1)l\,l_{z}\,*}_{Xl_{1}}(k,\mu)=(-1)^{l_{z}}\hat{\Gamma}^{(1)l\,l_{z}}_{Xl_{1}}(k,-\mu).$
(178)

Combining the properties of parity for the bipolar spherical harmonics, Eq.
(386), and for the propagator in redshift space, Eq. (101) with $n=1$, we have

$\hat{\Gamma}^{(1)l\,l_{z}}_{Xl_{1}}(k,\mu)\xrightarrow{\mathbb{P}}(-1)^{p_{X}+l+l_{z}+l_{1}}\hat{\Gamma}^{(1)l\,l_{z}}_{Xl_{1}}(k,\mu).$
(179)

If the Universe is statistically invariant under the parity, the function is
nonzero only when $p_{X}+l+l_{z}+l_{1}=\mathrm{even}$.

For the propagator of second order, $n=2$, we have

$\hat{\Gamma}^{(2)\,l_{z};l_{1}l_{2}}_{Xlm;m_{z};m_{1}m_{2}}(k_{1},k_{2};k,\mu)=\sum_{L}(-1)^{L}\sqrt{\\{L\\}}\left(l\,l_{z}\,L\right)_{mm_{z}M}\\\
\times\left(L\,l_{1}\,l_{2}\right)^{M}_{\phantom{M}m_{1}m_{2}}\hat{\Gamma}^{(2)\,l\,l_{z};L}_{Xl_{1}l_{2}}(k_{1},k_{2};k,\mu),$
(180)

with

$\hat{\Gamma}^{(2)\,l\,l_{z};L}_{Xl_{1}l_{2}}(k_{1},k_{2};k,\mu)\equiv(-1)^{l+l_{z}+l_{1}+l_{2}+L}\sqrt{\\{L\\}}\left(l\,l_{z}\,L\right)^{mm_{z}}_{\phantom{mm_{z}}M}\\\
\times\left(L\,l_{1}\,l_{2}\right)^{Mm_{1}m_{2}}\hat{\Gamma}^{(2)\,l\,l_{z};l_{1}l_{2}}_{Xmm_{z};m_{1}m_{2}}(k_{1},k_{2};k,\mu).$
(181)

The propagator is given by

$\hat{\Gamma}^{(2)}_{Xlm}(\bm{k}_{1},\bm{k}_{2};\hat{\bm{z}};k,\mu)\\\
=\sum_{l_{z},l_{1},l_{2},L}\hat{\Gamma}^{(2)\,l\,l_{z};L}_{Xl_{1}l_{2}}(k_{1},k_{2};k,\mu)X^{l_{z}l_{1}l_{2}}_{L;lm}(\hat{\bm{z}},\hat{\bm{k}}_{1},\hat{\bm{k}}_{2}).$
(182)

Combining the properties of complex conjugate for the tripolar spherical
harmonics, Eq. (394), and for the propagator in redshift space, Eq. (88) with
$n=2$, we have

$\hat{\Gamma}^{(2)\,l\,l_{z};L\,*}_{Xl_{1}l_{2}}(k_{1},k_{2};k,\mu)=(-1)^{l_{z}}\hat{\Gamma}^{(2)\,l\,l_{z};L}_{Xl_{1}l_{2}}(k_{1},k_{2};k,-\mu),$
(183)

Combining the properties of parity for the tripolar spherical harmonics, Eq.
(393), and for the propagator in redshift space, Eq. (101) with $n=2$, we have

$\hat{\Gamma}^{(2)\,l\,l_{z};L}_{Xl_{1}l_{2}}(k_{1},k_{2};k,\mu)\xrightarrow{\mathbb{P}}(-1)^{p_{X}+l+l_{1}+l_{2}}\hat{\Gamma}^{(2)\,l\,l_{z};L}_{Xl_{1}l_{2}}(k_{1},k_{2};k,\mu),$
(184)

If the Universe is statistically invariant under the parity, the function is
nonzero only when $p_{X}+l+l_{z}+l_{1}+l_{2}=\mathrm{even}$. The interchange
symmetry for the second-order propagator is given by

$\hat{\Gamma}^{(2)\,l\,l_{z};L}_{Xl_{2}l_{1}}(k_{2},k_{1};k,\mu)=(-1)^{l_{1}+l_{2}+L}\hat{\Gamma}^{(2)\,l\,l_{z};L}_{Xl_{1}l_{2}}(k_{1},k_{2};k,\mu).$
(185)

For the propagator of higher orders, $n\geq 3$, we have

$\hat{\Gamma}^{(n)\,l_{z};l_{1}\cdots l_{n}}_{Xlm;m_{z};m_{1}\cdots
m_{n}}(k_{1},\ldots,k_{n};k,\mu)\\\
=\sum_{L_{1},\ldots,L_{n-1}}(-1)^{L_{1}+\cdots+L_{n-1}}\sqrt{\\{L_{1}\\}\cdots\\{L_{n-1}\\}}\left(l\,l_{z}\,L_{1}\right)_{mm_{z}M_{1}}\\\
\times\left(L_{1}\,l_{1}\,L_{2}\right)^{M_{1}}_{\phantom{M_{1}}m_{1}M_{2}}\cdots\left(L_{n-2}\,l_{n-2}\,L_{n-1}\right)^{M_{n-2}}_{\phantom{M_{n-2}}m_{n-2}M_{n-1}}\\\
\times\left(L_{n-1}\,l_{n-1}\,l_{n}\right)^{M_{n-1}}_{\phantom{M_{n-1}}m_{n-1}m_{n}}\hat{\Gamma}^{(n)\,l\,l_{z};L_{1}\cdots
L_{n-1}}_{Xl_{1}\cdots l_{n}}(k_{1},\ldots,k_{n};k,\mu),$ (186)

with

$\hat{\Gamma}^{(n)\,l\,l_{z};L_{1}\cdots L_{n-1}}_{Xl_{1}\cdots
l_{n}}(k_{1},\ldots,k_{n};k,\mu)\equiv(-1)^{l+l_{z}+l_{1}+\cdots+l_{n}}\\\
\times(-1)^{L_{1}+\cdots+L_{n-1}}\sqrt{\\{L_{1}\\}\cdots\\{L_{n-1}\\}}\left(l\,l_{z}\,L_{1}\right)^{mm_{z}}_{\phantom{mm_{z}}M_{1}}\\\
\times\left(L_{1}\,l_{1}\,L_{2}\right)^{M_{1}m_{1}}_{\phantom{M_{1}m_{1}}M_{2}}\cdots\left(L_{n-2}\,l_{n-2}\,L_{n-1}\right)^{M_{n-2}m_{n-2}}_{\phantom{M_{n-2}m_{n-2}}M_{n-1}}\\\
\times\left(L_{n-1}\,l_{n-1}\,l_{n}\right)^{M_{n-1}m_{n-1}m_{n}}\hat{\Gamma}^{(n)\,l_{z};l_{1}\cdots
l_{n}}_{Xlm;m_{z};m_{1}\cdots m_{n}}(k_{1},\ldots,k_{n};k,\mu),$ (187)

The propagator is given by

$\hat{\Gamma}^{(n)}_{Xlm}(\bm{k}_{1},\cdots,\bm{k}_{n};\hat{\bm{z}};k,\mu)\\\
=\sum_{\begin{subarray}{c}l_{z},l_{1},\ldots,l_{n}\\\
L_{1},\ldots,L_{n-1}\end{subarray}}\hat{\Gamma}^{(n)\,l\,l_{z};L_{1}\cdots
L_{n-1}}_{Xl_{1}\cdots l_{n}}(k_{1},\cdots,k_{n};k,\mu)\\\ \times
X^{l_{z}l_{1}\cdots l_{n}}_{L_{1}\cdots
L_{n-1};lm}(\hat{\bm{z}},\hat{\bm{k}}_{1},\cdots,\hat{\bm{k}}_{n}).$ (188)

Combining the properties of complex conjugate for the polypolar spherical
harmonics, Eq. (400), and for the propagator in redshift space, Eq. (88), we
have

$\hat{\Gamma}^{(n)\,l\,l_{z};L_{1}\cdots L_{n-1}\,*}_{Xl_{1}\cdots
l_{n}}(k_{1},\ldots,k_{n};k,\mu)\\\
=(-1)^{l_{z}}\hat{\Gamma}^{(n)\,l\,l_{z};L_{1}\cdots L_{n-1}}_{Xl_{1}\cdots
l_{n}}(k_{1},\ldots,k_{n};k,-\mu).$ (189)

Combining the properties of parity for the polypolar spherical harmonics, Eq.
(399), and for the propagator in redshift space, Eq. (101), we have

$\hat{\Gamma}^{(n)\,l\,l_{z};L_{1}\cdots L_{n-1}}_{Xl_{1}\cdots
l_{n}}(k_{1},\ldots,k_{n};k,\mu)\\\
\xrightarrow{\mathbb{P}}(-1)^{p_{X}+l+l_{1}+\cdots+l_{n}}\hat{\Gamma}^{(n)\,l\,l_{z};L_{1}\cdots
L_{n-1}}_{Xl_{1}\cdots l_{n}}(k_{1},\ldots,k_{n};k,\mu).$ (190)

If the Universe is statistically invariant under the parity, the function is
nonzero only when $p_{X}+l+l_{z}+l_{1}+\cdots+l_{n}=\mathrm{even}$.
Corresponding to Eqs. (157) and (158) of the renormalized bias function, the
interchange symmetries for the higher-order propagators are given by

$\hat{\Gamma}^{(n)l\,l_{z};L_{1}\cdots L_{n-2}L_{n-1}}_{Xl_{l}\cdots
l_{n-2}l_{n}l_{n-1}}(k_{1},\ldots,k_{n-2},k_{n},k_{n-1};k,\mu)\\\
=(-1)^{l_{n-1}+l_{n}+L_{n-1}}\hat{\Gamma}^{(n)l\,l_{z};L_{1}\cdots
L_{n-1}}_{Xl_{l}\cdots l_{n}}(k_{1},\ldots,k_{n};k,\mu),$ (191)

and

$\hat{\Gamma}^{(n)\,l\,l_{z};L_{1}\cdots L_{i+1}L_{i}\cdots
L_{n-1}}_{Xl_{1}\cdots l_{i-1}l_{i+1}l_{i}l_{i+2}\cdots
l_{n}}(k_{1},\cdots,k_{i+1},k_{i},\cdots,k_{n};k,\mu)\\\
=(-1)^{l_{i}+l_{i+1}}\sum_{L^{\prime}}(2L^{\prime}+1)\begin{Bmatrix}l_{i}&L_{i}&L_{i+1}\\\
l_{i+1}&L_{i+2}&L^{\prime}\end{Bmatrix}\\\
\times\hat{\Gamma}^{(n)\,l\,l_{z};L_{1}\cdots L_{i}L^{\prime}L_{i+2}\cdots
L_{n-1}}_{Xl_{1}\cdots l_{n}}(k_{1},\ldots,k_{n};k,\mu),$ (192)

where $i\geq 2$. In the invariant propagator on the lhs of Eq. (192), the
order of only two arguments $k_{i}$ and $k_{i+1}$ is reversed.

###  III.5 Sample calculations of decomposed propagators

One can straightforwardly derive explicit forms of reduced propagators with
the formal decomposition as explained above. A simple way to derive the
invariant coefficients is to represent the angular dependence of the
propagators in terms of the spherical harmonics and polypolar spherical
harmonics, and read off the coefficients from the resulting expressions.

####  III.5.1 Tree-level first-order propagator in real space

The simplest example is the tree-level approximation of first-order propagator
of Eq. (75):

$\hat{\Gamma}_{Xlm}^{(1)}(\bm{k})=c_{Xlm}^{(1)}(\bm{k})+\left[\bm{k}\cdot\bm{L}_{1}(\bm{k})\right]c_{Xlm}^{(0)}.$
(193)

In real space, $\bm{L}_{1}(\bm{k})$ is given by Eq. (78), and we have
$\bm{k}\cdot\bm{L}_{1}(\bm{k})=1$. The relevant renormalized bias functions
are represented by invariant quantities through Eqs. (120) and (126).
Therefore the propagator of Eq. (193) is represented in a form of expansion
with spherical harmonics as

$\hat{\Gamma}_{Xlm}^{(1)}(\bm{k})=\frac{(-1)^{l}}{\sqrt{2l+1}}\left[c^{(1)}_{Xl}(k)+\delta_{l0}c^{(0)}_{X}\right]C_{lm}(\hat{\bm{k}}).$
(194)

Comparing this expression with Eq. (159), we have

$\hat{\Gamma}^{(1)}_{Xl}(k)=c^{(1)}_{Xl}(k)+\delta_{l0}c^{(0)}_{X}.$ (195)

In the Universe with parity symmetry, the rhs of Eq. (195) is nonzero only for
normal tensors, $p_{X}=0$, as discussed in Sec. III.3 and also the lhs
consistently vanishes as discussed in Sec. III.4.

####  III.5.2 Tree-level first-order propagator in redshift space with the
complete expansion

In redshift space, the tree-level approximation of the first-order propagator
is given by

$\hat{\Gamma}_{Xlm}^{(1)}(\bm{k};\hat{\bm{z}})=c_{Xlm}^{(1)}(\bm{k})+\left[\bm{k}\cdot\bm{L}^{\mathrm{s}}_{1}(\bm{k};\hat{\bm{z}})\right]c_{Xlm}^{(0)},$
(196)

where the displacement kernel is replaced by the one in redshift space as in
Eq. (82), and the relevant factor in the above equation is given by

$\bm{k}\cdot\bm{L}^{\mathrm{s}}_{1}(\bm{k};\hat{\bm{z}})=1+f(\hat{\bm{z}}\cdot\hat{\bm{k}})^{2}=1+\frac{f}{3}+\frac{2f}{3}\mathit{P}_{2}(\hat{\bm{z}}\cdot\hat{\bm{k}}),$
(197)

where $\mathit{P}_{2}(x)$ is one of the Legendre polynomials
$\mathit{P}_{l}(x)$ with $l=2$.

First we explain step-by-step calculations to obtain the invariant propagator
in the case that the dependence of $\mu$ is completely expanded into spherical
harmonics. The Legendre polynomials
$\mathit{P}_{l}(\hat{\bm{z}}\cdot\hat{\bm{k}})$ can be decomposed into a
superposition of products of spherical harmonics by the addition theorem of
Eq. (378), and represented by a special case of bipolar spherical harmonics as

$\mathit{P}_{l}(\hat{\bm{z}}\cdot\hat{\bm{k}})=(-1)^{l}\sqrt{2l+1}\,X^{ll}_{00}(\hat{\bm{z}},\hat{\bm{k}}),$
(198)

where the last function on the rhs is a special case of normalized bipolar
spherical harmonics given by Eq. (132).

Combining Eqs. (196)–(198) above, we derive

$\hat{\Gamma}^{(1)}_{Xlm}(\bm{k};\hat{\bm{z}})=\frac{(-1)^{l}}{\sqrt{2l+1}}c^{(1)}_{Xl}(k)C_{lm}(\hat{\bm{k}})\\\
+\delta_{l0}\delta_{m0}c^{(0)}_{X}\left[1+\frac{f}{3}+\frac{2f}{3}\sqrt{5}\,X^{22}_{00}(\hat{\bm{z}},\hat{\bm{k}})\right].$
(199)

Each term of the above equation can be represented by special cases of bipolar
spherical harmonics, because we have

$\displaystyle\delta_{l0}\delta_{m0}$
$\displaystyle=\delta_{l0}\,X^{00}_{lm}(\hat{\bm{z}},\hat{\bm{k}}),$ (200)
$\displaystyle C_{lm}(\hat{\bm{k}})$
$\displaystyle=(-1)^{l}\sqrt{2l+1}\,X^{0l}_{lm}(\hat{\bm{z}},\hat{\bm{k}}),$
(201)

as straightforwardly shown from the definitions. One can also rewrite the
bipolar spherical harmonics appeared in Eq. (199) as

$\delta_{l0}\delta_{m0}X^{22}_{00}(\hat{\bm{z}},\hat{\bm{k}})=\delta_{l0}X^{22}_{lm}(\hat{\bm{z}},\hat{\bm{k}}).$
(202)

Substituting the above identities into Eq. (199), one can readily read off the
coefficient of the bipolar spherical harmonics
$X^{l\,l_{z}}_{lm}(\hat{\bm{z}},\hat{\bm{k}})$ of Eq. (177) as

$\hat{\Gamma}^{(1)l\,l_{z}}_{Xl_{1}}(k)=\delta_{l_{z}0}\delta_{l_{1}l}\,c^{(1)}_{Xl}(k)\\\
+\delta_{l0}c^{(0)}_{X}\left[\left(1+\frac{f}{3}\right)\delta_{l_{z}0}\delta_{l_{1}0}+\frac{2\sqrt{5}f}{3}\delta_{l_{z}2}\delta_{l_{1}2}\right].$
(203)

Putting $f=0$, $l_{z}=0$ and $l_{1}=l$ in the above equation consistently
recovers the result in real space of Eq. (195). In the Universe with parity
symmetry, the rhs of Eq. (203) only nonzero when $p_{X}=0$ as discussed in
Sec. III.3, and apparently the rhs is nonzero only for
$l+l_{z}+l_{1}=\mathrm{even}$. Thus the lhs is nonzero only when
$p_{X}+l+l_{z}+l_{1}=\mathrm{even}$, which is consistent with the discussion
in Sec. III.4. The same result of Eq. (203) can also be derived more
straightforwardly by using an orthonormality relation of the bipolar spherical
harmonics given in Eq. (388). According to this orthonormality relation, the
Eq. (177) is inverted as

$\hat{\Gamma}^{(1)l\,l_{z}}_{Xl_{1}}(k)=(-1)^{l+l_{z}+l_{1}}\\{l\\}\\{l_{z}\\}\\{l_{1}\\}\int\frac{d^{2}\hat{z}}{4\pi}\frac{d^{2}\hat{k}}{4\pi}\,\hat{\Gamma}^{(1)}_{Xl0}(\bm{k};\hat{\bm{z}})X^{l_{z}l_{1}}_{l0}(\hat{\bm{z}},\hat{\bm{k}}).$
(204)

Substituting Eq. (199) into the above equation, using Eqs. (200)–(202) with
the orthonormal relation of Eq. (388), we derive the same result of Eq. (203)
as well.

####  III.5.3 Tree-level first-order propagator in redshift space with a
partial expansion

In the above result of Eq. (203) for the first-order invariant propagator in
redshift space is obtained by completely expanding the angular dependencies in
the propagator. As explained in Sec. III.4.2, angular dependencies which are
included in the external wavevector can optionally be kept unexpanded to make
the following calculations simpler. In the case of the tree-level first-order
propagator, the argument of the propagator in Eq. (196) is itself the external
wavevector.

Combining Eqs. (196) and (197), we trivially have

$\hat{\Gamma}_{Xlm}^{(1)}(\bm{k})=c_{Xlm}^{(1)}(\bm{k})+\left(1+f\mu^{2}\right)c_{Xlm}^{(0)}.$
(205)

Substituting Eqs. (120), (126), (200) (201) into the above equation, the
invariant propagator can be read off from Eq. (177), resulting in

$\hat{\Gamma}^{(1)l\,l_{z}}_{Xl_{1}}(k,\mu)=\delta_{l_{z}0}\delta_{l_{1}l}\left[c^{(1)}_{Xl}(k)+\delta_{l0}(1+f\mu^{2})c^{(0)}_{X}\right].$
(206)

The same result is also derived by substituting Eq. (205) into Eq. (204),
formally regarding as if the variable $\mu$ does not explicitly depend on
$\hat{\bm{k}}$ and $\hat{\bm{z}}$.

####  III.5.4 Tree-level second-order propagator in real space

For another example of calculating the invariant propagators, we consider the
second-order propagator $\hat{\Gamma}^{(2)}_{Xlm}$ in real space, Eq. (161).
The calculation in this case is quite similar to the derivation of Eq. (203).
With the tree-level approximation, the second-order propagator is given by Eq.
(76). Applying a similar procedure as in the above derivation using Eqs. (198)
and (200)–(202), the rhs of Eq. (76) is represented in a form of the rhs of
Eq. (161), and we can read off the invariant coefficient from the resulting
expression.

In the course of calculation, one needs to represent a product of spherical
harmonics of the same argument through the formula given by Eq. (382). The
coefficient of the spherical harmonics in this formula is known as the Gaunt
integral. Similarly to the simplified notation for the $3j$-symbols we
introduced above, it is useful to use a simplified symbol for the Gaunt
integral $[l_{1}\,l_{2}\,l_{3}]_{m_{1}m_{2}m_{3}}$ defined by Eq. (411) in our
notation. We consider the azimuthal indices in the Gaunt integral can be
raised and lowered by the spherical metric, such as

$\left[l_{1}\,l_{2}\,l_{3}\right]^{m_{1}m_{2}}_{\phantom{m_{1}m_{2}}m_{3}}=g_{(l_{1})}^{m_{1}m_{1}^{\prime}}g_{(l_{2})}^{m_{2}m_{2}^{\prime}}\left[l_{1}\,l_{2}\,l_{3}\right]_{m_{1}^{\prime}m_{2}^{\prime}m_{3}},$
(207)

and so forth. The Gaunt integral of Eq. (411) is nonzero only when
$l_{1}+l_{2}+l_{3}$ is an even number. With this notation, Eq. (382) is simply
represented by

$C_{l_{1}m_{1}}(\hat{\bm{k}})\,C_{l_{2}m_{2}}(\hat{\bm{k}})=\sum_{l}(2l+1)\left[l_{1}\,l_{2}\,l\right]_{m_{1}m_{2}}^{\phantom{m_{1}m_{2}}m}C_{lm}(\hat{\bm{k}}).$
(208)

Properties of the Gaunt integral, Eqs. (411)–(416) are useful in analytical
deductions of invariant propagators.

Following the procedures above to derive the invariant propagator of first
order, and additionally using the formula of Eq. (208), the invariant
propagator of second order can be straightforwardly derived. The result in
real space is given by

$\hat{\Gamma}^{(2)l}_{Xl_{1}l_{2}}(k_{1},k_{2})=c^{(2)l}_{Xl_{1}l_{2}}(k_{1},k_{2})\\\
+c^{(1)}_{Xl}(k_{1})\left[\delta_{l_{1}l}\delta_{l_{2}0}+\frac{k_{1}}{k_{2}}\frac{(-1)^{l}}{\sqrt{\\{l\\}}}\,\\{l_{1}\\}\begin{pmatrix}l&l_{1}&1\\\
0&0&0\end{pmatrix}\delta_{l_{2}1}\right]\\\
+c^{(1)}_{Xl}(k_{2})\left[\delta_{l_{1}0}\delta_{l_{2}l}+\frac{k_{2}}{k_{1}}\frac{(-1)^{l}}{\sqrt{\\{l\\}}}\,\\{l_{2}\\}\begin{pmatrix}l&l_{2}&1\\\
0&0&0\end{pmatrix}\delta_{l_{1}1}\right]\\\
+\delta_{l0}c^{(0)}_{X}\left[\frac{34}{21}\delta_{l_{1}0}\delta_{l_{2}0}-\sqrt{3}\left(\frac{k_{1}}{k_{2}}+\frac{k_{2}}{k_{1}}\right)\delta_{l_{1}1}\delta_{l_{2}1}\right.\\\
\left.+\frac{8\sqrt{5}}{21}\delta_{l_{1}2}\delta_{l_{2}2}\right].$ (209)

In the Universe with parity symmetry, the terms on the rhs of Eq. (209) except
the first term is nonzero only for normal tensors, $p_{X}=0$, as discussed in
Sec. III.3, and apparently nonzero only for $l+l_{1}+l_{2}=\mathrm{even}$. The
first term on the rhs is nonzero only for $p_{X}+l+l_{1}+l_{2}=\mathrm{even}$.
The lhs is nonzero only for $p_{X}+l+l_{1}+l_{2}=\mathrm{even}$, and therefore
the properties are consistent with each other. The corresponding result in
redshift space is explicitly given in Paper II [43].

Other propagators to arbitrary orders, both in real space and redshift space
can be calculated in a similar way as illustrated above, although the
calculations are more or less demanding in the case of calculation including
higher orders.

##  IV The spectra of tensor fields in perturbation theory

In this section, we consider how one can calculate the power spectrum and
higher-order spectra such as the bispectrum and so forth, for tensor fields in
general.

###  IV.1 The power spectrum

####  IV.1.1  The power spectrum in real space

The power spectrum $P^{(l_{1}l_{2})}_{X_{1}X_{2}\,m_{1}m_{2}}(\bm{k})$ of
irreducible components of tensor field $F_{Xlm}(\bm{k})$ can be defined by

$\left\langle
F_{X_{1}l_{1}m_{1}}(\bm{k}_{1})F_{X_{2}l_{2}m_{2}}(\bm{k}_{2})\right\rangle_{\mathrm{c}}\\\
=(2\pi)^{3}\delta_{\mathrm{D}}^{3}(\bm{k}_{1}+\bm{k}_{2})P^{(l_{1}l_{2})}_{X_{1}X_{2}m_{1}m_{2}}(\bm{k}_{1}),$
(210)

where $\langle\cdots\rangle_{\mathrm{c}}$ indicates the two-point cumulant, or
the connected part, and the appearance of the delta function is due to
translational symmetry. We consider generally the cross power spectra between
two fields, $X_{1}$ and $X_{2}$ with irreducible components $l_{1}$ and
$l_{2}$, respectively. The auto power spectrum is straightforwardly obtained
by putting $X_{1}=X_{2}=X$.

The power spectrum $P^{(l_{1}l_{2})}_{X_{1}X_{2}m_{1}m_{2}}(\bm{k})$ defined
above depends on the coordinates system, as it has azimuthal indices and
depends on the direction of wavevector $\bm{k}$ in the arguments. One can also
construct the reduced power spectrum which is invariant under the rotation, in
a similar manner to introducing the invariant propagators in Sec. III.4.

First, we consider the power spectrum in real space. The directional
dependence on the wavevector can be expanded by spherical harmonics as

$P^{(l_{1}l_{2})}_{X_{1}X_{2}m_{1}m_{2}}(\bm{k})=\sum_{l}P^{l_{1}l_{2};l}_{X_{1}X_{2}m_{1}m_{2};m}(k)\,C^{*}_{lm}(\hat{\bm{k}}).$
(211)

Due to the rotational transformations of Eqs. (92) and (22), and the
rotational invariance of the delta function in the definition of the power
spectrum, Eq. (210), the rotational transformation of the expansion
coefficient is given by

$P^{l_{1}l_{2};l}_{X_{1}X_{2}m_{1}m_{2};m}(k)\xrightarrow{\mathbb{R}}P^{l_{1}l_{2};l}_{X_{1}X_{2}m_{1}^{\prime}m_{2}^{\prime};m^{\prime}}(k)D^{m_{1}^{\prime}}_{(l_{1})m_{1}}(R)D^{m_{2}^{\prime}}_{(l_{2})m_{2}}(R)D^{m^{\prime}}_{(l)m}(R).$
(212)

For the statistically isotropic Universe, the functional form of the power
spectrum should not depend on the choice of coordinates system. Following the
same procedure as in Sec. III.3, the power spectrum in the statistically
isotropic Universe is shown to have a form,

$P^{(l_{1}l_{2})}_{X_{1}X_{2}m_{1}m_{2}}(\bm{k})=\sum_{l}\left(l_{1}\,l_{2}\,l\right)_{m_{1}m_{2}}^{\phantom{m_{1}m_{2}}m}C_{lm}(\hat{\bm{k}})P^{l_{1}l_{2};l}_{X_{1}X_{2}}(k),$
(213)

where

$P^{l_{1}l_{2};l}_{X_{1}X_{2}}(k)=(-1)^{l_{1}+l_{2}+l}\left(l_{1}\,l_{2}\,l\right)^{m_{1}m_{2}m}P^{l_{1}l_{2};l}_{X_{1}X_{2}m_{1}m_{2};m}(k)$
(214)

is the invariant power spectrum. The physical degrees of freedom of the power
spectrum are represented by the last invariant spectrum,
$P^{l_{1}l_{2};l}_{X_{1}X_{2}}(k)$. There are a finite number of independent
functions of power spectra which depend only on the modulus of the wavevector
$k=|\bm{k}|$ of Fourier modes, due to the triangle inequality of the
$3j$-symbol, $|l_{1}-l_{2}|\leq l\leq l_{1}+l_{2}$ for a fixed set of integers
$(l_{1},l_{2})$. For example, when $l_{1}=l_{2}=2$, only $l=0,1,\ldots,4$ are
allowed and there are five independent components in the power spectrum222In
Ref. [39], the power spectrum $P_{22}^{(q)}(k)$ of intrinsic alignment is
defined in a coordinate-dependent way, where $q=0,\pm 1,\pm 2$. The degree of
freedom of their power spectrum is five and is consistent with ours..

####  IV.1.2 Symmetries in real space

According to Eqs. (85) and (210), the complex conjugate of the power spectrum
is given by

$P^{(l_{1}l_{2})*}_{X_{1}X_{2}m_{1}m_{2}}(\bm{k})=(-1)^{l_{1}+l_{2}}P^{(l_{1}l_{2})m_{1}m_{2}}_{X_{1}X_{2}}(-\bm{k}).$
(215)

Using the property of the $3j$-symbol, Eq. (404), the symmetry of complex
conjugate is shown to be given by

$P^{l_{1}l_{2};l\,*}_{X_{1}X_{2}}(k)=P^{l_{1}l_{2};l}_{X_{1}X_{2}}(k),$ (216)

and thus the invariant power spectrum is a real function. Additionally, the
parity transformation of the power spectrum is given by

$P^{(l_{1}l_{2})}_{X_{1}X_{2}m_{1}m_{2}}(\bm{k})\xrightarrow{\mathbb{P}}(-1)^{p_{X_{1}}+p_{X_{2}}+l_{1}+l_{2}}P^{(l_{1}l_{2})}_{X_{1}X_{2}m_{1}m_{2}}(-\bm{k}),$
(217)

because of Eqs. (98) and (210). Accordingly, the parity transformation of the
invariant power spectrum is given by

$P^{l_{1}l_{2};l}_{X_{1}X_{2}}(k)\xrightarrow{\mathbb{P}}(-1)^{p_{X_{1}}+p_{X_{2}}+l_{1}+l_{2}+l}P^{l_{1}l_{2};l}_{X_{1}X_{2}}(k).$
(218)

In the Universe with parity symmetry, the power spectrum is invariant under
the parity transformation, we thus have

$p_{X_{1}}+p_{X_{2}}+l_{1}+l_{2}+l=\mathrm{even}.$ (219)

Together with the triangle inequality, $|l_{1}-l_{2}|\leq l\leq l_{1}+l_{2}$,
mentioned above, the numbers of independent components in the invariant power
spectrum $P^{l_{1}l_{2};l}_{X_{1}X_{2}}(k)$ are further reduced in the
presence of parity symmetry. For example, when $l_{1}=l_{2}=2$ and
$p_{X_{1}}=p_{X_{2}}=0$, only $l=0,2,4$ are allowed and there are three
independent components in the power spectrum333In Ref. [39], their power
spectrum of intrinsic alignment satisfies $P_{22}^{(q)}(k)=P_{22}^{(-q)}(k)$
in the presence of parity symmetry, and thus only three components, $q=0,1,2$,
are independent. Again, the degree of freedom of theirs is consistent with
ours.. The power spectrum in real space has an interchange symmetry,

$P^{(l_{2}l_{1})}_{X_{2}X_{1}m_{2}m_{1}}(\bm{k})=P^{(l_{1}l_{2})}_{X_{1}X_{2}m_{1}m_{2}}(-\bm{k}),$
(220)

and the corresponding symmetry for the invariant spectrum is given by

$P^{l_{2}l_{1};l}_{X_{2}X_{1}}(k)=(-1)^{l_{1}+l_{2}}P^{l_{1}l_{2};l}_{X_{1}X_{2}}(k).$
(221)

####  IV.1.3 The power spectrum in redshift space

Next, we consider the power spectrum in redshift space. The directional
dependence of the wavevector and the line of sight can be expanded by
spherical harmonics as

$P^{(l_{1}l_{2})}_{X_{1}X_{2}m_{1}m_{2}}(\bm{k};\hat{\bm{z}})=\sum_{l,l_{z}}P^{l_{1}l_{2};l\,l_{z}}_{X_{1}X_{2}m_{1}m_{2};mm_{z}}(k,\mu)C^{*}_{l_{z}m_{z}}(\hat{\bm{z}})C^{*}_{lm}(\hat{\bm{k}}).$
(222)

Corresponding to the discussion about the propagator in Sec. III.4.2, some
part of the dependence on directional cosine
$\mu=\hat{\bm{z}}\cdot\hat{\bm{k}}$ can be kept unexpanded, and we generally
include the argument $\mu$ in the coefficients in the above expansion. When
the dependence on $\mu$ is completely expanded into the spherical harmonics,
the argument $\mu$ is absent, and the expansion of Eq. (222) is unique.
Otherwise, unless which dependence on $\mu$ is specified and fixed, the
expansion is not unique.

Following similar considerations of rotational symmetry in the rest of this
paper, e.g., in deriving Eq. (177), the power spectrum in the statistically
isotropic Universe is shown to have a form,

$P^{(l_{1}l_{2})}_{X_{1}X_{2}m_{1}m_{2}}(\bm{k};\hat{\bm{z}})=\sum_{L}(-1)^{L}\sqrt{\\{L\\}}\left(l_{1}\,l_{2}\,L\right)_{m_{1}m_{2}}^{\phantom{m_{1}m_{2}}M}\\\
\times\sum_{l,l_{z}}X^{l_{z}l}_{LM}(\hat{\bm{z}},\hat{\bm{k}})P^{l_{1}l_{2};l\,l_{z};L}_{X_{1}X_{2}}(k,\mu),$
(223)

where

$P^{l_{1}l_{2};l\,l_{z};L}_{X_{1}X_{2}}(k,\mu)=(-1)^{l_{1}+l_{2}+l+l_{z}+L}\sqrt{\\{L\\}}\\\
\times\left(l_{1}\,l_{2}\,L\right)^{m_{1}m_{2}}_{\phantom{m_{1}m_{2}}M}\left(L\,l\,l_{z}\right)^{Mmm_{z}}P^{l_{1}l_{2};l\,l_{z}}_{X_{1}X_{2}m_{1}m_{2};mm_{z}}(k,\mu)$
(224)

is the invariant power spectrum in redshift space.

####  IV.1.4 Symmetries in redshift space

As we work in the distant-observer approximation, it is reasonable to assume
that the power spectrum is unchanged when the direction of the line of sight
is flipped, $\hat{\bm{z}}\rightarrow-\hat{\bm{z}}$, keeping the physical
system unchanged:

$P^{(l_{1}l_{2})}_{X_{1}X_{2}m_{1}m_{2}}(\bm{k};\hat{\bm{z}})=P^{(l_{1}l_{2})}_{X_{1}X_{2}m_{1}m_{2}}(\bm{k};-\hat{\bm{z}}).$
(225)

On one hand, when the directional dependence of
$\mu=\hat{\bm{z}}\cdot\hat{\bm{k}}$ is completely expanded in the expansion of
Eq. (223), the invariant spectrum satisfies
$P^{l_{1}l_{2};l\,l_{z};L}_{X_{1}X_{2}}(k)=(-1)^{l_{z}}P^{l_{1}l_{2};l\,l_{z};L}_{X_{1}X_{2}}(k)$,
and thus the integer $l_{z}$ in Eq. (223) should be even numbers. On the other
hand, when the dependence of $\mu$ is included in the invariant spectrum as
$P^{l_{1}l_{2};l\,l_{z};L}_{X_{1}X_{2}}(k,\mu)$, the integers $l_{z}$ are not
necessarily even numbers, and the invariant spectrum instead satisfies

$P^{l_{1}l_{2};l\,l_{z};L}_{X_{1}X_{2}}(k,\mu)=(-1)^{l_{z}}P^{l_{1}l_{2};l\,l_{z};L}_{X_{1}X_{2}}(k,-\mu).$
(226)

According to Eqs. (85), (210) and the symmetry of the $3j$-symbol, Eq. (404),
the complex conjugate of the invariant function is shown to be given by

$P^{l_{1}l_{2};l\,l_{z};L\,*}_{X_{1}X_{2}}(k,\mu)=(-1)^{l_{z}}P^{l_{1}l_{2};l\,l_{z};L}_{X_{1}X_{2}}(k,-\mu).$
(227)

Combining the above equation with the flip symmetry of Eq. (226), we derive

$P^{l_{1}l_{2};l\,l_{z};L\,*}_{X_{1}X_{2}}(k,\mu)=P^{l_{1}l_{2};l\,l_{z};L}_{X_{1}X_{2}}(k,\mu),$
(228)

i.e., the invariant spectrum are real functions. This conclusion holds
regardless of whether the dependence of $\mu$ is present or not in the
invariant spectrum.

The parity transformation of the power spectrum is given by

$P^{(l_{1}l_{2})}_{X_{1}X_{2}m_{1}m_{2}}(\bm{k};\hat{\bm{z}})\xrightarrow{\mathbb{P}}(-1)^{p_{X_{1}}+p_{X_{2}}+l_{1}+l_{2}}P^{(l_{1}l_{2})}_{X_{1}X_{2}m_{1}m_{2}}(-\bm{k};-\hat{\bm{z}}),$
(229)

because of Eqs. (98) and (210). Accordingly, the transformation of the
invariant power spectrum is given by

$P^{l_{1}l_{2};l\,l_{z};L}_{X_{1}X_{2}}(k,\mu)\xrightarrow{\mathbb{P}}(-1)^{p_{X_{1}}+p_{X_{2}}+l_{1}+l_{2}+l+l_{z}}P^{l_{1}l_{2};l\,l_{z};L}_{X_{1}X_{2}}(k,\mu).$
(230)

In the Universe with parity symmetry, the power spectrum in redshift space is
invariant under the parity transformation, we thus have

$p_{X_{1}}+p_{X_{2}}+l_{1}+l_{2}+l+l_{z}=\mathrm{even}.$ (231)

The power spectrum in redshift space has an interchange symmetry,

$P^{(l_{2}l_{1})}_{X_{2}X_{1}m_{2}m_{1}}(\bm{k};\hat{\bm{z}})=P^{(l_{1}l_{2})}_{X_{1}X_{2}m_{1}m_{2}}(-\bm{k};\hat{\bm{z}}),$
(232)

and the corresponding symmetry for the invariant spectrum is given by

$P^{l_{2}l_{1};l\,l_{z};L}_{X_{2}X_{1}}(k,\mu)=(-1)^{l_{1}+l_{2}+l+L}P^{l_{1}l_{2};l\,l_{z};L}_{X_{1}X_{2}}(k,-\mu).$
(233)

####  IV.1.5 The linear power spectrum in real space

We first consider the power spectrum with the lowest-order approximation, or
the linear power spectrum, which is the simplest application of our formalism.
The formal expression for the power spectrum of scalar fields in the
propagator formalism with Gaussian initial conditions [23, 53, 17] is
straightforwardly generalized to the case of tensor field. In the lowest
order, we simply have

$P^{(l_{1}l_{2})}_{X_{1}X_{2}m_{1}m_{2}}(\bm{k})=\Pi^{2}(\bm{k})\hat{\Gamma}^{(1)}_{X_{1}l_{1}m_{1}}(\bm{k})\hat{\Gamma}^{(1)}_{X_{2}l_{2}m_{2}}(-\bm{k})P_{\mathrm{L}}(k),$
(234)

where the parity symmetry for the vertex resummation factor,
$\Pi(-\bm{k})=\Pi(\bm{k})$ is taken into account.

The vertex resummation factor $\Pi(\bm{k})$ is given by Eq. (73). In the
lowest order of the perturbation theory, this factor can be approximated by
$\Pi(\bm{k})=1$. This factor strongly damps the power spectrum on small
scales, where the linear theory does not apply. However, sometimes it is
useful to keep this exponent even in the linear theory, especially in redshift
space where the peculiar velocities damp the power on small scales. When the
resummation factor is retained, the linear approximation for this factor is
derived from Eqs. (73), (78), (82), and we have [59, 17, 18]

$\Pi(k)=\exp\left[-\frac{k^{2}}{12\pi^{2}}\int dpP_{\mathrm{L}}(p)\right],$
(235)

in real space, and

$\Pi(k,\mu)=\exp\left\\{-\frac{k^{2}}{12\pi^{2}}\left[1+f(f+2)\mu^{2}\right]\int
dpP_{\mathrm{L}}(p)\right\\},$ (236)

in redshift space, where $\mu\equiv\hat{\bm{z}}\cdot\hat{\bm{k}}$ is the
direction cosine between the wavevector and the line of sight, as usual.
Substituting $f=0$ in the expression in redshift space reduces to the
expression in real space as a matter of course.

In real space, the first-order propagator in terms of the invariant propagator
is given by Eq. (159). Thus Eq. (234) reduces to

$P^{(l_{1}l_{2})}_{X_{1}X_{2}m_{1}m_{2}}(\bm{k})=\frac{(-1)^{l_{1}}}{\sqrt{\\{l_{1}\\}\\{l_{2}\\}}}\,\hat{\Gamma}^{(1)}_{X_{1}l_{1}}(k)\hat{\Gamma}^{(1)}_{X_{2}l_{2}}(k)\Pi^{2}(k)\,P_{\mathrm{L}}(k)\\\
\times\sum_{l}\\{l\\}\left[l_{1}\,l_{2}\,l\right]_{m_{1}m_{2}}^{\phantom{m_{1}m_{2}}m}C_{lm}(\hat{\bm{k}}).$
(237)

The invariant power spectrum is easily read off from the above expression and
Eq. (213), and we have

$P^{l_{1}l_{2};l}_{X_{1}X_{2}}(k)=\frac{(-1)^{l_{1}}\\{l\\}}{\sqrt{{\\{l_{1}\\}\\{l_{2}\\}}}}\begin{pmatrix}l_{1}&l_{2}&l\\\
0&0&0\end{pmatrix}\\\
\times\hat{\Gamma}^{(1)}_{X_{1}l_{1}}(k)\hat{\Gamma}^{(1)}_{X_{2}l_{2}}(k)\Pi^{2}(k)P_{\mathrm{L}}(k).$
(238)

Further substituting the lowest-order approximation of the first-order
propagator, Eq. (195), the above equation is explicitly given by

$P^{l_{1}l_{2};l}_{X_{1}X_{2}}(k)=\Pi^{2}(k)P_{\mathrm{L}}(k)\\\
\times\Biggl{\\{}\frac{(-1)^{l_{1}}\,\\{l\\}}{\sqrt{{\\{l_{1}\\}\\{l_{2}\\}}}}\begin{pmatrix}l_{1}&l_{2}&l\\\
0&0&0\end{pmatrix}c^{(1)}_{X_{1}l_{1}}(k)c^{(1)}_{X_{2}l_{2}}(k)\\\
+\left[\delta_{l_{1}l}\delta_{l_{2}0}c^{(1)}_{X_{1}l_{1}}(k)c^{(0)}_{X_{2}}+(-1)^{l}\delta_{l_{1}0}\delta_{l_{2}l}c^{(0)}_{X_{1}}c^{(1)}_{X_{2}l_{2}}(k)\right]\\\
+\delta_{l0}\delta_{l_{1}0}\delta_{l_{2}0}c^{(0)}_{X_{1}}c^{(0)}_{X_{2}}\Biggr{\\}}.$
(239)

This is one of the generic predictions of our theory.

As a consistency check, we consider the scalar case, $l_{1}=l_{2}=0$ of auto
power spectrum, $X_{1}=X_{2}=X$. In this case, Eq. (239) simply reduces to

$P^{00;0}_{XX}(k)=\left[b_{X}(k)\right]^{2}\Pi^{2}(k)P_{\mathrm{L}}(k),$ (240)

where

$b_{X}(k)\equiv c^{(1)}_{X0}(k)+c^{(0)}_{X}$ (241)

corresponds to the linear bias factor. Eq. (240) is consistent with the result
of the usual linear power spectrum of biased density field.

The Kronecker’s symbols in the Eq. (239) indicate that the resulting
expressions are quite different between scalar fields ($l_{1}=0$ or $l_{2}=0$)
and non-scalar fields ($l_{1}\geq 1$ or $l_{2}\geq 1$). When the fields
$X_{1}$ and $X_{2}$ are both non-scalar fields with $l_{1}\neq 0$ and
$l_{2}\neq 0$, Eq. (239) simply reduces to

$P^{l_{1}l_{2};l}_{X_{1}X_{2}}(k)=\frac{(-1)^{l_{1}}\,\\{l\\}}{\sqrt{{\\{l_{1}\\}\\{l_{2}\\}}}}\begin{pmatrix}l_{1}&l_{2}&l\\\
0&0&0\end{pmatrix}c^{(1)}_{X_{1}l_{1}}(k)c^{(1)}_{X_{2}l_{2}}(k)\\\
\times\Pi^{2}(k)P_{\mathrm{L}}(k),\quad(l_{1}\neq 0\mathrm{\ \ and\ \
}l_{2}\neq 0).$ (242)

When $X_{1}$ is a non-scalar field and $X_{2}$ is a scalar field with
$l_{1}\neq 0$ and $l_{2}=0$, Eq. (239) reduces to

$P^{l_{1}0;l}_{X_{1}X_{2}}(k)=\delta_{l_{1}l}c^{(1)}_{X_{1}l_{1}}(k)b_{X_{2}}(k)\Pi^{2}(k)P_{\mathrm{L}}(k),\quad(l_{1}\neq
0),$ (243)

where $b_{X}(k)$ is given by Eq. (241). The above equation survives only when
$l=l_{1}$. When $X_{1}$ and $X_{2}$ are both scalar fields with
$l_{1}=l_{2}=0$, Eq. (239) reduces to

$P^{00;l}_{X_{1}X_{2}}(k)=\delta_{l0}b_{X_{1}}(k)b_{X_{2}}(k)\Pi^{2}(k)P_{\mathrm{L}}(k),$
(244)

which survives only when $l=0$. The last equation coincides with Eq. (240)
when $X_{1}=X_{2}=X$ and $l=0$.

####  IV.1.6 The linear power spectrum in redshift space

In redshift space, the lowest-order propagator is given by Eq. (177), which
can be substituted into Eq. (234). As a result, there appears a product of
bipolar spherical harmonics, which can be represented by a superposition of a
single bipolar spherical harmonics using a formula with $9j$-symbols [46], as
given in Eq. (390) of Appendix B in our notation. The definition of the
$9j$-symbol is given by Eq. (421), and useful properties are summarized in
Eqs. (422)–(425). After some calculation along the line described above, we
have

$P^{(l_{1}l_{2})}_{X_{1}X_{2}m_{1}m_{2}}(\bm{k};\hat{\bm{z}})=\Pi^{2}(\bm{k})P_{\mathrm{L}}(k)\sum_{l_{z1},l_{z2},l_{1}^{\prime},l_{2}^{\prime}}(-1)^{l_{2}^{\prime}}\hat{\Gamma}^{(1)l_{1}l_{z1}}_{X_{1}l_{1}^{\prime}}(k,\mu)\\\
\times\hat{\Gamma}^{(1)l_{2}l_{z2}}_{X_{2}l_{2}^{\prime}}(k,-\mu)\sum_{l,l_{z},L}(-1)^{l+l_{z}+L}\\{l\\}\\{l_{z}\\}\\{L\\}\left(l_{1}\,l_{2}\,L\right)_{m_{1}m_{2}}^{\phantom{m_{1}m_{2}}M}\\\
\times\begin{pmatrix}l_{z1}&l_{z2}&l_{z}\\\
0&0&0\end{pmatrix}\begin{pmatrix}l_{1}^{\prime}&l_{2}^{\prime}&l\\\
0&0&0\end{pmatrix}\begin{Bmatrix}l_{z1}&l_{z2}&l_{z}\\\
l_{1}^{\prime}&l_{2}^{\prime}&l\\\
l_{1}&l_{2}&L\end{Bmatrix}X^{l_{z}l}_{LM}(\hat{\bm{z}},\hat{\bm{k}}).$ (245)

Substitution of Eq. (203) into the above gives an explicit result. The
invariant power spectrum is easily read off from the expression of Eq. (223),
including the case that the dependence of the directional cosine $\mu$ in
invariant propagators is present. We derive

$P^{l_{1}l_{2};l\,l_{z};L}_{X_{1}X_{2}}(k,\mu)=\Pi^{2}(k,\mu)P_{\mathrm{L}}(k)\,(-1)^{l+l_{z}}\\{l\\}\\{l_{z}\\}\sqrt{\\{L\\}}\\\
\times\sum_{l_{z1},l_{z2},l_{1}^{\prime},l_{2}^{\prime}}\begin{pmatrix}l_{z1}&l_{z2}&l_{z}\\\
0&0&0\end{pmatrix}\begin{pmatrix}l_{1}^{\prime}&l_{2}^{\prime}&l\\\
0&0&0\end{pmatrix}\begin{Bmatrix}l_{z1}&l_{z2}&l_{z}\\\
l_{1}^{\prime}&l_{2}^{\prime}&l\\\ l_{1}&l_{2}&L\end{Bmatrix}\\\
\times(-1)^{l_{2}^{\prime}}\hat{\Gamma}^{(1)l_{1}l_{z1}}_{X_{1}l_{1}^{\prime}}(k,\mu)\hat{\Gamma}^{(1)l_{2}l_{z2}}_{X_{2}l_{2}^{\prime}}(k,-\mu).$
(246)

This is one of the generic predictions of our theory.

As a consistency check, one can recover the result of the power spectrum in
real space, Eqs. (237) and (238). In fact, substituting $l_{z1}=l_{z2}=0$ into
Eqs. (245) and (246), using identities
$\hat{\Gamma}^{(1)l0}_{Xl^{\prime}}(k,\mu)=\delta_{l^{\prime}l}\hat{\Gamma}^{(1)}_{Xl}(k)$
when $f=0$, $(00l)_{000}=\delta_{l0}$, and a special case of $9j$-symbols
[46],

$\begin{Bmatrix}0&0&0\\\ l_{1}^{\prime}&l_{2}^{\prime}&l_{3}^{\prime}\\\