Deriving Kepler's Laws

Kepler's laws describe the motion of objects in the presence of a central inverse square force. For simplicity, we'll consider the motion of the planets in our solar system around the Sun, with gravity as the central force. Among other things, Kepler's laws allow one to predict the position and velocity of the planets at any given time, the time for a satellite to collapse into the surface of a planet, and the period of a planet's orbit as a function of its orbits' geometry. Though the laws were originally obtained by Kepler after careful analysis of empirical data, the complete understanding was missing until Newton derived each law as pieces of his orbital mechanics. In his footsteps we will obtain each law in turn, as we consider the orbit of a planet in the gravity of a massive star.

Kepler's laws of planetary motion state that

• A planet moves around the Sun in an elliptical path with the Sun as one of the focii.
• The line segment joining a planet and the Sun sweeps out equal areas during equal intervals of time, i.e. $\frac12\int r^2 d\theta \sim Lt$, where $L$ is a constant.
• The square of the orbital time period of a planet is proportional to the cube of the semi-major axis of its orbit, i.e. $T^2\propto r^3.$

Here, we list the basic assumptions underlying orbital mechanics:

• As the Sun, with mass $M_\odot$, is very large compared to any other object in the solar system, its motion is essentially unaffected by the gravity of the planets.
• The gravity of the Sun acts along the line between the Sun and a given planet ($F$ acts along $\hat{r}$). As a result, the motion of the planet is confined to a 2D plane.
• We assume that collisions with space dust and other methods of energy dissipation are negligible, so that the mechanical energy $E$ is a conserved quantity.

The central differential equation that describes planetary motion can be written as

$$m\ddot{\vec{r}} = - G\frac{M_\textrm{Sun}m}{r^2}\hat{r}.$$

This is a vector equation in two dimensions. The left hand side describes the kinematics of our object whose position relative to the Sun is given by $\vec{r}$, and the right hand side describes the force of gravity, which depends on the separation $\vec{r}$ only through the square of its magnitude.

Although this problem can be solved in a straightforward fashion in polar coordinates with unit vectors for the radius $\hat{r}$, and the angle about the Sun, $\hat{\theta}$, it requires us to keep track of some tricky infinitesimal quantities. To avoid this needless complication, we change over to Cartesian coordinates for the purposes of calculating our derivatives.

First, partially re-express the problem in the $\left(x,y\right)$ coordinate system. Notice that the $\cos$ and $\sin$ of the angle $\theta$ are given by $\frac{x}{r}$ and $\frac{y}{r},$ respectively. We recast the central equation in the form

$$\begin{aligned} m\ddot{x} &= -G\frac{M_\textrm{Sun}m}{r^2}\cos\theta \\ m\ddot{y} &= -G\frac{M_\textrm{Sun}m}{r^2} \sin\theta. \end{aligned}$$

We need to find the second derivative of the $x$ and $y$ coordinates in terms of the polar coordinates.

For $x$, we have

$$\begin{aligned} \dot{x} &= \dot{r}\cos\theta - r\dot{\theta}\sin\theta \\ \ddot{x} &= \ddot{r}\cos\theta -2 \dot{r}\dot{\theta}\sin\theta - r\dot{\theta}^2\cos\theta - r\ddot{\theta}\sin\theta. \end{aligned}$$

For $y$, we have

$$\begin{aligned} \dot{y} &= \dot{r}\sin\theta + r\dot{\theta}\cos\theta \\ \ddot{y} &= \ddot{r}\sin\theta + 2 \dot{r}\dot{\theta}\cos\theta - r\dot{\theta}^2\sin\theta + r\ddot{\theta}\cos\theta. \end{aligned}$$

We now obtain the orbital equations in polar coordinates by a trick applied in two different ways.

First, we multiply $\ddot{x}$ by $\cos\theta$, and $\ddot{y}$ by $\sin\theta$, and add them.

From the central equation, we have

$$\ddot{x}\cos\theta + \ddot{y}\sin\theta = -GM_\textrm{Sun}\left(\cos^2\theta + \sin^2\theta\right)\frac{1}{r^2} = -GM_\textrm{Sun}\frac{1}{r^2}.$$

From the derivative identities between Cartesian and polar coordinates, we have

