Project 2 - Stability of Planetary Orbits in a Binary Star System

Introduction

In this project, we will explore the effect of the gravitational forces at work within a binary star system on the stability of planetary orbits. This was an interesting theoretical question for a long time, but with the discovery of a planet in the Gamma Cephei system, it is clear that, under at least some circumstances, planets do form and survive in a binary system.

Project Calculations

The calculations for this project build upon the program you wrote for Exercise 14, with the main difference being a new, more complicated, force calculation due to the second star in the system. The coordinate system is three dimensional; we will actually set up the problem in such a way that the motion will be entirely in the X-Y plane. The initial conditions for our test planets should follow what was done for Exercises 14. Begin the calculation at t=0 with the planet at some location: x=r_test; y=0; z=0. Where r_test is the radial distance that we'd like to test in the calculation. The initial velocity for the test asteroid should be entirely in the y direction. We will vary this value to see what sorts of stable orbits might exist at different locations within the system.

Three Body Forces

The project was introduced in class with a lecture which may be found here. The equation of motion is:

where the Force includes the contributions from the gravitational attraction to the Primary and Secondary stars according to:

where the geometry is illustrated in Figure 1.

As for Exercise 14, I advise carrying out the calculation using Astronomical Units for the unit of length and Years for the unit of time. In this case:

Binary Star Position Model

We are going to simplify our binary system and assume that the stars move in a simple circles about the center of mass of the system. Set the distance between the two stars (D = r_s1+r_s2) so that the period of their mutual orbit is 50 years (approximately that of Gamma Cephei). Thus, we may describe the positions as follows:

where

and D and P are the distance between the stars and period of the mutual orbit respectively. What is D? Hint: Use Kepler's Third Law! Note that the initial conditions lead to the stars and the test planet aligned at t=0.

We will make our computations for a value of the secondary star's mass relative to the primary star of: M₂ = 0.25 M₁. Keep in mind that it might be fun to write your program so that you can try different relative values (see extra credit problems)!

Step by Step

Step 1: First we need to modify the Exercise 14 program so that it takes into account the forces due to two stars rather than one. Note that in this new system, the origin is the center of mass, which falls between the two stars.

Your new program should be capable of producing three graphs of the results:

1. a simple plot of the planet position in the inertial system (x and y). In this graph, the stars move on their circular path in orbit about one another.
2. a plot of the planet position in the coordinate system that rotates with the orbit of the stars. In this case the stars are at a fixed position in the plot and the test planet moves around.
3. [optional] a graph which shows that planet position plotted at time invtervals of exactly one stellar orbital period in the rotating coordinate system. This one is a bit tricky, but it is useful in assessing the behavior of the orbits.

It is also useful for your program to calculate the Jacobi Integral, which is a constant of the motion for this problem. This is useful to verify that the program is working and it will be useful for understanding the behavior of orbits later on.

Step 2: You will want to test this program against analytical solutions and be sure that you select a step size that gives adequate accuracy. What is adequate? You be the judge and be prepared to defend your decision in your report. You might find that there is no one good time step for all the orbits we consider on our problem!

HINT: In considering tests we need to compare our program to something that is predicted about the restricted three-body problem. Computing a two-body orbit solution, as in Exercise 14, allows a check of the orbital integration formulae, so a simple test would involve trying the program in a limiting case with secondary star mass MUCH smaller than the primary star mass to check the program. But that's only part of the story. Are there any tests that would really provide a check on the 3-body forces?

Step 3: Once we have a working program, we want to use it to investigate the stability of orbits. There are three obvious stable cases: (1) orbit around primary star; (2) orbit around secondary star; and (3) orbit outside the entire system at a great enough distance so that the difference between two stars and a single star at the center of mass is not large. Demonstrate that your program will produce a stable orbit for each of these cases.

In order to demonstrate stability, we are going to demand that the test planet remains in the system for 100 periods of the binary system, corresponding to 5,000 years. This is actually not a very tough standard for stability. Real calculations would attempt to show stability for 100's of millions of years! Even so, your program will take some time to run. Be sure to check it on short runs before your commit to doing a long one. Also, as discussed in class, be sure to avoid inefficient steps in your loop: Be especially careful to avoid allocation of memory in the loop!

Step 4: Now we can evaluate the stability of more interesting cases where the test planet orbits between the stars.

I think it is easiest to do this in a structured way, running the program for a fixed initial distance from the center of mass and then varying the initial velocity in some way. Essentially we are doing what we could have done with Exercise 14 to show how the shape of the orbit changes as the initial velocity is changed. If you have not tried this with your program from Exercise 14, take a few minutes to run the program with initial position= [1, 0, 0] and try different values for initial velocity to see the resulting shapes.

For our project, start the particle at a location along the x axis that is on the opposite side of star 1 from star 2 (see figure) with a negative y velocity, corresponding to the same sense of revolution as the stars (called the prograde direction). For a fixed initial position, try a range of velocities (which corresponds to a range of values of Jacobi Integral) and graph the results. As you start exploring, it is OK to run the program for a short time, say one period of the star orbit, to get an idea of what's going on. It is also useful to include some sort of test to tell you whether the test planet has run away in order to save time calculating the path of the escapee.

There are a variety of wonderful orbits to be explored. Try to explore the range of possibilities in a systematic way, as outlined above. It is easiest to see them in the plot of the planet's path in the rotating frame. Here are a few of my favorites:

All figures above have been started at the same position, but with different initial velocities. The red circle indicates the position of the primary star; blue circle indicates position of the secondary star; green circle indicates the initial position for the calculation. The title for each plot shows the value of the Jacobi Integral.

Questions:

1. Looking at the behavior of the orbit as the initial velocity is changed: how is the behavior of the test planet orbit in the binary system different from the behavior in a single star system?
2. There is a special case for a circular orbit. What's the value of the initial velocity required to make a circular orbit about the primary star for your chosen initial position? How does this value agree quantitatively with the value expected for an orbit around the primary star alone?
3. How does the distance of the planet from the primary star of the planet change during the orbit when you try cases with initial velocities just above and just below the critical circular value? Can you explain this behavior based on the behavior of two-body orbits?
4. Do you think there are likely to be stable orbits in this system other than the ones explored in step 3?

Extra Credit Problems: There is a lot we could explore about the restricted three-body problem with this program. The following are interesting things to try. They don't require major changes in the program, just more runs and a little thought to answer a question. I'll credit 1 bonus point for each, assuming correct answers of course.

• Investigate orbits within the stellar system with a different mass ratio for the binary stars.
• How close can an orbit exterior to the binary system come to the orbit of the secondary star before it becomes unstable?
• Run the program starting the test planet at the L4 point and show that this is stable. The result is kinda cool.
• Building on the previous run, try runs at L4 with different energies (e.g. a little extra velocity above the stable value). How big a "kick" can an object tolerate and still remain locked to the L4 point?
• Investigate the behavior of an object at the "unstable" Lagrange points: L1, L2, L3. In fact, in the solar system, L1 and L2 are used to "park" spacecraft in a fixed position with respect to Earth and Sun. What do you think of these locations as parking places? Which is the best?

Submit your Report

Write a lab report on the project, following the Project Report guidelines presented here. You may submit a paper copy or copy the report file to the directory: ~schloerb/ph281/username/Project2

Go to Ph281 Home

Go to Ph281 Topics