Celestial Mechanics I - Introduction to Dynamical Systems

The Kepler Problem

The Newtonian two body problem with masses $m_1$ and $m_2$ can be reduced to the canonical Kepler problem ODE

$$ \ddot{q} = - \frac{q}{|q|^3} $$ where $q = (x_1(t), x_2(t), \ldots, x_d(t))$ could be in an arbitrary dimension $d$ and $r = |q|$ is the 2-norm of $q$ . However, it is easy to show that for any central force problem of the form

$$ \ddot{q} = - q \ f(r) $$ (where $f(r)$ is a scalar function) the solution $q$ remains in a plane, so without loss of generality we can consider the problem in two dimensions, $q = (x(t), y(t))$ .

Once we restrict to a plane, it is sometimes convenient to use a complex coordinate $z = x + i y$ , especially in the polar form $z = r e^{i \theta}$ .

We calculate the first two time derivatives of $z$ :

$\dot{z} = \dot{r}e^{i \theta} + i r e^{i \theta} \dot{\theta}$

$\ddot{z} = \ddot{r}u + i 2 e^{i \theta} \dot{r} \dot{\theta} + i r e^{i \theta} \ddot{\theta} - r e^{i \theta} (\dot{\theta})^2$

From the Kepler ODE equation we now have

$\ddot{r}e^{i \theta} + i 2 e^{i \theta} \dot{r} \dot{\theta} + i r e^{i \theta} \ddot{\theta} - r e^{i \theta} (\dot{\theta})^2 = - e^{i \theta}/r^2$

Dividing by $e^{i \theta}$ we get:

$\ddot{r} + i 2 \dot{r} \dot{\theta} + i r \ddot{\theta} - r \dot{\theta}^2 = -1/r^2$

Separating out the real and imaginary parts we obtain two real ODEs:

$\ddot{r} - r \dot{\theta}^2 = -1/r^2$

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

Now note that only $\dot{\theta}$ and $\ddot{\theta}$ appear in these equations, so it is convenient to introduce a variable $w = \dot{\theta}$ for the angular velocity.

$\ddot{r} - r w^2 = -1/r^2$

$2 \dot{r} w + r \dot{w} = 0$

The second equation above can be immediately integrated (after multiplying by the integrating factor $r$ ) to obtain $$ r^2 w = L $$ where $L$ is a constant, the angular momentum. In $(x,y)$ coordinates

$L = x \dot{y} - y \dot{x}.$

The constancy of $L$ is equivalent to Kepler's second law of planetary motion, which states that the trajectory of $z$ will sweep out equal areas in equal times. It allows us to eliminate $w$ from the first ODE:

$\ddot{r} = L^2/r^3 -1/r^2$

This second-order ODE in one variable is not easy to solve directly; the solution is to usually parameterize $r$ by $\theta$ rather than the time $t$ , and then find an implicit relation between $t$ and $\theta$ .

The relation $r^2 w = L$ can be rearranged to $1/w = \frac{d t}{d \theta} = r^2/L$ , and then thinking of time as a function of $\theta$ we can write the derivatives of $r$ as functions of $\theta$ :

$\frac{d r}{d \theta} = \dot{r} r^2/L$

and

$\frac{d^2 r}{d \theta^2} = \ddot{r} r^4/L^2 + 2 \dot{r}^2 r^3/L^2$

That doesn't look very promising, but with a miraculous change of variables to $u = 1/r$ , we find

$\frac{d u}{d \theta} = - \frac{1}{r^2} \frac{d r}{d \theta} = - \frac{\dot{r}}{L}$

and

$\frac{d^2 u}{d \theta^2} = -\frac{\ddot{r}r^2}{L^2}$

Now we can write the ODE for $u$ in terms of $\theta$ :

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

So remarkably the solution for $u$ is just an offset harmonic motion:

$u = \frac{1}{L^2} + C_1 \cos{(\theta + \phi_0)}$

In terms of the radial variable $r$ , with some assumptions on the constants we get

$r = \frac{p}{1 + e \cos(f)} = \frac{a(1-e^2)}{1 + e \cos(f)}$

at least in the case of an elliptic (rather than parabolic or hyperbolic orbit). The middle expression is the standard form of an ellipse in polar coordinates with the origin at one of the foci. $p$ is called the parameter, or semi-latus rectum, of the ellipse, and it is the distance from the focus to the ellipse in the $y$ -direction. We have also chosen to assume that $\phi_0 = 0$ , rotating the coordinates so that the longest axis of the ellipse is on the $x$ -axis. This axis is called the Line of Apsides in celestial mechanics.

The parameters $a$ and $e$ in the last expression for $r$ are the semi-major axis length ( $a$ ) and the eccentricity ( $e$ ). It is also sometimes convenient with ellipses to use the semi-minor axis length $b$ , and the distance from the center to a focus $c$ . For a proper ellipse the eccentricity is in the interval $[0,1)$ .

If we use coordinates $\tilde{x}$ and $\tilde{y}$ centered on the center of the ellipse, which is probably more familiar, then the implicit formula for the ellipse is

$\frac{\tilde{x}^2}{a^2} + \frac{\tilde{y}^2}{b^2} = 1$

Some useful relationships between these parameters are

$e = \frac{c}{a} = \sqrt{1 - \frac{b^2}{a^2}}$

$c = \sqrt{a^2 + b^2}$

$p = \frac{b^2}{a} = a (1 - e^2)$

The sum of the distances from the foci to each point on the ellipse is equal to $2a$ .

