In Section 2, we derive a finite Hankel series representation of squared spherical Bessel functions, which is used to obtain closed analytic expressions for the above integrals suitable for high-precision calculations. The finite Hankel expansion admits term-by-term integration in terms of confluent hypergeometric functions. At integer power-law exponents μ$\mu$, the confluent functions become singular, so that the Hankel series requires epsilon expansion, see Section 3, where we perform the pole separation of the confluent functions appearing in the coefficients of the Hankel series. In Section 4, we obtain explicit formulas for the regularized Hankel series in closed form, arriving at a numerically efficient finite series representation of the Bessel integrals in terms of confluent hypergeometric functions.

In Section 5, we discuss a special case, where the averaging of the squared spherical Bessel function is carried out with an exponentially cut and modulated power-law density (a=0$a=0$ in the above integrand). In this case, the confluent functions in the Hankel coefficients are still singular at integer power-law exponents μ$\mu$ but become elementary, so that the regularized Hankel series are elementary functions as well. The Conclusion, Section 6, gives an overview of the results. In Appendix A, we use the ε$\epsilon$ pole extraction to derive combinatorial identities (sum rules) for the Hankel coefficients, by making use of the fact that the singular ε$\epsilon$ residuals of the series coefficients cancel one another if the integrals converge.

2. Integrals of squared spherical Bessel functions expressed as finite Hankel series

2.1. Bessel integrals of type${\int }_0^{\infty }{k}^{\mu +2}{\mathrm{e}}^{-{\mathit{ak}}^2-(b+\mathrm{i}\omega )k}{j}_l^2(\mathit{pk})\mathrm{d}k$

We study the Bessel average ${\int }_0^{\infty }g(k)j_l^2(\mathit{pk})k^2\mathrm{d}k$, where the distribution g(k)=k^μe^{-ak2-(b+iω)k}$g(k)=k^{\mu }{\mathrm{e}}^{-{\mathit{ak}}^2-(b+\mathrm{i}\omega )k}$ is a modulated Gaussian power law with real exponents μ,a>0,b$\mu \text{,}a>0\text{,}b$ and ω$\omega$ (or a=0$a=0$ and b>0)$b>0)$. The scale parameter p in the argument of the squared spherical Bessel (ssB) function $j_l^2(\mathit{pk})$ is positive. Spherical Bessel functions are rescaled ordinary Bessel functions of half-integer order, $j_l(x)=\sqrt{\pi /(2x)}{J}_{l+1/2}(x)$, the index l being a non-negative integer [3] and [4]. These integrals converge for μ+2+2l>-1$\mu +2+2l>-1$, as the ascending series of $j_l^2(x)$ starts with the power x^2l$x^{2l}$.

We will need an integral representation of a parabolic cylinder function [5],

equation(2.1)

$D_{\exp }(p\text{;}\mu \text{,}a\text{,}b\text{,}\omega )={\int }_0^{\infty }{k}^{\mu }{\mathrm{e}}^{-{\mathit{ak}}^2-(b+\mathrm{i}\omega )k}\text{exp}(2\mathrm{i}\mathit{pk})\mathrm{d}k=\frac{\mathrm{\Gamma }(\mu +1)}{{(2a)}^{(\mu +1)/2}}{\mathrm{e}}^{-y^2/4}{D}_{-(\mu +1)}(\mathrm{i}y)\text{,}$

Turn MathJax on

equation(2.2)

$y=\frac{\omega -2p-\mathrm{i}b}{\sqrt{2a}}\text{.}$

Turn MathJax on

equation(2.3)

$D_{\exp }(p\text{;}\mu \text{,}a\text{,}b\text{,}\omega )=\frac{\mathrm{\Gamma }(\mu +1)}{2^{\mu +1}{a}^{(\mu +1)/2}}U\left(\frac{\mu +1}{2}\text{,}\frac{1}{2}\text{,}\frac{{(\mathrm{i}y)}^2}{2}\right)\text{,}$

Turn MathJax on

equation(2.4)

