On dynamics of the Kepler problem

The two-body central-force problem is also known as the Kepler problem is a classical problem in the domain of classical mechanics. The textbooks pay much attention on the orbits, at least partly because
the orbits have closed-form. In this post, I put emphasis on the dynamic solutions, i.e. the time evolution of the motion, which are the direct solutions of the dynamic equation.

The notations

$r$: the radial coordinate, the distance to the center-of-mass;
$\theta$: the angular coordinate;

$m$: the reduced mass;
$l$: the module of the angular momentum ($l>0$);
$V(r) = - \frac{k m}{r}$: the potential energy. ($k > 0$);
$E$: total energy. ($E < 0$);
$\varepsilon$: the eccentricity. ($-1 < \varepsilon < 1$);
$a$: the semi-major axis;
$r_{\min} = a(1-\varepsilon), r_{\max} = a(1+\varepsilon)$: the periapsis and apoapsis distances.
$T=2\pi\sqrt{\frac{a^3}{k}}$: the period;

Here we have listed two sets of parameters: 1. the physical parameters $m, l, E, k$ (more precisely, $\frac{l}{m}, \frac{E}{m}, k$); 2. the orbit parameter (kinematical parameters) $a, \varepsilon, T$. Each set can determine the dynamic solution. We shall see their relation later.

The equations of motion

[ equation 1 ] $$\ddot r - \frac{l^2}{m^2 r^3} + \frac{1}{m}\frac{\partial }{\partial r} V(r) = 0; \\ \dot \theta = \frac{l}{m r^2}; \\ r(t=0) = r_{\min}, \quad \dot r(t=0) = 0, \quad \theta(t=0) = 0;$$ where $\dot{}$ stands for the time derivative.

The orbits

