The series that is subtracted

The two-body problem is finished. It was finished in 1687, it has a closed solution in terms of an ellipse and one transcendental equation, and nothing since has improved on it. That is a peculiar kind of luck, because the problem anybody actually wants to solve — a Sun with eight planets in it, or a satellite around an oblate Earth with a Moon nearby — has no closed solution at all and never will.

Celestial mechanics deals with this by subtraction. It does not attack the real problem; it takes the solved problem away from the real one and studies the difference. The difference has a name, a standard expansion, and a hierarchy among its terms that is the reason the method works rather than merely being a way of writing the equations down.

What is taken away, and what is left holding the equations #

Write the acceleration of a planet around the Sun with one other planet present. There are three attractions in it: the Sun on the planet, the other planet on the planet, and the Sun on the Sun — which sounds absurd until it is remembered that the coordinates are centred on the Sun, and a coordinate system anchored to an accelerating body is not an inertial frame.

Collecting the first of those on its own gives back Newton’s two-body equation exactly, and its solution is an ellipse. The rest is written as the gradient of a single scalar,

$$\mathcal{R} = Gm'\left(\frac{1}{\Delta} - \frac{\mathbf{r}\cdot\mathbf{r}'}{r'^3}\right), \qquad \Delta = |\mathbf{r} - \mathbf{r}'|,$$

which is called the disturbing function, and the whole of planetary theory is the study of that one object. The first term is the direct attraction of the perturber. The second is the indirect term, and it is the price of putting the origin on the Sun.

The indirect term is the part that gets forgotten, and forgetting it is not a small error. It is a pure cosine of the angle between the two bodies, so it is as large as the direct term at most separations and larger at many, and it averages to exactly zero while the direct term does not. Everything slow in a planetary system comes from that asymmetry: what survives averaging is not what dominates the instantaneous pull.

It is worth being concrete about where the indirect term comes from, because the algebra hides an ordinary fact. Jupiter pulls on Saturn, which is the direct term, and Jupiter also pulls on the Sun, which moves the Sun. A Saturn-centred astronomer would see the second effect as a displacement of the Sun and would have no difficulty with it. A Sun-centred astronomer sees it as a fictitious force on Saturn, because the ruler being used is attached to the thing that moved. The term is not an approximation and it is not optional: dropping it changes the perihelion precession of the inner planets by tens of arcseconds a century, which is an error of the same order as the total effect being computed.

The elements the perturbation acts on are the osculating elements: the six numbers of the ellipse the body would follow if every force but the primary’s were switched off at that instant. They are exactly constant in the two-body problem and they are what the disturbing function moves. Gauss’s equations say which component moves which one; the disturbing function says how large the components are and, crucially, how they depend on time.

Why an infinite series is not a defeat #

The trouble with the expression above is that $\Delta$ is a square root of a sum of squares and cosines, which is not a function anybody can integrate against time. So it is expanded — in powers of the ratio of the two semi-major axes, in powers of the eccentricities, and in powers of the sines of half the inclinations — and the result is an infinite sum of terms, each of the form

$$C_{jklm}\,\cos(j\lambda + k\lambda' + l\varpi + m\varpi'),$$

a coefficient multiplying the cosine of a linear combination of four angles. The coefficients involve the Laplace coefficients, which are integrals over the semi-major axis ratio; the angles are the mean longitudes and the longitudes of perihelion.

An infinite series looks like a step backwards from a closed solution, and it would be, except that the terms are not all the same size. They are stacked, and the stacking is a theorem.

D’Alembert’s rule is the theorem: a term whose argument contains $l\varpi + m\varpi'$ carries a factor $e^{|l|}e'^{|m|}$, and one containing the nodes carries the corresponding powers of the inclinations. The rule is not empirical. It follows from the fact that the physical quantities are analytic at $e = 0$, where the perihelion is undefined, and a function that must stay finite where an angle becomes meaningless can only depend on that angle multiplied by enough powers of the quantity that makes it meaningful.

That is the whole reason planetary perturbation theory exists as a practical subject. Earth’s eccentricity is 0.017, so a fourth-order term is $10^{-7}$ of a first-order one. Truncating after second order is not an act of faith; it is arithmetic with a known remainder.

The rule also explains a fact about the tables that would otherwise be a curiosity. Le Verrier’s theory of the inner planets ran to some four hundred terms; the theories of the outer planets, at comparable accuracy, needed fewer. The outer planets are further apart in the sense that matters — their ratios of semi-major axes are more favourable — but they are also on rounder orbits, and both of the two expansions are therefore working in their easy regime at once. Mars is the awkward case for the same two reasons run backwards: an eccentricity of 0.093 and a neighbour at a ratio of 0.66.

The other expansion, and the one that breaks #

The eccentricity expansion is not the only one. The direct term has to be expanded in the ratio of the two semi-major axes as well, and that expansion is the generating function of the Legendre polynomials:

$$\frac{1}{\Delta} = \frac{1}{a'}\sum_{n=0}^{\infty}\alpha^n P_n(\cos\psi), \qquad \alpha = \frac{a}{a'}.$$

At conjunction every Legendre polynomial equals one and the series is exactly geometric in $\alpha$. Its convergence is therefore completely transparent: each term is $\alpha$ times the last, and the number of terms needed for a given accuracy is set by nothing except how close the two orbits are.

Away from conjunction the Legendre polynomials are less than one and the convergence is faster, which is why conjunction is the case worth drawing. In the full theory these coefficients become the Laplace coefficients $b_s^{(j)}(\alpha)$, defined by an integral over one revolution and computed in practice by a recurrence that walks up in $j$. They are the only part of the whole construction that has to be evaluated numerically, and they are also the part that misbehaves: as $\alpha \to 1$ they diverge logarithmically, one after another, and the recurrence that generates them loses precision from the top down. A nineteenth-century theory’s practical accuracy limit was often set by how many Laplace coefficients could be tabulated by hand before the differencing became unreliable.

This is where the method has its boundary, and the boundary is sharp. As $\alpha \to 1$ the number of terms diverges, and if the orbits actually intersect there is a configuration in which $\Delta = 0$ and the disturbing function is infinite. No rearrangement helps. The expansion is not a bad approximation near an encounter; it is the wrong object.

That is the analytic statement of something the collection has already met numerically: a close approach has to be integrated. A step size that adapts to the encounter is the numerical version of the same admission, and an integrator that is right about the energy and wrong about the position is what one does instead. Perturbation theory and numerical integration are not competitors; they are the two halves of a subject divided along exactly this line.

Three kinds of term, and only one of them accumulates #

Sorting the terms by size is the first job. Sorting them by frequency is the second, and it is the one that decides what a planetary theory predicts.

Take one term, $C\cos(j\lambda + k\lambda' + \dots)$. Because $\lambda$ and $\lambda'$ advance at the two mean motions, the argument advances at $jn + kn'$, and the integral of the term over time carries a factor of $1/(jn + kn')$. Three cases follow, and they are qualitatively different:

• Short-period terms, where $j$ and $k$ are not both zero and the combination is not small. These oscillate at something near an orbital period, their integrals are divided by a large frequency, and they are correspondingly small. They make the ripple on an osculating element and nothing else.
• Secular terms, where $j = k = 0$. The argument does not contain the fast angles at all, so it does not oscillate on an orbital timescale; the integral has no small divisor because it has no divisor. These are what accumulate, and they are the whole of the theory that says no planet has an eccentricity of its own.
• Resonant terms, where $jn + kn'$ is nearly zero for small $j$ and $k$. The divisor is small, so a term that d’Alembert’s rule says is tiny gets multiplied by an enormous factor. This is the small-divisor problem, and it is the analytic origin of everything in resonance.

The practical consequence is the distinction between osculating and mean elements. An osculating semi-major axis is a real quantity, and it is the wrong one to publish: it jitters by parts in ten thousand over an orbit for reasons that have no long-term meaning. What a planetary theory publishes is the element with the short-period terms removed — a mean element — and the removal is done by the same expansion, term by term, because each short-period term’s contribution is known analytically once its frequency is known.

This is why the table that a modern ephemeris is does not consist of elements at all. Once the short-period structure is being modelled rather than averaged away, there is no advantage in elements over positions, and positions have the enormous merit of not requiring anybody to say which of the two conventions was used.

The distinction is not academic, and confusing the two is a live source of error. A satellite catalogue quotes mean elements, and the mean-element theory used to produce them has to be the same one used to propagate them: feeding an osculating state vector into a propagator expecting mean elements produces an orbit wrong by kilometres within a day, and wrong in a way that looks like a plausible perturbation rather than like a blunder. The published elements of a spacecraft and the published elements of an asteroid are, in general, not the same kind of object.

Where the subtraction is checked #

Nothing above is a measurement. The disturbing function is a rearrangement, and the honest question is where the rearrangement has been tested against something it could have failed.

The classical test is the one that made the method’s reputation. In the 1840s the residuals of Uranus’s position against a theory that included every known planet stood at about two arcminutes — small in absolute terms, enormous compared with the arcsecond precision of the observations, and, crucially, systematic in time rather than scattered. Le Verrier and Adams treated those residuals as a disturbing function with an unknown perturber in it and inverted for the perturber’s elements. Neptune was found within a degree of the predicted place. The subtraction had been carried out correctly enough that what was left over was a planet.

There is a second, less famous instance of the same inversion, and it went the other way. The same residuals that gave Neptune to Le Verrier gave him, twelve years later, an intra-Mercurial planet: the extra precession of Mercury’s perihelion was a leftover of exactly the same kind, treated with exactly the same method, and the answer was a body of about a tenth of Mercury’s mass inside its orbit. The method was not misapplied. What was wrong was the assumption that the only thing a residual can contain is more of the same force — and no amount of care with the subtraction could have revealed that, because the subtraction was correct.

The modern test is quieter and much sharper. Every planetary ephemeris is a fit whose model is a numerical integration, not a series — but the analysis of what that integration does is still done with the expansion. The secular frequencies of the solar system, the eight $g$ and seven $s$ eigenvalues that the linear theory produces, can be read out of a long numerical integration by Fourier analysis and compared against the values the second-order expansion predicts. They agree to a few parts in a thousand for the outer planets. Where they do not agree is more interesting: the inner planets’ frequencies drift, because $g_1 - g_5$ and $s_1 - s_2$ come close enough to each other for a small divisor to appear where the linear theory says there is none.

Where the picture stops #

Three limits stand out, and the third is the one that ended a two-hundred-year programme.

It assumes the orbits do not cross. Everything above is an expansion in $\alpha$, and a crossing orbit has no $\alpha$. Comets, planet-crossing asteroids and anything undergoing a close encounter are outside the theory and always were.

It assumes the divisors are not small. A resonance breaks the ordering the whole method rests on: the term the hierarchy says to discard is multiplied by a divisor the hierarchy knows nothing about. Resonant theories exist and are built by not expanding in the resonant angle — treating it as a slow variable in its own right, which is what turns a divergent series into a pendulum with a libration width.

And the series does not converge, in the end. This is Poincaré’s result of 1890 and it is not a technicality. The small divisors are dense: however far a system is from one resonance, there is another resonance with larger $j$ and $k$ arbitrarily nearby, and its divisor is arbitrarily small. So the coefficients of a general perturbation series do not stay bounded, and the series is asymptotic rather than convergent — it improves for a while as terms are added and then gets worse forever. A truncated planetary theory is accurate over a bounded span for a reason that has nothing to do with how many terms were kept.

That result is what makes the expiry date on a solar-system prediction a theorem rather than a limitation of computers, and it is why the question of what “no solution” means has an answer sharper than “it is hard”.

What the subtraction bought #

Be clear about what the method delivers, because the answer is not “the positions of the planets”. Numerical integration does that, and does it better.

What the expansion delivers is structure. It says which quantities are constant, which oscillate and which drift, and it says so with the constants of the system attached rather than as an observation about one particular integration. The secular eigenvalues are a property of the eight planets’ masses and spacings; the resonant arguments are a property of the ratios of the mean motions; d’Alembert’s rule is a property of analyticity. None of those can be read off a table of coordinates, and every one of them survives the failure of the series that produced it.

This is a pattern the collection meets repeatedly. An expansion that fails as a computation can still be right as an anatomy, and the anatomy is usually what was wanted. The forty-three arcseconds Mercury has left over were a residual against exactly this machinery — a number obtained by subtracting every term the theory could name, and finding the remainder neither zero nor random.

There is one more thing the subtraction buys, and it is the thing that makes the whole approach worth teaching after computers made it unnecessary for prediction. A term in the disturbing function is a physical statement in a way that a column of coordinates is not. The great inequality of Jupiter and Saturn — an oscillation of about twenty arcminutes in Jupiter’s longitude and fifty in Saturn’s, with a period of nine hundred years — is not a feature anybody would notice in an integration; it is a single resonant term with $j = 5$, $k = -2$, and a divisor of about a two-hundredth of the mean motion. Once it is named, the fact that it must exist follows from the ratio of two periods, its size follows from d’Alembert’s rule and the divisor, and its period follows from arithmetic. Laplace named it in 1785 and thereby closed a problem that had been open since Kepler noticed the two planets drifting from their tables in opposite directions.

Where the ladder goes next #

The next rung is the resonant case treated properly: what happens when the ordering fails, and the term the hierarchy discarded turns out to be the only one that matters. Further up sits the one this essay has quietly assumed — that the elements being perturbed are a good coordinate system in the first place, which stops being true when the eccentricity or the inclination goes to zero and the angles the expansion is written in cease to exist.