While it is common knowledge that the two-body Kepler problem can be solved exactly as an orbit equation, it is perhaps less commonly known that there is also a complete solution as a function of time in terms of Kapteyn series. These are series of the form

$∑k=0 ∞ ckJk (kz)$

where the index of the Bessel function appears as a factor in its argument. Such series with appropriate coefficients will converge as long as the quantity

$|z exp(1 -z2) 1+1 -z2 |$

remains finite, which is true even for extremely eccentric orbits in the Kepler problem. This presentation will evaluate the explicit temporal expansions for dynamic variables of interest and their combinations. It expands on the cursory treatment given in §17.21 of Watson’s venerable Bessel Functions.

The equation of an orbit with semimajor axis a and eccentricity e, along with its parametrization in terms of the eccentric anomaly u, is

$r=a( 1-e2) 1+ecosφ =a(1 -ecosu)$

The relation between the two angles in the invariant plane can be rearranged to give

$1-e2 1+ecosφ =1-ecosu ⇒ cosφ =cosu-e 1-ecosu sinφ =1 -e2sinu 1-ecosu$

from which one immediately has

$x=rcosφ =a(cosu -e) y=rsinφ =a1 -e2sinu$

Note that the x-coordinate is related linearly to the radial variable. Integrating the energy equation in terms of the eccentric anomaly produces Kepler’s equation,

$ωt =u-esinu ω=k ma3$

where the frequency of the orbit is dependent on the reduced mass m and the gravitational interaction constant k as well as the semimajor axis. Differentiating both sides of Kepler’s equation

$d(ωt) ≡dτ =(1 -ecosu) du$

provides a differential relationship between the variables in the equation. The symbol τ as defined here will be used as shorthand to simplify integral expressions.

The Kapteyn series representing dynamic variables arise from Fourier sine and cosine expansions. Given a periodic function of the eccentric anomaly, a sine series is determined by

$f(u) =∑k=1 ∞ sksinkωt sk=2π ∫0π f(u) sinkτ dτ$

while a cosine series is determined by

$f(u) =∑k=0 ∞ ckcoskωt c0=1π ∫0π f(u) dτ ck=2π ∫0π f(u) coskτ dτ$

The choice between the two expansions depends on the symmetry of the periodic function: an odd function will be expanded as a sine series and an even function as a cosine series.

Evaluation of the coefficients necessitates an integral representation of Bessel functions. While there are various ways to write such a representation, the one that is most useful in this context is

$Jn(z) =1π ∫0π cos(nθ -zsinθ) dθ$

which holds for n integral. As an aside, the corresponding integral with a sine instead of a cosine is related to the Bessel function ${Y}_{n}\left(z\right)$ .

Differentiation with respect to the argument z under the integral provides an integral representation of the Bessel function derivative:

$Jn′ (z) =1π ∫0π sinθ sin(nθ -zsinθ) dθ$

Application of product identities for circular functions gives

$1π ∫0π cosmθ cos(nθ -zsinθ) dθ =12 [Jn-m (z) +Jn+m (z)] 1π ∫0π sinmθ sin(nθ -zsinθ) dθ =12 [Jn-m (z) -Jn+m (z)]$

Two recursion relations for Bessel functions

$Jn-1 (z) +Jn+1 (z) =2nz Jn(z) Jn-1 (z) -Jn+1 (z) =2Jn′ (z)$

will be useful in what follows. The second can be seen immediately in the integral representations above, while the first implies

$1π ∫0π cosθ cos(nθ -zsinθ) dθ =nz Jn(z)$

Keep in mind that a prime will always indicate differentiation with respect to the entire argument of the Bessel function. A second derivative can always be replaced as

$Jn″ (z) =−1z Jn′ (z) +[n2 z2-1] Jn(z)$

using the defining differential equation for Bessel functions.

Inversion of Kepler’s equation requires an expansion for the sine of the eccentric anomaly. The twist to the integral involved and those that follow is that the integration is with respect to the mean anomaly τ, not the eccentric anomaly. This will add a factor to the integrand of the differential relationship between the variables in the equation.

Integrals can always be evaluated by expanding products of circular functions and applying recursion relations to the resulting Bessel functions, but often an intermediate integration by parts is more efficient. This is particularly true with a single differential factor in the integrand.

The coefficients for a sine expansion of the sine of the eccentric anomaly are

$2π ∫0π sinusinkτ dτ =2π ∫0π sinusin[k(u -esinu)] (1-ecosu) du =2πk ∫0π cosucos[k(u -esinu)] du =2ke Jk(ke)$

where the endpoint contributions to the integral are zero. The expansion for the sine of the eccentric anomaly is thus

$sinu=2e ∑k=1 ∞ Jk (ke)k sinkωt$

and that for the eccentric anomaly itself is

$u=ωt +2∑k=1 ∞ Jk (ke)k sinkωt$

