Escape Velocity / Newtonian Black Holes

Hopefully you completed the Slingshot with Gravity! coding activity. That activity explores the behavior of objects like planets or comets interacting with a massive object like the sun. An interesting part of that activity is showing that the gravitational force between two objects is inversely proportional to the distance squared. We found that this must be true: $$F_{grav} = \frac{G M m }{r^2}$$

and we have this formula instead of something proportional to $1/r^n$ because the planets are observed to have elliptical orbits and as the coding activity shows you only get elliptical orbits if $n = 2$.

In this activity we are going to keep using this formula and we are going to explore some of the consequenses of it. In particular, we want to know how fast you have to be going to escape the earth's gravitational field. We are also going to think about energy conservation and what it means for an object to be gravitationally "bound" by a massive orbit or not. These are all topics we dind't get to in the Slingshot with Gravity activity.

Step 0. Check out this modified version of the Slingshot with Gravity code

Click here to open up a modified version of the Slingshot with Gravity code

Here is a link to just run the code if that's all you want to do

Step 1. Check out the new features of this code

Step 1a. Click and drag to make the object drift far away from the earth

Unlike the Slingshot with gravity code, this program will zoom out if the object moves out of the field of view. To do this, click and drag to set the velocity vector and see what happens. It helps to have a nice long velocity vector so the object is moving fast. The speed is proportional to the length of the velocity vector.

Step 1b. Notice that you can die if you run into the Earth

Another feature of this code is that if the object runs into the earth then you will see a message saying "You're dead!" and the program will end. Click and drag to set a velocity vector towards the earth and show that this happens.

Step 1c. Notice that instead of x_sun and y_sun the code has x_massive and y_massive

We're not really thinking about the sun and planets like we were for the slingshot with gravity activity. So instead of using x_sun and y_sun we'll just use the variables x_massive and y_massive to indicate the position of the massive object. At least initially this massive object will be the earth. Later we will change it to a black hole. The important thing is that it is much more massive than the object that is orbiting it.

Step 1d. Notice that the picture of the earth is a famous one

We used a famous picture of the Earth from space called "The Blue Marble". This picture was one of the first photos of the entire earth and many credit this photo with launching the modern environmental movement in the 1970s.

Step 2. Configure the code so that the object launches from the surface of the earth

You can control the starting position of the object by changing this code:

x = 525;
y = 250;

to this:

x = x_massive + planet_radius;
y = y_massive;

What happens when you make this change? Does the object still orbit or do you die?

Step 3. Change the initial velocity of the object

Step 3a. Does having a velocity that is in the +x direction instead of +y direction help the object make it into orbit?

Change this:

vx = 0;
vy = 30;
to this:
vx = 30;
vy = 0;

What happens now? Does this allow the object to move far away?

Step 3b. Increase the velocity until the object escapes from the earth's gravitational field

We want to know how fast the object needs to be going to escape the earth's gravitational field, like we're going on an interplanetary road trip.

Note: This is NOT the same thing as the speed required to orbit the earth. If the object is escaping the gravitational field this means that it is moving farther and farther away, like a space probe launched from Earth heading to some other part of the solar system, rather than a satellite launched to orbit the earth indefinitely.

Modify the program to make the initial velocity larger and larger until this happens. Try this with the initial velocity in the +x direction and in the +y direction. Approximately what number do you get for when the object has enough velocity to escape as it does in this example?

Disclaimer #1: It is important that you find the smallest velocity that still allows the object to escape.Don't make the velocity any larger than it has to be! We will come back to this later.

Disclaimer #2: It is important to wait to see if the object falls back to the earth. It make take a little while for this to happen. Watch the dynamics for at least 20 seconds.

Step 4. Energy Conservation

Step 4a. Think about Energy Conservation

We are going to use energy conservtion to find an expression for what this minimum initial velocity for escape should be. The kinetic energy plus potential energy of an object moving in a gravitational field should add up to a constant: $$ KE(t) + PE(t) = E_{\rm tot} $$

