# On a General Method for Resolving Integrals of Multiple Spherical Bessel
Functions Against Power Laws into Distributions

Kiersten Meigs  Zachary Slepian  Department of Astronomy, University of
Florida, 211 Bryant Space Science Center, Gainesville, FL 32611, USA Physics
Division, Lawrence Berkeley National Laboratory, 1 Cyclotron Road, Berkeley,
CA 94709, USA

###### Abstract

We here present a method of performing integrals of products of spherical
Bessel functions (SBFs) weighted by a power-law. Our method, which begins with
double-SBF integrals, exploits a differential operator $\hat{D}$ defined via
Bessel’s differential equation. Application of this operator raises the power-
law in steps of two. We also here display a suitable base integral expression
to which this operator can be applied for both even and odd cases. We test our
method by showing that it reproduces previously-known solutions. Importantly,
it also goes beyond them, offering solutions in terms of singular
distributions, Heaviside functions, and Gauss’s hypergeometric,$\;{}_{2}{\rm
F}_{1}$ for all double-SBF integrals with positive semi-definite integer
power-law weight. We then show how our method for double-SBF integrals enables
evaluating arbitrary triple-SBF overlap integrals, going beyond the cases
currently in the literature. This in turn enables reduction of arbitrary
quadruple, quintuple, and sextuple-SBF integrals and beyond into tractable
forms.

##  1 Introduction

###  1.1 Motivation

Integrals of spherical Bessel functions (SBFs) are ubiquitous in cosmology.
Integrals against $j_{0}$ are used to transform the power spectrum into the
correlation function, or the bispectrum into the 3-Point Correlation Function
(e.g. [1], [2], [3]), and integrals against $j_{1}$ and $j_{2}$ respectively
enter models involving either a vector or a tensor dependence of the model
(e.g. [4]). They also appear in exact consideration of fully spherical
redshift-space distortions (e.g. [5], [6]). Integrals of two SBFs enter both
the covariance matrix of the correlation function (e.g. [7]) and the
computation of the CMB lensing angular power spectrum (often done using the
Limber approximation). Double SBF integrals, and even triple SBF integrals,
also enter computation of the galaxy 3PCF and 4PCF covariance matrices ([8],
[9], [10]) and indeed it can be shown that the Gaussian Random Field
approximation of the covariance for any N-Point Correlation Function is
expressible in terms of double and triple SBF integrals against the power
spectrum [10]. Famously, separating the double-SBF integrals against the
source function into a “geometric” and “numerical” part is at the heart of the
acceleration of the CMBFAST algorithm for CMB anisotropies [11]. Double and
triple SBF integrals also enter the computation of loop corrections to the
power spectrum [12, 13, 14] in perturbation theory and the redshift-space
power spectrum in Effective Field Theory of Large-Scale Structure [15].
Recently, there has been a movement to expand the power spectrum into complex
power laws and then perform the SBF-against-power-law integrals analytically
and sum them with the appropriate expansion-coefficient weights just once at
the end. This approach has been used both to accelerate computing perturbation
theory [16] and angular statistics models (both primary CMB [17], lensing [18]
and projected [19], [20]) and also in computing the relevant covariance
matrices [21]. SBF overlap integrals even appear in such diverse cosmological
contexts as a toy model for the effect of a survey boundary [22], analytic
treatment of the development of the Baryon Acoustic Oscillation feature in the
matter transfer function [23], and faster computation of the root-mean-square
mass fluctuations on a given smoothing scale [24].

With these many reasons for interest in such integrals, it is worthwhile to
expand our arsenal of techniques for performing them. In particular, a lacuna
in the literature has been resolution of singular versions of such integrals,
i.e. ones where the power-law weight is sufficient to make the integral
diverge in the UV. Of course, a separate matter of interest might be IR-
divergent such integrals, but we leave this unconsidered here. This lacuna has
developed because, while Dirac delta functions are of great use to physicists,
mathematicians tend to avoid them as defective functions. Yet, the types of
singular integrals we treat in the present work can arise in cosmological
Perturbation Theory (PT) (e.g. Hou & Slepian in prep.). These integrals are
also of intrinsic interest.

###  1.2 Prior work: convergent double and triple sBF integrals

Prior work on such integrals has been extensive and we summarize it here. We
organize our discussion by similar integrals; we first focus on the most
recent results and then proceed chronologically. Though little has been done
on indefinite integrals and they are not our concern here, we note for
completeness the work in this regard of [25] and [26].

Regarding definite integrals, the series of papers by Mehrem has been of
significant use and generalizes beyond early results e.g. in [27] or in
Watson’s Treatise [28]. Mehrem shows that a triple-SBF integral satisfying
certain conditions on the orders and the power can be written as a sum over
angular momentum symbols weighting Legendre polynomials in, essentially, the
cosine of the angle enclosed by the triangle formed by the three free
arguments of the integrand (e.g. $r_{1},r_{2}$ and $r_{3}$ in $\int
k^{2}dk\;j_{\ell_{1}}(kr_{1})j_{\ell_{2}}(kr_{2})j_{\ell_{3}}(kr_{3})$) [29],
[30]. Mehrem also included in double-SBF integrals [31] additional factors in
the integrand beyond power-laws, such as an exponential and polynomial [32]
and a Gaussian [33]. Importantly, Mehrem also displayed useful relationships
between the plane-wave expansion and such integrals [29]. Though his work
exploits that a Dirac delta function can be expressed as a Fourier Transform
of unity, and hence plane-wave expansions can be invoked, his end results
focus on only convergent integrals.

Much earlier work exploiting the plane-wave expansion is [34], which obtains
integrals of three spherical Bessel functions with $k^{2}$ weight and subject
to the condition that the sum of the orders is even and that the orders
satisfy triangular inequalities; if these two conditions hold, then the
integrals vanish unless the three free arguments also can form a triangle.
[35] obtains solutions for integrals of two Bessel functions of the first kind
$(J_{n})$ times a Hankel function of the first kind $(H_{n}^{(1)})$ when there
are constraints between the three orders and the power-law weight that enable
factorizing the Appell F4 function of the free arguments (appearing in
Bailey’s formula, [27]) into a product of two hypergeometrics; interestingly,
this work finds a triangular inequality between the three free arguments, as
often appears in these types of integrals. [36] extends the approach of the
previous work to modified Bessel functions of the first and second kind
(respectively $I_{n}$ and $K_{n}$) and to non-triangular combinations of free
arguments. [37] computes a three-SBF integral with $k^{2}$ weight and subject
to the restriction that the orders satisfy a triangle inequality. He also
discusses a factor of $1/2$ that must be introduced at the endpoints of the
integration domain that stem from use of a Delta function; we will encounter a
similar such factor in this work.111We note a typo in his result, equation 9;
his $\theta(k,k^{\prime},n)$ there should be $\theta(k,k^{\prime},\rho)$.

[38] uses an approach first outlined in [39] to recursively obtain triple-sBF
integrals with a $k^{2}$ weight subject to the conditions that the sum of the
orders is even and that the three free arguments form a triangle. The
recursion they find (their equation 47) is able to raise the order of one sBF
in the integrand at a time. [40] computes three-SBF integrals with $k^{2}$
weight (their equation A18 and following), using contour integration and
proving that the integral vanishes if the three free arguments cannot form a
triangle. They also obtain a four-SBF integral with the same weight (their
equation B11) and subject to the restrictions that its orders have an even
sum; they show there are also several inequalities on the free arguments.
Their result involves an infinite sum. [41] presents solutions for integrals
of three SBFs (see their equation 66 and onward) weighted by powers of $k^{2}$
and $k^{-1}$, representing the SBFs in terms of complex exponentials, giving a
result in terms of a finite series (but with many terms, which aspect the
authors note poses challenges to convergence when used in numerical work).
Their procedure is a development of one originally presented in [42], which
latter work they note did not handle completely the analysis of singularities
in the intermediate steps. There is no restriction on the sum of the orders in
[41].

###  1.3 Most Related Works, Other Works of Interest, & Plan of the Paper

Perhaps the two most closely related works to the present one are, first,
[43], which explores integrals of the form $\int
k^{2}dk\;j_{\ell_{1}}(kr_{1})j_{\ell_{2}}(kr_{2})$ where $\ell_{1}$ and
$\ell_{2}$ need not be equal. That work finds that for even-parity
$\ell_{1}+\ell_{2}$, the integrals can be expressed in terms of Dirac delta
functions, but for odd-parity $\ell_{1}+\ell_{2}$, the results never involve
Delta functions, but only Heavisides and kernels coupling $r_{1}$ and $r_{2}$
that are related to the Legendre functions. In the present work, we expand
this result to encompass all positive-semi-definite powers of $k$ in the
integrand.

Second, [44] uses an expansion of a spherical harmonic evaluated on the
gradient operator, as well as a translation formula, to show that double SBF
integrals can be obtained as derivatives of a Dirac delta function, subject to
the condition that the sum of the power-law and the two SBF orders be even
(that work’s equation (26), where $\sigma$ denotes a spherically-symmetric
Delta function). This is a restriction the present work will relax. Here we
also treat integrals of three or more SBFs weighted by an arbitrary (integer,
positive-semi-definite) power, extending Mehrem’s work in that regard. Indeed,
we show how the results we obtain allow one to evaluate an integral of an
arbitrary number of SBFs weighted by an arbitrary power.

Before turning to the main text, we mention briefly several other less-
directly related but interesting works. [45] finds an exact integral-to-sum
identity for double SBF overlap integrals, and we obtain a simpler proof of
this in [46] as well as extending it to products of an arbitrary number of
SBFs. [47] relates certain such integrals (e.g. one over a product of
$j_{0}$’s) to the behavior of random walkers, leading to interesting and
simple means of proving some otherwise unintuitive results. [48] adopts a
parametric-differentiation approach similar to ours (but using the Rayleigh
formula) to evaluate triple and quadruple SBF integrals, yet does not treat
divergent cases. [19] suggests using the spherical Bessel differential
equation to raise the power of the integrand, as we do here. [49] derives a
number of useful Fourier Transform pairs using single and double sBF integrals
as well as delta functions and expansions into spherical harmonics.

The paper is structured as follows. In §2, we summarize the properties of SBFs
we will find most germane. In §3 we outline our method and the base results
needed. §4 presents our detailed calculations and final formula. We then
display a few test cases in §5. §6 extends the work to triple-SBF integrals,
and then §7 shows how to treat an arbitrary number. We conclude in §8.

##  2 Brief Review of Relevant Properties

The spherical Bessel functions are solutions to Bessel’s differential
equation:

$\displaystyle{x^{2}}{\frac{d^{2}f}{dx^{2}}}+{2x}{\frac{df}{dx}}+{[{x^{2}}-{\ell}{(\ell+1)}]f}=0,$
(1)

