# Deriving Kepler's Laws Using Quaternions

Christopher J. Abel

###### Abstract

In the past, Kepler painstakingly derived laws of planetary motion using
difficult to understand and hard to follow techniques. In 1843 William
Hamilton created and described the quaternions, which extend the complex
numbers and can easily describe rotations in three dimensional space. In this
article, we will harness this system to provide a new and intuitive way to
derive Kepler’s laws. This will include using a quaternionic version of the
spatial Kepler problem differential equation, and using the general solution
to describe the motion of planets orbiting a central body. We use the standard
method for regularizing celestial mechanics, but this article will be solely
focused on showing the validity of Kepler’s laws.

##  I The Star Conjugate

In order to describe the KS transformation in the form of quaternions we
require a new operation that will help remove the fourth dimension for any
general quaternion. This setup was suggested by Waldvogel (2006). We define
the star conjugate of a quaternion $\boldsymbol{q}\in\mathbb{H}$ to be

$\boldsymbol{q}^{*}=q_{0}+q_{1}i+q_{2}j-q_{3}k.$ (1)

Note that all the star conjugate does is negate the $k$ term. We can easily
represent this in terms of the quaternion conjugate,

$\boldsymbol{q}^{*}=(-k)\overline{\boldsymbol{q}}k.$ (2)

One may easily show that for any $\boldsymbol{q},\boldsymbol{v}\in\mathbb{H}$,

$(\boldsymbol{q}^{*})^{*}=\boldsymbol{q},\qquad|\boldsymbol{q}^{*}|=|\boldsymbol{q}|,\qquad(\boldsymbol{q}\boldsymbol{v})^{*}=\boldsymbol{v}^{*}\boldsymbol{q}^{*}.$
(3)

We now define a mapping between the quaternion algebra to take advantage of
the star conjugate. Consider the mapping,

$\boldsymbol{q}\mapsto\boldsymbol{x}=\boldsymbol{q}\boldsymbol{q}^{*}.$ (4)

Immediately we notice,

$\boldsymbol{x}^{*}=(\boldsymbol{q}\boldsymbol{q}^{*})^{*}=(\boldsymbol{q}^{*})^{*}\boldsymbol{q}^{*}=\boldsymbol{q}\boldsymbol{q}^{*}=\boldsymbol{x},$
(5)

which implies $\boldsymbol{x}$ has the form

$\boldsymbol{x}=x_{0}+x_{1}i+x_{2}j.$ (6)

To find the values of $\\{x_{i}\\}_{0}^{2}$, we do the quaternion
multiplication. The result yields

$\displaystyle x_{0}$ $\displaystyle=q_{0}^{2}-q_{1}^{2}-q_{2}^{2}+q_{3}^{2}$
$\displaystyle x_{1}$ $\displaystyle=2(q_{0}q_{1}-q_{2}q_{3})$ (7)
$\displaystyle x_{2}$ $\displaystyle=2(q_{0}q_{2}+q_{1}q_{3}).$

This is the KS transformation in its classical form. Thus the KS
transformation $q=(q_{0},q_{1},q_{2},q_{3})\in\mathbb{R}^{4}\mapsto
x=(x_{0},x_{1},x_{2})\in\mathbb{R}^{3}$ is given by the quaternion relation

$\boldsymbol{x}=\boldsymbol{q}\boldsymbol{q}^{*},$ (8)

where $\boldsymbol{x}=x_{0}+x_{1}i+x_{2}j$,
$\boldsymbol{q}=q_{0}+q_{1}i+q_{2}j+q_{3}k$, and $\boldsymbol{q}^{*}$ is
defined as (1).

##  II Differentiation

Since the mapping (4) involving the star conjugate maps from four dimensions
to three, we are left with one degree of freedom. We take advantage of this by
simplifying differentiation. To do this we impose a bilinear relation between
a quaternion $\boldsymbol{q}$ and its differential, $d\boldsymbol{q}$, being

$\boldsymbol{q}\,d\boldsymbol{q}^{*}=d\boldsymbol{q}\,\boldsymbol{q}^{*}.$ (9)

Note (9) has the form of a commutator. The differential of the star conjugate
mapping (4) then becomes

$d\boldsymbol{x}=d\boldsymbol{q}\,\boldsymbol{q}^{*}+\boldsymbol{q}\,d\boldsymbol{q}^{*}=2d\boldsymbol{q}\,\boldsymbol{q}^{*}.$
(10)

