We consider a new identity involving integrals and sums of Bessel functions. The identity provides new ways to evaluate integrals of products of two Bessel functions. The identity is remarkably simple and powerful since the summand and the integrand are of exactly the same form and the sum converges to the integral relatively fast for most cases. A proof and numerical examples of the identity are discussed.
The Newtonian kernel 1.1is ubiquitous in mathematical physics and is essential to an understanding of both gravitation and electrostatics (Kellogg 1967). It is central in classical mechanics (Goldstein 1980), but plays an equally important role in quantum mechanics (Landau & Lifshitz 1965), where it mediates the dominant two-particle interaction in electronic Schrödinger equations of atoms and molecules.
Although the Newtonian kernel has many beautiful mathematical properties, the fact that it is both singular and long-ranged is awkward and expensive from a computational point of view (Hockney & Eastwood 1981) and this has led to a great deal of research into effective methods for its treatment. Of the many schemes that have been developed, Ewald partitioning (Ewald 1921), multi-pole methods (Greengard 1987) and Fourier transform techniques (Payne et al. 1992) are particularly popular and have enabled the simulation of large-scale particulate and continuous systems, even on relatively inexpensive computers.
A recent alternative (Gilbert 1996; Varganov et al. 2008; Gill & Gilbert 2009; Limpanuparb & Gill 2009, 2011; Limpanuparb 2011; Limpanuparb et al. 2011, 2012) to these conventional techniques is to resolve (1.1), a non-separable kernel, into a sum of products of one-body functions 1.2where Y lm(r) is a spherical harmonic (Olver et al. 2010) of the angular part of three-dimensional vector r, 1.3Jl(z) is a Bessel function of the first kind (Olver et al. 2010), and r=|r|. The resolution (1.2) is computationally useful because it decouples the coordinates r and r′ and allows the two-body interaction integral 1.4between densities ρa(r) and ρb(r) to be recast as 1.5where Anlm is a one-body integral of the product of ρa(r) and ϕnlm(r). If the one-body integrals can be evaluated efficiently and the sum converges rapidly, (1.5) may offer a more efficient route to E[ρa,ρb] than (1.4).
The key question is how best to obtain the Kl resolution Previous attempts (Gilbert 1996; Varganov et al. 2008; Gill & Gilbert 2009; Limpanuparb & Gill 2009) yielded complicated Knl whose practical utility is questionable but, recently, we have discovered the remarkable identity 1.6where εn is defined by 1.7a,b∈[0,π], and we take the appropriate limit for the n=0 term in the sum (1.6). This yields (Limpanuparb et al. 2011) the functions and these provide a resolution that is valid, provided that r<π. We note that ϕnlm vanish for n=0 unless l=0.
If we write (1.6) as an integral from to , 1.8the summand and integrand are of exactly the same form.
There have been a number of studies of this kind of sum-integral equality by various groups, for example, Krishnan & Bhatia in the 1940s (Bhatia & Krishnan 1948; Krishnan 1948a,b; Simon 2002) and Boas, Pollard & Shisha in the 1970s (Boas & Stutz 1971; Pollard & Shisha 1972; Boas & Pollard 1973). Their discovery was inspired by a practice to ‘approximate’ an intractable sum that arises in physics by an integral. They realized that the ‘approximation’ was in fact exact for a number of cases.
In this paper, however, our goal is the opposite. We originally aimed to use the sum to approximate the integral but later found that it was exact. Our identity (1.6) is also considerably different from theirs, but we may regard theirs (Boas & Pollard 1973) when c=0 as a special case of ours when ν=1/2.
The aim of this study was to prove an extended version of the identity (1.6) and demonstrate its viability in approximating the integral of Bessel functions.
Gauss' hypergeometric function is defined by Andrews et al. (1999) 2.2where (u)k is the Pochhammer symbol (or rising factorial), given by The series (2.2) converges absolutely for |z|<1 (Andrews et al. 1999). If Re(c−a−b)>0, we have (Andrews et al. 1999) 2.3
Many special functions can be defined in terms of the hypergeometric function. In particular, the Gegenbauer (or ultraspherical) polynomials are defined by Erdélyi et al. (1981) 2.4with and The even Gegenbauer polynomials can also be expressed in terms of the hypergeometric function by Erdélyi et al. (1981) 2.5
3. Main result
The discontinuous integral was investigated by Weber (1873), Sonine (1880) and Schafheitlin (1887). They proved that (Watson 1995) 3.1for 3.2and 0<b<a. The corresponding expression for the case when 0<a<b is obtained from (3.1) by interchanging a,b and also μ,ν. When a=b, we have (Watson 1995) provided that Re(μ+ν+1)>Re(λ)>0. This result also follows from Gauss' summation formula (2.3) and (3.1).
We will now prove our main theorem.
If 0<b<a<π,Re(μ+ν−2k)>−1, and then 3.3
From Watson (1995), the Bessel function Jν(z) has the integral representation 3.4valid when and . Replacing ν by μ−2k, l by 2k, z by an in (3.4), and changing the integration variable to , we obtain Since (2.5) implies that is an even function, we conclude that 3.5valid when , where 3.6Replacing μ by ν, k by 0 in (3.4), and changing the integration variable from x to y, we have, because , that 3.7
Making the change of integration variable in (3.10), and setting ω=b/a<1, we obtain
The special case of theorem 3.1 in which k=0 was derived by Cooke (1928), as part of his work on Schlömilch series.
From (3.3), we have
If 0<a,b<π, and with then
4. Alternative approach
After posting our manuscript on ArXiv,1 we received several comments on the paper. Tom H. Koornwinder,2 observed that a more conceptual proof of corollary 3.2 is obtained by the following proposition. The same method was used by Boas and Pollard in Boas & Pollard (1973) to find integral representations of some series.
Suppose that f,g∈L2[α,β], with 4.1Then, where and
Suppose that f,g∈L2[α,β]. Then, Parseval's relation for complex Fourier series gives (Zayed 1996)
We can use the above proposition to prove the result in corollary 3.2. Using the formula (Gradshteyn & Ryzhik 2000) valid for we obtain, after some substitutions, where
If we define then Using the proposition with α=π,β=−π, we conclude that or 4.2with and 0<a,b≤π.
Corollary 3.2 follows from (4.2), but the disadvantage of this method is that it does not yield an expression for the series, as we obtained in theorem 3.1.
5. Numerical results
In previous sections, we have proved the equality of integrals and sums of Bessel functions. However, before the identity is used in practice, we need to consider its convergence behaviour.
Firstly, the Weierstrass M-test shows that the series converges uniformly as long as μ+ν−2k>1 and μ,ν>0 because the numerator in corollary 3.2 is bounded |Jμ(at)Jν(bt)|≤1 (Watson 1995).
To explore the rate of convergence of the sum in corollary 3.3, we choose μ=ν=5/2. The exact value of the integral is 5.1and truncation of the infinite series yields the finite sum 5.2and its error 5.3The integral in (5.1) and sum in (5.2) are illustrated in figure 1 and their difference (5.3) is shown in figure 2 and table 1.
The excellent agreement region in figure 1, an apparently flat plateau in figure 2, and the decaying error from the 2nd to the 7th column of table 1 strongly suggest that corollary 3.3 is true over the larger domain |a|+|b|<2π for ν=5/2. Additional numerical experiments not shown here suggest that it is true for all ν. Furthermore, we believe that this larger domain conjecture also applies to corollary 3.2. The extension of the 0<a,b≤π proof to other three quadrants is trivial but we have not yet managed to find a proof for the diamond domain |a|+|b|<2π.
Table 1 reveals that the rate of convergence of the Bessel sum is strongly dependent on the value of a and b. If a and b are in the domain, where the integral equals the infinite sum, the difference (5.3) for general μ,ν,k is only due to truncation and can be recast as The dependence of the truncation error follows directly from the asymptotic expansion for the Bessel function. From Watson (1995), the first two terms in that expansion are given by 5.4
For the parameters considered here, namely, , we then have, for the moment keeping only the first term in the expansion, With a=π/2, we then have in which Choosing N+1 to be odd, the truncation error is then For b=π/4, we have The alternate terms cancel, resulting in terms of order N−4. The sum is thus of order N−3.
In similar fashion, for b=3π/4, we have again resulting in terms of order N−4 after cancellation and a sum of order N−3. On the other hand, for b=a=π/2, We then have and thus the sum is of order N−1.
For b=π, the first term in the asymptotic expansion of Jν(πn) is zero, since for all n. The second term in the expansion then gives from which, with , Again, successive terms cancel, leaving terms of order N−4 and the sum is of order N−3. These explain the N−1 and N−3 orders found in table 1. In this case (μ=ν=5/2), we have found that the decay rate is also of order N−p, with 1≤p≤3, for other values of a,b not shown here. For general, parameters a,b,μ,ν and k, the truncation error may be obtained directly using (5.4).
6. Concluding remark
We provide a proof of a new identity involving integral and sum of Bessel functions and discuss the rate of convergence of the sum in the identity. Generalization of the identity and detailed investigation of its convergence rate should be explored in the future work.
The work of D. Dominici was supported by a Humboldt Research Fellowship for Experienced Researchers from the Alexander von Humboldt Foundation. P.M.W. Gill thanks the Australian Research Council (grants DP0984806 and DP1094170) for funding and the NCI National Facility for a generous grant of supercomputer time. T. Limpanuparb acknowledges the financial support from the Royal Thai Government (DPST programme). We wish to express our gratitude to Lawrence Glasser, Tom H. Koornwinder, Christophe Vignat and two anonymous referees who provided us with invaluable suggestions and comments that greatly improved this paper.
- Received November 6, 2011.
- Accepted March 22, 2012.
- This journal is © 2012 The Royal Society