# Black hole orbits

• I

## Main Question or Discussion Point

Now my Kerr orbit simulator is pretty much feature complete, I have started to look at producing videos . . .

I have just started a channel on YouTube to accumulate some of the more interesting examples. Aside from creating the simulation, the most difficult part was to generate useful initial conditions (trial and error is of little use!), and as a result of this I can create and check my results against the literature, so I can assure you the results are accurate.

The basic equations are here*, and I have tested the output against (amongst other sources) http://staff.science.nus.edu.sg/~phylyk/downloads/reports/sp2172_report.pdf [Broken] and this one.

If anyone can think of a particular scenario they would like illustrated, ask and I'll see whether I can do it. If you want to have a go yourself, the sources are on GitHub as usual.

* this used to be on Arxiv (or similar) but seems to have been removed :(
Wilkins, D.C., "Bound Geodesics in the Kerr Metric", 1973, (some paywalled journal)

Last edited by a moderator:
• ProfuselyQuarky and edguy99

Related Special and General Relativity News on Phys.org
Thanks for the post! This is an automated courtesy bump. Sorry you aren't generating responses at the moment. Do you have any further information, come to any new conclusions or is it possible to reword the post?

why did you decide to do this project in python?

why did you decide to do this project in python?
Well, it originally started in Java (with commons-math minimization), but I wanted to compare the algorithms with Python (with Scipy for minimization) for speed / ease of programming (turns out standard Python is very slow, I used PyPy and performance became comparable). Python was the also obvious choice for the 3D plotting and graphs so it went all Python (but I was getting annoyed with all the "self" references that made the code hard to read).
Then I discovered Vala (using the gsl bindings for root solving) and converted all the computational stuff to that (it is very similar to Java/C# and compiles to native executables via C).
Note also that shell scripts play a central role in operation, build and testing.
I still have no better solution for the graphics/plotting, so that is likely to stay Python. Any suggestions welcome, although it is nice to be in maintenance mode for a while so no rush!

Now my Kerr orbit simulator is pretty much feature complete, I have started to look at producing videos . . .

I have just started a channel on YouTube to accumulate some of the more interesting examples. Aside from creating the simulation, the most difficult part was to generate useful initial conditions (trial and error is of little use!), and as a result of this I can create and check my results against the literature, so I can assure you the results are accurate.

The basic equations are here*, and I have tested the output against (amongst other sources) http://staff.science.nus.edu.sg/~phylyk/downloads/reports/sp2172_report.pdf [Broken] and this one.

If anyone can think of a particular scenario they would like illustrated, ask and I'll see whether I can do it. If you want to have a go yourself, the sources are on GitHub as usual.

* this used to be on Arxiv (or similar) but seems to have been removed :(
Wilkins, D.C., "Bound Geodesics in the Kerr Metric", 1973, (some paywalled journal)
Couldn't edit the OP, but tracked down a Wilkins paper link. Turns out it was 1972.

Last edited by a moderator:
edguy99
Gold Member
The current position of the stars around the central black hole in the milky way (sagatarius A*). Not sure if they have stable orbits or are on their way to a collision. Is it easy to add a time clock on the bottom to see the passage of time and know what time frame the animation is running at?

The current position of the stars around the central black hole in the milky way (sagatarius A*). Not sure if they have stable orbits or are on their way to a collision. Is it easy to add a time clock on the bottom to see the passage of time and know what time frame the animation is running at?
If you are asking what I think, then this simulation is not what you would use to analyze Sag A*, because there are no closed form N-Body solutions in GR.

As far as I can tell, the system has been studied using N-Body Newtonian simulators, or even separate 2-body ones. You have to get things really close for GR effects to kick in.

Now, it just so happens that by sheer coincidence, there is a Newtonian N-body simulator included in this package, as an independent test of the integrator ;) The only data file I have for it is our Solar System. I would love to get hold of suitable initial conditions data for Sag A* and try it out!

edguy99
Gold Member
Some initial data (although only a 2d view) of actual observations from 1995 to 2002. It would be easy to estimate the x/y position from this, but not sure where or how you would get the z axis position.

Some initial data (although only a 2d view) of actual observations from 1995 to 2002. It would be easy to estimate the x/y position from this, but not sure where or how you would get the z axis position.
Yes, in theory I could do this, and you could pan/zoom as you like while stuff is whizzing around. But unfortunately positions are not enough. I need 3D position data, as well as 3D velocity or momentum, and all the masses too.

If anyone can point me to a suitable source of data, it would be really quite cool (but not GR). Something like this, for example, which I used to generate initial conditions for the solar system simulations.

Since it's a quiet day, two more added, 12 in all now.

Both involve "capture" (or time-reversed escape) into strictly unstable orbits around a non-spinning black hole, the first is for a light pulse onto the photon sphere at r = 3M, marked as an orange ring. The second is for a particle into the region between the photon sphere and the Innermost Stable Circular Orbit (ISCO) marked as a magenta ring at r = 6M.

Some initial data (although only a 2d view) of actual observations from 1995 to 2002. It would be easy to estimate the x/y position from this, but not sure where or how you would get the z axis position.
BTW found this thread which might interest you.

edguy99
Gold Member
Just saw one more suggestion for you coming out of todays AAS press conference on black hole mergers. Not sure if you can do these, but they would be fun to watch. :)  Just saw one more suggestion for you coming out of todays AAS press conference on black hole mergers. Not sure if you can do these, but they would be fun to watch. :)
I'm afraid that is firmly in the realm of Numerical Relativity, which is well above my pay grade! If you want to watch something here is a typical simulation. This is on-topic because these are the sort of simulations that LIGO use to generate template signals for matching against their observations.

