Rotation method for accelerating multiple-spherical Bessel function integrals against a numerical source function

Understand

A common problem in cosmology is to integrate the product of two or more spherical Bessel functions (sBFs) with different configuration-space arguments against the power spectrum or its square, weighted by powers of wavenumber.

• Naively computing them scales as $N_{\rm g}^{p+1}$ with $p$ the number of configuration space arguments and $N_{\rm g}$ the grid size, and they cannot be done with Fast Fourier Transforms (FFTs).
• Here we show that by rewriting the sBFs as sums of products of sine and cosine and then using the product to sum identities, these integrals can then be performed using 1-D FFTs with $N_{\rm g} \log N_{\rm g}$ scaling.
• This "rotation" method has the potential to accelerate significantly a number of calculations in cosmology, such as perturbation theory predictions of loop integrals, higher order correlation functions, and analytic templates for correlation function covariance matrices.