$U\left(\frac{\mu +1}{2}\text{,}\frac{1}{2}\text{,}\frac{{(\mathrm{i}y)}^2}{2}\right)={\mathrm{e}}^{-y^2/2}\left[\frac{\sqrt{\pi }}{\mathrm{\Gamma }(\mu /2+1)}{}_1{F}_1\left(\frac{-\mu }{2}\text{,}\frac{1}{2}\text{,}\frac{y^2}{2}\right)-\frac{\sqrt{2\pi }\mathrm{i}y}{\mathrm{\Gamma }((\mu +1)/2)}{}_1{F}_1\left(\frac{1-\mu }{2}\text{,}\frac{3}{2}\text{,}\frac{y^2}{2}\right)\right]\text{,}$

Turn MathJax on

equation(2.5)

$D_{\exp }(p\text{;}\mu \text{,}a\text{,}b\text{,}\omega )=\frac{1}{2a^{(\mu +1)/2}}{\mathrm{e}}^{-y^2/2}\left[\mathrm{\Gamma }\left(\frac{\mu +1}{2}\right){}_1{F}_1\left(\frac{-\mu }{2}\text{,}\frac{1}{2}\text{,}\frac{y^2}{2}\right)-\mathrm{\Gamma }\left(\frac{\mu +2}{2}\right)\sqrt{2}\mathrm{i}y{}_1{F}_1\left(\frac{1-\mu }{2}\text{,}\frac{3}{2}\text{,}\frac{y^2}{2}\right)\right]=\frac{1}{2a^{(\mu +1)/2}}\left[\mathrm{\Gamma }\left(\frac{\mu +1}{2}\right){}_1{F}_1\left(\frac{1+\mu }{2}\text{,}\frac{1}{2}\text{,}-\frac{y^2}{2}\right)-\mathrm{\Gamma }\left(\frac{\mu +2}{2}\right)\sqrt{2}\mathrm{i}y{}_1{F}_1\left(1+\frac{\mu }{2}\text{,}\frac{3}{2}\text{,}-\frac{y^2}{2}\right)\right]\text{.}$

Turn MathJax on

We split the exponential exp(2ipk)$\text{exp}(2\mathrm{i}\mathit{pk})$ in integral (2.1) into real and imaginary parts, defining

equation(2.6)

$D_{\cos }(p\text{;}\mu \text{,}a\text{,}b\text{,}\omega )={\int }_0^{\infty }{k}^{\mu }{\mathrm{e}}^{-{\mathit{ak}}^2-(b+\mathrm{i}\omega )k}\cos (2\mathit{pk})\mathrm{d}k=\frac{1}{2}(D_{\exp }(p\text{;}\mu \text{,}a\text{,}b\text{,}\omega )+D_{\exp }(-p\text{;}\mu \text{,}a\text{,}b\text{,}\omega ))\text{,}$

Turn MathJax on

equation(2.7)

$D_{\sin }(p\text{;}\mu \text{,}a\text{,}b\text{,}\omega )={\int }_0^{\infty }{k}^{\mu }{\mathrm{e}}^{-{\mathit{ak}}^2-(b+\mathrm{i}\omega )k}\sin (2\mathit{pk})\mathrm{d}k=\frac{1}{2\mathrm{i}}(D_{\exp }(p\text{;}\mu \text{,}a\text{,}b\text{,}\omega )-D_{\exp }(-p\text{;}\mu \text{,}a\text{,}b\text{,}\omega ))\text{,}$

Turn MathJax on

equation(2.8)

$D_0(\mu \text{,}a\text{,}b\text{,}\omega )={\int }_0^{\infty }{k}^{\mu }{\mathrm{e}}^{-{\mathit{ak}}^2-(b+\mathrm{i}\omega )k}\mathrm{d}k=D_{\exp }(0\text{;}\mu \text{,}a\text{,}b\text{,}\omega )\text{.}$

Turn MathJax on

We will employ a finite Hankel representation of the squared spherical Bessel function,

equation(2.9)

$j_l^2(x)=\frac{1}{2x^2}\left({(-1)}^{l+1}{\hat{A}}_l(x)\cos 2x+{(-1)}^l{\hat{B}}_l(x)\sin 2x+{\hat{C}}_l(x)\right)\text{,}$

Turn MathJax on

${\hat{A}}_l(x)\text{,}{\hat{B}}_l(x)$

${\hat{C}}_l(x)$

equation(2.10)

$[l\text{,}k]=\frac{1}{\mathrm{\Gamma }(1+k)}\frac{\mathrm{\Gamma }(l+1+k)}{\mathrm{\Gamma }(l+1-k)}\text{.}$

Turn MathJax on

${\hat{A}}_l(x)$