Using these techniques, Waldvogel (2008) reduces the spatial Kepler problem to

$2\boldsymbol{u}^{\prime\prime}+h\boldsymbol{u}=0,$ (11)

with the parameter

$dt=r\cdot d\tau,\quad\frac{d}{d\tau}()=()^{\prime},$ (12)

where $t$ is time, $r$ is the radius, and $u=q^{*}$. For later importance, we
note the energy integral:

$-2\left|\boldsymbol{q}^{\prime}\right|^{2}=-\mu+rh.$ (13)

##  III Deriving Kepler’s laws

In this section, we will derive Kepler’s first law in the standard method for
regularizing celestial mechanics, then after which we will prove Kepler’s
second and third laws using new techniques based upon differentiation and
integration.

###  III.1 Kepler’s first law

We first solve (11) to find the general solution,

$\boldsymbol{u}=A\cos(\omega\tau)+B\sin(\omega\tau),\
\omega=\sqrt{\frac{h}{2}},\ A,B\in\mathbb{H}.$ (14)

There are two equations we note before continuing. The first is the $\cos$ and
$\sin$ functions are linear in $\mathbb{H}$, and thus are two dimensional. The
second is since the equation for motion of $\boldsymbol{u}$ is two
dimensional, we may rotate it into the complex plane without loss of
generality, to make solving the fine details a bit easier. So after rotating
to the complex plane, we rotate once more, such that the apocenter is a real
number greater than zero.

Note that we are not scaling, translating, or making any changes to the
equation at all, but merely changing its direction. Once we have solutions to
the equation of motion in this form, we may rotate it back to the original
equation of motion as needed.

After doing our rotations, (14) becomes,

$\boldsymbol{u}=A\cos(\omega\tau)+iB\sin(\omega\tau),\ A,B\in\mathbb{R}.$ (15)

Earlier we made the substitution
$\boldsymbol{x}=\boldsymbol{q}\boldsymbol{q}^{*}=\boldsymbol{u}^{*}\boldsymbol{u}$,
but with our rotations we have $\boldsymbol{u}\in\mathbb{C}$, not
$\boldsymbol{u}\in\mathbb{H}$. Thus for our substitution, we need not include
the star conjugate, so we have $\boldsymbol{x}=\boldsymbol{u}^{2}$. Therefore
our equation for motion of $\boldsymbol{x}$ is,

$\displaystyle\boldsymbol{x}$
$\displaystyle=\left(A\cos(\omega\tau)+iB\sin(\omega\tau)\right)^{2}$
$\displaystyle=A^{2}\cos^{2}(\omega\tau)-B^{2}\sin^{2}(\omega\tau)+2AB\cos(\omega\tau)\sin(\omega\tau)$
$=\frac{A^{2}-B^{2}}{2}+\frac{A^{2}+B^{2}}{2}\cos(2\omega\tau)+iAB\sin(2\omega\tau).$
(16)

Since $2\omega\tau$ defines the angle of the planet from the center, it must
be the eccentric anomaly. Thus the eccentric anomaly, which we denote $E$,
must be,

$E=2\omega\tau=\sqrt{2h}\tau.$ (17)

We substitute (17) into (16) to get,

$\boldsymbol{x}=\frac{A^{2}-B^{2}}{2}+\frac{A^{2}+B^{2}}{2}\cos(E)+iAB\sin(E).$
(18)

From equation (18) we see the center, which we denote $c$, must be,

$c=\frac{A^{2}-B^{2}}{2}.$ (19)

Since we know $(A^{2}+B^{2})/2\geq AB$ for all $A$ and $B$, we have,

$a=\frac{A^{2}+B^{2}}{2},$ (20)

where $a$ is the semi-major axis, and,

$b=AB,$ (21)

where $b$ is the semi-minor axis. We now may prove Kepler’s first law.

######  Theorem 1 (The Law of Orbits).

All planets move in elliptical orbits, with the sun at one focus.

###### Proof.

Note that the equation for the planets motion has the sun placed at the
origin. Also note that the equation for motion is obviously either an ellipse,
or a perfect circle. When using the semi-major axis $a$, the semi-minor axis
$b$, and the center $c$, we find,

$a^{2}-b^{2}-c^{2}=0.$ (22)

