Introduction

Our goal is efficient and robust numerical evaluation of Fourier integrals of the form $f(\vec{r}) = N_n(a,b) \int\, d^n k\, e^{+b i (\vec{k}\cdot\vec{r})}\, \tilde{f}(\vec{k}) \quad , \quad \tilde{f}(\vec{k}) = \tilde{N}_n(a,b) \int\, d^n r\, e^{-b i (\vec{k}\cdot\vec{r})}\,f(\vec{r}) \; , \label{eqn:fourier}$ with normalization factors $N_n(a,b) = |b|^{n/2} (2\pi)^{-n(1+a)/2} \quad , \quad \tilde{N}_n(a,b) = |b|^{n/2} (2\pi)^{-n(1-a)/2} \; ,$ where the constants $$a$$ and $$b$$ establish our choice of Fourier convention1. We focus on two- and three-dimensional ($$n=2,3$$) transforms of functions that can be adequately represented with a small number of (not necessarily low order) multipoles. Specifically, for $$n=2$$, we expand $f(r,\varphi_r) = \sum_{m=-\infty}^{+\infty} f_m(r)\,\Phi_m(\varphi_r) \quad , \quad \tilde{f}(k,\varphi_k) = \sum_{m=-\infty}^{+\infty} \tilde{f_m}(k)\,\Phi_m(\varphi_k) \; , \label{eqn:multipole2}$ using the polar basis functions $\Phi_m(\varphi) \equiv \frac{1}{\sqrt{2\pi}}\, e^{i m \varphi}$ with orthonormality2 ($$\delta_D$$ and $$\delta$$ are the Dirac and Kronecker delta functions, respectively): $\sum_{m=-\infty}^{+\infty} \Phi_m(\varphi)\Phi^\ast_m(\varphi') = \delta_D(\varphi-\varphi') \label{eqn:dirac2}$ and $\int_0^{2\pi} d\varphi\, \Phi_m(\varphi) \Phi_{m'}^\ast(\varphi) = \delta_{m m'} \; . \label{eqn:kronecker2}$ Similarly, for $$n = 3$$, we expand $f(r,\theta_r,\varphi_r) = \sum_{\ell=0}^{\infty}\sum_{m=-\ell}^{+\ell} f_{\ell m}(r)\, Y_{\ell m}(\theta_r,\varphi_r) \quad, \quad \tilde{f}(k,\theta_k,\varphi_k) = \sum_{\ell=0}^{\infty}\sum_{m=-\ell}^{+\ell} \tilde{f_{\ell m}}(k)\, Y_{\ell m}(\theta_k,\varphi_k) \; , \label{eqn:multipole3}$ using the spherical-harmonic basis functions3 (with associated Legendre polynomials $$P_{\ell}^m$$) $Y_{\ell m}(\theta,\varphi) \equiv \sqrt{\frac{2\ell+1}{2}\frac{(\ell-m)!}{(\ell+m)!}}\, P_{\ell}^m(\cos\theta) \Phi_m(\varphi)$ with orthonormality4 $\sum_{\ell=0}^{\infty}\sum_{m=-\ell}^{+\ell} Y_{\ell m}(\theta,\varphi)Y_{\ell m}^\ast(\theta',\varphi') = \delta_D(\cos\theta-\cos\theta') \delta_D(\varphi-\varphi') \label{eqn:dirac3}$ and5 $\int d\Omega\, Y_{\ell m}(\theta,\varphi) Y_{\ell' m'}^\ast(\theta,\varphi) = \delta_{\ell\ell'}\delta_{m m'} \; . \label{eqn:kronecker3}$ In the special case of a three-dimensional $$f(\vec{r})$$ that is cylindrically symmetric, i.e. has no $$\phi_r$$ dependence, only $$m = 0$$ terms contribute to the multipole expansion equation \ref{eqn:multipole3}. Since $Y_{\ell 0}(\theta,\varphi) = \sqrt{\frac{2\ell+1}{4\pi}}\,L_{\ell}(\mu) \label{eqn:sphericalmu}$ with $$L_{\ell}$$ the Legendre polynomial and $$\mu \equiv \cos\theta$$, it is then convenient to replace equation \ref{eqn:multipole3} with the equivalent expansion $f(r,\mu_r) = \sum_{\ell=0}^\infty f^{(\mu)}_{\ell}(r) L_{\ell}(\mu_r) \quad , \quad \tilde{f}(k,\mu_k) = \sum_{\ell=0}^\infty \tilde{f}^{(\mu)}_{\ell}(k) L_{\ell}(\mu_k) \; , \label{eqn:multipole3mu}$ in which the coefficient functions are simply rescaled $f^{(\mu)}_{\ell}(r) = \sqrt{\frac{2\ell+1}{4\pi}}\,f_{\ell 0}(r) \quad , \quad \tilde{f}^{(\mu)}_{\ell}(k) = \sqrt{\frac{2\ell+1}{4\pi}}\,\tilde{f_{\ell 0}}(k) \; . \label{eqn:coefsmu}$