equation(2.11)

${\hat{A}}_l(x)={\displaystyle \sum _{k=0}^l}\frac{a_{2k}(l)}{x^{2k}}\text{,}$

Turn MathJax on

equation(2.12)

$\begin{array}{cc}\hfill & a_{2k⩽l}(l)=-\frac{{(-1)}^k}{2^{2k}}{[l\text{,}k]}^2+\frac{{(-1)}^k}{2^{2k-1}}{\displaystyle \sum _{m=0}^k}[l\text{,}m][l\text{,}2k-m]\text{,}\hfill \\ \hfill & a_{2k⩾l}(l)=-\frac{{(-1)}^k}{2^{2k}}{[l\text{,}k]}^2+\frac{{(-1)}^k}{2^{2k-1}}{\displaystyle \sum _{m=0}^{l-k}}[l\text{,}2k-l+m][l\text{,}l-m]\text{,}\hfill \end{array}$

Turn MathJax on

As for the polynomials ${\hat{B}}_l(x)$ in (2.9), we define ${\hat{B}}_0(x)=0$ and

equation(2.13)

${\hat{B}}_l(x)={\displaystyle \sum _{k=0}^{l-1}}\frac{b_{2k+1}(l)}{x^{2k+1}}\text{,}$

Turn MathJax on

equation(2.14)

$\begin{array}{cc}\hfill & b_{2k+1⩽l}(l)=\frac{{(-1)}^k}{2^{2k}}{\displaystyle \sum _{m=0}^k}[l\text{,}m][l\text{,}2k+1-m]\text{,}\hfill \\ \hfill & b_{2k+1⩾l}(l)=\frac{{(-1)}^k}{2^{2k}}{\displaystyle \sum _{m=0}^{l-k-1}}[l\text{,}2k+1-l+m][l\text{,}l-m]\text{,}\hfill \end{array}$

Turn MathJax on

${\hat{C}}_l(x)$

equation(2.15)

${\hat{C}}_l(x)={\displaystyle \sum _{k=0}^l}\frac{c_{2k}(l)}{x^{2k}}\text{,}{c}_{2k}(l)=\frac{[l\text{,}k]}{2^{2k}}\frac{\mathrm{\Gamma }(2k+1)}{\mathrm{\Gamma }(k+1)}\text{,}$

Turn MathJax on

We substitute the Hankel representation (2.9) of the squared Bessel function into integral

equation(2.16)

$D_{\text{ssB}}(l\text{,}p\text{;}\mu \text{,}a\text{,}b\text{,}\omega )={\int }_0^{\infty }{k}^{\mu +2}{\mathrm{e}}^{-{\mathit{ak}}^2-(b+\mathrm{i}\omega )k}{j}_l^2(\mathit{pk})\mathrm{d}k$

Turn MathJax on

${\hat{A}}_l(x)\text{,}{\hat{B}}_l(x)$

${\hat{C}}_l(x)$

equation(2.17)

$D_{\text{ssB}}(l\text{,}p\text{;}\mu \text{,}a\text{,}b\text{,}\omega )=\frac{1}{2p^2}\left({(-1)}^{l+1}{D}_{\text{ssB}}^{(1)}+{(-1)}^l{D}_{\text{ssB}}^{(2)}+D_{\text{ssB}}^{(3)}\right)\text{,}$

Turn MathJax on

$D_{\text{ssB}}^{(1)}(l\text{,}p\text{;}\mu \text{,}a\text{,}b\text{,}\omega )={\displaystyle \sum _{k=0}^l}\frac{a_{2k}(l)}{p^{2k}}{D}_{\cos }(p\text{;}\mu -2k\text{,}a\text{,}b\text{,}\omega )\text{,}$

Turn MathJax on

$D_{\text{ssB}}^{(2)}(l\text{,}p\text{;}\mu \text{,}a\text{,}b\text{,}\omega )={\displaystyle \sum _{k=0}^{l-1}}\frac{b_{2k+1}(l)}{p^{2k+1}}{D}_{\sin }(p\text{;}\mu -2k-1\text{,}a\text{,}b\text{,}\omega )\text{,}$

Turn MathJax on

equation(2.18)

$D_{\text{ssB}}^{(3)}(l\text{,}p\text{;}\mu \text{,}a\text{,}b\text{,}\omega )={\displaystyle \sum _{k=0}^l}\frac{c_{2k}(l)}{p^{2k}}{D}_0(\mu -2k\text{,}a\text{,}b\text{,}\omega )\text{.}$

