# On a Generating Function for the Isotropic Basis Functions and Other
Connected Results

Zachary Slepian  Department of Astronomy, University of Florida, 211 Bryant
Space Science Center, Gainesville, FL 32611, USA Jessica Chellino  Jiamin Hou
Department of Astronomy, University of Florida, 211 Bryant Space Science
Center, Gainesville, FL 32611, USA

###### Abstract

Recently isotropic basis functions of $N$ unit vector arguments were
presented; these are of significant use in measuring the N-Point Correlation
Functions (NPCFs) of galaxy clustering. Here we develop the generating
function for these basis functions—i.e. that function which, expanded in a
power series, has as its angular part the isotropic functions. We show that
this can be developed using basic properties of the plane wave. A main use of
the generating function is as an efficient route to obtaining the Cartesian
basis expressions for the isotropic functions. We show that the methods here
enable computing difficult overlap integrals of multiple spherical Bessel
functions, and we also give related expansions of the Dirac Delta function
into the isotropic basis. Finally, we outline how the Cartesian expressions
for the isotropic basis functions might be used to enable a faster NPCF
algorithm on the CPU.

##  1 Introduction

Many situations in physics invoke rotational invariance; particularly, in
cosmology, the distribution of galaxies around a given center should appear
isotropic, modulo distortions induced by conversion of redshift into position.
Hence, functions that are invariant under simultaneous rotation of all of
their unit-vector arguments are particularly useful. In [1], these isotropic
basis functions were developed, and their importance and other related
previous work is more fully outlined there. The idea of using them to
parameterize N-Point Correlation Functions, quantifying the spatial clustering
of the distribution of galaxies, was presented in [2], with a detection of the
even-parity, “connected” (due to nonlinear gravitational evolution) 4-Point
Correlation Function (4PCF) in [3]. The idea of using the decomposition of the
4PCF into even and odd-parity 3-argument basis functions to search for parity
violation first presented in [4], and the covariance matrix required for this
search computed in the isotropic basis in [5]. The first detection of parity
violation in 3D large-scale structure was made in [6], with $7\sigma$ evidence
in one sample and $3\sigma$ evidence in another smaller sample, both of the
Sloan Digital Sky Survey Baryon Oscillation Spectroscopic Survey (SDSS BOSS).
Another work [7] also found evidence for parity violation in CMASS, using the
same covariance matrix template.111This work should unfortunately not be
viewed as an independent confirmatory study, since it was done based on ideas
its author took from work with our group, and the covariance matrix template
was computed by a member of our group and not intended for publication. Thus
this template was not fully calibrated to the final BOSS data used. This work
also uses incorrect systematics weights on the galaxies. It selected a CMASS-
like sample from the BOSS combined sample (LOWZ+CMASS) by imposing a redshift
cut, and thus ends up with the combined sample weights. However, those weights
correct for systematics in both CMASS and LOWZ, meaning that they will impose
some degree of essentially random fluctuations on the sample used (CMASS-like)
from any systematic corrections required by LOWZ. This latter work also used
the isotropic basis functions to perform calculations of an inflationary
model, and subsequent works have also done so, e.g. [8]. The isotropic basis
functions have also been used to compute a model of the even-parity, connected
4PCF (i.e. that due to gravity, beyond the Gaussian Random Field contribution)
in [9].

Here, our goal is to provide generating functions for these isotropic basis
functions, as well as present a few useful associated results and
applications.

First, we briefly recall what a generating function is. A frequent example in
physics is the multipole expansion, wherein one has

$\displaystyle\frac{1}{|\boldsymbol{r}_{1}-\boldsymbol{r}_{2}|}=\sum_{\ell=0}^{\infty}\frac{r_{<}^{\ell}}{r_{>}^{\ell+1}}\mathcal{L}_{\ell}(\hat{r}_{1}\cdot\hat{r}_{2}).$
(1)

Here, we have an expansion of the potential between two points at
$\boldsymbol{r}_{1}$ and $\boldsymbol{r}_{2}$, and the left-hand side is the
“generating function” of the Legendre polynomials as the $\ell^{\rm th}$
Legendre polynomial, $\mathcal{L}_{\ell}$, can be picked out if we focus on
the right-hand side term proportional to the $\ell^{th}$ power of $r_{<}$. We
may pull out $r_{<}$ from both sides and work in terms of
$t\equiv(r_{<}/r_{>})\leq 1$, such that the series’ convergence is manifest
and by expanding the vector sum magnitude on the left-hand side we may write

$\displaystyle\frac{1}{\sqrt{x^{2}-2xt+1}}=\sum_{\ell=0}^{\infty}t^{\ell}\mathcal{L}_{\ell}(x),$
(2)

where we have defined $x\equiv\hat{r}_{1}\cdot\hat{r}_{2}$. The left-hand side
above is the most common form for the “generating function” of the Legendre
polynomials, which themselves form an isotropic basis set for functions of two
unit vectors (i.e. $\hat{r}_{1}$ and $\hat{r}_{2}$; since $x$ is the dot
product it is invariant under simultaneous rotation of $\hat{r}_{1}$ and
$\hat{r}_{2}$). Since the isotropic basis functions of [1] are essentially
generalizations of the Legendre polynomials to more than two unit vectors, it
seems desirable to have a generating function for them as well. Motivated by
work proving the addition theorem for regular solid harmonics using the plane
wave expansion [10, 11], we do so here.

This work is structured as follows. First, we present the basic expansion we
will use and show how it works straightforwardly for $N=2$ and $3$ basis
functions (§2). We then discuss $N=4$ basis functions, where one needs to have
intermediate angular momenta to fully fix the coupling structure (§3). We next
move to $N=5$ functions, which have the same challenge (§4). We then (§5)
outline how to obtain generating functions for the odd-parity isotropic basis
functions, which is a more complicated matter and requires additional work
beyond what is needed for the even-parity case. We then turn to expanding the
Dirac Delta function in this basis (§6), and subsequently present connections
of the generating function to spherical Bessel function overlap integrals
(§7). In §8 we suggest one possible application of the Cartesian forms of the
isotropic functions, outlining a new algorithm to compute the N-Point
Correlation Function (NPCF) of a given density field; this approach may offer
a speed-up over the spherical-harmonic-based encore code [2] when run on a
CPU. §9 concludes.

##  2 Basic Expansion and Generalized Gaunt Integrals

Consider a vector

$\displaystyle\boldsymbol{R}\equiv\sum_{i=1}^{N}\boldsymbol{r}_{i},$ (3)

and a momentum $\boldsymbol{p}$ dual to $\boldsymbol{R}$. We note that, for
what follows, it is not essential that $\boldsymbol{p}$ be a momentum; it can
be simply considered as a dummy parameter if one prefers. Let us now examine
two different ways of writing the plane wave expansion into spherical
harmonics and spherical Bessel functions. First, we have

$\displaystyle e^{i\boldsymbol{p}\cdot\boldsymbol{R}}=4\pi\sum_{\ell
m}i^{\ell}j_{\ell}(pR)Y_{\ell m}^{*}(\hat{p})Y_{\ell m}(\hat{R}).$ (4)

We note that here and in all that follows, if the limits of an orbital
angular-momentum-like index (e.g. $\ell$ above) are not indicated explicitly,
they are from zero to infinity, and if the limits of a spin-like index (e.g.
$m$ above) are not indicated they are from $-\ell$ to $+\ell$, inclusive.
Throughout, we will refer to $m$-like indices as “spins”, though in detail
they are simply the $z$-component of the orbital angular momentum and need not
have the physical meaning of spin.

Now let us write equation (4) out by factoring the plane wave as
$\exp[i\boldsymbol{p}\cdot\boldsymbol{r}_{1}]\times\cdots\times\exp[i\boldsymbol{p}\cdot\boldsymbol{r}_{N}]$.
We then have

$\displaystyle
e^{i\boldsymbol{p}\cdot\boldsymbol{R}}=(4\pi)^{N}\prod_{{i=1}}^{{N}}\left\\{\sum_{\ell_{i}}i^{\ell_{i}}j_{\ell_{i}}(pr_{i})\sum_{m_{i}}Y_{\ell_{i}m_{i}}^{*}(\hat{p})Y_{\ell_{i}m_{i}}(\hat{r}_{i})\right\\}.$
(5)

We now integrate each of the above equations against $Y_{\ell m}(\hat{p})$,
and set them equal. We have

$\displaystyle\boxed{4\pi i^{\ell}j_{\ell}(pR)Y_{\ell
m}(\hat{R})=(4\pi)^{N}\sum_{\ell_{i=1,\ldots,N}}i^{\sum_{{i=1}}^{{N}}\ell_{i}}\prod_{{i=1}}^{{N}}\left\\{j_{\ell_{i}}(pr_{i})\sum_{m_{i}}Y_{\ell_{i}m_{i}}(\hat{r}_{i})\right\\}\int
d\Omega_{p}\;Y_{\ell
m}(\hat{p})\prod_{i=1}^{N}Y_{\ell_{i}m_{i}}^{*}(\hat{p}).}$ (6)

We note that each $\ell_{i}$ assumes a full range of values, over which we are
summing. We now convert $Y^{{*}}_{\ell
m}(\hat{p})=(-1)^{m}Y_{\ell-m}(\hat{p})$ so that our $d\Omega_{p}$ integral is
over all conjugated spherical harmonics, in which case it is then equal to
what it would be had they been all unconjugated (this equality follows from
the fact that the result must be real).222$d\Omega_{p}$ is the integration
measure for integrating over all directions of $\hat{p}$. We then may rewrite
this integral as

$\displaystyle\int d\Omega_{p}\;Y_{\ell
m}(\hat{p})\prod_{i=1}^{N}Y_{\ell_{i}m_{i}}^{*}(\hat{p})\equiv(-1)^{m}\mathcal{G}_{\ell_{1}\cdots\ell_{N}\ell}^{m_{1}\cdots
m_{N}-m},$ (7)

where $\mathcal{G}$ is the generalized Gaunt integral over $(N+1)$ all
unconjugated spherical harmonics with the indicated angular momenta and spins.
The last spin in the superscript is $-m$, not a subtraction of $m$ from
$m_{N}$. The generalized Gaunt can be considered to be defined by this
relation.

###  2.1 Generalized Gaunt Integrals

We will explicitly write down this integral for $N=2,3,4$ and $5$. We note the
$N$ here is the number of arguments of the basis functions we are studying,
i.e. the number of $\hat{r}_{i}$; the number of harmonics involved in the
integral (7) is then $N+1$, as one may see because there is one additional
harmonic beyond the product indexed by $i$.

####  2.1.1 $N=2$ Argument Basis Functions $\rightarrow$ 3 Harmonics

For $N=2$, we thus have three spherical harmonics, so:

$\displaystyle\mathcal{G}_{\ell_{1}\ell_{2}\ell}^{m_{1}m_{2}-m}=\mathcal{C}_{\ell_{1}\ell_{2}\ell}\begin{pmatrix}\ell_{1}&\ell_{2}&\ell\\\
0&0&0\end{pmatrix}\begin{pmatrix}\ell_{1}&\ell_{2}&\ell\\\
m_{1}&m_{2}&-m\end{pmatrix},$ (8)

with

$\displaystyle\mathcal{C}_{\ell_{1}\ell_{2}\ell}\equiv\sqrt{\frac{(2\ell_{1}+1)(2\ell_{2}+1)(2\ell+1)}{4\pi}}.$
(9)

####  2.1.2 $N=3$ Argument Basis Functions $\rightarrow$ 4 Harmonics

For $N=3$ we have, choosing the same coupling approach (“pyramidal”) as in
[1], and projecting the product of
$Y_{\ell_{1}m_{1}}(\hat{p})Y_{\ell_{2}m_{2}}(\hat{p})$ onto a sum over
$Y_{\ell_{12}m_{12}}(\hat{p})$ (where this last is unconjugated) using [12]
34.3.20, we find

$\displaystyle\mathcal{G}_{\ell_{1}\ell_{2}\ell_{3}\ell}^{m_{1}m_{2}m_{3}-m}=\sum_{\ell_{12}}(-1)^{m_{12}}\mathcal{G}_{\ell_{1}\ell_{2}\ell_{12}}^{m_{1}m_{2}-m_{12}}\mathcal{G}_{\ell_{12}\ell_{3}\ell}^{m_{12}m_{3}-m}.$
(10)

We do not need to sum over $m_{12}$ on the right-hand side because
$m_{1}+m_{2}-m_{12}=0$ fixes $m_{12}$. This constraint stems from the fact
that the usual Gaunt integrals contain a 3-$j$ symbol whose bottom row is
these spins.

####  2.1.3 $N=4$ Argument Basis Functions $\rightarrow$ 5 Harmonics

For $N=4$, we have

$\displaystyle\mathcal{G}_{\ell_{1}\ell_{2}\ell_{3}\ell_{4}\ell}^{m_{1}m_{2}m_{3}m_{4}-m}=\sum_{\ell_{12},\ell_{123}}(-1)^{m_{12}}\mathcal{G}_{\ell_{1}\ell_{2}\ell_{12}}^{m_{1}m_{2}-m_{12}}(-1)^{m_{123}}\mathcal{G}_{\ell_{12}\ell_{3}\ell_{123}}^{m_{12}m_{3}-m_{123}}\mathcal{G}_{\ell_{123}\ell_{4}\ell}^{m_{123}m_{4}-m}.$
(11)

Again we note there is no need for a sum over any multi-index spins (e.g.
$m_{12}$ or similar).

####  2.1.4 $N=5$ Argument Basis Functions $\rightarrow$ 6 Harmonics

For $N=5$ we have

$\displaystyle\mathcal{G}_{\ell_{1}\ell_{2}\ell_{3}\ell_{4}\ell_{5}\ell}^{m_{1}m_{2}m_{3}m_{4}m_{5}-m}=$
$\displaystyle\sum_{\ell_{12},\ell_{123},\ell_{1234}}\qquad(-1)^{m_{12}}\mathcal{G}_{\ell_{1}\ell_{2}\ell_{12}}^{m_{1}m_{2}-m_{12}}(-1)^{m_{123}}\mathcal{G}_{\ell_{12}\ell_{3}\ell_{123}}^{m_{12}m_{3}-m_{123}}(-1)^{m_{1234}}$
$\displaystyle\qquad\times\mathcal{G}_{\ell_{123}\ell_{4}\ell_{1234}}^{m_{123}m_{4}-m_{1234}}\mathcal{G}_{\ell_{1234}\ell_{5}\ell}^{m_{1234}m_{5}-m}.$
(12)

####  2.1.5 Use in Fundamental Relation

With these generalized Gaunt integrals in hand, we now return to our main
thread, the fundamental, boxed relation (6). With the integration over
$d\Omega_{p}$ now performed (we recall equation 7), it becomes:

$\displaystyle i^{\ell}j_{\ell}(pR)Y_{\ell
m}(\hat{R})=(4\pi)^{N-1}\sum_{\ell_{i=1,\ldots,N}}i^{\sum_{i=1}^{N}}{\ell_{i}}\prod_{{i=1}}^{{N}}\left\\{j_{\ell_{i}}(pr_{i})\sum_{m_{i}}Y_{\ell_{i}m_{i}}(\hat{r}_{i})\right\\}(-1)^{m}\mathcal{G}_{\ell_{1}\cdots\ell_{N}\ell}^{m_{1}\cdots
m_{N}-m}.$ (13)

####  2.1.6 Reduction when $\ell=0$

Let us now consider the case where $\ell\to 0$, meaning $m=0$ as well. Doing
so is motivated by the fact that we are seeking isotropic basis functions of
the collection of unit vectors $\hat{r}_{i}$; hence, the way we have set our
relations up, they must have coupled to an isotropic mode in $\hat{p}$ prior
to integration over it. We may impose this condition by taking it that
$\ell=0$ in the harmonic in $\hat{p}$ in equation (6).333We might have reached
this point more directly by simply integrating over $d\Omega_{p}$ in equation
(6), which is of course equivalent to integrating against
$1=\sqrt{4\pi}\;Y_{00}^{*}(\hat{p})$. The approach used above, however, is
more general. In particular, it can actually offer a set of “anisotropic”
basis functions as well if $\ell$ is allowed to be non-zero (e.g. as noted in
[13]; see also [14]). We obtain

$\displaystyle\boxed{\frac{1}{\sqrt{4\pi}}j_{0}(pR)=(4\pi)^{N-1}\sum_{{\ell_{i=1,\ldots,N}}}i^{\sum_{{i=1}}^{{N}}\ell_{i}}\prod_{{i=1}}^{{N}}\left\\{j_{\ell_{i}}(pr_{i})\sum_{m_{i}}Y_{\ell_{i}m_{i}}(\hat{r}_{i})\right\\}\mathcal{G}_{\ell_{1}\cdots\ell_{N}0}^{m_{1}\cdots
m_{N}0}.}$ (14)

###  2.2 Functions of Two Arguments: $N=2$

The generalized Gaunt integral becomes, for $N=2$ (recalling equation 8):

$\displaystyle\mathcal{G}_{\ell_{1}\ell_{2}0}^{m_{1}m_{2}0}=\frac{1}{\sqrt{4\pi}}\;(-1)^{m_{1}}\delta^{\rm
K}_{\ell_{1}\ell_{2}}\delta^{\rm K}_{m_{1}-m_{2}},$ (15)

obtained by applying NIST DLMF 34.3.1 twice. In particular, we used that a
3-$j$ symbol with zeros in the third column forces equality between the first
and second angular momenta, and that the associated two spins sum to zero.

Hence our $N=2$ result becomes