As usual, $ KE = (1/2) m v^2 $, and because we are dealing with a gravitational field we use $ PE = - G M m / r $. $$ \frac{1}{2} m v^2 - \frac{G M m}{r} = E_{\rm tot} $$

Step 4b. Think about whether $E_{\rm tot}$ is >, < or equal to zero

Depending on the situation, $E_{\rm tot}$ may be positive, negative or zero. Imagine three different scenarios: (1) the object is launched so fast that it moves very far away from the planet (like you found in step 3b), (2) the object is launched but it crashes into Earth and (3) the object is launched and it doesn't crash into earth but it doesn't move very far away either. The object just orbits over and over again.

For each of these three situations, is $E_{\rm tot}$ positive, negative or zero? Use the equation for energy conservation to explain why.

Advice: Remember that $m > 0$, $G > 0$, $r > 0$ and $v^2 > 0$. The only thing that is negative is the minus sign in the middle.

Step 5. Escape velocity

It turns out that if the object is launched very fast and is moving far away from the planet then $E_{\rm tot} > 0$, but if the object is launched much more slowly then $E_{\rm tot} < 0$ and the object is in some kind of orbit. So if $E_{\rm tot} = 0$ this means that the object has just barely enough speed to escape the system.

Use $E_{\rm tot} = 0$ to find a formula for $v$ when this is true. This $v$ is called the "escape velocity". How does this result compare to the number you got for the velocity in Step 3b?

Step 6. Fractional Error

Take the number you calculated for the escape velocity from Step 5 and compare it quantitatively to what you found from the program. Treat the calculation as "true" and the program as "measured". This is important because the program is breaking up time into little steps, so we do not expect the results to always be accurate (In general NEVER assume that the computer is always right!)

Here is the formula for fractional error: $$ {\rm fractional \, error}\, (\%) = 100 \cdot \frac{|v_{\rm meas} - v_{\rm true}|}{v_{\rm true}} $$

You should get a fractional error that is less than 10%. If this doesn't happen, maybe go back to Step 3b and see if you can make that initial velocity a little bit smaller and still find that the object moves very far away.

Step 7. Finally the black hole part!

Ok, instead of the earth let's have a black hole with the same mass as the earth (imagine a mad scientist experiment went terribly wrong).

Replace this:

drawPlanet(x_massive,y_massive,planet_radius,0,0,0,0);
 

with this:

drawBlackHole(x_massive,y_massive,0,0,0,0);	

Step 8. Not so fast! The initial speed can't be faster than the speed of light!

Let's assume that we can't hurl our object any faster than the speed of light. In the next few sub-steps we're going to make sure the object can't move any faster than the speed of light.

Step 8a. Set the value of c at the beginning of the code

The first step to making sure objects don't move faster than the speed of light is to put a variable near the beginning of the code that specifies the speed of light. We are going to set c = 50; in our code. If this seems a little weird, bear in mind that we are assuming that the gravitational constant is G = 100; instead of $G = 6.67 \cdot 10^{-11}$. We're just trying to get the essential behavior to play out in our little screen on a timescale that isn't ridiculously slow or fast. It turns out that c = 50; is in this goldilocks zone that is neither too fast for us to see or too slow for us to have the patience to watch!

Go ahead and add this to near the beginning of the code, and before the part where you set vx and vy!

  c = 50;

So your initial velocity section should look like this:

c = 50;

vx = 0; // or whatever value you have
vy = 30; // or whatever value you have

Step 8b. Limit the speed of the object to c

We need to add one or more if statements that will check if the object is moving faster than the speed of light and if it is moving faster the if statement will need to fix this. This is going to be more tricky than it sounds! Here is a group of four if statements that mostly fixes this problem:

  if ( vx > c ) {
       vx = c;
  }
  if ( vx < -c ) {
       vx = -c;
  }
  if ( vy > c ) {
       vy = c;
  }
  if ( vy < -c ) {
       vy = -c;
  }      

After you put this in your code, the program should behave like this

We will use this in what follows, but do you notice any potential or actual problems with this code? There is more than one! In what situations would this code fail or produce unwanted behavior? (Hint: The velocity and speed are related by this formula, $ {\rm speed} = \sqrt{v_x^2 + v_y^2} $).

