The formula that exists, and is not used

An earlier rung of this ladder said that the position of a body on its orbit has no formula and is computed anyway. That is the sentence every textbook uses, and it is wrong in a way worth an essay. Kepler’s equation

$$M = E - e \sin E$$

does have a closed-form solution for $E$ in terms of $M$ and $e$. It was written down by Lagrange in 1770, it is exact, it is an infinite series in the eccentricity, and every term of it can be produced by a rule. What it does not have is a radius of convergence large enough to be useful, and the boundary is a specific number that no orbit knows about.

What the series is #

The equation is transcendental in $E$, which is what stops it being inverted by algebra. But it is analytic in $e$, and that is a different property with different consequences. At $e = 0$ the answer is $E = M$; the question is how the answer moves as the eccentricity is turned up from there.

Lagrange’s inversion theorem answers exactly that question for any equation of the form $E = M + e\,\phi(E)$, and Kepler’s is one with $\phi = \sin$. It gives

$$E = M + \sum_{m=1}^{\infty} \frac{e^m}{m!}\,\frac{d^{m-1}}{dM^{m-1}}\bigl(\sin^m M\bigr),$$

and the first two terms are worth writing out because they are recognisable:

$$E = M + e\sin M + \tfrac{1}{2}e^2 \sin 2M + \cdots$$

The first correction is the one everybody meets — a body runs ahead of the mean near perihelion and behind it near aphelion, by an amount proportional to the eccentricity. The second is the first hint that the correction is not a sine wave. Keep going and the coefficients become trigonometric polynomials of rising order, each an exact rational combination of sines of multiples of $M$.

The coefficients, computed rather than differentiated #

Repeatedly differentiating $\sin^m M$ is exact and unpleasant. There is a shorter route through a second classical result: the Fourier expansion of the same function,

$$E = M + 2\sum_{n=1}^{\infty}\frac{1}{n}J_n(ne)\sin nM,$$

where $J_n$ is the Bessel function of the first kind. Expanding each $J_n(ne)$ in its own power series and collecting powers of $e$ gives the Lagrange coefficients directly: for $m = n + 2k$,

$$c_m(M) = 2\!\!\sum_{\substack{n\le m\\ n\equiv m\ (2)}}\!\! \frac{(-1)^{(m-n)/2}\,n^{m-1}\sin nM}{2^m\,\bigl(\tfrac{m-n}{2}\bigr)!\,\bigl(\tfrac{m+n}{2}\bigr)!}.$$

That expression gives $\sin M$ at $m = 1$ and $\tfrac12\sin 2M$ at $m = 2$, which is the check the generator behind these figures makes before it draws anything. It is also, incidentally, why Bessel functions exist at all: Bessel was working on this problem, and the functions now used for drumheads and waveguides arrived as the coefficients in a planetary series.

Where it stops #

A power series in $e$ has a radius of convergence, and the radius is set by the nearest singularity of the function in the complex $e$ plane — a place the eccentricity of a real orbit can never go. For $E(e, M)$ that singularity is where $\partial M/\partial E = 1 - e\cos E$ vanishes at the same time as the map ceases to be locally invertible, and working the condition out gives an equation with no trigonometry left in it:

$$e\,\exp\sqrt{1+e^2} = 1 + \sqrt{1+e^2}.$$

Its root is $e_{\mathrm L} = 0.6627434193\ldots$, the Laplace limit. Below it the series converges for every $M$; above it, for no $M$ at all except the two fixed points at $0$ and $\pi$.

The geometric rate is the most useful fact here. At the Earth’s eccentricity of 0.0167 the ratio is 0.025, so each term buys about a factor of forty: three terms reach one part in $10^5$ and six reach machine precision. At Mercury’s 0.2056 the ratio is 0.31 and about thirty terms are needed for the same. At 0.6 the ratio is 0.905, and the number of terms required for ten significant figures is around two hundred and thirty. The series does not fail suddenly at the limit; it becomes uneconomic long before, and then fails.

The number is not about orbits #

It is worth being explicit about what $e_{\mathrm L}$ is a property of. It is not a property of ellipses: an orbit at 0.66 and an orbit at 0.67 are geometrically indistinguishable to any measurement, have the same dynamics, obey the same three laws, and sweep the same equal areas. It is not a property of the physical problem: nothing observable changes at that eccentricity, and no comet has ever done anything peculiar on crossing it.

