⮜ Articles

The Kepler Problem

This article provides a detailed mathematical exploration of the Kepler Problem, which describes the motion of two bodies influenced by gravitational forces. It systematically derives and explains the principles underlying the motion, including the conservation of energy and angular momentum, and culminates in the establishment of Kepler's three laws of planetary motion. Key topics include:

• Setting up the problem with center-of-mass reference frames and vector notations.
• Deriving expressions for kinetic and potential energy in the context of celestial mechanics.
• Formulating the Lagrangian and its applications to conserved quantities.
• Mathematical proof and derivation of Kepler's First, Second, and Third Laws.

The content is rigorous and technical, aimed at readers with a background in physics or mathematics who seek an in-depth understanding of celestial mechanics.

Notation

On this page, I use the following notational conventions.

• Unbolded vectors: if $\mathbf a$ is a vector, then $a$ shall denote the magnitude $|\mathbf a|$.
• Overdot notation: if $w(t)$ and $\dot w(t)$ are functions of time, then $\dot w = \frac{dw}{dt}$.
• Gradient: if $f$ is a scalar function of a vector $\mathbf a$, then $\frac{\partial f}{\partial \mathbf a}$ denotes the gradient of $f$ with respect to $\mathbf a$.

The Situation

Consider two moving planets with properties represented through the symbols summarized in the table below. $$ \begin{array}{c|cc} & \text{planet 1} & \text{planet 2} \\ \hline \text{mass} & m_1 & m_2 \\ \text{position} & \mathbf x_1 & \mathbf x_2 \\ \text{velocity} & \dot{\mathbf x}_1 & \dot{\mathbf x}_2 \end{array} $$ Let us presume that we are observing these two planets from a reference frame wherein the center of mass is stationary.

\begin{equation} \label{eq:center-of-mass} \frac{m_1\dot{\mathbf x}_1 + m_2 \dot{\mathbf x}_2}{m_1+m_2} = \mathbf{0} \end{equation} Note that a consequence of this restriction is that the velocities $\dot{\mathbf x}_1$ and $\dot{\mathbf x}_2$ are parallel.

Kinetic Energy

The total kinetic energy $K$ is the sum of the individual kinetic energies of the planets.

\begin{equation} K = \frac{1}{2}m_1 \dot x_1^2 + \frac{1}{2}m_2 \dot x_2^2 \end{equation}

This expression can be generally rewritten in the following manner.

Show intermediate steps

\begin{gather} K = \frac{1}{2}\frac{(m_1 + m_2) m_1 \dot x_1^2 + (m_1 + m_2) m_2 \dot x_2^2}{m_1+m_2} \\[1em] K = \frac{1}{2}\frac{m_1^2 \dot x_1^2 + m_2^2 \dot x_2^2 + 2m_1m_2 \dot{\mathbf x}_1 \cdot \dot{\mathbf x}_2}{m_1+m_2} + \frac{1}{2} \frac{m_1 m_2\dot{x}_1^2 + m_1 m_2 \dot x_2^2 - 2m_1m_2 \dot{\mathbf x}_1 \cdot \dot{\mathbf x}_2}{m_1 + m_2} \end{gather}

\begin{equation} K = \frac{1}{2}(m_1+m_2) \left| \frac{m_1\dot{\mathbf x}_1 + m_2 \dot{\mathbf x}_2}{m_1+m_2} \right|^2 + \frac{1}{2}\left( \frac{m_1 m_2}{m_1 + m_2} \right) |\dot{\mathbf x}_1 - \dot{\mathbf x}_2|^2 \end{equation}

The first term vanishes in the frame of reference of the center of the mass according to equation \eqref{eq:center-of-mass}.

\begin{equation} K = \frac{1}{2}\left( \frac{m_1 m_2}{m_1 + m_2} \right) |\dot{\mathbf x}_1 - \dot{\mathbf x}_2|^2 \end{equation}

