Black hole orbits

[m4r35n357]

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)

 

[serp777]

why did you decide to do this project in python?

 

[m4r35n357]

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!

 

[m4r35n357]

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.

 

[edguy99]

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?

 

[m4r35n357]

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]

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.

 

[m4r35n357]

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.

 

[m4r35n357]

Just updated YouTube videos with some more light pulse orbits, fleshing out the range of photon spheres for a maximally spinning black hole:

https://www.youtube.com/playlist?list=PLvGnzGhIWTGRuIsHzDLdeWS4GSZ9CIY-b

Still not happy with video quality, but they are better than nothing for now.

 

[m4r35n357]

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.

 

[m4r35n357]

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]

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. :)

 

[m4r35n357]

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.

 

[Yukterez]

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 energy 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

 

[m4r35n357]

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

BTW your code link appears to be broken.

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

        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.

 

[m4r35n357]

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:

 

[Yukterez]

BTW your code link appears to be broken.

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:

 

[Yukterez]

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.

 

[m4r35n357]

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]

...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.

 

[m4r35n357]

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".

 

[Yukterez]

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.

 

[m4r35n357]

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 ;)

 

[Yukterez]

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
but when I ask
when the blue particle is at r=0.1, at which r is the red particle at that time in the blue particle's plane of simultaneity?
I think I really have to switch to another coordinate system than Schwarzschild, but maybe I also just don't see the obvious. I'll have to take a nap over that.

 

[PeterDonis]

I think I really have to switch to another coordinate system than Schwarzschild

Yes, you do. The Schwarzschild coordinate lines of simultaneity outside the horizon do not extend inside the horizon. So the question you are asking is meaningless in Schwarzschild coordinates.

 

[Yukterez]

I finished the simulator for rotating black holes using Kerr metric.

It takes intial positions in terms of {r0, θ0, ψ0} for radial distance and polar and spin angle, and initial velocities in terms of {v0, δ0, φ0} for total velcity and horizontal and vertical launch angle, or alternatively {vr, vθ, vψ} for the separate components.

Initial conditions:
Spin parameter: a=0.998, r0=4GM/c², initial position of the test projectile: θ0=90° (equatorial plane), horizontal launching angle: δ0=60° in south pole direction, vertical launch angle: φ0=0° (horizontal), local inital velocity: v0=0.6c

Top, side and front view relative to the spin axis of the black hole:

http://yukterez.net/org/kerr4.gif

It still needs some updates to have the same displays as the Schwarzschild-simulator which I will do when I have the time (hopefully soon).

 

[m4r35n357]

It still needs some updates to have the same displays as the Schwarzschild-simulator which I will do when I have the time (hopefully soon).

Nice. If you were to keep track of and display the max and min r values, and the maximum latitude, we could compare notes ;)

You can also compare against http://staff.science.nus.edu.sg/~phylyk/downloads/reports/sp2172_report.pdf [Broken] for particle orbits, and this for light orbits.

 

[Yukterez]

Nice. If you were to keep track of and display the max and min r values

That's next on the list, and also the local and delayed velocities as well as the time dilation factor.

and the maximum latitude

The latitude goes over all values, as seen in the right image (as well as the longitude, which is seen in the left image)

 

[m4r35n357]

That's next on the list, and also the local and delayed velocities as well as the time dilation factor.


The latitude goes over all values, as seen in the right image (as well as the longitude, which is seen in the left image)

You might want to reconsider your last statement, in the light of the top view on the left ;)
Anyway, this will become more explicit when you implement the latitude detection . . .

BTW this and this is what I mean by a maximum latitude/polar orbit.

 

[Yukterez]

You might want to reconsider your last statement, in the light of the top view on the left ;)

That's why I went online again, after thinking about it it became clear that this was an error.

 

[m4r35n357]

That's why I went online again, after thinking about it it became clear that this was an error.

Heh, I realize you know your stuff, just proving that I am awake!

 

[Yukterez]

Nice. If you were to keep track of and display the max and min r values, and the maximum latitude, we could compare notes

The most important values are now displayed below the animation (I call it "the butterfly", because with a little fantasy you can see one in the right image with the polar perspective):

http://i.imgur.com/5LuDiKl.gif

Unfortunately I did the index in german, but the latitude is the "Breitengrad θ" with 0° for the north and 180° for the south pole.
Edit: the picture embedding does not seem to work (maybe the gif is too large, the running numbers cost an extra megabyte) but if you click on the icon you should see it.

 

[m4r35n357]

why did you decide to do this project in python?

Just pulled in the Python implementation and made the algorithms identical. The Vala version is ~3.5 times as fast as the PyPy version, and at least 40 times the speed of bare Python. I scaled the number of timesteps in the run for Vala and PyPy (I didn't have the patience to wait for basic Python), and the execution times scaled by 10, so I think I've eliminated any VM startup delays from the calculation.

The difference is very significant when running on a Raspberry Pi.

 

[edguy99]

I finished the simulator for rotating black holes using Kerr metric...

Great animations, best I have seen. I do step by step animations using Newtons law calculation an any number of particles between each step (ie. force in 3d calculated, force added to velocity, velocity added to position, then re-display). Would love to convert to Kerr metric but I have a little trouble following your linked equations. Do you have any suggestions to assist in understanding how to calculate the x/y/z force on each particle for each frame using your linked equations?

Sample: http://www.animatedphysics.com/planets/planet_edit.htm