On Decoupling the Integrals of Cosmological Perturbation Theory

TL;DR

This paper introduces a method to decouple complex integrals in cosmological perturbation theory, reducing high-dimensional integrals to sums of independent 1-D integrals, significantly improving computational efficiency for large-scale structure analysis.

Contribution

The authors adapt a technique from Quantum Field Theory to efficiently compute higher-order corrections in cosmological perturbation theory by decoupling integrals into products of 1-D integrals.

Findings

Reduces 9-D integrals to sums over 1-D integrals

Enables faster computation with $N_g \, \log N_g$ scaling

Facilitates large-scale structure modeling for upcoming surveys

Abstract

Perturbation theory (PT) is often used to model statistical observables capturing the translation and rotation-invariant information in cosmological density fields. PT produces higher-order corrections by integration over linear statistics of the density fields weighted by kernels resulting from recursive solution of the fluid equations. These integrals quickly become high-dimensional and naively require increasing computational resources the higher the order of the corrections. Here we show how to decouple the integrands that often produce this issue, enabling PT corrections to be computed as a sum of products of independent 1-D integrals. Our approach is related to a commonly used method for calculating multi-loop Feynman integrals in Quantum Field Theory, the Gegenbauer Polynomial $x$-Space Technique (GPxT). We explicitly reduce the three terms entering the 2-loop power spectrum, formally requiring 9-D integrations, to sums over successive 1-D radial integrals. These 1-D integrals can further be performed as convolutions, rendering the scaling of this method $N_{\rm g} \log N_{\rm g}$ with $N_{\rm g}$ the number of grid points used for each Fast Fourier Transform. This method should be highly enabling for upcoming large-scale structure redshift surveys where model predictions at an enormous number of cosmological parameter combinations will be required by Monte Carlo Markov Chain searches for the best-fit values.

Equations

$${\mathcal{I}}_{ij,\mathrm{SV}}(k,\vec{\alpha },\vec{\beta })=$$

$$\int \frac{d{\Omega }_k}{4\pi }\int \frac{d^3{\vec{q}}_1}{(2\pi {)}^3}\frac{d^3{\vec{q}}_2}{(2\pi {)}^3}\frac{e^{i\vec{\alpha }\cdot {\vec{q}}_1}{e}^{i\vec{\beta }\cdot {\vec{q}}_2}}{q_1^{2n_1}|\vec{k}+{\vec{q}}_1{|}^{2n_1^{\prime }}{q}_2^{2n_2}|\vec{k}+{\vec{q}}_2{|}^{2n_2^{\prime }}}$$

$$\times \frac{P_{\mathrm{lin}}(q_1)P_{\mathrm{lin}}(q_2)}{|{\vec{q}}_1+{\vec{q}}_2{|}^{2n_3}|\vec{k}+{\vec{q}}_1+{\vec{q}}_2{|}^{2n_3^{\prime }}}{P}_{\mathrm{lin}}(|{\vec{w}}_{ij}|)$$

$${\vec{w}}_{15}=\vec{k},{\vec{w}}_{24}=\vec{k}+{\vec{q}}_2,{\vec{w}}_{33}=\vec{k}+{\vec{q}}_1+{\vec{q}}_2.$$

$$N({\vec{q}}_1,{\vec{q}}_2)=\sum _{\ell =0}^{{\ell }_{\max }}{N}_{\ell }(q_1,q_2){\mathcal{L}}_{\ell }({\hat{q}}_1\cdot {\hat{q}}_2)$$

$$=\sum _{\ell =0}^{{\ell }_{\max }}{N}_{\ell }^{[1]}(q_1)N_{\ell }^{[2]}(q_2){\mathcal{L}}_{\ell }({\hat{q}}_1\cdot {\hat{q}}_2)$$

$$=\sum _{\ell =0}^{{\ell }_{\max }}{N}_{\ell }^{[1]}(q_1)N_{\ell }^{[2]}(q_2)\frac{4\pi }{(2\ell +1)}\sum _{m=-\ell }^{\ell }{Y}_{\ell m}({\hat{q}}_1)Y_{\ell m}^{\ast }({\hat{q}}_1).$$

$${\mathcal{I}}_{ij}(k)=\sum _{\ell =0}^{{\ell }_{\max }}\sum _{m=-\ell }^{\ell }\frac{4\pi }{2\ell +1}\int \frac{d{\Omega }_k}{4\pi }\int \frac{d^3{\vec{q}}_1}{(2\pi {)}^3}\frac{d^3{\vec{q}}_2}{(2\pi {)}^3}$$

$$\times \frac{N_{\ell }^{[1]}(q_1)N_{\ell }^{[2]}(q_2)Y_{\ell m}({\hat{q}}_1)Y_{\ell m}^{\ast }({\hat{q}}_1)}{q_1^{2n_1}|\vec{k}+{\vec{q}}_1{|}^{2n_1^{\prime }}{q}_2^{2n_2}|\vec{k}+{\vec{q}}_2{|}^{2n_2^{\prime }}}$$

$$\times \frac{P_{\mathrm{lin}}(q_1)P_{\mathrm{lin}}(q_2)}{|{\vec{q}}_1+{\vec{q}}_2{|}^{2n_3}|\vec{k}+{\vec{q}}_1+{\vec{q}}_2{|}^{2n_3^{\prime }}}{P}_{\mathrm{lin}}(|{\vec{w}}_{ij}|)$$

$${\mathcal{I}}_{ij}(k)=\sum _{\ell =0}^{{\ell }_{\max }}\sum _{m=-\ell }^{\ell }\frac{4\pi }{2\ell +1}{\mathcal{I}}_{ij,\ell m}(k).$$

$${\mathcal{I}}_{ij,\ell m}(k)=$$

$$\int \frac{d{\Omega }_k}{4\pi }\int \frac{d^3{\vec{q}}_1}{(2\pi {)}^3}\frac{d^3{\vec{q}}_2}{(2\pi {)}^3}K(\vec{k},{\vec{q}}_1,{\vec{q}}_2)P_{\mathrm{lin}}(|{\vec{w}}_{ij}|).$$

$$K=K_{ij}^{\mathrm{c}}{K}_{ij}^{\mathrm{d}}$$

$$K_{15}^{\mathrm{c}}$$

$$={\mathcal{R}}_{\ell }^{n_1}(q_1){\mathcal{P}}^{n_1^{\prime }}(\vec{k}+{\vec{q}}_1){\mathcal{R}}_{\ell }^{n_2}(q_2){\mathcal{P}}^{n_1^{\prime }}((\vec{k}+{\vec{q}}_1)+{\vec{q}}_2)$$

$${\mathcal{R}}_{\ell }^n(q_t)=N_{\ell }^{[t]}(q_t)q_t^{-2n}{P}_{\mathrm{lin}}(q_t),$$

$${\mathcal{P}}^n(\vec{q})=|\vec{q}{|}^{-2n},$$

$$K_{15}^{\mathrm{d}}=\frac{1}{|{\vec{q}}_1+{\vec{q}}_2{|}^{2n_3}|\vec{k}+{\vec{q}}_2{|}^{2n_2^{\prime }}}.$$

$$K_{24}^{\mathrm{c}}={\mathcal{R}}_{\ell }^{n_1}(q_1){\mathcal{P}}^{n_1^{\prime }}(\vec{k}+{\vec{q}}_1){\mathcal{R}}_{\ell }^{n_2}(q_2){\mathcal{P}}^{n_1^{\prime }}(\vec{k}+{\vec{q}}_2)$$

$$K_{24}^{\mathrm{d}}=\frac{1}{|{\vec{q}}_1+{\vec{q}}_2{|}^{2n_3}|{\vec{q}}_1+(\vec{k}+{\vec{q}}_2){|}^{2n_3^{\prime }}}.$$

$$K_{33}^{\mathrm{c}}=K_{15}^{\mathrm{c}},K_{33}^{\mathrm{d}}=K_{15}^{\mathrm{d}}.$$

$$\frac{1}{|{\vec{p}}_1+{\vec{p}}_2{|}^2},\frac{1}{|{\vec{p}}_1+{\vec{p}}_2{|}^4},$$

$$\frac{1}{|{\vec{p}}_1+{\vec{p}}_2|}=\sum _{L=0}^{\infty }\frac{p_{<}^L}{p_{>}^{L+1}}{\mathcal{L}}_L({\hat{p}}_1\cdot {\hat{p}}_2).$$

$$\frac{2}{\pi }(2L+1)\int dx{j}_L(p_{1}x)j_L(p_{2}x)$$

$$=\frac{p_1^L}{p_2^{L+1}},p_1<p_2,\frac{p_2^L}{p_1^{L+1}},p_2<p_2.$$

$$\int dx{j}_L(p_{1}x)j_L(p_{2}x)=\sum _{n=0}^{\infty }n{ϵ}_n{j}_L(n{p}_1)j_L(n{p}_1)$$

$${ϵ}_n=1,n>0;=1/2,n=0.$$

$$\frac{1}{|{\vec{p}}_1+{\vec{p}}_2|}=\sum _{L=0}^{\infty }(2L+1)\sum _{n=0}^{\infty }{\varphi }_{nLM}({\vec{p}}_1){\varphi }_{nLM}^{\ast }({\vec{p}}_2),$$

$${\varphi }_{nLM}(\vec{p})=\sqrt{n{ϵ}_n}{j}_L(np)Y_{LM}(\hat{p}),$$

$$\frac{1}{|{\vec{p}}_1+{\vec{p}}_2{|}^{2\lambda }}=\sum _{L=0}^{\infty }\frac{p_{<}^L}{p_{>}^{L+2\lambda }}{C}_L^{(\lambda )}({\hat{p}}_1\cdot {\hat{p}}_2)$$

$$C_L^{(\lambda )}(x)=\sum _{J=0}^L{w}_J^{L,\lambda }\sum _{S=-J}^J{Y}_{JS}(\hat{a})Y_{JS}^{\ast }(\hat{b}),$$

$$\frac{2}{\pi }{\int }_0^{\infty }dxx\left[j_{L+1}(x{p}_1)j_L(x{p}_2)+j_{L+1}(x{p}_2)j_L(x{p}_1)\right]$$

$$=\frac{p_2^L}{p_1^{L+2}},p_1>p_2,\frac{p_1^L}{p_2^{L+2}},p_2>p_1.$$

$$\frac{2}{\pi }{\int }_0^{\infty }dxx\left[j_{L+1}(x{p}_1)j_L(x{p}_2)+j_{L+1}(x{p}_2)j_L(x{p}_1)\right]$$

$$=\frac{4}{{\pi }^2}\sum _{n=0}^{\infty }n{ϵ}_n\sqrt{p_1{p}_2}\left[j_{L+1}(n{p}_1)j_L(n{p}_2)+j_{L+1}(n{p}_2)j_L(n{p}_1)\right].$$

