Triple-spherical Bessel function integrals with exponential and Gaussian damping: towards an analytic N-point correlation function covariance model

Spherical Bessel functions (sBFs) appear commonly in many areas of physics wherein there is both translation and rotation invariance, and often integrals over products of several arise. Thus, analytic evaluation of such integrals with different weighting functions (which appear as toy models of a given physical observable, such as the galaxy power spectrum) is useful. Here, we present a generalization of a recursion-based method for evaluating such integrals. It gives relatively simple closed-form results in terms of Legendre functions (for the exponentially damped case) and Gamma, incomplete Gamma, and hypergeometric functions (for the Gaussian-damped case). We also present a new, non-recursive method to evaluate integrals of products of sBFs with Gaussian damping in terms of incomplete Gamma functions and hypergeometric functions.


Introduction
Spherical Bessel functions (sBFs) are ubiquitous in cosmology.The three-dimensional plane wave offers an eigenbasis for the gradient operators that appear in e.g.cosmological perturbation theory [1], and in general the plane waves are the eigenstates of momentum, which is the conserved quantity associated by Noether's theorem with the translation symmetry implied by the cosmological assumption of statistical homogeneity of the large-scale structure (LSS).However, cosmology also assumes isotropy (rotation symmetry) about any point, so spherical coordinates are often apt.They do not diagonalize the fluid equations as well (though efforts to solve these in the spherical basis have been made; e.g.[2]), but they are especially useful for numerical computation both of CMB anisotropies [3] and lensing [4] and LSS perturbative corrections [5,6,7,8,9].Efficient algorithms exist for performing sBF transforms numerically [10,11,12,13,14,15,16]. sBFs have also been explored as the "correct" basis for wide-angle redshift-space distortion (RSD) work [17] and, in a modified form, for LSS surveys with finite redshift width (leading to a spherical shell survey geometry) [18].
sBFs also commonly appear in other areas of physics.Mie scattering, in which electromagnetic waves are scattered off of a sphere, has scattering coefficients in terms of Riccati-Bessel functions [19,20,21].It is necessary to evaluate integrals of the product of sBFs in many scattering calculations (e.g.[22,23,24,25,26]).Such integrals may include Gaussian damping to ensure convergence at infinity.Modified Bessel functions often arise in the Ewald summation, which allows long-range interactions to be computed [27,28].Additionally, sBF integrals are relevant to spherical linear canonical transforms, which can be applied to studies of acoustics, optics and heat [29].
The plane wave expansion links the plane wave to the sBFs and spherical harmonics, and provides a powerful "Rosetta stone" for translating between situations where translation invariance is the more important symmetry and situations where rotation invariance is more germane.
While efficient numerical algorithms exist for transforms of one, two, and even three [7] sBFs against a source function given as a numerical lookup table, there are occasions where the integrand is given in closed form and thus an analytic evaluation of sBF integrals is of interest.For instance, in computations of the covariance matrix for pair correlations (2-Point Correlation Function, 2PCF), triplet correlations (3PCF, e.g.[30,31]), 4PCF (e.g.[32,33,34]), or beyond, a power law model of the galaxy power spectrum often provides insight into the gross properties of these large matrices.Such models for the 2PCF were explored in [35], for the 3PCF briefly in [36], §6.3, and at more length in [37], appendix D (the general integrals used for covariance of 3PCF, 4PCF and beyond).At the signal level, triple-sBF integrals were used prominently in the computation of the redshift-space galaxy power spectrum of [38], see in particular §3.2 and appendix B. This work used a method of parametric differentiation developed by [39], based on Rayleigh's formula for the sBFs.
There are numerous works on the integral of a known function against a single-sBF (a recent, incomplete list with fuller referencing within is [40,41]) and on the double-sBF case (more fully reviewed in [42]), but far fewer regarding triple-sBF integrals.Early work is in Watson's Treatise, [43], with further results in terms of Appell F4 functions in [44] (which simplify in some cases) and in [45] (reduction to products of associated Legendre polynomials of the external angles of a triangle formed by the free arguments) and [46,47], in terms of Legendre polynomials in the cosine of the internal enclosed angle of that triangle.These results are for a somewhat general power law ( [44], leading to the Appell F4, and simplifying if the power is a particular combination of the sBF orders) or a k 2 power law, i.e. the usual spherical coordinates Jacobian.Gervois and Navelet [48] pursues an integer power law and is subject to some restrictions on the sBF orders.There have been recent works in cosmology pursuing a variety of methods: [49,50], and most relevant for the current work, [51].This last work developed a recursion for triple-sBF integrals with a square power law weight, using the three-function recursion relation for the sBFs.
In the current work, we extend this recursion method to the case of a triple-sBF integral with an exponential damping and then with a Gaussian damping.Damped exponentials and Gaussians are both useful toy models for the galaxy power spectrum, and have the property that they lead to simple, finite 2PCFs, and that the integrals associated with the analytic covariance matrix formalism of [37] are all convergent.In a second paper, we explore an analytic covariance matrix model for the 2PCF-4PCF using these toy power spectra more fully.In the course of that work we found it necessary to develop methods to evaluate these integrals analytically, which we present here.

