We will consider the motion of two stars bounded by their mutual gravitational attraction. The unbounded case corresponding to parabolic and hyperbolic motion will not be considered as they are not relevant for the discussion which is about binary systems of stars. We will describe the relative motion of one star with respect to the other, and also the orbital motion of each star relative to the center of mass of the binary system. The stars interact gravitationally only. If the separation between the stars is much larger than the sum of their radii, which is equivalent to having large periods (more than one hundred days) according to the third Kepler’s law, we can consider them as point masses. We will study this simpler case and later we will address the problem of close binary systems.
Equations of motion
Let’s consider two stars \(s_1\) and \(s_2\), with respective masses \(m_1\) and \(m_2\). In a inertial reference system, their corresponding vector position are \(\vec{r}_1\) and \(\vec{r}_2\). The equations of motion in this reference system are
\begin{equation}
m_1\ddot{\vec{r}}_1=\frac{Gm_1m_2}{|\vec{r}_2-\vec{r}_1|^2}\frac{\vec{r}_2-\vec{r}_1}{|\vec{r}_2-\vec{r}_1|}
\end{equation}
\begin{equation}
m_2\ddot{\vec{r}}_2=-\frac{Gm_1m_2}{|\vec{r}_2-\vec{r}_1|^2}\frac{\vec{r}_2-\vec{r}_1}{|\vec{r}_2-\vec{r}_1|^3}\ .
\end{equation}
Equation for the relative motion
We define
\begin{equation}
\vec{r}=\vec{r}_2-\vec{r}_1\quad ;\quad |\vec{r}|=r\ .
\end{equation}
we can rewrite Eqs.(1,2) in the form
\begin{equation}
\ddot{\vec{r}}_1=G\frac{m_2}{r^2}\hat{r}
\end{equation}
\begin{equation}
\ddot{\vec{r}}_2=-G\frac{m_1}{r^2}\hat{r}\ ,
\end{equation}
subtracting we get
\begin{equation}
\ddot{\vec{r}}=-G\frac{(m_1+m_2)}{r^2}\hat{r}
\end{equation}
This equation describes the evolution of the position of one star relative to the other.
or
\begin{equation}
\ddot{\vec{r}}=-\frac{\alpha_{rel}}{r^2}\hat{r}
\end{equation}
where
\begin{equation}
\alpha_{rel}=G(m_1+m_2)
\end{equation}
If we multiply both sides of Eq.(6) by \(m_1m_2\) we can rewrite it in the form
\begin{equation}
\mu\ddot{\vec{r}}=-G\frac{(m_1+m_2)\mu}{r^2}\hat{r}\ .
\end{equation}
This equation implies that the orbital motion of one star from the point of view of the other star is
or
\begin{equation}
\mu\ddot{\vec{r}}=-G\frac{m_1m_2}{r^2}\hat{r}
\end{equation}
where
\begin{equation}
\mu=\frac{m_1m_2}{m_1+m_2}
\end{equation}
is the reduced mass.
The meaning of Eq.(8) is that from the point of view of \(m_1\), the mass \(m_2\) is attracted toward \(m_1\) wit the Newtonian force, but in its motion it behaves as if its mass were the reduced mass.
Equations of motion relative to the center of mass
The center of mass is located at
\begin{equation}
\vec{r}_{cm}=\frac{m_1\vec{r}_1+m_2\vec{r}_2}{m_1+m_2}
\end{equation}
the positions of the masses with respect to the center of mass are
\begin{align}
\vec{r}_1^{\ \prime}&=\vec{r}_1-\vec{r}_{cm}\\
\vec{r}_2^{\ \prime}&=\vec{r}_2-\vec{r}_{cm}
\end{align}
or
\begin{align}
\vec{r}_1^{\ \prime}&=-\frac{m_2\vec{r}}{m_1+m_2}\\
\vec{r}_2^{\ \prime}&=\frac{m_1\vec{r}}{m_1+m_2}
\end{align}
or
\begin{align}
\vec{r}&=-\frac{(m_1+m_2)\vec{r}_1^{\ \prime}}{m_2}\\
\vec{r}&=\frac{(m_1+m_2)\vec{r}_2^{\ \prime}}{m_1}
\end{align}
Substituting in Eq.(7) we get
\begin{align}
\ddot{\vec{r}}_1^{\ \prime}&=-\frac{\alpha_1^{\prime}}{r_1^{\prime\ 2}}\hat{r}_1^{\prime}\\
\ddot{\vec{r}}_2^{\ \prime}&=-\frac{\alpha_2^{\prime}}{r_2^{\prime\ 2}}\hat{r}_2^{\prime}
\end{align}
where
\begin{align}
\alpha_1^{\prime}&=\frac{Gm_2^3}{(m_1+m_2)^2}\\
\alpha_2^{\prime}&=\frac{Gm_1^3}{(m_1+m_2)^2}
\end{align}
Conservation of angular momentum
We take the vector product on both sides
\begin{equation}
\vec{r}\times\mu\ddot{\vec{r}}=-\frac{\alpha\mu}{r^2}\vec{r}\times\hat{r}
\end{equation}
the right side is zero because \(\vec{r}\) and \(\hat{r}\) are parallel. The left side can be written in a different way
\begin{align}
\frac{d}{dt}\left(\vec{r}\times\mu\dot{\vec{r}}\right)&=\dot{\vec{r}}\times\mu\dot{\vec{r}}+\vec{r}\times\mu\ddot{\vec{r}}\\
&=\vec{r}\times\mu\ddot{\vec{r}}
\end{align}
The angular momentum \(\vec{l}\) is
\begin{equation}
\vec{l}=\vec{r}\times\mu\dot{\vec{r}}
\end{equation}
we see that
\begin{equation}
\frac{d\vec{l}}{dt}=0\ .
\end{equation}
the magnitude \(l=|\vec{l}|\) of the angular momentum is
\begin{equation}
l=\mu r^2\dot{\theta}
\end{equation}
Geometrically we have
\begin{equation}
\dot{\vec{A}}=\frac{1}{2}\vec{r}\times\dot{\vec{r}}
\end{equation}
and therefore
\begin{equation}
\dot{\vec{A}}=\frac{\vec{l}}{2\mu}
\end{equation}
Consequences of the conservation of the angular momentum
The angular momentum \(\vec{l}\) is a constant vector. This implies that the vector \(\vec{r}\) is always in the plane whose normal is parallel to \(\vec{l}\), or in other words, the two masses move on a fixed plane up to a parallel translation of it.
Energy conservation
We take the dot product on both sides
\begin{equation}
\dot{\vec{r}}\cdot\mu\ddot{\vec{r}}=-\frac{\alpha\mu}{r^2}\dot{\vec{r}}\cdot\hat{r}
\end{equation}
remember that
\begin{equation}
\dot{\vec{r}}=\dot{r}\hat{r}+r\dot{\theta}\hat{\theta}
\end{equation}
and therefore
\begin{equation}
\dot{\vec{r}}\cdot\hat{r}=\dot{r}\ .
\end{equation}
the right side is
\begin{equation}
-\frac{\alpha\mu}{r^2}\dot{r}=\frac{d}{dt}\left(\frac{\alpha\mu}{r}\right)
\end{equation}
on the left side use the identity
\begin{align}
\frac{d}{dt}\left(\dot{\vec{r}}\cdot\mu\dot{\vec{r}}\right)&=\ddot{\vec{r}}\cdot\mu\dot{\vec{r}}+\dot{\vec{r}}\cdot\mu\ddot{\vec{r}}\\
&=2\mu\dot{\vec{r}}\cdot\dot{\vec{r}}\\
&=2\mu|\dot{\vec{r}}|^2
\end{align}
The equation can be written in the form
\begin{equation}
\frac{d}{dt}\left(\frac{1}{2}\mu|\dot{\vec{r}}|^2-\frac{\alpha\mu}{r}\right)=0
\end{equation}
The quantity between parenthesis is the sum of the kinetic and potential energies of the system (in the rest frame of \(m_1\)). This is the total energy $E$:
\begin{equation}
E=\frac{1}{2}\mu|\dot{\vec{r}}|^2-\frac{\alpha\mu}{r}
\end{equation}
Effective potential
The total energy can be written in the form
\begin{equation}
E=\frac{1}{2}\mu(\dot{r}^2+r^2\dot{\theta}^2)-\frac{\alpha\mu}{r}
\end{equation}
From the conservation of the angular momentum, we have \(\dot{\theta}=l/(\mu r^2)\), so that the total energy can be rewritten as
\begin{equation}
E=\frac{1}{2}\mu\dot{r}^2+\frac{l^2}{2\mu r^2}-\frac{\alpha\mu}{r}
\end{equation}
The effective potential is defined as
\begin{equation}
U_{eff}=\frac{l^2}{2\mu r^2}-\frac{\alpha\mu}{r}
\end{equation}
Integration of the equations of motion
The equation defining the total energy \(E\) is actually a differential equation for the radial distance \(r\) as a function of time.
\begin{align}
\dot{r}&=\dot{\theta}\frac{dr}{d\theta}\\
&=\frac{l}{\mu r^2}\frac{dr}{d\theta}\\
&=-\frac{l}{\mu}\frac{d}{d\theta}\left(\frac{1}{r}\right)
\end{align}
we can rewrite the equation for the total energy in the form
\begin{equation}
E=\frac{l^2}{2\mu}\left[\frac{d}{d\theta}\left(\frac{1}{r}\right)\right]^2+\frac{l^2}{2\mu r^2}-\frac{\alpha\mu}{r}
\end{equation}
we introduce the new variable \(\xi=1/r\) and the equation for the energy becomes
\begin{equation}
\frac{2\mu E}{l^2}=\left(\frac{d\xi}{d\theta}\right)^2+\xi^2-\left(\frac{2\mu^2\alpha}{l^2}\right)\xi
\end{equation}
or
\begin{equation}
\frac{2\mu E}{l^2}=\left(\frac{d\xi}{d\theta}\right)^2+\left(\xi-\frac{\mu^2\alpha}{l^2}\right)^2-\left(\frac{\mu^2\alpha}{l^2}\right)^2
\end{equation}
that can be rearranged in the form
\begin{equation}
\left(\frac{d\xi}{d\theta}\right)^2=\frac{2\mu E}{l^2}+\left(\frac{\mu^2\alpha}{l^2}\right)^2-\left(\xi-\frac{\mu^2\alpha}{l^2}\right)^2
\end{equation}
By inspection, one can check that
\begin{equation}
\xi-\frac{\mu^2\alpha}{l^2}=\sqrt{\frac{2\mu E}{l^2}+\left(\frac{\mu^2\alpha}{l^2}\right)^2}\ \cos(\theta-\theta_0)
\end{equation}
is the solution. Returning back to the variable \(r=1/\xi\), the solution becomes
\begin{equation}
r=\frac{1}{\frac{\mu^2\alpha}{l^2}+\sqrt{\frac{2\mu E}{l^2}+\left(\frac{\mu^2\alpha}{l^2}\right)^2}\ \cos(\theta-\theta_0)}\ .
\end{equation}
It is useful to write the solution in the form
\begin{equation}
r=\frac{l^2/(\mu^2\alpha)}{1+\sqrt{1+\frac{2l^2E}{\mu^3\alpha^2}}\ \cos(\theta-\theta_0)}\ .
\end{equation}
Types of orbits
We recognize from the solution the excentricity $\varepsilon$:
\begin{equation}
\varepsilon=\sqrt{1+\frac{2l^2E}{\mu^3\alpha^2}}
\end{equation}
The type of orbit is determined by \(\varepsilon\).
Circular orbits
The simplest of all is the circular orbit obtained when $\varepsilon=0$. In this case, the solution takes the simple form
\begin{equation}
r=l^2/(\mu\alpha)
\end{equation}
At this radius, the effective potential reaches its minimum:
\begin{align}
\frac{dU_{eff}}{dr}&=-\frac{l^2}{\mu r^3}+\frac{\alpha}{r^2}\\
\frac{dU_{eff}}{dr}=0&\Longrightarrow r=l^2/(\mu\alpha)\ .
\end{align}
The angular velocity $\dot{\theta}$ can be calculated from angular momentum conservation Eq.(19) together with Eq.(32):
\begin{align}
\dot{\theta}&=\frac{\sqrt{\mu\alpha r}}{\mu r^2}\\
&=\sqrt{\frac{\alpha}{\mu r^3}}\ .
\end{align}
The orbital velocity is just
\begin{align}
v&=r\dot{\theta}\\
&=\sqrt{\frac{\alpha}{\mu r}}\ .
\end{align}
Integrating Eq.(34) for one revolution gives the relation between the radius and the period $T$ of the orbit:
\begin{equation}
2\pi=\sqrt{\frac{\alpha}{\mu r^3}}\ T
\end{equation}
which is equivalent to the third Kepler’s law:
\begin{equation}
T^2=\left(\frac{4\pi^2\mu}{\alpha}\right)r^3
\end{equation}
The total energy in this case is
\begin{align}
E&=-\frac{\mu\alpha^2}{2l^2}\\
&=-\frac{\alpha}{2r}\ .
\end{align}
Elliptic orbit
Elliptic orbits are obtained when $0<\varepsilon<1$.