Volume 225, 1 December 2013, Pages 228–241

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

• Department of Physics, Hiroshima University, 1-3-1 Kagami-yama, Higashi-Hiroshima 739-8526, Japan

## Highlights

The high-index evaluation of integrals containing squared spherical Bessel functions is studied.

The integrals arise as spectral averages in multipole expansions of spherical Gaussian random fields.

A high-precision integration technique based on finite Hankel series in epsilon regularization is developed.

An Airy approximation of the integrals is derived using uniform Nicholson asymptotics.

The finite series evaluation is compared with the Airy approximation over an extended range of Bessel indices.

## Abstract

Bessel integrals of type are studied, where the squared spherical Bessel function 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.

## Keywords

• Squared spherical Bessel functions;
• Regularization of Hankel series;
• Gaussian power-law densities;
• Kummer distributions;
• Airy approximation of Bessel integrals

## 1. Introduction

We investigate the Bessel integrals arising in the multipole expansion of isotropic Gaussian random fields on the unit sphere [1] and [2]. is a squared spherical Bessel (ssB) function of integer index l⩾0,p is a positive scale parameter, a>0,b and ω are real constants, and μ is an integer power-law exponent.

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 μ, 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 in the above integrand). In this case, the confluent functions in the Hankel coefficients are still singular at integer power-law exponents μ 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 ε pole extraction to derive combinatorial identities (sum rules) for the Hankel coefficients, by making use of the fact that the singular ε 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<img height="23" border="0" style="vertical-align:bottom" width="201" alt="View the MathML source" title="View the MathML source" src="http://origin-ars.els-cdn.com/content/image/1-s2.0-S0096300313010072-si118.gif">$∫0∞kμ+2e-ak2-(b+iω)kjl2(pk)dk$

We study the Bessel average , where the distribution g(k)=kμe-ak2-(b+iω)k is a modulated Gaussian power law with real exponents μ,a>0,b and ω (or a=0 and b>0). The scale parameter p in the argument of the squared spherical Bessel (ssB) function is positive. Spherical Bessel functions are rescaled ordinary Bessel functions of half-integer order, , the index l being a non-negative integer [3] and [4]. These integrals converge for μ+2+2l>-1, as the ascending series of starts with the power x2l.

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

where D-(μ+1) denotes the cylinder function and
equation2.2
A complex exponent a is possible as well, with Re a>0 or Re a=0 and b>0, in which case the principal value of the root in (2.2) is assumed. We will also consider the limit case a=0 of integral (2.1) with positive b, cf. (2.19). Integral (2.1) also represents a Kummer function [5],where U is related to the confluent hypergeometric function 1F1 by [4] and [5]so that
equation2.5
In the second equality, we made use of the transformation e-x1F1(a,b,x)=1F1(b-a,b,-x). The representation (2.4) of the Kummer function is unambiguous at the branch cut of U(a,1/2,z), that is for real y, so that (2.5) is suitable for numerical calculations.

We split the exponential exp(2ipk) in integral (2.1) into real and imaginary parts, defining

The functions Dexp and Dcos,sin,0 will be considered as analytic continuations in μ of the respective integrals, which converge for Re μ>-1. This continuation is effected by the Kummer function in (2.4) or by (2.5). In the following sections, we study negative integer exponents μ, where the gamma function in (2.3) becomes singular, which can be dealt with by epsilon expansion. In this section, we assume a non-integer power-law exponent μ.

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

where and are polynomials in 1/x, explicitly stated in (2.10), (2.11), (2.12), (2.13), (2.14) and (2.15). In the following, we use the Hankel symbolThe polynomial in (2.9) readswith coefficients
equation2.12
where 0⩽kl,l⩾0, and [i,j] denotes Hankel’s symbol (2.10).

As for the polynomials in (2.9), we define and

with coefficientswhere 0⩽kl-1 and l⩾1. Finally, the polynomials in (2.9) readwhere l⩾0.

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

and interchange the integration with the summation of series and . In this way, we obtain integral (2.16) as linear combination of three finite Hankel series:whereThe series coefficients a2k(l),b2k+1(l), and c2k(l) are stated in (2.12) and (2.14), and (2.15), and the functions Dcos,sin,0 are defined in (2.6), (2.7) and (2.8), with Dexp in (2.5) substituted. We also note that ; a sum is regarded as void if the lower summation boundary exceeds the upper one. We thus obtain the Bessel integral DssB(l,p;μ,a,b,ω) in (2.16) as a finite linear combination of confluent hypergeometric functions, cf. (2.5).

