Central forces and the two-body problem

Two bodies that attract or repel each other along the line joining them form the most important problem in classical mechanics. It describes a planet and its star, a moon and its planet, a binary star, the electron and nucleus of the old Bohr atom, and an alpha particle passing a gold nucleus. Newton solved it in the Principia, and the method he and his successors developed, reducing two bodies to one and then using conservation laws to reduce three dimensions to one, is a model for how symmetry simplifies physics. This lesson develops the general theory for any central force; lesson 8 specializes to the inverse-square law and the orbits of Kepler, and lesson 9 to scattering.

Reducing two bodies to one

Bodies of mass \(m_1\) and \(m_2\) at \(\mathbf r_1\) and \(\mathbf r_2\) interact through a force that depends only on their separation \(\mathbf r = \mathbf r_1 - \mathbf r_2\) and points along it. With no external forces, lesson 2 tells us the centre of mass \(\mathbf R = (m_1\mathbf r_1 + m_2\mathbf r_2)/(m_1 + m_2)\) moves at constant velocity. Subtracting the two equations of motion, \(\ddot{\mathbf r}_1 = \mathbf F/m_1\) and \(\ddot{\mathbf r}_2 = -\mathbf F/m_2\), gives

The relative coordinate moves exactly like a single particle of reduced mass \(\mu\) attracted to a fixed centre. Once \(\mathbf r(t)\) is known, the individual positions are \(\mathbf r_1 = \mathbf R + (m_2/M)\mathbf r\) and \(\mathbf r_2 = \mathbf R - (m_1/M)\mathbf r\): both bodies trace similar orbits about the centre of mass, scaled in inverse proportion to their masses. When one mass is much larger, as for the Sun and a planet, \(\mu\) is close to the smaller mass and the large body barely moves. The Sun does wobble, by about its own radius under Jupiter's pull, and measuring such wobbles in other stars by the Doppler shift of their light was how the first planet around a Sun-like star, 51 Pegasi b, was found in 1995.

Conservation of angular momentum and Kepler's second law

A central force exerts no torque about the centre, since \(\mathbf r\times\mathbf F = 0\) when \(\mathbf F\) is parallel to \(\mathbf r\). The angular momentum \(\mathbf L = \mu\mathbf r\times\dot{\mathbf r}\) is therefore constant. Two consequences follow at once. First, since \(\mathbf r\) is always perpendicular to the fixed vector \(\mathbf L\), the motion lies in a fixed plane through the centre. We can use polar coordinates \((r, \phi)\) in that plane. Second, the magnitude

is constant. In a short time \(dt\) the line from the centre to the body sweeps a thin triangle of area \(dA = \tfrac12r^2\,d\phi\), so the areal velocity \(dA/dt = \ell/2\mu\) is constant. Equal areas are swept in equal times. This is Kepler's second law, and we have proved it for any central force, not only gravity. A planet moves fastest at perihelion, its closest approach, and slowest at aphelion. The Earth reaches perihelion in early January, which is why the northern winter half of the year, from the September equinox to the March equinox, is about a week shorter than the summer half.

The effective potential

With the motion confined to a plane and \(\dot\phi = \ell/\mu r^2\), the energy

involves only \(r\) and \(\dot r\). The radial motion is a one-dimensional problem in the effective potential \(U_{\text{eff}}\), and all the tools of lesson 4 apply. The extra term \(\ell^2/2\mu r^2\) is the centrifugal barrier: it is really the kinetic energy of the angular motion, which grows as the body comes closer to the centre because angular momentum is conserved. It pushes the body away from the centre and prevents it from falling in unless the attraction grows faster than \(1/r^2\) at short distance.

For gravity, \(U = -k/r\) with \(k = GMm\) (or \(Gm_1m_2\)). The effective potential has a minimum at \(r_0 = \ell^2/\mu k\), where a circular orbit is possible. At energies between that minimum and zero the radius oscillates between a perihelion \(r_{\min}\) and an aphelion \(r_{\max}\), the turning points where \(E = U_{\text{eff}}\). At \(E \ge 0\) there is only one turning point and the body escapes. Setting \(\dot r = 0\) gives a quadratic for the turning points, \(Er^2 + kr - \ell^2/2\mu = 0\). The program solves it for a sample orbit, then integrates the full two-dimensional motion to check.