Recursion Relation Approach
A recursion relation has been used to evaluate triple-sBF integrals [51].Here we extend this recursion to evaluate such integrals including an exponential damping term.We will find that the result involves Legendre functions of the second kind.
We wish to evaluate integrals of the form for {p, ℓ i , r i } ∈ R, ℓ i ∈ Z, ℓ i ≥ 0, r 1 > 0, r i = 0, and p = 0, as will be explained below.We begin with the recursion relation [52] for spherical Bessel functions: where ℓ must be a positive integer.This recursion can be extended [51] to a recursion for triple-sBF integrals as To use this recursion to evaluate all cases, we require three base cases to anchor it: f ℓ,0,0 , f ℓ,−1,0 , and f ℓ,−1,−1 .The first base case is from setting ℓ 2 = 0 = ℓ 3 in equation (1).By representing the first spherical Bessel function in the integrand (order ℓ, argument kr 1 ) as in [46], we obtain which requires r 1 ≥ 0. The first base case becomes Each j 0 in equation ( 5) has been rewritten explicitly (with k in the denominator of each) to cancel the k 2 in the Jacobian; we require r 2 , r 3 = 0.After rewriting the sines in equation ( 6) as complex exponentials, we find As long as p = 0, we can use [53] equation 3.310 to evaluate the k integral, obtaining We then use [53] equation 7.224.1 to integrate over µ, resulting in where r 1 = 0, Q ℓ is the order-ℓ Legendre function of the second kind and R ±± is defined as We note that [53] equation 7.224.1 has a discontinuity when the argument of the Q ℓ is along the real line from −1 to +1, but the arguments R ±± we have of the Q ℓ above are manifestly complex ({p, r i } ∈ R) and hence do not lie along this line.
The other base cases can be obtained in the same way, resulting in and The recursion relation (3) may then be used to construct from these base cases the integral (1) for any other non-negative integer values of the ℓ i .We note that equations ( 9), (11), and (12) depend on the same four Legendre functions, although the signs in front of each term differ.The sign differences as well as the different complex pre-factors arise from writing sines and cosines as complex exponentials: equation ( 9) depends on the product of two sines, equation (11) the product of a sine and cosine, and equation ( 12) the product of two cosines.

Proof That the Base Cases Are Real
We now show that the results obtained for the three base cases are real; they should be real since the beginning integrands and domain of integration for each of them are real.We note that for each base case (equations 9, 11, and 12), the argument of the fourth Q ℓ is negative the complex conjugate of that of the first, and similarly for the third and the second.For f ℓ,0,0 and f ℓ,−1,−1 we therefore must show that ] is real (this identity applies to each pair of Q ℓ as above).For f ℓ,−1,0 , we need to demonstrate that ] is real (and again, this identity then applies to each pair of Q ℓ noted above).
We write z = re iθ and expand the Legendre function as i.e. as an infinite series.The more familiar representation is as a finite polynomial in z multiplied by ln[(1 + z)/(1 − z)], and it is this singular second factor that when expanded renders the power series infinite.For f ℓ,0,0 and f ℓ,−1,−1 we then have which for odd ℓ, by taking the even part of the exponential becomes and for even ℓ, by taking the odd part of the exponential becomes Similarly, for f ℓ,−1,0 we have which by taking the even part of the exponential (for odd ℓ) becomes and from the odd part of the exponential (for even ℓ) is These results are manifestly real, which shows that all three base cases in §2.1 must be real, as expected.

