A Review of Procedures for Summing Kapteyn Series in Mathematical Physics

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.


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 that is, Kapteyn series of the first kind involve summations over terms containing one Bessel function of the form J n 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.

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 and the force, F, acting on the dipole in the z direction is 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/c sin θ with s being the offset distance of the dipole from the spin axis, and where Note that a ∈ −1, 1 is required so that the series S 1 and S 2 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 O 1 so that na 1 is almost nowhere valid.To investigate such situations one needs closed form expressions for the two series S 1 and S 2 .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.Consider then S 1 .Use the formula 15, Section 5.43

Manipulations with the
in the form e.g., 16, Section 6.681 Expression 2.11 shows that S 1 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 where a m is arbitrary but given, is known in closed form, then the series is given by two simple integrations because 14 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 Reversing the argument: if F z is known, then f z is given directly by differentiation of F z in 2.14 .

Reduction of S 1 to Closed Form
From 2.11 -2.14 we have that with a m 0 if m is odd, and a m m 4 if m is even where b a cos ψ.But it is well known that 15, Section 17.33 ∓1 m J m mz , 2.17 where z ∈ −1, 1 , so that, also for b ∈ −1, 1 ,

2.19
Carrying out the differentiations in 2.19 yields with b a cos ψ.

Advances in Mathematical Physics
Using a partial fraction expansion and where ξ a 2 / 2 − a 2 , and using Thus, this procedure shows that both S 1 and S 2 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 10 16 .The prototype of such Kapteyn series of the second kind was first given in closed form by Schott 17 , who evaluated

2.25
Note that, in 2.25 , there appears a factor 1−a 2 −7/2 which was missing in Lerche and Tautz 3 .
The basic procedure for evaluating Kapteyn series of the generic form where m is either an integer or half integer positive or negative and a n is either unity or −1 n , then follows the same recipe as given here, although the expressions rapidly become unwieldy as m becomes large.
Advances in Mathematical Physics 7

Calculation of L and F
The closed form expressions for S 1 and S 2 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

2.27
Note that the zero-order terms L 0 and F 0 agree with the expansion given by Harrison and Tademaru 2 in their Equations 5 and 7 .
Next, the expression for L is separated as and the functions L ρ , L φ , and L z 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 .
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.008 s/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.

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 has N charges equally spaced around a circle and all Advances in Mathematical Physics In panels 1, 3, and 5, the solid lines show the numerically calculated exact functions L ρ , L φ , and L z , respectively, and the dashed and dash-dot lines show the approximations F 0 and F 0 F 1 , respectively.All function values are normalized to Ω 4 /c 3 .In panels 2, 4, and 6, the relative deviation in percent from the exact function values is shown for the approximations L 0 solid lines and L 0 L 1 dashed lines , respectively.The two dotted lines mark deviations of 10% and 50%.In the upper panel, the solid line shows the numerically calculated exact function F, and the dashed and dash-dot lines show the approximations F 0 and F 0 F 1 , respectively.In the lower panel, the relative deviation in percent from the exact function values is shown for the approximations F 0 solid line and F 0 F 1 dashed line , respectively.The two dotted lines mark deviations of 10% and 50%.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 Advances in Mathematical Physics 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.

The Sets of Series
The twelve series in question are given by 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 with m 0 in the cases of S 1 , S 3 , and S 6 and the set of series with m −1 in the cases of S 2 , S 4 , and S 5 .
The reason for these reductions is as follows.One can write so that it is sufficient to obtain S 1 , S 2 , S 3 , and S 4 .

Advances in Mathematical Physics 11
Note also that But, because of Bessel's equation see 2.6 with a changed to b , one has Thus it is sufficient to obtain S 1 and S 2 .
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.

3.7
But the series is precisely of the form required in Watson's theorem, with a n 0 if n is odd and a n exp inπ/2 ln λ n 2m if n is even, so that where the differential operator L from 2.15 has been used.Hence, for m > 0 all series of the type S 1 can be reduced to the determination of h 0 b by differentiation.Equally, for m < 0 one can use Watson's theorem in the converse sense to note that so that, by two integrations, one has a recursive relation leading directly to h 0 .Thus, all twelve of the basic series needed can be written in terms of four fundamental series 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.

The Two Series Represented by
Now, in F , replace the Bessel functions using again 2.10 while in 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.16b
Numerical investigation by direct summation of F 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 10 4 ; this limit on resolution being caused by numerical round-off error.Figures 4 and 5  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.
Throughout this review, the numerical evaluation of infinite sums is carried out as follows: First, a number of terms usually 1000 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 F M ≡ 1/2 ∂F /∂b, written in the form J n nb J n nb , 3.17 could be represented by a single elliptic integral his 2.22 as

3.18
Figure 6 shows plots as a function of b of both the sum F M and the elliptic integral representation, G M 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.

3.19b
The series G has been known in closed form since the time of Schott 17 .Use the well-known fact 1 that J n nb cos n φ − b sin φ .

3.20
Integrate 3.20 over 0 φ π, thereby obtaining which is just Schott's 17 formula.The series G − is considerably more complicated to evaluate.Write J 0 2nb cos ψ cos 2nψ.

3.22
Now use the Schl ömilch 24 formula, which states that any function which is given through an arbitrary but known function Γ, can be rewritten as

3.24
With the identifications w 2ψ and x 2b cos ψ, 3.22 then yields where the upper integration limit is implicitly given by ψ b cos ψ , or b is given explicitly by b ψ secψ .One can then write 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 10 5 , 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.
Advances in Mathematical Physics

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 in 10 4 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.

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 .

19
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 where 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 ∈ N then k 2r with r ∈ N , while if p 2n 1 then k 2r 1 with r ∈ N .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.