$$\displaystyle\frac{1}{|\vec{p}_{1}+\vec{p}_{2}|^{2}}=\frac{4}{\pi^{2}}\sum_{L=0}^{\infty}\sum_{n=0}^{\infty}n\epsilon_{n}\;\sqrt{p_{1}p_{2}}\bigg{[}j_{L+1}(np_{1})j_{L}(np_{2})$$

Full text

††pubyear: 2018††pagerange: On Decoupling the Integrals of Cosmological Perturbation Theory–F

1 Introduction

In the consensus picture of cosmological structure formation, perturbations that were Gaussian-distributed in amplitude and uniform-randomly-distributed in phase (when written in Fourier space) were generated at the end of inflation (Starobinsky 1982, Bardeen et al. 1983). Baryon Acoustic Oscillations (BAO) in the ionized plasma present prior to decoupling111When an electron was last scattered by a photon; slightly later than the last scattering of a photon by an electron, which occurred at recombination z$$\sim$$1100. (z$$\sim$$1020) (Sakharov 1966, Silk 1968, Sunyaev & Zeldovich 1970, Peebles & Yu 1970, Bond & Szalay 1983, Holtzman 1989, Hu & Sugiyama 1996, Eisenstein & Hu 1998, Weinberg 2002, Eisenstein et al. 2005, Cole et al. 2005, Eisenstein et al. 2007, Slepian & Eisenstein 2016b) then imprinted additional spatial correlations on the perturbations, which were then amplified by gravitational instability (e.g. Goroff et al. 1986, Jain & Bertschinger 1994) into the large-scale structure observed in low-redshift galaxy surveys (Strauss & Willick 1995, Coil 2013, for reviews). It is typically taken that the dark matter density field was dominantly shaped by gravitational evolution, that halos form in this matter field according to their local environment, characterized by bias models (e.g. Kaiser 1984, Bardeen et al. 1986, Fry & Gaztanaga 1993, Mo & White 1996; Desjacques et al. 2018 for a recent review) or Effective Field Theory (EFT) parameters to be measured from simulations (e.g. McDonald & Roy 2009, Carrasco et al. 2012, Senatore 2015), and that galaxies form in the dark matter halos according to both environment and detailed gas and astrophysical processes. For purposes of large-scale structure these latter effects are also often compressed into the bias coefficients or EFT parameters.

The evolution of the dark matter field under gravity is a numerically straightforward problem in the sense that the gravitational force law on the relevant scales is well-understood. Nonetheless it is computationally demanding to evolve the matter field in large cosmological volumes for direct comparison with observational proxies for the true density field, given that we do not know the initial 3-D distribution of matter in the Universe (though recent work has made some progress in this direction, e.g. Kitaura & Enßlin 2008, Jasche & Kitaura 2010, Jasche & Wandelt 2013, Ata et al. 2015, Leclercq et al. 2015, Seljak et al. 2017, Schmittfull et al. 2017, Jasche & Lavaux 2018, Schmidt et al. 2018). Thus, as an alternative, translation and rotation-invariant clustering statistics such as the 2-Point Correlation Function (2PCF), power spectrum, 3-Point Correlation Function, and bispectrum, which measure excess clustering over random of galaxy pairs or triplets in respectively configuration and Fourier space, are used (Bernardeau et al. 2002, for a review). These statistics can be predicted by averaging many sub-regions of numerical simulations (though Pontzen et al. 2016 and Angulo & Pontzen 2016 suggest this can be improved by phase-matching in the initial conditions), but can also be, less accurately but far more quickly, obtained from analytic solution of the approximate equations of motion for the dark matter. These equations of motion are approximate for a number of reasons: they neglect possible velocity dispersion in the CDM trajectories, late-time couplings to other fluids (e.g. neutrinos), and they do not extend beyond shell-crossing, where CDM trajectories cross each other.

The fluid equations are solved in Eulerian coordinates by recursion relations (Goroff et al. 1986, Jain & Bertschinger 1994) giving kernels that systematically generate higher-order density corrections as integrals of lower-order fields against the kernels.222Recursive solutions also exist for Lagrangian PT, although we do not focus on them here; see e.g. Zheligovsky & Frisch (2014) for Einstein-de Sitter, Rampf et al. (2015) and Matsubara (2015) for $\Lambda$CDM. The higher-order corrections to the density can then be used to compute corrections to the clustering statistics. In particular, the corrections to the clustering statistics end up as integrals of products of the kernels against linear power spectra. The advantage is that in the initial field, modulo a small amount of possible primordial non-Gaussianity (PNG), the linear power spectrum contains all of the information. Further, the linear power spectrum can be quickly obtained numerically from linear Boltzmann solver codes such as CMBFAST (Seljak & Zaldarriaga, 1996), CAMB (Lewis et al., 2000), or CLASS (Lesgourgues, 2011). However, since the linear theory power spectrum is not known in closed form (save for approximately), the relevant integrals must be done numerically, and redone every time a new set of cosmological parameters is desired as these produce a new (numerical) linear power spectrum.

These integrals quickly become high-dimensional, making it numerically cumbersome to compute the corrections. This has not been an entirely limiting factor in previous analyses as often the cosmology was not varied, varied in a way that could be adjusted after this computation (e.g. BAO analyses varying the ratio $\alpha$ of the model to fiducial sound horizon) or varied over only a few fiducial models requiring recomputation. However, for future analyses of upcoming datasets such as DESI (DESI Collaboration et al., 2016) and LSST (Laureijs et al., 2011), it will be desirable to vary the cosmological parameters an enormous number of times, much as is already done with Monte Carlo Markov Chain (MCMC) analyses of the Cosmic Microwave Background (CMB; e.g. Lewis & Bridle 2002, Planck Collaboration et al. 2018). Thus greater speed in computation of the corrections to clustering statistics would be useful.

Several recent approaches to this problem have been proposed. Appendix A of Slepian & Eisenstein (2015a) shows that using decoupling into separated 3-D integrals and then factorization into radial and angular pieces, certain fourth order terms in the 2PCF (1-loop corrections) can be computed as 1-D Hankel transforms by performing the angular integrations analytically.333A Hankel transform is just an integral of a given function against a spherical Bessel function with a particular power-law weight. Ferraro et al. (2012) also exploited this idea even earlier (personal communication) to obtain their equations 13-15 although the details of the derivation do not appear in the paper. Angular-radial factorizations also substantially simplify the tree-level 3PCF into products of 1-D transforms of the power spectrum (e.g. Slepian & Eisenstein 2015a equations 51-56), as might be expected since the tree-level 3PCF is also fourth order and, like the 1-loop 2PCF, just involves 1-loop kernels (some of these transforms become 2-D in redshift space; Slepian & Eisenstein 2017 §2.2).

Schmittfull et al. (2016) develops this idea much more fully and demonstrates all 1-loop corrections to the power spectrum can be computed in this way. Schmittfull & Vlah (2016) extends this analysis to 2-loop corrections, but finds that several of the 2-loop terms are not amenable to the methods of separation employed in the previous works above and thus introduces a dummy free variable to be integrated over at the end subject to a Dirac Delta function. This trick allows formal separation. Thus at 2-loop one requires 2-D integrals for a number of terms; indeed one has an infinite sum over such integrals (though it converges quickly in practice). These 2-D integrals dominate the computational cost of the full 2-loop calculation (Schmittfull & Vlah 2016 §V.D).

McEwen et al. (2016) and Fang et al. (2017) combine these ideas with expressing the power spectrum as a sum of complex power laws to further accelerate the 1-loop power spectrum computations; the second of the above extends this to non-isotropic quantities relevant for working in redshift space. In this latter case, one introduces velocity as well as density kernels, but they are the same in form with different coefficients so the same techniques apply. Another complication in redshift space is having a preferred direction (the line of sight) and projecting onto it and transverse to it. So doing simply introduces additional angular dependences to the integrals, but these are known in closed form and can be factorized and dealt with analytically.

Simonović et al. (2018) takes a related approach, showing that the integrals required for power laws can all be done analytically, meaning that the challenging part, that of the coupled-multi-dimensional integrals, need only be done once, and then assembled in the linear combination dictated by the power law weights of the initial decomposition of the linear power spectrum. Assassi et al. (2017) uses this approach for angular statistics that come from integrating the power spectrum against spherical Bessel functions (sBFs), as does Grasshorn Gebhardt & Jeong (2018). Both of these works express the double spherical Bessel function integrals involved in the angular statistics in terms of hypergeometric functions, and the latter particularly focuses on developing a stable, fast, and accurate recursive approach to their computation.

Fonseca de la Bella et al. (2017) uses factorization and decoupling to reduce the 1-loop power spectrum computation to 1-D integrals in a slightly different way than Schmittfull et al. (2016) and McEwen et al. (2016). Writing the Fourier space integral for the loop correction with a Dirac Delta function to enforce momentum conservation rather than as an explicit convolution, they then rewrite the Delta function as the inverse FT of unity. This trick is also used in Schmittfull & Vlah (2016) for the 2-loop integrals. The inverse FT’s plane waves can be expanded into sBFs and spherical harmonics (their §3.2), and the integrals over sBFs done analytically (their Appendix B).444Their Appendix is also a rather complete reference for previous work on analytic evaluation of multiple sBF integrals, although attention should also be drawn to a method for evaluating triple-sBF integrals by recursion in the Appendix of Wang & Kamionkowski (2000), which they do not mention. Adkins (2013) also has a number of interesting results on singular integrals of sBFs. This resolution of the Delta function was employed earlier in Slepian & Eisenstein (2015b) (their equation 58) to reduce the disconnected piece of the 3PCF covariance to tractable 2-D integrals. As already mentioned earlier, Schmittfull & Vlah (2016) also used a Delta function for the 2-loop calculation, although with a different expansion of it than Fonseca de la Bella et al. (2017) use for the 1-loop calculation. Rewriting the Dirac delta function is also used in Böhm et al. (2016) (Appendix B footnote 6). This latter work also employs angular-radial decoupling (equations B1 and B2) and products in real space to reduce 3-D convolutions to integrals of 1-D product integrals (equation B10). Factorization of the Delta function by rewriting as an inverse FT is also used to accelerate the bispectrum algorithm of Scoccimarro (2015); we will further discuss algorithms taking advantage of factorization slightly later. Schmidt et al. (2018) shows that all of the perturbation theory corrections and biasing terms at 1-loop can be written in terms of a set of 28 independent integrals; their Appendix H discusses the fast evaluation of these integrals exploiting the same approach of using configuration space to perform products.555Another, less-well-known but useful case where a 3-D computation in configuration space can be simplified with analytic angular integration is Zehavi et al. (2005) equation 5 for $\sigma_{R}$, the standard deviation of the overdensity smoothed to a given scale $R$. $\sigma_{R}$ can be written as a 1-D integral in configuration space, despite that it naively appears to require a convolution there. The trick is to convolve geometrically as the overlap lens between two spheres. Of course $\sigma_{R}$ is standardly expressed as a 1-D integral in Fourier space, but working in configuration space is useful if one wishes to compute the variance of data. One then avoids needing to grid the data for an FT. A similar geometric trick is used in Hand et al. (2017) §3.3. to examine a toy-model of boundary effects on the power spectrum multipoles.

Taruya et al. (2018) exploits the fact that, given a realization of the initial, Gaussian, linear density and velocity fields, computation of higher-order corrections to the fields themselves (as opposed to their statistics) is a series of 3-D forward and inverse FTs. This work thus achieves $N_{\rm g}\log N_{\rm g}$ evaluation of these corrections, where $N_{\rm g}$ is the number of grid points used for the FT. However the price of this method is cosmic variance: many realizations must be averaged over to obtain precise density field statistics. McDonald (personal communication) also shows that the 1-loop density and velocity field corrections and their statistics can be expressed with configuration space products rather than the Fourier-space convolutions the PT recursion relations imply. The method presented in that work is analogous to the split-step Fourier method in spectral solution of partial differential equations, where one takes diffusion steps in Fourier space but time-steps in configuration space (e.g. Agrawal 2001). McDonald suggests taking products in configuration space but applying derivative operators in Fourier space.

The purpose of the present work is evaluating integrals in cosmological perturbation theory; of course, this is only worthwhile if PT offers a reasonable description of the clustering of matter. Full discussion of the convergence of PT is beyond the scope of this work. In 3-D, the corrections rapidly become complicated, with a large number of terms. Zheligovsky & Frisch (2014) (Einstein-de Sitter) and Rampf et al. (2015) ($\Lambda$CDM) show the convergence in 3-D up to a certain time away from the initialization. In 1-D corrections of arbitrarily higher order are more tractable, rendering it easier to compare the fully-summed PT calculation to simulations. Several works have investigated convergence from this perspective: Novikov 1969, McQuinn & White 2016, Taruya & Colombi 2017, Pietroni 2018. We also point to Rampf & Frisch 2017 which works in quasi-1-D. McQuinn & White (2016) finds that it converges but not to the correct answer, although Pietroni (2018) suggest adding non-perturbative terms is able to greatly improve the agreement, a point also implicit in the discussion of Rampf & Frisch (2017). Pajer & van der Woude (2018) extend the work of McQuinn & White (2016) to non-Gaussian initial conditions and the bispectrum, as well as showing that divergences occur in the configuration-space 2PCF. Rampf (2017) analyses the quasi-spherical collapse of a perturbed tophat, finding that generically collapse occurs earlier and the linear density at collapse is reduced relative to the spherical case. Saga et al. 2018 investigates Lagrangian PT in a semi-3-D limit that might be described as $1+1+1$, i.e. halos are modeled as being seeded by three crossed sine waves along the three Cartesian axes. This model allows derivation of recursion relations that are simple enough to evaluate to high (tenth) order and can be compared with simulations; reasonable agreement is found.

McDonald & Vlah (2018) presents a new method in principle applicable in 3-D although only tested numerically in 1-D thus far. The method performs well in 1-D, going beyond the point where PT traditionally breaks down (shell-crossing, where the determinant of the Jacobian between Eulerian and Lagrangian space becomes singular). It is inspired by the Hubbard-Stratonovich transformation typically used in Quantum Field Theory (QFT), which is in McDonald & Vlah (2018) used to carry an exact exponential of the displacement field farther into the calculation, allowing increased accuracy. Motivated by this work we investigated use of a Hubbard-Stratonovich transformation to simplify the integrals of perturbation theory treated here, but without success. Overall, application of techniques from QFT to cosmological perturbation theory is a rich subject of future investigation.666Suggestively, Simonović et al. (2018) notes that the 2-loop standard PT integrals are equivalent to those in a massless scalar QFT with a cubic interaction. The method presented in the present work is in fact similar, up to a point, to a well-known method for computing multi-loop Feynman integrals called the Gegenbauer Polynomial $x$-Space Technique (GPXT; Chetyrkin et al. 1980). We return to this and other connections to QFT in §7.2.

The idea of factorization has also been more generally exploited to develop fast algorithms for measuring the clustering of large-scale structure. Slepian & Eisenstein (2015b, 2016a, 2018)777Implementations of these latter two are presented in Portillo et al. (2018) and Friesen et al. (2017)., and Sugiyama et al. (2018) all use splitting of Legendre polynomials into spherical harmonics by the addition theorem (100) to accelerate measurement of the 3PCF as well as the anisotropic 2PCF, and Hand et al. (2017) also exploits this factorization to develop a fast algorithm for measuring the anisotropic power spectrum to high multipole.888Publicly available through nbodykit; see Hand et al. 2018. Bianchi et al. (2015) also uses an angular-radial splitting to measure the anisotropic power spectrum, although their basis is not orthogonal and so ends up requiring more FTs than Hand et al. (2017).

In this work we take a somewhat different tack from the PT works discussed above although conceptually the motives are the same. Much as these earlier works largely do, we seek to exploit the isotropy of the power spectrum in combination with angular-radial factorization to perform as many integrals as possible analytically just once. We also require a general way of decomposing the types of integration kernel that generically lead to the coupled nature of the loop integrals, and thus their high dimensionality and numerical costliness.

All of the 2-loop contributions are over two linear momenta. If certain coupled terms could be decoupled, each of these two momentum integrals can be written as a convolution, one nested inside the other. This idea is similar to that in Taruya et al. (2018) except at the level of the density field statistics rather than the raw field. The nested double convolution gives one four different arguments in the integrand “for free,” two for each convolution, for instance $\vec{q}_{1},\vec{q}_{1}+\vec{a}$ and $\vec{q}_{2},\vec{q}_{2}+\vec{b}$ if one is integrating over dummy momenta $\vec{q}_{1}$ and $\vec{q}_{2}$. The convolutions can then be factored into radial and angular pieces and reduced to 1-D radial integrals against the linear power spectrum.

Our main problem is to decouple the coupled terms in the 2-loop integrands that, left unaddressed, spoil the double convolutionality. We need both to decouple them and then factorize them into angular and radial pieces to be able to incorporate them into the convolutional terms in the integrand. There are six arguments in the 2-loop integrals, and as noted above, convolutions deal with four “for free”; hence we will need to decouple at most two terms. We note that throughout this work, we use the word “decouple” to denote pulling apart two 3-D vectors so they may be integrated over separately, and the word “factorize” to denote separating radial and angular dependences within functions of a single 3-D vector. We also use the word decouple to indicate removing a constraint on the magnitude of one 3-D vector relative to that of another, i.e. that it be less than the other; our meaning should be clear from context.

This work is laid out as follows. In §2 we present the general integrals we treat and show that if certain parts of the integrand could be decoupled the rest of the calculation becomes a set of nested 3-D convolutions. In §3 we show how to decouple these parts of the integrands by deriving a decoupled and factorized eigenfunction expansion for them. In §4 we insert these expansions and show the whole calculation is now a set of 3-D convolutions. In §5 we evaluate these convolutions, ultimately reducing them to an infinite sum over 1-D integrals that can be done quickly using FFTs. §6 discusses possible extensions of our results, connections to QFT, and outlines how a numerical implementation might proceed. We conclude in §8.

2 Casting the 2-Loop Integrals as Double Convolutions

Throughout, our convention will be that forward 3-D FTs will have a positive $i$ in the plane wave, and inverse FTs a negative $i$ and be normalized by $(2\pi)^{-3}$. Schmittfull & Vlah (2016) express the three 2-loop integrals as

[TABLE]

with

[TABLE]

Relative to Schmittfull & Vlah (2016) we have added averaging over $d\Omega_{k}$, which they presumably intended since the power spectrum and its loop corrections are isotropic (ignoring for the moment redshift-space distortions (RSD)). The vectors $\vec{\alpha}$ and $\vec{\beta}$ are parameters with respect to which one can differentiate to generate more complicated numerators, and then take the limit $\vec{\alpha},\vec{\beta}\to 0$. We note that any dependence on $\vec{q}_{1}$ and $\vec{q}_{2}$ in the numerator must be rotation-invariant as the PT kernels cannot depend on absolute direction. We also pause to note that these integrals are not divergent: for instance, by momentum conservation, if $\vec{q}_{1}=-\vec{q}_{2}$, then $k\to 0$ and $P_{\rm lin}(0)=0$. An analogous point applies to the other denominators that could potentially vanish.

We notice that all three 2-loop integrals have largely the same structure, with just a different argument of the third linear power spectrum. We choose to drop the exponentials in favor of writing out a Legendre series for the terms in the numerator they are meant to capture. We have

[TABLE]

Noticing that any parametric differentiation of the exponentials in equation (1) would give rise to a factorizable function of $q_{1}$ and $q_{2}$ allows us to assert the second equality above, and in the third line we have factored the angular dependence using the spherical harmonic addition theorem (100). We do note that the numerator may need an additional sum at each $\ell$ to be factored radially; a simpler example of this type of behavior is in the dipole term of the $F^{(2)}$ kernel, which is $(1/2)(q_{1}/q_{2}+q_{2}/q_{1})(\hat{q}_{1}\cdot\hat{q}_{2})$. This extra layer will not alter any of our results as it can be incorporated as a trivial sum at the end of our calculations.

The fundamental set of integrals we will pursue here thus becomes

[TABLE]

and for compactness of notation we define all terms but the outer sums and the $\ell$-dependent pre-factor as $\mathcal{I}_{ij,\ell m}(k)$ so that

[TABLE]

We will focus on the $\mathcal{I}_{ij,\ell m}$ as then summing over $\ell$ and $m$ at the end is a trivial additional step. We now introduce a kernel $K$ such that

[TABLE]

In other words, we may take advantage of the fact that the three integrals (15, 24, and 33) differ only with respect to the argument of the final power spectrum to focus on integrals of that final power spectrum against the same kernel $K$ for all three.

We seek to cast each of the $\mathcal{I}_{ij,\ell m}$ as double 3-D convolutions, with the inner integral over $\vec{q}_{2}$ a 3-D convolution the result of which is then further convolved with additional factors as we perform the outer integral over $\vec{q}_{1}$. In short, we wish to render the computation “doubly convolutional.”

Suppose we could find a way to split $K$, perhaps differently for each of the three 2-loop terms, into a piece that would be explicitly doubly convolutional and another piece that, if we could factor it into separated functions of $\vec{k},\vec{q}_{1},$ and $\vec{q}_{2}$, could also be incorporated into the double convolutions. Pursuing such a splitting, we write

[TABLE]

where ${\rm c}$ is for “convolutional” and ${\rm d}$ is for “decoupled.” We have added subscripts to the righthand side to indicate that the required splitting may be different for each of the three 2-loop contributions.

For $\mathcal{I}_{15}$, we set

[TABLE]

where we have defined

[TABLE]

with $\mathcal{R}$ for “radial” and $\mathcal{P}$ for “power law”. We have included a subscript in the argument of $\mathcal{R}$ becauase there it should match the superscript of $N$. We set

[TABLE]

We have rearranged the factors in equation (8) relative to those in equation (4) so as to group the terms in $\vec{q}_{1}$ together and the same for those in $\vec{q}_{2}.$ $K_{15}^{\rm c}$ is clearly convolutional in that we could integrate over $\vec{q}_{2}$ at an offset of $\vec{k}+\vec{q}_{1}$. We could then integrate that result against $\vec{q}_{1}$, leading to a convolution over $\vec{q}_{1}$ at an offset $\vec{k}$. We note that the third power spectrum in equation (4), which for $\mathcal{I}_{15}$ becomes $P_{\rm lin}(k)$, can be pulled outside these integrals.

For $\mathcal{I}_{24}$ we are motivated by the argument of the third power spectrum (i.e. $\vec{w}_{24}$ in equation 2) to define

[TABLE]

and

[TABLE]

$K_{24}^{\rm c}$ is convolutional: we may first convolve just the $\vec{q}_{2}$-dependent factors over $\vec{q}_{2}$ at an offset of $\vec{k}$, obtaining a result dependent only on $\vec{k}$. We note that the power spectrum has argument $\vec{k}+\vec{q}_{2}$ (see equation 4) and so could also be included in this inner convolution. We may then perform a separate convolution of the $\vec{q}_{1}$-dependent factors over $\vec{q}_{1}$. We would then finally integrate the product of the inner and outer convolution results over $\vec{k}$.

Finally, motivated by $\vec{w}_{33}=\vec{k}+\vec{q}_{1}+\vec{q}_{2}$ in the third power spectrum in $I_{33}$, we set

[TABLE]

We note that the power spectrum in $\vec{k}+\vec{q}_{1}+\vec{q}_{2}$ can be included in the inner convolution.

3 Decoupling

We now have two tasks. First, we need to obtain explicit representations of the $K^{\rm d}$ that are separated into functions of one of its three arguments each. Second, we need to factor both the convolutional ($K^{\rm c}$) and decoupled ($K^{\rm d}$) kernels into radial and angular pieces. We may then exploit the isotropy of the power spectrum, the only component that is not available in closed form, to evaluate the angular dependence analytically. Doing so will reduce the numerical integrals we must do to 1-D radial transforms that can be done efficiently.

We discuss each task in turn, first considering the decoupling. Schmittfull & Vlah (2016) have the powers $n_{t}$, $t\in\{1,1^{\prime},2,2^{\prime},3,3^{\prime}\}$, take on values of $1$ and $2$ (see equation 1), so we need to decouple denominators of the form

[TABLE]

where $\vec{p}_{1}$ and $\vec{p}_{2}$ stand in for any of $\vec{k},\;\vec{q}_{1}$, and $\vec{q}_{2}$.

We take motivation from Limpanuparb (2011) and Dominici et al. (2012). These works use a multipole expansion to write

[TABLE]

However, in equation (15) the two momentum magnitudes are only ${\it formally}$ decoupled; in reality, one must be constrained by the other when we integrate to insure that the lesser-greater restriction on the right-hand side is met. Thus one cannot do the radial integrals over $p_{1}$ and $p_{2}$ in a truly separated way, so the computational cost remains $N_{\rm g}^{2}$ rather than $2N_{\rm g}$, with $N_{\rm g}$ the number of grid points in the $p_{t}$, $t=1,2$ used for the integration. We require a means of truly decoupling the radial part of the multipole expansion (15).

Limpanuparb (2011) and Dominici et al. (2012) do this via the integral

[TABLE]

We prove this integral in Appendix C. Rewriting the radial piece of the multipole series (15) via the integral (16) thus offers a true rather than merely formal decoupling: but at the price of an additional integral over $x$ at the end.

Limpanuparb (2011) and Dominici et al. (2012) circumvent this issue by finding that the integral is exactly equal to an infinite sum. This is essentially a sampling rule for performing the integral numerically; one can take this identity to mean that taking the integrand at integer values only of $x$ is the best integration scheme. In particular, Dominici et al. (2012) show the integral-sum identity that

[TABLE]

where we have taken their Corollary 3.3 and rewritten in terms of spherical Bessel functions using that $j_{L}(x)=\sqrt{\pi/(2x)}J_{L+1/2}(x)$. Following Dominici et al. (2012) we have defined

[TABLE]

The identity (17) leads them to the decomposition

[TABLE]

which has the additional advantage of being factorized into radial and angular pieces. We term the integral (16) a “decoupling integral.” We note that the integral-to-sum identity is only valid in the range $0<p_{1},\;p_{2}<\pi$, but since the power spectrum has effectively compact support due to the damping from non-linear structure formation, we can always rescale our integration domain of momentum magnitudes appropriately.

To adopt this approach to decoupling here, we need to solve three problems. First is to find an appropriate series expansion of $1/|\vec{p}_{1}+\vec{p}_{2}|^{2}$ and $1/|\vec{p}_{1}+\vec{p}_{2}|^{4}$ to be used in the kernels $K^{\rm d}$, in essence a generalization of the multipole expansion for powers other than the inverse first power. Second is to find a decoupling integral giving the right ratio of powers of the momentum magnitudes for the series developed in the first step. Third is to hope that this decoupling integral has an exact integral-to-sum identity like that Limpanuparb (2011) and Dominici et al. (2012) exploit.

These three problems are surmountable and in the following subsection we discuss how.

3.1 Decoupling $K^{\rm d}$

Motivated by the discussion above, we adopt the generalization of Legendre polynomials, Gegenbauer polynomials $C^{(\lambda)}_{L}$, and use the expansion

[TABLE]

(e.g. Arfken et al. 2013 or Sack 1964 equation 4). It is clear that this form will cover the required cases for powers of the vector sums we wish to decompose, but we also need to decouple the angular dependence. This decoupling can be done using the addition theorem for Gegenbauer polynomials (e.g. Koornwinder 1977 equation 3.1). This theorem is a special case of the addition theorem for Jacobi polynomials $P^{\nu,\mu}_{n}$, which reduce to Gegenbauer polynomials for $\nu=\mu$. However, the Gegenbauer polynomial addition theorem expresses the Gegenbauer polynomial of a dot product as a sum over products of Gegenbauers of the individual angles as well as powers of sines of the individual angles.

For our case it will be more convenient to use a “mixed” addition theorem we prove in Appendix B, which is to our knowledge novel to this work. We decompose a Gegenbauer polynomial of a dot product $x\equiv\hat{a}\cdot\hat{b}$ into a finite sum of products of spherical harmonics each of only one unit vector. We have

[TABLE]

where $w_{J}^{L,\lambda}$ is a constant coefficient defined in equation (130).

We note that a seeming alternative would have been to use parametric differentiation with respect to $\cos\theta_{12}\equiv\hat{p}_{1}\cdot\hat{p}_{2}$ on the Legendre series (15), having written $|\vec{p}_{1}+\vec{p}_{1}|=\sqrt{p_{1}^{2}+p_{2}^{2}+2p_{1}p_{2}\cos\theta_{12}}$ in the denominator. However this will only raise the denominator’s power by even steps, so we could generate $|\vec{p}_{1}+\vec{p}_{2}|^{-3},|\vec{p}_{1}+\vec{p}_{2}|^{-5}$, etc. but never the even powers (2 and 4) we need.

We now seek the decoupling integral to rewrite the radial part of the expansion (20), concentrating first on $\lambda=1$. After some experimentation, we find that the integral of the symmetrized sum

[TABLE]

The second term integrates to zero when $p_{1}>p_{2}$, and the first term integrates to zero when $p_{2}>p_{1}$, so in combination they always provide the desired ratio. We prove this integral in Appendix D.

We now need a suitable integral-to-sum identity to convert the decoupling integral (22) to a sum. Dominici et al. (2012) supplies one, though with the restriction that one free argument (e.g. $p_{1}$) be less than the other (e.g. $p_{2}$). Of course this is exactly the restriction we sought to relax in moving to a decoupling integral in the first place. However it can be shown that the restriction under which Dominici et al. (2012) proves the integral can be relaxed to cover the cases we need; details are given in Appendix E.

The integral-to-sum identity discussed there, with our extension to cover both $p_{1}<p_{2}$ and $p_{2}<p_{1}$, gives

[TABLE]

Consequently we see that

[TABLE]

with

[TABLE]

$\phi^{2\pm}_{nL}$ represents the two radial eigenfunctions. Superscript $2$ indicates an inverse-square expansion, $+$ a spherical Bessel function of greater index $(L+1)$ and $-$ one of lesser index $(L)$.

We can generate all of the other even powers by parametric differentiation with respect to $\cos\theta_{12}\equiv\hat{p}_{1}\cdot\hat{p}_{2}$. Expanding the magnitude in the denominator using the binomial theorem we notice that

[TABLE]

so that we may solve for the desired inverse fourth power as

[TABLE]

We may then insert our representation for $1/|\vec{p}_{1}+\vec{p}_{2}|^{2}$ in terms of $p_{<}/p_{>}$ and Gegenbauer polynomials (equation 20 with $\lambda=1$) on the righthand side of equation (28) and apply the differential operator to just the Gegenbauer polynomials as they are the only factor dependent on $\cos\theta_{12}$.

We prefer to use parametric differentiation rather than directly changing $\lambda$ from unity to two in equation (20), as this latter would alter the radial structure of the expansion and so render our particular decoupling integral (22) inapt. In contrast, as noted above, parametric differentiation with respect to $\cos\theta_{12}$ acts only on the Gegenbauer polynomial piece of the expansion, preserving the radial structure and hence the utility of our particular decoupling integral. Handled in this way, the radial structure of the expansion for $1/|\vec{p}_{1}+\vec{p}_{2}|^{4}$ will be the same as for that of $1/|\vec{p}_{1}+\vec{p}_{2}|^{2}$ save for multiplying by $1/(2p_{1}p_{2})$.

Explicitly, we have

[TABLE]

We can use the relation that

[TABLE]

to conveniently evaluate this parametric differentiation while remaining in the basis of Gegenbauer polynomials. This latter is desirable so that we may still apply our “mixed” addition theorem (130) splitting the Gegenbauer polynomials into spherical harmonics to factorize the angular dependence. We now simply alter the coefficients $w_{J}^{L,\lambda}$ to $2w_{J}^{L-1,\lambda+1}$, and the $2$ will cancel with the $1/2$ from equation (28). We also note that we do not want the sum (29) to now involve Gegenbauer polynomials of negative order. Of course, the derivative of the zero-order Gegenbauer polynomial (a constant) vanishes, and so it does not enter the sum. So we leave the indexing of the sums the same but take it that $w_{J}^{-1,2}=0$.

We find the factorized eigenfunction expansion for an inverse fourth power as

[TABLE]

Parallel to our notation in equation (24), here superscript $4$ indicates an inverse-fourth-power expansion, $+$ the spherical Bessel function of greater index, and $-$ that of lesser index.

Examining equations (24) and (3.1) we see that the eigenfunctions may all be conveniently written as

[TABLE]

with $\alpha=2$ and $4$ for the inverse square and inverse fourth power expansions respectively. This notation enables writing our expansions for both the inverse square and inverse fourth powers as

[TABLE]

where we also rewrote the superscripts on the weight $w$ so that they would correctly handle either value of $\alpha$.

Finally, we now write out the full decoupled kernels $K^{\rm d}$ by inserting the expansion (33) into equation (10). We have

[TABLE]

Relative to equation (10) we have chosen to use $\alpha$ in place of $2n_{3}$ and $\alpha^{\prime}$ in place of $2n_{2}^{\prime}$; $\alpha$ and $\alpha^{\prime}$ can each take on the values $2$ and $4$ as needed and will point to the correct eigenfunctions (24) and (3.1). We have also rearranged factors within each term relative to equation (10) to keep the $\vec{q}_{2}$-dependent factors together. Finally, the sums over $n,\;n^{\prime},L,$ and $L^{\prime}$ run from zero to infinity, but the sums over $S$ and $S^{\prime}$ are bounded by $J$ and $J^{\prime}$, which are in turn bounded by $L$ and $L^{\prime}$. The basic structure is the same as that of a spherical harmonic expansion of a Legendre series but then with two additional indices $(n$ and $n^{\prime})$. In what follows we will not write out the explicit bounds of the sums to save space.

In equation (34) we have also defined

[TABLE]

To keep the notation compact moving forward, we also define the product

[TABLE]

Upper indices correspond to the first argument, lower to the second. We will write $\Pi$ with just one argument when both factors on the righthand side have the same argument, as will be the case for $\Pi$ from the $\vec{q}_{2}$ terms in equation (34). With these definitions, we find

[TABLE]

Above we have suppressed the $nLJS,\;n^{\prime}L^{\prime}J^{\prime}S^{\prime}$ lower and upper indices on $\Pi$ since they are the same for all of the factors in each of the four terms above. $*$ denotes conjugate and should be taken to apply to the entire $\Pi$, not just the part indexed by $\alpha$. As discussed near equation (13), equation (34) also covers the 33 decoupling, as $K^{\rm d}_{33}=K^{\rm d}_{15}$.

For completeness, we now write out the expression for $K_{24}^{\rm d}$. By comparing equations (10) and (12), we see that leaving the $\alpha$ terms unchanged and replacing $\vec{k}$ with $\vec{q}_{1}$ and $\vec{q}_{2}$ with $\vec{k}+\vec{q}_{2}$ in the terms in $\alpha^{\prime}$ in equation (37) will give the desired result. It is

[TABLE]

We note that with these replacements, in contrast to equation (37), the first factor of $\Pi$ in each term now just has one argument, $\vec{q}_{1}$, while the second factor in each term now has two, $\vec{q}_{2}+\vec{k}$ and $\vec{q}_{2}$. With these factorizations of the $K^{\rm d}$ in hand, we are now ready to return to the full problem of evaluating $\mathcal{I}_{ij,\ell m}$. It is acceptable to have $\vec{q}_{2}+\vec{k}$ in equation (38) as this is a convolutional variable for $\mathcal{I}_{24}$.

4 2-Loop Terms as 3-D Convolutions

We now have

[TABLE]

We have defined

[TABLE]

where $\mathcal{R}$ and $\mathcal{P}$ are defined in equation (9). We have omitted the parameters and arguments of $T$ in favor $\cdots$ in all of the $T$ in equation (39) save for the first because they are exactly the same as those in the first. We note that $T$ has all of the power law indices involved in $I_{15}$ in its argument save for $n_{2}^{\prime}$ and $n_{3}$; these instead enter through respectively $\alpha^{\prime}$ and $\alpha$. $JJ$ and $SS^{\prime}$ are indices on the $\Pi$ but have been suppressed.

Inserting the more detailed forms for $\Pi$ of equation (36), we can write $T$ in a form that will facilitate analytic evaluation of the angular pieces:

[TABLE]

where we have defined the inner integral

[TABLE]

As we see from equations (40) through (42), while equation (39) looks involved, each of the four double integrals given by the $T$ is of the same form, and just comes from multiplying the four terms of equation (37) with the factors in equation (8) for $K^{\rm c}_{15}$. Equations (39), and (41) with the definition (42) are main results of this paper, and display the $15$ contribution to the 2-loop power spectrum as a sum of nested 3-D convolutions.

We now write down the analogous results for the $24$ and $33$ contributions. We treat the $33$ contribution first as it is more similar to the $15$ contribution than is the $24$ one. The only change in going from the $15$ to the $33$ contribution is that what was formerly the $P(k)$ pre-factor now has a different argument (see equation 2) and so must come inside the integrals. We have

[TABLE]

with

[TABLE]

and

[TABLE]

We now turn to the 24 contribution. Comparing equations (10) and (12), we see that $n_{2}^{\prime}$ becomes $n_{3}^{\prime}$ in the definition of $\alpha^{\prime}$, and consequently any $n_{3}^{\prime}$ that was written explicitly must now be replaced with $n_{2}^{\prime}$. From this we also see that the unprimed $\alpha$ terms are the same, but the primed ones now have $\vec{k}$ becoming $\vec{q}_{1}$ and also $\vec{q}_{2}$ becoming $\vec{k}+\vec{q}_{2}$. And finally from equation (2) the linear power spectrum with argument $\vec{k}+\vec{q}_{1}+\vec{q}_{2}$ in the $33$ integral now has argument $\vec{k}+\vec{q}_{2}$. We thus have

[TABLE]

with

[TABLE]

and

[TABLE]

5 Evaluating the Convolutions

We now evaluate the $T_{ij}$. For each of the cases (15, 33, and 24) we split the work into two steps, the inner convolution first and then the outer. We give detailed derivations for the 15 case and simply quote results for the other two as the steps to follow are analogous.

5.1 15: Inner convolution

We first focus on the inner integral over $\vec{q}_{2}$, defined in equation (42). We see that

[TABLE]

where star denotes convolution. The convolution is an integration over $\vec{q}_{2}$ but we have suppressed these arguments inside the square brackets and will do so throughout.

To evaluate the convolution, we use the Convolution Theorem to write

[TABLE]

We have again suppressed $\vec{q}_{2}$, the arguments of the functions being inverse Fourier-Transformed, but for clarity have included $\vec{r}$, the argument the result of the inverse FT will have in configuration space, as well as $\vec{k}+\vec{q}_{1}$, the final argument once we return to Fourier space.

We now evaluate the required inverse and forward FTs. We first treat the more complicated one, the first convolvand in equation (50). We have

[TABLE]

where we have defined

[TABLE]

To obtain the second equality we used the plane wave expansion (104) to expand the complex exponential into spherical Bessel functions and spherical harmonics, leading to the $Y_{L_{2}M_{2}}$ and the $j_{L_{2}}$ inside the radial integral in equation (51).

We see that the angular integral in equation (51) is just the overlap integral of four spherical harmonics, which can be evaluated in terms of Gaunt integrals, as we show in Appendix A. We notice that since the result is real, we may take the conjugate of the whole angular integral without altering it, and so we can rewrite it as the integral over four unconjugated spherical harmonics.

Meanwhile, the radial integral over $q_{2}$ is just a triple spherical Bessel transform of the linear power spectrum weighted by powers of $q_{2}$ as well as whatever $q_{2}$ dependence entered the numerator as in equation (3). We define it as

[TABLE]

where we take this as the definition of the $f$-tensor on the lefthand side. We note that the minimum power-law weight applied to the linear power spectrum in this tensor is from setting $\alpha$ and $\alpha^{\prime}$ to their maxima, $4$, and setting $n_{2}$ to its maximum as well, $2$. In this case we find a weight of $q_{2}^{-13/2}$, to which we add $2$ to account for the radial Jacobian in spherical coordinates. The power spectrum scales as $q_{2}$ in the infrared limit, and the spherical Bessel functions scale as $q_{2}^{L+L^{\prime}+L_{2}}$ in this limit. We thus have a net power of $[2(L+L^{\prime}+L_{2})-7]/2$: the $f$-tensor is infrared-divergent for low orders of the spherical Bessel functions. However, in practice the numerical integrations will be over a finite range, and as long as this is done carefully, these divergences should not affect the end result. Indeed, FFTLog (Hamilton, 2000), one acceleration method for performing these integrals, biases the integrands by power-law weights in any case and so will not see the divergence numerically at all. We will discuss this in more detail in §6.

Now, the maximum power-law weight in the integrand occurs when $\alpha$ and $\alpha^{\prime}$ are at their minima, $2$, as is $n_{2}$, $n_{2}=1$, leading to an overall power law weight of $[2(L+L^{\prime}+L_{2})+1]/2$; this is not divergent in the infrared. We do not concern ourselves with high-$q_{2}$ (ultraviolet) divergences as we assume that to avoid artifacts from a finite and discretized integration grid, the linear power spectrum will always be multiplied by a smoothing Gaussian, i.e. $\exp[-q_{2}^{2}\sigma^{2}]$ with \sigma$$\simeq$$1 Mpc/$h$.

With the notation outlined above in hand, our full result for the inverse FT of the first, more complicated convolvand in equation (50) is

[TABLE]

$\mathcal{H}$ is the integral of four spherical harmonics and is defined in equation (113). Noticing that the $f$-tensor is spin-independent, as is the weight $\mathcal{W}$, we see that at the end of the calculation we will be able to sum over the spins $S$ and $S^{\prime}$ analytically (we recall equation 39). This will also prove true of $M_{2}$ in the end. We will also be able to sum over $J$ and $J^{\prime}$ analytically because in equation (39) the only other dependence on these is through the weight $\mathcal{W}$ (itself defined in equation 35).

We now compute the inverse FT of the second convolvand in equation (50). For $n_{3}^{\prime}=2$ it scales as $1/r$; for $n_{3}^{\prime}=4$ it is formally divergent so we have regularized with an infrared cutoff of the integral at $\epsilon$; we may take the limit $\epsilon\to 0$ at the end of the full calculation. Details are discussed in Appendix F, as well as an alternative regularization scheme. The choice of regularization does not affect the structure of the calculation and so either may be used as dicated by numerical work. We find

[TABLE]

where ${\rm Si}(x)=\int_{x}^{\infty}dq\;\sin q/q$ is the sine integral. A key aspect of the expressions in equation (55) is that neither has any angular dependence, so they will not alter the angular momentum structure of the product in configuration space. As $r\to 0$, $g$ diverges as $1/r$ for $n_{3}^{\prime}=1$ but tends to zero for $n_{3}^{\prime}=2$.

We now write down the product of the two inverse FTs for the inner integral equation (50) as

[TABLE]

We defined

[TABLE]

and also set

[TABLE]

to simplify notation and facilitate the outer convolution we will perform in the next section (as dicated by equation 39). We note that the IR divergence structure of $h$ is given by combing that of $f$ (equation 53) and $g$ (equation 55) with the small-argument behavior of the sBF in $h$. We find that $h$ scales as $q_{2}^{L_{2}-1+[2(L+L^{\prime}+L_{2})-7]}$ as $q_{2}\to 0$, in the worst-case scenario where $f$ and $g$ diverge most strongly.

The other inner integrals, i.e. equation (49) but with the superscripted signs changed to $+-$, $-+$, and $++$, are generated by raising either the first index $L$, second index $L^{\prime}$, or both indices of the $h$-tensor by unity. This is the only change, but for completeness we write out the full expressions below.

[TABLE]

[TABLE]

and

[TABLE]

5.2 15: Outer convolution

5.2.1 Core calculation

We now turn to performing the outer convolution. From equation (39) we have

[TABLE]

We have suppressed the arguments of the convolvands. Using the Convolution Theorem we may write

[TABLE]

We use $\vec{x}$ as our configuration-space variable above rather than $\vec{r}$ to distinguish from our work on the inner convolution. Evaluating the first inverse FT in equation (5.2.1) we obtain

[TABLE]

$\mathcal{G}$ is a Gaunt integral, defined in equation (109). We have also defined

[TABLE]

we recall from equation (9) that $\mathcal{R}$ contains a linear power spectrum and so the integal above cannot be performed analytically. This $f$-tensor is a limit of the more general three-index one in equation (53) with $\alpha^{\prime}=0$, $L^{\prime}=0=n^{\prime}$, $n_{2}\to n_{1}$, $L\to L+1$, $L_{2}\to L_{1}$, $r\to x$, $q_{2}\to q_{1}$. Thus the analysis of the divergences presented there can be used here as well.

We need only to treat one additional case, that with a negative sign in the superscript of $\phi$. The result for that case is

[TABLE]

i.e. only the order of the spherical Bessel function with argument $nq_{1}$ changes, from $L+1$ to $L$.

We now obtain the inverse FT of the second term in equation (5.2.1). Noticing that the multiplication of $I_{15}$ by $\mathcal{P}^{n_{3}^{\prime}}$ does not alter the angular structure of the expansion, we see we will have an inverse FT at the same order sBF as is in $I_{15}$. We find

[TABLE]

we have suppressed the momentum argument of $I_{15}$, consistent with our earlier convention for inverse FTs, but retain the parameters $n,n^{\prime},$ etc. for clarity. We defined

[TABLE]

The IR behavior of $H$ is as that of $h$ with additional weight $q_{2}^{L_{2}+2(1-n_{1}^{\prime})}$.

Thus the configuration-space product of the two inverse FTs (equations 64 and 67) required by equation (5.2.1) is:

[TABLE]

We notice that $L_{1}$ and $L_{2}$ are controlled by $\ell,J$, and $J^{\prime}$, and $J$ and $J^{\prime}$ are in turn respectively controlled by $L$ and $L^{\prime}$. These latter two are just the angular momentum indices in our eigenfunction expansions of the two denominators that were expanded into Gegenbauer polynomials. We do not quote the intermediate results for the other combinations of $\phi^{\alpha\pm}$ and $I_{15}^{\alpha\alpha^{\prime}\pm\pm}$ as these can easily be generated by raising or lowering appropriate angular momentum indices; we will state all results for $T$ at the end of the calculation.

We now take the FT from $\vec{x}$ to $\vec{k}$ as dictated by equation (5.2.1), and observe that we will get a third angular momentum and spin, $\mathcal{L}$ and $\mathcal{M}$, from expanding the plane wave entering this FT. However, since we will integrate $T$ over $d\Omega_{k}$, we will see that $\mathcal{L}\mathcal{M}\to 00$ and so $L_{2}\to L_{2},\;M_{2}\to M_{1}$. But we do not apply this simplication yet. We have for $T$

[TABLE]

The IR behavior of $T$ is as that of $f$ (equation 65) combined with $H$ (equation 68) and weighted by an additional $x^{\mathcal{L}+2}$.

We now need the three other terms, corresponding to the other options for the plus and minus signs in the superscript of $T$. Changes to the innermost two signs affect only $H$, while the outermost signs affect respectively $f_{L+1,L_{1}}^{\alpha\ell}(n;n_{1};x)$ and $\phi_{n^{\prime}L^{\prime}}^{\alpha^{\prime}+}(k)$. Writing out the additional terms entering equation (39) we find

[TABLE]

[TABLE]

and

[TABLE]

We now integrate against $d\Omega_{k}/(4\pi)$ (angle-average) as required by equation (39). This cancels the leading $4\pi$ above and sets $\mathcal{L}=0=\mathcal{M}$, which in turn requires that $L_{1}=L_{2}$ and $M_{1}=M_{2}$. This brings the relevant Gaunt integral to $(-1)^{M_{1}}/\sqrt{4\pi}$ (see equation 109) using Olver et al. (2010) §34.3.1. to evaluate the two 3j-symbols that enter it.999https://dlmf.nist.gov/34.3 Denoting the angle-averaged $T$ as $\bar{T}(k)$, we find

[TABLE]

We notice that in equation (39), the only dependence on $S$ and $S^{\prime}$ is in the $T$, and so we may now perform the sums over these spins. We then may make the definition

[TABLE]

where we used the definition of $\mathcal{H}$ equation (113) to obtain the second line; defining $\omega$ is desirable to simplify notation moving forward. We thus have

[TABLE]

So we have expressed $T$ in terms of successive transforms against spherical Bessel functions of linear power spectra weighted by power-laws. Returning to equation (39) we see that this means that the full 15 contribution can be written in this way. The expressions for $T$ above are a second major result of this work: they mean that the 15 contribution to the 2-loop power spectrum can be expressed purely as a sum of successive 1-D integral transforms.

5.2.2 Control of Angular Momenta

Before going on to address the 33 and 24 contributions, which will have a similar structure, we briefly discuss the structure of our results thus far. We have sums over $n,n^{\prime},L,L^{\prime},J^{\prime},L_{1},L_{2},$ and $\mathcal{L}$, with $\ell$ as an additional angular momentum. We discuss their provenance in the order they were produced. $\mathcal{\ell}$ came from expanding an arbitrary numerator in our original integral $I_{15}$; in practice these numerators will have very compact support in $\ell$, so the range of $\ell$ is controlled. The sums over $L$ and $L^{\prime}$ came from our expansion of the coupled denominators into a sum over Gegenbauer polynomials (the radial piece of which was then decoupled using the decoupling integrals). The sum over $J^{\prime}$ came from the fact that at each $L$, the Gegenbauer polynomial was expanded into a finite number of spherical harmonics. The range of $J^{\prime}$ at fixed $L$ is thus finite. The sums over $n$ and $n^{\prime}$ came from our expansion of the decoupling integral into a product of sBFs.

We now notice that all angular momenta in the problem are controlled by $\ell,L$, and $L^{\prime}$. Due to the triangle inequality on the 3j-symbol, $L_{1}$ is controlled to be $|J^{\prime}-\ell|\leq L_{1}\leq|J^{\prime}+L_{1}|$. Returning to equation (54), we see that $\mathcal{H}$ involves $L_{2}$. The three other angular momenta in $\mathcal{H}$ are controlled at fixed $L$: $\ell$ comes from the numerators, as just discussed, $J=0$ and $J^{\prime}$ is finite and ranges from [math] or $1$ up to $L$ or $L^{\prime}$ depending on $\alpha$ (i.e. whether we expanded an inverse square or inverse fourth power). So $L_{2}$ is controlled at fixed $L$. Since $L_{1}$ and $L_{2}$ are controlled, $\mathcal{L}$ is by the triangle inequality on the 3j-symbols. So we see that deciding at which $L,L^{\prime},n,$ and $n^{\prime}$ to truncate our expansion of the denominators we expanded will fully fix the extent of the rest of the sums. Fortunately, Dominici et al. (2012) find that these expansions converge for a small number of terms, of order less than 10.

We make three further comments. First, the $f$-tensors are just multiple sBF transforms of the linear power spectrum weighted by power laws; therefore they should be amenable to the Limber approximation as the order of the sBFs grows large.101010The Limber approximation takes it that as their order grows the sBFs become similar to Dirac Delta functions about the point where they first begin to oscillate; Limber (1953). Thus, if going to a large number of terms in the expansion of the denominators were required (it is likely not), it could possibly be done without the computational cost of more radial integrals. This also holds true of our final integral over $j_{\mathcal{L}}(kr)$.

Second, $\ell$ is fixed by the numerators in the original integral, and should be small, as the 2-loop PT kernels do not have complicated angular structure in their numerators. Meanwhile $J^{\prime}$ can range up to $L^{\prime}$, which itself is controlled by our truncation. If we needed to continue to larger $L^{\prime}$ before truncating for convergence, the Wigner 3-j symbol would become very squeezed, with J^{\prime}$$\sim$$L_{1}\gg\ell; these symbols become very small in this limit, scaling as $1/\sqrt{J^{\prime}+L+1}$, meaning that terms at high $L$ likely contribute rather less to the result. Since $\mathcal{H}$ couples $\ell,J^{\prime}$ and $L_{2}$ in the same way as $J^{\prime}$, $\ell$, and $L_{1}$ are coupled (since $J=0$, $\mathcal{H}$ reduces to just 3j-symbols), this conclusion holds for $L_{2}$ as well.

Third, the $f$-tensors may be evaluated using FFTLogs (see Hamilton 2000), as we discuss further in §6, so the scaling of the integrals becomes $N_{\rm g}\log N_{\rm g}$ rather than $N_{\rm g}^{2}$, where $N_{\rm g}$ is the number of points used for the $k$ and $r$ integrations each.

5.3 33: Inner convolution

We now turn to the 33 contribution to the 2-loop power spectrum. From equation (45), we read off that we need

[TABLE]

The first convolvand has the same form as the first convolvand in equation (49), so we can simply adopt equation (54) with the definition (53) of $f_{LL^{\prime}L_{2}}^{\alpha\alpha^{\prime}\ell}(n,n^{\prime};n_{2};r)$. The second convolvand is easily evaluated as

[TABLE]

with the last integral serving to define $\xi_{\rm lin}(n_{3}^{\prime};r;\epsilon)$ and the infrared cut-off $\epsilon$ only genuinely needed for $n_{3}^{\prime}=2$. For $n_{3}^{\prime}=1$, the integral converges with $\epsilon=0$, as $P_{\rm lin}(q_{2})\propto q_{2}$ and $j_{0}\to 1$ as $q_{2}\to 0$, so overall one has $q^{2}$ as the leading infrared behavior of the integrand. We have termed this result $\xi_{\rm lin}$ becuase it is simply the same transform as used for the linear correlation function, but here weighted by an additional power law specified by $n_{3}^{\prime}$.

Using these results we now have

[TABLE]

with

[TABLE]

The above are the analogs of equations (56) and (57), used for the inner convolution in $\mathcal{I}_{15}$.

5.4 33: Outer convolution

We now need the outer convolution for $\mathcal{I}_{33}$, the analog of equation (62). From equation (4) we have

[TABLE]

The first and second convolvands are both exactly the same as those in equation (62). For the first convolvand, we may therefore just adopt equation (64); for the second, equation (67). We make a new definition $H_{\xi,LL^{\prime}L_{2}}^{\alpha\alpha^{\prime}\ell}$ as our analog here of equation (68): it is just equation (68) but with the replacement $h_{LL^{\prime}L_{2}}^{\alpha\alpha^{\prime}\ell}\to h_{\xi,LL^{\prime}L_{2}}^{\alpha\alpha^{\prime}\ell}$ on the righthand side. Our results for the various required $T_{33}$ and $\bar{T}_{33}$ are then exactly as equations (70)-(73) but with $H\to H_{\xi}$. For the sake of clarity, we write out explicitly the analog of equation (70):

[TABLE]

with

[TABLE]

with $h_{\xi}$ defined in equation (80).

Overall, the 33 integrations shared the same angular momentum structure as the 15. The only difference in our work was the addition of one more linear power spectrum inside the integrals; the 15 contribution had two inside in total, but the 33 has three. This of course means that summing over spins and the simplifications leading to equation (76) proceed exactly in the same manner here, and we find

[TABLE]

This shows that $T$ for the 33 contributions to the two-loop power spectrum can be expressed in exactly the same form as those for the 15 contribution, namely as successive transforms against spherical Bessel functions of linear power spectra weighted by power-laws.

5.5 24: Inner convolution

We now obtain the inner convolution for the 24 contribution to the two-loop power spectrum. The angular structure is the same as for the other two contributions, so we may simply adopt equation (76) but with $H\to H_{\xi}$. We find

[TABLE]

This convolution is not the same as equation (49) but as already noted does have the same angular structure.

For the inverse FT of the first convolvand, we find

[TABLE]

with

[TABLE]

This result is the analog of equation (54) for the 24 contribution; the IR behavior can be read off from the discussion below equation (65).

For the second convolvand in equation (85), we have one more factor, $\phi_{n^{\prime}L^{\prime}}^{\alpha^{\prime}-}$, than in the analogous 33 piece (compare equations 85 and 78). We find

[TABLE]

The IR behavior of $\xi_{\phi}^{\alpha^{\prime}}$ is as $q_{2}^{9/2+L^{\prime}-2n_{2}^{\prime}-\alpha^{\prime}/2}$.

By comparing with our work for the 33 contribution, we see that defining $h_{\xi\phi}$ analogously to $h_{\xi}$ of equation (80) but with $\xi\to\xi_{\phi L^{\prime}}$ and with $f_{LL^{\prime}L_{2}}\to f_{LL_{2}}$ on the righthand side will allow us to adapt equation (79) to give the result here required. We find

[TABLE]

5.6 24: Outer convolution

We now present the outer convolution and full result for the 24 contribution. The 24 contribution is particularly easy to deal with as the outer and inner convolutions are independent of each other, because the inner convolution is independent of $\vec{q}_{1}$. From equation (4) the outer convolution is

[TABLE]

To obtain this result, we noticed several parallels with our previous work. First, the first convolvand is just the first term of equation (49) with $n_{2}\to n_{1}$, $--\to++$, and no conjugates on the spherical harmonics. The lack of conjugates here does not matter as the first convolvand of equation (49) and of the above are both real, meaning one can take the conjugate of all spherical harmoincs (the only complex part of either) wihtout altering the result. Second, the second convolvand above is just the second convolvand of equation (49) but with $n_{3}^{\prime}\to n_{1}^{\prime}$. Thus we could adopt equations (54) and (55) with these alterations for the two inverse FTs of the two convolvands above, and then use equation (56) with a conjugate to obtain our ultimate result (90). We simply multiply equation (90) by our result (89) for $I_{24,\ell m}$ to obtain $T$. We obtain

[TABLE]

Integrating over $d\Omega_{k}/(4\pi)$ sets $L_{2}M_{2}=L_{1}(m+S)$ by orthogonality of the spherical harmonics, eliminating the sums over $L_{2}$ and $M_{2}$ (and canceling a factor of $4\pi$). We represent this integrated $T$ as $\bar{T}$. Now summing over spins $S$ and $S^{\prime}$ we find

[TABLE]

where we defined the angular weight

[TABLE]

by analogy with $\omega_{\ell JJ^{\prime}L_{1}}^{m}$ of equation (75). The difference is that $\omega$ was a sum (over spins) of the product of an integral of four spherical harmonics with one of three harmonics, while here we have a sum (over spins) of the product of two integrals over four spherical harmonics. This is a direct result of the structural difference of the 24 contribution from the 15 and 33 contributions; it stems from that fact that for 24, the inner and outer convolutions are actually completely independent.

6 Discussion

Here we discuss a number of possible extensions of this work and other applications of the work. We then make some comments on how a numerical implementation might proceed.

6.1 Extensions

6.1.1 Incorporating Anisotropy

We now outline how this work could be extended to redshift space. Redshift-space distortions (RSD) generate a preferred direction (the line of sight), breaking isotropy, and typically multipoles of the power spectrum with respect to the angle to the line of sight, or wedges (bins) in angle, are used to quantify the anisotropy. Since RSD stem from peculiar velocities, the velocity kernels $G^{(i)}$ also enter the calculation. These have the same form as the density kernels $F^{(i)}$ but with different numerical pre-factors (see e.g. Bernardeau et al. 2002, equations 43-46). Thus there would be two changes to the calculations presented here.

First, one would replace the density kernels with the appropriate mix of velocity and density kernels; this would affect only the numerator $N$ of equation (4), and that trivially as the constant coefficients would come outside the integrals. Second, one would average over the orientation of the wave-vector $\vec{k}$ at which the power spectrum is measured with some weights, i.e. in equation (4) replace $d\Omega_{k}\to\mathcal{L}_{\ell}(\hat{k}\cdot\hat{n})$ with $\hat{n}$ the line of sight. This weight can be factored using the spherical harmonic addition theorem (100) to decompose the Legendre polynomial, and the integration over $d\Omega_{k}$ easily performed.

We note in particular that due to the azimuthal symmetry about the line of sight that remains even with RSD, the spins do not couple to this change, only the total angular momentum. Thus the simplification of summing $T$ over spins $S$ and $S^{\prime}$ (e.g. equation 76) still holds. We simply first promote $\bar{T}$ of equation (74) to $T_{0}$, the monopole moment of $T$ with respect to the line of sight, and also seek the higher multipoles $T_{\ell}$ with respect to the line of sight. Alternatively, it is also trivial to obtain wedges by averaging $T$ over bins in $\mu\equiv\hat{k}\cdot\hat{n}$; one just adds a binning function $\Phi(\mu;\bar{\mu}_{i})$ in equation (4) as an integration weight rather than the Legendre polynomials, with $\bar{\mu}_{i}$ denoting the $i^{th}$ angular wedge.

6.1.2 Lagrangian Perturbation Theory

We now briefly turn to consider extending this work to Lagrangian Perturbation Theory (LPT). There is no closed-form expression of order-by-order recursion relations in LPT (Bernardeau et al. 2002 §2.7); however recent work (Matsubara 2015) obtained approximate LPT kernels, and it would be worth considering whether evaluation using them benefits from the technique presented here.

We note that in the Zeldovich Approximation (ZA) within LPT, the recursion for the density at all orders takes on a particularly simple form where there are only linear momenta in the denominators, and the only appearance of a sum of momenta is in the numerator (Bernardeau et al., 2002). Thus the technique presented here is likely not needed for the ZA. Nonetheless we observe that the ZA integrals can likely be done as 1-D transforms by performing all of the angular integrations analytically and exploiting the isotropy of the power spectrum, one of the enabling ideas in this work as well.

6.1.3 Polyspectra

The techniques of this work could also be used to enable loop calculations of the bispectrum and trispectrum. At a given, fixed order in SPT, the kernels will be the same, just the combination of density fields will be different. But the key ideas here, factorization, decoupling to reduce to convolutions, and then evaluating the angular pieces analytically, should still enable bispectrum and trispectrum to be reduced to 1-D integrals. For instance, see Simonović et al. (2018) equation 3.6. for the one-loop bispectrum and Bertolini & Solon (2016) for the trispectrum. Furthermore, in many modified gravity (MG) models, the kernels are similar to those of SPT (Vernizzi et al. 2018, in prep.); these can thus also likely be dealt with in the same way.

6.1.4 Further Accelerating Cosmological Parameter Searches

Derivatives of the power spectrum with respect to the cosmological parameters might be used to further accelerate recomputation of the loop corrections for many different input cosmologies. Clearly, if one Taylor-expanded the linear power spectrum about some fiducial case with respect to variation in the cosmological parameters, one would have fundamental integrals of the derivatives evaluated at the fiducial case. These fundamental integrals would have the same kernels in the integrand as the integrals over the power spectrum computed in this work, as the kernels are cosmology-independent.111111At least for variations close enough to the fiducial cosmology that assuming Einstein-de Sitter remains valid; see Hivon et al. (1995). The variations in the cosmological parameters would sit outside the integrals. One could thus use the techniques outlined in this work to compute the additional required integrals. Similar ideas are discussed in Fendt & Wandelt (2007), Cataneo et al. (2017), and Lewandowski & Senatore (2018)).

