On the Bessel solution of Kepler’s equation

Since its introduction in 1650, Kepler’s equation has never ceased to fascinate mathematicians, scientists, and engineers. Over the course of five centuries, a large number of different solution strategies have been devised and implemented. Among them, the one originally proposed by J. L. Lagrange and later by F. W. Bessel still continue to be a source of mathematical treasures. Here, the Bessel solution of the elliptic Kepler equation is explored from a new perspective offered by the theory of the Stieltjes series. In particular, it has been proven that a complex Kapteyn series obtained directly by the Bessel expansion is a Stieltjes series. This mathematical result, to the best of our knowledge, is a new integral representation of the KE solution. Some considerations on possible extensions of our results to more general classes of the Kapteyn series are also presented.


I. INTRODUCTION
There is no other equation in the history of science that boasts more solution strategies than the famous Kepler equation (hereafter KE).Since 1650, a huge number of different strategies and algorithms have been devised and implemented to deal with KE.The classic textbook by Colwell [1] lists and briefly describes most of these methods.However, since 1993 (the year of publication of Colwell's book), this list has has grown so rapidly that it is difficult to account for all of the strategies that have been developed.Just to gain an idea, quoting some of the papers that have been published on the subject in the last five years would be sufficient [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17].
In this paper, we focus on a semi-analytic approach to the solution of KE, originally conceived by J. L. Lagrange and later finalised by F. W. Bessel via a Fourier series expansion, in which his famous functions appeared.The whole third chapter of [1] is devoted to an account of the history of the infinite series solutions of KE, starting with an important theorem proven by Lagrange in 1770, to the series expansion proposed by T. Levi-Civita in 1904.Today, the Bessel solution of KE has been abandoned for practical, computational purposes, although it continues to arouse some interest, especially regarding its mathematical properties.The Bessel solution is the first example of a class of expansions called the Kapteyn series (hereafter KS) [18], which has gained considerable importance in mathematics and theoretical physics, even in recent times [19][20][21][22][23][24][25][26][27][28][29][30].
In the present paper, the Bessel solution of KE is studied from the perspective of the theory of so-called Stieltjes functions [31][32][33].In particular, two main results are established: (i) a Kapteyn series obtained by a simple "complexification" of the Bessel series is proven to be Stieltjes and (ii) thanks to such a (constructive) proof, to our knowledge, a new integral representation of the KE solution will be obtained.We also hope that, under this perspective, the availability of a well-established convergence theory for Stieltjes series based on the Padé approximant [32,34], could provide new life to the Bessel solution, even as a computational tool.

II. BESSEL' SOLUTION OF ELLIPTIC KEPLER'S EQUATION
Consider the motion of a planet around the star S along an elliptic orbit with eccentricity ϵ, as sketched in Figure 1 [35].
The geometry of Kepler's Equation.

arXiv:2312.01437v2 [math-ph] 4 Jan 2024
Suppose P is the position of the planet at time t, which is measured by assuming the pericenter A as the initial (i.e., at t = 0) position (motion is counter-clockwise).Upon introducing the mean angular speed ω = 2π/T , with T being the orbital period, the so-called mean anomaly, say M , is defined by M = ωt.To retrieve the position of the planet along the orbit at time t > 0, the angle ψ is first evaluated by solving the following trascendental equation: which provides the position of point Q along the circumscribed circle (dotted in the figure).The planet's position P is then immediately obtained through the geometrical construction sketched in Figure 1.Equation ( 1) is called the elliptic KE.The Bessel solution of Equation ( 1) is provided by [1] (Ch. 3) where function S(ϵ; M ) represents the Fourier series Although the series converges for any ϵ ∈ [0, 1), its convergence turns out to be extremely slow, especially when ϵ → 1.This fact, together with the considerable number of Bessel function evaluations to be implemented, unavoidably resulted in Equations ( 2) and ( 3) being abandoned as far as the practical use of KE was concerned.Nevertheless, the Bessel series expansion (3) remained a subject of considerable interest in both mathematics and theoretical physics.The scope of the present paper is to explore why this occurred.To this end, we introduce the complex function S(ϵ; M ), defined through the Kapteyn series in such a way that S = Im{S}.The convergence of the series (4) can be proven, for example, by invoking the following Bessel function asymptotics, valid for ϵ < 1 [36] (Sec.8.4): where χ = √ 1 − ϵ 2 denotes the ellipse's aspect ratio and λ is a negative parameter, provided by The parameter λ will play a key role in the rest of the paper.After substituting Equation (5) into Equation (4) and upon introducing the complex parameter z = exp(λ + iM ), it can immediately be seen that the single terms of the series (4) asymptotically approach those of the paradigmatic model series which, for |z| < 1, is nothing but the Dirichlet series representation of the so-called polylogarithm function L 3/2 (z), where For |z| ≥ 1, the series in Equation ( 8) can be analytically continuated to the whole complex plane, besides the half-axis Re{z} > 1.
More importantly (for the scopes of the present paper), functions L ν (z) are examples of Stieltjes function.For the reader's convenience, the main definitions of Stieltjes functions and of Stieltjes series will now be briefly summarized.Interested readers can found some more useful references for instance in Ref. [37].Moreover, the formalism of Allen et al. [38] will be adopted.Consider then an increasing real-valued function µ(t) defined for t ∈ [0, ∞], with infinitely many points of increase.The measure dµ is then positive on [0, ∞].Assume all moments to be finite.Then, the formal power series is called a Stieltjes series.Such series turns out to be asymptotic, for z → 0, to the function F (z) defined by which is called Stieltjes integral.To give a significative example, the Stieltjes nature of functions L ν will now be proved.To this end, it is sufficient to start from the following integral representation [39] (25.12.11): which, on setting t = exp(−x), gives thus proving that the function L ν (z)/z is a Stieltjes integral exactly of the form (11), with the measure where symbols Γ(•) and Γ(•, •) denote the gamma and the incomplete gamma functions, respectively [39].The above described "asymptotic connection" between series (4) and (7) via Equation ( 5) would at first sight suggest the function S(ϵ; M ) could be Stieltjes too.In the next section such a conjecture will definitely be proved.

III. A CONSTRUCTIVE PROOF THAT S(ϵ; M ) IS A STIELTJES SERIES
We are now going to prove the following theorem: Function S(ϵ; M ) defined through the series (4) is a Stieltjes function.
First of all, the definition of S(ϵ; M ) given in Equation ( 4) must be recast as follows: with z = exp(λ + iM ).What we are going to show is that, similarly as for Equation ( 13), S can be written as where the measure dµ = ρ(t)dt represents the proof goal.Accordingly, the following infinite system must necessarily hold: The key tool of our proof is Watson's integral representation of J n (nϵ), namely [40] J where Function F (θ; ϵ) has features which reveal very important for our scopes.One of them is which suggests Equation ( 18) to be recast as where function Another property of F (θ; ϵ) is that ∂F (θ; ϵ)/∂θ > 0 for θ ∈ (0, π] [40].Then, in order for Equation (21) to be "tuned" with Equation ( 17), the new integration variable t = exp[− G χ (θ)] is introduced, in such a way that for all values of θ except for θ = 0 and θ = π, where t ′ (θ) = 0. Accordingly, function t = t(θ) can be uniquely inverted into the function θ = θ(t), formally defined by where the symbol F −1 denotes the inverse of F (θ; 1 − χ 2 ) with respect to θ.Although its explicit analytical expression cannot be found, θ(t) can be easily numerically computed for any values of t and χ, for instance on using Mathematica's native command InverseFunction.In Figure 2, the behaviour of θ(t), evaluated according to Equation (23), is plotted for different values of the aspect ratio χ.Now, we are going to prove that the function θ(t) actually coincides (apart some constant proportional factors), with the density ρ(t) into Equation (16).To this end, it is sufficient to change the integration variable in Equation ( 21) from θ to t, which yields On performing the t-integration by parts we have or, equivalently, which, once compared to Equation ( 17), leads to ρ(t) = 2θ(t)/π.The uniqueness of the function ρ(t) can further be proved by invoking Calerman's theorem for the Hausdorff problem, i.e., by showing that In order to prove the divergence of the above series, it is worth starting from the following inequality, proved by Siegel [41]: according to which we have and then which leads to Equation (27).Our proof is now complete.S(ϵ; M ) is a Stieltjes function.

IV. A NEW INTEGRAL REPRESENTATION OF KE'S SOLUTION
The proof given in Section III is the main result of the present paper.Here, an interesting byproduct of such proof is now illustrated.In particular, the "Stieltjesness" of S(ϵ; M ) will be employed to conceive a new, up to our knowledge, integral representation of the KE solution.To this end, Equation ( 16) is first recast as follows: so that formal partial integration gives The change of variable t = exp[−G χ (θ)] then yields which can be rearranged into a simpler form with a few work.First of all, we have or, equivalently, To retrieve the KE solution, it is sufficient to recall that S(ϵ; M ) = Im{S(ϵ; M )}, which gives at once that can also be recast as where the term 2π in Equation ( 36) must be removed, due to the phase wrapping of the principal branch of the natural logarithm, −π < arg{log z} ≤ π.Equation ( 37) is the other main result of the present paper.Compared with other, formally similar integral representations of the KE solution, for instance that in [24], Equation ( 37) is expected to provide computational advantages due to the monotonic decreasing behaviour of the function exp[−F (θ; ϵ)] within the integration interval θ ∈ [0, π].
Just to provide an idea, a single but significant numerical experiment will now be carried out by using only native commands of Mathematica.A case of particular interest is that of nearly parabolic orbits around the periapsis, which correspond to ϵ → 1. Accurate numerical solutions for such a critical regime can be found, for instance, in [12,42].In particular, in [42], it was emphasized how solving the elliptic KE for a given couple (ϵ, M ) resulted in a bad conditioned problem in the neighbourhood of (1, 0).To provide an example of the computational use of the new integral representation in Equation ( 37), we considered the sole case of a parabolic orbit, i.e., ϵ = 1.
In Figure 3, the relative error, evaluated with respect to the "exact KE solution", is plotted against M ∈ (0, π).In the above experiments, such an "exact value" was evaluated by solving Equation ( 1) via Mathematica's native command FindRoot with the parameter WorkingPrecision set to 50 and an initial guess of ψ = π.Function S(ϵ; M ) was evaluated by implementing Equation ( 36

V. DISCUSSIONS
Here, the classical solution of Kepler's equation was revisited from a new perspective, offered by the beautiful theory of the Stieltjes series.In particular, after introducing a "complexified" version of the original Kapteyn series (3), its Stieltjes nature was mathematically proven.Our principal tool was Watson's integral representation of Bessel functions.As a potentially interesting by-product of our analysis, to our knowledge, a new integral representation of the KE solution was achieved.
Exploring the Stieltjes nature of Equation ( 4) could also reveal new interesting aspects concerning the more general framework of the Kapteyn series.The need to develop new methods to achieve the summability of the Kapteyn series of both the first and second kind have already been invoked in the recent past [25].In this respect, we believe the results presented here could be of some help.
Just to gain an idea, consider again the Kapteyn series m≥1 z m m J m (mϵ), but suppose that |z| > exp(−λ).In that case, the series diverges.However, the analysis carried out in Section III is also still valid for z belonging to the whole complex plane, provided that |arg(z)| < π.In other words, our Kapteyn series can be continued to a Stieltjes function according to which lives, for ϵ ∈ [0, 1), on the whole complex plane, besides the infinite interval (1, ∞).
Proving that a divergent power series is a Stieltjes series also implies that such a series would be Padé summable.To quote an important example, in [43], numerical evidence for the factorially divergent perturbation expansion for energy eigenvalues of the PT-symmetric Hamiltonian H(λ) = p 2 + 1/4x 2 + iλx 3 being a Stieltjes series and thus Padé summable, were provided.Such conjecture was later rigorously proven in [44].The Padé summability of our Kapteyn series also implies that suitable resummation algorithms, the most well-known being Wynn's ϵ algorithm, can successfully be employed to retrieve the correct value provided in Equation (38).Padé approximants remain the most favourite mathematical tool to resume divergent and slowly convergent series in theoretical physics.This is because a solid convergence theory has been established for them, especially in the case of Stieltjes asymptotic series [31][32][33].In particular, if the input data are the partial sums of a Stieltjes series, it can be rigorously proven that certain subsequences of the Padé table converge to the value of the corresponding Stieltjes function [32,34] (Chapter 5).A typical feature of Wynn's ϵ algorithm is that only the input of the numerical values of a finite substring of the partial sequence of the series is required to achieve a meaningful resummation.
In the case of the Stieltjes series, however, important additional structural information about the behaviour of the series terms is available.For example, it is known that the magnitudes of the truncation errors are bounded by the first term neglected in the partial sum, and they also have the same sign patterns.In the seminal 1973 paper by Levin [45], a new class of sequence transformations was introduced that used such valuable a priori structural information to improve the efficiency and convergence speed of the transformed sequences.For a general discussion on the construction principles of Levin-type sequence transformations, readers are encouraged to go through the paper by Weniger [46], as well as through Ref. [47], which contains a slightly upgraded and formally improved version of the former.To provide a simple example, Table I shows the performance of a particular Levin-type transformation, the so-called Weniger δ-transformation [46], in the resummation of the Kapteyn series on the left side of Equation (38) for ϵ = 9/10 and z = 10 exp(iπ/3).
The convergence to the value obtained by numerically computing the integral on the right side of Equation ( 38) is evident.Finally, before concluding, it is worth illustrating at least one experiment of the use of the Weniger transformation to solve Kepler's equation.Once again, the parabolic case will be analysed, but now only "small values" of M , namely less than 1/10, will be considered.In Figure 4, the relative error is still plotted versus M ∈ 1 1000 , Again, the "exact solution" was evaluated by solving Equation (1) via Mathematica's command FindRoot with the parameter WorkingPrecision set to 50 and an initial guess of ψ = π, but now the function S(ϵ; M ) is thought of as the imaginary part of S(ϵ; M ).The latter is computed, as done in

VI. CONCLUSIONS
In common with almost any scientific problem which achieves a certain longevity and whose literature exceeds a certain critical mass, the Kepler problem has acquired luster and allure for the modern practitioner.Any new technique for the treatment of trascendental equations should be applied to this illustrious case; any new insight, however slight, lets its conceiver join an eminent list of contributors.
Nothing can illustrate the fascination that Kepler's equation has held for nearly 400 years better than the above paragraph, which appears in the introduction to Colwell's book.Kepler's problem continues to inspire mathematicians, scientists, and engineers to search for further contributions to the subject.In the present paper, the Bessel solution of KE has been studied from the perspective of Stieltjes function theory.In particular, it has been proven that the Kapteyn series (3) obtained by the "complexification" of the Bessel series solution (2) is a Stieltjes series.Moreover, thanks to our constructive proof, to our knowledge, a new integral representation of the KE solution has been obtained.More importantly, the limit of validity of such a representation extends far beyond that fixed by the convergence circle of the Kapteyn series itself.
Stieltjes series are strictly connected with Padé summability via a well established convergence theory.Divergent Stieltjes series can be resummed to the corresponding Stieltjes integral values using Padé approximants.In the last 40 years, Levin-type transformations have been shown to possess retrieving capabilities considerably superior to those of Padé, as far as the resummation of the Stieltjes series is concerned.At present, the main theoretical drawback of Levin-type transformations is the lack of a rigorous convergence theory, like that available for Padé.Thus, we also hope that what is contained in the present paper could provide a further contribution to the construction of such a convergence theory, also in order to provide Levin-type transformations the scientific recognition they certainly deserve, but that, unfortunately, in our opinion, still seems too far away.
) using the native Mathematica command NIntegrate with different degrees of accuracy, measured by the parameter WorkingPrecision, which was set to 10 (a), 15 (b), 20 (c), and 25 (d).

FIG. 3 :
FIG.3: Behaviour of the relative error against M ∈ (0, π).In the above experiments, the "exact value" of the KE solution was evaluated by solving Equation (1) via Mathematica's native command FindRoot with the parameter WorkingPrecision set to 50 and the initial guess of ψ = π.Function S(ϵ; M ) was evaluated by implementing Equation (36) through the native Mathematica command NIntegrate with different degrees of accuracy, measured by the parameter WorkingPrecision, which was set to 10 (a), 15 (b), 20 (c), and 25 (d).
Table I, via the Weniger δ-transformation with different orders.
FIG. 4: The same as in Figure 3, but for M ∈ 1 1000 , 1 10 .Note that now the function S(ϵ; M ) is thought of as the imaginary part of S(ϵ; M ), which is computed, similarly to that in Table I, via the Weniger δ-transformation with an order of 20 (black circles), 30 (open circles), 40 (open squares), and 50 (black squares).