I have updated my code so it now also displays the local shell velocity and the shapiro delayed coordinate velocity. The conversion factors are $\nu_{\parallel} = v_{\parallel} \cdot (1-r_s/r)$ for the radial, $\nu_{\perp}= v_{\perp} \cdot \sqrt{1-r_s/r}$ for the transversal component and $\nu = \sqrt{\nu_{\parallel}^2+\nu_{\perp}^2}$.

The input parameters are initial radial distance, local or coordinate velocity, plane angle and the local launch angle of the test particle, or alternatively the initital cartesian coordinates and velocity components. The output is given either in proper or coordinate time steps, while the latter takes more CPU ressources. The syntax is Mathematica.

Comparison of newtonian and einsteinian orbits with the same initial conditions: +2% to the circular orbital velocity ar r=2.5rs gives an eye

http://org.yukterez.net/2.5rs,1.02vo,newton,E.gif http://org.yukterez.net/2.5rs,1.02vo,einstein,E.gif

while an increase of the circular orbital velocity by 111‰ at 2.8rs gives a flower

http://yukterez.net/org/2.8rs,1.111vo,newton,E.gif http://org.yukterez.net/2.8rs,1.111vo,einstein,E.gif

as one can see on the icons, Newton is on the left and Einstein on the right. It might also be of interest to see a high enery particle fall towards a black hole in coordinate time steps where it delecerates and freezes at the horizon at rs=2GM/c², while with newton it reaches infinite velocity at r=0 (the v0 is 0.9999c at r=5rs)

http://yukterez.net/org/c,newton.vs.einstein.gif

On the right it actually never reaches c in external time, the 1.0 at the end of the animation in the local velocity line is just displayed because there are so many 9s behind the 0. that it is rounded up to 1.0

Last edited:
On the right it actually never reaches c in external time, the 1.0 at the end of the animation in the local velocity line is just displayed because there are so many 9s behind the 0. that it is rounded up to 1.0

First a couple of questions; are you using proper time as your integration variable, and are you using a variable step integrator?

I suspect what you are seeing might not be "freezing" at the horizon, but your integrator reducing the time step towards zero as the (coordinate time) error increases towards the horizon.

I just verified that if I don't stop my HTML simulator at the horizon, it goes all the way to $r = 0$. You can try changing my models.js file line 275 to