where the solution $f(x)=j_{\ell}(x)$ is the spherical Bessel of the first
kind, and the solution $f(x)=y_{\ell}(x)$ is the spherical Bessel of the
second kind.

Spherical Bessel functions of the first kind are related to Bessel functions
by

$\displaystyle j_{\ell}(x)=\sqrt{\frac{\pi}{2x}}J_{\ell+1/2}(x).$ (2)

The first few spherical Bessel functions of the first kind are:

$\displaystyle j_{0}(x)$
$\displaystyle={\frac{\sin{x}}{x}}=\mathrm{sinc{(x)}},$ $\displaystyle
j_{1}(x)$ $\displaystyle={\frac{\sin{x}}{x^{2}}}-{\frac{\cos{x}}{x}},$
$\displaystyle j_{2}(x)$
$\displaystyle=\left(\frac{3}{x^{2}}-1\right)\frac{\sin{x}}{x}-\frac{3\cos{x}}{x^{2}}.$
(3)

Rayleigh’s formula generates the order-$\ell$ SBF from the order-zero one, as

$\displaystyle
j_{\ell}(x)={(-x)^{\ell}}{\left(\frac{1}{x}\frac{d}{dx}\right)}^{\ell}{\frac{\sin{x}}{x}}.$
(4)

The behavior for small values of $x$ is

$\displaystyle j_{\ell}(x)\approx{\frac{x^{\ell}}{(2{\ell}+1)!!}},$ (5)

where $!!$ denotes the double factorial, which is the product of all integers
from one up to $n$ (inclusive) that share the parity of $n$. The behavior for
large values of $x$ is

$\displaystyle
j_{\ell}(x)\approx{\frac{1}{x}}{\sin{\left(x-{\frac{{\ell}{\pi}}{2}}\right)}}.$
(6)

The closure relation is:

$\displaystyle\int_{0}^{\infty}x^{2}dx\;j_{\ell}(kx)j_{\ell}(k^{\prime}x)={\frac{\pi}{2kk^{\prime}}}{\delta_{\rm
D}^{[1]}}(k-k^{\prime}),$ (7)

where $\delta_{\rm D}^{[1]}$ is a 1D Dirac delta function.

The orthogonality relation is:

$\displaystyle\int_{-\infty}^{\infty}{j_{\ell}(x)}{j_{\ell^{\prime}}(x)}{dx}={\frac{\pi}{2{\ell}+1}}{\delta^{\rm
K}_{\ell\ell^{\prime}}},$ (8)

where $\delta^{\rm K}_{\ell\ell^{\prime}}$ is the Kronecker delta, which is
unity when $\ell=\ell^{\prime}$ and zero otherwise.

##  3 Method and Fundamental Results Used

###  3.1 Integral Form and Operator

We are concerned with integrals of the form:

$\displaystyle
I_{\ell\ell^{\prime}}^{[n]}(r,r^{\prime})=\int_{0}^{\infty}dkk^{n}j_{\ell}(kr)j_{\ell^{\prime}}(kr^{\prime}).$
(9)

In what follows, we suppress the upper and lower bounds as we will always be
considering half-infinite integrals.

We may use Bessel’s differential equation to define an operator that will act
on the SBF and by parametric differentiation introduce an additional power of
$k^{2}$:

$\displaystyle\hat{D}_{\ell,\chi}j_{\ell}(k\chi)\equiv-\left[\frac{\partial^{2}}{\partial\chi^{2}}+\frac{2}{\chi}\frac{\partial}{\partial\chi}-\frac{\ell(\ell+1)}{\chi^{2}}\right]\;j_{\ell}(k\chi)=k^{2}j_{\ell}(k\chi).$
(10)

This method was suggested and used to write a formal result for even powers of
$k$ against two SBFs of the same order in [19].

To use this operator, we need a base integral for equation (9) for some value
of $n$. There are certain $n$ for which the integral is convergent; we begin
with them. We will need a formula that has both an odd and an even allowed
starting $n$ since applying our operator cannot alter the parity of the power
of $k$.

We observe that

$\displaystyle\hat{D}_{\ell,r}I_{\ell\ell^{\prime}}^{[n]}(r,r^{\prime})=I_{\ell\ell^{\prime}}^{[n+2]}(r,r^{\prime})=\int
dk\,k^{n+2}j_{\ell}(kr)j_{\ell^{\prime}}(kr^{\prime}).$ (11)

The $\hat{D}_{\ell,r}$ operator can be repeatedly applied to form a “ladder”
of integrals with $n$ ascending in steps of two. This means that an integral
of the form (9) can be evaluated for arbitrary $n$ by applying the $\hat{D}$
operator a suitable number of times, as long as we have a result for some
starting value of $n$; in particular we will need both even and odd parity-$n$
starting cases as the operator cannot change the parity of the power.

###  3.2 Base Case

We find in [50] integrals that offer the required base cases for even and odd
$n$. We set variables in their work to those in ours as follows: $\alpha\to
r$, $\beta\to r^{\prime}$, $\nu\to\ell+1/2$, $\mu\to\ell^{\prime}+1/2$, $t\to
k$, $\lambda\to 1-n$. We use colors to separate the cases in [50]; blue is for
$r^{\prime}>r$ and red is for $r>r^{\prime}$. We now present the base results
we will require, showing both the full form and then a shorthand one in terms
of coefficients which in each case we will use in what follows. We use
equation (2) to convert the spherical Bessel functions into Bessel functions;
the [50] formulae we use are in terms of these latter.

  * •

Case where $r^{\prime}>r$ ([50] equation 6.574.1):

$\displaystyle{\color[rgb]{0,0,1}I_{r^{\prime}>}^{[n]}=\frac{\pi}{2^{2-n}}\frac{r^{\ell}}{r^{\prime\ell+1+n}}\frac{\Gamma\left(\frac{1}{2}(\ell+\ell^{\prime}+1+n)\right)}{\Gamma\left(\frac{1}{2}(\ell^{\prime}-\ell+2-n)\right)\Gamma\left(\ell+\frac{3}{2}\right)}}$
$\displaystyle\qquad{\color[rgb]{0,0,1}\times\;_{2}F_{1}\left(\frac{1}{2}(\ell+\ell^{\prime}+n+1),\frac{1}{2}(\ell-\ell^{\prime}+n);\ell+\frac{3}{2};\left({\frac{r}{r^{\prime}}}\right)^{2}\right).}$
(12)

We define both a constant for the term in the first line that involves Gamma
functions, as

$\displaystyle{\color[rgb]{0,0,1}c_{\Gamma
r^{\prime}>}=\frac{\Gamma\left(\frac{1}{2}(\ell+\ell^{\prime}+1+n)\right)}{\Gamma\left(\frac{1}{2}(\ell^{\prime}-\ell+2-n)\right)\Gamma\left(\ell+\frac{3}{2}\right)},}$
(13)

and shorthand notation for the hypergeometric as:

$\displaystyle{\color[rgb]{0,0,1}F_{{}^{\prime}xyz}=\;_{2}F_{1}\left(\frac{1}{2}(\ell+\ell^{\prime}+n+x),\frac{1}{2}(\ell-\ell^{\prime}+n+y);\ell+\frac{z}{2};\left({\frac{r}{r^{\prime}}}\right)^{2}\right)}.$
(14)

Using our definitions to shorten the coefficients and the hypergeometric, this
becomes

$\displaystyle{\color[rgb]{0,0,1}\boxed{I_{r^{\prime}>}^{[n]}=\frac{\pi}{2^{2-n}}c_{\Gamma
r^{\prime}>}\frac{r^{\ell}}{r^{\prime\ell+1+n}}F_{{}^{\prime}103}}}$ (15)

where $\ell+\ell^{\prime}+n+1>0$, $n<2$, and $0<r<r^{\prime}$.

  * •

Case where $r^{\prime}<r$ ([50] equation 6.574.3)

$\displaystyle{\color[rgb]{1,0,0}I_{r>}^{[n]}=\frac{\pi}{2^{2-n}}\frac{r^{\prime\ell}}{r^{\ell^{\prime}+1+n}}\frac{\Gamma\left(\frac{1}{2}(\ell+\ell^{\prime}+1+n)\right)}{\Gamma\left(\frac{1}{2}(\ell-\ell^{\prime}+2-n)\right)\Gamma\left(\ell+\frac{3}{2}\right)}}$
$\displaystyle\qquad{\color[rgb]{1,0,0}\times\;_{2}F_{1}\left(\frac{1}{2}(\ell+\ell^{\prime}+n+1),\frac{1}{2}(\ell^{\prime}-\ell+n);\ell^{\prime}+\frac{3}{2};\left({\frac{r^{\prime}}{r}}\right)^{2}\right).}$

We again define a constant for the term in the first line involving the Gamma
functions:

$\displaystyle{\color[rgb]{1,0,0}c_{\Gamma
r>}=\frac{\Gamma\left(\frac{1}{2}(\ell+\ell^{\prime}+1+n)\right)}{\Gamma\left(\frac{1}{2}(\ell-\ell^{\prime}+2-n)\right)\Gamma\left(\ell+\frac{3}{2}\right)},}$
(16)

and a shorthand for the hypergeometric:

$\displaystyle{\color[rgb]{1,0,0}F_{xyz}=\;_{2}F_{1}\left(\frac{1}{2}(\ell+\ell^{\prime}+n+x),\frac{1}{2}(\ell^{\prime}-\ell+n+y);\ell^{\prime}+\frac{z}{2};\left({\frac{r^{\prime}}{r}}\right)^{2}\right)}.$
(17)

Using our definitions to shorten the coefficients and the hypergeometric, this
becomes

$\displaystyle{\color[rgb]{1,0,0}\boxed{I_{r>}^{[n]}=\frac{\pi}{2^{2-n}}c_{\Gamma
r>}\frac{r^{\prime\ell}}{r^{\ell^{\prime}+1+n}}F_{103}}}$ (18)

where $\ell+\ell^{\prime}+n+1>0$, $n<2$, and $0<r^{\prime}<r$.

Instead of applying the operator to each of these cases separately, we combine
them into one expression for the integral to which the operator may then be
applied. We do so using Heaviside functions, as

$\displaystyle
I_{\ell\ell^{\prime}}^{[n]}(r,r^{\prime})={\color[rgb]{0,0,1}H(r^{\prime}-r)I_{r^{\prime}>}}+{\color[rgb]{1,0,0}H(r-r^{\prime})I_{r>}}.$
(19)