1. Use $$a = 1$$ and $$b = 1$$ for the convention of references (Dodelson 2003, Weinberg 2013).

2. http://dlmf.nist.gov/1.17#E12

3. http://dlmf.nist.gov/14.30#E1

4. http://dlmf.nist.gov/1.17#E25

5. http://dlmf.nist.gov/14.30#E8

Our general approach is to translate the multi-dimensional Fourier transform, equation \ref{eqn:fourier}, into a small number (one per multipole) of one-dimensional Hankel ($$n=2$$) or spherical Bessel ($$n=3$$) transforms between the weight functions $$f_m(r)$$ and $$\tilde{f_m}(k)$$ or $$f_{\ell m}(r)$$ and $$\tilde{f_{\ell m}}(k)$$. We derive these transforms below then, in the following sections, describe our numerical algorithm for calculating them. Note that we rely on the fact that the Fourier transform does not mix the multipoles, so that there is a one-to-one correspondence between modes in r-space and k-space.

Using equation \ref{eqn:kronecker2}, we find for $$n=2$$ $f_m(r) = \int_0^{2\pi} d\varphi_r\, f(r,\varphi_r)\,\Phi_m^\ast(\varphi_r) \quad , \quad \tilde{f_m}(k) = \int_0^{2\pi} d\varphi_k\, \tilde{f}(k,\varphi_k)\,\Phi_m^\ast(\varphi_k)$ and similarly for $$n=3$$ using equation \ref{eqn:kronecker3} $f_{\ell m}(r) = \int d\Omega_r\, f(r,\theta_r,\varphi_r)\, Y_{\ell m}^\ast(\theta_r,\varphi_r) \quad , \quad \tilde{f_{\ell m}}(k) = \int d\Omega_k\, \tilde{f}(k,\theta_k,\varphi_k)\, Y_{\ell m}^\ast(\theta_k,\varphi_k) \; .$ In the case of $$n = 3$$ cylindrical symmetric, equations \ref{eqn:sphericalmu} and \ref{eqn:coefsmu} yield the equivalent $f^{(\mu)}_{\ell}(r) = \frac{2\ell+1}{2} \int_{-1}^{+1} d\mu_r\, f(r,\mu_r)\, L_{\ell}(\mu_r) \quad , \quad \tilde{f}^{(\mu)}_{\ell}(k) = \frac{2\ell+1}{2} \int_{-1}^{+1} d\mu_k\, \tilde{f}(k,\mu_k)\, L_{\ell}(\mu_k) \; .$ Next, we replace $$f(r,\varphi_r)$$ or $$f(r,\theta_r,\varphi_r)$$ in these equations with their Fourier expansions of equation \ref{eqn:fourier}, obtaining $f_m(r) = N_n(a,b) \int_0^\infty k dk\, \int_0^{2\pi} d\varphi_r\,\int_0^{2\pi} d\varphi_k\, e^{+b i (\vec{k}\cdot\vec{r})}\, \tilde{f}(k,\varphi_k)\,\Phi_m^\ast(\varphi_r)$ and $f_{\ell m}(r) = N_n(a,b) \int_0^\infty k^2 dk\, \int d\Omega_r\, \int d\Omega_k\, e^{+b i (\vec{k}\cdot\vec{r})}\, \tilde{f}(k,\theta_k,\varphi_k)\,Y_{\ell m}^\ast(\theta_r,\varphi_r) \; .$ Substituting the multipole expansions of $$\tilde{f}$$ from equations \ref{eqn:multipole2} and \ref{eqn:multipole3}, we find $f_m(r) = N_n(a,b) \int_0^\infty k dk \int_0^{2\pi} d\varphi_r \int_0^{2\pi} d\varphi_k\, e^{+b i (\vec{k}\cdot\vec{r})}\, \sum_{m'=-\infty}^{+\infty} \tilde{f_{m'}}(k)\,\Phi_{m'}(\varphi_k)\, \Phi_m^\ast(\varphi_r)$ and $f_{\ell m}(r) = N_n(a,b) \int_0^\infty k^2 dk \int d\Omega_r \int d\Omega_k\, e^{+b i (\vec{k}\cdot\vec{r})}\, \sum_{\ell'=0}^\infty \sum_{m'=-\ell'}^{+\ell'} \tilde{f_{\ell'm'}}(k)\,Y_{\ell'm'}(\theta_k,\varphi_k)\, Y_{\ell m}^\ast(\theta_r,\varphi_r) \; .$ Finally we expand the plane wave factor in terms of polar ($$n=2$$) basis functions $e^{+b i (\vec{k}\cdot\vec{r})} = 2\pi \sum_{m=-\infty}^{+\infty}(b i)^{m}\, J_{m}(kr)\,\Phi_{m}^\ast(\varphi_r)\,\Phi_{m}(\phi_k)$ or spherical harmonic ($$n=3$$) basis functions $e^{+b i (\vec{k}\cdot\vec{r})} = 4\pi \sum_{\ell=0}^{\infty}\sum_{m=-\ell}^{+\ell}\, (b i)^{\ell} j_{\ell}(kr)\, Y_{\ell m}^\ast(\theta_r,\varphi_r)\, Y_{\ell m}(\theta_k,\varphi_k) \; ,$ where $$J_m$$ and $$j_\ell$$ are the Bessel and spherical Bessel functions, respectively. Using equations \ref{eqn:dirac2} or \ref{eqn:dirac3} to simplify the resulting expressions, we find that the moments of $$f$$ and $$\tilde{f}$$ are related by a scaled Hankel transform when $$n = 2$$, $f_{m}(r) = c \int_0^\infty k dk\, J_m(kr)\, \tilde{f_m}(k) \quad , \quad c = 2\pi N_n(a,b) (b i)^m \; , \label{eqn:hankelr}$ and by a scaled spherical Bessel transform when $$n = 3$$, $f_{\ell m}(r) = c \int_0^\infty k^2 dk\, j_{\ell}(kr)\, \tilde{f_{\ell m}}(k) \quad , \quad c = 4\pi N_n(a,b) (b i)^{\ell} \; . \label{eqn:sbtr}$ Substituting $$f \leftrightarrow \tilde{f}$$, $$k \leftrightarrow r$$, $$N_n \rightarrow \tilde{N_n}$$, and $$b \rightarrow -b$$, we find similar results for the inverse transforms $\tilde{f_m}(k) = c \int_0^\infty r dr\, J_m(kr)\, f_m(r) \quad , \quad c = 2\pi \tilde{N_n}(a,b) (-b i)^m \; , \label{eqn:hankelk}$ and $\tilde{f_{\ell m}}(k) = c \int_0^\infty r^2 dr\, j_{\ell}(kr)\, f_{\ell m}(r) \quad , \quad c = 4\pi \tilde{N_n}(a,b) (-b i)^{\ell} \; . \label{eqn:sbtk}$ The $$n = 3$$ cylindrically symmetric $$f'_{\ell}(r)$$ and $$\tilde{f}'_{\ell}(k)$$ transform identically to $$f_{\ell m}(r)$$ and $$\tilde{f_{\ell m}}(k)$$, as given in equations \ref{eqn:sbtr} and \ref{eqn:sbtk}, since they are simply related by an $$\ell$$-dependent rescaling factor.