# A Review of  Procedures for Summing Kapteyn Series in Mathematical Physics

Authors: R. C. Tautz; I. Lerche
Journal: Advances in Mathematical Physics
 (2009)

---

## Abstract

Since the discussion of Kapteyn series occurrences in astronomical problems the wealth of
mathematical physics problems in which such series play dominant roles has burgeoned massively. One of the major concerns is the ability to sum such series in closed form so that one can better understand the structural and functional behavior of the basic physics
problems. The purpose of this review article is to present some of the recent methods for providing such series in closed form with applications to: (i) the summation of Kapteyn series for radiation from pulsars; (ii) the summation of other Kapteyn series in radiation problems; (iii) Kapteyn series arising in terahertz sideband spectra of quantum systems modulated by an alternating electromagnetic field; and (iv) some plasma problems involving sums of Bessel functions and their closed form summation using variations of the techniques developed for Kapteyn series. In addition, a short review is given of some other Kapteyn series to illustrate the ongoing deep interest and involvement of scientists in such problems and to provide further techniques for attempting to sum divers Kapteyn series.


---

## Body

## 1. Introduction

This review article is concerned with exhibiting techniques leading to either closed form expressions for Kapteyn series or integral representations that cannot be further  reduced.In general there are two sorts of Kapteyn series [1]. Kapteyn series of the first kind are infinite sums of Bessel functions of the form(1.1)F(x)=∑n=1∞fn(x)Jn(nx);
that is, Kapteyn series of the first kind involve summations over terms containing one Bessel function of the form Jn(nx), while Kapteyn series of the second kind involve terms each of which is proportional to a product of two such Bessel functions. Note that the index of summation n appears both in the order and in the argument of the Bessel functions.Kapteyn series arise in a host of mathematical physics problems. The range extends from pulsar physics [2, 3] through radiation from rings of discrete charges [4, 5] through quantum modulated systems [6, 7] through traffic queuing problems [8, 9] and on to plasma physics problems in ambient magnetic fields [10, 11] to name but a few such disciplines. Therefore, it seems appropriate to spell out a variety of techniques that can be used separately or in combination to sum such series efficiently.While some procedures for summation of selected Kapteyn series in mathematical physics have been known for over a century, the purpose here is to provide more general methods of broad use for many categories of such series. This purpose is based on many physical applications that have arisen over the last half century where, to date, either only asymptotic representations of the relevant Kapteyn series have been given or where recourse to direct numerical investigations have been given without considering whether closed form expressions exist at all for the series.In the latter situation it is difficult to determine whether the numerical methods provide accurate results because one has no basic template in closed form (or at worst in integral form) against which a comparison can be made. In the former case, while it is often that one can compare a known asymptotic representation of a Kapteyn series against numerical results, often one does not know the domain of validity of the asymptotic expansion nor does one know the functional behavior of the Kapteyn series in regions removed from the asymptotic result nor, indeed, does one have available the general domain of convergence of the desired Kapteyn series.For all of these reasons it is appropriate to review some general methods that can be used to sum a large array of Kapteyn series in mathematical physics.


## 2. Kapteyn Series in Pulsar Radiation Problems

In discussing radiation in vacuum from a rotating magnetic dipole, which is off-center with respect to a rotating pulsar, but which is “frozen” in the pulsar body, Harrison and Tademaru [2] showed that the total power radiated, L, is given by(2.1)L=Ω4c3∫0πdθsinθ[(μρ2+μϕ2cos2θ)a-2S1(a)+(μρ2cos2θ+μϕ2)S2(a)+μz2sin2θS1(a)],
and the force, F, acting on the dipole in the z direction is(2.2)F=2Ω4c3μzμϕ∫0πdθsin2θcos2θa-1S1(a),
where (μρ,μϕ,μz) are the magnetic dipole components in a cylindrical (ρ,ϕ,z) coordinate system (see Figure 1), Ω is the angular velocity of the pulsar, a=Ωs/csinθ with s being the offset distance of the dipole from the spin axis, and whereFigure 1
Sketch of the spin and dipole coordinates from Harrison and Tademaru [2].(2.3)S1(a)=∑n=1∞n4Jn2(na),(2.4)S2(a)=∑n=1∞n4Jn′2(na).Note thata∈(-1,1) is required so that the series S1 and S2 are convergent.Harrison and Tademaru [2] argued that for values of na≪1 one could approximate the power, L, and the force F as given in their Equations (5′) and (7′). However, the fact that n is in the range 1≤n≤∞ means that it is not easy to justify their expansion procedure. Further, in situations where a pulsar has a high spin rate and where the offset distance can approach the radius, R, of the pulsar, the factor Ωs/c is 𝒪(1) so that na≪1 is almost nowhere valid. To investigate such situations one needs closed form expressions for the two series S1 and S2. Watson [12] refers to these series as Kapteyn [1] series of the second kind, which series have been investigated to some extent by Nielsen [13].This section provides the general procedure for evaluating the series (2.3) and (2.4), although the method is of much greater generality as it will become clear in the course of its development.


### 2.1. Manipulations with the SeriesS1 and S2

First, differentiateS2 with respect to a to obtain(2.5)dS2da=2∑n=1∞n5Jn′(na)Jn′′(na).
Use Bessel's equation (e.g., [14, Section 9.1]) for Jn(na) in the form(2.6)Jn′′(na)=1a2[-anJn′(na)+(1-a2)Jn(na)]
in  (2.5) to obtain(2.7)a2dS2da+2aS2(a)=(1-a2)dS1da
so that(2.8)S2(a)=a-2[(1-a2)S1(a)+2∫0adxxS1(x)],
because the integration constant in  (2.8) is zero by evaluation as a→0.Thus, it is sufficient to evaluateS1 in closed form and to perform the integration in  (2.8) to obtain S2.Consider thenS1. Use the formula [15, Section 5.43](2.9)Jν(z)Jμ(z)=2π∫0π/2dψJν+μ(2zcosψ)cos[(μ-ν)θ]
in the form (e.g., [16, Section 6.681])(2.10)Jn2(na)=2π∫0π/2dψJ2n(2nacosψ)
so that(2.11)S1(a)=18π∫0π/2dψ∑n=1∞(2n)4J2n(2nacosψ).Expression (2.11) shows that S1 is expressed as an integral over a Kapteyn series of the first kind, for which several theorems are available as expressed in [15]. The most important result needed is the following.If the Kapteyn series(2.12)f(z)=∑m=1∞amJm(mz),
where am is arbitrary but given, is known in closed form, then the series(2.13)F(z)=∑m=1∞amm2Jm(mz)
is given by two simple integrations because(2.14)ℒzF(z)=f(z),
by direct differentiation of  (2.13). Again, z∈(-1,1) is required so that the series is convergent. Furthermore, the differential operator in  (2.14) with respect to an arbitrary variable is introduced as(2.15)ℒx=11-x2(xddx)2≡x1-x2ddx(xddx).Reversing the argument: ifF(z) is known, then f(z) is given directly by differentiation of F(z) in  (2.14).


### 2.2. Reduction ofS1 to Closed Form

From (2.11)–(2.14) we have that (with am=0 if m is odd, and am=m4 if m is even)(2.16)S1(a)=18π∫0π/2dψ(ℒb∘ℒb)[∑n=1∞J2n(2nb)],
where b=acosψ. But it is well known that [15, Section 17.33](2.17)11±z=1+2∑m=1∞(∓1)mJm(mz),
where z∈(-1,1), so that, also for b∈(-1,1),(2.18)∑n=1∞J2n(2nb)=b22(1-b2).Hence(2.19)S1(a)=18π∫0π/2dψ(ℒb∘ℒb)(b21-b2).Carrying out the differentiations in  (2.19) yields(2.20)S1(a)=1π∫0π/2dψb2(1-b2)7[1+14b2+21b4+4b6],
with b=acosψ.Using a partial fraction expansion and(2.21)∫0π/2dψ(1-b2)n=14(1-a2/2)n∫02πdu(1-ξcosψ)n≡Rn,
where ξ=a2/(2-a2), and using(2.22)R1=2π1-b2,Rn+1=Rn+ξn∂Rn∂ξ,
the integral in  (2.20) can be completed in closed form yielding(2.23)S1(a)=a2(64+592a2+472a4+27a6)256(1-a2)13/2.Inserting this expression forS1 into  (2.8) and performing the integral leads to(2.24)S2(a)=64+624a2+632a4+45a6256(1-a2)11/2.Thus, this procedure shows that bothS1 and S2 are available analytically. Numerical comparison of direct series evaluation term by term with the closed form analytical expressions confirms agreement to at least one part in 1016. The prototype of such Kapteyn series of the second kind was first given in closed form by Schott [17], who evaluated(2.25)∑n=1∞n2Jn2(na)=a2(4+a2)16(1-a2)7/2.
(Note that, in  (2.25), there appears a factor (1-a2)-7/2 which was missing in Lerche and Tautz [3].)The basic procedure for evaluating Kapteyn series of the generic form(2.26)∑n=1∞ann2mJn2(na),∑n=1∞ann2mJn'2(na),
where m is either an integer or half integer (positive or negative) and an is either unity or (-1)n, then follows the same recipe as given here, although the expressions rapidly become unwieldy as m becomes large.


### 2.3. Calculation ofL and F

The closed form expressions forS1 and S2 can then be used in (2.1) and (2.2) to evaluate the radiated power and force on the dipole in terms of ɛ≡Ωs/c. The result leads to elliptic integrals which cannot be solved analytically. But an expansion of the result in powers of ɛ yields L=L0+L1+⋯ and F=F0+F1+⋯, where(2.27)L0=2Ω43c3(μρ2+μϕ2+2Ω2s25c2μz2),L1=Ω6s215c5(94μρ2+92μϕ2+54Ω2s2c2μz2),F0=2Ω5s15c5μzμϕ,F1=6Ω7s35c7μzμϕ.Note that the zero-order termsL0 and F0 agree with the expansion given by Harrison and Tademaru [2] in their Equations (5′) and (7′).Next, the expression forL is separated as(2.28)L=Ω4c3(Lρμρ2+Lϕμϕ2+Lzμz2),
and the functions Lρ, Lϕ, and Lz as well as F are calculated numerically. Comparing the exact function values to the approximations from (2.27), drastic deviations are revealed even from the first-order approximations, as illustrated in Figures 2 and 3. The relative deviations are shown to reach 10% even for ɛ as low as 0.1 (zero-order approximation) and ~0.4 (first-order approximation).Comparison of the exact and approximate values for the three components of the functionL. In panels 1, 3, and 5, the solid lines show the (numerically calculated) exact functions Lρ, Lϕ, and Lz, respectively, and the dashed and dash-dot lines show the approximations F0 and F0+F1, respectively. All function values are normalized to Ω4/c3. In panels 2, 4, and 6, the relative deviation (in percent) from the exact function values is shown for the approximations L0 (solid lines) and L0+L1 (dashed lines), respectively. The two dotted lines mark deviations of 10% and 50%.
(a)(b)(c)(d)(e)(f)Comparison of the exact and approximate values for the functionF, normalized to Ω4c-4μzμϕ. In the upper panel, the solid line shows the (numerically calculated) exact function F, and the dashed and dash-dot lines show the approximations F0 and F0+F1, respectively. In the lower panel, the relative deviation (in percent) from the exact function values is shown for the approximations F0 (solid line) and F0+F1 (dashed line), respectively. The two dotted lines mark deviations of 10% and 50%.
(a)(b)The expansion parameterɛ, however, is normally very small as will be illustrated by two examples: (i) the fastest rotating pulsar [18] (PSR J1748–244ad) has a rotation period of 1/716 s with a radius of <8 km and, therefore, the offset of the dipole from the spin axis, s, would have to be as large as 6.7 km in order to have ɛ=0.1, with these parameters one would have obtained a deviation of 10% resulting from the approximations of Harrison and Tademaru [2]; (ii) for the Crab pulsar [19] (PSR B0531+21) the parameter yields ɛ=0.008s/R with R the pulsar radius, which is small even for large offsets s.If the surface velocity approached the speed of light, the expansion parameter would be given byɛ~s/R (without taking into account any relativistic effects); thus, ɛ can, at least in principle, attain values where both the zero-order and the first-order approximations from (2.27) become invalid.


## 2.1. Manipulations with the SeriesS1 and S2

First, differentiateS2 with respect to a to obtain(2.5)dS2da=2∑n=1∞n5Jn′(na)Jn′′(na).
Use Bessel's equation (e.g., [14, Section 9.1]) for Jn(na) in the form(2.6)Jn′′(na)=1a2[-anJn′(na)+(1-a2)Jn(na)]
in  (2.5) to obtain(2.7)a2dS2da+2aS2(a)=(1-a2)dS1da
so that(2.8)S2(a)=a-2[(1-a2)S1(a)+2∫0adxxS1(x)],
because the integration constant in  (2.8) is zero by evaluation as a→0.Thus, it is sufficient to evaluateS1 in closed form and to perform the integration in  (2.8) to obtain S2.Consider thenS1. Use the formula [15, Section 5.43](2.9)Jν(z)Jμ(z)=2π∫0π/2dψJν+μ(2zcosψ)cos[(μ-ν)θ]
in the form (e.g., [16, Section 6.681])(2.10)Jn2(na)=2π∫0π/2dψJ2n(2nacosψ)
so that(2.11)S1(a)=18π∫0π/2dψ∑n=1∞(2n)4J2n(2nacosψ).Expression (2.11) shows that S1 is expressed as an integral over a Kapteyn series of the first kind, for which several theorems are available as expressed in [15]. The most important result needed is the following.If the Kapteyn series(2.12)f(z)=∑m=1∞amJm(mz),
where am is arbitrary but given, is known in closed form, then the series(2.13)F(z)=∑m=1∞amm2Jm(mz)
is given by two simple integrations because(2.14)ℒzF(z)=f(z),
by direct differentiation of  (2.13). Again, z∈(-1,1) is required so that the series is convergent. Furthermore, the differential operator in  (2.14) with respect to an arbitrary variable is introduced as(2.15)ℒx=11-x2(xddx)2≡x1-x2ddx(xddx).Reversing the argument: ifF(z) is known, then f(z) is given directly by differentiation of F(z) in  (2.14).