The corresponding cosine expansion for the cosine of the eccentric anomaly has a constant term

$1π ∫0π cosudτ =1π ∫0π cosu(1 -ecosu)du =−e2$

and the remaining coefficients are

$2π ∫0π cosucoskτ dτ =2π ∫0π cosucos[k(u -esinu)] (1-ecosu) du =2πk ∫0π sinusin[k(u -esinu)] du =2k Jk′ (ke)$

where endpoint contributions to the integral are again zero. The complete expansion of the cosine of the eccentric anomaly is

$cosu=−e2 +2∑ k=1∞ Jk′ (ke)k coskωt$

Having expansions of both the sine and cosine of the eccentric anomaly immediately gives expansions for the radial variable and its Cartesian components. The expansion of the radial variable is

$ra =1+e22 -2e∑ k=1∞ Jk′ (ke)k coskωt$

while those of its two components are

$xa =−3e2 +2∑ k=1∞ Jk′ (ke)k coskωt$

$ya =21 -e2e ∑k=1 ∞ Jk (ke)k sinkωt$

An expansion of the inverse of the radial variable is particularly easy because of cancellations of linear factors. The constant term in its cosine expansion is

$1π ∫0π ar dτ =1π ∫0π du =1$

while the remaining coefficients are

$2π ∫0π ar coskτ dτ =2π ∫0π cos[k(u -esinu)] du =2Jk (ke)$

so that the complete expansion is

$ar =1 +2∑ k=1∞ Jk(ke) coskωt$

The same cancellations occur in expansions for the sine and cosine of the polar angle in the invariant plane. The constant term in the cosine expansion is

$1π ∫0π cosφdτ =1π ∫0π (cosu-e) du =−e$

and the remaining coefficients are

$2π ∫0π cosφ coskτ dτ =2π ∫0π (cosu-e) cos[k(u -esinu)] du =2(1 -e2)e Jk(ke)$

so that the complete expansion is

$cosφ=−e +2(1 -e2)e ∑k=1 ∞ Jk(ke) coskωt$

This expression can also be found using the relationship

$cosφ=−1e +1-e2 e ar$

between dynamic variables and the previous result. The coefficients for the corresponding sine expansion are

$2π ∫0π sinφ sinkτ dτ =2π ∫0π 1-e2sinu sin[k(u -esinu)] du =21-e2 Jk′ (ke)$

with the result

$sinφ=21 -e2 ∑k=1 ∞ Jk′ (ke) sinkωt$

An explicit expansion of the polar angle φ itself is possible, but it is much more complicated than either of the last two expressions and will not be included.

With the basic dynamic variables in place, one can begin to consider expansions for some of their combinations. First consider a cosine expansion of the square of the sine of the eccentric anomaly. The constant term is

$1π ∫0π sin2u dτ =1π ∫0π sin2u (1-ecosu) du =12$

and the remaining coefficients are

$2π ∫0π sin2ucoskτ dτ =2π ∫0π sin2u cos[k(u -esinu)] (1-ecosu) du =−2πk ∫0π sin2usin[k(u -esinu)] du =1k [Jk+2 (ke) -Jk-2 (ke)] =1k [2(k +1)ke Jk+1 (ke) -2(k -1)ke Jk-1 (ke)] =4k2 e2 Jk (ke) -4ke Jk′ (ke)$

where ${J}_{0}\left(ke\right)$ is added and subtracted on the third line. The complete expansion for this square is thus

$sin2u =12 +4e2 ∑k=1 ∞ Jk (ke) k2 coskωt -4e ∑k=1 ∞ Jk′ (ke)k coskωt$

The expansion for the square of the cosine of the eccentric anomaly follows from the simple relationship between the two circular functions:

$cos2u =12 -4e2 ∑k=1 ∞ Jk (ke) k2 coskωt +4e ∑k=1 ∞ Jk′ (ke)k coskωt$

Another combination that follows from the integal in the last evaluation is a sine series for the product of the sine of the eccentricity and the cosine of the angular variable in the invariant plane. The coefficients are

$2π ∫0π sinucosφ sinkτ dτ =2π ∫0π sinu(cosu -e) sin[k(u -esinu)] du =−2e Jk′ (ke) +1π ∫0π sin2u sin[k(u -esinu)] du =−2 ke2 Jk (ke) +2(1 -e2)e Jk′ (ke)$

with the result

$sinucosφ =−2e2 ∑k=1 ∞ Jk (ke)k sinkωt +2(1 -e2)e ∑k=1 ∞ Jk′ (ke) sinkωt$

The product of the sine of the eccentricity and the sine of the angular variable is a cosine series whose constant term is

$1π ∫0π sinusinφ dτ =1π ∫0π 1-e2 sin2u du =12 1-e2$

The remaining coefficients are

