Schwarzschild Orbits in Cartesian coordinates

**CarlB** · Jul24-06, 03:43 AM

My Java applet gravity simulator http://www.gaugegravity.com/testappl...etGravity.html
draws beautiful orbits, however the GR simulation is very badly broken as one can tell when comparing it with Newton at long distances. The source code is at:
http://www.gaugegravity.com/testapplet/SwGrav_Top.java

Test bodies are required to all lie in a plane. The Newtonian physics is produced by writing the x and y accelerations in terms of the phase space parameters of the test body, that is, in terms of its (x,y) position and its (Vx,Vy) velocity. With Newton, the acceleration has no dependence on velocity of course.

The GR "physics" is produced by a complicated calculation that follows the artificial potential that is discussed with respect to the Schwarzschild orbits in MTW and in various places all over the web. These methods use energy and angular momentum as constants of the motion to describe the orbits but they are less than optimal when used in an integration form.

I've seen a computer method for computing GR orbits that should work better, but I am hesitating to implement it because it is based on polar coordinates. This means that the calculation has to use trig functions to convert back and forth and these are slow.

So I'm wondering how hard it is to put the Schwarzschild orbits into phase space form in Cartesian coordinates. Here's the basic plan:

(1) Write the Schwarzschild metric in Cartesian coordinates.

(2) Write the proper length of a path as an integral over coordinate time.

(3) Vary the path and use the Euler-Lagarange equation to determine a pair of 2nd order differential equations that the orbits solve.

(4) Find

and

from the two 2nd order DEs.

Now for the poll: Do you think that this can be done reasonably easily? Part of the reason I am working on it is that I want to get a better idea on how these gravity theories differ one from another. I will post calculations as I make them, if convenient.

Carl

**CarlB** · Jul24-06, 03:51 AM

Step (1), write the Schwarzschild metric in Cartesian coordinates.

$LaTeX Code: x = r \\cos(\\phi)$
$LaTeX Code: y = r \\sin(\\phi)$

$LaTeX Code: dx = \\cos(\\phi) dr - r\\sin(\\phi) d\\phi$
$LaTeX Code: dy = \\sin(\\phi) dr + r\\cos(\\phi) d\\phi$

$LaTeX Code: dr = \\cos(\\phi) dx + \\sin(\\phi) dy$
$LaTeX Code: rd\\phi = -\\sin(\\phi) dx + \\cos(\\phi) dy$

$LaTeX Code: dr^2 = \\cos^2(\\phi) dx^2 + \\sin^2(\\phi) dy^2 + \\sin(2\\phi) dx\\;dy$
$LaTeX Code: r^2d\\phi^2 = \\sin^2(\\phi) dx^2 + \\cos^2(\\phi) dy^2 - \\sin(2\\phi) dx\\;dy$

$LaTeX Code: dr^2 = (x^2dx^2 + y^2dy^2 + 2xy dx\\;dy)/r^2$
$LaTeX Code: r^2d\\phi^2 = (y^2dx^2 + x^2dy^2 - 2xy dx\\;dy)/r^2$

$LaTeX Code: ds^2 = -\\frac{(r-2)}{r}dt^2 + \\frac{r}{r-2}dr^2 + r^2d\\phi^2$

$LaTeX Code: ds^2 = -\\frac{(r-2)}{r}dt^2 + \\frac{1}{r(r-2)}(x^2dx^2 + y^2dy^2 + 2xy dx\\;dy) + \\frac{1}{r^2}(y^2dx^2 + x^2dy^2 - 2xy dx\\;dy)$

$LaTeX Code: ds^2 = -\\frac{(r-2)}{r}dt^2 + \\frac{1}{r^2(r-2)}((rx^2+ ry^2 -2y^2)dx^2 + (ry^2+rx^2-2x^2)dy^2 + (2rxy -2rxy +4rxy) dx\\;dy)$

$LaTeX Code: ds^2 = -\\frac{(r-2)}{r}dt^2 + \\frac{1}{r^2(r-2)}((r^3 -2y^2)dx^2 + (r^3-2x^2)dy^2 + (4rxy) dx\\;dy)$