Evaluation of the Kapteyn Series
For p 2n and so k 2r , that is, for the even-order sideband spectra, one has to evaluate with a ε 1 / 2Δ , for all nonnegative integers n.
For p 2n 1 with n ∈ N and so k 2r 1 with r ∈ N , that is, for the odd-order side spectra, one has to evaluate with a ε 1 / 2Δ , for all integers n including n 0. It is the closed form evaluation of the Kapteyn series S E n and S O n that is of concern here.Thus, 4.4 and 4.5 may be regarded as the starting point of our study.

The Even-Order Side Spectra Summation
Consider first the even-order sideband spectrum summation written in the form with

4.14
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 argument b.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 Again use 4.8 with μ − ν 2b and μ ν 2n to obtain Advances in Mathematical Physics 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 for a 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 .

The Odd-Order Side Spectra Summation
Consider 4.5 written in the form By a procedure similar to that followed for the even-order series, one replaces the product of the Bessel functions in 4. 18  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 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 if ab 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 that b 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.

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.

Even-Order Numerical Results
Start with the even-order representations.As shown in Figure 8 for the case of n 1, b 0.5, the agreement between direct computation of the value of K E from the series solid and the analytic closed-form expression for K E 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.Consider now the value of K E 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 10 16 throughout most of the range of b < 1 and even at b 1 the inaccuracy is still only a part in 10 14 , thus showing the appropriateness of the closed form analytic expression.Many other values of n have been checked and all curves show marked agreement.For instance, the case of n 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 10 16 throughout most of the range of b and rising only to about a part in 10 2 at the end of range b 1.

Odd-Order Numerical Results
Similar to the even-order spectra, here we present some illustrative examples of the oddorder spectra results for direct summation of the series in comparison to the analytical results.Figure 11 shows the behavior of K O 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.
This point can be further extended by considering the case of n 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 K O out as far as a 1.This point is further underscored by the inset where agreement to better than a part in 10 10 obtains for a < 0.5 and, even at a 1, the disagreement is still only about a part in 10 2 .Rugged stability is again seen.Considering the variation of K O 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 10 15 throughout this range of 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 all n, a, and b values is not known at the present time but is likely worth exploring at some future date.

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 of ab.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.

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, Advances in Mathematical Physics 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 onedimensional 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 the 3 × 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 both 3 × 3 dispersion relations for plasmas.

Plasma Dispersion Relation Series
In the 3 × 3 dispersion relation for plasma waves the basic series that occurs is 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 P 1 either by integration or direct differentiation so that it is both sufficient and necessary to perform P 1 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 series P 1 was evaluated in closed form, apparently for the first time, by Lerche 31 .The closed form result he gives for the sum is Advances in Mathematical Physics

29
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 K E and K O 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 series P 1 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 P 1 and confirmed the correctness of the sum.In a recent article, Qin et al. 34 have used the summation of P 1 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 P 1 .
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 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 .

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.

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 where, in the specific context, ε ∈ 0, 1 is understood as the ellipse eccentricity.Expanding ϕ − ψ in a Fourier sine series as C n ε sin nψ .

6.2
An easy calculation allows one to specialize the formal expressions for the coefficients C n ε as A n t z n , 6.8 with the coefficients A n t given by 6.9 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 C n 6.10 was solved.These few additional examples underline the importance of Kapteyn series not only in plasma physics but in a variety of different problems.

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.

Advances in Mathematical Physics
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.

Figure 1 :
Figure 1: Sketch of the spin and dipole coordinates from Harrison and Tademaru 2 .

Figure 2 :
Figure 2: Comparison of the exact and approximate values for the three components of the function L.In panels 1, 3, and 5, the solid lines show the numerically calculated exact functions L ρ , L φ , and L z , respectively, and the dashed and dash-dot lines show the approximations F 0 and F 0 F 1 , respectively.All function values are normalized to Ω 4 /c 3 .In panels 2, 4, and 6, the relative deviation in percent from the exact function values is shown for the approximations L 0 solid lines and L 0 L 1 dashed lines , respectively.The two dotted lines mark deviations of 10% and 50%.

Figure 3 :
Figure 3: Comparison of the exact and approximate values for the function F, normalized to Ω 4 c −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 F 0 and F 0 F 1 , respectively.In the lower panel, the relative deviation in percent from the exact function values is shown for the approximations F 0 solid line and F 0 F 1 dashed line , respectively.The two dotted lines mark deviations of 10% and 50%.

Figure 4 :
Figure 4:The series F from 3.12a with the relative error when compared to the integral from 3.16a .

Figure 5 :
Figure 5: The series F − from 3.12b with the relative error when compared to the integral from 3.16b .

Figure 6 :
Figure 6:The series F M from 3.17 solid line compared to the integral representation G M 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.

,Figure 7 :
Figure 7: 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 .

Figure 8 :
Figure 8: Comparison of the direct series evaluation solid and the analytic representation dashed for the even-order spectra K E 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.

Figure 9 :
Figure 9: Comparison of the direct series evaluation solid and the analytic representation dashed for the even-order spectra K E with n 1, a π/6 as b is varied.The inset shows the relative difference between the direct series evaluation and the analytic representation.

Figure 10 :
Figure 10: Comparison between the direct series evaluation solid and the analytic representation dashed for the even-order spectra K E 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.

Figure 11 :Figure 12 :Figure 13 :
Figure 11: Comparison of the direct series evaluation solid and the analytic representation dashed for the odd-order spectra K O 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.
Series S 1 and S 2 because the integration constant in 2.8 is zero by evaluation as a → 0.Thus, it is sufficient to evaluate S 1 in closed form and to perform the integration in 2.8 to obtain S 2 .
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 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 C and t ∈ R, can be written in terms of a power series in z in the form ∞ n 1