## 2.2. Reduction ofS1 to Closed Form

From (2.11)–(2.14) we have that (with am=0 if m is odd, and am=m4 if m is even)(2.16)S1(a)=18π∫0π/2dψ(ℒb∘ℒb)[∑n=1∞J2n(2nb)],
where b=acosψ. But it is well known that [15, Section 17.33](2.17)11±z=1+2∑m=1∞(∓1)mJm(mz),
where z∈(-1,1), so that, also for b∈(-1,1),(2.18)∑n=1∞J2n(2nb)=b22(1-b2).Hence(2.19)S1(a)=18π∫0π/2dψ(ℒb∘ℒb)(b21-b2).Carrying out the differentiations in  (2.19) yields(2.20)S1(a)=1π∫0π/2dψb2(1-b2)7[1+14b2+21b4+4b6],
with b=acosψ.Using a partial fraction expansion and(2.21)∫0π/2dψ(1-b2)n=14(1-a2/2)n∫02πdu(1-ξcosψ)n≡Rn,
where ξ=a2/(2-a2), and using(2.22)R1=2π1-b2,Rn+1=Rn+ξn∂Rn∂ξ,
the integral in  (2.20) can be completed in closed form yielding(2.23)S1(a)=a2(64+592a2+472a4+27a6)256(1-a2)13/2.Inserting this expression forS1 into  (2.8) and performing the integral leads to(2.24)S2(a)=64+624a2+632a4+45a6256(1-a2)11/2.Thus, this procedure shows that bothS1 and S2 are available analytically. Numerical comparison of direct series evaluation term by term with the closed form analytical expressions confirms agreement to at least one part in 1016. The prototype of such Kapteyn series of the second kind was first given in closed form by Schott [17], who evaluated(2.25)∑n=1∞n2Jn2(na)=a2(4+a2)16(1-a2)7/2.
(Note that, in  (2.25), there appears a factor (1-a2)-7/2 which was missing in Lerche and Tautz [3].)The basic procedure for evaluating Kapteyn series of the generic form(2.26)∑n=1∞ann2mJn2(na),∑n=1∞ann2mJn'2(na),
where m is either an integer or half integer (positive or negative) and an is either unity or (-1)n, then follows the same recipe as given here, although the expressions rapidly become unwieldy as m becomes large.


## 2.3. Calculation ofL and F

The closed form expressions forS1 and S2 can then be used in (2.1) and (2.2) to evaluate the radiated power and force on the dipole in terms of ɛ≡Ωs/c. The result leads to elliptic integrals which cannot be solved analytically. But an expansion of the result in powers of ɛ yields L=L0+L1+⋯ and F=F0+F1+⋯, where(2.27)L0=2Ω43c3(μρ2+μϕ2+2Ω2s25c2μz2),L1=Ω6s215c5(94μρ2+92μϕ2+54Ω2s2c2μz2),F0=2Ω5s15c5μzμϕ,F1=6Ω7s35c7μzμϕ.Note that the zero-order termsL0 and F0 agree with the expansion given by Harrison and Tademaru [2] in their Equations (5′) and (7′).Next, the expression forL is separated as(2.28)L=Ω4c3(Lρμρ2+Lϕμϕ2+Lzμz2),
and the functions Lρ, Lϕ, and Lz as well as F are calculated numerically. Comparing the exact function values to the approximations from (2.27), drastic deviations are revealed even from the first-order approximations, as illustrated in Figures 2 and 3. The relative deviations are shown to reach 10% even for ɛ as low as 0.1 (zero-order approximation) and ~0.4 (first-order approximation).Comparison of the exact and approximate values for the three components of the functionL. In panels 1, 3, and 5, the solid lines show the (numerically calculated) exact functions Lρ, Lϕ, and Lz, respectively, and the dashed and dash-dot lines show the approximations F0 and F0+F1, respectively. All function values are normalized to Ω4/c3. In panels 2, 4, and 6, the relative deviation (in percent) from the exact function values is shown for the approximations L0 (solid lines) and L0+L1 (dashed lines), respectively. The two dotted lines mark deviations of 10% and 50%.
(a)(b)(c)(d)(e)(f)Comparison of the exact and approximate values for the functionF, normalized to Ω4c-4μzμϕ. In the upper panel, the solid line shows the (numerically calculated) exact function F, and the dashed and dash-dot lines show the approximations F0 and F0+F1, respectively. In the lower panel, the relative deviation (in percent) from the exact function values is shown for the approximations F0 (solid line) and F0+F1 (dashed line), respectively. The two dotted lines mark deviations of 10% and 50%.
(a)(b)The expansion parameterɛ, however, is normally very small as will be illustrated by two examples: (i) the fastest rotating pulsar [18] (PSR J1748–244ad) has a rotation period of 1/716 s with a radius of <8 km and, therefore, the offset of the dipole from the spin axis, s, would have to be as large as 6.7 km in order to have ɛ=0.1, with these parameters one would have obtained a deviation of 10% resulting from the approximations of Harrison and Tademaru [2]; (ii) for the Crab pulsar [19] (PSR B0531+21) the parameter yields ɛ=0.008s/R with R the pulsar radius, which is small even for large offsets s.If the surface velocity approached the speed of light, the expansion parameter would be given byɛ~s/R (without taking into account any relativistic effects); thus, ɛ can, at least in principle, attain values where both the zero-order and the first-order approximations from (2.27) become invalid.


## 3. Kapteyn Series in Other Radiation Problems

One problem in radiation that was considered of great interest at the beginning of the 20th Century is the following. It is well known that a single point charge, moving uniformly in a circle, radiates. Suppose then that one hasN charges equally spaced around a circle and all moving at the same circular speed. Then they, too, radiate. Now as the number N of charges is increased, all other conditions being held fixed, then the spacing between charges decreases proportional to 1/N. The limit of this process is a continuous uniform charge distribution moving with constant circular motion, that is, a steady-state ring current. But it is also well known that such a current formation does not radiate. Then the question is as N→∞ how does the radiation diminish so that, finally, there is no radiation from a continuous ring current?Investigations of this basic problem immediately encountered Kapteyn series of the second kind (see, e.g., [1, 15]) in a variety of forms and guises. While the formula describing the radiation output was expressible as a set of terms involving sums of Kapteyn series, at first only approximations to the series could be obtained for arbitrary N [4]. The work of Budden [20] provided a systematic determination of the Kapteyn series involved and evaluated the radiation field of the N like particles in terms of factors summed to N/2-1. The advantage was that, along the way, Budden managed to effect solutions in closed analytical form to some of the Kapteyn series involved. The upshot was that, as N→∞, one could show how the radiation field diminished to zero.Since that time there has been, and continues to be, interest in a variety of such radiation types of problems. Alternating positive and negative point charges spread uniformly around a ring, each of which moves at constant circular speed, is one such problem [17]. As the number of charges increases without limit the spacing between successive charges tends to zero so that, in the limit, there is a charge neutral ring that does not radiate. The approach of the radiation field to zero as the number of charges tends to infinity is the problem of interest. Fortunately this problem is just a variant of the problem solved by Budden [20] because it represents two rings of opposite charges with twice the spacing. Budden's solution is then immediately appropriate by superposition and charge reversal.Radiation from a magnetic dipole, off-center from a pulsar that spins, is another such problem, as we have seen earlier in this review [2, 3], as is the radiation field from a charged particle undergoing elliptical motion [21].In all such problems there have arisen, to date, twelve basic Kapteyn series of the second kind, some of which have been known in closed form for a while while others are often referred to as “solved” but seem to be not readily available, if at all.The next section of the review provides the basic methodology to handle all twelve of the series and shows which are expressible in closed analytic form, and which are only expressible only as integrals that cannot be reduced to analytic form.


### 3.1. Manipulations with Basic Sets of Kapteyn Series

#### 3.1.1. The Sets of Series

The twelve series in question are given by(3.1)S1(λ,m,b)=∑n=1∞λnn2mJn2(nb),S2(λ,m,b)=∑n=1∞λnn2m+1Jn2(nb),S3(λ,m,b)=∑n=1∞λnn2mJn′2(nb),S4(λ,m,b)=∑n=1∞λnn2m+1Jn′2(nb),S5(λ,m,b)=∑n=1∞λnn2mJn(nb)Jn′(nb),S6(λ,m,b)=∑n=1∞λnn2m+1Jn(nb)Jn′(nb),
where λ∈{±1} and m∈Z.Determination of the sets of series can be reduced to the simpler problem of determining only the set of series withm=0 (in the cases of S1, S3, and S6) and the set of series with m=-1 (in the cases of S2, S4, and S5).The reason for these reductions is as follows. One can write(3.2)2S6(λ,m,b)=∂S1∂b,2S5(λ,m,b)=∂S2∂b
so that it is sufficient to obtain S1, S2, S3, and S4.Note also that(3.3)∂S3∂b=2∑n=1∞λnn2m+1Jn′(nb)Jn′′(nb).But, because of Bessel's equation (see  (2.6) with a changed to b), one has(3.4)b2∂S3∂b+2bS3=(1-b2)∂S1∂b
so that(3.5)S3(λ,m,b)=1b2[(1-b2)S1(λ,m,b)+2∫0bdxS1(λ,m,x)].
Equally(3.6)S4(λ,m,b)=1b2[(1-b2)S2(λ,m,b)+2∫0bdxS2(λ,m,b)].Thus it is sufficient to obtainS1 and S2.One can also use the theorem due to Watson [12] of  (2.14), which was derived in Section 2.1, and which yields f(b) if g(b) is known. Alternatively, if f(b) is known then g(b) is given by direct differentiation.Consider thenS1. Use  (2.10) so that(3.7)S1(λ,m,b)=2π∫0π/2dψ∑n=1∞λnn2mJ2n(2nbcosψ).But the series(3.8)hm(b)=∑n=1∞λnn2mJ2n(2nbcosψ)≡122m∑n=1∞λn(2n)2mJ2n(2nbcosψ)
is precisely of the form required in Watson's theorem, with an=0 if n is odd and an=exp[inπ/2lnλ]n2m if n is even, so that(3.9)hm(b)=ℒbhm-2(b),
where the differential operator ℒ from  (2.15) has been used. Hence, for m>0 all series of the type S1 can be reduced to the determination of h0(b) by differentiation. Equally, for m<0 one can use Watson's theorem in the converse sense to note that(3.10)h-|m|(b)=ℒbh-|m|+2(b)
so that, by two integrations, one has a recursive relation leading directly to h0.Thus, all twelve of the basic series needed can be written in terms of four fundamental series(3.11)F(λ,b)=∑n=1∞λnnJn2(nb),G(λ,b)=∑n=1∞λnJn2(nb)
for λ∈{±1}. All other series (with m≠0, or m≠-1, resp.) are directly given as simple differentials or simple integrals with respect to b of one or the other of the four fundamental series. It is, therefore, both necessary and sufficient to consider F and G.


#### 3.1.2. The Two Series Represented byF