$\displaystyle
j_{0}(pR)=4\pi\sum_{\ell_{1}}(-1)^{\ell_{1}}j_{\ell_{1}}(pr_{1})j_{\ell_{1}}(pr_{2})\sum_{m_{1}}Y_{\ell_{1}m_{1}}(\hat{r}_{1})Y^{*}_{\ell_{1}m_{1}}(\hat{r}_{2}).$
(16)

We used the conjugate identity to absorb the spin-dependent phase into the
second harmonic above, and we note that we set $m=0$ as implied by $\ell=0$.
We note that the left-hand side depends on the angles between $\hat{r}_{1}$
and $\hat{r}_{2}$ because $R=|\bf{r}_{1}+\bf{r}_{2}|$.

Comparing with [1] equation (3), we have a fundamental relation for generating
the $N=2$ isotropic basis functions, as

$\displaystyle\boxed{j_{0}(pR)=4\pi\sum_{\ell_{1}}\sqrt{2\ell_{1}+1}\;j_{\ell_{1}}(pr_{1})j_{\ell_{1}}(pr_{2})\mathcal{P}_{\ell_{1}\ell_{1}}(\hat{r}_{1},\hat{r}_{2}).}$
(17)

We may now take the limit that $p\to 0$, and using that $j_{0}(x)=\sin x/x$,
we have the series for the left-hand side as

$\displaystyle\boxed{j_{0}(pR)\approx
1-\frac{(pR)^{2}}{6}+\frac{(pR)^{4}}{120}-\frac{(pR)^{6}}{5,040}+\cdots.}$
(18)

We box this relation because it will be of much use in what follows. On the
right-hand side, we have the series for each spherical Bessel Function (sBF)
at lowest order as $(pr_{i})^{\ell_{1}}/(2\ell_{1}+1)!!$. The double factorial
$n!!$ is the product of all the integers from 1 up to $n$ (inclusive) that
have the same parity as $n$.444E.g.
https://mathworld.wolfram.com/DoubleFactorial.html.

We may now demand term-by-term equality for the coefficients of powers of $p$
on each side of equation (17); hence the above relation shows that $j_{0}(pR)$
with $R\equiv|\boldsymbol{r}_{1}+\boldsymbol{r}_{2}|$ is a generating function
for the $N=2$ isotropic basis functions.555There is no approximation involved
in this, even though we used a series for the sBFs. It is a power-matching
argument, just like that typically employed in solving differential equations
via the method of Frobenius. Fundamentally, it is because, roughly, the powers
of $pR$ may be considered an orthogonal basis, which is, roughly, the result
of the Stone-Weierstrass theorem.

###  2.3 Examples for $N=2$

We may easily read off results for the first few
$\mathcal{P}_{\ell_{1}\ell_{1}}$. Throughout this section, we will make use of
equation (18) for the left-hand side of equation (17).

For $\ell_{1}=0$, we have after a line of algebra that
$\mathcal{P}_{00}=1/(4\pi)$. For $\ell_{1}=1$, we expand
$|\boldsymbol{r}_{1}+\boldsymbol{r}_{2}|^{2}$ on the left-hand side and retain
only the part proportional to $r_{1}r_{2}$, to match the right-hand side. We
find
$\mathcal{P}_{11}(\hat{r}_{1},\hat{r}_{2})=-\sqrt{3}/(4\pi)\;(\hat{r}_{1}\cdot\hat{r}_{2})$.

For $\ell_{1}=2$, there is a subtlety. The sub-leading expansion of the
$\ell_{1}=0$ term on the right-hand side can contribute a term in $p^{4}$ (the
relevant power for $\ell_{1}=2$) as well. This contribution also will enter
proportional to $r_{1}^{2}r_{2}^{2}$, so we cannot eliminate it by selecting
based on the $r_{i}$-dependence.666The $\ell_{1}=1$ terms can contribute as
$p^{4}$ but, given the parity of each $j_{1}$, never with dependence on
$r_{1}^{2}r_{2}^{2}$; an even-parity term in $p$ will have to have mismatched
powers of $r_{1}$ and $r_{2}$. We also note that no $\ell_{1}>2$ can
contribute such a low power, given the leading-order behavior of the SBFs at
small argument.

However, we must indeed count the contribution of the sub-leading
$\ell_{1}=0$; it must be subtracted off to isolate the contribution to the
right-hand side sum due only to $\mathcal{P}_{22}$. We have, focusing only on
terms proportional to $p^{4}$ and $r_{1}^{2}r_{2}^{2}$,

$\displaystyle\frac{R^{4}}{120}\bigg{|}_{r_{1}^{2}r_{2}^{2}}=4\pi\left[\frac{1}{36}\mathcal{P}_{00}(\hat{r}_{1},\hat{r}_{2})+\frac{\sqrt{5}}{15^{2}}\mathcal{P}_{22}(\hat{r}_{1},\hat{r}_{2})\right].$
(19)

Our notation on the left-hand side means that when $R^{4}$ is expanded only
terms in $r_{1}^{2}r_{2}^{2}$ should be retained (and we that we have divided
the $r_{1}^{2}r_{2}^{2}$ out from both sides). We find, working in terms of
$\mu_{12}\equiv\hat{r}_{1}\cdot\hat{r}_{2}$, that

$\displaystyle\frac{2+4\mu_{12}^{2}}{120}=\frac{1}{36}+4\pi\frac{\sqrt{5}}{15^{2}}\mathcal{P}_{2,2}(\hat{r}_{1},\hat{r}_{2});$
(20)

solving for $\mathcal{P}_{2,2}$ yields

$\displaystyle\mathcal{P}_{2,2}(\hat{r}_{1},\hat{r}_{2})=\frac{3\sqrt{5}}{8\pi}\left[(\hat{r}_{1}\cdot\hat{r}_{2})^{2}-\frac{1}{3}\right]$
(21)

in agreement with [1] equation A.1., last line.

We observe that due to the power-matching in the $r_{i}$, we need only to
consider one term from each series for lower $\ell$ than the $\ell_{i}$ under
consideration, and we also need only to consider expansions of lower $\ell$
sBFs of the same parity here.

Working through $\ell_{1}=3$, we find the fundamental equation

$\displaystyle-\frac{R^{6}}{5,040}\bigg{|}_{r_{1}^{3}r_{2}^{3}}=4\pi\left[\frac{\sqrt{3}}{30^{2}}\mathcal{P}_{1,1}(\hat{r}_{1},\hat{r}_{2})+\frac{\sqrt{7}}{105^{2}}\mathcal{P}_{3,3}(\hat{r}_{1},\hat{r}_{2})\right].$
(22)

Expanding the left-hand side gives

$\displaystyle-\frac{R^{6}}{5,040}\bigg{|}_{r_{1}^{3}r_{2}^{3}}=-\frac{1}{5,040}\left[8\mu_{12}^{3}+12\mu_{12}\right].$
(23)

Solving for $\mathcal{P}_{3,3}$ we obtain

$\displaystyle\mathcal{P}_{3,3}(\hat{r}_{1},\hat{r}_{2})=-\frac{\sqrt{7}}{4\pi}\mathcal{L}_{3}(\hat{r}_{1}\cdot\hat{r}_{2}),$
(24)

where $\mathcal{L}_{3}$ is the Legendre polynomial of order three. This result
is in agreement with [1] equation 3 with $\Lambda=3$ there.

###  2.4 Functions of Three Arguments: $N=3$

We now seek to extend our work to $N>2$ functions. For $N=3$, the generalized
Gaunt integral becomes

$\displaystyle\mathcal{G}_{\ell_{1}\ell_{2}\ell_{3}0}^{m_{1}m_{2}m_{3}0}=\frac{1}{\sqrt{4\pi}}\mathcal{G}_{\ell_{1}\ell_{2}\ell_{3}}^{m_{1}m_{2}m_{3}}.$
(25)

Generally, we note that setting the last angular momentum and spin to $``00"$
always introduces a factor of $1/\sqrt{4\pi}$.

Using this in our fundamental relation (14), we find after some simplification

$\displaystyle j_{0}(pR)$
$\displaystyle=(4\pi)^{2}\sum_{\ell_{1}\ell_{2}\ell_{3}}i^{\ell_{1}+\ell_{2}+\ell_{3}}j_{\ell_{1}}(pr_{1})j_{\ell_{2}}(pr_{2})j_{\ell_{3}}(pr_{3})\begin{pmatrix}\ell_{1}&\ell_{2}&\ell_{3}\\\
0&0&0\end{pmatrix}\mathcal{C}_{\ell_{1}\ell_{2}\ell_{3}}$
$\displaystyle\qquad\times\sum_{m_{1}m_{2}m_{3}}\begin{pmatrix}\ell_{1}&\ell_{2}&\ell_{3}\\\
m_{1}&m_{2}&m_{3}\end{pmatrix}Y_{\ell_{1}m_{1}}(\hat{r}_{1})Y_{\ell_{2}m_{2}}(\hat{r}_{2})Y_{\ell_{3}m_{3}}(\hat{r}_{3}).$
(26)

Using [1] equations 16 and 11 this becomes

$\displaystyle\boxed{j_{0}(pR)=(4\pi)^{2}\sum_{\ell_{1}\ell_{2}\ell_{3}}i^{\ell_{1}+\ell_{2}+\ell_{3}}j_{\ell_{1}}(pr_{1})j_{\ell_{2}}(pr_{2})j_{\ell_{3}}(pr_{3})\begin{pmatrix}\ell_{1}&\ell_{2}&\ell_{3}\\\
0&0&0\end{pmatrix}\mathcal{C}_{\ell_{1}\ell_{2}\ell_{3}}\mathcal{P}_{\ell_{1}\ell_{2}\ell_{3}}(\hat{r}_{1},\hat{r}_{2},\hat{r}_{3}).}$
(27)

We note that the right-hand side will always be real, since the zero-spin
3-$j$ symbol guarantees that $\ell_{1}+\ell_{2}+\ell_{3}$ is even. We recall
that $\mathcal{C}_{\ell_{1}\ell_{2}\ell_{3}}$ is defined by equation (9).

Now, let us consider the structure above further. A given power of $p$ on the
left-hand side (when we expand in a series in $p$) must be equal to
$\ell_{1}+\ell_{2}+\ell_{3}$. One can think of this as analogous to an overall
energy level (as in the usual model of an atom, which has overall energy
levels and then shells within them, s p d f etc.) for the set of isotropic
functions satisfying the constraint. Then, within each level, we may select
out one or another isotropic function by requiring a given behavior in
$r_{1},r_{2},r_{3}$.

This latter aspect is analogous to how the Legendre polynomials are organized
when written using their generating function as in equation (1). There, one
requires a given behavior in $r_{<}/r_{>}$ at each order. For Legendres, since
they have only two arguments, they have only one $\ell$ (as we indeed saw in
our $N=2$ isotropic functions above), so there is no distinction, and hence no
need for ordering, within a given “energy level”—there is only one “shell” in
each. Here, when we go to $N=3$, we have multiple “shells” within each “energy
level”, and we can use $p$ to indicate the energy level and then the power of
each of the $r_{i}$ to indicate the shell.

We make one final remark. $j_{0}(pR)$ is even-parity, and hence only even
powers of $p$ will enter the expansion for $p\ll 1$. Consequently, we must
have that $\ell_{1}+\ell_{2}+\ell_{3}$ is even and so we can only have even-
parity isotropic functions on the right-hand side. This is also enforced by
the even-sum rule on the zero-spin 3-$j$ symbol that appears. We turn to the
issue of odd-parity basis functions in §5.

###  2.5 Examples for $N=3$

Here we work out several examples of the three-argument basis functions using
this work’s formalism. We begin with $\mathcal{P}_{110}$. Here, we have

$\displaystyle-\frac{(pR)^{2}}{6}\bigg{|}_{r_{1}r_{2}}=(4\pi)^{2}\left[-\frac{1}{9}\left(-\frac{1}{\sqrt{3}}\right)\frac{3}{\sqrt{4\pi}}\right]\mathcal{P}_{110}(\hat{r}_{1},\hat{r}_{2},\hat{r}_{3}),$
(28)

where we used equation (18) for the left-hand side and equation (27) for the
right-hand side, evaluating the phase and 3-$j$ symbol with
$\ell_{1}=\ell_{2}=1$ and $\ell_{3}=0$, and using the definition of
$\mathcal{C}_{\ell_{1}\ell_{2}\ell_{3}}$ equation (9). We recall that our
notation $\bigg{|}_{r_{1}r_{2}}$ means that we have taken only terms involving
this set of powers of the $r_{i}$ and that we have then cancelled any powers
of $p$ and the $r_{i}$ from both sides.

Now expanding $R^{2}=|\boldsymbol{r}_{1}+\boldsymbol{r}_{2}|^{2}$ and
retaining only the terms in $r_{1}r_{2}$, a few lines of algebra shows that

$\displaystyle\mathcal{P}_{110}(\hat{r}_{1},\hat{r}_{2},\hat{r}_{3})=-\frac{1}{4\pi}\sqrt{\frac{3}{4\pi}}\;\left(\hat{r}_{1}\cdot\hat{r}_{2}\right),$
(29)

in agreement with [1] Appendix A.2., first line.

We now pursue a case requiring subtraction of sub-leading terms in the series
of lower-order sBFs on the right-hand side, as we already saw was required in
the $N=2$ functions. We consider $\mathcal{P}_{112}$. We have, using equation
(18) for the left-hand side and equation (27) for the right-hand side, that

$\displaystyle\frac{(pR)^{4}}{120}\bigg{|}_{r_{1}r_{2}r_{3}^{2}}=(4\pi)^{2}\left[-\frac{1}{54}\sqrt{\frac{3}{4\pi}}\;\mathcal{P}_{110}+\frac{1}{45}\sqrt{\frac{2}{3(4\pi)}}\;\mathcal{P}_{{112}}\right],$
(30)

where the first term on the right-hand side came from noticing that expanding
$j_{1}(pr_{1})j_{1}(pr_{2})$ to leading order each, but $j_{0}(pr_{3})$ to
second order, will also give an overall $p^{4}$. The second term is simply
from taking the lowest-order terms in the series for, respectively,
$j_{1}(pr_{1}),j_{1}(pr_{2})$ and $j_{2}(pr_{3})$.

Computing the expansion of $R^{4}$ and retaining only the term indicated by
the evaluation bar on the left-hand side, we find

$\displaystyle\frac{(pR)^{4}}{120}\bigg{|}_{r_{1}r_{2}r_{3}^{2}}=\frac{1}{30}\left[\hat{r}_{1}\cdot\hat{r}_{2}+2(\hat{r}_{2}\cdot\hat{r}_{3})(\hat{r}_{3}\cdot\hat{r}_{1})\right].$
(31)