We intend the Heaviside functions as left-continuous $(H(0)=0)$ here. In our
representation (19), we omit the value of $I_{\ell\ell^{\prime}}^{[n]}$ at
$r=r^{\prime}$, which is given by [50] equation 6.574.2, and is not a limit of
the other two cases (as when $r\to r^{\prime}$ the argument of the
hypergeometric for them becomes unity, which can cause the hypergeometric to
diverge). We omit this value because it is only relevant on a set of measure
zero. Hence it is not needed if the results obtained in this work are
integrated over, as we expect would be the case in applications. Indeed, the
use of Dirac delta functions is justified only if they are ultimately
integrated over. In contrast to the $r=r^{\prime}$ value, we note that Delta
functions, defined for instance as the limit of a Gaussian as its width goes
to zero, are always supported on a set of nonzero measure even as this limit
is taken; hence they must be retained in the calculation.

##  4 Double-SBF Overlap Integral General Formula

The $\hat{D}$ operator outlined in equation (10) can now be applied to the RHS
of equation (19). We preserve our organization by colors by maintaining a
term’s color when a derivative is applied.

Beginning with the second term of the operator in equation (10) and for the
moment focusing only on the derivative and leaving the factor of $2/r$
unapplied (as we will use this raw first derivative to obtain the second
derivative shortly), we differentiate equation (19). We use the product rule,
the chain rule, and the relation $dH(x)/dx=\delta_{\rm D}^{[1]}(x)$, where
$\delta_{\rm D}^{[1]}(x)$ is the 1D Dirac delta function:

$\displaystyle\frac{\partial}{\partial
r}[I_{\ell\ell^{\prime}}^{[n]}(r,r^{\prime})]={\color[rgb]{0,0,1}-\frac{1}{2}\delta_{\rm
D}^{[1]}(r^{\prime}-r)I_{r^{\prime}>}+H(r^{\prime}-r)\frac{\partial}{\partial
r}I_{r^{\prime}>}}$ $\displaystyle+{\color[rgb]{1,0,0}\frac{1}{2}\delta_{\rm
D}^{[1]}(r-r^{\prime})I_{r>}+H(r-r^{\prime})\frac{\partial}{\partial
r}I_{r>}}.$ (20)

We retain the color of signs when they originate from applying the chain rule.

We now pause to outline an important subtlety: a careful reader will have
noted the pre-factors of $1/2$ introduced by hand on the Dirac Delta
functions. This is because in fact we wish to take only left-handed or right-
handed derivatives of the Heaviside functions; the blue region derivative may
only be taken as $r$ approaches $r^{\prime}$ from the left, whereas the red
region derivative may only be taken as $r$ approaches $r^{\prime}$ from the
right. The symmetry of the Dirac delta function implies that each of these
“handed” derivatives must be half of the total derivative. This argument is
along the same lines as that in [37] where he argues that a Dirac delta
function evaluated through the endpoint of an interval should only be counted
with $1/2$ weight (comment below his equation 8).

We now differentiate equation (4) to yield the first term of the operator in
equation (19):

$\displaystyle\frac{\partial^{2}}{\partial
r^{2}}I_{\ell\ell^{\prime}}^{[n]}(r,r^{\prime})$
$\displaystyle={\color[rgb]{0,0,1}-\frac{\partial}{\partial
r}\left[\frac{1}{2}\delta_{\rm
D}^{[1]}(r^{\prime}-r)\right]I_{r^{\prime}>}-\frac{1}{2}\delta_{\rm
D}^{[1]}(r^{\prime}-r)\frac{\partial}{\partial r}I_{r^{\prime}>}}$
$\displaystyle{\color[rgb]{0,0,1}-\frac{1}{2}\delta_{\rm
D}^{[1]}(r^{\prime}-r)\frac{\partial}{\partial
r}I_{r^{\prime}>}+H(r^{\prime}-r)\frac{\partial^{2}}{\partial
r^{2}}I_{r^{\prime}>}}$
$\displaystyle+{\color[rgb]{1,0,0}\frac{\partial}{\partial
r}\left[\frac{1}{2}\delta_{\rm
D}^{[1]}(r-r^{\prime})\right]I_{r>}+\frac{1}{2}\delta_{\rm
D}^{[1]}(r-r^{\prime})\frac{\partial}{\partial r}I_{r>}}$
$\displaystyle{\color[rgb]{1,0,0}+\frac{1}{2}\delta_{\rm
D}^{[1]}(r-r^{\prime})\frac{\partial}{\partial
r}I_{r>}+H(r-r^{\prime})\frac{\partial^{2}}{\partial r^{2}}I_{r>}}.$ (21)

Now substituting equations (4) and (4) into equation (19) for the operator, we
find:

$\displaystyle
I_{\ell\ell^{\prime}}^{[n+2]}(r,r^{\prime})=-\left[\frac{\partial^{2}}{\partial
r^{2}}+\frac{2}{r}\frac{\partial}{\partial
r}-\frac{\ell(\ell+1)}{r^{2}}\right]I_{\ell\ell^{\prime}}^{[n]}(r,r^{\prime})$
$\displaystyle\qquad=-\bigg{\\{}\bigg{[}{\color[rgb]{0,0,1}-\frac{\partial}{\partial
r}\left[\frac{1}{2}\delta_{\rm
D}^{[1]}(r^{\prime}-r)\right]I_{r^{\prime}>}-\delta_{\rm
D}^{[1]}(r^{\prime}-r)\frac{\partial}{\partial
r}I_{r^{\prime}>}+H(r^{\prime}-r)\frac{\partial^{2}}{\partial
r^{2}}I_{r^{\prime}>}}$
$\displaystyle\qquad+{\color[rgb]{1,0,0}\frac{\partial}{\partial
r}\left[\frac{1}{2}\delta_{\rm D}^{[1]}(r-r^{\prime})\right]I_{r>}+\delta_{\rm
D}^{[1]}(r-r^{\prime})\frac{\partial}{\partial
r}I_{r>}+H(r-r^{\prime})\frac{\partial^{2}}{\partial r^{2}}I_{r>}}\bigg{]}$
$\displaystyle\qquad+\frac{2}{r}\bigg{[}-{\color[rgb]{0,0,1}\frac{1}{2}\delta_{\rm
D}^{[1]}(r^{\prime}-r)I_{r^{\prime}>}+H(r^{\prime}-r)\frac{\partial}{\partial
r}I_{r^{\prime}>}}+{\color[rgb]{1,0,0}\frac{1}{2}\delta_{\rm
D}^{[1]}(r-r^{\prime})I_{r>}+H(r-r^{\prime})\frac{\partial}{\partial
r}I_{r>}}\bigg{]}$
$\displaystyle\qquad-\bigg{[}\frac{\ell(\ell+1)}{r^{2}}\Big{(}{\color[rgb]{0,0,1}H(r^{\prime}-r)I_{r^{\prime}>}}+{\color[rgb]{1,0,0}H(r-r^{\prime})I_{r>}}\Big{)}\bigg{]}\bigg{\\}}.$
(22)

We now group terms by the region in $(r,r^{\prime})$ space on which they will
apply, while retaining their original color. We introduce a green bracket to
group terms that will apply when $r=r^{\prime}$. We also factor out the Delta
functions and their derivatives by using their even parity, i.e. that
$\delta_{\rm D}^{[1]}(r-r^{\prime})=\delta_{\rm D}^{[1]}(r^{\prime}-r)$, and
also the relation ${\partial}\,\delta_{\rm D}^{[1]}(r-r^{\prime})/{\partial
r}={\partial}\,\delta_{\rm D}^{[1]}(r^{\prime}-r)/{\partial r}$. We obtain:

$\displaystyle
I_{\ell\ell^{\prime}}^{[n+2]}(r,r^{\prime})=-\bigg{\\{}{\color[rgb]{0,0,1}H(r^{\prime}-r)\bigg{(}\frac{\partial^{2}}{\partial
r^{2}}I_{r^{\prime}>}+\frac{2}{r}\frac{\partial}{\partial
r}I_{r^{\prime}>}-\frac{\ell(\ell+1)}{r^{2}}I_{r^{\prime}>}\bigg{)}}$
$\displaystyle+{\color[rgb]{1,0,0}H(r-r^{\prime})\bigg{(}\frac{\partial^{2}}{\partial
r^{2}}I_{r>}+\frac{2}{r}\frac{\partial}{\partial
r}I_{r>}-\frac{\ell(\ell+1)}{r^{2}}I_{r>}\bigg{)}}$
$\displaystyle+{\color[rgb]{0,0.88,0}\bigg{[}}\delta_{\rm
D}^{[1]}(r-r^{\prime})\bigg{(}{\color[rgb]{0,0,1}-\frac{\partial}{\partial
r}I_{r^{\prime}>}-\frac{I_{r^{\prime}>}}{r}}+{\color[rgb]{1,0,0}\frac{\partial}{\partial
r}I_{r>}+\frac{I_{r>}}{r}}\bigg{)}+\frac{\partial}{\partial
r}\left[\frac{1}{2}\delta_{\rm
D}^{[1]}(r-r^{\prime})\right]\bigg{(}{\color[rgb]{0,0,1}-I_{r^{\prime}>}}+{\color[rgb]{1,0,0}I_{r>}}\bigg{)}{\color[rgb]{0,0.88,0}\bigg{]}}\bigg{\\}}.$
(23)

We now must substitute our expressions for
${\color[rgb]{0,0,1}I_{r^{\prime}>}}$ and ${\color[rgb]{1,0,0}I_{r>}}$
(respectively, equations (15) and (18)) and their first and second derivatives
into equation (4). We display these latter here:

