l_{1}&l_{2}&l_{3}\end{Bmatrix}=\frac{\delta_{l_{1}l_{1}^{\prime}}\delta_{l_{2}l_{2}^{\prime}}\delta_{l_{3}l_{3}^{\prime}}}{\sqrt{\\{l_{1}\\}\\{l_{2}\\}\\{l_{3}\\}}},$
(247)

the corresponding equations in real space are recovered.

As another consistency check, one can also see if Eq. (245) reproduces the
well-known result for the scalar case, i.e., the Kaiser’s formula [62]. In the
case that the directional dependence is completely expanded, substituting Eq.
(203) into Eq. (245), putting $l_{1}=l_{2}=m_{1}=m_{2}=0$, $X_{1}=X_{2}=X$,
and using Eqs. (410) and (247), we derive

$P^{(00)}_{XX00}(\bm{k};\hat{\bm{z}})=\Pi^{2}(k,\mu)\left[b_{X}(k)\right]^{2}P_{\mathrm{L}}(k)\\\
\times\left[\left(1+\frac{2}{3}\beta+\frac{1}{5}\beta^{2}\right)\mathit{P}_{0}(\mu)+\left(\frac{4}{3}\beta+\frac{4}{7}\beta^{2}\right)\mathit{P}_{2}(\mu)\right.\\\
\left.+\frac{8}{35}\beta^{2}\mathit{P}_{4}(\mu)\right],$ (248)

where $\beta\equiv c^{(0)}_{X}f/b_{X}(k)$ corresponds to the redshift-space
distortion parameter. When the scalar field is given by density field of
biased objects, this expression is equivalent to a known result of the linear
power spectrum of a biased (scalar) field in redshift space [63].

In the lowest-order approximation of the first-order propagator in redshift
space, Eq. (203), either $l$ or $l_{z}$ is zero in each term. Therefore, at
least two of the indices $l_{1}$, $l_{2}$, $l_{z1}$, $l_{z2}$ are zero in the
$9j$-symbol of Eq. (246). When two of the indices in the $9j$-symbol are zero,
there are formulas [46],

$\displaystyle\begin{Bmatrix}l_{1}&l_{2}&l_{3}\\\ l_{4}&l_{5}&l_{6}\\\
l_{7}&0&0\end{Bmatrix}$
$\displaystyle=\delta_{l_{7}0}\frac{\delta_{l_{1}l_{4}}\delta_{l_{2}l_{5}}\delta_{l_{3}l_{6}}}{\sqrt{\\{l_{1}\\}\\{l_{2}\\}\\{l_{3}\\}}},$
(249) $\displaystyle\begin{Bmatrix}l_{1}&l_{2}&l_{3}\\\ l_{4}&0&l_{6}\\\
l_{7}&l_{8}&0\end{Bmatrix}$
$\displaystyle=(-1)^{l_{1}-l_{2}-l_{3}}\frac{\delta_{l_{4}l_{6}}\delta_{l_{2}l_{8}}\delta_{l_{3}l_{6}}\delta_{l_{7}l_{8}}}{\sqrt{\\{l_{2}\\}\\{l_{3}\\}}},$
(250)

A $9j$-symbol, in which two of the indices are zero, can always be represented
in a form of the lhs of the above equations by means of the symmetric
properties of the $9j$-symbol, Eqs. (422)–(424). Thus, with the lowest-order
approximation of the first-order propagator given by Eq. (203), the summation
of Eq. (246) is explicitly expanded to give

$\displaystyle
P^{l_{1}l_{2};l\,l_{z};L}_{X_{1}X_{2}}(k,\mu)=\Pi^{2}(k,\mu)P_{\mathrm{L}}(k)$
$\displaystyle\quad\times\Biggl{\\{}\delta_{l_{z}0}\delta_{lL}\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)$
$\displaystyle\hskip
36.0pt+\left(1+\frac{f}{3}\right)\delta_{l_{z}0}\delta_{lL}\Biggl{[}\delta_{l_{1}l}\delta_{l_{2}0}c^{(1)}_{X_{1}l_{1}}(k)c^{(0)}_{X_{2}}$
$\displaystyle\hskip
132.0pt+(-1)^{l}\delta_{l_{1}0}\delta_{l_{2}l}c^{(0)}_{X_{1}}c^{(1)}_{X_{2}l_{2}}(k)\Biggr{]}$
$\displaystyle\hskip
36.0pt+\frac{2f}{3}\delta_{l_{z}2}\\{l\\}\Biggl{[}(-1)^{l}\delta_{l_{1}L}\delta_{l_{2}0}\frac{1}{\sqrt{\\{l_{1}\\}}}\begin{pmatrix}2&l_{1}&l\\\
0&0&0\end{pmatrix}c^{(1)}_{X_{1}l_{1}}(k)c^{(0)}_{X_{2}}$ $\displaystyle\hskip
90.0pt+\delta_{l_{1}0}\delta_{l_{2}L}\frac{1}{\sqrt{\\{l_{2}\\}}}\begin{pmatrix}2&l_{2}&l\\\
0&0&0\end{pmatrix}c^{(0)}_{X_{1}}c^{(1)}_{X_{2}l_{2}}(k)\Biggl{]}$
$\displaystyle\hskip
36.0pt+\delta_{L0}\delta_{l_{1}0}\delta_{l_{2}0}\delta_{l\,l_{z}}\Biggl{[}\delta_{l_{z}0}\left(1+\frac{2f}{3}+\frac{f^{2}}{5}\right)$
$\displaystyle\hskip
66.0pt+\sqrt{5}\delta_{l_{z}2}\left(\frac{4f}{3}+\frac{4f^{2}}{7}\right)+3\,\delta_{l_{z}4}\frac{8f^{2}}{35}\Biggr{]}c^{(0)}_{X_{1}}c^{(0)}_{X_{2}}\Biggr{\\}},$
(251)

where an identity

$\begin{pmatrix}2&2&l\\\
0&0&0\end{pmatrix}^{2}=\frac{1}{5}\delta_{l0}+\frac{2}{35}\delta_{l2}+\frac{2}{35}\delta_{l4}$
(252)

is used in deriving the above result. The Eq. (251) is one of the generic
predictions of our theory.

The above Eq. (251) is nonzero only when $l_{z}=0,2,4$, and the redshift-space
distortions in linear theory contains only monopole ($l_{z}=0$), quadrupole
($l_{z}=2$) and hexadecapole ($l_{z}=4$) components, as is well known in the
scalar perturbation theory, Eq. (248), while higher-order corrections contain
higher-order multipoles for the redshift-space distortions.

In linear theory besides the resummation factor $\Pi(k,\mu)$, higher-rank
tensors are not affected by the redshift-space distortions as the first-order
propagator of Eq. (203) does not neither. Thus, if neither $l_{1}$ nor $l_{2}$
is zero, Eq. (251) is simply given by

$P^{l_{1}l_{2};l\,l_{z};L}_{X_{1}X_{2}}(k,\mu)=\delta_{l_{z}0}\delta_{lL}\Pi^{2}(k,\mu)P_{\mathrm{L}}(k)\\\
\times\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),\\\
(l_{1}\neq 0\mathrm{\ \ and\ \ }l_{2}\neq 0),$ (253)

which survives only when $l_{z}=0$ and $l=L$. The above equation without the
resummation factor exactly coincides with the result in real space, Eq. (242).
Namely, redshift-space distortions on higher-rank tensors are nonlinear
effects. In fact, nonlinear loop corrections in the power spectrum do
introduce the redshift-space distortion effects, as explicitly shown in Paper
II [43].

As in the case of real space, Kronecker’s symbols in the Eq. (251) indicate
that the resulting expression is 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$, the expression reduces to Eq. (253) above. 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. (251) reduces to

$P^{l_{1}0;l\,l_{z};L}_{X_{1}X_{2}}(k,\mu)=\delta_{Ll_{1}}\Pi^{2}(k,\mu)P_{\mathrm{L}}(k)b_{X_{2}}(k)c^{(1)}_{X_{1}l_{1}}(k)\\\
\times\left[\delta_{l_{z}0}\delta_{ll_{1}}\left(1+\frac{\beta_{2}}{3}\right)\right.\\\
\left.+(-1)^{l}\delta_{l_{z}2}\frac{2\beta_{2}}{3}\frac{\\{l\\}}{\sqrt{\\{l_{1}\\}}}\begin{pmatrix}2&l_{1}&l\\\
0&0&0\end{pmatrix}\right],$ (254)

where $\beta_{2}\equiv c^{(0)}_{X_{2}}f/b_{X_{2}}(k)$ and $b_{X}(k)$ is the
linear bias factor of Eq. (241). When $X_{1}$ and $X_{2}$ are both scalar
fields with $l_{1}=l_{2}=0$, Eq. (251) reduces to

$P^{00;l\,l_{z};L}_{X_{1}X_{2}}(k,\mu)=\delta_{L0}\delta_{ll_{z}}\Pi^{2}(k,\mu)P_{\mathrm{L}}(k)b_{X_{1}}(k)b_{X_{2}}(k)\\\
\times\Biggl{\\{}\delta_{l_{z}0}\left[1+\frac{1}{3}\left(\beta_{1}+\beta_{2}\right)+\frac{1}{5}\beta_{1}\beta_{2}\right]\\\
+\sqrt{5}\delta_{l_{z}2}\left[\frac{2}{3}\left(\beta_{1}+\beta_{2}\right)+\frac{4}{7}\beta_{1}\beta_{2}\right]+3\,\delta_{l_{z}4}\frac{8}{35}\beta_{1}\beta_{2}\Biggr{\\}},$
(255)

where $\beta_{1}\equiv c^{(0)}_{X_{1}}f/b_{X_{1}}(k)$. The last form is
consistent with Eq. (248).

####  IV.1.7 Alternative representation of the linear power spectrum in
redshift space

In the above examples, we consider the case that the dependence on the
directional cosine $\mu$ of the first-order normalized propagator in redshift
space is completely expanded into spherical harmonics, and the propagator is
given by Eq. (203). We now consider the other case that the dependence is not
completely expanded and the propagator is given by Eq. (206). In this case, we
substitute Eq. (206) into Eq. (246). Because of the simple form of Eq. (206),
the summations over $l_{z1}$, $l_{z2}$, $l_{1}^{\prime}$, $l_{2}^{\prime}$ are
trivially evaluated, resulting in,

$P^{l_{1}l_{2};l\,l_{z};L}_{X_{1}X_{2}}(k,\mu)=\delta_{l_{z}0}\delta_{Ll}\Pi^{2}(k,\mu)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(1+f\mu^{2}\right)\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]\\\
+\left(1+f\mu^{2}\right)^{2}\delta_{l0}\delta_{l_{1}0}\delta_{l_{2}0}c^{(0)}_{X_{1}}c^{(0)}_{X_{2}}\Biggr{\\}}.$
(256)

As a consistency check, one immediately see that this equation reduces to Eq.
(239) when $f=0$.

When the fields $X_{1}$ and $X_{2}$ are both non-scalar fields with $l_{1}\neq
0$ and $l_{2}\neq 0$, Eq. (256) simply reduces to

$P^{l_{1}l_{2};l\,l_{z};L}_{X_{1}X_{2}}(k,\mu)=\delta_{l_{z}0}\delta_{Ll}\Pi^{2}(k,\mu)P_{\mathrm{L}}(k)\\\
\times\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),\\\
(l_{1}\neq 0\mathrm{\ \ and\ \ }l_{2}\neq 0).$ (257)

which survives only when $l_{z}=0$ and $L=l$. This equation is exactly the
same as Eq. (253). 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. (256) reduces to

$P^{l_{1}0;l\,l_{z};L}_{X_{1}X_{2}}(k)=\delta_{l_{z}0}\delta_{lL}\delta_{l_{1}l}\left(1+f\mu^{2}\right)\\\
\times
c^{(1)}_{X_{1}l_{1}}(k)b_{X_{2}}(k)\Pi^{2}(k)P_{\mathrm{L}}(k),\quad(l_{1}\neq
0),$ (258)

which survives only when $l_{z}=0$ and $L=l=l_{1}$. When $X_{1}$ and $X_{2}$
are both scalar fields with $l_{1}=l_{2}=0$, Eq. (256) reduces to

$P^{00;l\,l_{z};L}_{X_{1}X_{2}}(k)=\delta_{l0}\delta_{l_{z}0}\delta_{L0}\left(1+f\mu^{2}\right)^{2}b_{X_{1}}(k)b_{X_{2}}(k)\Pi^{2}(k)P_{\mathrm{L}}(k),$
(259)

which survives only when $l=l_{z}=L=0$. The last equation corresponds to the
Kaiser’s formula for the density power spectrum in redshift space [62].

###  IV.2 The correlation function

While the power spectrum is defined in Fourier space, the counterpart in
configuration space is the correlation function,
$\xi^{(l_{1}l_{2})}_{X_{1}X_{2}m_{1}m_{2}}(\bm{x})$, which is defined by

$\left\langle
F_{X_{1}l_{1}m_{1}}(\bm{x}_{1})F_{X_{2}l_{2}m_{2}}(\bm{x}_{2})\right\rangle_{\mathrm{c}}=\xi^{(l_{1}l_{2})}_{X_{1}X_{2}m_{1}m_{2}}(\bm{x}_{1}-\bm{x}_{2}),$
(260)

where the tensor field $F_{Xlm}(\bm{x})$ on the lhs now corresponds to the one
in configuration space. On the rhs, the correlation function is a function of
the relative vector between the two positions, $\bm{x}_{1}-\bm{x}_{2}$,
because of the statistical homogeneity of the Universe. It is standard that
the correlation function and the power spectrum are related by a three-
dimensional Fourier transform as

$\xi^{(l_{1}l_{2})}_{X_{1}X_{2}m_{1}m_{2}}(\bm{x})=\int\frac{d^{3}k}{(2\pi)^{3}}e^{i\bm{k}\cdot\bm{x}}P^{(l_{1}l_{2})}_{X_{1}X_{2}m_{1}m_{2}}(\bm{k}).$
(261)

First, we consider the correlation function in real space. In exactly the same
way of deriving Eq. (213), the rotational symmetry requires that the
correlation function should have a form,

$\xi^{(l_{1}l_{2})}_{X_{1}X_{2}m_{1}m_{2}}(\bm{x})=\sum_{l}i^{l}\left(l_{1}\,l_{2}\,l\right)_{m_{1}m_{2}}^{\phantom{m_{1}m_{2}}m}C_{lm}(\hat{\bm{x}})\,\xi^{l_{1}l_{2};l}_{X_{1}X_{2}}(x),$
(262)

where the last factor $\xi^{l_{1}l_{2};l}_{X_{1}X_{2}}(x)$ corresponds to the
invariant correlation function. Comparing with the definition of the invariant
power spectrum of Eq. (213), we should note that the phase factor $i^{l}$ is
additionally present. One can show that the invariant correlation function of
$\xi^{l_{1}l_{2};l}_{X_{1}X_{2}}(x)$ defined above with the phase factor is a
real function.

Substituting Eq. (213) into Eq. (261), and using a formula of Eq. (381) in
Appendix B, we have

$\xi^{l_{1}l_{2};l}_{X_{1}X_{2}}(x)=\int\frac{k^{2}dk}{2\pi^{2}}j_{l}(kx)P^{l_{1}l_{2};l}_{X_{1}X_{2}}(k),$
(263)

i.e., the invariant correlation function in real space is given by a Hankel
transform of the invariant power spectrum. Using the completeness relation of
the spherical Bessel functions, the inverse relation of the above transform is
given by

$P^{l_{1}l_{2};l}_{X_{1}X_{2}}(k)=4\pi\int
r^{2}drj_{l}(kx)\xi^{l_{1}l_{2};l}_{X_{1}X_{2}}(x).$ (264)

The correlation function in redshift space is also represented by the
invariant spectrum in redshift space. The relation between the correlation
function and the power spectrum is just given by Eq. (261) as well in redshift
space, provided that both explicitly depend on the direction of the line of
sight, $\hat{\bm{z}}$,

$\xi^{(l_{1}l_{2})}_{X_{1}X_{2}m_{1}m_{2}}(\bm{x};\hat{\bm{z}})=\int\frac{d^{3}k}{(2\pi)^{3}}e^{i\bm{k}\cdot\bm{x}}P^{(l_{1}l_{2})}_{X_{1}X_{2}m_{1}m_{2}}(\bm{k};\hat{\bm{z}}).$
(265)

Due to the rotational symmetry, we have

$\xi^{(l_{1}l_{2})}_{X_{1}X_{2}m_{1}m_{2}}(\bm{x};\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}}i^{l}X^{l_{z}l}_{LM}(\hat{\bm{z}},\hat{\bm{x}})\xi^{l_{1}l_{2};l\,l_{z};L}_{X_{1}X_{2}}(x),$
(266)

just as in the case of the power spectrum of Eq. (223), and the last factor
$\xi^{l_{1}l_{2};l\,l_{z};L}_{X_{1}X_{2}}(x)$ corresponds to the invariant
correlation function. When the invariant power spectrum in redshift space
$P^{l_{1}l_{2};l\,l_{z};L}_{X_{1}X_{2}}(k)$ does not contain any dependence on
the directional cosine of the wavevector, $\mu=\cos\theta_{k}$, the invariant
correlation function $\xi^{l_{1}l_{2};l\,l_{z};L}_{X_{1}X_{2}}(x)$ is
straightforwardly related to the invariant spectrum via the Hankel transform
just in the similar way in the case of real space above,

$\xi^{l_{1}l_{2};l\,l_{z};L}_{X_{1}X_{2}}(x)=\int\frac{k^{2}dk}{2\pi^{2}}j_{l}(kx)P^{l_{1}l_{2};l\,l_{z};L}_{X_{1}X_{2}}(k),$
(267)

and

$P^{l_{1}l_{2};l\,l_{z};L}_{X_{1}X_{2}}(k)=4\pi\int
r^{2}drj_{l}(kx)\xi^{l_{1}l_{2};l\,l_{z};L}_{X_{1}X_{2}}(x).$ (268)

However, when the invariant spectrum in redshift space
$P^{l_{1}l_{2};l\,l_{z};L}_{X_{1}X_{2}}(k,\mu)$ does depend on the directional
cosine, the above equations no longer hold. In this case, the invariant
correlation function in can still be represented by the invariant power
spectrum. In order to derive the relation, we first note that Eq. (266) can be
inverted by applying orthogonality relations of the $3j$-symbols, Eq. (407),
and of the bipolar spherical harmonics, Eq. (388), and we have

$\xi^{l_{1}l_{2};l\,l_{z};L}_{X_{1}X_{2}}(x)=i^{l}(-1)^{l_{1}+l_{2}+l_{z}+L}\\{l\\}\\{l_{z}\\}\sqrt{\\{L\\}}\\\
\times\int\frac{d^{2}\hat{z}}{4\pi}\frac{d^{2}\hat{x}}{4\pi}X^{l_{z}l}_{LM}(\hat{\bm{z}},\hat{\bm{x}})\xi^{(l_{1}l_{2})}_{X_{1}X_{2}m_{1}m_{2}}(\bm{x};\hat{\bm{z}}).$
(269)

We substitute the expression of the power spectrum, Eq. (223), into Eqs.
(265), and then the result is substituted into (269). The resulting equation
is calculated by using orthogonality relations of the $3j$-symbols, Eq. (407),
the Rayleigh expansion formula, Eq. (379), sum rules of spherical harmonics
and $6j$-symbols, Eqs. (378) and (426), and the formulas of Eqs. (382) and
(383). After these manipulations, the expression of the invariant correlation
function in redshift space finally reduces to