Set(3.12a)F+(b)=∑n=1∞Jn2(nb)n,(3.12b)F-(b)=∑n=1∞(-1)nJn2(nb)n.Now, inF+, replace the Bessel functions using again  (2.10) while in F- replace (-1)nJn2(nb)=Jn(nb)J-n(nb) and(3.13)Jn(nb)J-n(nb)=2π∫0π/2dψJ0(2nbcosψ)cos2nψ.Then write(3.14)J2n(2nbcosψ)=2π∫0π/2dθcos(2nbcosψsinθ)cos2nθ,(3.15)J0(2nbcosψ)=2π∫0π/2dθcos(2nbcosψsinθ)
(see [15]).In principle, one could also use a representation of the Bessel function in exponential form [16] see and then carry out the summation. However, because (3.12a) and (3.12b) are a product of two Bessel functions, this ansatz would be even more difficult than the approach followed here.Now, inserting (2.10) and (3.14) into expression (3.12a) for F+ and inserting (3.13) and (3.15) into expression (3.12b) for F- and then performing directly the infinite sums lead, after some tedious but elementary algebra, to(3.16a)F+(b)=-1π2∫0π/2dϕ∫0π/2dθln[sin2(θ-bcosϕsinθ)sin2(θ+bcosϕsinθ)sin4θ],(3.16b)F-(b)=-1π2∫0π/2dϕ∫0π/2dθln[cos2(θ-bcosϕsinθ)cos2(θ+bcosϕsinθ)cos4θ],=-1π2∫0π/2dϕ∫0π/2dθln[sin2(θ-bcosϕcosθ)sin2(θ+bcosϕcosθ)sin4θ].Numerical investigation by direct summation ofF+ and F- as given in (3.12a), (3.12b) and comparison with the simple integral formulations given in (3.16a), (3.16b) shows that the series are indeed given by (3.16a), (3.16b) to better than a part in 104; this limit on resolution being caused by numerical round-off error. Figures 4 and 5 show the comparison between the integrals and direct summation as a function of increasing b∈(0,1) for both F+ and F-, respectively, with the relative error (in percent) also being plotted. (Note that, for numerical reasons, the relative error increases above 10-4 percent as b→1 (Figure 4) and as b→0 (Figure 5), respectively. Such depends heavily on the numerical summation and integration methods as well as on the computer time. By expansion of the integrals around b=1 and b=0, however, one can get almost exact agreement of the series and the integral.)The seriesF+ from  (3.12a) with the relative error when compared to the integral from  (3.16a).
(a)(b)The seriesF- from  (3.12b) with the relative error when compared to the integral from  (3.16b).
(a)(b)Throughout this review, the numerical evaluation of infinite sums is carried out as follows: First, a number of terms (usually1000) is summed directly; to accelerate the convergence of the sum, then Wynn's epsilon method (see, e.g., [22, 23]) is used, which samples a number of additional terms (usually 100) in the sum, and then tries to fit them to a polynomial multiplied by a decaying exponential. Thus, the series are well approximated and the required computer time is kept moderate. The convergence of the sums, in addition, is guaranteed by analytical considerations. Furthermore, numerical integrations are carried out using standard techniques such as adaptive grids. However, some care has to be taken of the square-root singularity (e.g., at ϕ=θ=0 in (3.16a) and (3.16b)). Since we used Mathematica version 6.0, this problem is dealt with automatically. Using other packages, however, appropriate measures would have to be taken manually.Marshall [21] suggested that the sum FM≡(1/2)∂F+/∂b, written in the form(3.17)FM(b)=∑n=1∞Jn(nb)Jn′(nb),
could be represented by a single elliptic integral (his  (2.22)) as(3.18)GM(b)=1πb∫1∞du(uu2-b2sin2u-1).Figure6 shows plots (as a function of b) of both the sum FM and the elliptic integral representation, GM from  (3.18), suggested in [21]. There is no agreement even at the crudest level of approximation that indicating the elliptic integral is not appropriate.The seriesFM from  (3.17) (solid line) compared to the integral representation GM from  (3.18), as given in Marshall [21] (dashed line). In the lower panel, the relative error with respect to the direct summation of the series is shown.
(a)(b)


#### 3.1.3. The Two Series Represented byG

Set(3.19a)G+(b)=∑n=1∞Jn2(nb),(3.19b)G-(b)=∑n=1∞(-1)nJn2(nb).The seriesG+ has been known in closed form since the time of Schott [17]. Use the well-known fact [1] that(3.20)11-bcosϕ=1+2∑n=1∞Jn(nb)cos[n(ϕ-bsinϕ)].Integrate  (3.20) over 0⩽ϕ⩽π, thereby obtaining(3.21)∑n=1∞Jn(nb)2=12(11-b2-1),
which is just Schott's [17] formula.The seriesG- is considerably more complicated to evaluate. Write(3.22)G-(b)≡∑n=1∞Jn(nb)J-n(nb)=2π∫0π/2dψ∑n=1∞J0(2nbcosψ)cos2nψ.Now use the Schlömilch [24] formula, which states that any function(3.23a)γ(x)=2π∫0π2dϕΓ(xsinϕ),
which is given through an arbitrary (but known) function Γ, can be rewritten as(3.23b)γ(x)=1π∫0πduΓ(u)+2π∫0πduΓ(u)∑n=1∞J0(nx)cosnu.
Set Γ(u)=δ(u-w) with 0⩽w⩽π so that(3.24)∑n=1∞J0(nx)cosnw=12[πγ(x)-1].With the identificationsw=2ψ and x=2bcosψ,  (3.22) then yields(3.25a)G-(b)=-12+1π∫0ψ⋆dψb2cos2ψ-ψ2,
where the upper integration limit is implicitly given by ψ⋆=bcosψ⋆, or b is given explicitly by b=ψ⋆secψ⋆. One can then write(3.25b)G-(b)=-12+cosψ⋆π∫01dzcos2(ψ⋆z)-z2cos2ψ⋆,
which might be more amenable when numerical integration is required. Figure 7 compares G- given by  (3.25a) with direct term by term summation of the series in  (3.19b), showing that, to within about 1 part in 105, the two are identical in the interval 0<b<1 (cf. footnote). Note also that the integral representation of G- is convergent for all values of b, including b>1.The values forψ⋆ as a function of b (a) and the series G- from  (3.19b) together with the relative error when compared to the integral from  (3.25a) ((b), (c)).
(a)(b)(c)


### 3.2. Discussion

A general method has been presented for the evaluation of twelve Kapteyn series of the second kind. Such series are important for the analytic description of radiation processes in various astrophysical applications such as the radiation from off-centered dipoles in neutron stars. Originally, the Kapteyn series described here arose when the attempt was made to describe the radiation from a distribution of a finite number of discrete point charges, all moving at uniform spacing at constant speed in a circle.Previously, most of the Kapteyn series have not been evaluated or, in the case of one of the series, were written in terms of a single elliptic integral, which turned out to be invalid when evaluated numerically (see  (3.18)).  Equation (3.25a) is more appropriate because it represents the series G- in terms of a different, but also elliptic, integral.As has been shown here by recurrence relations, there are only four basic series that need to be calculated, one of which was already known in closed algebraic form. All other of the twelve series can be obtained from direct differentiation or integration of one or other of the four basic series. The series can be evaluated in terms of closed analytic expressions or in terms of integrals that cannot be further reduced. Numerical calculations were carried out to compare the values obtained by direct summation to those obtained from the integral representations, and the relative errors (less than a part in104) were shown to be limited by numerical round-off errors that are responsible for the differences occurring between direct series representations and integral representations of the series.Furthermore, the method presented here may be useful when one has other Kapteyn series of the second kind to consider, thereby providing an additional reason to consider such series anew.


## 3.1. Manipulations with Basic Sets of Kapteyn Series

### 3.1.1. The Sets of Series

The twelve series in question are given by(3.1)S1(λ,m,b)=∑n=1∞λnn2mJn2(nb),S2(λ,m,b)=∑n=1∞λnn2m+1Jn2(nb),S3(λ,m,b)=∑n=1∞λnn2mJn′2(nb),S4(λ,m,b)=∑n=1∞λnn2m+1Jn′2(nb),S5(λ,m,b)=∑n=1∞λnn2mJn(nb)Jn′(nb),S6(λ,m,b)=∑n=1∞λnn2m+1Jn(nb)Jn′(nb),
where λ∈{±1} and m∈Z.Determination of the sets of series can be reduced to the simpler problem of determining only the set of series withm=0 (in the cases of S1, S3, and S6) and the set of series with m=-1 (in the cases of S2, S4, and S5).The reason for these reductions is as follows. One can write(3.2)2S6(λ,m,b)=∂S1∂b,2S5(λ,m,b)=∂S2∂b
so that it is sufficient to obtain S1, S2, S3, and S4.Note also that(3.3)∂S3∂b=2∑n=1∞λnn2m+1Jn′(nb)Jn′′(nb).But, because of Bessel's equation (see  (2.6) with a changed to b), one has(3.4)b2∂S3∂b+2bS3=(1-b2)∂S1∂b
so that(3.5)S3(λ,m,b)=1b2[(1-b2)S1(λ,m,b)+2∫0bdxS1(λ,m,x)].
Equally(3.6)S4(λ,m,b)=1b2[(1-b2)S2(λ,m,b)+2∫0bdxS2(λ,m,b)].Thus it is sufficient to obtainS1 and S2.One can also use the theorem due to Watson [12] of  (2.14), which was derived in Section 2.1, and which yields f(b) if g(b) is known. Alternatively, if f(b) is known then g(b) is given by direct differentiation.Consider thenS1. Use  (2.10) so that(3.7)S1(λ,m,b)=2π∫0π/2dψ∑n=1∞λnn2mJ2n(2nbcosψ).But the series(3.8)hm(b)=∑n=1∞λnn2mJ2n(2nbcosψ)≡122m∑n=1∞λn(2n)2mJ2n(2nbcosψ)
is precisely of the form required in Watson's theorem, with an=0 if n is odd and an=exp[inπ/2lnλ]n2m if n is even, so that(3.9)hm(b)=ℒbhm-2(b),
where the differential operator ℒ from  (2.15) has been used. Hence, for m>0 all series of the type S1 can be reduced to the determination of h0(b) by differentiation. Equally, for m<0 one can use Watson's theorem in the converse sense to note that(3.10)h-|m|(b)=ℒbh-|m|+2(b)
so that, by two integrations, one has a recursive relation leading directly to h0.Thus, all twelve of the basic series needed can be written in terms of four fundamental series(3.11)F(λ,b)=∑n=1∞λnnJn2(nb),G(λ,b)=∑n=1∞λnJn2(nb)
for λ∈{±1}. All other series (with m≠0, or m≠-1, resp.) are directly given as simple differentials or simple integrals with respect to b of one or the other of the four fundamental series. It is, therefore, both necessary and sufficient to consider F and G.


### 3.1.2. The Two Series Represented byF

Set(3.12a)F+(b)=∑n=1∞Jn2(nb)n,(3.12b)F-(b)=∑n=1∞(-1)nJn2(nb)n.Now, inF+, replace the Bessel functions using again  (2.10) while in F- replace (-1)nJn2(nb)=Jn(nb)J-n(nb) and(3.13)Jn(nb)J-n(nb)=2π∫0π/2dψJ0(2nbcosψ)cos2nψ.Then write(3.14)J2n(2nbcosψ)=2π∫0π/2dθcos(2nbcosψsinθ)cos2nθ,(3.15)J0(2nbcosψ)=2π∫0π/2dθcos(2nbcosψsinθ)
(see [15]).In principle, one could also use a representation of the Bessel function in exponential form [16] see and then carry out the summation. However, because (3.12a) and (3.12b) are a product of two Bessel functions, this ansatz would be even more difficult than the approach followed here.Now, inserting (2.10) and (3.14) into expression (3.12a) for F+ and inserting (3.13) and (3.15) into expression (3.12b) for F- and then performing directly the infinite sums lead, after some tedious but elementary algebra, to(3.16a)F+(b)=-1π2∫0π/2dϕ∫0π/2dθln[sin2(θ-bcosϕsinθ)sin2(θ+bcosϕsinθ)sin4θ],(3.16b)F-(b)=-1π2∫0π/2dϕ∫0π/2dθln[cos2(θ-bcosϕsinθ)cos2(θ+bcosϕsinθ)cos4θ],=-1π2∫0π/2dϕ∫0π/2dθln[sin2(θ-bcosϕcosθ)sin2(θ+bcosϕcosθ)sin4θ].Numerical investigation by direct summation ofF+ and F- as given in (3.12a), (3.12b) and comparison with the simple integral formulations given in (3.16a), (3.16b) shows that the series are indeed given by (3.16a), (3.16b) to better than a part in 104; this limit on resolution being caused by numerical round-off error. Figures 4 and 5 show the comparison between the integrals and direct summation as a function of increasing b∈(0,1) for both F+ and F-, respectively, with the relative error (in percent) also being plotted. (Note that, for numerical reasons, the relative error increases above 10-4 percent as b→1 (Figure 4) and as b→0 (Figure 5), respectively. Such depends heavily on the numerical summation and integration methods as well as on the computer time. By expansion of the integrals around b=1 and b=0, however, one can get almost exact agreement of the series and the integral.)The seriesF+ from  (3.12a) with the relative error when compared to the integral from  (3.16a).
(a)(b)The seriesF- from  (3.12b) with the relative error when compared to the integral from  (3.16b).
(a)(b)Throughout this review, the numerical evaluation of infinite sums is carried out as follows: First, a number of terms (usually1000) is summed directly; to accelerate the convergence of the sum, then Wynn's epsilon method (see, e.g., [22, 23]) is used, which samples a number of additional terms (usually 100) in the sum, and then tries to fit them to a polynomial multiplied by a decaying exponential. Thus, the series are well approximated and the required computer time is kept moderate. The convergence of the sums, in addition, is guaranteed by analytical considerations. Furthermore, numerical integrations are carried out using standard techniques such as adaptive grids. However, some care has to be taken of the square-root singularity (e.g., at ϕ=θ=0 in (3.16a) and (3.16b)). Since we used Mathematica version 6.0, this problem is dealt with automatically. Using other packages, however, appropriate measures would have to be taken manually.Marshall [21] suggested that the sum FM≡(1/2)∂F+/∂b, written in the form(3.17)FM(b)=∑n=1∞Jn(nb)Jn′(nb),
could be represented by a single elliptic integral (his  (2.22)) as(3.18)GM(b)=1πb∫1∞du(uu2-b2sin2u-1).Figure6 shows plots (as a function of b) of both the sum FM and the elliptic integral representation, GM from  (3.18), suggested in [21]. There is no agreement even at the crudest level of approximation that indicating the elliptic integral is not appropriate.The seriesFM from  (3.17) (solid line) compared to the integral representation GM from  (3.18), as given in Marshall [21] (dashed line). In the lower panel, the relative error with respect to the direct summation of the series is shown.
(a)(b)