6.2 Other Applications

6.2.1 Redshift-Space 3PCF, Hierarchical 3PCF Model

The techniques outlined here also have several other LSS applications. First, even in the tree-level redshift-space 3PCF (Slepian & Eisenstein 2017, equation 39) one has a term involving $1/k_{3}^{2}$, with $k_{3}=|\vec{k}_{1}+\vec{k}_{2}|$. There, the integration that occurs when transforming this to the multipole basis in configuration space (of $r_{1},r_{2}$, and Legendre polynomials in the cosine of their enclosed angle $\hat{r}_{1}\cdot\hat{r}_{2}$) was done by introducing a Delta function, rewriting it as the Fourier Transform (FT) of unity, and performing all of the angular integrals.121212See Szapudi (2004), Slepian & Eisenstein (2015a), and Slepian & Eisenstein (2017) for development of the multipole basis for the 3PCF. This procedure resulted in an infinite sum over terms each involving 2-D radial integrals in the end (see equations 16-18). However, using the full decoupling technique presented here, the denominator $1/k_{3}^{2}$ could have been decomposed into an eigenfunction expansion, leaving as an end result an infinite sum of 1-D integrals instead. Numerically, the $1/k_{3}^{2}$ contribution to the 3PCF multipoles turns out to be small, because of near-exact cancellation of the 2-D integrals in the result of Slepian & Eisenstein (2017) with each other; see discussion below equation 19 of that work. Nonetheless, at the level of precision DESI (DESI Collaboration et al., 2016) will likely offer on the 3PCF, even small systematic errors in the theoretical modeling should be avoided if possible. More importantly, we suspect that the expansion proposed here for these terms may shed light on why the nearly-exact cancellation making it so small occurs, offering additional insight on the physics of the 3PCF.