$LaTeX Code: ds^2 = -\\frac{(r-2)}{r}dt^2 + \\frac{(r^3-2y^2)}{(r^3-2r^2)}dx^2 + \\frac{(r^3-2x^2)}{(r^3-2r^2)}dy^2 +\\frac{4rxy}{(r^3-2r^2)}dx\\;dy$

**Jorrie** · Jul24-06, 05:02 AM

Originally Posted by CarlB

...
Now for the poll: Do you think that this can be done reasonably easily? Part of the reason I am working on it is that I want to get a better idea on how these gravity theories differ one from another. I will post calculations as I make them, if convenient.
Carl

Interesting excercise Carl, but how do you think having the orbital equation in Cartesian form (as against polar form) will help you get a "better idea on how these gravity theories differ one from another."?

**pervect** · Jul24-06, 12:19 PM

There isn't any such thing as "Cartesian coordinates" for a black hole. Cartesian coordinates work only for flat space-times. The space-time around a black hole is not flat. In particular, r is not a radial coordinate.

The closest you can come to "cartesian" coordiantes are isotropic coordinates.

Howver, from experience, I can tell you it will be far far easier to get your program working for Schwarzschild coordinates than isotropic cooridinates.

As I mentioned before, you can just use the geodesic equations to express your set of equations as a second order system that look pretty much like the second order set you have already for Newton's equations.

This is the "geodesic" equation

http://en.wikipedia.org/wiki/Solving...esic_equations

(Wiki somewhat confusingly uses "t" instead of "tau" for the variables.

The necessary Christoffel symbols are simplest for the schwarzschild metric, they are

Code:


                        r                 m
                      CC   [r r] = - -----------
                                     (r - 2 m) r


                         theta
                       CC       [theta r] = 1/r


                           phi
                         CC     [phi r] = 1/r


                         t               m
                       CC   [t r] = -----------
                                    (r - 2 m) r


                         theta
                       CC       [r theta] = 1/r


                      r
                    CC   [theta theta] = -r + 2 m


                     phi                cos(theta)
                   CC     [phi theta] = ----------
                                        sin(theta)


                           phi
                         CC     [r phi] = 1/r


                     phi                cos(theta)
                   CC     [theta phi] = ----------
                                        sin(theta)


                 r                                   2
               CC   [phi phi] = -(r - 2 m) sin(theta)


               theta
             CC       [phi phi] = -sin(theta) cos(theta)


                         t               m
                       CC   [r t] = -----------
                                    (r - 2 m) r


                         r          (r - 2 m) m
                       CC   [t t] = -----------
                                         3
                                        r

As I also mentioned, you should be able to obtain the same set of geodesic equations from the conserved quantites, as well, by replacing "conserve quantity = constant" with "d/dtau (conserved quantity) = 0"

**pervect** · Jul24-06, 01:00 PM

Let's try to make things really, really, really simple.

theta will always be 90 degrees. Therfore d theta / dtau will be zero.

We then have for our variables

r(tau), t(tau), phi(tau).

The second order geodesic equations of motion (the geodesic equations) will therfore be

$LaTeX Code: \\frac{d^2 r}{d tau^2} + \\Gamma^r{}_{rr} \\left( \\frac{dr}{dtau} \\right) ^2 + \\Gamma^r{}_{tt} \\left( \\frac{dt}{d tau} \\right)^2 + \\Gamma^r{}_{\\phi\\phi} \\left( \\frac{d \\phi}{d tau} \\right)^2 = 0 $

$LaTeX Code: \\frac{d^2 t}{d tau^2} + \\left( \\Gamma^t{}_{tr} + \\Gamma^t{}_{rt} \\right) \\left( \\frac{dr}{d tau} \\right) \\left( \\frac{dt}{dtau} \\right) = 0 $

$LaTeX Code: \\frac{d^2 \\phi}{d tau^2} + \\left( \\Gamma^{\\phi}{}_{r\\phi} + \\Gamma^{\\phi}{}_{\\phi r} \\right) \\left( \\frac{dr}{d tau} \\right) \\left( \\frac{d \\phi} {d tau} \\right) = 0 $

See the previous post for the values for the Christoffel symbols.

Note for instance that since r^2 $LaTeX Code: {d \\phi}/ {d tau}$ = constant, when differentiated, gives by the chain rule

$LaTeX Code: 2 r \\frac{dr}{d\\tau} \\frac{d \\phi}{d tau} + r^2 \\frac{d^2 \\phi} {d tau^2} = 0 $

and that this is equivalent to the last equation above when you substitute in the correct value for the Christoffel symbol.

The conseved energy when differentiated will give one of the other geodesic equations (the one for time).

The last geodesic equation is a consequence of the others plus the equation for the metric.

**Jorrie** · Jul24-06, 01:38 PM

Originally Posted by pervect

There isn't any such thing as "Cartesian coordinates" for a black hole. Cartesian coordinates work only for flat space-times. The space-time around a black hole is not flat. In particular, r is not a radial coordinate.

The closest you can come to "Cartesian" coordinates are isotropic coordinates.

Are you not a bit harsh on CarlB? I understood that he simply wants to express Schwarzschild coordinates in terms of x,y,z rather than in terms of r,phi,theta. Nothing fundamentally wrong with that, I think! Maybe his (Cartesian) choice of words was wrong, but let’s not get too semantic!

Nevertheless, I agree that we should rather stick to the tried and tested procedure of treating black holes.

**CarlB** · Jul24-06, 01:48 PM

pervect, one of the things I'm trying to get away from is differential equations written in proper time. It makes for a more complicated simulation and also makes it a lot easier for the simulation to fail in various ways. So I'm continuing with the process of writing the equations in Cartesian coordinates. I will allow you to continue calling them isotropic coordinates if that's what you want to do.

Carl

**pervect** · Jul24-06, 03:12 PM

Carl (and Jorrie, too) I think you've totally missed the point.

The 'r' in the Scwarzschild coordinate system is not the same 'r' in the isotropic coordinate system.

In fact, one can see that r = r'(1+2m/r')^2 from the webpage, where r' is the isotropic radial coordinate and r is the Schwarzschild r coordinate.