$$\begin{aligned} \ddot{x}\cos\theta + \ddot{y}\sin\theta = &\ddot{r}\cos^2\theta -2 \dot{r}\dot{\theta}\cos\theta\sin\theta - r\dot{\theta}^2\cos^2\theta - r\ddot{\theta}\cos\theta\sin\theta \\ &+ \ddot{r}\sin^2\theta + 2 \dot{r}\dot{\theta}\sin\theta\cos\theta - r\dot{\theta}^2\sin^2\theta + r\ddot{\theta}\sin\theta\cos\theta. \end{aligned}$$

We see that every term with a $\sin\theta\cos\theta$ cancels so that we're left with

$$\ddot{x}\cos\theta + \ddot{y}\sin\theta = \ddot{r} - r\dot{\theta}^2.$$

Using the result from the central equations, we have

$$\boxed{\left(\ddot{r} - r\dot{\theta}^2\right) = -GM_\textrm{Sun}\frac{1}{r^2}}.$$

We notice that if we insist that $r$ is constant and $\ddot{r}$ is zero, then this is just the equation for a circular orbit. Allowing for $r$ to vary opens our problem up to more general orbits like ellipses and hyperbolas.

If we had instead multiplied $\ddot{x}$ by $\sin\theta$, $\ddot{y}$ by $\cos\theta$, and subtract the equations, we would find

$$r\ddot{\theta} + 2\dot{r}\dot{\theta} = 0.$$

If we multiply this equation by $r$, we find $r^2\ddot{\theta} + 2r\dot{r}\dot{\theta}=0$. However, this is just the time derivative of $r^2\dot{\theta}$, and thus we have shown

$$\frac{d}{dt}r^2\dot{\theta} = 0.$$

But $mr^2\dot{\theta}$ is the angular momentum of the planet. Thus, the angular momentum of the planet is conserved.

This result is somewhat anti-climactic. For one thing, gravity acts along the displacement vector between the Sun and the planet, and thus there is no torque on the system, and the angular momentum must be conserved. On a deeper level, if we wrote down the Hamiltonian for the system, we'd see it has no dependence on $\theta$, and thus the momentum associated with $\theta$ must be a constant of the motion.

If we integrate this equation with respect to time, we find that $mr^2\dot{\theta} = L,$ where $L$ is a constant. Integrating once more in time, we find

$$\begin{aligned} \int L dt &= m\int r^2 \frac{d\theta}{dt} dt\\ &= m\int_{\theta_i}^{\theta_f} r^2 d\theta. \end{aligned}$$

The integral $\frac12 \int r^2 d\theta$ is the area swept out by the radial vector from the Sun to the planet in moving from $\theta_i$ to $\theta_f$. However, the result is independent of $\theta_i$ and $\theta_f$, but it only depends on $t$ since the angular momentum is constant.

Thus, we have derived Kepler's second law, i.e. segments of orbits sweep out equal areas in equal intervals of time:

$$\boxed{\displaystyle\frac12\int_{\text{any path}} r^2 d\theta = \frac{Lt}{2m}}.$$

The speed of a certain planet at the perihelion is $v_p$ and, at this position, the distance of the sun from the planet is $r_p$. Relate $\{r_p,v_p\}$ to the corresponding quantities at the aphelion $\{r_a,v_a\}$.

---

The magnitude of the angular momentum at perihelion is $L_p=m_pr_pv_p$ because $r_p$ and $v_p$ are mutually perpendicular. Similarly, $L_a=m_pr_av_a$. Using the conservation of angular momentum,

$$m_pr_pv_p=m_pr_av_a\Rightarrow r_pv_p=r_av_a \Rightarrow \dfrac{v_p}{v_a} = \dfrac{r_a}{r_p}.$$

To make progress, we need to solve our central equation for $r$. We have

$$\ddot{r} - r\dot{\theta}^2 = -GM\frac{1}{r^2}.$$

The differential equation becomes easy to solve if we make the substitution $r \rightarrow u^{-1}$. We have

$$\begin{aligned} \frac{dr}{dt} &= -u^{-2}\frac{du}{dt} \\ &= -u^{-2}\frac{du}{d\theta} \frac{Lu^2}{m} \\ &= -\frac{L}{m}\frac{du}{d\theta}. \end{aligned}$$