Define the reduced mass $M = \frac{m_1 m_2}{m_1 + m_2}$ and the relative position vector $\mathbf x = \mathbf x_1 - \mathbf x_2$ in order to express the kinetic energy in a more compact form.

\begin{equation} \label{eq:K} K = \frac{1}{2} M\dot{x}^2 \end{equation}

Potential Energy

The gravitational potential energy $U$ is given through Newton's Law of Universal Gravitation, an empirical law.

\begin{equation} U = -\frac{Gm_1 m_2}{|\mathbf x_1 - \mathbf x_2|} \text{, where } G \approx 6.67\times 10^{-11} \frac{\text{N} \cdot \text{m}^2}{\text{kg}^2} \end{equation}

For simplicity, let $\alpha = G m_1 m_2$ and also adopt the relative position vector $\mathbf x = \mathbf x_1 - \mathbf x_2$ from the previous section.

\begin{equation} \label{eq:U} U = -\frac{\alpha}{x} \end{equation}

The Lagrangian

The Lagrangian $L$ can be written from its definition. \begin{equation} \label{eq:L} L \overset{\text{def}}{=} K - U \overset{\eqref{eq:K}\eqref{eq:U}}{=} \frac{1}{2} M\dot{x}^2 + \frac{\alpha}{x} \end{equation} The gradient of the Lagrangian with respect to the velocity vector $\frac{\partial L}{\partial \dot{\mathbf x}}$ will be important for the section that follows.

Show intermediate steps

\begin{equation} L = \frac{1}{2} M \dot{\mathbf x} \cdot \dot{\mathbf x} + \frac{\alpha}{x} \end{equation}

The gradient $\frac{\partial L}{\partial \dot{\mathbf x}}$ is defined as the unique vector field such that for any vector $\mathbf v$ the dot product $\frac{\partial L}{\partial \dot{\mathbf x}} \cdot \mathbf v$ represents the directional derivative of $L$ along $\mathbf v$.

\begin{gather} \frac{\partial L}{\partial \dot{\mathbf x}} \cdot \dot{\mathbf x} \overset{\text{def}}{=} \lim_{\epsilon \to 0} \frac{\left( \frac{1}{2} M \left( \dot{\mathbf x} + \epsilon \dot{\mathbf x} \right) \cdot \left( \dot{\mathbf x} + \epsilon \dot{\mathbf x} \right) + \frac{\alpha}{x} \right) - \left( \frac{1}{2} M \dot{\mathbf x} \cdot \dot{\mathbf x} + \frac{\alpha}{x} \right) }{\epsilon} \\[0.5em] \frac{\partial L}{\partial \dot{\mathbf x}} \cdot \dot{\mathbf x} = \lim_{\epsilon \to 0} \frac{\left( \epsilon + \frac{1}{2}\epsilon^2 \right) M \dot{\mathbf x} \cdot \dot{\mathbf x}}{\epsilon} \\[0.5em] \frac{\partial L}{\partial \dot{\mathbf x}} \cdot \dot{\mathbf x} = M \dot{\mathbf x} \cdot \dot{\mathbf x} \end{gather}

\begin{equation} \label{eq:gradient} \frac{\partial L}{\partial \dot{\mathbf x}} = M \dot{\mathbf x} \end{equation}

Conserved Quantities

Angular Momentum

The Lagrangian $L$ is invariant under spatial rotations as it is only dependent on the vector magnitudes $x$ and $\dot x$. Hence, angular momentum \begin{equation} \label{eq:angular-momentum} \mathbf l \overset{\text{def}}{=} \mathbf x \times \frac{\partial L}{\partial \dot{\mathbf x}} \overset{\eqref{eq:gradient}}{=} M \mathbf x \times \dot{\mathbf x} \end{equation} is constant, which reflects the familiar law of conservation of angular momentum.

Energy