Inserting this result in equation (30), using our result for
$\mathcal{P}_{110}$ from equation (29, and solving for $\mathcal{P}_{112}$
yields

$\displaystyle\mathcal{P}_{112}(\hat{r}_{1},\hat{r}_{2},\hat{r}_{3})=\sqrt{\frac{27}{2(4\pi)^{3}}}\left[(\hat{r}_{1}\cdot\hat{r}_{3})(\hat{r}_{2}\cdot\hat{r}_{3})-\frac{1}{3}(\hat{r}_{1}\cdot\hat{r}_{2})\right],$
(32)

in agreement with [1] Appendix A.2., third line.

##  3 Dealing with Intermediate Angular Momenta: $N=4$

We here discuss how to deal with the intermediate angular momenta required as
one goes to isotropic functions of four arguments or more. We first show that
the most obvious extension of the approach of the previous sections fails, and
then outline the correct method.

###  3.1 Failure of the Naive Approach

We now consider the generalized Gaunt for $N=4$. Again making use of NIST DLMF
34.3.1 and after some simplification, we have

$\displaystyle\mathcal{G}_{\ell_{1}\ell_{2}\ell_{3}\ell_{4}0}^{m_{1}m_{2}m_{3}m_{4}0}=\frac{1}{\sqrt{4\pi}}\sum_{\ell_{12}}(-1)^{m_{12}}\mathcal{G}_{\ell_{1}\ell_{2}\ell_{12}}^{m_{1}m_{2}-m_{12}}\mathcal{G}_{\ell_{12}\ell_{3}\ell_{4}}^{m_{12}m_{3}m_{4}}.$
(33)

Inserting this in our fundamental relation equation (14), we find

$\displaystyle
j_{0}(pR)=(4\pi)^{3}\sum_{\ell_{1},\ldots,\ell_{4}}i^{\ell_{1}+\ell_{2}+\ell_{3}+\ell_{4}}\prod_{i=1}^{4}\left\\{j_{\ell_{i}}(pr_{i})\sum_{m_{i}}Y_{\ell_{i}m_{i}}(\hat{r}_{i})\right\\}\mathcal{G}_{\ell_{1}\ell_{2}\ell_{3}\ell_{4}0}^{m_{1}m_{2}m_{3}m_{4}0}.$
(34)

Using [1] equation (12), we have

$\displaystyle j_{0}(pR)=$
$\displaystyle(4\pi)^{2}\sum_{\ell_{1},\ldots,\ell_{4}}i^{\ell_{1}+\ell_{2}+\ell_{3}+\ell_{4}}\left[\prod_{i=1}^{4}j_{\ell_{i}}(pr_{i})\right](-1)^{\ell_{1}+\ell_{2}+\ell_{3}+\ell_{4}}\mathcal{D}^{\rm
P}$
$\displaystyle\qquad\qquad\qquad\times\sum_{\ell_{12}}\sqrt{2\ell_{12}+1}\begin{pmatrix}\ell_{1}&\ell_{2}&\ell_{12}\\\
0&0&0\end{pmatrix}\begin{pmatrix}\ell_{12}&\ell_{3}&\ell_{4}\\\
0&0&0\end{pmatrix}(-1)^{\ell_{12}}$
$\displaystyle\qquad\qquad\qquad\times\mathcal{P}_{\ell_{1}\ell_{2}(\ell_{12})\ell_{3}\ell_{4}}(\hat{r}_{1},\hat{r}_{2},\hat{r}_{3},\hat{r}_{4}).$
(35)

The additional phase $(-1)^{\ell_{1}+\ell_{2}+\ell_{3}+\ell_{4}}$ comes from
the definition of $\mathcal{P}_{\Lambda}$, [1] equation 11, in particular
$\epsilon(\Lambda)$ as defined in equation 9 there;
$\Lambda\equiv\\{\ell_{1},\ell_{2},\ell_{12},\ldots\\}$, i.e. it is a
specification of both the primary and intermediate angular momenta. The symbol
$\mathcal{D}^{\rm P}$ is a product of $\sqrt{2\ell_{i}+1}$ over all primary
angular momenta (defined in [1] equation 35). Briefly, a primary angular
momentum is one that corresponds to a spherical harmonic entering the
expansion; we always denote them with a single numerical subscript, e.g.
$\ell_{1}$. An intermediate angular momentum is one representing the coupling
of either two primaries, or a primary and an intermediate, and always will
appear with a subscript involving multiple numbers, e.g.
$\ell_{12},\ell_{123}$, etc. The intermediates do not have spherical harmonics
associated with them. The need for intermediates is discussed more fully in
[1]; see also [15] and [16].

We now observe that expanding both sides in $p$ and requiring power-by-power
equality as we have in the previous cases, we will have no way to pick out any
particular $\ell_{12}$, but will be left with the sum over all of those
allowed. This is because the sBFs on the right-hand side have no index
$\ell_{12}$.

###  3.2 Picking Out the Intermediates

The failure of the naive approach above means we require another way to
proceed.

Since each of our angular momenta before were associated with a unit vector
$\hat{p}$ that is integrated out, we now consider also associating the
intermediate angular momentum with a vector, which we denote
$\boldsymbol{u}_{12}$. Let us examine the expansion

$\displaystyle
e^{i\boldsymbol{p}\cdot(\boldsymbol{r}_{1}+\boldsymbol{r}_{2})}e^{i\boldsymbol{p}\cdot\boldsymbol{u}_{12}}e^{i\boldsymbol{p}^{\prime}\cdot(\boldsymbol{r}_{3}+\boldsymbol{r}_{4})}e^{i\boldsymbol{p}^{\prime}\cdot\boldsymbol{u}_{12}}=(4\pi)^{4}\sum_{\ell\ell^{\prime}}\;\sum_{\ell_{12}\ell_{12}^{\prime}}\;\sum_{{\rm
all}\;m}i^{\ell+\ell^{\prime}+\ell_{12}+\ell_{12}^{\prime}}j_{\ell}(pR)j_{\ell^{\prime}}(p^{\prime}R^{\prime})Y_{\ell
m}(\hat{p})Y_{\ell m}^{*}(\hat{R})$ $\displaystyle\times
Y_{\ell^{\prime}m^{\prime}}(\hat{p}^{\prime})Y_{\ell^{\prime}m^{\prime}}^{*}(\hat{R}^{\prime})j_{\ell_{12}}(pu_{12})j_{\ell^{\prime}_{12}}(p^{\prime}u_{12})Y_{\ell_{12}m_{12}}(\hat{p})Y_{\ell_{12}m_{12}}^{*}(\hat{u}_{12})Y_{\ell_{12}^{\prime}m_{12}^{\prime}}(\hat{p}^{\prime})Y_{\ell_{12}^{\prime}m_{12}^{\prime}}^{*}(\hat{u}_{12}),$
(36)

where $\boldsymbol{R}\equiv\boldsymbol{r}_{1}+\boldsymbol{r}_{2}$ and
$\boldsymbol{R^{\prime}}\equiv\boldsymbol{r}_{3}+\boldsymbol{r}_{4}$.
Integrating equation (36) over $d\Omega_{p}d\Omega_{p^{\prime}}$ gives
Kronecker deltas on the right-hand side as

$\displaystyle\int d\Omega_{p}d\Omega_{p^{\prime}}Y_{\ell
m}(\hat{p})Y_{\ell_{12}m_{12}}(\hat{p})Y_{\ell^{\prime}m^{\prime}}(\hat{p}^{\prime})Y_{\ell_{12}^{\prime}m_{12}^{\prime}}(\hat{p}^{\prime})=\delta^{\rm
K}_{\ell\ell_{12}}\delta^{\rm K}_{m-m_{12}}\delta^{\rm
K}_{\ell^{\prime}\ell^{\prime}_{12}}\delta^{\rm
K}_{m^{\prime}-m_{12}^{\prime}}(-1)^{m_{12}+m^{\prime}_{12}},$ (37)

which will eliminate several sums. We used that777https://dlmf.nist.gov/14.30
$Y_{\ell m}=(-1)^{m}Y_{\ell-m}^{*}$ and then the spherical harmonics’
orthogonality relation.

We then integrate our result again, over $d\Omega_{12}$. After simplifying
what results and then applying the spherical harmonic addition theorem to
convert it from a product of harmonics into a Legendre polynomial, we finally
find

$\displaystyle\int
d\Omega_{p}d\Omega_{p^{\prime}}d\Omega_{12}\;e^{i\boldsymbol{p}\cdot(\boldsymbol{r}_{1}+\boldsymbol{r}_{2})}e^{i\boldsymbol{p}\cdot\boldsymbol{u}_{12}}e^{i\boldsymbol{p}^{\prime}\cdot(\boldsymbol{r}_{3}+\boldsymbol{r}_{4})}e^{i\boldsymbol{p}^{\prime}\cdot\boldsymbol{u}_{12}}$
$\displaystyle=(4\pi)^{3}\sum_{\ell}(2\ell+1)\mathcal{L}_{\ell}(\hat{R}\cdot\hat{R}^{\prime})j_{\ell}(pR)j_{\ell}(p^{\prime}R^{\prime})j_{\ell}(pu_{12})j_{\ell}(p^{\prime}u_{12}).$
(38)

We note that the phase from equation (36) cancelled because
$\ell+\ell_{12}+\ell^{\prime}+\ell_{12}^{\prime}=2(\ell_{12}+\ell^{\prime}_{12})$
due to equation (37) and then the $d\Omega_{12}$ integration set
$\ell_{12}^{\prime}=\ell_{12}$, so the power of $i$ in equation (36) becomes
an even number divisible by two.

We now seek to relate this set of plane waves to the $N=4$ isotropic basis
functions. We have, by rewriting each plane wave of a sum as the product of
plane waves prior to expanding, that

$\displaystyle
e^{i\boldsymbol{p}\cdot(\boldsymbol{r}_{1}+\boldsymbol{r}_{2})}e^{i\boldsymbol{p}\cdot\boldsymbol{u}_{12}}e^{i\boldsymbol{p}^{\prime}\cdot(\boldsymbol{r}_{3}+\boldsymbol{r}_{4})}e^{i\boldsymbol{p}^{\prime}\cdot\boldsymbol{u}_{12}}$
$\displaystyle=(4\pi)^{6}\sum_{\ell_{1}\ell_{2}\ell_{12}}\;\sum_{\ell_{3}\ell_{4}\ell_{12}^{\prime}}\;\sum_{{\rm
all\;}m}i^{\ell_{1}+\ell_{2}+\ell_{12}+\ell_{3}+\ell_{4}+\ell_{12}^{\prime}}j_{\ell_{1}}(pr_{1})j_{\ell_{2}}(pr_{2})j_{\ell_{12}}(pu_{12})j_{\ell_{3}}(p^{\prime}r_{3})j_{\ell_{4}}(p^{\prime}r_{4})j_{\ell_{12}^{\prime}}(p^{\prime}u_{12})$
$\displaystyle\qquad\times
Y_{\ell_{1}m_{1}}(\hat{p})Y_{\ell_{1}m_{1}}^{*}(\hat{r}_{1})Y_{\ell_{2}m_{2}}(\hat{p})Y_{\ell_{2}m_{2}}^{*}(\hat{r}_{2})Y_{\ell_{12}m_{12}}(\hat{p})Y_{\ell_{12}m_{12}}^{*}(\hat{u}_{12})Y_{\ell_{3}m_{3}}(\hat{p}^{\prime})Y_{\ell_{3}m_{3}}^{*}(\hat{r}_{3})$
$\displaystyle\qquad\times
Y_{\ell_{4}m_{4}}(\hat{p}^{\prime})Y_{\ell_{4}m_{4}}^{*}(\hat{r}_{4})Y_{\ell_{12}^{\prime}m_{12}^{\prime}}(\hat{p}^{\prime})Y_{\ell_{12}^{\prime}m_{12}^{\prime}}^{*}(\hat{u}_{12}).$
(39)

We now perform the same manipulations on this relation as we did on equation
(36), to wit, integrating over $d\Omega_{p}d\Omega_{p^{\prime}}d\Omega_{12}$.
Integrating over $d\Omega_{p}d\Omega_{p^{\prime}}$ replaces the product of
spherical harmonics in $\ell_{1},\ell_{2},\ell_{12}$ with a single Gaunt
integral, and likewise for those in $\ell_{3},\ell_{4},\ell_{12}^{\prime}$.
Integrating what results over $d\Omega_{12}$, gives a product of Kronecker
deltas on the right-hand side as $\delta_{\ell_{12}\ell_{12}^{\prime}}^{\rm
K}\delta_{m_{12}-m_{12}^{\prime}}^{\rm K}(-1)^{m_{12}}$.

Inserting these results and simplifying, we obtain

$\displaystyle\int
d\Omega_{p}d\Omega_{p^{\prime}}d\Omega_{12}\;e^{i\boldsymbol{p}\cdot(\boldsymbol{r}_{1}+\boldsymbol{r}_{2})}e^{i\boldsymbol{p}\cdot\boldsymbol{u}_{12}}e^{i\boldsymbol{p}^{\prime}\cdot(\boldsymbol{r}_{3}+\boldsymbol{r}_{4})}e^{i\boldsymbol{p}^{\prime}\cdot\boldsymbol{u}_{12}}=(4\pi)^{6}$
$\displaystyle\times\sum_{\ell_{1}\ell_{2}\ell_{12}}\;\sum_{\ell_{3}\ell_{4}}\;\sum_{m_{1}m_{2}m_{3}}(-1)^{\ell_{12}}i^{\ell_{1}+\ell_{2}+\ell_{3}+\ell_{4}}j_{\ell_{1}}(pr_{1})j_{\ell_{2}}(pr_{2})j_{\ell_{3}}(p^{\prime}r_{3})j_{\ell_{4}}(p^{\prime}r_{4})j_{\ell_{12}}(pu_{12})j_{\ell_{12}}(p^{\prime}u_{12})$
$\displaystyle\qquad\times\mathcal{G}_{\ell_{1}\ell_{2}\ell_{12}}^{m_{1}m_{2}m_{12}}\mathcal{G}_{\ell_{12}\ell_{3}\ell_{4}}^{-m_{12}m_{3}m_{4}}(-1)^{m_{12}}Y_{\ell_{1}m_{1}}^{*}(\hat{r}_{1})Y_{\ell_{2}m_{2}}^{*}(\hat{r}_{2})Y_{\ell_{3}m_{3}}^{*}(\hat{r}_{3})Y_{\ell_{4}m_{4}}^{*}(\hat{r}_{4}).$
(40)

There is no sum over $m_{12}$ because the Gaunt integrals fix it; this is also
the case for $m_{4}$.

Inserting the definition of the $N=4$ isotropic functions ([1] equations 4, 5,
and 12), we find

$\displaystyle\int
d\Omega_{p}d\Omega_{p^{\prime}}d\Omega_{12}\;e^{i\boldsymbol{p}\cdot(\boldsymbol{r}_{1}+\boldsymbol{r}_{2})}e^{i\boldsymbol{p}\cdot\boldsymbol{u}_{12}}e^{i\boldsymbol{p}^{\prime}\cdot(\boldsymbol{r}_{3}+\boldsymbol{r}_{4})}e^{i\boldsymbol{p}^{\prime}\cdot\boldsymbol{u}_{12}}$
(41)
$\displaystyle=(4\pi)^{5}\sum_{\ell_{1}\ell_{2}\ell_{12}}\sum_{\ell_{3}\ell_{4}}(-1)^{(\ell_{1}+\ell_{2}+\ell_{3}+\ell_{4})/2}j_{\ell_{1}}(pr_{1})j_{\ell_{2}}(pr_{2})j_{\ell_{3}}(p^{\prime}r_{3})j_{\ell_{4}}(p^{\prime}r_{4})j_{\ell_{12}}(pu_{12})j_{\ell_{12}}(p^{\prime}u_{12})\mathcal{D}^{\rm
P}$
$\displaystyle\qquad\times(-1)^{\ell_{1}+\ell_{2}+\ell_{3}+\ell_{4}}\sqrt{2\ell_{12}+1}\begin{pmatrix}\ell_{1}&\ell_{2}&\ell_{12}\\\
0&0&0\end{pmatrix}\begin{pmatrix}\ell_{12}&\ell_{3}&\ell_{4}\\\
0&0&0\end{pmatrix}\mathcal{P}_{\ell_{1}\ell_{2}(\ell_{12})\ell_{3}{\ell_{4}}}(\hat{r}_{1},\hat{r}_{2},\hat{r}_{3},\hat{r}_{4}).$

We used that $\ell_{1}+\ell_{2}+\ell_{12}$ and $\ell_{12}+\ell_{3}+\ell_{4}$
are even (due to the two zero-spin 3-$j$ symbols) to conclude that the parity
of $\ell_{1}+\ell_{2}$ is that of $\ell_{3}+\ell_{4}$, and hence their sum is
even and we may rewrite

$\displaystyle
i^{-(\ell_{1}+\ell_{2}+\ell_{3}+\ell_{4})}=i^{2\times-(\ell_{1}+\ell_{2}+\ell_{3}+\ell_{4})/2}=(-1)^{-(\ell_{1}+\ell_{2}+\ell_{3}+\ell_{4})/2}=(-1)^{(\ell_{1}+\ell_{2}+\ell_{3}+\ell_{4})/2}.$
(42)

The phase in the last line of equation (41) comes from the isotropic
functions’ behavior under conjugation ([1] equation 11), as the harmonics in
equation (40) are conjugated and thus give rise to a conjugated isotropic
function ([1] equation 10), which we then transform to a non-conjugated one.
The phases due to the intermediate angular momentum $\ell_{12}$ and its spin
$m_{12}$ are absorbed in the definition of the isotropic functions ([1]
equations 12 and 13, the factor $\kappa$).

Now, equating our two different ways of expanding the plane wave, Eqs. (38)
and (41), and eliminating appropriate factors of $4\pi$ on each side, we find

$\displaystyle
H_{4}\equiv\sum_{\ell}(2\ell+1)\mathcal{L}_{\ell}(\hat{R}\cdot\hat{R}^{\prime})j_{\ell}(pR)j_{\ell}(p^{\prime}R^{\prime})j_{\ell}(pu_{12})j_{\ell}(p^{\prime}u_{12})=$
$\displaystyle(4\pi)^{2}\sum_{\ell_{1}\ell_{2}\ell_{12}}\sum_{\ell_{3}\ell_{4}}(-1)^{(\ell_{1}+\ell_{2}+\ell_{3}+\ell_{4})/2}j_{\ell_{1}}(pr_{1})j_{\ell_{2}}(pr_{2})j_{\ell_{3}}(p^{\prime}r_{3})j_{\ell_{4}}(p^{\prime}r_{4})j_{\ell_{12}}(pu_{12})j_{\ell_{12}}(p^{\prime}u_{12})\mathcal{D}^{\rm
P}\sqrt{2\ell_{12}+1}$
$\displaystyle\qquad\times(-1)^{\ell_{1}+\ell_{2}+\ell_{3}+\ell_{4}}\begin{pmatrix}\ell_{1}&\ell_{2}&\ell_{12}\\\
0&0&0\end{pmatrix}\begin{pmatrix}\ell_{12}&\ell_{3}&\ell_{4}\\\
0&0&0\end{pmatrix}\mathcal{P}_{\ell_{1}\ell_{2}(\ell_{12})\ell_{3}{\ell_{4}}}(\hat{r}_{1},\hat{r}_{2},\hat{r}_{3},\hat{r}_{4}).$
(43)

The zero-spin 3-$j$ symbols imply that the parity of $\ell_{12}$ must equal
that of $\ell_{1}+\ell_{2}$ and also that of $\ell_{3}+\ell_{4}$; hence the
parity of $\ell_{1}+\ell_{2}+\ell_{3}+\ell_{4}$ must be even. This sum gives
the parity of $\mathcal{P}_{\ell_{1}\ell_{2}(\ell_{12})\ell_{3}\ell_{4}}$.
This means that, as was the case for $N=3$ as well, our approach thus far only
allows us to explore even-parity isotropic basis functions. Furthermore, the
parity of the intermediate must be equal to that of the first two primary
angular momenta and of the last two. Within the category of overall even basis
functions, this gives a further restriction on what intermediates can be
generated. We discuss how to generate odd-parity isotropic basis functions in
§5.

We now inspect equation (43) to confirm that it will give us a way to control
both the primary angular momenta and the intermediate. As usual, we must,
after expanding each side in a power series (and not necessarily a priori
restricting $\ell=0$ on the left-hand side) match the overall power of $p$ on
each side; here we demand that the power of $p^{\prime}$ match as well. Then,
comparing powers of the $r_{i}$ on each side and demanding a match fixes the
$\ell_{i}$, up to the need to account for sub-leading terms in the sBF
expansions and to subtract them and the basis functions that go with them off,
as we did before. The only addition beyond previous sections is that now, the
powers of $u_{12}$ may also be matched on each side, and this enables
controlling the intermediate, $\ell_{12}$. We now show two examples of how
this proceeds in practice.

###  3.3 $N=4$ Examples

####  3.3.1 $\mathcal{P}_{10(1)01}$

Here, looking at equation (43), we have $\ell_{1}=\ell_{4}=1$, and this means
we need combinations on the left- and right-hand sides that behave as
$r_{1}r_{4}$. We also observe that, since $\ell_{12}=1$, to match the powers
of $u_{12}$ on the right-hand side that we will have from expanding the two
$j_{1}$’s to leading order, we need $u_{12}^{2}$ on the left-hand side. This
implies $\ell=1$ on the left-hand side, if we look at the leading order of the
sBf series that would result. In particular, since the series for $j_{0}$ is
even parity, the sub-leading terms if we set $\ell=0$ on the left-hand side
would give $u_{12}^{4}$ and higher powers, so $\ell=0$ cannot be relevant.

Setting $\ell=1$, we thus have the left-hand side of equation (43) as

$\displaystyle H_{4,10(1)01}\to
3\hat{R}\cdot\hat{R}^{\prime}\,j_{1}(pR)j_{1}(p^{\prime}R^{\prime})j_{1}(pu_{12})j_{1}(p^{\prime}u_{12}),$
(44)

where our additional subscript on $H_{4}$ indicates that it is just the terms
in the sum of equation (43) relevant for $\mathcal{P}_{10(1)01}$. We adopt
this notation here for clarity now that we have an intermediate as well,
rather than the vertical bar and subscripted $r_{i}$ of previous sections. In
the vertical bar notation, we would have
$H_{4}\big{|}_{r_{1}r_{4}u_{12}^{2}}$.

Now expanding each sBF to leading order, we find from equation (44) that

$\displaystyle H_{4,10(1)01}\to
3\hat{R}\cdot\hat{R}^{\prime}\,\frac{pR}{3}\,\frac{p^{\prime}R^{\prime}}{3}\,\frac{pu_{12}}{3}\,\frac{p^{\prime}u_{12}}{3}=\frac{1}{27}\,{\bf
R}\cdot{\bf R^{\prime}}\,p^{2}p^{\prime 2}u_{12}^{2}.$ (45)

Next, we expand the dot product of ${\bf R}\cdot{\bf R^{\prime}}$ and retain
only the term in $r_{1}r_{4}$. We find that this term has coefficient unity,
so we have

$\displaystyle H_{4,10(1)01}\to
3\hat{R}\cdot\hat{R}^{\prime}\,\frac{pR}{3}\frac{p^{\prime}R^{\prime}}{3}\frac{pu_{12}}{3}\frac{p^{\prime}u_{12}}{3}=\frac{1}{27}\,r_{1}r_{4}\,(\hat{r}_{1}\cdot{\hat{r}_{4}})p^{2}p^{\prime
2}u_{12}^{2}.$ (46)

We now set this result equal to the right-hand side of equation (43). We first
expand each sBF in a series, and note that demanding equality of the powers of
$p$, $p^{\prime}$, $u_{12}$, and the $r_{i}$ removes the sums over all
$\ell_{i}$ and $\ell_{12}$. We then cancel the powers $p$, $p^{\prime}$,
$u_{12}$, and the $r_{i}$, and evaluate the factor $\mathcal{D}^{\rm P}$,
$\sqrt{2\ell_{12}+1}$, and the 3-$j$ symbols. We find

$\displaystyle\frac{1}{27}\,\hat{r}_{1}\cdot\hat{r}_{4}=(4\pi)^{2}(-1)\frac{1}{3\cdot
27}\sqrt{3}^{3}\left(-\frac{1}{\sqrt{3}}\right)^{2}\mathcal{P}_{10(1)01}.$
(47)

Simplifying and solving for $\mathcal{P}_{10(1)01}$, we obtain

$\displaystyle\mathcal{P}_{10(1)01}(\hat{r}_{1},\hat{r}_{2},\hat{r}_{3},\hat{r}_{4})=-\frac{\sqrt{3}}{(4\pi)^{2}}\,\hat{r}_{1}\cdot\hat{r}_{4},$
(48)

in agreement with [1] Appendix A.3., second line.

####  3.3.2 $\mathcal{P}_{11(0)11}$

Here, we use the same approach as in the previous subsection. Now we see that
on the right-hand side of equation (43), we will need $p^{2}p^{\prime
2}u_{12}^{0}$ and combinations behaving as $r_{1}r_{2}r_{3}r_{4}$. We might at
first think $\ell=1$ on the left-hand side is the correct choice, as this
gives the right powers of $p$ and $p^{\prime}$ if the sBFs there are expanded
to leading order. However, this would give the wrong power of $u_{12}$; we
would have $u_{12}^{2}$ at leading order, but we need $u_{12}^{0}$. If we set
$\ell=0$ and take the the second-order terms from expanding
$j_{0}(pR)j_{0}(p^{\prime}R^{\prime})$ but the leading-order terms from
expanding $j_{0}(pu_{12})j_{0}(p^{\prime}u^{\prime}_{12})$, we will then have
the desired powers of $p$, $p^{\prime}$, and $u_{12}$. We have

$\displaystyle
H_{4,11(0)11}\to\frac{1}{36}(pR)^{2}(p^{\prime}R^{\prime})^{2}.$ (49)

In the vertical bar notation, we would have
$H_{4}\big{|}_{r_{1}r_{2}r_{3}r_{4}}$. We now expand the powers of $R$ and
$R^{\prime}$ and retain only the term in $r_{1}r_{2}r_{3}r_{4}$, which has a
coefficient of $4$. We find

$\displaystyle H_{4,11(0)11}\to\frac{1}{9}p^{2}p^{\prime
2}r_{1}r_{2}r_{3}r_{4}\,(\hat{r}_{1}\cdot\hat{r}_{2})(\hat{r}_{3}\cdot\hat{r}_{4}).$
(50)

Returning to equation (43 and demanding equality of all powers, we see this
eliminates the sums over the $\ell_{i}$ and over $\ell_{12}$, and we have,
cancelling the factor of $p^{2}p^{\prime 2}r_{1}r_{2}r_{3}$, that

$\displaystyle\frac{1}{9}\,(\hat{r}_{1}\cdot\hat{r}_{2})(\hat{r}_{3}\cdot\hat{r}_{4})=(4\pi)^{2}\frac{1}{3^{4}}\sqrt{3}^{4}\left(-\frac{1}{\sqrt{3}}\right)^{2}\mathcal{P}_{11(0)11}(\hat{r}_{1},\hat{r}_{2},\hat{r}_{3},\hat{r}_{4}).$
(51)

Solving for $\mathcal{P}_{11(0)11}$ and simplifyng, we obtain

$\displaystyle\mathcal{P}_{11(0)11}(\hat{r}_{1},\hat{r}_{2},\hat{r}_{3},\hat{r}_{4})=\frac{3}{(4\pi)^{2}}\,(\hat{r}_{1}\cdot\hat{r}_{2})(\hat{r}_{3}\cdot\hat{r}_{4}),$
(52)

in agreement with [1] Appendix A.3., first line.

##  4 Generating Function for $N=5$

Let us now apply the same approach to deriving the generating function for the
$N=5$ isotropic basis functions. We consider

$\displaystyle P_{5}\equiv
e^{i\boldsymbol{p}\cdot(\boldsymbol{r}_{1}+\boldsymbol{r}_{2})}e^{i\boldsymbol{p}^{\prime}\cdot\boldsymbol{r}_{3}}e^{i\boldsymbol{p}^{\prime\prime}\cdot(\boldsymbol{r}_{4}+\boldsymbol{r}_{5})}e^{i\boldsymbol{p}\cdot\boldsymbol{u}_{12}}e^{i\boldsymbol{p}^{\prime}\cdot\boldsymbol{u}_{12}}e^{i\boldsymbol{p}^{\prime}\cdot\boldsymbol{u}_{123}}e^{i\boldsymbol{p}^{\prime\prime}\cdot\boldsymbol{u}_{123}}.$
(53)

Defining $R\equiv|\boldsymbol{r}_{1}+\boldsymbol{r}_{2}|$,
$R^{\prime}\equiv|\boldsymbol{r}_{3}|$, and
$R^{\prime\prime}\equiv|\boldsymbol{r}_{4}+\boldsymbol{r}_{5}|$, we have

$\displaystyle P_{5}=$
$\displaystyle(4\pi)^{7}\sum_{\ell\ell^{\prime}\ell^{\prime\prime}}\sum_{\ell_{12}\ell_{12}^{\prime}}\sum_{\ell_{123}^{\prime}\ell_{123}^{\prime\prime}}\sum_{{\rm
all\;m}}i^{\ell+\ell^{\prime}+\ell^{\prime\prime}+\ell_{12}+\ell_{12}^{\prime}+\ell_{123}^{\prime}+\ell_{123}^{\prime\prime}}j_{\ell}(pR)j_{\ell^{\prime}}(p^{\prime}R^{\prime})j_{\ell^{\prime\prime}}(p^{\prime\prime}R^{\prime\prime})j_{\ell_{12}}(pu_{12})j_{\ell_{12}^{\prime}}(p^{\prime}u_{12})$
$\displaystyle\times
j_{\ell_{123}^{\prime}}(p^{\prime}u_{123})j_{\ell_{123}^{\prime\prime}}(p^{\prime\prime}u_{123})Y_{\ell
m}(\hat{R})Y_{\ell
m}^{*}(\hat{p})Y_{\ell^{\prime}m^{\prime}}(\hat{R}^{\prime})Y_{\ell^{\prime}m^{\prime}}^{*}(\hat{p}^{\prime})Y_{\ell^{\prime\prime}m^{\prime\prime}}(\hat{R}^{\prime\prime})Y_{\ell^{\prime\prime}m^{\prime\prime}}^{*}(\hat{p}^{\prime\prime})$
(54) $\displaystyle\times
Y_{\ell_{12}m_{12}}(\hat{u}_{12})Y_{\ell_{12}m_{12}}^{*}(\hat{p})Y_{\ell_{12}^{\prime}m_{12}^{\prime}}(\hat{u}_{12})Y_{\ell_{12}^{\prime}m_{12}^{\prime}}^{*}(\hat{p}^{\prime})Y_{\ell_{123}^{\prime}m_{123}^{\prime}}(\hat{u}_{123})Y_{\ell_{123}^{\prime}m_{123}^{\prime}}^{*}(\hat{p}^{\prime})$
$\displaystyle\times
Y_{\ell_{123}^{\prime\prime}m_{123}^{\prime\prime}}(\hat{u}_{123})Y_{\ell_{123}^{\prime\prime}m_{123}^{\prime\prime}}^{*}(\hat{p}^{\prime\prime}).$

Now performing an integral over
$d\Omega_{p}d\Omega_{p^{\prime}}d\Omega_{p^{\prime\prime}}$ we obtain

$\displaystyle P_{5}$ $\displaystyle=(4\pi)^{7}\sum_{\rm
all}i^{(\cdots)}j_{\ell}(pR)\cdots
j_{\ell_{123}^{\prime\prime}}(p^{\prime\prime}u_{123})Y_{\ell
m}(\hat{R})Y_{\ell^{\prime}m^{\prime}}(\hat{R}^{\prime})Y_{\ell_{12}m_{12}}(\hat{u}_{12})Y_{\ell_{12}^{\prime}m_{12}^{\prime}}(\hat{u}_{12})$
(55) $\displaystyle\times
Y_{\ell_{123}^{\prime}m_{123}^{\prime}}(\hat{u}_{123})Y_{\ell_{123}^{\prime\prime}m_{123}^{\prime\prime}}(\hat{u}_{123})Y_{\ell^{\prime\prime}m^{\prime\prime}}(\hat{R}^{\prime\prime})\delta^{\rm
K}_{\ell\ell_{12}}(-1)^{m_{12}}\delta^{\rm
K}_{m-m_{12}}\mathcal{G}_{\ell^{\prime}\ell_{12}^{\prime}\ell_{123}^{\prime}}^{m^{\prime}m_{12}^{\prime}m_{123}^{\prime}}\delta^{\rm
K}_{\ell^{\prime\prime}\ell_{123}^{\prime\prime}}(-1)^{m_{123}^{\prime\prime}}\delta^{\rm
K}_{m^{\prime\prime}-m_{123}^{\prime\prime}},$

where the subscript “${\rm all}$” on the sum indicates to sum over all the
subscripts on the right-hand side, and the $\cdots$ in the phase indicate the
same phase as in equation (4). The ellipsis $(\cdots)$ in between the sBfs
indicate they are all copied over from equation (4). We note that there were
just two harmonics involving $\hat{p}$ and also just two involving
$\hat{p^{\prime\prime}}$, hence the Kronecker deltas for those angular
momenta, as well as the phases and flipped signs on the spins due to using the
conjugate identity. We also used that the Gaunt integral is real, so taking
its conjugate, which is equivalent to integrating over three conjugated
spherical harmonics, leaves it unchanged.

Simplifying equation (55), we obtain

$\displaystyle P_{5}$ $\displaystyle=(4\pi)^{7}\sum_{\rm
all}i^{2(\ell+\ell^{\prime\prime})+\ell^{\prime}+\ell_{12}^{\prime}+\ell_{123}^{\prime}}j_{\ell}(pR)\cdots
j_{\ell^{\prime\prime}}(p^{\prime\prime}u_{123})Y_{\ell
m}(\hat{R})Y_{\ell^{\prime}m^{\prime}}(\hat{R}^{\prime})Y_{\ell^{\prime\prime}m^{\prime\prime}}(\hat{R}^{\prime\prime})\mathcal{G}_{\ell^{\prime}\ell_{12}^{\prime}\ell_{123}^{\prime}}^{m^{\prime}m_{12}^{\prime}m_{123}^{\prime}}$
$\displaystyle\qquad\times Y_{\ell
m}^{*}(\hat{u}_{12})Y_{\ell_{12}^{\prime}m_{12}^{\prime}}(\hat{u}_{12})Y_{\ell_{123}^{\prime}m_{123}^{\prime}}(\hat{u}_{123})Y_{\ell^{\prime\prime}m^{\prime\prime}}^{*}(\hat{u}_{123}).$
(56)

We now further integrate over $d\Omega_{12}d\Omega_{123}$ and obtain

$\displaystyle P_{5}$ $\displaystyle=(4\pi)^{7}\sum_{\rm
all}(-1)^{\ell+\ell^{\prime}}i^{\ell+\ell^{\prime}+\ell^{\prime\prime}}j_{\ell}(pR)j_{\ell^{\prime}}(p^{\prime}R^{\prime})j_{\ell^{\prime\prime}}(p^{\prime\prime}R^{\prime\prime})$
$\displaystyle\times
j_{\ell}(pu_{12})j_{\ell}(p^{\prime}u_{12})j_{\ell^{\prime\prime}}(p^{\prime}u_{123})j_{\ell^{\prime\prime}}(p^{\prime\prime}u_{123})Y_{\ell
m}(\hat{R})Y_{\ell^{\prime}m^{\prime}}(\hat{R}^{\prime})Y_{\ell^{\prime\prime}m^{\prime\prime}}(\hat{R}^{\prime\prime})\mathcal{G}_{\ell^{\prime}\ell\ell^{\prime\prime}}^{m^{\prime}mm^{\prime\prime}}.$
(57)

We note that the Gaunt integral above contains a zero-spin 3-$j$ symbol, which
in turn guarantees that $\ell+\ell^{\prime}+\ell^{\prime\prime}$ is even and
hence that $i$ to this power is real. For this reason we may also reorder the
Gaunt to be in order of ascending primes without picking up any overall phase.
We observe that the generating function for the $N=5$ isotropic basis
functions thus turns out to be a tripolar spherical harmonic in
$\hat{R},\,\hat{R}^{\prime},$ and $\hat{R}^{\prime\prime}$. It now remains to
obtain the right-hand side of the expansion of the plane wave here, which will
explicitly display the relationship of the above to the $N=5$ basis functions.
This is extensive but conceptually much simpler than was the expansion for the
left-hand side. We simply write

$\displaystyle P_{5}$
$\displaystyle=e^{i\boldsymbol{p}\cdot\boldsymbol{r}_{1}}e^{i\boldsymbol{p}\cdot\boldsymbol{r}_{2}}e^{i\boldsymbol{p}^{\prime}\cdot\boldsymbol{r}_{3}}e^{i\boldsymbol{p}^{\prime\prime}\cdot\boldsymbol{r}_{4}}e^{i\boldsymbol{p}^{\prime\prime}\cdot\boldsymbol{r}_{5}}e^{i\boldsymbol{p}^{\prime}\cdot\boldsymbol{u}_{12}}e^{i\boldsymbol{p}\cdot\boldsymbol{u}_{12}}e^{i\boldsymbol{p}^{\prime}\cdot\boldsymbol{u}_{123}}e^{i\boldsymbol{p}^{\prime\prime}\cdot\boldsymbol{u}_{123}}=$
$\displaystyle(4\pi)^{9}\sum_{\rm
all}i^{\ell_{1}+\ell_{2}+\ell_{3}+\ell_{4}+\ell_{5}+\ell_{12}+\ell_{12}^{\prime}+\ell_{123}^{\prime}+\ell_{123}^{\prime\prime}}$
$\displaystyle\times
j_{\ell_{1}}(pr_{1})j_{\ell_{2}}(pr_{2})j_{\ell_{3}}(p^{\prime}r_{3})j_{\ell_{4}}(p^{\prime\prime}r_{4})j_{\ell_{5}}(p^{\prime\prime}r_{5})j_{\ell_{12}}(pu_{12})j_{\ell^{\prime}_{12}}(p^{\prime}u_{12})j_{\ell_{123}^{\prime}}(p^{\prime}u_{123})j_{\ell_{123}^{\prime\prime}}(p^{\prime\prime}u_{123})$
$\displaystyle\times\mathcal{Y}_{\ell_{1}\ell_{2}\ell_{12}}^{m_{1}m_{2}m_{12}*}(\hat{p})\mathcal{Y}_{\ell_{3}\ell^{\prime}_{12}\ell_{123}^{\prime}}^{m_{3}m_{12}^{\prime}m_{123}^{\prime}*}(\hat{p}^{\prime})\mathcal{Y}_{\ell_{4}\ell_{5}\ell_{123}^{\prime\prime}}^{m_{4}m_{5}m_{123}^{\prime\prime}*}(\hat{p}^{\prime\prime})Y_{\ell_{1}m_{1}}(\hat{r}_{1})Y_{\ell_{2}m_{2}}(\hat{r}_{2})Y_{\ell_{3}m_{3}}(\hat{r}_{3})Y_{\ell_{4}m_{4}}(\hat{r}_{4})Y_{\ell_{5}m_{5}}(\hat{r}_{5})$
$\displaystyle\times
Y_{\ell_{12}m_{12}}(\hat{u}_{12})Y_{\ell_{12}^{\prime}m_{12}^{\prime}}(\hat{u}_{12})Y_{\ell_{123}^{\prime}m_{123}^{\prime}}(\hat{u}_{123})Y_{\ell_{123}^{\prime\prime}m_{123}^{\prime\prime}}(\hat{u}_{123}),$
(58)

with

$\displaystyle\mathcal{Y}_{abc}^{def*}(\hat{x})\equiv
Y_{ad}^{*}(\hat{x})Y_{be}^{*}(\hat{x})Y_{cf}^{*}(\hat{x}).$ (59)

Now integrating equation (4) over
$d\Omega_{p}d\Omega_{p^{\prime}}d\Omega_{p^{\prime\prime}}$ removes the three
$\mathcal{Y}$ in favor of three Gaunt integrals, as

$\displaystyle\int
d\Omega_{p}d\Omega_{p^{\prime}}d\Omega_{p^{\prime\prime}}\,\mathcal{Y}_{\ell_{1}\ell_{2}\ell_{12}}^{m_{1}m_{2}m_{12}*}(\hat{p})\mathcal{Y}_{\ell_{3}\ell^{\prime}_{12}\ell_{123}^{\prime}}^{m_{3}m_{12}^{\prime}m_{123}^{\prime}*}(\hat{p}^{\prime})\mathcal{Y}_{\ell_{4}\ell_{5}\ell_{123}^{\prime\prime}}^{m_{4}m_{5}m_{123}^{\prime\prime}*}(\hat{p}^{\prime\prime})$
(60)
$\displaystyle\qquad=\mathcal{G}_{\ell_{1}\ell_{2}\ell_{12}}^{m_{1}m_{2}m_{12}}\mathcal{G}_{\ell_{3}\ell_{12}^{\prime}\ell_{123}^{\prime}}^{m_{3}m_{12}^{\prime}m_{123}^{\prime}}\mathcal{G}_{\ell_{4}\ell_{5}\ell_{123}^{\prime\prime}}^{m_{4}m_{5}m_{123}^{\prime\prime}},$

where the second and third Gaunt integrals above may be non-cyclically
permuted without acquiring any net phase, as the phases from permutation will
vanish because each involves an even power of $(-1)$, as guaranteed by the
presence of a zero-spin 3-$j$ symbol within each Gaunt.

We now integrate over $d\Omega_{12}d\Omega_{123}$, finding that

$\displaystyle\int
d\Omega_{12}d\Omega_{123}\,Y_{\ell_{12}m_{12}}(\hat{u}_{12})Y_{\ell_{12}^{\prime}m_{12}^{\prime}}(\hat{u}_{12})Y_{\ell_{123}^{\prime}m_{123}^{\prime}}(\hat{u}_{123})Y_{\ell_{123}^{\prime\prime}m_{123}^{\prime\prime}}(\hat{u}_{123})$
$\displaystyle\qquad=\delta^{\rm
K}_{\ell_{12}^{\prime}\ell_{12}}(-1)^{m_{12}}\delta^{\rm
K}_{m_{12}^{\prime}-m_{12}}\delta^{\rm
K}_{\ell_{123}^{\prime\prime}\ell_{123}^{\prime}}(-1)^{m_{123}^{\prime}}\delta^{\rm
K}_{m_{123}^{\prime\prime}-m_{123}^{\prime}}.$ (61)

Putting the results above together and rewriting in terms of the $N=5$
isotropic basis functions using [1] equations (4) and (12), we find that

$\displaystyle P_{5}$ $\displaystyle=(4\pi)^{9}\sum_{\rm
all}(-1)^{(\ell_{1}+\ell_{2}+\ell_{3}+\ell_{4}+\ell_{5})/2}j_{\ell_{1}}(pr_{1})j_{\ell_{2}}(pr_{2})j_{\ell_{3}}(p^{\prime}r_{3})j_{\ell_{4}}(p^{\prime\prime}r_{4})j_{\ell_{5}}(p^{\prime\prime}r_{5})j_{\ell_{12}}(pu_{12})j_{\ell_{12}}(p^{\prime}u_{12})$
$\displaystyle\times
j_{\ell_{123}}(p^{\prime}u_{123})j_{\ell_{123}}(p^{\prime\prime}u_{123})\begin{pmatrix}\ell_{1}&\ell_{2}&\ell_{12}\\\
0&0&0\end{pmatrix}\begin{pmatrix}\ell_{12}&\ell_{3}&\ell_{123}\\\
0&0&0\end{pmatrix}\begin{pmatrix}\ell_{123}&\ell_{4}&\ell_{5}\\\
0&0&0\end{pmatrix}\mathcal{D}^{\rm P}$
$\displaystyle\times\sqrt{2\ell_{12}+1}\sqrt{2\ell_{123}+1}\,\mathcal{P}_{\ell_{1}\ell_{2}(\ell_{12})\ell_{3}(\ell_{123})\ell_{4}\ell_{5}}(\hat{r}_{1},\hat{r}_{2},\hat{r}_{3},\hat{r}_{4},\hat{r}_{5}).$
(62)

We note that as a last step we have renamed $\ell_{123}^{\prime}$ as
$\ell_{123}$ without loss of generality and to be consistent with the
convention in [1]. By comparing the equation above with equation (57) we see
that we have obtained the generating function for the $N=5$ isotropic basis
functions. We also see that expanding the above in small $p$, $p^{\prime}$,
and $p^{\prime\prime}$ allows us to select a total of the primary $\ell_{i}$,
and then the power of each $r_{i}$ to control the individual $\ell_{i}$. Then,
the powers of respectively $u_{12}$ and $u_{123}$ may be used to control the
intermediate angular momenta $\ell_{12}$ and $\ell_{123}$.

##  5 Generating Function for Odd-Parity Isotropic Functions

As noted earlier, the generating functions thus far presented produce only
even isotropic basis functions. However, in 3D, one can have odd parity basis
functions of three or more arguments, i.e. for $N\geq 3$. Below we outline how
the generating function may be obtained for $N=3$; we defer the $N>3$ cases to
future work. The $N=3$ odd-parity basis functions are already sufficient to
perform a parity violation search with the galaxy 4-Point Correlation
Function, as first proposed in [4] and carried out in [6].

###  5.1 For $N=3$

We notice that $\mathcal{P}_{111}$ is the simplest odd-parity function (recall
there are none for $N=2$). We consider if we can use it as a building block to
produce the other odd parity $N=3$ functions. Let us consider its product with
an arbitrary even-parity $N=3$ function.

The linearization identity for two isotropic functions of the same arguments
is

$\displaystyle\mathcal{P}_{\Lambda}\mathcal{P}_{\Lambda^{\prime}}=\sum_{\Lambda^{\prime\prime}}\mathcal{G}_{\Lambda\Lambda^{\prime}\Lambda^{\prime\prime}}\mathcal{P}_{\Lambda^{\prime\prime}},$
(63)

where $\mathcal{G}$ is a generalized Gaunt integral as defined in [1]; it is
explicitly given for $N=3$ in their equation (62). We set
$\Lambda=(\ell_{1},\ell_{2},\ell_{3})$ and $\Lambda^{\prime}=(1,1,1)$, since
using the lowest-order odd isotropic basis function will convert even
functions into their nearest-lying odd neighbors. We can view
$\mathcal{P}_{111}$ almost as an operator, which raises/lowers an even state
to the proximate odd states. We find

$\displaystyle\mathcal{G}_{\Lambda\Lambda^{\prime}\Lambda^{\prime\prime}}=(4\pi/3)^{-3/2}\bigg{[}\prod_{i=1}^{3}(2\ell_{i}+1)\begin{pmatrix}\ell_{i}&1&\ell_{i}^{\prime\prime}\\\
0&0&0\end{pmatrix}\bigg{]}\left\\{\begin{matrix}\ell_{1}&1&\ell_{1}^{\prime\prime}\\\
\ell_{2}&1&\ell_{2}^{\prime\prime}\\\ \ell_{3}&1&\ell_{3}^{\prime\prime}\\\
\end{matrix}\right\\}.$ (64)

The $3\times 3$ matrix is a Wigner 9-$j$ symbol. There will be two allowed
values of each $\ell^{\prime\prime}_{i}$, namely
$\ell^{\prime\prime}_{i}=\ell_{1}+1,\ell_{1}-1$; thus the right-hand side of
equation (63) will contain eight odd-parity isotropic functions. The
orthogonality relation for 9-$j$ symbols, NIST DLMF
34.7.2,888https://dlmf.nist.gov/34.7 We note that we can use the reflection
symmetry of the 9-$j$ to bring it to the form required by the orthogonality
relation, where the two non-matching $j$’s are the bottom left two elements.
Then orthogonality will work on the coefficients where the bottom right
element is matched, yet some symbols will not have that as the case and they
will give the “anomalous” contribution. shows that if we sum over two of the
$j$’s then we can pick out a desired isotropic function plus an anomalous
contribution that involves other isotropic functions. Doing this with all the
possible choices of the two $j$’s to sum over will result in eight independent
linear equations for the eight odd isotropic functions; these may then be
solved. This will give an expression for each of them in terms of the products
$\mathcal{P}_{\ell_{1}\ell_{2}\ell_{3}}\mathcal{P}_{111}$ weighted by angular
momentum symbols; replacing $\mathcal{P}_{\ell_{1}\ell_{2}\ell_{3}}$ by its
Cartesian form, which can be obtained from the foregoing sections of the
present work, then gives a Cartesian form for the desired odd isotropic
functions.

##  6 Expansion of the Dirac Delta Function

We now outline how to quickly obtain the expansion of the Dirac Delta function
of multiple arguments in terms of the appropriate isotropic basis functions.
We recall that the Delta function has spherical symmetry and is also the 3D
inverse FT of unity; hence

$\displaystyle\delta_{\rm
D}^{[3]}(R)=\int\frac{p^{2}dp}{2\pi^{2}}\;j_{0}(pR).$ (65)

$\delta_{\rm D}^{[3]}$ is a 3D Dirac Delta function, and we note it is
spherically symmetric and hence depends only on $R$, the magnitude of
$\bf{R}$. This may be seen easily by recalling that a 3D Delta function is
just the zero-width limit of a 3D Gaussian, which is spherically symmetric.

The integral above may be shown using the plane wave expansion into spherical
harmonics and spherical Bessel functions along with the spherical symmetry of
the Delta function’s Fourier-space counterpart, unity. Thus, if we integrate
both sides of our generating function relation against $p^{2}dp/(2\pi^{2})$,
we obtain immediately the expansion of the Delta function in the isotropic
basis. Schematically, it is

$\displaystyle\delta_{\rm D}^{[3]}(R)\to\sum_{{\rm
all}\;\ell}\int\frac{p^{2}dp}{2\pi^{2}}\,j_{\ell_{1}}(pr_{1})\cdots
j_{\ell_{N}}(pr_{N})\mathcal{P}^{(N)}_{\ell_{1}\cdots\ell_{N}},$ (66)

where $R\equiv|\sum_{i}^{N}{\bf r}_{i}|$ (we recall equation 3) and the
arguments of the $N$-argument isotropic basis function
$\mathcal{P}^{N}_{\ell_{1}\cdots\ell_{N}}$ are the $\hat{r}_{i}$.

###  6.1 Delta Function for $N=2$ and $N=3$

This is straightforward for $N=2$ and $N=3$, when there are no intermediate
angular momenta; we simply write down the results below. For $N=2$, we have

$\displaystyle\int\frac{p^{2}dp}{2\pi^{2}}\,j_{0}(pR)=(4\pi)^{3/2}\sum_{\ell}\sqrt{2\ell+1}\;\mathcal{P}_{\ell\ell}(\hat{r}_{1},\hat{r}_{2})\int\frac{p^{2}dp}{2\pi^{2}}\,j_{\ell}(pr_{1})j_{\ell}(pr_{2}),$
(67)

and the overlap integral on the right-hand side may be simplified to $1/(4\pi
r_{1}r_{2})\,\delta_{\rm D}^{[1]}(r_{1}-r_{2})$ if desired.

For $N=3$, we have

$\displaystyle\int\frac{p^{2}dp}{2\pi^{2}}\,j_{0}(pR)=$ (68)
$\displaystyle(4\pi)^{2}\sum_{\ell_{1}\ell_{2}\ell_{3}}(-1)^{(\ell_{1}+\ell_{2}+\ell_{3})/2}\mathcal{C}_{\ell_{1}\ell_{2}\ell_{3}}\begin{pmatrix}\ell_{1}&\ell_{2}&\ell_{3}\\\
0&0&0\end{pmatrix}\mathcal{P}_{\ell_{1}\ell_{2}\ell_{3}}(\hat{r}_{1},\hat{r}_{2},\hat{r}_{3})\int\frac{p^{2}dp}{2\pi^{2}}\,j_{\ell_{1}}(pr_{1})j_{\ell_{2}}(pr_{2})j_{\ell_{3}}(pr_{3}).$

The sum of the $\ell_{i}$ is guaranteed to be even by the zero-spin 3-$j$
symbol, thus our simplification of the phase we obtained from the plane wave
expansions.

The triple-sBF integral on the right-hand side can be done using a result in
[17], where they denote it $I$ (more precisely our integral above is
$I/(2\pi^{2})$ under their definition for $I$). We have

$\displaystyle
I(\ell_{1},\ell_{2},\ell_{3};r_{1},r_{2},r_{3})=\frac{\pi\beta(\Delta)}{4r_{1}r_{2}r_{3}}i^{\ell_{1}+\ell_{2}+\ell_{3}}\sqrt{2\ell_{3}+1}\left(\frac{r_{1}}{r_{3}}\right)^{\ell_{3}}\begin{pmatrix}\ell_{1}&\ell_{2}&\ell_{3}\\\
0&0&0\end{pmatrix}^{-1}\sum_{\mathcal{L}}^{\ell_{3}}{{2\ell_{3}}\choose{2\mathcal{L}}}^{1/2}\left(\frac{r_{2}}{r_{1}}\right)^{\mathcal{L}}$
$\displaystyle\times\sum_{\ell}(2\ell+1)\begin{pmatrix}\ell_{1}&\ell_{3}-\mathcal{L}&\ell\\\
0&0&0\end{pmatrix}\begin{pmatrix}\ell_{2}&\mathcal{L}&\ell\\\
0&0&0\end{pmatrix}\left\\{\begin{matrix}\ell_{1}&\ell_{2}&\ell_{3}\\\
\mathcal{L}&\ell_{3}-\mathcal{L}&\ell\end{matrix}\right\\}\mathcal{L}_{\ell}(\Delta),$
(69)

with

$\displaystyle\Delta\equiv\frac{r_{1}^{2}+r_{2}^{2}-r_{3}^{2}}{2r_{1}r_{2}},\;\;\;\beta(\Delta)\equiv\theta(1-\Delta)\theta(1+\Delta),$
$\displaystyle\theta(y)=0,\;y<0;\;\;\;1/2,\;y=0;\;\;\;1,\;y>0.$ (70)

We note that the Legendre polynomial on the right-hand side does not
explicitly depend upon the individual unit vectors $\hat{r}_{1},\hat{r}_{2},$
or $\hat{r}_{3}$; its argument is actually just the cosine of the angle
between $r_{1}$ and $r_{2}$ in a triangle with sides $r_{1},r_{2}$, and
$r_{3}$. We also note that the $2\times 3$ matrix in curly brackets is a
Wigner 6-$j$ symbol (e.g. NIST DLMF 34.4999https://dlmf.nist.gov/34.4),
defined as a sum over spins of a specified product of 3-$j$ symbols.

###  6.2 Delta Function for $N=4$: Naive Approach

For $N=4$ and $N=5$, the treatment of the intermediates requires additional
consideration. However, an interesting observation is that, for the Delta
function derivation, we may begin with the “naive” (incorrect) $N=4$
generating function equation (35) and still obtain a correct expansion. We
first do this, and then begin with the proper generating function, showing
that, interestingly, using these two approaches enables expanding slightly
different versions of the Delta function. Proceeding first from equation (35),
we have, with
$\boldsymbol{R}\equiv\boldsymbol{r}_{1}+\boldsymbol{r}_{2}+\boldsymbol{r}_{3}+\boldsymbol{r}_{4}$
and $R\equiv|\boldsymbol{R}|$ there and here,

$\displaystyle\int\frac{p^{2}dp}{2\pi^{2}}\,j_{0}(pR)=(4\pi)^{2}\sum_{\ell_{1},\ldots,\ell_{4}}(-1)^{(\ell_{1}+\ell_{2}+\ell_{3}+\ell_{4})/2}\,\mathcal{D}^{\rm
P}\sum_{\ell_{12}}\sqrt{2\ell_{12}+1}\begin{pmatrix}\ell_{1}&\ell_{2}&\ell_{12}\\\
0&0&0\end{pmatrix}\begin{pmatrix}\ell_{12}&\ell_{3}&\ell_{4}\\\
0&0&0\end{pmatrix}$
$\displaystyle\times(-1)^{\ell_{12}}\mathcal{P}_{\ell_{1}\ell_{2}(\ell_{12})\ell_{3}\ell_{4}}(\hat{r}_{1},\hat{r}_{2},\hat{r}_{3},\hat{r}_{4})\int\frac{p^{2}dp}{2\pi^{2}}\,j_{\ell_{1}}(pr_{1})j_{\ell_{2}}(pr_{2})j_{\ell_{3}}(pr_{3})j_{\ell_{4}}(pr_{4}).$
(71)

We simplified the phase since the two 3-$j$ symbols imply that the sum of the
primary $\ell_{i}$ is even. If we denote the quadruple-sBF overlap integral on
the right-hand side (including the factor of $1/(2\pi^{2})$) as $I_{4}$, it
may be evaluated using a technique presented by [18], §6. One rewrites the
integral as a double integral where the outer integrand is a pair of triple-
sBF integrals; this approach relies on the fact that a double-sBF overlap
integral with the appropriate weight and at matching orders gives a Dirac
Delta function, as

$\displaystyle\int
x^{2}dx\;j_{\ell}(ax)j_{\ell}(bx)=\frac{\pi}{2ab}\,\delta_{\rm D}^{[1]}(a-b),$
(72)

e.g. NIST DLMF 1.17.14.101010https://dlmf.nist.gov/1.17

We have

$\displaystyle I_{4}=\frac{2}{\pi}\cdot\frac{1}{2\pi^{2}}\int
q^{2}dq\left[\int
p^{2}dp\;j_{\ell_{1}}(pr_{1})j_{\ell_{2}}(pr_{2})j_{L}(qp)\int p^{\prime
2}dp^{\prime}\;j_{\ell_{3}}(p^{\prime}r_{3})j_{\ell_{4}}(p^{\prime}r_{4})j_{L}(qp^{\prime})\right].$
(73)

The question now becomes whether we can always evaluate the two triple-sBF
integrals required above. Examining the result of [17] for each of the three-
sBf integrals $I$ we require, we see that we may as long as $L$ is chosen so
as to be in a triangular relation with both $\ell_{1},\ell_{2}$ and
$\ell_{3},\ell_{4}$. If we choose that $L\equiv\ell_{12}$, this will always be
guaranteed, and indeed provides an interesting connection to the idea of
coupling two angular momenta into an intermediate; this now happens on the
orders of the sBFs. Making this choice, we find that

$\displaystyle I_{4}=\frac{2}{\pi}\cdot\frac{1}{2\pi^{2}}\int
q^{2}dq\left[I(\ell_{1},\ell_{2},\ell_{12};r_{1},r_{2},q)\;I(\ell_{12},\ell_{3},\ell_{4};q,r_{3},r_{4})\right]\equiv\pi^{-3}\,\bar{I}_{4}.$
(74)

Hence, the Dirac Delta function may be expanded in the isotropic basis as

$\displaystyle\delta_{\rm D}^{[3]}(|{\bf r}_{1}+{\bf r}_{2}+{\bf r}_{3}+{\bf
r}_{4}|)=\sum_{\ell_{i}}c_{\ell_{1}\ell_{2}\ell_{3}\ell_{4}}(r_{1},r_{2},r_{3},r_{4})$
$\displaystyle\qquad\qquad\times\sum_{\ell_{12}}(-1)^{\ell_{12}}\sqrt{2\ell_{12}+1}\begin{pmatrix}\ell_{1}&\ell_{2}&\ell_{12}\\\
0&0&0\end{pmatrix}\begin{pmatrix}\ell_{12}&\ell_{3}&\ell_{4}\\\
0&0&0\end{pmatrix}\mathcal{P}_{\ell_{1}\ell_{2}(\ell_{12})\ell_{3}\ell_{4}}(\hat{r}_{1},\hat{r}_{2},\hat{r}_{3},\hat{r}_{4})$
(75)

with

$\displaystyle
c_{\ell_{1}\ell_{2}\ell_{3}\ell_{4}}(r_{1},r_{2},r_{3},r_{4})\equiv\frac{16}{\pi}\,(-1)^{(\ell_{1}+\ell_{2}+\ell_{3}+\ell_{4})/2}\,\mathcal{D}^{\rm
P}\,\bar{I}_{4}(r_{1},r_{2},r_{3},r_{4}),$ (76)

where $\mathcal{D}^{\rm P}$ is a product of $\sqrt{2\ell_{i}+1}$ over all
primary angular momenta.

It is notable that the radial expansion coefficients do not depend on the
intermediate angular momentum, $\ell_{12}$, which is summed over at fixed
$\ell_{i}$ in our expansion. The symmetry of the Delta function with respect
to its four vector arguments $\bf{r}_{i}$ suggests that any expansion for it
should be “symmetrized”, or summed, over all possible couplings of the four
primary angular momenta that pertain to them, and this is indeed what occurs.
This summation also means that we may define “recoupling-symmetrized”
isotropic basis functions by the sum in the second line of equation (75);
these are a subset of the full isotropic basis which is suitable for expanding
functions that have global interchange symmetry among all arguments, as the
Dirac Delta function does. Ultimately, the irrelevance of the intermediate is
why we can obtain a correct expansion of the Delta function even using the
naive $N=4$ generating function approach (which does not allow control of the
intermediate).

###  6.3 Delta Function for $N=4$: Approach with the Correct Generating
Function

We now show that we obtain the same result as in the previous section if we
begin with the correct generating function at $N=4$, which is given by
equation (43).

Here, it is informative to begin with the Dirac delta function itself. We
first write the desired 3D Dirac delta function as the inverse Fourier
Transform of unity:

$\displaystyle\delta_{\rm D}^{[3]}({\bf R}+{\bf
R^{\prime}})=\int\frac{d^{3}\boldsymbol{p}}{(2\pi)^{3}}\;e^{-i\boldsymbol{p}\cdot\boldsymbol{R}}e^{-i\boldsymbol{p}\cdot\boldsymbol{R}^{\prime}}.$
(77)

Expanding the plane waves into spherical harmonics and spherical Bessel
functions, integrating over $d\Omega_{p}$, invoking orthogonality, and summing
over spins, we obtain

$\displaystyle\delta_{\rm D}^{[3]}({\bf R}+{\bf
R^{\prime}})=\sum_{\ell}(-1)^{\ell}(2\ell+1)\mathcal{L}_{\ell}(\hat{R}\cdot\hat{R}^{\prime})\int\frac{p^{2}dp}{2\pi^{2}}\;j_{\ell}(pR)j_{\ell}(pR^{\prime}),$
(78)

If we integrate equation (43) over $u_{12}^{2}du_{12}$ we will obtain a result
resembling our desired $p$ integrand in equation (78). In particular, we get
from the $u_{12}$ integration:

$\displaystyle\int
u_{12}^{2}du_{12}\;H_{4}=\sum_{\ell}(-1)^{\ell}(2\ell+1)\mathcal{L}_{\ell}(\hat{R}\cdot\hat{R}^{\prime})\,j_{\ell}(pR)j_{\ell}(p^{\prime}R^{\prime})\,\frac{\pi}{2p^{2}}\delta_{\rm
D}^{[1]}(p-p^{\prime}).$ (79)

The Delta function will set $p=p^{\prime}$, so if we integrate over $p$ we can
come close to equation (78); however we do not have the right pre-factor. We
thus need to multiply by $p^{4}/\pi^{3}$ before integrating. We then have

$\displaystyle\pi^{-3}\int
p^{4}dp\,u_{12}^{2}du_{12}\;H_{4}=\sum_{\ell}(-1)^{\ell}(2\ell+1)\mathcal{L}_{\ell}(\hat{R}\cdot\hat{R}^{\prime})\,\int\frac{p^{2}dp}{2\pi^{2}}\;j_{\ell}(pR)j_{\ell}(p^{\prime}R^{\prime})=\delta_{\rm
D}^{[3]}({\bf R}+{\bf R^{\prime}}),$ (80)

where the last equality follows from equation (78).

Now, using equation (43) to replace $H_{4}$ with the expansion we have in four
spherical Bessel functions, we find that

$\displaystyle\delta_{\rm D}^{[3]}({\bf R}+{\bf R^{\prime}})$
$\displaystyle=\frac{16}{\pi}\sum_{\ell_{1}\ell_{2}\ell_{12}}\sum_{\ell_{3}\ell_{4}}(-1)^{(\ell_{1}+\ell_{2}+\ell_{3}+\ell_{4})/2}\,\mathcal{D}^{\rm
P}\sqrt{2\ell_{12}+1}$ $\displaystyle\times\int
p^{4}dp\,u_{12}^{2}du_{12}\;j_{\ell_{1}}(pr_{1})j_{\ell_{2}}(pr_{2})j_{\ell_{3}}(p^{\prime}r_{3})j_{\ell_{4}}(p^{\prime}r_{4})j_{\ell_{12}}(pu_{12})j_{\ell_{12}}(p^{\prime}u_{12})$
$\displaystyle\qquad\times\begin{pmatrix}\ell_{1}&\ell_{2}&\ell_{12}\\\
0&0&0\end{pmatrix}\begin{pmatrix}\ell_{12}&\ell_{3}&\ell_{4}\\\
0&0&0\end{pmatrix}\mathcal{P}_{\ell_{1}\ell_{2}(\ell_{12})\ell_{3}{\ell_{4}}}(\hat{r}_{1},\hat{r}_{2},\hat{r}_{3},\hat{r}_{4}).$
(81)

We have used the remark below equation (43), that the sum
$\ell_{1}+\ell_{2}+\ell_{3}+\ell_{4}$ is even, to remove that phase factor,
relative to equation (43).

Now, performing the $u_{12}$ integration we find

$\displaystyle\delta_{\rm D}^{[3]}({\bf R}+{\bf R^{\prime}})$
$\displaystyle=\frac{16}{\pi}\sum_{\ell_{1}\ell_{2}}\sum_{\ell_{3}\ell_{4}}(-1)^{(\ell_{1}+\ell_{2}+\ell_{3}+\ell_{4})/2}\,\mathcal{D}^{\rm
P}\int
p^{2}dp\;j_{\ell_{1}}(pr_{1})j_{\ell_{2}}(pr_{2})j_{\ell_{3}}(pr_{3})j_{\ell_{4}}(pr_{4})$
$\displaystyle\qquad\times\sum_{\ell_{12}}\sqrt{2\ell_{12}+1}\begin{pmatrix}\ell_{1}&\ell_{2}&\ell_{12}\\\
0&0&0\end{pmatrix}\begin{pmatrix}\ell_{12}&\ell_{3}&\ell_{4}\\\
0&0&0\end{pmatrix}\mathcal{P}_{\ell_{1}\ell_{2}(\ell_{12})\ell_{3}{\ell_{4}}}(\hat{r}_{1},\hat{r}_{2},\hat{r}_{3},\hat{r}_{4}).$
(82)

$\mathcal{D}^{\rm P}$ is defined in the paragraph below equation (35).

It is again notable that the intermediate $\ell_{12}$ here is simply summed
over, with no effect on the spherical Bessel overlap integral on the right-
hand side; this is because the Delta function is totally symmetric under
interchange of its four arguments.

###  6.4 Delta Function for $N=5$

With the two subsections treating the $N=4$ case, we have already displayed
all the important conceptual components involved in isotropic basis expansions
of the Delta function. For completeness, we now present the results for $N=5$.
We first integrate the left-hand side of equation (57) against
$u_{123}^{2}du_{123}$, obtaining a spherically-symmetric Delta function
setting $p^{\prime\prime}=p^{\prime}$. We then integrate over
$u_{12}^{2}du_{12}$, which gives a spherically symmetric Delta function
setting $p^{\prime}=p$. Overall we find that

$\displaystyle\int u_{12}^{2}du_{12}\int u_{123}^{2}du_{123}\;P_{5}$
$\displaystyle=(4\pi)^{7}\left(\frac{\pi}{2}\right)^{2}\frac{1}{p^{4}}\sum_{\rm
all}(-1)^{\ell+\ell^{\prime}+(\ell+\ell^{\prime}+\ell^{\prime\prime})/2}j_{\ell}(pR)j_{\ell^{\prime}}(pR^{\prime})j_{\ell^{\prime\prime}}(pR^{\prime\prime\prime})$
$\displaystyle\qquad\times Y_{\ell
m}(\hat{R})Y_{\ell^{\prime}m^{\prime}}(\hat{R}^{\prime})Y_{\ell^{\prime\prime}m^{\prime\prime}}(\hat{R}^{\prime\prime})\mathcal{G}_{\ell^{\prime}\ell\ell^{\prime\prime}}^{m^{\prime}mm^{\prime\prime}}.$
(83)

We used the remark below equation (57), that the sum
$\ell+\ell^{\prime}+\ell^{\prime\prime}$ is even, to simplify a power of $i$
into a power of $(-1)$ above. Now, we note that if the parity of
$\ell+\ell^{\prime}$ is equal to that of
$(\ell+\ell^{\prime}+\ell^{\prime\prime})/2$, the phase above will vanish. If
$\ell+\ell^{\prime}$ is odd, then so is $\ell^{\prime\prime}$ (by the zero-
spin 3-$j$ symbol) and hence $(\ell+\ell^{\prime}+\ell^{\prime\prime})/2$ is
odd, so the phase vanishes. Now, if $\ell+\ell^{\prime}$ is even, so is
$\ell^{\prime}$ (again by the zero-spin 3-$j$ symbol), and hence
$(\ell+\ell^{\prime}+\ell^{\prime\prime})/2$ is even, so again the phase
vanishes. With this simplification in hand, and summing over the spins in the
spherical harmonics times the Gaunt integral, we can see that the above
relation is simply equation (68) (the relation for the $N=3$ delta function)
prior to the integration of that latter over $p^{2}dp/(2\pi^{2})$ and also
multiplied by $(4\pi)^{5}(\pi/2)^{2}p^{-4}$ and with the replacements
$\hat{r}_{1}\to\hat{R}$, $\hat{r}_{2}\to\hat{R}^{\prime}$, and
$\hat{r}_{3}\to\hat{R}^{\prime\prime}$. Hence, we see that integrating
equation (83) against $p^{6}dp/(2\pi^{2})$ and dividing by
$(4\pi)^{5}(\pi/2)^{2}$ will give a Dirac Delta function $\delta_{\rm
D}^{[3]}(\boldsymbol{R}+\boldsymbol{R}^{\prime}+\boldsymbol{R}^{\prime\prime})$.

Hence, we may write that

$\displaystyle\delta_{\rm
D}^{[3]}(\boldsymbol{R}+\boldsymbol{R}^{\prime}+\boldsymbol{R}^{\prime\prime})=\frac{1}{(4\pi)^{5}(\pi/2)^{2}}\int\frac{p^{6}dp}{2\pi^{2}}\int
u_{12}^{2}du_{12}\int u_{123}^{2}du_{123}\;P_{5}.$ (84)

Having thus shown that applying the operations outlined above to the
generating function for the $N=5$ isotropic basis functions indeed yields the
desired Delta function, we now need to apply the same operations to the right-
hand side of equation (62) to obtain the Delta function’s expansion in the
isotropic basis. We see that the $u_{123}$ integration will produce a Dirac
Delta function setting $p^{\prime\prime}=p^{\prime}$, and then the $u_{12}$
integration will set $p^{\prime}=p$. So we have overall a factor of $p^{-4}$
to deal with. Integrating against $p^{6}dp/(2\pi^{2})$, we find

$\displaystyle\delta_{\rm
D}^{[3]}(\boldsymbol{R}+\boldsymbol{R}^{\prime}+\boldsymbol{R}^{\prime\prime})=(4\pi)^{4}\sum_{\rm
all}(-1)^{(\ell_{1}+\ell_{2}+\ell_{3}+\ell_{4}+\ell_{5})/2}\int\frac{p^{2}dp}{2\pi^{2}}\;j_{\ell_{1}}(pr_{1})j_{\ell_{2}}(pr_{2})j_{\ell_{3}}(pr_{3})j_{\ell_{4}}(pr_{4})j_{\ell_{5}}(pr_{5})$
$\displaystyle\qquad\times\begin{pmatrix}\ell_{1}&\ell_{2}&\ell_{12}\\\
0&0&0\end{pmatrix}\begin{pmatrix}\ell_{12}&\ell_{3}&\ell_{123}\\\
0&0&0\end{pmatrix}\begin{pmatrix}\ell_{123}&\ell_{4}&\ell_{5}\\\
0&0&0\end{pmatrix}\mathcal{D}^{\rm P}\sqrt{2\ell_{12}+1}\sqrt{2\ell_{123}+1}$
$\displaystyle\qquad\times\mathcal{P}_{\ell_{1}\ell_{2}(\ell_{12})\ell_{3}(\ell_{123})\ell_{4}\ell_{5}}(\hat{r}_{1},\hat{r}_{2},\hat{r}_{3},\hat{r}_{4},\hat{r}_{5}).$
(85)

We notice that the parity of $\ell_{1}+\ell_{2}+\ell_{3}+\ell_{4}+\ell_{5}$ is
even, as can be shown by manipulating the rules implied by the zero-spin 3-$j$
symbols and canceling the parities of the intermediates.

##  7 Connection to Overlap Integrals of Multiple Spherical Bessel Functions

The foregoing work with the Dirac Delta function already suggests that there
is a connection between the isotropic basis functions and overlap integrals of
multiple sBFs; here we explore this connection further. A complete discussion
of previous works on overlap integrals of spherical Bessel functions is
presented in [19]; here we focus only on the results of immediate use in the
present work.

###  7.1 $N=2$

We consider integrating both sides of our fundamental relation for the $N=2$
functions. We have

$\displaystyle\int
dp\;j_{0}(pR)=\frac{\pi}{2}\frac{1}{R}=(4\pi)^{3/2}\sum_{\ell}\sqrt{2\ell+1}\;\mathcal{P}_{\ell\ell}(\hat{r}_{1},\hat{r}_{2})\int
dp\;j_{\ell}(pr_{1})j_{\ell}(pr_{2}).$ (86)

The first equality is simply from changing variables and using the result for
the Sine integral.111111https://mathworld.wolfram.com/SineIntegral.html By
comparing $(\pi/2)(1/R)$ to the usual expansion of a two-argument denominator
in Legendre polynomials, we may conclude that the two-sBF overlap integral on
the right-hand side is $(\pi/2)r_{<}^{\ell}/r_{>}^{\ell+1}$, where $r_{<}$ is
the lesser of $r_{1}$ and $r_{2}$ and $r_{>}$ the greater. This recovers the
standard result for this sBF integral (e.g. [20] equation C8 and references
below).

###  7.2 $N=3$

We now consider the same method applied to the $N=3$ functions. We have

$\displaystyle\int dp\;j_{0}(pR)=\frac{\pi}{2}\frac{1}{R}$
$\displaystyle=(4\pi)^{2}\sum_{\ell_{1}\ell_{2}\ell_{3}}\begin{pmatrix}\ell_{1}&\ell_{2}&\ell_{3}\\\
0&0&0\end{pmatrix}\mathcal{C}_{\ell_{1}\ell_{2}\ell_{3}}(-1)^{(\ell_{1}+\ell_{2}+\ell_{3})/2}\,\mathcal{P}_{\ell_{1}\ell_{2}\ell_{3}}(\hat{r}_{1},\hat{r}_{2},\hat{r}_{3})$
$\displaystyle\qquad\times\int
dp\;j_{\ell_{1}}(pr_{1})j_{\ell_{2}}(pr_{2})j_{\ell_{3}}(pr_{3}),$ (87)

where we again used the Sine integral to evaluate the left-hand side. The
right-hand side above is simply from equation (27), and we recall that
$\mathcal{C}_{\ell_{1}\ell_{2}\ell_{3}}$ is defined below equation (8).

We now integrate both sides against
$\mathcal{P}^{*}_{\ell_{1}\ell_{2}\ell_{3}}$. By orthonormality of the
$\mathcal{P}_{\Lambda}$, this reduces the right-hand side to simply the
triple-sBF overlap integral (with a constant pre-factor).

Now, to obtain an explicit result in terms of $r_{1},r_{2}$ and $r_{3}$ we
consider expanding the $1/R$ on the left-hand side. We have

$\displaystyle\frac{1}{R}$
$\displaystyle=\frac{1}{|\boldsymbol{r}_{1}+\boldsymbol{r}_{2}+\boldsymbol{r}_{3}|}=I_{0}^{0}(\boldsymbol{r}_{12}+\boldsymbol{r}_{3})$
$\displaystyle=\sum_{\lambda=0}^{\infty}{{2\lambda+1}\choose{2\lambda}}^{1/2}\sum_{\mu=-\lambda}^{\lambda}R_{\lambda}^{\mu}(\boldsymbol{r}_{12})I_{\lambda}^{-\mu}(\boldsymbol{r}_{3})\left<\lambda,\mu;\lambda,-\mu|00\right>.$
(88)

We have defined
$\boldsymbol{r}_{12}\equiv\boldsymbol{r}_{1}+\boldsymbol{r}_{2}$.
$I_{\ell}^{m}$ is the irregular solid harmonic,

$\displaystyle I_{\ell}^{m}(\boldsymbol{R})$
$\displaystyle\equiv\sqrt{\frac{4\pi}{2\ell+1}}\,\frac{Y_{\ell
m}(\hat{R})}{R^{\ell+1}};$ (89) $\displaystyle I_{0}^{0}(\boldsymbol{R})$
$\displaystyle=\sqrt{4\pi}\,\frac{Y_{00}}{R}=\frac{1}{R}.$

and $R_{\lambda}^{\mu}$ is the regular solid harmonic.

On the second line of equation (7.2), the first factor within the sum is a
binomial coefficient, and the final factor, in angle brackets, is a Clebsch-
Gordan coefficient, standardly related to the 3-$j$ symbol121212E.g.
https://mathworld.wolfram.com/Clebsch-GordanCoefficient.html (but more
convenient for our intermediate steps involving the shift theorem). For
clarity, we use a semicolon to separate the orbital angular momentum-spin pair
for the first and second of the two being added into the third, which appears
after the bar.

To obtain the second line of equation (7.2), we used the addition theorem for
irregular solid harmonics. The result we have written requires that
$r_{12}\leq r_{3}$. If $r_{12}>r_{3}$, one would interchange ${\bf
r}_{12}\leftrightarrow{\bf r}_{3}$ in the second line. We now use the addition
theorem for regular solid harmonics to expand
$R_{\lambda}^{\mu}(\boldsymbol{r}_{12})$; this imposes no restriction on the
magnitudes of $r_{1}$ and $r_{2}$. We have

$\displaystyle\frac{1}{R}=\sum_{\lambda\mu}(-1)^{\lambda-\mu}I_{\lambda}^{-\mu}(\boldsymbol{r}_{3})\sum_{\lambda^{\prime}=0}^{\lambda}{2\lambda\choose
2\lambda^{\prime}}^{1/2}\sum_{\mu^{\prime}=-\lambda^{\prime}}^{\lambda^{\prime}}R_{\lambda^{\prime}}^{\mu^{\prime}}(\boldsymbol{r}_{1})R_{\lambda-\lambda^{\prime}}^{\mu-\mu^{\prime}}(\boldsymbol{r}_{2})\left<\lambda^{\prime},\mu^{\prime};\lambda-\lambda^{\prime},\mu-\mu^{\prime}|\lambda\mu\right>.$
(90)

We used the relationship between Clebsch-Gordan coefficients and 3-$j$
symbols,131313E.g. https://mathworld.wolfram.com/Clebsch-
GordanCoefficient.html equation 7. and then NIST DLMF 34.3.1, to simplify the
zero-total angular momentum Clebsch-Gordan coefficient, and we also simplified
the binomial coefficient.

Rewriting the solid harmonics explicitly using spherical harmonics and powers
of the $r_{i}$, we find

$\displaystyle\frac{1}{R}$
$\displaystyle=4\pi\sum_{\lambda\mu,\lambda^{\prime}\mu^{\prime}}\mathcal{C}_{\lambda,\lambda^{\prime},\lambda-\lambda^{\prime}}^{-1}{{2\lambda}\choose{2\lambda^{\prime}}}^{1/2}\begin{pmatrix}\lambda^{\prime}&\lambda-\lambda^{\prime}&\lambda\\\
\mu^{\prime}&\mu-\mu^{\prime}&-\mu\end{pmatrix}r_{3}^{-\lambda-1}\,r_{1}^{\lambda^{\prime}}\,r_{2}^{\lambda-\lambda^{\prime}}$
(91) $\displaystyle\qquad\qquad\qquad\qquad\times
Y_{\lambda,-\mu}(\hat{r}_{3})Y_{\lambda^{\prime}\mu^{\prime}}(\hat{r}_{1})Y_{\lambda-\lambda^{\prime},\mu-\mu^{\prime}}(\hat{r}_{2}).$

We recall that $r_{12}\leq r_{3}$ here. We note that our original integral has
no preference for $\ell_{1}$, $\ell_{2}$, or $\ell_{3}$. We have intentionally
not simplified the above or canceled constant factors so that it is
transparent whence each piece comes. Of course these simplifications can be
performed at the final step.

Now, following our prescription below equation (7.2), we integrate equation
(91) against $\mathcal{P}^{*}_{\ell_{1}\ell_{2}\ell_{3}}$ will set
$\lambda^{\prime},\mu^{\prime}=\ell_{1},m_{1}$,
$\lambda-\lambda^{\prime},\mu-\mu^{\prime}=\ell_{2},m_{2}$, and
$\lambda,-\mu=\ell_{3},m_{3}$ and eliminate the sums.

We may now sum over the spins using NIST DLMF 34.3.18, and find

$\displaystyle\int
d\hat{r}_{1}d\hat{r}_{2}d\hat{r}_{3}\;\frac{1}{R}\;\mathcal{P}_{\ell_{1}\ell_{2}\ell_{3}}^{*}(\hat{r}_{1},\hat{r}_{2},\hat{r}_{3})=4\pi\,\mathcal{C}_{\ell_{1}\ell_{2}\ell_{3}}^{-1}{2\ell_{3}\choose
2\ell_{1}}\,r_{3}^{-\ell_{3}-1}\,r_{1}^{\ell_{1}}\,r_{2}^{\ell_{2}}.$ (92)

We highlight that the integral result we have obtained above is simply the
expansion coefficient of $1/R$ in the isotropic basis functions of three
arguments.

We also note that one could pursue the same steps if one had split the problem
as $\boldsymbol{r}_{13}\equiv\boldsymbol{r}_{1}+\boldsymbol{r}_{3}$ and
$\boldsymbol{r}_{2}$, or even
$\boldsymbol{r}_{23}\equiv\boldsymbol{r}_{2}+\boldsymbol{r}_{3}$ and
$\boldsymbol{r}_{1}$ and one would find equivalent answers. They will pertain
respectively to regions on which $r_{3}>r_{12}$, $r_{2}>r_{13}$ and
$r_{1}>r_{23}$. Since there is full interchange symmetry among the vectors
$\boldsymbol{r}_{i}$, we can simply choose an ordering and so will only need
one result; we may therefore adopt the one explicitly obtained above.

Setting equation (92) (times $\pi/2$) equal to equation (7.2) and solving for
our desired integral overlap integral of the three sBFs, which we denote
$I_{\ell_{1}\ell_{2}\ell_{3}}^{[0]}$, where the zero indicates the power of
$p$ weighting the integrand, we find

$\displaystyle I_{\ell_{1}\ell_{2}\ell_{3}}^{[0]}(r_{1},r_{2},r_{3})$
$\displaystyle=\frac{1}{8}\,\mathcal{C}_{\ell_{1}\ell_{2}\ell_{3}}^{-2}\begin{pmatrix}\ell_{1}&\ell_{2}&\ell_{3}\\\
0&0&0\end{pmatrix}^{-1}{2\ell_{3}\choose
2\ell_{1}}\,(-1)^{(\ell_{1}+\ell_{2}+\ell_{3})/2}\,r_{3}^{-\ell_{3}-1}\,r_{1}^{\ell_{1}}\,r_{2}^{\ell_{2}},\;\;\;r_{12}<r_{3},$
(93)

and where there is no constraint between $r_{1}$ and $r_{2}$. This result is
analogous to that of [17] but with a weight in the integral of $p^{0}$ rather
than his case which would have $p^{2}$ in our notation; discussion of this
latter type of integral from a different perspective is also in [21].

###  7.3 $N=4$

We focus on equation (43). We first compute the integral of the left-hand side
over $u_{12}^{2}du_{12}$; this integration yields a factor of
$\pi/(2pp^{\prime})\delta_{\rm D}^{[1]}(p-p^{\prime})$. Integrating over
$p^{\prime 2}dp^{\prime}$ what remains on the left-hand side against this
Dirac Delta function factor, we then find

$\displaystyle\mathcal{J}_{4}\equiv\int p^{\prime 2}dp^{\prime}\;\int
u_{12}^{2}du_{12}\;H_{4}=\frac{\pi}{2}\sum_{\ell}(2\ell+1)\mathcal{L}_{\ell}(\hat{R}\cdot\hat{R}^{\prime})j_{\ell}(pR)j_{\ell}(pR^{\prime}).$
(94)

Finally, integrating over $p\,dp$, we obtain

$\displaystyle\bar{\mathcal{J}}_{4}\equiv\int p\,dp\;\mathcal{J}_{4}$
$\displaystyle=\frac{\pi}{2}\sum_{\ell}(2\ell+1)\mathcal{L}_{\ell}(\hat{R}\cdot\hat{R}^{\prime})\int
dp\;j_{\ell}(pR)j_{\ell}(pR^{\prime})=\left(\frac{\pi}{2}\right)^{2}\sum_{\ell}\mathcal{L}_{\ell}(\hat{R}\cdot\hat{R}^{\prime})\frac{R_{<}^{\ell}}{R_{>}^{\ell+1}}$
$\displaystyle=\left(\frac{\pi}{2}\right)^{2}\frac{1}{|\boldsymbol{r}_{1}+\boldsymbol{r}_{2}+\boldsymbol{r}_{3}+\boldsymbol{r}_{4}|}.$
(95)

For the integral over $p\,dp$, we used the result equation (C8) from
[20].141414The integral there has a typo; each $L_{1}$ inside the integral
should be $L$. We now observe that the fraction
$1/|\boldsymbol{r}_{1}+\cdots+\boldsymbol{r}_{4}|$ is simply the irregular
solid harmonic of order and degree zero,
$I_{0}^{0}(\boldsymbol{R}+\boldsymbol{R}^{\prime})$. We may then apply to it
the shift theorem for the irregular solid harmonic; with $R\leq R^{\prime}$,
we find

$\displaystyle\bar{\mathcal{J}}_{4}=\left(\frac{\pi}{2}\right)^{2}\sum_{\lambda\mu}\mathcal{D}^{+}_{0,\lambda,0,\mu}\;R_{\lambda}^{\mu}(\hat{R})I_{\lambda}^{-\mu}(\hat{R}^{\prime}).$
(96)

We have defined

$\displaystyle\mathcal{D}_{\lambda_{1},\lambda_{2},\mu_{1},\mu_{2}}^{+}$
$\displaystyle\equiv{{2(\lambda_{1}+\lambda_{2})+1}\choose{2\lambda_{2}}}^{1/2}\left<\lambda_{2},\mu_{2};\lambda_{1}+\lambda_{2},\mu_{1}-\mu_{2}|\lambda_{1}\mu_{1}\right>,$
(97)

as it will be of use in what follows; in equation (96) we have the limit where
$\lambda_{1}=0,\lambda_{2}=\lambda,\mu_{1}=0,\mu_{2}=\mu$.

Each solid harmonic on the right-hand side of equation (96) can be expanded
using the appropriate shift theorems (there is one for the regular solid
harmonic, $R_{\lambda}^{\mu}$, and another for the irregular solid harmonic,
of which we already once made use). We require $r_{3}\leq r_{4}$ for this
latter expansion. We find

$\displaystyle\bar{\mathcal{J}}_{4}=$
$\displaystyle\left(\frac{\pi}{2}\right)^{2}\sum_{\lambda\mu}\mathcal{D}_{0,\lambda,0,\mu}^{+}\sum_{\lambda^{\prime}\mu^{\prime}}\mathcal{D}_{\lambda,\lambda^{\prime},\mu,\mu^{\prime}}\sum_{\lambda^{\prime\prime}\mu^{\prime\prime}}\mathcal{D}^{+}_{\lambda,\lambda^{\prime\prime},-\mu,\mu^{\prime\prime}}$
$\displaystyle\times(4\pi)^{2}\left[(2\lambda^{\prime}+1)(2(\lambda-\lambda^{\prime})+1)(2\lambda^{\prime\prime}+1)(2(\lambda^{\prime}+\lambda^{\prime\prime})+1)\right]^{1/2}r_{1}^{\lambda^{\prime}}r_{2}^{\lambda-\lambda^{\prime}}r_{3}^{\lambda^{\prime\prime}}r_{4}^{-\lambda-\lambda^{\prime\prime}-1}$
$\displaystyle\times
Y_{\lambda^{\prime}\mu^{\prime}}(\hat{r}_{1})Y_{\lambda-\lambda^{\prime},\mu-\mu^{\prime}}(\hat{r}_{2})Y_{\lambda^{\prime\prime},\mu^{\prime\prime}}(\hat{r}_{3})Y_{\lambda+\lambda^{\prime\prime},-\mu-\mu^{\prime\prime}}(\hat{r}_{4}).$
(98)

The constants in the middle line stem from the normalizations of the four
solid harmonics that were ultimately involved. We have also defined

$\displaystyle\mathcal{D}_{\lambda_{1},\lambda_{2},\mu_{1},\mu_{2}}$
$\displaystyle\equiv{{2\lambda_{1}}\choose{2\lambda_{2}}}^{1/2}\left<\lambda_{2},\mu_{2};\lambda_{1}-\lambda_{2},\mu_{1}-\mu_{2}|\lambda_{1}\mu_{1}\right>;$
(99)

we note it is different from $\mathcal{D}^{+}$ as defined in equation (96).

We now integrate the right-hand side of equation (43) against
$u_{12}^{2}du_{12}$ and then against $p^{\prime 2}dp^{\prime}$, just as we did
for the left-hand side above. We then integrate over $p\,dp$, finding

$\displaystyle\bar{\mathcal{J}}_{4}$
$\displaystyle=(4\pi)^{2}\left(\frac{\pi}{2}\right)\int
p\,dp\;j_{\ell_{1}}(pr_{1})j_{\ell_{2}}(pr_{2})j_{\ell_{3}}(pr_{3})j_{\ell_{4}}(pr_{4})\mathcal{D}^{\rm
P}\sqrt{2\ell_{12}+1}$
$\displaystyle\qquad\times(-1)^{(\ell_{1}+\ell_{2}+\ell_{3}+\ell_{4})/2}\begin{pmatrix}\ell_{1}&\ell_{2}&\ell_{12}\\\
0&0&0\end{pmatrix}\begin{pmatrix}\ell_{12}&\ell_{3}&\ell_{4}\\\
0&0&0\end{pmatrix}\mathcal{P}_{\ell_{1}\ell_{2}(\ell_{12})\ell_{3}\ell_{4}}(\hat{r}_{1},\hat{r}_{2},\hat{r}_{3},\hat{r}_{4}).$
(100)

Setting equations (100) and (98) equal and then integrating both sides against
$\mathcal{P}_{\ell_{1}\ell_{2}(\ell_{12})\ell_{3}\ell_{4}}^{*}$ sets
$\lambda^{\prime}=\ell_{1},\;\lambda-\lambda^{\prime}=\ell_{2},\;\lambda^{\prime\prime}=\ell_{3}$
and $\lambda+\lambda^{\prime\prime}=\ell_{4}$, which, taken together, imply
that $\ell_{1}+\ell_{2}+\ell_{3}=\ell_{4}$. This gives an expression for the
overlap integral of four spherical Bessel functions, weighted by $p$, in terms
of a sum over powers of the free arguments $r_{i}$. Explicitly, we obtain

$\displaystyle\int
p\,dp\;j_{\ell_{1}}(pr_{1})j_{\ell_{2}}(pr_{2})j_{\ell_{3}}(pr_{3})j_{\ell_{4}}(pr_{4})=\left[\mathcal{D}^{\rm
P}\sqrt{2\ell_{12}+1}\right]^{-1}\begin{pmatrix}\ell_{1}&\ell_{2}&\ell_{12}\\\
0&0&0\end{pmatrix}^{-1}\begin{pmatrix}\ell_{12}&\ell_{3}&\ell_{4}\\\
0&0&0\end{pmatrix}^{-1}$
$\displaystyle\qquad\times\frac{\pi}{2}\sum_{\lambda\mu}\sum_{\lambda^{\prime}\mu^{\prime}}\sum_{\lambda^{\prime\prime}\mu^{\prime\prime}}\mathcal{D}_{0,\lambda,0,\mu}^{+}\mathcal{D}_{\lambda,\lambda^{\prime},\mu,\mu^{\prime}}\mathcal{D}^{+}_{\lambda,\lambda^{\prime\prime},-\mu,\mu^{\prime\prime}}$
$\displaystyle\qquad\times\left[(2\lambda^{\prime}+1)(2(\lambda-\lambda^{\prime})+1)(2\lambda^{\prime\prime}+1)(2(\lambda^{\prime}+\lambda^{\prime\prime})+1)\right]^{1/2}r_{1}^{\lambda^{\prime}}r_{2}^{\lambda-\lambda^{\prime}}r_{3}^{\lambda^{\prime\prime}}r_{4}^{-\lambda-\lambda^{\prime\prime}-1}$
$\displaystyle\qquad\times\int
d\Omega_{R}\;Y_{\lambda^{\prime},\mu^{\prime}}(\hat{r}_{1})\cdots
Y_{\lambda+\lambda^{\prime\prime},-\mu-\mu^{\prime\prime}}(\hat{r}_{4})\,\mathcal{P}^{*}_{\ell_{1}\ell_{2}(\ell_{12})\ell_{3}\ell_{4}}(\hat{r}_{1},\hat{r}_{2},\hat{r}_{3},\hat{r}_{4}).$
(101)

The last line may be evaluated using the definition of the $N=4$ isotropic
basis function, and the Kronecker deltas resulting from the orthogonality of
the spherical harmonics involved will eliminate the sums over
$\lambda,\lambda^{\prime},\lambda^{\prime\prime},\mu,\mu^{\prime}$ and
$\mu^{\prime\prime}$. Importantly, there will then be no infinite sums.

Notably, comparing equations (95) and (100) reveals an expansion of the
inverse of a sum of four vectors in terms of the $N=4$ isotropic functions;
this is another illustration that these latter really are generalizations of
the Legendre polynomials.

##  8 Use in Enabling a Faster NPCF Algorithm for $N\geq 4$

Recently, an algorithm to compute N-Point Correlation Functions (NPCFs),
$N\geq 4$ has been proposed [2], based on the spherical-harmonic 3PCF
algorithm of [22, 14, 23, 24]. This algorithm has already been used to make
the first measurement of the connected 4PCF [25]. This algorithm projects the
multi-point correlation functions onto the isotropic basis functions of one
less index; i.e. an N-Point Correlation Function (NPCF) would require the
$(N-1)$-argument basis functions. This is because one point can be placed at
the origin of coordinates. The algorithm proceeds by sitting on each galaxy in
a survey, one at a time (the “primary”), and then taking its neighbors out to
some maximum distance $R_{\rm max}$ away and placing them in spherical shells
(“bins”). On each shell, the angular dependence of the pattern (about the
primary) of galaxies on that shell is expanded into spherical harmonics, with
coefficients $a_{\ell m}$. These coefficients may then be combined, weighted
by a 3-$j$ symbol or an appropriate combination of 3-$j$ symbols, to form the
projection of the NPCF about that primary onto the isotropic basis. This
process is repeated on every primary and the results are then averaged.

This algorithm is quite fast because it never explicitly needs to consider
$(N-1)$-tuplets about a primary; obtaining the harmonic expansion coefficients
about each primary scales as $nV_{\rm max}$ where $n$ is the number density
and $V_{\rm max}$ is the spherical volume out to $R_{\rm max}$ about the
primary, and must be done per primary, so the total scaling is $N_{\rm
g}nV_{\rm max}$ if there are $N_{\rm g}$ galaxies in the survey. However,
about each primary, the appropriately-weighted combination of $(N-1)$
coefficients $a_{\ell m}$ must be formed, and this must be done for all bin
combinations. Put more abstractly, if we consider the $a_{\ell m}$ at each
fixed $\ell$ and $m$ a vector of length $N_{\rm bins}$, then we form the outer
product on bins of the $a_{\ell m}$. Furthermore, at each fixed combination of
bins, there are many different combinations of $\ell$ that are permitted by
isotropy. And the combinations of $a_{\ell m}$ must be summed over spins (this
comes from the fact that the clustering should be agnostic as to the $z$-axis,
yet the spins know about the $z$-axis—so they must be summed out).

Overall, this forming many $\ell$ combinations and then summing over spins
leaves to a complicated, multi-level loop. In practice, in the encore code,
one has a nested 9-level loop. Fundamentally there are selection rules on the
$\ell$ combinations, and so one cannot unroll and parallelize the loops any
further. Consequently on the CPU, this step takes the dominant amount of time
for 4PCF and higher.

Here, we suggest an algorithm that is less efficient, judged in terms of
FLOPS, but may be faster on the CPU, and takes advantage of the generating
functions presented above. If each isotropic basis function could be written
directly in Cartesian coordinates, as the generating functions indeed make
straightforward, then the orthogonality integral used to project the full NPCF
onto the basis could be done directly at each combination of radial bins, with
no need for any loops over the $\ell$ or the $m$. This approach would be less
efficient in terms of FLOPS, because the Cartesian elements from which basis
function is built are not orthogonal to those for other basis functions. Thus,
one would repeat computation. The advantage of the spherical-harmonic-based
approach is that one never repeats any computation—each $a_{\ell m}$ is
independent and computed once. However, computation is not the limiting step
in the 4PCF and higher algorithm, so there is presumably room to repeat some
of it to enable avoiding the slow nested loops.

However, there is still the problem of factorization: naively, the Cartesian
isotropic basis functions involve dot and cross products of pairs of unit
vectors, and so might require forming an $(N-1)$-tuplet about a primary to be
able to perform the projection integral; this would scale catastrophically.
Here, however, we may use the trick originally developed by [26], where one
multiplies out the components of the Cartesian dot and cross products, and
thus has a sum of integrals, each of which can be factorized. For instance,
$\mathcal{P}_{112}$ involves
$(\hat{r}_{1}\cdot\hat{r}_{3})(\hat{r}_{2}\cdot\hat{r}_{3})$ (see [1],
Appendix A, equation A.2). Multiplying this out, we find

$\displaystyle\mathcal{P}_{112}\supset(\hat{r}_{1}\cdot\hat{r}_{3})(\hat{r}_{2}\cdot\hat{r}_{3})=r_{1x}^{2}r_{2x}r_{3x}+r_{1x}r_{1y}r_{2y}r_{3x}+\cdots+r_{1z}^{2}r_{2z}r_{3z};$
(102)

$r_{ix}$ is the $x$-component of the vector $\boldsymbol{r}_{i}$, etc.

If we then wish to integrate the product of binned densities,
$\bar{\delta}(R_{i};\hat{r}_{i})$ (where the bin is designated by $R_{i}$ and
the bar denotes “binned”), against $\mathcal{P}^{*}_{112}$, we may form
integrals like

$\displaystyle\int d\Omega_{1}\;\bar{\delta}(R_{1};\hat{r}_{1})r_{1x}^{2}\int
d\Omega_{2}\;\bar{\delta}(R_{2};\hat{r}_{2})r_{2x}\int
d\Omega_{3}\;\bar{\delta}(R_{3};\hat{r}_{3})r_{3x}+\cdots.$ (103)

For a discrete density field such as galaxies, these integrals simply become
sums over each galaxy in the shell, weighted by the appropriate power of
$r_{ix}$, etc. Since they are factorized, they allow one to avoid examining an
$(N-1)$-tuplet around each primary. Indeed, obtaining each integral simply
scales as $nV_{\rm max}$ just as in the spherical-harmonic-based algorithm.

Now, in forthcoming work (Slepian, Warner, Hou & Cahn in prep.), we outline a
GPU-based implementation of encore called cadenza. One might wonder if the
above scheme would offer improvements for cadenza in the same way we suspect
it will for encore. However, the GPU code never loops over the $\ell$ and $m$.
GPUs are able to run a large number of concurrent threads, of the order of 1.4
million for typical hardware (NVIDIA A100), so each combination of the $\ell$
and $m$ is unrolled to a separate thread and already computed in parallel on
the GPU. Thus, there is no further advantage to be gained using the scheme
outlined above, whose primary goal was to avoid the loops used on the CPU.

##  9 Concluding Discussion

The isotropic basis functions [1] offer a natural basis for expanding
functions of $N$ unit vectors that are invariant under joint rotations of all
the arguments. Here, we have shown systematically how to generate isotropic
basis functions’ representations as sums of simple dot products (for the even-
parity sector) by expanding a plane wave in two different ways, and then
taking Taylor series and matching power by power. For the odd-parity sector
(possible for basis functions with three arguments or more), the presence of a
scalar triple product complicates this. However we have shown (§5) that given
the lowest-lying parity-odd basis function, one can solve for the other
parity-odd ones given the even ones (as obtained from our generating function
procedure) and the linearization identity for products of isotropic basis
functions (in [1]).

The introduction discusses a number of applications that have already been
made of the isotropic basis functions, specifically in measuring galaxy
correlations. The main utility of the present work is that it enables easily
obtaining simple, dot product (“Cartesian”) expressions for the isotropic
basis functions. These expressions are not easily obtained from the spherical-
harmonic-based definition of the isotropic basis functions presented in [1].
In turn, the Cartesian expressions can enable a speedier CPU-based algorithm
for measuring correlations in the basis, which avoids the many-layered,
highly-restricted loop required by the selection rules of the 3-$j$ symbols
involved in the current approach to computing galaxy clustering in the
isotropic basis (§8).

We have also given useful expansions of both multi-argument Dirac delta
functions (§6), and inverse powers of vector sums (§7), in the isotropic
basis, and shown the intimate link between these latter and spherical Bessel
overlap integrals. We expect these expressions will find use in inflationary
calculations of multi-point correlations, as well as cosmological perturbation
theory for the late-time evolution of galaxy clustering (e.g. the type of
integrals discussed in [20])

## Acknowledgments

We thank RN Cahn for many useful conversations, and the Slepian research group
members for helpful discussions over the past years.

## References

  * [1] Robert N. Cahn and Zachary Slepian.  Isotropic N-point basis functions and their properties.  Journal of Physics A Mathematical General, 56(32):325204, August 2023.
  * [2] Oliver H. E. Philcox, Zachary Slepian, Jiamin Hou, Craig Warner, Robert N. Cahn, and Daniel J. Eisenstein.  ENCORE: An {O}(N_g(2)) Estimator for Galaxy N-Point Correlation Functions.  MNRAS, October 2021.
  * [3] Oliver H. E. Philcox, Jiamin Hou, and Zachary Slepian.  A First Detection of the Connected 4-Point Correlation Function of Galaxies Using the BOSS CMASS Sample.  arXiv e-prints, page arXiv:2108.01670, August 2021.
  * [4] Robert N. Cahn, Zachary Slepian, and Jiamin Hou.  Test for cosmological parity violation using the 3d distribution of galaxies.  Phys. Rev. Lett., 130:201002, May 2023.
  * [5] Jiamin Hou, Robert N. Cahn, Oliver H. E. Philcox, and Zachary Slepian.  Analytic Gaussian covariance matrices for galaxy N -point correlation functions.  Phys. Rev. D, 106(4):043515, August 2022.
  * [6] Jiamin Hou, Zachary Slepian, and Robert N. Cahn.  Measurement of parity-odd modes in the large-scale 4-point correlation function of Sloan Digital Sky Survey Baryon Oscillation Spectroscopic Survey twelfth data release CMASS and LOWZ galaxies.  MNRAS, 522(4):5701–5739, May 2023.
  * [7] Oliver H. E. Philcox.  Probing parity violation with the four-point correlation function of boss galaxies.  Phys. Rev. D, 106:063501, Sep 2022.
  * [8] Giovanni Cabass, Mikhail M. Ivanov, and Oliver H. E. Philcox.  Colliders and ghosts: Constraining inflation with the parity-odd galaxy four-point function.  Phys. Rev. D, 107(2):023523, January 2023.
  * [9] William Ortolá Leonard, Zachary Slepian, and Jiamin Hou.  A Model for the Redshift-Space Galaxy 4-Point Correlation Function.  arXiv e-prints, page arXiv:2402.15510, February 2024.
  * [10] Sil, NC Deb, NC.  A simple approach to prove the addition theorem for regular solid harmonics.  Indian Journal of Physics, pages (3), 267 – 270, May 1994.
  * [11] S Chakrabarti and D P Dewangan.  Addition theorems for solid harmonics and the second born amplitudes.  Journal of Physics B: Atomic, Molecular and Optical Physics, 28(24):L769, dec 1995.
  * [12] Boisvert, R Olver, F, Lozier, D and Clark, C.  NIST Handbook of Mathematical Functions.  NIST and Cambridge University Press, 2010.
  * [13] Oliver H. E. Philcox and Zachary Slepian.  Efficient computation of N-point correlation functions in D dimensions.  Proceedings of the National Academy of Science, 119(33):e2111366119, August 2022.
  * [14] Zachary Slepian and Daniel J. Eisenstein.  A practical computational method for the anisotropic redshift-space three-point correlation function.  MNRAS, 478(2):1468–1483, August 2018.
  * [15] Zachary Slepian.  Algorithm to Produce a Density Field with Given Two, Three, and Four-Point Correlation Functions.  arXiv e-prints, page arXiv:2306.05383, May 2023.
  * [16] Wayne Hu.  Angular trispectrum of the cosmic microwave background.  PRD, 64(8):083005, October 2001.
  * [17] R. Mehrem, J. T. Londergan, and M. H. Macfarlane.  Analytic expressions for integrals of products of spherical Bessel functions.  Journal of Physics A Mathematical General, 24(7):1435–1453, April 1991.
  * [18] Rami Mehrem.  The Plane Wave Expansion, Infinite Integrals and Identities involving Spherical Bessel Functions.  arXiv e-prints, page arXiv:0909.0494, September 2009.
  * [19] Kiersten Meigs and Zachary Slepian.  On a general method for resolving integrals of multiple spherical bessel functions against power laws into distributions, 2021.
  * [20] Zachary Slepian.  On decoupling the integrals of cosmological perturbation theory.  MNRAS, 507(1):1337–1360, October 2021.
  * [21] Andrew D. Jackson and Leonard C. Maximon.  Integrals of products of bessel functions.  Siam Journal on Mathematical Analysis, 3:446–460, 1972.
  * [22] Zachary Slepian and Daniel J. Eisenstein.  Computing the three-point correlation function of galaxies in O(N2̂) time.  MNRAS, 454(4):4142–4158, December 2015.
  * [23] Zachary Slepian and Daniel J. Eisenstein.  Accelerating the two-point and three-point galaxy correlation functions using Fourier transforms.  MNRAS, 455(1):L31–L35, January 2016.
  * [24] Stephen K. N. Portillo, Zachary Slepian, Blakesley Burkhart, Sule Kahraman, and Douglas P. Finkbeiner.  Developing the 3-point Correlation Function for the Turbulent Interstellar Medium.  ApJ, 862(2):119, August 2018.
  * [25] Oliver H. E. Philcox, Jiamin Hou, and Zachary Slepian.  A First Detection of the Connected 4-Point Correlation Function of Galaxies Using the BOSS CMASS Sample.  arXiv e-prints, page arXiv:2108.01670, August 2021.
  * [26] Davide Bianchi, Héctor Gil-Marín, Rossana Ruggeri, and Will J. Percival.  Measuring line-of-sight-dependent Fourier-space clustering using FFTs.  MNRAS, 453(1):L11–L15, October 2015.