Computing trajectories in Schwarzschild spacetime

[bcrowell]

I want to produce some realistic figures showing the spatial trajectories of test particles in a Schwarzschild spacetime. For instance, I'd like to start a massive test particle at aponegricon (how often do you get to use that word!?) in an orbit that Kepler and Newton would have predicted to be elliptical, but that in fact crosses the event horizon. Basically if I can produce a list of (r,phi) pairs, I can produce such a plot.

The strategy I tried initially was to use the equation of motion expressed in terms of r and phi, as described here: http://en.wikipedia.org/wiki/Schwarzschild_geodesics#Orbits_of_test_particles . We have an equation of the form (dr/dphi)^2=f(r), which I tried to integrate numerically. However, numerical integration of this equation is kind of a hassle. Since it gives you the square of dr/dphi, you have to decide on the sign, and the code has to figure out when you've reached a turning point and its time to flip the sign. In the case of a circular or near-circular orbit, I had an especially hard time making my code smart enough to get this right. Also, because of numerical errors you can have f(r)<0, which you then have to deal with, and you have to distinguish this from the case where the solution actually doesn't continue to larger phi, e.g., a hyperbolic orbit.

The WP article also describes an exact solution using elliptic functions. This is a possible approach, although I would have to figure out if I can find an implementation of these functions in a programming language that's convenient. (I don't own a copy of Mathematica, but I do use Maxima sometimes.)

A third approach would simply be to use (t,x,y) coordinates and apply the geodesic equation. I would think that this would be more numerically well-behaved than the (r,phi) version.

Does anyone have any suggestions as to the easiest way to go about this?

 

[Mentz114]

I want to produce some realistic figures showing the spatial trajectories of test particles in a Schwarzschild spacetime. For instance, I'd like to start a massive test particle at aponegricon (how often do you get to use that word!?) in an orbit that Kepler and Newton would have predicted to be elliptical, but that in fact crosses the event horizon. Basically if I can produce a list of (r,phi) pairs, I can produce such a plot.

The strategy I tried initially was to use the equation of motion expressed in terms of r and phi, as described here: http://en.wikipedia.org/wiki/Schwarzschild_geodesics#Orbits_of_test_particles . We have an equation of the form (dr/dphi)^2=f(r), which I tried to integrate numerically. However, numerical integration of this equation is kind of a hassle. Since it gives you the square of dr/dphi, you have to decide on the sign, and the code has to figure out when you've reached a turning point and its time to flip the sign. In the case of a circular or near-circular orbit, I had an especially hard time making my code smart enough to get this right. Also, because of numerical errors you can have f(r)<0, which you then have to deal with, and you have to distinguish this from the case where the solution actually doesn't continue to larger phi, e.g., a hyperbolic orbit.

The WP article also describes an exact solution using elliptic functions. This is a possible approach, although I would have to figure out if I can find an implementation of these functions in a programming language that's convenient. (I don't own a copy of Mathematica, but I do use Maxima sometimes.)

A third approach would simply be to use (t,x,y) coordinates and apply the geodesic equation. I would think that this would be more numerically well-behaved than the (r,phi) version.

Does anyone have any suggestions as to the easiest way to go about this?

I think this does it all. Hope it helps

http://arxiv.org/pdf/1201.5611v1.pdf

 

[bcrowell]

Thanks, Lut, that's very helpful!

 

[George Jones]

See posts #16 and #31 of

https://www.physicsforums.com/threa...-schwarzschild-blackhole.347915/#post-2420613

When I wrote a computer program to to produce figures, I didn't combine the equations into an r \theta equation, but maybe something similar will work for this equation.

 

[bcrowell]

See posts #16 and #31 of

https://www.physicsforums.com/threa...-schwarzschild-blackhole.347915/#post-2420613

When I wrote a computer program to to produce figures, I didn't combine the equations into an r \theta equation, but maybe something similar will work for this equation.

I think your #16 and 31 and my #21 in that thread are all calculations for a massless test particle, but I'm interested in a massive particle now.

 

[George Jones]

The same technique applies to the case of timelike geodesics.

 

[bcrowell]

The same technique applies to the case of timelike geodesics.

