BRS: Using Maxima for gtr computations

[Chris Hillman]

Many of you are probably already familiar with the following phenomena:

• Many gtr-related questions in PF are often best answered using the result of some computation; for example, someone often incorrectly believes that the Einstein tensor of his manifold vanishes),
• Computations "by hand" in gtr are error prone and time consuming,
• SA/Ms typically don't have a lot of free time.

What to do?

Use a computer!

I think it is fairly well known among SA/Ms that a free package which runs under Maple (which is definitely not free), GRTensorII, is extremely powerful and very well suited for computations with explicit spacetimes. (Since we can define these using undetermined metric functions, explicit spacetimes is not the same thing as a specific spacetime!)

In another BRS thread, I plan to give some hints about using GRTensorII with Maple to do all kinds of powerful stuff. Here I want to explain how to use a free package, Maxima, an open source cousin of Maple which is currently, sad to say, much more primitive than Maple, but which does contain packages which you can use for some simple index gymnastics (itensor) and for componentwise computations for explicit spacetimes (ctensor). Best of all, ctensor supports frame field computations, although there are few trouble spots.

I should caution that the Maple/GRTensorII combo is much more powerful than Maxima/ctensor: GRTensorII offers

• much better input/output, e.g. easy raising and lowering of indices,
• useful derivative operators, e.g. covariant derivative, Lie derivative,
• automatic computation of any other tensorial quantities needed to compute the requested components,
• more flexible spacetime definitions (can use frame or coframe, can use constraint equations),
• the ability to define and compute the components of (almost) arbitrary tensorial objects,
• the ability to repeatedly apply constraint equations (very useful in working with families like the Weyl family of all static axisymmetric vacuum solutions!),
• much better defined quantities, e.g. Ricci spinors in addition to Weyl spinors,
• some sophisticated trickery which improves performance in gtr computations,
• while data passing between GRTensorII and Maple is not always as simple as one might desire, one can easily take advantage of Maple's powerful symbolic DE solvers, symbolic integrators, and the powerful casesplit command.

But for many SA/Ms, Maxima/ctensor offers an advantage which trumps all these: it's available free of charge.

Actually, there are at least two aspects in which I prefer Maxima to Maple:

• .mpl files are not as convenient to comment as .mac files,
• Maple is not always very good at letting the user declare dependent and independent variables, so symbolic expressions involving partial derivatives in Maple often appear unneccessarily fussy, at times even unreadable, although this can be suppressed in the output, this is not so easy to suppress in the input.)

Before I forget, let me say that I would be happy to hear from any programmers amongst you who might be interested in improving the capabilities of ctensor or even Maxima itself, although I have as yet no connection to these open source projects.

Maxima (which comes with itensor and ctensor as standard packages) is available for Windows, Mac, and *nix users. I myself use Debian Linux, and I can't help anyone with installation issues in Windows or Mac, I can probably help with Linux installation, so shoot me a PM if you try to intall it and have problems. If you use Debian Lenny (the current "stable" branch of Debian), you can easily install Maxima using your package manager (which will pull in any further required packages). I recommend also installing wxmaxima, a GUI which seems to work better than the alternative GUI on offer in the Debian repositories, but caution that wxmaxima has suffered from a security hole (which has been patched in the current version of the package). (To be fair, Maple uses Java, which is a kinda scary thought.)

I think the simplest way to get started computing some simple quantities useful for gtr, once you have installed Maxima and wxmaxima, is to use some sample batch files which I will provide. If you run these under wxmaxima in batch mode, Maxima will simply load ctensor and execute some commands in sequence. Then, roughly speaking, simply copying a template and editing it to define the frame you want is all that is required to use Maxima to compute components using the frame you have defined.

 

[Chris Hillman]

BRS: Using Maxima for GTR: First Example

Let's start by deriving the Schwarzschild vacuum. Directions for the demo

• cut and paste the file below into your favorite text editor, to make a file called something like Schwarz_vac_deriv.mac
• put this the directory where you want to store your Maxima files,
• start wxmaxima
• in the wxmaxima menu, choose File -> Batch,
• navigate to the appropriate directory
• open Schwarz_vac_deriv.mac

Maxima will load ctensor and execute some commands. The action is rather fast, but you can scroll back to review.

/* 
Ansatz for Schwarzschild vac; Painleve chart; Lemaitre coframe

*/
load(ctensor);
cframe_flag: true;
ratchristof: true;
ctrgsimp: true;
/* define the dimension */
dim: 4;
/* list the coordinates */
ct_coords: [t,r,theta,phi];
/* define background metric */
lfg: ident(4);
lfg[1,1]: -1;
/* Declare the dependent and independent variables */
depends(f,r);
/* Define the coframe covectors */
/* only need enter the nonzero components */
fri: zeromatrix(4,4);
fri[1,1]: -1;
fri[2,1]:  f;
fri[2,2]:  1;
fri[3,3]:  r;
fri[4,4]:  r*sin(theta);
/* setup the spacetime definition */
cmetric();
/* compute a matrix whose rows give frame vectors */
fr;
/* metric tensor g_(ab) */
lg;
/* compute g^(ab) */
ug: trigsimp(invert(lg));
christof(false);
/* Compute fully covariant Riemann components R_(mijk) = riem[i,k,j,m] */
lriemann(true);
/* Compute R^(mijk) */
uriemann(false);
/* Compute Ricci componets R_(jk) */
ricci(false);
/* Compute trace of Ricci tensor */
tracer;
/* Compute R^(jk) */
uricci(false);
/* Compute and display MIXED Einstein tensor G^a_b */
/* For (-1,1,1,1) sig Flip sign of top row to get G^(ab) */
einstein(false);
cdisplay(ein);
/* WARNING! leinstein(false) only works for metric basis! */
/* Solve the Einstein field equation */
ein[1,1];
ode2(%,f,r);
radcan(%);
constant(m); subst( sqrt(2*m/r), f, ein[1,1]); ev(%,diff); ratsimp(%);
subst( sqrt(2*m/r), f, ein[3,3]); ev(%,diff); ratsimp(%);

Here's what that file is supposed to illustrate:

