Planetary Motion and Central Forces: From Conservation of Angular Momentum to Kepler's Three Laws

0. Key points

• The two-body problem separates completely into the motion of the centre of mass and the relative motion. The relative motion obeys the same equation as that of a single particle carrying the reduced mass $\mu = m_1 m_2/(m_1+m_2)$.
• From the sole assumption that the force is central (that it acts along the line joining the two bodies), the angular momentum $\boldsymbol{L}$ is conserved. As a consequence the motion is confined to a plane and the areal velocity is constant. This is Kepler’s second law, and it holds no matter how the magnitude of the force depends on distance.
• Using conservation of angular momentum, the two-dimensional motion reduces to one-dimensional motion in the effective potential $U_{\mathrm{eff}}(r) = U(r) + L^2/(2\mu r^2)$. The qualitative classification of orbits (bound, unbound, circular) can be read off from that single picture.
• The shape of the orbit is determined by Binet’s equation. Only for the inverse-square force $f(r) = -k/r^2$ does the equation become linear, and its solutions are the conic sections $r = \ell/(1 + e\cos\theta)$. This is Kepler’s first law. The eccentricity is expressed through the energy and the angular momentum as $e^2 = 1 + 2EL^2/(\mu k^2)$.
• Combining the constancy of the areal velocity with the area of an ellipse gives $T^2 = 4\pi^2 a^3 / \bigl(G(m_1+m_2)\bigr)$. This is Kepler’s third law, carrying the correction $(m_1+m_2)$ that was absent from Kepler’s own statement.
• Besides the angular momentum and the energy, the inverse-square force possesses one further conserved quantity, the Laplace–Runge–Lenz vector, and this is the reason the orbit closes (the perihelion does not move).

1. Motivation: where do the three empirical rules come from?

Johannes Kepler spent more than a decade analysing the observational record of Mars left by Tycho Brahe, and extracted three rules from it.

1. A planet traces an ellipse with the Sun at one focus.
2. The line segment joining the Sun to the planet sweeps out equal areas in equal times.
3. The square of the orbital period is proportional to the cube of the semi-major axis.

These are summaries of observation, not explanations. Why an ellipse? Why areas? Why a square and a cube? Kepler himself had no means of answering.

Newton’s answer was that the three rules are not independent. The second law follows from nothing more than the fact that the force points towards the Sun; it does not depend at all on how the strength of the force varies with distance. The first and third laws emerge once one adds the further condition that the force falls off as the inverse square of the distance. The three empirical rules therefore collapse, essentially, into the single assumption of an inverse-square central force.

In this article we carry that derivation through to the end. The only tools we need are the equation of motion treated in The foundations of Newtonian mechanics (the second law（Axiom 3.3）［Foundations of Newtonian Mechanics］) and analysis at roughly the level of Differentiation of functions of several variables and partial derivatives. The Lagrangian formalism of a later chapter (the same central-force motion is treated from the Lagrangian in polar coordinates in Example 5.3［Lagrangian Mechanics］) and Noether’s theorem will teach us that the conservation laws we obtain here by hand are in fact consequences of symmetries independent of any choice of coordinates. So that the article can be reread with that viewpoint in mind, we state each time which conservation law comes from which assumption.

Note

A word on notation. We write $m_1, m_2$ for the masses of the two bodies, $\boldsymbol{r}_1, \boldsymbol{r}_2$ for their position vectors, $\boldsymbol{r} = \boldsymbol{r}_1 - \boldsymbol{r}_2$ for the relative position, $r = |\boldsymbol{r}|$ for its magnitude and $\hat{\boldsymbol{r}} = \boldsymbol{r}/r$ for the corresponding unit vector. Time derivatives are denoted by dots.

2. Preliminaries: reducing the two-body problem to a one-body problem

2.1. Definition of a central force

Definition 2.1（Central force）

The force $\boldsymbol{F}$ exerted on particle 1 by particle 2 is called a central force if, for some real-valued function $f$, it can be written as

$$\boldsymbol{F} = f(r)\,\hat{\boldsymbol{r}},\qquad \boldsymbol{r} = \boldsymbol{r}_1 - \boldsymbol{r}_2,\quad r = |\boldsymbol{r}|,\quad \hat{\boldsymbol{r}} = \boldsymbol{r}/r .$$

That is, the force is directed along the line joining the two particles and its magnitude is determined by their separation alone. The force is attractive when $f(r) < 0$ and repulsive when $f(r) > 0$.

A central force is necessarily conservative. Indeed, setting $U(r) = -\int^r f(s)\,ds$ gives $f(r) = -U'(r)$, and the gradient of the spherically symmetric function $U(r)$ is

$$\nabla U(r) = U'(r)\,\nabla r = U'(r)\,\hat{\boldsymbol{r}},$$

so that $\boldsymbol{F} = -\nabla U$ holds (the identity $\nabla r = \hat{\boldsymbol{r}}$ is verified at once by differentiating $r = \sqrt{x^2+y^2+z^2}$ componentwise; for instance $\partial r/\partial x = x/r$). The mechanical energy is therefore conserved (Theorem 7.5［Foundations of Newtonian Mechanics］).

Newtonian gravitation has $f(r) = -Gm_1m_2/r^2$, the Coulomb force between point charges has $f(r) = q_1q_2/(4\pi\varepsilon_0 r^2)$, and the isotropic harmonic oscillator has $f(r) = -\mu\omega^2 r$; all of these are central forces. In what follows we write the inverse-square attraction as

$$f(r) = -\frac{k}{r^2},\qquad U(r) = -\frac{k}{r},\qquad k > 0 ,$$

so that $k = Gm_1m_2$ for gravitation.

2.2. Separating the centre-of-mass motion from the relative motion

Theorem 2.2（Separation of the two-body problem）

Let two particles of masses $m_1, m_2$ exert on each other a central force in the sense of Definition 2.1, with no external force present. Put $M = m_1+m_2$ for the total mass, $\boldsymbol{R} = (m_1\boldsymbol{r}_1+m_2\boldsymbol{r}_2)/M$ for the centre of mass, and

$$\mu = \frac{m_1 m_2}{m_1 + m_2}$$

for the reduced mass. Then the following hold.

1. $\ddot{\boldsymbol{R}} = \boldsymbol{0}$; that is, the centre of mass moves uniformly along a straight line.
2. The relative position $\boldsymbol{r} = \boldsymbol{r}_1 - \boldsymbol{r}_2$ satisfies $\mu\ddot{\boldsymbol{r}} = f(r)\hat{\boldsymbol{r}}$.
3. The total kinetic energy separates as $T = \tfrac12 M|\dot{\boldsymbol{R}}|^2 + \tfrac12\mu|\dot{\boldsymbol{r}}|^2$.

Proof（Theorem 2.2）

By the law of action and reaction (Axiom 3.4［Foundations of Newtonian Mechanics］), particle 1 feels the force $\boldsymbol{F} = f(r)\hat{\boldsymbol{r}}$ and particle 2 the force $-\boldsymbol{F}$. The equations of motion are

$$m_1\ddot{\boldsymbol{r}}_1 = \boldsymbol{F},\qquad m_2\ddot{\boldsymbol{r}}_2 = -\boldsymbol{F}.$$

(1) Adding the two equations gives $m_1\ddot{\boldsymbol{r}}_1 + m_2\ddot{\boldsymbol{r}}_2 = \boldsymbol{0}$. The left-hand side equals $M\ddot{\boldsymbol{R}}$, so $\ddot{\boldsymbol{R}} = \boldsymbol{0}$.

(2) Dividing the first equation by $m_1$, the second by $m_2$, and subtracting,

$$\ddot{\boldsymbol{r}} = \ddot{\boldsymbol{r}}_1 - \ddot{\boldsymbol{r}}_2 = \frac{\boldsymbol{F}}{m_1} + \frac{\boldsymbol{F}}{m_2} = \left(\frac{1}{m_1}+\frac{1}{m_2}\right)\boldsymbol{F} = \frac{\boldsymbol{F}}{\mu} .$$