### 3.1.3. The Two Series Represented byG

Set(3.19a)G+(b)=∑n=1∞Jn2(nb),(3.19b)G-(b)=∑n=1∞(-1)nJn2(nb).The seriesG+ has been known in closed form since the time of Schott [17]. Use the well-known fact [1] that(3.20)11-bcosϕ=1+2∑n=1∞Jn(nb)cos[n(ϕ-bsinϕ)].Integrate  (3.20) over 0⩽ϕ⩽π, thereby obtaining(3.21)∑n=1∞Jn(nb)2=12(11-b2-1),
which is just Schott's [17] formula.The seriesG- is considerably more complicated to evaluate. Write(3.22)G-(b)≡∑n=1∞Jn(nb)J-n(nb)=2π∫0π/2dψ∑n=1∞J0(2nbcosψ)cos2nψ.Now use the Schlömilch [24] formula, which states that any function(3.23a)γ(x)=2π∫0π2dϕΓ(xsinϕ),
which is given through an arbitrary (but known) function Γ, can be rewritten as(3.23b)γ(x)=1π∫0πduΓ(u)+2π∫0πduΓ(u)∑n=1∞J0(nx)cosnu.
Set Γ(u)=δ(u-w) with 0⩽w⩽π so that(3.24)∑n=1∞J0(nx)cosnw=12[πγ(x)-1].With the identificationsw=2ψ and x=2bcosψ,  (3.22) then yields(3.25a)G-(b)=-12+1π∫0ψ⋆dψb2cos2ψ-ψ2,
where the upper integration limit is implicitly given by ψ⋆=bcosψ⋆, or b is given explicitly by b=ψ⋆secψ⋆. One can then write(3.25b)G-(b)=-12+cosψ⋆π∫01dzcos2(ψ⋆z)-z2cos2ψ⋆,
which might be more amenable when numerical integration is required. Figure 7 compares G- given by  (3.25a) with direct term by term summation of the series in  (3.19b), showing that, to within about 1 part in 105, the two are identical in the interval 0<b<1 (cf. footnote). Note also that the integral representation of G- is convergent for all values of b, including b>1.The values forψ⋆ as a function of b (a) and the series G- from  (3.19b) together with the relative error when compared to the integral from  (3.25a) ((b), (c)).
(a)(b)(c)


## 3.1.1. The Sets of Series

The twelve series in question are given by(3.1)S1(λ,m,b)=∑n=1∞λnn2mJn2(nb),S2(λ,m,b)=∑n=1∞λnn2m+1Jn2(nb),S3(λ,m,b)=∑n=1∞λnn2mJn′2(nb),S4(λ,m,b)=∑n=1∞λnn2m+1Jn′2(nb),S5(λ,m,b)=∑n=1∞λnn2mJn(nb)Jn′(nb),S6(λ,m,b)=∑n=1∞λnn2m+1Jn(nb)Jn′(nb),
where λ∈{±1} and m∈Z.Determination of the sets of series can be reduced to the simpler problem of determining only the set of series withm=0 (in the cases of S1, S3, and S6) and the set of series with m=-1 (in the cases of S2, S4, and S5).The reason for these reductions is as follows. One can write(3.2)2S6(λ,m,b)=∂S1∂b,2S5(λ,m,b)=∂S2∂b
so that it is sufficient to obtain S1, S2, S3, and S4.Note also that(3.3)∂S3∂b=2∑n=1∞λnn2m+1Jn′(nb)Jn′′(nb).But, because of Bessel's equation (see  (2.6) with a changed to b), one has(3.4)b2∂S3∂b+2bS3=(1-b2)∂S1∂b
so that(3.5)S3(λ,m,b)=1b2[(1-b2)S1(λ,m,b)+2∫0bdxS1(λ,m,x)].
Equally(3.6)S4(λ,m,b)=1b2[(1-b2)S2(λ,m,b)+2∫0bdxS2(λ,m,b)].Thus it is sufficient to obtainS1 and S2.One can also use the theorem due to Watson [12] of  (2.14), which was derived in Section 2.1, and which yields f(b) if g(b) is known. Alternatively, if f(b) is known then g(b) is given by direct differentiation.Consider thenS1. Use  (2.10) so that(3.7)S1(λ,m,b)=2π∫0π/2dψ∑n=1∞λnn2mJ2n(2nbcosψ).But the series(3.8)hm(b)=∑n=1∞λnn2mJ2n(2nbcosψ)≡122m∑n=1∞λn(2n)2mJ2n(2nbcosψ)
is precisely of the form required in Watson's theorem, with an=0 if n is odd and an=exp[inπ/2lnλ]n2m if n is even, so that(3.9)hm(b)=ℒbhm-2(b),
where the differential operator ℒ from  (2.15) has been used. Hence, for m>0 all series of the type S1 can be reduced to the determination of h0(b) by differentiation. Equally, for m<0 one can use Watson's theorem in the converse sense to note that(3.10)h-|m|(b)=ℒbh-|m|+2(b)
so that, by two integrations, one has a recursive relation leading directly to h0.Thus, all twelve of the basic series needed can be written in terms of four fundamental series(3.11)F(λ,b)=∑n=1∞λnnJn2(nb),G(λ,b)=∑n=1∞λnJn2(nb)
for λ∈{±1}. All other series (with m≠0, or m≠-1, resp.) are directly given as simple differentials or simple integrals with respect to b of one or the other of the four fundamental series. It is, therefore, both necessary and sufficient to consider F and G.


## 3.1.2. The Two Series Represented byF

Set(3.12a)F+(b)=∑n=1∞Jn2(nb)n,(3.12b)F-(b)=∑n=1∞(-1)nJn2(nb)n.Now, inF+, replace the Bessel functions using again  (2.10) while in F- replace (-1)nJn2(nb)=Jn(nb)J-n(nb) and(3.13)Jn(nb)J-n(nb)=2π∫0π/2dψJ0(2nbcosψ)cos2nψ.Then write(3.14)J2n(2nbcosψ)=2π∫0π/2dθcos(2nbcosψsinθ)cos2nθ,(3.15)J0(2nbcosψ)=2π∫0π/2dθcos(2nbcosψsinθ)
(see [15]).In principle, one could also use a representation of the Bessel function in exponential form [16] see and then carry out the summation. However, because (3.12a) and (3.12b) are a product of two Bessel functions, this ansatz would be even more difficult than the approach followed here.Now, inserting (2.10) and (3.14) into expression (3.12a) for F+ and inserting (3.13) and (3.15) into expression (3.12b) for F- and then performing directly the infinite sums lead, after some tedious but elementary algebra, to(3.16a)F+(b)=-1π2∫0π/2dϕ∫0π/2dθln[sin2(θ-bcosϕsinθ)sin2(θ+bcosϕsinθ)sin4θ],(3.16b)F-(b)=-1π2∫0π/2dϕ∫0π/2dθln[cos2(θ-bcosϕsinθ)cos2(θ+bcosϕsinθ)cos4θ],=-1π2∫0π/2dϕ∫0π/2dθln[sin2(θ-bcosϕcosθ)sin2(θ+bcosϕcosθ)sin4θ].Numerical investigation by direct summation ofF+ and F- as given in (3.12a), (3.12b) and comparison with the simple integral formulations given in (3.16a), (3.16b) shows that the series are indeed given by (3.16a), (3.16b) to better than a part in 104; this limit on resolution being caused by numerical round-off error. Figures 4 and 5 show the comparison between the integrals and direct summation as a function of increasing b∈(0,1) for both F+ and F-, respectively, with the relative error (in percent) also being plotted. (Note that, for numerical reasons, the relative error increases above 10-4 percent as b→1 (Figure 4) and as b→0 (Figure 5), respectively. Such depends heavily on the numerical summation and integration methods as well as on the computer time. By expansion of the integrals around b=1 and b=0, however, one can get almost exact agreement of the series and the integral.)The seriesF+ from  (3.12a) with the relative error when compared to the integral from  (3.16a).
(a)(b)The seriesF- from  (3.12b) with the relative error when compared to the integral from  (3.16b).
(a)(b)Throughout this review, the numerical evaluation of infinite sums is carried out as follows: First, a number of terms (usually1000) is summed directly; to accelerate the convergence of the sum, then Wynn's epsilon method (see, e.g., [22, 23]) is used, which samples a number of additional terms (usually 100) in the sum, and then tries to fit them to a polynomial multiplied by a decaying exponential. Thus, the series are well approximated and the required computer time is kept moderate. The convergence of the sums, in addition, is guaranteed by analytical considerations. Furthermore, numerical integrations are carried out using standard techniques such as adaptive grids. However, some care has to be taken of the square-root singularity (e.g., at ϕ=θ=0 in (3.16a) and (3.16b)). Since we used Mathematica version 6.0, this problem is dealt with automatically. Using other packages, however, appropriate measures would have to be taken manually.Marshall [21] suggested that the sum FM≡(1/2)∂F+/∂b, written in the form(3.17)FM(b)=∑n=1∞Jn(nb)Jn′(nb),
could be represented by a single elliptic integral (his  (2.22)) as(3.18)GM(b)=1πb∫1∞du(uu2-b2sin2u-1).Figure6 shows plots (as a function of b) of both the sum FM and the elliptic integral representation, GM from  (3.18), suggested in [21]. There is no agreement even at the crudest level of approximation that indicating the elliptic integral is not appropriate.The seriesFM from  (3.17) (solid line) compared to the integral representation GM from  (3.18), as given in Marshall [21] (dashed line). In the lower panel, the relative error with respect to the direct summation of the series is shown.
(a)(b)


## 3.1.3. The Two Series Represented byG

Set(3.19a)G+(b)=∑n=1∞Jn2(nb),(3.19b)G-(b)=∑n=1∞(-1)nJn2(nb).The seriesG+ has been known in closed form since the time of Schott [17]. Use the well-known fact [1] that(3.20)11-bcosϕ=1+2∑n=1∞Jn(nb)cos[n(ϕ-bsinϕ)].Integrate  (3.20) over 0⩽ϕ⩽π, thereby obtaining(3.21)∑n=1∞Jn(nb)2=12(11-b2-1),
which is just Schott's [17] formula.The seriesG- is considerably more complicated to evaluate. Write(3.22)G-(b)≡∑n=1∞Jn(nb)J-n(nb)=2π∫0π/2dψ∑n=1∞J0(2nbcosψ)cos2nψ.Now use the Schlömilch [24] formula, which states that any function(3.23a)γ(x)=2π∫0π2dϕΓ(xsinϕ),
which is given through an arbitrary (but known) function Γ, can be rewritten as(3.23b)γ(x)=1π∫0πduΓ(u)+2π∫0πduΓ(u)∑n=1∞J0(nx)cosnu.
Set Γ(u)=δ(u-w) with 0⩽w⩽π so that(3.24)∑n=1∞J0(nx)cosnw=12[πγ(x)-1].With the identificationsw=2ψ and x=2bcosψ,  (3.22) then yields(3.25a)G-(b)=-12+1π∫0ψ⋆dψb2cos2ψ-ψ2,
where the upper integration limit is implicitly given by ψ⋆=bcosψ⋆, or b is given explicitly by b=ψ⋆secψ⋆. One can then write(3.25b)G-(b)=-12+cosψ⋆π∫01dzcos2(ψ⋆z)-z2cos2ψ⋆,
which might be more amenable when numerical integration is required. Figure 7 compares G- given by  (3.25a) with direct term by term summation of the series in  (3.19b), showing that, to within about 1 part in 105, the two are identical in the interval 0<b<1 (cf. footnote). Note also that the integral representation of G- is convergent for all values of b, including b>1.The values forψ⋆ as a function of b (a) and the series G- from  (3.19b) together with the relative error when compared to the integral from  (3.25a) ((b), (c)).
(a)(b)(c)


## 3.2. Discussion

A general method has been presented for the evaluation of twelve Kapteyn series of the second kind. Such series are important for the analytic description of radiation processes in various astrophysical applications such as the radiation from off-centered dipoles in neutron stars. Originally, the Kapteyn series described here arose when the attempt was made to describe the radiation from a distribution of a finite number of discrete point charges, all moving at uniform spacing at constant speed in a circle.Previously, most of the Kapteyn series have not been evaluated or, in the case of one of the series, were written in terms of a single elliptic integral, which turned out to be invalid when evaluated numerically (see  (3.18)).  Equation (3.25a) is more appropriate because it represents the series G- in terms of a different, but also elliptic, integral.As has been shown here by recurrence relations, there are only four basic series that need to be calculated, one of which was already known in closed algebraic form. All other of the twelve series can be obtained from direct differentiation or integration of one or other of the four basic series. The series can be evaluated in terms of closed analytic expressions or in terms of integrals that cannot be further reduced. Numerical calculations were carried out to compare the values obtained by direct summation to those obtained from the integral representations, and the relative errors (less than a part in104) were shown to be limited by numerical round-off errors that are responsible for the differences occurring between direct series representations and integral representations of the series.Furthermore, the method presented here may be useful when one has other Kapteyn series of the second kind to consider, thereby providing an additional reason to consider such series anew.