This is the definition of a focus of an ellipse, thus the planet orbits the
origin. Also if the orbit is a perfect circle, the center is at zero, and thus
orbits the sun. Therefore since the sun is at the origin, the planet orbits
the sun at one focus. ∎

###  III.2 Kepler’s second law

Since we have derived Kepler’s first law, we know that

$e=\frac{c}{a},$ (23)

where $e$ is the eccentricity of the orbit. We then use (19), (20), and (23)
to solve for $A$ and $B$. We find,

$A=\sqrt{a(1+e)},\qquad B=\sqrt{a(1-e)}.$ (24)

Plugging (24) into (18) yields

$\boldsymbol{x}=a(e+\cos(E))+ia\sqrt{1-e^{2}}\sin(E).$ (25)

Using (25) we calculate the magnitude of $\boldsymbol{x}$ to find,

$r=a(1+e\cos(E)).$ (26)

We are now prepared to derive Kepler’s second law.

######  Theorem 2 (The Law of Areas).

A line that connects a planet to the sun sweeps out equal areas in equal
times.

###### Proof.

Since $\boldsymbol{x}=a(e+\cos(E))+ia\sqrt{1-e^{2}}\sin(E)$, we have,

$\boldsymbol{x}_{0}=a(e+\cos(E)),\qquad\boldsymbol{x}_{1}=a\sqrt{1-e^{2}}\sin(E).$
(27)

We then let,

$\varphi=\tan^{-1}\left(\frac{\boldsymbol{x}_{1}}{\boldsymbol{x}_{0}}\right),$
(28)

and let the area swept out by the line segment be $\Sigma$.

$r$$rd\varphi$$\varphi$ Figure 1: Change in Area

From figure (1), we can define the change of area being swept out as,

$\frac{d\Sigma}{dt}=\frac{1}{2}r^{2}\frac{d\varphi}{dt}=\frac{1}{2}r^{2}\frac{d\varphi}{dE}\frac{dE}{dt}=\frac{1}{2}r^{2}\frac{d\varphi}{dE}\frac{dE}{d\tau}\frac{d\tau}{dt}.$
(29)

We take the derivative of (17) to get,

$\frac{dE}{d\tau}=\sqrt{2h}.$ (30)

We then substitute (30), and (121) into (29), which yields

$\frac{d\Sigma}{dt}=\frac{1}{2}r^{2}\left(\frac{d}{dE}\tan^{-1}\frac{\boldsymbol{x}_{1}}{\boldsymbol{x}_{0}}\right)\sqrt{2h}\frac{1}{r}.$
(31)

We then take the derivative, cancel one $r$, substitute $r=a(1+e\cos(E))$, and
move the constants to the left, which transforms (31) into

$a\sqrt{\frac{h}{2}(1-e^{2})}(1+e\cos(E))\left(\frac{f(E)}{g(E)}\right),$ (32)

where

$f(E)=\frac{\cos(E)}{\cos(E)+e}+\frac{\sin^{2}(E)}{(\cos(E)+e)^{2}},$ (33)

and

$g(E)=1-\frac{(e^{2}-1)\sin^{2}(E)}{(\cos(E)+e)^{2}}.$ (34)

We now reduce the fraction. First we look at the numerator, $f(E)$. We
multiply the top and bottom of the first fraction in the numerator by its
denominator, and combine with the second fraction, yielding,

$\frac{e\cos(E)+1}{(\cos(E)+e)^{2}}.$ (35)

We then multiply the $1$ in $g(E)$ on a common denominator, and combining with
the second fraction to get,

$\frac{(\cos(E)+e)^{2}-(e^{2}-1)\sin^{2}(E)}{(\cos(E)+e)^{2}}.$ (36)

Since the $f(E)$ and $g(E)$ have the same denominator, we cancel it, which
reduces the fraction $f(E)/g(E)$ to

$\frac{e\cos(E)+1}{(\cos(E)+e)^{2}-(e^{2}-1)\sin^{2}(E)}.$ (37)

Note that the numerator is identical to $r/a$. We multiply in our $r$ term to
get,

$\frac{(e\cos(E)+1)^{2}}{(\cos(E)+e)^{2}-(e^{2}-1)\sin^{2}(E)}.$ (38)

We then expand the numerator and denominator to get,

$\frac{e^{2}\cos^{2}(E)+2e\cos(E)+1}{\cos^{2}(E)\\!+\\!2e\cos(E)+e^{2}-e^{2}\sin^{2}(E)\\!+\\!\sin^{2}(E)}.$
(39)