Turn MathJax on

$D_{\text{ssB}}^{(2)}(l=0)=0$

If a=0$a=0$ and b>0$b>0$ in integral (2.16), then integral D_exp$D_{\exp }$ in (2.1) is elementary, and the confluent hypergeometric D_exp(p;μ,a,b,ω)$D_{\exp }(p\text{;}\mu \text{,}a\text{,}b\text{,}\omega )$ in (2.3) and (2.5) is replaced by

equation(2.19)

$D_{\exp }(p\text{;}\mu \text{,}0\text{,}b\text{,}\omega )=\frac{\mathrm{\Gamma }(\mu +1)}{{(b+\mathrm{i}(\omega -2p))}^{\mu +1}}\text{.}$

Turn MathJax on

equation(2.20)

$D_{\exp }(p\text{;}\mu \text{,}a\text{,}b\text{,}\omega )\sim \frac{\mathrm{\Gamma }((\mu +1)/2)}{2a^{(\mu +1)/2}}\text{,}$

Turn MathJax on

2.2. A limit case: Kummer averages${\int }_0^{\infty }{k}^{\mu +2}{\mathrm{e}}^{-(b+\mathrm{i}\omega )k}{j}_l^2(\mathit{pk})\mathrm{d}k$

We consider the modulated and exponentially cut power law k^μe^-(b+iω)k$k^{\mu }{\mathrm{e}}^{-(b+\mathrm{i}\omega )k}$, with b>0$b>0$ and real ω$\omega$. The power-law index μ$\mu$ is real and non-integer. (Integer indices will be studied in Section 5.) This is the limit case a=0$a=0$ in (2.16) and (2.19). We use the shortcut β=b+iω,Reβ>0$\beta =b+\mathrm{i}\omega \text{,Re}\beta >0$, so that the Bessel integral (2.16) reads

equation(2.21)

$D_{\text{ssB}}(l\text{,}p\text{;}\mu \text{,}0\text{,}b\text{,}\omega )={\int }_0^{\infty }{k}^{\mu +2}{\mathrm{e}}^{-\beta k}{j}_l^2(\mathit{pk})\mathrm{d}k\text{.}$

Turn MathJax on

equation(2.22)

$D_{\exp }(p\text{;}\mu \text{,}0\text{,}b\text{,}\omega )={\int }_0^{\infty }{k}^{\mu }{\mathrm{e}}^{-\beta k}\text{exp}(2\text{i}\mathit{pk})\mathrm{d}k=\frac{\mathrm{\Gamma }(\mu +1)}{{(\beta -2\mathrm{i}p)}^{\mu +1}}\text{.}$

Turn MathJax on

equation(2.23)

$D_{\cos }(p\text{;}\mu \text{,}0\text{,}b\text{,}\omega )={\int }_0^{\infty }{k}^{\mu }{\mathrm{e}}^{-\beta k}\cos (2\mathit{pk})\mathrm{d}k=\frac{1}{2}\mathrm{\Gamma }(\mu +1)\left({(\beta -2\mathrm{i}p)}^{-\mu -1}+{(\beta +2\mathrm{i}p)}^{-\mu -1}\right)\text{,}$

Turn MathJax on

equation(2.24)

$D_{\sin }(p\text{;}\mu \text{,}0\text{,}b\text{,}\omega )={\int }_0^{\infty }{k}^{\mu }{\mathrm{e}}^{-\beta k}\sin (2\mathit{pk})\mathrm{d}k=\frac{1}{2\mathrm{i}}\mathrm{\Gamma }(\mu +1)\left({(\beta -2\mathrm{i}p)}^{-\mu -1}-{(\beta +2\mathrm{i}p)}^{-\mu -1}\right)$

Turn MathJax on

equation(2.25)

$D_0(\mu \text{,}0\text{,}b\text{,}\omega )={\int }_0^{\infty }{k}^{\mu }{\mathrm{e}}^{-\beta k}\mathrm{d}k=\mathrm{\Gamma }(\mu +1){\beta }^{-\mu -1}\text{.}$

Turn MathJax on

