Schwarzschild Orbits in Cartesian coordinates

[CarlB]

Looks like my previous effort dropped a 2 on the r squared term in the denominator, so as usual, let me begin by repairing the previous damage.

$$I = 1-\frac{2}{r}-(\dot{x}^2+\dot{y}^2) -2\frac{(x\dot{x}+y\dot{y})^2}{r^3-2r^2}$$

The first order partial derivatives are:

$$\left(\frac{\partial I}{\partial x}\right) = \frac{2x}{r^3}-\frac{4\dot{x}(x\dot{x}+y\dot{y})}{r^3-2r^2} +2\frac{(x\dot{x}+y\dot{y})^2 (3rx-4x)}{(r^3-2r^2)^2}$$

$$\left(\frac{\partial I}{\partial y}\right) = \frac{2y}{r^3}-\frac{4\dot{y}(x\dot{x}+y\dot{y})}{r^3-2r^2} +2\frac{(x\dot{x}+y\dot{y})^2 (3ry-4y)}{(r^3-2r^2)^2}$$

$$\left(\frac{\partial I}{\partial \dot{x}}\right) = -2\dot{x}-\frac{4x(x\dot{x}+y\dot{y})}{r^3-2r^2}$$

$$\left(\frac{\partial I}{\partial \dot{y}}\right) = -2\dot{y}-\frac{4y(x\dot{x}+y\dot{y})}{r^3-2r^2}$$


2nd order partials with respect to velocity:

$$\left(\frac{\partial^2 I}{\partial \dot{x}^2}\right) = -2-\frac{4x^2}{r^3-2r^2}$$

$$\left(\frac{\partial^2 I}{\partial \dot{y}^2}\right) = -2-\frac{4y^2}{r^3-2r^2}$$

$$\left(\frac{\partial^2 I}{\partial \dot{x}\partial \dot{y}}\right) = -\frac{4xy}{r^3-2r^2}$$

Mixed 2nd order partials:

$$\left(\frac{\partial^2 I}{\partial \dot{x}\partial x}\right) = -\frac{8x\dot{x}+4y\dot{y}}{r^3-2r^2} +\frac{(x\dot{x}+y\dot{y}) (12r-16)x^2}{(r^3-2r^2)^2}$$

$$\left(\frac{\partial^2 I}{\partial \dot{y}\partial y}\right) = -\frac{8y\dot{y}+4x\dot{x}}{r^3-2r^2} +\frac{(x\dot{x}+y\dot{y}) (12r-16)y^2}{(r^3-2r^2)^2}$$

$$\left(\frac{\partial^2 I}{\partial \dot{x}\partial y}\right) = -\frac{4x\dot{y}}{r^3-2r^2} +\frac{(x\dot{x}+y\dot{y}) (12r-16)xy}{(r^3-2r^2)^2}$$

$$\left(\frac{\partial^2 I}{\partial \dot{y}\partial x}\right) = -\frac{4y\dot{x}}{r^3-2r^2} +\frac{(x\dot{x}+y\dot{y}) (12r-16)xy}{(r^3-2r^2)^2}$$

It remains to substitute these into post #26, and then simplify. In doing this, it is clear that I'll want to cancel out the factors of $$r^3-2r^2,$$ which was the cause of the problem at r=2 that I wanted to avoid by using coordinate time instead of proper time.

Of the terms A, B, C, D, E, and F of post #26, the number of times $$r^3-2r^2$$ appears in the denominator is 2 for A, B, D and E, and 3 for C and F. Thus there will be Consequently, in the formulas for $$\ddot{x}, \ddot{y},$$ there will be 4 cancellations of $$1/(r^3-2r^2)$$ between the numerator and denominator, and the numerator will have a left over copy of $$1/(r^3-2r^2).$$ This means that, at first glance, the solution blows up at r=2. Alternatively, we should be able to factor 2 copies of (r-2) out of the numerator.