Suppose we want to find a vacuum solution of the EFE such that test particles fall spherically inward from r=infinity. The simplest frame Ansatz which could possibly work is
$$\begin{array}{rcl} \vec{e}_1 & = & \partial_t - f \; \partial_r \\ \vec{e}_2 & = & \partial_r \\ \vec{e}_3 & = & \frac{1}{r} \; \partial_\theta \\ \vec{e}_4 & = & \frac{1}{r \; \sin(\theta)} \; \partial_\phi$$
(look for fr in the ctensor output, or type "fr;" in the wxmaxima command pane). Here, f is some undetermined (positive) real valued function of r. This Ansatz says that we assume the spatial vectors behave just like the spatial vectors for the obvious frame for polar trig chart in euclidean three-space, while the timelike unit vector leans toward the origin by an amount depending only on r. (If we had been less bold and had assumed instead that f depends on t as well as r, we'd wind up with the same result, although this is not obvious!)

The dual coframe (look for fri in the output) is
$$\begin{array}{rcl} \sigma^1 & = & -dt \\ \sigma^2 & = & dr + f \, dt\\ \sigma^3 & = & r \; d\theta \\ \sigma^4 & = & r \; \sin(\theta) \; d\phi$$
By definition (look for lfg), the metric tensor is
$$g = -\sigma^1 \otimes \sigma^1 + \sigma^2 \otimes \sigma^2 + \sigma^3 \otimes \sigma^3 + \sigma^4 \otimes \sigma^4$$
which expanded out (look for lg in the output) is
$$ds^2 = -(1-f^2) \; dt^2 + 2 \, f \, dt \, dr + dr^2 + r^2 \, \left( d\theta^2 + \sin(\theta)^2 \; d\phi^2 \right)$$
Computing the Einstein tensor, we find this has only two distinct nonzero components. Killing the simplest gives a solution $f = c/\sqrt{r}$ which, when substituted into the remaining components turns out to kill them too. So we found a vacuum solution. With a bit of foresight, let's put $c = \sqrt{2m}$.

Notice that unlike GRTensorII, ctensor requires us to define spacetimes in terms of coframes, rather than frames. In practice this can be slightly annoying, but it is not generally a huge drawback.

Next, plug our solution into a new file:

/* 
Schwarzschild vacuum; Painleve chart; Lemaitre nsi coframe 
*/
load(ctensor);
cframe_flag: true;
ratchristof: true;
ctrgsimp: true;
/* define the dimension */
dim: 4;
/* list the coordinates */
ct_coords: [tau,r,theta,phi];
/* define background metric */
lfg: ident(4);
lfg[1,1]: -1;
/* Declare the dependent and independent variables */
constant(m);
/* Define the coframe covectors */
/* only need enter the nonzero components */
fri: zeromatrix(4,4);
fri[1,1]: -1;
fri[2,1]:  sqrt(2*m/r);
fri[2,2]:  1;
fri[3,3]:  r;
fri[4,4]:  r*sin(theta);
/* setup the spacetime definition */
cmetric();
/* compute a matrix whose rows give frame vectors */
fr;
/* metric tensor g_(ab) */
lg;
/* compute g^(ab) */
ug: trigsimp(invert(lg));
christof(false);
/* Compute fully covariant Riemann components R_(mijk) = riem[i,k,j,m] */
lriemann(true);
/* Compute R^(mijk) */
uriemann(false);
/* Compute Ricci componets R_(jk) */
ricci(true);
/* Compute trace of Ricci tensor */
tracer;
/* Compute R^(jk) */
uricci(false);
/* Compute and display MIXED Einstein tensor G^a_b */
/* For (-1,1,1,1) sig Flip sign of top row to get G^(ab) */
einstein(false);
cdisplay(ein);
/* WARNING! leinstein(false) only works for metric basis! */
/* Compute tidal tensor */
EX[2,2] = lriem[2,2,1,1];
EX[2,3] = lriem[2,3,1,1];
EX[2,4] = lriem[2,4,1,1];
EX[3,2] = lriem[3,2,1,1];
EX[3,3] = lriem[3,3,1,1];
EX[3,4] = lriem[3,4,1,1];
EX[4,2] = lriem[4,2,1,1];
EX[4,3] = lriem[4,3,1,1];
EX[4,4] = lriem[4,4,1,1];
/* Construct NP tetrad for our frame, compute Weyl spinors and Petrov type */
nptetrad(true);
weyl(false);
psi(true);
factor(realpart(psi[2]));
factor(imagpart(psi[2]));
petrov();
/* Compute Kretschmann scalar */
rinvariant();
expand(factor(%));

Running this, you find that we computed the tidal tensor, verifying that we made the correct choice of integration constant, and we also computed an NP tetrad, the Weyl spinors (the only one which is nonzero turns out to be $\Psi_2 = m/r^3$, and verified that we found a "degenerate vacuum solution" which is Petrov type D. And just for fun, we computed the Kretschmann scalar
$$R_{abcd} \, R^{abcd} = \frac{48 \, m^2}{r^6}$$

If this doesn't look like the Schwarzschild vacuum to you, that's probably because not everyone is familiar with the Painleve chart (1921) which the earliest and still the best of all charts adapted for studying radial infall. What we did was to derive the Lemaitre frame
$$\begin{array}{rcl} \vec{e}_1 & = & \partial_t - \sqrt{2m/r} \; \partial_r \\ \vec{e}_2 & = & \partial_r \\ \vec{e}_3 & = & \frac{1}{r} \; \partial_\theta \\ \vec{e}_4 & = & \frac{1}{r \; \sin(\theta)} \; \partial_\phi$$
which also defines a dual coframe and thus defines the metric tensor. So this illustrates a very physical approach to solving the EFE: we set up a thought experiment involving some family of test particles whose motion will implicitly define the spacetime, where we assume for convenience that this motion exhibits various symmetries, and finally, we can use a Newtonian limit or other physical reasoning (perhaps evaluating a Komar mass) to deduce appropriate values for any parameters which arose as integration constants when we solved an ODE. As illustrated in the BRS:Weyl vacuums thread, we won't always obtain an ODE; in general we obtain rather nasty system of PDEs which can usually be greatly simplified (triangularized) using powerful Groebner-like computational techniques from the theory of differential rings. Maple can do that stuff; Maxima cannot, not yet.

Still not convinced? It's a vacuum solution, and its spherically symmetric, indeed static, so it must be the Schwarzschild vacuum. And it's not very hard to find a coordinate transformation from Painleve to Schwarzschild (or Droste-Weyl) coordinates.

 

[Chris Hillman]

BRS: Using Maxima for gtr: Second example

Here is a more impressive example

/* 
Kerr vacuum; Boyer-Lindquist chart; B-L coframe

This does give the correct Weyl scalars wrt standard B-L NP tetrad!
*/
load(ctensor);
cframe_flag: true;
ratchristof: true;
ctrgsimp: true;
/* define the dimension */
dim: 4;
/* list the coordinates */
ct_coords: [t,r,theta,phi];
/* define background metric */
lfg: ident(4);
lfg[1,1]: -1;
/* define the frame */
fri: zeromatrix(4,4);
fri[1,1]: -sqrt(r^2-2*m*r+a^2)/sqrt(r^2+a^2*cos(theta)^2);
fri[1,4]:  sqrt(r^2-2*m*r+a^2)/sqrt(r^2+a^2*cos(theta)^2)*a*sin(theta)^2;
fri[2,2]:  sqrt(r^2+a^2*cos(theta)^2)/sqrt(r^2-2*m*r+a^2);
fri[3,3]:  sqrt(r^2+a^2*cos(theta)^2);
fri[4,1]: -a*sin(theta)/sqrt(r^2+a^2*cos(theta)^2);
fri[4,4]:  (r^2+a^2)*sin(theta)/sqrt(r^2+a^2*cos(theta)^2);
/* setup the spacetime definition */
cmetric();
/* display matrix whose rows give coframe covectors */
fri;
/* compute a matrix whose rows give frame vectors */
fr;
/* metric tensor g_(ab) */
lg;
/* compute g^(ab) */
ug: invert(lg);
christof(false);
/* Compute fully covariant Riemann components R_(mijk) = lriem[i,k,j,m] */
lriemann(true);
/* Compute R^(mijk) */
uriemann(false);
/* Compute Ricci componets R_(jk) */
ricci(true);
/* Compute trace of Ricci tensor */
tracer;
/* Compute R^(jk) */
uricci(false);
/* Compute and display MIXED Einstein tensor G^a_b */
/* For (-1,1,1,1) sig Flip sign of top row to get G^(ab) */
einstein(false);
cdisplay(ein);
/* WARNING! leinstein(false) only works for coordinate basis! */
/* Compute tidal tensor */
EX[2,2] = lriem[2,2,1,1];
EX[2,3] = lriem[2,3,1,1];
EX[2,4] = lriem[2,4,1,1];
EX[3,2] = lriem[3,2,1,1];
EX[3,3] = lriem[3,3,1,1];
EX[3,4] = lriem[3,4,1,1];
EX[4,2] = lriem[4,2,1,1];
EX[4,3] = lriem[4,3,1,1];
EX[4,4] = lriem[4,4,1,1];
/* Compute magneto-Riemann tensor */
BX[2,2] = lriem[2,4,3,1];
BX[2,3] = lriem[2,2,4,1];
BX[2,4] = lriem[2,3,2,1];
BX[3,2] = lriem[3,4,3,1];
BX[3,3] = lriem[3,2,4,1];
BX[3,4] = lriem[3,3,2,1];
BX[4,2] = lriem[4,4,3,1];
BX[4,3] = lriem[4,2,4,1];
BX[4,4] = lriem[4,3,2,1];
/* Construct NP tetrad for our frame, compute Weyl spinors and Petrov type */
nptetrad(true);
weyl(false);
psi(true);
factor(realpart(psi[2]));
factor(imagpart(psi[2]));
petrov();

This file computes

• the electro-Riemann or tidal tensor and the magneto-Riemann tensor for the Boyer-Lindquist observers, defined by the Boyer-Lindquist coframe in Boyer-Lindquist chart,
• the NP tetrad corresponding to the BL coframe
• the Weyl spinors
• the Petrov type

Only one spinor is nonvanishing; its real part corresponds to the tidal tensor, while its imaginary part corresponds to the magneto-Riemann tensor (or "magnetic part of the Riemann tensor"), which decays faster as r grows.
$$\begin{array}{rcl} \Re \Psi_2 & = & \frac{m\,r\,\left( 3\,{a}^{2}\,{\cos\left( \theta\right) }^{2}-{r}^{2}\right) }{{\left( {a}^{2}\,{\cos\left( \theta\right) }^{2}+{r}^{2}\right) }^{3}} \\ \Im \Psi_2 & = & -\frac{a\,m\,\cos\left( \theta\right) \,\left( {a}^{2}\,{\cos\left( \theta\right) }^{2}-3\,{r}^{2}\right) }{{\left( {a}^{2}\,{\cos\left( \theta\right) }^{2}+{r}^{2}\right) }^{3}}$$
(I used the "copy Tex" feature of wxmaxima to copy the output from maxima in latex form into my browser, and then introduced a few minor changes to make the latex markup which the PF server used to run tex to produce the equations just displayed.) The magneto-Riemann tensor controls things like spin-spin forces and is generally much smaller than the tidal tensor (the major exception: pure gravitational radiation).

The Weyl spinor $\Psi_2$ is associated with a Coulomb or nonradiative contribution to the curvature. If we can find an NP tetrad (essentially the same thing as a frame, just expressed a bit differently) in which this is real, we have a family of observers who measure no spin-spin forces, which is indicative of a nonrotating source. It turns out that the Kerr vacuum does not admit such observers; the source is rotating, as shown the scalar $R_{abcd} {^{\ast} R}^{abcd}$ which ctensor unfortunately does not compute (yet). See Wheeler and Ciufolini, Gravitation and Inertia.

The fact that only $\Psi_2$ is nonvanishing when we compute the Weyl spinors using the BL tetrad (these quantities are tetrad dependent) shows that the BL tetrad is adapted to the symmetries of the Kerr vacuum. Here is another frame on the Kerr vacuum, the Doran frame, which is not so well adapted to the Weyl tensor, but which has the virtue that it defines inertial observers who are analogous to the Painleve observers in the nonrotating case (Schwarzschild vacuum):

  /* 
Kerr vacuum; Doran influx chart; Doran coframe

The Doran frame is inertial but spinning about e_4
The timelike geodesic congruence obtained from e_1 is vorticity-free
(hypersurface-orthogonal), so can consider riemann tensor of hyperslices
Differences in Doran time coordinate between two events on world line of
a Doran observer corresponds to lapse proper time measured by that obsvr.

Electroriemann tensor is:
factor(trigsimp(lriem[2,2,1,1]));
factor(trigsimp(lriem[2,3,1,1]));
factor(trigsimp(lriem[3,3,1,1]));
factor(trigsimp(lriem[4,4,1,1]));

WARMING! Maxima is unable to reduce certain expressions to zero, 
so ctensor thinks this is Petrov I, but Kerr vacuum is actually Petrov D!
*/
load(ctensor);
cframe_flag: true;
ratchristof: true;
ctrgsimp: true;
/* define the dimension */
dim: 4;
/* list the coordinates */
ct_coords: [tau,r,theta,phi];
/* define background metric */
lfg: ident(4);
lfg[1,1]: -1;
/* define the coframe */
fri: zeromatrix(4,4);
fri[1,1]: -1;
fri[2,1]:  sqrt(2*m*r)/sqrt(r^2+a^2*cos(theta)^2);
fri[2,2]:  sqrt(r^2+a^2*cos(theta)^2)/sqrt(r^2+a^2);
fri[2,4]:  a*sin(theta)^2*sqrt(2*m*r)/sqrt(r^2+a^2*cos(theta)^2);
fri[3,3]:  sqrt(r^2+a^2*cos(theta)^2);
fri[4,4]:  sqrt(r^2+a^2)*sin(theta);
/* setup the spacetime definition */
cmetric();
/* display matrix whose rows give coframe covectors */
fri;
/* compute a matrix whose rows give frame vectors */
fr;
/* metric tensor g_(ab) */
lg;
/* compute g^(ab) */
ug: invert(lg);
christof(false);
/* Compute fully covariant Riemann components R_(mijk) = lriem[i,k,j,m] */
lriemann(true);
/* Compute R^(mijk) */
uriemann(false);
/* Compute Ricci componets R_(jk) */
ricci(true);
/* Compute trace of Ricci tensor */
tracer;
/* Compute R^(jk) */
uricci(false);
/* Compute and display MIXED Einstein tensor G^a_b */
/* For (-1,1,1,1) sig Flip sign of top row to get G^(ab) */
einstein(false);
cdisplay(ein);
/* WARNING! leinstein(false) only works for coordinate basis! */
/* Compute tidal tensor */
EX[2,2] = lriem[2,2,1,1];
EX[2,3] = lriem[2,3,1,1];
EX[2,4] = lriem[2,4,1,1];
EX[3,2] = lriem[3,2,1,1];
EX[3,3] = lriem[3,3,1,1];
EX[3,4] = lriem[3,4,1,1];
EX[4,2] = lriem[4,2,1,1];
EX[4,3] = lriem[4,3,1,1];
EX[4,4] = lriem[4,4,1,1];
/* Compute magneto-Riemann tensor */
BX[2,2] = lriem[2,4,3,1];
BX[2,3] = lriem[2,2,4,1];
BX[2,4] = lriem[2,3,2,1];
BX[3,2] = lriem[3,4,3,1];
BX[3,3] = lriem[3,2,4,1];
BX[3,4] = lriem[3,3,2,1];
BX[4,2] = lriem[4,4,3,1];
BX[4,3] = lriem[4,2,4,1];
BX[4,4] = lriem[4,3,2,1];
/* Compute and display NP tetrad and Weyl scalars */
nptetrad(true);
weyl(false);
psi(false);
factor(realpart(psi[0]));
factor(imagpart(psi[0]));
factor(realpart(psi[1]));
factor(imagpart(psi[1]));
factor(realpart(psi[2]));
factor(imagpart(psi[2]));
factor(realpart(psi[3]));
factor(imagpart(psi[3]));
factor(realpart(psi[4]));
factor(imagpart(psi[4]));
petrov();

As you can see, in the new tetrad, the Weyl spinors look quite different. We are using a new chart, to be sure, so the coordinate functions are a bit different, but even if we transformed the coordinate expressions for the Doran frame vector fields from Doran to BL coordinates, the Weyl spinors computed using the Doran chart would still look different from those computed using the BL frame. However, the Petrov type is a true invariant, if computed correctly. But sad to say, ctensor fails to simply certain quantities far enough to recognize that the Petrov type is D, and incorrectly reports it as algebraically general (type I)! GRTensorII will do better, although in more complicated cases it can make the same kind of mistake, which is why GRTensorII wisely reports "Petrov D, or simpler"!

However, probably the single biggest lack in the current version of Ctensor is that it does not compute as a defined quantity the kinematic decomposition of a timelike congruence (acceleration vector, expansion tensor, vorticity tensor). The Maxima help pages claims there is a (very crude) way to coax Maxima into computing any tensorial quantity you can define using itensor, but I have not been able to make that work very well yet.

 

[Chris Hillman]

BRS: Using Maxima for gtr computations: Editing a .mac template

Here are some quick hints for what you can change and what you cannot change, without getting into trouble:

• change the header comments appropriately
• don't touch the lines
  
  load(ctensor);
  cframe_flag: true;
  which load ctensor and tell it we will be using a coframe and its dual frame
• you can modify these flags which set up certain default simplifications (see the wxmaxima help for some other possibilities you may wish to explore):
  
  ratchristof: true;
  ctrgsimp: true;
• you probably won't want to touch the line
  
  dim: 4;
  unless you want to change the dimension of the manifold!
• modify the list of coordinate names as appropriate:
  
  ct_coords: [t,r,theta,phi];
  (note that order matters in this list!),
• you probably won't want to touch the lines
  
  lfg: ident(4);
  lfg[1,1]: -1;
  which ensures that the frame will be understood as an orthonormal basis, first timelike, next three timelike,
• modify appropriately the lines
  
  fri: zeromatrix(4,4);
  fri[1,1]: -sqrt(r^2-2*m*r+a^2)/sqrt(r^2+a^2*cos(theta)^2);
  fri[1,4]:  sqrt(r^2-2*m*r+a^2)/sqrt(r^2+a^2*cos(theta)^2)*a*sin(theta)^2;
  fri[2,2]:  sqrt(r^2+a^2*cos(theta)^2)/sqrt(r^2-2*m*r+a^2);
  fri[3,3]:  sqrt(r^2+a^2*cos(theta)^2);
  fri[4,1]: -a*sin(theta)/sqrt(r^2+a^2*cos(theta)^2);
  fri[4,4]:  (r^2+a^2)*sin(theta)/sqrt(r^2+a^2*cos(theta)^2);
  which defines the coframe by expanding each covector wrt the coordinate basis (it is convenient to define all the components to be zero and then just enter the ones which are nonzero),
• don't touch the lines
  
  /* setup the spacetime definition */
  cmetric();
  /* display matrix whose rows give coframe covectors */
  fri;
  /* compute a matrix whose rows give frame vectors */
  fr;
  /* metric tensor g_(ab) */
  lg;
  /* compute g^(ab) */
  ug: invert(lg);
  christof(false);
  /* Compute fully covariant Riemann components R_(mijk) = lriem[i,k,j,m] */
  lriemann(true);
  /* Compute R^(mijk) */
  uriemann(false);
  /* Compute Ricci componets R_(jk) */
  ricci(true);
  /* Compute trace of Ricci tensor */
  tracer;
  /* Compute R^(jk) */
  uricci(false);
  /* Compute and display MIXED Einstein tensor G^a_b */
  einstein(false);
  cdisplay(ein);
  because ctensor must compute these quantities in a certain order in order to obtain the Einstein tensor
• you probably will want to keep the lines
  
  /* Compute tidal tensor */
  EX[2,2] = lriem[2,2,1,1];
  EX[2,3] = lriem[2,3,1,1];
  EX[2,4] = lriem[2,4,1,1];
  EX[3,2] = lriem[3,2,1,1];
  EX[3,3] = lriem[3,3,1,1];
  EX[3,4] = lriem[3,4,1,1];
  EX[4,2] = lriem[4,2,1,1];
  EX[4,3] = lriem[4,3,1,1];
  EX[4,4] = lriem[4,4,1,1];
  /* Compute magneto-Riemann tensor */
  BX[2,2] = lriem[2,4,3,1];
  BX[2,3] = lriem[2,2,4,1];
  BX[2,4] = lriem[2,3,2,1];
  BX[3,2] = lriem[3,4,3,1];
  BX[3,3] = lriem[3,2,4,1];
  BX[3,4] = lriem[3,3,2,1];
  BX[4,2] = lriem[4,4,3,1];
  BX[4,3] = lriem[4,2,4,1];
  BX[4,4] = lriem[4,3,2,1];
  which reports the "electric" and "magnetic" components of the Riemann tensor
• if you want to compute the Weyl spinors wrt the NP tetrad corresponding to your frame field, include the lines
  
  nptetrad(true);
  weyl(false);
  psi(true);

 

[Chris Hillman]

BRS: Using Maxima for gtr computations: some pitfalls to avoid

Ctensor does not behave well in some respects:

• when the cframe flag is set, "leinstein" does not correctly give $G_{ab}$, but fortunately, "einstein" gives ${G^a}_b$ and then we need only flip the sign of the top row to obtain $G^{ab}$ (this simple rule would not be true if we were using a coordinate basis, of course!)
• "lriem" gives $R_{abcd}$, but ctensor uses the opposite sign of the Riemann tensor, compared to the Landau Lifschitz sign conventions,
• ctensor uses a nonstandard indexing of the Riemann tensor
• you can compensate for both the sign flip and the nonstandard indexing using the rule
  $$R_{mijk} = \operatorname{lriem} \left[ i,k,j,m \right]$$
• GRTensorII observes a convention, awkward but universal in the literature, in which the +--- signature is always employed when using NP formalism; ctensor is happy to use -+++ signature (the LLSC sign convention, which is used by slightly more than half of the most popular gtr textbooks), which means that even though "npi" reports the same NP tetrad as GRTensorII's "grdisplay(nullt(up))" command, the signs of the Weyl spinors are flipped wrt GRTensorII.

 

[Chris Hillman]

Computing the electroriemann and magnetoriemann tensors

A minor improvement upon my Post #3:

Given a coframe such as
$$\begin{array}{rcl} \sigma^1 & = & -1/\sqrt{2} \, (1-x \, (v-\exp(x))/2) \, du - 1/\sqrt{2} \, dv \\ \sigma^2 & = & -1/\sqrt{2} \, (1+x \, (v-\exp(x))/2) \, du + 1/\sqrt{2} \, dv \\ \sigma^3 & = & \exp(x/2) \, dx \\ \sigma^4 & = & \exp(x/2-u) \, dy \end{array}$$
it is useful to compute and display second rank tensors, particularly ones living in three dimensional spatial hyperplane elements, as matrices.

In particular, consider the electroriemann and magnetoriemann tensors defined by the timelike congruence (not neccessarily geodesic or irrotational) defined by the timelike basis vector field $\vec{U} = \vec{e}_1$ in the dual frame field--- recall that this vector field is a particular timelike unit vector field, so it defines a congruence of timelike curves in our spacetime. These are simply a rearrangement of certain frame components of the Riemann tensor:
$${E\left[\vec{e}_1\right]}_{ab} = R_{ambn} \, U^m \, U^n, \; \; {B\left[\vec{e}_1\right]}_{ab} = {R^\ast}_{ambn} \, U^m \, U^n$$
(right Hodge dual), defined analogously to the way the EM field tensor (second rank antisymmetric) is decomposed (wrt a timelike congruence) into electric and magnetic components. (We could also use left Hodge dual, which annoyingly gives a slightly different result except in the vacuum case, when the Riemann tensor reduces to the Weyl tensor.) Here, the right Hodge dual of the Riemann tensor is
$${R^\ast}_{abcd} = {R_{ab}}^{mn} \, \varepsilon_{mncd}$$
(Levi-Civita tensor, the same one we use to write the volume form). If we write the Riemann tensor in terms of curvature two-forms, this really is what we obtain by taking the Hodge duals of those two-forms, using the formalism of exterior calculus (see Flanders), and reassembling them into a fourth rank tensor with the symmetries of the Riemann tensor. See also MTW for several other ways of understanding the electric and magnetic components of the Riemann tensor.

(Strictly speaking, if we are using a timelike unit vector other than $\vec{e}_1$, we should project the rank two tensors to ensure they live in hyperplane elements orthogonal to $\vec{U}$, but this extra effort is not needed if we are computing the Bel decomposition associated with a coframe in terms of the same coframe.)

We can display the electroriemann tensor and magnetoriemann tensors as matrices, like this:

/* electroriemann tensor */
matrix([lriem[2,2,1,1], lriem[2,3,1,1],lriem[2,4,1,1]],
[lriem[3,2,1,1],lriem[3,3,1,1],lriem[3,4,1,1]],
[lriem[4,2,1,1],lriem[4,3,1,1],lriem[4,4,1,1]]);
/* magnetoriemann tensor */
matrix([lriem[2,4,3,1],lriem[2,2,4,1],lriem[2,3,2,1]],
[lriem[3,4,3,1],lriem[3,2,4,1],lriem[3,3,2,1]],
[lriem[4,4,3,1],lriem[4,2,4,1],lriem[4,3,2,1]]);

Here is a batch file you can play with, using the same example I started with (the coframe defining Petrov's example of a Petrov III Kundt wave, a generalization of the pp wave spacetimes to include the longitudinal shearing terms which are characteristic of Petrov III radiation, a phenomenon unknown in linearized gravitational wave theory):

/* 
Kundt vacuum plane fronted wave; Petrov's example;  Petrov chart 

Chart defined on
        -infty < u, v,x,y < infty

This is a Kundt plane fronted vacuum wave with Weyl tensor type Petrov III
The wave vector is null geodesic congruence @_v
which has vanishing optical scalars, but is not covariantly constant

Curvature clearly displays Petrov III nature:
Electroriemann tensor (aka tidal tensor) shows 
both 
	plus polarized transverse/spintwo tidal acc Petrov N component
	(polarized) longitudinal/shearing tidal acc Petrov III component
of the wave:
	E_(23) = R_(1213) = lriem[2,3,1,1] = -exp(-x/2)/2/sqrt(2)
	E_(33) = R_(1313) = lriem[3,3,1,1] =  (3+x)/8+v*exp(-x)/8
	E_(44) = R_(1414) = lriem[4,4,1,1] = -(3+x)/8-v*exp(-x)/8
These correspond to Psi4 and Psi3 respectively (both real).
*/
load(ctensor);
cframe_flag: true;
ctrgsimp: true;
/* define the dimension */
dim: 4;
/* list the coordinates */
ct_coords: [u,v,x,y];
/* define background metric */
lfg: ident(4);
lfg[1,1]: -1;
/* define the dependent and independent variables */
depends([a,b,h,k],[u,x,y]); 
/* define the coframe; only need enter the nonzero components */
fri: zeromatrix(4,4);
fri[1,1]: -1/sqrt(2)*(1-x*(v-exp(x))/2);
fri[1,2]: -1/sqrt(2);
fri[2,1]: -1/sqrt(2)*(1+x*(v-exp(x))/2);
fri[2,2]:  1/sqrt(2);
fri[3,3]:  exp(x/2);
fri[4,4]:  exp(x/2-u);
/* setup the spacetime definition */
cmetric();
/* compute a matrix whose rows give frame vectors */
fr;
/* metric tensor g_(ab) */
lg;
/* compute g^(ab) */
ug: invert(lg);
/* compute geodesic equations */
christof(false);
cgeodesic(true);
/* Compute fully covariant Riemann components R_(mijk) = riem[i,k,j,m] */
lriemann(true);
/* Compute R^(mijk) */
uriemann(false);
/* Compute Ricci componets R_(jk) */
ricci(true);
/* Compute trace of Ricci tensor */
tracer;
/* Compute R^(jk) */
uricci(false);
/* Compute and display MIXED Einstein tensor G^a_b */
/* For (-1,1,1,1) sig Flip sign of top row to get G^(ab) */
einstein(false);
cdisplay(ein);
/* WARNING! leinstein(false) only works for metric basis! */
/* electroriemann tensor */
print("electroriemann tensor");
matrix([lriem[2,2,1,1], lriem[2,3,1,1],lriem[2,4,1,1]],
[lriem[3,2,1,1],lriem[3,3,1,1],lriem[3,4,1,1]],
[lriem[4,2,1,1],lriem[4,3,1,1],lriem[4,4,1,1]]);
/* magnetoriemann tensor */
print("magnetoriemann tensor");
matrix([lriem[2,4,3,1],lriem[2,2,4,1],lriem[2,3,2,1]],
[lriem[3,4,3,1],lriem[3,2,4,1],lriem[3,3,2,1]],
[lriem[4,4,3,1],lriem[4,2,4,1],lriem[4,3,2,1]]);
print("spin two transverse (Type N) contrib to tidal tensor components look like this:");
expand(lriem[3,3,1,1]);
wxplot3d(%, [x,-5,5], [v,-5,5]);
print("shearing longitudinal (Type III) contrib to tidal tensor components look like this:");
lriem[2,3,1,1];
wxplot2d([%], [x,-5,5]);
/* Compute Kretschmann scalar */
rinvariant();
expand(factor(%));
/* Compute NP tetrad, Weyl spinors, and Petrov type */
weyl(false);
nptetrad(true);
psi(true);
petrov();

The third piece of the Bel decomposition, the toporiemann tensor, is obtained from the left dual of the right dual, aka the double dual, of the Riemann tensor
$${L\left[\vec{e}_1\right]}_{ab} = {{}^\ast \! R^\ast}_{ambn} \, U^m \, U^n$$
In general, E and L are symmetric (6 algebraically independent components) while B is traceless (8 algebraically independent components), making up 20 algebraically independent components in all. As this suggests, the Riemann tensor is completely characterized by its Bel decomposition. As happens with the decomposition of the EM field tensor into electric and magnetic pieces, the Bel decomposition depends on what congruence you are using, but offers great physical insight.

In a vacuum spacetime, B, E,L are traceless symmetric and L is algebraically equivalent to E, so we have only 10 algebraic independent components in all. This is of course the Bel decomposition of the Weyl tensor, which is completely characterized by its Bel decomposition.

In higher dimensions, a fourth tensor is needed (and unfortunately Bel's first initial is L). Therefore higher dimensions are deprecated on notational grounds

Another example:

/* 
Reissner-Nordstrom nnevac; exterior psph; static coframe 

This models exterior of spherically symmetric electrically charged object
Given coframe corresponds to static observers who use their rocket engines
to hover over this object.  The coordinate r is Schwarzschild coordinate,
defined such that area of spheres r=r0 is %pi*r0^2

From Killing vectors we obtain three first integrals for t1,theta1,phi1.
Plugging these into epsilon = ds^2 we obtain one for r1.  This gives
general geodesic up to integration over proper time.
WLOG, consider null geodesic in the equatorial plane theta=Pi/2:
	t1 = E/(1-2*m/r+q^2/r^2)
	r1 = +/-sqrt(E^2-(L/r)^2*(1-2*m/r+q^2/r^2)
	theta1 = 0
	phi1 = L/r^2
where variable is affine parameter s (unique up to constant scalar multiple).
Note that r_min is positive real root of
	E^2-(L/r)^2*(1-2*m/r+q^2/r^2
Thus the unparameterized world line in "graphing form" is
	dr/dt = sqrt(1-(L/E/r)^2*(1-2*m/r+q^2/r^2))*(1-2*m/r+q^2/r^2)
	dphi/dt = L/E/r^2*(1-2*m/r+q^2/r^2)
Note this depends only on ratio L/E, so the "world line of a wave packet"
is indpt of frequency!  Also
	dphi/dr = L/E/r^2/sqrt(1-(L/E/r)^2*(1-2*m/r+q^2/r^2))
Integrating from r=-infty to r=r_min and doubling gives exact light bending
as an elliptic integral.

*/
load(ctensor);
cframe_flag: true;
ratchristof: true;
ctrgsimp: true;
/* define the dimension */
dim: 4;
/* list the coordinates */
ct_coords: [t,r,theta,phi];
/* define background metric */
lfg: ident(4);
lfg[1,1]: -1;
/* rows of this matrix give the coframe covectors */
/* only need enter the nonzero components */
fri: zeromatrix(4,4);
fri[1,1]: -sqrt(1-2*m/r+q^2/r^2);
fri[2,2]:  1/sqrt(1-2*m/r+q^2/r^2);
fri[3,3]:  r;
fri[4,4]:  r*sin(theta);
/* setup the spacetime definition */
cmetric();
/* compute a matrix whose rows give frame vectors */
fr;
/* metric tensor g_(ab) */
lg;
/* compute g^(ab) */
ug: invert(lg);
christof(false);
/* Compute fully covariant Riemann components R_(mijk) = riem[i,k,j,m] */
/* Note: for each j,m, cdisplay(riem) displays riem[i,k,-,-] as matrix */
lriemann(true);
/* Compute R^(mijk) */
uriemann(false);
/* Compute Ricci componets R_(jk) */
ricci(true);
/* Compute trace of Ricci tensor */
tracer;
/* Compute R^(jk) */
uricci(false);
/* Compute and display MIXED Einstein tensor G^a_b */
/* For (-1,1,1,1) sig Flip sign of top row to get G^(ab) */
einstein(false);
cdisplay(ein);
/* WARNING! leinstein(false) only works for metric basis! */
/* electroriemann tensor */
print("electroriemann tensor");
expand(matrix([lriem[2,2,1,1], lriem[2,3,1,1],lriem[2,4,1,1]],
[lriem[3,2,1,1],lriem[3,3,1,1],lriem[3,4,1,1]],
[lriem[4,2,1,1],lriem[4,3,1,1],lriem[4,4,1,1]]));
/* magnetoriemann tensor */
print("magnetoriemann tensor");
matrix([lriem[2,4,3,1],lriem[2,2,4,1],lriem[2,3,2,1]],
[lriem[3,4,3,1],lriem[3,2,4,1],lriem[3,3,2,1]],
[lriem[4,4,3,1],lriem[4,2,4,1],lriem[4,3,2,1]]);
/* Compute Kretschmann scalar */
rinvariant();
expand(factor(%));
/* Compute NP tetrad, Weyl spinors, and Petrov type */
weyl(false);
nptetrad(true);
psi(true);
petrov();

It turns out that the electroweyl tensor does not capture the gravitational effect of the "immediate presence of EM field energy and momentum" at each event in the exterior of a static spherically symmetric electrically charged object. For an electrovacuum, the electroriemann tensor includes a just enough from the Einstein tensor to account for this, so it is the correct "tidal tensor" to use in a nonvacuum. To be sure, the difference is generally very small; in this example
$${E\left[\vec{e}_1\right]}_{22} = \frac{2m}{r^3} + \frac{3 \, q^2}{r^4}, \; \; {EW\left[\vec{e}_1\right]}_{22} = \frac{2m}{r^3} + \frac{2 \, q^2}{r^4}$$
where q is generally quite small, compared with m, for an astrophysical object.

At the time of writing, IMO the major obstacle to regarding Maxima and Itensor/Ctensor as tools worthy of attention from serious students is the lack of any easy way to (roughly in decreasing order of importance):

• compute the acceleration vector, expansion tensor, and vorticity tensor wrt a timelike congruence,
• compute optical scalars (expansion, shear, twist) for a null geodesic congruence (this is should "secretly" be an intermediate step in computing the Weyl spinor components in the NP formalism, however!),
• define and apply constraint equations,
• apply simplifications to an entire tensor, rather than component by component,
• solve nontrivial ODEs,
• solve PDEs by separation of variables,
• use differential rings to triangularize systems of possibly nonlinear PDEs,
• compute and solve the Killing equations (easiest route to solving the geodesic equations, in many cases),
• compute multivariable Taylor series, asymptotic expansions, componentwise,
• define and compute useful scalar invariants,
• use "exotic" metrics on tangent spaces (e.g. for Coll canonical charts).

all things GRTensorII under Maple does very well. I'd mention

• decide equivalence up to local isometry

but that wouldn't be fair, because GRTensorII doesn't (yet) do this either.

 

[Chris Hillman]

How to master the world of exact solutions of the EFE

I have finally been able to borrow a copy of the new "field guide to exact solutions of the EFE", which joins the current generation of gtr textbooks from the CUP stable:

• Jerry B. Griffiths and Jiri Podolsky,
  Exact Space-Times in Einstein's General Relativity,
  Cambridge University Press, 2009
• Jerzy Plebanski and Andrzej Krasinski,
  An Introduction to General Relativity and Cosmology,
  Cambridge University Press, 2006
• Eric Poisson,
  A Relativist's Toolkit,
  Cambridge University Press, 2004

(CUP: the last honest publisher, IMO--- may they never fail, or go over to dark side of the universal population surveillance industry, the arms trade, or blood diamonds.)

Since I love Griffith's delightful and very readable monograph on colliding plane waves, which he has generously made available for free download at his academic website (see "Links for SA/Ms"), and know both authors as the authors of research papers I like, I had high expectations, and after my first pass through this book, I am delighted to say my expectations have been exceeded!

I love books, obviously, and it seems I never met a gtr book I didn't love for one reason or another. But even though I know I can't simply recommend that every SA/M buy all those published since 1995 plus the indispensable MTW, I think the new book by Griffiths and Podolsky really is an essential desk reference which you can think of as a "field guide to the important known exact solutions of the EFE". With the not unreasonable exception, given length constraints, of static spherically symmetric perfect fluids, LTB dusts, and other important classes which are well treated in review papers or other recent textbooks (Chapter 14 of P&K is the essential reference for LTB dusts, a core topic for any serious student of cosmology).

(MTW plus D'Inverno plus the above three would make a fine bookshelf for anyone truly devoted to gtr. I love several other modern textbooks such as Carroll, which has perhaps the most modern approach of the recent "elementary gtr textbooks", but today I'll say D'Inverno because of its breadth of coverage, which should help modern students appreciate the variety of topics and techniques they can seek to master. P&K is however the first textbook which attempts to present the way capable researchers in classical gravitation actually work--- routinely employing frame fields, null tetrads, kinematic decomposition of geodesic congruences, and other geometric techniques--- so today I'll say D'Inverno followd by P&K is your best route into the current research literature. The fun part is, papers on exact solutions are easy to read once you know this previously hard to find background material!)

So what will you find in G&P? A superbly organized field guide introducing all the most important classes of solutions and clearly explaining their causal structure and the rather complicated interrelationships between them. You can expect to find a clear discussion of a good selection of the most important coordinate charts (including transformations between them) for the most important solutions, but you won't find detailed discussion of the physical experience of various classes of observers in the manner I favor, based upon using frame fields and computing geometrically/physically significant quantities such as optical expansion scalars of important null geodesic congruences, the Ricci and Weyl spinor components, etc.--- the "canonical approach" for studying and interpreting specific solutions which is well illustrated by P&K in their Chapter 18 on the LTB dusts. Since even this introductory sketch and overview would run over 500 pages had the important fluid solutions been included, I think the author's choice of goals is wise: their book offers a perfect springboard for anyone who wishes to master of the world of exact solutions, by giving you enough information to dig in using the computational engines, GRTensorII under Maple, or if that is not available, Ctensor under Maxima.

So here's my suggestion for ambitious SA/Ms: unsystematically, as time permits, create Maxima Ctensor files using appropriate coframe fields (or even better, Maple GRTensorII files) for each solution and class of solution in P&K (supplemented by the important classes of fluid solutions), and follow the model of my posts here or P&K to study the physical experience of the observers corresponding to your coframe. Keep good notes, of course!

I can't resist adding a few more comments explaining why I am so delighted by the new "field guide":

• the authors's choices about what are the essential solutions and families of solutions, how much weight to give to each, and key points should be made in each case, are virtually identical to my own (the differences mainly amount to my desire to say more or to throw fluid solutions and more dust models into the mix),
• random example one: they take care to mention the important fact that Minkowski vacuum has two representations using the Weyl canonical chart (the importance of which I hinted at in my BRS on the Weyl vacuums),
• random example two: they devote entire chapters to the Taub-NUT vacuum, the Kundt class, and some of the other essential but generally underappreciated solutions/families, and in every case I enthusiastically agree that the classes they chose to devote entire chapters too amply merit such devotion,
• I've always tried to stress the point of view generally adopted by researchers like Bonnor who know the literature really well and understand what the most important insights are, that physical interpretation has received far too little attention in the research literature, even though (IMO) this is by far the most fun portion of grappling with gtr--- my background is entirely in math, and I am naturally inclined towards fascination with techniques for solving systems of DEs, but it seem very striking that so many physics grad students don't seem to recognize that simply finding a solution is only the first step in doing useful work in this area!--- so I am delighted that G&P also stress the importance of physical interpretation,
• In particular, at the expense of neglecting the stuff a good student can do himself (use a computer tool to compute frame components and interpret results), the authors emphasize the crucial role of using Carter-Penrose diagrams to understand the global structure of suitable extensions of what we might call "partial solutions" (expressions for a metric tensor and any physical fields, written in terms of some coordinate chart which is extensible).

The next BRS I had planned was in fact to survey Carter-Penrose diagrams, and if I can muster the energy I still plan to write that, as a kind of invitation to study G&P.

Oh, one other thing: I wrote above

GRTensorII observes a convention, awkward but universal in the literature, in which the +--- signature is always employed when using NP formalism; ctensor is happy to use -+++ signature (the LLSC sign convention, which is used by slightly more than half of the most popular gtr textbooks), which means that even though "npi" reports the same NP tetrad as GRTensorII's "grdisplay(nullt(up))" command, the signs of the Weyl spinors are flipped wrt GRTensorII.

G&K break this schizophrenic approach, employing -+++ signature consistently for both NP formalism and tensor gymnastics. This is the best choice for gtr, but few previous authors had the courage to adopt a -+++ signature for NP formalism, which requires careful changes to some formulas from NP formalism, as applied in the vast majority of the research literature which appeals to null tetrad techniques. So when computing Ricci and Weyl spinors from NP tetrads using GRTensorII under Maple, or Ctensor under Maxima, bear in mind that GRTensorII uses +--- while Ctensor uses -+++.