## 4. Kapteyn Series in Quantum-Modulated Systems

Kapteyn series of the second kind also appear in models of even- and odd-order sideband spectra in the optical regime of a quantum system modulated by a high-frequency (e.g., terahertz) electromagnetic field [6] and in certain time-periodic transport problems in superlattices [25, 26]. This section shows that both the even- and the odd-order Kapteyn series that appear can be summed in closed form, thereby allowing more transparent insight into the structural dependence of the sideband spectra and also providing an analytic control for the accuracy of numerical procedures designed to evaluate the series (see also [7]).In discussing an optical analogue for phase-sensitive measurements in quantum transport through a quantum dot whose energy levels are modulated periodically in time, Citrin [6] has considered optical propagation of a monochromatic optical beam at frequency ω (known as the fundamental frequency) transmitted through or reflected from a quantum well modulated by a high-frequency field (henceforth called the terahertz field) at frequency Ω. The transmitted and reflected optical beams are shown to contain new frequencies ω+pΩ where p is an integer, known as terahertz sidebands [6]. The amplitude of such signals as a function of ω is known as terahertz sideband spectra. In the limit that only one modulated energy level (at time-averaged energy ω0) is relevant and the periodic modulation of that energy level is sinusoidal, a simple and useful model can be obtained that permits considerable analytic progress to be made before numerical methods need to be brought to bear on the problem. Such a model then permits one to study in a straightforward fashion how the terahertz sidebands scale with various parameters such as Ω and the modulation strength (the degree to which the energy level varies with respect to its time average ω0).A formally similar analytic model also arises in connection with miniband transport in a superlattice subjected to a strong terahertz field [25, 26]. The phases of the reflected and transmitted complex electromagnetic amplitudes for each sideband (with respect to the initial optical beam at angular frequency ω) provide information on the quantum system. The detailed development given by Citrin [6] has its basic underpinning from the calculation of the amplitude of the transmitted optical electric field, T(ω′,ω), at frequency ω′.  Equation (2.1) of Citrin [6] provides(4.1)T(ω',ω)=2πζ[ω-ɛ0ω-ɛ0′δω',ω+Kp(ω)δω′-ω,pζ],
with(4.2)Kp(ω)=2iΓΔe-ipαS,
where(4.3)S=∑k=1∞1Δ2-(kζ/2)2J(k+p)/2(kɛ12Δ)J(k-p)/2(kɛ12Δ).
The series S is the Kapteyn series of the second kind of interest here. The notation in (4.1) through (4.3) is that given by Citrin [6]. In particular, the prime on the summation indicates that only terms where the parity of k is that of p are retained and Δ=ω-μζ/2-ω0 is the (sideband order μ-dependent) detuning between the average energy ω-μζ/2 of the fundamental and relevant sideband and the time-average energy of the modulated level ω0. The first term in  (4.1) gives the transmitted beam at the input frequency ω′=ω in the absence of the modulation field, while the second contains the terahertz sidebands at ω′=ω+pζ. The cardinal point for this section is the requirement that the sum in  (4.3) is the sum over integers with the same parity as p. Thus if p=2n (with n∈ℕ) then k=2r (with r∈ℕ), while if p=2n+1 then k=2r+1 (with r∈ℕ). Note that due to the form of  (4.3), there is no need to consider negative values of p. Citrin [6] notes that by expanding  (4.3) in powers of ɛ1/ζ1/2 one can identify the various multiphoton processes contributing to each sideband, and he provides the appropriate expansion. Numerical evaluation at this stage is required and has the consequence that convergence of an infinite product inside an infinite sum must be proven, a less than trivial task.The purpose here is to show that the Kapteyn series represented in  (4.3) can indeed be summed in closed form, thereby facilitating not only the general understanding of the sideband spectra but also obviating the need to prove convergence of an infinite product inside an infinite sum—a serendipitous result that is definitely a welcome blessing. Moreover, the closed-form expressions found as well as the approach by which they are obtained are likely to be of interest for other areas of physics and applied mathematics.


### 4.1. Evaluation of the Kapteyn Series

Forp=2n (and so k=2r), that is, for the even-order sideband spectra, one has to evaluate(4.4)SE(n)=∑r=1∞1Δ2-(rζ)2Jr+n(ar)Jr-n(ar),
with a=ɛ1/(2Δ), for all nonnegative integers n.Forp=2n+1 (with n∈ℕ) and so k=2r+1 (with r∈ℕ), that is, for the odd-order side spectra, one has to evaluate(4.5)SO(n)=∑r=0∞1Δ2-(r+1/2)2ζ2Jr+n+1(a(r+12))Jr-n(a(r+12)),
with a=ɛ1/(2Δ), for all integers n including n=0.It is the closed form evaluation of the Kapteyn seriesSE(n) and SO(n) that is of concern here. Thus, (4.4) and (4.5) may be regarded as the starting point of our study.


#### 4.1.1. The Even-Order Side Spectra Summation

Consider first the even-order sideband spectrum summation written in the form(4.6)SE=-(1ζ)2KE(a,b),
with(4.7)KE(a,b)=∑r=1∞1r2-b2Jr+n(ar)Jr-n(ar)
and b=Δ/ζ. Closed-form evaluation of KE(a,b) proceeds as follows. From Watson [15, (1) in Chpter 5.43, page 150] one has(4.8)Jμ(z)Jν(z)=2π∫0π/2dθJμ+ν(2zcosθ)cos(μ-ν)θ,
which is  valid in general when μ and ν are arbitrary integers, and is otherwise valid so long as ℜ𝔢(μ+ν)>-1. One also has Jr-n(ar)=(-1)r-nJn-r(ar). Thus with μ=n+r, ν=n-r, and z=ar it follows that(4.9)Jn+r(z)Jr-n(z)=(-1)r-n∫0π/2dθJ2n(2zrcosθ)cos(2rθ).Now use the representation from  (3.14) in  (4.9) and substitute the result into  (4.7) to obtain(4.10)KE(a,b)=(-1)n(2π)2∫0π/2dψcos(2nψ)∫0π/2dθA(θ),
with(4.11)A(θ)=∑r=1∞(-1)rr2-b2cos(2rθ)cos(2arcosθsinψ).
Use the fact [16, Equation IIII 545 in Chpter 1.445, page 47] that(4.12)∑r=1∞(-1)rr2-b2cos(rf)=12b2-π2bcsc(πb)cos(bf)
valid in the range -π⩽f⩽π. In fact, as is readily obtained from  (4.12), one shows(4.13)∑r=1∞(-1)rr2-b2cos(rf)cos(rg)=12b2-π2bcsc(πb)cos(bf)cos(bg),
which holds for f,g∈[-π,π]. Consequently, we obtain(4.14)KE(a,b)=(-1)n+12πcsc(πb)1b∫0π/2dψcos(2nψ)∫0π/2dθcos(2bθ)cos(2abcosθsinψ).
Care must be exercised that the relevant ranges of the cosine arguments in  (4.14) lie in the appropriate range of modulo (2π) to ensure that one handles the integrals in the correct domain. The bookkeeping associated with values of the cosine arguments outside the range (0,2π) is cumbersome but the general sense of evaluation of the double integral in  (4.14) remains unaltered. For ease of exposition here we treat solely the case where the cosine arguments are restricted to the range (0,2π); all other ranges can be dealt with accordingly, mutatis mutandis.There is also a slight restriction on the argumentb. As Citrin [6] has noted, neglect of any imaginary component of b allows one to obtain an optical theorem [27]. To the same extent, neglect of the imaginary part of b in  (4.14) is equally justified. Then use  (3.14) to write(4.15)KE(a,b)=(-1)n+12πcsc(πb)1b∫0π/2dθcos(2bθ)Jn(2abcosθ).
Again use  (4.8) with μ-ν=2b and μ+ν=2n to obtain(4.16)KE(a,b)=(-1)n+1π2bcsc(πb)Jn+b(ab)Jn-b(ab),
which is the summation required and is valid for n an integer and n⩾1, with 0<a<1 and 0<b<1.Outside of these ranges fora and b one must proceed with the evaluation using the argument given above for validation of the cosine integrals with considerably more bookkeeping as a and b increase systematically. In principle there is no difficulty in completing the evaluations because the method is precisely as given above but the resulting expressions become increasingly unwieldy compared to  (4.16).


#### 4.1.2. The Odd-Order Side Spectra Summation

Consider  (4.5) written in the form(4.17)SO(n)=-(1ζ)2KO(a,b),
with(4.18)KO(a,b)=∑r=0∞(-1)n-r(r+1/2)2-b2Jn+r+1(a(r+12))Jn-r(a(r+12)).By a procedure similar to that followed for the even-order series, one replaces the product of the Bessel functions in  (4.18) by an integral over one Bessel function using  (4.8), then one replaces the single Bessel function occurring under the integral by (see [15, Section 2.2])(4.19)J2n+1(a(r+12)cosθ)=2π∫0π/2dψsin((2n+1)ψ)sin(a(r+12)cosθsinψ),
and then finally one performs the summation over r from r=0 to ∞. Then the reversal of the integral representations is undertaken, just as for the even-order spectra, with the result that one finds(4.20)2KO(a,b)=(-1)nπ2bsec(πb)Jn+1/2+b(ab)Jn+1/2-b(ab),
which is the summation sought, and is valid in 0<b<1/2 and 0<a<1. For values of a and b outside these ranges one has to ensure that the arguments of the various cosine and sine terms in the relevant integrals sit in the appropriate ranges—just as is required for the even-order series.It is noteworthy that the forms of the results for both the even- and odd-order sideband spectra are similar. It is also immediately evident that the given sideband spectrum will vanish ifab is chosen such that it is a zero of the relevant Bessel function.Fortunately, as Citrin [6] has discussed, the parameter b is directly proportional to the detuning frequency and so is considered in some sense as small, in which smallness allowed Citrin [6] to expand the Kapteyn sums in ascending powers of b.The suggestion then is thatb≪1 so that there will be little need to include the higher argument ranges. However, the evaluation of the Kapteyn series for such higher range values for a and b is not complicated, rather fraught with bookkeeping and so is tedious. For this reason only the outline of the procedure has been given here for such ranges. For the ranges most appropriate for the quantum optics and transport experiments discussed by [6], the closed-form detailed evaluations have been given here of the even- and odd-order Kapteyn series.


### 4.2. Numerical Comparison

To illustrate the degree of agreement between the analytical closed form solutions and direct evaluation of the Kapteyn series summations within the ranges chosen, this section of the review provides a few illuminating cases for both the even- and the odd-order summations.


#### 4.2.1. Even-Order Numerical Results

Start with the even-order representations. As shown in Figure8 for the case of n=1, b=0.5, the agreement between direct computation of the value of KE from the series (solid) and the analytic closed-form expression for KE (dashed) is so close that there is no discernable difference between the two curves when plotted as a function of the parameter a in a<1, as is evident also from the inset, which shows the relative deviation between the two curves.Comparison of the direct series evaluation (solid) and the analytic representation (dashed) for the even-order spectraKE(a,b) with n=1, b=0.5 as a is varied. The inset shows the relative difference between the direct series evaluation and the analytic representation.
(a)(b)Consider now the value ofKE at the fixed parameter value a=π/6 as the parameter b varies, again for n=1, the lowest even-order sideband, in Figure 9. The inset clearly shows that there is no discernable difference between the series and the closed form expression for b<1. Indeed, the inset indicates an accuracy of about a part in 1016 throughout most of the range of b<1 and even at b=1 the inaccuracy is still only a part in 1014, thus showing the appropriateness of the closed form analytic expression.Comparison of the direct series evaluation (solid) and the analytic representation (dashed) for the even-order spectraKE with n=1, a=π/6 as b is varied. The inset shows the relative difference between the direct series evaluation and the analytic representation.
(a)(b)Many other values of n have been checked and all curves show marked agreement. For instance, the case ofn=3, shown in Figure 10, indicates no discernable difference between the analytic and series representations in the range a<1, for b=0.95. As seen in the inset, the mismatch is around a part in 1016 throughout most of the range of b and rising only to about a part in 102 at the end of range b=1.Comparison between the direct series evaluation (solid) and the analytic representation (dashed) for the even-order spectraKE with n=3, b=0.95 as a is varied. The inset shows the relative difference between the direct series evaluation and the analytic representation.
(a)(b)


#### 4.2.2. Odd-Order Numerical Results

