The elements that stop existing

The six classical orbital elements are one of the most successful sets of coordinates ever devised. Five of them are constants of a two-body motion and the sixth is a clock, so an orbit becomes a point that does not move, and a perturbation becomes a slow drift of a point that would otherwise stay put. Nearly all of celestial mechanics since Lagrange is written in them.

Three of them are also undefined for orbits that are entirely ordinary, and the failure is not a rare edge case: it happens to satellites in commercial use, to planets in the solar system, and to every orbit that is nearly circular or nearly equatorial — which is to say, to most of them.

Which three, and why #

The classical set is $(a, e, i, \Omega, \omega, M)$: size, shape, tilt, the direction of the line where the orbit plane crosses the reference plane, the angle from that line round to pericentre, and the position along the orbit. Each of the last three is an angle measured from a feature, and each feature can be absent.

• $\omega$, the argument of pericentre, is measured from the ascending node to the point of closest approach. A circular orbit has no closest approach, so there is no direction to measure to.
• $\Omega$, the longitude of the ascending node, is measured to the line in which the orbit plane cuts the reference plane. An orbit in the reference plane does not cut it, so there is no line.
• $M$, the mean anomaly, is measured from pericentre, and inherits $\omega$’s problem.

It is worth being precise about the kind of failure. This is not a physical singularity like the collision the regularised two-body problem removes, where something in the equations of motion really does become infinite. It is a coordinate singularity, of the same species as the one at the pole of a latitude–longitude grid: longitude is undefined at the North Pole, an aircraft flying over it passes through every longitude in an instant, and nothing has happened to the aircraft.

The rate is exactly one over the eccentricity #

The blow-up in the opening figure has a closed form, which is what makes it a property of the coordinates rather than a numerical artefact.

Write the eccentricity vector in components $k = e\cos\varpi$ and $h = e\sin\varpi$. Then $\varpi = \operatorname{atan2}(h, k)$ and

$$\dot\varpi = \frac{k\dot h - h\dot k}{e^2}.$$

The numerator is a perfectly finite quantity — it is the rate at which the eccentricity vector sweeps area about the origin, and a secular perturbation makes it constant. The denominator goes to zero. So the pericentre’s angular rate goes as $1/e^2$ times a fixed area rate, which is $1/e$ times a fixed linear speed: as the vector passes at distance $e_{\min}$ from the origin, the angle swings at $Ag/e_{\min}$ for a free amplitude $A$ turning at rate $g$. In the figure that is 3,500° per cycle against a free rate of 360°, and the ratio is the ratio of the amplitudes: $0.034/0.003$.

This is not an exotic configuration. The forced-plus-free decomposition is the standard first-order description of what the planets do to each other, and every planet’s eccentricity vector is a small circle offset from the origin by a forced term. Venus’s eccentricity varies between about 0.0000 and 0.0700 over roughly 400,000 years, and at the low end its pericentre is doing precisely what the top panel of the opening figure does.

The replacement, and what it costs #

The cure is to stop measuring from features that can vanish. The equinoctial elements do it by combining the offending angles with the eccentricity and inclination they belong to:

$$h = e\sin(\Omega+\omega),\quad k = e\cos(\Omega+\omega),\quad p = \tan\tfrac{i}{2}\sin\Omega,\quad q = \tan\tfrac{i}{2}\cos\Omega,$$

with $a$ kept as it is and the position given by the mean longitude $\lambda = M + \omega + \Omega$. Each of the four new quantities goes smoothly to zero as its own singularity is approached, because each is a product of the small quantity and a bounded trigonometric function. There is no angle left to blow up: the information that was in $\omega$ is still there, distributed between $h$ and $k$, and it simply becomes less and less as the orbit becomes rounder, which is the honest accounting.

The price is legibility. An orbital element ought to mean something, and $a$, $e$ and $i$ all do — a size, a shape, a tilt. $h$ and $k$ do not. An engineer reading a state vector in equinoctial elements cannot see at a glance where pericentre is, and the conversion back is a two-line calculation with its own edge cases at $i = 180°$, where $\tan(i/2)$ diverges and a retrograde equatorial orbit becomes the singular case instead.

What was actually measured #

The clearest evidence that this is a real problem rather than a mathematical fussiness is in the published element sets themselves.