Regarding the physics of the 3PCF in the multipole basis, we further note that the Gegenbauer polynomial expansion and decoupling presented here also enables compactly evaluating the “hierarchical ansatz” for the 3PCF (Groth & Peebles 1977) in this basis. This ansatz takes it that the 3PCF $\zeta$ scales as a cyclic sum of products of 2PCFs $\xi$. At large separations $r$ these in turn scale roughly as power laws, $\xi\to 1/r^{2}$. We thus have, with $\xi_{i}\equiv\xi(r_{i})$,

[TABLE]

The phase factor of $(-1)^{J}$ comes from using parity to convert the spherical harmonic that would have argument $-\hat{r}_{2}$ into one with argument $\hat{r}_{2}$ (see equation 103).

We now project $\zeta_{\rm hier.}$ onto the multipole basis by integrating against $(2\ell+1)/(16\pi^{2})$. The factor of $(2\ell+1)/2$ corrects for the fact that the Legendre polynomials are orthogonal but not orthonormal, and we have an additional $1/(8\pi^{2})$ to deal with the redundancy of integrating over $d\Omega_{1}d\Omega_{2}$ when the Legendres only depend on the cosine of the relative angle, $\hat{r}_{1}\cdot\hat{r}_{2}$. To perform the integral we expand the Legendre polynomial using the spherical harmonic addition theorem (100) and invoke orthogonality. We thus find the Legendre expansion of the hierarchical-model 3PCF $\zeta_{\rm hier.}$ as