Recursion Relation Approach for Triple-sBF Integrals with Gaussian Damping
The integral of three sBFs with Gaussian damping can be evaluated by a method similar to that described in §2.1.To evaluate integrals of the form for {p, ℓ i , r i } ∈ R, ℓ i ∈ Z, ℓ i ≥ 0, r 1 > 0, r i = 0, and p = 0, the same recursion relation given in equation (3) [51] can be used with base cases f ℓ,0,0 , f ℓ,−1,0 , and f ℓ,−1,−1 .We will find that the base cases can be expressed using regularized hypergeometric and incomplete Gamma functions.The first base case is from setting ℓ 2 = 0 = ℓ 3 in equation (20): By using equation ( 5) to rewrite the j ℓ (kr 1 ), and writing the two factors of j 0 explicitly to cancel the k 2 in the Jacobian, we obtain where {r 2 , r 3 } = 0.After rewriting the sines as complex exponentials (obtaining the same set of exponentials as are in the square brackets in the integrand of equation 7) and completing the square on what results, this becomes where we define We note that the ordering of signs on r 2 and r 3 in the integrand matches that outside the integral in each line because the factor outside is what needs to be subtracted off after completing the square in each integrand.The inner k integrals can be evaluated using a change of variable x = k − (iR µ±± )/(2p 2 ) and [53] equation 2.33.16, as long as p = 0: We now must integrate over µ in equation (23), where in the first line we have a Legendre polynomial in µ.We expand this Legendre polynomial as1 this is simply a falling power series composed of powers of µ ranging from ℓ down to either one (if ℓ is odd) or zero (if ℓ is even) in even steps.We note that ⌊ℓ/2⌋ is the floor function, which gives the largest integer less than or equal to its argument; if ℓ is even, ℓ/2, but if ℓ is odd, (ℓ − 1)/2.
The first base case then becomes The µ integrals in equation ( 27) have the form After changing the integration variable from µ to R µ±± , then using the binomial expansion on µ ℓ−2n , we find for r 1 = 0: Multiplying the two terms in the square bracket by the Gaussian damping in the integrand of equation ( 29) and then, in what results, focusing on the term proportional to unity, we have from [53] equation 2.33.10 where Γ is the incomplete Gamma function and we have defined and Multiplying the two terms in the square bracket of equation ( 29) by the Gaussian damping in the integrand and then focusing on the term proportional to the error function, we obtain where Γ is the (complete) Gamma function and 2 F2 is the regularized hypergeometric function.We have defined The first base case can be written in its final form as where the arguments of G ±±± and χ (m) ±±± have been suppressed for legibility.The other two base cases are obtained by following the same steps outlined throughout this section, resulting in 2 https://functions.wolfram.com/GammaBetaErf/Erf/21/01/02/04/01/0007/and We note that the powers of the pairs of r i serving as pre-factors in each line above are always greater than or equal to zero; m ≤ ℓ − 2n, and ℓ − 2n ≥ 0 due to the floor function in the outer sum (which ultimately stems from expanding the Legendre polynomial in a series as in equation 26).The sign on each pair of r i stems from writing the j 0 (kr i ) and j −1 (kr i ) in terms of sine or cosine, then rewriting as complex exponentials.
Given these three base cases, the recursion relation ( 3) may now be used to construct the integral (20) for any non-negative integer values of the ℓ i .As in §2.1, the pre-factors of equations ( 35), (36), and (37) as well as the signs on each line differ because of the multiplication of sines and cosines as complex exponentials.