The last equality is precisely the definition of the reduced mass, $1/\mu = 1/m_1 + 1/m_2$. Multiplying both sides by $\mu$ gives the assertion.

(3) From the definition of the centre of mass, $\boldsymbol{r}_1 = \boldsymbol{R} + (m_2/M)\boldsymbol{r}$ and $\boldsymbol{r}_2 = \boldsymbol{R} - (m_1/M)\boldsymbol{r}$. Substituting these,

$$\begin{aligned} T &= \tfrac12 m_1\left|\dot{\boldsymbol{R}} + \tfrac{m_2}{M}\dot{\boldsymbol{r}}\right|^2 + \tfrac12 m_2\left|\dot{\boldsymbol{R}} - \tfrac{m_1}{M}\dot{\boldsymbol{r}}\right|^2 \\ &= \tfrac12 (m_1+m_2)|\dot{\boldsymbol{R}}|^2 + \left(\tfrac{m_1m_2}{M} - \tfrac{m_2m_1}{M}\right)\langle \dot{\boldsymbol{R}}, \dot{\boldsymbol{r}}\rangle + \tfrac12\frac{m_1m_2^2 + m_2m_1^2}{M^2}|\dot{\boldsymbol{r}}|^2 . \end{aligned}$$

The coefficient of the cross term is $0$. The coefficient of the last term is $m_1m_2(m_2+m_1)/M^2 = m_1m_2/M = \mu$, so $T = \tfrac12 M|\dot{\boldsymbol{R}}|^2 + \tfrac12\mu|\dot{\boldsymbol{r}}|^2$.

∎

Thanks to this theorem it suffices, from now on, to consider only the one-body problem $\mu\ddot{\boldsymbol{r}} = f(r)\hat{\boldsymbol{r}}$. Working in the centre-of-mass frame (the inertial frame in which $\dot{\boldsymbol{R}} = \boldsymbol{0}$), the actual orbits of the two bodies are the relative orbit $\boldsymbol{r}(t)$ scaled by $m_2/M$ and by $-m_1/M$. In the solar system $m_2 \gg m_1$, so $\mu \approx m_1$ and the Sun may be regarded as essentially at rest. In a binary system, however, where the two masses are comparable, both stars trace similar ellipses about their common centre of mass.

3. Conservation of angular momentum and the areal velocity

Here the main argument begins. Let us first see what follows from the assumption that the force is central, and from that assumption alone.

Theorem 3.1（Conservation of angular momentum under a central force）

For motion obeying $\mu\ddot{\boldsymbol{r}} = f(r)\hat{\boldsymbol{r}}$, the angular momentum

$$\boldsymbol{L} = \boldsymbol{r}\times\boldsymbol{p} = \mu\,\boldsymbol{r}\times\dot{\boldsymbol{r}}$$

is independent of time. Here $f$ may be an arbitrary (continuous) function.

Proof（Theorem 3.1）

Apply the product rule to the cross product:

$$\frac{d\boldsymbol{L}}{dt} = \mu\frac{d}{dt}(\boldsymbol{r}\times\dot{\boldsymbol{r}}) = \mu\,(\dot{\boldsymbol{r}}\times\dot{\boldsymbol{r}}) + \mu\,(\boldsymbol{r}\times\ddot{\boldsymbol{r}}).$$

The first term is the cross product of a vector with itself and hence $\boldsymbol{0}$. Substituting the equation of motion $\mu\ddot{\boldsymbol{r}} = f(r)\hat{\boldsymbol{r}}$ from Theorem 2.2 (2) into the second term,

$$\boldsymbol{r}\times(\mu\ddot{\boldsymbol{r}}) = f(r)\,\boldsymbol{r}\times\hat{\boldsymbol{r}} = \frac{f(r)}{r}\,\boldsymbol{r}\times\boldsymbol{r} = \boldsymbol{0}.$$

All that was used here is that the force is parallel to $\boldsymbol{r}$, that is, the definition of a central force (Definition 2.1). Hence $d\boldsymbol{L}/dt = \boldsymbol{0}$.

∎

In one sentence, the reason for the conservation is that a central force exerts no torque $\boldsymbol{N} = \boldsymbol{r}\times\boldsymbol{F}$ about the origin. As we shall see in a later chapter, this is a consequence of the rotational symmetry of the relative coordinate system, and it is the most basic instance (Example 4.4［対称性と保存則］) of Symmetry and conservation laws (Noether’s theorem).

Corollary 3.2（Planarity of the motion）

If $\boldsymbol{L}\neq\boldsymbol{0}$, the entire motion takes place in a single plane through the origin perpendicular to $\boldsymbol{L}$. If $\boldsymbol{L}=\boldsymbol{0}$, the motion is confined to a single straight line through the origin.

Proof（Corollary 3.2）

Since $\boldsymbol{L} = \mu\,\boldsymbol{r}\times\dot{\boldsymbol{r}}$ is a cross product, $\langle \boldsymbol{L}, \boldsymbol{r}\rangle = 0$ holds at all times. By Theorem 3.1 the vector $\boldsymbol{L}$ is constant, so when $\boldsymbol{L}\neq\boldsymbol{0}$ the relation $\langle \boldsymbol{L}, \boldsymbol{r}(t)\rangle = 0$ says that $\boldsymbol{r}(t)$ lies in the plane through the origin with normal $\boldsymbol{L}$. This holds for every $t$, so the motion is planar.

If $\boldsymbol{L}=\boldsymbol{0}$ then $\boldsymbol{r}\times\dot{\boldsymbol{r}} = \boldsymbol{0}$, that is, $\dot{\boldsymbol{r}}$ is always parallel to $\boldsymbol{r}$. Computing the derivative of $\hat{\boldsymbol{r}}$ on an interval where $\boldsymbol{r}\neq\boldsymbol{0}$,

$$\frac{d\hat{\boldsymbol{r}}}{dt} = \frac{\dot{\boldsymbol{r}}}{r} - \frac{\dot r}{r^2}\boldsymbol{r} = \frac{1}{r}\left(\dot{\boldsymbol{r}} - \dot r\,\hat{\boldsymbol{r}}\right),$$

and writing $\dot{\boldsymbol{r}} = \lambda\boldsymbol{r}$ gives $\dot r = \lambda r$, whence $\dot{\boldsymbol{r}} - \dot r\hat{\boldsymbol{r}} = \lambda\boldsymbol{r} - \lambda r\hat{\boldsymbol{r}} = \boldsymbol{0}$. Thus the direction $\hat{\boldsymbol{r}}$ is constant and the motion is one-dimensional, along a line through the origin.

∎

From now on we assume $\boldsymbol{L}\neq\boldsymbol{0}$ and introduce polar coordinates $(r,\theta)$ in the plane of motion. The position is $\boldsymbol{r} = r\,\hat{\boldsymbol{e}}_r$, and the unit vectors satisfy

$$\frac{d\hat{\boldsymbol{e}}_r}{dt} = \dot\theta\,\hat{\boldsymbol{e}}_\theta,\qquad \frac{d\hat{\boldsymbol{e}}_\theta}{dt} = -\dot\theta\,\hat{\boldsymbol{e}}_r$$

(these follow immediately upon differentiating $\hat{\boldsymbol{e}}_r = (\cos\theta,\sin\theta)$ and $\hat{\boldsymbol{e}}_\theta = (-\sin\theta,\cos\theta)$ with respect to $t$). Writing out the velocity and the acceleration with their help,

$$\dot{\boldsymbol{r}} = \dot r\,\hat{\boldsymbol{e}}_r + r\dot\theta\,\hat{\boldsymbol{e}}_\theta,\qquad \ddot{\boldsymbol{r}} = (\ddot r - r\dot\theta^2)\,\hat{\boldsymbol{e}}_r + (r\ddot\theta + 2\dot r\dot\theta)\,\hat{\boldsymbol{e}}_\theta .$$