[TABLE]

with coefficients

[TABLE]

$\delta^{\rm K}$ is the Kronecker delta, unity when its subscripts are equal and zero otherwise.

Equation (96) is thus a compact expression for a toy-model of the 3PCF in the multipole basis based on the hierarchical ansatz plus the large-scale behavior of the galaxy correlation function. This toy model could be of use in understanding the scale dependence of $\ell>0$ multipoles of the 3PCF in the space of $r_{1}$ and $r_{2}$, such as are shown in e.g. Figure 9 of Slepian & Eisenstein (2015a) or Figure 7 of Slepian & Eisenstein (2017).

6.2.2 An $\mathcal{O}(N)$ $N$-body Integrator

Another application of the decoupling approach presented here is to $N$-body calculations. Naive evaluation of the gravitational forces a given set of $N$ particles (i.e. a discretization of the matter field) exert on each other scales as $N^{2}$. Most currently standard algorithms use approximate methods to avoid this scaling. For instance, one might compute forces from nearby particles directly using a pair count but farther away particles using a multipole expansion of their spatial coordinates (e.g. as Abacus does; Garrison et al. 2018, in prep.). However one must recompute the multipole expansion about every particle on which one wants the force (although Abacus uses the spherical harmonic shift theorem to mitigate this cost). Alternatively, one might compute the far-field forces using fast Fourier transforms (Ewald summation). This latter scales as $N_{\rm g}\log N_{\rm g}$ with $N_{\rm g}$ the number of grid points used for the FFT.