$\xi^{l_{1}l_{2};l\,l_{z};L}_{X_{1}X_{2}}(x)=(-1)^{L}\\{l\\}\\{l_{z}\\}\sum_{l^{\prime},l_{z}^{\prime},L^{\prime}}(-1)^{L^{\prime}}\\{L^{\prime}\\}\\\
\times\begin{pmatrix}l_{z}&l_{z}^{\prime}&L^{\prime}\\\
0&0&0\end{pmatrix}\begin{pmatrix}l&l^{\prime}&L^{\prime}\\\
0&0&0\end{pmatrix}\begin{Bmatrix}l_{z}&l_{z}^{\prime}&L^{\prime}\\\
l^{\prime}&l&L\end{Bmatrix}\\\
\times\int\frac{d\mu}{2}P_{L^{\prime}}(\mu)\int\frac{k^{2}dk}{2\pi^{2}}j_{l}(kx)P^{l_{1}l_{2};l^{\prime}l_{z}^{\prime};L}_{X_{1}X_{2}}(k,\mu),$
(270)

where $P_{L^{\prime}}(\mu)$ is the Legendre polynomial of order $L^{\prime}$.
When the invariant spectrum does not depend on $\mu$, it is straightforwardly
shown that the above equation reduces to the previous result of Eq. (267),
noting a special case of the $6j$-symbol [46],

$\begin{Bmatrix}l_{z}&l_{z}^{\prime}&0\\\
l^{\prime}&l&L\end{Bmatrix}=(-1)^{l+l_{z}+L}\frac{\delta_{l_{z}^{\prime}l_{z}}\delta_{l^{\prime}l}}{\sqrt{\\{l_{z}\\}\\{l\\}}}.$
(271)

###  IV.3  Scale-dependent bias in the power spectrum with non-Gaussian
initial conditions

In the presence of bias, the local-type non-Gaussianity is known to produce
the scale-dependent bias in the power spectrum on very large scales [3]. The
same effect appears even in the shape correlations [64, 65], which is
considered to be a second-rank ($l=2$) tensor field in our context. Higher
moments of the shape correlations are also investigated [42].

In our formalism, these previous findings are elegantly reproduced as shown
below. The derivation of scale-dependent bias from the primordial non-
Gaussianity in the iPT formalism of scalar fields is already given in Ref.
[19]. We can simply generalize the last method to the case of higher-rank
tensor fields. We consider the clustering in real space below, while the
generalization to that in redshift space is straightforward.

The primordial bispectrum $B_{\mathrm{L}}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$
is defined by three-point correlations of linear density contrast in Fourier
space,

$\left\langle\delta_{\mathrm{L}}(\bm{k}_{1})\delta_{\mathrm{L}}(\bm{k}_{2})\delta_{\mathrm{L}}(\bm{k}_{3})\right\rangle_{\mathrm{c}}=(2\pi)^{3}\delta_{\mathrm{D}}^{3}(\bm{k}_{1}+\bm{k}_{2}+\bm{k}_{3})B_{\mathrm{L}}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3}),$
(272)

where $\langle\cdots\rangle_{\mathrm{c}}$ denotes the three-point cumulant, or
the connected part. Applying a straightforward generalization of the
corresponding result of iPT [19] to tensor fields, the lowest-order
contribution of the primordial bispectrum to the power spectrum in the
formalism of iPT is given by

$P^{\mathrm{NG}(l_{1}l_{2})}_{X_{1}X_{2}m_{1}m_{2}}(\bm{k})=\frac{1}{2}\Gamma^{(1)}_{X_{1}l_{1}m_{1}}(-\bm{k})\\\
\times\int\frac{d^{3}p}{(2\pi)^{3}}\Gamma^{(2)}_{X_{2}l_{2}m_{2}}(\bm{p},\bm{k}-\bm{p})B_{\mathrm{L}}(\bm{k},-\bm{p},\bm{p}-\bm{k})\\\
+[(X_{1}l_{1}m_{1})\leftrightarrow(X_{2}l_{2}m_{2})],$ (273)

where we implicitly assume the parity symmetry of the Universe. The linear
power spectrum with the Gaussian initial condition, Eq. (234), is added to the
above in order to obtain the total power spectrum up to the linear order of
the initial power spectrum and bispectrum. However, the Gaussian contribution
does not generate scale-dependent bias in the limit of large scales,
$k\rightarrow 0$, and thus the non-Gaussian contribution of Eq. (273)
dominates in that limit.

The prefactor of the integral in Eq. (273) is given by Eq. (194):

$\Gamma^{(1)}_{X_{1}l_{1}m_{1}}(\bm{k})=\frac{(-1)^{l_{1}}}{\sqrt{2l_{1}+1}}\left[c^{(1)}_{X_{1}l_{1}}(k)+\delta_{l_{1}0}c^{(0)}_{X_{1}}\right]C_{l_{1}m_{1}}(\hat{\bm{k}}),$
(274)

where we consider the lowest-order approximation and the resummation factor
$\Pi(k)$ is just replaced by unity. In order to evaluate the integral in Eq.
(273) with the spherical basis, one can insert a unity $\int
d^{3}p^{\prime}\delta_{\mathrm{D}}^{3}(\bm{p}+\bm{p}^{\prime}-\bm{k})=1$ in
the integral, and re-express the delta function by a Fourier integral, and
primordial bispectrum can be also expanded into spherical harmonics by means
of the plane-wave expansion of Eq. (379). In this way, the angular
integrations are analytically performed, leaving only the radial integrals
over $p$ and $p^{\prime}$.

We are interested in the scale-dependent bias from the primordial non-
Gaussianity, and instead of deriving a general expression according to the
above procedure, we directly consider a limiting case of $k\rightarrow 0$ in
Eq. (273). In this case, the integral is approximately given by

$\int\frac{d^{3}p}{(2\pi)^{3}}\Gamma^{(2)}_{X_{2}l_{2}m_{2}}(\bm{p},-\bm{p})B_{\mathrm{L}}(\bm{k},-\bm{p},\bm{p}).$
(275)

The second-order propagator in the integrand is given by Eq. (76) with
$\bm{k}_{12}=\bm{0}$, $\bm{k}_{1}=\bm{p}$ and $\bm{k}_{2}=-\bm{p}$, and thus
we have

$\displaystyle\Gamma^{(2)}_{X_{2}l_{2}m_{2}}(\bm{p},-\bm{p})$
$\displaystyle=c^{(2)}_{X_{2}l_{2}m_{2}}(\bm{p},-\bm{p})$
$\displaystyle=\sum_{l_{2}^{\prime},l_{2}^{\prime\prime}}(-1)^{l_{2}^{\prime\prime}}c^{(2)l_{2}}_{X_{2}l_{2}^{\prime}l_{2}^{\prime\prime}}(p,p)X^{l_{2}^{\prime}l_{2}^{\prime\prime}}_{l_{2}m_{2}}(\hat{\bm{p}},\hat{\bm{p}}),$
(276)

where Eq. (134) is substituted. According to Eqs. (132), (383), (208) and
(407), we have

$X^{l_{2}^{\prime}l_{2}^{\prime\prime}}_{l_{2}m_{2}}(\hat{\bm{p}},\hat{\bm{p}})=\begin{pmatrix}l_{2}^{\prime}&l_{2}^{\prime\prime}&l_{2}\\\
0&0&0\end{pmatrix}C_{l_{2}m_{2}}(\hat{\bm{p}}).$ (277)

We decompose the primordial bispectrum with the spherical harmonics as

$B_{\mathrm{L}}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})=\sum_{l_{1},l_{2},l_{3}}B_{\mathrm{L}}^{l_{1}l_{2}l_{3}}(k_{1},k_{2},k_{3})\left(l_{1}\,l_{2}\,l_{3}\right)^{m_{1}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}).$
(278)

One can always decompose the primordial bispectrum in the above form, which is
shown in the similar of way of deriving Eq. (141) when $l=m=0$. The appearance
of the $3j$-symbols is due to the rotational symmetry. The inverse relation of
the above is given by

$B_{\mathrm{L}}^{l_{1}l_{2}l_{3}}(k_{1},k_{2},k_{3})=\\{l_{1}\\}\\{l_{2}\\}\\{l_{3}\\}\left(l_{1}\,l_{2}\,l_{3}\right)^{m_{1}m_{2}m_{3}}\\\
\times\int\frac{d^{2}\hat{k}_{1}}{4\pi}\frac{d^{2}\hat{k}_{2}}{4\pi}\frac{d^{2}\hat{k}_{3}}{4\pi}B_{\mathrm{L}}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{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})$
(279)

The particular bispectrum in the integrand of Eq. (275) is given by

$B_{\mathrm{L}}(\bm{k},-\bm{p},\bm{p})=\sum_{l}C_{lm}(\hat{\bm{k}})C_{lm}^{*}(\hat{\bm{p}})\sum_{l^{\prime},l^{\prime\prime}}(-1)^{l^{\prime}}\\\
\times\begin{pmatrix}l^{\prime}&l^{\prime\prime}&l\\\
0&0&0\end{pmatrix}B_{\mathrm{L}}^{l\,l^{\prime}l^{\prime\prime}}(k,p,p).$
(280)

Using the above equations together for the integral of Eq. (275), the angular
integration over $\hat{\bm{p}}$ can be analytically performed by using the
orthonormality relation of the spherical harmonics, Eq. (380). Comparing the
resulting equation with Eq. (213), the corresponding invariant power spectrum
is derived, and the result is given by

$P^{\mathrm{NG}\,l_{1}l_{2};l}_{X_{1}X_{2}}(k)=\frac{(-1)^{l}}{2}\frac{\\{l\\}}{\sqrt{\\{l_{1}\\}}\,\\{l_{2}\\}}\begin{pmatrix}l_{1}&l_{2}&l\\\
0&0&0\end{pmatrix}\left[c^{(1)}_{X_{1}l_{1}}(k)+\delta_{l_{1}0}c^{(0)}_{X_{1}}\right]\\\
\times\sum_{l_{1}^{\prime},l_{1}^{\prime\prime},l_{2}^{\prime},l_{2}^{\prime\prime}}(-1)^{l_{1}^{\prime\prime}+l_{2}^{\prime\prime}}\begin{pmatrix}l_{1}^{\prime}&l_{1}^{\prime\prime}&l_{2}\\\
0&0&0\end{pmatrix}\begin{pmatrix}l_{2}^{\prime}&l_{2}^{\prime\prime}&l_{2}\\\
0&0&0\end{pmatrix}\\\
\times\int\frac{p^{2}dp}{2\pi^{2}}c^{(2)l_{2}}_{X_{2}l_{2}^{\prime}l_{2}^{\prime\prime}}(p,p)B_{\mathrm{L}}^{l_{2}l_{1}^{\prime}l_{1}^{\prime\prime}}(k,p,p)\\\
+[(X_{1}l_{1})\leftrightarrow(X_{2}l_{2})].$ (281)

This is one of the generic predictions of our theory.

The scale-dependent bias from the primordial non-Gaussianity is conveniently
derived from the cross power spectrum between the linear density field and
biased objects. In particular, we define the lowest-order of scale-dependent
bias factor by

$\Delta b_{Xlm}(\bm{k})\equiv\frac{P^{\mathrm{NG}(0l)}_{\delta
X\,0m}(\bm{k})}{P_{\mathrm{L}}(k)},$ (282)

where $X_{1}=\delta$ indicates that the biased field
$F_{X_{1}l_{1}m_{1}}(\bm{k})$ is replaced by the mass density field
$\delta(\bm{k})$ with $l_{1}=m_{1}=0$, and thus $c^{(0)}_{X_{1}}=1$ and
$c^{(n)}_{X_{1}l_{1}m_{1}}=0$ for $n\geq 1$. Due to the rotational symmetry,
and as explicitly seen from Eq. (213), the directional dependence on the
wavevector on the rhs of Eq. (282) should be proportional to
$C_{lm}(\hat{\bm{k}})$, thus we naturally define an invariant bias factor
$\Delta b_{Xl}(k)$ by

$\Delta b_{Xlm}(\bm{k})=\frac{(-1)^{l}}{\sqrt{\\{l\\}}}\Delta
b_{Xl}(k)C_{lm}(\hat{\bm{k}}).$ (283)

Substituting Eq. (213) into Eq. (282), the above definition is equivalent to

$\Delta b_{Xl}(k)=\frac{P^{\mathrm{NG}\,0l;l}_{\delta
X}(k)}{P_{\mathrm{L}}(k)}.$ (284)

Substituting $l_{1}=l$, $l_{2}=0$, $c^{(0)}_{X_{2}}=1$, and
$c^{(1)}_{X_{2}l_{2}}=c^{(2)l_{2}}_{X_{2}l_{2}^{\prime}l_{2}^{\prime\prime}}=0$
into Eq. (281), the invariant scale dependent bias of Eq. (284) reduces to

$\Delta
b_{Xl}(k)=\frac{1}{2P_{\mathrm{L}}(k)}\sum_{l_{1}^{\prime},l_{1}^{\prime\prime},l_{2}^{\prime},l_{2}^{\prime\prime}}(-1)^{l_{1}^{\prime\prime}+l_{2}^{\prime\prime}}\begin{pmatrix}l_{1}^{\prime}&l_{1}^{\prime\prime}&l\\\
0&0&0\end{pmatrix}\begin{pmatrix}l_{2}^{\prime}&l_{2}^{\prime\prime}&l\\\
0&0&0\end{pmatrix}\\\
\times\int\frac{p^{2}dp}{2\pi^{2}}c^{(2)l}_{Xl_{2}^{\prime}l_{2}^{\prime\prime}}(p,p)B_{\mathrm{L}}^{l\,l_{1}^{\prime}l_{1}^{\prime\prime}}(k,p,p).$
(285)

This is one of the generic predictions of our theory.

The above formula is applicable for arbitrary models of primordial non-
Gaussianity with a given bispectrum at the lowest order. We consider a
particular model of non-Gaussianity, which is a generalization of the so-
called local-type model [66, 64]:

$B_{\mathrm{L}}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})=\frac{2\mathcal{M}(k_{3})}{\mathcal{M}(k_{1})\mathcal{M}(k_{2})}\\\
\times\sum_{l}f_{\mathrm{NL}}^{(l)}\mathit{P}_{l}(\hat{\bm{k}}_{1}\cdot\hat{\bm{k}}_{2})P_{\mathrm{L}}(k_{1})P_{\mathrm{L}}(k_{2})+\mathrm{cyc.}$
(286)

The function $\mathcal{M}(k)$ in the prefactor is given by

$\mathcal{M}(k)=\frac{2}{3}\frac{D(z)}{(1+z_{*})D(z_{*})}\frac{k^{2}T(k)}{{H_{0}}^{2}\Omega_{\mathrm{m0}}},$
(287)

where $D(z)$ is the linear growth factor with an arbitrary normalization,
$z_{*}$ is an arbitrary redshift at the matter-dominated epoch, $T(k)$ is the
transfer function, $H_{0}$ is the Hubble constant, and $\Omega_{\mathrm{m0}}$
is the density parameter of total matter. The factor $(1+z_{*})D(z_{*})$ does
not depend on the choice of $z_{*}$ as long as $z_{*}$ is deep in the matter-
dominated epoch. Some authors choose a normalization of the linear growth
factor by $(1+z_{*})D(z_{*})=1$ in which case the expression of Eq. (287) has
the simplest form. Substituting Eq. (286) into Eq. (279), we have

$B_{\mathrm{L}}^{l_{1}l_{2}l_{3}}(k_{1},k_{2},k_{3})=\frac{2\mathcal{M}(k_{3})}{\mathcal{M}(k_{1})\mathcal{M}(k_{2})}P_{\mathrm{L}}(k_{1})P_{\mathrm{L}}(k_{2})\\\
\times\delta_{l_{1}l_{2}}\delta_{l_{3}0}(-1)^{l_{1}}\sqrt{\\{l_{1}\\}}\,f_{\mathrm{NL}}^{(l_{1})}+\mathrm{cyc.}$
(288)

We need a special form of the bispectrum to further evaluate Eqs. (281) and
(285). With the above model of Eq. (286), we have

$B_{\mathrm{L}}^{l\,l_{1}^{\prime}l_{1}^{\prime\prime}}(k,p,p)\simeq\frac{2P_{\mathrm{L}}(k)}{\mathcal{M}(k)}P_{\mathrm{L}}(p)\left(\delta_{l_{1}^{\prime}l}\delta_{l_{1}^{\prime\prime}0}+\delta_{l_{1}^{\prime}0}\delta_{l_{1}^{\prime\prime}l}\right)\\\
\times(-1)^{l}\sqrt{\\{l\\}}\,f_{\mathrm{NL}}^{(l)},$ (289)

where only two terms in which $\mathcal{M}(k)$ appeared in the denominator are
retained because the other term is negligible in the limit of $k\rightarrow
0$.

Substituting Eq. (289) into Eq. (281), and using Eqs. (410) and (416), we have

$P^{\mathrm{NG}\,l_{1}l_{2};l}_{X_{1}X_{2}}(k)=\left[1+(-1)^{l_{2}}\right]f_{\mathrm{NL}}^{(l_{2})}\frac{P_{\mathrm{L}}(k)}{\mathcal{M}(k)}\frac{\\{l\\}}{\sqrt{\\{l_{1}\\}}}\begin{pmatrix}l_{1}&l_{2}&l\\\
0&0&0\end{pmatrix}\\\
\times\left[c^{(1)}_{X_{1}l_{1}}(k)+\delta_{l_{1}0}c^{(0)}_{X_{1}}\right]\sum_{l^{\prime},l^{\prime\prime}}(-1)^{l^{\prime\prime}}\begin{pmatrix}l^{\prime}&l^{\prime\prime}&l_{2}\\\
0&0&0\end{pmatrix}\\\
\times\int\frac{p^{2}dp}{2\pi^{2}}c^{(2)l_{2}}_{X_{2}l^{\prime}l^{\prime\prime}}(p,p)P_{\mathrm{L}}(p)\\\
+[(X_{1}l_{1})\leftrightarrow(X_{2}l_{2})].$ (290)

Substituting $l_{1}=0$, $l_{2}=l$, $c^{(0)}_{X_{1}}=1$, and
$c^{(1)}_{X_{1}l_{1}}=c^{(2)l_{1}}_{X_{1}l_{1}^{\prime}l_{1}^{\prime\prime}}=0$
into the above equation, and using Eq. (284), we derive

$\Delta
b_{Xl}(k)=\left[1+(-1)^{l}\right]\frac{f_{\mathrm{NL}}^{(l)}}{\mathcal{M}(k)}\sqrt{\\{l\\}}\sum_{l^{\prime},l^{\prime\prime}}(-1)^{l^{\prime\prime}}\begin{pmatrix}l^{\prime}&l^{\prime\prime}&l\\\
0&0&0\end{pmatrix}\\\
\times\int\frac{p^{2}dp}{2\pi^{2}}c^{(2)l}_{Xl^{\prime}l^{\prime\prime}}(p,p)P_{\mathrm{L}}(p).$
(291)

Apparently, the rhs survives only in the case of $l=\mathrm{even}$. This
proves that the scale-dependent bias of the rank-$l$ tensor field is sensitive
to the multipole moment $l$ of the Legendre polynomials
$P_{l}(\hat{\bm{k}}_{1}\cdot\hat{\bm{k}}_{2})$ in the model of Eq. (286) for
the primordial non-Gaussianity. That is previously shown in Ref. [42] by quite
different considerations and a method.

As a consistency check, we can see if the known results of scalar
perturbations are derived from the general result above. In the scalar case of
$l=0$, Eq. (291) reduces to

$\Delta
b_{X0}(k)=\frac{2f_{\mathrm{NL}}}{\mathcal{M}(k)}\sum_{l}\frac{1}{\sqrt{\\{l\\}}}\int\frac{p^{2}dp}{2\pi^{2}}c^{(2)0}_{Xll}(p,p)P_{\mathrm{L}}(p),$
(292)

where $f_{\mathrm{NL}}\equiv f_{\mathrm{NL}}^{(0)}$ is the original parameter
of local-type non-Gaussianity. One can confirm that the last equation is
consistent with the previously known results for the scale-dependent bias in
the halo models. In fact, the second-order renormalized bias function of a
simple halo model is given by [19]