$\displaystyle{\color[rgb]{0,0,1}\frac{\partial}{\partial
r}I_{r^{\prime}>}=\frac{\pi}{2^{2-n}}c_{\Gamma
r^{\prime}>}\bigg{(}\frac{2r^{\ell+1}}{r^{\prime\ell+3+n}}d_{{}^{\prime}103}F_{{}^{\prime}325}+\frac{\ell
r^{\ell-1}}{r^{\prime\ell+1+n}}F_{{}^{\prime}103}\bigg{)}},$
$\displaystyle{\color[rgb]{0,0,1}\frac{\partial^{2}}{\partial
r^{2}}I_{r^{\prime}>}=\frac{\pi}{2^{2-n}}c_{\Gamma
r^{\prime}>}\bigg{(}\frac{4r^{\ell+2}}{r^{\prime\ell+5+n}}d_{{}^{\prime}103}d_{{}^{\prime}325}F_{547}+\frac{(4\ell+2)r^{\ell}}{r^{\prime\ell+3+n}}d_{{}^{\prime}103}F_{{}^{\prime}325}+\frac{(\ell^{2}-\ell)r^{\ell-2}}{r^{\prime\ell+1+n}}F_{{}^{\prime}103}\bigg{)}},$
$\displaystyle{\color[rgb]{1,0,0}\frac{\partial}{\partial
r}I_{r>}=-\frac{\pi}{2^{2-n}}c_{\Gamma
r>}\bigg{(}\frac{(\ell^{\prime}+1+n)r^{\prime\ell}}{r^{\ell^{\prime}+2+n}}F_{103}+\frac{2r^{\prime\ell+2}}{r^{\ell^{\prime}+4+n}}d_{103}F_{325}\bigg{)}},$
$\displaystyle{\color[rgb]{1,0,0}\frac{\partial^{2}}{\partial
r^{2}}I_{r>}=\frac{\pi}{2^{2-n}}c_{\Gamma
r>}\bigg{(}\frac{(\ell^{\prime}+2+n)(\ell^{\prime}+1+n)r^{\prime\ell}}{r^{\ell^{\prime}+3+n}}F_{103}+\frac{(4\ell^{\prime}+10+4n)r^{\prime\ell+2}}{r^{\ell^{\prime}+5+n}}d_{103}F_{325}}$
$\displaystyle\qquad\qquad\;{\color[rgb]{1,0,0}+\frac{4r^{\prime\ell+4}}{r^{\ell^{\prime}+7+n}}d_{103}d_{325}F_{{}^{\prime}547}\bigg{)}}.$
(24)

We have used the following relation for the derivative of a hypergeometric
function:222https://functions.wolfram.com/HypergeometricFunctions/Hypergeometric2F1/20/01/05/,
first formula.

$\displaystyle\frac{\partial}{\partial
u}\;_{2}F_{1}(a,b;c;u)=\frac{ab}{c}\;_{2}F_{1}(a+1,b+1;c+1;u).$ (25)

Given the factors of $1/2$ on each parameter in our shorthand $F_{xyz}$
(equations 14 and 17), raising each parameter by unity raises each lower index
by two, so the hypergeometric factor in the first derivative $F_{103}$ becomes
$F_{325}$.

We have also defined a shorthand for the constant $(ab)/c$ that stems from
differentiating the hypergeometric function:

$\displaystyle{\color[rgb]{0,0,1}d_{{}^{\prime}xyz}=\frac{(\ell+\ell^{\prime}+n+x)(\ell-\ell^{\prime}+n+y)}{4(\ell+{z}/{2})}},$
(26)

and

$\displaystyle{\color[rgb]{1,0,0}d_{xyz}=\frac{(\ell+\ell^{\prime}+n+x)(\ell^{\prime}-\ell+n+y)}{4(\ell^{\prime}+{z}/{2})}};$
(27)

the $4$ in the denominators comes from the fact that the first two parameters
of the hypergeometric have overall factors of $1/2$ that may be multiplied
together.

Finally, substituting these expressions into equation (4), simplifying and
factoring out constants, we have:

$\displaystyle
I_{\ell\ell^{\prime}}^{[n+2]}(r,r^{\prime})=-\frac{\pi}{2^{2-n}}\Bigg{\\{}{\color[rgb]{0,0,1}H(r^{\prime}-r)c_{\Gamma
r^{\prime}>}\bigg{(}\frac{4r^{\ell+2}}{r^{\prime\ell+5+n}}d_{{}^{\prime}103}d_{{}^{\prime}325}F_{{}^{\prime}547}+\frac{(4\ell+2)r^{\ell}}{r^{\prime\ell+3+n}}d_{{}^{\prime}103}F_{{}^{\prime}325}}$
$\displaystyle\;\;\;{\color[rgb]{0,0,1}+\frac{(\ell^{2}-\ell)r^{\ell-2}}{r^{\prime\ell+1+n}}F_{{}^{\prime}103}+\left(\frac{4r^{\ell}}{r^{\prime\ell+3+n}}d_{{}^{\prime}103}F_{{}^{\prime}325}+\frac{2\ell
r^{\ell-2}}{r^{\prime\ell+1+n}}F_{{}^{\prime}103}\bigg{)}-\frac{\ell(\ell+1)r^{\ell-2}}{r^{\prime\ell+1+n}}F_{{}^{\prime}103}\right)}$
$\displaystyle\;\;\;+{\color[rgb]{1,0,0}H(r-r^{\prime})c_{\Gamma
r>}\bigg{(}\frac{(\ell^{\prime}+2+n)(\ell^{\prime}+1+n)r^{\prime\ell}}{r^{\ell^{\prime}+3+n}}F_{103}+\frac{(4\ell^{\prime}+10+4n)r^{\prime\ell+2}}{r^{\ell^{\prime}+5+n}}d_{103}F_{325}}$
$\displaystyle\;\;\;{\color[rgb]{1,0,0}+\frac{4r^{\prime\ell+4}}{r^{\ell^{\prime}+7+n}}d_{103}d_{325}F_{547}+\left(-\frac{2(\ell^{\prime}+1+n)r^{\prime\ell}}{r^{\ell^{\prime}+3+n}}F_{103}-\frac{4r^{\prime\ell+2}}{r^{\ell^{\prime}+5+n}}d_{103}F_{325}\right)}{\color[rgb]{1,0,0}-\frac{\ell(\ell+1)r^{\prime\ell}}{r^{\ell+3+n}}F_{103}\bigg{)}}$
$\displaystyle\;\;\;+{\color[rgb]{0,0.88,0}\Bigg{[}}\delta_{\rm
D}^{[1]}(r-r^{\prime})\bigg{(}{\color[rgb]{0,0,1}c_{\Gamma
r^{\prime}>}\left(-\frac{\ell
r^{\ell-1}}{r^{\prime\ell+1+n}}F_{103}-\frac{2r^{\ell+1}}{r^{\prime\ell+3+n}}d_{{}^{\prime}103}F_{{}^{\prime}325}-\frac{r^{\ell-1}}{r^{\prime\ell+1+n}}F_{{}^{\prime}103}\right)}$
$\displaystyle\;+{\color[rgb]{1,0,0}c_{\Gamma
r>}\left(-\frac{(\ell^{\prime}+1+n)r^{\prime\ell}}{r^{\ell^{\prime}+2+n}}F_{103}-\frac{2r^{\prime\ell+2}}{r^{\ell^{\prime}+4+n}}d_{103}F_{325}+\frac{r^{\prime\ell}}{r^{\ell^{\prime}+2+n}}F_{103}\right)}\bigg{)}$
$\displaystyle\;\;\;+\frac{\partial}{\partial r}\left[\frac{1}{2}\delta_{\rm
D}^{[1]}(r-r^{\prime})\right]\left({\color[rgb]{0,0,1}-\frac{c_{\Gamma
r>}r^{\ell}}{r^{\prime\ell+1+n}}F_{{}^{\prime}103}}+{\color[rgb]{1,0,0}\frac{c_{\Gamma
r^{\prime}>}r^{\prime\ell}}{r^{\ell^{\prime}+1+n}}F_{103}}\right){\color[rgb]{0,0.88,0}\bigg{]}}\Bigg{\\}}.$
(28)

Combining terms we find:

$\displaystyle
I_{\ell\ell^{\prime}}^{[n+2]}(r,r^{\prime})=-\frac{\pi}{2^{2-n}}\Bigg{\\{}{\color[rgb]{0,0,1}H(r^{\prime}-r)c_{\Gamma
r^{\prime}>}\bigg{(}\frac{4r^{\ell+2}}{r^{\prime\ell+5+n}}d_{{}^{\prime}103}d_{{}^{\prime}325}F_{{}^{\prime}547}+\frac{(4\ell+6)r^{\ell}}{r^{\prime\ell+3+n}}d_{{}^{\prime}103}F_{{}^{\prime}325}\bigg{)}}$
$\displaystyle\;+{\color[rgb]{1,0,0}H(r-r^{\prime})c_{\Gamma
r>}\bigg{(}\frac{\left(2n\ell^{\prime}+\ell^{\prime}+n^{2}+\ell^{\prime
2}-\ell^{2}-\ell+n\right)r^{\prime\ell}}{r^{\ell^{\prime}+3+n}}F_{103}+\frac{(4\ell^{\prime}+6+4n)r^{\prime\ell+2}}{r^{\ell^{\prime}+5+n}}d_{103}F_{325}}$
$\displaystyle\;{\color[rgb]{1,0,0}+\frac{4r^{\prime\ell+4}}{r^{\ell^{\prime}+7+n}}d_{103}d_{325}F_{547}\bigg{)}}$
$\displaystyle\;+{\color[rgb]{0,0.88,0}\Bigg{[}}\delta_{\rm
D}^{[1]}(r-r^{\prime})\bigg{(}{\color[rgb]{0,0,1}c_{\Gamma
r^{\prime}>}\left(-\frac{2r^{\ell+1}}{r^{\prime\ell+3+n}}d_{{}^{\prime}103}F_{{}^{\prime}325}-\frac{(\ell+1)r^{\ell-1}}{r^{\prime\ell+1+n}}F_{{}^{\prime}103}\right)}$
$\displaystyle\;+{\color[rgb]{1,0,0}c_{\Gamma
r>}\left(-\frac{2r^{\prime\ell+2}}{r^{\ell^{\prime}+4+n}}d_{103}F_{325}-\frac{(\ell^{\prime}+n)r^{\prime\ell}}{r^{\ell^{\prime}+2+n}}F_{103}\right)}\bigg{)}$
$\displaystyle\;+\frac{\partial}{\partial r}\left[\frac{1}{2}\delta_{\rm
D}^{[1]}(r-r^{\prime})\right]{}\left({\color[rgb]{0,0,1}-\frac{c_{\Gamma
r^{\prime}>}r^{\ell}}{r^{\prime\ell+1+n}}F_{{}^{\prime}103}}+{\color[rgb]{1,0,0}\frac{c_{\Gamma
r>}r^{\prime\ell}}{r^{\ell^{\prime}+1+n}}F_{103}}\right){\color[rgb]{0,0.88,0}\Bigg{]}}\Bigg{\\}}.$
(29)