Code:
        if (this.r > 0.0) {
And see for yourself with these parameters: If you watch the times, you should see the coordinate time increase up to the horizon and decrease slightly inside, whilst the proper time just ticks away merrily until the end.

Last edited:
Comparison of newtonian and einsteinian orbits with the same initial conditions: +2% to the circular orbital velocity ar r=2.5rs gives an eye

while an increase of the circular orbital velocity by 111‰ at 2.8rs gives a flower
Just for fun, since it's another quiet day . . .

Here's an eye: and here's a flower: I have no idea why that is and can't edit it anymore, but luckily I have put the url as a watermark on the images (bottom right). Here it is again: yukterez.net/org/schwarzschild,code.txt

First a couple of questions; are you using proper time as your integration variable, and are you using a variable step integrator?
The equation of motion differentiates coordinate distances by proper time (so internally the rapidity, which is infinite when v=c and imaginary when v>c, not the velocity is used). When I plot orbits that don't get near the event horizon I use a standard Runge Kutta or even automatic variable step size at MachinePrecision because it is much faster but accurate enough when the trajectory doesn't tipple the horizon and not too many revolutions where the small errors could accumulate are plotted. When the particle gets near the event horizon I include a StiffnessSwitcher and when plotting in coordinate time steps also an increased WorkingPrecision for the last micrometers before the event horizon.

The same scenario in proper (left) and coordinate (right) time steps: (The local initial launch angle of 85.5° seems flatter since the radial shell distances are contracted relative to the circumference, which is a strong effect at r0=2.002GM/c²)

I suspect what you are seeing might not be "freezing" at the horizon,
But that's what is expected to happen in external coordinate time. When I switch to proper time steps in the animation frames the particle penetrates the horizon and ends in the singularity at the center.

but your integrator reducing the time step towards zero as the (coordinate time) error increases towards the horizon.
If I plot in coordinate time steps the coordinate time increases while the proper time freezes asymptotically. If I plot in proper time steps the coordinate time goes to infinity when approaching r=rs. That's ok, the increasing error near the horizon in the recursive numerical solving from t(τ) to τ(t) can be reduced with increasing precision if needed.

I just verified that if I don't stop my HTML simulator at the horizon, it goes all the way to r = 0.
In proper time this is ok, in coordinate time this shouldn't happen in a finite time.

If you watch the times, you should see the coordinate time increase up to the horizon and decrease slightly inside, whilst the proper time just ticks away merrily until the end.
If you trace the test particles world line this are the numbers you shoud get for the coordinate time, but inside the horizon they are not meaningful physically for an outside observer: Last edited:
• edguy99
To understand why the local launch angle of almost 90° in my last posting seems like almost flat in Schwarzschild coordinates one has to look at Flamm's Paraboloid: where the x-Axis represents the radial distance as it is seen by a coordinate observer, and the true number of rulers is the one that would fit on the blue graph; so between the coordinate r1=2GM/c² and r2=24>GM/c² there fits a ruler not of length r2-r1=22GM/c² but of length ∫(1/√(1-rs/r), r=r1..r2)=26.8GM/c². The same goes for the velocity components, so a grid that would look like this: with Newton would look like this in Schwarzschild coordinates: On the shell ar r0=1.001rs where the test particle is lauched the radial component is contracted by a factor of more than 30 relative to the transversal component. That's why an almost straight up launching angle seems like an almost flat angle for an external observer.

I have no idea why that is and can't edit it anymore, but luckily I have put the url as a watermark on the images (bottom right). Here it is again: yukterez.net/org/schwarzschild,code.txt
Don't see the watermarks but got it now anyway, thanks. I don't normally even look at closed source stuff, but at least I can see what you are doing and run this on a Raspberry Pi. If I were to actually do this myself, I would use GNU Octave.

If I plot in coordinate time steps the coordinate time increases while the proper time freezes asymptotically. If I plot in proper time steps the coordinate time goes to infinity when approaching r=rs. That's ok, the increasing error near the horizon in the recursive numerical solving from t(τ) to τ(t) can be reduced with increasing precision if needed.

In proper time this is ok, in coordinate time this shouldn't happen in a finite time.
OK just thought I'd mention it as I think I suffered from this a couple of years back when I last used the geodesic equations for this stuff. All water under the bridge now . . .

edguy99
Gold Member
...I just verified that if I don't stop my HTML simulator at the horizon, it goes all the way to $r = 0$. You can try changing my models.js file line 275 ...
Are videos displayed in coordinate time or proper time? From the video, I had assumed they we are in coordinate time (ie. what you would see from earth looking at the black hole in the center of the milky way), but now I see they must be proper time (from the particles point of view that is falling into the black hole). Still a little confused.

... In proper time this is ok, in coordinate time this shouldn't happen in a finite time ...
I agree.

Are videos displayed in coordinate time or proper time? From the video, I had assumed they we are in coordinate time (ie. what you would see from earth looking at the black hole in the center of the milky way), but now I see they must be proper time (from the particles point of view that is falling into the black hole). Still a little confused.
All my simulations are done in particle proper time, but can be plotted in coordinate time (well, against any of the coordinates really if you want to . . . ). The main difference when plotting against coordinate time is that particles travel slower near to the mass and faster far away. I think it's more dramatic to have them "whoosh" by, and requires no post processing ;)

Of course, if I were to simulate n test particles, it would only make sense to use coordinate time to show where they "are".

Of course, if I were to simulate n test particles, it would only make sense to use coordinate time to show where they "are".
If I had to simulate more than 1 particle free falling I would also use coordinate time steps: But i was also wondering where the red particle would be in terms of the blue particles system after the blue particle has already disappeared behind the horizon. I'll see if I can figure that out too, but right now I have no idea how that would look like.

Last edited:
• edguy99
But i was also wondering where the red particle would be in terms of the blue particles system after the blue particle has already disappeared behind the horizon. I'll see if I can figure that out too, but right now I have no idea how that would look like.
Off the top of my head, you might want to look at "rain coordinates" (Gullstrand-Painleve), but to be honest I'm still wondering precisely what the question means ;)

In the above example I can say
in coordinate time, when the blue particle is at r=4.277, the red particle is ar r=11.61