It is a property of one method applied to the problem — the method of expanding in a small parameter and hoping the parameter stays small. The Laplace limit is where that hope runs out, and it runs out at a place decided by a singularity in a plane that contains no orbits.

What is done instead #

Every ephemeris in use solves the equation by iteration, and the reason is now stateable precisely: iteration has no radius of convergence to run out of.

The comparison is not close. The series wants two hundred terms at $e = 0.6$ and infinitely many at $e = 0.7$; Newton wants four steps at either, each costing one sine and one cosine. The series is the older answer and the better one to have thought of, and it is not the one anybody runs.

There is a second reason, and it matters more for the objects the series is worst at. Nearly all the bodies with eccentricities above the Laplace limit are comets, and a comet’s orbit is not merely eccentric but sometimes not an ellipse at all. A series in $e$ built around the ellipse has nothing to say about $e = 1.0001$; the universal formulation has the same thing to say about it as about everything else. The series is not just slow at the interesting end, it stops being about the right object.

The eccentricities that actually occur #

It is fair to ask how much of the solar system the series would have handled, since somebody had to make that judgement in 1770 without a computer.

All eight planets are far below the limit; the largest is Mercury at 0.2056, where about thirty terms give ten digits, and that was within reach of a person with a table of sines and a winter. The Moon is at 0.055. The asteroids are mostly below 0.3. So for the bodies whose positions anybody needed in the eighteenth century, the series was a genuine method and was used as one — Lagrange’s and Laplace’s planetary theories are built out of expansions of exactly this kind, and the equation of the centre, which is $\nu - M$ rather than $E - M$ but has the same structure, is tabulated in series form in every nautical almanac of the period. What defeated the method was the comets, which is also what made it worth defeating. Encke’s comet is at 0.85, Halley’s at 0.967, and the long-period comets crowd against 1. Those are the objects whose returns are worth predicting and whose orbits the series cannot touch.

What a series is still for #

Iteration wins on speed and on range, and it loses one thing that matters in a narrower place than it sounds: a series is an expression, and an expression can be differentiated.

Orbit determination and perturbation theory both need the derivatives of position with respect to the six elements, because that is what a least-squares correction is built out of and what a first-order perturbation is expanded in. A Newton iteration returns a number and no derivative; getting $\partial E/\partial e$ out of it means either differentiating the equation implicitly — which is easy here, $\partial E/\partial e = \sin E/(1 - e\cos E)$ — or differencing, which loses digits. A series returns the whole function of $e$ at once, and every derivative with respect to $e$ is another series with the same convergence. That is why nineteenth-century celestial mechanics is written in series and not in iterations, and why the literature on the Laplace limit is largely nineteenth-century. When the object of study is not “where is this body tonight” but “what does a small change in Jupiter’s mass do to Saturn’s longitude over a thousand years”, the thing needed is an expression in the parameters, and the eccentricities involved are all comfortably small. The elements that drift are analysed by exactly that route.

The modern division of labour is therefore not that the series lost. It is that the two questions separated: numerical ephemerides iterate, and analytic theories expand, and the Laplace limit is a constraint on the second which the first never meets.

Both of the points above are about the same asymmetry: a closed form is judged by what it costs to evaluate and by what can be done with it, and an iteration is judged by whether it can be relied on to start. Neither judgement is about whether a solution exists, which is the question the title of this essay is about and the one that turns out to matter least.

What the ladder has established #

Three rungs of this anchor have now said three different things about one curve. The first is that the orbit is an ellipse with the primary at a focus, which fixes the shape. The second is that the timing along that shape is set by a transcendental equation. The third is that an orbit can be indistinguishable from a circle and still not be one, which is a statement about how little eccentricity it takes to matter.

This one adds the correction to the second. The position is not uncomputable in closed form; it is computable in closed form over part of the range, and the part is bounded by a number that belongs to complex analysis rather than to celestial mechanics. That distinction is not pedantry — it is the difference between “no formula exists” and “the formula that exists has a domain”, and the second is a statement somebody can act on. Acting on it produced the universal variables, which are what happens when the problem is re-posed so that no expansion parameter is required at all.

The figure this essay cannot draw is the one that would settle it visually: a picture of the singularity, which lives at a complex eccentricity of about $0.66 + 0.68i$ and is not on any page that also shows an orbit. The two panels of the opening figure are the closest available, and what they show is the consequence rather than the cause — a sum that closes and a sum that does not, from an equation that looks the same in both.