$2π ∫0π sinusinφ coskτ dτ =2π ∫0π 1-e2 sin2u cos[k(u -esinu)] du =1-e2 Jk (ke) -1π ∫0π cos2u cos[k(u -esinu)] du =1-e2 [Jk (ke) -12 Jk+2 (ke) -12 Jk-2 (ke)] =1-e2 [2Jk (ke) -(k+1) ke Jk+1 (ke) -(k-1) ke Jk-1 (ke)] =1-e2 [−2(1 -e2) e2 Jk (ke) +2ke Jk′ (ke)]$

so that the complete result is

$sinusinφ =1-e2 [12 -2(1 -e2) e2 ∑k=1 ∞ Jk (ke) coskωt +2e ∑k=1 ∞ Jk′ (ke)k coskωt]$

Finally, since the equations of motion in Cartesian coordinates are

$d2x dt2 =−km xr3 =−ω2 a3 cosφr2 d2y dt2 =−km yr3 =−ω2 a3 sinφr2$

one can differentiate existing expansions for the immediate results

$a2r2 cosφ =2∑ k=1∞ kJk′ (ke) coskωt$

$a2r2 sinφ =21 -e2e ∑k=1 ∞ kJk (ke) sinkωt$

The direct evaluation of these expansions is most easily accomplished with an intermediate integration by parts using the integrals

$∫du cosu-e (1-ecosu )2 =sinu 1-ecosu ∫dusinu (1-ecosu )2 =−1e (1-ecosu)$

where endpoint contributions to integrals are zero in both cases. The coefficients for the cosine expansion are

$2π ∫0π a2r2 cosφcoskτ dτ =2π ∫0π cosu-e (1-ecosu )2 cos[k(u -esinu)] du =2π ∫0π ksinu sin[k(u -esinu)] du =2kJk′ (ke)$

while the coefficients for the sine series are

$2π ∫0π a2r2 sinφsinkτ dτ =2π ∫0π 1-e2 sinu (1-ecosu )2 sin[k(u -esinu)] du =2π 1-e2e ∫0π kcos[k(u -esinu)] du =2k1 -e2e Jk (ke)$

Both are identical as expected with the results from the equations of motion. A third way to get to the cosine expansion is via partial differentiation of the expansion of the inverse radial variable:

where defining differential equation for Bessel functions has been used to replace the second derivative. This last somewhat surprising evaluation raises the question of what relationships hold among the various Kapteyn series appearing in this presentation and the extent to which such relationships can produce further expansions.

To facilitate comparisons of expressions, introduce the notation for the series

$Sn=∑ k=1∞ knJk (ke) sinkτ Cn=∑ k=1∞ knJk (ke) coskτ Sn′ =∑k=1 ∞ kn Jk′ (ke) sinkτ Cn′ =∑k=1 ∞ kn Jk′ (ke) coskτ$

where a prime indicates the series is constructed from derivatives of Bessel functions and $\tau =\omega t$ as before. Derivatives with respect to this mean anomaly increase the index while changing the circular function:

$∂Sn ∂τ =Cn+1 ∂Cn ∂τ =−Sn+1 ∂Sn′ ∂τ =Cn+1 ′ ∂Cn′ ∂τ =−S n+1′$

Derivatives with respect to eccentricity increase the index and alternate whether the series is constructed from derivatives or not:

$∂Sn ∂e =Sn+1 ′ ∂Cn ∂τ =Cn+1 ′ ∂(e Sn′) ∂e =1-e2 eSn+1 ∂(e Cn′) ∂e =1-e2 eCn+1$

The defining differential equation for Bessel functions has been used here to replace second derivatives.

In terms of this notation, simple dynamical variables are

and their combinations are

Algebraic relationships among the series can be determined from the algebraic relationships of these variables. For example, multiplying the radial variable and its inverse gives

Equating the expressions for $racosφ$ and $xa$ produces the same result. Using expressions for the y-coordinate one has

This relationship arises from adding the squares of circular functions,

but there does not appear to be a simple way to determine each of these squared series separately. This can be seen by beginning to evaluate the expansion for the squared inverse of the radial variable, which would be needed for a separate expression for ${C}_{0}^{2}$ . The constant term is a known integral,

$1π ∫0π a2r2 dτ =1π ∫0π du1 -ecosu =11 -e2$

but the remaining coefficients are not so simple:

$2π ∫0π a2r2 coskτdτ =2π ∫0π cos[k(u -esinu)] 1-ecosu du$

The denominator can be expanded as an infinite geometric series and rearranged to terms of the form $cosmu$ for explicit evaluation, but that means that each of these coefficients is itself an infinite series. In fact this same infinite series appears for the coefficients of the expansion of φ mentioned above but not evaluated.

There are clearly some simplifications among products of these interesting series that are akin to the product formulae of circular functions. To be continued...

Uploaded 2018.08.30 — Updated 2020.11.15 analyticphysics.com