How to find orbital elements of a binary system

Question

I'm trying to find the rate of precession of a binary system consists of two similar mass stars A and B, so I need to find the value of the argument of periapsis $\omega$ first, with velocity and position of both stars are given.

I have done the calculation of the value of $\omega$ for Mercury and the method used can be seen here How to programmatically calculate orbital elements using position/velocity vectors?, or as the person in that post suggested, you can refer to the book Fundamentals of Astrodynamics and Applications, by Vallado, 2007.

However, for the Sun-Mercury case, Sun is much heavier and basically Mercury just rotates around and the Sun almost stays still, so it makes sense to just put the Sun(which is also the center of mass for this system) as the origin and the computation just follows, but for the A-B binary system which has similar mass and the center of mass is approximately half way between A and B, what should I put as the origin for my computation(and what is the mass of this origin that A and B rotate around)?

Update: So now I'm using center of mass as my origin and use the effective mass $\frac{m_{B}^{3}}{(m_{A}+m_{B})^{2}}$ as the mass. I plotted the $\omega$(of star A) vs. $t$ but it doesn't make sense to me because the $\omega$ is not zero which it should be since there are no other planets to perturb and I didn't include GR term, any suggestions?

(The mass of the two stars mA = 0.29097 Solar mass, mB = 0.356 Solar mass. This system is PSR J0737-3039 if details needed.)

Answer 1

Make the stationary center of your coordinate system the stationary center of mass of the two bodies (called the "barycenter"). Then by conservation of momentum, the total momentum of the system will always be zero. If the net momentum of the two bodies is not zero, then first subtract the velocity of their center of mass from both to make it stationary.

What you normally see as the orbit solution is in reality the vector difference between the positions of the two bodies. If one of the bodies is effectively not moving because it is so much larger than the other, then you simply have that solution as the motion of the smaller body.

For the general case of the two-body problem, you solve for the vector difference of the two positions using the sum of the two masses as the central mass, getting the usual Kepler solution for that vector difference. Then to get the motions of the two actual bodies, you solve for two equations in two unknowns, where the difference of the two positions is your solution, and the sum of the momentums is zero.

Try this simple example in two dimensions:

$$\vec{r_1}=\left(2,0\right),\,\,\vec{r_2}=\left(-1,0\right)$$ $$\vec{v_1}=\left(0,2\right),\,\,\vec{v_2}=\left(0,-1\right)$$ $$\mu_1=6,\,\,\mu_2=12$$

where $\mu$ is the notation for the mass times Newton's gravitational constant $G$. Verify that the CG position and velocity are at the origin. Then determine the orbit for:

$$\vec{r}=\vec{r_1}-\vec{r_2}$$ $$\vec{v}=\vec{v_1}-\vec{v_2}$$ $$\mu=\mu_1+\mu_2$$

where $\vec{r}$ and $\vec{v}$ are the initial position and velocity of a massless object orbiting a central body with mass times $G$ of $\mu$.

You get $a=6$, $e=1/2$, $\omega=0$, and $\nu=0$. We then have:

$$\vec{r}=\left(6\left(-{1\over 2}+\cos\tau\right),6{\sqrt{3}\over 2}\sin\tau\right)$$

where $\tau$ is the eccentric anomaly. The time equation is the same for both the zero mass and general two-body cases:

$$t=2\sqrt{3}\left(\tau-{1\over 2}\sin\tau\right)$$

Solving the simultaneous equations, we get in general:

$$\vec{r_1}={\mu_2\over\mu}\vec{r},\,\,\vec{r_2}=-{\mu_1\over\mu}\vec{r}$$

That gives:

$$\vec{r_1}=\left(-2+4\cos\tau,\,\,2\sqrt{3}\sin\tau\right)$$ $$\vec{r_2}=\left(1-2\cos\tau,\,\,-\sqrt{3}\sin\tau\right)$$

Now we can plot it:

Update for updated question:

Here is the same Newtonian solution for PSR J0737-3039, where this presentation provides the masses, period, and eccentricity:

The scale on the plot is km, and the period is 2.45 hours (the animated gif runs a little faster than that).

The eccentricity is small, 0.0878, and the masses are close to equal, where A is 1.338 and B is 1.249 solar masses. As a result, they are in close to a common circular orbit about their center of mass. (A is blue and B is orange.) The positions in km and time in seconds are:

$$\vec{r_A}=\left(424000\,(-0.0878+\cos\tau),\,\,422000\sin\tau\right)$$ $$\vec{r_B}=\left(-454000\,(-0.0878+\cos\tau),\,\,-452000\sin\tau\right)$$ $$t=1404\,(\tau-0.0878\sin\tau)$$

where $\tau$ goes from $0$ to $2\pi$.

Their speed in their orbits is an impressive 0.1% the speed of light.

Answer 2

When solving the Kepler problem is often solved for the case in which the central body is stationary, such that the body orbiting it has a negligible mass.

For a two body problem this solution can still be used. For this the center of mass of the two bodies, also called the barycenter, can be chosen as inertial reference frame, if there are not external forces acting on this system. In this reference frame the following will then always be true,

$$ \mu_1 \vec{r}_1 = -\mu_2 \vec{r}_2, \tag{1} $$

where $\mu_i$ and $\vec{r}_i$ respectively are the gravitational parameter and position vector of body with index $i$. Equation $(1)$ can be used to simplify the acceleration experienced by body with index $i$,

$$ \ddot{\vec{r}}_i = \frac{\mu_j}{\|\vec{r}_j - \vec{r}_i\|^3} (\vec{r}_j - \vec{r}_i) = \frac{-\mu_j^3}{(\mu_i+\mu_j)^2 \|\vec{r}_i\|^3} \vec{r}_i. \tag{2} $$

The final expression of equation $(2)$ looks similar to that of the case in which the central body is so massive that it can be assumed to be stationary. The only difference is that the constant/gravitational parameter is different. This different constant can be used as the effective gravitational parameter, $\mu_{eff}$, for a fictional body positioned at the barycenter for the body with index $i$,

$$ \mu_{eff,i} = \frac{\mu_j^3}{(\mu_i+\mu_j)^2}, \tag{3} $$

such that this fictional body be assumed to be so massive, relative to the body with index $i$, that it can be assumed to be stationary. Finding the orbital elements will now be the same as when the central body is so massive that it can be assumed to be stationary.

It can also be noted that you only have to find the orbital elements of one of the orbits, because the orbital elements of the second orbit can be derived from the orbital elements of the first orbit. Starting with both orbital planes are the same, thus the longitude of ascending node, $\Omega$, and inclination, $i$, have to be the same. From equation $(1)$ one can be derive that: the eccentricities, $e$, have to be the same; the semi-major axes have to satisfy $\mu_1a_1=\mu_2a_2$; the arguments of periapsis, $\omega$, are separated by 180° (or $\pi$ radians); the true anomalies, $\nu$, have to be the same.