The orbit is an equation of the coordinates $F(r,\theta) =0$. It is convenient to express $r$ as a function of $\theta$: $r = r(\theta)$. The classical way to get the orbit is first to use the chain rule: $\frac{d}{dt} = \dot{\theta}\frac{d}{d\theta} = \frac{l}{mr^2}\frac{d}{d\theta}$. Then the equation of motion becomes
[ equation 2 ] $$r''_\theta - 2\frac{r^{'2}_\theta}{r} - r + \frac{m}{l^2}r^4\frac{\partial}{\partial r}V(r) = 0.$$ By doing a change of variable $u \equiv \frac{1}{r}$, the orbit equation becomes,
[ equation 3 ] $$u''_\theta + u + \frac{m}{l^2}\frac{\partial}{\partial u}V(1/u) = 0$$
For $V = -\frac{k m}{r}$, it becomes $u''_\theta + u = \frac{km^2}{l^2}$. And the solution is,
[ solution 1 ] $$u = \frac{1}{r} = \frac{km^2}{l^2}(1+\varepsilon \cos(\theta))$$ The orbit is an ellipse. The semi-major axis $a (1-\varepsilon^2) = \frac{l^2}{km^2}$.

The dynamic solutions

To solve for the dynamics $\vec{r}(t)$ (in polar coordinate $\vec{r}(t) = (r(t), \theta(t))$), we shall return to equation 1. Note that equation 1 does not depend on $t$ explicitly. This type of 2nd order differential equation can always be solved by using the chain rule: $\frac{d}{dt} = \frac{dr}{dt}\frac{d}{dr}$. Equation 1 becomes:
[ equation 4 ] $$\dot r \frac{d \dot r}{dr} + \frac{d}{dr} \left( \frac{1}{m}V(r) + \frac{l^2}{2m^2 r^2} \right) = 0 \\ \implies \frac{1}{2} m \dot r^2 + V(r) + \frac{l^2}{2mr^2} = E$$ Substitute the orbit equation, we get two relations $a(1-\varepsilon^2) = \frac{l^2}{km^2}, E = -\frac{km}{2a}$ thus $$a = -\frac{km}{2E}, \quad \varepsilon = \sqrt{1 + 2\frac{l^2E}{k^2m^3}}, \quad T = 2\pi \sqrt{-\frac{k^2m^3}{8E^3}},$$ where the expression about $T$ comes from Kepler's third law. We'll derive it later. The second one in equation 4 is nothing but the energy conservation. To solve equation 4, we note $\dot r = \frac{dr}{dt}$. Thus, $$t = \int_C \frac{dr}{\pm\sqrt{\frac{2}{m} \left(E- \frac{l^2}{2m r^2}-V(r) \right)}}$$ where $C$ is an integral path. The path $C$ also encodes the initial condition. In our case, $r(t=0) = a(1-\varepsilon)$. For $V = -\frac{km}{r}$, $E - \frac{l^2}{2m r^2} - V(r)$ in the denominator of equation 4 is a rational polynomial and can be rewritten as $E (r-r_+)(r-r_-)/r^2$, where $r_\pm$ are the roots of the polynomial. On the other hand, $E - \frac{l^2}{2m r^2_\pm} - V(r_\pm) = 0 \Rightarrow \dot r = 0$. So $r_\pm$ are the two extrema, $r_\pm = a(1\pm \varepsilon)$. Equation 5 can be written as
[ equation 5 ] $$t = \sqrt{-\frac{m}{2E}} \int_C \frac{r dr}{\pm\sqrt{(r- r_-)(r_+-r)}}$$ Note that $\sqrt{-\frac{m}{2E}} = \sqrt{\frac{a}{k}}$.

Fig. The plot of the imaginary part of the integrand with $r_- = 1, r_+ = 2$. $r_\pm$ are the two brach points of the integrand.

To evaluate this integral, we need proper parametrization $r = r(\xi)$. There are several viable parametrization schemes.

scheme 1. If we integrate from $r_{\min}$ to $r< r_{\max}$, then we can evaluate the integral directly (numerically for example):
[ equation 6-1 ] $$t = \sqrt{\frac{a}{k}} \int_{a(1-\varepsilon)}^r \frac{\rho d\rho}{\sqrt{(\rho-a+a\varepsilon)(a+a\varepsilon-\rho)}}$$

scheme 2. $r = a(1-\varepsilon)\cos^2\varphi + a(1+\varepsilon)\sin^2\varphi = a(1-\varepsilon\cos2\varphi), \quad (0 \le \varphi \le \frac{\pi}{2})$. So the integral becomes,
[ equation 6-2 ] $$t = \sqrt{\frac{a}{k}}\int_0^{\frac{1}{2}\arccos\frac{a-r}{a\varepsilon}} 2a(1-\varepsilon\cos2\varphi) d \varphi \\ = \sqrt{\frac{a^3}{k}} \left[ \arccos\frac{a-r}{a\varepsilon}-\sqrt{\varepsilon^2-\left(\frac{a-r}{a}\right)^2}\right],$$ were $\arccos$ is defined on $[0, \pi]$. $t = \frac{1}{2}T$ when $r = a(1+\varepsilon)$. Then $T = 2\pi \sqrt{\frac{a^3}{k}}$, as we know from Kepler's third law. In order to get the expression $r = r(t)$, we need to invert the function, especially $f_\varepsilon(z) \equiv \arccos(z) - \varepsilon\sqrt{1-z^2}$. $f_\varepsilon(z)$ is monotone and its inverse function exists. Then $r(t) = f^{-1}_\varepsilon\left(2\pi\frac{t}{T}\right)$. Unfortunately, $f^{-1}_\varepsilon$ can not be expressed by elementary function. In numerical calculation, however, this is not as bad as one may think. The inverse function can be computed via Newton's algorithm.

scheme 3. We can also parametrize $r$ with $\theta$ directly, $r = a(1-\varepsilon^2)/(1+\varepsilon\cos\theta)$.
[ equation 6-3 ] $$t = \sqrt{\frac{a^3(1-\varepsilon^2)^3}{k}} \int_0^\theta \frac{\mathrm d \phi}{(1+\varepsilon\cos\phi)^2}.$$ We can do a change of variable $s = \tan\frac{\phi}{2}$. Then $d\phi = \frac{2}{1+s^2}ds, \cos\phi = \frac{1-s^2}{1+s^2}$ The integral part becomes, $$\int  \frac{\mathrm d \phi}{(1+\varepsilon\cos\phi)^2} = \int ds \frac{2(1+s^2)}{\left(1+s^2+\varepsilon(1-s^2)\right)^2} = \frac{1}{(1-\varepsilon)^2}\int \frac{ds}{s^2+\frac{1+\varepsilon}{1-\varepsilon}} - \frac{2\varepsilon}{(1-\varepsilon)^3}\int \frac{ds}{\left(s^2+\frac{1+\varepsilon}{1-\varepsilon}\right)^2}$$ Note that $$\int \frac{ds}{s^2+a^2} = \frac{1}{a}\arctan \frac{s}{a} \\ \int \frac{ds}{(s^2+a^2)^2} = \frac{s}{2a^2(s^2+a^2)} + \frac{1}{2a^3}\arctan\frac{s}{a}$$
Therefore,
$$t = \sqrt{\frac{a^3}{k}} \left( 2 \arctan \left(\sqrt{\frac{1-\varepsilon}{1+\varepsilon}} \tan\frac{\theta}{2}\right) - \sqrt{1-\varepsilon^2} \frac{\varepsilon\sin\theta}{1+\varepsilon\cos\theta} \right)$$ (note that the $\arctan$ function returns values on $[0, 2\pi]$. Sometimes, the so-called two-argument inverse tangent function $\arctan(x,y)$ is used.)

Again, in this case, we need to invert the function and obtain an expression for $\theta$.

With either of the above method, we can solve for $\vec{r}(t) = (r(t), \theta(t))$. Finally, we can put the orbit and dynamic solution together, in an animation.

Applications

1. the quasi-satellite

A quasi-satellite of a planet is an object that orbiting the same center and with the same period as the planet (co-orbital configuration). Usually, the quasi-satellite has a different eccentricity. The relative motion of the object with respect to the planet form a closed orbit. So, the object looks as if it is orbiting around the planet with irregular orbits (quasi-orbit). Historically, quasi-satellites resulted several false claims of the existence of other moons of the earth.

Fig. 1a, the relative orbit (yellow), the orbit of the planet (blue), the orbit of the object (purple).

Fig. 2, Blue: the relative orbit (yellow), the orbit of the planet (blue), the orbit of the object (purple).

cf. The quantum Hamilton-Jacobi theory of a central-force problem