Sure, it's the same geodesic equation, but the practical issues may be different. For example, in the implementation in your #16 in the other thread, you have an equation for (dr/dlambda)^2. This leads to a hassle similar to the one I described in my OP. A typical timelike geodesic will repeatedly reach turning points. You have to figure out when the sign of dr/dlambda is flipping, and you want an algorithm that doesn't misbehave for circular or near-circular orbits, and that doesn't mess up due to numerical errors. I suppose the same problem could also crop up for null geodesics near a black hole, but the typical case for a null geodesic near a less compact body would be that it only has one turning point.

Note that the paper linked to from Lut's #2 only deals with timelike geodesics.

 

[George Jones]

I have run up against this problem. The procedure in post #16 works for the timelike case; I have used it. Differentiating

$$ \left( \frac{d r}{d\lambda} \right)^2 = f \left( r \left( \lambda \right) \right)$$

leads to, after dividing both sides by $2dr / d \lambda$,

$$ \left( \frac{d^2 r}{d\lambda^2} \right) = \frac{1}{2} \frac{df}{dr}$$

Set $p = dr / d \lambda$, so that

$$ \left( \frac{d p}{d\lambda} \right) = \frac{1}{2} \frac{df}{dr}$$, so you have a coupled set of three first-order equations,

$$\begin{align}
\frac{dr}{d \lambda} &= p \\
\frac{dp}{d \lambda} &= \frac{1}{2} \frac{df}{dr} \\
\frac{d\phi}{d \lambda} &=
\end{align}$$

Compute $df/dr$ by hand, and solve the resulting system numerically. No needs to worry about the turning points; they are handled automatically.

 

[Mentz114]

I used Kostic's analysis to plot ellipsoidal bound orbits from these equations

The 4-velocity that satisfies the geodesic equation in the equatorial plane
$\begin{align} \vec{u} &= \frac{Er}{r-2m}\ \partial_t + \frac{\sqrt{\left( 2m-r\right) {L}^{2}+(E^2-1){r}^{3}+2m{r}^{2}}}{{r}^{\frac{3}{2}}} \ \partial_r + \frac{L}{{r}^{2}}\ \partial_\phi \end{align}$

from which
$\frac{dr}{d\phi} =\sqrt{ \frac{{r}^{4}\left( {E}^{2}-1\right) }{{L}^{2}}+\frac{2m{r}^{3}}{{L}^{2}}-{r}^{2}+2mr }$
It was very difficult because the parameter space $(E,L)$ is almost chaotic and all the problems that Ben mentions crop up in buckets.