Neither of these 'r' coordinates can be interpreted as 'radial distance'. The isotropic coordinates, though, define a set of coordinates where light cones will apear to be circles.

Thus assuming that r^2 = x^2 + y^2 does not make any physical sense. This is assumed implicitly when one writes x = r cos(phi) y = r sin(phi), one assumes that x^2 + y^2 = r^2. One can think of it as a purely formal mathematical expression, I suppose.

As such a formal exercise, unless Carl has made an algebraic error he's got a new and totally non-standard metric that doesn't even have a name in terms of new and toatally non-standard coordiantes that don't have names either (except what Carl has christened them). They definitely shouldn't be confused with Cartesian coordinates, which is what Carl has chosen to call his variables.

One can easily see that Carl's metric is not the same as the isotropic metric either by comparing the line elements.

Aside from the fact that it is non-standard, and thus will be hard to check solutions for, even deriving the correct solutions is going to be a giant pain. Re-writing a simple metric in a more complicated form is not going to get anywyere as far as generating a solution for the metric.

If the problem is with the factor of 'tau' in the standard equations, just eliminate that tau. This shouldn't be too hard since

(1-2m/r) dt/dtau = E

The most annoying thing is going to be this

If we write r as a function of t : r = rr(t) then

r(tau) = rr(t(tau))

and dr/dtau = drr/dt dt/dtau

but by the chain rule

$LaTeX Code: \\frac{d^2 r} {d tau^2} = \\frac{d^2 rr}{d t^2} \\left( \\frac{dt}{d tau} \\right) ^2 + \\frac{dr}{dt} \\frac{d^2 t}{d tau^2} $

Still, a little bit of algebra should eliminate tau as a variable if this is really what is desired.