This is our final expression for the integral form outlined in equation (9),
and is a key result of this work. We have used Mathematica to confirm the
reduction starting directly from equation (19).333Relative to the work in
Mathematica, we note three points. First, when both $r$ and $r^{\prime}$ are
written purely symbolically (i.e. without any values inserted for either,
especially $r^{\prime}$), Mathematica writes the Heaviside function’s second
derivative simply as ${\rm DiracDelta^{\prime}[r-r^{\prime}]}$ or ${\rm
DiracDelta^{\prime}[r^{\prime}-r]}$, where prime is ambiguous and this result
does not capture the negative sign produced on differentiating the Heaviside
$H(r^{\prime}-r)$ for the first time (by chain rule). In contrast, Mathematica
does perform the differentiation correctly if a numerical value is specified
for $r^{\prime}$. Mathematica also performs a first differentiation of the
Heaviside correctly even with both variables’ remaining symbolic; the issue
only appears when the second derivative is requested ab initio. Thus, the sign
of the first blue term in the last line of equation (4) requires some thought
to compare with the Mathematica result; the negative sign we have above is
intended. Second, the factors of $1/2$ on the Delta functions and their
derivatives (as discussed around equation (4)) are not present in the
mathematica notebook as it does not handle correctly the question of left and
right-handed derivatives. Third, mathematica often uses the regularized
hypergeometric$\;{}_{2}\tilde{{\rm F}}_{1}$; a helpful relation is
that$\;{}_{2}\tilde{{\rm F}}_{1}(a,b,c;u)=[1/\Gamma(c)]\;_{2}{\rm
F}_{1}(a,b,c;u)$.

##  5 Tests

Let us begin by observing the behavior of some of the constants that we have
defined. The constants $d_{xyz}$ and $d_{{}^{\prime}xyz}$ that arose from
differentiating the hypergeometric function defined in equations (26) and (27)
will both become zero when $n+y=0$ and $\ell=\ell^{\prime}$ for any $\ell$ and
$\ell^{\prime}$. Additionally, $d_{xyz}=d_{{}^{\prime}xyz}$ when
$\ell=\ell^{\prime}$. Finally, we note that when $\ell=\ell^{\prime}$,
$c_{\Gamma r>}=c_{\Gamma r^{\prime}>}$.

###  5.1 $\ell=0=\ell^{\prime},\;n=0$

This test will ultimately yield an integral $\int
k^{2}dk\;j_{0}(kr)j_{0}(kr^{\prime})$, which we know gives a spherically-
symmetric Dirac delta function by the closure relation (7).

We will need to substitute our test parameters into equation (4). We first
observe from equations (27) and (26) that $d_{103}=0=d_{{}^{\prime}103}$ in
this case, meaning all terms in equation (4) containing either of these
coefficients will immediately drop out. Any terms proportional solely to sums
of powers of $\ell$, $\ell^{\prime}$, and $n$, will also drop out. We find

$\displaystyle
I_{00}^{[0+2]}(r,r^{\prime})=-\frac{\pi}{4}\Bigg{\\{}{\color[rgb]{0,0.88,0}\Bigg{[}}\delta_{\rm
D}^{[1]}(r-r^{\prime})\bigg{(}{\color[rgb]{0,0,1}c_{\Gamma
r^{\prime}>}\left(-\frac{1}{rr^{\prime}}F_{{}^{\prime}103}\right)}$
$\displaystyle\;+\frac{\partial}{\partial r}\left[\frac{1}{2}\delta_{\rm
D}^{[1]}(r-r^{\prime})\right]\left({\color[rgb]{0,0,1}-\frac{c_{\Gamma
r^{\prime}>}}{r^{\prime}}F_{{}^{\prime}103}}+{\color[rgb]{1,0,0}\frac{c_{\Gamma
r>}}{r}F_{103}}\right){\color[rgb]{0,0.88,0}\Bigg{]}}\Bigg{\\}}.$ (30)

We see that for this case, from equations (16) and (13) we have $c_{\Gamma
r>}=2=c_{\Gamma r^{\prime}>}$. Using this, we find

$\displaystyle
I_{00}^{[0+2]}(r,r^{\prime})=-\frac{\pi}{4}\Bigg{\\{}{\color[rgb]{0,0.88,0}\Bigg{[}}\delta_{\rm
D}^{[1]}(r-r^{\prime})\bigg{(}{\color[rgb]{0,0,1}\left(-\frac{2}{rr^{\prime}}F_{{}^{\prime}103}\right)}\bigg{)}$
$\displaystyle\;+\frac{\partial}{\partial r}\left[\frac{1}{2}\delta_{\rm
D}^{[1]}(r-r^{\prime})\right]\left({\color[rgb]{0,0,1}-\frac{2}{r^{\prime}}F_{{}^{\prime}103}}+{\color[rgb]{1,0,0}\frac{2}{r}F_{103}}\right){\color[rgb]{0,0.88,0}\Bigg{]}}\Bigg{\\}}.$
(31)

Using that $\partial\;\delta_{\rm D}^{[1]}(r-r^{\prime})/\partial
r=-(r-r^{\prime})\delta_{\rm D}^{[1]}(r-r^{\prime})$,444By applying
https://mathworld.wolfram.com/DeltaFunction.html equation 14 and chain rule.
we see that in both lines, we must have $r=r^{\prime}$. In view of equations
(14) and (17), this enables us to evaluate $F_{103}=1=F_{{}^{\prime}103}$,
which in turn renders the second line zero. Simplifying, we find

$\displaystyle I_{00}^{[2]}(r,r^{\prime})=\frac{\pi}{2rr^{\prime}}\delta_{\rm
D}^{[1]}(r-r^{\prime}).$ (32)

This recovers the closure relation (7) as desired.

###  5.2 $\ell=\ell^{\prime},\;n=0$

Just as with the above case, since we have $\ell=\ell^{\prime}$ and $n=0$, all
terms in equation (4) containing $d_{103}$ or $d_{{}^{\prime}103}$ will
vanish. We also set $n=0$. We find:

$\displaystyle
I_{00}^{[0+2]}(r,r^{\prime})=-\frac{\pi}{4}\Bigg{\\{}{\color[rgb]{1,0,0}H(r-r^{\prime})c_{\Gamma
r^{\prime}>}\bigg{(}\frac{\left(\ell^{\prime}+\ell^{\prime
2}-\ell^{2}-\ell\right)r^{\prime\ell}}{r^{\ell^{\prime}+3}}F_{103}\bigg{)}}$
$\displaystyle\;+{\color[rgb]{0,0.88,0}\Bigg{[}}\delta_{\rm
D}^{[1]}(r-r^{\prime})\bigg{(}{\color[rgb]{0,0,1}c_{\Gamma
r^{\prime}>}\left(-\frac{(\ell+1)r^{\ell-1}}{r^{\prime\ell+1}}F_{{}^{\prime}103}\right)}+{\color[rgb]{1,0,0}c_{\Gamma
r>}\left(-\frac{(\ell^{\prime})r^{\prime\ell}}{r^{\ell^{\prime}+2}}F_{103}\right)}\bigg{)}$
$\displaystyle\;+\frac{\partial}{\partial r}\left[\frac{1}{2}\delta_{\rm
D}^{[1]}(r-r^{\prime})\right]\left({\color[rgb]{0,0,1}-\frac{c_{\Gamma
r^{\prime}>}r^{\ell}}{r^{\prime\ell+1}}F_{{}^{\prime}103}}+{\color[rgb]{1,0,0}\frac{c_{\Gamma
r>}r^{\prime\ell}}{r^{\ell^{\prime}+1}}F_{103}}\right){\color[rgb]{0,0.88,0}\Bigg{]}}\Bigg{\\}}.$
(33)

Setting $\ell=\ell^{\prime}$ in the first line makes it vanish. In the second
line, we use that the Delta function requires $r=r^{\prime}$ and hence
$F_{{}^{\prime}103}=F_{103}=1$, and $c_{\Gamma r^{\prime}>}=c_{\Gamma r>}$. In
the third line, we again use $c_{\Gamma r>}=c_{\Gamma r^{\prime}>}$ and also
the comment about the Delta function’s derivative made in §5.1, which means
$r=r^{\prime}$ and $F_{103}=F_{{}^{\prime}103}$. With these substitutions the
third line drops out. We have

$\displaystyle I_{00}^{[0+2]}(r,r^{\prime})=\frac{\pi}{4}\delta_{\rm
D}^{[1]}(r-r^{\prime})\bigg{(}\frac{(2\ell+1)c_{\Gamma r>}}{r^{2}}\bigg{)}.$
(34)

Using now that $c_{\Gamma r>}=c_{\Gamma
r^{\prime}>}=\Gamma(\ell+1/2)/\Gamma(\ell+3/2)=2/(2\ell+1)$, and replacing one
$r$ with $r^{\prime}$ for symmetry, we find

$\displaystyle
I_{00}^{[0+2]}(r,r^{\prime})=\frac{\pi}{2rr^{\prime}}\delta_{\rm
D}^{[1]}(r-r^{\prime})$ (35)

as desired.

###  5.3 $\ell=0,\;\ell^{\prime}=1,\;n=0$

For this case, the integral resulting after raising the power is captured by
equation (19) (a splicing together of two formulae from [50]), as the result
is non-singular. However, the reduction is involved, and rather than showing
it, we simply allow the accompanying Mathematica notebook to perform it, and
plot the results, comparing to the analytic expression in [50]. We note that
none of the subtleties discussed in the footnote below equation (4) regarding
Mathematica’s treatment of the Heavisides or Delta functions arise, as these
subtleties only enter the singular terms and there are none here. We display
results for this case in Figure 1.

Figure 1: Comparison of using our formula (left) vs. using Gradshteyn & Ryzhik
([50]) directly (right). The axes are $r$ (horizontal) and $r^{\prime}$
(vertical) and we have taken the logarithm of the absolute value to show the
dynamic range of these integrals. Overall, one can see that our approach
agrees quite well with the known Gradshteyn & Ryzhik [50] result (given by our
equation (19)); we do not show a ratio to avoid dividing by numbers near zero.
It is worth observing the lack of reflection symmetry about the diagonal; this
is because there is no interchange symmetry between $r$ and $r^{\prime}$, as
we expect since they are the free arguments of sBFs of different orders
(respectively, zero and one).

##  6 Triple-SBF Integral Result

We first outline our overall approach, then present a result from the
literature that we will need, and end with several examples showing the
utility of our method.

###  6.1 Overall Approach

Consider an integral of three SBFs. We have

$\displaystyle I_{\ell_{1}\ell_{2}\ell_{3}}^{[n]}(r_{1},r_{2},r_{3})\equiv\int
k^{n}dk\;j_{\ell_{1}}(kr_{1})j_{\ell_{2}}(kr_{2})j_{\ell_{3}}(kr_{3}).$ (36)

If $n=2$ and the sum $\ell_{1}+\ell_{2}+\ell_{3}$ is even, then this integral
may be done and results in a polynomial, as can be seen in examples in [15]
and as we prove in the next subsection. We seek a fully general formula.
First, let us rewrite

$\displaystyle k^{n-2}j_{\ell_{3}}(kr_{3})=\int dk^{\prime}\;k^{\prime
n-2}j_{\ell_{3}}(k^{\prime}r_{3})\delta_{\rm
D}^{[1]}(k-k^{\prime})=\frac{2}{\pi}\int k^{\prime
n}j_{\ell_{3}}(k^{\prime}r_{3})\int u^{2}du\;j_{L}(uk)j_{L}(uk^{\prime}).$
(37)

Inserting this result in equation (36) we have

$\displaystyle
I_{\ell_{1}\ell_{2}\ell_{3}}^{[n]}(r_{1},r_{2},r_{3})=\frac{2}{\pi}\int
u^{2}du\;\int k^{\prime
n}dk^{\prime}\;j_{\ell_{3}}(k^{\prime}r_{3})j_{L}(k^{\prime}u)\;\int
k^{2}dk\;j_{\ell_{1}}(kr_{1})j_{\ell_{2}}(kr_{2})j_{L}(ku).$ (38)

We now perform the inner integral over $k$. Since $L$ is free, we can choose
it so that $\ell_{1}+\ell_{2}+L$ is even and this inner integral can be done
(some examples are given in e.g. [15]), and results in a sum over terms each
of which just involves powers of $r_{1},r_{2}$ and $r_{3}$. We present and
discuss the detailed general result, obtained by [31], in the following
subsection. We then may use the methods outlined in earlier sections of the
present work to perform the integral over $k^{\prime}$ in closed form. It will
yield Heaviside functions, Delta functions and their derivatives, and
hypergeometrics.

We can easily integrate Heaviside functions and Delta functions and their
derivatives over $u$. Since any additional dependence on $u$ from our
innermost integral will be a sum of power-laws, we can also perform that
integral analytically using formulae for integrals of power-laws against
hypergeometrics such as NIST DLMF 15.4.1555https://dlmf.nist.gov/15.14 (the
integral is just its Mellin transform) or, if a smoothing factor is needed for
convergence to then have a limit taken at the end, the formula found in e.g.
Wolfram.666https://functions.wolfram.com/HypergeometricFunctions/Hypergeometric2F1/21/02/02/
This point will become more concrete when we move, shortly, to presenting
several examples.