The eigenfunction expansion developed in this work offers an exact expansion of the force law between two particles in absolute coordinates. This latter point is key. By using the decoupling integral and then the integral-to-sum identity to derive this eigenfunction expansion, we removed the restriction that one distance from the global origin be less than another. Thus one can work always in an absolute coordinate system and never need to compare the distances of two particles from the global origin. This avoids a pair count.

In more detail, here is how an $N$-body algorithm based on this decoupling would proceed. Setting Newton’s constant $G=1$ and assuming that all particles have the same mass, we have the acceleration of the $i^{th}$ particle as

[TABLE]

With an eigenfunction expansion of the denominator, in this case in Legendre polynomials using parametric differentiation as discussed in §3, around equations (19), (27) and (28) (for the Gegenbauer case, but the same idea can be used to obtain an inverse-cube from the Legendre series), we have

[TABLE]

$J$ and $J^{\prime}$ run from zero up to $L$ and $L^{\prime}$ respectively and come from re-expressing the derivative of a Legendre polynomial as a series in lower-order ones with weights $w_{J}^{L}$ and $w_{J^{\prime}}^{L^{\prime}}$. These Legendres are then expanded into spherical harmonics. We may now rewrite the sum that excluded $i$ as one including $i$ and then subtract the value of the summand at $i$.

However we assume that the truncation of our series in $n,n^{\prime},L$, and $L^{\prime}$ will offer a natural softening so that the force kernel does not diverge as $\vec{r}_{j}\to\vec{r}_{i}$. In particular, Figure 1 of Dominici et al. (2012) shows that the sum approximating the integral, when truncated, does not develop a singularity at equal values of the free frequencies ($a$ and $b$ in that work, $r_{i}$ and $r_{j}$ here), at least for the orders of sBF they tested. This point is akin to the idea that if one expands a Dirac delta function in Fourier modes but limits the bandwidth, one will recover an oscillatory, smoothed function with infinite support. Given Parseval’s theorem, that the integral of a function is conserved from configuration space to Fourier space, the broadening implies that the singularity of the Delta function at the origin must now be finite. Consequently we suspect that truncating the sums in the eigenfunction expansion, as is required to compute them numerically in any case, will offer a “natural” softening.131313Given that softening has recently been found to cause problems in $N$-body simulations (van den Bosch & Ogiya, 2018), it may be of interest in and of itself to investigate whether the “natural” softening proposed here offers improvements over standard softenings. If this holds, then the value of the summand at $j=i$ will tend to zero, so that we can drop that term. We then find

[TABLE]

We now notice that, at fixed $i$, the two sums depending on $i$ are just numerical values, and can be pre-computed for every particle in the simulation box (i.e. stored for each $i$). The sums over $j$ are both over all particles in the box, and can therefore just be computed once. Furthermore, once the eigenfunctions $\phi$ and $Y_{JS}$ are found at all particle positions, these can be used to quickly perform both of these computations. The acceleration of the $i^{th}$ particle in the box at a given time-step is then just found using the simple linear combination of products these results as given above.

6.3 Comments on Numerical Implementation

We defer a full numerical implementation to future work, although it should be straightforward save for requiring some care to track all indices without error. However, we outline a few ideas on how an efficient numerical implementation might proceed.

We first sketch how even the most naive scheme could be rather efficient. Each successive transform required for our final results is just an integral of the power spectrum weighted by power laws against one or more spherical Bessel functions. This is exactly the type of integral done in CMBFAST (Seljak & Zaldarriaga 1996) to compute predictions for the temperature and polarization anisotropies of the Cosmic Microwave Background (CMB). There, the integration is split into a source function and geometric part, the source function being the perturbation theory potentials (for instance $h$ and $\eta$ in synchronous gauge, see e.g. Ma & Bertschinger 1995). The geometric part is the spherical Bessel functions, which can be pre-computed into a lookup table, meaning they need only be obtained once. Since the source functions are smooth but the spherical Bessel functions are oscillatory, these latter drive the cost of sampling and evaluation. However since the geometric part is known and does not change with changing source functions, one can optimize and precompute the sampling points for the integrations. This point is discussed in the CMBFAST implementation (Seljak & Zaldarriaga 1996), and there only about 50 points are required for the integrations.