2.3. Laplace/Fourier asymptotics of integral ${\int }_0^{\infty }{k}^{\mu +2}{\mathrm{e}}^{-{\mathit{ak}}^2-(b+\mathrm{i}\omega )k}{j}_l^2(\mathit{pk})\mathrm{d}k$

We calculate integral D_ssB(l,p;μ,a,b,ω)$D_{\text{ssB}}(l\text{,}p\text{;}\mu \text{,}a\text{,}b\text{,}\omega )$ in (2.16) by making use of the ascending series expansion of the squared spherical Bessel function $j_l^2(x)$[6],

equation(2.26)

$\begin{array}{cc}\hfill & j_l^2(x)={\displaystyle \sum _{k=0}^{\infty }}{\alpha }_{2k}(l)x^{2k+2l}\text{,}\hfill \\ \hfill & {\alpha }_{2k}(l)=\frac{{(-1)}^k}{\mathrm{\Gamma }(k+1)}\frac{2^{2l+2k}{\mathrm{\Gamma }}^2(l+k+1)}{\mathrm{\Gamma }(2l+k+2)\mathrm{\Gamma }(2l+2k+2)}\text{.}\hfill \end{array}$

Turn MathJax on

equation(2.27)

${\alpha }_{2k}(l)=\frac{{(-1)}^k}{4\sqrt{\pi }}\frac{{\mathrm{e}}^{-2k(\log k-1)}}{k^{2l+5/2}}\left(1+\mathrm{O}\left(\frac{1}{k}\right)\right)\text{,}$

Turn MathJax on

equation(2.28)

$\mathrm{\Gamma }(n+\alpha )=\sqrt{2\pi }\frac{{\mathrm{e}}^{n(\log n-1)}}{n^{1/2-\alpha }}\left(1+\mathrm{O}\left(\frac{1}{n}\right)\right)\text{.}$

Turn MathJax on

equation(2.29)

$D_{\text{ssB}}(l\text{,}p\text{;}\mu \text{,}a\text{,}b\text{,}\omega )={\displaystyle \sum _{k=0}^{\infty }}{\alpha }_{2k}(l)p^{2k+2l}{D}_0(\mu +2k+2l+2\text{,}a\text{,}b\text{,}\omega )\text{,}$

Turn MathJax on

$y=(\omega -\mathrm{i}b)/\sqrt{2a}$

In series (2.29), we consider the limit a=0$a=0$ and b>0$b>0$, so that D₀(μ,0,b,ω)=Γ(μ+1)/(b+iω)^μ+1$D_0(\mu \text{,}0\text{,}b\text{,}\omega )=\mathrm{\Gamma }(\mu +1)/{(b+\mathrm{i}\omega )}^{\mu +1}$, cf. (2.19). Series (2.29) is thus majorized by a polylogarithm ${\sum }_{k=1}^{\infty }{k}^{\mu }{z}^k$, with z=4p²/(b+iω)²$z=4p^2/{(b+\mathrm{i}\omega )}^2$, which converges in the unit disk |z|<1$\left|z\right|<1$ for arbitrary μ$\mu$[7]. This can be seen from the asymptotic series coefficients (2.27) and the limit (2.28) of the gamma function. Rapid convergence apparently requires a large parameter b (Laplace asymptotics) or ω$\omega$ (Fourier asymptotics).

In the opposite limit a→∞$a\to \infty$, we can use (2.20) as estimate for D₀(μ,a,b,ω)$D_0(\mu \text{,}a\text{,}b\text{,}\omega )$ in series (2.29). In this case, series (2.29) is convergent due to a factor e^-k(logk-1)${\mathrm{e}}^{-k(\log k-1)}$ emerging in the series coefficients for large k, cf. (2.27) and (2.28). However, this factor only emerges for k≫l$k\gg l$ (since the asymptotic limit in (2.27) only applies in this regime), which requires to sum a large number of terms in series (2.29) if the Bessel index l is large. Rapid convergence requires a large exponent a or a large |b+iω|$\left|b+\mathrm{i}\omega \right|$ in integral (2.16). This expansion also works if both a and |b+iω|$\left|b+\mathrm{i}\omega \right|$ are large and y is moderate, cf. (2.1) and (2.2). In brief, if the exponents a and/or |b+iω|$\left|b+\mathrm{i}\omega \right|$ in the integrand (2.16) are large numbers and the index l is low or moderate, the Laplace/Fourier series expansion (2.29) is efficient. The parameter p in series (2.18) and (2.29) can be put equal to one, as it can be scaled into the exponents a,b$a\text{,}b$ and ω$\omega$ of the integrand in (2.16).