Finally, we note that this trick of introducing a Delta function is related to
the method of [51], as that work’s equation (2) is another form of it though
one that avoids explicit appeal to the Delta function. [51] exploits this
trick to prove a number of interesting and complicated integrals. [52] obtains
further results using this idea.

###  6.2 Result for Even Sum of the Orders and Square Weight

[31] equation (3.21) gives the result for an even sum of the three SBF orders
($\ell_{1}+\ell_{2}+L=$ even) and a weight of $k^{2}dk$ in the integral, as we
require to perform the innermost integral in equation (38). We duplicate the
result here as we will need to discuss some of its properties. We have

$\displaystyle I_{\ell_{1}\ell_{2}L,\;{\rm even}}^{[2]}$
$\displaystyle=\frac{\pi\beta(\Delta)}{4r_{1}r_{2}u}i^{\ell_{1}+\ell_{2}-L}\sqrt{2L+1}\left(\frac{r_{1}}{u}\right)^{L}\begin{pmatrix}\ell_{1}&\ell_{2}&L\\\
0&0&0\end{pmatrix}^{-1}\sum_{\mathcal{L}=0}^{L}{{2L}\choose{2\mathcal{L}}}^{1/2}\left(\frac{r_{2}}{r_{1}}\right)^{\mathcal{L}}$
(39)
$\displaystyle\qquad\times\sum_{\ell}(2\ell+1)\begin{pmatrix}\ell_{1}&L-\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}&L\\\
\mathcal{L}&L-\mathcal{L}&\ell\\\ \end{matrix}\right\\}P_{\ell}(\Delta).$

We may see from the inverse of the 3-$j$ symbol in the $\ell_{i}$ that the
result of [31] is only valid when their sum is even; otherwise this 3-$j$
symbol vanishes and hence its inverse is ill-defined. We note that the
$2\times 3$ matrix in curly brackets is a Wigner 6-$j$ symbol (see e.g. NIST
DLMF §34.4777https://dlmf.nist.gov/34.4).

Equation (3.13) of [31] defines $\beta$ as

$\displaystyle\beta(\Delta)\equiv\bar{\theta}(1-\Delta)\bar{\theta}(1+\Delta)$
(40)

and equation (3.14) of [31] gives the modified step function $\bar{\theta}$ as

$\displaystyle\bar{\theta}(y)=\left\\{0,\;y<0;\;\;\;1/2,\;y=0;\;\;\;1,\;y>0\right\\}.$
(41)

$\Delta$ is ([31] equation 3.9)

$\displaystyle\Delta\equiv\frac{r_{1}^{2}+r_{2}^{2}-u^{2}}{2r_{1}r_{2}}.$ (42)

We may see from equation (39) that the ratios $(r_{1}/u)$ and $(r_{2}/r_{2})$
will only ever give rise to terms that behave like powers of each $r_{i}$. We
may see from the explicit representation of the Legendre polynomial $P_{\ell}$
as a finite series of powers of $\Delta$ descending in same-parity steps from
$\ell$ (e.g. Wolfram equation
32888https://mathworld.wolfram.com/LegendrePolynomial.html) that, upon
expanding these powers of $\Delta$, we will simply have a sum over products of
powers of the $r_{i}$. Consequently, our integral $I_{\ell_{1}\ell_{2}L,\;{\rm
even}}^{[2]}$ will always be a sum of terms each of which is a product of
powers of the $r_{i}$, as claimed. The only additional complication is
$\beta$, which involves a modified step function constraining the relationship
between the $r_{i}$.

We might worry that we could not integrate the hypergeometrics that may result
from the $k^{\prime}$ integral in equation (38) against the result (39) for
the $k^{2}dk$ integral, which will involve $\beta$. However, we have results
for indefinite integrals of the hypergeometric against power-laws, so this
integral may always be performed. In particular, we have

$\displaystyle\int z^{\alpha-1}\;_{2}{\rm
F}_{1}(a,b,c;z)dz=\frac{z^{\alpha}}{\alpha}\;_{3}{\rm
F}_{2}(a,b,\alpha;c,\alpha+1;z)$ (43)

from e.g.
Wolfram.999https://functions.wolfram.com/HypergeometricFunctions/Hypergeometric2F1/21/01/02/01/,
second formula. As this is an indefinite integral, there is no restriction on
$\alpha$ so we can also deal with negative powers if required. For some terms
in the result (39) this may be necessary as both the ratios and the Legendre
polynomial in $\Delta$ bring in negative powers of $r_{1},r_{2}$, or $u$.

Finally, we note that [48], though it claims to be very general as to the
three-SBF cases it treats, actually does not indicate how to do a triple-SBF
integral where the sum of the three SBF orders and the power is even; there
the cosine in [48]’s equation (9) vanishes and it is not well-described how to
proceed. Hence, though the [48] might in principle give a simpler result than
[31] if it could be applied, in practice it is in fact the [31] result that is
most useful for this case.

###  6.3 Examples

To demonstrate the utility of our approach, we here give a few examples that
would not be obtainable by the other methods previously in the literature.

####  6.3.1 $n=2,\;\ell_{1}=0=\ell_{2},\;\ell_{3}=1$

We seek to evaluate the integral

$\displaystyle I_{001}^{[2]}(r_{1},r_{2},r_{3})=\int
k^{2}dk\;j_{0}(kr_{1})j_{0}(kr_{2})j_{1}(kr_{3}).$ (44)

We apply equation (38). We first set $L=0$ in the innermost integral of that
relation, and to evaluate it then use [30] Appendix B, first unnumbered result
below their equation B3:

$\displaystyle\int
k^{2}dk\;j_{0}(kr_{1})j_{0}(kr_{2})j_{0}(ku)=\frac{1}{4}\frac{\pi\beta(\Delta)}{r_{1}r_{2}u},$
(45)
$\displaystyle\qquad\qquad\Delta=\frac{r_{1}^{2}+r_{2}^{2}-u^{2}}{2r_{1}r_{2}}.$

Our middle integral, over $k^{\prime}$, is then

$\displaystyle\int k^{\prime
2}dk^{\prime}\;j_{1}(k^{\prime}r_{3})j_{0}(k^{\prime}u)$
$\displaystyle=H(u-r_{3})\left[\frac{1}{r_{3}^{3}-r_{3}u^{2}}+\frac{{\rm
arctanh}[r_{3}/u]}{r_{3}^{2}u}\right]$
$\displaystyle+H(r_{3}-u)\left[\frac{2u}{3(r_{3}^{3}-r_{3}u^{2})}+\frac{2\;{\rm
arctanh}[u/r_{3}]}{3r_{3}^{2}}\right]$ (46)

We note that this integral is a case that is marginally convergent as the out-
of-phase behavior of the SBFs in the UV (e.g. see equation (6)) converts the
otherwise-constant asymptotic (as $k^{2}\times(1/k)^{2}$) into an oscillating
function. Hence, though it can be evaluated with our operator method, it also
can be taken directly from [50]. Nonetheless this example is still a non-
trivial case where our method of reducing a triple SBF integral is required,
relative to what [31] would be able to evaluate.

The product of equations (45) and (46) must then be integrated against
$u^{2}du$, as equation (38) dictates, for the final result. Importantly, the
factor of $1/u$ in equation (45) means that, effectively, we are integrating
all of equation (46) against a weight of $u\,du$. We also note that due to the
factor of $\beta(\Delta)$ in equation (45), this final integral will be over a
finite rather than a half-infinite interval, but since results for the
indefinite integrals for even the most complicated piece in the integrand,
${\rm arctanh}$ times a power-law, are available, this is not problematic.
Indeed, since ${\rm arctanh}$ may be represented as a
hypergeometric$\;{}_{2}{\rm F}_{1}$, equation (43) gives the required result.

We may solve the restriction on $\Delta$ to find the bounds of integration,
denoted $u_{-}=r_{2}-r_{1}$ and $u_{+}=r_{2}+r_{1}$, and overall we have

$\displaystyle
I_{001}^{[2]}(r_{1},r_{2},r_{3})=\frac{2}{\pi}\left\\{H(u-r_{3})I_{u>}+H(r_{3}-u)I_{u<}\right\\}\bigg{|}_{u_{-}}^{u_{+}},$
(47)

with

$\displaystyle
I_{u>}=-\frac{1}{2}r_{3}^{-1}\ln[r_{3}^{2}-u^{2}]+\frac{1}{2r_{3}}r_{3}^{-2}\left[u\;{\rm
arctanh}[r_{3}/u]+\frac{1}{2}r_{3}\ln[r_{3}^{2}-u^{2}]\right],$ $\displaystyle
I_{u<}=2\left[-\frac{u}{3r_{3}}+\frac{1}{6}\ln\left[\frac{u+r_{3}}{u-r_{3}}\right]\right]+\frac{1}{6}r_{3}^{-2}\left[2u^{2}\;{\rm
arctanh}[r_{3}/u]+r_{3}\left(2u+r_{3}\ln\left[\frac{u-r_{3}}{u+r_{3}}\right]\right)\right].$
(48)

We have dropped the contribution of measure zero when $u=r_{3}$ and we also do
not need to track the $1/2$ in $\bar{\theta}$ (equation 41) within $\beta$
(equation 40) because it is a set of measure zero with respect to our $u$
integration. These results may now be straightforwardly evaluated at $u_{+}$
and $u_{-}$.

####  6.3.2 $n=4,\;\ell_{1}=0=\ell_{2},\;\ell_{3}=2$

We now seek to evaluate the integral

$\displaystyle I_{002}^{[4]}(r_{1},r_{2},r_{3})=\int
k^{4}dk\;j_{0}(kr_{1})j_{0}(kr_{2})j_{2}(kr_{3}).$ (49)

This integral is even-parity in the integrand, so the 3-$j$ symbol of [31]
would allow it, yet the power of $k$ is too high for that work’s formula to
apply; the work of [48] is not applicable because the integral is even-parity
and also not convergent.

First, we return to equation (38) and notice that we may use the same $L$ and
innermost integral as we did in §6.3.1. In this case our middle integral, over
$k^{\prime}$, is

$\displaystyle I_{\rm mid}\equiv\int k^{\prime
4}dk^{\prime}\;j_{2}(k^{\prime}r_{3})j_{0}(k^{\prime}u)$
$\displaystyle=\hat{D}_{0,u}\int k^{\prime
2}dk^{\prime}\;j_{2}(k^{\prime}r_{3})j_{0}(k^{\prime}u)$
$\displaystyle=\frac{3\pi}{2r_{3}^{3}}\left[\frac{1}{u}\delta_{\rm
D}^{[1]}(r_{3}-u)-\frac{1}{2}\frac{\partial}{\partial u}\delta_{\rm
D}^{[1]}(r_{3}-u)\right].$ (50)

We note that the operator with respect to $u$ must be applied only to this
integral, and not act on the innermost integral, the result of which contains
$u$-dependence; by the same token the operator must be inside the outermost
integral over $u$. Multiplying this result with equation (45) for our
innermost integral, we have

$\displaystyle I_{002}^{[4]}(r_{1},r_{2},r_{3})$
$\displaystyle=\frac{1}{2}\int
u^{2}du\;\frac{\beta(\Delta)}{r_{1}r_{2}u}I_{\rm mid}$
$\displaystyle=\frac{3\pi}{4r_{1}r_{2}r_{3}^{3}}\int_{u_{-}}^{u_{+}}udu\;\left[\frac{1}{u}\delta_{\rm
D}^{[1]}(r_{3}-u)-\frac{1}{2}\frac{\partial}{\partial u}\delta_{\rm
D}^{[1]}(r_{3}-u)\right].$ (51)

The integral above is trivial to perform, and we note that it will vanish
unless $r_{3}$ satisfies a triangular inequality with $r_{1}$ and $r_{2}$;
otherwise the value of $u$ picked out by the Delta functions, $u=r_{3}$, will
not be within the domain of integration.

####  6.3.3 $n=4,\;\ell_{1}=0=\ell_{2},\;\ell_{3}=1$

We now seek to evaluate the integral

$\displaystyle I_{001}^{[4]}(r_{1},r_{2},r_{3})=\int
k^{4}dk\;j_{0}(kr_{1})j_{0}(kr_{2})j_{1}(kr_{3}).$ (52)

This integral is odd-parity in the integrand and also not convergent, and
again, cannot be performed using previous techniques in the literature.

In view of equation (38), we observe that we may use the same innermost
integral as in the previous two subsubsections, and our middle integral is

$\displaystyle I_{\rm mid}=\int k^{\prime
4}dk^{\prime}\;j_{1}(k^{\prime}r_{3})j_{0}(k^{\prime}u).$ (53)

It may be computed by applying our operator to the known integral below on the
righthand side, as

$\displaystyle I_{\rm mid}=\hat{D}_{0,u}\int k^{\prime
2}dk^{\prime}\;j_{1}(k^{\prime}r_{3})j_{0}(k^{\prime}u).$ (54)

##  7 Dealing with More than Three SBFs

To resolve all integrals of more than three SBFs we may use a technique
suggested by [53]. Consider an integral

$\displaystyle
I_{\ell_{1}\ell_{2}\ell_{3}\ell_{4}}^{[n]}(r_{1},r_{2},r_{3},r_{4})\equiv\int
k^{n}dk\;j_{\ell_{1}}(kr_{1})j_{\ell_{2}}(kr_{2})j_{\ell_{3}}(kr_{3})j_{\ell_{4}}(kr_{4}).$
(55)

Using the same trick as in §6 we see that

$\displaystyle I_{\ell_{1}\ell_{2}\ell_{3}\ell_{4}}^{[n]}=\frac{2}{\pi}\int
u^{2}du\;\int k^{\prime
n}dk^{\prime}\;j_{\ell_{3}}(k^{\prime}r_{3})j_{\ell_{4}}(k^{\prime}r_{4})j_{L}(k^{\prime}u)\int
k^{2}dk\;j_{\ell_{1}}(kr_{1})j_{\ell_{2}}(kr_{2})j_{L}(ku).$ (56)

We make two observations. First, these integrals cannot just be done using the
techniques of [15] because i) one of them (that over $k^{\prime}$) has a
power-law other than $2$, and ii) we are not guaranteed that both sums
$\ell_{1}+\ell_{2}+L$ and $\ell_{3}+\ell_{4}+L$ will be even. One might hope
that at least ii) could be dealt with, but in fact one cannot choose an $L$
that makes both of these sums even unless the sums $\ell_{1}+\ell_{2}$ and
$\ell_{3}+\ell_{4}$ themselves have the same parity (see also discussion in
Slepian and Hou in prep.). This is not guaranteed. Thus, for both reasons i)
and ii) the technique outlined in §6 is required to enable performing an
arbitrary integral of four SBFs.

We now investigate how to split an integral of five SBFs. We have

$\displaystyle
I_{\ell_{1}\ell_{2}\ell_{3}\ell_{4}\ell_{5}}^{[n]}(r_{1},r_{2},r_{3},r_{4}),r_{5}\equiv\int
k^{n}dk\;j_{\ell_{1}}(kr_{1})j_{\ell_{2}}(kr_{2})j_{\ell_{3}}(kr_{3})j_{\ell_{4}}(kr_{4})j_{\ell_{5}}(kr_{5}).$
(57)

We use the trick of §6 twice to find

$\displaystyle I_{\ell_{1}\ell_{2}\ell_{3}\ell_{4}\ell_{5}}^{[n]}$
$\displaystyle=\left(\frac{2}{\pi}\right)^{2}\int u^{\prime
2}du^{\prime}\;\int u^{2}du\;\int
k^{n}dk\;j_{\ell_{1}}(kr_{1})j_{L}(ku)j_{L^{\prime}}(ku^{\prime})$
$\displaystyle\times\int
dk^{\prime}\;j_{\ell_{2}}(k^{\prime}r_{2})j_{\ell_{3}}(k^{\prime}r_{3})j_{L}(k^{\prime}u)\;\int
dk^{\prime\prime}\;j_{\ell_{4}}(k^{\prime\prime}r_{4})j_{\ell_{5}}(k^{\prime\prime}r_{5})j_{L^{\prime}}(k^{\prime\prime}u^{\prime}).$
(58)

We now see three triple-SBF integrals that we can perform using the results of
§6. Furthermore, they will involve integrals of at most hypergeometrics, and
these integrals, by change of variable, can be seen to result in powers of
e.g. $u$ and $u^{\prime}$; hence the final outer integrals can be done in
closed form too.

We note that if one wished to perform an integral of six SBFs, one could use
our trick of §6 just once; this would give the analog of equation (56) but
with four SBFs each in the integrals over $k^{\prime}$ and $k$ (imagine adding
$j_{\ell_{5}}(k^{\prime}r_{5})$ to the product within the first and
$j_{\ell_{6}}(kr_{6})$ to that within the second, for instance). But we have
already shown that any integral over four SBFs can be reduced as in equation
(56); hence the result of the reduction on the six-SBF integral can be
evaluated. Similar logic shows that an arbitrary number of SBFs in the
integrand may be dealt with in this fashion.

##  8 Concluding Discussion

In this work, we have presented a solution to the overlap integral of any pair
of spherical Bessel functions with any arbitrary (integer and positive semi-
definite) power-law weight. This work has applications in cosmology as more
fully outlined in §1. We have also shown how the ability to do any such double
integral enables performing any triple-sBF integral, exploiting a technique
originally due to [53] but which requires the ability to do an arbitrary
double-sBF integral to enable an arbitrary (positive semi-definite integer
power-law weight) triple one. Once any triple sBF integral can be done,
further use of this Mehrem technique enables doing any quadruple, quintuple,
etc., as we outline. An intriguing direction of future work would be to extend
these approaches to integrals with negative power-law weights, i.e. $k^{-1}$,
etc. We believe this could be done beginning from the integrals in the present
work using a technique of [54], wherein one performs “parametric integration”,
to lower the power-law. Also of interest would be further investigating the
reason for the singularity of some of these overlap integrals but not others;
complementary perspectives may come from considering contour integration as an
approach, as in [40], and also from analyzing the asymptotic behavior of the
phase (e.g. as in equation (6)) at large argument [Cahn & Slepian in prep.].

## Acknowledgments

We thank J. Chellino, F. Kamalinejad, A. Krolewski, O.H.E. Philcox, M.
Simonovic, and M. White for many useful conversations, and especially R.N.
Cahn for much insight into singular sBF integrals over several years, and J.
Hou for help in the final stages of various reductions.

