Multipole Transforms


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).





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.

Spherical Bessel Transforms

Our algorithm for evaluating the scaled spherical Bessel transform equations \ref{eqn:sbtr} and \ref{eqn:sbtk}1 follows earlier work in references (Talman 2009, Hamilton 2000).

With the change to dimensionless variables \[\kappa = \log (k/k_0) \quad , \quad \rho = \log (r/r_0) \; ,\] we can write equation \ref{eqn:sbtr} as a convolution suitable for FFT evaluation \[(r/r_0)^\alpha\cdot f_{\ell m}(r) = \int_{-\infty}^{+\infty} G(\rho + \kappa) F(\kappa) d\kappa \label{eqn:conv}\] with dimensionless functions \[G(s) = e^{\alpha s}\, j_{\ell}(k_0 r_0 e^s) \quad , \quad F(s) = c\, e^{(3-\alpha)s}\,k_0^3\,\tilde{f_{\ell m}}(k_0 e^s) \; ,\] where the \(r\)-weighting parameter \(\alpha\) is arbitrary at this point. The function \(F\) depends on \(f\), but \(G(s)\) depends only on the values of \(k_0 r_0\) and \(\alpha\) (see Fig. \ref{fig:fplot}), and has asymptotes \[G_{+}(s) = \frac{1}{k_0 r_0}\, e^{(\alpha-1)s} \quad , \quad G_{-}(s) = \frac{2^{-(\ell+1)}\sqrt{\pi}}{\Gamma(\ell+3/2)}\,(k_0 r_0)^\ell\,e^{(\alpha+\ell)s} \; ,\] where \(G_+\) is the envelope of oscillations with a rapidly decreasing period \(2\pi e^{-s}/(k_0 r_0)\).

In order to simplify the calculations, we use our freedom to chose \(k_0 r_0\) and \(\alpha\) to symmetrize \(G(s)\) via \[\alpha = \frac{1-\ell}{2} \quad , \quad \left( \frac{k_0 r_0}{2}\right)^{\ell+1} = \frac{\Gamma(\ell+3/2)}{\sqrt{\pi}} \;, \label{eqn:symm}\] which yields \[G_+(s) = G_-(-s) = \frac{1}{2}\left[ \frac{\sqrt{\pi}}{\Gamma(\ell+3/2)}\right]^{1/(\ell+1)}\cdot \exp\left(-(\ell+1)s/2\right) \; .\] We can write this using \[s_0 = \frac{2}{\ell + 1} \quad , \quad G_0 = \frac{1}{k_0 r_0} \;,\] as \[G_+(s) = G_-(-s) = G_0 e^{-s/s_0} \;.\]

  1. Since these are simply related by the substitution \(r \leftrightarrow k\), we use the notation of equation \ref{eqn:sbtr} in the following, without any loss of generality.