$c^{(2)}_{\mathrm{h}}(\bm{k}_{1},\bm{k}_{2})=\frac{\delta_{\mathrm{c}}b^{\mathrm{L}}_{1}}{\sigma^{2}}W(k_{1}R)W(k_{2}R),$
(293)

where

$\sigma^{2}=\int\frac{p^{2}dp}{2\pi^{2}}W^{2}(pR)P_{\mathrm{L}}(p),$ (294)

and $b^{\mathrm{L}}_{1}$ is a Lagrangian bias parameter of linear order. Eq.
(293) is a result for the simplest case with a high-peaks (or high-mass) limit
of halos in the Press-Schechter mass function (see Ref. [19] for results of
other extended halo models). The Lagrangian bias parameter
$b^{\mathrm{L}}_{1}$ is related to the Eulerian bias parameter $b_{1}$ by
$b^{\mathrm{L}}_{1}=b_{1}-1$. Comparing the above function with Eq. (131) or
Eq. (134), and noting
$c^{(2)0}_{\mathrm{h}00}(k_{1},k_{2})=c^{(2)}_{\mathrm{h}}(k_{1},k_{2})$, the
invariant coefficient in this case is given by

$c^{(2)0}_{\mathrm{h}\,l_{1}l_{2}}(k_{1},k_{2})=\delta_{l_{1}0}\delta_{l_{2}0}\frac{\delta_{\mathrm{c}}b^{\mathrm{L}}_{1}}{\sigma^{2}}W(k_{1}R)W(k_{2}R).$
(295)

Substituting the above expression into Eq. (292), we have

$\Delta
b_{\mathrm{h}0}(k)=\frac{2f_{\mathrm{NL}}\delta_{\mathrm{c}}b^{\mathrm{L}}_{1}}{\mathcal{M}(k)}.$
(296)

Due to the normalization of the spherical basis for the scalar component, Eqs.
(7) and (31), the scale-dependent bias of the halo number density is given by
$\Delta b_{\mathrm{h}}(k)=\Delta b_{\mathrm{h}00}(k)\mathsf{Y}^{(0)}=\Delta
b_{\mathrm{h}0}$, and thus we have

$\Delta
b_{\mathrm{h}}(k)=\frac{2(b_{1}-1)f_{\mathrm{NL}}\delta_{\mathrm{c}}}{\mathcal{M}(k)}.$
(297)

This exactly reproduces a well-known result of the scale-dependent bias in the
simplest halo model, which has been first derived by a quite different method
[3].

###  IV.4  The bispectrum

####  IV.4.1  The invariant bispectrum

For yet another example of applications of our formalism, we consider the
bispectrum of the tensor field. For an illustrative purpose, we only consider
the bispectrum in real space below, while generalizing the results to those in
redshift space is straightforward with similar methods employed for the power
spectrum in Sec. IV.1.

The bispectrum
$B^{(l_{1}l_{2}l_{3})}_{X_{1}X_{2}X_{3}\,m_{1}m_{2}m_{3}}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$
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})F_{X_{3}l_{3}m_{3}}(\bm{k}_{3})\right\rangle_{\mathrm{c}}\\\
=(2\pi)^{3}\delta_{\mathrm{D}}^{3}(\bm{k}_{1}+\bm{k}_{2}+\bm{k}_{3})B^{(l_{1}l_{2}l_{3})}_{X_{1}X_{2}X_{3}\,m_{1}m_{2}m_{3}}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3}),$
(298)

as a simple generalization of the definition of the power spectrum of Eq.
(210). We consider generally the cross bispectrum among three fields, $X_{1}$,
$X_{2}$ and $X_{3}$ with irreducible components $l_{1}$, $l_{2}$ and $l_{3}$,
respectively. The auto bispectrum is straightforwardly obtained by putting
$X_{1}=X_{2}=X_{3}$.

The bispectrum
$B^{(l_{1}l_{2}l_{3})}_{X_{1}X_{2}X_{3}\,m_{1}m_{2}m_{3}}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$
defined above depends on the coordinates system. Generalizing the method of
constructing the invariant power spectrum, one can also construct the
invariant bispectrum as explicitly shown below. Due to the presence of the
delta function on the lhs of Eq. (298), the third argument $\bm{k}_{3}$ in the
bispectrum can be replaced by $-\bm{k}_{1}-\bm{k}_{2}$, and regarded as a
function of only $\bm{k}_{1}$ and $\bm{k}_{2}$. However, the resulting
expression loses apparent symmetry of permutating the subscripts $(1,2,3)$ of
$\bm{k}_{i}$, $l_{i}$ and $m_{i}$. As is usually the case in the analytic
predictions in perturbation theory, we can decompose the symmetric bispectrum
into a sum of asymmetric components as

$B^{(l_{1}l_{2}l_{3})}_{X_{1}X_{2}X_{3}\,m_{1}m_{2}m_{3}}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})=\tilde{B}^{(l_{1}l_{2}l_{3})}_{X_{1}X_{2}X_{3}\,m_{1}m_{2}m_{3}}(\bm{k}_{1},\bm{k}_{2})+\mathrm{cyc.}$
(299)

Although the decomposition into asymmetric components is not unique, it is
always possible to specify the originally symmetric bispectrum by giving the
asymmetric bispectra, and the predictions of perturbation theory are usually
given in the above form.

The directional dependence of the asymmetric bispectrum on the wavevectors can
be expanded by spherical harmonics as

$\tilde{B}^{(l_{1}l_{2}l_{3})}_{X_{1}X_{2}X_{3}m_{1}m_{2}m_{3}}(\bm{k}_{1},\bm{k}_{2})=\sum_{l_{1}^{\prime},l_{2}^{\prime}}\tilde{B}^{l_{1}l_{2}l_{3};l_{1}^{\prime}l_{2}^{\prime}}_{X_{1}X_{2}X_{3}m_{1}m_{2}m_{3};m_{1}^{\prime}m_{2}^{\prime}}(k_{1},k_{2})\\\
\times
C^{*}_{l_{1}^{\prime}m_{1}^{\prime}}(\hat{\bm{k}}_{1})C^{*}_{l_{2}^{\prime}m_{2}^{\prime}}(\hat{\bm{k}}_{2}).$
(300)

Due to the rotational transformations of Eqs. (92) and (22), and rotational
invariance of the delta function in the definition of the bispectrum, Eq.
(298), the rotational transformation of the expansion coefficients is given by

$\tilde{B}^{l_{1}l_{2}l_{3};l_{1}^{\prime}l_{2}^{\prime}}_{X_{1}X_{2}X_{3}m_{1}m_{2}m_{3};m_{1}^{\prime}m_{2}^{\prime}}(k_{1},k_{2})\\\
\xrightarrow{\mathbb{R}}\tilde{B}^{l_{1}l_{2}l_{3};l_{1}^{\prime}l_{2}^{\prime}}_{X_{1}X_{2}X_{3}m_{1}^{\prime\prime}m_{2}^{\prime\prime}m_{3}^{\prime\prime};m_{1}^{\prime\prime\prime}m_{2}^{\prime\prime\prime}}(k_{1},k_{2})D^{m_{1}^{\prime\prime}}_{(l_{1})m_{1}}(R)\\\
\times
D^{m_{2}^{\prime\prime}}_{(l_{2})m_{2}}(R)D^{m_{3}^{\prime\prime}}_{(l_{2})m_{3}}(R)D^{m_{1}^{\prime\prime\prime}}_{(l)m_{1}^{\prime}}(R)D^{m_{2}^{\prime\prime\prime}}_{(l)m_{2}^{\prime}}(R).$
(301)

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

$\tilde{B}^{(l_{1}l_{2}l_{3})}_{X_{1}X_{2}X_{3}m_{1}m_{2}m_{3}}(\bm{k}_{1},\bm{k}_{2})=\sum_{\begin{subarray}{c}L,L^{\prime}\\\
l_{1}^{\prime},l_{2}^{\prime}\end{subarray}}(-1)^{L+L^{\prime}}\sqrt{\\{l_{1}^{\prime}\\}\\{l_{2}^{\prime}\\}\\{L\\}\\{L^{\prime}\\}}\\\
\times\left(l_{1}\,l_{2}\,L\right)_{m_{1}m_{2}}^{\phantom{m_{1}m_{2}}M}\left(L\,l_{3}\,L^{\prime}\right)_{Mm_{3}}^{\phantom{Mm_{3}}M^{\prime}}\\\
\times
X^{l_{1}^{\prime}l_{2}^{\prime}}_{L^{\prime}M^{\prime}}(\hat{\bm{k}}_{1},\hat{\bm{k}}_{2})\tilde{B}^{l_{1}l_{2}l_{3};l_{1}^{\prime}l_{2}^{\prime};LL^{\prime}}_{X_{1}X_{2}X_{3}}(k_{1},k_{2}),$
(302)

where

$\tilde{B}^{l_{1}l_{2}l_{3};l_{1}^{\prime}l_{2}^{\prime};LL^{\prime}}_{X_{1}X_{2}X_{3}}(k_{1},k_{2})=(-1)^{l_{1}+l_{2}+l_{3}+l_{1}^{\prime}+l_{2}^{\prime}+L+L^{\prime}}\sqrt{\frac{\\{L\\}\\{L^{\prime}\\}}{\\{l_{1}^{\prime}\\}\\{l_{2}^{\prime}\\}}}\\\
\times\left(l_{1}\,l_{2}\,L\right)^{m_{1}m_{2}}_{\phantom{m_{1}m_{2}}M}\left(L\,l_{3}\,L^{\prime}\right)^{Mm_{3}}_{\phantom{Mm_{3}}M^{\prime}}\left(L^{\prime}\,l_{1}^{\prime}\,l_{2}^{\prime}\right)^{M^{\prime}m_{1}^{\prime}m_{2}^{\prime}}\\\
\times\tilde{B}^{l_{1}l_{2}l_{3};l_{1}^{\prime}l_{2}^{\prime}}_{X_{1}X_{2}X_{3}m_{1}m_{2}m_{3};m_{1}^{\prime}m_{2}^{\prime}}(k_{1},k_{2})$
(303)

is the invariant bispectrum.

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

$B^{l_{1}l_{2}l_{3};l_{1}^{\prime}l_{2}^{\prime};LL^{\prime}\,*}_{X_{1}X_{2}X_{3}}(k)=B^{l_{1}l_{2}l_{3};l_{1}^{\prime}l_{2}^{\prime};LL^{\prime}}_{X_{1}X_{2}X_{3}}(k),$
(304)

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

$\tilde{B}^{(l_{1}l_{2}l_{3})}_{X_{1}X_{2}X_{3}m_{1}m_{2}m_{3}}(\bm{k}_{1},\bm{k}_{2})\xrightarrow{\mathbb{P}}(-1)^{p_{X_{1}}+p_{X_{2}}+p_{X_{3}}+l_{1}+l_{2}+l_{3}}\\\
\times\tilde{B}^{(l_{1}l_{2}l_{3})}_{X_{1}X_{2}X_{3}m_{1}m_{2}m_{3}}(-\bm{k}_{1},-\bm{k}_{2}),$
(305)

because of Eqs. (98) and (298). Accordingly, the parity transformation of the
invariant bispectrum is given by

$\tilde{B}^{l_{1}l_{2}l_{3};l_{1}^{\prime}l_{2}^{\prime}}_{X_{1}X_{2}X_{3}}(k_{1},k_{2})\xrightarrow{\mathbb{P}}(-1)^{p_{X_{1}}+p_{X_{2}}+p_{X_{3}}+l_{1}+l_{2}+l_{3}+l_{1}^{\prime}+l_{2}^{\prime}}\\\
\times\tilde{B}^{l_{1}l_{2}l_{3};l_{1}^{\prime}l_{2}^{\prime}}_{X_{1}X_{2}X_{3}}(k_{1},k_{2}).$
(306)

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

$p_{X_{1}}+p_{X_{2}}+p_{X_{3}}+l_{1}+l_{2}+l_{3}+l_{1}^{\prime}+l_{2}^{\prime}=\mathrm{even}.$
(307)

####  IV.4.2  The lowest-order bispectrum

We consider the bispectrum of the lowest-order approximation or the tree-level
bispectrum. The formal expression of the lowest-order bispectrum in terms of
propagators is given by

$\tilde{B}^{(l_{1}l_{2}l_{3})}_{X_{1}X_{2}X_{3}\,m_{1}m_{2}m_{3}}(\bm{k}_{1},\bm{k}_{2})=\Gamma^{(1)}_{X_{1}l_{1}m_{1}}(\bm{k}_{1})\Gamma^{(1)}_{X_{2}l_{2}m_{2}}(\bm{k}_{2})\\\
\times\Gamma^{(2)}_{X_{3}l_{3}m_{3}}(-\bm{k}_{1},-\bm{k}_{2})P_{\mathrm{L}}(k_{1})P_{\mathrm{L}}(k_{2}).$
(308)

We do not use the higher-order resummation factor and substitute $\Pi(k)=1$.
Substituting Eqs. (159) and (161), and applying Eq. (208), the above equation
is represented in terms of invariant functions as

$\tilde{B}^{(l_{1}l_{2}l_{3})}_{X_{1}X_{2}X_{3}\,m_{1}m_{2}m_{3}}(\bm{k}_{1},\bm{k}_{2})=\frac{(-1)^{l_{3}}}{\sqrt{\\{l_{1}\\}\\{l_{2}\\}\\{l_{3}\\}}}P_{\mathrm{L}}(k_{1})P_{\mathrm{L}}(k_{2})\\\
\times\hat{\Gamma}^{(1)}_{X_{1}l_{1}}(k_{1})\hat{\Gamma}^{(1)}_{X_{2}l_{2}}(k_{2})\sum_{l_{1}^{\prime},l_{2}^{\prime},l_{1}^{\prime\prime},l_{2}^{\prime\prime}}(-1)^{l_{1}^{\prime}+l_{2}^{\prime}}\\{l_{1}^{\prime}\\}\\{l_{2}^{\prime}\\}\\\
\times\left[l_{1}\,l_{1}^{\prime}\,l_{1}^{\prime\prime}\right]_{m_{1}}^{\phantom{m_{1}}m_{1}^{\prime}m_{1}^{\prime\prime}}\left[l_{2}\,l_{2}^{\prime}\,l_{2}^{\prime\prime}\right]_{m_{2}}^{\phantom{m_{2}}m_{2}^{\prime}m_{2}^{\prime\prime}}\left(l_{1}^{\prime\prime}\,l_{2}^{\prime\prime}\,l_{3}\right)_{m_{1}^{\prime\prime}m_{2}^{\prime\prime}m_{3}}\\\
\times
C_{l_{1}^{\prime}m_{1}^{\prime}}(\hat{\bm{k}}_{1})C_{l_{2}^{\prime}m_{2}^{\prime}}(\hat{\bm{k}}_{2})\hat{\Gamma}^{(2)\,l_{3}}_{X_{3}l_{1}^{\prime\prime}l_{2}^{\prime\prime}}(k_{1},k_{2}).$
(309)

Substituting Eqs. (194) and (209), the above expression is given by invariant
forms of renormalized bias functions, $c^{(0)}_{X}$, $c^{(1)}_{X\,l}(k)$ and
$c^{(2)\,l}_{Xl_{1}l_{2}}(k_{1},k_{2})$.

The sum of the products of three $3j$-symbols appeared in the above expression
is represented by a $9j$-symbol using the formulas of Eqs. (411) and (427) in
Appendix C. Comparing the resulting equation with Eq. (302), the invariant
bispectrum is derived, and the result is given by

$\tilde{B}^{l_{1}l_{2}l_{3};l_{1}^{\prime}l_{2}^{\prime};LL^{\prime}}_{X_{1}X_{2}X_{3}}(k_{1},k_{2})=(-1)^{L+L^{\prime}}\sqrt{\frac{\\{l_{1}^{\prime}\\}\\{l_{2}^{\prime}\\}\\{L\\}\\{L^{\prime}\\}}{\\{l_{1}\\}\\{l_{2}\\}\\{l_{3}\\}}}\\\
\times\hat{\Gamma}^{(1)}_{X_{1}l_{1}}(k_{1})\hat{\Gamma}^{(1)}_{X_{2}l_{2}}(k_{2})P_{\mathrm{L}}(k_{1})P_{\mathrm{L}}(k_{2})\sum_{l_{1}^{\prime\prime},l_{2}^{\prime\prime}}(-1)^{l_{1}^{\prime\prime}+l_{2}^{\prime\prime}}\\\
\times\begin{pmatrix}l_{1}&l_{1}^{\prime}&l_{1}^{\prime\prime}\\\
0&0&0\end{pmatrix}\begin{pmatrix}l_{2}&l_{2}^{\prime}&l_{2}^{\prime\prime}\\\
0&0&0\end{pmatrix}\begin{Bmatrix}l_{1}&l_{1}^{\prime}&l_{1}^{\prime\prime}\\\
l_{2}&l_{2}^{\prime}&l_{2}^{\prime\prime}\\\
L&L^{\prime}&l_{3}\end{Bmatrix}\hat{\Gamma}^{(2)\,l_{3}}_{X_{3}l_{1}^{\prime\prime}l_{2}^{\prime\prime}}(k_{1},k_{2}).$
(310)

This is one of the generic predictions of our theory.

As a consistency check, let us see if the bispectrum of the scalarly biased
field in the lowest-order perturbation theory can be correctly reproduced.
Considering the case of $l_{1}=l_{2}=l_{3}=m_{1}=m_{2}=m_{3}=0$ and
$X_{1}=X_{2}=X_{3}\equiv X$ in Eq. (309), and using Eqs. (410), (416) and
(378), we straightforwardly derive

$\tilde{B}^{(000)}_{XXX\,000}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})=P_{\mathrm{L}}(k_{1})P_{\mathrm{L}}(k_{2})\sum_{l}\frac{(-1)^{l}}{\sqrt{\\{l\\}}}\mathit{P}_{l}(\hat{\bm{k}}_{1}\cdot\hat{\bm{k}}_{2})\\\
\times\hat{\Gamma}^{(1)}_{X0}(k_{1})\hat{\Gamma}^{(1)}_{X0}(k_{2})\hat{\Gamma}^{(2)\,0}_{Xll}(k_{1},k_{2}).$
(311)

The same result is also derived from Eqs. (302) and (310) with Eq. (247). We
substitute Eqs. (195) and (209) in the above equation, and straightforward
calculations derive the scalar bispectrum as

$B_{X}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})=B^{(000)}_{XXX\,000}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})\\\
=b_{1}(k_{1})b_{1}(k_{2})P_{\mathrm{L}}(k_{1})P_{\mathrm{L}}(k_{2})\Biggl{\\{}b^{\mathrm{L}}_{2}(\bm{k}_{1},\bm{k}_{2})+b_{1}(k_{1})+b_{1}(k_{2})\\\
-\frac{4}{7}+\left[\frac{k_{2}}{k_{1}}b_{1}(k_{1})+\frac{k_{1}}{k_{2}}b_{1}(k_{2})\right]\frac{\bm{k}_{1}\cdot\bm{k}_{2}}{k_{1}k_{2}}+\frac{4}{7}\left(\frac{\bm{k}_{1}\cdot\bm{k}_{2}}{k_{1}k_{2}}\right)^{2}\Biggr{\\}}\\\
+\mathrm{cyc.},$ (312)

where $b_{1}(k)\equiv c^{(0)}_{X}+c^{(1)}_{X0}(k)$ corresponds to the Eulerian
bias factor and is the same as $b_{X}(k)$ defined in Eq. (241), and

$b^{\mathrm{L}}_{2}(\bm{k}_{1},\bm{k}_{2})\equiv
c^{(2)}_{X00}(\bm{k}_{1},\bm{k}_{2})=\sum_{l}\frac{(-1)^{l}}{\sqrt{\\{l\\}}}\mathit{P}_{l}(\hat{\bm{k}}_{1}\cdot\hat{\bm{k}}_{2})c^{(2)0}_{Xll}(k_{1},k_{2})$
(313)