The magnitude of the angular momentum is $L = |\boldsymbol{L}| = \mu|r\hat{\boldsymbol{e}}_r \times (\dot r\hat{\boldsymbol{e}}_r + r\dot\theta\hat{\boldsymbol{e}}_\theta)| = \mu r^2|\dot\theta|$. Choosing the sense of $\theta$ so that $\dot\theta > 0$, we write

$$L = \mu r^2\dot\theta .$$

This is the relation we shall use again and again below.

Theorem 3.3（Kepler's second law (constancy of the areal velocity)）

For planar motion under a central force, let $A(t)$ be the area swept out by the segment joining the origin to the particle between the times $t_0$ and $t$. Then

$$\frac{dA}{dt} = \frac{1}{2}r^2\dot\theta = \frac{L}{2\mu} = \text{constant}.$$

In particular, equal areas are swept out in equal times. Here again $f$ may be an arbitrary function; the inverse-square law is not used.

Proof（Theorem 3.3）

Suppose that during an infinitesimal time $dt$ the radius vector turns from $\theta$ to $\theta + d\theta$ while its length changes from $r$ to $r + dr$. The region swept out is a thin sector with sides $r$ and $r+dr$ and central angle $d\theta$, whose area is the polar area element $dA = \tfrac12 r^2\,d\theta$ plus terms of higher order $O(dr\,d\theta)$. Rigorously, the area swept out is given by the double integral

$$A = \int_{\theta_0}^{\theta_1}\!\!\int_0^{r(\theta)} \rho \,d\rho\,d\theta = \int_{\theta_0}^{\theta_1}\frac{r(\theta)^2}{2}\,d\theta$$

(we used that the Jacobian of polar coordinates is $\rho$; see Example 6.6［重積分と累次積分］ in Multiple integrals and iterated integrals). Differentiating both sides with respect to $t$ and using the chain rule gives $dA/dt = \tfrac12 r(\theta)^2\,\dot\theta$.

Substituting $L = \mu r^2\dot\theta$, obtained from Theorem 3.1, so that $r^2\dot\theta = L/\mu$, we get

$$\frac{dA}{dt} = \frac{L}{2\mu},$$

and since $L$ and $\mu$ are constants the areal velocity is constant.

∎

Note that the second law, which Kepler extracted from observations of Mars, has come out without any use whatever of how the force depends on distance. The second law is no more than a geometric restatement of the fact that the Sun attracts, and it is no evidence for the inverse-square law. If the second law appeared to be violated, that would signal that the force is not central — that there are perturbations from bodies other than the Sun, or relativistic effects.

flowchart TD
A["force is central: F = f(r) r̂"] --> B["torque r × F = 0"]
B --> C["angular momentum L conserved"]
C --> D["motion is planar (L ≠ 0)"]
C --> E["areal velocity L/2μ constant = Kepler's 2nd law"]
A --> F["F is conservative: U(r) exists"]
F --> G["mechanical energy E conserved"]
C --> H["1-D motion in effective potential U + L²/2μr²"]
G --> H
H --> I["add f(r) = -k/r²"]
I --> J["conic-section orbits = Kepler's 1st and 3rd laws"]

The conserved quantities of central-force motion and the assumptions each one rests on

4. The effective potential and the radial motion

Having secured two conservation laws, we now reduce the number of degrees of freedom.

Definition 4.1（Effective potential）

With the magnitude $L$ of the angular momentum held fixed,

$$U_{\mathrm{eff}}(r) = U(r) + \frac{L^2}{2\mu r^2}$$

is called the effective potential. The second term $L^2/(2\mu r^2)$ is called the centrifugal barrier.

Proposition 4.2（Reduction to a one-dimensional radial problem）

For planar motion under a central force, the mechanical energy

$$E = \frac{1}{2}\mu|\dot{\boldsymbol{r}}|^2 + U(r)$$

is conserved, and moreover can be written as

$$E = \frac{1}{2}\mu\dot r^2 + U_{\mathrm{eff}}(r).$$

That is, $r(t)$ obeys the same equation as the one-dimensional motion of a particle of mass $\mu$ moving in the potential $U_{\mathrm{eff}}$.

Proof（Proposition 4.2）

Let us first verify conservation of energy. Taking the inner product of both sides of $\mu\ddot{\boldsymbol{r}} = -\nabla U$ with $\dot{\boldsymbol{r}}$,

$$\mu\langle\ddot{\boldsymbol{r}},\dot{\boldsymbol{r}}\rangle = -\langle\nabla U,\dot{\boldsymbol{r}}\rangle \;\Longleftrightarrow\; \frac{d}{dt}\left(\frac{1}{2}\mu|\dot{\boldsymbol{r}}|^2\right) = -\frac{d}{dt}U(r(t))$$

(the right-hand side by the chain rule), so the time derivative of $E$ vanishes. What was used here is the fact established in §2.1 that a central force is conservative.

Next we use the decomposition of the velocity, $\dot{\boldsymbol{r}} = \dot r\,\hat{\boldsymbol{e}}_r + r\dot\theta\,\hat{\boldsymbol{e}}_\theta$. Since $\hat{\boldsymbol{e}}_r$ and $\hat{\boldsymbol{e}}_\theta$ are mutually orthogonal unit vectors,

$$|\dot{\boldsymbol{r}}|^2 = \dot r^2 + r^2\dot\theta^2 .$$

Substituting the relation $\dot\theta = L/(\mu r^2)$ from Theorem 3.1 gives $r^2\dot\theta^2 = r^2\cdot L^2/(\mu^2r^4) = L^2/(\mu^2r^2)$. Hence

$$E = \frac{1}{2}\mu\dot r^2 + \frac{L^2}{2\mu r^2} + U(r) = \frac{1}{2}\mu\dot r^2 + U_{\mathrm{eff}}(r),$$

using the definition Definition 4.1.

∎

This reduction is powerful. Since $\dot r^2 \ge 0$, the motion is confined to the range of $r$ satisfying $U_{\mathrm{eff}}(r)\le E$, and at a point where equality holds (a turning point) we have $\dot r = 0$. For the inverse-square attraction $U = -k/r$,

$$U_{\mathrm{eff}}(r) = -\frac{k}{r} + \frac{L^2}{2\mu r^2}$$

tends to $+\infty$ as $r\to 0^+$ (the centrifugal barrier wins) and to $0^-$ as $r\to\infty$, with a single minimum in between. The position of the minimum follows from $U_{\mathrm{eff}}'(r) = k/r^2 - L^2/(\mu r^3) = 0$, namely

$$r_c = \frac{L^2}{\mu k},\qquad U_{\mathrm{eff}}(r_c) = -\frac{\mu k^2}{2L^2} .$$

The effective potential of the inverse-square attraction. The height of the energy fixes the type of the orbit

Let us collect what the figure tells us. When $E = U_{\mathrm{eff}}(r_c) = -\mu k^2/(2L^2)$, the coordinate $r$ cannot move at all and the orbit is circular. When $U_{\mathrm{eff}}(r_c) < E < 0$, the coordinate $r$ oscillates between two turning points $r_{\min}, r_{\max}$ and the orbit stays within a finite distance of the origin (bound motion). When $E\ge 0$ there is only one turning point, and the particle departs to infinity after its closest approach (unbound motion). In §5 we confirm that this classification corresponds exactly to ellipse, parabola and hyperbola.

Example 4.3（Stability of circular orbits for power-law central forces）

Consider the attraction $f(r) = -k/r^n$ (with $k>0$ and $n\neq 1$). The corresponding potential is $U(r) = -k/\bigl((n-1)r^{n-1}\bigr)$, and the effective potential is

$$U_{\mathrm{eff}}(r) = -\frac{k}{(n-1)r^{n-1}} + \frac{L^2}{2\mu r^2} .$$

Circular orbits correspond to stationary points of $U_{\mathrm{eff}}$. From

$$U_{\mathrm{eff}}'(r) = \frac{k}{r^n} - \frac{L^2}{\mu r^3} = 0 \iff r_c^{\,3-n} = \frac{L^2}{\mu k}$$