A geostationary satellite’s inclination is held below about 0.1° by north–south station-keeping, and its eccentricity below about 0.0003. At those values the reported $\Omega$ and $\omega$ in a two-line element set jump by tens of degrees between consecutive epochs a day apart, while the satellite is doing nothing but sitting still over its longitude. Neither number is an error; each is the argument of a vector shorter than its own uncertainty. What is stable, and what a controller can act on, is the pair $(k, h)$ and the pair $(p, q)$, and those are what appear in the manoeuvre planning.

The same shows up in the planetary theories. The secular solutions for the solar system — Laplace’s, Le Verrier’s, and the modern numerical ones — are all written in $(k, h, p, q)$ rather than in $(e, \varpi, i, \Omega)$, and the reason given in every one of them is the same: the equations in the classical variables have $1/e$ and $1/\sin i$ on the right-hand side, and Mercury is the only planet whose eccentricity is comfortably far from zero.

Where the model stops #

Equinoctial elements remove two singularities and do not remove the third. The mean longitude $\lambda$ still requires a rate, and that rate involves $a^{-3/2}$; an orbit whose semi-major axis passes through infinity — a bound orbit perturbed onto an escape trajectory — is outside the set entirely, as it is outside the classical one. That case is handled by yet another parameterisation, and the pattern by now is familiar: each element set is adapted to a family of orbits and fails at that family’s boundary.

Nor do elements of any kind help with the case where the two-body approximation itself is the problem. Osculating elements are defined by pretending that the current position and velocity belong to an unperturbed Kepler orbit; where the perturbation is not small — inside a close encounter, or for a satellite deep in an atmosphere — that pretence produces elements that gyrate wildly and describe nothing.

The history, which is shorter than the problem #

Lagrange wrote the planetary equations in 1808, and the $1/e$ and $1/\sin i$ terms are there in his own derivation. He knew what they meant, and the workaround for most of the nineteenth century was to avoid the cases: perturbation theories were written for the planets, whose eccentricities are small but not that small, and where a variable did approach zero the theorist changed variables by hand for that one problem.

The general fix arrived with computers, because computers do not change variables by hand. The equinoctial set in its modern form is due to Broucke and Cefola in 1972, working on satellite theory for exactly the orbits — near-circular, near-equatorial — that the classical elements are worst at. The paper’s argument is not that the classical elements are wrong; it is that a program has to run on every input it is given, and a set of variables that requires the user to know in advance which special case applies is not a set of variables, it is a family of them.

That is a recognisable shape. The move from a chart that works for the typical case to one that works for every case is what turns a technique into a piece of software, and it happened to attitude representation at almost exactly the same time and for the same reason.

The sixth element has the same problem #

The three angles are the famous cases, but the mean anomaly inherits the difficulty and is the one that bites in propagation rather than in bookkeeping.

$M$ is measured from pericentre. If pericentre is undefined, so is $M$, and the two are undefined together in a way that partly cancels: the sum $\varpi + M$ — the mean longitude — is measured from the fixed reference direction and is perfectly well behaved. That is why $\lambda = M + \omega + \Omega$ is what the equinoctial set carries, and it is a small piece of good luck rather than a design: the two bad angles enter the position with opposite signs, so their difference is unstable and their sum is not. Three further points belong with the argument rather than after it: what a genuine singularity looks like by comparison, what happens at the other end of the family of orbits, and how a published catalogue copes with any of it.

The singularity that is not a coordinate #

The essay has insisted that this failure is a property of the description, and the contrast is worth drawing against a case where it is not.

Two point masses on a collision course reach zero separation in finite time, and at that instant the acceleration is infinite and the equations of motion have no solution to continue. That is a genuine singularity of the problem, not of the chart, and no relabelling removes the collision.

What can be removed is the numerical difficulty it causes. Long before a collision the separation is small, the acceleration is large, and any integrator taking uniform steps in time is either taking absurdly small steps everywhere or losing accuracy near the closest approach.

The remedy is a change of both the dependent and the independent variable. Replacing the time by a fictitious time whose increments are the true time divided by the separation makes the steps automatically shrink where the motion is fast; and replacing the position by a set of variables whose squares give the position turns the equation of motion into that of a harmonic oscillator, which has no singularity at all.