To find the position as a function of time, we need to know $\theta$ in terms of $t$ . It seems impossible to do this explicitly, but Kepler (later simplified by Gauss) found that if we introduce a new angle $E$ , the "eccentric anomaly", then from the two relations

$\tan(\frac{\theta}{2}) = \sqrt{\frac{1 + e}{1 - e}} \tan(\frac{E}{2})$

and Kepler's equation

$\frac{2 \pi}{P} t = E - e \sin(E)$

we can find first $E$ from $t$ and then $\theta$ from $E$ .

There are many other geometric facts about ellipses on the Wikipedia entry for ellipses.

In terms of the Kepler problem, some more relations between the ellipse and the constants $h$ (energy), $L$ (angular momentum), $P$ (the period), and $m = m_1 + m_2$ (total mass of the two-body problem) are:

$\frac{L^2}{m} = a |1-e^2|$

(Here I am including the total mass $m = m_1 + m_2$ because it is not always obvious where it would go.)

$a = \frac{m}{2 |h|}$

$P = \frac{2 \pi}{\sqrt{m a^3}} = \frac{2 \pi a b}{L}$

The extra symmetry of the Kepler problem manifests itself in a variety of ways. One of them is the conserved "Laplace-Runge-Lenz" or LRL vector:

$A = |v|^2 q - (q \cdot v) v - \frac{m q}{r}$

In the plane if $A=(\alpha, \beta)$ then

$e = \frac{|A|}{m}$

$\alpha = L \dot{y} - \frac{m x}{r}$

$\beta = -L \dot{x} - \frac{m y}{r}$

$mr = L^2 - \alpha x - \beta y$

This last relation is another way to see how the shape of the orbit is a conic, since it defines a plane section of the cone $r^2 = x^2 + y^2$ in $(x,y,r)$ space.

The Restricted Three-Body Problem

The full three-body problem in the Newtonian model of gravitation has been intensively studied, and while many interesting and useful results are known the dynamics are extremely rich and no explicit solutions are really possible in general. For many practical purposes it is useful to consider the much simpler case in which one of the masses is much smaller than the other two, and its influence can be safely neglected. This is a good approximation for many objects in the solar system such as human-made satellites and spacecraft, asteroids, comets, small moons, and even small planets.

A further simplication can be made to the circular restricted three-body problem, in which the two massive (primary) bodies are assumed to be on circular (rather than elliptical) orbits. This is a good approximation for the orbit of the Earth-Moon barycenter and the Sun, and for the Sun and Jupiter which make up the vast majority of the mass in our solar system. We can choose coordinates so that the primary masses are in the plane of the first two coordinates.

Some standard choices for the circular restricted three-body problem is to choose units with the gravitational constant $G=1$ , the total mass $m_1 + m_2 = 1$ , and the diameter of the circular orbit to be $1$ . The period of the two massive bodies is then $2 \pi$ . Usually the mass $m_2$ is assumed to be smaller (sometimes much smaller) than $m_1$ , and $\mu = m_2$ is used as the single mass parameter (so that $m_1 = 1 - \mu$ ).

We can also choose to have $t=0$ correspond to the primary masses being on the x-axis, so $q_1(0) = (-\mu, 0, 0)$ and $q_2(0) = (1-\mu,0,0)$ . Then $q_1(t) = -\mu (\cos(t), \sin(t),0 )$ and $q_2(t) = (1-\mu)(\cos(t), \sin(t),0)$ .

Next, to get a time-independent dynamical system we will work in rotating coordinates so that the primaries are not moving, but rather fixed at their initial positions. The new coordinates for the third mass will be

$q = R \ q_3 = \left ( \begin{array}{ccc} \cos(t) & \sin(t) & 0 \\ -\sin(t) & \cos(t) & 0 \\ 0 & 0 & 1 \end{array} \right ) q_3$

Differentiating the above relationship we obtain the ODEs for $q(t) = (x(t), y(t), z(t))$ . The time dependence of the coordinates means that several new types of terms get introduced (sometimes called Coriolis and centripedal "forces"; they are not true forces, just artifacts of the coordinate system).

For the first derivative:

$\dot{q} = \dot{R} q_3 + R \dot{q}_3 = R (R^{-1}\dot{R} q_3 + \dot{q}_3) = R (J q_3 + \dot{q}_3)$

where $J$ is a constant matrix,

$J = \left ( \begin{array}{ccc} 0 & 1 & 0 \\ -1 & 0 & 0 \\ 0 & 0 & 0 \end{array} \right )$

Then the second derivative:

$\ddot{q} = \dot{R} (J q_3 + \dot{q}_3) + R (J \dot{q}_3 + \ddot{q}_3)$

Writing everything in terms of $q$ , $\dot{q}$ , and $\ddot{q}$ we get (after some algebra)

$\ddot{x} = x + 2\dot{y} + \frac{\partial V}{\partial x}$

$\ddot{y} = y - 2\dot{x} + \frac{\partial V}{\partial y}$

$\ddot{z} = \frac{\partial V}{\partial z}$

where $V$ is the potential in the new coordinates

$V = \frac{1-\mu}{\sqrt{(x+\mu)^2 + y^2 + z^2}} + \frac{\mu}{\sqrt{(x-1+\mu)^2+y^2+z^2}}$

We can also write these equations as a classical mechanical system with the energy-like function

$H = \frac{1}{2}|\dot{q}|^2 - V - \frac{1}{2}(x^2 + y^2)$

The fact that $H$ is conserved (i.e. time-invariant) was noticed by Jacobi, and the quantity $-2H$ is sometimes called the Jacobi constant.

Below is a contour plot of the function $V+\frac{1}{2}(x^2 + y^2)$ .

[an error occurred while processing this directive]