It's getting late and I'm going to start making mistakes (and possibly typing with a Texas accent). Let me rewrite out the results of this page in form XXX/(r^n(r-2)^m, probably drop a 2, and then go to bed. [Nope, a man's got to know his limitations.]

Carl

 

[CarlB]

Factoring out r and r-2 gives:
$$\begin{array}{rcl} I &=& (r(r-2)^2-(\dot{x}^2+\dot{y}^2)r^2(r-2)-2(x\dot{x}+y\dot{y})^2) /r^2(r-2)\\ \left(\frac{\partial I}{\partial x}\right) &=& (2xr(r-2)^2 - 4\dot{x}(x\dot{x}+y\dot{y})r^2(r-2) +2x(x\dot{x}+y\dot{y})^2(3r-4) ) /r^4(r-2)^2\\ \left(\frac{\partial I}{\partial y}\right) &=& (2yr(r-2)^2 - 4\dot{y}(x\dot{x}+y\dot{y})r^2(r-2) +2y(x\dot{x}+y\dot{y})^2(3r-4) ) /r^4(r-2)^2\\ \left(\frac{\partial I}{\partial \dot{x}}\right) &=& (-2\dot{x}r^2(r-2) -4x(x\dot{x}+y\dot{y}))/r^2(r-2)\\ \left(\frac{\partial I}{\partial \dot{y}}\right) &=& (-2\dot{y}r^2(r-2) -4y(x\dot{x}+y\dot{y}))/r^2(r-2)\\ \left(\frac{\partial^2 I}{\partial \dot{x}^2}\right) &=& (-2r^2(r-2)-4x^2)/r^2(r-2)\\ \left(\frac{\partial^2 I}{\partial \dot{y}^2}\right) &=& (-2r^2(r-2)-4y^2)/r^2(r-2)\\ \left(\frac{\partial^2 I}{\partial \dot{x}\partial \dot{y}}\right) &=& (-4xy)/r^2(r-2)\\ \left(\frac{\partial^2 I}{\partial \dot{x}\partial x}\right) &=& (-(8x\dot{x}+4y\dot{y})r^2(r-2) + (x\dot{x}+y\dot{y}) (12r-16)x^2)/r^4(r-2)^2\\ \left(\frac{\partial^2 I}{\partial \dot{y}\partial y}\right) &=& (-(8y\dot{y}+4x\dot{x})r^2(r-2) + (x\dot{x}+y\dot{y}) (12r-16)y^2)/r^4(r-2)^2\\ \left(\frac{\partial^2 I}{\partial \dot{x}\partial y}\right) &=& (-4x\dot{y}r^2(r-2) + ( x\dot{x}+y\dot{y}) (12r-16)xy)/r^4(r-2)^2\\ \left(\frac{\partial^2 I}{\partial \dot{y}\partial x}\right) &=& (-4y\dot{x}r^2(r-2) + ( x\dot{x}+y\dot{y}) (12r-16)xy)/r^4(r-2)^2 \end{array}$$

Let me first just see what happens to all those 1/r and 1/(r-2) factors. Treating these alone, I have

$$\begin{array}{rcl} A &=& (XXX)/r^4(r-2)^2\\ B &=& (XXX)/r^4(r-2)^2\\ C &=& (XXX)/r^6(r-2)^3\\ D &=& (XXX)/r^4(r-2)^2\\ E &=& (XXX)/r^4(r-2)^2\\ F &=& (XXX)/r^6(r-2)^3 \end{array}$$

where "XXX" are various polynomials in the phase space variables. The equations of motion are then of the form:

$$\begin{array}{rcl} \ddot{x} &=& (BF-CE)/(BD-AE)\\ &=&\frac{(XXX)/r^{10}(r-2)^5}{(YYY)/r^8(r-2)^4}\\ &=&\frac{XXX}{YYY}\frac{1}{r^2(r-2)}.\end{array}$$

So I don't need to keep around all the r and r-2 factors dividing I and the various partial derivatives. Instead, I just keep around the numerators. Let me rewrite these various terms as primed:

$$\begin{array}{rcl} I' &=& r(r-2)^2-(\dot{x}^2+\dot{y}^2)r^2(r-2)-2(x\dot{x}+y\dot{y})^2\\ \left(\frac{\partial I}{\partial x}\right)' &=& 2xr(r-2)^2 - 4\dot{x}(x\dot{x}+y\dot{y})r^2(r-2) +2x(x\dot{x}+y\dot{y})^2(3r-4)\\ \left(\frac{\partial I}{\partial y}\right)' &=& 2yr(r-2)^2 - 4\dot{y}(x\dot{x}+y\dot{y})r^2(r-2) +2y(x\dot{x}+y\dot{y})^2(3r-4)\\ \left(\frac{\partial I}{\partial \dot{x}}\right)' &=& -2\dot{x}r^2(r-2) -4x(x\dot{x}+y\dot{y})\\ \left(\frac{\partial I}{\partial \dot{y}}\right)' &=& -2\dot{y}r^2(r-2) -4y(x\dot{x}+y\dot{y})\\ \left(\frac{\partial^2 I}{\partial \dot{x}^2}\right)' &=& -2r^2(r-2)-4x^2\\ \left(\frac{\partial^2 I}{\partial \dot{y}^2}\right)' &=& -2r^2(r-2)-4y^2\\ \left(\frac{\partial^2 I}{\partial \dot{x}\partial \dot{y}}\right)' &=& -4xy\\ \left(\frac{\partial^2 I}{\partial \dot{x}\partial x}\right)' &=& -(8x\dot{x}+4y\dot{y})r^2(r-2) + (x\dot{x}+y\dot{y}) (12r-16)x^2\\ \left(\frac{\partial^2 I}{\partial \dot{y}\partial y}\right)' &=& -(8y\dot{y}+4x\dot{x})r^2(r-2) + (x\dot{x}+y\dot{y}) (12r-16)y^2\\ \left(\frac{\partial^2 I}{\partial \dot{x}\partial y}\right)' &=& -4x\dot{y}r^2(r-2) + ( x\dot{x}+y\dot{y}) (12r-16)xy\\ \left(\frac{\partial^2 I}{\partial \dot{y}\partial x}\right)' &=& -4y\dot{x}r^2(r-2) + ( x\dot{x}+y\dot{y}) (12r-16)xy \end{array}$$

The above puts the answer into rational form, that is, the answer is the ratio of two polynomials in $$x, y, \dot{x}, \dot{y}, r.$$ The only problem is that the answer is going to have a LOT of terms. Fortunately, we have MAXIMA, a free symbolic computing engine.

 

[CarlB]

In putting the above into MAXIMA, I discovered that defining r = sqrt(x^2+y^2) is a bad idea, because it prevents keeping the answer in r, which is what you want. Furthermore, I found that reasonably small results could be obtained by using factorout: factorout(expression, x, y, r); This keeps the terms that are going to be reduced by x^2+y^2 = r^2 together. The input program for the partial derivatives is:

ii:expand(r*(r-2)^2-(xx^2+yy^2)*r^2*(r-2)-2*(x*xx+y*yy)^2);
didx:expand(2*x*r*(r-2)^2-4*xx*(x*xx+y*yy)*r^2*(r-2)+2*x*(x*xx+y*yy)^2*(3*r-4));
didy:expand(2*y*r*(r-2)^2-4*yy*(x*xx+y*yy)*r^2*(r-2)+2*y*(x*xx+y*yy)^2*(3*r-4));
didxx:expand(-2*xx*r^2*(r-2)-4*x*(x*xx+y*yy));
didyy:expand(-2*yy*r^2*(r-2)-4*y*(x*xx+y*yy));
d2idxx2:expand(-2*r^2*(r-2)-4*x^2);
d2idyy2:expand(-2*r^2*(r-2)-4*y^2);
d2idxxdyy:expand(-4*x*y);
d2idxxdx:expand(-(8*x*xx+4*y*yy)*r^2*(r-2)+(x*xx+y*yy)*(12*r-16)*x^2);
d2idyydy:expand(-(8*y*yy+4*x*xx)*r^2*(r-2)+(x*xx+y*yy)*(12*r-16)*y^2);
d2idxxdy:expand(-4*x*yy*r^2*(r-2)+(x*xx+y*yy)*(12*r-16)*x*y);
d2idyydx:expand(-4*y*xx*r^2*(r-2)+(x*xx+y*yy)*(12*r-16)*x*y);

Then the A, B, C, D, E, F variables are:

aa:expand(ii*d2idxx2-0.5*didxx^2);
bb:expand(ii*d2idxxdyy-0.5*didxx*didyy);
cc:expand((0.5*didxx*didx-ii*d2idxxdx)*xx+(0.5*didxx*didy-ii*d2idxxdy)*yy+ii*didx);
dd:expand(ii*d2idxxdyy-0.5*didyy*didxx);
ee:expand(ii*d2idyy2-0.5*didyy^2);
ff:expand((0.5*didyy*didy-ii*d2idyydy)*yy+(0.5*didyy*didx-ii*d2idyydx)*xx+ii*didy);

Finally, the x and y numerators and the shared denominator are defined as:

xxxn:expand(bb*ff-cc*ee);
yyyn:expand(cc*dd-aa*ff);
denom:expand(bb*dd-aa*ee);

To get the MAXIMA program to reduce these to sizes that an old forklift driver feels like crunching by hand on, I used:

grind(factorout(denom,x,y,r));
grind(factorout(xxxn,x,y,r));

and did the denominator first because it's smaller. Running the MAXIMA program, we find that the denominator is particularly easy and reduces to:

$$B'*D'-A'*E' = 8r^6(r-2)^3(x\dot{x}+y\dot{y})^2 +4r^7(r-2)^5+4r^8(r-2)^4(\dot{x}^2+\dot{y})^2$$

The Schwarzschild coordinates halt at the event horizon (r=2), and we already have an (r-2) factor in the denominator besides the three above, so we expect to factor at least one more (r-2) from the numerator, for a total of at least 5.

Ouch! The numerator is hard to reduce. Time to get some sleep.

 

[CarlB]

The Schwarzschild coordinates halt at the event horizon (r=2), and we already have an (r-2) factor in the denominator besides the three above, so we expect to factor at least one more (r-2) from the numerator, for a total of at least 5. Ouch! The numerator is hard to reduce. Time to get some sleep.

Well no matter what i do I can't factor out 3 copies of (r-2) from the numerator. I went back and checked the equations going in. No problems. I abused MAXIMA by rewriting all the equations with r replaced with s = r+2, and then derived the equation for small s. Finally I solved the differential equation for small s and discovered that one doesn't need to suppose that there are 5 copies of (r-2) in order to make particles stall on the event horizon. 3 copies is enough, and 3 is what I found. This means that the final equation for the acceleration is going to have one copy of (r-2) in the denominator.

Most of the terms in the numerator will have 2 copies of (r-2) and so will be finite across the event horizon. The remaining terms will tend to infinity but will be kept finite (and in fact will tend to zero) because they are proportional to 3 or 4 copies of the velocity.

Technically, none of this is much of an issue for the Schwarzschild metric as orbits never reach the event horizon. But it is an issue for the Painleve metric. So my guess is that the Painleve metric will have two more factors of (r-2) in the numerator and so will be a bit more elegant. And it could be a problem that would cause a division by zero in a simulation of the Schwarzschild metric.

To make the division by zero explicit, I will factor the terms according to powers of (r-2). That might have some use in the computer program, if it turns out that orbits that crash onto the event horizon cause divisions by zero.

By the way, the division by zero is not a consequence of doing the problem in x-y coordinates. I know this because I obtained the differential equation in (q,t) by substituting y=0, q = x-2 > 0, so that r-2 = q. These are equivalent to spherical coordinates for a particle on a ray through the origin. I guess these things are obvious to people who mess around with gravity all the time, but I've certainly learned a few things. But if it was obvious to them I would have appreciated their informing me not to expect 5 copies of (r-2) in the numerator. On the other hand, since some of them had voted that I wasn't going to get this completed, not mentioning anything sort of makes sense.

And now, let me finish off the Schwarzschild numerator, cancel what terms I can, and then test the equations in my gravity simulator.

Carl

 

[CarlB]

As usual, my earlier formula for the numerator is off a bit. The corrected formula is:

$$B'D'-A'E' = 4 r^6 (r-2)^3 (2(x\dot{x}+y\dot{y})^2 + r^2(r-2)(\dot{x}^2+\dot{y}^2)-r(r-2)^2)$$

For the numerator, MAXIMA gives me:

+4*r^5*s^7*x

+4*r^6*s^6*(x*yy^2-2*xx*y*yy-x*xx^2)

-4*r^4*s^5*(2*r^3*x*yy^4-2*r^3*xx*y*yy^3+3*r*x*y^2*yy^2-2*x*y^2*yy^2+2*r^3*x*xx^2*yy^2+4*xx*y^3*yy-2*r^3*xx^3*y*yy+6*r*x^2*xx*y*yy+4*x*xx^2*y^2+3*r*x^3*xx^2+2*x^3*xx^2)

+4*r^5*s^4*(y*yy+x*xx)*(3*r*x*y*yy^3-8*x*y*yy^3+8*xx*y^2*yy^2+3*r*x^2*xx*yy^2-4*x^2*xx*yy^2+3*r*x*xx^2*y*yy+4*xx^3*y^2+3*r*x^2*xx^3)

+8*r^3*s^3*(y*yy+x*xx)^3*(3*r*x*y*yy-4*x*y*yy+4*xx*y^2+3*r*x^2*xx)

In the above $$s = r-2, xx = \dot{x}, yy =\dot{y}.$$ To get the final result, I need to divide by $$r^2(r-2)$$ times the denominator calculated two posts ago:

$$r^2(r-2)(B'*D'-A'*E') = 4 r^8 (r-2)^4 (2(x\dot{x}+y\dot{y})^2 + r^2(r-2)(\dot{x}^2+\dot{y}^2)-r(r-2)^2)$$

Some of the factors of r and r-2 will cancel between the numerators and denominators. The resulting equation for the accelerations of the Schwarzschild metric in "flat" coordinates is:

$$\begin{array}{rcl} \ddot{x} &=& (+r(r-2)^3x +r^2(r-2)^2(x\dot{y}^2-2\dot{x}y\dot{y}-x\dot{x}^2)\\ &&\;\;-(r-2)(2r^3x\dot{y}^4-2r^3\dot{x}y\dot{y}^3+3rxy^2\dot{y}^2-2xy^2\dot{y}^2+2r^3x\dot{x}^2\dot{y}^2+4\dot{x}y^3\dot{y}-2r^3\dot{x}^3y\dot{y}+6rx^2\dot{x}y\dot{y}+4x\dot{x}^2y^2+3rx^3\dot{x}^2+2x^3\dot{x}^2)\\&&+r(y\dot{y}+x\dot{x})(3rxy\dot{y}^3-8xy\dot{y}^3+8\dot{x}y^2\dot{y}^2+3rx^2\dot{x}\dot{y}^2-4x^2\dot{x}\dot{y}^2+3rx\dot{x}^2y\dot{y}+4\dot{x}^3y^2+3rx^2\dot{x}^3))\\ &&/(4 r^4 (2(x\dot{x}+y\dot{y})^2 + r^2(r-2)(\dot{x}^2+\dot{y}^2)-r(r-2)^2))\\ &+&(2(y\dot{y}+x\dot{x})^3(3rxy\dot{y}-4xy\dot{y}+4\dot{x}y^2+3rx^2\dot{x}))\\ &&/(4 r^5 (r-2) (2(x\dot{x}+y\dot{y})^2 + r^2(r-2)(\dot{x}^2+\dot{y}^2)-r(r-2)^2)) \end{array}$$

$$\begin{array}{rcl} \ddot{y} &=& (+r(r-2)^3y +r^2(r-2)^2(y\dot{x}^2-2\dot{y}x\dot{x}-y\dot{y}^2)\\ &&\;\;-(r-2)(2r^3y\dot{x}^4-2r^3\dot{y}x\dot{x}^3+3ryx^2\dot{x}^2-2yx^2\dot{x}^2+2r^3y\dot{y}^2\dot{x}^2+4\dot{y}x^3\dot{x}-2r^3\dot{y}^3x\dot{x}+6ry^2\dot{y}x\dot{x}+4y\dot{y}^2x^2+3ry^3\dot{y}^2+2y^3\dot{y}^2)\\&&+r(x\dot{x}+y\dot{y})(3ryx\dot{x}^3-8yx\dot{x}^3+8\dot{y}x^2\dot{x}^2+3ry^2\dot{y}\dot{x}^2-4y^2\dot{y}\dot{x}^2+3ry\dot{y}^2x\dot{x}+4\dot{y}^3x^2+3ry^2\dot{y}^3))\\ &&/(4 r^4 (2(y\dot{y}+x\dot{x})^2 + r^2(r-2)(\dot{y}^2+\dot{x}^2)-r(r-2)^2))\\ &+&(2(x\dot{x}+y\dot{y})^3(3ryx\dot{x}-4yx\dot{x}+4\dot{y}x^2+3ry^2\dot{y}))\\ &&/(4 r^5 (r-2) (2(y\dot{y}+x\dot{x})^2 + r^2(r-2)(\dot{y}^2+\dot{x}^2)-r(r-2)^2)) \end{array}$$

Let me check this to see if it gives the right answer in the Newtonian limit (i.e. $$\dot{x}=\dot{y}=0$$ ). It checks out, that is, I get $$x r^4/(-r^7) = -x/r^3.$$ The next step is to program this into my all gravity simulator and see if it looks okay. I think that's enough for one night. I can hardly wait to see the Painleve coordinates.

Carl

 

[CarlB]

Some of the factors of r and r-2 will cancel between the numerators and denominators. The resulting equation for the accelerations of the Schwarzschild metric in "flat" coordinates is:

It turns out that the three things that are most useful in MAXIMA are "grind", "factorout", and "factor". I applied "factor" to the above equations of motion and eventually got them into a fairly simple form. I managed to cancel out everything in the denominator except $$r^5(r-2),$$ and put the numerators into a form where I can do the x and y coordinates by a combined calculation. The simplified equations of motion are now:

$$\begin{array}{rcl} \ddot{x} &=& (x(3r(x\dot{x}+y\dot{y})^2-r(r-2)^2)\\ &&+4y(y\dot{x}-x\dot{y})(x\dot{x}+y\dot{y})\\ &&+2\dot{y}(y\dot{x}-x\dot{y})r^2(r-2)) /(r^5(r-2)) \\ \ddot{y} &=& \left(y(3r(x\dot{x}+y\dot{y})^2-r(r-2)^2)\\ &&-4x(y\dot{x}-x\dot{y})(x\dot{x}+y\dot{y})\\ &&-2\dot{x}(y\dot{x}-x\dot{y})r^2(r-2)\right) /(r^5(r-2)) \end{array}$$

For comparison, here are the equations in polar coordinates, from pervect's post #17, with m=1:

$$\begin{array}{rcl} \ddot{r} &=& \frac {3 \dot{r}^2}{r(r-2)} + (r-2)(\dot{\phi}^{2}-{\frac {1}{r^3}}) \\ \ddot{\phi} &=& -\frac {2 \dot{r}\dot{\phi}(r-3)}{r(r-2)} \end{array}$$

Are these the same? One quick test is to see what accelerations they give for a test mass that is stationary. So set all velocities to zero and get:

$$\begin{array}{rcl} \ddot{x} &=& -xr(r-2)^2 /(r^5(r-2)) = -x(r-2)/r^4 \\ \ddot{y} &=& -y(r-2)/r^4\\ \ddot{r} &=& -(r-2)/r^3 \\ \ddot{\phi} &=& 0 \end{array}$$

Putting x=r, we see that they are identical. In any case, I've programmed the (x,y) version into my computer and sure enough it does draw what seem like reasonable orbits. I've not yet updated the Java applet because I want to include the Painleve coordinates too. These will be a great check of the equations because they should send particles down the same orbits but with different times (and should spiral into the singularity at the center of the black hole rather than getting stuck on the unphysical event horizon). I've also slightly changed the simulation so that it draws the gray "sun" out to the event horizon, r=2, rather than out to r=1. This makes the Schwarzschild tracks get stuck on the edge of the gray region, which is what you expect sort of.

I'm note sure if these are different. Let me run a simulation...

I still can't get pervect's equation to run correctly. This could be because I am converting from polar to rectangular coordinates wrong, but the first parts of the orbits match, so I think it's probably okay.

What happens is that when the orbit goes down near the event horizon (i.e. so d phi /d t is maximum), the test mass ends up picking up huge amounts of energy and then zooms off to infinity. Even when the orbit doesn't get near the event horizon, they tend to pick up energy. It's like the calculation for acceleration in phi is too large.

I've just tested my version against the impact parameter, $$3\sqrt{3}/2$$ which is near 5.2 and it passes with flying colors. That is, when you throw relativistic particles at it from (+1000,+5.2) with velocity (-1,0), they do interesting things. Also, for my programming of both my version and pervect's, the scattering is approximately twice the Newtonian scattering.

Another test is to compare with:
http://www.fourmilab.ch/gravitation/orbits/
but it's late and I haven't figured out how to convert from his proper time coordinates to my system time. His simulation has a particle whose maximum radius is 12 has interesting orbits if its angular momentum L=3.57, I get similar orbits for x=12 when I set dy/dt = 0.26084. Converting this into angular momentum, one first gets dphi/dt = 0.26084/12. This gives:

$$3.57= L = r^2 \frac{d\phi}{d\tau} = 12^2 \frac{0.26084}{12}\frac{dt}{d\tau}$$

which gives me that $$dt/d\tau = 3.57/(12 X .26084) = 1.14055$$, which seems reasonable, but I don't know if it's close to exact. To check this, I have to substitute into the Schwarzschild metric:

$$\begin{array}{rcl} \left(\frac{d\tau}{dt}\right)^2 &=& \frac{r-2}{r} - \frac{r}{r-2}\left(\frac{dr}{dt}\right)^2 - r^2\left(\frac{d\phi}{dt}\right)^2\\ &=& \frac{10}{12} - 0 - \left(\frac{dy}{dt}\right)^2\\ &=& (1/1.1431)^2\\ \end{array}$$

and the first couple digits of dt/ds match. Therefore, I conclude that I've probably got the correct Cartesian differential equation for the Schwarzschild metric, as desired. Next, on to the Painleve metric, which may or may not be more difficult.

Carl

 

[Fjern]

What happens is that when the orbit goes down near the event horizon (i.e. so d phi /d t is maximum), the test mass ends up picking up huge amounts of energy and then zooms off to infinity. Even when the orbit doesn't get near the event horizon, they tend to pick up energy. It's like the calculation for acceleration in phi is too large.
Carl

I have also encountered this phenoma. Total energy of the particle suddenly starts for infinity, while total angular momentum falls a bit less dramatically.
I am working with spherical coordinates in 3D, with no simlifications of the specific metric. I am solving the geodesic equation for all for coordinates (t,r,theta,phi).

and by the way, if you write metric in cartesian coordinates you get something like
g_tt = -1 + 2m/r
g_xx = g_yy=g_zz = 1 +2m/r
for covariant version of metric tensor.
If you put them in geodesic equation, work out the christoffel symbols and there you go...nice, clean and easy.

 

[pervect]

I also had numerical stability problems back when I did this. The reason I believed my results is that the problems were lessened when I tightened the tolerance on the integration algorithm. So I put the problem down to numerical integration issues, rather than an incorrect equation.

 

[Chris Hillman]

Some possibly useful citations

This old thread seems to continue another thread, so at this late date it is rather hard to follow. For one thing, in this thread, Carl Brannen didn't define what he meant by "cartesian-like coordinates", although presumably he did so in the original posts. I conjecture he might have been trying to reinvent one of two classic types of chart: Kerr-Schild charts, or harmonic charts. If you search the arXiv using the search bar at the site in my sig on ""Kerr-Schild numerical gr-qc" and "harmonic coordinates" numerical gr-qc" quite a few results turn up, which might help in tracking down the source of the problem. BTW, see Chandrasekhar, The Mathematical Theory of Black Holes for the exact solution of the geodesic equations in the Kerr vacuum, including diagrams, which you can use as a "benchmark" for testing numerical integration schemes.

Incidently, gr-qc/0704.2508 is a new eprint extending a classical result of Mueller zum Hagen to the effect that a static vacuum solution (i.e. with timelike Killing vector $X$) can be written in a harmonic chart such that the metric components are analytic in this chart. In the eprint, Tod (coauthor of the LMS student series textbook on gtr) extends this to electrovacuum solutions, irrespective of whether or not the EMF shares the timelike Killing vector (i.e. $L_X F$ might be nonzero). This shows one way in which harmonic charts can be useful for theoretical work, but as the discussion in Chandrasekhar shows, they may not really be the best for searching for closed form solutions of the geodesic equations. (I guess everyone here knows that the fact that such can be found for the Kerr vacuum is in fact one of the more remarkable properties of this solution, which depends upon the existence of a Killing tensor; see D'Inverno's textbook.)

 

[pervect]

Carl's version of "cartesian" coordinates are defined back in post #2

x = r cos phi
y = r sin phi

and I believe it is implied that theta =0, as one is interesting in orbital simulations.

Here r and phi are the Schwarzschild r and phi coordinates.

I never understood the motivation for doing this, you'd have to ask Carl if he's still around. I had assumed that Carl simply confused r with radial distance, so I complained a bit (which didn't go over very well with him if you read the thread). But possibly he has some other motivation, who knows.

The thing that is interesting is that Maple gave pages and pages of output for the differential equations for x and y using these coordinates. If Carl is correct, Maple didn't simplify the equations completely (which is easy to believe).

I don't think this has much to do with the references cited by Chris from my quick scan.

 

[CarlB]

Chris, the whole thing is here. The "cartesian coordinates" are defined in the first few posts, and should be obvious to anyone in the field. All they amount to are the usual (r,theta, phi,t) coordinates converted to (x,y,z,t), which is commonly done in the industry.

I managed to get that part of the program to work and the simulation runs beautifully. It is here:
http://www.gaugegravity.com/testapplet/SweetGravity.html

What is incomplete is the work on Painleve coordinates. After I get those to run, I'll add a whole bunch of initial conditions that the student can select. These will illustrate all the various cool facts about classical and GR gravity.

The reason for including Painleve coordinates (which frankly, I thought would be simpler than the usual Cartesian GR), is so that test particles can hit the singularity at zero, while in the usual Cartesian they get stuck on the event horizon.

Carl

 

[CarlB]

I never understood the motivation for doing this, you'd have to ask Carl if he's still around. I had assumed that Carl simply confused r with radial distance, so I complained a bit (which didn't go over very well with him if you read the thread). But possibly he has some other motivation, who knows.

For now my secrets will remain my own.

Also, I would like to have equations that allow theta non zero. I don't allow that in the simulation now, but it would be nice to have the ability to model Earth orbit crossing comets, etc.

The thing that is interesting is that Maple gave pages and pages of output for the differential equations for x and y using these coordinates. If Carl is correct, Maple didn't simplify the equations completely (which is easy to believe).

I was unfamiliar with MAXIMA before starting this. From what I know now, I've got doubts about it. For example, I can't seem to solve the following six simple quadratic equations in six unknowns:

$$\begin{equation} \begin{array}{rcl} F_I&=&F_IF_I+F_JF_K+F_KF_J+F_RF_R+F_GF_G+F_BF_B,\\ F_J&=&F_IF_J+F_JF_I+F_KF_K+F_RF_G+F_GF_B+F_BF_R,\\ F_K&=&F_IF_K+F_JF_J+F_KF_I+F_RF_B+F_GF_R+F_BF_G,\\ F_R&=&F_IF_R+F_JF_G+F_KF_B+F_RF_I+F_GF_K+F_BF_J,\\ F_G&=&F_IF_G+F_JF_B+F_KF_R+F_RF_J+F_GF_I+F_BF_K,\\ F_B&=&F_IF_B+F_JF_R+F_KF_G+F_RF_K+F_GF_J+F_BF_I, \end{array} \end{equation}$$

I had to solve it by hand. See setion (5) of
http://brannenworks.com/dmfound.pdf

I'm curious to see what a higher dollar program can do with them. Part of the problem is the messy topology. There are about a dozen finite solutions (i.e. 0-manifolds or point solutions), but there are also a couple of 2-manifolds that do not get near the point solutions. If I ever teach college algebra again, I'll assign problems like this as "take home extra credit ".

An example point solution is $$(F_I,F_J,F_K,F_R,F_G,F_B) = (1,0,0,0,0,0)$$. For the 2-manifold solutions, see the above link.

Carl

 

[Chris Hillman]

Points to ponder

Hi, Carl,

Chris, the whole thing is here. The "cartesian coordinates" are defined in the first few posts, and should be obvious to anyone in the field.

I guess you might be saying that you started with the standard "polar spherical" Schwarzschild chart and employed
$$x = \sin (\theta) \, \cos(\phi), \; y = \sin (\theta) \, \sin(\phi), \; z = \cos(\theta)$$
to obtain a new chart.

(Additional: I see now that in a post from some months ago, Carl did say that at least at that point in time this is exactly what he was doing.)

Fine, but trust me, that is not neccessarily the first thing an expert might guess "pseudo-cartesian coordinates" should mean in this context. (My first guess would have been something similar, but starting with the standard "spatially isotropic" chart.) Such a "pseudo-cartesian" version of the usual Schwarzschild chart would not neccessarily be a particularly good idea if you want "simple" geodesic equations; depending upon how you define "simple" there may be better choices.

I still have no idea what you were trying to do, BTW. Since this is a very old thread which someone revived, could you indulge me and remind us what you meant by "Cartesian DE form"? If you can express this precisely I guess you will turn out to be placing some demands upon the form of the geodesic equations in the chart you seek.

All they amount to are the usual (r,theta, phi,t) coordinates converted to (x,y,z,t), which is commonly done in the industry.

Unless I misunderstand what "industry" you have in mind (generic computer graphics? generic engineering?), these applications probably don't involve curved manifolds. For example, if you aren't familiar with the ways of curved manifolds, you are probably at least partially misinterpreting your radial coordinate. Roughly, "spherical-type charts" can preserve various properties of the radial coordinate familar from the usual polar spherical chart on $E^3$, but never all of them.

I managed to get that part of the program to work and the simulation runs beautifully.

Did you compare your plots with analytical solutions transformed into your coordinate chart?

What is incomplete is the work on Painleve coordinates. After I get those to run, I'll add a whole bunch of initial conditions that the student can select. These will illustrate all the various cool facts about classical and GR gravity.

Sounds great--- just make sure your plots are correct!

The reason for including Painleve coordinates (which frankly, I thought would be simpler than the usual Cartesian GR), is so that test particles can hit the singularity at zero, while in the usual Cartesian they get stuck on the event horizon.

Right, that's one reason why I suggested them. Another is that we have the Killing vector which suggests a quotient space, and we have the slicing into constant Painleve time. These each suggest different ways of reducing to a kind of "three-dimensional space". Of course, projections of geodesics will not be geodesics in the slice metric! (Otherwise there would be no light bending, since the slice metric is euclidean!) However, the Painleve chart is neither orthogonal nor harmonic.

Did you draw local light cones in the Painleve vs Schwarzschild charts? See http://physics.syr.edu/courses/modules/LIGHTCONE/ for local light cones in Minkowski vacuum. Plotting local light cones in your chart is a very good idea if you want to understand what your plots of geodesics as represented in that chart are trying to tell you! Try plotting in both a slice t=t0, and also plotting equatorial plane geodesics in an E^(1,2) picture of Painleve.

Did you look up harmonic charts and Kerr-Schild charts yet? These are very well worth looking into. See http://casa.colorado.edu/~ajsh/schw.shtml for some more useful charts, including Eddington and Kruskal-Szekeres.

 

[CarlB]

Painleve Coordinates, step 1, convert to Cartesian

Okay, the program of writing Newton and Schwarzschild in Cartesian (see post #11 for references to literature) coordinates without trig has been completed and the results test out nicely. It is time to do the same with Painleve coordinates:

$$ds^2 = -dt^2 + (dr-(2/r)^{0.5} dt)^2 + r^2 d\phi^2.$$

First step is to convert them to "Cartesian" form. Following post #2, we find:

$$\begin{array}{rcl} ds^2 &=& -dt^2 + ((xdx+ydy)r^{-1}-\sqrt{2}\;r^{-0.5} dt)^2 + (y^2dx^2 + x^2dy^2 -2xy\;dx\;dy)r^{-2},\\ &=&-dt^2 + dx^2 + dy^2 + 2r^{-1}dt^2 -2\sqrt{2}\;r^{-1.5}(x\;dx\;dt+y\;dy\;dt),\\ I = \left(\frac{ds}{dt}\right)^2 &=&\frac{2}{r}-1 + \dot{x}^2 + \dot{y}^2-2\sqrt{2}r^{-1.5}(x\dot{x}+y\dot{y}) \end{array}$$

The idea is to apply the Euler Lagrange equations to the integral over t for s:

$$s = \int \sqrt{\left(\frac{ds}{dt}\right)^2}\;dt = \int \sqrt{I}\;dt$$

to obtain the equations of motion in a Cartesian coordinate system. As supplied by pervect in post #22, the resulting equations of motion are:

$$\frac{d}{d t} \left[ \frac{1}{2 \sqrt{I}} \left( \frac{\partial I}{\partial \dot{x}} \right) \right] = \frac{1}{2 \sqrt{I}} \frac{\partial I}{\partial x}$$

 

[Chris Hillman]

Hi, Carl,

Okay, the program of writing Newton and Schwarzschild in Cartesian (see post #11 for references to literature) coordinates without trig has been completed and the results test out nicely... Chris, I really appreciate your help on this problem, but I don't have time to bring you up to speed when you can do it for yourself with little effort. If you are willing to swear that you've carefully read each post in the rest of the thread, and at least looked at the articles that were linked to it, then I will answer whatever questions you have left, and we won't clutter up the thread going over the same old same old. Actually, after looking through the thread it appears that I really should rewrite the secret reasons for working on this sort of thing.

Well, I don't think you are giving me enough credit for "reasonable effort" (I can see that for your part, you don't think I am giving you sufficient credit for knowledge and careful work, and you may be right!), but it seems that you have now recognized some of the expository deficiencies of your previous posts, and I applaud your attempt to rectify the lacunae you noticed.

By the way, if you said previously that you have a background in geometric algebra, I must have either missed or forgotten that comment. As it happens, many years ago I knew a fair amount about Cayley-Dickson algebras (see for example my posts on noneuclidean trigonometries, which will probably be familiar to you but not to many others here), and I do know the papers by Doran, so I think I can say with some confidence that it was not nearly so apparent as you thought that this was part of the background for your computations. (By the way, there is an important point about the Doran chart for the Kerr vacuum which we can discuss if you like, which illustrates my point about unanticipated subtleties which can easily trip up newbies.)

Also, if you had seen as many inaccurate plots drawn by enthusiastic computer jocks for fun, or as many accurate but terribly misinterpreted plots, as I have done, you would probably appreciate why it is in your own best interests to quickly convince sceptics that your plots do not fall into these categories!
But even if you don't appreciate why I am skeptical by default, please not that following up on my suggestions should be fun and if you have indeed made no errors of math or interpretation should greatly strengthen the presention of your plots.

OK, I hope this clears the air!

Let me try to continue by annotating Carl's post, saying some things which I thing should help others to follow along.

Painleve coordinates:
$$ds^2 = -dt^2 + (dr-(2/r)^{0.5} dt)^2 + r^2 d\phi^2.$$

The Painleve chart is one of the most important charts on the Schwarzschild vacuum solution, and I've been a big fan of it ever since I independently "discovered" it some time around 1984. (Much later I learned it was introduced by Painleve in 1921.)

Like the Eddington chart, it comes in ingoing and outgoing flavors and the ingoing (outgoing) Painleve chart covers the same region as the ingoing (outgoing) Eddington chart. We are interested in ingoing Painleve chart here because Carl is (I take it) planning to study the motion of freely and radially infalling observers who pass from the "right exterior" into the "future interior" region. The Painleve chart is in fact ideally adapted for this purpose. To obtain it from the exterior Schwarzschild chart
$$ds^2 = -(1-2m/r) \, dt^2 + \frac{dr^2}{1-2m/r} + r^2 \, \left( d\theta^2 + \sin(\theta)^2 \, d\phi^2 \right),$$
$$-\infty < t < \infty, \; 2 m < r < \infty, \; 0 < \theta < \pi, \; -\pi < \phi <\pi$$
put
$$T = t + \int \frac{\sqrt{2m/r}}{1-2m/r} \, dr = t + \sqrt{8 m r} - 4m \operatorname{arctanh} \left( \sqrt{\frac{r}{2 m}}\right)$$
which gives
$$ds^2 = -(1-2m/r) \, dT^2 + 2 \, \sqrt{2m/r} \, dT \, dr + dr^2 + r^2 \, \left( d\theta^2 + \sin(\theta)^2 \, d\phi^2 \right),$$
$$-\infty < T < \infty, \; 2m < r < \infty, \; 0 < \theta < \pi, \; -\pi < \phi <\pi$$

I'm too lazy to draw pictures, and I've told this story before, but let me try to give some indication of how I discovered this transformation before I knew any calculus (!). Looking at the pictures in MTW, my idea was to "pull down" the constant time slices by defining a new time coordinate such that the new constant time slices would be everywhere orthogonal to the world lines of the Lemaitre observers (whom I learned about from the little textbook by Dirac), in order to combine the best features of the ingoing Lemaitre and ingoing Eddington charts. You can actually do this if you know only the appropriate trignometry and trust a table of integrals :-/ and that's what I did.

Now we immediately notice two things. First, the domain can be extended from the right exterior into the "future interior" by increasing the range of the radial coordinate $0 < r < \infty$ (we'll need to interpret this carefully inside the horizon, but it still deserves the name of "Schwarzschild radial coordinate"). Second, the constant Painleve coordinate time slices $T=T_0$ have the metric of euclidean three-space! This was my first noteworthy "discovery" in gtr, so you can see why I am so fond of the Painleve chart!

With a bit of insight (or practice!) we can read off an "obvious" coframe field:
$$\sigma^0 = -dT, \; \sigma^1 = \sqrt{\frac{2m}{r}} \, dT + dr, \; \sigma^2 = r \, d\theta, \; \sigma^3 = r \, \sin(\theta) \, d\phi$$
whose dual frame field consists of the unit vector fields (first timelike, last three spacelike)
$$\vec{e}_0 = \partial_T - \sqrt{\frac{2m}{r}} \, \partial_r, \; \vec{e}_1 = \partial_r, \; \vec{e}_2 = \frac{1}{r} \, \partial_\theta, \; \vec{e}_3 = \frac{1}{r \, \sin(\theta)} \, \partial_\phi$$
This frame field is suitable for studying the physical experience of observers who fall in freely and radially "from rest at $r=\infty$". That is, the world lines of this class of observers are precisely the integral curves of the timelike unit vector field from our frame field, and the three spacelike vector fields then give the "spatial frame vectors". IOW, our frame field attaches the "local Lorentz frame" of these observers to each event. The point of introducing this frame here is to point out that in the Schwarzschild vacuum geometry, the timelike vectors $\vec{e}_0$ are everywhere orthogonal to the hyperslices $T=T_0$. IOW, in a sense, for these particular observers, we can consider their "space at a time" to be flat, but of course the spacetime itself is not flat. (A pedantic aside: the spatial frame vectors have vanishing Fermi derivatives, after projection into the orthogonal hyperplane elements, so this frame is in fact inertial nonspinning, which is the closest we can come in a curved Lorentzian manifold to a global Lorentz frame as in the standard formalism of str.

We can simultaneously employ alternative frame fields to study other families of observers, and this is in fact the best way to clear up the kind of confusion so often evident in discussions of alleged "paradoxes". I have done that in great detail elsewhere and it's off topic here, so moving on:

By the spherical symmetry, we can suppress one angular coordinate by restricting attention to the equatorial plane $\theta=\pi/2$. Also, the mass parameter can be set to one by rescaling coordinates, so let's do that too. This gives the line element Carl wrote down, after renaming $T \mapsto t$.

Going one step further, we can temporarily suppress $\phi$ also and draw the local light cones at various events using the frame vectors
$$\vec{e}_0 = \partial_T - \sqrt{\frac{2m}{r}} \, \partial_r, \; \vec{e}_1 = \partial_r$$
Switching to a more computer graphics friendly notation,
$$\vec{A}(T,r) = (1, -\sqrt{\frac{2m}{r}} ), \; \vec{B}(T,r) = (0,1)$$
Fixing event $P = (T_0,r_0)$, and pretending we are working in Minkowski spacetime (really, we are working in the tangent space at our event),
draw a segment to the event $P+\varepsilon \, (\vec{A}+\var{B})$,
then from that event to the event $P+\varepsilon \, (\vec{A}-\var{B})$,
then back to $P$. That's the top half of our local light cone. Draw the frame vectors to the same scale (i.e. scalar multiplied by $\varepsilon$) to see how our "local Lorentz frame" relates to the shape of our "local light cones". Now draw some more of these gadgets, to scale. Draw some in the exterior, some in the interior, and some on the horizon itself. Notice how the light cones progressively shear inward as the Schwarzschild radius decreases. This is really a way of visualizing how our local light cones, which are can be considered as a feature of the local "conformal structure" although perceptive readers will note that I am treating them as a "metrical structure", relate to the Killing vector field $\partial_t$, which is timelike on the right exterior (hence our solution is "static" there), but null on the horizon, and spacelike in the future interior.

Next, if we want to study general geodesics, it makes sense to write down the geodesic equations. The easiest way to do this is to read off the geodesic Lagrangian and then compute the Euler-Lagrange equations; putting the result in "monic" form we obtain the geodesic equations in their standard form, from which we can read off the Christoffel symbols, if we so desire. Alternatively, we can use the methods of Cartan (see MTW). I'll spare you the gory mess: suffice it to say that we get something quite a bit messier than what we obtain for the Schwarzschild exterior chart (see the discussion of effective potentials in MTW for a convenient method of analyzing the geodesics in the interior, in terms of the Schwarzschild exterior chart).

Now, anyone who did the "coordinate tutorial" I posted long ago to sci.physics.relativity knows very well that adopting the "right" coordinate chart can greatly simplify the analysis of the geodesics! As the active reader has now seen, while the Painleve chart has many virtues, it does not succeed in simplifying the geodesic equations if that is our goal; indeed, it makes them appear somewhat more complicated.

Another reason for exploring multiple charts is that Lie's symmetry methods can be much easier to crank through, which can enable us to discover "especially symmetrical" geodesics. In particular, if we didn't already know about the circular photon orbits at $r=3m$, we would discover them by a systematic study of this kind. Actually, Lie's methods are probably even more relevant for the problem of solving the wave equation, Klein-Gordon equation, Dirac equation, etc., on curved spacetimes, where changing charts can again have drastically beneficial effects, and where again this can be detected (even sometimes predicted) using a Lie theoretic symmetry analysis of the equation to be solved. This game is important for physics for various reasons; one which should sound familiar to many readers here is the study of perturbations, e.g. we might want to know how radiation of various kinds is "scattered" by a Schwarzschild object.

At this point, I'd highly recommend that anyone who hasn't much experience comparing say the wave equation in various charts should pause and collect a bunch of familiar charts for the Schwarzschild vacuum (paying careful attention to noting the coordinate ranges), such as Schwarzschild, static and spatially isotropic, Lemaitre, Eddington, Painleve, Kruskal-Szekeres, and figure out the transformations between them. Next review how to compute the wave operator using the elegant methods of Cartan (see MTW, or my long ago expository post to sci.math.relativity called "The Joy of Forms"), and compute for your charts. Check your results by applying the transformations you found earlier.

OK, so now let's confront the question: how can we find charts in which the geodesic equations, or the wave equation, or whatever, look "as simple as possible"? I don't think there's an algorithmic answer to that, but there are some basic things to try.

One approach is to try to "rationalize" our coordinates. For example, the usual polar spherical line element
$$ds^2 = dr^2 + r^2 \; \left( d\theta^2 + \sin(\theta)^2 \, d\phi^2 \right)$$
contains trig functions. We can rationalize this by setting $\eta=\cos(\theta)$, so that the line element becomes
$$ds^2 = dr^2 + r^2 \; \left( \frac{d\eta^2}{\eta^2} + (1-\eta^2) \, d\phi^2 \right)$$
which contains only rational expressions in the coordinates.

A particularly dramatic example of this arises in a much more general context which includes the solution under discussion: in studying the Ernst vacuums, a major breakthrough was the recognition that the "rational prolate spheroidal chart" results in a dramatic simplication of the Ernst equation (that is, the "master equation" in the set of equations to which we can reduce the vacuum EFE in studying the Weyl canonical chart for a stationary axisymmetric spacetime). If you're curious, the rational prolate spheroidal chart for E^3 has the line element
$$\frac{ds^2}{A^2} = (\chi^2-\xi^2) \, \left( \frac{d\chi^2}{\chi^2-1} + \frac{d\xi^2}{1-\xi^2} \right) + (\chi^2-1) \, (1-\xi^2) \, d\phi^2,$$
$$1 < \chi < \infty, \; -1 < \xi < 1, \; -\pi < \phi < \pi$$
where $A > 0$ is a scale factor. The transformation from the standard cylindrical chart is
$$z = A \, \chi \, \xi, \; r = A \, \sqrt{\chi^2-1} \, \sqrt{1-\xi^2}$$
Using the elegant exterior calculus methods of Cartan, you can verify with surprisingly little effort that the Laplace operator becomes, for an axisymmetric function,
$$\Delta = \frac{1}{A^2 \, (\chi^2-\xi^2)} \left( \partial_\chi (\chi^2-1) \partial_\chi + \partial_\xi (1-\xi^2) \partial_\xi \right)$$
Quite remarkably, this is symmetric under interchanging $\chi, \xi$, which Chandrasekhar took good advantage of in his work on the Ernst vacuums! Even better, you can readily attack this by the elementary method of separation of variables, obtaining a useful series expansion for asymptotically vanishing harmonic functions in terms of Legrendre polynomials/functions. This looks completely different from the analogous expression we'd obtain by analyzing the same equation in a cylindrical or polar spherical chart! (This is far more impressive for the Ernst equation, where as I said it turns out that the rational prolate spheroidal Weyl canonical chart on a stationary axisymmetric spacetime achieves several substantial simplifications of several interesting cases of the EFE for this situation, including vacuum, electrovacuum, minimally coupled massless scalar fields, dusts, and so on.)

Another approach is to try harmonic coordinates. These are charts in which the coordinates, considered as monotonic functions on the domain of our chart, are harmonic functions (in the "Laplace-Beltrami sense": they satisfy the curved spacetime wave equation).

Still another is try geometric algebra, which is apparently what Carl is working on, ultimately in the context of trying to simplify some equation other than the geodesic equation. Geometric algebra is a kind of elaboration of the formalisms of vector calculus, exterior calculus, and quaternionic calculus. So far it has been studied mainly in the UK, especially Cambridge, and I think Carl will agree that, like Hamilton's pursuit of his program of rewriting nineteenth century physics in terms of quaternions, geometric algebra is either ignored or considered by many observers to be of dubious value, given the intricacy of the notation. However, its enthusiastic proponents (and back in the day, I knew someone studying this stuff at Cambridge pretty well) can point to a few notable successes, including the Doran chart for the Kerr vacuum, which generalizes the Painleve chart (aha!) to the Kerr vacuum, which was obtained using the methods of geometric algebra. As a simpler example, it often happens when using frame fields that one wishes to "despin" a spinning frame by rotating by just the right angle about just the right axis of rotation (said angle and axis defined at each event) to kill the projected Fermi derivatives, and this kind of operation would indeed probalby be more conveniently carried out in geometric algebra, except in special cases. For example, the "obvious" coframe read off the Doran chart is in fact dual to a spinning frame, i.e. its timelike unit vector is indeed what we want, but the spacelike unit vectors need to be "despun" to obtain a nonspinning inertial frame analogous to the Lemaitre frame discussed above. However, this Doran-Lemaitre frame can in this case be obtained easily by elementary methods, once one has the Doran chart in hand, and I have little doubt that the Doran chart itself can be obtained without recourse to geometric algebra.

First step is to convert them to "Cartesian" form.

Looking back, I see that Carl originally referred to "Cartesian coordinates" for the Schwarzschild vacuum solution, which in the context of this forum, or elementary gtr generally, is terribly confusing because those in the know will fear that someone is trying to treat a curved spacetime using coordinates which can exist only in a flat spacetime! However, eventually it turned out that Carl just meant that he is applying a transformation given by a formula familiar from vector calculus:
$$x = r \, \sin(\theta) \, \cos(\phi), \; y = r \, \sin(\theta) \, \sin(\phi), \; z = r \, \cos(\theta)$$
Although of course we are working in part of the Schwarzschild vacuum, not Euclidean three-space, this makes perfect sense, so long as we are careful to restrict the ranges of the coordinates to ensure that we obtain a well defined diffeomorphism on the overlap of the domains (in this case, they are almost the same--- as an exercise, what is the "bad locus" which must be excised from the domain of the original Painleve chart but which we can reasonably hope to include in the variant Carl is switching to?)

Following post #2, we find:
$$\begin{array}{rcl} ds^2 &=& -dt^2 + ((xdx+ydy)r^{-1}-\sqrt{2}\;r^{-0.5} dt)^2 + (y^2dx^2 + x^2dy^2 -2xy\;dx\;dy)r^{-2},\\ &=&-dt^2 + dx^2 + dy^2 + 2r^{-1}dt^2 -2\sqrt{2}\;r^{-1.5}(x\;dx\;dt+y\;dy\;dt),\\ I = \left(\frac{ds}{dt}\right)^2 &=&\frac{2}{r}-1 + \dot{x}^2 + \dot{y}^2-2\sqrt{2}r^{-1.5}(x\dot{x}+y\dot{y}) \end{array}$$

The idea is to apply the Euler Lagrange equations to the integral over t for s:

$$s = \int \sqrt{\left(\frac{ds}{dt}\right)^2}\;dt = \int \sqrt{I}\;dt$$

to obtain the equations of motion in a Cartesian coordinate system.

In the standard "geodesic Lagrangian", the dots refer to differentiation with respect to an affine parameter, which for timelike geodesics we can take to be an arc length paramter. Here, however, Carl's dots refer to differentiation with respect to the Painleve time coordinate, which I am denoting by $T$. But recall that in the Painleve chart, the Painleve time coordinate happens to give the arc length parameterization of the Lemaitre congruence. The expressions Carl quotes (attributing them to pervect)
$$\frac{d}{d t} \left[ \frac{1}{2 \sqrt{I}} \left( \frac{\partial I}{\partial \dot{x}} \right) \right] = \frac{1}{2 \sqrt{I}} \frac{\partial I}{\partial x}$$
do agree with the Euler-Lagrange equations obtained via elementary variational calculus from the standard geodesic Lagrangian, except for what would be an inessential scale factor if I didn't depend upon time. Pervect?

 

[Chris Hillman]

Continuing (my post grew too long!):

the Dirac equation in a gravitational field is particularly simple in Painleve-Cartesian coordinates. This may not make a hill of beans to a GR specialist,

For other readers: I think Carl was being disingenuous (or sarcastic) here, but for the record, simplifying the Dirac equation on the Schwarzschild or Kerr vacuums, or other "favorite gtr models", is in fact a "minor industry" (meaning that these topics are of considerable interest to researchers working in the field of gravitation physics!)--- maybe that's the "industry" Carl had in mind when he mentioned this term in an earlier post.

But nevertheless, this is a subject of interest in physics at this time, and this is precisely the language that is used to describe it. The reason is sort of subtle and has to do with the fact that GR can be derived from a graviton theory on a flat space.

My first guess was that Carl was referring here to Deser's approach to obtaining a "local mimic" of gtr by starting with a naive quantization formulated on a flat spacetime background, then fixing up inconsistencies, noticing new inconsistencies, fixing those, and so on until, miraculously, one winds up with gtr. (This is discussed briefly in MTW, including a citation to the original paper.)

However,

An example of a recent paper in the literature using "Gullstrand-Painleve-Cartesian" coordinates is:
http://arxiv.org/abs/gr-qc/0411060
The above has 12 citations and gives easily understood reasons for looking at things this way that can be understood at a high-school level.

the Hamilton/Lisle eprint (I seem to have overlooked this--- thanks for bringing it to my attention!) seems, at a glance, to refer to the flat metric inherited by the constant Painleve time hyperslices. They write "picture space as flowing like a river into the Schwarzschild black hole", which suggests to me they are thinking of projecting $\vec{e}_0$ into $T=0$ to form a flow (in the sense of differential equations on manifolds).

(The world lines of the Doran-Lemaitre observers in the Kerr vacuum form an irrotational congruence, hence a hypersurface orthogonal congruence, but the orthogonal hyperslices no longer have vanishing three-dimensional Riemann tensor, so this presumably would not apply to the Kerr generalization.)

Continuing with Carl's program, his transformation yields the line element
$$ds^2 = -\left(1-\frac{2m}{\sqrt{x^2+y^2+z^2}} \right) \, dT^2 + 2 \, \frac{\sqrt{2m}}{(x^2+y^2+z^2)^{3/4}} \, \left( x \, dx + y \, dy + z \, dz \right) \; dT + dx^2 + dy^2 + dz^2,$$
$$-\infty < T, \, x, \, y, \, z < \infty, \; x^2+y^2+z^2 \neq 0$$
The Lemaitre frame field becomes:
$$\vec{e}_0 = \partial_T - \frac{\sqrt{2m}}{(x^2+y^2+z^2)^{3/4}} \; \left( x \, \partial_x + y \, \partial_y + z \, \partial_z \right)$$
$$\vec{e}_1 = \frac{1}{\sqrt{x^2+y^2+z^2}} \; \left( x \, \partial_x + y \, \partial_y + z \, \partial_z \right)$$
$$\vec{e}_2 = \frac{z}{\sqrt{x^2+y^2+z^2} \, \sqrt{x^2+y^2}} \; \left( x \, \partial_x + y \, \partial_y \right) - \frac{\sqrt{x^2+y^2}}{\sqrt{x^2+y^2+z^2}} \; \partial_z$$
$$\vec{e}_3 = \frac{1}{\sqrt{x^2+y^2}} \; \left( -y \partial_x + x \partial_y \right)$$
To see this, we simply write down the coordinate vectors of the old chart in terms of the new one and plug these into our previous expressions for the frame field. The easiest way to transform a coordinate vector is to simply apply it to the new coordinates, which are simply functions on our manifold! For example,
$$\partial_r z = \partial_r (r \cos(\theta)) = \cos(\theta) = \frac{z}{\sqrt{x^2+y^2+z^2}}$$
and so forth, so that
$$\partial_r = \frac{1}{\sqrt{x^2+y^2+z^2}} \; \left( x \, \partial_x + y \, \partial_y + z \, \partial_z \right)$$
and so on.

The geodesic equations turn out to be a bit messy. It's much easier to state what the four standard Killing vectors look like:
$$\partial_T, \; -y \partial_x + x \partial_y, \; -z \partial_y + y \partial_z, \; -x \partial_z + z \partial_x$$
We should in fact take this simple result as the defining characteristic of "the cartesian-type Painleve chart" for our static spherically symmetric vacuum solution-- maybe that was Carl's point in one of his remarks above.

The first two geodesic equations are
$$0 = \ddot{T} + \frac{\sqrt{2m^3}}{(x^2+y^2+z^2)^{5/4}} \, \dot{T}^2 + \frac{2 m \, \dot{T}}{(x^2+y^2+z^2)^{3/2}} \; \left( x \dot{x} + y \dot{y} + z \dot{z} \right)$$
$$\hspace{0.5in} -\frac{\sqrt{m/2}}{(x^2+y^2+z^2)^{7/4}} \; \left( (x^2-2\,(y^2+z^2)) \, \dot{x}^2 + (y^2-2\,(z^2+x^2)) \, \dot{y}^2 + (z^2-2\,(x^2+y^2)) \, \dot{z}^2 \right)$$
$$\hspace{0.5in} + \frac{3 \, \sqrt{2m}}{(x^2+y^2+z^2)^{7/4}} \; \left( x y \, \dot{x} \dot{y} + y z \, \dot{y} \dot{z} + z x \, \dot{z} \dot{x} \right)$$
and
$$0 = \ddot{x} - m x \, \frac{2m - \sqrt{x^2+y^2+z^2}}{(x^2+y^2+z^2)^2} \, \dot{T}^2 - \frac{2 \sqrt{2 m^3} \, x \, \dot{T}}{(x^2+y^2+z^2)^{9/4}} \; \left( x \dot{x} + y \dot{y} + z \dot{z} \right)$$
$$\hspace{0.5in} -\frac{m x}{(x^2+y^2+z^2)^{5/2}} \; \left( (x^2-2 \, (y^2+z^2)) \, \dot{x} + (y^2-2 \, (z^2+x^2)) \, \dot{y} + (z^2-2 \, (x^2+y^2)) \, \dot{z} \right)$$
$$\hspace{0.5in} - \frac{6 m x}{(x^2+y^2+z^2)^{5/2}} \; \left( x y \, \dot{x}\dot{y} + y z \, \dot{y}\dot{z} + z x \, \dot{z}\dot{x} \right)$$
(similarly for $\ddot{y}, \, \ddot{z}$), where the dots denote differentiation wrt the affine parameter. These expressions are not particularly simple, but they are at least symmetric between the coordinates x,y,z (as must be the case). As a check, note that in relativistic units, the right hand sides have the dimensions of inverse length.

If I understand correctly, Carl's intention is not simply to write down the geodesic equations and numerically integrate them and plot results in illuminating ways (I enthusiastically second his recommendation of Andrew Hamilton's website, BTW; see the link at RelWWW below), but to express affine parameterized null geodesics (or their "tracks" in a constant time slice) in closed form. Most textbooks show how to do this (using the Schwarzschild exterior chart) in the exterior region, but to my knowledge none discuss how to do it in the future interior region (using for example ingoing Eddington or Painleve charts).

As a check, hardy readers can compute the connection one-forms and curvature two-forms, read off the Riemann components with respect to the Lemaitre frame (don't forget that this is an anholonomic frame!), and compute the characteristic polynomial, which is of course an invariant. If you do this correctly, you will obtain
$$\lambda^6 - \frac{6 m^2}{(x^2+y^2+z^2)^3} \, \lambda^4 + \frac{4 m^3}{(x^2+y^2+z^2)^{9/2}} \, \lambda^3 + \frac{9 m^4}{(x^2+y^2+z^2)^6} \, \lambda^2$$
$$\hspace{0.5in} - \frac{12 m^5}{(x^2+y^2+z^2)^{15/2}} \, \lambda + \frac{4 m^6}{(x^2+y^2+z^2)^9}$$
Of course this should agree with the result of simplying plugging Carl's transformation into the characteristic computed in the original Schwarzschild chart, and it does! Similarly, the principle Lorentz invariants of the Riemann tensor are
$$R_{abcd} \, R^{abcd} = \frac{48 m^2}{(x^2+y^2+z^2)^3}, \; \; R_{abcd} \, (\star R)^{abcd} = 0$$
(no intrinsic gravitomagnetism in the sense of Ciufolini and Wheeler).

By the way, it is easy to rationalize the original Painleve chart. Set
$$\tau=T, \; \rho = \sqrt{r}, \; \eta = \sin(\theta)$$
Writing $M=\sqrt{2m}$, the line element becomes
$$ds^2 = -(1-M^2/\rho^2) \, d\tau^2 + 4 M \, d\tau d\rho + 4\rho^2 \, d\rho^2 + \rho^4 \; \left( \frac{d\eta^2}{1-\eta^2} + \eta^2 \, d\phi^2 \right),$$
$$-\infty < \tau < \infty, \; 0 < \rho < \infty, \; 0 < \eta < 1, \; -\pi < \phi < \pi$$
The geodesic equations become
$$\ddot{\tau} + \frac{2 M^2}{\rho^3} \, \dot{\tau} \dot{\rho} + \frac{M^3}{2 \rho^5} \, \dot{\tau}^2 + \frac{2M}{\rho} \, \dot{\rho}^2 - M \rho \; \left( \frac{\dot{\eta}^2}{1-\eta^2} + \eta^2 \dot{\phi}^2 \right) = 0$$
$$\ddot{\rho} - \frac{M^3}{\rho^5} \, \dot{\tau} \dot{\rho} + (\rho^2-M^2) \; \left( \frac{M^2 \, \dot{\tau}^2}{4 \rho^7} + \frac{\dot{\rho}^2}{\rho^3} + \frac{1}{2 \rho} \; \left( \frac{\dot{\eta}^2}{1-\eta^2} + \frac{\dot{\phi}^2}{\eta^2} \right) \right) = 0$$
$$\ddot{\eta} + \frac{4}{\rho} \, \dot{\rho} \dot{\eta} + \eta \; \left( \frac{\dot{\eta}^2}{1-\eta^2} - (1-\eta)^2 \, \dot{\phi}^2 \right) = 0$$
$$\ddot{\phi} + 2 \, \dot{\phi} \; \left( \frac{2 \dot{\rho}}{\rho} + \frac{\dot{\eta}}{\eta} \right) = 0$$
In the second equation, note that the horizon is located at $\rho=M$, so this equation simplifies greatly there!

In the equatorial plane $\eta=1$, the geodesic equations become
$$\ddot{\tau} + \frac{2 M^2}{\rho^3} \, \dot{\tau} \dot{\rho} + \frac{M^3}{2 \rho^5} \, \dot{\tau}^2 + \frac{2M}{\rho} \, \dot{\rho}^2 - M \rho \, \dot{\phi}^2 \right) = 0$$
$$\ddot{\rho} - \frac{M^3}{\rho^5} \, \dot{\tau} \dot{\rho} + (\rho^2-M^2) \; \left( \frac{M^2 \, \dot{\tau}^2}{4 \rho^7} + \frac{\dot{\rho}^2}{\rho^3} + \frac{\dot{\phi}^2}{2 \rho} \right) \right) = 0$$
$$\ddot{\phi} + \frac{4}{\rho} \, \dot{\rho} \dot{\phi} = 0$$
where the last immediately yields the first integral
$$\dot{\phi} = \frac{L}{\rho^4}$$
Now plug this into the first two equations. We obtain a coupled system which is first order in $\dot{\tau}, \dot{\rho}$ provided that we pretend to forget the relationship between $\rho$ and $\dot{\rho}$. Thus, near some fixed locus of constant "rootradius", we can apply a phase plane analysis.

As this shows, for analyzing or plotting geodesics, it can pay to use several coordinate charts, e.g. using this one to numerically integrate the geodesic equations, then transforming to another chart, such as the "cartesian-type Painleve chart", to plot the results.

Hope this helps others to follow along--- speak up, lurking students, if you want me to continue trying to provide background as above!

 

[CarlB]

Now where was I? The equations of motion for Painleve-Cartesian coordinates. As usual, be warned I make LOTS of mistakes in calculation, and corrections are appreciated (and the reason I do this publicly).

$$\frac{1}{2\sqrt{I}} \frac{\partial I}{\partial x} =\frac{d}{d t} \left[ \frac{1}{2\sqrt{I}} \left( \frac{\partial I}{\partial \dot{x}} \right) \right],$$ or

$$\frac{I^{-0.5}}{2} \frac{\partial I}{\partial x} =$$
$$\frac{I^{-0.5}}{2}\left( \frac{\partial^2 I}{\partial\dot{x}\partial x}\dot{x}+ \frac{\partial^2 I}{\partial\dot{x}\partial y}\dot{y}+ \frac{\partial^2 I}{\partial\dot{x}^2}\ddot{x}+ \frac{\partial^2 I}{\partial\dot{x}\partial\dot{y}}\ddot{y}\right)$$
$$-\frac{I^{-1.5}}{4}\left( \frac{\partial I}{\partial x}\frac{\partial I}{\partial\dot{x}}\dot{x}+ \frac{\partial I}{\partial y}\frac{\partial I}{\partial\dot{x}}\dot{y}+ \frac{\partial I}{\partial \dot{x}}\frac{\partial I}{\partial\dot{x}}\ddot{x}+ \frac{\partial I}{\partial \dot{y}}\frac{\partial I}{\partial\dot{x}}\ddot{y} \right)$$

I want to get rid of the sqrt(I) and the 4s, and rearrange the terms because we need to solve for the 2nd derivatives in terms of the 1st derivatives and constants. Moving the 2nd derivatives to the right hand side we get:

$$2I\frac{\partial I}{\partial x} + \left(\frac{\partial I}{\partial x}\frac{\partial I}{\partial\dot{x}}- 2I\frac{\partial^2 I}{\partial\dot{x}\partial x}\right)\dot{x} + \left(\frac{\partial I}{\partial y}\frac{\partial I}{\partial\dot{x}}- 2I\frac{\partial^2 I}{\partial\dot{x}\partial y}\right)\dot{y} =$$

$$\left(2I\frac{\partial^2 I}{\partial\dot{x}^2}- \frac{\partial I}{\partial \dot{x}}\frac{\partial I}{\partial\dot{x}} \right)\ddot{x}+ \left(2I\frac{\partial^2 I}{\partial\dot{x}\partial\dot{y}}- \frac{\partial I}{\partial \dot{y}}\frac{\partial I}{\partial\dot{x}} \right)\ddot{y}$$

Similarly, the Euler Lagrange equation in y gives:

$$2I\frac{\partial I}{\partial y} + \left(\frac{\partial I}{\partial y}\frac{\partial I}{\partial\dot{y}}- 2I\frac{\partial^2 I}{\partial\dot{y}\partial y}\right)\dot{y} + \left(\frac{\partial I}{\partial x}\frac{\partial I}{\partial\dot{y}}- 2I\frac{\partial^2 I}{\partial\dot{y}\partial x}\right)\dot{x} =$$

$$\left(2I\frac{\partial^2 I}{\partial\dot{y}^2}- \frac{\partial I}{\partial \dot{y}}\frac{\partial I}{\partial\dot{y}} \right)\ddot{y}+ \left(2I\frac{\partial^2 I}{\partial\dot{y}\partial\dot{x}}- \frac{\partial I}{\partial \dot{x}}\frac{\partial I}{\partial\dot{y}} \right)\ddot{x}$$

The above equations mix the 2nd deriviatives. To use them, we need to solve for $$\ddot{x}$$ and $$\ddot{y}$$. The equations are linear so they can be written in 2x2 matrix form. To solve the equations requires that we invert the matrix, which can be done provided its determinant is not zero. The determinant is:

$$\left(2I\frac{\partial^2 I}{\partial\dot{x}^2}-\frac{\partial I}{\partial \dot{x}}\frac{\partial I}{\partial\dot{x}}\right)\left(2I\frac{\partial^2 I}{\partial\dot{y}^2}-\frac{\partial I}{\partial \dot{y}}\frac{\partial I}{\partial\dot{y}}\right) -\left(2I\frac{\partial^2 I}{\partial\dot{x}\partial\dot{y}}-\frac{\partial I}{\partial \dot{y}}\frac{\partial I}{\partial\dot{x}}\right)\left(2I\frac{\partial^2 I}{\partial\dot{y}\partial\dot{x}}-\frac{\partial I}{\partial \dot{x}}\frac{\partial I}{\partial\dot{y}}\right)$$

Next to calculate the various derivatives for the Painleve coordinates.

 

[CarlB]

The first order partial derivatives of I are as follows:
$$\frac{\partial I}{\partial x} = -2r^{-3}x - \sqrt{8}r^{-1.5}\dot{x} +\sqrt{18}r^{-3.5}x(x\dot{x}+y\dot{y}),$$
$$\frac{\partial I}{\partial y} = -2r^{-3}y - \sqrt{8}r^{-1.5}\dot{y} +\sqrt{18}r^{-3.5}y(x\dot{x}+y\dot{y}),$$
$$\frac{\partial I}{\partial \dot{x}} = 2\dot{x}-\sqrt{8}r^{-1.5}x,$$
$$\frac{\partial I}{\partial \dot{y}} = 2\dot{y}-\sqrt{8}r^{-1.5}y$$

The required second order partial derivatives are:

$$\frac{\partial^2 I}{\partial x\partial\dot{x}} = -\sqrt{8}r^{-1.5} + \sqrt{18}r^{-3.5}x^2$$
$$\frac{\partial^2 I}{\partial x\partial\dot{y}} = \sqrt{18}r^{-3.5}xy$$
$$\frac{\partial^2 I}{\partial y\partial\dot{x}} = \sqrt{18}r^{-3.5}xy$$
$$\frac{\partial^2 I}{\partial y\partial\dot{y}} = -\sqrt{8}r^{-1.5} + \sqrt{18}r^{-3.5}y^2$$

$$\frac{\partial^2 I}{\partial \dot{x}^2} = 2$$
$$\frac{\partial^2 I}{\partial \dot{y}^2} = 2$$
$$\frac{\partial^2 I}{\partial \dot{x}\partial\dot{y}} = 0$$

And the value of I was:
$$I =\frac{2}{r}-1 + \dot{x}^2 + \dot{y}^2-\sqrt{8}r^{-1.5}(x\dot{x}+y\dot{y})$$

The 2nd partials are fairly simple. Substituting them into the formula for the determinant gives:

$$\left(4I-\frac{\partial I}{\partial \dot{x}}\frac{\partial I}{\partial\dot{x}}\right)\left(4I-\frac{\partial I}{\partial \dot{y}}\frac{\partial I}{\partial\dot{y}}\right) - \left(\frac{\partial I}{\partial \dot{y}}\frac{\partial I}{\partial\dot{x}}\right)\left(\frac{\partial I}{\partial \dot{x}}\frac{\partial I}{\partial\dot{y}}\right)$$

Dividing by 4, this factors to:

$$I\left(4I - \left(\frac{\partial I}{\partial \dot{x}}\frac{\partial I}{\partial \dot{y}}\right)^2\right)$$

which is zero when I=0 or when

$$I = \frac{1}{4}\left(\frac{\partial I}{\partial \dot{x}}\frac{\partial I}{\partial \dot{y}}\right)^2$$

This is nice and simple, and I suspect it's correct. I should mention that when I last messed around with this, I ended up with very ugly Painleve calculations. MAXIMA barfed on it. This time I'm doing it by hand and it seems to be going better. Also, I was fixated at breaking out the singularity at the origin and I don't think that helped the calculation.

Putting (x,y) = (2,0) so that r=2 to see what happens on the event horizon eventually seems to lead to a result that I only end up with a zero determinant there for a light ray exiting the event horizon. I hope that I can divide these zeroes out of the equations of motion so that I'll have equations of motion that work for all particles, massless or massive (and perhaps also tachyonic), everywhere.

But before doing that, I'll go ahead and code up the above equations into the Schwarzschild applet to see if I can get motion that looks correct.

What I would really like would be motion that follows the Schwarzschild orbits but with different t so that the particles can penetrate the event horizon. However, to do this, I have to adjust the initial conditions between Painleve and Schwarzschild. I believe that initial positions are compatible with this hope, but that initial velocities depend on the difference in time. As usual, velocity here means dx/dt where t is parameter time, not proper time.

 

[CarlB]

As usual, be warned I make LOTS of mistakes in calculation, and corrections are appreciated (and the reason I do this publicly).

As usual, I find at least one mistake in the above equations. The square root in the Painleve coordinates had the wrong sign. The correct sign for the Painleve coordinates is:

$$ds^2 =(\frac{2}{r}-1)dt^2 + dx^2 + dy^2 + \sqrt{8}r^{-1.5}(x\;dx\;dt+y\;dy\;dt)$$

This changes the sign of the square root here:

$$I =\frac{2}{r}-1 + \dot{x}^2 + \dot{y}^2+\sqrt{8}r^{-1.5}(x\dot{x}+y\dot{y})$$

and everywhere else the sqrt(8) or sqrt(18) appears.

But on correcting this, and putting the equations into my Java differential equation integrator, the behavior is fairly crazy, so I suspect that there is still another problem.

I've received some PMs asking about the theoretical justification for this method of calculating geodesics. The calculation follows a problem from Kip Thorne's 2006-2007 class:

Exercise 24.5 Problem: Action principle for geodesic motion
Show, by introducing a specific but arbitrary coordinate system, that among all timelike world lines that a particle could take to get from event $$P_0$$ to $$P_1$$, the one or ones whose proper time lapse is stationary under small variations of path are the free-fall geodesics. In other words, an action principle for a timelike geodesic $$P(\lambda)$$ [i.e., $$x^\alpha(\lambda)$$ in any coordinate system $$x^\alpha$$ ] is

$$\delta\int_{P_0}^{P_1}d\tau = \int_0^1\left(g_{\alpha\beta}\frac{dx^\alpha}{d\lambda}\frac{dx^\beta}{d\lambda}\ \right)^{1/2} d\lambda = 0,$$

where $$\lambda$$ is an arbitrary parameter which, by definition, ranges from 0 at $$P_0$$ to 1 at $$P_1$$. [Note: unless, after the variation, you choose the arbitrary parameter $$\lambda$$ to be ``affine'' ($$\lambda = a\tau+b$$ where $$a$$ and $$b$$ are constants), your equation for $$d^2x^\alpha/d\lambda^2$$ will not look quite like (24.26).]

http://www.pma.caltech.edu/Courses/ph136/yr2004/0424.1.K.pdf
http://www.pma.caltech.edu/Courses/ph136/yr2006/

In my exercise, $$\lambda$$ is coordinate time, and the difference in time between the beginning and end of coordinate time is 1.

$$I = g_{\alpha\beta}\frac{dx^\alpha}{d\lambda}\frac{dx^\beta}{d\lambda} = \left(\frac{ds}{dt}\right)^2$$

Carl

 

[Chris Hillman]

Just trying to help, as you requested...

Hi all, I am the PF member whom CarlB mentioned, although I don't think he understood my PMs very well.

As usual, I find at least one mistake in the above equations...wrong sign... putting the equations into my Java differential equation integrator, the behavior is fairly crazy, so I suspect that there is still another problem.

Well, you did ask in the thread for help checking your work, and it seems you could use some. Do you still want my input?

I feel that this series of threads has been unneccessarily confusing, due to the fact that Carl

1. has declined every opportunity to outline his motivation for seeming to do something simple in a peculiar and difficult way,

2. frequently uses standard terms to mean something different, without bothering to warn the readers whose help he requested,

3. at times appears to use two letters to mean the same quantity, and also to use one letter to mean two different quantities,

4. in various other ways writes with minimal clarity.

At this point, I am tiring of trying to help out, but for the benefit of other PF members, let me try to add a bit more clarification before I duck out.

Carl says he is trying to compute "the geodesic equations". Used without qualification, this phrase always means what all the textbooks say it means, indeed what Kip Thorne's notes which he just cited say it means. See (24.26) of Thorne's notes.

Carl also says he is trying to compute them "the Euler-Lagrange way". To those who have studied variational calculus, this can only mean "write down a Lagrangian and apply the Euler-Lagrange operator". This E-L operator has the schematic form
$$\frac{d}{d\lambda} \frac{\partial}{\partial \dot{q}} - \frac{\partial}{\partial q}$$
(See Peter J. Olver, Applications of Lie Groups to Differential Equations, 2nd ed. Springer-Verlag, 1993 for a discussion of the Euler-Lagrange operator in connection with the general Noether symmetry principle). To those who have studied gtr, the Lagrangian one expects is what MTW call in Box 13.3 the "dynamical Lagrangian", which is read right off the line element by replacing $dx$ by $\dot{x}$ and so on. Applying the E-L operator then yields a system of second order ODEs which can be written in the form
$$\ddot{x^i} + \left( \operatorname{quadratic} \, \operatorname{combination} \dot{x^j}, \dot{x^k} \right) = 0$$
This is very easy and efficient; see Box 14.4 for a worked example.

It is crucial to understand that in this standard approach, "dot" refers to differentiation with respect to an affine parameter for the unknown curve.

Unfortunately, it turns out that Carl is not seeking the geodesic equations at all. (Too bad, because I took the trouble to write them out for his "Cartesian-type" Painleve chart in the middle of post #51 in this thread.) Furthermore, he is not using the standard "E-L" way, nor is he using some alternative methods some of us guessed he might be working with. After fifty posts, one really ought to be to expect some sense of what he is actually trying to do, but unfortunately, I still am not very confident I understand what he is up to, still less whether he is likely to eventually shed a little light, after blowing so much smoke over what he called, in one of the posts above, his "secrets".

FWIW, it now comes to light that Carl is apparently seeking not "the geodesic equations" but something I'll call "shape equations", which can often be obtained fairly easily from the geodesic equations. It's probably best to explain by example:

Consider the line element of the hyperbolic plane in the upper half space chart familiar from elementary complex analysis:
$$ds^2 = \frac{dx^2 + dy^2}{y^2}, \; -\infty < x < \infty, \; 0 < y < \infty$$
Now the unqualified phrase "the geodesic equations" can only mean one thing:
$$\ddot{x} - \frac{ 2 \, \dot{x} \, \dot{y} } {y} = 0, \; \ddot{y} + \frac{ \dot{x}^2 - \dot{y}^2 }{y} = 0$$
These are easily obtained by reading off what MTW call the dynamic Lagrangian from the line element
$$L = \frac{\dot{x}^2 + \dot{y}^2}{y^2}$$
and applying the operator
$$\left( \frac{d}{d\lambda} \frac{\partial}{\partial \dot{x}} - \frac{\partial}{\partial x}, \, \frac{d}{d\lambda} \frac{\partial}{\partial \dot{y}} - \frac{\partial}{\partial y} \right)$$
to L, which you can check gives the equations I wrote above, after making the second order derivatives "monic". The solution to these equations gives the geodesics of H^2 as affine parameterized curves in H^2.

Sometimes one isn't interested in any parameterization, but just the "shape" of the curves. In our example, we can obtain a differential equation for the "shape" y(x) from the geodesic equations as follows:

First, note that we immediately obtain a first integral from the "logarithmic" form of the first equation: $\dot{x} = A \, y^2$. Next, we can "force" an arc length parameter by writing
$$1 = \frac{\dot{x}^2 + \dot{y}^2}{y^2}$$
whence
$$\dot{y}^2 = y^2 \, \left( 1-A^2 \, y^2 \right)$$
Now we have

(EDIT: removed pseudolatex tags due to apparent cache exhaustion)

\frac{dx}{dy} = \frac{A \, y}{\sqrt{1-A^2 y^2}
Solutions of this equation give the "shape" of the geodesics, and in this case the shape has a nice description: the geodesics look like semicircles orthogonal to the locus y=0.

(EDIT: subsequent posters in this thread look out!--- we might have reached the maximal number of pseudo-latex formatted expressions in this thread.)

The point is that it seems that Carl seeks something like "shape equations", not what are called "the geodesic equations" in all the textbooks, including the course notes he himself cited. Indeed, he apparently seeks to differentiate with respect to the Painleve coordinate time, which I am writing T to avoid confusing with the Schwarzschild time t (the other coordinates are the very same nonconstant monotonic functions on our Lorentzian manifold, so we can use the same letters without any possibility of confusion). Just to add to the opportunity for confusion, the Painleve time just so happens to be an arc length parameter for certain timelike geodesics, namely the integral curves of the Lemaitre congruence, i.e. the world lines of observers who fall in freely and radially "from rest at infinity". The point of Box 13.3 in MTW is that solutions of the geodesic equations (obtained as above from the dynamical Lagrangian) are always affine parameterized curves. A slightly tricky point in Lorentzian manifolds: affine parameters make perfect sense for null geodesics, but we can't turn these into arc length parameters, as we can for timelike or spacelike geodesics parameterized by an affine parameter! At times, Carl seems to think that affine parameters are useless for null geodesics, but this is quite wrong: the geodesic equations work perfectly well for describing all geodesics, including null geodesics.

Another possible complication is that Carl is suppressing one spatial coordinate (reasonable because of the spherical symmetry, and a good idea if you want to make a spacetime plot in the end), but at times he also seems to be projecting into a constant Painleve time slice. The eprint by Hamilton and Lisle which he cited visualizes Kerr geometry by studying how the local light cones (see my post #50 above) are sheared from Minkowski background to the Painleve cones, and also how a gyrostabilized spatial frame carried by Doran observers (the analog of the Lemaitre observers) appear to spin in the Doran chart (the analog of the Painleve chart) as the observers fall in. They express this as a combination of Galilei transformation (for the shearing) plus a twist (for the spinning). This is a just a novel way of expressing, for the Kerr vacuum (and so far only for the Kerr vacuum!), the common phenomenon I often mention: gyrostabilized frames often appear to spin in a given coordinate chart. In particular, I pointed out long ago that the obvious frame in the original paper presenting the Doran chart is in fact spinning, but is easily "despun" because the frame is actually spinning about one of the spatial vectors, the one aligned with \partial_\phi.

Let me finish up by saying that new ways to visualize Schwarzschild and Kerr are constantly appearing, and the more the merrier, say I, subject to one important injunction: make no errors. I have repeatedly tried to suggest to Carl that he tackle some simpler examples first and compare results with those obtained using standard methods. Being systematic is more important, not less so, if one is prone to making errors. Published papers on visualizations so far tend to be OK, in my experience (but then I don't read trashcan journals), but there are quite a few websites out there which give offer alleged visualizations which are in fact erroneous, and we don't want to add to that plague.

 

[Kea]

Um, Chris, regarding "secrets" ... note the use of inverted commas here. You are quite free to read more of Carl's writing, widely available on the crackpot web. And of course we know about Baez's website.

Just trying to help.

 

[CarlB]

Carl says he is trying to compute "the geodesic equations". Used without qualification, this phrase always means what all the textbooks say it means, indeed what Kip Thorne's notes which he just cited say it means.

"Geodesic equations" does carry a meaning in physics as Chris states. However, if he will bother to read the thread, he will find that I never use the phrase "geodesic equations" to describe what I am looking for.

Most of the rest of the post is based on this failure to carefully read what's been written.

Carl also says he is trying to compute them "the Euler-Lagrange way". To those who have studied variational calculus, this can only mean "write down a Lagrangian and apply the Euler-Lagrange operator".

It is crucial to understand that in this standard approach, "dot" refers to differentiation with respect to an affine parameter for the unknown curve.

Anyone who wants to search for the "Euler-Lagrange" method will find that it does not imply the use of a Lagrangian. Anyone who reads through my posts will find that I have not claimed that I am using a Lagrangian, and have specifically noted otherwise in several places.

Here's a link to Mathworld's definition of Euler-Lagrange:
http://mathworld.wolfram.com/Euler-LagrangeDifferentialEquation.html

If you read the above link, you will find that (a) they use "L" to refer to the function being minimized, and (b) they never mention the word "Lagrangian". I don't know how Chris could have got the impression that Euler-Lagrange is a technique that only works with Lagrangians, and I certainly don't see how he can justify this comment. Furthermore, you will find that they use the "dot" notation to refer to the derivatives with respect to the integration variable, which they label as "t". I've followed these conventions precisely.

After fifty posts, one really ought to be to expect some sense of what he is actually trying to do ...

Really? You mean that someone really couldn't understand the first post I made on this? Let me quote myself:

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.

What part of "vary the path" is confusing here? Most of the rest of Chris' post just continues to boringly beat the same drum over and over. Yes, we all know how to calculate geodesic equations the standard way. So what.

Before I forget, this concept of writing GR in Cartesian coordinates is not a heresy that I invented. Here's a handout for students at U. Colorado:

ASTR 3740. Spring 2004.
Using the River Model to Draw Geodesics around Black Holes

Fancy creating a computer program, maybe a Java applet, to draw the orbits of particles or photons around black holes? The river model provides a nice way to implement this.

...

Rewrite the river metric (1) in Cartesian coordinates $$x^\mu \equiv (x^0,x^1,x^2,x^3) \equiv (t_{ff},x,y,z)$$ with origin at the center of the black hole:

$$ds^2 = \eta_{\mu\nu}(dx^\mu-v^\mu dt_{ff}(dx^\nu - v^\nu dt_{ff})$$

where $$\eta_{\mu\nu} \equiv$$ diag(-1,1,1,1) is the Minkowski metric and $$v^\mu$$ are the components of the river velocity

$$v^\mu = v(0,-x/r,-y/r,-z/r)$$

http://casa.colorado.edu/~ajsh/astr3740_04/rivergeodesics.ps

Another possible complication is that Carl is suppressing one spatial coordinate (reasonable because of the spherical symmetry, and a good idea if you want to make a spacetime plot in the end), but at times he also seems to be projecting into a constant Painleve time slice.

Cool, a possible error. Could you point it out? By the way, when I'm done, symmetry will indicate how to modify the equations for three spatial dimensions, so there's no reason to carry them around now.

Carl

Update: Looking for the problem. I started plotting the acceleration as a function of position and velocity for the Newton, Schwarzschild, and Painleve coordinates. I was just a little surprised that the calculated acceleration (i.e. Cartesian acceleration with respect to coordinate time) is identical between the Schwarzschild and Painleve metrics for particles that have no velocity (again, in the Cartesian coordinates). The plot does cool things inside the event horizon.

Where the accelerations differ is when the test particle is moving. When I check phase space points with slight velocities, the accelerations diverge quite a bit. Looking at particle speeds from -0.009 to +0.007. Turns out Einstein acceleration changes by 0.000436, but my version of Painleve acceleration changes by 0.2366, probably waaaaay too big. Very likely I'm missing some divisions by r in terms that depend on velocity.

 

[Chris Hillman]

Well, I tried...

Before I forget, this concept of writing GR in Cartesian coordinates is not a heresy that I invented. Here's a handout for students at U. Colorado:
http://casa.colorado.edu/~ajsh/astr3740_04/rivergeodesics.ps

Carl, maybe you should ask yourself why I understand and enjoy reading AH's papers, but despite my good faith effort, I have been unable to avoid an exponential growth of confusion and misery in this thread. Apparently I am not the only one who has similar experiences trying to help you "check your work" or whatever it is you are here for.

Kea, I didn't get the joke, but I'll leave you and Carl to it. Good luck.

 

[CarlB]

VICTORY!

Painleve coordinates are working. The last problem was that my Java was missing part of a term. The orbits are gorgeous, and match the Schwarzschild orbits I earlier computed.

[edit] The new simulation is here:
http://www.gaugegravity.com/testapplet/SweetGravity.html [/edit]

The only thing I don't like about them is that the Schwarzschild particles don't follow down the same paths as the Painleve particles even though they have the same initial conditions. This is just due to the difference in time for the two coordinate systems. To fix it, I can use the differences in proper time / coordinate time between the two metrics to make their velocities measure up the same. I might add a switch box to do that.

The other thing that is ugly is that I haven't tried to reduce the equations to lower terms. Instead, the equations of motion are:

$$\left(\frac{ds}{dt}\right)^2 = I =\frac{2}{r}-1 + \dot{x}^2 + \dot{y}^2+\sqrt{8}r^{-1.5}(x\dot{x}+y\dot{y})$$

$$\frac{\partial I}{\partial x} = -2r^{-3}x + \sqrt{8}r^{-1.5}\dot{x} -\sqrt{18}r^{-3.5}x(x\dot{x}+y\dot{y}),$$
$$\frac{\partial I}{\partial y} = -2r^{-3}y + \sqrt{8}r^{-1.5}\dot{y} -\sqrt{18}r^{-3.5}y(x\dot{x}+y\dot{y}),$$
$$\frac{\partial I}{\partial \dot{x}} = 2\dot{x}+\sqrt{8}r^{-1.5}x,$$
$$\frac{\partial I}{\partial \dot{y}} = 2\dot{y}+\sqrt{8}r^{-1.5}y$$

The required second order partial derivatives are:

$$\frac{\partial^2 I}{\partial x\partial\dot{x}} = +\sqrt{8}r^{-1.5} - \sqrt{18}r^{-3.5}x^2$$
$$\frac{\partial^2 I}{\partial x\partial\dot{y}} = -\sqrt{18}r^{-3.5}xy$$
$$\frac{\partial^2 I}{\partial y\partial\dot{x}} = -\sqrt{18}r^{-3.5}xy$$
$$\frac{\partial^2 I}{\partial y\partial\dot{y}} = +\sqrt{8}r^{-1.5} - \sqrt{18}r^{-3.5}y^2$$

$$\frac{\partial^2 I}{\partial \dot{x}^2} = 2$$
$$\frac{\partial^2 I}{\partial \dot{y}^2} = 2$$
$$\frac{\partial^2 I}{\partial \dot{x}\partial\dot{y}} = 0$$

The equations of motion are then defined by computing the following six functions of phase space:

$$A = \left(2I\frac{\partial^2 I}{\partial\dot{x}^2}- \frac{\partial I}{\partial \dot{x}}\frac{\partial I}{\partial\dot{x}} \right)$$

$$B = \left(2I\frac{\partial^2 I}{\partial\dot{x}\partial\dot{y}}- \frac{\partial I}{\partial \dot{y}}\frac{\partial I}{\partial\dot{x}} \right)$$

$$C = 2I\frac{\partial I}{\partial x} + \left(\frac{\partial I}{\partial x}\frac{\partial I}{\partial\dot{x}}- 2I\frac{\partial^2 I}{\partial\dot{x}\partial x}\right)\dot{x} + \left(\frac{\partial I}{\partial y}\frac{\partial I}{\partial\dot{x}}- 2I\frac{\partial^2 I}{\partial\dot{x}\partial y}\right)\dot{y}$$

$$D = \left(2I\frac{\partial^2 I}{\partial\dot{y}\partial\dot{x}}- \frac{\partial I}{\partial \dot{x}}\frac{\partial I}{\partial\dot{y}} \right)$$

$$E = \left(2I\frac{\partial^2 I}{\partial\dot{y}^2}- \frac{\partial I}{\partial \dot{y}}\frac{\partial I}{\partial\dot{y}} \right)$$

$$F = 2I\frac{\partial I}{\partial y} + \left(\frac{\partial I}{\partial y}\frac{\partial I}{\partial\dot{y}}- 2I\frac{\partial^2 I}{\partial\dot{y}\partial y}\right)\dot{y} + \left(\frac{\partial I}{\partial x}\frac{\partial I}{\partial\dot{y}}- 2I\frac{\partial^2 I}{\partial\dot{y}\partial x}\right)\dot{x}$$

and combining them as follows:

$$\frac{d^2x}{dt^2} = (EC-BF)/(AE-BD),$$

$$\frac{d^2y}{dt^2} = (AF-DC)/(AE-BD).$$

This eliminates the need to keep track of proper time of particles when computing simultaneously multiple orbits of relativistic particles in the Painleve metric. I expect that I will be able to simplify it with MAXIMA, but the important thing is that it works.

As an aside, I almost wasted my time programming up the usual geodesic equations. In fact, I actually sat down to do it. What stopped me was the realization of how much more difficult the usual geodesic equations are to use than the simple equations of motion I've found here.

The geodesic equations require that you keep track of the particle's "proper time". But this won't work for photons, so you have to have a different algorithm for massless particles (and your algorithm for massive particles probably blows up when faced with very relativistic massive particles). But the equations I've got cleanly do the job for all particles, massive, massless and tachyonic.

The phase space I'm using (after you add back in the z dimension), has 6 dimensions, x,y,z and their velocities. If you do it with the geodesic equations, you will have 8 dimensions, x, y, z, and t and their velocities with respect to s. This means that in the case of the usual geodesic equations, you have 33% more variables to keep track of, and in addition you need to properly initialize the dt/ds variable. In my equations of motion (note I don't call them "geodesic equations"), all you do is plug and go.

So when I remembered that I'd have to deal with initialization hassle, (which I once did have in my Schwarzschild simulation), plus the hassle of having to a given number of position iterations not correspond to always the same amount of time, I decided to put more effort into debugging the Painleve case.

What I've done, essentially, is integerated out the extra dimension normally used in the geodesic equations.

Carl

 

[MeJennifer]

Exiting! When can we see the working applet?

 

[CarlB]

Exiting! When can we see the working applet?

Here it is:
http://www.gaugegravity.com/testapplet/SweetGravity.html

I'd have put it up yesterday, but it still had the debug prints turned on and I didn't have protection against division by zero at the singularity.

The Painleve orbits are labeled as "gauge" because the purpose of the applet is to show off the gauge gravity equations. The Schwarzschild coordinates are labeled "Einstein".

The Schwarzschild and Painleve orbits are similar but not quite identical. I can make them identical by correcting the initial conditions. That is the initial conditions for (coordinate) velocity:

$$\left. \frac{dx}{dt}\right|_0$$

means different things for Schwarzschild and Painleve coordinates because t is different between them. I'm contemplating adding a button to the applet that will change the initial conditions for Schwarzschild to match the Painleve or vice versa.

By the way, if you look around in the literature, you will find formulas for approximate relativistic corrections to Newton's equations of motion. What I've computed here are the exact relativistic corrections to Newton's equations of motion. After I do the rotating black hole in Doran coordinates I will publish this.

Here are some papers talking about 1st order corrections to Newton:

http://www.numdam.org/numdam-bin/fitem?id=AIHPA_1985__43_1_107_0
http://www.obs-azur.fr/gemini/pagesperso/pireaux/proc/SCRMI_integrator.pdf
http://syrte.obspm.fr/journees2004/PDF/Pireaux.pdf
http://adsabs.harvard.edu/abs/1986IAUS..114..105N

[Latest Update]: The equations take the form of a ratio of two quantities. In post #53 I factored the denominator factors into:

$$4 \;I \;\left(4I - \left(\frac{\partial I}{\partial \dot{x}}\frac{\partial I}{\partial \dot{y}}\right)^2\right)$$

As usual, this is wrong. The correct factorization is:

$$4 \;I \;\left(4I - \left(\frac{\partial I}{\partial \dot{x}}\right)^2 - \left(\frac{\partial I}{\partial \dot{y}}\right)^2\right)$$

It turns out that the complicated part of the factorization simplifies to -4. Thus the denominator is:

$$-16 I$$

I've now factored the $$I$$ out of the numerator. Consequently I have a set of equations that have no singularities other than at the origin. The numerator is still messy, but probably can be simplified (which I'm doing).
[/Latest Update]

Carl

 

[CarlB]

From reducing the equations of motion, it's pretty clear that your life is easier if you write your Painleve metric in the following fashion:

$$ds^2 = \left(dx + \sqrt{2}xr^{-1.5}\;dt\right)^2 + \left(dy + \sqrt{2}yr^{-1.5}\;dt\right)^2 + \left(dz + \sqrt{2}zr^{-1.5}\;dt\right)^2 - dt^2$$

Then the calculations end up using terms like $$\dot{x}+\sqrt{2}xr^{-1.5}$$, which in the river vernacular is the velocity relative to the river. For example, if $$x$$ and $$\dot{x}$$ are both positive, then you are climbing out of the black hole and swimming against the stream. Therefore the relative velocity is increased.

I've got the terms fairly well reduced and well behaved, but I suspect that if I play with them for a little longer I'll get them much better.

Carl

 

[ZapperZ]

Based on all the reports that I have received, this thread has been going in the wrong path for a long time.

If you wish to continue, please do so in the IR forum.

Zz.