[NOTES/CM-05008] Keplar Problem --- Solving Differential Equation

Differential equation of the orbit 

Here we obtain a second order differential equation and solve it for the motion of the planets around the Sun. The energy conservation for a a particle in a spherically symmetric potential takes the form
\begin{equation}\label{EQ21}
E = \frac{1}{2}\mu \Big(\dd[r]{t}\Big)^2 + V(r) +\frac{L^2}{2\mu r^2}
\end{equation}
where \(\mu \) is the mass and \(L\) is the angular momentum of the particle. Differentiating \eqref{EQ21} w.r.t. time we get the EOM for radial coordinate \(r\): 
\begin{equation}\label{EQ22}
\mu \DD[r]{t} = - \dd[V(r)]{r} + \frac{L^2}{\mu r^3}.
\end{equation}
The conservation of angular momentum in plane polar coordinate is of the form:
\begin{equation} \label{EQ04}
\mu r^{2} \dot{\phi} = \text{constant, say} L \end{equation}
Let $f(r) = - \frac{\partial V}{ \partial r}$ be the force law, then \eqref{EQ22} and \eqref{EQ04} give
%
\begin{equation}\label{EQ05}
\mu\, \ddot{r} = f(r) + \frac{L^{2}}{\mu r^{3}}
\end{equation}
Next use \eqref{EQ04} to convert $\ddot{r}$ into $ \frac{d^{2} r}{d \phi^{2}}$
\begin{eqnarray}
\frac{dr}{dt}
&=& \frac{dr}{d \phi} \left( \frac{d \phi}{dt} \right)\label{EQ06}\\
&=& \frac{L}{\mu r^{2}} \frac{dr}{d \phi}\label{EQ06A}
\end{eqnarray}
and
\begin{equation}\label{EQ07}
\frac{d^{2} r}{dt^{2}} = \frac{L}{\mu r^{2}} \frac{d}{d \phi} \frac{L}{\mu
r^{2}}
\frac{dr}{d \phi}= \frac{L^2}{\mu^2r^2}\frac{d}{d \phi} \frac{1}{r^{2}}
\frac{dr}{d \phi}
\end{equation}
Changing variables from \(r\) to $u = \frac{1}{r}$, we get
\begin{eqnarray}\label{EQ09A}
\frac{du}{d \phi} &=& - \frac{1}{r^{2}} \frac{dr}{d \phi}
\end{eqnarray}
and
\begin{eqnarray}\label{EQ09}
\frac{d^{2} u}{d \phi^{2}} &=& - \frac{d}{d \phi} \frac{1}{r^{2}} \frac{dr}{d
\phi^{2}}.
\end{eqnarray}
Using \eqref{EQ07} and \eqref{EQ09} in \eqref{EQ05}, we get
\begin{eqnarray}
\frac{L^2}{ \mu r^{2}} \frac{d}{d \phi} \frac{1}{r^{2}} \frac{dr}{d \phi}
&=&
f(r) + \frac{L^{2}}{ \mu r^{3}}\\
\therefore \qquad - \frac{L^{2}}{r^{2}} \frac{1}{\mu} \frac{d^{2} u}{d
\phi^{2}} &=& f \left(
\frac{1}{u} \right) + \frac{L^{2}}{ \mu} u^{3}
\end{eqnarray}
\begin{equation}\label{EQ12}
\frac{L^{2} u^{2}}{ \mu} \frac{d^{2}u}{d \phi^{2}} = - f \left( \frac{1}{u}
\right) - \frac {L^{2}}{\mu} u^{3}
\end{equation}
This is the required differential equation of the orbit. If the equation of the orbit is known the force law can be found.

Solution of differential equation
We shall now solve the differential \eqref{EQ12} for the attractive inverse square law of force.

\begin{equation}
f = - \frac{k}{r^{2}} = - k u^{2}
\end{equation}
The differential equation of the orbit, \eqref{EQ12} becomes
\begin{equation}
\frac{L^{2}}{\mu} \frac{d^{2}u}{d \phi^{2}} = k - \frac{L^{2}}{\mu} u.
\end{equation}
or
\begin{equation}
\frac{d^{2}u}{d \phi^{2}} + u - \frac{\mu k}{L^{2}}=0
\end{equation}
We now change the variable to $w = u - \frac{mk}{L^{2}}$ to get
\begin{equation}
\frac{d^{2}w}{d \phi^{2}} + w = 0
\end{equation}
This equation can be solved immediately to give
\begin{equation}
w = b \cos(\phi - \phi_{0})
\end{equation}
\begin{equation}\label{EQ19}
u = \frac{\mu k }{L^{2}} + b \cos(\phi - \phi_{0})
\end{equation}
Hence the equation of the orbit can be written in the form
\begin{equation}
\frac{1}{r} = \frac{\mu k}{L^{2}} ( 1 + \epsilon \cos(\phi - \phi_{0}))
\end{equation}

Orbit parameters
Recall that equation of the orbit is given by
\begin{equation}
\frac{\ell_0}{r} = \frac{\mu k}{L^{2}} 1 + \epsilon \cos(\phi - \phi_{0})
\end{equation}
where \(\ell_0=\frac{L^{2}}{\mu k}\), and \(L\) is angular momentum, $\epsilon$, and $ \phi_{0}$ are constants of integration. We shall now relate maximum and minimum values of \(r\), \(r_1,r_2,\), with the eccentricity \(\epsilon\) and energy \(E\).

If $ r_{1}$ and $ r_{2}$ are the minimum and maximum values of \(r\), then from \eqref{EQ19}, we have
\begin{equation}\label{EQ20}
\frac{1}{r_{1,2}} = \frac{\mu k}{L^{2}} ( 1 \pm \epsilon)
\end{equation}

We relate these values \(r_{1,2}\) in terms to \(E\). At the maximum of \(r\), and also at minimum, $ \dot{r} = 0$. The energy conservation equation,
\begin{equation}
E = \frac{1}{2} \mu \dot{r}^{2} + \frac{L^{2}}{2 \mu r^{2}} + V(r),
\end{equation}
with \(\dot{r}=0\),becomes
\begin{equation}
E = \frac{L^{2}}{2 \mu r^{2}} - \frac{k}{r} \qquad \qquad \ \ \ \ \ (\because \dot{r}
= 0),
\end{equation}
or
\begin{equation}
\frac{1}{r^{2}} - \frac{2 \mu k}{L^{2} r} - \frac{2 \mu E}{L^{2}} = 0.
\end{equation}
The values $ r_{1}$ and $r_{2}$ are roots of this equation, hence
\begin{eqnarray}
\frac{1}{r_{1}} + \frac{1}{r_{2}} &=& \frac{2 \mu u}{L^{2}}\\
\frac{1}{r_{1}} \frac{1}{r_{2}} &=& \frac{- 2 \mu E}{L^{2}}\label{EQ25}
\end{eqnarray}
\eqref{EQ20} and \eqref{EQ25} give us the desired expressions for $\epsilon$ in terms of energy
\begin{equation}
\frac{\mu^{2} k^{2}}{L^{4}} ( 1 - \epsilon^{2}) = \frac{-2 \mu E}{L^{2}}
\end{equation}
or
\begin{eqnarray}
\epsilon^{2} - 1 &=&\frac{2 E L^{2}}{\mu k^{2}}\\
\epsilon &=& \sqrt{1+\frac{2 E L^{2}}{\mu k^{2}}}
\end{eqnarray}
Therefore, the final form of the equation of the orbit is
\begin{equation}
\frac{1}{r} = \frac{\mu k }{L^{2}} ( 1 + \epsilon \cos(\phi - \phi_{0}))
\end{equation}
$\phi_{0}$ is the value of $\phi$ at the turning point $r = r_{min}$. This equation represents a conic section with eccentricity = $\epsilon$. Types of orbits traversed, corresponding to different values of \(\epsilon\), are summarised below.

Kepler's Laws

1. The motion in $\frac{1}{r}$ gravitational fields is in elliptic orbits.
2. Areal velocity is constant follows from the conservation laws of angular momentum. \[\mu r^{2} \dot{\phi} = \text{constant, say} L \Longrightarrow r^{2} \dot{\phi} = \frac{L}{\mu} = \text{constant}\] \begin{eqnarray} \text{Areal velocity} &=& \lim_{\Delta t\to 0}\frac{\text{area swept in time} \Delta t} {\Delta t} \\ &=& \lim_{\Delta t\to 0} \frac{1}{2} r^{2} \frac{\Delta \phi}{\Delta t} = \frac{1}{2} r^{2} \dot{\phi} &=& \frac{L}{2 \mu} = \text{constant} \end{eqnarray} \FigBelow{30,10}{60}{90}{arealvelocity}{}
3. To calculate the time period and verify Kepler's third law, since areal velocity is constant