3. Bessel integrals with integer power-law exponent: epsilon expansion and Hermite residuals

In Sections 2.1 and 2.2, we considered non-integer power-law exponents μ$\mu$ in integral ${\int }_0^{\infty }{k}^{\mu +2}{\mathrm{e}}^{-{\mathit{ak}}^2-(b+\mathrm{i}\omega )k}{j}_l^2(\mathit{pk})\mathrm{d}k$, cf. (2.16), where $j_l^2(\mathit{pk})$ is a squared spherical Bessel function of integer index l⩾0$l⩾0$, and p a positive scale parameter. Here and in Sections 4 and 5, we study integer power-law exponents μ=m$\mu =m$, so that $j_l^2(\mathit{pk})$ is averaged with the distribution k^me^{-ak2-(b+iω)k}$k^m{\mathrm{e}}^{-{\mathit{ak}}^2-(b+\mathrm{i}\omega )k}$, with a>0$a>0$, real b and ω$\omega$. The k²$k^2$ factor in the integrand typically stems from a 3D volume integration in polar coordinates [2]. The exponent m can be negative, provided that m+2⩾-2l$m+2⩾-2l$, so that the integral is convergent. We will calculate integral (2.16) by performing ε$\epsilon$ expansion for arbitrary integer power-law exponents μ=m⩾-2l-2$\mu =m⩾-2l-2$.

The analytic μ$\mu$ continuations D_cos,sin(p;μ-n,a,b,ω)$D_{\cos \text{,}\sin }(p\text{;}\mu -n\text{,}a\text{,}b\text{,}\omega )$ and D₀(μ-n,a,b,ω)$D_0(\mu -n\text{,}a\text{,}b\text{,}\omega )$ of the integrals in (2.6), (2.7) and (2.8), which are required in the Hankel series (2.18), become singular due to poles in the gamma function in (2.3). We put μ=m+ε$\mu =m+\epsilon$, where m is an arbitrary integer, and n=m+1+j$n=m+1+j$, so that 1+μ-n=-j+ε$1+\mu -n=-j+\epsilon$. Poles in D_exp(p;μ-n,a,b,ω)$D_{\exp }(p\text{;}\mu -n\text{,}a\text{,}b\text{,}\omega )$ (defined by (2.3) or (2.5)) and in its trigonometric components D_cos,sin,0$D_{\cos \text{,}\sin \text{,}0}$ in (2.6), (2.7) and (2.8) emerge only at integer μ-n=-j-1$\mu -n=-j-1$ with j⩾0$j⩾0$.

We start with the ε$\epsilon$ expansion of integral D_exp(p;μ-n,a,b,ω)$D_{\exp }(p\text{;}\mu -n\text{,}a\text{,}b\text{,}\omega )$ in (2.3), and put μ=m+ε$\mu =m+\epsilon$ and n=m+1+j$n=m+1+j$, where the pole index j is a non-negative integer; no singularities arise for negative integer j. Thus, cf. (2.3),

equation(3.1)

$D_{\exp }(p\text{;}\mu -n\text{,}a\text{,}b\text{,}\omega )=\frac{\mathrm{\Gamma }(-j+\epsilon )}{2^{-j+\epsilon }{a}^{(-j+\epsilon )/2}}U\left(\frac{-j+\epsilon }{2}\text{,}\frac{1}{2}\text{,}\frac{{(\mathrm{i}y)}^2}{2}\right)\text{,}$

Turn MathJax on

$y=(\omega -2p-\mathrm{i}b)/\sqrt{2a}$

equation(3.2)

$\mathrm{\Gamma }(-j+\epsilon )=\frac{{(-1)}^j}{\mathrm{\Gamma }(j+1)}\frac{1}{\epsilon }(1+\epsilon \psi (j+1)+\mathrm{O}({\epsilon }^2))$

Turn MathJax on

equation(3.3)

$D_{\exp }(p\text{;}-j-1+\epsilon \text{,}a\text{,}b\text{,}\omega )=\frac{1}{\epsilon }{D}_{\exp }^{\text{sing}}(j\text{;}p\text{,}a\text{,}b\text{,}\omega )+D_{\exp }^{\text{reg}}(j\text{;}p\text{,}a\text{,}b\text{,}\omega \text{;}\epsilon )\text{.}$