Since $\frac{\partial L}{\partial t}=0$, we have that the energy \begin{equation} \label{eq:energy} E \overset{\text{def}}{=} \frac{\partial L}{\partial \dot{\mathbf x}} \cdot \dot{\mathbf x} - L \overset{\eqref{eq:L}}{=} \frac{1}{2}M\dot x^2 - \frac{\alpha}{x} \end{equation} is constant, which reflects the familiar law of conservation of energy.

Kepler's 1st Law: Elliptical Orbits

Polar Coordinates

From equation \eqref{eq:angular-momentum}, it follows that $\mathbf x, \dot{\mathbf x} \perp \mathbf l$, so since $\mathbf l$ is constant, we must have that $\mathbf x$ and $\dot{\mathbf x}$ remain in a single plane as time evolves. Given that the vectors of the system are restricted to a plane, let's set up a system of polar coordinates $(r, \theta)$ in that plane. The origin of this coordinate system will correspond to $\mathbf x = 0$. First, let's establish some key facts about this polar coordinate system. \begin{align} \mathbf x &= \begin{bmatrix} r \vphantom{\dot \theta} \\ 0 \vphantom{\dot \theta} \\ 0 \vphantom{\dot \theta} \end{bmatrix}, & \dot{\mathbf x} &= \begin{bmatrix} \dot r \vphantom{\dot \theta} \\ r \dot \theta \\ 0 \vphantom{\dot \theta} \end{bmatrix}, & x &= r, & \dot x^2 &= \dot r^2 + r^2 \dot \theta^2, & |\mathbf x \times \dot{\mathbf x}| &= r^2 \dot \theta \end{align} From this, we can write the magnitude of the angular momentum $l$ from equation \eqref{eq:angular-momentum} and energy $E$ from equation \eqref{eq:energy} each in polar coordinate representations. \begin{equation} \label{eq:angular-momentum2} l = M r^2 \dot \theta \end{equation} \begin{equation} E = \frac{1}{2} M \dot{r}^2 + \frac{1}{2}Mr^2 \dot{\theta}^2 - \frac{\alpha}{r} \end{equation} Combining these equations to eliminate $\dot \theta$ yields a differential equation for $r = r(t)$. \begin{equation} \label{eq:r-only} E = \frac{1}{2} M \dot{r}^2 + \frac{1}{2} \frac{l^2}{Mr^2} - \frac{\alpha}{r} \end{equation}

Solving the Differential Equation

Rearrange equation \eqref{eq:r-only} to provide an expression for $\dot r$ in terms of $r$.

\begin{equation} \dot r = \sqrt{\frac{2E}{M} - \frac{l^2}{M^2r^2} + \frac{2\alpha}{Mr}} \end{equation}

Dividing this result by $\dot \theta = \frac{l}{Mr^2}$ from equation \eqref{eq:angular-momentum2} yields an expression for $\frac{dr}{d\theta}$, which yields a separable differential equation for the function $r = r(\theta)$.

\begin{equation} \frac{dr}{d\theta} = \frac{\dot r}{\dot \theta} = \frac{\sqrt{\frac{2E}{M} - \frac{l^2}{M^2r^2} + \frac{2\alpha}{Mr}}}{\frac{l}{Mr^2}} \end{equation}

Through some remarkable and nontrivial algebraic rearrangement, this can be reexpressed as follows.

\begin{equation} d\theta = -\frac{\displaystyle -\frac{\frac{l^2}{M\alpha} \frac{1}{r^2}}{\sqrt{1 + \frac{2El^2}{M\alpha^2}}}}{\displaystyle \sqrt{1 - \left( \frac{\frac{l^2}{M\alpha} \frac{1}{r}-1}{\sqrt{1 + \frac{2El^2}{M\alpha^2}}} \right)^2}}\ dr \end{equation}

Let $p = \frac{l^2}{M\alpha}$ and $e = \sqrt{1 + \frac{2El^2}{M\alpha^2}}$. Integrate the result.

