# Computing trajectories in Schwarzschild spacetime

1. Apr 13, 2015

### bcrowell

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

2. Apr 13, 2015

### Mentz114

I think this does it all. Hope it helps

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

3. Apr 13, 2015

### bcrowell

Staff Emeritus

4. Apr 13, 2015

### George Jones

Staff Emeritus
5. Apr 13, 2015

### bcrowell

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

6. Apr 13, 2015

### George Jones

Staff Emeritus
The same technique applies to the case of timelike geodesics.

7. Apr 13, 2015

### bcrowell

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

8. Apr 13, 2015

### George Jones

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

9. Apr 13, 2015

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

10. Apr 14, 2015

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

Last edited: Apr 14, 2015
11. Apr 14, 2015

### George Jones

Staff Emeritus
12. Apr 14, 2015

### pervect

Staff Emeritus
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 :(.

Last edited: Apr 14, 2015
13. Apr 14, 2015

### Staff: Mentor

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.

14. Apr 14, 2015

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

Last edited: Apr 14, 2015
15. Apr 14, 2015

### George Jones

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

16. Apr 14, 2015

### Mentz114

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.

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

Last edited: Apr 14, 2015
17. Apr 15, 2015

### bcrowell

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

18. Apr 15, 2015

### pervect

Staff Emeritus
As long as the mass is sufficiently small, the particle will follow a geodesic
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.

19. Apr 15, 2015

### m4r35n357

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?

Last edited: Apr 15, 2015
20. Apr 15, 2015

### Mentz114

That is great. Nice work.

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

21. Apr 15, 2015

### Staff: Mentor

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.

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.

22. Apr 15, 2015

### 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!

23. Apr 16, 2015

### rollingstein

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

24. Apr 16, 2015

### bcrowell

Staff Emeritus
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).

25. Apr 16, 2015

### George Jones

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