Turn MathJax on

equation(3.4)

$D_{\exp }^{\text{reg}}(j\text{;}p\text{,}a\text{,}b\text{,}\omega \text{;}0)=a^{j/2}\frac{{(-1)}^j{2}^j}{\mathrm{\Gamma }(j+1)}\left[U\left(\frac{-j}{2}\text{,}\frac{1}{2}\text{,}\frac{{(\mathrm{i}y)}^2}{2}\right)(\psi (j+1)-\frac{1}{2}\log (4a))+\frac{\mathrm{d}}{\mathrm{d}\epsilon }{\left.U\left(\frac{-j+\epsilon }{2}\text{,}\frac{1}{2}\text{,}\frac{{(\mathrm{i}y)}^2}{2}\right)\right|}_{\epsilon =0}\right]\text{,}$

Turn MathJax on

$D_{\exp }^{\text{reg}}(j\text{;}p\text{,}a\text{,}b\text{,}\omega \text{;}\epsilon )=D_{\exp }^{\text{reg}}(j\text{;}p\text{,}a\text{,}b\text{,}\omega \text{;}0)+\mathrm{O}(\epsilon )$

$D_{\exp }^{\text{sing}}$

equation(3.5)

$D_{\exp }^{\text{sing}}(j\text{;}p\text{,}a\text{,}b\text{,}\omega )=a^{j/2}\frac{{(-1)}^j}{\mathrm{\Gamma }(j+1)}{2}^{j}U\left(\frac{-j}{2}\text{,}\frac{1}{2}\text{,}\frac{{(\mathrm{i}y)}^2}{2}\right)\text{,}$

Turn MathJax on

equation(3.6)

$2^{j}U\left(\frac{-j}{2}\text{,}\frac{1}{2}\text{,}\frac{{(\mathrm{i}y)}^2}{2}\right)=H_j\left(\frac{\mathrm{i}y}{\sqrt{2}}\right)\text{.}$

Turn MathJax on

$D_{\exp }^{\text{sing}}$

equation(3.7)

$U\left(-n\text{,}\frac{1}{2}\text{,}\frac{{(\mathrm{i}y)}^2}{2}\right)=\frac{{(-1)}^n}{\sqrt{\pi }}\mathrm{\Gamma }(n+1/2){}_1{F}_1\left(-n\text{,}\frac{1}{2}\text{,}\frac{{(\mathrm{i}y)}^2}{2}\right)=2^{-2n}{H}_{2n}\left(\frac{\mathrm{i}y}{\sqrt{2}}\right)={(-1)}^n\frac{\mathrm{\Gamma }(2n+1)}{2^{2n}}{\displaystyle \sum _{k=0}^n}\frac{2^k{y}^{2k}}{\mathrm{\Gamma }(2k+1)\mathrm{\Gamma }(1+n-k)}\text{.}$

Turn MathJax on

$D_{\exp }^{\text{sing}}(2n\text{;}p\text{,}a\text{,}b\text{,}\omega )$

For odd integer pole index, we put j=2n+1,n⩾0$j=2n+1\text{,}n⩾0$, in the residue (3.5) to find

equation(3.8)

$U\left(\frac{-2n-1}{2}\text{,}\frac{1}{2}\text{,}\frac{{(\mathrm{i}y)}^2}{2}\right)={(-1)}^n\sqrt{\frac{2}{\pi }}\mathrm{\Gamma }(n+3/2)\mathrm{i}y{}_1{F}_1\left(-n\text{,}\frac{3}{2}\text{,}\frac{{(\mathrm{i}y)}^2}{2}\right)=2^{-2n-1}{H}_{2n+1}\left(\frac{\mathrm{i}y}{\sqrt{2}}\right)={(-1)}^n\frac{\mathrm{\Gamma }(2n+2)}{2^{2n+1}}\sqrt{2}\mathrm{i}y{\displaystyle \sum _{k=0}^n}\frac{2^k{y}^{2k}}{\mathrm{\Gamma }(2k+2)\mathrm{\Gamma }(1+n-k)}\text{.}$

Turn MathJax on

$D_{\exp }^{\text{sing}}(2n+1\text{;}p\text{,}a\text{,}b\text{,}\omega )$

At very small but finite ε$\epsilon$, the regular part in (3.3) can be calculated with high precision as