If a=0 and b>0 in integral (2.16), then integral Dexp in (2.1) is elementary, and the confluent hypergeometric Dexp(p;μ,a,b,ω) in (2.3) and (2.5) is replaced by

We will further discuss this limit in Section 2.2. In the opposite limit, a,y→0, cf. (2.2), we findwhich can readily be seen from (2.5) since 1F1∼1. This limit is best dealt with in Laplace asymptotics, cf. Section 2.3.

### 2.2. A limit case: Kummer averages<img height="22" border="0" style="vertical-align:bottom" width="175" alt="View the MathML source" title="View the MathML source" src="http://origin-ars.els-cdn.com/content/image/1-s2.0-S0096300313010072-si193.gif">$∫0∞kμ+2e-(b+iω)kjl2(pk)dk$

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

Integral Dexp in (2.1) is now elementary,The trigonometric components (2.6), (2.7) and (2.8) of this integral are elementary as well,andWe note that β=b+iω can be complex, with b>0, and the power-law exponent μ in integral (2.21) is usually real, although the functions Dcos,sin,0 will be regarded as analytic μ continuations of the respective integrals. The Hankel series in (2.18) are compiled with these elementary functions, and integral DssB(l,p;μ,0,b,ω) in (2.21) is assembled as stated in (2.17). If the power-law exponent μ is integer, the coefficients Dcos,sin,0 of the Hankel series (2.18) become singular, which requires ε expansion, performed in Section 5 for the Kummer averages (2.21).

### 2.3. Laplace/Fourier asymptotics of integral <img height="23" border="0" style="vertical-align:bottom" width="201" alt="View the MathML source" title="View the MathML source" src="http://origin-ars.els-cdn.com/content/image/1-s2.0-S0096300313010072-si215.gif">$∫0∞kμ+2e-ak2-(b+iω)kjl2(pk)dk$

We calculate integral DssB(l,p;μ,a,b,ω) in (2.16) by making use of the ascending series expansion of the squared spherical Bessel function [6],

We note the asymptotic expansion of the series coefficients α2k(l) for large k and fixed l,obtained by means of Stirling’s estimate [5]The Laplace/Fourier series of the Bessel integral (2.16) is found by substituting the ascending series (2.26), and by interchanging integration and summation,where the integrals D0 are calculated in (2.5) and (2.8), with , cf. (2.2). The series (2.29) can also be used for integer exponents μ, as no singularities arise, provided that integral (2.16) converges at the lower integration boundary, which is the case if μ+2+2l>-1.