## References

  * [1] Zachary Slepian and Daniel J. Eisenstein.  Modelling the large-scale redshift-space 3-point correlation function of galaxies.  MNRAS, 469(2):2059–2076, August 2017.
  * [2] Obinna Umeh.  Optimal computation of anisotropic galaxy three point correlation function multipoles using 2DFFTLOG formalism.  JCAP, 2021(5):035, May 2021.
  * [3] Alejandro Aviles, Arka Banerjee, Gustavo Niz, and Zachary Slepian.  Clustering in massive neutrino cosmologies via Eulerian Perturbation Theory.  JCAP, 2021(11):028, November 2021.
  * [4] Zachary Slepian and Daniel J. Eisenstein.  On the signature of the baryon-dark matter relative velocity in the two- and three-point galaxy correlation functions.  MNRAS, 448(1):9–26, March 2015.
  * [5] A. J. S. Hamilton.  Uncorrelated modes of the non-linear power spectrum.  MNRAS, 312(2):257–284, February 2000.
  * [6] Emanuele Castorina and Martin White.  Beyond the plane-parallel approximation for redshift surveys.  MNRAS, 476(4):4403–4417, June 2018.
  * [7] Xiaoying Xu, Nikhil Padmanabhan, Daniel J. Eisenstein, Kushal T. Mehta, and Antonio J. Cuesta.  A 2 per cent distance to z = 0.35 by reconstructing baryon acoustic oscillations - II. Fitting techniques.  MNRAS, 427(3):2146–2167, December 2012.
  * [8] 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.
  * [9] 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.
  * [10] Jiamin Hou, Robert N. Cahn, Oliver H. E. Philcox, and Zachary Slepian.  Analytic Gaussian Covariance Matrices for Galaxy $N$-Point Correlation Functions.  arXiv e-prints, page arXiv:2108.01714, August 2021.
  * [11] Uros Seljak and Matias Zaldarriaga.  A Line-of-Sight Integration Approach to Cosmic Microwave Background Anisotropies.  ApJ, 469:437, October 1996.
  * [12] Marcel Schmittfull and Zvonimir Vlah.  Reducing the two-loop large-scale structure power spectrum to low-dimensional, radial integrals.  PRD, 94(10):103530, November 2016.
  * [13] Marcel Schmittfull, Zvonimir Vlah, and Patrick McDonald.  Fast large scale structure perturbation theory using one-dimensional fast Fourier transforms.  PRD, 93(10):103528, May 2016.
  * [14] Zachary Slepian.  On decoupling the integrals of cosmological perturbation theory.  MNRAS, 507(1):1337–1360, October 2021.
  * [15] Lucía Fonseca de la Bella, Donough Regan, David Seery, and Shaun Hotchkiss.  The matter power spectrum in redshift space using effective field theory.  JCAP, 2017(11):039, November 2017.
  * [16] Marko Simonović, Tobias Baldauf, Matias Zaldarriaga, John Joseph Carrasco, and Juna A. Kollmeier.  Cosmological perturbation theory using the FFTLog: formalism and connection to QFT loop integrals.  JCAP, 2018(4):030, April 2018.
  * [17] Nils Schöneberg, Marko Simonović, Julien Lesgourgues, and Matias Zaldarriaga.  Beyond the traditional line-of-sight approach of cosmological angular statistics.  JCAP, 2018(10):047, October 2018.
  * [18] Xiao Fang, Elisabeth Krause, Tim Eifler, and Niall MacCrann.  Beyond Limber: efficient computation of angular power spectra for galaxy clustering and weak lensing.  JCAP, 2020(5):010, May 2020.
  * [19] Valentin Assassi, Marko Simonović, and Matias Zaldarriaga.  Efficient evaluation of angular power spectra and bispectra.  JCAP, 2017(11):054, November 2017.
  * [20] Henry S. Grasshorn Gebhardt and Donghui Jeong.  Fast and accurate computation of projected two-point functions.  PRD, 97(2):023504, January 2018.
  * [21] Xiao Fang, Tim Eifler, and Elisabeth Krause.  2D-FFTLog: efficient computation of real-space covariance matrices for galaxy clustering and weak lensing.  MNRAS, 497(3):2699–2714, September 2020.
  * [22] Nick Hand, Yin Li, Zachary Slepian, and Uroš Seljak.  An optimal FFT-based anisotropic power spectrum estimator.  JCAP, 2017(7):002, July 2017.
  * [23] Zachary Slepian and Daniel J. Eisenstein.  A simple analytic treatment of linear growth of structure with baryon acoustic oscillations.  MNRAS, 457(1):24–37, March 2016.
  * [24] Alex Krolewski and Zachary Slepian.  Accelerating computation of the density-field filtering scale $\sigma$(R) and non-linear mass by an order of magnitude.  Monthly Notices of the Royal Astronomical Society, 500(4):4439–4447, 11 2020.
  * [25] Jolyon K. Bloomfield, Stephen H. P. Face, and Zander Moss.  Indefinite Integrals of Spherical Bessel Functions.  arXiv e-prints, page arXiv:1703.06428, March 2017.
  * [26] Dr Teboho, Moloi.  Spherical Bessel functions.  arXiv e-prints, page arXiv:2102.02634, February 2021.
  * [27] W. N. Bailey.  Some infinite integrals involving bessel functions.  Proceedings of the London Mathematical Society, s2-40(1):37–48, 1936.
  * [28] G. N. Watson.  A Treatise on the Theory of Bessel Functions.  1966\.
  * [29] Rami Mehrem.  The Plane Wave Expansion, Infinite Integrals and Identities involving Spherical Bessel Functions.  arXiv e-prints, page arXiv:0909.0494, September 2009.
  * [30] R. Mehrem and A. Hohenegger.  A generalization for the infinite integral over three spherical Bessel functions.  Journal of Physics A Mathematical General, 43(45):455204, November 2010.
  * [31] 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.
  * [32] R. Mehrem.  An Analytical Evaluation For The Integral Of Two Spherical Bessel Functions With An Additional Exponential And Polynomial Factor.  arXiv e-prints, page arXiv:1110.6147, October 2011.
  * [33] Rami Mehrem.  The Integral Over 2 Spherical Bessel Functions Multiplied by a Gaussian.  arXiv e-prints, page arXiv:1908.07374, August 2019.
  * [34] AD Jackson and LC Maximon.  Integrals of products of bessel functions.  SIAM Journal on Mathematical Analysis, 3(3):446–460, 1972.
  * [35] A Gervois and H Navelet.  Differentiation method: a way to get new integrals with three bessel functions; application to spherical bessel functions.  Technical report, 1985.
  * [36] Annie Gervois and Henri Navelet.  Integrals involving three bessel or modified bessel functions.  Technical report, Commissariat a l’energie atomique et aux energies alternatives-CEA, 1986.
  * [37] Colm T Whelan.  On the evaluation of integrals over three spherical bessel functions.  Journal of Physics B: Atomic, Molecular and Optical Physics, 26(23):L823, 1993.
  * [38] Limin Wang and Marc Kamionkowski.  Cosmic microwave background bispectrum and inflation.  Phys. Rev. D, 61(6):063504, March 2000.
  * [39] R. Friedberg and O. Martin.  Fast evaluation of lattice Green functions.  Journal of Physics A Mathematical General, 20(15):5095–5111, October 1987.
  * [40] Jian-Xing Chen and Jun-Chen Su.  Glueball spectrum based on a rigorous three-dimensional relativistic equation for two-gluon bound states. II. Calculation of the glueball spectrum.  Phys. Rev. D, 69(7):076003, April 2004.
  * [41] J. R. Fergusson and E. P. S. Shellard.  Primordial non-Gaussianity and the CMB bispectrum.  Phys. Rev. D, 76(8):083523, October 2007.
  * [42] A Gervois and H Navelet.  Infinite integrals involving three spherical bessel functions.  SIAM journal on mathematical analysis, 20(4):1006–1018, 1989.
  * [43] Leonard C. Maximon.  On the evaluation of the integral over the product of two spherical Bessel functions.  Journal of Mathematical Physics, 32(3):642–648, March 1991.
  * [44] S.N. Stuart.  Non-classical integrals of bessel functions.  J. Austral. Math. Soc. (Series B), 1981.
  * [45] D. E. Dominici, P. M. W. Gill, and T. Limpanuparb.  A remarkable identity involving Bessel functions.  Proceedings of the Royal Society of London Series A, 468(2145):2667–2681, September 2012.
  * [46] Oliver H. E. Philcox and Zachary Slepian.  An exact integral-to-sum relation for products of Bessel functions.  Proceedings of the Royal Society of London Series A, 477(2253):20210376, September 2021.
  * [47] Satya N. Majumdar and Emmanuel Trizac.  When Random Walkers Help Solving Intriguing Integrals.  Phys. Rev. Lett., 123(2):020201, July 2019.
  * [48] V.I. Fabrikant.  Computation of infinite integrals involving three bessel functions by introduction of new formalism.  ZAMM - Journal of Applied Mathematics and Mechanics / Zeitschrift für Angewandte Mathematik und Mechanik, 83(6):363–374, 2003.
  * [49] Gregory S. Adkins.  Three-dimensional Fourier transforms, integrals of spherical Bessel functions, and novel delta function identities.  arXiv e-prints, page arXiv:1302.1830, February 2013.
  * [50] I. S. Gradshteyn, I. M. Ryzhik, Alan Jeffrey, and Daniel Zwillinger.  Table of Integrals, Series, and Products.  2007\.
  * [51] Howard S. Cohl.  Definite integrals using orthogonality and integral transforms.  Symmetry, Integrability and Geometry: Methods and Applications, Oct 2012.
  * [52] Howard S. Cohl, Sean J. Nair, and Rebekah M. Palmer.  Some dual definite integrals for bessel functions, 2015.
  * [53] R. Mehrem.  Analytical Evaluation Of An Infinite Integral Over Four Spherical Bessel Functions.  arXiv e-prints, page arXiv:1012.6013, December 2010.
  * [54] Roman Tomaschitz.  Weber and Beltrami integrals of squared spherical Bessel functions: finite series evaluation and high-index asymptotics.  Mathematical Methods in the Applied Sciences, 37(9):1249–1272, June 2014.


# High speed shadowgrpah of a L/D cavity at Mach 0.7 and 1.5

Ryan Schmit, Frank Semmelmayer, Mitch Haverkamp and James Grove

Air Force Research Laboratory,
Wright-Patterson Air Force Base, OH 45433, USA

###### Abstract

This is article highlights the fluid dynamics video of a rectangular cavity
with an L/D of 5.67 at Mach 0.7 and 1.5.

A rectangular cavity with an L/D of 5.67 at Mach 0.7 and 1.5, Reynolds number
$2x10^{6}$ and $2.3x10^{6}$ respectively, was examined using high speed
shadowgraph imaging. The cavity motion is shown in the video. The three movies
clips presented were sampled at at 75kHz and played back at 20 Hz. The camera
shutter-speed was 0.37$\mu$sec. The first clip shows the side and top view of
the cavity at Mach 0.7. Note that these two views are not in sync. The second
clip shows the side and top view of the cavity at Mach 1.5. Again the views
are not in sync. The third clip shows a zoomed out side view of the cavity at
Mach 1.5.

For more information please refer to Schmit, Semmelmayer, Haverkamp and Grove,
”Fourier Analysis of High Speed Shadowgraph Images around a Mach 1.5 Cavity
Flow Field”, 29th AIAA Applied Aerodyanmic Conference, Honolulu, Hi, pp 1-24,
2011, AIAA 2011-3961