## Abstract

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.

## 1. Introduction

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

**, 1.3**

*r**J*

_{l}(

*z*) is a Bessel function of the first kind (Olver

*et al.*2010), and

*r*=|

**|. The resolution (1.2) is computationally useful because it decouples the coordinates**

*r***and**

*r*

*r*^{′}and allows the two-body interaction integral 1.4between densities

*ρ*

_{a}(

**) and**

*r**ρ*

_{b}(

**) to be recast as 1.5where**

*r**A*

_{nlm}is a one-body integral of the product of

*ρ*

_{a}(

**) and**

*r**ϕ*

_{nlm}(

**). If the one-body integrals can be evaluated efficiently and the sum converges rapidly, (1.5) may offer a more efficient route to**

*r**E*[

*ρ*

_{a},

*ρ*

_{b}] than (1.4).

The key question is how best to obtain the *K*_{l} resolution
Previous attempts (Gilbert 1996; Varganov *et al.* 2008; Gill & Gilbert 2009; Limpanuparb & Gill 2009) yielded complicated *K*_{nl} whose practical utility is questionable but, recently, we have discovered the remarkable identity
1.6where *ε*_{n} is defined by
1.7*a*,*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 1948*a*,*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.

## 2. Preliminaries

The Bessel function of the first kind *J*_{ν}(*z*) is defined by Watson (1995)
2.1It follows from (2.1) that
is an entire function of *z* and we have

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.

### Theorem 3.1

*If 0<b<a<π,Re(μ+ν−2k)>−1, and* *then
*3.3

### Proof.

From Watson (1995), the Bessel function *J*_{ν}(*z*) has the integral representation
3.4valid when and . Replacing *ν* by *μ*−2*k*, *l* by 2*k*, *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

It follows from (3.5) and (3.7) that
and
where
Using the identity (Olver *et al.* 2010),
where *δ*(*z*) is the Dirac delta function, we obtain
Hence, since 0<*x*,*y*<*a*<*π*, we get
3.8

Using Euler's transformation (Olver *et al.* 2010)
in (2.5), we have
3.9

Using (3.9) in (3.8), we get 3.10

Making the change of integration variable in (3.10), and setting *ω*=*b*/*a*<1, we obtain

Using the formula (Olver *et al.* 2010),
valid for Re(*c*)>Re(*s*)>0, and , with and *c*=*ν*+1, we get
Using (3.6), we have
and (3.3) follows. ■

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.

### Corollary 3.2

*If* *then*

### Proof.

From (3.3), we have

Taking *λ*=*μ*+*ν*−2*k* in (3.1), we get
and the result follows. Note that since for all
and
the conditions (3.2) are satisfied. ■

### Corollary 3.3

*If* 0<*a*,*b*<*π*, *and* *with* *then*

### Proof.

The result is a consequence of corollary 3.2 (taking and a special case of the integral (3.1) (Watson 1995) ■

## 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.

### Proposition 4.1

*Suppose that* *f*,*g*∈*L*^{2}[*α*,*β*], *with*
4.1*Then*,
*where*
*and*

### Proof.

Suppose that *f*,*g*∈*L*^{2}[*α*,*β*]. Then, Parseval's relation for complex Fourier series gives (Zayed 1996)

On the other hand, we know from Parseval's relation for the Fourier transform (Zayed 1996) that Using (4.1), the result follows. ■

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 *μ*+*ν*−2*k*>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.

## Acknowledgements

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.

## Footnotes

- Received November 6, 2011.
- Accepted March 22, 2012.

- This journal is © 2012 The Royal Society