import math

# Units: G M = 1, reduced mass mu = 1. Effective potential U_eff = l^2 / (2 r^2) - 1/r.
def turning_points(E, l):
    # roots of E r^2 + r - l^2 / 2 = 0
    disc = 1 + 2 * E * l * l
    return sorted(((-1 + s * math.sqrt(disc)) / (2 * E) for s in (1, -1)))

E, l = -0.3, 1.0
rmin, rmax = turning_points(E, l)
print(f"E = {E}, l = {l}: perihelion r = {rmin:.4f}, aphelion r = {rmax:.4f}")
a = (rmin + rmax) / 2; e = (rmax - rmin) / (rmax + rmin)
print(f"semi-major axis a = {a:.4f} (= -1/(2E) = {-1 / (2 * E):.4f}), eccentricity e = {e:.4f}")

# Integrate the orbit in Cartesian coordinates with velocity Verlet, adding a small 1/r^3 force
def orbit(eps, turns=6, dt=1e-4):
    x, y = rmin, 0.0
    vx, vy = 0.0, l / rmin
    def acc(x, y):
        r = math.hypot(x, y)
        f = -1 / r ** 2 - eps / r ** 3          # radial force per unit mass
        return f * x / r, f * y / r
    ax, ay = acc(x, y); peri, r_prev, r_now, t = [], None, rmin, 0.0
    E0, L0 = 0.5 * (vx * vx + vy * vy) - 1 / rmin - eps / (2 * rmin * rmin), x * vy - y * vx
    while len(peri) < turns:
        x += vx * dt + 0.5 * ax * dt * dt; y += vy * dt + 0.5 * ay * dt * dt
        nax, nay = acc(x, y)
        vx += 0.5 * (ax + nax) * dt; vy += 0.5 * (ay + nay) * dt; ax, ay = nax, nay
        r_new = math.hypot(x, y)
        if r_prev is not None and r_now < r_prev and r_now < r_new:   # passed a perihelion
            peri.append(math.degrees(math.atan2(y, x)))
        r_prev, r_now, t = r_now, r_new, t + dt
    E1 = 0.5 * (vx * vx + vy * vy) - 1 / r_now - eps / (2 * r_now * r_now)
    return peri, abs(E1 - E0) / abs(E0), abs((x * vy - y * vx) - L0)

for eps in (0.0, 0.02):
    peri, dE, dL = orbit(eps)
    steps = [(b - a) % 360 for a, b in zip(peri, peri[1:])]
    shift = sum(s if s < 180 else s - 360 for s in steps) / len(steps)
    shift = 0.0 if abs(shift) < 5e-4 else shift
    print(f"eps = {eps}: perihelion advance per orbit = {shift:+.3f} deg, relative energy error {dE:.1e}, |change in l| {dL:.1e}")
print(f"prediction for the 1/r^3 force: 360 (1/sqrt(1 - eps/l^2) - 1) = {360 * (1 / math.sqrt(1 - 0.02 / l ** 2) - 1):.3f} deg")

E = -0.3, l = 1.0: perihelion r = 0.6126, aphelion r = 2.7208
semi-major axis a = 1.6667 (= -1/(2E) = 1.6667), eccentricity e = 0.6325
eps = 0.0: perihelion advance per orbit = +0.000 deg, relative energy error 7.0e-14, |change in l| 7.5e-14
eps = 0.02: perihelion advance per orbit = +3.655 deg, relative energy error 2.2e-13, |change in l| 4.4e-14
prediction for the 1/r^3 force: 360 (1/sqrt(1 - eps/l^2) - 1) = 3.655 deg

