I Explore Black Hole Orbits with My Kerr Orbit Simulator on YouTube

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

 

[Yukterez]

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?

If you use Mathematica you can automatically differentiate r'[t], θ'[t] and ψ'[t] (which are defined in the differential equation) with

Evaluate[r''[t] /. sol][[1]]
Evaluate[θ''[t] /. sol][[1]]
Evaluate[ψ''[t] /. sol][[1]]

to get the second (proper) time derivative of the motion. The transformation from r, θ, ψ components to x, y, z components is

x = r Sin[θ] Cos[ψ])
y = r Sin[θ] Sin[ψ]
z = r Cos[θ]

 

[pervect]

If you use Mathematica you can automatically differentiate r'[t], θ'[t] and ψ'[t] (which are defined in the differential equation) with

Evaluate[r''[t] /. sol][[1]]
Evaluate[θ''[t] /. sol][[1]]
Evaluate[ψ''[t] /. sol][[1]]

to get the second (proper) time derivative of the motion. The transformation from r, θ, ψ components to x, y, z components is

x = r Sin[θ] Cos[ψ])
y = r Sin[θ] Sin[ψ]
z = r Cos[θ]

I assume r, theta, and psi are Schwarzschild coordinates? The difficulty with constructing x,y, and z in the manner you did above from r, theta, and psi is that the speed of light is not isotropic in Schwarzschild coordinates, i.e. it depends on direction. Thus with these relationships, you have the speed of light being different in the x direction, the y direction, and the z direction. This means that the distance scales are not the same - if you consider the distance light travels in one "tick" of proper time (or, if you prefer, one "tick" of coordinate time), the distance the light moves depends on the direction, i.e. the distance is different in the x, y, and z directions.

This is rather incompatible with the standard SI definition of the meter if you attempt to interpret your coordinate values as having physical significance (i.e. representing distances). The proportionality between a change in the coordinate value and the change in distance depends on which coordinate you pick.

The usual solution to this is to use an alternate coordinate system, the isotropic coordinate system, which you can find the metric for in The "Alternate coordinate" section of https://en.wikipedia.org/wiki/Schwarzschild_metric. It's rather painful to compute things in isotropic coordinates, though it can be done, the equations are are lot more complicated. It would probably be easiest to compute the answer in Schwarzschild coordinates, and convert the answer into isotropic coordinates. Then you could convert that answer into x,y, and z, and the results would be much more physically significant.

I don't recall the transformation equations offhand, however, and I didn't notice them in the wiki article I quoted.

 

[Yukterez]

Thus with these relationships, you have the speed of light being different in the x direction, the y direction, and the z direction.

In Schwarzschild the factor for the transversal direction is √(1-rs/r) and for the radial direction (1-rs/r) without the square (because you have time dilation and radial length expansion of the same magnitude, also see this link). In Kerr you replace the Schwarzschild radius rs=2 with rk=1+√(1-a²Cos²[θ]) where a is the spin parameter.

 

[Yukterez]

I have a little trouble following your linked equations

I see I have a bug in the .txt-File code and mixed up some r0 and r(t) in the formula for the digit velocity display (the simulation is ok, but the digit output below the animation needs a repair). I'll fix that tomorrow, until then you can use the equations of motion from my animation code, but if you want to display the numbers transform from r', θ' & ψ' to vr, vθ & vψ like in my last posting above, where all r are r(t).

 

[Yukterez]

The bug in the code is fixed now, the transformation from proper to coordinate velocity should give the right numers now. If you use the kerr.txt force refresh (ctrl f5) to make sure to get the updated version.

I'll put the transformation in Latex (with G, M and c set to 1):
to transform from the derivatives in the differential equation to the local velocity components use

$\dot{r} = \frac{v_r (1-r_k/r)^{(3/2)}}{\sqrt{1-v^2}}$

$\dot{\theta} = \frac{v_{\theta}}{r \sqrt{1-v^2}}$

$\dot{\psi} = \frac{v_{\psi} }{r \sqrt{1-v^2}}$

The $v_{\psi}$ component is observed faster by

$d_{\psi} = \pm \frac{2 \alpha r^2 }{ \Sigma}$

with

$r_k = \sqrt{1-a^2 \cos^2 (\theta )}+1$

because of frame dragging, where

$\Sigma =\left(\alpha ^2+r^2\right)^2-\alpha ^2 \sin ^2(\theta ) \Delta$

and

$\Delta =r^{2}-2r+\alpha ^{2}$

After solving for $v_r$, $v_{\theta}$, $v_{\psi}$ and $d_{v_{\psi}}$ you can use Pythagoras for the total velocity $v$.

 

[Drakkith]

I'll post again and ask a moderator to delete the previous post:

Done.

 

[Yukterez]

I do step by step animations using Newtons law calculation an any number of particles between each step
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?