equation(3.9)

$D_{\exp }^{\text{reg}}(j\text{;}p\text{,}a\text{,}b\text{,}\omega \text{;}\epsilon )=D_{\exp }(p\text{;}-j-1+\epsilon \text{,}a\text{,}b\text{,}\omega )-\frac{1}{\epsilon }{D}_{\exp }^{\text{sing}}(j\text{;}p\text{,}a\text{,}b\text{,}\omega )\text{,}$

Turn MathJax on

$D_{\exp }^{\text{sing}}$

4. Regularized Hankel series

4.1. Pole subtraction

In Section 3, we performed the ε$\epsilon$ expansion of integral D_exp$D_{\exp }$ in (2.1), and calculated the pole term $D_{\exp }^{\text{sing}}/\epsilon$, cf. (3.3). The pole separation and ε$\epsilon$ expansion of the trigonometric integrals D_cos,sin,0$D_{\cos \text{,}\sin \text{,}0}$ defined in (2.6), (2.7) and (2.8) can be assembled from the coefficients $D_{\exp }^{\text{sing}}$ and $D_{\exp }^{\text{reg}}$ in (3.3).

As in (3.1), we put μ=m+ε$\mu =m+\epsilon$ and n=m+1+j$n=m+1+j$, where j is a non-negative integer, and write the integrals D_cos,sin,0(p;μ-n,a,b,ω)$D_{\cos \text{,}\sin \text{,}0}(p\text{;}\mu -n\text{,}a\text{,}b\text{,}\omega )$ in (2.6), (2.7) and (2.8) as

equation(4.1)

$D_{\cos \text{,}\sin }(p\text{;}-j-1+\epsilon \text{,}a\text{,}b\text{,}\omega )={\int }_0^{\infty }{k}^{-j-1+\epsilon }{\mathrm{e}}^{-{\mathit{ak}}^2-(b+\mathrm{i}\omega )k}(\cos (2\mathit{pk})\text{,}\sin (2\mathit{pk}))\mathrm{d}k=\frac{1}{\epsilon }{D}_{\cos \text{,}\sin }^{\text{sing}}(j\text{;}p\text{,}a\text{,}b\text{,}\omega )+D_{\cos \text{,}\sin }^{\text{reg}}(j\text{;}p\text{,}a\text{,}b\text{,}\omega \text{;}\epsilon )$

Turn MathJax on

equation(4.2)

$D_0(-j-1+\epsilon \text{,}a\text{,}b\text{,}\omega )={\int }_0^{\infty }{k}^{-j-1+\epsilon }{\mathrm{e}}^{-{\mathit{ak}}^2-(b+\mathrm{i}\omega )k}\mathrm{d}k=\frac{1}{\epsilon }{D}_0^{\text{sing}}(j\text{;}a\text{,}b\text{,}\omega )+D_0^{\text{reg}}(j\text{;}a\text{,}b\text{,}\omega \text{;}\epsilon )\text{.}$

Turn MathJax on

$D_{\cos }^{\text{reg}}(j\text{;}p\text{,}a\text{,}b\text{,}\omega \text{;}\epsilon )=\frac{1}{2}(D_{\exp }^{\text{reg}}(j\text{;}p\text{,}a\text{,}b\text{,}\omega \text{;}\epsilon )+D_{\exp }^{\text{reg}}(j\text{;}-p\text{,}a\text{,}b\text{,}\omega \text{;}\epsilon ))\text{,}$

Turn MathJax on

$D_{\sin }^{\text{reg}}(j\text{;}p\text{,}a\text{,}b\text{,}\omega \text{;}\epsilon )=\frac{1}{2\mathrm{i}}(D_{\exp }^{\text{reg}}(j\text{;}p\text{,}a\text{,}b\text{,}\omega \text{;}\epsilon )-D_{\exp }^{\text{reg}}(j\text{;}-p\text{,}a\text{,}b\text{,}\omega \text{;}\epsilon ))\text{,}$

Turn MathJax on

equation(4.3)

$D_0^{\text{reg}}(j\text{;}a\text{,}b\text{,}\omega \text{;}\epsilon )=D_{\exp }^{\text{reg}}(j\text{;}p=0\text{,}a\text{,}b\text{,}\omega \text{;}\epsilon )\text{,}$

Turn MathJax on

$D_{\exp }^{\text{reg}}$