Lawrence Berkeley National Laboratory

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.