Stability of the Recursion Relations
For many applications within cosmology, such as covariance matrix calculations, the orders of two of the sBFs within a triple-sBF integral each have a maximum value ≤ 10 (e.g.[31,34,36,37,54]).The maximum value of the order of the third sBF is set by the sum of the orders of the other two sBFs; this is because of the triangle inequalities that stem from Wigner 3-j symbols coupling the sBF orders.Thus, it is desirable for the recursion relations of §2.1 and §3 to be stable for ℓ ≤ 20.
We duplicate the recursion for triple-sBF integrals (equation 3) below: f ℓ,0,0 , f ℓ,−1,0 , and f ℓ,−1,−1 are the three base cases needed to anchor this recursion.We will use an example to show how the integrals with higher-order sBFs are sums of integrals with lower-order sBFs.We begin by evaluating one of the lowest-order f -integrals that is not one of the base cases: f 1,1,0 .We use the recursion (equation 38) to write f 1,1,0 as The f 0,0,0 , f 2,0,0 , and f 1,−1,0 in equation ( 39) are all base cases.Thus, we see that f 1,1,0 can be evaluated by first determining only three base cases.We now consider f 1,1,1 , an f -integral with higher sBF orders than f 1,1,0 .Using the recursion In equation (40), f 0,0,1 is a base case; f 2,0,1 and f 1,−1,1 are not base cases and will each need to be evaluated in terms of base cases by using the recursion.f 2,0,1 is the same as f 2,1,0 ; switching the ordering of the sBFs within a triple-sBF integral does not affect the integral.We thus have The f 1,0,0 , f 3,0,0 , and f 2,−1,0 in equation ( 41) are all base cases.The factor (r 1 /r 2 ) in equation (38) has been replaced by (r 1 /r 3 ) in equation ( 41) since we have switched the ordering of the sBFs with arguments (kr 2 ) and (kr 3 ) when evaluating f 2,1,0 instead of f 2,0,1 .Finally, we must use the recursion to evaluate the f 1,−1,1 (equal to f 1,1,−1 ) in equation ( 40): In equation ( 42), the f 0,0,−1 , f 2,0,−1 , and f 1,−1,−1 are all base cases.Again, the factor (r 1 /r 3 ) in equation ( 42) arises from swapping the ordering of two sBFs within the desired integral.We now see that f 1,1,1 (equation 40) depends on seven base cases: one is f 0,0,1 (equation 40), three are from f 2,0,1 (equation 41), and three are from f 1,−1,1 (equation 42).This example shows that non-base case integrals are simply sums over some number of base case integrals.For example, the lowest order sBF to be used in the recursion is ℓ 1 = −1.If ℓ 2 and ℓ 3 are 0 or -1, the recursion is simply composed of base cases.By using the recursion along with these base cases, any triple sBF integral in the plane ℓ 1 = −1 with arbitrary {ℓ 2 , ℓ 3 } can be evaluated.This same strategy can be used for any other ℓ 1 > −1 plane: first evaluate the base cases when ℓ 2 and ℓ 3 are 0 or 1, then use the recursion to build up to integrals with higher-order ℓ 2 and ℓ 3 .Since swapping the order of any of the sBFs within a triple-sBF integral does not alter the integral, the same recursive technique can be used to evaluate integrals in the planes ℓ 2 ≥ −1 and ℓ 3 ≥ −1.
The stability of the recursion depends on whether a small error present in the base cases will propagate with each iteration of the recursion.As we use the recursion to evaluate higher-order triple-sBF integrals, we are simply adding more terms.Each new term may have a random error ǫ, which is independent of previous terms as it represents numerical error from evaluating the special functions that are required for the base cases.When adding terms in the recursion, the error from each will add up in quadrature, as is typical for a random walk.Thus, the stability of the recursion can be viewed as the expected error after N iterations of the recursion, which is √ N ǫ.To practically use the recursive work of this paper, we must address two issues: accurately evaluating each base case and co-adding as many base cases as needed without loss of numerical precision.This second issue is widespread in numerical work with many remedies for it.For instance, the mp.math library in Python offers arbitrary precision as a custom type.Balanced sum algorithms also exist to handle addition of values that significantly differ in order of magnitude, as could occur when co-adding a number of the base cases.Although it is possible to accurately evalulate the base cases, underflow may be present in the computed values.This underflow will build up throughout the recursion and can make the recursion unstable if high-precision results are needed.Additionally, sBF integrals solved using the spherical harmonic shift theorem (as in [55]) have a known numerical issue when one of the free arguments of the sBFs is much less than the other two free arguments.However, our approach does not require this shift theorem, so we do not expect to see any numerical issues with our results, similar to the explanation in [55] that the integrals within [45] do not have any numerical issues.We now turn to the first issue of accurately evaluating the base cases.

Stability of Recursion with Exponential Damping
The base cases (equations 9, 11, and 12) for the exponentially-damped triple-sBF integral ( §2.1) depend on Legendre functions of the second kind Q ℓ (R ±± ), with R ±± ≡ (−ip 2 ±r 2 ±r 3 )/r 1 previously defined in equation (10).These Legendre functions can be written in terms of the hyperbolic arctangent and powers of the argument R ±± ([53] equation 8.827), which will diverge when R ±± is along the real line from -1 to +1.Since R ±± is manifestly complex, Q ℓ (R ±± ) has no divergence.
We display the five lowest-order Legendre functions below: The hyperbolic arctangent is standardly implemented in numerical libraries such as scipy.Thus, the base cases for exponentially-damped triple-sBF integrals are stable.

