RIEMANNIAN GEOMETRY & TENSOR CALCULUS
A set of functions for tensor calculations in Riemannian Geometry and General Relativity

version 3.8.2 - June 2011

S. Bonanos

Institute of Nuclear Physics,  NCSR "Demokritos", 
Aghia Paraskevi 15310, Greece. 

Please send comments, suggestions, bug reports to:  sbonano@inp.demokritos.gr

Δωρϵ(Overscript[α, `])ν ϵλ(Overscript[α, '])βϵτϵ, δωρϵ(Overscript[α, `])ν δ(Overscript[o, '])τϵ.

1  Instructions

1.0  Introduction

This package allows calculations in Riemannian Geometry by defining all basic tensors and tensor operations. Tensors are represented as nested lists of their components in a particular frame (default = coordinate frame). Functions are defined for Raising/Lowering indices, Contracting (multiple) indices on one (or between two) tensors, for Covariant and Lie Differentiation and for classifying the Weyl and Ricci tensors. One can also use standard Mathematica functions like Dot or Outer (giving tensor products), and Transpose (index rearrangement) to create new (lists representing) tensors.

Given a metric and a set of coordinates, the main routine in the package computes explicit expressions for all common Riemannian Geometry tensors (Riemann, Ricci, Einstein, Weyl) and tests if the space belongs to any of the following categories: Flat, Conformally Flat, Ricci Flat, Einstein Space or Space of Constant Curvature. Each tensor is stored as a nested list under an appropriate global name. One can optionally specify a frame (if different from the coordinate frame) with respect to which tensor components are to be computed.

The Riemann and Ricci tensor sign conventions used are:  (Subscript[V, b;c;d]) - (Subscript[V, b;d;c]) = (Subscript[V, a]) (Subscript[(R^a), bcd]) and (Subscript[R, ab]) = (Subscript[(R^s), asb]), so that Einstein's equations are (Subscript[R, ab])-(1/2)R (Subscript[g, ab])=+κ (Subscript[T, ab]) (κ=8 π G/(c^4)).

1.1  Basic Input /Output & Notation

To define the geometry, begin by computing all basic Riemannian Geometry tensors for a particular metric. To do this:
	(0) In a new notebook, evaluate the initialization code: <<EDCRGTCcode.m
	 (1) Define a list symbols (coordinate names) and give it a name -- say xIN={t,r,θ,ϕ};
	 (2) Define a symmetric matrix (list of lists) of functions of the coordinates (the metric) and give it a name -- say 
gIN={{1-2M/r,0,0,0},{0,-1/(1-2M/r),0,0},{0,0,-r^2,0},{0,0,0,-(r^2 Sin[θ]^2)}};
	(3) Assign appropriate simplification rules to the global list simpRules as needed (see Examples and Usage Tips). For example: simpRules={Cos[θ]^2+Sin[θ]^2->1};
	(4) Evaluate RGtensors[gIN,xIN]

RGtensors  computes all common Riemannian Geometry tensors, printing messages along the way. It also tests if the space is Flat, Conformally Flat, Ricci Flat, Einstein Space or Space of Constant Curvature and, if yes, prints corresponding messages. When it is finished, the following "tensors" - nested lists - are stored in memory and can be used in further calculations (the letters U, d  in the tensor names serve as mnemonic devices for the position - Up or down - of the relevant indices) :

gdd (covariant form of the metric tensor), with components gdd[[a,b]] = (Subscript[g, ab])

gUU (contravariant form of the metric tensor), with components gUU[[a,b]] = (g^ab)

GUdd (Christoffel symbols of the 2nd kind), with components GUdd[[a,b,c]] = (Subscript[Γ, bc]^a)

Rdddd (covariant form of the Riemann tensor), with components Rdddd[[a,b,c,d]] = (Subscript[R, abcd])

RUddd (mixed form of the Riemann tensor), with components RUddd[[a,b,c,d]] = (Subscript[(R^a), bcd])

Rdd (covariant form of the Ricci tensor), with components Rdd[[a,b]] = (Subscript[R, ab]) ≡ (Subscript[(R^s), asb])

EUd (mixed form of the Einstein tensor), with components EUd[[a,b]] = (Subscript[(E^a), b]) ≡ (Subscript[(R^a), b])- R/2 (Subscript[(δ^a), b])

R is the curvature scalar    (Important: Do not use the symbol R in coordinates or metric!)

Wdddd (covariant form of the Weyl tensor), with components Wdddd[[a,b,c,d]] = (Subscript[W, abcd])

Finally, RGtensors defines the following global variables: (i) coordList contains the coordinates defined in (1) above, (ii) detg is the determinant of the metric, (iii) Dim stores the dimension of the space and (iv) LineElement stores, of course, the line element! Tensor components refer to the coordinate frame and are numbered 1,...,Dim, rather than 0,...,Dim-1, so that the {i,j...k} component of a tensor is the [[i,j...k]] part of the corresponding nested list. See section 5, "Calculating Tensor Components in an Arbitrary Frame", on how to do calculations in an arbitrary frame.

Note: One can choose to skip the most time consuming computations: RUddd and/or Wdddd. This can be achieved by calling RGtensors  with an optional extra argument -- a list of two numbers {i,j}. When the first number is 0, RUddd is not computed and when the second number is 0 Wdddd is not computed. Thus RGtensors[gIN,xIN,{1,0}]  will compute RUddd but not Wdddd. (See Examples 7 and 8)

1.1.1  Totally antisymmetric symbol / tensor

Definitions for the totally antisymmetric symbol (Subscript[ε, a1...an]) and tensor  (Subscript[η, a1...an]) are provided, but the tensor is not evaluated unless called by the user. As the number of indices on these objects equals the dimension of the mannifold, the U/d notation cannot be used. Instead, the following global names are used:

epsilon is a function that takes Dim integers as arguments and returns one of { 0, 1, -1} - the sign of the permutation of the arguments. For example, in 5 dimensions,

epsilon[1,2,3,4,5] returns 1,  epsilon[1,3,2,4,5] returns -1, and  epsilon[1,2,2,4,5] returns 0.

eta[] can take 0 to Dim arguments and returns the Dim-rank antisymmetric tensor with raised indices at the position(s) indicated by the argument(s). Thus eta[] is the totally covariant antisymmetric tensor, with components

eta[][[a1,a2,...,an]]=(((|detg|))^(1/2))epsilon[a1,a2,...,an],

while

eta[2,3][[a1,a2,...,an]]=(((|detg|))^(1/2))(g^(a2 b2))(g^(a3 b3))epsilon[a1,b2,b3,..,an].

Warning: The evaluation of eta[] tensors (Dim-fold nested lists) is very time consuming in high dimensions as the number of components grows as (Dim^Dim). It also requires large amounts of memory. For example, for the simple case of a conformally flat metric, the {Time, Storage} requirements for computing each eta[] grow as follows: {~1 sec, <1 MB} for Dim=6, {~20 sec, ~7 MB} for Dim=7, and {> 10 min, >120 MB} for Dim=8! Thus, calculations with eta[] should not be attempted if Dim>7. Fortunately, for Dim=4 the requirements are quite modest, even for complicated metrics.

Once the basic tensors have been computed, one can use the functions defined in the following sections to perform further calculations.

1.2  Manipulation of tensor indices & tensor products

The first or last index of a tensor can be raised / lowered by inner (Dot) multiplication with gUU / gdd:

Edd=gdd.EUd has components Edd[[a,b]] = (Subscript[E, ab])

WdddU=Wdddd.gUU has components WdddU[[a,b,c,d]] = ((Subscript[W, abc])^d).

To manipulate other indices the following functions have been defined:

Lower[x_,ni__]:: lowers the index (indices) at position(s) ni (assumed U) of the nested list x using the matrix  gdd.

Raise[x_,ni__]:: raises the index (indices) at position(s) ni (assumed d) of the nested list x using the matrix  gUU.

For example, the mixed Riemann tensor with  components (Subscript[((Subscript[(R^a), b])^c), d])  is  RUdUd=Raise[RUddd,3], while  RUUUd=Raise[RUddd,2,3]  has components (Subscript[(R^abc), d]).

Note: The functions Raise / Lower apply the rules in simpRules and factor the expressions obtained. Thus Raise[RUddd,4] is equivalent but not identical to RUddd.gUU.

Tensor products are obtained by using the Mathematica function Outer. For example,

gg=Outer[Times,gdd,gdd]; is a 4th rank (d-type) tensor with components gg[[a,b,c,d]]= (Subscript[g, ab])(Subscript[g, cd]).

Using the Mathematica function Transpose,  one can then define the tensor: gxg=gg-Transpose[gg,{1,3,2,4}], which has components  gxg[[a,b,c,d]] = (Subscript[g, ab])(Subscript[g, cd ])- (Subscript[g, ac])(Subscript[g, bd]).

1.3  Contraction of tensor indices

The function  Contract[x_,{i1_,j1_},{i2_,j2_},..]  takes as arguments an nth rank tensor x, and one or more lists of  pairs of numbers, indicating the positions of the indices to be contracted. It returns a tensor of rank n-2k, where k is the number of lists. It is assumed that the two numbers in each list refer to one U and one d type index, in any order.  Examples (the Einstein summation convention holds):

(a) Contract[GUdd,{1,3}]  gives a d-type vector with components (Subscript[Γ, ba]^a).

(b) Contract[RUddd,{1,3}]  gives Rdd.

(c) Contract[Outer[Times,RUddd,RUddd],{1,7},{4,5}]  produces a dddd tensor with components (Subscript[(R^s), abt])(Subscript[(R^t), csd]).

Note1: Even though the last result can be obtained by calling Contract twice, each time with a single pair of indices, it is much more efficient to perform multiple contractions simultaneously, as above.

Note2: The result of (c) can also be obtained as  Contract[RUddd.RUddd,{1,5}], since the Mathematica function Dot performs the {4,5} contraction of the outer product (last index of first tensor with first index of second tensor) giving a 6th rank tensor.

The function multiDot[x_,y_,{i1_,j1_},{i2_,j2_},..] is similar to Contract, and generalizes the built-in function Dot.  It takes as arguments two tensors -- x, y -- and one or more lists of  pairs of numbers.  For each pair of numbers {i,j}, it contracts the ith index of the first tensor (x) with the jth index of the second tensor (y),  returning a tensor of rank m+n-2k, where m, n are the ranks of x, y, respectively, and k is the number of index pairs. Again, it is assumed that the indices at {i,j} are of opposite (U/d) types. Using multiDot, the result of example (c) above can also be obtained as  multiDot[RUddd,RUddd,{1,3},{4,1}].

Note1: When contracting a product of two tensors, as in example (c), one should always use multiDot to perform the contractions, thus avoiding the outer product calculation, which is very time / memory consuming. See also Example 1a in section 3 and Example 4b in Section 5.1

Note2: One can use Contract and multiDot to contract nested lists with unequal index dimensions (provided, of course, that the indices to be contracted have the same dimension); thus, if A has dimensions {2, 3, 3, 4} and B {3, 2, 4}, Contract[A, {2, 3}] will return a tensor with dimensions {2, 4} and multiDot[A, B, {1, 2}, {3, 1}] will return a tensor with dimensions {3, 4, 4}.

1.4  Covariant differentiation of tensors

The following functions, involving differentiations of tensors, are defined by RGTC:

(i) The function  covD[x_,{i___}]  takes two arguments -- an nth rank tensor, x, and a list of numbers indicating the positions of the contravariant  (U) components of x -- and produces an n+1 rank tensor: the covariant derivative of x. The last index is the differentiation index (d type). If x has only d-type components, the second argument (empty list) can be omitted. Examples:

covD[Rdd] is a 3rd rank (ddd type) tensor with components (Subscript[R, ab;c]).

covD[RUddd,{1}] is a 5th rank (Udddd type) tensor with components (Subscript[(R^a), bcd;e]).

(ii) The function covDiv[x_,{i___,{s_}}] acting on an nth rank tensor x with contravariant indices at positions {i..., s} returns an n-1 rank tensor: the covariant divergence of x with respect to the contravariant index at position s. If there is only one contravariant index, the second argument can be written simply {s} rather than {{s}}. Examples (the Einstein summation convention holds):

covDiv[RUddd,{1}] is a 3rd rank (ddd type) tensor with components (Subscript[(R^s), abc;s]).

covDiv[RUdUd,{1,{3}}] is a 3rd rank (Udd type) tensor with components (Subscript[((Subscript[(R^a), b])^s), c;s]).

Note1: It is much more efficient to use covDiv  than to use covD  followed by Contract.

Note2: The functions covD  and covDiv  can be very slow on high-rank tensors in high dimensions. However, they remember their values, so calling them again with the same arguments takes no time.

(iii)  The function  LieD[h_,x_,{i___}] computes the Lie derivative of the tensor x, having contravariant indices at index positions i__, with respect to the contravariant (U type) vector h. For totally covariant tensors x (or scalars), the last argument (empty list) can be omitted. The result of LieD is a tensor of the same rank (and type) as x. As with covD  and covDiv,  LieD remembers its values. Examples:

LieD[{h1,h2,..,hn},gdd] is a 2nd rank (dd type) tensor with components (h^s)(Subscript[g, ab,s])+(Subscript[(h^s), ,a])(Subscript[g, sb])+(Subscript[(h^s), ,b])(Subscript[g, as]).

LieD[{h1,h2,..,hn},RUddd,{1}] is a 4th rank (Uddd type) tensor with components (h^s)(Subscript[(R^a), bcd,s])+(Subscript[(h^s), ,b])(Subscript[(R^a), scd])+(Subscript[(h^s), ,c])(Subscript[(R^a), bsd])+(Subscript[(h^s), ,d])(Subscript[(R^a), bcs])- (Subscript[(h^a), ,s])(Subscript[(R^s), bcd]).

(iv) The function  Laplacian[x_,{i___}] computes the Laplacian ((Subscript[=(x^(;s)), ;s])) of the tensor x, having contravariant indices at index positions i__. Again, for totally covariant tensors x, the last argument (empty list) can be omitted. The result of Laplacian is a tensor of the same rank (and type) as x.

(v) The function  Grad2Norm[x__] takes either 1 or 2  scalar functions as arguments and computes the square norm of the gradient of x (inner product): Grad2Norm[x]= (g^ab)(Subscript[x, ,a])(Subscript[x, ,b]) and  Grad2Norm[x,y]= (g^ab)(Subscript[x, ,a])(Subscript[y, ,b]).

(vi) Using covD and covDiv,  the function Bianchi[n_], for n = 0, 1, 2, computes the full, once, twice contracted Bianchi identities and returns a list containing any non-vanishing terms. This list is normally the empty list, except when the simplification rules used imply relations among the 3rd derivatives of the metric functions (see Example 5).

Evaluation of the Bianchi Identities takes long: for the metrics in the examples, the time needed to evaluate Bianchi[1] is ~ 1 - 15 times the time taken by RGtensors to compute all basic tensors. Bianchi[0] takes considerably longer. (For the Bondi metric the times are ~ 12 min, 6 min, 7 sec).

1.5  Index of global names / variables

It is essential that, in choosing coordinates and metric functions, one avoids using the symbols listed in this subsection, as they are defined by (or have special meanings in) RGTC. In particular, the  single-letter symbols R, d, e must be avoided-- see subsection 4.0.

(i) Riemannian Geometry tensors / scalars

gdd, gUU, GUdd, Rdddd, RUddd, Rdd, EUd, R, Wdddd, detg		(section 1.1)
	epsilon, eta								(section 1.1.1)
	ωUd 									(section 5.1)
	α,β,...,  𝔇[_],Δ[_]..., Ψ[i], Φ[i, j], Λ				(section 5.1.3)

(ii) Tensor operations / functions

RGtensors							     (section 1.1, 5.1, 5.2)
	Raise, Lower, Outer, Transpose					(section 1.2) 
	Contract, Dot, multiDot						(section 1.3) 
	covD, covDiv, LieD, Laplacian, Grad2Norm, Bianchi		(section 1.4)
	HStar									(section 5.1.2)
	Classify, Plebanski							(section 7)

(iii) Other global variables

Dim, coordList, LineElement					(section 1.1)
	simpRules, TrigRules, HypTrigRules, varList			(section 4.3)
	e[], d[]									(section 5.1)
	deList, dxRuleList							(section 5.2)

(iv) Auxiliary functions

zeroQ, nonZeroN, nonZeroL, indepTerms, indepCoefs, FacSimp	(section 4.2)
	FuncRepRules							    	    (section 4.4)
	eTO$dx, Clear$dx								(section 5.1.1)
	GenCoef, metric								(section 4.2)

(v) Functions / variables defined in EDC

e.g., Wedge, d[], DeclareForms, Bar[], DirectSum etc.		see  EDCmanual.nb

Important Note: As the behavior of RGtensors depends on the values of global variables like simpRules, deList, dxRuleList, one should normally restart the kernel to do calculations with a new metric.

1.6  Input Palettes

The notebooks OperatorPLT.nb and NPsymbolPLT.nb are input palettes. The first contains the more common functions/variables used in calculations, so they can be entered easily without misspellings, while the second contains the symbols for the Newman-Penrose variables (see section 5.1.3).

It is assumed that the user has a working knowledge of the basic Mathematica structures used in symbolic calculations -- in particular that he/she understands the properties and use of Lists, Patterns, Functions, Transformation Rules and Definitions.  Although not necessary for elementary evaluations, familiarity with these structures is needed to make full use of all features of RG&TC.

The executable cells in this notebook must be evaluated in the order they appear, after placing the file EDCRGTCcode.m in an appropriate directory. Also read subsection 4.0 before doing calculations with a new metric.

2  RG&TC-Code

<<EDCRGTCcode.m

The file EDCRGTCcode.m  must be placed in one of the directories that Mathematica looks at when searching for extra packages (e.g., /Applications/Mathematica/AddOns/AutoLoad). Such directories can be found by evaluating the command $Path:

$Path

EDCRGTCcode.m  is the combined EDC and RGTC code, allowing calculations with differential forms, symbolic matrices, etc. (See the  EDCmanual.nb).

3  Examples (~ 1.5 min) 
(All timings refer to an Apple iBook G3/900 computer running Mathematica 3.0 under Mac OS 9.2.2)

In each case, after RGtensors is finished, examples of further calculations, using the other functions defined in the package, are given.

stTim=AbsoluteTime[];

Example 0: Flat spacetime in rotating prolate spheroidal coordinates

Transform the flat line element from cylindrical coordinates (ρ, z, phi) to (λ, μ, ϕ), where (ρ^2)=((λ^2)-1) (1-(μ^2)), z=λ μ, phi=ϕ+Ω t:

xCoord={t,φ,λ,μ};

gFLAT={{1-Ω^2 (λ^2-1)(1-μ^2),-Ω (λ^2-1)(1-μ^2),0,0},{-Ω (λ^2-1)(1-μ^2),-(λ^2-1)(1-μ^2),0,0},{0,0,-(λ^2-μ^2)/(λ^2-1),0},{0,0,0,-(λ^2-μ^2)/(1-μ^2)}};

RGtensors[gFLAT,xCoord]

GUdd[[3,1,2]]

Check covariant derivatives of metric tensor:

covD[gdd]

covD[gUU,{1,2}]

Example 1: Schwarzschild

xCoord={t,r,θ,φ};

g=DiagonalMatrix[{1-2M/r,-1/(1-2M/r),-r^2,-r^2Sin[θ]^2}];

simpRules=TrigRules;

TrigRules are predefined rules to implement the identities ((Cos[θ])^2) + ((Sin[θ])^2)=1 and ((Csc[θ])^2) - ((Cot[θ])^2)=1.

RGtensors[g,xCoord]

RUddd[[1,2,1,2]]

Check the full Bianchi identities:

zeroQ[Bianchi[0]]//Timing

Example 1a: Schwarzschild with Cosmological Constant

xCoord={t,r,θ,φ};

g=DiagonalMatrix[{1-2M/r-Λ r^2/3,-1/(1-2M/r-Λ r^2/3),-r^2,-r^2Sin[θ]^2}];

simpRules=TrigRules;

RGtensors[g,xCoord]

Rdd[[3,3]]

R

Factor[Rdd-R/Dim gdd]

Check the once-contracted Bianchi identities:

zeroQ[Bianchi[1]]

Use multiDot and Contract to compute the quadratic and cubic Weyl tensors and invariants:

Off[General::"spell",General::"spell1"];

WUUdd=Raise[Wdddd,1,2];

(** W2=Weyl squared **)W2UUdd=multiDot[WUUdd,WUUdd,{3,1},{4,2}];
(** W3=Weyl cubed **)W3UUdd=multiDot[W2UUdd,WUUdd,{3,1},{4,2}];

quadraticInv=Contract[W2UUdd,{1,3},{2,4}];
cubicInv=Contract[W3UUdd,{1,3},{2,4}];

{quadraticInv,cubicInv,quadraticInv^3/cubicInv^2}

The invariants can be computed directly using multiDot:

{multiDot[WUUdd,WUUdd,{1,3},{2,4},{3,1},{4,2}],multiDot[W2UUdd,WUUdd,{1,3},{2,4},{3,1},{4,2}]}

On[General::"spell",General::"spell1"];

Example 2: Robertson-Walker

xCoord={t,r,θ,φ};

g={{1,0,0,0},{0,-a[t]^2/(1-k r^2),0,0},{0,0,-a[t]^2r^2,0},{0,0,0,-r^2a[t]^2Sin[θ]^2}};

simpRules=TrigRules;

RGtensors[g,xCoord]

Rdd[[1,1]]

EUd[[1,1]]

Check the twice-contracted Bianchi identities:

zeroQ[Bianchi[2]]

Example 3: Conformally Flat 5-dim metrics

xCoord={x1,x2,x3,x4,x5};

g=Ω[x1,x2,x3,x4,x5]^2IdentityMatrix[5];

simpRules={};

(1) General Case:

RGtensors[g,xCoord]

GUdd[[2,2,2]]

EUd[[1,1]]

Of the 5^4 components in Rdddd

nonZeroN[Rdddd]

are non-zero, of which

Length[nonZeroL[Rdddd]]

are different. Of these

Length[indepTerms[Rdddd]]

are independent (modulo numerical factors).

Useful Hint: To find which components of Rdddd equal a particular expression in nonZeroL, use the Mathematica function Position:

Position[Rdddd,nonZeroL[Rdddd][[1]]]

Rdddd[[1,4,1,5]]===nonZeroL[Rdddd][[1]]

Rdddd[[2,4,2,5]]===nonZeroL[Rdddd][[1]]

(2) Constant Curvature Case:

RGtensors[g/.{Ω[x1,x2,x3,x4,x5]->1/(1+K/4(x1^2+x2^2+x3^2+x4^2+x5^2))},xCoord]

RUddd[[1,2,1,2]]

R

Of the 5^4 components in RUddd

nonZeroN[RUddd]

are non-zero, of which the following are independent (modulo numerical factors):

indepTerms[RUddd]

Example 4: Kerr in Boyer-Lindquist coordinates (Cos[θ] -> x): Finding Killing Vectors

Examples 4a, 4b, 4c are in section 5.1

Clear[x];xCoord={t,r,x,φ};

Sigma=r^2+a^2x^2;Delta=r^2+a^2-2M r;

g={{1-2M r/Sigma,0,0,2a M r(1-x^2)/Sigma},{0,-Sigma/Delta,0,0},{0,0,-Sigma/(1-x^2),0},{2a M r(1-x^2)/Sigma,0,0,-(1-x^2)(r^2+a^2+2a^2M r(1-x^2)/Sigma)}};

simpRules={};

RGtensors[g,xCoord]

GUdd[[2,2,3]]

RUddd[[1,2,1,3]]

Check the once-contracted Bianchi identities:

zeroQ[Bianchi[1]]//Timing

Bianchi[1] calls covD and covDiv which remember their values; therefore, calling it again takes no time:

zeroQ[Bianchi[1]]//Timing

Finding Killing Vectors

Assume (r, x) dependence for the contravariant components of ξ:

ξU={ξ[1][r,x],ξ[2][r,x],ξ[3][r,x],ξ[4][r,x]};

covDξdd=covD[gdd.ξU];//Timing

KillingEquations=Collect[indepTerms[Flatten[FacSimp[covDξdd+Transpose[covDξdd]]]],{Derivative[__][ξ[_]][r,x],ξ[_][r,x]},Factor]

Length[%]

Solve[Take[KillingEquations,2]==0,{ξ[2][r,x],ξ[3][r,x]}]

The function FuncRepRules generates replacement rules for replacing a function (and its derivatives) by its values; see section 4.4.

Collect[indepTerms[FacSimp[KillingEquations/.FuncRepRules[ξ[2][r,x],0]/.FuncRepRules[ξ[3][r,x],0]]],{Derivative[__][ξ[_]][r,x],ξ[_][r,x]},Factor]

Solve[%==0,{Derivative[1,0][ξ[1]][r,x],Derivative[1,0][ξ[4]][r,x],Derivative[0,1][ξ[1]][r,x],Derivative[0,1][ξ[4]][r,x]}]

FacSimp[KillingEquations/.FuncRepRules[ξ[1][r,x],C1]/.FuncRepRules[ξ[2][r,x],0]/.FuncRepRules[ξ[3][r,x],0]/.FuncRepRules[ξ[4][r,x],C4]]

Lie Derivatives

LieD[h_,x__] ccomputes the Lie derivative of the tensor x with respect to the contravariant vector h:

TimeLikeKV={1,0,0,0}

AxialKV={0,0,0,1}

GeneralKV=C1 TimeLikeKV+C2 AxialKV

LieD[TimeLikeKV,gdd]

LieD[AxialKV,RUddd,{1}]

LieD[GeneralKV,gUU,{1,2}]

Example 5: Floratos-Kehagias 8-dim Ricci-flat metric (hep_th/9802107)

Clear[u,g,h,a,b,c];

xCoord={x1,x2,x3,x4,x5,x6,x7,x8};

gFloKeh={{u^2 g^2 h^2,0,0,0,0,0,0,0},{0,u^2 g^2 h^2,0,0,0,0,0,0},{0,0,1/u^2,0,0,0,0,2a/u^2},{0,0,0,1/g^2,2b/g^2,0,0,0},{0,0,0,2b/g^2,g^2+4b^2/g^2,0,0,0},{0,0,0,0,0,h^2+4c^2/h^2,2c/h^2,0},{0,0,0,0,0,2c/h^2,1/h^2,0},{0,0,2a/u^2,0,0,0,0,4a^2/u^2+u^2}};

The symbols u, g, h, a, b, c are functions of  x1, x2 :

u=U[x1,x2];g=G[x1,x2];h=H[x1,x2];a=A[x1,x2];b=B[x1,x2];c=C[x1,x2];

The complex functions u^2 +2 i a, g^2 +2 i b, h^2 +2 i c  are assumed to be analytic functions of x1+i x2. Therefore, use simpRules to implement the Cauchy-Riemann & harmonic function conditions:

simpRules={D[a,x1]->-u D[u,x2],D[a,x2]->u D[u,x1],D[b,x1]->-g D[g,x2],D[b,x2]->g D[g,x1],D[c,x1]->-h D[h,x2],D[c,x2]->h D[h,x1],Z_.Derivative[0,2][G_][x1,x2]+Z_.Derivative[2,0][G_][x1,x2]:>-Z(Derivative[1,0][G][x1,x2]^2+Derivative[0,1][G][x1,x2]^2)/G[x1,x2]}

RGtensors[gFloKeh,xCoord]

zeroQ[Bianchi[1]] in this case gives False (after ~30 sec)! The reason is that the simpRules used imply relations among the 3rd derivatives of the functions U[x1,x2], G[x1,x2], H[x1,x2].

Exercise: Find the additional two rules needed to satisfy the Bianchi identities.

Solution: See AdditionalExamples.html

Clear[u,g,h,a,b,c,U,G,H]

Example 6: Harrison metric II-A-7 (Phys. Rev. 116, 1285)

See the "Exact Solutions..." book: section 15.4 (1st Edition) or section 17.3 (2nd Edition)

Clear[x];xCoord={t,x,y,z};

rt3=Sqrt[3];F1=x/rt3-t^2/12 ;F2'[z]=F2[z]v[z]/z;v'[z]=(v[z]^2-1)/(4z)(2  v[z]/(z-1)+4/rt3);

n1[t]=2+rt3;n1[x]=1+rt3;n1[y]=-rt3;n1[z]=3+rt3;

n2[t]=1+2/rt3;n2[x]=1+2/rt3;n2[y]=-1;n2[z]=1+2/rt3;

n3[t]=1/rt3;n3[x]=(1+1/rt3);n3[y]=-1/rt3;n3[z]=(2+1/rt3);

g=DiagonalMatrix[{F1^n1[t]F2[z]^n2[t](-1/z)^n3[t],-F1^n1[x]F2[z]^n2[x](-1/z)^n3[x],-F1^n1[y]F2[z]^n2[y](-1/z)^n3[y],-F1^n1[z]F2[z]^n2[z](-1/z)^n3[z](v[z]^2-1)/(z-1)}];

simpRules={};

RGtensors[g,xCoord]

detg

Mathematica 4 (or higher) is able to factor the long polynomial ((t^12)...)  appearing in detg  as  ((((t^2)-4 (Sqrt[3]) x))^6) (and also the other polynomials in gUU). As a result the calculation is completed in 1/3 the time!

Example 7: Bondi metric (Proc. Roy. Soc. A269, 21)

xCoord={u,r,θ,φ};

g={{Exp[2 β[u,r,θ]]V[u,r,θ]/r-r^2Exp[2γ[u,r,θ]]U[u,r,θ]^2,Exp[2 β[u,r,θ]],r^2Exp[2γ[u,r,θ]]U[u,r,θ],0},{Exp[2 β[u,r,θ]],0,0,0},{r^2Exp[2γ[u,r,θ]]U[u,r,θ],0,-r^2Exp[2γ[u,r,θ]],0},{0,0,0,-r^2Exp[-2γ[u,r,θ]]Sin[θ]^2}};

simpRules=TrigRules;

RGtensors[g,xCoord]

GUdd[[2,2,3]]

Rdd[[2,2]]

EUd[[2,3]]

covD[gdd]

LieD[{0,0,0,1},Rdd]

Verification of the Bianchi identities for this metric takes ~ 19 min (12.5 min + 6 min + 7 sec).

One can choose to skip the time consuming RUddd and/or Wdddd computations by calling RGtensors wirh a 3rd argument:

RGtensors[g,xCoord,{0,0}]

Example 8: Kerr with Cos[θ],  Sin[θ]

Transform Boyer-Lindquist t, φBL coordinates:  d[t]->d[u]+(r^2+a^2) d[r]/(r^2-2M*r+a^2), d[φBL]->d[φ]+a d[r]/(r^2-2M*r+a^2);

xCoord={u,r,θ,φ};

g={{1-((2 M r)/((r^2)+(a^2) ((Cos[θ])^2))),1,0,((2a  M r Sin[θ]^2)/((r^2)+(a^2) ((Cos[θ])^2)))},{1,0,0,-a  Sin[θ]^2},{0,0,-(r^2+a^2Cos[θ]^2),0},{((2a  M r Sin[θ]^2)/((r^2)+(a^2) ((Cos[θ])^2))),-a  Sin[θ]^2,0,-Sin[θ]^2(2 (a^2) M r+(a^2) (r^2)+(r^4)+(a^2) ((a^2)-2 M r+(r^2)) ((Cos[θ])^2))/((r^2)+(a^2) ((Cos[θ])^2))}};

(1)  Use TrigRules in simpRules and do not compute theWeyl tensor:

simpRules=TrigRules;

RGtensors[g,xCoord,{1,0}]

Note that, using TrigRules, RGtensors is unable to verify that the metric is Ricci-Flat. Rdd has two independent non-vanishing components:

indepTerms[Rdd]

Simplify[%]

RUddd[[1,2,1,3]]

(2)  To simplify denominators in gUU and elsewhere, use extra rule in simpRules:

simpRules=Join[TrigRules,{((Sin[θ])^(n_)):>((Sin[θ])^(n-2))(1-((Cos[θ])^2))/;n>1}]

RGtensors[g,xCoord]

GUdd[[2,2,3]]

RUddd[[1,2,1,3]]

Note the speed improvement with the extra simplification rule!

Print["Total time = ",Round[AbsoluteTime[]-stTim]," sec"]

Sections 5 and 6 contain examples illustrating more advanced features of RG&TC. The examples in each subsection must be evaluated in sequence as some definitions given in the first one are also used later.
	
	More complicated examples can be found in AdditionalExamples.html.

4  Usage Tips & Auxiliary Functions

4.0  Do not use symbols defined in RGtensors

The most common cause of unexpected behavior or meaningless results is the inappropriate use, in the coordinates or the metric, of a symbol that is defined by RGtensors. For example, if one's metric contains the function R[x,y...], RGtensors will go in an infinite loop trying to define the scalar curvature, denoted by R, in terms of the metric function R[x,y...]! The same will happen if the symbol R is used as a coordinate. The situation is worse when an NP frame is used, where several more single-letter symbols are defined by RGtensors (the spin coefficients and Λ, Φ  and Ψ) or have special meanings (d, e). One must, therefore, make sure that the symbols used in one's coordinates, metric or frame are distinct from (or are consistent with) the global variables listed in section 1.5.

4.1  Avoid displaying long expressions & keep intermediate results

The most important thing to keep in mind when working with nested lists is that one can easily generate huge outputs that can tax the graphics capabilities of any system and may even cause Mathematica to crash. Besides, not much is gained by viewing on screen deeply nested lists containing mostly zeros! Thus, as a general rule, never evaluate an input line that does not end in a semi-colon(;). As a safeguard, the printed form of any expression that exceeds 200 KB is shortened, occasionally resulting in the output <<1>> (which means "too big to be displayed"). One can still use this quantity in calculations, however, and display parts of it.

One must, therefore, get used to working with symbols representing tensors and only occasionally view individual components or subparts. Also, because operations on lists are time-consuming, one must save intermediate results in case they are needed later on. For these reasons, a consistent system of naming tensors is necessary to keep track of what has been calculated.

The U/d naming convention used here is reasonably short and visually conveys the properties (rank, index type) of the tensor represented. For contractions, the convention can be extended by using the same letter in upper/lower case (S, s) for pairs of contracted indices. Finally, for covariant derivatives, one can use the letter j (resembling ;) to separate covariant differentiation indices. For example, the symbol RUdddjd is an obvious candidate for covD[RUddd,{1}], while RSdddjs can be used for covDiv[RUddd,{1}].

Another notation, similar to the one used for eta (Section 1.1.1), is to define an appropriately named function for each tensor, with arguments the positions of the contravariant indices. For example, for the Riemann tensor, the definition

Riem[]=Rdddd;Riem[a__]:=Riem[a]=Raise[Rdddd,a];

would enable one to use the notation Riem[1] for the tensor RUddd, Riem[3,4] for the tensor RddUU, etc. Of course, all these are just names for referring to nested lists stored in memory; any name that is meaningful to the user can be used!

4.2  Auxiliary functions for examining and simplifying lists

The following auxiliary functions have been defined and are useful in examining high rank tensors:

(1) zeroQ[x_] tests if the tensor x is the zero tensor, returning True or False.

(2) nonZeroN[x_] returns the Number of non-zero components of the tensor x.

(3) nonZeroL[x_] returns a List containing the distinct non-zero components of the tensor x (if Z is any expression, Z,-Z,2Z  are considered distinct). To find which components of x equal a particular element of nonZeroL[x], use the Mathematica function Position. Thus, Position[Rdddd,nonZeroL[Rdddd][[1]]] will produce a list of the indices -- components of Rdddd -- which are equal to the second argument of Position. See Example 3 (1).

(4) indepTerms[x_] returns a list of the "independent" components of the tensor x. Independent here means not equal to  ±  a numerical multiple of another component. Note that indepTerms does not perform any reductions:

indepTerms[{2a+2b,a+b}]

indepTerms[{2(a+b),a+b}]

nor does it test for linear dependence:

indepTerms[{a-b,b-c,c-a}]

Plus@@%

indepTerms can take an optional second argument - a list of symbols or expressions which, when appearing as overall factors, can be ignored in testing equivalence.  Consider, for example, the following list S (Z can be any expression):

S={Z,(a^2)Z,2a*M*Z,-2Sin[th]Cos[th]Z,((Sin[th])^2)Z}

indepTerms[S]

indepTerms[S,{a,Sin[th]}]

indepTerms[S,{a,M,Sin[th],Cos[th]}]

This option can be useful in reducing the number of independent components of a tensor.

(5) The function indepCoefs[x_]  returns a list of the "independent" coefficients of all symbolic differential forms or wedge-products in the elements of the list x. Independent here means not equal to  ±  a numerical multiple of another coefficient.

The symbols e[i] have been defined as symbolic 1-forms in RGTC.

indepCoefs[a[1] e[1]+a[2] e[2]+a[3] e[3]]

indepCoefs[(a[2]+a[3]) e[2]⋀e[3]+(a[3]+a[1]) e[3]⋀e[1]+(a[1]-a[2]) e[1]⋀e[2]]

Like indepTerms, indepCoefs does not test for linear dependence:

Plus@@%

indepCoefs has no effect on expressions that are not differential forms:

indepCoefs[a[1] b[1]+a[2] b[2]+a[3] b[3]]

indepCoefs can be applied to lists:

sym33mat=Table[(i+j)(a[i]-a[j])e[i]⋀e[j],{i,3},{j,3}]

indepCoefs[%]

Note that indepTerms does not recognize that (-a[1]+a[2])=-(a[1]-a[2]), as it does not carry out any simplifications:

indepTerms[sym33mat]

indepTerms[FacSimp[sym33mat]]

(6) FacSimp[x_] applies the simplification rules defined in simpRules and factors each component of the tensor x. All basic tensors calculated initially by RGtensors as well as those returned by the functions Raise, Lower, Contract, multiDot, covD, covDiv , LieD, Laplacian , Grad2Norm and Bianchi have been simplified using FacSimp. Even though it can be very slow on high-rank tensors and /or when simpRules contains many rules, FacSimp should always be used to simplify sums and products (Dot or Outer) of these tensors. (see section 4.5 "Taking advantage of symmetries")

If zeroQ[expr] on an expression that is expected to vanish gives False,  try zeroQ[FacSimp[expr]] or,  for high-rank tensors with symmetries,  zeroQ[FacSimp[indepTerms[expr]]].

(7) GenCoef[x_,y_] ("Generalized Coefficient") is very similar to the Mathematica function Coefficient[x_, y_]. Unlike Coefficient however, GenCoef returns a series when x is a series and, when both x and y are Lists, returns a nested List:

seriesExpr=Series[Sum[(a+b)^k*x^k,{k,0,3}],{x,0,3}]

Coefficient[seriesExpr,b]

Coefficient[seriesExpr,b,2]

GenCoef[seriesExpr,b]

GenCoef[seriesExpr,b,2]

Coefficient[{a[1]x+a[2]y,b[1]x+b[2]y},{x,y}]

GenCoef[{a[1]x+a[2]y,b[1]x+b[2]y},{x,y}]

Thus GenCoef[LinearEquations_List,unknowns_List] gives directly the matrix of coefficients of the unknowns in the given equations.

(8) metric[ds2_,xIN_] constructs the metric tensor from the line element (ds2) and the list of coordinates (xIN):

metric [2d[u]d[r]+(1-2M/r)d[u]^2-r^2d[θ]^2-r^2Sin[θ]^2d[ϕ]^2,{u,r,θ,ϕ}]

4.3  Choosing simpRules and using varList

It is impossible to overemphasize the importance of choosing appropriate simplification rules in the global variable simpRules to be applied by FacSimp. As FacSimp is used repeatedly by all functions, minor changes in simpRules can have dramatic effects.

RGTC  has two predefined sets of rules for basic simplification of expressions depending on trigonometric (ordinary or hyperbolic) functions: TrigRules and HypTrigRules. For metrics depending on such functions, one should set simpRules=TrigRules (or HypTrigRules).

When applying RGtensors to a new metric, it may not be clear what, if any, simpRules should be used. In such cases, start with no rules and abort RGtensors after Gamma has been computed. Look at the expressions for gUU and GUdd: apart from being factored, are they as simple as they can be? Could a simple substitution rule further simplify these expressions? If yes, include it in simpRules. As gUU and GUdd are used repeatedly by other functions, even minor improvements can have significant effects (see Example 8, and Example 4b in section 5.1).

Exercise: Find the simpRules that must be used in order that RGtensors prints Ricci Flat when applied to the Kruskal form of the Schwarzschild metric:

gKRUSKAL={{0,-8M^2(1-2M/r[u,v])/(u*v),0,0},{-8M^2(1-2M/r[u,v])/(u*v),0,0,0},{0,0,-r[u,v]^2,0},{0,0,0,-r[u,v]^2Sin[θ]^2}};

xCoord={u,v,θ,φ};

Where r[u,v] is related to the coordinates u,v by the equation (2M-r[u,v])Exp[r[u,v]/(2M)]-2M*u*v=0.

Solution: See AdditionalExamples.html

In formulating substitution rules it is better to use many specific rules rather than few general ones. The reason is that a general pattern will match more instances (even if it fails) and thus slow down the computation. For rules and patterns see the Mathematica Book, chapters 2.3 and 2.4 (=2.5 in version 5).

The global variable varList is a list of variables/patterns that, when non-empty, is used by FacSimp to group terms. It can be useful in cases like (see also Example 4c after section 5.1.3):

varList

expr=FacSimp[(1+x y)^2(x-y)^3+Sqrt[(1-x^2)]Sum[y^n,{n,5}]]

Simplify[expr]//Timing

FullSimplify[expr]//Timing

varList={(Sqrt[1-(x^2)])};FacSimp[expr]//Timing

varList={};

FacSimp[expr]

Warning: Inappropriate choice of varList variables can prevent cancellations and slow down RGtensors. It should be used with caution. A good check, before calling RGtensors, is to see the effect of particular varList variable(s) on the expression for the inverse metric FacSimp[Inverse[gIN]].

Note: The global variables simpRules and varList are assumed to be lists. Therefore, one must never use Clear[...] for these variables. One should instead set them equal to the empty list {}.

4.4  Substituting new functions

The function FuncRepRules[f[x__],g_,n_] ("Function Replacement Rules") generates a list of rules for replacing the function f[x] and all its derivatives up to order n by g and the corresponding derivatives of g, where g can be any expression. The third argument (n) is optional. When it is omitted, rules up to the second derivatives are generated. For example:

FuncRepRules[f[x],u[x]/x]

FuncRepRules[f[x],u[x]/x,3]

FuncRepRules[U[t,r],(f1[r-t]+f2[r+t])/r]

FuncRepRules[H[u,x,y],a[u](x^2-y^2)+2b[u]x y]

One can thus compute the Ricci tensor for a metric that depends on arbitrary functions of the coordinates and then use FuncRepRules to substitute specific functions / functional forms to search for solutions.

4.5  Taking advantage of symmetries

The greatest drawback of the functions defined in this package is that they do not take advantage of the symmetry properties of tensors to reduce the amount (and time) of a computation. (This is not true of the main procedure - RGtensors - but this is called only once!). Thus, FacSimp[Rdddd] factors every component of the Riemann tensor! As a result, high rank tensor calculations tend to be very slow, especially in high dimensions.

In such cases one should define a temporary function for the components of the tensor to be computed, having the correct symmetries, and then use Array to construct the tensor (see how Rdddd, RUddd, Rdd, Wdddd are computed by RGtensors). FacSimp should always be applied to the independent components only!

4.6  When the results appear wrong: Can I trust the results of RGTC?

The basic functions defined in RGTC have been tested on a large number of metrics using different frames and always give correct results. Successive versions have improved the performance and have corrected even some very rare bugs (like the one affecting HStar when acting on some SeriesData objects).  Of course, one can never claim that a piece of software is completely bug-free. However, if one encounters unexpected results, there are several checks one can make before concluding that there is a bug involved.

(i) Make sure no symbols defined by RGtensors are used in the coordinates or metric --  see section 4.0 above. Repeat the calculation using different notation and a fresh Kernel.

(ii) Common mistakes that can lead to erroneous results are: (a) in covD, covDiv, LieD, Laplacian, etc. -- not giving the correct positions of contravariant indices as second arguments, and (b) in multiDot or Contract -- giving wrong index positions or index pairs of the same type {U,U} or {d,d}. Double check such inputs.

(iii) Unexpected results can often be due to inadequate simplification: it is possible that RGtensors does not print Ricci Flat (or Conformally Flat) for a metric that is known to have this property; or that it fails to recognize that a given set of 2-forms deList and coordinate differential rules dxRuleList define a closed implicit coframe basis (see section 5.2) and thus aborts evaluation. This can happen when simpRules are used to simplify calculations in complicated metrics, particularly ones involving trigonometric functions. 

The case when Rdd or Wdddd does not get simplified sufficiently to give a zero tensor is easy to deal with: just use FacSimp or Simplify repeatedly to check if the expressions returned by RGtensors  actually do vanish (or repeat the calculation adding extra simplification rules as in example 8 of section 3). This is illustrated with the metrics on n-dimensional spheres, where RGtensors fails to print Conformally Flat, Einstein Space, and Space of Constant Curvature for n>4 -- see the notebook nDimSpheres.nb  in AdditionalExamples.html. 

In the case when RGtensors does not recognize that the given deList is closed, and therefore stops the calculation, it is best if one changes coordinates to eliminate the trigonometric functions. For an ad-hoc example where this happens see the notebook EliminateTrigFunctions.nb  in AdditionalExamples.html.

(iv) Differential Geometry provides a number of identities one can use to check the validity of the results returned by RGtensors and the other functions defined in RGTC (covD, Contract, etc.). Such identities follow from (1) the vanishing of the covariant derivatives of the metric tensor and the Bianchi identities, (2) properties of special vector fields (geodesics, Killing vectors) that may be present in particular metrics, and (3) agreement of the results obtained using different coordinate systems or frames or different ways of computing the same tensors / invariants. For an example of the latter see the notebook KerrInvariants.nb  in AdditionalExamples.html.

5  Tensor Components in an Arbitrary Frame

Warning: It is best to study the Examples in Section 3 and become familiar with calculations in the default coordinate-frame before reading this section. Otherwise, it may appear very confusing!

The computation of tensor components in an arbitrary frame involves some elementary operations of exterior differential calculus. We use the notation and conventions of the package Exterior Differential Calculus (EDC). Familiarity with that package is not required, unless one wants to do further calculations with the connection 1-form tensor ωUd which is computed by RGtensors.

The new intialization file EDCRGTCcode.m  evaluated in Section 2 above contains the combined EDC and RGTC code, allowing the complete functionality of the EDC package to be used.

5.1  Explicit Frame Calculations

In this case, one must

(i) Define, in addition to the coordinates (xIN) and the metric (gIN),  a list of 1-forms (eIN) -- the coframe -- written as linear combinations of the coordinate differentials d[xa]. The components of gIN must, of course, be given with respect to the coframe eIN.

(ii) Evaluate RGtensors with the coframe as a third argument: RGtensors[gIN,xIN,eIN].

RGtensors now computes all tensor components with respect to the given frame. (The same holds true, of course, for all other tensor functions - covD, covDiv, etc.). When it is finished, in addition to the other global tensors (Rdddd, Rdd, etc.), the connection 1-form tensor ωUd with components  ωUd[[a,b]] = (Subscript[(ω^a), b]) is now stored in memory. ωUd satisfies the equation  d[eIN]+ωUd⋀eIN= 0, or, in components, d[ea]+(Subscript[(ω^a), b])⋀eb= 0, where ea=eIN[[a]]. (Wedge multiplication of lists acts like Dot - it sums the last index of the first list with the first index of the second list). The notation e[a] is used for ea and all differential forms are written as linear combinations of the bases e[a], e[a]⋀e[b], etc.

Note: The exterior derivative operator d[] treats all symbols as variables. Therefore, constants must be explicitly defined to have zero exterior derivative, e.g., d[a]=0 and d[M]=0 for the Kerr metric. Otherwise, warning messages are generated and calculations take longer.

◼   Example 4a: Kerr (Example 4)  @ "rotating coordinate frame"

xCoord={t,r,x,φ};

coFrame={d[t]+2a M r(1-x^2)/(Sigma-2M r)d[φ],d[r],d[x],d[φ]};

Sigma=r^2+a^2x^2;Delta=r^2+a^2-2M r;

g=DiagonalMatrix[{1-2M r/Sigma,-Sigma/Delta,-Sigma/(1-x^2),-Delta Sigma(1-x^2)/(Sigma-2M r)}];

d[a]=0;d[M]=0;

simpRules={};

RGtensors[g,xCoord,coFrame]

detg

RUddd[[1,2,1,2]]

ωUd[[1,2]]

It is still possible to skip computation of RUddd and/or Wdddd by calling RGtensors with an extra (fourth!) argument of the form {i,j} (see note at the end of section 1.1). Thus RGtensors can be called with 2, 3 or 4 arguments, in the following order:

RGtensors[gIN,xIN] - computes all tensors in the coordinate frame

RGtensors[gIN,xIN,{1,0}] - computes all tensors except Wdddd in the coordinate frame

RGtensors[gIN,xIN,eIN] - computes all tensors in the co-frame specified by eIN

RGtensors[gIN,xIN,eIN,{0,0}] - computes all tensors except RUddd and Wdddd in the co-frame specified by eIN

5.1.1   Important note on repeated frame calculations:   eTO$dx  and  Clear$dx

After a calculation in a particular frame, the coordinate differentials are stored as linear combinations of the coframe basis e[a]. Thus, after Example 4a:

{d[t],d[t]⋀d[r]}

The global variable eTO$dx is a list of substitution rules, defined by RGtensors, that remembers the original frame definitions in terms of the coordinate differentials -- with 𝒹  used in place of d. It is useful in converting differential forms from the e[i]⋀e[j] basis to the d[xi]⋀d[xj] basis. For example

FacSimp[%/.eTO$dx]

Setting the d[x] equal to their expressions in terms of the e[a] makes calculations involving differentiations in this frame more efficient. However, if a subsequent calculation in a different frame is to be carried out, these values must be cleared. This can be done by evaluating:

Clear$dx

This procedure clears the coordinate differentials and eTO$dx. Now

{d[t],d[r]}

It is always safer, however, to reevaluate the initialization code EDCRGTCcode.m  -- which clears ALL variables -- before starting a calculation with a new frame (or metric).

5.1.2   Hodge Duality

A new function that can be used now is HStar[x_]:: applied to any p-form scalar or tensor expression x, it gives the "Hodge Dual" of x, i.e., the result of replacing each p-form basis e[a1]⋀e[a2]⋀...⋀e[ap] by (the Einstein summation convention holds)

(((|detg|))^(1/2))(g^(a1 b1))(g^(a2 b2))..(g^(ap bp))epsilon[b1,b2,..bp,bp+1,..,Dim]e[bp+1]⋀..⋀e[Dim]/(Dim-p)!

Note: The default volume Dim-form used as the definition of HStar[1] is  (((|detg|))^(1/2))e[1]⋀e[2]⋀...⋀e[Dim], except in the case of a Newman-Penrose null frame (see section 5.1.3) where the (real) volume 4-form  I e[1]⋀e[2]⋀e[3]⋀e[4] is used.  If in a particular case (to agree with the user's conventions) a different volume form must be used as the definition of HStar[1], one must define a new Mathematica function for Hodge-Dual, as a simple multiple of HStar: newHstar[x_]:= α HStar[x], where α is the ratio newVolumeForm / defaultVolumeForm.

◼   Example 4b: Kerr (Example 4)  @  orthonormal frame

Clear$dx

ONframe=FacSimp[DiagonalMatrix[PowerExpand[{Sqrt[(Sigma-2M r)/Sigma],Sqrt[Sigma/Delta],Sqrt[Sigma/(1-x^2)],Sqrt[Delta Sigma(1-x^2)/(Sigma-2M r)]}]].coFrame]

Note: PowerExpand must always be used when the argument of Sqrt is a product of factors.  Otherwise simplifications cannot be carried out.

gON=DiagonalMatrix[{1,-1,-1,-1}];

simpRules={(-1+x)(1+x)/(Sqrt[1-(x^2)])->-(Sqrt[1-(x^2)]),((Sqrt[(r^2)+(a^2) (x^2)])-(Sqrt[-2 M r+(r^2)+(a^2) (x^2)]))((Sqrt[(r^2)+(a^2) (x^2)])+(Sqrt[-2 M r+(r^2)+(a^2) (x^2)]))->2 M r};

This example illustrates the use of particular simpRules, chosen by looking at inadequately simplified expressions in GUdd when evaluated with no simpRules, as suggested in section 4.3.

RGtensors[gON,xCoord,ONframe,{0,0}]

detg

Rdddd[[1,2,1,2]]

ωUd[[1,2]]

◼   Sample calculations with tensor-valued differential forms and use of HStar:

• Construct the (rank-2) 2-form curvature tensor ΩUd with components ΩUd[[a,b]] = (Subscript[(Ω^a), b])=(1/2)(Subscript[(R^a), bcd])(e^c)⋀(e^d)  (eUU is the 2-form-basis tensor: eUU[[a,b]]=e[a]⋀e[b]).

eUU=Outer[Wedge,e/@{1,2,3,4},e/@{1,2,3,4}];Off[General::"spell",General::"spell1"];

ΩUd=multiDot[gUU.Rdddd/2,eUU,{3,1},{4,2}];//Timing

• Verify that  (Subscript[(Ω^a), b]) equals  d[(Subscript[(ω^a), b])]+(Subscript[(ω^a), s])⋀(Subscript[(ω^s), b]). (Wedge multiplication of tensors acts like Dot - contracts last index of first tensor with first index of second tensor).

zeroQ[FacSimp[d[ωUd]+ωUd⋀ωUd-ΩUd]]//Timing

• Verify Bianchi identities using ΩUd:

zeroQ[FacSimp[d[ΩUd]+ωUd⋀ΩUd-ΩUd⋀ωUd]]//Timing

• If both tensors to be contracted with multiDot are differential form-valued, then their components are multiplied using Wedge:

multiDot[ΩUd,ΩUd,{1,2},{2,1}]//Timing

• HStar gives the Hodge-dual:

pseudoScalar=HStar[%[[2]]]

• This invariant equals -(1/4)(Subscript[(R^ab), cd])(Subscript[(R^cd), ef])(Subscript[(η^ef), ab])  and can also be obtained (using eta[1,2]) as follows:

RUUdd=Raise[Rdddd,1,2];//Timing

RstarUUdd=multiDot[RUUdd/2,eta[1,2],{3,1},{4,2}];//Timing

multiDot[-RUUdd/2,RstarUUdd,{1,3},{2,4},{3,1},{4,2}]

Factor[%-pseudoScalar]

• The invariant  (1/2)(Subscript[(R^ab), cd])(Subscript[(R^cd), ab])  can be obtained either as:

multiDot[RUUdd,RUUdd,{1,3},{2,4},{3,1},{4,2}]/2

or, using ΩUd and HStar, as:

HStar[multiDot[ΩUd,HStar[ΩUd],{1,2},{2,1}]]//Timing

Factor[%[[2]]-%%]

On[General::"spell",General::"spell1"];

5.1.3   Newman-Penrose Null Frame

When the metric equals {{0,1,0,0},{1,0,0,0},{0,0,0,-1},{0,0,-1,0}}, RGtensors assigns values to a number of new symbols (global variables):
	(i)   the spin coefficents:  τ, κ, ρ, σ, γ, ϵ, α, β, ν, $π, λ, μ
	(ii)  the Ricci and Weyl scalars:  Φ[i,j], Λ, and Ψ[i] for the Newman-Penrose scalars (Subscript[Φ, ij]), Λ, and (Subscript[Ψ, i]),
and defines the differential operators:  {𝔇[_], Δ[_], δ[_], δB[_]}, 
assuming that the null basis given in the coframe is {𝕝,𝕟,𝕞,(OverBar[𝕞])} (standard NP notation). 

	RGtensors prints messages after each set of NP symbols is defined. Differential forms are still written in the e[i] notation, where e[1]=𝕝,e[2]=𝕟, etc. The only other differences from the standard notation are the use of $π in place of π (which in Mathematica is the constant 3.14159...), the use of the Gothic D - 𝔇 - for the NP differential operator D (to distinguish it from the built-in operator D) and the use of  δB instead of OverBar[δ] (easier to write and to  see!). The palette NPsymbolPLT.nb  facilitates entering these NP symbols/operators.

The function Bar now knows the reality conditions of the frame, e.g., that Bar[e[1]]=e[1], Bar[e[3]]=e[4],  etc., but assumes that all other symbols used are real (the spin coefficients are names that refer to expressions in terms of the coordinates which may contain the imaginary unit I). Thus, if complex symbols are to be used in an NP frame, their Bar[] relations must be defined before calling RGtensors.

◼   Example 4c: Kerr (Example 4)  @ NP[(1,2), (3,4)] frame

Clear$dx

NPframe=FacSimp[{(ONframe[[1]]-ONframe[[2]])/(Sqrt[2]),(ONframe[[1]]+ONframe[[2]])/(Sqrt[2]),(ONframe[[3]]+I ONframe[[4]])/(Sqrt[2]),(ONframe[[3]]-I ONframe[[4]])/(Sqrt[2])}]

gNP={{0,1,0,0},{1,0,0,0},{0,0,0,-1},{0,0,-1,0}};

varList={(Sqrt[(a^2)-2 M r+(r^2)])}

This is an example where using varList is essential (see section 4.3). With the default empty list, RGtensors seems to get stuck (more than 10 minutes with Mathematica 4) in computing Rdddd (Mathematica 3 actually gets stuck -- cannot be aborted!).

RGtensors[gNP,xCoord,NPframe,{0,0}]

Use the input palette NPsymbolPLT.nb to enter the spin coefficients.

κ

Bar[κ]

phi00=𝔇[ρ]-δB[κ]-(ρ^2+σ Bar[σ]+(ϵ+Bar[ϵ])ρ-Bar[κ]τ-κ(3α+Bar[β]-$π));//Timing

This takes long because the expressions for ρ and κ  are being differentiated.

ByteCount[phi00]

FacSimp[phi00]//Timing

Bar works as expected with exterior and partial differentiation:

drho=FacSimp[d[ρ]];Imrho=FacSimp[ρ-Bar[ρ]];

FacSimp[drho-Bar[drho]-d[Imrho]]

i.e., the imaginary part of the 1-form d[ρ] equals the exterior derivative of the imaginary part of the scalar ρ. Similarly:

rhor=D[ρ,r];

FacSimp[rhor-Bar[rhor]-D[Imrho,r]]

Rdddd[[1,2,1,2]]

ωUd[[1,3]]

One can also define other NP quantities, say the spin connection 1-form

SpinOmega01=FacSimp[γ e[1]+ϵ e[2]-α e[3]-β e[4]]

5.2  Implicit Frame Calculations

For some problems, giving the explicit form of the coframe (in terms of the coordinate differentials) is not the most efficient way of doing the calculations. This is particularly true if the exterior derivatives of the coframe satisfy simple expressions, not involving the coordinates, as in the Bianchi cosmological models, for which d[ea]= -(1/2)(Subscript[(C^a), bc])eb⋀ec  with (Subscript[(C^a), bc])=constants.

To do the calculations in such cases, RGtensors makes use of two global variables - deList and dxRuleList - which must be defined by the user. The first (deList) must be set equal to the expressions for d[ea] (in terms of  e[a]⋀e[b]),  and the second (dxRuleList) must be set equal to a list of rules for replacing the coordinate differentials in terms of the e[a] basis. Actually, only the differentials of  those coordinates that appear in the metric and in the expressions for d[e[a]] need to be included in dxRuleList. Finally, one must call RGtensors with {e[]} as third argument: RGtensors[gIN,xIN,{e[]}]. The coframe thus remains undefined in terms of the coordinate differentials, but its exterior derivative is known.

For the last example (4c), the lists dxRuleList and deList can be obtained as follows:

Clear$dx

varList={(Sqrt[(a^2)-2 M r+(r^2)])};

NPframe

d[xCoord]

dxRuleList=FacSimp[Solve[NPframe==e/@{1,2,3,4},d[xCoord]][[1]]]

deList=FacSimp[d[NPframe]/.dxRuleList]

Now RGtensors[gNP,xCoord,{e[]},{0,0}] gives the same results faster, except for the line element, which is just 2 e[1] e[2]-2 e[3] e[4].

RGtensors[gNP,xCoord,{e[]},{0,0}]

κ

Clear$dx

d[a]=.;d[M]=.;varList={};

Warning: When the differentials of only some of the coordinates are included in dxRuleList, covD and covDiv behave as if differentiations with respect to the remaining coordinates give zero.

◼   Example 9: Diagonal, Class-A, Bianchi Cosmological Models

Clear$dx

xCoord={x1,x2,x3,t};

g={{-A[t],0,0,0},{0,-B[t],0,0},{0,0,-C[t],0},{0,0,0,F[t]}};

The components of g refer to an undefined coframe {e[1],e[2],....,e[Dim]}.

dxRuleList={d[t]->e[4]};

deList={n[1]e[2]⋀e[3],n[2]e[3]⋀e[1],n[3]e[1]⋀e[2],0}

The elements of deList are, by definition, {d[e[1]],d[e[2]],....,d[e[Dim]]}.

(** d[n[_]]=0; **)

d[n[_]] is left undefined intentionally to generate a warning message.

simpRules={};

RGtensors[g,xCoord,{e[]}]

Clear$dx

Clear$dx clears deList and dxRuleList also, so they must be redefined 

	Choosing n[i_] and A[t], B[t], C[t], F[t] appropriately, particular models can be selected:

(** Taub (Bianchi type IX) **)dxRuleList={d[t]->e[4]};deList={n[1]e[2]⋀e[3],n[2]e[3]⋀e[1],n[3]e[1]⋀e[2],0}/.n[_]->-1;

RGtensors[g/.{A[t]->1+t^2,B[t]->1+t^2,C[t]->4(1+2a t-t^2)/(1+t^2),F[t]->(1+t^2)/(1+2a t-t^2)},xCoord,{e[]}]

zeroQ[Bianchi[0]]//Timing

Clear$dx

(** Kasner (Bianchi type I) - u is a constant parameter **)d[u]=0;dxRuleList={d[t]->e[4]};
deList={n[1]e[2]⋀e[3],n[2]e[3]⋀e[1],n[3]e[1]⋀e[2],0}/.n[_]->0;

RGtensors[g/.{F[t]->1,A[t]->t^(-2u/(1+u+u^2)),B[t]->t^(2(1+u)/(1+u+u^2)),C[t]->t^(2u(1+u)/(1+u+u^2))},xCoord,{e[]}]

zeroQ[Bianchi[1]]

Clear$dx;Clear[n,g];d[u]=.;

See also Example 12 in Section 6.

6  Approximations - Series Expansions

The built-in Mathematica function Series[f,{x,x0,n}], which gives the series expansion of the function f around x0 up to order n, can be applied to lists, and thus to tensors. RGtensors and the other functions defined in this package (version 2.3 or later) can handle SeriesData objects as input data -- gIN (metric) and eIN (coframe) --  and thus allow approximate calculations (perturbations of a given solution, asymptotic expansion of the metric or the coframe, etc.).

	There are several points, however, that must be kept in mind in doing calculations with series:

6.1  Be aware of how your version of Mathematica handles series

There are significant differences in the way different versions of Mathematica handle series. Under Mathematica 3.0, algebraically equivalent expressions do not have identical series approximations at infinity (this problem does not arise in version 4.0 or higher):

Series[1-2m/r+a^2/r^2,{r,#,2}]&/@{0,Infinity}

Series[(r-2m+a^2/r)/r,{r,#,2}]&/@{0,Infinity}

Series[(r^2-2m r+a^2)/r^2,{r,#,2}]&/@{0,Infinity}

As a consequence, with Mathematica 3.0, calculations with such expansions may not have the correct accuracy -- see Example 12 below.

RGTC corrects some of the notorious bugs of SeriesData objects  in Mathematica 3 and 4. For example, before evaluating EDCRGTCcode, multiplication of a series by zero gives a non-zero result

0*Series[1+a x+b x^2,{x,0,2}]

Mathematica versions 5 or higher correct this and other problems, but cannot handle substitutions of series expressions in the coefficients of a series, a situation that arises in local frame transformations. For example, given a series

ser1=Series[Sum[c[k]x^k,{k,4}],{x,0,4}]

suppose one needs to substitute new series expansions for the coefficients:

serRules={c[k_]:>Series[Sum[b[m,k]/x^m,{m,-2,2}],{x,0,4}]}

c[1]/.serRules

Mathematica leaves such expressions unevaluated:

Simplify[ser1/.serRules]

RGTC enables Series to handle such expressions. Thus, after evaluating EDCRGTCcode, one gets:

Simplify[ser1/.serRules]

Also Mathematica versions 5 or higher are considerably slower in handling series expressions. In some cases, they can be more than 10 times slower compared to Mathematica 4 (substitution of series expansions in complicated expressions -- see alternate ways to verify the Tomimatsu-Sato solutions,  TS123AltEV.nb,  in AdditionalExamples.html.). The time taken to complete all calculations in this section differ significantly between versions:

({{MathematicaVersion, 3, 4, 5, 6, 7, 8}, {Time(sec), 11, 13, 21, 18, 18, 21}});On[General::"spell1"];

The slower times in later versions is probably due to the fact that newer Mathematica versions keep terms to higher order than Mathematica 3 or 4. For the examples in the next two subsections, the results (Rdd, LineElement) are identical between versions 3 and 4, and between all versions 5 or higher. However, the version 5 or higher results are given to higher order (at least for some tensor components) than those of version 4. Thus, for a given required accuracy in the result (Rdd), lower order expansions of the metric are needed in later versions.

6.2  Series expansion of the metric and /or coframe

Let us consider a metric depending on a parameter b, which reduces to Schwarzschild when b=0:

xCoord={t,r,θ,φ};simpRules=TrigRules;

A[r]=1-2m  r/(r^2-b^2);B[r]=Exp[2b/r];

gPERT=DiagonalMatrix[{A[r],-1/A[r],-r^2/B[r],-r^2B[r]Sin[θ]^2}]

We can use the Mathematica function Series to find its series expansion around b=0 to any order:

gSERb0[n_]:=FacSimp[Series[gPERT,{b,0,n}]]

gSERb0[3]

Using this series expansion for the metric in RGtensors results in all tensor components being computed as series expansions to the appropriate order:

◼   Example 10a: Expansion of a "perturbed Schwarzschild" metric around b=0

RGtensors[gSERb0[4],xCoord]

Rdd

%/.b->0

The parameter b has been set to a constant:

d[b]

One can also expand the metric with respect to one of the coordinates to examine the geometry near a hypersurface. In this case the notation d[x] for the differential of that coordinate cannot be used because Series will consider d[x] a function of x and formally expand it (using d[0]=0):

Series[d[x]/(1-x),{x,0,2}]

When such an expansion of the metric w.r.t. a coordinate is used in the default coordinate frame, RGtensors will automatically introduce the notation dx for d[x]. The symbol dx will then be treated correctly as a closed 1-form. The user can still input d[x].

◼   Example 10b: Expansion of a "perturbed Schwarzschild" metric as r→∞

gSERrInf[n_]:=FacSimp[Series[gPERT,{r,Infinity,n}]]

gSERrInf[3]

RGtensors introduces the symbol dr  for d[r]. Note its use in LineElement:

RGtensors[gSERrInf[4],xCoord]

Rdd

When b=0 the metric is Ricci-flat:

RGtensors[gSERrInf[4]/.b->0,xCoord]

The "order n" refers to the components of the mixed tensor RUd:

Rdd

gUU.Rdd

Of course, one can use any coframe to do this calculation. In that case the coframe must be expanded in a series also. But since the coframe depends on the differentials of the coordinates, when expanding with respect to the coordinate x, the user must introduce the symbol dx  for d[x] before expanding the coframe. This is achieved with the following commands (where x is the symbol for the expansion coordinate) :

DeclareForms[1,dx];d[x]=dx;d[dx]=0;

Then Series does not expand d[x] (= dx):

Series[d[x]/(1-x),{x,0,2}]

And dx is treated correctly as a closed 1-form:

d[%]

%%⋀%%

◼   Example 11: Example 10b in an orthonormal frame

Clear$dx

d[m]=0;d[a]=0;d[b]=0;

ONframe=PowerExpand[{Sqrt[A[r]]d[t],d[r]/Sqrt[A[r]],r d[θ]/Sqrt[B[r]],r Sin[θ]Sqrt[B[r]]d[φ]}];

gON=DiagonalMatrix[{1,-1,-1,-1}];

To expand with respect to the coordinate r, the notation dr  for d[r] must be used:

DeclareForms[1,dr];d[r]=dr;d[dr]=0;

ONfrSerInf[n_]:=FacSimp[Series[ONframe,{r,Infinity,n}]]

ONfrSerInf[4]

RGtensors[gON,xCoord,ONfrSerInf[4]]

Rdd

%/.b->0

After a frame calculation, the coordinate differentials are expressed in terms of the coframe:

d[r]

If repeated calculations (e.g., for different orders) are to be performed, the differentials of the coordinates (and of the symbol dx) must be cleared before each new calculation.This is accomplished with the command Clear$dx  - see section 5.1.1 "Note on repeated frame calculations".

Clear$dx;d[r]

(** Order 5 **)RGtensors[gON,xCoord,ONfrSerInf[5]]

Rdd

Remember that the coordinate differentials have been set equal to their expressions in terms of the e[i] everywhere. To recover the original expressions use eTO$dx (see section 5.1.1):

LineElement

FacSimp[LineElement]

FacSimp[LineElement/.eTO$dx]

Coefficient[%,𝒹[r]^2]

The function GenCoef[x_,y_] ("GeneralizedCoefficient") returns a series for the coefficient of y in x when x is a series:

GenCoef[%%,𝒹[r]^2]

It is also possible to do the calculation with both the metric and the coframe expanded in series:

◼   Example 11a: Examples 10a,b in a "special" frame

Clear$dx

SPECframe=PowerExpand[{d[t],d[r],r d[θ]/Sqrt[B[r]],r Sin[θ]Sqrt[B[r]]d[φ]}];

gSPEC=DiagonalMatrix[{A[r],-1/A[r],-1,-1}];

(i) Expand both frame and metric around b=0:

SPECframeb0[n_]:=FacSimp[Series[SPECframe,{b,0,n}]]

gSPECb0[n_]:=FacSimp[Series[gSPEC,{b,0,n}]]

SPECframeb0[4]

gSPECb0[4]

RGtensors[gSPECb0[4],xCoord,SPECframeb0[4]]

Rdd

%/.b->0

(ii) Expand both frame and metric for r→∞:

Again the notation dr  for d[r] must be introduced before expanding the coframe w.r.t. r:

Clear$dx

DeclareForms[1,dr];d[r]=dr;d[dr]=0;

SPECframeINF[n_]:=FacSimp[Series[SPECframe,{r,Infinity,n}]]

gSPECinf[n_]:=FacSimp[Series[gSPEC,{r,Infinity,n}]]

SPECframeINF[4]

gSPECinf[4]

RGtensors[gSPECinf[4],xCoord,SPECframeINF[4]]

Rdd

%/.b->0

6.3  Expanding deList and dxRuleList in implicit frame calculations

If an implicit frame is used (see Section 5.2), then, in addition to the metric, the expressions in deList and dxRuleList must be expanded in a series, especially if they are complicated. The reason is that before a normal expression can be multiplied by a series, it must be expanded itself in a series. Thus, if deList and dxRuleList are given as normal expressions, Mathematica will keep expanding these expressions again and again every time it performs a differentiation, causing the calculations to proceed very slowly.

◼   Example 12: Example 10b in an implicit NP null frame

Clear$dx

NPframe=FacSimp[{(ONframe[[1]]-ONframe[[2]])/(Sqrt[2]),(ONframe[[1]]+ONframe[[2]])/(Sqrt[2]),(ONframe[[3]]+I ONframe[[4]])/(Sqrt[2]),(ONframe[[3]]-I ONframe[[4]])/(Sqrt[2])}]

NPframe=NPframe/.Power[1-((2 m r)/(-(b^2)+(r^2))),n_]:> Power[-(b^2)-2 m r+(r^2),n]Power[-(b^2)+(r^2),-n]

dxTOe=FacSimp[Solve[NPframe=={e[1],e[2],e[3],e[4]},{d[t],d[r],d[θ],d[φ]}][[1]]]

dNPframeTOee=FacSimp[d[NPframe]/.dxTOe]

For implicit frame calculations, the user must define the global variables deList and dxRuleList. For this problem, only the θ, r differentials need to be included in dxRuleList :

deList=dNPframeTOee;

dxTOe2=Drop[dxTOe,{2,3}];

dxRuleList=dxTOe2

One should check that deList is a set of closed 2-forms, giving the d[e[i]]:

checkClosure=FacSimp[d[deList]/.d[e[i_]]:>deList[[i]]/.dxRuleList]//Timing

(1)  Exact calculation

RGtensors[{{0,1,0,0},{1,0,0,0},{0,0,0,-1},{0,0,-1,0}},xCoord,{e[]}]

ExactRdd=Rdd

{Φ[0,0],Φ[0,2],Φ[1,2]}

ExactPSI=Ψ/@{0,1,2,3,4}

(1)  Approximate calculation

Clear$dx

After an implicit frame calculation, Clear$dx clears deList and dxRuleList also.

deList=FacSimp[Series[dNPframeTOee,{r,Infinity,3}]]

dxRuleList=Map[Rule[#[[1]],FacSimp[Series[#[[2]],{r,Infinity,3}]]]&,dxTOe2]

checkClosure=FacSimp[d[deList]/.d[e[i_]]:>deList[[i]]/.dxRuleList]//Timing

GenCoef[x_,y_] returns the correct series coefficient of y in x, when x is a series - otherwise acts as FacSimp[Coefficient[x_,y_]]

Coefficient[deList,e[1]⋀e[3]]

GenCoef[deList,e[1]⋀e[3]]

FacSimp[%%]

GenCoef[Normal[deList],e[1]⋀e[3]]

RGtensors[{{0,1,0,0},{1,0,0,0},{0,0,0,-1},{0,0,-1,0}},xCoord,{e[]}]

Clear$dx

ExactRdd-Rdd

The function zeroQ[] ignores ((O[(1/r)])^n) terms:

zeroQ[%]

ExactPSI-Ψ/@{0,1,2,3,4}

Mathematica 4.0 (or higher) has better handling of SeriesData objects. As a result it returns ((O[(1/r)])^5) for all elements in the last two lists.

If only deList is expanded, RGtensors takes ~20% longer, while if only dxRuleList is expanded, it takes ~100% longer!

7  Classification of Weyl and Ricci Tensors

Beginning with RGTC version 3.6.5, the new function Classify[x_] carries out the Petrov-Penrose classification when x is a 4th-rank tensor with the symmetries of the Weyl tensor (Wdddd), and the Plebanski-Segre classification when x is a 2nd-rank symmetric tensor (Rdd). Classify is defined for 4-dimensional spaces only. Classify tests that the given tensor has the correct symmetry for its rank before classifying it. The tensors must be given with all indices covariant (down).

Classify[Wdddd] prints the Petrov-Penrose type (I, II, IID, III, N, O) and returns a list containing one of the numbers (1111, 211, 22, 31, 4, 0) describing the alignment of principal null directions.

Classify[Rdd] prints the correponding Segre type (in the notation of the Exact Solutions Book) and additional information if Rdd is traceless, has zero as an eigenvalue, or corresponds to a physical energy momentum tensor. As the complete Plebanski-Segre classification requires, in some cases, deciding whether an expression is positive or negative (real or complex eigenvalues, timelike or spacelike character of eigenvectors), it cannot be carried out completely on a computer. Thus Classify[Rdd] will print multiple Segre types in such cases.

Classify[Rdd] returns a symbol indicating the highest root multiplicity and the degree and structure of the minimal polynomial. The first part of the symbol is a roman numeral (as in the Petrov-Penrose case) indicating the root multiplicity. Next is the letter M - minimal - followed by an integer, n, indicating the degree of the minimal polynomial (except when the minimal polynomial coincides with the characteristic polynomial in which case the letter G - general - is used). Finally, in parenthesis, there are n digits (0's or 1's), which correspond to the coefficients (Subscript[c, k]) of different powers in the minimal polynomial: (x^n)-(Subscript[c, n-1])(x^(n-1))-(Subscript[c, n-2])(x^(n-2))-...-(Subscript[c, 0]) (x^0)=0. A zero digit signifies that the corresponding coefficient vanishes. The  following examples explain the notation:

The symbol II-M3(010) means there is one double root (II), the minimal polynomial is cubic (M3) and has the form (x^3)-(Subscript[c, 1])x=0, i.e., the coefficients (Subscript[c, 2]) and (Subscript[c, 0]) vanish (010).

The symbol IID-M3(011) means there are two double roots (IID), the minimal polynomial is cubic (M3) and has the form (x^3)-(Subscript[c, 1])x-(Subscript[c, 0])=0, i.e., the coefficient (Subscript[c, 2]) vanishes (011).

The symbol II-G(0100) means there is one double root (II), the minimal polynomial is quartic (G) and has the form (x^4)-(Subscript[c, 2])(x^2)=0, i.e., the coefficients (Subscript[c, 3, ])(Subscript[c, 1 ])and (Subscript[c, 0]) vanish (0100).

The symbol IV-M2(01) means there is a quadruple root (IV), the minimal polynomial is quadratic (M2) and has the form (x^2)-(Subscript[c, 0])=0, i.e., the coefficient (Subscript[c, 1]) vanishes (01).

The different combinations of non-vanishing coefficients (for the same root multiplicity and polynomial degree) define invariant subclasses for each Segre type: they imply algebraic relations between the eigenvalues which do not change the multiplicity. Note that, unlike the root multiplicity and the degree of the minimal polynomial, the digits in parenthesis can change when a term proportional to the metric is added to Rdd. Thus, to classify the energy-momentum tensor, one must use Classify[Rdd-1/2R gdd] or Classify[gdd.EUd]. And, of course, the minimal polynomial for a tensor Tdd translates to a matrix identity relating different powers of the 4x4 matrix TUd=gUU.Tdd.

Another new (in version 3.6.5) function is Plebanski[x_], which constructs the Plebanski tensor (4th-rank covariant tensor quadratic in x with the symmetries and trace properties of the Weyl tensor) associated to the symmetric covariant 2nd-rank tensor x. Plebanski is defined for any space dimension. Plebanski is also defined with two symmetric tensors as arguments: Plebanski[x_,y_] is the unique 4th-rank covariant tensor bilinear and symmetric in x and y with the symmetries and trace properties of the Weyl tensor. The definition of P(X,Y) in terms of components is:

(Subscript[P, abcd])(X,Y)=(Subscript[X, ac])(Subscript[Y, bd])-(Subscript[X, bc])(Subscript[Y, ad])-(Subscript[X, ad])(Subscript[Y, bc])+(Subscript[X, bd])(Subscript[Y, ac])+(1/(n-2))((Subscript[g, ac])(Subscript[Z, bd])-(Subscript[g, bc])(Subscript[Z, ad])-(Subscript[g, ad])(Subscript[Z, bc])+(Subscript[g, bd])(Subscript[Z, ac]))+(2/((n-1)(n-2)))(X Y-(Subscript[(X^e), f])(Subscript[(Y^f), e]))((Subscript[g, ac])(Subscript[g, bd])-(Subscript[g, bc])(Subscript[g, ad])),
where (Subscript[Z, ab])=(Subscript[X, af])(Subscript[(Y^f), b])+(Subscript[Y, af])(Subscript[(X^f), b])-X (Subscript[Y, ab])-Y (Subscript[X, ab]), X=(Subscript[(X^a), a]), Y=(Subscript[(Y^a), a]),. When x=y  Plebanski satisfies the identity Plebanski[x,x]=2*Plebanski[x], while when x (or y) is the metric tensor Plebanski[x,g] vanishes.

In 4 dimensions, if y is the square of the traceless part of x, then Plebanski[x,y] is identically zero as follows from the 3rd order tensor identity discussed in this artticle. For example:

Dim=4;gdd=DiagonalMatrix[{-1,1,1,1}];gUU=gdd;

s[i_,j_]:=s[j,i]/;i>j;SS=Array[s,{4,4}]

SS0=SS-tr[SS.gUU]/4gdd;

zeroQ[Plebanski[SS,SS0.gUU.SS0]]

Both Classify and Plebanski use the metric gdd, gUU (and eta[..] for the Weyl classification) defined by RGtensors to raise indices, construct duals or form tensor products. Thus, these functions must be defined separately if one wishes to classify a particular Weyl- or Ricci-like tensor without calling RGtensors (see the preliminary definitions in the "Canonical Forms" section of the examples).

In the examples below, use is made of a new auxiliary function defined in version 3.6.5 (see Section 4.2), metric[ds2_, xIN_], which constructs the metric tensor from the line element (ds2) and the list of coordinates (xIN). All calculations are carried out in the coordinate frame, as the classification algorithm does not depend on relations between components in a particular frame (e.g., NP frame). Use of an appropriate frame in each case would speed up the calculations considerably!

Examples from the Exact Solutions Book (Second Edition)

The examples are referred to by their equation number in the book: H Stephani, D. Kramer, M. A. H. MacCallum, C. Hoenselaers, E. Herlt, Exact Solutions of Einstein's Field Equations, Second Edition, Cambridge University Press, 2003.

7.1  Canonical Forms of the Ricci Tensor (Eq. 5.3a - 5.3d)

Preliminary definitions

To classify a given tensor we must first define the global variables gdd, gUU and Dim (the space dimension). As the canonical forms are given with respect to unit vectors in an orthonormal frame, we define

Dim=4;gdd=DiagonalMatrix[{1,1,1,-1}]

gUU=gdd;

We also need the definition of eta[..], to be able to classify the Weyl tensor.

eta[]:=eta[]=Array[Signature[{##}]&,Table[4,{4}]];eta[a__]:=eta[a]=Raise[eta[],a];

To construct the canonical forms, we need to define unit (or special null) vectors along the frame axes and the operation of symmetric tensor product of 2 vectors (SymProd[x_, y_]).

Xd={1,0,0,0};Yd={0,1,0,0};Zd={0,0,1,0};Ud={0,0,0,1};

Kd=(Ud-Zd)/Sqrt[2];Ld=(Ud+Zd)/Sqrt[2];

SymProd[x_,y_]:=(Outer[Times,x,y]+Outer[Times,y,x])/2

Clear[A1,A2,A3,B,λ];

A1 class (5.3a)

A1=λ[1]*SymProd[Xd,Xd]+λ[2]*SymProd[Yd,Yd]+λ[3]*SymProd[Zd,Zd]-λ[4]*SymProd[Ud,Ud]

Classify[A1]

Classify[A1/.λ[1]->λ[2]]

Classify[A1/.λ[1]->λ[2]/.λ[3]->λ[2]]

Classify[A1/.λ[1]->λ[2]/.λ[3]->λ[4]]

Classify[A1/.λ[1]->λ[2]/.λ[4]->-λ[2]/.λ[3]->-λ[2]]

Classify[A1/.λ[_]->λ[2]]

Classify[A1/.λ[_]->0]

A2 class (5.3b)

A2=λ[1]*SymProd[Xd,Xd]+λ[2]*SymProd[Yd,Yd]+λ[3]*SymProd[Kd,Ld]+λ[4]*(SymProd[Kd,Kd]-SymProd[Ld,Ld])

Classify[A2]

Classify[A2/.λ[1]->λ[2]]

Classify[A2/.λ[1]->λ[2]/.λ[3]->λ[2]]

Classify[A2/.λ[1]->0/.λ[2]->0/.λ[3]->0]

A3 class (5.3c)

A3=λ[1]*SymProd[Xd,Xd]+λ[2]*SymProd[Yd,Yd]-2λ[3]*SymProd[Kd,Ld]+(SymProd[Kd,Kd])

Classify[A3]

Classify[A3/.λ[1]->λ[2]]

Classify[A3/.λ[3]->λ[2]]

Classify[A3/.λ[3]->λ[1]/.λ[2]->λ[1]]

tr[A3.gUU]

Classify[A3/.λ[3]->-(λ[1]+λ[2])/2]

Classify[A3/.λ[_]->0]

B class (5.3d)

B=λ[1]*SymProd[Xd,Xd]+λ[2]*SymProd[Yd,Yd]+λ[3]*SymProd[Kd,Ld]+(SymProd[Kd,Xd])

Classify[B]

Classify[B/.λ[3]->0]

Classify[B/.λ[1]->λ[2]]

Classify[B/.λ[_]->λ[1]]

tr[B.gUU]

Classify[B/.λ[3]->λ[1]+λ[2]]

Classify[B/.λ[1]->-λ[2]/.λ[3]->0]

Classify[B/.λ[1]->0/.λ[2]->0/.λ[3]->1]

Classify[B/.λ[1]->λ[2]/.λ[3]->0]

checkMatrixIdentity=Factor[MM.MM.MM-c[2]MM.MM/.MM->B.gUU/.λ[1]->λ[2]/.λ[3]->0]

%/.c[2]->λ[2]

Classify[B/.λ[_]->0]

7.2  Perfect Fluids

Petrov-Type I (Eq. 21.59 -- Corrected)

The metric is incorrect as printed in the book. The correct form is given here.

ds2=(d[t]-2b y d[ph]/a)^2/x-x(1-al)^2Exp[-y](d[x]^2/(a x^2-b^2)+d[y]^2/a/(y+al(Exp[-y]-1)))-4x(y+al(Exp[-y]-1))d[ph]^2/a

xIN={t,x,y,ph}

gIN=metric[ds2,xIN]

RGtensors[gIN,xIN]

EUd

Classify[gdd.EUd]

Classify[Wdddd](** this takes long! **)

Petrov-Type II (Eq. 33.14)

xIN={u,w,r,ph}

ds2=(1-r^2)/(2ec)(d[w]^2/(4w)+w Exp[w] d[ph]^2)-2d[u]d[r]-ec(r^2+2w+3)d[u]^2

gIN=metric[ds2,xIN]

RGtensors[gIN,xIN]

Classify[Wdddd]

Classify[gdd.EUd]

Classify[Plebanski[Rdd]]

Petrov-Type D (Eq. 21.61 -- sign correction)

Clear[MM,NN,G]

ds2=1/NN^2(d[z]^2+d[x]^2/G+G*d[ϕ]^2-x^2(d[t]+nn/x^2*d[ϕ])^2)

NN=p*Exp[z]+q*Exp[-z];G=-x^2+b*Log[x]+cc+nn^2/x^2;(** Note minus sign! **)

xIN={t,x,z,ϕ}

gIN=FacSimp[metric[ds2,xIN]]

RGtensors[gIN,xIN]

Classify[Wdddd]

Classify[gdd.EUd]

Petrov-Type III (Eq. 33.45)

xIN={u,t,x,y}

ds2=2 Exp[3t+4x]d[x]^2+(3Exp[t+2x-2y]+Exp[3t+6x]/a)d[u]^2+a*Exp[t]/(Exp[2y]-1)d[y]^2-2Exp[2t+3x-y]d[u]d[t+3x+y]

gIN=metric[ds2,xIN]

RGtensors[gIN,xIN]

Classify[Wdddd]

Note: Inadequate simplification of the exponentials in some versions of Mathematica, may prevent Classify from verifying that Wdddd and gdd.EUd have the correct symmetries!

Classify[gdd.EUd]

Petrov-Type N (Eq. 33.48 - Class I)

xIN={u,t,x,y}

rtt=Sqrt[t]

Special simplificaion rules, implied by the differential equation satisfied by the metric functions, must be defined.

ghRules=Join[FuncRepRules[D[g1[x,u],x,x],-a[u]^2g1[x,u]],FuncRepRules[D[h1[y,u],y,y],-b^2h1[y,u]]];

G1rules=FacSimp[FuncRepRules[G1[x,y,u],g1[x,u]h1[y,u],3]/.ghRules];

simpRules=Join[G1rules,ghRules];

H1=-2(a[u]^2rtt+b^2rtt^3)

ds2=Collect[2G1[x,y,u] d[t]d[u]+2G1[x,y,u]^2H1 d[u]^2-rtt^3(d[x]-2/rtt*D[G1[x,y,u],x]d[u])^2-rtt(d[y]+2rtt*D[G1[x,y,u],y]d[u])^2,{d[_]^2,d[i_]d[j_]},FacSimp]

gIN=FacSimp[metric[ds2,xIN]]

RGtensors[gIN,xIN]

Classify[Wdddd]

Classify[Rdd]

simpRules={};

7.3  Electromagnetic Fields

Kerr-Newman (Eq. 21.24)

xIN={t,r,x,ϕ} (* x = Cos[θ] *)

ds2=(r^2+a^2x^2)(d[r]^2/(r^2-2M*r+a^2+e^2)+d[x]^2/(1-x^2))-(1-(2M*r-e^2)/(r^2+a^2x^2))d[t]^2+(1-x^2)(r^2+a^2+a^2(1-x^2)(2M*r-e^2)/(r^2+a^2x^2))d[ϕ]^2-2a(1-x^2)(2M*r-e^2)/(r^2+a^2x^2)d[ϕ]d[t]

gIN=metric[ds2,xIN]

RGtensors[gIN,xIN]

Classify[Wdddd]

Classify[gdd.EUd]

Classify[Plebanski[Rdd]]

Type D with G2 on Null Orbits (Eq. 24.21)

xIN={u,v,x,y}

ds2=(x^2+y^2)/FF*d[x]^2+FF/(x^2+y^2)(d[u]-y^2d[v])^2+2d[y](d[u]+x^2d[v])

FF=2n*x-(e^2+g^2);

gIN=metric[ds2,xIN]

RGtensors[gIN,xIN]

Classify[Wdddd]

Classify[Rdd]

Classify[Plebanski[Rdd]]

McLenaghan-Tariq-Tupper (Eq. 12.21)

xIN={t,x,y,ϕ}

ds2=a^2/x^2(d[x]^2+d[y]^2)+x^2d[ϕ]^2-(d[t]-2y*d[ϕ])^2

gIN=metric[ds2,xIN]

RGtensors[gIN,xIN]

Classify[Wdddd]

Classify[Rdd]

Classify[Plebanski[Rdd]]

PP wave (Eq. 24.40 -- real coordinates)

xIN={u,v,x,y}

ds2=d[x]^2+d[y]^2-2d[u]d[v]-2HH[u,x,y]d[u]^2

gIN=metric[ds2,xIN]

RGtensors[gIN,xIN]

Classify[Wdddd]

Classify[Rdd]

Conformally Flat wave (Eq. 24.51)

xIN={u,v,x,y}

ds2=(1-Sin[ω*u])d[x]^2+(1+Sin[ω*u])d[y]^2-2d[u]d[v]

gIN=metric[ds2,xIN]

simpRules=TrigRules;

RGtensors[gIN,xIN]

Classify[Wdddd]

Classify[Rdd]

simpRules={};

7.4  Miscellaneous test metrics

Curzon (Eq. 20.5 @ 20.3)

xIN={t,r,θ,ϕ}

ds2=Exp[2M/r](Exp[-M^2/r^2Sin[θ]^2](d[r]^2+r^2d[θ]^2)+r^2Sin[θ]^2d[ϕ]^2)-Exp[-2M/r]d[t]^2

gIN=metric[ds2,xIN]

simpRules=TrigRules;

RGtensors[gIN,xIN]

Classify[Wdddd]

Classify[Rdd]

simpRules={};

Robinson-Trautman Type II (Eq. 28.78 @ 28.26)

xIN={u,r,z,zB}

Bar[z]=zB;Bar[zB]=z;

ds2=2r^2/(z+zB)^3d[z]d[zB]-2d[u]d[r]+(3(z+zB)+2M/r+Λ/3r^2)d[u]^2

gIN=metric[ds2,xIN]

RGtensors[gIN,xIN]

Classify[Wdddd]

Classify[gdd.EUd]

Classify[Plebanski[Rdd]]

Khan-Penrose colliding waves (Eq. 25.24)

xIN={η,μ,x,y}

W=Sqrt[1-η^2]Sqrt[1-μ^2]

ds2=(1-(η^2))^(3/4) /(1-(μ^2))^(1/4)(d[μ]^2/(1-μ^2)-d[η]^2/(1-η^2))+W((1-η)/(1+η)d[x]^2+(1+η)/(1-η)d[y]^2)

gIN=Factor[metric[ds2,xIN]]

RGtensors[gIN,xIN]

Classify[Wdddd]

This message sometimes appears due to Mathematica's failure to recognize that an expression vanishes. This happens when the determinant of the metric tensor is not a perfect square, as in this case, and Mathematica writes (1-η) as -(-1+η):

detg

The expression that fails to vanish is:

Simplify[1-(((Sqrt[-1+η]) (Sqrt[1+η]) (Sqrt[1-(μ^2)]))/((Sqrt[1-(η^2)]) (Sqrt[-1+μ]) (Sqrt[1+μ])))]

In this case, a coordinate transformation solves the problem:

Factor[W/.{η->2h/(1+h^2),μ->2k/(1+k^2)}]

Simplify[%/.Power[1-((4 ((x_)^2))/(((1+((x_)^2)))^2)),Rational[m_,n_]]:>Power[((1- (x^2))/((1+(x^2)))),Rational[2m,n]]]

The line element now becomes:

ds2=Collect[PowerExpand[ds2/.{η->2h/(1+h^2),μ->2k/(1+k^2)}/.Power[1-((4 ((x_)^2))/(((1+((x_)^2)))^2)),Rational[m_,n_]]:>Power[((1- (x^2))/((1+(x^2)))),Rational[2m,n]]],_d,Factor]

xIN={h,k,x,y}

gIN=Factor[metric[ds2,xIN]]

RGtensors[gIN,xIN]

detg

Classify[Wdddd]

Warning: Even though Classify takes very little time in these examples, it can take a very long time on algebraically complicated metrics, as it evaluates expressions that are up to order 12 in the components of the input tensors. Classify[Wdddd] usually takes a few times as long as RGtensors takes to compute all tensors; but for complicated metrics it can take hundreds of times as long!

Classify[Rdd] finds the root multiplicity first and then tests the possible minimal polynomials for that multiplicity, starting with the one of lowest degree. Lastly, for the particular minimal polynomial found, it tests whether a matrix identity with some terms missing is satisfied. In some cases Mathematica 3 gives erroneous results for this last step, because it cannot solve the system of nonlinear equations obtained.

It is possible, especially with Mathematica versions <5.0, that due to inadequate simplification, Classify[Wdddd] will print the message: "Cannot perform classification; try different coordinates or frame!". This happens in the last example in section 7.4, where a change of coordinate solves the problem. For the same reason, it can happen that Classify will be unable to verify that Wdddd (or gdd.EUd) has the correct symmetries and will not proceed with the classification!