we see that for $n\neq 3$ each value of $L$ determines exactly one circular radius $r_c$. Stability is decided by the sign of the second derivative,

$$U_{\mathrm{eff}}''(r) = -\frac{nk}{r^{n+1}} + \frac{3L^2}{\mu r^4} .$$

Substituting the stationarity condition $L^2/\mu = k\,r_c^{\,3-n}$ into the second term gives $3L^2/(\mu r_c^4) = 3k\,r_c^{\,-1-n}$, so that

$$U_{\mathrm{eff}}''(r_c) = \frac{(3-n)k}{r_c^{\,n+1}} .$$

Since $k>0$ and $r_c>0$, this is positive only when $n < 3$. In other words, for attractions steeper than the inverse cube the circular orbit is unstable, and the slightest perturbation makes the particle fall into the centre or fly off to infinity. The inverse-square force has $n=2$, so $U_{\mathrm{eff}}''(r_c) = k/r_c^3 > 0$ and the circular orbit is stable. This stability is one of the reasons a planetary system can persist.

5. The shape of the orbit: Binet’s equation and Kepler’s first law

So far we have followed the time evolution of $r(t)$. If what we want is the shape of the orbit, however, it is quicker to eliminate the time and find $r$ as a function of $\theta$. The key is the substitution $u = 1/r$.

Lemma 5.1（Binet's orbit equation）

Let $\boldsymbol{L}\neq\boldsymbol{0}$ and put $u(\theta) = 1/r(\theta)$. Then the orbit under the central force $f(r)$ satisfies

$$\frac{d^2u}{d\theta^2} + u = -\frac{\mu}{L^2u^2}\,f\!\left(\frac{1}{u}\right).$$

Proof（Lemma 5.1）

From $L = \mu r^2\dot\theta$ we have $\dot\theta = Lu^2/\mu$. Since $L\neq 0$, the sign of $\dot\theta$ never changes and we may take $\theta$ as the independent variable.

First rewrite $\dot r$ in terms of derivatives with respect to $\theta$. Since $r = 1/u$, the chain rule gives

$$\dot r = \frac{d}{dt}\left(\frac{1}{u}\right) = -\frac{1}{u^2}\frac{du}{d\theta}\dot\theta = -\frac{1}{u^2}\frac{du}{d\theta}\cdot\frac{Lu^2}{\mu} = -\frac{L}{\mu}\frac{du}{d\theta}.$$

The point of the substitution is precisely that the factors $u^2$ cancel cleanly. Differentiating once more,

$$\ddot r = -\frac{L}{\mu}\frac{d^2u}{d\theta^2}\dot\theta = -\frac{L}{\mu}\frac{d^2u}{d\theta^2}\cdot\frac{Lu^2}{\mu} = -\frac{L^2u^2}{\mu^2}\frac{d^2u}{d\theta^2}.$$

On the other hand, from the $\hat{\boldsymbol{e}}_r$ component of the acceleration obtained in §3, the radial component of the equation of motion is

$$\mu(\ddot r - r\dot\theta^2) = f(r).$$

Here

$$r\dot\theta^2 = \frac{1}{u}\left(\frac{Lu^2}{\mu}\right)^2 = \frac{L^2u^3}{\mu^2},$$

so substituting gives

$$\mu\left(-\frac{L^2u^2}{\mu^2}\frac{d^2u}{d\theta^2} - \frac{L^2u^3}{\mu^2}\right) = f(1/u) \;\Longleftrightarrow\; -\frac{L^2u^2}{\mu}\left(\frac{d^2u}{d\theta^2} + u\right) = f(1/u).$$

Since $u\neq 0$ ($r$ being finite), we may divide both sides by $-L^2u^2/\mu$ to obtain the stated equation.

∎

The right-hand side of Lemma 5.1 is in general a nonlinear function of $u$. But precisely when $f(r) = -k/r^2$, that is $f(1/u) = -ku^2$, the factor $u^2$ cancels and the right-hand side becomes a constant. This single point is what makes the inverse-square law special.

Theorem 5.2（Kepler's first law (orbits are conic sections)）

Under the inverse-square attraction $f(r) = -k/r^2$ (with $k>0$), the orbit of a motion with $\boldsymbol{L}\neq\boldsymbol{0}$ is given by

$$r(\theta) = \frac{\ell}{1 + e\cos(\theta-\theta_0)},\qquad \ell = \frac{L^2}{\mu k},\quad e \ge 0 .$$

This is a conic section with the centre of force (the origin) at one focus: an ellipse if $e<1$, a parabola if $e=1$, and one branch of a hyperbola if $e>1$. In particular, for bound motion ($e<1$) the orbit is an ellipse with the centre of force at one of its foci.

Proof（Theorem 5.2）

Substitute $f(1/u) = -ku^2$ into Lemma 5.1:

$$\frac{d^2u}{d\theta^2} + u = -\frac{\mu}{L^2u^2}\cdot(-ku^2) = \frac{\mu k}{L^2} = \frac{1}{\ell}.$$

This is a second-order linear inhomogeneous ordinary differential equation with constant coefficients. A particular solution is the constant function $u_p = 1/\ell$, and the general solution of the homogeneous equation $u'' + u = 0$ is $C\cos(\theta-\theta_0)$ (with $C\ge 0$ and $\theta_0$ constant; this is Theorem 5.2［Foundations of Newtonian Mechanics］ with the angular frequency set to $1$ and the time variable read as $\theta$). Hence the general solution is

$$u(\theta) = \frac{1}{\ell} + C\cos(\theta-\theta_0).$$

Setting $e = C\ell \ge 0$ and taking reciprocals,

$$r(\theta) = \frac{1}{u(\theta)} = \frac{\ell}{1 + e\cos(\theta-\theta_0)} .$$

From now on we fix the reference direction of $\theta$ so that $\theta_0 = 0$.

Let us check that this is a conic section. Writing $r + er\cos\theta = \ell$ in Cartesian coordinates $x = r\cos\theta$, $y = r\sin\theta$ gives $\sqrt{x^2+y^2} = \ell - ex$, and squaring both sides,

$$x^2 + y^2 = \ell^2 - 2e\ell x + e^2x^2 \iff (1-e^2)x^2 + 2e\ell x + y^2 = \ell^2 .$$

When $e<1$, completing the square in $x$ gives

$$(1-e^2)\left(x + \frac{e\ell}{1-e^2}\right)^2 + y^2 = \ell^2 + \frac{e^2\ell^2}{1-e^2} = \frac{\ell^2}{1-e^2},$$

and dividing both sides by the right-hand side yields the standard form of an ellipse,

$$\frac{\bigl(x + a e\bigr)^2}{a^2} + \frac{y^2}{b^2} = 1,\qquad a = \frac{\ell}{1-e^2},\quad b = \frac{\ell}{\sqrt{1-e^2}} = a\sqrt{1-e^2}.$$

Its centre is at $(-ae, 0)$, at distance $ae$ from the origin. Since $a^2 - b^2 = a^2e^2$, this distance is exactly the focal distance, so the origin is one of the foci of the ellipse. When $e=1$ the $x^2$ term disappears and we get the parabola $y^2 = \ell^2 - 2\ell x$; when $e>1$ we have $1-e^2<0$ and the same manipulation produces the standard form of a hyperbola.

∎

Definition 5.3（Orbital elements）

The quantities appearing in Theorem 5.2 are named as follows. We call $e$ the eccentricity, $\ell = L^2/(\mu k)$ the semi-latus rectum, and, for an elliptical orbit, $a = \ell/(1-e^2)$ the semi-major axis and $b = a\sqrt{1-e^2}$ the semi-minor axis. The point at $\theta=0$, where $r$ is smallest, is the periapsis (the perihelion for motion about the Sun), and the point at $\theta=\pi$, where $r$ is largest, is the apoapsis; thus