In our case, the power spectrum plays the role of source function, and even despite BAO it will be rather smooth relative to the highly oscillatory sBFs. Thus in our case as well the sBFs will drive the integration points. A further advantage of our decomposition is that many of the spherical Bessel functions appearing in our expansions have integer frequencies; this makes locating the nodes especially easy, likely enabling further optimization of the sampling for the numerical integration. We note that the most naive integration scheme would scale as $N_{\rm g}^{2}$ with $N_{\rm g}$ grid points for the momenta $q_{i}$ and $k$ and the same for the configuration space-variables $r$ and $x$. At each momentum $q_{i}$ or $k$, one would need to evaluate an integral over all $r$ or $x$.

However, these integrals can be done scaling as $N_{\rm g}\log N_{\rm g}$ using a generalization of the FFTLog algorithm of Hamilton (2000), as mentioned briefly at the end of §5.2. The change of variable outlined there renders integrals against a single spherical Bessel convolutional, and this change of variable is actually agnostic as to how many sBFs one has, so can also be used to render integrals against multiple sBFs so. These ideas are further discussed in Li et al. (2018) and Li and Slepian 2018 (in prep.).141414https://github.com/eelregit/mcfit There are two advantages of this approach, in addition to the speed increase. First, it naturally handles logarithmically gridded input kernels, such as one would typically have for a power spectrum formed from transfer functions from e.g. CMBFAST (Seljak & Zaldarriaga, 1996), CAMB (Lewis et al., 2000), or CLASS (Lesgourgues, 2011). Second, the FFTLog transformation “biases” the integrands by power laws to take out any smooth, secular fall-off or rise and thereby increase the numerical precision of the integration. Hamilton (2000) shows that this biasing can then be undone analytically in the final space. So FFTLog will not actually “see” divergences in the integrals caused by power-law weights, which are exactly the divergences that arise in the approach of this work.

7 Connections to Other Work

Here we discuss connections to other work, in particular highlighting the power-law-basis approach to these integrals and also highlighting similar techniques used in QFT.

7.1 Complex Power-Law Basis

We first discuss a possible connection with the power-law work of Simonović et al. (2018). There, the power spectra are expanded as a sum of complex power laws (the imaginary part of the power laws enables capturing oscillations, such as BAO) and the integrals of the full PT kernels against complex power laws are done analytically. Thus when one alters the input power spectrum one only needs to update the expansion coefficients into the complex-power-law basis, and can then sum up the integrals (done already in closed form) with these new weights. The same approach could be applied here: writing $P_{\rm lin}$ as a sum of complex power laws would allow analytic evaluation of the 1-D integrals in our final results.

The advantage of this “hybrid” approach over that presented in Simonović et al. (2018) is that one would likely have easier results to evaluate than the hypergeometric functions required by their analytic integrals. The disadvantage is that some of these might diverge and some care might be required in treatment so as to obtain a convergent sum in all cases, and also that one would have several infinite (but countably so) sums. Details may merit attention in future work. However, taking that aside, we note that the formulae derived here must, by comparison with the results of Simonović et al. (2018), offer a new expansion of the hypergeometric functions into infinite sums of likely relatively simple terms. In particular, as noted above, mating our approach with a complex-power law expansion of the power spectrum, one would end up with integrals of complex power laws against pairs or triplets of spherical Bessel functions in place of the hypergeometrics appearing in Simonović et al. (2018). One can do these former in simple form, including resolution of the singular terms, as shown in Slepian (2018, in prep.) and Cahn and Slepian (2018, in prep.).

7.2 Multi-Loop Calculations in Quantum Field Theory

The present work also is related to techniques for performing multi-loop calculations in QFT. In particular, the first part of our approach here—expansion of certain denominators into a Gegenbauer series—is exactly what is done in the Gegenbauer Polynomial $x$-space Technique (GPXT) developed in Chetyrkin et al. (1980) and originally used to compute the first 3-loop results for $e^{+}e^{-}\to$ hadrons. GPXT was further developed in Chetyrkin et al. (1979), Celmaster & Gonsalves (1980), Terrano (1980), with some early calculations presented in Chetyrkin & Smirnov (1984), Lampe & Kramer (1983), Smirnov (1991), Chetyrkin & Tarasov (1994), and Kotikov (1996). However in QFT the Feynman integrals have no source functions in the integrand that are known solely numerically; all the functions in the integrand are written down in closed form. Thus the momentum-magnitude integrals may reduce to a set of nested integrals on a simplex (e.g. equations 2.8 and 2.9 of Chetyrkin et al. 1980), which can then be done just once analytically. In our case, where one may want to use many different linear power spectra, this would be inefficient. Avoiding the inefficiency of working over simplexes is where the utility of our decoupling integral and integral-to-sum identity lies.

We note that our approach does likely simplify some QFT integrals even though the integrand is completely known in closed form. For instance, in Bekavac (2006) GPXT is applied to a 5-loop calculation, leading to a result that has 120 different terms corresponding to the $5!$ different simplexes one needs to consider for the different orderings of the position-space lengths (Fourier dual of the momentum magnitudes) being integrated over. Although due to symmetries only 30 terms are genuinely different, this is still a substantial number. The eigenfunction expansion presented here, by breaking the simplicial coupling, enables avoiding individual treatment of each simplex, potentially reducing the computational complexity.

There are several other works in QFT that seem conceptually related to the present one. Bollini & Giambiagi (1996) exploits the spherical symmetry of the Feynman propagators in Fourier space to rewrite them as $j_{0}$ transforms of configuration space functions, and then performs the convolution (e.g. of two) as a product in configuration space. This product gives two sBFs. Exploiting the spherical symmetry of the result and transforming back to Fourier space yields a third. The transform to Fourier space also gives an integration over the configuration space variable: this can be performed analytically as the overlap integral of the three sBFs. One then has integrals left over for the two initial $j_{0}$ transforms now weighted by a coupling kernel from the triple-sBF overlap integral. This coupling kernel has a geometric interpretation: it is proportional to the area of a triangle whose sides are the external momentum and these two internal momenta. This work is related insofar as it exploits spherical symmetry and therefore works in terms of sBFs; however the result, of a coupled integral over momenta, is exactly what we here wished to avoid as this would couple the integrals over power spectra that must be done numerically, leading to a 2-D or more grid for numerical integration.

An interesting connection to the technique above is work presented in Davydychev & Delbourgo (1998) showing that 1-loop Feynman integrals of $N$ propagators can be related to the volumes of $N$-dimensional simplexes. Ortner & Wagner (1995) also relates this type of integral to simplexes. The connection to the technique above is the appearance there of the area of the triangle between the external and two internal momenta: in essence the 2-D volume of a simplex. It is not surprising that simplexes show up in these works: as we have shown here, expansion of the integrands into Gegenbauer series enforces “simplicial” constraints between the momentum magnitudes.

Our work here is also realted to that of Easther et al. (2000), which expands the Feynman propagator in a sum of generalized sinc functions. This expansion converges quickly, and also renders all of the integrals one needs to perform Gaussian, allowing most of the calculation to be done analytically just once. This is conceptually similar to Simonović et al. (2018), which also exploits the opportunity to perform integrals over a set of basis functions analytically.

Finally, regarding QFT, we note that our technique, as a method of handling simplexes, can also be used as an alternative representation of the time-ordering operator $\hat{T}$ that appears in Dyson series. In particular, our decoupling integral (reference equation) for the even powers, when not symmetrized, is identically zero when the frequency of the smaller-order sBF exceeds that of the larger-order sBF. This can be used to rewrite the time-ordered product that enters the integrals giving the Dyson series. Thus one obtains a purely algebraic representation of $\hat{T}$. Whether this has utility in performing QFT calculations with Dyson series, or is simply of conceptual interest, may be worth exploring in future work.

8 Conclusion

In this paper, we have shown that the integrals required for the 2-loop power spectrum in cosmological Eulerian Standard Perturbation Theory can be written as nested double 3-D convolution. The kernels one integrates against linear power spectra to generate the 2-loop corrections involve six different denominators each containing dummy momenta and possibly the external momentum. Four of these denominators can be automatically handled as convolutions, with the particular choice of four dictated by the arguments of the linear power spectra. The remaining two need to be decoupled so that they can be grouped into the convolvands entering each of the nested convolutions.

We accomplish this decoupling in four steps. First, we write each of the denominators as a Gegenbauer series; Gegenbauer polynomials in the dot product of the relevant unit vectors appear because the denominators all involve even powers (two and four). This factorizes the radial and angular dependences from each other. Second, we decouple the dot product in the Gegenabauers’ arguments by expanding each Gegenbauer into a sum of spherical harmonics depending on one unit vector each. Third, we decouple the radial piece, which is formally expanded in powers of the ratio of the greater to the lesser magnitude vector, by rewriting it as an overlap integral of two spherical Bessel functions. This integral automatically selects the correct (smaller) magnitude to be in the numerator of the ratio. However, to avoid paying the price of performing an integral at the end of our convolutions, which would require the convolutions themselves to be of higher dimension, we use an integral-to-sum identity to rewrite the overlap integral as an infinite series, though one that quickly converges. These steps give a complete decoupling and factorization of the 2-loop integrals into nested, double 3-D convolution. Since all of the angular dependence is known in closed form (as the power spectra are isotropic), we perform all of the angular integrations in each 3-D convolution to obtain 1-D radial integrals. We then have expressions for the 2-loop corrections that simply involve a sequence of 1-D integrations.

This approach is related to the long-used QFT multi-loop calculation aid Gegenbauer Polynomial $x$-Space Technique (GPXT), although this latter does not decouple the radial piece appearing in the Gegenbauer series. This work solves the same problem, but in a different way, as recent work by Simonović et al. (2018). That work expands each linear power spectrum as a series in complex power laws, and performs the integrals of the SPT kernels against these basis functions analytically, resulting in hypergeometric functions. The computational cost of that method is therefore front-loaded; initial calculation of the hypergeometric functions to high precision and with appropriate analytic continuations is challenging, but once this is done, only one FFT is required to obtain the complex power law coefficients of the linear power spectrum, and this vector can then be matrix-multiplied by the basis integrals.

In contrast, in the approach of the present work, no pre-computation is done, and one simply has a series of forward and inverse FFTs of the linear power spectra. The exact number of FFTs required will depend on the order to which both the Gegenbauer series and the series for the integral-to-sum identity used to handle the sBF overlap integral. The number of integrals could likely be somewhat reduced using recursion relations—taking a linear combination of two sbF integrals with a given power law weight could provide another sBF integral with a different power-law weight. More detailed comparison of the computational cost of the present method with that of Simonović et al. (2018), in particular the exact balance of less initial cost but more run-time calculation, is probably best done once the present method has been numerically implemented.

We have also shown that a few desirable extensions, such as to redshift-space and to polyspectra (bispectrum, trispectrum), should be straightforward. Further, we have presented several additional applications of the technique outlined here, such as its possible utility for simplifying multi-loop integrals in QFT and its potential to serve as a basis for a novel $N$-body integration scheme.

More generally, we note that the importance of factorization into radial and angular pieces as well as the appearance of sBFs and spherical harmonics is not surprising. Given the isotropy of the power spectrum, this factorization makes sense, and many other works exploit this framework as well. For instance, Dai et al. (2012) and Dai et al. (2013) work in a total angular momentum basis where they expand the observables using spherical harmonics and their generalizations motivated by the spherical symmetry in cosmology. Earlier work by Zaroubi & Hoffman (1993) displayed cosmological SPT to second-order in the basis of spherical harmonics.

Looking forward, we expect that techniques such as that outlined here will be useful in the context of large-volume, high-precision surveys such as DESI or LSST. DESI especially has BAO as one major focus, and BAO scales are large enough that perturbative methods are valid. However, DESI will achieve percent precision on of order ten redshift bins, for sub-percent precision overall, and so accurately including higher-order PT corrections will be worthwhile. Furthermore, given the high precision of DESI, it is likely that MCMC exploration of variations in the cosmological parameters will also be rewarded. This will require recomputation of the loop corrections a large number of times for different input linear power spectra. Thus the speed-up in loop calculations offered by the method presented here (or that of Simonović et al. 2018) will be particularly enabling.

Acknowledgements

ZS thanks Emanuele Castorina, Jolyon K. Bloomfield, Joanne D. Cohn, Alex Krolewski, Daniel J. Eisenstein, Douglas P. Finkbeiner, Yin Li, Stephen Portillo, Marcel Schmittfull, Marko Simonović, Masahiro Takada, and Martin White for useful conversations, as well as Cornelius Rampf and Marko Simonović for comments on a draft. ZS especially expresses gratitude to Bob Cahn for numerous discussions and encouragement over the long period of this work as well as comments on many calculations. ZS also thanks Pat McDonald and Jolyon Bloomfield for each sharing unpublished notes. ZS particularly thanks Donghui Jeong, Sarah Shandera, Matias Zaldarriaga, and Matt McQuinn for hospitality at respectively the Pennsylvania State University, the Institute for Advanced Study, and the University of Washington during portions of this work, as well as Joshua Silver, Roberta Cohen, Hillary Child, and Salman Habib for the same in Chicago and at Argonne National Laboratory. Support for portions of this work was provided by the National Aeronautics and Space Administration through Einstein Postdoctoral Fellowship Award Number PF7-180167 issued by the Chandra X-ray Observatory Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of the National Aeronautics Space Administration under contract NAS8-03060. ZS also received support during portions of this work from the Berkeley Center for Cosmological Physics as well as affiliate resources from Lawrence Berkeley National Laboratory.

Appendix A Relevant Mathematical Theorems and Identities