Solving for phi shouldn't be terriby difficult, since

r^2 dphi/dtau =L

therfore

r^2 dphi/dt * dt/dtau = L

And we've already noted that dt/dtau is easy to eliminate

**pervect** · Jul24-06, 06:28 PM

Unless I've made a mistake (quite possible, even with computer assistance) the equations with proper time eliminated for the Schwarzschld metric are just

r(t), phi(t)

d phi/dt = K*(1-2*m/r) / r^2

K is a constant of motion, equal to "L/E" in the old notation. Knowing dphi/dt and r at one point will determine the value of K

The second order equation for r(t) is

$LaTeX Code: \\frac{d^2 r}{dt^2} = \\frac{2\\,m}{r^2 \\left( 1 - \\frac{2m}{r} \\right) } \\left( \\frac{dr}{dt} \\right) ^2 - \\left( 1-{\\frac {2 m}{r}} \\right) \\frac{m}{r^2}+ \\left( 1-{\\frac {2m}{r}} \\right) ^{3}\\frac{{K}^{2}}{r^3} $

It looks sensible at first glance, dr/dt is forced to zero, however I have not tested it.

Compare this to the newtonian equations in polar coordinates (r, phi)

r^2 * d phi/dt = K

where K = L/m

d^2 r / dt^2 = -m/r^2 + r (dphi/dt)^2 = -m/r^2 + K^2 / r^3

we can also see that the Schwarzschild equations reduce to the Newtonian equations when dr/dt << 1 and r >> 2m as they should, another sign that they are probably right.

**pervect** · Jul24-06, 06:38 PM

The Christoffel symbols for the Painleve metric are:

Code:

                                   1/2      1/2
                       T          2    (M/r)    M
                     CC   [T T] = ---------------
                                         2
                                        r


                         r          (r - 2 M) M
                       CC   [T T] = -----------
                                         3
                                        r


                            T           M
                          CC   [r T] = ----
                                         2
                                        r


                                    1/2      1/2
                      r            2    (M/r)    M
                    CC   [r T] = - ---------------
                                          2
                                         r


                            T           M
                          CC   [T r] = ----
                                         2
                                        r


                                    1/2      1/2
                      r            2    (M/r)    M
                    CC   [T r] = - ---------------
                                          2
                                         r


                                       1/2
                        T             2    M
                      CC   [r r] = -------------
                                          1/2  2
                                   2 (M/r)    r


                           r             M
                         CC   [r r] = - ----
                                          2
                                         r


                         theta
                       CC       [theta r] = 1/r


                           phi
                         CC     [phi r] = 1/r


                         theta
                       CC       [r theta] = 1/r


                  T                    1/2      1/2
                CC   [theta theta] = -2    (M/r)    r


                      r
                    CC   [theta theta] = -r + 2 M


                     phi                cos(theta)
                   CC     [phi theta] = ----------
                                        sin(theta)


                           phi
                         CC     [r phi] = 1/r


                     phi                cos(theta)
                   CC     [theta phi] = ----------
                                        sin(theta)


              T                1/2      1/2             2
            CC   [phi phi] = -2    (M/r)    r sin(theta)


                 r                                   2
               CC   [phi phi] = -(r - 2 M) sin(theta)


               theta
             CC       [phi phi] = -sin(theta) cos(theta)

a similar analysis using the geodesic equations and eliminating tau should be possible.

**CarlB** · Jul24-06, 10:15 PM

Originally Posted by pervect

Carl (and Jorrie, too) I think you've totally missed the point.

The 'r' in the Scwarzschild coordinate system is not the same 'r' in the isotropic coordinate system.

In fact, one can see that r = r'(1+2m/r')^2 from the webpage, where r' is the isotropic radial coordinate and r is the Schwarzschild r coordinate.

Neither of these 'r' coordinates can be interpreted as 'radial distance'. The isotropic coordinates, though, define a set of coordinates where light cones will apear to be circles.