$$r_{\min} = \frac{\ell}{1+e} = a(1-e),\qquad r_{\max} = \frac{\ell}{1-e} = a(1+e).$$

In particular $\ell = a(1-e^2) = b^2/a$.

Geometry of an elliptical orbit. The Sun sits at a focus, not at the centre of the ellipse

Corollary 5.4（Relation between eccentricity and energy）

For the orbit of Theorem 5.2, the mechanical energy $E$ and the magnitude $L$ of the angular momentum are related by

$$e^2 = 1 + \frac{2EL^2}{\mu k^2}.$$

Hence $E<0 \iff e<1$ (ellipse), $E=0\iff e=1$ (parabola) and $E>0\iff e>1$ (hyperbola). Moreover, for an elliptical orbit,

$$a = -\frac{k}{2E},$$

so that the semi-major axis is determined by the energy alone.

Proof（Corollary 5.4）

Write the general solution from the proof of Theorem 5.2 as $u = 1/\ell + C\cos\theta$ (with $C = e/\ell$). From $\dot r = -(L/\mu)\,du/d\theta$, obtained in the proof of Lemma 5.1,

$$\dot r = -\frac{L}{\mu}\cdot(-C\sin\theta) = \frac{LC}{\mu}\sin\theta .$$

Insert this into the energy expression of Proposition 4.2. Noting that $U = -k/r = -ku$,

$$\begin{aligned} E &= \frac{1}{2}\mu\dot r^2 + \frac{L^2u^2}{2\mu} - ku \\ &= \frac{L^2C^2}{2\mu}\sin^2\theta + \frac{L^2}{2\mu}\left(\frac{1}{\ell}+C\cos\theta\right)^2 - k\left(\frac{1}{\ell}+C\cos\theta\right) \\ &= \frac{L^2C^2}{2\mu}\sin^2\theta + \frac{L^2}{2\mu}\left(\frac{1}{\ell^2} + \frac{2C\cos\theta}{\ell} + C^2\cos^2\theta\right) - \frac{k}{\ell} - kC\cos\theta . \end{aligned}$$

Using $\ell = L^2/(\mu k)$, that is $L^2/(\mu\ell) = k$, the terms proportional to $\cos\theta$ cancel: $\dfrac{L^2}{\mu}\dfrac{C\cos\theta}{\ell} - kC\cos\theta = kC\cos\theta - kC\cos\theta = 0$. Also, $\sin^2\theta + \cos^2\theta = 1$ lets the $C^2$ terms combine, while the constant terms give $\dfrac{L^2}{2\mu\ell^2} - \dfrac{k}{\ell} = \dfrac{k}{2\ell} - \dfrac{k}{\ell} = -\dfrac{k}{2\ell} = -\dfrac{\mu k^2}{2L^2}$. Altogether

$$E = \frac{L^2C^2}{2\mu} - \frac{\mu k^2}{2L^2}.$$

Substituting $C = e/\ell = e\mu k/L^2$ gives $L^2C^2/(2\mu) = \mu k^2e^2/(2L^2)$, so

$$E = \frac{\mu k^2}{2L^2}\left(e^2 - 1\right) \;\Longleftrightarrow\; e^2 = 1 + \frac{2EL^2}{\mu k^2} .$$

Since $\mu, k, L^2$ are all positive, the signs of $E$ and of $e^2-1$ agree, and the classification of types follows.

In the elliptical case, Definition 5.3 gives $a = \ell/(1-e^2)$, and since $1-e^2 = -2EL^2/(\mu k^2)$,

$$a = \frac{L^2}{\mu k}\cdot\frac{\mu k^2}{-2EL^2} = -\frac{k}{2E}$$

(with $a>0$ because $E<0$).

∎

This corollary is extremely convenient in practice. The angular momentum fixes only the “slenderness” of the orbit (the eccentricity), and the energy only its “size” (the semi-major axis). A circular orbit is the case $e=0$, that is $E = -\mu k^2/(2L^2)$, which coincides with the minimum of $U_{\mathrm{eff}}$ found in §4. Two independent derivations have agreed.

Example 5.5（The hyperbolic orbit of an interstellar object）

The object 1I/ʻOumuamua, discovered in 2017, passed through the solar system on an orbit with eccentricity $e \approx 1.20$ and perihelion distance $q \approx 0.255\ \mathrm{au}$. Since $e>1$, Corollary 5.4 gives $E>0$; the object is not bound to the Sun.

We have $r\to\infty$ when the denominator vanishes, that is when $\cos\theta_\infty = -1/e = -0.833$, so $\theta_\infty = 146.4^\circ$. The angle between the incoming and outgoing directions (the deflection of the orbit) is $2\theta_\infty - 180^\circ = 112.8^\circ$: the Sun’s gravity bent the direction of travel by about $113$ degrees.

Let us find the speed $v_\infty$ at infinity. From $q = a'(e-1)$ (where $a' = |a| = k/(2E)$ is the real semi-axis of the hyperbola) we get $a' = 0.255/0.20 = 1.275\ \mathrm{au} = 1.907\times10^{11}\ \mathrm{m}$. From $E = \tfrac12\mu v_\infty^2$ and $E = k/(2a')$ (the sign-reversed version of Corollary 5.4 valid for $e>1$), with $\mu\approx m$ (the mass of the object being negligible compared with the Sun’s),