Similar to the even-order spectra, here we present some illustrative examples of the odd-order spectra results for direct summation of the series in comparison to the analytical results. Figure11 shows the behavior of KO as a function of the parameter a for n=1 and b=0.25. Note that while the analytic result is justified for a<1/2 and b<1/2, the numerical evaluation shows that there is a high degree of overlap beyond the limit for a. Indeed, from Figure 11 one sees that the two results are in agreement out to about a=0.9—perhaps indicative of the larger domain of correctness of the analytic result than is derivable from the arguments given above.Comparison of the direct series evaluation (solid) and the analytic representation (dashed) for the odd-order spectraKO with n=1, b=0.25 as a is varied. The inset shows the relative difference between the direct series evaluation and the analytic representation.
(a)(b)This point can be further extended by considering the case ofn=0 and b=0.95, as shown in Figure 12—well beyond the range where analytic justification can be given without inclusion of the additional terms arising from the arduous bookkeeping. Note that as a function of a there is virtually no difference between the direct series evaluation and the analytic results for the odd-order KO out as far as a=1. This point is further underscored by the inset where agreement to better than a part in 1010 obtains for a<0.5 and, even at a=1, the disagreement is still only about a part in 102. Rugged stability is again seen.Comparison of the direct series evaluation (solid) and the analytic representation (dashed) for the odd-order spectraKO with n=0, b=0.95 as a is varied. The inset shows the relative difference between the direct series evaluation and the analytic representation.
(a)(b)Considering the variation ofKO as a function of the parameter b for fixed values of a (exhibited in Figure 13 is the case of a=0.2 and n=1, although all other cases yield similar results) one sees that there is almost perfect agreement out to b=0.5 and the relative degree of mismatch shown in the inset indicates agreement to better than a part in about 1015 throughout this range of b.Comparison of the direct series evaluation (solid) and the analytic representation (dashed) for the odd-order spectraKO with n=0, b=0.95 as the parameter b is varied. The inset shows the relative difference between the direct series evaluation and the analytic representation.
(a)(b)In short, the analytic evaluations of both the even- and odd-order spectra presented are numerically accurate outside of the ranges where one needs to include extra terms from the phase variation of the various cosine and sine factors appearing under the integral signs. To what extent this general pattern persists for alln, a, and b values is not known at the present time but is likely worth exploring at some future date.


### 4.3. Discussion

In Feise and Citrin [26] the approximate vanishing of the transport at zeros of Bessel functions was noted based on approximate formula. Here we see that this result is exact in the model. If there is loss (b complex) one can also find ab (for a real) that are (complex) zeros of the relevant Bessel functions but then the order of the Bessel functions is also complex so that one would have to work through the problem de novo allowing for a complex value of b from the onset. Such is, regrettably, necessary in order to take into account the complex values of b in summing the series to obtain analytic representations. However, it is also possible to obtain approximate representation of such values because one notes that the original series are even functions of b, whether b is real or complex. As such the series are meromorphic functions of b and so can be analytically continued into the complex plane. Such a detailed investigation is beyond the scope of the present section but the point is noted here for future investigation.It is also opportune to note here that the infinite product representation given by Citrin [6] (obtained by expansion of the infinite Kapteyn series in powers of b) enables one to obtain simply useful expressions for the infinite products by expansion of the analytic results presented for the even- and odd-order series also in powers of b. Such a development is a bit trickier than it appears at first glance because of the presence of the parameter b in both the order and argument of the analytic representations of the series but poses no fundamental difficulties in principle. Such a development would, however, make for a very long paper indeed and so is deferred to the future.Perhaps one of the more interesting points to be made is that the Kapteyn series arising in the sideband spectra can be given in closed form, enabling more insight to be gained into the response of such quantum systems, as illustrated by the vanishing of the sideband spectrum at selected values ofab. The other point to make is that the ability to produce closed form expressions for the Kapteyn series is of considerable benefit when one attempts to perform numerical computations because such closed form expressions act as strong controls on the accuracy determination of any numerical scheme. In addition, the general procedure for summing such Kapteyn series may be of use in other problems where similar Kapteyn series arise.


## 4.1. Evaluation of the Kapteyn Series

Forp=2n (and so k=2r), that is, for the even-order sideband spectra, one has to evaluate(4.4)SE(n)=∑r=1∞1Δ2-(rζ)2Jr+n(ar)Jr-n(ar),
with a=ɛ1/(2Δ), for all nonnegative integers n.Forp=2n+1 (with n∈ℕ) and so k=2r+1 (with r∈ℕ), that is, for the odd-order side spectra, one has to evaluate(4.5)SO(n)=∑r=0∞1Δ2-(r+1/2)2ζ2Jr+n+1(a(r+12))Jr-n(a(r+12)),
with a=ɛ1/(2Δ), for all integers n including n=0.It is the closed form evaluation of the Kapteyn seriesSE(n) and SO(n) that is of concern here. Thus, (4.4) and (4.5) may be regarded as the starting point of our study.


### 4.1.1. The Even-Order Side Spectra Summation

Consider first the even-order sideband spectrum summation written in the form(4.6)SE=-(1ζ)2KE(a,b),
with(4.7)KE(a,b)=∑r=1∞1r2-b2Jr+n(ar)Jr-n(ar)
and b=Δ/ζ. Closed-form evaluation of KE(a,b) proceeds as follows. From Watson [15, (1) in Chpter 5.43, page 150] one has(4.8)Jμ(z)Jν(z)=2π∫0π/2dθJμ+ν(2zcosθ)cos(μ-ν)θ,
which is  valid in general when μ and ν are arbitrary integers, and is otherwise valid so long as ℜ𝔢(μ+ν)>-1. One also has Jr-n(ar)=(-1)r-nJn-r(ar). Thus with μ=n+r, ν=n-r, and z=ar it follows that(4.9)Jn+r(z)Jr-n(z)=(-1)r-n∫0π/2dθJ2n(2zrcosθ)cos(2rθ).Now use the representation from  (3.14) in  (4.9) and substitute the result into  (4.7) to obtain(4.10)KE(a,b)=(-1)n(2π)2∫0π/2dψcos(2nψ)∫0π/2dθA(θ),
with(4.11)A(θ)=∑r=1∞(-1)rr2-b2cos(2rθ)cos(2arcosθsinψ).
Use the fact [16, Equation IIII 545 in Chpter 1.445, page 47] that(4.12)∑r=1∞(-1)rr2-b2cos(rf)=12b2-π2bcsc(πb)cos(bf)
valid in the range -π⩽f⩽π. In fact, as is readily obtained from  (4.12), one shows(4.13)∑r=1∞(-1)rr2-b2cos(rf)cos(rg)=12b2-π2bcsc(πb)cos(bf)cos(bg),
which holds for f,g∈[-π,π]. Consequently, we obtain(4.14)KE(a,b)=(-1)n+12πcsc(πb)1b∫0π/2dψcos(2nψ)∫0π/2dθcos(2bθ)cos(2abcosθsinψ).
Care must be exercised that the relevant ranges of the cosine arguments in  (4.14) lie in the appropriate range of modulo (2π) to ensure that one handles the integrals in the correct domain. The bookkeeping associated with values of the cosine arguments outside the range (0,2π) is cumbersome but the general sense of evaluation of the double integral in  (4.14) remains unaltered. For ease of exposition here we treat solely the case where the cosine arguments are restricted to the range (0,2π); all other ranges can be dealt with accordingly, mutatis mutandis.There is also a slight restriction on the argumentb. As Citrin [6] has noted, neglect of any imaginary component of b allows one to obtain an optical theorem [27]. To the same extent, neglect of the imaginary part of b in  (4.14) is equally justified. Then use  (3.14) to write(4.15)KE(a,b)=(-1)n+12πcsc(πb)1b∫0π/2dθcos(2bθ)Jn(2abcosθ).
Again use  (4.8) with μ-ν=2b and μ+ν=2n to obtain(4.16)KE(a,b)=(-1)n+1π2bcsc(πb)Jn+b(ab)Jn-b(ab),
which is the summation required and is valid for n an integer and n⩾1, with 0<a<1 and 0<b<1.Outside of these ranges fora and b one must proceed with the evaluation using the argument given above for validation of the cosine integrals with considerably more bookkeeping as a and b increase systematically. In principle there is no difficulty in completing the evaluations because the method is precisely as given above but the resulting expressions become increasingly unwieldy compared to  (4.16).


### 4.1.2. The Odd-Order Side Spectra Summation

Consider  (4.5) written in the form(4.17)SO(n)=-(1ζ)2KO(a,b),
with(4.18)KO(a,b)=∑r=0∞(-1)n-r(r+1/2)2-b2Jn+r+1(a(r+12))Jn-r(a(r+12)).By a procedure similar to that followed for the even-order series, one replaces the product of the Bessel functions in  (4.18) by an integral over one Bessel function using  (4.8), then one replaces the single Bessel function occurring under the integral by (see [15, Section 2.2])(4.19)J2n+1(a(r+12)cosθ)=2π∫0π/2dψsin((2n+1)ψ)sin(a(r+12)cosθsinψ),
and then finally one performs the summation over r from r=0 to ∞. Then the reversal of the integral representations is undertaken, just as for the even-order spectra, with the result that one finds(4.20)2KO(a,b)=(-1)nπ2bsec(πb)Jn+1/2+b(ab)Jn+1/2-b(ab),
which is the summation sought, and is valid in 0<b<1/2 and 0<a<1. For values of a and b outside these ranges one has to ensure that the arguments of the various cosine and sine terms in the relevant integrals sit in the appropriate ranges—just as is required for the even-order series.It is noteworthy that the forms of the results for both the even- and odd-order sideband spectra are similar. It is also immediately evident that the given sideband spectrum will vanish ifab is chosen such that it is a zero of the relevant Bessel function.Fortunately, as Citrin [6] has discussed, the parameter b is directly proportional to the detuning frequency and so is considered in some sense as small, in which smallness allowed Citrin [6] to expand the Kapteyn sums in ascending powers of b.The suggestion then is thatb≪1 so that there will be little need to include the higher argument ranges. However, the evaluation of the Kapteyn series for such higher range values for a and b is not complicated, rather fraught with bookkeeping and so is tedious. For this reason only the outline of the procedure has been given here for such ranges. For the ranges most appropriate for the quantum optics and transport experiments discussed by [6], the closed-form detailed evaluations have been given here of the even- and odd-order Kapteyn series.


## 4.1.1. The Even-Order Side Spectra Summation

Consider first the even-order sideband spectrum summation written in the form(4.6)SE=-(1ζ)2KE(a,b),
with(4.7)KE(a,b)=∑r=1∞1r2-b2Jr+n(ar)Jr-n(ar)
and b=Δ/ζ. Closed-form evaluation of KE(a,b) proceeds as follows. From Watson [15, (1) in Chpter 5.43, page 150] one has(4.8)Jμ(z)Jν(z)=2π∫0π/2dθJμ+ν(2zcosθ)cos(μ-ν)θ,
which is  valid in general when μ and ν are arbitrary integers, and is otherwise valid so long as ℜ𝔢(μ+ν)>-1. One also has Jr-n(ar)=(-1)r-nJn-r(ar). Thus with μ=n+r, ν=n-r, and z=ar it follows that(4.9)Jn+r(z)Jr-n(z)=(-1)r-n∫0π/2dθJ2n(2zrcosθ)cos(2rθ).Now use the representation from  (3.14) in  (4.9) and substitute the result into  (4.7) to obtain(4.10)KE(a,b)=(-1)n(2π)2∫0π/2dψcos(2nψ)∫0π/2dθA(θ),
with(4.11)A(θ)=∑r=1∞(-1)rr2-b2cos(2rθ)cos(2arcosθsinψ).
Use the fact [16, Equation IIII 545 in Chpter 1.445, page 47] that(4.12)∑r=1∞(-1)rr2-b2cos(rf)=12b2-π2bcsc(πb)cos(bf)
valid in the range -π⩽f⩽π. In fact, as is readily obtained from  (4.12), one shows(4.13)∑r=1∞(-1)rr2-b2cos(rf)cos(rg)=12b2-π2bcsc(πb)cos(bf)cos(bg),
which holds for f,g∈[-π,π]. Consequently, we obtain(4.14)KE(a,b)=(-1)n+12πcsc(πb)1b∫0π/2dψcos(2nψ)∫0π/2dθcos(2bθ)cos(2abcosθsinψ).
Care must be exercised that the relevant ranges of the cosine arguments in  (4.14) lie in the appropriate range of modulo (2π) to ensure that one handles the integrals in the correct domain. The bookkeeping associated with values of the cosine arguments outside the range (0,2π) is cumbersome but the general sense of evaluation of the double integral in  (4.14) remains unaltered. For ease of exposition here we treat solely the case where the cosine arguments are restricted to the range (0,2π); all other ranges can be dealt with accordingly, mutatis mutandis.There is also a slight restriction on the argumentb. As Citrin [6] has noted, neglect of any imaginary component of b allows one to obtain an optical theorem [27]. To the same extent, neglect of the imaginary part of b in  (4.14) is equally justified. Then use  (3.14) to write(4.15)KE(a,b)=(-1)n+12πcsc(πb)1b∫0π/2dθcos(2bθ)Jn(2abcosθ).
Again use  (4.8) with μ-ν=2b and μ+ν=2n to obtain(4.16)KE(a,b)=(-1)n+1π2bcsc(πb)Jn+b(ab)Jn-b(ab),
which is the summation required and is valid for n an integer and n⩾1, with 0<a<1 and 0<b<1.Outside of these ranges fora and b one must proceed with the evaluation using the argument given above for validation of the cosine integrals with considerably more bookkeeping as a and b increase systematically. In principle there is no difficulty in completing the evaluations because the method is precisely as given above but the resulting expressions become increasingly unwieldy compared to  (4.16).


## 4.1.2. The Odd-Order Side Spectra Summation