In units with \(GM = \mu = 1\), the orbit with \(E = -0.3\) and \(\ell = 1\) swings between \(r = 0.613\) and \(r = 2.721\). The average of these, 1.667, equals \(-k/2E\), a result lesson 8 derives: for the inverse-square law the size of a bound orbit depends only on its energy. The numerical orbit returns to exactly the same perihelion after each revolution, and the integrator keeps energy and angular momentum constant to about \(10^{-13}\).

Stability of circular orbits

A circular orbit exists at any radius where \(U_{\text{eff}}'(r_0) = 0\), that is, where the attraction supplies exactly the centripetal force: \(F(r_0) = -\ell^2/\mu r_0^3\). It is stable if \(U_{\text{eff}}''(r_0) > 0\). For a power-law force \(F = -Cr^n\), the test reduces to \(n > -3\). Inverse-square gravity (\(n = -2\)) and a spring (\(n = 1\)) give stable circular orbits. A force falling faster than \(1/r^3\) does not: a slight push inward and the body spirals in, a slight push outward and it spirals away. Small oscillations about a stable circular orbit have radial frequency \(\omega_r = \sqrt{U_{\text{eff}}''(r_0)/\mu}\), and comparing it with the orbital frequency \(\omega_\phi = \dot\phi\) gives

For gravity the ratio is exactly 1: the radius completes one oscillation per orbit, so the perihelion returns to the same direction and the orbit closes. For the spring, the ratio is 2, giving a closed ellipse centred on the origin. For any other power the ratio is irrational in general, and the perihelion advances or regresses each orbit.

The Binet equation

To find the shape of the orbit rather than its time dependence, eliminate time. With \(u = 1/r\) and \(\dot\phi = \ell u^2/\mu\), the chain rule gives \(\dot r = -(\ell/\mu)\,du/d\phi\) and \(\ddot r = -(\ell^2u^2/\mu^2)\,d^2u/d\phi^2\). Substituting into the radial equation \(\mu(\ddot r - r\dot\phi^2) = F(r)\) yields the Binet equation, named for the French mathematician Jacques Binet:

For any central force, this is a second-order equation for the orbit's shape \(u(\phi)\). For the inverse-square force \(F = -ku^2\), the right side is the constant \(\mu k/\ell^2\), and the equation is that of a harmonic oscillator in \(\phi\) with a constant push. Its solution, \(u = (\mu k/\ell^2)(1 + e\cos\phi)\), is a conic section, as lesson 8 shows. The Binet equation also works backward: given an orbit, it tells you the force that produces it. Newton used essentially this reasoning to show that an elliptical orbit with the centre of force at a focus requires an inverse-square law.

Precessing orbits and Bertrand's theorem

Add a small inverse-cube term to gravity: \(F = -k/r^2 - \varepsilon/r^3\). In the Binet equation the new term is proportional to \(u\), and it can be moved to the left side: \(u'' + (1 - \mu\varepsilon/\ell^2)u = \mu k/\ell^2\). The orbit is still an ellipse in shape, but the radial oscillation now has angular period \(2\pi/\sqrt{1 - \mu\varepsilon/\ell^2}\), slightly more than \(2\pi\). The perihelion advances each orbit by

With \(\varepsilon = 0.02\) in the program's units the predicted advance is \(3.655^\circ\) per orbit, and the numerical integration gives \(3.655^\circ\), agreeing to every printed digit. The figure shows the resulting rosette for a larger perturbation.

This is not just a mathematical curiosity. General relativity modifies the effective potential for a planet by a term proportional to \(1/r^3\), producing a perihelion advance of Mercury's orbit of 43 seconds of arc per century beyond what the other planets' pulls explain. Urbain Le Verrier found the discrepancy in 1859, and Einstein's explanation in 1915 was the first triumph of his theory.

Joseph Bertrand proved in 1873 that only two central forces make every bounded orbit close on itself: the inverse-square attraction and the linear spring. For every other force law, almost all bound orbits precess and never exactly repeat, filling an annulus between the turning points over time. That the planets move on nearly closed ellipses is therefore strong evidence that gravity falls off as exactly the inverse square of distance.