corresponds to the Lagrangian (nonlocal) bias function of second order [17].
More commonly known form of the bispectrum in the lowest-order perturbation
theory is represented by the Eulerian bias parameters $b_{1}$, $b_{2}$,
$b_{s^{2}}$ [67, 68, 69, 70]. It is natural to take a normalization $\langle
F^{\mathrm{L}}_{X}\rangle=1$ [c.f., Eq. (49)] and we have $c^{(0)}_{X}=1$ in
this scalar case. On one hand, in the special case of local bias, in which
$b_{1}$ and $b^{\mathrm{L}}_{2}$ are constants and do not depend on
wavevectors, Eq. (312) reduces to

$B_{X}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})={b_{1}}^{2}P_{\mathrm{L}}(k_{1})P_{\mathrm{L}}(k_{2})\\\
\times\left[b^{\mathrm{L}}_{2}+2b_{1}-\frac{4}{7}+b_{1}\left(\frac{k_{2}}{k_{1}}+\frac{k_{1}}{k_{2}}\right)\frac{\bm{k}_{1}\cdot\bm{k}_{2}}{k_{1}k_{2}}+\frac{4}{7}\left(\frac{\bm{k}_{1}\cdot\bm{k}_{2}}{k_{1}k_{2}}\right)^{2}\right]\\\
+\mathrm{cyc.}$ (314)

On the other hand, the prediction of scalar bispectrum in the standard
(Eulerian) perturbation theory is given by [70],

$B_{X}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})={b_{1}}^{2}P_{\mathrm{L}}(k_{1})P_{\mathrm{L}}(k_{2})\\\
\times\left\\{b_{2}+2b_{s^{2}}\left[\left(\frac{\bm{k}_{1}\cdot\bm{k}_{2}}{k_{1}k_{2}}\right)^{2}-\frac{1}{3}\right]+b_{1}F_{2}(\bm{k}_{1},\bm{k}_{2})\right\\}\\\
+\mathrm{cyc.},$ (315)

where

$F_{2}(\bm{k}_{1},\bm{k}_{2})=\frac{10}{7}+\left(\frac{k_{2}}{k_{1}}+\frac{k_{1}}{k_{2}}\right)\frac{\bm{k}_{1}\cdot\bm{k}_{2}}{k_{1}k_{2}}+\frac{4}{7}\left(\frac{\bm{k}_{1}\cdot\bm{k}_{2}}{k_{1}k_{2}}\right)^{2}$
(316)

is the second-order kernel444Our choice of the normalizations for kernel
functions $F_{n}$ is different from those in many previous references, where
our $F_{n}$ corresponds to $n!F_{n}$ in the latter. We believe our choice
simplifies the expressions of equations in general [17]. of the Eulerian
perturbation theory [8]. Comparing the coefficients of Eqs. (314) and (315),
we find the relations between Lagrangian local bias parameters and Eulerian
local bias parameters as

$b_{1}=1+b^{\mathrm{L}}_{1},\quad
b_{2}=b^{\mathrm{L}}_{2}-\frac{8}{21}b^{\mathrm{L}}_{1},\quad
b_{s^{2}}=-\frac{2}{7}b^{\mathrm{L}}_{1}.$ (317)

The above relations exactly agree with those which are found in the literature
[69, 70]. It is also possible to include the nonlocal Lagrangian bias
$b^{\mathrm{L}}_{s^{2}}$ in $b^{\mathrm{L}}_{2}(\bm{k}_{1},\bm{k}_{2})$, and
derive consistent results with those in the literature [69, 70].

##  V  Semi-local models of bias for tensor fields

So far the formulation is general enough, and we have not assumed any model of
bias for the tensor fields so far. The bias model is naturally specified in
Lagrangian space through the renormalized bias functions $c^{(n)}_{Xlm}$
defined by Eq. (70), which is evaluated once a model of bias
$F^{\mathrm{L}}_{Xlm}(\bm{q})$ is analytically given as a functional of the
linear density field $\delta_{\mathrm{L}}(\bm{q})$. In this section, we
consider a category of bias models which are commonly adopted in models of
cosmological structure formation. We call this category “semi-local models” of
bias.

###  V.1 The semi-local models of bias

####  V.1.1 Defining the semi-local models

The concept of the semi-local models is somehow similar to the EFTofLSS
approach in perturbation theory of biased tracers [71, 4], where the biased
field is given by a finite set of semi-local operators made of spatial and
temporal derivatives of the gravitational potential in the perturbative
expansions order by order. However, while the EFTofLSS approach is based on
phenomenologically perturbative expansions of the biased field, our iPT
approach does not assume the underlying operators are perturbative, because we
essentially use orthogonal expansions instead of Taylor expansions to include
the fully nonlinear biasing into the cosmological perturbation theory through
the renormalized bias functions.

This category of semi-local models of bias is defined so that the field
$F^{\mathrm{L}}_{Xlm}(\bm{q})$ is given by functions (instead of functionals)
of the spatial derivatives of the gravitational potential at the same position
$\bm{q}$, smoothed with (generally, multiple numbers of) window functions
$W^{(a)}(k)$:

$\chi^{(a)}_{i_{1}i_{2}\cdots
i_{L_{a}}}(\bm{q})=\partial_{i_{1}}\partial_{i_{2}}\cdots\partial_{i_{L_{a}}}\psi^{(a)}(\bm{q}),$
(318)

where

$\psi^{(a)}(\bm{q})=\int\frac{d^{3}k}{(2\pi)^{3}}e^{i\bm{k}\cdot\bm{q}}\delta_{\mathrm{L}}(\bm{k})\,k^{-L_{a}}W^{(a)}(k)$
(319)

is the smoothed linear density field with an isotropic window function
$k^{-L_{a}}W^{(a)}(k)$. The label “$a$” distinguishes different kinds of the
linear tensor field with a particular rank $L_{a}$, i.e., the label $a$
uniquely specifies the rank $L_{a}$ and the form of the window function
$W^{(a)}(k)$. The window function usually contains a parameter of smoothing
radius $R_{a}$, which can take different values for each window function
$W^{(a)}$. While typical window functions include the Gaussian window,
$W^{(a)}(k)=\exp(-k^{2}{R_{a}}^{2}/2)$, and the top-hat window,
$W^{(a)}(k)=3j_{1}(kR_{a})/(kR_{a})$, we do not assume any special form of the
window function here. When only a single linear tensor field of a fixed window
function is considered, one can omit the label $a$ in the above. We include
the possibility of using multiple numbers of fields and window functions in
the general formulation below.

For a concrete example of biasing using multiple numbers of fields and window
functions, see Ref. [20]. For example, the second derivatives of the
(normalized and smoothed) linear gravitational potential $\varphi$ correspond
to $\varphi=\mathop{}\\!\mathbin{\bigtriangleup}^{-1}\delta_{R}=-\psi^{(a)}$
with a rank of $L_{a}=2$, and
$\chi^{(a)}_{ij}=-\partial_{i}\partial_{j}\mathop{}\\!\mathbin{\bigtriangleup}^{-1}\delta_{R}$
where $\delta_{R}$ is a smoothed linear density field with a smoothing window
function $W^{(a)}(k)=W(kR)$.

The linear density contrast $\delta_{\mathrm{L}}$ is a real function in
configuration space, and thus we have
$\delta_{\mathrm{L}}^{*}(\bm{k})=\delta_{\mathrm{L}}(-\bm{k})$ as shown from
Eq. (61). Thus we see the function $\psi^{(a)}(\bm{q})$ of Eq. (319) is a real
function, and so as the function $\chi^{(a)}_{i_{1}\cdots i_{L_{a}}}(\bm{q})$
of Eq. (318).

In the semi-local models of bias, the field value
$F^{\mathrm{L}}_{Xlm}(\bm{q})$ at a particular position is given by a function
of the field values $\chi^{(a)}_{i_{1}\cdots i_{l}}(\bm{q})$ at the same
position in Lagrangian space. This function is common among all the positions.
Therefore, without loss of generality, one can consider the particular
position at the origin, $\bm{q}=\bm{0}$, to describe the relation in the semi-
local model of bias, and thus the field $F^{\mathrm{L}}_{Xlm}$ is a function
of

$\chi^{(a)}_{i_{1}\cdots
i_{L_{a}}}=i^{L_{a}}\int\frac{d^{3}k}{(2\pi)^{3}}\,\hat{k}_{i_{1}}\cdots\hat{k}_{i_{L_{a}}}\,\delta_{\mathrm{L}}(\bm{k})W^{(a)}(k).$
(320)

The tensor field $F^{\mathrm{L}}_{Xlm}$ can depend on multiple of local
tensors of various ranks, $\chi^{(a_{0})}$, $\chi^{(a_{1})}_{i}$,
$\chi^{(a_{2})}_{ij}$, $\chi^{(a_{3})}_{ijk}$, etc.

####  V.1.2 Irreducible decomposition of semi-local variables

A local tensor is decomposed by irreducible spherical basis according to the
the procedure explained in Sec. II. The tensor of Eq. (320) is decomposed into
traceless tensor of rank $L_{a},L_{a}-2,\ldots$, and each component of
traceless tensor is decomposed by a spherical basis, as in the Eqs. (35)–(42).
The decomposition is given by

$\chi^{(a)}_{i_{1}\cdots i_{l}}=\chi^{(a,l)}_{i_{1}i_{2}\cdots
i_{l}}+\frac{l(l-1)}{2(2l-1)}\delta_{(i_{1}i_{2}}\chi^{(a,l-2)}_{i_{3}\cdots
i_{l})}+\cdots,$ (321)

where we simply denote $L_{a}=l$ in the above, and
$\chi^{(a,l-2)}_{i_{3}\cdots i_{l}}$ is the traceless part of the first trace
$\chi^{(a)}_{jji_{3}\cdots i_{l}}$ and so forth.

The decomposed traceless tensors are further decomposed by spherical basis.
Similarly to Eq. (41), we have

$\displaystyle\chi^{(a,l)}_{i_{1}i_{2}\cdots i_{l}}$
$\displaystyle=i^{l}A_{l}\,\chi^{(a[l])}_{lm}\mathsf{Y}^{(m)}_{i_{1}i_{2}\cdots
i_{l}},$ (322) $\displaystyle\chi^{(a,l-2)}_{i_{3}\cdots i_{l}}$
$\displaystyle=i^{l-2}A_{l-2}\,\chi^{(a[l])}_{l-2,m}\mathsf{Y}^{(m)}_{i_{3}\cdots
i_{l}},$ (323)

and so forth. In our simplified notation with $l=L_{a}$, the rank of the
linear tensor field $a$ is not obvious if we write, e.g., $\chi^{(a)}_{l-2,m}$
for the first trace of the original tensor, and thus we instead employ a
notations $\chi^{(a[l])}_{lm}$, $\chi^{(a[l])}_{l-2,m}$, etc. to remind the
rank of the original field $a$ is $l$ in the lower-rank trace parts of the
tensor. Namely,

$\chi^{(a[l])}_{lm}=\left.\chi^{(a)}_{L_{a}m}\right|_{L_{a}=l},\quad\chi^{(a[l])}_{l-2,m}=\left.\chi^{(a)}_{L_{a}-2,m}\right|_{L_{a}=l},$
(324)

and so forth.

The spherical tensors are given by

$\displaystyle\chi^{(a[l])}_{lm}$
$\displaystyle=\frac{(-i)^{l}}{A_{l}}\,\chi^{(a)}_{i_{1}\cdots
i_{l}}\mathsf{Y}^{(m)*}_{i_{1}\cdots i_{l}}$
$\displaystyle=(-i)^{l}\int\frac{d^{3}k}{(2\pi)^{3}}\delta_{\mathrm{L}}(\bm{k})C_{lm}(\hat{\bm{k}})W^{(a)}(k),$
(325) $\displaystyle\chi^{(a[l])}_{l-2,m}$
$\displaystyle=\frac{(-i)^{l-2}}{A_{l-2}}\,\chi^{(a)}_{jji_{3}\cdots
i_{l}}\mathsf{Y}^{(m)*}_{i_{3}\cdots i_{l}}$
$\displaystyle=(-i)^{l-2}\int\frac{d^{3}k}{(2\pi)^{3}}\delta_{\mathrm{L}}(\bm{k})C_{l-2,m}(\hat{\bm{k}})W^{(a)}(k),$
(326)

where $A_{l}$ is given by Eq. (19), and $C_{lm}(\hat{\bm{k}})$ is the
spherical harmonics with Racah’s normalization. In general, we have

$\displaystyle\chi^{(a[l])}_{l-2p,m}$
$\displaystyle=\frac{(-i)^{l-2p}}{A_{l-2p}}\,\chi^{(a)}_{j_{1}j_{1}\cdots
j_{p}j_{p}i_{2p+1}\cdots i_{l}}\mathsf{Y}^{(m)*}_{i_{2p+1}\cdots i_{l}}$
$\displaystyle=(-i)^{l-2p}\int\frac{d^{3}k}{(2\pi)^{3}}\delta_{\mathrm{L}}(\bm{k})C_{l-2p,m}(\hat{\bm{k}})W^{(a)}(k),$
(327)

where $p=0,1,\cdots,[l/2]$. Changing the labels of ranks by $l\rightarrow L$
and $l-2p\rightarrow l$, Eq. (327) is equivalently represented by

$\chi^{(a[L])}_{lm}=(-i)^{l}\int\frac{d^{3}k}{(2\pi)^{3}}\delta_{\mathrm{L}}(\bm{k})C_{lm}(\hat{\bm{k}})W^{(a)}(k),$
(328)

where $l=L,L-2$, …, (0 or 1) for a linear tensor field $a$ of rank $L$. The
smallest value of $l$ in the last equation is 0 if $L$ is even, and 1 if $L$
is odd. When the rank $L_{a}$ of the original tensor $\chi^{(a)}_{i_{1}\cdots
i_{L_{a}}}$ is obvious, one can employ a simplified notation $\chi^{(a)}_{lm}$
instead of $\chi^{(a[L])}_{lm}$, where $l=L,L-2,\ldots$

####  V.1.3 Simple examples

For a simple example, we consider the second derivatives of the normalized
potential $\varphi$, smoothed with a window function $W(kR)$. In this case, we
have

$\chi^{(\varphi)}_{ij}=-\partial_{i}\partial_{j}\varphi=-\int\frac{d^{3}k}{(2\pi)^{3}}\,\hat{k}_{i}\hat{k}_{j}\delta_{\mathrm{L}}(\bm{k})W(kR).$
(329)

The general variables in Eqs. (319) and (320) correspond to
$\psi^{(a)}=-\varphi$ and $W^{(\varphi)}(k)=W(kR)$. The second-rank tensor is
decomposed into irreducible components as

$\displaystyle\chi^{(\varphi)}_{ij}$
$\displaystyle=-\left(\partial_{i}\partial_{j}-\frac{\delta_{ij}}{3}\mathop{}\\!\mathbin{\bigtriangleup}\right)\varphi-\frac{\delta_{ij}}{3}\mathop{}\\!\mathbin{\bigtriangleup}\varphi$
$\displaystyle=-\sqrt{\frac{2}{3}}\,\chi^{(\varphi)}_{2m}\,\mathsf{Y}^{(m)}_{ij}-\frac{\delta_{ij}}{3}\,\chi^{(\varphi)}_{00}\,\mathsf{Y}^{(0)},$
(330)

where, corresponding to Eqs. (325) and (326), we have

$\displaystyle\chi^{(\varphi[2])}_{2m}$
$\displaystyle=-\int\frac{d^{3}k}{(2\pi)^{3}}\delta_{\mathrm{L}}(\bm{k})C_{2m}(\hat{\bm{k}})W(kR),$
(331) $\displaystyle\chi^{(\varphi[2])}_{00}$
$\displaystyle=\int\frac{d^{3}k}{(2\pi)^{3}}\delta_{\mathrm{L}}(\bm{k})C_{00}(\hat{\bm{k}})W(kR).$
(332)

The last variable is apparently given by

$\chi^{(\varphi[2])}_{00}=\int\frac{d^{3}k}{(2\pi)^{3}}\delta_{\mathrm{L}}(\bm{k})W(kR)\equiv\delta_{R},$
(333)

where $\delta_{R}$ is the smoothed density field with a smoothing function
$W(kR)$.

Similarly, the second derivatives of the smoothed density field is given by

$\chi^{(\delta)}_{ij}=\partial_{i}\partial_{j}\delta_{R}=-\int\frac{d^{3}k}{(2\pi)^{3}}\,k^{2}\,\hat{k}_{i}\hat{k}_{j}\delta_{\mathrm{L}}(\bm{k})W(kR).$
(334)

The general variables in Eqs. (319) and (320) correspond to
$\psi^{(a)}=\delta_{R}$ and $W^{(\delta)}(k)=k^{2}W(kR)$. The second-rank
tensor is decomposed into irreducible components as

$\displaystyle\chi^{(\delta)}_{ij}=\partial_{i}\partial_{j}\delta_{R}$
$\displaystyle=-\left(\partial_{i}\partial_{j}-\frac{\delta_{ij}}{3}\mathop{}\\!\mathbin{\bigtriangleup}\right)\delta_{R}-\frac{\delta_{ij}}{3}\mathop{}\\!\mathbin{\bigtriangleup}\delta_{R}$
$\displaystyle=-\sqrt{\frac{2}{3}}\chi^{(\delta[2])}_{2m}\mathsf{Y}^{(m)}_{ij}-\frac{\delta_{ij}}{3}\chi^{(\delta[2])}_{00}\mathsf{Y}^{(0)},$
(335)

where

$\displaystyle\chi^{(\delta[2])}_{2m}$
$\displaystyle=-\int\frac{d^{3}k}{(2\pi)^{3}}\delta_{\mathrm{L}}(\bm{k})C_{2m}(\hat{\bm{k}})\,k^{2}W(kR),$
(336) $\displaystyle\chi^{(\delta[2])}_{00}$
$\displaystyle=\int\frac{d^{3}k}{(2\pi)^{3}}\delta_{\mathrm{L}}(\bm{k})C_{00}(\hat{\bm{k}})\,k^{2}W(kR),$
(337)

and the last variable is apparently given by

$\chi^{(\delta[2])}_{00}=\int\frac{d^{3}k}{(2\pi)^{3}}\delta_{\mathrm{L}}(\bm{k})\,k^{2}W(kR)=-\mathop{}\\!\mathbin{\bigtriangleup}\delta_{R},$
(338)

where $\mathop{}\\!\mathbin{\bigtriangleup}\delta_{R}$ is the Laplacian of the
smoothed density field.

###  V.2 Renormalized bias parameters

####  V.2.1 Renormalized bias functions in the semi-local models

In the semi-local models of bias, the tensor field $F^{\mathrm{L}}_{Xlm}$ is
generally considered as a function of various irreducible tensors

$F^{\mathrm{L}}_{Xlm}=F^{\mathrm{L}}_{Xlm}\left(\left\\{\chi^{(a)}_{l^{\prime}m^{\prime}}\right\\}\right)$
(339)

at every position in (Lagrangian) configuration space. In this case, the
functional derivative in the definition of renormalized bias functions, Eq.
(70), can be replaced by

$(2\pi)^{3}\frac{\delta}{\delta\delta_{\mathrm{L}}(\bm{k})}\rightarrow\sum_{a}W^{(a)}(k)\sum_{l=L_{a},L_{a}-2,\ldots}C_{lm}(\hat{\bm{k}})\frac{\partial}{\partial\chi^{(a)}_{lm}},$
(340)

where the integer $L_{a}$ is the original rank of the linear tensor
$\chi^{(a)}_{i_{1}\cdots i_{L_{a}}}$. The summation over $m$ in Eq. (340) is
implicitly assumed according to the Einstein summation convention, just as in
the rest of this paper. Therefore, the renormalized bias functions of Eq. (70)
are given by