The spherical harmonic addition theorem decouples Legendre polynomials $\mathcal{L}_{\ell}$ whose argument is a dot product of two unit vectors into a sum of products of spherical harmonics $Y_{\ell m}$ each in just one unit vector, as (e.g. Arfken et al. 2013)

[TABLE]

The spherical harmonics we use are complex:

[TABLE]

where $\theta$ and $\phi$ are respectively the angle of $\hat{r}$ with the $z$-axis and the azimuthal angle of $\hat{r}$ with respect to the $x$-axis (i.e. we take the $x$-axis as our zero point for $\phi$). $\phi$ runs from [math] to $2\pi$ and $\theta$ from [math] to $\pi$. $P_{\ell}^{m}$ is an associated Legendre polynomial.

The conjugate identity for spherical harmonics changes the sign of $m$ and gives a pre-factor of $(-1)^{m}$, i.e.

[TABLE]

This identity follows from noting that the complex part of the spherical harmonic comes solely from $\exp[im\phi]$ and flipping $m$ is equivalent to conjugation. However then the associated Legendre polynomial part of the spherical harmonic changes from $P_{\ell}^{m}$ to $P_{\ell}^{-m}$ under a flip of $m$, necessitating the compensating phase factor $(-1)^{m}$ above.

Under parity $\hat{P}$, spherical harmonics behave as

[TABLE]

The plane wave expansion writes a complex exponential as an infinite sum over products of spherical harmonics and spherical Bessel functions (e.g. Arfken et al. 2013):

[TABLE]

It can be understood as a convenient “Rosetta stone” enabling translation between problems with translation invariance (homogeneity) and rotation invariance (isotropy). Cosmological settings frequently involve both, rendering this identity highly useful.

The Gaunt integral is defined

[TABLE]

with

[TABLE]

the $2\times 3$ matrices are Wigner 3j-symbols (e.g. Olver et al. 2010 §34.2, with basic properties in §34.3; more detailed discussion is in Varshalov (2013)).

The integral of four spherical harmonics can be obtained by first linearizing two spherical harmonics into a sum over single spherical harmonics using the Gaunt integral. We begin with

[TABLE]

where the coefficients $c_{LM}(l_{1},l_{2};m_{1},m_{2})$ are given by integrating both sides against $Y^{*}_{LM}(\hat{r})$ and invoking orthogonality, so that

[TABLE]

We then have

[TABLE]

by inserting equation (111) with the definition (112) into the first line above and then integrating. We note that there is no sum over $M$ because it is set by the zero-sum rule on the spins enforced by the Gaunt integrals: $m_{1}+m_{2}=M$. We further note that the sum over $L$ has compact support because of the triangle rules on total angular momenta: $|l_{1}-l_{2}|\leq L\leq l_{1}+l_{2}$ and the same constraint holds replacing $l_{1}\to l_{3}$ and $l_{2}\to l_{4}$. Finally, we chose to couple the first two angular momenta to $L$ and the last two to $L$, but of course choosing other combinations to recouple would lead to the same result.

The Gegenbauer polynomials can be written explicitly as

[TABLE]

e.g. San Kim et al. (2012) or Durand et al. (1976). The rising factorial is defined as $a^{(k)}=a(a+1)(a+2)\cdots(a+k-1)$.151515http://mathworld.wolfram.com/RisingFactorial.html.

Appendix B Factorization of Gegenbauer polynomials into spherical harmonics

Here we show that the Gegenbauer polynomials can be written as a sum of products of spherical harmonics each with one unit vector as their argument. The Gegenbauer polynomials have an explicit, finite representation in powers of their argument, as

[TABLE]

where the superscripts $\ell,\lambda$ identify the Gegenbauer polynomial being expanded. The second line comes from expanding $[(x-1)/2]^{k}$ in equation (118) using the binomial theorem and comparing with the first line above.

Powers can be represented using Legendre polynomials as

[TABLE]

where superscript $n$ identifies the power being expanded.161616http://mathworld.wolfram.com/LegendrePolynomial.html, equation (15), originally from Schmied (2005) personal communication to E. Weisstein. The double factorial has a different definition for even versus odd arguments; see e.g. equation (1) of Weisstein (2018).171717http://mathworld.wolfram.com/DoubleFactorial.html

Inserting equation (127) into equation (126) we find

[TABLE]

Summarizing, the Gegenbauer polynomials can be written as a finite series of powers and each power written as a finite series of Legendre polynomials, leading to a double sum. The Legendre expansion coefficients of a given Gegenbauer polynomial can then be obtained by holding the order in the second series fixed, so that we take a sum over all powers that contribute to a particular Legendre polynomial as it enters the Gegenbauer expansion. This approach can be conveniently understood as writing down a matrix showing the coupling between powers and Legendre polynomials in the series for the Gegenbauer polynomial, with the row being the power and the column being the Legendre. The way we constructed the representation was by staying in a fixed row (power) and reading across columns to obtain the Legendre coefficients for that power, but to obtain the desired Legendre coefficients we now instead read across rows at fixed column.

Our final step is to apply the spherical harmonic addition theorem (100) to write the Gegenbauer as

[TABLE]

if $x\equiv\hat{a}\cdot\hat{b}$. For more compact notation in the main text, we define $w_{j}^{\ell,\lambda}=4\pi v_{j}^{\ell,\lambda}/(2j+1)$, so that

[TABLE]

though we emphasize that $w_{j}^{\ell,\lambda}$ is spin-independent.

Appendix C Proof of decoupling integral for $|\vec{p}_{1}+\vec{p}_{2}|^{-1}$

Here we show that the radial part of the Legendre series can be represented in terms of the overlap integral of two spherical Bessel functions, i.e.

[TABLE]

We begin with the identity

[TABLE]

This identity is equation (46) of Slepian & Eisenstein (2017) and was proven there using two different ways to write the plane wave, expanding each side into spherical harmonics and spherical Bessel functions with the plane wave expansion (104), and integrating over solid angle. The $L=0$ case of this identity is the version of Gegenbauer’s addition theorem applying to spherical Bessel functions (see Sack 1964 equation 59). $C_{L_{1}L_{2}L}$ is defined in equation (110); the $2\times 3$ matrices are Wigner 3j-symbols.

Setting $L=0$ means that the Wigner 3j-symbols force $L_{1}=L_{2}$, and $L=0$ means $M=0$ so that $M_{2}=-M_{1}$. We find

[TABLE]

where we used that $Y_{00}=1/\sqrt{4\pi}$ and evaluated the zero-spin 3j-symbol as $1/\sqrt{2L_{1}+1}$ and the 3j-symbol with spins (with $M_{2}=-M_{1}$) as $(-1)^{M_{1}}/\sqrt{2L_{1}+1}$ using Olver et al. (2010) §34.3.1. We then used the conjugate identity to rewrite the second spherical harmonic from the previous simplification as unconjugated and cancel the $(-1)^{M_{1}}$ phase factor. Simplifying using that $\mathcal{C}_{L_{1}L_{1}0}=(2L_{1}+1)/\sqrt{4\pi}$, we obtain

[TABLE]

Canceling the $1/\sqrt{4\pi}$ and then resumming the right-hand side using the spherical harmonic addition theorem (100), we find

[TABLE]

We now integrate both sides with respect to $k$ as

[TABLE]

To obtain the first line above, we simply changed variables and used that the Sine integral has value $\pi/2$. We now observe that the middle expression above may be written as a Legendre series, so we find

[TABLE]

Integrating both sides against $[(2L+1)/2]\mathcal{L}_{L}(\mu_{12})$ with $\mu_{12}\equiv\hat{r}_{1}\cdot\hat{r}_{2}$, we invoke the orthogonality of the Legendre polynomials (the pre-factor is to compensate for the fact that they are not orthonormal) to find

[TABLE]

as desired. This result also appears in Bloomfield et al. (2017) equation (14), Mehrem (2009) equation (4.13), and in Watson (1944). An alternative derivation is to use the plane-wave expansion into sBFs and Legendre polynomials from the beginning, and then use the spherical harmonic addition theorem to integrate over the Legendres that will appear on the righthand side with arguments $\hat{k}\cdot\hat{r}_{1}$ and $\hat{k}\cdot\hat{r}_{2}$.

Appendix D Proof of decoupling integral for $|\vec{p}_{1}+\vec{p}_{2}|^{-2}$

We now prove that

[TABLE]

We may prove this second decoupling integral using the first, proven in the previous section, and the recursion relation for spherical Bessel functions. We first assume that $a>b$. We invoke the recursion relation (Olver et al. 2010 §10.51.2)181818https://dlmf.nist.gov/10.51

[TABLE]

where prime denotes $d/du$, to find that

[TABLE]

We set $u=ax$, so that $d/du=(1/x)\partial/\partial a$. Replacing $j_{\ell+1}(ax)$ in our desired integral (144) and simplifying, we find the equality

[TABLE]

where the last line followed by using our result (143) from the previous section with $a>b$ (which we assumed above) and simplifying.

We now address the case $a<b$. The first equality of equation (147) still holds, but we now use the $a<b$ result from the previous section to replace the matched-order integral, finding

[TABLE]

We note that our eigenfunction expansion, based on this integral, offers a factorized representation of a Heaviside step function centered at $a$. We have

[TABLE]

The above (save for the value at $b=a$, which is not $1/2$) is a Heaviside function $H$ shifted so that it switches on at $b=a$ rather than zero. We used the notation $\simeq$ to indicate the difference from the conventional definition of $H$ at $a=b$.

Appendix E Integral-to-sum identity for $|\vec{p}_{1}+\vec{p}_{2}|^{-2}$ decoupling integral

We now convert the spherical Bessel functions to Bessel functions (the relation is given below equation 17). We then use Dominici et al. (2012) Corollary 3.2, which is

[TABLE]

with $0<b<a<\pi$, ${\rm Re}(\mu+\nu-2k)>-1$ and $k$ a natural number. We need this identity to correctly reproduce the two integrals in our expression (22) both for $p_{2}<p_{1}$ and $p_{1}<p_{2}$ however. In particular we want the vanishing of one of the two integrals at all times to be correctly reproduced by the sums. However, the corollary above does have sufficient freedom to do this, because there is no constraint between the magnitudes of sBFs’ orders.

Consider our first integral in (22), which turns on when $p_{1}>p_{2}$ and is zero otherwise. We consider the two cases and show that the sum identity covers both. First, suppose $p_{1}>p_{2}$, so we set $a=p_{1}$, $b=p_{2}$, and $\mu=\ell+3/2$, $\nu=\ell+1/2$. $\mu+\nu=2\ell+2$ so to get the right power in the denominator (zero, as the factor of $x$ in the integrand is canceled when we convert to Bessel functions, each of which brings $1/\sqrt{x}$), we set $k=\ell+1$.

We then have

[TABLE]

We now consider the case where $p_{1}<p_{2}$ (the integral vanishes here). In equation (150) we now set $a=p_{2}$ and $b=p_{1}$; but we can also set $\nu=\ell+3/2$ and $\mu=\ell+1/2$. We then have

[TABLE]

The same approach shows that the second integral in equation (22) is also faithfully represented by the sum on both regions; the proof is exactly as above but replacing everywhere $p_{1}$ with $p_{2}$ and $p_{2}$ with $p_{1}$, so that $p_{2}$ is paired with the larger-index Bessel function.

Finally, we note that Dominici et al. (2012) give an alternative proof of Corollary 3.2 in §4 of their paper (suggested to them by T. Koornwinder) that uses Parseval’s relation for the Fourier transform and for Fourier series to show that the analog of the corollary where the integral runs from $-\infty$ to $\infty$ and the sum from $-\infty$ to $\infty$ as well holds for $0<a,b<\pi$, i.e. with no restriction on $b$ relative to $a$ (see Dominici et al. 2012 equation 4.2). Since our particular case has an even integrand (the sum of the Bessel function indices is even) and an even summand, one can divide both sides by two. Properly treating the case where $n=0$ using the definition of $\epsilon_{n}$ as $1/2$ at $n=0$ in Dominici et al. (2012), this argument shows that the $a<b$ restriction is not relevant for our case.

Appendix F Inverse Fourier Transform of $1/|\vec{q}|^{n},\;N=2,4$ and regularization

Here we obtain the two inverse FTs required to perform the radial (over momentum magnitudes) convolution integrals as products in configuration space. See for instance equation (55). The $1/q^{2}$ case is simple; the $1/q^{4}$ case is formally divergent, but the multiplication by other functions, in particular the correlation function, in the subsequent configuration space radial integrals will effectively regularize this divergence. Thus we will compute the inverse transform of a softened version of this kernel and take the softening to zero after the configuration-space radial integrals have been performed.

For the first kernel, we have

[TABLE]

where the 3D inverse FT reduced to a 1D $j_{0}$ transform due to the spherical symmetry of the kernel. We used the change of variables $u=qr$ and then that the sine integral $\int du\;\sin u/u=\pi/2$.

For the second kernel, we have

[TABLE]

We note that we chose to regularize the divergence by softening the kernel here; we could also have chosen to cut off the integral above zero. We present the result of that as well, on the view that one or the other may prove more suitable for numerical work.

To evaluate the integral (154), we notice that $(1/2)\partial^{2}/\partial\epsilon^{2}$ of our previous integral (153) will generate the integral (154). To perform this shifted-denominator version of our previous integral, we change variables so that $u=q+\epsilon$. We then in the numerator have $\sin[(u-\epsilon)r](u-\epsilon)^{2}$ and can expand this using respectively the sum-to-product identity for trigonometric functions and the binomial theorem.

We find

[TABLE]

where ${\rm Ci}$ is the cosine integral and ${\rm Si}$ the sine integral. Performing the parametric differentiation on equation (155), we have

[TABLE]

As noted above, a second regularization would have been just to truncate the integral at $\epsilon$ rather than integrate all the way down to zero:

[TABLE]

We obtained this result using integration by parts.