Step 9. Where is the boundary between bound and unbound?

Earlier we increased the velocity until the object was clearly not gravitationally "bound" to the planet. Now we're going to keep the velocity equal to c in either the +x or +y direction, but we're going to change the initial position to see if the object is bound or not.

Step 9a. Set the initial velocity to c

Change this:

c = 50;
    
vx = 0; // or whatever value you have
vy = 30; // or whatever value you have

to this:

c = 50;
   
vx = 0;
vy = c;

Step 9b. Change the initial position until the object can escape to far away

In Step 2 you changed the inital position of the object to this:

x = x_massive + planet_radius;
y = y_massive;

Replace planet_radius with its numerical value (which happens to be 50) so you end up with this:

x = x_massive + 50.0;
y = y_massive;

Try this out and notice that the object is still bound to the black hole. Increase 50.0 to larger numbers until the object is unbound!

If you increase this value large enough, eventually the object will drift very far away with no sign of returning as in this example

What value do you get? Does it depend on whether you assume vy = c; and vx = 0; or do you get the same result from vx = c; and vy = 0; ?

Step 10. Does it depend on the mass of the object?

So far we've been discussing the behavior of an object of non-zero mass that is experiencing a gravitational force from the Earth or a black hole. For black holes we are really more interested in how light behaves near a black hole. A key property of light particles is that they have no mass. So does our program model light near a black hole?

Decrease the mass of the object in the program (the default value is mass = 3.0;) but do not set it to zero. Why is it important to NOT set the mass to zero even though that's what we're interested in? (Hint: it's a code reason)

Does the mass of the object affect the result? If it doesn't then our program is telling us something about how light would behave.

Step 11. Use energy conservation to explain the number you got in Step 9

Here again is the equation for energy conservation: $$ \frac{1}{2} m v^2 - \frac{G M m}{r} = E_{\rm tot} $$

Assume that the object is moving at the speed of light (so $v^2 = c^2$) and use this to solve for $r$ if the object is just above the threshold for being unbound ($E_{\rm tot} = 0$).

What do you get for $r$? How does it compare to the number you found in Step 9? What is the fractional error?

Step 12. Draw this radius on the screen

Go ahead and draw a circle on the screen with the radius you found in Step 11. You can use this function:

drawCircle(x_massive,y_massive,0.0);

Make sure to replace 0.0 with the number you came up with in Step 11! When you are done your code should behave like this

Double check that speed-of-light moving objects that start just inside this circle are trapped in orbit while objects that start just outside of it are free to travel very far distances away.

Step 13. Think about what this means

Imagine we are observing this system from very far away. If a particle of light were to cross through this circle, would we ever see it in our telescopes?

Note: The radius you found is often called the Schwarzchild radius. By another name it is the "Event Horizon" of the black hole, which is a kind of point of no return.

Step 14. Newtonian Black Hole

Our program is now a description of a "Newtonian" black hole. There is still quite a bit more work to be done to make our code into something that would model a black hole more realistically (including relativistic effects). What kind of effects still need to be added? (Hint: the movie Interstellar will help you figure out one of the things that we haven't included.)

Challenge: Getting the launch angle right while keeping the speed = c

In Step 8b. we gave you some code to limit the velocity of the object to the speed of light, but we said that there were some problems with this code. One of the problems is that if you drag the velocity arrow so that it exceeds the speed of light in one direction (say x) but you are not exceeding the speed of light in the other direction (say y) then when those if statements correct the velocity, it may accidentally change the angle you wanted the object to move (go ahead and try this out and see if you notice).

We want to get rid of those if statements and replace it with something that preserves the launch angle that the user specified. Here is some code to get you part of the way there:

   // if object is moving faster than speed of light
   if (sqrt(vx*vx + vy*vy) > c) { 
      vel_angle = atan2(vy,vx);
      vx = ????   // what goes here?
      vy = ????   // what goes here?
   }

Put this into your code and figure out "what goes here?" In the end your code should behave like this