$c^{(n)}_{Xlm}(\bm{k}_{1},\ldots\bm{k}_{n})=\sum_{a_{1},\ldots,a_{n}}W^{(a_{1})}(k_{1})\cdots
W^{(a_{n})}(k_{n})\\\
\times\sum_{l_{1},\ldots,l_{n}}C_{l_{1}m_{1}}(\hat{\bm{k}}_{1})\cdots
C_{l_{n},m_{n}}(\hat{\bm{k}}_{n})\left\langle\frac{\partial^{n}F^{\mathrm{L}}_{Xlm}}{\partial\chi^{(a_{1})}_{l_{1}m_{1}}\cdots\partial\chi^{(a_{n})}_{l_{n}m_{n}}}\right\rangle,$
(341)

where integers $l_{i}$ ($i=1,\ldots,n$) run over
$l_{i}=L_{a_{i}},L_{a_{i}}-2,\ldots$, (0 or 1). Comparing the above equation
with Eq. (107), we have

$c^{(n)\,l;l_{1}\cdots l_{n}}_{Xm;m_{1}\cdots
m_{n}}(k_{1},\ldots,k_{n})=\sum_{a_{1},\ldots,a_{n}}W^{(a_{1})}(k_{1})\cdots
W^{(a_{n})}(k_{n})\\\ \times g^{(l_{1})}_{m_{1}m_{1}^{\prime}}\cdots
g^{(l_{n})}_{m_{n}m_{n}^{\prime}}\left\langle\frac{\partial^{n}F^{\mathrm{L}}_{Xlm}}{\partial\chi^{(a_{1})}_{l_{1}m_{1}^{\prime}}\cdots\partial\chi^{(a_{n})}_{l_{n}m_{n}^{\prime}}}\right\rangle.$
(342)

The renormalized bias functions in the form of Eq. (341) is evaluated for a
given model of the biased tensor field, once the underlying statistics of the
fields $\chi^{(a)}_{lm}$ are specified. These fields just linearly depend on
the linear density contrast $\delta_{\mathrm{L}}$ through Eqs. (325)–(328) and
therefore the statistics are straightforwardly given by those of linear
density field.

Using Eqs. (74) and (328), and the orthonormality relation of the spherical
harmonics, Eq. (380), the covariance of the fields (at the same position, the
same applies hereafter) can be straightforwardly calculated and is given by

$\left\langle\chi^{(a)}_{lm}\chi^{(b)}_{l^{\prime}m^{\prime}}\right\rangle=\delta_{ll^{\prime}}\gamma^{ab}_{(l)}g^{(l)}_{mm^{\prime}}.$
(343)

where

$\gamma^{ab}_{(l)}\equiv\frac{1}{2l+1}\int\frac{k^{2}dk}{2\pi^{2}}W^{(a)}(k)W^{(b)}(k)P_{\mathrm{L}}(k).$
(344)

One can regard the set of parameters $\gamma^{ab}_{(l)}$ as matrix elements of
a matrix $\bm{\gamma}_{(l)}$, which components are given by
$[\bm{\gamma}_{(l)}]_{ab}=\gamma^{ab}_{(l)}$. We denote the matrix elements of
the inverse matrix as

$\gamma^{(l)}_{ab}\equiv\left[\bm{\gamma}_{(l)}^{-1}\right]_{ab}.$ (345)

When one considers a set of variables $\chi^{(a)}_{lm}$ as a vector with the
set of indices $(a,l,m)$, the inverse of the covariance matrix of Eq. (343) is
given by $\delta_{ll^{\prime}}\gamma^{(l)}_{ab}g_{(l)}^{mm^{\prime}}$ with our
notations.

When the initial density field $\delta_{\mathrm{L}}$ is a Gaussian random
field, the two-point covariance of Eq. (343) contains all the information for
the statistical distribution of the variables. In this case, the distribution
function is given by a multivariate Gaussian distribution function,

$P_{\mathrm{G}}\left(\left\\{\chi^{(a)}_{lm}\right\\}\right)=\frac{1}{\sqrt{(2\pi)^{N}\det\bm{\gamma}_{(l)}}}\\\
\times\exp\left[-\frac{1}{2}\sum_{l}\gamma^{(l)}_{ab}\,g_{(l)}^{mm^{\prime}}\chi^{(a)}_{lm}\chi^{(b)}_{lm^{\prime}}\right],$
(346)

where the Einstein summation convention is applied also to the indices $a$,
$b$ as well as azimuthal indices $m$, $m^{\prime}$ and thus the summation over
these indices are implicitly assumed in the exponent. While the variables
$\chi^{(a)}_{lm}$ are complex numbers, the exponent of Eq. (346) is a real
number. This can be readily shown by noting a property
$\chi^{(a)\,*}_{lm}=(-1)^{l}g_{(l)}^{mm^{\prime}}\chi^{(a)}_{lm^{\prime}}$,
which is shown from Eq. (328) with
$\delta_{\mathrm{L}}^{*}(\bm{k})=\delta_{\mathrm{L}}(-\bm{k})$ and Eqs. (376)
and (377). We have

$\gamma^{(l)}_{ab}\,g_{(l)}^{mm^{\prime}}\chi^{(a)}_{lm}\chi^{(b)}_{lm^{\prime}}=(-1)^{l}\gamma^{(l)}_{ab}\,\chi^{(a)}_{lm}\chi^{(b)\,*}_{lm},$
(347)

which is a real number because the coefficients $\gamma^{(l)}_{ab}$ are real
numbers and symmetric with respect to the indices $a$ and $b$. Moreover, the
above factor is positive as long as the matrix $\bm{\gamma}_{(l)}$ is positive
definite.

####  V.2.2 Defining the renormalized bias parameters

The expectation value of Eq. (341) is calculated by multivariate Gaussian
integrals with the distribution function of Eq. (346), which evaluation is
analytically possible in many cases when the semi-local bias function of Eq.
(339) is given by an analytic function of a finite number of variables. When
there is a small non-Gaussianity in the initial condition, and the linear
density field is not exactly a random Gaussian field, one can evaluate non-
Gaussian corrections to the Gaussian distribution function of Eq. (346). The
procedure of deriving the non-Gaussian corrections are found in Ref. [23, 72],
which is straightforward to apply in this case. However, for the illustrative
examples below of this paper, we do not need to include the non-Gaussian
corrections in the evaluations of renormalized bias functions.

The expectation value on the rhs of Eqs. (341) and (342) should be represented
by rotationally invariant variables, just in the case of the renormalized bias
function $c^{(n)l;l_{1}\cdots l_{n}}_{Xm;m_{1}\cdots
m_{n}}(k_{1},\ldots,k_{n})$, while the arguments of $k_{i}$’s are not present
here. Similarly to Eqs. (125), (129), (139), (150), we have

$\displaystyle\left\langle\frac{\partial
F^{\mathrm{L}}_{Xlm}}{\partial\chi^{(a_{1})}_{l_{1}m_{1}}}\right\rangle$
$\displaystyle=\frac{(-1)^{l}}{\sqrt{2l+1}}\delta_{ll_{1}}\delta_{m}^{m_{1}}b^{(1:a_{1})}_{Xl},$
(348)
$\displaystyle\left\langle\frac{\partial^{2}F^{\mathrm{L}}_{Xlm}}{\partial\chi^{(a_{1})}_{l_{1}m_{1}}\partial\chi^{(a_{2})}_{l_{2}m_{2}}}\right\rangle$
$\displaystyle=\left(l\,l_{1}\,l_{2}\right)_{m}^{\phantom{m}m_{1}m_{2}}b^{(2:a_{1}a_{2})}_{Xl;l_{1}l_{2}},$
(349)
$\displaystyle\left\langle\frac{\partial^{3}F^{\mathrm{L}}_{Xlm}}{\partial\chi^{(a_{1})}_{l_{1}m_{1}}\partial\chi^{(a_{2})}_{l_{2}m_{2}}\partial\chi^{(a_{3})}_{l_{3}m_{3}}}\right\rangle$
$\displaystyle=\sum_{L}(-1)^{L}\sqrt{\\{L\\}}\left(l\,l_{1}\,L\right)_{m}^{\phantom{m}m_{1}M}$
$\displaystyle\hskip
12.0pt\times\left(L\,l_{2}\,l_{3}\right)_{M}^{\phantom{M}m_{2}m_{3}}b^{(3:a_{1}a_{2}a_{3})L}_{Xl;l_{1}l_{2}l_{3}},$
(350)

and

$\left\langle\frac{\partial^{3}F^{\mathrm{L}}_{Xlm}}{\partial\chi^{(a_{1})}_{l_{1}m_{1}}\cdots\partial\chi^{(a_{n})}_{l_{n}m_{n}}}\right\rangle=\sum_{L_{2},\ldots,L_{n-1}}(-1)^{L_{2}+\cdots
L_{n-1}}\\\
\times\sqrt{\\{L_{2}\\}\cdots\\{L_{n-1}\\}}\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_{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}}b^{(n:a_{1}\cdots
a_{n})L_{2}\cdots L_{n-1}}_{Xl;l_{1}\cdots l_{n}}.$ (351)

The invariant constants $b^{(1;a_{1})}_{Xl}$,
$b^{(2;a_{1}a_{2})}_{Xl;l_{1}l_{2}}$,
$b^{(3;a_{1}a_{3})L}_{Xl;l_{1}l_{2}l_{3}}$ etc. are parameters which
characterize the semi-local models of bias.

Therefore, comparing Eq. (342) with Eqs. (125), (129), (134) and (150),
invariant functions of the renormalized bias functions are given by

$\displaystyle c^{(1)}_{Xl}(k)$
$\displaystyle=\sum_{a}b^{(1:a)}_{Xl}W^{(a)}(k),$ (352) $\displaystyle
c^{(2)\,l}_{Xl_{1}l_{2}}(k_{1},k_{2})$
$\displaystyle=\sum_{a_{1},a_{2}}b^{(2:a_{1}a_{2})}_{Xl;l_{1}l_{2}}W^{(a_{1})}(k_{1})W^{(a_{2})}(k_{2}),$
(353) $\displaystyle c^{(3)\,l;L}_{Xl;l_{1}l_{2}l_{3}}(k_{1},k_{2},k_{3})$
$\displaystyle=\sum_{a_{1},a_{2},a_{3}}b^{(3:a_{1}a_{2}a_{3})L}_{Xl;l_{1}l_{2}l_{3}}$
$\displaystyle\qquad\times
W^{(a_{1})}(k_{1})W^{(a_{2})}(k_{2})W^{(a_{3})}(k_{3}),$ (354)

and

$c^{(n)\,l;L_{2}\cdots L_{n-1}}_{Xl_{1}\cdots
l_{n}}(k_{1},\ldots,k_{n})=\sum_{a_{1},\ldots,a_{n}}b^{(n:a_{1}\cdots
a_{n})L_{2}\cdots L_{n-1}}_{Xl;l_{1}\cdots l_{n}}\\\ \times
W^{(a_{1})}(k_{1})\cdots W^{(a_{n})}(k_{n}).$ (355)

Thus, the renormalized bias functions are given by superpositions of products
of window functions with constant parameters which are determined by a given
semi-local model of bias. We naturally call these parameters as the
“renormalized bias parameters.” The scale dependencies in the invariant
renormalized functions are all contained in the window functions which are
fixed by a construction of the semi-local model with a finite number of scale-
independent parameters. These properties simplify the calculations of loop
corrections in our applications.

Corresponding to the interchange symmetries of the renormalized bias functions
with orders greater than two, Eqs. (138), (148), (149), (157) and (158), the
same symmetries for the constant coefficients of bias are given by

$b^{(2:a_{2}a_{1})}_{Xl;l_{2}l_{1}}=(-1)^{l+l_{1}+l_{2}}b^{(2:a_{1}a_{2})}_{Xl;l_{1}l_{2}}$
(356)

for the second-order bias,

$\displaystyle b^{(3:a_{1}a_{3}a_{2})L}_{Xl;l_{l}l_{3}l_{2}}$
$\displaystyle=(-1)^{l_{2}+l_{3}+L}b^{(3:a_{1}a_{2}a_{3})L}_{Xl;l_{l}l_{2}l_{3}},$
(357) $\displaystyle b^{(3:a_{2}a_{1}a_{3})L}_{Xl;l_{2}l_{1}l_{3}}$
$\displaystyle=(-1)^{l_{1}+l_{2}}\sum_{L^{\prime}}(2L^{\prime}+1)\begin{Bmatrix}l_{1}&l&L\\\
l_{2}&l_{3}&L^{\prime}\end{Bmatrix}b^{(3:a_{1}a_{2}a_{3})L^{\prime}}_{Xl;l_{l}l_{2}l_{3}}$
(358)

for the third-order bias, and

$\displaystyle b^{(n:a_{1}\cdots a_{n-2}a_{n}a_{n-1})L_{2}\cdots
L_{n-2}L_{n-1}}_{Xl;l_{l}\cdots
l_{n-2}l_{n}l_{n-1}}=(-1)^{l_{n-1}+l_{n}+L_{n-1}}b^{(n)L_{2}\cdots
L_{n-1}}_{Xl;l_{l}\cdots l_{n}},$ (359) $\displaystyle b^{(n:a_{1}\cdots
a_{i-1}a_{i+1}a_{i}a_{i+2}\cdots a_{n});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}}$
$\displaystyle\qquad=(-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}$ $\displaystyle\hskip 96.0pt\times
b^{(n:a_{1}\cdots a_{n});L_{2}\cdots L_{i}L^{\prime}L_{i+2}\cdots
L_{n-1}}_{Xl;l_{1}\cdots l_{n}}$ (360)

for higher-order bias with $n>3$, where $i=1,\cdots,n-2$.

####  V.2.3 An example: the scale-dependent bias in tensor fields

As a specific example, there is a loop integral in the scale-dependent bias of
Eq. (291). Applying an expression Eq. (353) of semi-local models, the integral
is given by

$\int\frac{p^{2}dp}{2\pi^{2}}c^{(2)l}_{Xl_{1}l_{2}}(p,p)P_{\mathrm{L}}(p)=\sum_{a_{1},a_{2}}b^{(2:a_{1}a_{2})}_{Xl;l_{1}l_{2}}\sigma^{(a_{1}a_{2})},$
(361)

where

$\sigma^{(a_{1}a_{2})}\equiv\int\frac{k^{2}dk}{2\pi^{2}}P_{\mathrm{L}}(k)W^{(a_{1})}(k)W^{(a_{2})}(k).$
(362)

Thus Eq. (291) reduces to

$\Delta
b_{Xl}(k)=\frac{2f_{\mathrm{NL}}^{(l)}}{\mathcal{M}(k)}\sqrt{\\{l\\}}\sum_{l_{1},l_{2}}(-1)^{l_{2}}\begin{pmatrix}l_{1}&l_{2}&l\\\
0&0&0\end{pmatrix}\sum_{a_{1},a_{2}}b^{(2:a_{1}a_{2})}_{Xl;l_{1}l_{2}}\sigma^{(a_{1}a_{2})},$
(363)

for $l=\mathrm{even}$ and $\Delta b_{Xl}(k)=0$ for $l=\mathrm{odd}$.

For a simple case of scalar bias with the high-mass limit of the halo model,
comparison of Eq. (295) with Eq. (353) shows that only a term with
$b^{(2)\,0}_{\mathrm{h}\delta\delta;00}=\sigma^{-2}\delta_{\mathrm{c}}b^{\mathrm{L}}_{1}$
survives, and one can readily see the result of Eq. (296) is reproduced as a
consistency check.

If we consider a simple model that the tensor bias is a local function of only
the second-order derivatives of the gravitational potential
$\partial_{i}\partial_{j}\varphi$, the window function is given by
$W^{(\varphi)}(k)=W(kR)$ as shown in Eq. (329). We can omit the label
$(\varphi)$ in this case, because the linear tensor field consists of only a
single tensor. Thus Eq. (362) simply reduces to

$\sigma^{2}=\int\frac{k^{2}dk}{2\pi^{2}}P_{\mathrm{L}}(k)W^{2}(kR),$ (364)

which is the variance of the smoothed linear density field. The ranks $l_{1}$
and $l_{2}$ take only values of 0 and 2 in the summation of Eq. (363) in this
case, and due to the $3j$-symbol in Eq. (363), only the cases of $l=0,2,4$ are
nonzero. Substituting concrete numbers of the $3j$-symbols, Eq. (363) is
explicitly expanded in a finite number of terms and the result is given by

$\Delta
b_{Xl}(k)=\frac{2\sigma^{2}}{\mathcal{M}(k)}\Biggl{[}\delta_{l0}f^{(0)}_{\mathrm{NL}}\left(b^{(2)}_{0;00}+\frac{1}{\sqrt{5}}b^{(2)}_{0;22}\right)\\\
+\delta_{l2}f^{(2)}_{\mathrm{NL}}\left(2b^{(2)}_{2;02}-\sqrt{\frac{2}{7}}b^{(2)}_{2;22}\right)\\\
+\delta_{l4}f^{(4)}_{\mathrm{NL}}\,\frac{3\sqrt{2}}{\sqrt{35}}\,b^{(2)}_{4;22}\Biggr{]},$
(365)

where we use an interchange symmetry $b^{(2)}_{2;20}=b^{(2)}_{2;02}$ derived
from Eq. (356). Therefore, if the tensor bias is modeled as a local function
of only a second-rank linear tensor field, the scale-dependent bias of the
tensor field is nonzero only when the rank of the biased tensor is 0, 2, or 4.
Due to the rotational symmetry, only a small number of bias parameters,
$b^{(2)}_{0;00}$, $b^{(2)}_{0;22}$, $b^{(2)}_{2;02}$, $b^{(2)}_{2;22}$, and
$b^{(2)}_{4;22}$ appear in this particular semi-local model.

The above examples probably do not sufficiently exhibit the practical merits
of semi-local models of bias. Their advantages are more apparently shown in
the calculation of loop corrections in the higher-order perturbation theory.
The present formalism with the full use of the spherical basis is quite
compatible with the FFT-PT framework [73, 74] or the FAST-PT framework [75,
76], which dramatically reduces the dimensionality of multidimensional
integrals in the calculation of loop corrections in the higher-order
perturbation theory. Technical details of calculating the loop corrections in
the present formalism will be considered in a separate paper of the series,
Paper II [43].

##  VI  Conclusions

In this paper, the formalism iPT is generalized to calculate statistics of
generally tensor-valued objects, using the nonlinear perturbation theory.
Higher-rank tensors are conveniently decomposed into spherical tensors which
are irreducible representations of the three-dimensional rotation group SO(3),
and mathematical techniques developed in the theory of angular momentum in
quantum mechanics are effectively applied.

The fundamental formalism of iPT can be recycled without essential
modification, and the only difference is that the biased field is colored by
tensor values. Because the original formalism of iPT is supposed to be applied
to scalar fields, the rotational symmetry in this original version of the
theory is relatively trivial. When the theory is generalized to include
tensor-valued fields, the statistical properties of biasing are largely
constrained by symmetries. In particular, the renormalized bias functions,
which are ones of the important ingredients in the iPT, can be represented by
a set of invariant functions under rotations of the coordinates system.
Explicit constructions of these invariants are derived and presented in this
paper. These invariants can also be seen as coefficients of angular expansions
by polypolar spherical harmonics.

The iPT offers a systematic way of calculating the propagators of nonlinear
perturbation theory, both in real space and in redshift space. The methodology
of deriving propagators of tensor-valued fields is generally explained in
detail and several examples of lower-order propagators are explicitly derived
in this paper. For an illustrative purpose, these examples of propagators are
applied to simple cases of predicting the correlation statistics, including
the linear power spectra and correlation functions, the lowest-order power
spectrum with primordial non-Gaussianity, the bispectrum in the tree-level
approximation, for generally tensor fields. So as to check consistencies in
the examples, we confirm that each derived formula reproduces the known result
for scalar fields, just substituting $l=m=0$ in the general formulas.