Damped Triple-sBF Integrals Weighted by Higher Powers of k
In the preceding sections, we evaluated recursion relations for damped triple-sBF integrals involving the Jacobian k 2 dk.It was necessary to begin with k 2 so that the sBFs with orders 0 and -1 in the base cases, when rewritten as sines and cosines divided by their arguments kr i , would cancel the Jacobian.Recursions for triple-sBF integrals involving higher powers of k cannot be evaluated using the methods described in §2.1 and §3; however, they can be obtained by parametric differentiation.
We begin with the triple-sBF integral with exponential damping given by equation ( 1): for {p, ℓ i , r i } ∈ R, ℓ i ∈ Z, ℓ i ≥ 0, r 1 > 0, r i = 0 and p = 0.The partial derivative of f ℓ1,ℓ2,ℓ3 (r 1 , r 2 , r 3 ; p) with respect to −p 2 can be moved inside the integral since the derivative is with respect to −p 2 and the integral is with respect to k: where the partial derivative has resulted in one extra power of k in the integral since the only −p 2 dependence is in the exponential.Each successive derivative will add one to the power of k; therefore, this can be generalized to We can therefore evaluate the recursion for triple-sBF integrals with exponential damping for k to any power ≥ 2 by first evaluating the recursion for k 2 as done in §2.1 and then taking a partial derivative of the result with respect to −p 2 for each additional power of k desired.This parametric differentiation method can also be used to evaluate triple-sBF integrals with Gaussian damping involving k to any even power ≥ 2. We begin with equation ( 20): e −(pk) 2 j ℓ1 (kr 1 )j ℓ2 (kr 2 )j ℓ3 (kr 3 ), (50) for {p, ℓ i , r i } ∈ R, ℓ i ∈ Z, ℓ i ≥ 0, r 1 > 0, r i = 0, and p = 0.The partial derivative of f ℓ1,ℓ2,ℓ3 (r 1 , r 2 , r 3 ; p) with respect to −p 2 can be moved inside the integral since the derivative is with respect to −p 2 while the integral is with respect to k.Since the only −p 2 dependence is in the Gaussian, taking this derivative will increase the power of k by two: k 2 dk e −(pk) 2 j ℓ1 (kr 1 )j ℓ2 (kr 2 )j ℓ3 (kr 3 ) = ∞ 0 k 4 dk e −(pk) 2 j ℓ1 (kr 1 )j ℓ2 (kr 2 )j ℓ3 (kr 3 ).
This can be generalized to since each successive derivative with respect to −p 2 will increase the power of k by two.By evaluating the recursion for triple-sBF integrals with Gaussian damping involving the Jacobian k 2 dk ( §3) and then taking the necessary amount of partial derivatives with respect to −p 2 of the result, we can evaluate triple-sBF integrals with Gaussian damping weighted by k to any even power ≥ 2.