Consider  (4.5) written in the form(4.17)SO(n)=-(1ζ)2KO(a,b),
with(4.18)KO(a,b)=∑r=0∞(-1)n-r(r+1/2)2-b2Jn+r+1(a(r+12))Jn-r(a(r+12)).By a procedure similar to that followed for the even-order series, one replaces the product of the Bessel functions in  (4.18) by an integral over one Bessel function using  (4.8), then one replaces the single Bessel function occurring under the integral by (see [15, Section 2.2])(4.19)J2n+1(a(r+12)cosθ)=2π∫0π/2dψsin((2n+1)ψ)sin(a(r+12)cosθsinψ),
and then finally one performs the summation over r from r=0 to ∞. Then the reversal of the integral representations is undertaken, just as for the even-order spectra, with the result that one finds(4.20)2KO(a,b)=(-1)nπ2bsec(πb)Jn+1/2+b(ab)Jn+1/2-b(ab),
which is the summation sought, and is valid in 0<b<1/2 and 0<a<1. For values of a and b outside these ranges one has to ensure that the arguments of the various cosine and sine terms in the relevant integrals sit in the appropriate ranges—just as is required for the even-order series.It is noteworthy that the forms of the results for both the even- and odd-order sideband spectra are similar. It is also immediately evident that the given sideband spectrum will vanish ifab is chosen such that it is a zero of the relevant Bessel function.Fortunately, as Citrin [6] has discussed, the parameter b is directly proportional to the detuning frequency and so is considered in some sense as small, in which smallness allowed Citrin [6] to expand the Kapteyn sums in ascending powers of b.The suggestion then is thatb≪1 so that there will be little need to include the higher argument ranges. However, the evaluation of the Kapteyn series for such higher range values for a and b is not complicated, rather fraught with bookkeeping and so is tedious. For this reason only the outline of the procedure has been given here for such ranges. For the ranges most appropriate for the quantum optics and transport experiments discussed by [6], the closed-form detailed evaluations have been given here of the even- and odd-order Kapteyn series.


## 4.2. Numerical Comparison

To illustrate the degree of agreement between the analytical closed form solutions and direct evaluation of the Kapteyn series summations within the ranges chosen, this section of the review provides a few illuminating cases for both the even- and the odd-order summations.


### 4.2.1. Even-Order Numerical Results

Start with the even-order representations. As shown in Figure8 for the case of n=1, b=0.5, the agreement between direct computation of the value of KE from the series (solid) and the analytic closed-form expression for KE (dashed) is so close that there is no discernable difference between the two curves when plotted as a function of the parameter a in a<1, as is evident also from the inset, which shows the relative deviation between the two curves.Comparison of the direct series evaluation (solid) and the analytic representation (dashed) for the even-order spectraKE(a,b) with n=1, b=0.5 as a is varied. The inset shows the relative difference between the direct series evaluation and the analytic representation.
(a)(b)Consider now the value ofKE at the fixed parameter value a=π/6 as the parameter b varies, again for n=1, the lowest even-order sideband, in Figure 9. The inset clearly shows that there is no discernable difference between the series and the closed form expression for b<1. Indeed, the inset indicates an accuracy of about a part in 1016 throughout most of the range of b<1 and even at b=1 the inaccuracy is still only a part in 1014, thus showing the appropriateness of the closed form analytic expression.Comparison of the direct series evaluation (solid) and the analytic representation (dashed) for the even-order spectraKE with n=1, a=π/6 as b is varied. The inset shows the relative difference between the direct series evaluation and the analytic representation.
(a)(b)Many other values of n have been checked and all curves show marked agreement. For instance, the case ofn=3, shown in Figure 10, indicates no discernable difference between the analytic and series representations in the range a<1, for b=0.95. As seen in the inset, the mismatch is around a part in 1016 throughout most of the range of b and rising only to about a part in 102 at the end of range b=1.Comparison between the direct series evaluation (solid) and the analytic representation (dashed) for the even-order spectraKE with n=3, b=0.95 as a is varied. The inset shows the relative difference between the direct series evaluation and the analytic representation.
(a)(b)


### 4.2.2. Odd-Order Numerical Results

Similar to the even-order spectra, here we present some illustrative examples of the odd-order spectra results for direct summation of the series in comparison to the analytical results. Figure11 shows the behavior of KO as a function of the parameter a for n=1 and b=0.25. Note that while the analytic result is justified for a<1/2 and b<1/2, the numerical evaluation shows that there is a high degree of overlap beyond the limit for a. Indeed, from Figure 11 one sees that the two results are in agreement out to about a=0.9—perhaps indicative of the larger domain of correctness of the analytic result than is derivable from the arguments given above.Comparison of the direct series evaluation (solid) and the analytic representation (dashed) for the odd-order spectraKO with n=1, b=0.25 as a is varied. The inset shows the relative difference between the direct series evaluation and the analytic representation.
(a)(b)This point can be further extended by considering the case ofn=0 and b=0.95, as shown in Figure 12—well beyond the range where analytic justification can be given without inclusion of the additional terms arising from the arduous bookkeeping. Note that as a function of a there is virtually no difference between the direct series evaluation and the analytic results for the odd-order KO out as far as a=1. This point is further underscored by the inset where agreement to better than a part in 1010 obtains for a<0.5 and, even at a=1, the disagreement is still only about a part in 102. Rugged stability is again seen.Comparison of the direct series evaluation (solid) and the analytic representation (dashed) for the odd-order spectraKO with n=0, b=0.95 as a is varied. The inset shows the relative difference between the direct series evaluation and the analytic representation.
(a)(b)Considering the variation ofKO as a function of the parameter b for fixed values of a (exhibited in Figure 13 is the case of a=0.2 and n=1, although all other cases yield similar results) one sees that there is almost perfect agreement out to b=0.5 and the relative degree of mismatch shown in the inset indicates agreement to better than a part in about 1015 throughout this range of b.Comparison of the direct series evaluation (solid) and the analytic representation (dashed) for the odd-order spectraKO with n=0, b=0.95 as the parameter b is varied. The inset shows the relative difference between the direct series evaluation and the analytic representation.
(a)(b)In short, the analytic evaluations of both the even- and odd-order spectra presented are numerically accurate outside of the ranges where one needs to include extra terms from the phase variation of the various cosine and sine factors appearing under the integral signs. To what extent this general pattern persists for alln, a, and b values is not known at the present time but is likely worth exploring at some future date.


## 4.2.1. Even-Order Numerical Results

Start with the even-order representations. As shown in Figure8 for the case of n=1, b=0.5, the agreement between direct computation of the value of KE from the series (solid) and the analytic closed-form expression for KE (dashed) is so close that there is no discernable difference between the two curves when plotted as a function of the parameter a in a<1, as is evident also from the inset, which shows the relative deviation between the two curves.Comparison of the direct series evaluation (solid) and the analytic representation (dashed) for the even-order spectraKE(a,b) with n=1, b=0.5 as a is varied. The inset shows the relative difference between the direct series evaluation and the analytic representation.
(a)(b)Consider now the value ofKE at the fixed parameter value a=π/6 as the parameter b varies, again for n=1, the lowest even-order sideband, in Figure 9. The inset clearly shows that there is no discernable difference between the series and the closed form expression for b<1. Indeed, the inset indicates an accuracy of about a part in 1016 throughout most of the range of b<1 and even at b=1 the inaccuracy is still only a part in 1014, thus showing the appropriateness of the closed form analytic expression.Comparison of the direct series evaluation (solid) and the analytic representation (dashed) for the even-order spectraKE with n=1, a=π/6 as b is varied. The inset shows the relative difference between the direct series evaluation and the analytic representation.
(a)(b)Many other values of n have been checked and all curves show marked agreement. For instance, the case ofn=3, shown in Figure 10, indicates no discernable difference between the analytic and series representations in the range a<1, for b=0.95. As seen in the inset, the mismatch is around a part in 1016 throughout most of the range of b and rising only to about a part in 102 at the end of range b=1.Comparison between the direct series evaluation (solid) and the analytic representation (dashed) for the even-order spectraKE with n=3, b=0.95 as a is varied. The inset shows the relative difference between the direct series evaluation and the analytic representation.
(a)(b)


## 4.2.2. Odd-Order Numerical Results

Similar to the even-order spectra, here we present some illustrative examples of the odd-order spectra results for direct summation of the series in comparison to the analytical results. Figure11 shows the behavior of KO as a function of the parameter a for n=1 and b=0.25. Note that while the analytic result is justified for a<1/2 and b<1/2, the numerical evaluation shows that there is a high degree of overlap beyond the limit for a. Indeed, from Figure 11 one sees that the two results are in agreement out to about a=0.9—perhaps indicative of the larger domain of correctness of the analytic result than is derivable from the arguments given above.Comparison of the direct series evaluation (solid) and the analytic representation (dashed) for the odd-order spectraKO with n=1, b=0.25 as a is varied. The inset shows the relative difference between the direct series evaluation and the analytic representation.
(a)(b)This point can be further extended by considering the case ofn=0 and b=0.95, as shown in Figure 12—well beyond the range where analytic justification can be given without inclusion of the additional terms arising from the arduous bookkeeping. Note that as a function of a there is virtually no difference between the direct series evaluation and the analytic results for the odd-order KO out as far as a=1. This point is further underscored by the inset where agreement to better than a part in 1010 obtains for a<0.5 and, even at a=1, the disagreement is still only about a part in 102. Rugged stability is again seen.Comparison of the direct series evaluation (solid) and the analytic representation (dashed) for the odd-order spectraKO with n=0, b=0.95 as a is varied. The inset shows the relative difference between the direct series evaluation and the analytic representation.
(a)(b)Considering the variation ofKO as a function of the parameter b for fixed values of a (exhibited in Figure 13 is the case of a=0.2 and n=1, although all other cases yield similar results) one sees that there is almost perfect agreement out to b=0.5 and the relative degree of mismatch shown in the inset indicates agreement to better than a part in about 1015 throughout this range of b.Comparison of the direct series evaluation (solid) and the analytic representation (dashed) for the odd-order spectraKO with n=0, b=0.95 as the parameter b is varied. The inset shows the relative difference between the direct series evaluation and the analytic representation.
(a)(b)In short, the analytic evaluations of both the even- and odd-order spectra presented are numerically accurate outside of the ranges where one needs to include extra terms from the phase variation of the various cosine and sine factors appearing under the integral signs. To what extent this general pattern persists for alln, a, and b values is not known at the present time but is likely worth exploring at some future date.


## 4.3. Discussion

In Feise and Citrin [26] the approximate vanishing of the transport at zeros of Bessel functions was noted based on approximate formula. Here we see that this result is exact in the model. If there is loss (b complex) one can also find ab (for a real) that are (complex) zeros of the relevant Bessel functions but then the order of the Bessel functions is also complex so that one would have to work through the problem de novo allowing for a complex value of b from the onset. Such is, regrettably, necessary in order to take into account the complex values of b in summing the series to obtain analytic representations. However, it is also possible to obtain approximate representation of such values because one notes that the original series are even functions of b, whether b is real or complex. As such the series are meromorphic functions of b and so can be analytically continued into the complex plane. Such a detailed investigation is beyond the scope of the present section but the point is noted here for future investigation.It is also opportune to note here that the infinite product representation given by Citrin [6] (obtained by expansion of the infinite Kapteyn series in powers of b) enables one to obtain simply useful expressions for the infinite products by expansion of the analytic results presented for the even- and odd-order series also in powers of b. Such a development is a bit trickier than it appears at first glance because of the presence of the parameter b in both the order and argument of the analytic representations of the series but poses no fundamental difficulties in principle. Such a development would, however, make for a very long paper indeed and so is deferred to the future.Perhaps one of the more interesting points to be made is that the Kapteyn series arising in the sideband spectra can be given in closed form, enabling more insight to be gained into the response of such quantum systems, as illustrated by the vanishing of the sideband spectrum at selected values ofab. The other point to make is that the ability to produce closed form expressions for the Kapteyn series is of considerable benefit when one attempts to perform numerical computations because such closed form expressions act as strong controls on the accuracy determination of any numerical scheme. In addition, the general procedure for summing such Kapteyn series may be of use in other problems where similar Kapteyn series arise.


## 5. Bessel Function Summations in Plasma Instability Problems

The history of small amplitude plasma wave instabilities, treated from a plasma kinetic point of view, has had a long history since the Second World War. In the 1950's there was a major effort made to understand the unstable behavior of one component and multicomponent kinetic plasmas with a large number of instabilities being uncovered, ranging from the Bunemann [28] two-stream instability through to the transverse Weibel [29] instability. Now most of these early investigations were undertaken either as one-dimensional plasma problems or as three dimensional problems involving wave dispersion relations with coupling of the transverse and longitudinal components of the perturbing fields. Considerable progress in cataloguing such instabilities was made in one, two, and three dimensions. In most cases the early investigations were undertaken assuming there was no background magnetic field because the complexity of such an ambient magnetic field was formidable.The main problem that arose in plasma investigations involving an ambient magnetic field was the fact that each entry in the3×3 determinant describing the wave dispersion characteristics involved an infinite sum of Bessel functions. (This problem does not occur in the absence of the ambient magnetic field when the entries in the 3×3 dispersion relation are surprisingly simple by comparison). Thus, in evaluating the general 3×3 determinant one was left with triply infinite sums of Bessel functions that also involved integrals over the particles distribution functions. Except for special directions of wave propagation (e.g., parallel or exactly perpendicular to the ambient magnetic field) this complexity effectively stymied a general investigation of the dispersion characteristics for arbitrary plasma distribution functions. In addition, the speed of computers was severely limited in comparison to modern computers so that term by term evaluation of the triply infinite sums and performance of the distribution function integrals was seriously impaired. What was needed was a procedure to sum, in closed form if possible or as an integral representation if not analytically possible, the Bessel function series. If such were to be attainable then one could make significant progress with the general 3×3 dispersion relation for arbitrary particle distribution functions.This section of the review is concerned with expressing the plasma dispersion relation series in closed form. Because the technical details are somewhat involved, the purpose here is restricted to bringing to the reader's attention the facts that such sums can be done with references to the original literature where such series and the methods can be found, together with the influence such summation procedures have on simplifying both3×3 dispersion relations for plasmas.