In the last section, the concept of semi-local models of bias is introduced.
The formalism of iPT does not have to assume any model of bias, and all the
uncertainties of bias is included in the renormalized bias functions, which
are determined by a nonlocal functional of the biased field in terms of the
linear density field, and have infinite degrees of freedom in general. The
semi-local models of bias are defined so that the bias functional should be a
function of a finite number of variables derived from the linear density field
in Lagrangian space. As a result, the biasing is modeled by a finite number of
parameters in this category of models. Almost all the existing models of bias
introduced in the literature fall into this category, including the halo
model, the peak theory of bias, excursion set peaks, and so forth. In this
paper, only the formal definition of the semi-local models of bias is
described, and rotationally invariant parameters in the models are identified.
Their applications to individual models of concrete targets will be addressed
in future work.

In this paper, the power spectrum and bispectrum are only evaluated in a
lowest-order, or tree-level approximation. Some further techniques are
required in calculating higher-order loop corrections, and they are described
in a subsequent paper of the series, Paper II [43]. In realistic observations,
we can only observe projected components of tensors onto the two-dimensional
sky. Predictions of the statistics of three-dimensional tensors given in this
paper can be transformed to those of projected tensors. Explicit formulations
of this transformation are given in another subsequent paper of the series,
Paper III [44]. We hope that the present formalism offers a way to carve out
the future of applying the cosmological perturbation theory to the analysis of
observations in the era of precision cosmology.

###### Acknowledgements.

I thank Y. Urakawa, K. Kogai, K. Akitsu, and A. Taruya for useful discussions.
This work was supported by JSPS KAKENHI Grants No. JP19K03835 and No.
21H03403.

##  Appendix A  Spherical tensor basis

The spherical tensors $\mathsf{Y}^{(m)}_{i_{1}\cdots i_{l}}$ are the basis to
expand symmetric tensors in the Cartesian coordinates system to the
irreducible tensors on the spherical basis. In this Appendix, the expressions
of the spherical tensor basis for ranks $l=0,1,2,3,4$ are explicitly given in
terms of the spherical bases $\mathbf{e}^{0}$ and $\mathbf{e}^{\pm}$, whose
components are denoted by ${\mathrm{e}^{m}}_{i}=[\mathbf{e}^{m}]_{i}$ with
$m=0,\pm$. The Cartesian indices of the spherical tensor basis are totally
symmetric, $\mathsf{Y}^{(m)}_{i_{1}\cdots
i_{l}}=\mathsf{Y}^{(m)}_{(i_{1}\cdots i_{l})}$, where the round brackets
indicate the symmetrization with respect to the indices inside the brackets,
i.e.,

$A_{(i_{1}\cdots
i_{n})}\equiv\frac{1}{n!}\sum_{\sigma\in\mathcal{S}_{n}}A_{i_{\sigma(1)}\cdots
i_{\sigma(n)}},$ (366)

and $\mathcal{S}_{n}$ is the symmetric group of order $n$. The following
results are straightforwardly derived from the explicit form of spherical
harmonics $Y_{lm}(\theta,\phi)$, which uniquely determine the form of
spherical basis tensors $\mathsf{Y}^{(m)}_{i_{1}\cdots i_{l}}$ to satisfy Eq.
(16), i.e.,

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

The procedure to derive the results are given as follows: First, we represent
the complex conjugate of the spherical harmonics $Y_{lm}^{*}(\theta,\phi)$ in
terms of the components of
$\bm{n}=(\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta)=(n_{1},n_{2},n_{3})$,
and the resulting expressions only contain factors of $n_{1}\mp in_{2}$ and
$n_{3}$, because of the structure of spherical harmonics. Second, we
substitute $n_{1}\mp
in_{2}\rightarrow\mp\sqrt{2}\,{\mathbf{e}^{\pm}}_{i}n_{i}$ and
$n_{3}\rightarrow{\mathbf{e}^{0}}_{i}n_{i}$ [c.f., Eq. (2)] in the expression.
Comparing the results with Eq. (367), the explicit forms of spherical basis
tensors are uniquely and obviously determined. We give the explicit
expressions up to $l\leq 4$ below for the reader’s convenience.

Spherical basis tensors of rank 0 and 1:

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

Spherical basis tensors of rank 2:

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

Spherical tensor basis of rank 3:

$\displaystyle\mathsf{Y}^{(0)}_{ijk}$
$\displaystyle=\sqrt{\frac{5}{2}}\left[{\mathrm{e}^{0}}_{i}\,{\mathrm{e}^{0}}_{j}\,{\mathrm{e}^{0}}_{k}-\frac{3}{5}\delta_{(ij}{\mathrm{e}^{0}}_{k)}\right],$
$\displaystyle\mathsf{Y}^{(\pm 1)}_{ijk}$
$\displaystyle=\frac{\sqrt{15}}{2}\left[{\mathrm{e}^{0}}_{(i}\,{\mathrm{e}^{0}}_{j}\,{\mathrm{e}^{\pm}}_{k)}-\frac{1}{5}\delta_{(ij}{\mathrm{e}^{\pm}}_{k)}\right],$
(370) $\displaystyle\mathsf{Y}^{(\pm 2)}_{ijk}$
$\displaystyle=\sqrt{3}\,{\mathrm{e}^{0}}_{(i}\,{\mathrm{e}^{\pm}}_{j}\,{\mathrm{e}^{\pm}}_{k)},$
$\displaystyle\mathsf{Y}^{(\pm 3)}_{ijk}$
$\displaystyle={\mathrm{e}^{\pm}}_{i}\,{\mathrm{e}^{\pm}}_{j}\,{\mathrm{e}^{\pm}}_{k}.$
(371)

Spherical tensor basis of rank 4:

$\displaystyle\mathsf{Y}^{(0)}_{ijkl}$
$\displaystyle=\frac{1}{2}\sqrt{\frac{35}{2}}\left[{\mathrm{e}^{0}}_{i}\,{\mathrm{e}^{0}}_{j}\,{\mathrm{e}^{0}}_{k}\,{\mathrm{e}^{0}}_{l}-\frac{6}{7}\delta_{(ij}{\mathrm{e}^{0}}_{k}{\mathrm{e}^{0}}_{l)}+\frac{3}{35}\delta_{(ij}\delta_{kl)}\right],$
$\displaystyle\mathsf{Y}^{(\pm 1)}_{ijkl}$
$\displaystyle=\sqrt{7}\left[{\mathrm{e}^{0}}_{(i}\,{\mathrm{e}^{0}}_{j}\,{\mathrm{e}^{0}}_{k}\,{\mathrm{e}^{\pm}}_{l)}-\frac{3}{7}\delta_{(ij}{\mathrm{e}^{0}}_{k}\,{\mathrm{e}^{\pm}}_{l)}\right],$
(372) $\displaystyle\mathsf{Y}^{(\pm 2)}_{ijkl}$
$\displaystyle=\sqrt{7}\left[{\mathrm{e}^{0}}_{(i}\,{\mathrm{e}^{0}}_{j}\,{\mathrm{e}^{\pm}}_{k}\,{\mathrm{e}^{\pm}}_{l)}-\frac{1}{7}\delta_{(ij}{\mathrm{e}^{\pm}}_{k}\,{\mathrm{e}^{\pm}}_{l)}\right],$
$\displaystyle\mathsf{Y}^{(\pm 3)}_{ijkl}$
$\displaystyle=2\,{\mathrm{e}^{0}}_{(i}\,{\mathrm{e}^{\pm}}_{j}\,{\mathrm{e}^{\pm}}_{k}\,{\mathrm{e}^{\pm}}_{l)},$
(373) $\displaystyle\mathsf{Y}^{(\pm 4)}_{ijkl}$
$\displaystyle={\mathrm{e}^{\pm}}_{i}\,{\mathrm{e}^{\pm}}_{j}\,{\mathrm{e}^{\pm}}_{k}\,{\mathrm{e}^{\pm}}_{l}.$
(374)

##  Appendix B  Formulas for spherical harmonics and polypolar spherical
harmonics

In this Appendix, various formulas regarding the spherical harmonics are
summarized, which are repeatedly used in the main text. As we mainly use
Racah’s normalization of the spherical harmonics in this paper, the
corresponding formulas look differently from those found in standard
textbooks. Accordingly bipolar, tripolar, and generally polypolar spherical
harmonics are also defined in convenient normalizations used in this paper.
Most of the formulas in this Appendix are derived from those found in standard
literatures [45, 46].

While the spherical harmonics are frequently represented by
$Y_{l}^{m}(\theta,\phi)$ in the literature, this function transforms as a
covariant tensor with respect to the coordinate rotation in the spherical
basis introduced in Sec. II, thus we instead prefer the notation
$Y_{lm}(\theta,\phi)$ with the lower position of the azimuthal index $m$. For
simplicity, the spherical harmonics are frequently represented by
$Y_{lm}(\bm{n})=Y_{lm}(\theta,\phi)$, where the angular coordinates are
represented by a unit vector $\bm{n}$ whose spherical coordinates are
$(\theta,\phi)$. We also use the simplified notation that the integration over
the angular coordinates $\int\sin\theta\,d\theta\,d\phi\cdots$ is represented
by $\int d^{2}n\cdots$. The standard normalization of the spherical harmonics
$Y_{lm}(\bm{n})$ is given by Eqs. (14) and (15). However, as described above,
we intensively use the Racah’s normalization of the harmonics, defined by Eq.
(17), i.e.,

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

where $P_{l}^{m}(x)$ is the associated Legendre polynomial.

The parity symmetry of the spherical harmonics is given by

$C_{lm}(-\bm{n})=(-1)^{l}C_{lm}(\bm{n}).$ (376)

The complex conjugate is given by

$C_{lm}^{*}(\bm{n})=g_{(l)}^{mm^{\prime}}C_{lm^{\prime}}(\bm{n})=(-1)^{m}C_{l,-m}(\bm{n}),$
(377)

where $g^{mm^{\prime}}_{(l)}=(-1)^{m}\delta_{m,-m^{\prime}}$ is the metric
tensor in the spherical basis and the summation over $m^{\prime}$ are
implicitly assumed when the same azimuthal index appears as a pair of lower
and upper indices, just as in the Einstein summation convention. An addition
theorem is given by

$g_{(l)}^{mm^{\prime}}C_{lm}(\bm{n}_{1})C_{lm^{\prime}}(\bm{n}_{2})=\mathit{P}_{l}(\bm{n}_{1}\cdot\bm{n}_{2}),$
(378)

where $P_{l}(x)$ is the Legendre polynomial of order $l$. The Rayleigh
expansion or plane-wave expansion is given by

$e^{\pm i\bm{k}\cdot\bm{x}}=\sum_{l=0}^{\infty}(\pm
i)^{l}(2l+1)j_{l}(kx)g_{(l)}^{mm^{\prime}}C_{lm}(\hat{\bm{k}})C_{lm^{\prime}}(\hat{\bm{x}}),$
(379)

where $k=|\bm{k}|$, $\hat{\bm{k}}=\bm{k}/k$, $r=|\bm{x}|$,
$\hat{\bm{x}}=\bm{x}/r$. The orthonormality relation is given by

$\int\frac{d^{2}n}{4\pi}\,C_{lm}(\bm{n})C_{l^{\prime}m^{\prime}}(\bm{n})=\frac{\delta_{ll^{\prime}}}{2l+1}g^{(l)}_{mm^{\prime}},$
(380)

where $g^{(l)}_{mm^{\prime}}=(-1)^{m}\delta_{m,-m^{\prime}}$ is the inverse of
$g_{(l)}^{mm^{\prime}}$, which matrix elements are the same as
$g_{(l)}^{mm^{\prime}}$. Combining Eqs. (379) and (380), we have

$\int\frac{d^{2}\hat{k}}{4\pi}\,e^{\pm
i\bm{k}\cdot\bm{x}}C_{lm}(\hat{\bm{k}})=(\pm
i)^{l}j_{l}(kx)C_{lm}(\hat{\bm{x}}).$ (381)

A product of two spherical harmonics with the same direction is given by

$C_{l_{1}m_{1}}(\bm{n})C_{l_{2}m_{2}}(\bm{n})=\sum_{l}(2l+1)\begin{pmatrix}l_{1}&l_{2}&l\\\
0&0&0\end{pmatrix}\begin{pmatrix}l_{1}&l_{2}&l\\\
m_{1}&m_{2}&m\end{pmatrix}C_{lm}^{*}(\bm{n}),$ (382)

where the last two factors but one are Wigner’s $3j$-symbols.

The bipolar spherical harmonics in our notation are defined by

$X^{l_{1}l_{2}}_{lm}(\bm{n}_{1},\bm{n}_{2})=\left(l\,l_{1}\,l_{2}\right)_{m}^{\phantom{m}m_{1}m_{2}}C_{l_{1}m_{1}}(\bm{n}_{1})C_{l_{2}m_{2}}(\bm{n}_{2}),$
(383)

where

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

is a $3j$-symbol with two azimuthal indices raised by the metric tensor of the
spherical basis. The relation to the bipolar spherical harmonics introduced in
standard textbooks [46], $\\{Y_{l_{1}}(\bm{n}_{1})\otimes
Y_{l_{2}}(\bm{n}_{2})\\}_{lm}$, is given by

$\left\\{Y_{l_{1}}(\bm{n}_{1})\otimes
Y_{l_{2}}(\bm{n}_{2})\right\\}_{lm}=\frac{(-1)^{l}}{4\pi}\sqrt{(2l+1)(2l_{1}+1)(2l_{2}+1)}\,X^{l_{1}l_{2}}_{lm}(\bm{n}_{1},\bm{n}_{2}).$
(385)

The parity symmetry of the bipolar spherical harmonics is given by

$X^{l_{1}l_{2}}_{lm}(-\bm{n}_{1},-\bm{n}_{2})=(-1)^{l_{1}+l_{2}}X^{l_{1}l_{2}}_{lm}(\bm{n}_{1},\bm{n}_{2}).$
(386)

The complex conjugate of the bipolar spherical harmonics is given by

$X^{l_{1}l_{2}*}_{lm}(\bm{n}_{1},\bm{n}_{2})=(-1)^{l+l_{1}+l_{2}}g_{(l)}^{mm^{\prime}}X^{l_{1}l_{2}}_{lm^{\prime}}(\bm{n}_{1},\bm{n}_{2}).$
(387)

The orthonormality relation for the bipolar spherical harmonics is given by

$\int\frac{d^{2}n_{1}}{4\pi}\frac{d^{2}n_{2}}{4\pi}X^{l_{1}l_{2}}_{lm}(\bm{n}_{1},\bm{n}_{2})X^{l_{1}^{\prime}l_{2}^{\prime}}_{l^{\prime}m^{\prime}}(\bm{n}_{1},\bm{n}_{2})=\frac{(-1)^{l+l_{1}+l_{2}}\delta_{ll^{\prime}}\delta_{l_{1}l_{1}^{\prime}}\delta_{l_{2}l_{2}^{\prime}}\delta^{\triangle}_{l_{1}l_{2}l}}{(2l+1)(2l_{1}+1)(2l_{2}+1)}g^{(l)}_{mm^{\prime}},$
(388)

where

$\delta^{\triangle}_{l\,l_{1}l_{2}}=\begin{cases}1,&\mathrm{if\ \
}|l_{1}-l_{2}|\leq l\leq l_{1}+l_{2},\\\ 0,&\mathrm{otherwise},\end{cases}$
(389)

i.e., $\delta^{\triangle}_{l_{1}l_{2}l_{3}}$ is unity when the set of three
numbers $(l_{1},l_{2},l_{3})$ satisfies the triangular condition, and is zero
otherwise. A product of two bipolar spherical harmonics with the same set of
directions is given by

$X^{l_{1}l_{2}}_{lm}(\bm{n}_{1},\bm{n}_{2})X^{l_{1}^{\prime}l_{2}^{\prime}}_{l^{\prime}m^{\prime}}(\bm{n}_{1},\bm{n}_{2})=\sum_{l^{\prime\prime}}(-1)^{l^{\prime\prime}}(2l^{\prime\prime}+1)\left(l\,l^{\prime}\,l^{\prime\prime}\right)_{mm^{\prime}}^{\phantom{mm^{\prime}}m^{\prime\prime}}\\\
\times\sum_{l_{1}^{\prime\prime},l_{2}^{\prime\prime}}(-1)^{l_{1}^{\prime\prime}+l_{2}^{\prime\prime}}(2l_{1}^{\prime\prime}+1)(2l_{2}^{\prime\prime}+1)\begin{pmatrix}l_{1}&l_{1}^{\prime}&l_{1}^{\prime\prime}\\\
0&0&0\end{pmatrix}\begin{pmatrix}l_{2}&l_{2}^{\prime}&l_{2}^{\prime\prime}\\\
0&0&0\end{pmatrix}\begin{Bmatrix}l&l^{\prime}&l^{\prime\prime}\\\
l_{1}&l_{1}^{\prime}&l_{1}^{\prime\prime}\\\
l_{2}&l_{2}^{\prime}&l_{2}^{\prime\prime}\end{Bmatrix}X^{l_{1}^{\prime\prime}l_{2}^{\prime\prime}}_{l^{\prime\prime}m^{\prime\prime}}(\bm{n}_{1},\bm{n}_{2}),$
(390)

where the last factor but one is the Wigner’s $9j$-symbol.

The tripolar spherical harmonics are defined by

$X^{l_{1}l_{2}l_{3}}_{L;lm}(\bm{n}_{1},\bm{n}_{2},\bm{n}_{3})=(-1)^{L}\sqrt{2L+1}\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}}C_{l_{1}m_{1}}(\bm{n}_{1})C_{l_{2}m_{2}}(\bm{n}_{2})C_{l_{3}m_{3}}(\bm{n}_{3}).$
(391)

The relation to the tripolar spherical harmonics introduced in standard
textbooks [46] is given by

$\left\\{Y_{l_{1}}(\bm{n}_{1})\otimes\left\\{Y_{l_{2}}(\bm{n}_{2})\otimes
Y_{l_{3}}(\bm{n}_{3})\right\\}_{L}\right\\}_{lm}=\frac{(-1)^{l}}{(4\pi)^{3/2}}\sqrt{(2l+1)(2l_{1}+1)(2l_{2}+1)(2l_{3}+1)}\,X^{l_{1}l_{2}l_{3}}_{L;lm}(\bm{n}_{1},\bm{n}_{2},\bm{n}_{3}).$
(392)

The parity symmetry of the tripolar spherical harmonics is given by

$X^{l_{1}l_{2}l_{3}}_{lm}(-\bm{n}_{1},-\bm{n}_{2},-\bm{n}_{3})=(-1)^{l_{1}+l_{2}+l_{3}}X^{l_{1}l_{2}l_{3}}_{lm}(\bm{n}_{1},\bm{n}_{2},\bm{n}_{3}).$
(393)

The complex conjugate of the bipolar spherical harmonics is given by

$X^{l_{1}l_{2}l_{3}*}_{lm}(\bm{n}_{1},\bm{n}_{2},\bm{n}_{3})=(-1)^{l+l_{1}+l_{2}+l_{3}}g_{(l)}^{mm^{\prime}}X^{l_{1}l_{2}l_{3}}_{lm^{\prime}}(\bm{n}_{1},\bm{n}_{2},\bm{n}_{3}).$
(394)

The orthonormality relation for the tripolar spherical harmonics is given by

$\int\frac{d^{2}n_{1}}{4\pi}\frac{d^{2}n_{2}}{4\pi}\frac{d^{2}n_{3}}{4\pi}X^{l_{1}l_{2}l_{3}}_{L;lm}(\bm{n}_{1},\bm{n}_{2},\bm{n}_{3})X^{l_{1}^{\prime}l_{2}^{\prime}l_{3}^{\prime}}_{L^{\prime};l^{\prime}m^{\prime}}(\bm{n}_{1},\bm{n}_{2},\bm{n}_{3})=\frac{(-1)^{l+l_{1}+l_{2}+l_{3}}\delta_{ll^{\prime}}\delta_{l_{1}l_{1}^{\prime}}\delta_{l_{2}l_{2}^{\prime}}\delta_{l_{3}l_{3}^{\prime}}\delta_{LL^{\prime}}\delta^{\triangle}_{l\,l_{1}L}\delta^{\triangle}_{Ll_{2}l_{3}}}{(2l+1)(2l_{1}+1)(2l_{2}+1)(2l_{3}+1)}g^{(l)}_{mm^{\prime}}.$
(395)