A Second Method to Evaluate Triple-sBF Integrals with Gaussian Damping
We now show an alternate, non-recursive method to evaluate triple-sBF integrals with Gaussian damping weighted by any non-negative integer power of k.The result will depend upon incomplete Gamma functions and hypergeometric functions.We term this the "Hankel-Bowman" method since it uses a result from Bowman [61] derived from Hankel's contour integral.We begin with for {p, ℓ i , r i } ∈ R, {ℓ i , n} ∈ Z, ℓ i ≥ 0, r i = 0, p = 0 and n ≥ 0. The spherical Bessel functions are related to cylindrical Bessel functions J ℓ (x) by Hankel's contour integral [61] can be used to rewrite the cylindrical Bessel functions as for ℓ > −(1/2).We first use equation ( 54) to write each sBF in equation ( 53) as a cylindrical Bessel function, then use equation ( 55) to rewrite each cylindrical Bessel function as an integral over q i , and change the order of integration so that the k integral will be performed first: where we have defined L ≡ ℓ 1 + ℓ 2 + ℓ 3 and The inner k integral, which will be defined as I k below, will be performed first: We first complete the square in the Gaussian (to obtain the second line of equation 58), then use a change of variable x = k − [(is 123 )/(2p 2 )] (to obtain the third line), and finally use a binomial expansion on the polynomial within the integrand (to obtain the fourth line).The integral in the last line of equation ( 58) is known in closed form in terms of an incomplete gamma function; we thus obtain3 We have defined α ≡ n + L − ζ.Replacing the innermost, k integral in equation ( 56) using equation (59), equation ( 56) becomes The upper incomplete gamma function is the difference between the complete gamma function and the lower incomplete gamma function, γ: The lower incomplete gamma function can be rewritten in terms of the confluent hypergeometric function of the first kind so that equation (61) becomes5 After replacing the upper incomplete gamma function in equation ( 60) using equation (62), we obtain where we have defined two new integrals and I pl depends on a power law and I hg depends on a hypergeometric function.To obtain equation (65), we used ζ + α + 1 = n + L + 1.
We begin with the integral I pl .We need to recall that s 123 depends on q 3 (equation 57).To evaluate the q 3 integral in equation (64), we make the change of variable x = −[s 123 /(2p)] 2 , then use 3 ) ℓ3 in terms of x and use the binomial expansion on it: where we have defined We evaluate the innermost integral with respect to x using [53] equation 2.33.10.Equation (66) then becomes To evaluate the innermost, q 2 integral, we make the change of variable x = [(s 12 ± r 3 )/(2p)] 2 , where the ± sign represents the sign of r 3 in the arguments of the incomplete gamma functions in equation ( 68).After making this change of variable, we use a binomial expansion on the resulting factor of s c 12 = (2p 2 ) ℓ2 , then use another binomial expansion on this factor to obtain where we have defined ǫ ≡ [c + 2d − g − h − 1]/2 and We now integrate over x, resulting in 6 where we have defined (with arguments suppressed) 6 https://functions.wolfram.com/GammaBetaErf/Gamma2/21/01/02/01/0001/and To evaluate the q 1 integral of equation ( 71), we make the change of variable x = [r q±± /(2p)] 2 , which corresponds to the arguments of the incomplete gamma functions in equation (71).We use ) ℓ1 , then use a binomial expansion on it.We use another binomial expansion on the factor of (−r 1 q 1 ∓ r 3 ) g in the third and fourth lines of equation ( 71): where we have defined η ≡ [g + 2m − t − u − v − 1]/2 and By evaluating the integrals of equation ( 74), we obtain a final result for I pl :7 where we have defined and We have suppressed the arguments of Z ±±± and Y ±±± .We now turn to I hg , given in equation ( 65) and duplicated below: To evaluate the innermost, q 3 integral, we must recall that s 123 depends on q 3 (equation 57).We make the change of variable x = −[s 123 /(2p)] 2 , which is the argument of the hypergeometric function in equation (79).Recalling the definition of s 12 from equation (67), we use q 3 = [2pi √ x − s 12 ]/r 3 to rewrite the (1 − q 2 3 ) ℓ3 factor in equation (79) before using a binomial expansion on it: where we have defined κ ≡ [n + L + 2b − c]/2.We now integrate over x, resulting in 8 To evaluate the innermost, q 2 integral of equation ( 81), we use the change of variable x = −[(s 12 ± r 3 )/(2p)] 2 , where the ± sign represents the sign of r 3 in the arguments of the hypergeometric functions in equation (81).Since we can rewrite the s c 12 factor in equation (81) as s c 12 = (2ip √ x ∓ r 3 ) c , then use a binomial expansion on it.We also use a binomial expansion on the (1 − q 2 2 ) ℓ2 factor of equation (81) after writing it in terms of x: where we recall r q±± ≡ r 1 q 1 ± r 2 ± r 3 and define where we have defined (with arguments suppressed) To evaluate the q 1 integral, we use the change of variable x = −[r q±± /(2p)] 2 .We use ℓ1 factor in equation ( 83), then use a binomial expansion on it.We also rewrite the (−r 1 q 1 ∓ r 3 ) g factor of equation ( 83) in terms of x, then use another binomial expansion on this factor: where we recall r ±±± ≡ ±r 1 ± r 2 ± r 3 and define φ ≡ where we have defined The arguments of F ±±± have been suppressed.The integral of three sBFs with a Gaussian was given in terms of I pl and I hg in equation ( 63); we duplicate this below: We now use the results for I pl (equation 76) and I hg (equation 86) to evaluate equation (88): where we have defined Equation ( 89) is our final result for the triple-sBF integral with Gaussian damping weighted by any non-negative integer power of k, given by equation (53).This final result is in terms of incomplete Gamma functions and hypergeometric functions; the base cases of the recursion relation for the triple-sBF integral with Gaussian damping in §3 also depend on the same types of functions.Additionally, both these base cases and equation (89) have a similar structure in front of the incomplete Gamma and hypergeometric functions: powers of (∓r 2 ∓ r 3 ).
The recursion relation in §3 was limited to triple-sBF integrals with Gaussian damping and the Jacobian k 2 dk.In §5, it was shown that parametric differentiation can be used to extend this recursion to integrals weighted by k to any even power ≥ 2. The result given in equation ( 89) allows more freedom in the power of k; with this result, we are able to evaluate triple-sBF integrals with Gaussian damping weighted by k to any non-negative integer power ≥ 0.
The steps outlined in this section may be used to evaluate integrals of any number of sBFs with Gaussian damping.The sBFs can be rewritten as integrals over q i , as we have done using equation (55).Then, as in equations ( 58) and ( 59), the Gaussian integral with respect to k can be evaluated.Finally, the remaining q i integrals will require many changes of variables and binomial expansions, following the work in this section to give a result in terms of incomplete Gamma and hypergeometric functions.