The second equation can also be found in ( with $m'=1$ and $r_s=2m$ )

Gunter Scharf,
Schwarzschild Geodesics in Terms of Elliptic Functions and the Related Red Shift
Journal of Modern Physics, 2011, 2, 274-283

 

[m4r35n357]

@bcrowell, I did a side by side comparison of Newton/GR(Kerr) orbits in Javascript a couple of years ago, using Stormer-Verlet integration of the effective potentials. SV sidesteps the square root and is a symplectic integrator (excellent long-term accuracy).

BTW this all started when I used the python geodesic equations in your GR book so I'd like to give something back. Let me know if you if you are interested in my sources of info, there is a demo online at http://m4r35n357.github.io/orbits/
The sources are also on GitHub . . .

 

[George Jones]

using stormer-verlet intergration of the effective potentials. SV sidesteps the square root and is a symplectic integrator (excellent long-term accuracy).

This is because Stormer-Verlet also uses the second derivative.

 

[pervect]

I think your #16 and 31 and my #21 in that thread are all calculations for a massless test particle, but I'm interested in a massive particle now.

I am interpreting this to mean that you'd like to include the effects of gravitational radiation. I believe one has to go beyond the "test particle" idea to perform this sort of calculation, though. That's one reason I'm not sure if it's what you want - you say "test particle" but ask about a "massive" test particle, and my recollection is that test particles are presumed massless. Assuming that I've interpreted the question correctly - I'm not sure how to do this, but I suspect it gets very complicated :(. Maybe there's an easier way to do it than the full numerical GR simulations like they do for BH merger simulations, but I don't know what it is. I have a hazy recollection that the amount of gravitational radiation produced depends on the mass ratio of the infalling object to the central mass (presumed to be much larger than the mass of the infalling object), I know I read about this on a post at PF, but I can't find the post in question. It might not answer your question totally, but it would be of some help to at least estimate the accuracy of the "test particle" approximation by having some idea of the amount of perturbing gravity waves emitted.

[add]There's always the approximate relation for gravitational wave emission of $$\frac{1}{5} \left( \frac{d^3 Q} {dt^3} \right) ^2 \frac{c^5}{G}$$ where Q is the magnitude of the quadrapole moment, but that's not what I was thinking of :(.

 

[PeterDonis]

test particles are presumed massless

Not in the sense bcrowell is using "massless"--test particles don't have to be lightlike. What is assumed is that they do not affect the spacetime geometry, i.e., their contribution to the stress-energy tensor is negligible.

 

[Mentz114]

I did the simplest thing and it seems to work. Integrating $\dot{r}$ and $\dot{\phi}$ from the 4-velocity in my earlier post using this Maxima script. The values of $E$ and $L$ are guesses based on Kostic's analysis.

The script is

/********************************
Schwarzschild geodesic plot
***************************/
kill(all);
drdphi:sqrt( r^4*(E^2-1)/L^2+2*m*r^3/L^2-r^2+2*m*r)$
rdot:sqrt( ((2*m-r)*L^2+ (E^2-1)*r^3+2*m*r^2)/r^3)$
pdot:L/r^2$

rdotval( pr,pm,pE,pL) := block( [res],res:subst(pr,r,rdot),
res:(subst(pm,m,res)),
res:(subst(pL,L,res)),
res:(subst(pE,E,res)),
res:ev(res,nouns),
return(res))$

Lval:3.0$
numer:true$
numits:100;
rn:500$
phin:0$
dtau:6;
for p thru numits do [
phin:phin+(Lval/rn^2)*dtau,
rn:rn-rdotval(rn,1,1.2,Lval)*dtau,
print( phin,",",rn,",",rn*cos(phin),",",rn*sin(phin))
]$

[Edut] I deleted the picture because it is wrong - sorry. I'll repost it when it's fixed

 

[George Jones]

I did the simplest thing and it seems to work. Integrating $\dot{r}$ and $\dot{\phi}$ from the 4-velocity in my earlier post using this Maxima script. The values of $E$ and $L$ are guesses based on Kostic's analysis.

But you you don't seem to have any turning points (i.e., dr/d phi doesn't seem to change sign on your plot), which is the situation about which Ben is concerned.

 

[Mentz114]

But you you don't seem to have any turning points (i.e., dr/d phi doesn't seem to change sign on your plot), which is the situation about which Ben is concerned.

This is true. dr/d\phi falls monotonically in the orbit I plotted. But I'm integrating velocities to get displacement so ( I hope) avoiding the vagaries of dr/d\phi.

[Edit] I made an error in the plotting ( surprise, surprise) and I'll repost the picture when its fixed.

 

[bcrowell]

@bcrowell, I did a side by side comparison of Newton/GR(Kerr) orbits in Javascript a couple of years ago, using Stormer-Verlet integration of the effective potentials. SV sidesteps the square root and is a symplectic integrator (excellent long-term accuracy).

BTW this all started when I used the python geodesic equations in your GR book so I'd like to give something back. Let me know if you if you are interested in my sources of info, there is a demo online at http://m4r35n357.github.io/orbits/
The sources are also on GitHub . . .

That is just super super cool! It's kind of hypnotic watching the orbits. Can you point us to the github url?

 

[pervect]

Not in the sense bcrowell is using "massless"--test particles don't have to be lightlike. What is assumed is that they do not affect the spacetime geometry, i.e., their contribution to the stress-energy tensor is negligible.

As long as the mass is sufficiently small, the particle will follow a geodesic

Not in the sense bcrowell is using "massless"--test particles don't have to be lightlike. What is assumed is that they do not affect the spacetime geometry, i.e., their contribution to the stress-energy tensor is negligible.

Not what I was concerned about at all. I was hoping Ben would clarify. Geodesics are fine for particles of negligible mass (not zero, as you point out). i.e test particles. But Ben seemed concerned about the effects of particle mass on the orbits. Since there isn't any effect if the mass is negligible, I thought maybe he might be asking about the non-test particle case, to determine the effects of non-negligible masses when you go beyond the test particle approximation. But I'm not sure if that's what he was asking, really.

 

[m4r35n357]

That is just super super cool! It's kind of hypnotic watching the orbits. Can you point us to the github url?

Thanks for that!
https://github.com/m4r35n357/JsBlackHole should get you there. The physics is supposed to be encapsulated in models.js. Since you've replied, I'll dig out the references for the effective potential equations that I used . . .
http://eagle.phys.utk.edu/guidry/astro421/lectures/lecture490_ch13.pdf

Incidentally, I check all my work against this program: http://stuleja.org/grorbits/, are you aware of it?

 

[Mentz114]

@bcrowell, I did a side by side comparison of Newton/GR(Kerr) orbits in Javascript a couple of years ago, using Stormer-Verlet integration of the effective potentials. SV sidesteps the square root and is a symplectic integrator (excellent long-term accuracy).

BTW this all started when I used the python geodesic equations in your GR book so I'd like to give something back. Let me know if you if you are interested in my sources of info, there is a demo online at http://m4r35n357.github.io/orbits/
The sources are also on GitHub . . .

That is great. Nice work.

I never got past the periapsis turning point problem ...

 

[PeterDonis]

As long as the mass is sufficiently small, the particle will follow a geodesic

Yes, but it could be a timelike geodesic or a lightlike geodesic, so the term "mass" isn't really appropriate here since it invites confusion about what kind of "mass" you are talking about. It can't be invariant mass because photons have zero invariant mass but they can still potentially affect the spacetime geometry; but if you mean relativistic mass, that isn't what determines whether spacetime geometry is affected either (and that term isn't usually used to describe photons anyway). What determines that is the particle's stress-energy; that is what must be assumed to be negligible for a test particle.