Now I see what you mean, but that needs some further explanation. You can compute the force acting on the test particle from the velocity derivative (after substracting the frame drag), but you don't use the force to compute the velocity. Normally with Newton or Schwarzschild you plug in your function for the accelerations r''(t), ψ''(t) & θ''(t) to equate the motion and get your velocities r'(t), ψ'(t) & θ'(t) and positions r(t), ψ(t) & θ(t). But in Kerr even when a particle has 0 velocity, the velocity in the equatorial direction - ψ'(t) - is not zero from begin with. So your test particle already starts with a velocity even when it has none (because of the frame dragging). Therefore under Kerr the equation of motion is not second order but first order, and therefore a little more complicated because you have to write the equation in terms of conserved quantities that go into the definions of r'(t), ψ'(t) & θ'(t). Therefore it is not possible to plug in Kerr corrections into a Newtonian simulator (with Schwarzschild on the other hand this is no problem, the equation of motion only differs by a small additional term). If you want to use Kerr you have to switch from second order derivative to first order derivative, the force is only on the output side but not on the input.

 

[Yukterez]

I had a numeric error in the last posting, because it was too late to edit I deleted and repost

This trajectory is a useful example to explain the numeric display below the animation, with focus on the last three numbers on the bottom right:

At radial coordinate r=3GM/c² we launch a test particle to a prograde orbit spiraling into the ergosphere of a rotating Kerr black hole.

When the test particle approaches the horizon at 1+√[1-a²],
it has a local velocity of almost c relative to a local probe which is locally at rest,
a delayed velocity of almost 0 relative to the same probe when described from the perspective of the observer at infinity, and
an observed velocity of 0.45c (that is the frame dragging velocity at the latitude where the particle nears the horizon, so we observe the locally stationary probe and the locally almost light speed fast particle to corotate with almost the same 0.45c).

First in proper time steps:

Now the same situation in coordinate time steps (the interval is as above 1/8GM/c³):

Comparison of proper speed versus coordinate velocity as viewed from different observers (left: particle, right: coordinate observer):

Initial settings: a=0.9, r0=3GM/c², v0=√(1/r0)=1/√(3)c, ψ0=0, θ0=π/2=90°, φ0=0, δ0=π/5=36°

The red tail length is 1/4 GM/c³ of proper time (so before the particle enters the horizon it gets infinitely long because it whirls up while freezing and corotating from the perspective of the coordinate observer). As one can see, the test particle plunges into the horizon in a finite proper time. But it also makes an infinite amount of revolutions with the black hole before it does so.

I heard that there is a coordinate system where one can transform this infinty away and follow the particle on its way behind the horizon, but as far as I am with Kerr right now I can only simulate until there and then I get an infinity in the spin rate when differentiating coordinates by proper time at the inner horizon (the particle proper time freezes while it is frame dragged around the horizon and corotating with a constant velocity, in this case 0.45c)

Code.txt

 

[edguy99]

At radial coordinate r=3GM/c² we launch a test particle to a prograde orbit spiraling into the ergosphere of a rotating Kerr black hole.

Thanks for the post. I would like to see the simulation of a pair of black holes orbiting one another (or a star orbiting a small black hole) where the mass of one object is affecting the position of the other object and you can expand it to more objects if you want (ie. http://www.animatedphysics.com/planets/moon_orbit.htm this animation calculates both the effect of Newtonian gravity of the Earth's pull on the moon as well as the moons pull on the earth, where both the Earth's position and the moons position changes between each frame due to the gravity pull of the moon). Is this doable with your equations? Specifically, I don't see anything that would represent the mass of the particle around the black hole? I could be not reading German properly of just significantly confused.

The only way I see to attempt to model a pair of black holes is through x/y/z timesteps, with calculations and movement of each object between timesteps. Perhaps there is an easier way to do this?

 

[Yukterez]

Is this doable with your equations?

With my Newtonian equations it is doable, but only in the nonrelativistic limit, see here. With Kerr's metric it is not doable, since

The Kerr metric or Kerr geometry describes the geometry of empty spacetime around a rotating uncharged axially-symmetric black hole with a spherical event horizon.

So the mass of the test particle is assumed to be small in Kerr metric, just like in Schwarzschild metric. If you want to simulate full relativistic n-body you need to go for the BSSN formalism but I doubt that you can do that on a home computer, as far as I know supercomputers are busy for weeks if not months with that.

 

[Yukterez]

The old simulator had a bug, as one can see in the animations above the ergosphere was not displayed correctly in it's pumpkin-shape but as an ellipsoid. That was fixed in an update; also the numerical display was extended and the plot now shows not only the outer but also the inner horizon and ergosphere. So use the new code instead of the old one. Example: plunge orbit with v0=vz=0.975c from r0=3GM/c², a=0.998:

The initial conditions can be fed in in terms of the local 3-velocity (relative to a ZAMO) either vx, vy, vz or equivalently vr, vφ, vθ or in terms of the total Energy and angular momentum components E, Lz, pθ0 or alternatively Q (the carter constant). Index for the numerical display: