Developing Kepler's first law from the two-body problem

Solution 1:

Starting point

Newtons law of gravity gives for a 2-body system \begin{align} \ddot X_1&=-\frac{m_2G(X_1-X_2)}{\|X_1-X_2\|^3} \\ \ddot X_2&=-\frac{m_1G(X_2-X_1)}{\|X_2-X_1\|^3} \end{align}

From 2-body to 1-body around center-of-mass

Note that the right side only depends on the difference of the positions, and that the left side does not change if one adds linear functions to the positions. In total, the system is invariant if one adds the same linear function to both. Now consider that for the linear combination $m_1X_1+m_2X_2$ you get $$ m_1\ddot X_1+m_2\ddot X_2=0 $$ so the center-of-mass moves linearly, this linear motion is computable from the initial values.

To complement this, for the difference of positions one now gets $$ \ddot X_1-\ddot X_2=-\frac{(m_1+m_2)G(X_1-X_2)}{\|X_2-X_1\|^3} $$ Set $m=m_1+m_2$. With $m_1X_1+m_2X_2=m(X_0+tV_0)$ and $D=X_1-X_2$ the trjectories of the two bodies are obtained as solutions of the linear system, so $$ X_1 = (X_0+tV_0) + \frac{m_2}{m}D\\ X_2 = (X_0+tV_0) - \frac{m_1}{m}D\\ $$

Reduction to polar coordinates

Represent $D=X_1-X_2=rS$ in a polar decomposition where $r$ is a positive scalar and $S$ a unit vector. Then per the Leibniz product rule $$ \ddot rS+2\dot r\dot S+r\ddot S=-\frac{mGS}{r^2} $$ A first conclusion is that $\ddot S$ lies in the plane of $S$ and $\dot S$, so that this plane is invariant. Let $N$ be a normalized unit vector orthogonal to that plane. The radius vector of a curve on a sphere is always orthogonal to its tangent vector. Together this means that $$ \dot S=\omega(N\times S)\implies \ddot S =\dot\omega(N\times S)+\omega(N\times \dot S)=\dot\omega(N\times S)-\omega^2 S $$ with some function $\omega$ which is the angular speed of the orbit around the origin.

The equation of motion thus splits into two components $$ \left(\ddot r+\frac{mG}{r^2}-r\omega^2\right) S +(2\dot r\omega+r\dot\omega)(N\times S)=0 $$

Keplers second law

The second component has an integrating factor $r$ to get $$ 2\dot r\omega+r\dot\omega=0\implies r^2\omega=c $$ where $c/2$ is the constant "area velocity".

Keplers first law

The remaining part of the equation of motion can now be formulated in $r$ alone $$ \ddot r+\frac{mG}{r^2}-\frac{c^2}{r^3}=0 $$ Multiply with $2\dot r$ and integrate $$ \dot r^2-\frac{2mG}r+\frac{c^2}{r^2}=const. $$ The last two terms can be completed to a square by adding a constant, so finalize that equation as $$ \dot r^2+\left(\frac{c}{r}-\frac{mG}{c}\right)^2=\rho^2. $$ This can be read as a circle equation of constant radius $\rho$, so parametrize it as such $$ \dot r=\rho\sin\phi,~~~ \frac{c}{r}-\frac{mG}{c}=\rho\cos\phi. $$ At this point, $\phi$ is just some time-dependent function without a direct relation to the orbit position. But take the derivative of the second equation and insert the first to get $$ -\rho\dot \phi\sin\phi=-\frac{c\dot r}{r^2}=-\frac{c\rho\sin\phi}{r^2} \implies \dot\phi=\frac{c}{r^2}=\omega $$ which means that $\phi$ plus some offset is the angle of the orbit.

The second circle component equation can be transformed into the usual radius parametrization $$ r=\frac{c}{\frac{mG}{c}+\rho\cos\phi}=\frac{\frac{c^2}{mG}}{1+\frac{c\rho}{mG}\cos\phi}=\frac{R}{1+e\cos\phi} $$ To see that this is indeed a conic section/quadratic curve in coordinates $x=r\cos\phi$, $y=r\sin\phi$ take the last equation to get $$ R=\sqrt{x^2+y^2}+ex\\ (R-ex)^2=x^2+y^2\\ R^2=(1-e^2)x^2+2eRx+y^2\\ (e\ne1):~~R^2=\left((1-e^2)x+eR\right)^2+(1-e^2)y^2 $$ Depending on the value of the eccentricity $e$ this gives [$e<1$] an ellipse, [$e=1$] a parabola or [$e>1$] a hyperbola.

Now do some clean-up of the ad-hoc constants. Comparing the parametrizations one has $c^2=mGR$, $\rho=e\sqrt{\frac{mG}{R}}$.