\begin{equation} \theta = \int -\frac{-\frac{p}{er^2}}{\sqrt{1 - \left( \frac{ \frac{p}{r}-1}{e} \right)^2}}\ dr = \arccos \left( \frac{ \frac{p}{r}-1}{e} \right) + C \end{equation}

Without loss of generality, set $C = 0$ and rearrange for $r = r(\theta)$.

\begin{equation} r(\theta) = \frac{p}{1 + e\cos \theta} \text{ where } \begin{cases} p = \frac{l^2}{M\alpha} \vphantom{\sqrt{\frac{a^2}{b^2}}} & \text{semi-latus rectum} \\[0.5em] e = \sqrt{1 + \frac{2El^2}{M\alpha^2}} & \text{eccentricity} \\[0.5em] M = \frac{m_1 m_2}{m_1 + m_2} \vphantom{\sqrt{\frac{a^2}{b^2}}} & \text{reduced mass} \\[0.5em] \alpha = G m_1 m_2 \vphantom{\sqrt{\frac{a^2}{b^2}}} \\[0.5em] G \approx 6.67\times 10^{-11} \frac{\text{N} \cdot \text{m}^2}{\text{kg}^2} & \text{gravitational constant} \\[0.5em] E = \frac{1}{2} M \dot{r}^2 + \frac{1}{2} \frac{l^2}{Mr^2} - \frac{\alpha}{r} & \text{energy (conserved)} \vphantom{\sqrt{\frac{a^2}{b^2}}} \\[0.5em] l = Mr^2 \dot \theta & \text{angular momentum (conserved)} \vphantom{\sqrt{\frac{a^2}{b^2}}} \end{cases} \end{equation}

Therefore, the curve swept out by $\mathbf x = \mathbf x_1 - \mathbf x_2$ is an ellipse with latus rectum $p$, eccentricity $e$, and one focus at the origin of the polar coordinate system. This establishes Kepler's First Law.

Kepler's 2nd Law: Equal Areas in Equal Time

Note that polar area $A$ is swept out throughout the orbit at the rate \begin{equation} \frac{dA}{dt} = \frac{1}{2}r^2 \dot{\theta} = \frac{l}{2M}. \end{equation} Since angular momentum $l$ is constant, then the rate at which area is swept out by the orbit $\frac{dA}{dt}$ must also be constant. In other words, equal polar areas are swept out in equal times, which is Kepler's Second Law.

Kepler's 3rd Law: A Proportionality Law for Period and Orbital Size

For some constants $\xi$ and $\eta$, perform the space and time scalings $\mathbf x \mapsto \xi \mathbf x$ and $t \mapsto \eta t$ on the Lagrangian as shown in equation \eqref{eq:L}. \begin{equation} L = \frac{\xi^2}{\eta^2} \frac{1}{2} M \dot{x}^2 + \frac{1}{\xi} \frac{\alpha}{x} \end{equation} This new Lagrangian $L$ must be a constant multiple of the Lagrangian shown in equation \eqref{eq:L} in order to yield the same dynamics. Hence, we must conclude that $\frac{\xi^2}{\eta^2} = \frac{1}{\xi}$, which implies \begin{equation} \xi^3 = \eta^2. \end{equation} The constant $\xi$ can be interpreted as a scale factor of orbital size and $\eta$ as a scale factor of orbital period for this particular transformation. Hence, this relationship establishes Kepler's Third Law. Commonly in the literature, one writes $\xi = \frac{r_1}{r_2}$ and $\eta = \frac{T_1}{T_2}$, where $r_1$ and $r_2$ are corresponding length measurements (such as orbital radii, in the case of a circular orbit) and $T_1$ and $T_2$ are corresponding time measurements (such as orbital periods) in two geometrically similar orbits. \begin{equation} \left( \frac{r_1}{r_2} \right)^3 = \left( \frac{T_1}{T_2} \right)^2 \end{equation} This is Kepler's 3rd Law.