Blast.jl: Differentiable Non-Limber Power Spectra for Joint Clustering, Shear, and CMB lensing Analyses

Upcoming large-scale structure surveys require accurate theoretical predictions for angular power spectra on the largest angular scales, where the commonly used Limber approximation breaks down and multiple relativistic and observational effects become relevant. At the same time, modern cosmological analyses increasingly rely on gradient-based inference techniques, motivating the development of fast, fully differentiable algorithms. In this work, we present a major extension of Blast.jl, a numerical toolkit for the efficient computation of non-Limber angular power spectra. Building on the Chebyshev-polynomial decomposition of the original framework, the updated algorithm incorporates redshift-space distortions, magnification bias, primordial non-Gaussianity, intrinsic alignments, CMB lensing, and the integrated Sachs-Wolfe effect, while remaining fully differentiable through custom automatic differentiation rules. Despite the increased physical complexity, the algorithm retains favorable scaling by precomputing all cosmology-independent quantities. We validate the expanded framework against brute-force integration and established cosmological codes, and provide a detailed analysis of the trade-off between computational speed and precision as a function of the algorithm hyperparameters. We demonstrate the readiness of Blast.jl for gradient-based inference through a simulated likelihood analysis with LSST Y10-like mock data, recovering unbiased cosmological and nuisance parameter constraints from a 35-parameter model using gradient-based samplers.