Thus assuming that r^2 = x^2 + y^2 does not make any physical sense. This is assumed implicitly when one writes x = r cos(phi) y = r sin(phi), one assumes that x^2 + y^2 = r^2. One can think of it as a purely formal mathematical expression, I suppose.

As such a formal exercise, unless Carl has made an algebraic error he's got a new and totally non-standard metric that doesn't even have a name in terms of new and toatally non-standard coordiantes that don't have names either (except what Carl has christened them). They definitely shouldn't be confused with Cartesian coordinates, which is what Carl has chosen to call his variables.

Aside from the fact that it is non-standard, and thus will be hard to check solutions for, even deriving the correct solutions is going to be a giant pain. Re-writing a simple metric in a more complicated form is not going to get anywyere as far as generating a solution for the metric.

It is my intention to work only in the Schwarzschild metric, but with the polar coordinates replaced by Cartesian coordinates. I realize that you think this is unphysical, but here in the United States computers are built with screens that have rectangular arrays of pixels, and so any computation I make must be eventually so converted. Now that you've explained what they are, I don't mean to use "isotropic" coordinates at all. It's just extra confusion you've brought in, not me. I thought it was you calling my coordinates "isotropic". Recall, it was I who "allowed" you to call them that if you like.

Nor am I the only person in the world who calls these things I'm using "Cartesian". The United States is a very primitive place, filled with horses, saloons and cowboys. Other obscure, amateurish authors here use the same confusing and incorrect terminology. For example see:

Numerical Relativity for Inspiraling Binaries in Co-Rotating Coordinates: Test Bed for Lapse and Shift Equations
Kip S. Thorne, 1998
http://arxiv.org/abs/gr-qc/9808024

By the way, I should mention that while I only took one graduate level General Relativity class at the University of California, Irvine, I was the top student, the only student in the class receiving an "A+". I don't think that taking one class makes you an expert in a subject, but you can be sure that I understand the principles of general relativity at least a little better than some.

[edit] But I see that it is not only on this side of the Atlantic that one finds the use of Cartesian coordinates in describing black holes:

Shortcut Method of Solution of Geodesic Equations for Schwarzschild Black Hole
Jean-Alain Marck,
D´epartement d’Astrophysique Relativiste et de Cosmologie
Observatoire de Paris - Section de Meudon
F-92195 Meudon Cedex, France

Classical and Quantum Gravity 13 (1996) 393-402

$LaTeX Code: r^2 = x^2 + y^2 + z^2\\;\\;\\;\\;\\;\\;\\;(6)$

http://arxiv.org/abs/gr-qc/9505010
[/edit]

Carl

**CarlB** · Jul24-06, 11:02 PM

Originally Posted by pervect

Unless I've made a mistake (quite possible, even with computer assistance) the equations with proper time eliminated for the Schwarzschld metric are just:

This is exactly what I need and much thanks. To convert it into Cartesian form, I need only split the radial acceleration into x and y accelerations. I should have a new simulation running later today, but Seattle is sooooooo hot (about 90 F) that everyone has left for the beach and very little work is getting done. Maybe I'll wait until it begins to cool down.

I was in mathematics grad school at the University of Washington 10 years ago and got a student priced copy of Mathematica. I played around with it, but never worked at it long enough to get around to messing with the symbolic calculus. Since then, most of my work has been in Clifford algebra and for that I wrote Java programs. But after this fiasco I'm convinced that I need to get some sort of a symbolic math package.

Carl

**pervect** · Jul24-06, 11:42 PM

Symbolic math really rocks. Couldn't live without it.

I've also added horizon crossing solutions to the Painleve thread (but only in terms of conseved quantites) I haven't eliminated dt/dtau, it appears to be a royal pain to do so. dt/dtau will range in values from 1 to E for the horizon crossing solutions.

I note that non-horizon crossing alternate solutions exist, but I haven't explored them very far.

I've also computed d^2 r / dtau^2 only in the special case where dr/dtau = 0.

**CarlB** · Jul25-06, 05:52 AM

Originally Posted by pervect