$$v_\infty = \sqrt{\frac{GM_\odot}{a'}} = \sqrt{\frac{1.327\times10^{20}}{1.907\times10^{11}}} = \sqrt{6.958\times10^{8}} \approx 2.64\times10^{4}\ \mathrm{m/s},$$

that is about $26\ \mathrm{km/s}$. This agrees well with the observed value and was the ground for concluding that the object came from outside the solar system. The speed at perihelion follows from $v_q = \sqrt{GM_\odot(1+e)/q}$:

$$v_q = \sqrt{\frac{1.327\times10^{20}\times 2.20}{3.815\times10^{10}}} = \sqrt{7.65\times10^{9}} \approx 8.75\times10^{4}\ \mathrm{m/s},$$

reaching about $87\ \mathrm{km/s}$.

6. Kepler’s third law

Theorem 6.1（Kepler's third law (the harmonic law)）

Under the inverse-square attraction $f(r) = -k/r^2$, the orbital period $T$ and the semi-major axis $a$ of an elliptical orbit ($e<1$) satisfy

$$T^2 = \frac{4\pi^2\mu}{k}\,a^3 .$$

In particular, for Newtonian gravitation with $k = Gm_1m_2$ and $\mu = m_1m_2/(m_1+m_2)$,

$$T^2 = \frac{4\pi^2}{G(m_1+m_2)}\,a^3 .$$

The constant of proportionality does not depend on the eccentricity; it is fixed by the sum of the two masses alone.

Proof（Theorem 6.1）

By Theorem 3.3 the areal velocity is the constant $L/(2\mu)$. The area swept out in one revolution is the whole area $\pi ab$ of the ellipse, so

$$T = \frac{\pi ab}{L/(2\mu)} = \frac{2\pi\mu\,ab}{L}.$$

Now express $L$ through the orbital elements. Equating $\ell = L^2/(\mu k)$ and $\ell = b^2/a$ from Definition 5.3,

$$L^2 = \mu k \ell = \frac{\mu k b^2}{a} \;\Longrightarrow\; L = b\sqrt{\frac{\mu k}{a}}$$

(taking $L>0$). Substituting this into the expression for $T$ above, the factor $b$ cancels:

$$T = \frac{2\pi\mu\,ab}{b}\sqrt{\frac{a}{\mu k}} = 2\pi\mu a\sqrt{\frac{a}{\mu k}} = 2\pi\sqrt{\frac{\mu a^3}{k}} .$$

Squaring both sides gives $T^2 = 4\pi^2\mu a^3/k$.

For Newtonian gravitation, $\mu/k = \dfrac{m_1m_2/(m_1+m_2)}{Gm_1m_2} = \dfrac{1}{G(m_1+m_2)}$, and substituting this yields the second formula.

∎

Remark 6.2

Kepler’s own third law asserted that ”$T^2/a^3$ is the same for every planet”. Theorem 6.1 corrects this. Exactly, $T^2/a^3 = 4\pi^2/\bigl(G(M_\odot + m)\bigr)$, so the value differs from planet to planet by the amount of the planetary mass $m$. Even for Jupiter, the heaviest in the solar system, $m/M_\odot \approx 9.5\times10^{-4}$, so the ratio changes only by about $0.1\%$ — undetectable at Kepler’s observational precision. In binary systems, on the other hand, where the two stellar masses are comparable, this correction term is precisely what makes it possible to measure stellar masses: measuring the period and the semi-major axis gives $m_1+m_2$ directly.

Example 6.3（Computing the orbital period of the Earth）

The Sun’s gravitational parameter is $GM_\odot = 1.32712\times10^{20}\ \mathrm{m^3/s^2}$ and the Earth’s semi-major axis is $a = 1.49598\times10^{11}\ \mathrm{m}$. The Earth’s mass is $3\times10^{-6}$ times the Sun’s, so we neglect the correction of Remark 6.2 and set $G(m_1+m_2)\approx GM_\odot$. By Theorem 6.1,

$$T = 2\pi\sqrt{\frac{a^3}{GM_\odot}} .$$

Working through the arithmetic,

$$a^3 = (1.49598\times10^{11})^3 = 3.3479\times10^{33}\ \mathrm{m^3},$$$$\frac{a^3}{GM_\odot} = \frac{3.3479\times10^{33}}{1.32712\times10^{20}} = 2.5227\times10^{13}\ \mathrm{s^2},$$$$\sqrt{2.5227\times10^{13}} = 5.0226\times10^{6}\ \mathrm{s},\qquad T = 2\pi\times5.0226\times10^{6} = 3.1559\times10^{7}\ \mathrm{s}.$$

Converting to days, $3.1559\times10^{7}/86400 = 365.3$ days, which agrees with the actual sidereal year of $365.256$ days to four significant figures. Note that the observed value is reproduced to this accuracy from the single assumption of the inverse-square law.

Example 6.4（The isotropic harmonic oscillator: the other force with closed orbits）

Consider the central force $f(r) = -\mu\omega^2 r$ (with potential $U = \tfrac12\mu\omega^2r^2$). In Cartesian coordinates the equations of motion separate completely into $\ddot x = -\omega^2 x$ and $\ddot y = -\omega^2 y$, so the general solution is

$$x(t) = a\cos\omega t,\qquad y(t) = b\sin\omega t$$

(with the phases adjusted by the initial conditions), and the orbit is $(x/a)^2 + (y/b)^2 = 1$, that is, an ellipse whose centre is the centre of force. Contrast this with the inverse-square case, where the centre of force is a focus.

Let us check the angular momentum:

$$L = \mu(x\dot y - y\dot x) = \mu\bigl(a\cos\omega t\cdot b\omega\cos\omega t + b\sin\omega t\cdot a\omega\sin\omega t\bigr) = \mu ab\omega ,$$

which is indeed constant (consistent with Theorem 3.1). The period is $T = 2\pi/\omega$, entirely independent of the amplitude — a different dependence from the $T\propto a^{3/2}$ of the Kepler problem.

It is known that the only central forces for which every bounded orbit closes are the inverse-square force $-k/r^2$ and the harmonic force $-\mu\omega^2 r$ (Bertrand’s theorem). This fact is two sides of one coin with the existence of an extra conserved quantity for exactly these two forces; we examine one of them in the next section.

7. A hidden symmetry: the Laplace–Runge–Lenz vector

Motion in three dimensions has three degrees of freedom, so phase space is six-dimensional. The energy $E$ and the angular momentum $\boldsymbol{L}$ (three components) give four conserved quantities, but the inverse-square force possesses one more independent conserved quantity.

Proposition 7.1（Conservation of the Laplace–Runge–Lenz vector）

Under the inverse-square attraction $\mu\ddot{\boldsymbol{r}} = -\dfrac{k}{r^2}\hat{\boldsymbol{r}}$, setting $\boldsymbol{p} = \mu\dot{\boldsymbol{r}}$, the vector

$$\boldsymbol{A} = \boldsymbol{p}\times\boldsymbol{L} - \mu k\,\hat{\boldsymbol{r}}$$

is conserved. Moreover $\boldsymbol{A}$ lies in the plane of motion, points towards the periapsis, and has magnitude $|\boldsymbol{A}| = \mu k e$.

Proof（Proposition 7.1）

By Theorem 3.1 the vector $\boldsymbol{L}$ is constant, so

$$\frac{d}{dt}(\boldsymbol{p}\times\boldsymbol{L}) = \dot{\boldsymbol{p}}\times\boldsymbol{L} = \left(-\frac{k}{r^3}\boldsymbol{r}\right)\times\left(\mu\,\boldsymbol{r}\times\dot{\boldsymbol{r}}\right)$$

(using $\hat{\boldsymbol{r}} = \boldsymbol{r}/r$). Applying the vector triple product identity $\boldsymbol{a}\times(\boldsymbol{b}\times\boldsymbol{c}) = \boldsymbol{b}\langle\boldsymbol{a},\boldsymbol{c}\rangle - \boldsymbol{c}\langle\boldsymbol{a},\boldsymbol{b}\rangle$,

$$\boldsymbol{r}\times(\boldsymbol{r}\times\dot{\boldsymbol{r}}) = \boldsymbol{r}\langle\boldsymbol{r},\dot{\boldsymbol{r}}\rangle - \dot{\boldsymbol{r}}\,|\boldsymbol{r}|^2 = r\dot r\,\boldsymbol{r} - r^2\dot{\boldsymbol{r}},$$

where the relation $\langle\boldsymbol{r},\dot{\boldsymbol{r}}\rangle = r\dot r$ is obtained by differentiating both sides of $r^2 = \langle\boldsymbol{r},\boldsymbol{r}\rangle$. Hence

$$\frac{d}{dt}(\boldsymbol{p}\times\boldsymbol{L}) = -\frac{\mu k}{r^3}\left(r\dot r\,\boldsymbol{r} - r^2\dot{\boldsymbol{r}}\right) = \mu k\left(\frac{\dot{\boldsymbol{r}}}{r} - \frac{\dot r}{r^2}\boldsymbol{r}\right) = \mu k\,\frac{d\hat{\boldsymbol{r}}}{dt}$$

(the last equality being the formula for the derivative of $\hat{\boldsymbol{r}}$ already used in the proof of Corollary 3.2). Therefore $d\boldsymbol{A}/dt = \boldsymbol{0}$.

Next we determine the direction and magnitude of $\boldsymbol{A}$. Both $\boldsymbol{p}\times\boldsymbol{L}$ and $\hat{\boldsymbol{r}}$ are orthogonal to $\boldsymbol{L}$, so $\boldsymbol{A}$ is a vector in the plane of motion. Taking the inner product with $\boldsymbol{r}$ and using the cyclic property of the scalar triple product,

$$\langle\boldsymbol{r}, \boldsymbol{p}\times\boldsymbol{L}\rangle = \langle\boldsymbol{L}, \boldsymbol{r}\times\boldsymbol{p}\rangle = \langle\boldsymbol{L},\boldsymbol{L}\rangle = L^2 ,$$

so that

$$\langle\boldsymbol{A},\boldsymbol{r}\rangle = L^2 - \mu k r .$$

If $\theta$ denotes the angle between $\boldsymbol{A}$ and $\boldsymbol{r}$, the left-hand side is $|\boldsymbol{A}|\,r\cos\theta$; solving for $r$ gives

$$r = \frac{L^2}{\mu k + |\boldsymbol{A}|\cos\theta} = \frac{L^2/(\mu k)}{1 + \bigl(|\boldsymbol{A}|/\mu k\bigr)\cos\theta} .$$

Comparing with the orbit equation of Theorem 5.2 we read off $\ell = L^2/(\mu k)$ and $e = |\boldsymbol{A}|/(\mu k)$, that is $|\boldsymbol{A}| = \mu k e$. Furthermore $r$ is smallest at $\theta=0$, in the direction of $\boldsymbol{A}$, so $\boldsymbol{A}$ points towards the periapsis.

∎

This computation deserves attention. Without solving a differential equation, merely by taking an inner product of conserved quantities, we obtained the orbit equation. Its physical meaning is equally clear: that $\boldsymbol{A}$ is a constant vector means the direction of the periapsis does not move. The reason the orbit closes is thus explained by a conserved quantity.

Conversely, when the force departs from the exact inverse square, $\boldsymbol{A}$ is not conserved and the periapsis rotates slowly. The advance of Mercury’s perihelion (the unexplained part of about 43 arcseconds per century) was due to precisely such a departure from $-k/r^2$, as predicted by general relativity. In Exercise 8.3 we compute explicitly how far the periapsis moves when a $1/r^2$ correction is added to the potential.

Remark 7.2

The conservation of $\boldsymbol{A}$ cannot be explained by ordinary spatial rotations, and is therefore called a hidden symmetry. The bound motion of the Kepler problem is known to possess the symmetry of the rotation group $SO(4)$ of four-dimensional space (Fock showed this in 1935 in the quantum theory of the hydrogen atom). The language for this symmetry is set up in Canonical transformations and Poisson brackets and Symmetry and conservation laws (Noether’s theorem). Note also that only one of the three components of $\boldsymbol{A}$ is independent, since they are constrained by the two relations $\langle\boldsymbol{A},\boldsymbol{L}\rangle = 0$ and $|\boldsymbol{A}|^2 = \mu^2k^2 + 2\mu EL^2$ (Exercise 8.4).

8. Exercises

Exercise 8.1Easy

Let $v_p$ be the speed of a planet on an elliptical orbit at periapsis and $v_a$ its speed at apoapsis. Using the periapsis distance $r_{\min}=a(1-e)$ and the apoapsis distance $r_{\max}=a(1+e)$, express the ratio $v_p/v_a$ in terms of the eccentricity $e$ alone. Then compute the ratio numerically for the Earth ($e = 0.0167$).

Solution

At periapsis and apoapsis the coordinate $r$ attains an extremum, so the radial velocity is $\dot r = 0$. The velocity therefore points along $\hat{\boldsymbol{e}}_\theta$ only, and the speed is $v = r|\dot\theta|$. By Theorem 3.1, $L = \mu r^2\dot\theta = \mu r v$ holds at both points, and since $L$ is conserved,

$$r_{\min}v_p = r_{\max}v_a \;\Longrightarrow\; \frac{v_p}{v_a} = \frac{r_{\max}}{r_{\min}} = \frac{1+e}{1-e}.$$

For the Earth, $(1+0.0167)/(1-0.0167) = 1.0167/0.9833 = 1.0340$. At perihelion (in early January) the orbital speed is about $3.4\%$ greater than at aphelion; the actual values are $30.29\ \mathrm{km/s}$ and $29.29\ \mathrm{km/s}$. This is the most direct manifestation of Kepler’s second law.

Exercise 8.2Standard

From the Earth’s orbital period $T = 3.156\times10^{7}\ \mathrm{s}$, its semi-major axis $a = 1.496\times10^{11}\ \mathrm{m}$ and the gravitational constant $G = 6.674\times10^{-11}\ \mathrm{m^3\,kg^{-1}\,s^{-2}}$, determine the mass of the Sun. The mass of the Earth may be neglected.

Solution

Neglecting $m_\oplus$ in the gravitational form $T^2 = 4\pi^2a^3/\bigl(G(M_\odot+m_\oplus)\bigr)$ of Theorem 6.1 gives

$$M_\odot = \frac{4\pi^2 a^3}{G T^2}.$$

Compute the numerator. Since $a^3 = (1.496\times10^{11})^3 = 3.348\times10^{33}$ and $4\pi^2 = 39.478$,

$$4\pi^2a^3 = 39.478\times3.348\times10^{33} = 1.3218\times10^{35}.$$

For the denominator, $T^2 = (3.156\times10^{7})^2 = 9.960\times10^{14}$, so

$$GT^2 = 6.674\times10^{-11}\times9.960\times10^{14} = 6.647\times10^{4}.$$

Hence

$$M_\odot = \frac{1.3218\times10^{35}}{6.647\times10^{4}} = 1.988\times10^{30}\ \mathrm{kg},$$

in agreement with the published value $1.989\times10^{30}\ \mathrm{kg}$. The same procedure gives the mass of a planet from the period and orbital radius of one of its satellites. This is the standard way of weighing celestial bodies.

Exercise 8.3Hard

Suppose the potential is $U(r) = -\dfrac{k}{r} - \dfrac{\beta}{r^2}$, where $\beta$ is a small constant. Find the increase in $\theta$ between one periapsis passage and the next. Show that it differs from $2\pi$, and compute the periapsis shift to first order in $\beta$.

Solution

The force is $f(r) = -U'(r) = -\dfrac{k}{r^2} - \dfrac{2\beta}{r^3}$, that is $f(1/u) = -ku^2 - 2\beta u^3$. Substituting into Lemma 5.1,

$$\frac{d^2u}{d\theta^2} + u = -\frac{\mu}{L^2u^2}\left(-ku^2 - 2\beta u^3\right) = \frac{\mu k}{L^2} + \frac{2\mu\beta}{L^2}u .$$

Moving the term in $u$ to the left-hand side,

$$\frac{d^2u}{d\theta^2} + \left(1 - \frac{2\mu\beta}{L^2}\right)u = \frac{\mu k}{L^2}.$$

Setting $\gamma^2 = 1 - 2\mu\beta/L^2$ (positive, since $\beta$ is small), this is exactly the same form of equation as in the proof of Theorem 5.2, with general solution

$$u(\theta) = \frac{\mu k}{\gamma^2L^2} + C\cos(\gamma\theta).$$

Now $u$ is maximal ($r$ minimal, i.e. at periapsis) when $\gamma\theta = 2\pi n$, so the angle between two consecutive periapsis passages is

$$\Delta\theta = \frac{2\pi}{\gamma} = \frac{2\pi}{\sqrt{1 - 2\mu\beta/L^2}} .$$

Expanding to first order in $\beta$ with $(1-x)^{-1/2} \approx 1 + x/2$ (see Theorem 5.3［Mean Value Theorems and Taylor's Theorem］ in The mean value theorem and Taylor’s theorem),

$$\Delta\theta \approx 2\pi\left(1 + \frac{\mu\beta}{L^2}\right) = 2\pi + \frac{2\pi\mu\beta}{L^2}.$$

Thus the periapsis advances by $\delta = 2\pi\mu\beta/L^2$ per revolution (in the same sense as the orbital motion if $\beta>0$). If $\beta=0$ then $\gamma=1$ and $\Delta\theta=2\pi$: the orbit closes, consistently with the conservation of $\boldsymbol{A}$ in Proposition 7.1. The relativistic advance of Mercury’s perihelion is an effect of the same kind, arising from a $1/r^3$ term added to the effective potential.

Exercise 8.4Hard

For the vector $\boldsymbol{A} = \boldsymbol{p}\times\boldsymbol{L} - \mu k\hat{\boldsymbol{r}}$ of Proposition 7.1, show that

$$|\boldsymbol{A}|^2 = \mu^2k^2 + 2\mu E L^2 .$$

Using this together with $|\boldsymbol{A}| = \mu k e$, rederive Corollary 5.4.

Solution

Compute $|\boldsymbol{A}|^2 = |\boldsymbol{p}\times\boldsymbol{L}|^2 - 2\mu k\langle\boldsymbol{p}\times\boldsymbol{L}, \hat{\boldsymbol{r}}\rangle + \mu^2k^2$ term by term.

First term. Since $\boldsymbol{p}$ and $\boldsymbol{L}$ are orthogonal ($\langle\boldsymbol{L},\boldsymbol{p}\rangle = \langle\boldsymbol{r}\times\boldsymbol{p},\boldsymbol{p}\rangle = 0$), we have $|\boldsymbol{p}\times\boldsymbol{L}| = |\boldsymbol{p}||\boldsymbol{L}|$ and hence $|\boldsymbol{p}\times\boldsymbol{L}|^2 = p^2L^2$.

Second term. By the cyclic property of the scalar triple product,

$$\langle\boldsymbol{p}\times\boldsymbol{L}, \hat{\boldsymbol{r}}\rangle = \frac{1}{r}\langle\boldsymbol{p}\times\boldsymbol{L},\boldsymbol{r}\rangle = \frac{1}{r}\langle\boldsymbol{L},\boldsymbol{r}\times\boldsymbol{p}\rangle = \frac{L^2}{r}.$$

Altogether,

$$|\boldsymbol{A}|^2 = p^2L^2 - \frac{2\mu kL^2}{r} + \mu^2k^2 = 2\mu L^2\left(\frac{p^2}{2\mu} - \frac{k}{r}\right) + \mu^2k^2 .$$

The bracket is exactly the mechanical energy $E = p^2/(2\mu) - k/r$ (Proposition 4.2). Hence $|\boldsymbol{A}|^2 = \mu^2k^2 + 2\mu EL^2$.

Substituting $|\boldsymbol{A}| = \mu k e$ gives $\mu^2k^2e^2 = \mu^2k^2 + 2\mu EL^2$, and dividing both sides by $\mu^2k^2$,

$$e^2 = 1 + \frac{2EL^2}{\mu k^2},$$

which reproduces Corollary 5.4. The fact that the eccentricity is fixed by the energy and the angular momentum has emerged naturally as the norm of a conserved quantity.

References

• H. Goldstein, C. Poole, J. Safko, Classical Mechanics, 3rd ed., Addison-Wesley, 2002 — Chapter 3, “The Central Force Problem”. Treats the effective potential, Binet’s equation, the Laplace–Runge–Lenz vector and Bertrand’s theorem systematically.
• L. D. Landau, E. M. Lifshitz, Mechanics, 3rd ed., Butterworth-Heinemann, 1976 — Chapter III, “Integration of the equations of motion”. The classic account, reaching the reduction of the two-body problem and the Kepler problem by the shortest route.
• Harashima Akira, Rikigaku I (Mechanics I), Shokabo, 1972 (in Japanese) — the chapters on central forces and planetary motion. A standard introductory text available in Japanese.
• Yamamoto Yoshitaka, Koten Rikigaku no Keisei: Newton kara Lagrange e (The Formation of Classical Mechanics: From Newton to Lagrange), Nippon Hyoron Sha, 1997 (in Japanese) — follows, from the primary sources, the historical route by which the inverse-square law of gravitation was derived from Kepler’s three laws.
• V. I. Arnold, Mathematical Methods of Classical Mechanics, 2nd ed., Springer, 1989 — Chapter 2. Treats motion in a central field from a differential-geometric viewpoint and discusses when orbits close.

Appendix: Kepler’s equation, linking time to position

The problem that remains. In the main text we determined completely the shape $r(\theta)$ of the orbit, but not when the planet is where. It should suffice to integrate the constancy of the areal velocity of Theorem 3.3, but integrating $dt = (\mu/L)\,r(\theta)^2 d\theta$ directly over $\theta$ does not produce an inverse function in closed form. What helps here is an auxiliary variable, the eccentric anomaly.

Introducing the eccentric anomaly. For an elliptical orbit define the variable $\psi$ by

$$r = a(1 - e\cos\psi).$$

Here $\psi=0$ corresponds to periapsis ($r = a(1-e)$) and $\psi=\pi$ to apoapsis ($r = a(1+e)$), and $\psi$ covers the range of $r$ exactly once. Geometrically it is the central angle obtained by projecting the point of the orbit perpendicularly onto the circle of radius $a$ circumscribing the ellipse (the auxiliary circle).

Integrating the time. Solve the energy relation of Proposition 4.2 for $\dot r$ and substitute $E = -k/(2a)$ from Corollary 5.4 together with $L^2 = \mu k a(1-e^2)$ (which follows from $\ell = a(1-e^2)$ and $\ell = L^2/\mu k$ in Definition 5.3):

$$\dot r^2 = \frac{2}{\mu}\left(E + \frac{k}{r}\right) - \frac{L^2}{\mu^2r^2} = \frac{k}{\mu a r^2}\left(2ar - r^2 - a^2(1-e^2)\right) = \frac{k}{\mu a r^2}\left(a^2e^2 - (r-a)^2\right).$$

The last equality is the completion of the square $2ar - r^2 - a^2 + a^2e^2 = -(r-a)^2 + a^2e^2$. Inserting $r - a = -ae\cos\psi$ turns the bracket into $a^2e^2\sin^2\psi$. On the other hand $\dot r = ae\sin\psi\,\dot\psi$, so equating the two sides,

$$a^2e^2\sin^2\psi\,\dot\psi^2 = \frac{k}{\mu a r^2}\,a^2e^2\sin^2\psi \;\Longrightarrow\; \dot\psi = \frac{1}{r}\sqrt{\frac{k}{\mu a}} = \frac{1}{a(1-e\cos\psi)}\sqrt{\frac{k}{\mu a}} .$$

Separating variables and integrating from the time $t_p$ of periapsis passage,

$$\int_0^{\psi}(1-e\cos\psi')\,d\psi' = \sqrt{\frac{k}{\mu a^3}}\,(t-t_p),$$

that is,

$$\psi - e\sin\psi = n\,(t-t_p),\qquad n = \sqrt{\frac{k}{\mu a^3}} = \frac{2\pi}{T}$$

(the last equality is Theorem 6.1 itself). The right-hand side $M = n(t-t_p)$ is called the mean anomaly, and the relation

$$M = \psi - e\sin\psi$$

is Kepler’s equation.

Numerical solution. Kepler’s equation cannot be solved for $\psi$ in elementary functions. In practice one solves it by Newton’s method. Setting $g(\psi) = \psi - e\sin\psi - M$, we have $g'(\psi) = 1 - e\cos\psi \ge 1-e > 0$, so for $0\le e<1$ the function $g$ is strictly increasing, the solution is unique, and Newton’s method converges stably.

import numpy as np

def solve_kepler(M, e, tol=1e-12, max_iter=60):
    """Solve Kepler's equation M = psi - e*sin(psi) by Newton's method."""
    psi = M if e < 0.8 else np.pi
    for _ in range(max_iter):
        g = psi - e * np.sin(psi) - M
        dpsi = -g / (1.0 - e * np.cos(psi))
        psi += dpsi
        if abs(dpsi) < tol:
            break
    return psi

def position(t, a, e, n, t_p=0.0):
    """Return the angle theta from periapsis and the radius r at time t."""
    psi = solve_kepler(n * (t - t_p), e)
    r = a * (1.0 - e * np.cos(psi))
    theta = 2.0 * np.arctan2(np.sqrt(1 + e) * np.sin(psi / 2),
                             np.sqrt(1 - e) * np.cos(psi / 2))
    return theta, r

The final expression for theta is the relation $\tan(\theta/2) = \sqrt{(1+e)/(1-e)}\,\tan(\psi/2)$ between the true anomaly and the eccentric anomaly, written with arctan2 so that the quadrant is handled correctly. That relation is derived by equating $r=a(1-e\cos\psi)$ with $r = a(1-e^2)/(1+e\cos\theta)$, expressing $\cos\theta$ through $\cos\psi$, and applying the half-angle formulae.

Computing an ephemeris is this procedure repeated. One determines the orbital elements $(a, e, \ldots)$ from observation and then reproduces the position at any time from Kepler’s equation. We perform daily the exact inverse of the operation by which Kepler read his laws out of tables of observations.