BLAST: Beyond Limber Angular power Spectra Toolkit. A fast and efficient algorithm for 3x2 pt analysis

Technical Summary: Blast.jl – A Fast and Efficient Algorithm for 3 × 2 pt Analysis

Problem Statement
The advent of next-generation photometric and spectroscopic surveys (e.g., Euclid, LSST, Roman, DESI) promises high-precision measurements of Large Scale Structure (LSS). These surveys rely heavily on "3 × 2 pt" statistics, which combine galaxy clustering, weak lensing shear-shear correlations, and galaxy-galaxy lensing cross-correlations. The theoretical modeling of these statistics requires calculating angular power spectra ($C_\ell$) by integrating the product of the 3D matter power spectrum and highly oscillating spherical Bessel functions over three dimensions.

The standard approach to reduce computational cost is the Limber approximation, which simplifies the 3D integral to a 1D integral by assuming modes at different epochs do not contribute to the line-of-sight integration. However, this approximation introduces significant errors, particularly on large angular scales (low $\ell$). To fully exploit the constraining power of future surveys and control systematics, accurate non-Limber calculations are essential. Yet, performing these full 3D integrals is computationally prohibitive for parameter inference methods like Markov Chain Monte Carlo (MCMC), which require millions of evaluations. For instance, recent analyses using the Limber approximation required 0.5 million CPU hours to converge.

Methodology
The authors present Blast.jl, a novel algorithm designed to compute full non-Limber angular power spectra efficiently without assuming a scale-dependent growth rate. The core methodology relies on the decomposition of the matter power spectrum using Chebyshev polynomials.

1. Chebyshev Decomposition: The algorithm approximates the smooth matter power spectrum $P(k, \chi_1, \chi_2)$ on the interval $[k_{min}, k_{max}]$ using a Chebyshev series expansion. This allows the separation of the integral into geometric components and cosmological components.
2. Pre-computation of Geometric Terms: The integral is restructured such that the most computationally expensive part—the integration of the product of two spherical Bessel functions and Chebyshev polynomials—becomes independent of cosmological parameters. These terms ($\tilde{T}^{AB}_{n; \ell}$) are computed once and stored.
3. Efficient Evaluation: For any new cosmological model, the algorithm only needs to recalculate the Chebyshev coefficients of the power spectrum, a process that is computationally inexpensive. The remaining outer integrals over comoving distance ($\chi$) and the ratio $R = \chi_2/\chi_1$ are performed using standard quadrature rules (Simpson's rule for $\chi$ and Clenshaw-Curtis quadrature for $R$).
4. Coordinate Transformation: A change of variables from $(\chi_1, \chi_2)$ to $(\chi, R)$ is employed to better sample the region where $\chi_1 \approx \chi_2$, which contributes most significantly to the integral.
5. Non-linear Correction: The method treats the linear and non-linear power spectra separately. The non-linear contribution, which is significant only on small scales where the Limber approximation is valid, is corrected using a Limber-based approach, while the full non-Limber integral is applied to the linear component.

Key Contributions

• Algorithmic Innovation: Blast.jl introduces a decomposition-based approach using Chebyshev polynomials that avoids the logarithmic grid spacing constraints of FFTLog-based methods (like FKEM and matter) and the iterative differential equation solving of the Levin method.
• Pre-computation Strategy: By isolating the cosmology-independent geometric integrals, the algorithm drastically reduces the time required for repeated evaluations during parameter inference.
• Flexibility: Unlike some competing methods, Blast.jl does not assume a scale-independent growth rate, making it applicable to cosmologies with massive neutrinos or modified gravity.
• Open Source: The code is publicly available on GitHub, accompanied by a repository detailing its usage.

Results
The algorithm was evaluated against the "N5K" (Non-local No-Nonsense Non-Limber Numerical Knock-out) challenge, organized by the LSST Dark Energy Science Collaboration to benchmark non-Limber integration methods. The tests were conducted on a high-performance computing node (Narval) with 64 threads, comparing Blast.jl against the challenge winners: FKEM (CosmoLike), matter, and Levin.

• Speed: At the fiducial accuracy requirement ($\Delta\chi^2 < 0.2$ for $\ell < 200$), Blast.jl achieved a runtime of 0.00512 seconds using 64 cores. This is approximately 10 times faster than FKEM (the previous fastest method, 0.05871 s) and roughly 15 times faster when optimized to 32 cores (0.00360 s).
• Scaling: Blast.jl demonstrates excellent scaling with respect to the number of tomographic redshift bins. While other methods show runtime increases with the number of bins, Blast.jl's runtime remains effectively independent of the number of bins (ranging from 2 to 15).
• Accuracy: The method meets the challenge's accuracy requirements. The authors note that accuracy can be tuned by adjusting the number of Chebyshev polynomials and sampling points without significantly impacting runtime.
• Precision vs. Speed Trade-off: Reducing the floating-point precision from Float64 to Float32 improved runtime by a factor of 2 with negligible impact on accuracy ($\Delta\chi^2$ increased from 0.1054 to 0.1064).

Significance
The paper claims that Blast.jl represents the fastest known non-Limber integration algorithm currently available, outperforming the state-of-the-art methods by over an order of magnitude in the fiducial scenario. Its primary significance lies in making cosmological analyses using 3 × 2 pt statistics computationally feasible for next-generation surveys without relying on the error-prone Limber approximation.

The authors emphasize that this speed is critical for enabling MCMC-based parameter inference with the massive datasets expected from future surveys. Furthermore, the algorithm's design allows for future extensions, such as the inclusion of redshift space distortions, magnification bias, and relativistic corrections, as well as the integration of CMB lensing. The authors also highlight plans to incorporate automatic differentiation, which would enable the use of gradient-based inference methods, further accelerating cosmological analysis pipelines.