\begin{eqnarray}\label{EQ34} \text{Time Period } T &=&\frac{\text{ Total area of the orbit}}{\text{areal velocity}}\\ &=& \frac{\pi a b }{(\frac{L}{2\mu})} = \frac{2 \mu}{ L} ( \pi a b ) \label{EQ35} \end{eqnarray}

Now semi major axis \begin{equation}\label{EQ36} 2a = r_{1} + r_{2} = \frac{k}{|E|} \end{equation} where (24) and (25) has been used. Remember \(E\) is negative for elliptic orbits, so $E = - |E|$. Now \begin{equation}\label{EQ37} b = a \sqrt{1- \epsilon^{2}} = a \left( \frac{2 |E| L^{2}}{\mu k^{2}} \right)^{\frac{1}{2}}, \end{equation} Using \eqref{EQ35}-\eqref{EQ37} we get \begin{equation} T = \pi a b \frac{2 \mu}{L} = \frac{2 \pi a \mu}{L} a \left( \frac{2 |E| L^{2}}{\mu k^{2}} \right)^{\frac{1}{2}} \end{equation} use $ |E| = \frac{k}{2a}$, from \eqref{EQ36}, and eliminate \(E\) \begin{eqnarray} T &=& \frac{2 \pi a^{2} \mu}{L} \left( \frac{2L^{2}}{\mu k^{2}} \frac{k}{2a} \right) ^{\frac{1}{2}} = \frac{2 \pi a^{2} \mu}{L} \left( \frac{L^{2}}{k \mu a} \right)^{\frac{1}{2}},\\ &=& 2 \pi a^{\frac{3}{2}} \sqrt{\frac{\mu}{k}}, \qquad \end{eqnarray} Now use \(k=GMm\) and \(\mu= \frac{mM}{M+m}\) to get \begin{equation} \boxed{T^2 = {\frac{4 \pi^2a^3}{G(M+m)}}.} \end{equation} In this expression mass \(m\) of the planet can be neglected compared to the mass of the sun \(M\). This proves Kepler's law that square of time period is proportional to the cube of semi major axis. Note that the time period has a tiny dependence on the mass of planet which has been neglected.

Hyperbolic Orbits
We will now derive equation of orbit when \(E>0\) and the orbits are hyperbolic. The distance of closest approach corresponds to \(r =\) minimum and $\frac{1}{r}$ will be a maximum. This happens at $\cos(\phi - \phi_{0}) = 1$ 
$\frac{1}{r} = a (1 + \epsilon \cos( \phi - \phi_{0}))$
$r=r_{min} \Longrightarrow \phi = \phi_{0}$

As $ r\longrightarrow \infty, (note \ \epsilon > 0) \frac{1}{r} \longrightarrow 0 \ and \cos(\phi - \phi_{0}) \longrightarrow \frac{-1}{\epsilon}$
\begin{eqnarray}
\phi - \phi_{0} &=& \pi \pm \cos^{-1} \left( \frac{1}{\epsilon} \right)\\
\phi_{\pm} &=& \phi_{0} + \pi \pm \cos^{-1} \left( \frac{1}{\epsilon} \right).
\end{eqnarray}
The total deflection of the particle in hyperbolic orbit is
\begin{equation}
\theta = \phi_{+} - \phi_{-}=2\cos^{-1} \left( \frac{1}{\epsilon} \right)
\end{equation}
or
\begin{eqnarray}
\cos \left( \frac{\theta}{2} \right) &=& \frac{1}{\epsilon} = \left( 1 +
\frac{2E L^{2}}{\mu k^{2}} \right)^{-1} \\
\cos \left( \frac{\theta}{2} \right) &=& \left( 1 + \frac{2E L^{2}}{\mu k^{2}}
\right)^{-1} \\
&=& \frac{\mu k^{2}}{(\mu k^{2} + 2EL^{2})}
\end{eqnarray}