# Calculating Stable Circular Orbits for Binary Stars

Binary stars are cool. But if you have two stars interacting with a gravitational interaction, there are technically 6 coordinates you need to deal with (3 vector components for each star). It’s possible to reduce this to an equivalent 1 dimensional problem, but I’m not going to go over that today. Instead, I’m going to show you two things: how to calculate the initial conditions for stable orbits and how to model a binary system in python.

Let’s get started.

**Initial Conditions**

Suppose we have two stars (because it’s really hard to have a binary star system with only one star). For now, I will just say that one star is twice the mass of the other star and they are separated by a distance “r”. Here’s what I want:

- The center of mass of the binary star system to be at my arbitrary origin (just because).
- Momentum of the center of mass to be zero (the zero vector).
- Circular orbits.

First the center of mass. I’m going to put my 2M star on the negative x-axis and then put my 1M star on the positive x-axis. Like this:

The center of mass can be calculated as:

I want this center of mass to be at <0,0,0> (the origin). But since the two masses are already on the x-axis, I just need to solve for the x-positions of star 1 and 2. But don’t forget, I also know the distance between the two stars is r.

Note that x_1 has a negative value — right? OK, I have two equations with two unknowns (x_1 and x_2). I’m going to solve for x_2 from the r-equation and sub into the center of mass equation. I can multiple both sides by 3M to make it simpler. From that, I get:

If we have some arbitrary distance units, then the more massive star (star 1) would have an x-position of -1 unit. Star 2 would have a position of 2 units. That puts the center of mass for this system at the origin.

Now, what about the momentum of the center of mass? Because of the definition of the center of mass, the momentum of the center of mass is equal to the sum of the momentum of the two stars.

For this situation, I can use momentum as the product of mass and velocity. So, if the center of mass momentum is zero (vector), then the momentum of the two stars will be equal and opposite. This means I just need to find the velocity (and momentum) of one of the stars and I will get the other one too.

Of course a star could start with any particular velocity — but I don’t want just any velocity. I want a velocity that will produce circular (and stable) orbits. Here is a diagram showing two stars in circular orbits.

Yes, they both have circular orbits with the center of mass at the center. You can trust me — but you don’t have to because I’m going to show you this below.

Why do the stars move in a circle? It’s because of the gravitational force. The magnitude of this force can be calculated as:

This is just the magnitude of the force. It pulls towards the other star. Also, both stars have the exact same force on them — just in opposite directions. I’m using the scalar version of this force just to make things simpler. The “r” in this case is the distance from star 1 to star 2. But there you have it. Oh, G is the gravitational constant.

So, what do forces do to an object? A net force makes an object accelerate. Since acceleration is the rate of change of velocity and because velocity is a vector, an object moving in a circle (even with a constant speed) will have an acceleration. This is the centripetal acceleration. Let me write down the centripetal acceleration for star 1 (the one with a mass of 2M).

Again, this is the magnitude of the acceleration for star 1. The direction of the acceleration is towards the center of the circle — which just so happens to be the same direction as the gravitational force (isn’t that nice). However, there is something very important with both the gravitational force and the acceleration. The gravitational force depends on the distance between the two stars (r) and the acceleration depends on the radius of the orbit (r_1). These are not the same.

OK, let’s put these two together. Newton’s second law says that the net force is equal to the mass multiplied by the acceleration. I’m going to write this as a scalar equation since there’s only one force and it’s in the same direction as the acceleration. Here it is for star 1.

Solving for v_1:

But wait! I already found the value of r_1 in terms of r (from the center of mass). That means the velocity of star 1 will be:

That’s it. I just need to multiply by the mass of star 1 to get the momentum. Star 2 starts off with the same momentum but in the opposite direction. But does this really work?

Oh, I just realized something. I should have solved this for the generic binary star system and not my 2M to M. Oh well. You can do that for homework.

**Modeling a Binary Star System in Python**

I’m going to use GlowScript VPython since it has built in stuff to visualize motion.

So, we have these two planets moving in 3D. They have 3 dimensions for the momentum. They have 3 dimensions for the positions. The force is in 3D. There’s a lot of stuff going on. Oh, it gets even worse. The forces aren’t even constant. Even for the case of a star moving in a circle, the direction of the gravitational force changes. It’s just a tough problem.

This is a perfect time to use a numerical calculation. The basic idea is to break the problem up into small (where small is a relative term) time intervals. During each time interval, I can assume the forces are constant. I can use this constant force to find the momentum and position at the end of the time interval. After that, I just do the same stuff for the next time interval. Since no one wants to do a whole bunch of these calculations, we can just make a computer do it (computers rarely complain).

Here is the basic recipe I’m going to use:

- Start with initial conditions (good thing I just calculated that).
- Find the vector position from one star to the other (this is the vector r).
- Use r to calculate the vector gravitational force.
- Use the force to update the momentum of both stars.
- Use the momentum to update the position of both stars.
- Keep doing this until you want to stop.

Yes, I glossed over some of the details — but you get the idea. The main thing is that I want to show you that the calculation for circular orbits does indeed work.

Oh, what about the length of the time interval? If this was a ball tossed across the room, I might use a step interval of 0.01 seconds. But a star is moving ginourmous distances. Even with super fast orbital speeds, it would still be basically in the same place after 0.01 second. Then what would be an appropriate time interval? I’m just going to guess — if I’m wrong, I will pick a different value. Let’s go with 1000 seconds.

Since this is a numerical calculation, I need actual numbers. I can’t use “r” and “M”. I will just pick some values that seem legit for a star (like our Sun) — but you can change them if you want. Here is the actual code for you.

Let me go over some parts of this code. This is the first part.

Comments:

- Line 3–5 are just some constants and values.
- Lines 7–12 create the two stars (with masses) and the center of mass. Yes, I put the center of mass as a visible sphere.
- Line 11: here I actually calculate the vector location of the center of mass and then print it just to make sure I didn’t mess something up.
- Line 16 is the magic. That’s the calculation of the velocity of star 1 for a circular orbit. Velocity is a vector, so in line 17 use this value for the y-component of the star’s velocity.
- Line 20: star 2’s momentum is opposite of star 1.

Here is the rest of the code where the stuff is actually calculated.

Comments:

- Line 27 is the start of the loop. I mean, who wants to write out the all those calculations over and over? This loop just repeats those calculations until some particular time. Yes, that time is 10 million seconds.
- Line 29: this one is special for GlowScript. It tells the program to not do any more than 1000 loops every second. If you increase this number, the simulation will run faster. Of course you don’t want it to run in real time — ain’t nobody got time for that.
- Line 30–31 calculates the vector from star 1 to star 2 and then calculates the vector force.
- Lines 32–36: update momentum, update position, update time. Done.

OK, now check it out. You should of course really run the code (it’s free), but this is what it looks like.

That’s not physics, that’s art. But there you have it. It’s a stable and circular orbit using the initial velocity that I calculated. Both stars are orbiting the center of mass. If you want, you can change the initial velocity of star 1 and get some non-circular orbits. It’s fun.