$LaTeX Code: \\frac{d^2 r}{dt^2} = \\frac{2\\,m}{r^2 \\left( 1 - \\frac{2m}{r} \\right) } \\left( \\frac{dr}{dt} \\right) ^2 - \\left( 1-{\\frac {2 m}{r}} \\right) \\frac{m}{r^2}+ \\left( 1-{\\frac {2m}{r}} \\right) ^{3}\\frac{{K}^{2}}{r^3} $

Unfortunately, I couldn't get it to work. This could always be a programming error of course. It acts like there is too much energy in the angular kinetic energy term, that is the K^2 term. In fact, the simulation doesn't match Newton at large r because of this term. To get it to work, I have to divide by an extra factor of r.

Take a look at the units for large r, in the first and third terms of the above and see if there is a transcription error, most likely in the last term of the above, but at least conceivably in the formula for $LaTeX Code: K = \\dot{\\phi} r^2 / (1-2m/r)$ which can be used to simplify the third term.

Having a differential equation in r and a constant of the motion in phi is unsatisfying to me. I realize that I can put the Newtonian equation into the same form, but I am one of those crazy guys who thinks that the equations of motion are closer to the physics than the constants of the motion. Yes, I believe that the constants of the motion are nothing more than convenient methods we use for simplifying the original physical problem, which to me are always defined in terms of positions and velocities. So it looks like I will continue my process of deriving the Newton style equations of motion for the Schwarzschild metric.

Carl

**CarlB** · Jul25-06, 06:00 AM

Step (2), write the proper length of a path as an integral over coordinate time.

Here's where I begin:

Originally Posted by CarlB

$LaTeX Code: ds^2 = \\frac{(r-2)}{r}dt^2 - \\frac{(r^3-2y^2)}{(r^3-2r^2)}dx^2 - \\frac{(r^3-2x^2)}{(r^3-2r^2)}dy^2 -\\frac{4rxy}{(r^3-2r^2)}dx\\;dy$

Note that in the above I've swapped over to the other coast metric in order to make s time like instead of spacelike. Sorry for starting the problem with the spacelike signature, it's a habit from Euclidean relativity where s is interpreted as a little curled up spatial dimension. Dividing by (dt)^2 gives:

$LaTeX Code: \\left(\\frac{ds}{dt}\\right)^2 = \\frac{(r-2)}{r} - \\frac{(r^3-2y^2)}{(r^3-2r^2)}\\dot{x}^2 - \\frac{(r^3-2x^2)}{(r^3-2r^2)}\\dot{y}^2 -\\frac{4rxy}{(r^3-2r^2)}\\dot{x}\\dot{y}$

where the dots indicate derivatives with respect to coordinate time t rather than the more usual proper time. Thus the proper time along the path is:

$LaTeX Code: S = \\int_{t_1}^{t_2} \\frac{ds}{dt}\\; dt$

$LaTeX Code: = \\int_{t_1}^{t_2} \\sqrt{\\frac{(r-2)}{r} - \\frac{(r^3-2y^2)}{(r^3-2r^2)}\\dot{x}^2 - \\frac{(r^3-2x^2)}{(r^3-2r^2)}\\dot{y}^2 -\\frac{4rxy}{(r^3-2r^2)}\\dot{x}\\dot{y} }\\;\\; dt$

**CarlB** · Jul25-06, 06:12 AM

The next step is to vary the path by the Euler-Lagrange equations. Discretion being the better part of valor, I will wait until I have a symbolic mathematics engine before I attempt that slope. Basically, I need to calculate the partial derivatives with respect to the coordinates and velocities. The equations of motion are then the time derivatives of the partials with respect to velocities set equal to the partials with respect to positions. That is, for example,

$LaTeX Code: \\frac{d}{dt} \\left(\\frac{\\partial ds/dt}{\\partial x}\\right) = \\frac{\\partial ds/dt}{\\partial \\dot{x}}$

Carl

Okay, now I've got Maxima running on my laptop. Pretty cool. And cheap.