### 5.1. Plasma Dispersion Relation Series

In the3×3 dispersion relation for plasma waves the basic series that occurs is(5.1)P1(a,b)=∑n=-∞∞Jn2(b)n-a,
where a and b are related to the plasma wave frequency and to the wavenumber, with b=0 for waves propagating parallel to the magnetic field. Other sums that occur in the 3×3 determinant are related to P1 either by integration or direct differentiation so that it is both sufficient and necessary to perform P1 in closed form and then the other series fall into line. This point has been emphasized quite strongly in Schlickeiser [30].From a purely plasma physics point of view the seriesP1 was evaluated in closed form, apparently for the first time, by Lerche [31]. The closed form result he gives for the sum is(5.2)P1=-πcsc(aπ)Ja(b)J-a(b).
(The procedure is to replace the square of the Bessel function in  (5.1) using  (4.9) and to then perform the summation before reversing the argument—just as for the KE and KO Kapteyn summations in the previous section of the review).In the late 1960's and early 1970's the astrophysical community was starting to discuss the influence of an ambient interstellar field on plasma waves for both nonrelativistic and relativistic plasmas and ran directly up against the seriesP1 plus a couple of simple (sic!) other Bessel function series. In order to promulgate the simplifications that ensued given the summed series, it was appropriate to reiterate for the astrophysical community the basic summation obtained in 1966 and also to present the residual Bessel function sums. Such results can be found in Lerche [32].Some 16 years after the original summation was given by Lerche [31], Newberger [33] rediscovered the summation of P1 and confirmed the correctness of the sum. In a recent article, Qin et al. [34] have used the summation of P1 and, surprisingly, attribute its first derivation to Newberger, seemingly being completely unaware of the Bessel function summation procedure arising from the earlier investigations by Lerche in the fields of both plasma physics and astrophysics. Swanson [35] seems similarly unaware of the work of Lerche [31, 32, 36] on the Bessel functions summation procedure for he, too, quotes only the later work by Newberger [33]. This section of the review has then set the historical record straight in regard to the closed form expression for the Bessel function sum P1.What is also of interest is then the use of this summation procedure in astrophysical plasmas and of laboratory plasma physics. As was shown in Lerche [31, 36] and ably summarized in Schlickeiser [30] in complete form for relativistic and nonrelativistic plasmas, one can then write a much simpler form for the full 3×3 dispersion relation for plasma waves in an ambient magnetic field than would otherwise be possible. In particular, Equations (55a), (55b), and (55c) of Lerche [36] provide the closed form expressions for the determinant components; this result precedes the development given in Swanson [35] by approximately 15 years despite all protestations to the contrary. For instance, the work of San Martin [37] on plasma waves in oscillating electromagnetic fields was shown by Lerche [36] to be addressed more easily using the Bessel function sum rule and its generalizations given in Lerche [36]. In particular, (40a) and (42a) of Lerche [36] give one of the generalizations as(5.3)Q(a,b)=∑n=-∞∞(-1)nn-aJn+μ(b)J-(n+ν)(b)=-πcsc(aπ)J-(ν+a)(b)Jμ+a(b)
and many other generalizations are also provided in Lerche [36].This general result would seem to precede the same result given later by [33]. This original set of investigations has been confirmed in a recent article by Qin et al. [34] who, again surprisingly, make no broad-based mention of the original investigations by Lerche [31, 32, 36] and Schlickeiser [30], and would seem to be more concerned with noting the simplification. The nominal reason would appear to be that with modern computers one can cut the numerical computation time by orders of magnitudes with such summation methods—a point that was also made considerably earlier [36] after the original discovery of the summation technique by Schlickeiser [31].


### 5.2. Discussion

This section has spelled out the reduction of infinite sums of Bessel functions to closed form expressions for plasma problems. In addition, the use of such summations in simplifying discussion of plasma dispersion relations has been noted and emphasized following on from the original work on the summations of over forty years ago, most of which work seems to be underappreciated at the present day. It would be appropriate if modern and future articles were to be aware of the historical development so that further reinvention of existing information is not necessary.


## 5.1. Plasma Dispersion Relation Series

In the3×3 dispersion relation for plasma waves the basic series that occurs is(5.1)P1(a,b)=∑n=-∞∞Jn2(b)n-a,
where a and b are related to the plasma wave frequency and to the wavenumber, with b=0 for waves propagating parallel to the magnetic field. Other sums that occur in the 3×3 determinant are related to P1 either by integration or direct differentiation so that it is both sufficient and necessary to perform P1 in closed form and then the other series fall into line. This point has been emphasized quite strongly in Schlickeiser [30].From a purely plasma physics point of view the seriesP1 was evaluated in closed form, apparently for the first time, by Lerche [31]. The closed form result he gives for the sum is(5.2)P1=-πcsc(aπ)Ja(b)J-a(b).
(The procedure is to replace the square of the Bessel function in  (5.1) using  (4.9) and to then perform the summation before reversing the argument—just as for the KE and KO Kapteyn summations in the previous section of the review).In the late 1960's and early 1970's the astrophysical community was starting to discuss the influence of an ambient interstellar field on plasma waves for both nonrelativistic and relativistic plasmas and ran directly up against the seriesP1 plus a couple of simple (sic!) other Bessel function series. In order to promulgate the simplifications that ensued given the summed series, it was appropriate to reiterate for the astrophysical community the basic summation obtained in 1966 and also to present the residual Bessel function sums. Such results can be found in Lerche [32].Some 16 years after the original summation was given by Lerche [31], Newberger [33] rediscovered the summation of P1 and confirmed the correctness of the sum. In a recent article, Qin et al. [34] have used the summation of P1 and, surprisingly, attribute its first derivation to Newberger, seemingly being completely unaware of the Bessel function summation procedure arising from the earlier investigations by Lerche in the fields of both plasma physics and astrophysics. Swanson [35] seems similarly unaware of the work of Lerche [31, 32, 36] on the Bessel functions summation procedure for he, too, quotes only the later work by Newberger [33]. This section of the review has then set the historical record straight in regard to the closed form expression for the Bessel function sum P1.What is also of interest is then the use of this summation procedure in astrophysical plasmas and of laboratory plasma physics. As was shown in Lerche [31, 36] and ably summarized in Schlickeiser [30] in complete form for relativistic and nonrelativistic plasmas, one can then write a much simpler form for the full 3×3 dispersion relation for plasma waves in an ambient magnetic field than would otherwise be possible. In particular, Equations (55a), (55b), and (55c) of Lerche [36] provide the closed form expressions for the determinant components; this result precedes the development given in Swanson [35] by approximately 15 years despite all protestations to the contrary. For instance, the work of San Martin [37] on plasma waves in oscillating electromagnetic fields was shown by Lerche [36] to be addressed more easily using the Bessel function sum rule and its generalizations given in Lerche [36]. In particular, (40a) and (42a) of Lerche [36] give one of the generalizations as(5.3)Q(a,b)=∑n=-∞∞(-1)nn-aJn+μ(b)J-(n+ν)(b)=-πcsc(aπ)J-(ν+a)(b)Jμ+a(b)
and many other generalizations are also provided in Lerche [36].This general result would seem to precede the same result given later by [33]. This original set of investigations has been confirmed in a recent article by Qin et al. [34] who, again surprisingly, make no broad-based mention of the original investigations by Lerche [31, 32, 36] and Schlickeiser [30], and would seem to be more concerned with noting the simplification. The nominal reason would appear to be that with modern computers one can cut the numerical computation time by orders of magnitudes with such summation methods—a point that was also made considerably earlier [36] after the original discovery of the summation technique by Schlickeiser [31].


## 5.2. Discussion

This section has spelled out the reduction of infinite sums of Bessel functions to closed form expressions for plasma problems. In addition, the use of such summations in simplifying discussion of plasma dispersion relations has been noted and emphasized following on from the original work on the summations of over forty years ago, most of which work seems to be underappreciated at the present day. It would be appropriate if modern and future articles were to be aware of the historical development so that further reinvention of existing information is not necessary.


## 6. Kapteyn Series in Other Fields

There are many other fields in which Kapteyn series have been used. The first use of Kapteyn series had to do with Kepler's problem of a particle moving in an ellipse under the action of a center of force at the focus, attracting the particle according to an inverse square law (e.g., [38]). The problem is completely determined by the coordinates r, ϑ, and φ as functions of the time t or according to Kepler's Second Law as functions of the mean anomaly ψ. For example, the equation connecting φ and ψ can be written as(6.1)ψ=φ-ɛsinφ,
where, in the specific context, ɛ∈[0,1] is understood as the ellipse eccentricity. Expanding φ-ψ in a Fourier sine series as(6.2)φ-ψ=∑n=1∞Cn(ɛ)sin(nψ).
An easy calculation allows one to specialize the formal expressions for the coefficients Cn(ɛ) as(6.3)Cn(ɛ)=2π∫0πdψ(φ-ψ)sin(nψ)=2πn∫0πdϕcos(nφ-nɛsinφ)=2nJn(nɛ).
Thus, the solution to  (6.1) is written as [38](6.4)φ(ψ,ɛ)=ψ+2∑n=1∞sin(nψ)nJn(nɛ),
which was the first appearence of a Kapteyn series [1].Another example is the field neighboring to plasma instability theory, namely the theory of charged particles being scattered in turbulent electromagnetic fields. In Shalchi and Schlickeiser [10] it was shown that(6.5)∑n=1∞n-2Jn'2(nz)=14-z28
which was derived from the well-known Kapteyn series Watson [15](6.6)∑n=1∞n-2Jn2(nz)=z24.Furthermore, Dominici [8] has proven that(6.7)F(z,t)=∑n=1∞tnJn(nz),
where z∈C and t∈R, can be written in terms of a power series in z in the form(6.8)F(z,t)=∑n=1∞An(t)zn,
with the coefficients An(t) given by(6.9)An(t)=∑k=0ncos(π(n-k)/2)kn(n-k)!!(n+k)!!=(-1)nn!∑k=0⌊n/2⌋(-1)k(nk)(k-n2)ntn-2k.Other work focuses on specific examples that use Kapteyn-type series of generalized Bessel functions [38–41], or problems from queueing theory [9], where, in the latter case, a transcendental equation of the type(6.10)1=D2(D+1)∑n=-∞∞Jn(nD+1)Cn
was solved.These few additional examples underline the importance of Kapteyn series not only in plasma physics but in a variety of different problems.


## 7. Summary and Conclusion

The dominant purpose of this review has been to show, using specific examples, how one can bring some form of order to physics problems involving Kapteyn series. Rather than being content to discuss just approximations to such series, what has been shown here is that there are available techniques for summing such series in closed form or of reducing the series to simple integrals that cannot be further reduced. The method and procedures for handling such Kapteyn series are by no means exhausted with the presentation given here; indeed one would seriously doubt any such claim because the number and types of such Kapteyn series are legion.However, what has been shown is that the basic development of an extremely large number of Kapteyn series in mathematical physics problems can be resolved into a small subset of basic arguments. Basically one either reduces the desired Kapteyn series to differential or integral operators acting on a known Kapteyn series or one uses basic properties of Bessel functions to write the desired Kapteyn series in terms of an integral transform and reverses the order of integration and summation. Combinations of both procedures have proven powerful tools in allowing one to effect summations.A further advantage of such summation procedures is that one can determine to what extent approximations to the series are accurate and over what range of parameters accuracy is maintained.Perhaps one other advantage is that the basic procedures can be used for other forms of Bessel function summation, as was also demonstrated, and most likely for Kapteyn series occurring in physics problems other than the selected variety exhibited here.The ability to understand the underpinning physics is enhanced once one obtains closed for solutions to Kapteyn series, a point made in the review directly in connection with examples from the realms of both pulsars and quantum systems.Modern methods of Kapteyn series summations owe their development to the underpinning basic ideas expressed so many years ago at the turn of and in the first quarter of the 20th Century. While modern computers are sufficiently fast often to enable one to numerically sum Kapteyn series without needing to have available analytical representations of such series, there is not only a loss of elegance in such numerical procedures but there is also only a limited understanding brought to the fundamental physics in comparison to what can be gleaned from closed form expression for the desired Kapteyn series. This point is, perhaps, one of the major reasons one undertakes direct summation procedures.We would be highly interested to hear of other procedures for effecting Kapteyn series summations that either are based on the fundamental procedures detailed here or that are completely novel. Until such time it would appear that one should not be too hasty to approximate Kapteyn sums analytically or numerically for there is much that can be learned from the methods that is of great importance to mathematical physics problems that have as kernels some form of Kapteyn series—and many such problems are strewn throughout the scientific literature as the examples given here have been designed to illustrate.



---

*Source: 425164-2010-03-02.xml*