Ben seemed concerned about the effects of particle mass on the orbits.

I think when he said "massive", he meant "timelike", meaning he was looking at the behavior of timelike geodesics as opposed to lightlike geodesics--i.e., he was referring to invariant mass, not stress-energy. The effective potential and the behavior of the orbits are qualitatively different for the two cases. Ben can clarify, of course, if I've misinterpreted.

 

[m4r35n357]

Thanks guys, appreciated. BTW I used the same hammer (symplectic integration) in Python to get a working 4D Kerr geodesic generator (for light/particles) and a Newtonian N-Body simulator. Both are also on Github . . .
It's a damn good hammer!

 

[rollingstein]

The strategy I tried initially was to use the equation of motion expressed in terms of r and phi, as described here: http://en.wikipedia.org/wiki/Schwarzschild_geodesics#Orbits_of_test_particles . We have an equation of the form (dr/dphi)^2=f(r), which I tried to integrate numerically.

Naive question: What boundary condition do I use when solving this ODE? Or can one start from essentially any given (r, phi) combination?

 

[bcrowell]

Naive question: What boundary condition do I use when solving this ODE? Or can one start from essentially any given (r, phi) combination?

You fix an initial position (r,phi), and you also fix the values of the conserved quantities E and L (which have to be consistent with the initial position). You also have to fix the initial direction of motion (sign of dr/dphi, if it's not zero).

 

[George Jones]

Equivalently, in my simulation, the user launches a projectile from a hovering platform. The use inputs r for the hovering a platform and the initial velocity of the projectile with respect to (an orthonormal frame for) the platform. The user inputs speed and angle with respect to "horizontal" for the platform. The program takes the initial phi to be zero, and uses these inputs to calculate E and L.

 

[rollingstein]

You fix an initial position (r,phi), and you also fix the values of the conserved quantities E and L (which have to be consistent with the initial position). You also have to fix the initial direction of motion (sign of dr/dphi, if it's not zero).

This being a 1st order ODE you ought to need only a single Initial Condition, correct? So why are two independant conditions fixable here? Both (r,phi) & the sign of dr/dphi?

 

[Jimster41]

I would love to have some more explication of what we (or I as a layperson) can learn from this. It's very interesting/cool. I gather we are able to see a divergence from a Newtonian Model of a particle orbit and a GR model (or two GR models)? Can someone say more about how why the Schwarzschild (BH) case is special, why the modeler chose to look at it in particular in this case, how it does or doesn't generalize. And, maybe related, how the discussion of Gravitational Radiation bears. Does that mean gravitational waves?

 

[Jimster41]

[add]There's always the approximate relation for gravitational wave emission of

15(d3Qdt3)2c5G​

\frac{1}{5} \left( \frac{d^3 Q} {dt^3} \right) ^2 \frac{c^5}{G} where Q is the magnitude of the quadrapole moment, but that's not what I was thinking of :(.

Any translation, if it were possible, of what the "quadrupole moment" means, and how it's related to Gravitational Radiation would be very much appreciated. I'm studying the wiki, but it's over my head presently. Edit: Scratch that. I'm able to follow the formalism. Still could use (would love to have) some more context for the modelling and the relevance of Gravitational Radiation.

What does the plot at bottom left show?

 

[Jimster41]

Question about copyright.

@bcrowell, I did a side by side comparison of Newton/GR(Kerr) orbits in Javascript a couple of years ago, using Stormer-Verlet integration of the effective potentials. SV sidesteps the square root and is a symplectic integrator (excellent long-term accuracy).

BTW this all started when I used the python geodesic equations in your GR book so I'd like to give something back. Let me know if you if you are interested in my sources of info, there is a demo online at http://m4r35n357.github.io/orbits/
The sources are also on GitHub . . .

I downloaded the source for this simulator from github. Is this code I can play with? I assumed some copyright, yours, by default. But then I've never heard of github, and got the idea it was considered open source. I'm pretty good with C#, and I code with it often (in application/algorithm contexts, not as a pure developer), but I'm pretty new(100%noob) to java. I made some changes and was having no luck saving them. I assumed because I am a doof with java. But then I thought I better check to be sure that I wasn't misunderstanding the copyright context (and my saves weren't working because it was locked).

 

[bcrowell]

Question about copyright.
I downloaded the source for this simulator from github. Is this code I can play with? I assumed some copyright, yours, by default. But then I've never heard of github, and got the idea it was considered open source. I'm pretty good with C#, and I code with it often (in application/algorithm contexts, not as a pure developer), but I'm pretty new(100%noob) to java. I made some changes and was having no luck saving them. I assumed because I am a doof with java. But then I thought I better check to be sure that I wasn't misunderstanding the copyright context (and my saves weren't working because it was locked).

He has an open-source license in a comment in the source code: https://github.com/m4r35n357/JsBlackHole/blob/master/display.js

 

[Jimster41]

I saw that but wasn't sure if that was for the orbit model code or just required boilerplate for distributing anything running Java.

 

[bcrowell]

I saw that but wasn't sure if that was for the orbit model code or just required boilerplate for distributing anything running Java.

It's an open-source license. There is no such license required for java code in general.

 

[Mentz114]

Question about copyright.
I downloaded the source for this simulator from github. Is this code I can play with? I assumed some copyright, yours, by default. But then I've never heard of github, and got the idea it was considered open source. I'm pretty good with C#, and I code with it often (in application/algorithm contexts, not as a pure developer), but I'm pretty new(100%noob) to java. I made some changes and was having no luck saving them. I assumed because I am a doof with java. But then I thought I better check to be sure that I wasn't misunderstanding the copyright context (and my saves weren't working because it was locked).

Actually it is JavaScript, not Java. JavaScript can have objects see here for instance https://developer.mozilla.org/en-US/docs/Web/JavaScript/Introduction_to_Object-Oriented_JavaScript
I downloaded the zip file from github and the code is working locally in FireFox. I'm still impressed.

 

[Jimster41]

Actually it is JavaScript, not Java. JavaScript can have objects see here for instance https://developer.mozilla.org/en-US/docs/Web/JavaScript/Introduction_to_Object-Oriented_JavaScript
I downloaded the zip file from github abd the code is working locally in FireFox. I'm still impressed.

I was able to run it no problem after the download. Just trying to figure out how to monkey around with the pieces. Thank s for the ref. I've got a good book and I'm pretty familiar with OO architecture, just not HTML, or CSS, or .JS. I've been working in a proprietary environment for... Too long. Its all drag and drop, with bits of custom C#. Lost my skillz.

 

[Mentz114]

I was able to run it no problem after the download. Just trying to figure out how to monkey around with the pieces. Thank s for the ref. I've got a good book and I'm pretty familiar with OO architecture, just not HTML, or CSS, or .JS. I've been working in a proprietary environment for... Too long. Its all drag and drop, with bits of custom C#. Lost my skillz.

OK, you'll soon get familiar with Javascript. It's easy and pretty 'fiendly'.

I've changed some colours and I'm going to lay out the screen differently and add a trace list to log (say) the differences between succesive $\phi_{min}$ values to get the precession of the major axis.

All good fun ...