A product of two tripolar spherical harmonics with the same set of directions
is given by

$X^{l_{1}l_{2}l_{3}}_{L;lm}(\bm{n}_{1},\bm{n}_{2},\bm{n}_{3})X^{l_{1}^{\prime}l_{2}^{\prime}l_{3}^{\prime}}_{L^{\prime};l^{\prime}m^{\prime}}(\bm{n}_{1},\bm{n}_{2},\bm{n}_{3})=\sqrt{(2L+1)(2L^{\prime}+1)}\sum_{l^{\prime\prime}}(-1)^{l^{\prime\prime}}(2l^{\prime\prime}+1)\left(l\,l^{\prime}\,l^{\prime\prime}\right)_{mm^{\prime}}^{\phantom{mm^{\prime}}m^{\prime\prime}}\\\
\times\sum_{l_{1}^{\prime\prime},l_{2}^{\prime\prime},l_{3}^{\prime\prime}}(-1)^{l_{1}^{\prime\prime}+l_{2}^{\prime\prime}+l_{3}^{\prime\prime}}(2l_{1}^{\prime\prime}+1)(2l_{2}^{\prime\prime}+1)(2l_{3}^{\prime\prime}+1)\begin{pmatrix}l_{1}&l_{1}^{\prime}&l_{1}^{\prime\prime}\\\
0&0&0\end{pmatrix}\begin{pmatrix}l_{2}&l_{2}^{\prime}&l_{2}^{\prime\prime}\\\
0&0&0\end{pmatrix}\begin{pmatrix}l_{3}&l_{3}^{\prime}&l_{3}^{\prime\prime}\\\
0&0&0\end{pmatrix}\\\
\times\sum_{L^{\prime\prime}}\sqrt{2L^{\prime\prime}+1}\begin{Bmatrix}l&l^{\prime}&l^{\prime\prime}\\\
l_{1}&l_{1}^{\prime}&l_{1}^{\prime\prime}\\\
L&L^{\prime}&L^{\prime\prime}\end{Bmatrix}\begin{Bmatrix}L&L^{\prime}&L^{\prime\prime}\\\
l_{2}&l_{2}^{\prime}&l_{2}^{\prime\prime}\\\
l_{3}&l_{3}^{\prime}&l_{3}^{\prime\prime}\end{Bmatrix}X^{l_{1}^{\prime\prime}l_{2}^{\prime\prime}l_{3}^{\prime\prime}}_{L^{\prime\prime};l^{\prime\prime}m^{\prime\prime}}(\bm{n}_{1},\bm{n}_{2},\bm{n}_{3}).$
(396)

We generally define the polypolar spherical harmonics of arbitrary order by

$X^{l_{1}l_{2}\cdots l_{n}}_{L_{2}L_{3}\cdots
L_{n-1};lm}(\bm{n}_{1},\bm{n}_{2},\ldots,\bm{n}_{n})=(-1)^{L_{2}+\cdots+L_{n-1}}\sqrt{(2L_{2}+1)(2L_{3}+1)\cdots(2L_{n-1}+1)}\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}}\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}}(\bm{n}_{1})C_{l_{2}m_{2}}(\bm{n}_{2})\cdots
C_{l_{n}m_{n}}(\bm{n}_{n}).$ (397)

They satisfy recursive relations,

$X^{l_{1}l_{2}\cdots l_{n}}_{L_{2}L_{3}\cdots
L_{n-1};lm}(\bm{n}_{1},\bm{n}_{2},\ldots,\bm{n}_{n})=(-1)^{L_{2}}\sqrt{(2L_{2}+1)}\left(l\,l_{1}\,L_{2}\right)_{m}^{\phantom{m}m_{1}M_{2}}C_{l_{1}m_{1}}(\bm{n}_{1})X^{l_{2}\cdots
l_{n}}_{L_{3}\cdots L_{n-1}}(\bm{n}_{2},\bm{n}_{3},\ldots,\bm{n}_{n}).$ (398)

The parity symmetry of the tripolar spherical harmonics is given by

$X^{l_{1}l_{2}\cdots
l_{n}}_{lm}(-\bm{n}_{1},-\bm{n}_{2},\ldots,-\bm{n}_{n})=(-1)^{l_{1}+l_{2}+\cdots+l_{n}}X^{l_{1}l_{2}\cdots
l_{n}}_{lm}(\bm{n}_{1},\bm{n}_{2},\ldots,\bm{n}_{n}).$ (399)

The complex conjugate of the bipolar spherical harmonics is given by

$X^{l_{1}l_{2}\cdots
l_{n}*}_{lm}(\bm{n}_{1},\bm{n}_{2},\ldots,\bm{n}_{n})=(-1)^{l+l_{1}+l_{2}+\cdots+l_{n}}g_{(l)}^{mm^{\prime}}X^{l_{1}l_{2}\cdots
l_{n}}_{lm^{\prime}}(\bm{n}_{1},\bm{n}_{2},\ldots,\bm{n}_{n}).$ (400)

The orthonormality relation for the polypolar spherical harmonics in general
is straightforwardly derived from the orthonormality of the $3j$-symbols, and
the result is given by

$\int\frac{d^{2}n_{1}}{4\pi}\cdots\frac{d^{2}n_{n}}{4\pi}\,X^{l_{1}l_{2}\cdots
l_{n}}_{L_{2}L_{3}\cdots
L_{n-1};lm}(\bm{n}_{1},\bm{n}_{2},\ldots,\bm{n}_{n})X^{l_{1}^{\prime}l_{2}^{\prime}\cdots
l_{n}^{\prime}}_{L_{2}^{\prime}L_{3}^{\prime}\cdots
L_{n-1}^{\prime};l^{\prime}m^{\prime}}(\bm{n}_{1},\bm{n}_{2},\ldots,\bm{n}_{n})\\\
=\frac{(-1)^{l+l_{1}+\cdots+l_{n}}\delta_{ll^{\prime}}\delta_{l_{1}l_{1}^{\prime}}\cdots\delta_{l_{n}l_{n}^{\prime}}\delta_{L_{2}L_{2}^{\prime}}\cdots\delta_{L_{n-1}L_{n-1}^{\prime}}\delta^{\triangle}_{l\,l_{1}L_{2}}\delta^{\triangle}_{L_{2}l_{2}L_{3}}\cdots\delta^{\triangle}_{L_{n-2}l_{n-2}L_{n-1}}\delta^{\triangle}_{L_{n-1}l_{n-1}l_{n}}}{(2l+1)(2l_{1}+1)(2l_{2}+1)\cdots(2l_{n}+1)}g^{(l)}_{mm^{\prime}}.$
(401)

A product of polypolar harmonics with the same set of directions is given by

$X^{l_{1}l_{2}\cdots l_{n}}_{L_{2}L_{3}\cdots
L_{n-1};lm}(\bm{n}_{1},\bm{n}_{2},\ldots,\bm{n}_{n})X^{l_{1}^{\prime}l_{2}^{\prime}\cdots
l_{n}^{\prime}}_{L_{2}^{\prime}L_{3}^{\prime}\cdots
L_{n-1}^{\prime};l^{\prime}m^{\prime}}(\bm{n}_{1},\bm{n}_{2},\ldots,\bm{n}_{n})=\sqrt{(2L_{2}+1)\cdots(2L_{n-1}+1)(2L_{2}^{\prime}+1)\cdots(2L_{n-1}^{\prime}+1)}\\\
\times\sum_{l^{\prime\prime}}(-1)^{l^{\prime\prime}}(2l^{\prime\prime}+1)\left(l\,l^{\prime}\,l^{\prime\prime}\right)_{mm^{\prime}}^{\phantom{mm^{\prime}}m^{\prime\prime}}\sum_{l_{1}^{\prime\prime},\ldots,l_{n}^{\prime\prime}}(-1)^{l_{1}^{\prime\prime}+\cdots+l_{n}^{\prime\prime}}(2l_{1}^{\prime\prime}+1)\cdots(2l_{n}^{\prime\prime}+1)\begin{pmatrix}l_{1}&l_{1}^{\prime}&l_{1}^{\prime\prime}\\\
0&0&0\end{pmatrix}\cdots\begin{pmatrix}l_{n}&l_{n}^{\prime}&l_{n}^{\prime\prime}\\\
0&0&0\end{pmatrix}\\\
\times\sum_{L_{2}^{\prime\prime},\ldots,L_{n-1}^{\prime\prime}}\sqrt{(2L_{2}^{\prime\prime}+1)\cdots(2L_{n-1}^{\prime\prime}+1)}\begin{Bmatrix}l&l^{\prime}&l^{\prime\prime}\\\
l_{1}&l_{1}^{\prime}&l_{1}^{\prime\prime}\\\
L_{2}&L_{2}^{\prime}&L_{2}^{\prime\prime}\end{Bmatrix}\begin{Bmatrix}L_{2}&L_{2}^{\prime}&L_{2}^{\prime\prime}\\\
l_{2}&l_{2}^{\prime}&l_{2}^{\prime\prime}\\\
L_{3}&L_{3}^{\prime}&L_{3}^{\prime\prime}\end{Bmatrix}\cdots\begin{Bmatrix}L_{n-2}&L_{n-2}^{\prime}&L_{n-2}^{\prime\prime}\\\
l_{n-2}&l_{n-2}^{\prime}&l_{n-2}^{\prime\prime}\\\
L_{n-1}&L_{n-1}^{\prime}&L_{n-1}^{\prime\prime}\end{Bmatrix}\begin{Bmatrix}L_{n-1}&L_{n-1}^{\prime}&L_{n-1}^{\prime\prime}\\\
l_{n-1}&l_{n-1}^{\prime}&l_{n-1}^{\prime\prime}\\\
l_{n}&l_{n}^{\prime}&l_{n}^{\prime\prime}\end{Bmatrix}\\\ \times
X^{l_{1}^{\prime\prime}l_{2}^{\prime\prime}\cdots
l_{n}^{\prime\prime}}_{L_{2}^{\prime\prime}L_{3}^{\prime\prime}\cdots
L_{n-1}^{\prime\prime};l^{\prime\prime}m^{\prime\prime}}(\bm{n}_{1},\bm{n}_{2},\ldots,\bm{n}_{n}).$
(402)

The last result can be derived by applying the sum rules of Eqs. (427) and
(428) in Appendix C.

##  Appendix C  The $\bm{3j}$-, $\bm{6j}$\- and $\bm{9j}$-symbols

Wigner’s $3j$-, $6j$\- and $9j$-symbols are the coefficients in adding and
recoupling of angular momenta in quantum mechanics. These coefficients are
frequently used in the main text of this paper, because they are useful as
well in our formulation, regarding rotational symmetry of the Universe in a
statistical sense. We summarize some of the formulas in this Appendix for
readers’ convenience. All the formulas below can be found in, or derived from,
Refs. [45, 46]. In the following, we assume $l,l_{1},l_{2},\ldots$ are non-
negative integers (excluding the possibility of half-integers) as well as
azimuthal indices, and we assume $(-1)^{l}=(-1)^{-l}$, $(-1)^{m}=(-1)^{-m}$ in
the following.

In this paper, we employ a simplified notation for the $3j$-symbol,

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

where the $3j$-symbol is nonzero only when $m_{1}+m_{2}+m_{3}=0$ in the above,
and azimuthal indices $m$’s are raised or lowered by spherical metric tensors,
$g^{(l)}_{mm^{\prime}}=(-1)^{m}\delta_{m,-m^{\prime}}=g_{(l)}^{mm^{\prime}}$,
as exemplified in Eq. (384). In particular, the symmetries of the $3j$-symbol
indicates

$\left(l_{1}\,l_{2}\,l_{3}\right)^{m_{1}m_{2}m_{3}}=\begin{pmatrix}l_{1}&l_{2}&l_{3}\\\
-m_{1}&-m_{2}&-m_{3}\end{pmatrix}=(-1)^{l_{1}+l_{2}+l_{3}}\begin{pmatrix}l_{1}&l_{2}&l_{3}\\\
m_{1}&m_{2}&m_{3}\end{pmatrix}=(-1)^{l_{1}+l_{2}+l_{3}}\left(l_{1}\,l_{2}\,l_{3}\right)_{m_{1}m_{2}m_{3}},$
(404)

since $m_{1}+m_{2}+m_{3}=0$. Even permutations of the $3j$-symbols have the
same value,

$\left(l_{1}\,l_{2}\,l_{3}\right)_{m_{1}m_{2}m_{3}}=\left(l_{2}\,l_{3}\,l_{1}\right)_{m_{2}m_{3}m_{1}}=\left(l_{3}\,l_{1}\,l_{2}\right)_{m_{3}m_{1}m_{2}},$
(405)

while odd permutations are given by

$\left(l_{2}\,l_{1}\,l_{3}\right)_{m_{2}m_{1}m_{3}}=\left(l_{1}\,l_{3}\,l_{2}\right)_{m_{1}m_{3}m_{2}}=\left(l_{3}\,l_{2}\,l_{1}\right)_{m_{3}m_{2}m_{1}}=(-1)^{l_{1}+l_{2}+l_{3}}\left(l_{1}\,l_{2}\,l_{3}\right)_{m_{1}m_{2}m_{3}}.$
(406)

The above equalities hold even when some of the azimuthal indices are raised
on both sides of the equations.

The orthogonality relations of the $3j$-symbol are represented by

$\displaystyle\left(l\,l_{1}\,l_{2}\right)_{mm_{1}m_{2}}\left(l^{\prime}\,l_{1}\,l_{2}\right)_{m^{\prime}}^{\phantom{m^{\prime}}m_{1}m_{2}}=\frac{(-1)^{l+l_{1}+l_{2}}\delta_{ll^{\prime}}}{2l+1}g^{(l)}_{mm^{\prime}}\delta^{\triangle}_{l\,l_{1}l_{2}},$
(407)
$\displaystyle\sum_{l}(-1)^{l}(2l+1)\left(l\,l_{1}\,l_{2}\right)_{mm_{1}m_{2}}\left(l\,l_{1}\,l_{2}\right)^{m}_{\phantom{m}m_{1}^{\prime}m_{2}^{\prime}}=(-1)^{l_{1}+l_{2}}g^{(l_{1})}_{m_{1}m_{1}^{\prime}}g^{(l_{2})}_{m_{2}m_{2}^{\prime}},$
(408)

where $\delta^{\triangle}_{l_{1}l_{2}l_{3}}=1$ only when $(l_{1},l_{2},l_{3})$
satisfies the triangle inequality, and is zero otherwise. From Eq. (407), we
have a derived formula,

$\left(l_{1}\,l_{2}\,l_{3}\right)_{m_{1}m_{2}m_{3}}\left(l_{1}\,l_{2}\,l_{3}\right)^{m_{1}m_{2}m_{3}}=(-1)^{l_{1}+l_{2}+l_{3}}\delta^{\triangle}_{l_{1}l_{2}l_{3}}.$
(409)

When one or two of the $l$’s in the $3j$-symbols are zero, we have

$\left(l\,0\,0\right)_{m00}=\delta_{l0}\delta_{m0},\qquad\left(l\,l^{\prime}\,0\right)_{mm^{\prime}0}=\frac{(-1)^{l}\delta_{ll^{\prime}}}{\sqrt{2l+1}}g^{(l)}_{mm^{\prime}}.$
(410)

In this paper, the Gaunt integral with our notation is denoted by

$\left[l_{1}\,l_{2}\,l_{3}\right]_{m_{1}m_{2}m_{3}}\equiv\int\frac{d^{2}n}{4\pi}\,C_{l_{1}m_{1}}(\hat{\bm{n}})C_{l_{2}m_{2}}(\hat{\bm{n}})C_{l_{3}m_{3}}(\hat{\bm{n}})=\begin{pmatrix}l_{1}&l_{2}&l_{3}\\\
0&0&0\end{pmatrix}\begin{pmatrix}l_{1}&l_{2}&l_{3}\\\
m_{1}&m_{2}&m_{3}\end{pmatrix},$ (411)

and azimuthal indices are raised and lowered by the spherical metric
$g_{(l)}^{mm^{\prime}}$, $g^{(l)}_{mm^{\prime}}$ as in the case of
$3j$-symbol. A $3j$-symbol with vanishing $m$’s are nonzero only when the sum
of $l$’s are even, i.e.,

$\begin{pmatrix}l_{1}&l_{2}&l_{3}\\\ 0&0&0\end{pmatrix}=0,\qquad\mathrm{when}\
\ l_{1}+l_{2}+l_{3}=\mathrm{odd}.$ (412)

Correspondingly, the Gaunt integrals are nonzero only when
$l_{1}+l_{2}+l_{3}=\mathrm{even}$ in the above equation. The covariant version
of the Gaunt integral has the same components as the contravariant one:

$\left[l_{1}\,l_{2}\,l_{3}\right]_{m_{1}m_{2}m_{3}}=\left[l_{1}\,l_{2}\,l_{3}\right]^{m_{1}m_{2}m_{3}},$
(413)

because of a symmetry of the $3j$-symbols, Eq. (404). The Gaunt integrals are
totally symmetric under the permutation,

$\left[l_{1}\,l_{2}\,l_{3}\right]_{m_{1}m_{2}m_{3}}=\left[l_{2}\,l_{3}\,l_{1}\right]_{m_{2}m_{3}m_{1}}=\left[l_{3}\,l_{1}\,l_{2}\right]_{m_{3}m_{1}m_{2}}=\left[l_{2}\,l_{1}\,l_{3}\right]_{m_{2}m_{1}m_{3}}=\left[l_{1}\,l_{3}\,l_{2}\right]_{m_{1}m_{3}m_{2}}=\left[l_{3}\,l_{2}\,l_{1}\right]_{m_{3}m_{2}m_{1}}.$
(414)

An orthogonality relation of the $3j$-symbols, Eq. (407), indicates the
orthogonality relation of the Gaunt integral:

$\left[l\,l_{1}\,l_{2}\right]_{mm_{1}m_{2}}\left[l^{\prime}\,l_{1}\,l_{2}\right]_{m^{\prime}}^{\phantom{m^{\prime}}m_{1}m_{2}}=\frac{1}{2l+1}\begin{pmatrix}l&l_{1}&l_{2}\\\
0&0&0\end{pmatrix}^{2}\delta_{ll^{\prime}}g^{(l)}_{mm^{\prime}}.$ (415)

When one or two of the $l$’s in the Gaunt integral are zero, we have

$\left[l\,0\,0\right]_{m00}=\delta_{l0}\delta_{m0},\qquad\left[l\,l^{\prime}\,0\right]_{mm^{\prime}0}=\frac{\delta_{ll^{\prime}}}{2l+1}g^{(l)}_{mm^{\prime}}.$
(416)

The Wigner’s $6j$-symbol is associated with the recoupling of three angular
momenta, and defined by

$\begin{Bmatrix}l_{1}&l_{2}&l_{3}\\\
l_{4}&l_{5}&l_{6}\end{Bmatrix}=(-1)^{l_{1}+l_{2}+l_{3}+l_{4}+l_{5}+l_{6}}\left(l_{1}\,l_{2}\,l_{3}\right)^{m_{1}m_{2}m_{3}}\left(l_{1}\,l_{5}\,l_{6}\right)_{m_{1}}^{\phantom{m_{1}}m_{5}m_{6}}\left(l_{4}\,l_{2}\,l_{6}\right)^{m_{4}}_{\phantom{m_{4}}m_{2}m_{6}}\left(l_{4}\,l_{5}\,l_{3}\right)_{m_{4}m_{5}m_{3}}.$
(417)

or equivalently,

$\begin{Bmatrix}l_{1}&l_{2}&l_{3}\\\
l_{4}&l_{5}&l_{6}\end{Bmatrix}\left(l_{1}\,l_{2}\,l_{3}\right)_{m_{1}m_{2}m_{3}}=(-1)^{l_{1}+l_{2}+l_{3}}\left(l_{1}\,l_{5}\,l_{6}\right)_{m_{1}m_{5}m_{6}}\left(l_{4}\,l_{2}\,l_{6}\right)_{m_{4}m_{2}}^{\phantom{m_{4}m_{2}}m_{6}}\left(l_{4}\,l_{5}\,l_{3}\right)^{m_{4}m_{5}}_{\phantom{m_{4}m_{5}}m_{3}}.$
(418)

The $6j$-symbols are invariant against any permutation of the columns:

$\begin{Bmatrix}l_{1}&l_{2}&l_{3}\\\
l_{4}&l_{5}&l_{6}\end{Bmatrix}=\begin{Bmatrix}l_{2}&l_{3}&l_{1}\\\
l_{5}&l_{6}&l_{4}\end{Bmatrix}=\begin{Bmatrix}l_{3}&l_{1}&l_{2}\\\
l_{6}&l_{4}&l_{5}\end{Bmatrix}=\begin{Bmatrix}l_{2}&l_{1}&l_{3}\\\
l_{5}&l_{4}&l_{6}\end{Bmatrix}=\begin{Bmatrix}l_{1}&l_{3}&l_{2}\\\
l_{4}&l_{6}&l_{5}\end{Bmatrix}=\begin{Bmatrix}l_{3}&l_{2}&l_{1}\\\
l_{6}&l_{5}&l_{4}\end{Bmatrix},$ (419)

and also invariant against the interchange of the upper and lower arguments in
each of any two columns,

$\begin{Bmatrix}l_{1}&l_{2}&l_{3}\\\
l_{4}&l_{5}&l_{6}\end{Bmatrix}=\begin{Bmatrix}l_{1}&l_{5}&l_{6}\\\
l_{4}&l_{2}&l_{3}\end{Bmatrix}=\begin{Bmatrix}l_{4}&l_{2}&l_{6}\\\
l_{1}&l_{5}&l_{3}\end{Bmatrix}=\begin{Bmatrix}l_{4}&l_{5}&l_{3}\\\
l_{1}&l_{2}&l_{6}\end{Bmatrix}.$ (420)

The Wigner’s $9j$-symbol is associated with the recoupling of four angular
momenta, and defined by

$\begin{Bmatrix}l_{1}&l_{2}&l_{3}\\\ l_{4}&l_{5}&l_{6}\\\
l_{7}&l_{8}&l_{9}\end{Bmatrix}=\left(l_{1}\,l_{2}\,l_{3}\right)^{m_{1}m_{2}m_{3}}\left(l_{4}\,l_{5}\,l_{6}\right)^{m_{4}m_{5}m_{6}}\left(l_{7}\,l_{8}\,l_{9}\right)^{m_{7}m_{8}m_{9}}\left(l_{1}\,l_{4}\,l_{7}\right)_{m_{1}m_{4}m_{7}}\left(l_{2}\,l_{5}\,l_{8}\right)_{m_{2}m_{5}m_{8}}\left(l_{3}\,l_{6}\,l_{9}\right)_{m_{3}m_{6}m_{9}}.$
(421)

The $9j$-symbols are invariant under transposition,

$\begin{Bmatrix}l_{1}&l_{2}&l_{3}\\\ l_{4}&l_{5}&l_{6}\\\
l_{7}&l_{8}&l_{9}\end{Bmatrix}=\begin{Bmatrix}l_{1}&l_{4}&l_{7}\\\
l_{2}&l_{5}&l_{8}\\\ l_{3}&l_{6}&l_{9}\end{Bmatrix}.$ (422)

and also invariant against cyclic (even) permutations of the columns or lows:

$\begin{Bmatrix}l_{1}&l_{2}&l_{3}\\\ l_{4}&l_{5}&l_{6}\\\
l_{7}&l_{8}&l_{9}\end{Bmatrix}=\begin{Bmatrix}l_{2}&l_{3}&l_{1}\\\
l_{5}&l_{6}&l_{4}\\\
l_{8}&l_{9}&l_{7}\end{Bmatrix}=\begin{Bmatrix}l_{3}&l_{1}&l_{2}\\\
l_{6}&l_{4}&l_{5}\\\
l_{9}&l_{7}&l_{8}\end{Bmatrix}=\begin{Bmatrix}l_{7}&l_{8}&l_{9}\\\
l_{1}&l_{2}&l_{3}\\\
l_{4}&l_{5}&l_{6}\end{Bmatrix}=\begin{Bmatrix}l_{4}&l_{5}&l_{6}\\\
l_{7}&l_{8}&l_{9}\\\ l_{1}&l_{2}&l_{3}\end{Bmatrix}.$ (423)

They change the sign by non-cyclic (odd) permutations of the columns or lows,
if the sum of all the arguments is odd:

$(-1)^{R}\begin{Bmatrix}l_{1}&l_{2}&l_{3}\\\ l_{4}&l_{5}&l_{6}\\\
l_{7}&l_{8}&l_{9}\end{Bmatrix}=\begin{Bmatrix}l_{2}&l_{1}&l_{3}\\\
l_{5}&l_{4}&l_{6}\\\
l_{8}&l_{7}&l_{9}\end{Bmatrix}=\begin{Bmatrix}l_{3}&l_{2}&l_{1}\\\
l_{6}&l_{5}&l_{4}\\\
l_{9}&l_{8}&l_{7}\end{Bmatrix}=\begin{Bmatrix}l_{1}&l_{3}&l_{2}\\\
l_{4}&l_{6}&l_{5}\\\
l_{7}&l_{9}&l_{8}\end{Bmatrix}=\begin{Bmatrix}l_{4}&l_{5}&l_{6}\\\
l_{1}&l_{2}&l_{3}\\\
l_{7}&l_{8}&l_{9}\end{Bmatrix}=\begin{Bmatrix}l_{7}&l_{8}&l_{9}\\\
l_{4}&l_{5}&l_{6}\\\
l_{1}&l_{2}&l_{3}\end{Bmatrix}=\begin{Bmatrix}l_{1}&l_{2}&l_{3}\\\
l_{7}&l_{8}&l_{9}\\\ l_{4}&l_{5}&l_{6}\end{Bmatrix},$ (424)

where

$R=l_{1}+l_{2}+l_{3}+l_{4}+l_{5}+l_{6}+l_{7}+l_{8}+l_{9}.$ (425)

Many formulas among $3j$-, $6j$-, and $9j$-symbols, including orthogonality
relations and sum rules, can be found in the standard literature mentioned
above. In particular, we use the following sum rules in this paper for
products of $3j$-symbols:

$\left(l_{1}\,l_{2}\,L\right)_{m_{1}m_{2}}^{\phantom{m_{1}m_{2}}M}\left(L\,l_{3}\,l_{4}\right)_{Mm_{3}m_{4}}=(-1)^{l_{2}+l_{3}}\sum_{L^{\prime}}(2L^{\prime}+1)\left(l_{1}\,l_{3}\,L^{\prime}\right)_{m_{1}m_{3}}^{\phantom{m_{1}m_{3}}M^{\prime}}\left(L^{\prime}\,l_{2}\,l_{4}\right)_{M^{\prime}m_{2}m_{4}}\begin{Bmatrix}l_{1}&l_{2}&L\\\
l_{4}&l_{3}&L^{\prime}\end{Bmatrix},$ (426)
$\left(L\,l_{1}\,l_{2}\right)_{M}^{\phantom{M}m_{1}m_{2}}\left(L^{\prime}\,l_{1}^{\prime}\,l_{2}^{\prime}\right)_{M^{\prime}}^{\phantom{M^{\prime}}m_{1}^{\prime}m_{2}^{\prime}}\left(l_{2}^{\prime}\,l_{2}\,l_{2}^{\prime\prime}\right)_{m_{2}^{\prime}m_{2}}^{\phantom{m_{2}^{\prime}m_{2}}m_{2}^{\prime\prime}}=\sum_{l_{1}^{\prime\prime},L^{\prime\prime}}(2l_{1}^{\prime\prime}+1)(2L^{\prime\prime}+1)\left(L\,L^{\prime}\,L^{\prime\prime}\right)_{MM^{\prime}}^{\phantom{MM^{\prime}}M^{\prime\prime}}\left(L^{\prime\prime}\,l_{2}^{\prime\prime}\,l_{1}^{\prime\prime}\right)_{M^{\prime\prime}}^{\phantom{m^{\prime\prime}}m_{2}^{\prime\prime}m_{1}^{\prime\prime}}\left(l_{1}\,l_{1}^{\prime}\,l_{1}^{\prime\prime}\right)^{m_{1}m_{1}^{\prime}}_{\phantom{m_{1}m_{1}^{\prime}}m_{1}^{\prime\prime}}\begin{Bmatrix}l_{1}&l_{1}^{\prime}&l_{1}^{\prime\prime}\\\
l_{2}&l_{2}^{\prime}&l_{2}^{\prime\prime}\\\
L&L^{\prime}&L^{\prime\prime}\end{Bmatrix},$ (427)

and

$\left(L\,l_{1}\,l_{2}\right)_{M}^{\phantom{M}m_{1}m_{2}}\left(L^{\prime}\,l_{1}^{\prime}\,l_{2}^{\prime}\right)_{M^{\prime}}^{\phantom{M^{\prime}}m_{1}^{\prime}m_{2}^{\prime}}\left(l_{1}^{\prime}\,l_{1}\,l_{1}^{\prime\prime}\right)_{m_{1}^{\prime}m_{1}}^{\phantom{m_{1}^{\prime}m_{1}}m_{1}^{\prime\prime}}\left(l_{2}^{\prime}\,l_{2}\,l_{2}^{\prime\prime}\right)_{m_{2}^{\prime}m_{2}}^{\phantom{m_{2}^{\prime}m_{2}}m_{2}^{\prime\prime}}=\sum_{L^{\prime\prime}}(2L^{\prime\prime}+1)\left(L\,L^{\prime}\,L^{\prime\prime}\right)_{MM^{\prime}}^{\phantom{MM^{\prime}}M^{\prime\prime}}\left(L^{\prime\prime}\,l_{2}^{\prime\prime}\,l_{1}^{\prime\prime}\right)_{M^{\prime\prime}}^{\phantom{m^{\prime\prime}}m_{2}^{\prime\prime}m_{1}^{\prime\prime}}\begin{Bmatrix}l_{1}&l_{1}^{\prime}&l_{1}^{\prime\prime}\\\
l_{2}&l_{2}^{\prime}&l_{2}^{\prime\prime}\\\
L&L^{\prime}&L^{\prime\prime}\end{Bmatrix}.$ (428)

## References

  * [1] D. H. Lyth and A. R. Liddle, “The primordial density perturbation: Cosmology, inflation and the origin of structure,”
  * [2] D. J. Eisenstein et al. [SDSS], Astrophys. J. 633, 560-574 (2005)
  * [3] N. Dalal, O. Dore, D. Huterer and A. Shirokov, Phys. Rev. D 77, 123514 (2008)
  * [4] V. Desjacques, D. Jeong and F. Schmidt, Phys. Rept. 733, 1-193 (2018)
  * [5] C. Hunter, Astrophys. J. 139, 570 (1964)
  * [6] K. Tomita, Progr. Theor. Phys. 37, 831 (1967)
  * [7] P. J. E. Peebles, Large-Scale Structure of the Universe (Princeton University Press, Princeton, 1980)
  * [8] F. Bernardeau, S. Colombi, E. Gaztanaga and R. Scoccimarro, Phys. Rept. 367, 1-248 (2002)
  * [9] H. J. Mo and S. D. M. White, Mon. Not. Roy. Astron. Soc. 282, 347 (1996)
  * [10] H. J. Mo, Y. P. Jing and S. D. M. White, Mon. Not. Roy. Astron. Soc. 284, 189 (1997)
  * [11] A. Cooray and R. K. Sheth, Phys. Rept. 372, 1-129 (2002)
  * [12] W. H. Press and P. Schechter, Astrophys. J. 187, 425-438 (1974)
  * [13] J. R. Bond, S. Cole, G. Efstathiou and N. Kaiser, Astrophys. J. 379, 440 (1991)
  * [14] R. K. Sheth and G. Tormen, Mon. Not. Roy. Astron. Soc. 308, 119 (1999)
  * [15] N. Kaiser, Astrophys. J. Lett. 284, L9-L12 (1984)
  * [16] J. M. Bardeen, J. R. Bond, N. Kaiser and A. S. Szalay, Astrophys. J. 304, 15-61 (1986)
  * [17] T. Matsubara, Phys. Rev. D 83, 083518 (2011)
  * [18] T. Matsubara, Phys. Rev. D 90, 043537 (2014)
  * [19] T. Matsubara, Phys. Rev. D 86, 063518 (2012)
  * [20] T. Matsubara and V. Desjacques, Phys. Rev. D 93, no.12, 123522 (2016)
  * [21] T. Buchert, Astron. Astrophys. 223, 9-24 (1989)
  * [22] F. Moutarde, J. M. Alimi, F. R. Bouchet, R. Pellat and A. Ramani, Astrophys. J. 382, 377 (1991)
  * [23] T. Matsubara, Astrophys. J. Suppl. 101, 1-28 (1995)
  * [24] M. Bartelmann and P. Schneider, Phys. Rept. 340, 291-472 (2001)
  * [25] P. Catelan, M. Kamionkowski and R. D. Blandford, Mon. Not. Roy. Astron. Soc. 320, L7-L13 (2001)
  * [26] T. Okumura and Y. P. Jing, Astrophys. J. Lett. 694, L83-L86 (2009)
  * [27] B. Joachimi, M. Cacciato, T. D. Kitching, A. Leonard, R. Mandelbaum, B. M. Schäfer, C. Sifón, H. Hoekstra, A. Kiessling and D. Kirk, et al. Space Sci. Rev. 193, no.1-4, 1-65 (2015)
  * [28] K. Kogai, T. Matsubara, A. J. Nishizawa and Y. Urakawa, JCAP 08, 014 (2018)
  * [29] T. Okumura and A. Taruya, Mon. Not. Roy. Astron. Soc. 493, no.1, L124-L128 (2020)
  * [30] T. Okumura, A. Taruya and T. Nishimichi, Phys. Rev. D 100, no.10, 103507 (2019)
  * [31] T. Okumura, A. Taruya and T. Nishimichi, Mon. Not. Roy. Astron. Soc. 494, no.1, 694-702 (2020)
  * [32] A. Taruya and T. Okumura, Astrophys. J. Lett. 891, L42 (2020)
  * [33] T. Okumura and A. Taruya, Phys. Rev. D 106, no.4, 043523 (2022)
  * [34] N. E. Chisari and C. Dvorkin, JCAP 12, 029 (2013)
  * [35] N. Arkani-Hamed and J. Maldacena, [arXiv:1503.08043 [hep-th]].
  * [36] J. Blazek, Z. Vlah and U. Seljak, JCAP 08, 015 (2015)
  * [37] J. Blazek, N. MacCrann, M. A. Troxel and X. Fang, Phys. Rev. D 100, no.10, 103506 (2019)
  * [38] D. M. Schmitz, C. M. Hirata, J. Blazek and E. Krause, JCAP 07, 030 (2018)
  * [39] Z. Vlah, N. E. Chisari and F. Schmidt, JCAP 01, 025 (2020)
  * [40] Z. Vlah, N. E. Chisari and F. Schmidt, JCAP 05, 061 (2021)
  * [41] A. Taruya and K. Akitsu, JCAP 11, no.11, 061 (2021)
  * [42] K. Kogai, K. Akitsu, F. Schmidt and Y. Urakawa, JCAP 03, 060 (2021)
  * [43] T. Matsubara, [arXiv:2210.11085] (Paper II)
  * [44] T. Matsubara, [arXiv:2304.13304] (Paper III)
  * [45] A. R. Edmonds, “Angular momentum in quantum mechanics,” (Princeton University Press, Princeton, NJ, 1996)
  * [46] D. A. Varshalovich, A. N. Moskalev and V. K. Khersonskii, “Quantum Theory Of Angular Momentum,” (World Scientific Publishing, Singapore, 1988)
  * [47] J. J. Sakurai and J. Napolitano, “Modern Quantum Mechanics,” (Cambridge University Press, Cambridge, 2020).
  * [48] V. Desjacques, M. Crocce, R. Scoccimarro and R. K. Sheth, Phys. Rev. D 82, 103529 (2010)
  * [49] E. P. Wigner, On the Matrices Which Reduce the Kronecker Products of Representations of S. R. Groups In: A. S. Wightman (eds), The Collected Works of Eugene Paul Wigner, vol A / 1. (Springer, Berlin, Heidelberg, 1993)
  * [50] A. J. M. Spencer, Int. J. Eng. Sci., 8(6) 475-481 (1970)
  * [51] J. Applequist, J. Phys. A:Math. Gen. 22, 4303 (1989)
  * [52] T. Matsubara, Phys. Rev. D 92, no.2, 023534 (2015)
  * [53] F. Bernardeau, M. Crocce and R. Scoccimarro, Phys. Rev. D 78, 103521 (2008)
  * [54] F. Bernardeau, M. Crocce and E. Sefusatti, Phys. Rev. D 82, 083507 (2010)
  * [55] V. Zheligovsky and U. Frisch, J. Fluid Mech. 749, 404 (2014)
  * [56] P. Catelan, Mon. Not. Roy. Astron. Soc. 276, 115 (1995)
  * [57] P. Catelan and T. Theuns, Mon. Not. Roy. Astron. Soc. 282, 455 (1996)
  * [58] A. Aviles and J. L. Cervantes-Cota, Phys. Rev. D 96, no.12, 123526 (2017)
  * [59] T. Matsubara, Phys. Rev. D 77, 063530 (2008)
  * [60] P. J. E. Peebles, Astrophys. J. , 155, 393 (1969)
  * [61] W. Hu, Phys. Rev. D 64, 083005 (2001)
  * [62] N. Kaiser, Mon. Not. Roy. Astron. Soc. 227, 1-27 (1987)
  * [63] A. J. S. Hamilton, Astrophys. J. Lett. 385, L5-L8 (1992)
  * [64] F. Schmidt, N. E. Chisari and C. Dvorkin, JCAP 10, 032 (2015)
  * [65] K. Akitsu, T. Kurita, T. Nishimichi, M. Takada and S. Tanaka, Phys. Rev. D 103, no.8, 083508 (2021)
  * [66] M. Shiraishi, E. Komatsu, M. Peloso and N. Barnaby, JCAP 05, 002 (2013)
  * [67] S. Matarrese, L. Verde and A. F. Heavens, Mon. Not. Roy. Astron. Soc. 290, 651-662 (1997)
  * [68] P. McDonald and A. Roy, JCAP 08, 020 (2009)
  * [69] R. K. Sheth, K. C. Chan and R. Scoccimarro, Phys. Rev. D 87, no.8, 083002 (2013)
  * [70] T. Baldauf, U. Seljak, V. Desjacques and P. McDonald, Phys. Rev. D 86, 083540 (2012)
  * [71] M. Mirbabayi, F. Schmidt and M. Zaldarriaga, JCAP 07, 030 (2015)
  * [72] T. Matsubara, Phys. Rev. D 101, no.4, 043532 (2020)
  * [73] M. Schmittfull, Z. Vlah and P. McDonald, Phys. Rev. D 93, no.10, 103528 (2016)
  * [74] M. Schmittfull and Z. Vlah, Phys. Rev. D 94, no.10, 103530 (2016)
  * [75] J. E. McEwen, X. Fang, C. M. Hirata and J. A. Blazek, JCAP 09, 015 (2016)
  * [76] X. Fang, J. A. Blazek, J. E. McEwen and C. M. Hirata, JCAP 02, 030 (2017)