That is regularisation, and the transformation used for the three-dimensional case has an unexpected structure: it maps three dimensions into four, because the algebra that linearises the problem is the algebra of quaternions rather than of vectors. The extra coordinate is constrained rather than free, and the resulting equations are linear in the fictitious time.

The distinction between the two kinds of trouble is worth carrying because the responses are different. A coordinate singularity is fixed by choosing better variables and the physics is untouched. A physical singularity is not fixed at all — what regularisation buys is a computation that runs through it accurately rather than one that survives it, and whether the collision itself is meaningful is a question about the model rather than about the arithmetic.

The comparison also explains why the two are treated so differently in practice: one is a bug to be fixed once, and the other is a property of the problem that every integrator has to be built around.

Two of the drifts the essay describes are worth drawing at other settings, because each of them is a rate rather than a state and the rates depend on where the orbit is.

The chart for the other end of the family #

Equinoctial elements remove the singularities at zero eccentricity and zero inclination and leave one at the other extreme, and the fix for that is a different construction with the same motivation.

The difficulty is that the classical elements are a chart on the bound orbits. A parabolic orbit has infinite semi-major axis, a hyperbolic one has a negative one, and the mean anomaly’s definition changes form in each case — so a program that integrates a comet through a perturbation carrying it from bound to unbound has to switch parameterisations partway.

The universal-variable formulation removes the switch. Instead of the eccentric anomaly, which is defined differently for each conic, it uses a variable proportional to the accumulated path in a fictitious time, and expresses everything through a family of functions of a single argument that combines that variable with the reciprocal of the semi-major axis.

Those functions are analytic through zero, and their arguments’ sign is what distinguishes the conics — positive for an ellipse, zero for a parabola, negative for a hyperbola. So one set of expressions covers the whole family, evaluates without branching, and is accurate near the parabolic limit, where the classical forms are worst.

The practical importance is exactly the case the classical set handles worst and that occurs constantly: a long-period comet whose eccentricity is 0.99998, or a spacecraft on a trajectory that a gravity assist tips across the parabolic boundary. Both are ordinary, and both are near the place where a two-branch formulation loses digits.

A chart’s boundaries are where the interesting objects are, which is a fair summary of both halves of this essay.

What a catalogue does about it #

The theory says to use different variables. Published element sets mostly do not, and it is worth saying how the difficulty is handled in practice.

The two-line element sets for Earth satellites carry the classical angles, because the format was fixed in the 1960s and is consumed by a great deal of software that would break if it changed. So a geostationary satellite’s file contains an argument of pericentre that jumps by tens of degrees between epochs, and every consumer of that file has to know not to differentiate it.

The convention that mitigates the worst of it is to quote the longitude of pericentre, the sum of the node and the argument of pericentre, rather than the two separately. For a low-inclination orbit the individual angles are both ill-determined and their sum is not, because the ill-determination is a rotation of the reference direction within the orbital plane and it cancels in the sum. The same trick appears in planetary work for the same reason, which is why the planetary tables quote a longitude of perihelion rather than an argument of it.

Small-body catalogues have the opposite problem and do not need the trick: an asteroid’s eccentricity and inclination are rarely small enough for the angles to misbehave, so the classical set is used throughout and nobody notices the singularity exists.

The one place a catalogue does change representation is in the internal working of the fitting software. An orbit determination is a least-squares problem, and a least-squares problem in a variable that is unbounded when its partner is small is ill-conditioned — so the fits are done in equinoctial variables and the answers are converted to classical ones for publication.

The singularity is therefore hidden rather than removed, and the place it shows is in the difference between the numbers a catalogue publishes and the numbers the software that produced them was working in.

And the triaxial case drawn on its own, since it is the one singular configuration that is a physical equilibrium rather than a defect of the chart.

Where this ladder goes next #

The lesson generalises past orbits, and stating it that way is the point of the rung: a coordinate system is a hypothesis about what varies, and it fails where the hypothesis does. The classical elements assume there is a pericentre and a node, and they are superb wherever there is.

The next rungs on this ladder are about what happens when the perturbation stops being small enough to treat as a drift of the elements at all. A distant companion can trade an orbit’s inclination for its eccentricity on a timescale far longer than either period, driving $e$ from zero to near unity — a secular effect that is enormous, and that is naturally written in exactly the variables this essay has been arguing for, because it carries the eccentricity vector clean through the origin more than once.