The first of them is that the essay’s title is only true of one of the two series Lagrange’s contemporaries wrote down, and the other one has no such limitation.

The series that does converge #

The Laplace limit is a statement about one series, and the essay has already written down a second one without remarking that it behaves differently.

The Fourier form,

$$E = M + 2\sum_{n=1}^{\infty}\frac{1}{n}J_n(ne)\sin nM,$$

is not a power series in the eccentricity. It is an expansion in harmonics of the mean anomaly, whose coefficients happen to be functions of the eccentricity — and it converges for every eccentricity below one, at every mean anomaly, with no Laplace limit anywhere in it.

That is not a contradiction. The two series are different objects: expanding each Bessel function in powers of $e$ and collecting terms is a rearrangement, and rearranging a series can destroy its convergence. What the power series is asking is how the answer varies as $e$ is turned up from zero, and that question has a radius of convergence set by a singularity off the real axis. What the Fourier series is asking is how the answer decomposes into harmonics at a fixed eccentricity, and that question has no such obstruction.

So there is a convergent closed-form solution to Kepler’s equation for every ellipse, and it is a hundred and eighty years old.

It is also not used, for a reason that is arithmetic rather than principled. The Bessel coefficients decay slowly as the eccentricity approaches one — the decay is governed by how fast $J_n(ne)$ falls with $n$, and for $e$ near unity that is a slow algebraic fall rather than a geometric one. At $e = 0.9$ something like fifty terms are needed for six digits, and each term costs a Bessel function evaluation, which is far more expensive than the sine and cosine a Newton step needs.

The comparison is therefore not between a formula that fails and an iteration that works. It is between a formula that works everywhere and costs a hundred special-function evaluations, and an iteration that costs four trigonometric ones.

A closed form is worth having when it is cheap or when it can be differentiated, and this one is neither cheap nor needed, which is why a result of Bessel’s is a curiosity in the field it was invented for.

Two further points belong to the comparison between the two methods, and the first of them undercuts the essay’s own title slightly.

The guess, which is where the difficulty actually is #

Newton’s method converges quadratically once it is close, and nothing guarantees that it gets close. For Kepler’s equation at high eccentricity that is the whole practical problem, and it is worth setting out because it is where the engineering effort has gone.

The naive starting guess is $E_0 = M$, which is exact at zero eccentricity and is what the equation reduces to there. It works comfortably for a planet. At an eccentricity near one and a small mean anomaly it does not: the function’s derivative $1 - e\cos E$ is nearly zero there, so the first Newton step divides by a very small number and throws the iterate a long way — sometimes past $\pi$, occasionally onto a value from which the iteration converges to the answer for a different revolution.

The standard defences are of two kinds.

One is a better starter. Several closed-form approximations exist whose error is small enough everywhere that Newton catches immediately: the simplest useful one is $E_0 = M + e\sin M/(1 - \sin(M+e) + \sin M)$, and there are cubic starters that solve a reduced form of the equation exactly and are accurate to a few per cent over the whole range. With one of those, three iterations suffice at any eccentricity below one.

The other is a safeguarded iteration. Bracketing the root between zero and $\pi$ — which is always valid, since the function is monotonic there — and falling back to bisection whenever a Newton step leaves the bracket gives a method that cannot fail, at the cost of occasionally converging linearly.

Production ephemeris code does both, and the reason is not fastidiousness. A solver embedded in an integration is called millions of times with arguments nobody inspects, and a failure that happens once in ten million calls is a failure that happens.

The transcendental equation is not the hard part; the hard part is a starting value that behaves at both ends of a range, which is the usual shape of a numerical problem once the mathematics has been settled.

Both halves of that division are worth keeping in view, because the phrase “no closed-form solution” is used about a great many problems and rarely means the same thing twice.

Neither of the two properties that decide the question — cost per evaluation, and whether the result can be differentiated — has anything to do with existence.

Where the ladder goes next #

The rung after this one is the equation of the centre proper, $\nu - M$, which is the quantity an almanac tabulates and which has its own series with its own convergence — the same limit, arrived at through a different function. Past that is the question this essay has kept out of the way: what happens to all of it when the orbit is not quite the two-body orbit it was assumed to be, so that $e$ itself is a slowly changing number and the series is being expanded about a moving point.