Our Methods Compared to Direct Integration
Our analytic results for triple-sBF integrals with exponential (base cases of §2.1) and Gaussian damping (base cases of §3 as well as equation 89) can be computed more efficiently than a direct integration of damped triple-sBF integrals.For a direct integration, one would need a vector of N k k points to perform the integral at each triple {r 1 , r 2 , r 3 } on a three-dimensional grid of the latter.If there are N r sample points for each r i , the total cost of a direct integration is N k N 3 r .We now turn to the computational complexity of our methods.
The Legendre functions of the second kind in the analytic result for the triple-sBF integral with exponential damping (equations 9, 11, and 12) all depend on R ±± ≡ −ip 2 ± r 2 ± r 3 /r 1 (equation 10).Given values for r 1 , r 2 , and r 3 (the free arguments of the sBFs), a new variable R can be defined such that it ranges from the minimum of R ±± to the maximum of R ±± .We may then map all required R ±± to R. Once this is done, the Legendre functions of the second kind may be evaluated on a one-dimensional vector of R. The numerator of R ±± causes this R grid to be be 4N r in length.However, the r 1 in the denominator will shrink this grid, since in a typical covariance calculation, the minimum r 1 is larger than 1 Mpc and so dividing by it will reduce the range covered by the numerator.The length of the R grid stems from our method's results rotating r 1 , r 2 , and r 3 to 45 • lines coming out of the origin of that coordinate system.The integrals only depend on the values along these lines.This concept was previously used in a numerical algorithm for sBF integrals ( [7], Figure 1).The computational cost of our method thus scales as N r rather than N k N 3 r .The incomplete Gamma and hypergeometric functions in the results for the triple-sBF integral with Gaussian damping (equations 35, 36, 37, and 89) all depend on r ±±± ≡ ±r 1 ± r 2 ± r 3 .Similarly to R, a new variable r ranging from the minimum of r ±±± to the maximum of r ±±± can be defined.The incomplete gamma and hypergeometric functions can be evaluated on a one-dimensional r vector with length 6N r .As the computational cost of our method scales as N r , it will be orders of magnitude faster than a direct numerical integration, even when accounting for the prefactors that arise in our method due to summing over a number of incomplete gamma and hypergeometric functions.
Additionally, the triple-sBF integrals do not need to be calculated explicitly for every {ℓ 1 , ℓ 2 , ℓ 3 } and {r 1 , r 2 , r 3 } combination.The recursion allows triple-sBF integrals to be expressed in terms of three base cases.If any of these base cases have already been evaluated and the result stored, less computations are required to obtain the desired integral.If, for example, all three required base cases have previously been evaluated, the recursion method is simply a combination of the base cases rather than a method with computational complexity N r .Thus, the computational cost of the recursion method will scale at most as N r , but will often be less costly.

Comparison of Results for the Gaussian-Damped triple-sBF Integral
We have derived results for the Gaussian-damped triple sBF integral in two different ways: with a recursion relation ( §3) and by the Hankel-Bowman method ( §6).The recursion relation can be used to evaluate Gaussian-damped triple-sBF integrals with the Jacobian k 2 dk ( §3).By parametric differentiation ( §5), this recursion method can be extended to integrals involving k to any even power ≥ 2. The Hankel-Bowman method ( §6), however, can be used to evaluate Gaussian-damped triple-sBF integrals with the Jacobian k n dk, as long as n is a non-negative integer.Thus, a major advantage of the Hankel-Bowman method is that it can be used to evaluate a wider range of integrals.
The recursion method requires at least three base cases (equations 35, 36, and 37), each of which depend on eight incomplete gamma functions and eight regularized hypergeometric functions.As demonstrated in §4, using the recursion to evaluate integrals with higher-order sBFs will require more than three base cases; therefore, each recursion will require computation of at least 24 incomplete gamma functions and at least 24 regularized hypergeometric functions.
The result for the Hankel-Bowman method (equation 89) appears to depend upon 32 incomplete gamma functions and eight hypergeometric functions.However, the parameters of these functions depend on α (ranging from 0 to n + ℓ 1 + ℓ 2 + ℓ 3 ), as well as the values of the variables being summed over in equation (89).Thus, in most cases, more than 32 incomplete gamma functions and more than eight hypergeometric functions will need to be computed for the Hankel-Bowman method.The recursion method will therefore typically be more efficient than the Hankel-Bowman method for evaluating Gaussian-damped triple-sBF integrals.

