Bessel integrals of type integral(infinity)(0)k(mu+2)e(-ak2) (-(b+i omega)k)j(l)(2)(pk)dk are studied, where the squared spherical Bessel function j(l)(2) 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 mu, singularities emerge in the series coefficients, which requires epsilon 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 epsilon 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. (C) 2013 Elsevier Inc. All rights reserved.