Note we can now reduce the numerator and denominator to be the same value, so

$\frac{e^{2}(1-\sin^{2}(E))+2e\cos(E)+1}{e^{2}(1-\sin^{2}(E))+2e\cos(E)+1}=1.$
(40)

Thus substituting (40) into (32) yields,

$\frac{d\Sigma}{dt}=a\sqrt{\frac{h}{2}(1-e^{2})}.$ (41)

We now have $d\Sigma/dt$ equal to a constant, which proves Kepler’s second
law. ∎

###  III.3 Kepler’s third law

We now work to represent the energy $h$ with variables describing the motion
of the planet. We first note that (13) gives us one equation representing
${2\left|\boldsymbol{q}^{\prime}\right|^{2}=2\left|\boldsymbol{u}^{\prime}\right|^{2}}$.
We then take the derivative of (15), and find that,

$2\left|\boldsymbol{u}^{\prime}\right|^{2}=ah(1-e\cos(E)).$ (42)

We set (42) and the negative of (13) equal, which yields $2ah=\mu$, or,

$h=\frac{\mu}{2a}.$ (43)

We now are prepared to derive Kepler’s third law.

######  Theorem 3 (The Law of Periods).

The square of the period of any planet is proportional to the cube of the
semi-major axis of its orbit.

###### Proof.

Firstly we define,

$n:=\sqrt{\frac{\mu}{a^{3}}}.$ (44)

We square and divide (44) by two to show

$\frac{n^{2}}{2}=\frac{\mu}{2a^{3}}.$ (45)

We substitute (43) into (45) to get

$\frac{n^{2}}{2}=\frac{h}{a^{2}}.$ (46)

We then multiply by two, take the square root, and substitute (30) to get

$n=\frac{1}{a}\frac{dE}{d\tau}=\frac{1}{a}\frac{dE}{dt}\frac{dt}{d\tau}.$ (47)

From (121) we substitute $dt/d\tau$ and rearrange to get

$ndt=\frac{r}{a}dE.$ (48)

We then substitute (26) in and take the definite integral over the period $T$,
in which the eccentric anomaly goes from $0$ to $2\pi$,

$\displaystyle\int_{0}^{T}ndt$ $\displaystyle=$
$\displaystyle\int_{0}^{2\pi}\frac{a(1+e\cos(E))}{a}dE,$ (49) $\displaystyle
nT$ $\displaystyle=$ $\displaystyle 2\pi.$ (50)

Finally we substitute (44) into (50), square both sides, and rearrange to
yield

$\frac{T^{2}}{a^{3}}=\frac{4\pi^{2}}{\mu}.$ (51)

Since we have the ratio of the period squared to the semi-major axis cubed
being a constant, the derivation is complete. ∎

Now that Kepler’s third law is derived, we note that $n$ is actually the mean
motion of the planet, and rearranging (50) shows that,

$n=\frac{2\pi}{T}.$ (52)

We summarize by deriving the area of an ellipse.

######  Theorem 4 (Area of an Ellipse).

The area of the ellipse generated by the orbit of the planet, $A_{e}$, is
found to be

$A_{e}=\pi ab.$ (53)

###### Proof.

We substitute (24) into (21) to find

$b=a\sqrt{1-e^{2}}.$ (54)

We then substitute (54) into (41) to find

$\frac{d\Sigma}{dt}=b\sqrt{\frac{h}{2}}.$ (55)

We then multiply the $dt$ on both sides, and integrate. We want to integrate
over the whole area of the ellipse, $A_{e}$, which takes the period, $T$,
time, so we have

$\begin{split}\int_{0}^{A_{e}}d\Sigma=&\int_{0}^{T}b\sqrt{\frac{h}{2}}dt\\\
A_{e}=&bT\sqrt{\frac{h}{2}}.\end{split}$ (56)

We then substitute (43), and (50) into (56),

$A_{e}=b\frac{2\pi}{n}\sqrt{\frac{\mu}{4a}}.$ (57)

Thus after substituting in (44) and canceling terms, (57) yields

$A_{e}=\pi ab.$ (58)

∎

## References

  * Waldvogel (2006) J. Waldvogel, Celestial Mech. Dynam. Astronom. 95, 201 (2006).
  * Waldvogel (2008) J. Waldvogel, Celestial Mech. Dynam. Astronom. 102, 149 (2008).