Bessel integrals in epsilon expansion: Squared spherical Bessel functions averaged with Gaussian power-law distributions

Bessel integrals of type ∫0∞kμ+2e-ak2-(b+iω)kjl2(pk)dk are studied, where the squared spherical Bessel function jl2 is averaged with a modulated Gaussian power-law density. These integrals define the multipole moments of Gaussian random fields on the unit sphere, arising in multipole fits of temperature and polarization power spectra of the cosmic microwave background. The averages can be calculated in closed form as finite Hankel series, which allow high-precision evaluation. In the case of integer power-law exponents μ, singularities emerge in the series coefficients, which requires ε expansion. The pole extraction and regularization of singular Hankel series is performed, for integer Gaussian power-law densities as well as for the special case of Kummer averages (a=0 in the exponential of the integrand). The singular ε residuals are used to derive combinatorial identities (sum rules) for the rational Hankel coefficients, which serve as consistency checks in precision calculations of the integrals. Numerical examples are given, and the Hankel evaluation of Gaussian and Kummer averages is compared with their high-index Airy approximation over a wide range of integer Bessel indices l.