In series (2.29), we consider the limit a=0 and b>0, so that D0(μ,0,b,ω)=Γ(μ+1)/(b+iω)μ+1, cf. (2.19). Series (2.29) is thus majorized by a polylogarithm , with z=4p2/(b+iω)2, which converges in the unit disk |z|<1 for arbitrary μ[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 ω (Fourier asymptotics).

In the opposite limit a, we can use (2.20) as estimate for D0(μ,a,b,ω) in series (2.29). In this case, series (2.29) is convergent due to a factor e-k(logk-1) emerging in the series coefficients for large k, cf. (2.27) and (2.28). However, this factor only emerges for kl (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ω| in integral (2.16). This expansion also works if both a and |b+iω| are large and y is moderate, cf. (2.1) and (2.2). In brief, if the exponents a and/or |b+iω| 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 and ω 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 μ in integral , cf. (2.16), where is a squared spherical Bessel function of integer index l⩾0, and p a positive scale parameter. Here and in Sections 4 and 5, we study integer power-law exponents μ=m, so that is averaged with the distribution kme-ak2-(b+iω)k, with a>0, real b and ω. The k2 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, so that the integral is convergent. We will calculate integral (2.16) by performing ε expansion for arbitrary integer power-law exponents μ=m⩾-2l-2.

The analytic μ continuations Dcos,sin(p;μ-n,a,b,ω) and D0(μ-n,a,b,ω) 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+ε, where m is an arbitrary integer, and n=m+1+j, so that 1+μ-n=-j+ε. Poles in Dexp(p;μ-n,a,b,ω) (defined by (2.3) or (2.5)) and in its trigonometric components Dcos,sin,0 in (2.6), (2.7) and (2.8) emerge only at integer μ-n=-j-1 with j⩾0.

We start with the ε expansion of integral Dexp(p;μ-n,a,b,ω) in (2.3), and put μ=m+ε and 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),

where μ-n=-j-1+ε, and we use the shortcut , cf. (2.2). In (3.1), we substitute the ε expansion of the gamma function, assuming integer j⩾0,and split off the ε pole,The psi function in (3.2) is the logarithmic derivative of the gamma function. The constant term in the ε Laurent expansion of Dexp is
equation3.4
so that . The residue in (3.3) readswhere we can substitute the Hermite polynomial, cf. (3.7) and (3.8),In , we consider even integer pole indices j=2n,n⩾0, and note [5]
equation3.7
The singular term in (3.5) is thus real for real y.

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

equation3.8
Here, the Kummer function U at the branch cut (real y) and the Hermite polynomial are imaginary. This defines the singular term in (3.5) for odd pole index, which is imaginary at b=0, that is for real y, cf. (2.2).

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

with explicitly given in (3.5), (3.7) and (3.8). This can be more efficient than the numerical differentiation of the Kummer function with respect to the ε parameter in (3.4).

## 4. Regularized Hankel series

### 4.1. Pole subtraction

In Section 3, we performed the ε expansion of integral Dexp in (2.1), and calculated the pole term , cf. (3.3). The pole separation and ε expansion of the trigonometric integrals Dcos,sin,0 defined in (2.6), (2.7) and (2.8) can be assembled from the coefficients and in (3.3).

As in (3.1), we put μ=m+ε and n=m+1+j, where j is a non-negative integer, and write the integrals Dcos,sin,0(p;μ-n,a,b,ω) in (2.6), (2.7) and (2.8) as

equation4.1
andBy comparing to (2.1) and (3.3), we identify the series coefficients as, cf. (2.6), (2.7) and (2.8),with in (3.9). The same relations (4.3) also hold for the residues in (4.1) and (4.2), with replaced by in (3.5).

Returning to the Hankel series and in (2.18), we put μ=m+ε. The gamma function in the analytic ε continuations Dcos(p;μ-2k,a,b,ω) and D0(μ-2k,a,b,ω) (defined by (2.3) and (2.6) and (2.8)) stays regular for ε→0 if the non-negative integer k satisfies 2km. It is singular for 2km+1. Accordingly, we split the series defining and in (2.18) as

and use Dcos(p;m-2k,a,b,ω) and D0(m-2k,a,b,ω) in (2.6) and (2.8) in the first summation like in (2.18). In the second summation in (4.4), we use, instead of Dcos and D0, the regular coefficients and in (4.3), with j=2k-1-m (so that μ-2k=-j-1+ε, cf. (4.1) and (4.2). The brackets in the upper and lower summation boundaries in (4.4) denote the integer part (that is, the largest integer less than or equal to the indicated argument). We also use the convention that a sum is zero if the lower summation boundary exceeds the upper one.

As for series in (2.18), with μ=m+ε and Dsin in (2.7), we note that the gamma function in Dsin(p;μ-2k-1,a,b,ω) stays regular for ε→0 if k satisfies 2k+1⩽m, cf. (2.3) and (2.7), whereas poles emerge if 2k+1⩾m+1. We split the series defining in (2.18) as

In the first summation, we can use Dsin(p;m-2k-1,a,b,ω) in (2.7). In the second summation, we replace Dsin by in (4.3) and put j=2k-m (since μ-2k-1=-j-1+ε). The pole terms in (4.1) and (4.2) can be ignored in the calculation of integral (4.6) below, as they cancel one another, cf. Appendix A.

The ε decomposition of integral (2.16) and its series representation (2.17) with power-law exponent μ=m+ε,0<ε≪1, is thus

equation4.6
where are the regularized Hankel series (2.18):
equation4.7
with Dcos calculated via (2.5) and (2.6), and via (3.9) and (4.3). The brackets in the summation boundaries denote the integer part, cf. (4.4). The second regularized Hankel series in (4.6) reads
equation4.8
with Dsin in (2.5) and (2.7), and in (3.9) and (4.3); the summation is explained in (4.5). The third Hankel series in (4.6) reads
equation4.9
with D0 defined in (2.5) and (2.8), and in (3.9) and (4.3). The series coefficients a2k(l),b2k+1(l) and c2k(l) are listed in (2.12), (2.14) and (2.15). The regularized Hankel series (4.7), (4.8) and (4.9) can be used for integer power-law exponents m⩾-2l-2 to calculate integral (4.6); the latter inequality is required for this integral to converge for ε→0. The series (4.6), (4.7), (4.8) and (4.9) remain valid if we put ε=0 and define in (4.3) by substituting in (3.4). This is actually done in Section 5, where we consider the special case of Kummer averages, where the confluent hypergeometric functions in (3.4) become elementary. In contrast, if ε is kept small but finite in the Hankel series (4.6), (4.7), (4.8) and (4.9), we can calculate the regular coefficients in (4.3) by means of in (3.9). As mentioned, the pole subtraction in (3.9) is handier than the calculation of the ε derivative of the confluent hypergeometric function required in (3.4), and can readily be done in high precision. High-precision evaluation of the coefficients Dcos,sin,0 and is necessary in the Hankel series (4.7), (4.8) and (4.9), since the combinatorial coefficients a2k(l),b2k+1(l) and c2k(l) become very large rational numbers at high Bessel index l. Therefore, the finite summations in (4.7), (4.8) and (4.9) give very large numbers, which cancel one another almost entirely in the total sum (4.6) resulting in precision loss, so that high-precision evaluation of these series is required. In Table 1 and Table 2, we give numerical examples of the Hankel expansion (4.6), (4.7), (4.8) and (4.9), covering a wide range of Bessel indices, and compare with the Airy approximation discussed in Section 4.2.

Table 1.

Hankel expansion and Airy approximation of the integrals , cf. (2.16) (p=1,μ=ω=0), with exponents a=7.0×10-3 and b=-0.5, employed in a multipole fit of the cosmic microwave background (CMB) temperature power spectrum [2]. The integrals are labeled by their Bessel index l (first column). The Hankel evaluation in the second column is based on the finite series (4.6), (4.7), (4.8) and (4.9). The Airy approximation in column 3 is calculated with (4.21). The indicated decimal digits of the Hankel series evaluation are all significant, the numbers being truncated. The Airy approximation is accurate to 3 decimal places for moderate Bessel index, but substantially deteriorates when the integrals become negligible as compared to their peak value.

lDssB(l) Hankel seriesDssB(l) Airy approx.
07.991354106316505 ×1047.992338873 ×104
17.999250409735869 ×1048.000256164 ×104
58.113326663187338 ×1048.114670865 ×104
108.481208178016143 ×1048.483594389 ×104
201.042557080617444 ×1051.042268591 ×105
301.238911735222280 ×1051.236931157 ×105
501.511142597127369 ×1041.509695163 ×104
701.503953736747144 ×1011.180405467 ×101
1008.233597800611870 ×10-88.635323282 ×10-9
1509.739275886661490 ×10-304.028033041×10-36
Full-size table
Table 2.

Comparison of Hankel evaluation and Airy approximation of the integrals , cf. (2.16) and (4.6), with exponents a=1.067×10-4 and b=-0.11. The caption to Table 1 applies. The order l of the spherical Bessel function is listed in the first column, the finite Hankel series expansion (see Sections 4 and 6) in column 2, and the Airy approximation (4.21) in column 3. The indicated digits of the Hankel evaluation are significant, whereas 5 to 6 digit accuracy is reached by the Airy approximation. The Hankel expansion for Bessel indices above l∼400 requires very high precision evaluation of the confluent hypergeometric functions in (2.5), and is therefore not efficient at high l, also see the remarks following (4.9).

lDssB(l) Hankel seriesDssB(l) Airy approx.
01.761712987728264 ×10141.761713864 ×1014
11.761720004433907 ×10141.761720881 ×1014
51.761818247890288 ×10141.761819125 ×1014
101.762099042033407 ×10141.762099922 ×1014
501.770734121286065 ×10141.770735073 ×1014
1001.798359213865244 ×10141.798360424 ×1014
2001.924862478039058 ×10141.924865675 ×1014
3002.233775172960886 ×10142.233790888 ×1014
4003.075282836068884 ×10143.075217334 ×1014
6009.669919328 ×1013
8002.505261104 ×1010
10001.580371521 ×103
12002.090425596 ×10-8
Full-size table

### 4.2. Consistency checks for Hankel series: comparison with Weber integrals and Airy approximation

For m=0,ε=0 and b=ω=0, integral (4.6) degenerates into a Weber integral [6],e

convergent for Re a>0. A finite representation of the spherical Bessel function jl(x) valid for complex argument iswith the Lommel polynomials [5]The brackets of the upper summation boundaries indicate the integer part (largest integer less than or equal to the indicated argument), not to be confused with the Hankel symbol in the summands,More generally, we may compare integral DssB(l,p;m+ε,a,b,ω) in (4.6) to its Airy approximation. To this end, we consider an arbitrary kernel function g(k) and writewhere we may put g(k)=km+εe-ak2-(b+iω)k to recover integral DssB(l,p,m+ε,a,b,ω). We rescale the integration variable in (4.14) with (l+1/2)/p,and substitute the high-l Nicholson approximation of the spherical Bessel function [4],whereIn this transitional steepest-descent approximation, a squared Airy function appears in the integrand of (4.15), which admits the integral representation [8]also see [9], [10] and [11] for representations of Airy functions by trigonometric kernels. The indicated approximation in (4.18) holds for real argument and large z, so that the cubic term in the argument of the cosine can be dropped (Riemann–Lebesgue lemma). θ(z) is the Heaviside step function. Thus we put in the interval 0<x<1, where ξ(x) is positive, which means to replace the lower integration boundary in integral (4.15) by 1. In the range x⩾1, we introduce a new integration variable, , and substitute (4.18) into the squared Nicholson approximation (4.16), with , to obtainApplying this approximation and the indicated variable change y=x2-1 to integral (4.15), we obtain the Airy approximation of integral (4.14),On specifying the kernel function g as stated after (4.14), we find the Airy approximation of integral (4.6),
equation4.21
We may apply a further variable change, 1+y=(1+t)2,y=t(2+t), to remove the root in the exponential,
equation4.22
The Airy approximation (4.21) or (4.22) is numerically tame and turns out to be efficient even at very low Bessel index l, but its accuracy is limited as compared with the Hankel series evaluation in (4.6), cf. Table 1, Table 2 and Table 3. The Airy approximation of integrals containing products of Bessel derivatives is discussed in [12].

Table 3.

Hankel expansion and Airy approximation of the Kummer averages in (5.10) (p=1,m=0), with exponents b=2.3×10-3 and ω=2.15×10-2 defining the oscillatory multipole component of the CMB fluctuations in [2]. In columns 2 and 4, the Hankel evaluation of integral Re  and is recorded as a function of Bessel index l in column 1. The finite Hankel series expansion of Kummer averages is explained in Section 5 and compared with the Airy approximation (4.21) in columns 3 and 5. The indicated digits of the Hankel evaluation are significant. The Airy approximation is thus accurate to 3 or 4 decimal digits.

lRe DssB(l) Hankel seriesRe DssB(l) Airy approx.ImDssB(l) Hankel seriesImDssB(l) Airy approx.
02.4593949946914532.456975786−2.29953729177942 ×101−2.29993855 ×101
12.4390237403995052.436746712−2.30369798295024 ×101−2.30396971 ×101
52.1766235464024762.174578376−2.34117813788856 ×101−2.34127427 ×101
101.4695375772071731.467707567−2.41348458087903 ×101−2.41348612 ×101
50−1.49154918284621 ×101−1.49147674 ×101−2.66093614013266 ×101−2.66074707 ×101
100−3.54240464462342 ×101−3.54214934 ×101−3.37196089026216−3.37281990
3002.126889250211135 ×1012.127187079 ×101−3.07546677806975 ×101−3.07505100 ×101
5001.536265267664984 ×1011.535679165 ×1012.612978498620274 ×1012.613253539 ×101
8001.214846414087418 ×1011.214373489 ×1011.484542229293373 ×1011.484978720 ×101
1000−1.30411841949000 ×101−1.30432858 ×1013.5989465845212003.593693228
50004.215571594 ×10-4−3.02117450 ×10-3
104−2.26686563 ×10-8−3.73505755 ×10-8
2×104−5.96529699 ×10-182.147340549 ×10-18
Full-size table

## 5. Kummer averages<img height="22" border="0" style="vertical-align:bottom" width="177" alt="View the MathML source" title="View the MathML source" src="http://origin-ars.els-cdn.com/content/image/1-s2.0-S0096300313010072-si416.gif">$∫0∞km+2e-(b+iω)kjl2(pk)dk$ with integer power-law index

The Kummer averages (2.21) are a limit case of integral (2.16), where the quadratic term in the exponent of the integrand vanishes, a=0. We can proceed as in Section 4.1, starting with the ε expansion of integral Dcos(p;μ-n,a=0,b,ω) in (2.23), where μ=m+ε and n=m+1+j, with integer j⩾0. As in Section 2.2, we use the shortcut β=b+iω, and write, cf. (2.23) and (4.1),

equation5.1
The singular coefficient isIn contrast to (3.4), the regular constant term of the ε expansion in (5.1) is elementary,
equation5.3
Here, j is a non-negative integer; no singularities arise for negative integer j. Principal values are assumed for the logarithms. These coefficients can readily be derived by substituting the ε expansion (3.2) of the gamma function into (5.1). The constant coefficient (5.3) differs from in (5.1) by terms of O(ε).

Analogously, integral Dsin(p;μ-n,0,b,ω) in (2.24) admits the ε expansion

equation5.4
with residueand the regular constant coefficient
equation5.6
Finally, the ε expansion of integral D0(μ-n,0,b,ω) in (2.25) reads, cf. (4.2),whereand β=b+iω.

As mentioned, in contrast to the Gaussian power-law averages (4.6), the regular coefficients of Kummer averages are elementary functions at ε=0. Therefore, it is numerically efficient to perform the regularized Hankel series expansion for integer power-law index m directly at ε=0. The analogue to the series expansion (4.6) is thus

equation5.10
The Hankel series in (5.10) are assembled by substituting the regular constant coefficients calculated in (5.3), (5.6) and (5.9) into series (4.7), (4.8) and (4.9) (where we put ε=0). Also required in these series are the non-singular coefficients Dcos,sin,0 listed in (2.23), (2.24) and (2.25). Integral (5.10) converges if m⩾-2l-2, with integer Bessel index l⩾0 and exponent b>0, since jl(x)∝xl(1+O(x2)) and jl(x)=O(1/x) for large argument.

As a consistency check, we note, cf. (5.2) and (5.5),

which coincides with in (3.4) in the limit y (that is a→0 in (2.2)), cf. (A.4). A check of the regular constant coefficients in (5.3) and (5.6) is provided bywhich can be recovered from in (3.4) in the same limit y by substituting U(a,b,z)∼z-a.

A different consistency check of the Hankel expansion (5.10) is obtained for power-law index m=-1. In this case, integral (5.10) can be calculated as [6]

where Ql(x) is a Legendre function of the second kind with branch cut -1⩽x⩽1[4], which is elementary for integer degree l⩾0:where Pl(x) is the Legendre polynomialand γE is Euler’s constant. Other integer power-law indices m in (5.10) can be checked by considering multiple β derivatives or antiderivatives of the Beltrami integral (5.13), cf. [13]. The sum rules (A.10) and (A.11) as well as the closed expression (2.15) for the Hankel coefficients c2k(l) can be obtained by comparing the finite Legendre series in (5.13) to the Hankel expansion (5.10). Alternatively, the derivation of the combinatorial identities in Appendix A is based on the residues stated in (5.2), (5.5) and (5.8), and the Hankel coefficients c2k(l) in (2.15) are derived by making use of Lommel’s identity [12]. Table 3 illustrates the Hankel evaluation (5.10) of Kummer averages and their Airy approximation (4.21).

## 6. Results, discussion & conclusion

We studied Bessel integrals with squared spherical Bessel functions of integer index l⩾0 in the integrand,

averaged with the distribution kme-ak2-(b+iω)k. The constant p in the argument of the Bessel function is a positive scale parameter. The constants a,b and ω defining the modulated exponential are real, and a is positive, ensuring convergence at the upper integration boundary. (The limit case a=0 with b>0 is also convergent, defining Kummer averages worked out in Section 5. Here, we summarize the general case, a>0.) The power-law index m is integer and can be negative, provided that m+2⩾-2l, which is the condition for integral (6.1) to converge at the lower integration boundary. In multipole expansions of spherical random fields used to analyze temperature fluctuations of the cosmic background radiation, the exponents m,a,b,ω and the scale parameter p are kept constant, and the multipole index l varies over a wide range; presently multipoles up to l∼104 are resolvable [1] and [2]. The complex exponent b+iω generates modulations in the multipole spectrum.

We start by summarizing the precision evaluation of integral (6.1) based on regularized Hankel series, cf. Section 4. Integral (6.1) can be assembled as

with the finite Hankel series
The brackets in the summation boundaries denote the integer part. The series coefficients a2k(l),b2k+1(l) and c2k(l) are stated in (2.12), (2.14) and (2.15). These series are derived by making use of the fact that the asymptotic Hankel expansion of Bessel functions [5] in reciprocal powers of the argument terminates in the case of spherical Bessel functions, cf. Section 2. The coefficients Dcos,sin,0 and in (6.3) are obtained by making use of term-by-term integration combined with analytic continuation in the power-law index m treated as complex variable. Analytic continuation is required because of divergent integrals arising in the term-by-term integration, and can be performed with confluent hypergeometric functions, cf. Section 2: The coefficients Dcos,sin,0 in (6.3) are explicitly defined byand D0(z)=Dexp(0,z), whereA representation of the Kummer function U in terms of the confluent hypergeometric function 1F1 is given in (2.4), which can also be used at the branch cut of U, that is for real argument y in (6.5).

At negative integer z=-j-1,j⩾0, pole singularities emerge in Dexp and thus in the coefficient functions Dcos,sin,0(z). The regularized coefficients in (6.3) read, after pole subtraction,

and , wherewith the Kummer function U in (2.4). The argument y is the same as in (6.5). The regularized coefficients are defined for integer j⩾0. The pole subtraction is performed by means of epsilon expansion, cf. Section 3,with residue stated in (3.5). The ε residues of the coefficients Dcos,sin,0 are defined analogous to (6.8) and calculated as in (6.6) with replaced by .

The residues do not enter in the actual calculation of integral (6.1). If we replace the regular coefficients in the Hankel series (6.3) by the residues , then the second summation in each of the three series in (6.3) must vanish identically. (The poles cancel one another because integral (6.1) is finite.) This pole cancellation happens due to combinatorial identities, sum rules satisfied by the series coefficients a2k(l),b2k+1(l) and c2k(l) in (6.3), which are derived in Appendix A. These identities, cf. (A.10) and (A.11), provide useful consistency checks for the series coefficients in high-precision calculations, especially for large Bessel index l; they only involve rational numbers, and are independent of the parameters defining integral (6.1). Finally, the Bessel index l and the power-law exponent m in the Hankel series (6.3) satisfy l⩾0 and the convergence condition m⩾-2l-2 of integral (6.1). It is evident from the summation boundaries in (6.3) (second summation signs), that regularization by pole subtraction (6.8) is required if the Bessel index exceeds m/2. Since the power-law exponent m of integral (6.1) is usually a small integer, this already happens at very low Bessel index in multipole expansions.

The Hankel method summarized in (6.1), (6.2), (6.3), (6.4), (6.5), (6.6), (6.7) and (6.8) is based on finite series, and can be used for integer Bessel index l, allowing safe evaluation in any desired precision. We have also discussed asymptotic techniques and special cases of integral (6.1) providing cross-checks by numerical comparison. In Section 2.3, we studied the Laplace and Fourier asymptotics of integral (6.1), namely the limit a or |b+iω|→, keeping the Bessel index l fixed. This limit, however, is not applicable in multipole expansions, in which the exponents are moderate and the Bessel index becomes large. In Section 4.2, we pointed out a Weber integral defined by exponents m=0 and b=ω=0 in (6.1), which can be expressed by a spherical Bessel function with imaginary argument. Weber integrals with even integer power-law exponent m are obtained as multiple a derivatives or antiderivatives of identity (4.10), cf. [13]. In Section 4.2, we also derived the Airy approximation of integral (6.1), making use of Nicholson’s steepest-descent asymptotics of the squared Bessel function in the integrand of (6.1). In Table 1 and Table 2, we compare Hankel expansion and Airy approximation by listing integral (6.1) for various Bessel indices.

In Sections 2.2 and 5, we studied the special case of Kummer averages, that is integral (6.1) with a=0 and b>0. In this limit case, the confluent hypergeometric functions in (6.5) and (6.7) become elementary power laws; the Hankel series of these averages are explicitly compiled in Section 5. In this section, we also mentioned a Beltrami integral defined by exponents m=-1,a=0 in (6.1), which can be expressed as an elementary Legendre function of the second kind and integer degree l, cf. (5.14). Table 3 compares the Hankel evaluation of Kummer averages with their Airy approximation (4.21), illustrating the accuracy of the Airy approximation over an extended range of Bessel indices.

## Appendix A. Combinatorial identities for Hankel coefficients

First, we point out an identity relating to the residues and in (4.1) and (4.2). We start with integral (4.6), where m⩾-2l-2 (convergence condition). This integral admits the ε expansion (0<ε≪1)

and is finite for ε→0. Accordingly, the residue in (A.1) must vanish identically (and is therefore not indicated in (4.6)), and we can thus identify DssB(l,p;m,a,b,ω) with the constant term in (A.1). The residue is assembled like the regularized Hankel series in (4.6),where denote the finite seriesThese series are obtained by replacing the regular terms in the respective second series in (4.7), (4.8) and (4.9) by their singular counterparts , cf. (4.1) and (4.2). The latter are assembled by means of in (3.5) (with the Hermite polynomials in (3.7) and (3.8) substituted) and the singular counterpart to identities (4.3). The fact that identity (A.2) holds for arbitrary constants a,b and ω as well as integers l⩾0 and m⩾-2l-2 (convergence condition of integral (4.6) at the lower boundary) allows for significant consistency checks, in particular regarding numerical precision, since the Hankel coefficients a2k(l),b2k+1(l) and c2k(l) become very large rational numbers for high Bessel index l.

We consider identity (A.2) in the limit a→0 , b=0,ω→0. The left-hand side of this identity is a polynomial in a as well as ω, since is a polynomial, cf. (3.5) and (3.6). Thus this identity stays valid for a=ω=0, even if integral (4.6) does not converge at the upper integration boundary. (In this limit, integral (2.16) degenerates into a Schafheitlin integral [5], convergent in the strip -2l-3<Reμ<-1.) In (3.6), we perform the limit |y|→, replacing the Hermite polynomial by its asymptotic limit Hj(z)∼2jzj and the Kummer function by

This limit is realized by a→0 in , cf. (2.2). On substituting (A.4) into (3.5), we obtainand the trigonometric components are defined as in (4.3). In this way, we find series in (A.3) at a=b=0 asHere, we perform the limit ω→0. It is easy to see that vanishes if m is even or if m⩽-3, since all powers of ω in the series are positive. A non-zero limit is obtained only for odd m⩾-1, provided that l⩾(m+1)/2 (so that the sum is not void). In this case, the first term of the series is independent of ω. Thus, for m=-1,1,3,… and l⩾(m+1)/2, we obtainReturning to the trigonometric components in (4.3), we substitute in (A.5) and perform the limit ω→0 to findwhereas and vanish. Accordingly, the Hankel series and in (A.3) vanish for even m, like series discussed above. If m is odd, we find, by substituting (A.8) into (A.3),We can draw the following conclusions. The identity holds for odd m⩽-3, since the Hankel series vanishes in this case. Thus we find, for m=-3,-5,-7… and l⩾-(m+1)/2 [the latter being the convergence condition m⩾-2l-2 of integral (4.6) for odd power-law exponent m, required for identity (A.2) to hold], the sum rulesFor odd integer still vanishes if l<(m+1)/2, but the summations in (A.9) are void in this case, and no new identities are obtained.

It remains to consider the case of odd m⩾-1, with l⩾(m+1)/2, where does not vanish and is given by (A.7). In this case, identity holds according to (A.2), and we find, by substitution of the Hankel series (A.7) and (A.9),

equationA.11
This sum rule is valid for odd integer m⩾-1 and integer l⩾(m+1)/2, whereas identity (A.10) applies for odd m⩽-3 and l⩾-(m+1)/2. The rational identities (A.10) and (A.11) provide significant numerical consistency checks for the combinatorial coefficients a2k(l),b2k+1(l) and c2k(l) stated in (2.12), (2.14) and (2.15).

## References

• [1]
• Planck Collaboration et al., Planck 2013 results. I. Overview of products and scientific results, 2013. arXiv:1303.5062.
• [3]
• Scattering Theory of Waves and Particles

• Springer, New York (1982)

• [4]
• F.W.J. Olver, D.W. Lozier, R.F. Boisvert, C.W. Clark (Eds.), NIST Handbook of Mathematical FunctionsCambridge University Press, Cambridge (2010) <http://dlmf.nist.gov>

• [5]
• Formulas and Theorems for the Special Functions of Mathematical Physics

• Springer, New York (1966)

• [6]
• A Treatise on the Theory of Bessel Functions

• Cambridge University Press, Cambridge (1996)

• [7]
• Polylogarithms and Associated Functions

• North-Holland, Amsterdam (1981)

• [9]
• Electric-field effects on optical absorption near thresholds in solids

• Phys. Rev., 147 (1966), pp. 554–566

• [11]
• Airy Functions and Applications to Physics

• (second ed.)Imperial College Press, London (2010)

• [13]
• Weber and Beltrami integrals of squared spherical Bessel functions: Finite series evaluation and high-index asymptotics

• Math. Methods Appl. Sci. (2013) http://dx.doi.org/10.1002/mma.2882