Further, we have

$$\begin{aligned} \frac{d^2r}{dt^2} &= -\frac{L}{m}\frac{d}{dt}\frac{du}{d\theta} \\ &= -\frac{L}{m} \frac{d\theta}{dt}\frac{d}{d\theta}\frac{du}{d\theta} \\ &= -\left(\frac{L}{m}\right)^2u^2\frac{d^2u}{d\theta^2}. \end{aligned}$$

With this identity in hand, our central equation becomes

$$-\left(\frac{L}{m}\right)^2u^2\frac{d^2u}{d\theta^2} - \left(\frac{L}{m}\right)^2u^3 = -GMu^2,$$

which becomes

$$-\frac{d^2u}{d\theta^2} + \frac{GMm^2}{L^2} = u,$$

which has the simple solution $u = A\cos\left(\theta + \delta\right) + \frac{GMm^2}{L^2}$.

We can always define our coordinates so that $\delta=0$, and thus we set it to zero for the remainder of the discussion.

In terms of $r(\theta)$, we have

$$\begin{aligned} r &= \frac{1}{A\cos\theta+ \frac{GMm^2}{L^2}} \\ &= \frac{1}{\frac{GMm^2}{L^2}\left(e\cos\theta + 1\right)}. \end{aligned}$$

Note that in the line above we made the useful substitution $e \Rightarrow \dfrac{AL^2}{GMm^2}$.

From this form of $r$, it is clear that the aphelion and perihelion (points of furthest and closest distance, respectively, to the Sun) are given by

$$r_\text{max} =\dfrac{1}{\dfrac{GMm^2}{L^2}\left(1 - e\right)}, \quad r_\text{min} = \dfrac{1}{\dfrac{GMm^2}{L^2}\left(1 + e\right)}.$$

The semi-major axis is given by

$$\boxed{a=\dfrac12\left(r_\text{min} + r_\text{max}\right)=\dfrac{L^2}{GMm^2}\left(1-e^2\right)^{-1}}.$$

We see that the orbit is given by an ellipse as Kepler found from Brahe's dataset. Moreover, since $r_\text{min}$ and $r_\text{max}$ are distances from the Sun, we see that the Sun is at one focus of the orbit. Thus, we have derived Kepler's first law.

$r\sim\left(1+e\cos\theta\right)^{-1}$ is the general form of an ellipse in polar coordinates, with the origin placed at a focus. In the study of ellipses, the parameter $e$ is often called the eccentricity. When the eccentricity of a planet's orbit is zero, the orbit is perfectly circular. As $e$ approaches one, the orbit is stretched out into more elongated elliptical trajectories. To demonstrate this feature, we plot the orbit below for several values of the eccentricity, $e$.

$$\text{Plot of the orbit for increasing eccentricity}, e$$

What happens if $e \geq 1$?

If in the expression $\displaystyle\int_\text{any path} r^2d\theta = Lt/m$ we take the path to be one complete orbit of the Sun, $t=T$ and the area swept out by the radial vector is the area of the elliptical orbit, $A=2\pi ab$. Here $a$ and $b$ are the semi-major and semi-minor axes of the elliptical orbit.

From the expression for $a$ obtained above, we can see that the square of the angular momentum is equal to the semi-major axis of the elliptical orbit multiplied by some constants, $L^2 \propto a$. This is the crucial information we need in order to obtain the third law.

Squaring both sides in $\displaystyle\int_\text{path} r^2\theta = LT/m = 2\pi ab$, we have $L^2T^2\propto a^2b^2$. Now, $a$ and $b$ are simply related since they're both linear dimensions of the fixed elliptical orbit, so they are proportional and thus we have $L^2T^2 \propto a^4$.

Finally, we showed that $L^2 \propto a$, so we have $T^2\propto a^3$, which is Kepler's third law.

Were we to do more careful record keeping in the analysis above, we could obtain the factor of $\dfrac{4\pi^{2}}{GM}$, to get an exact statement of the third law:

$$\boxed{\displaystyle T^2 = \frac{4\pi^{2}}{GM}a^3}.$$

Note that this law holds for all elliptical orbits, regardless of their eccentricities.