Conclusion
We have shown that the spherical Bessel function recursion relation can be used to obtain a recursion for triple-sBF integrals with exponential or Gaussian damping, generalizing an approach pioneered by [51].This approach requires three base cases to anchor it; once those base cases are determined, any triple-sBF integral with damping can be evaluated as long as the (integer) orders of the sBFs are at least zero, the argument of one of the sBFs is positive, and the arguments of the other two sBFs are nonzero.The orders of the sBFs and their arguments must also be real.The recursion allows the triple-sBF integral with exponential damping to be written in a simple form dependent upon Legendre functions of the second kind.For three sBFs with Gaussian damping, the recursion shows that the integral can be evaluated in terms of regularized hypergeometric and incomplete Gamma functions.
For the case with exponential damping, the base cases depend on Legendre functions with complex arguments.Each base case requires four Legendre functions (two negative complex conjugate pairs).By expanding the Legendre functions in an infinite series, we demonstrated that the base cases are indeed real, as expected given the reality of the integrand and of the domain of integration.
The work shown here is new; previously, a recursion relation had only been used to evaluate integrals with no damping term [51]; this is the p = 0 limit of our work (which is the same for the exponentially-damped and Gaussian cases).Notably, the previous work also required that the free arguments of the three sBFs, r i , be able to form a closed triangle, and also that the sum of the three sBF orders, ℓ i , be even.The results of the present work require only that {p, ℓ i , r i } ∈ R, ℓ i ∈ Z, ℓ i ≥ 0, r 1 > 0, and r i = 0, but demand no constraints between the individual r i and ℓ i .
The stability of the recursion depends on two factors: whether the base cases can be evaluated accurately and whether as many base cases as needed can be co-added without loss of numerical precision.Accurate numerical algorithms can be used to evaluate the special functions within the base cases.Numerical precision can be handled with existing libraries or balanced sum algorithms, which are useful for adding values with order of magnitude differences.The recursions are thus stable for both the exponentially-and Gaussian-damped integrals.
Additionally, the recursion relations for triple-sBF integrals with the Jacobian k 2 dk and exponential or Gaussian damping can be used to obtain results for such integrals weighted by higher powers of k.The partial derivative with respect to −p 2 of the damped triple-sBF integral can be moved inside the integral so that it acts on the damping.Each partial derivative increases the power of k by one for exponential damping and by two for Gaussian damping.Thus, by first evaluating the damped triple-sBF integral using the recursion relations outlined in §2.1 and §3 and then taking the necessary amount of partial derivatives with respect to −p 2 , the damped triple-sBF integrals can be evaluated for k to any integer power ≥ 2 for exponential damping and k to any even power ≥ 2 for Gaussian damping.
Finally, we have found a new method (the Hankel-Bowman method) to evaluate triple-sBF integrals with Gaussian damping weighted by any non-negative integer power of k.This method requires that {p, ℓ i , r i } ∈ R, {ℓ i , n} ∈ Z, ℓ i ≥ 0, r i = 0, p = 0, and n ≥ 0. The sBFs are rewritten as integrals derived from Hankel's contour integral, then many changes of variable and binomial expansions are used to obtain a result in terms of incomplete Gamma functions and hypergeometric functions.This method can also be used to evaluate integrals of any number of sBFs with Gaussian damping.While the recursion relation for the triple-sBF integral with Gaussian damping only can be used for the Jacobian k 2 dk, and extended by parametric differentiation to even powers of k that are greater than or equal to two, the work outlined in §6 allows such integrals weighted by any non-negative integer power of k to be evaluated.
We have demonstrated that our methods for Gaussian-damped triple-sBF integrals are more computationally efficient than a direct integration: our methods scale linearly as N r , while a direct integration would scale as N k N 3 r .The recursion method for Gaussian-damped triple-sBF integrals ( §3) requires computation of a minimum of 24 incomplete Gamma functions and 24 hypergeometric functions.The Hankel-Bowman method ( §6) requires at least 32 incomplete Gamma functions and eight hypergeometric functions to be evaluated.However, due to the many finite sums in the final result of the Hankel-Bowman method (equation 89), typically many more incomplete Gamma and hypergeometric functions will need to be computed.Therefore, it is preferable to use the recursion method over the Hankel-Bowman method for Gaussian-damped triple-sBF integrals with the Jacobian k n dk, where n is an even power ≥ 2.
the integrals in equation (85), we obtain a final result for I hg :10