# Curvature Tensors

The defining operation for a tensor is the coordinate transformation, here effected on a tensor with one covariant (lower) index:
T'b'(x') = Ta(x') d xa(x') / d x'b'
As always, repeated indices are summed over, and the primed indices indicate the new coordinates. Note that the valence of the primed index is inverted: an upper index in the denominator is equivalent to a lower index in the numerator.

This transformation consists of two operations. The first involves a substitution of the original coordinates as functions of the new ones in the original tensor. The resulting tensor is then multiplied by the matrix whose elements are the derivatives of the old coordinates as functions of the new. For tensors with multiple covariant indices, the matrix multiplication is performed on each original index: the transformation is a multilinear operation, the tensor transforming linearly in each index.

Tofxpll = Tll /. Thread[ xu -> xfu];

dxofxpdxplu = Table[ Table[

D[ xfu[[j]], xpu[[i]]],

{j, dim}], {i, dim}];

Tpll = Table[ Table[ simpler[ Sum[ Sum[
Tofxpll[[k, l]] dxofxpdxplu[[i, k]] dxofxpdxplu[[j, l]],

{l, dim}], {k, dim}]],

{j, dim}], {i, dim}];
"Tofxpll" might be read "T of x' with two lower indices", while "dxofxpdxplu" might be read "the derivative of x as a function of x' with respect to x', whose first index is lower and whose second index is upper". The use of meaningful variable names is always encouraged as a way of making your code at least somewhat self-documenting.
coordxform::usage = "arguments are Tl...l(x), xu, xu(x'), x'u";

coordxform[ T_List, x_List, xf_List, xp_List] := Block[ {dim, dimp, rank, Tp, Tofxp, dxofxpdxp}, (

rank = TensorRank[T];

If[ rank > 5,(Print[ "Error - tensor has too many indices for this function"]; Return[];),];

dim = Length[x];

If[ Length[T] == dim,,(Print[ "Error - tensor is incorrect dimension"]; Return[];)];

If[ Length[xf] == dim,,(Print[ "Error - function list x(xp) is incorrect dimension"]; Return[];)];

dimp = Length[xp];

Tofxp = T /. Thread[ x -> xf];

dxofxpdxp = Table[ Table[

D[ xf[[jj]], xp[[ii]]],

{jj, dim}], {ii, dimp}];

If[ rank == 1,
Tp = Table[ simpler[ Sum[
dxofxpdxp[[ii, jj]] Tofxp[[jj]],

{jj, dim}]],

{ii, dimp}],];
If[ rank == 2,
Tp = Table[ Table[ simpler[ Sum[ Sum[
dxofxpdxp[[ii, kk]] dxofxpdxp[[jj, ll]] Tofxp[[kk, ll]],

{ll, dim}], {kk, dim}]],

{jj, dimp}], {ii, dimp}],];
If[ rank == 3,
Tp = Table[ Table[ Table[ simpler[ Sum[ Sum[ Sum[
dxofxpdxp[[ii, ll]] dxofxpdxp[[jj, mm]] dxofxpdxp[[kk, nn]] Tofxp[[ll, mm, nn]],

{nn, dim}], {mm, dim}], {ll, dim}]],

{kk, dimp}], {jj, dimp}], {ii, dimp}],];
If[ rank == 4,
Tp = Table[ Table[ Table[ Table[ simpler[ Sum[ Sum[ Sum[ Sum[
dxofxpdxp[[ii, mm]] dxofxpdxp[[jj, nn]] dxofxpdxp[[kk, pp]] dxofxpdxp[[ll, qq]] Tofxp[[mm, nn, pp, qq]],

{qq, dim}], {pp, dim}], {nn, dim}], {mm, dim}]],

{ll, dimp}], {kk, dimp}], {jj, dimp}], {ii, dimp}],];
If[ rank == 5,
Tp = Table[ Table[ Table[ Table[ Table[ simpler[ Sum[ Sum[ Sum[ Sum[ Sum[
dxofxpdxp[[ii, nn]] dxofxpdxp[[jj, pp]] dxofxpdxp[[kk, qq]] dxofxpdxp[[ll, vv]] dxofxpdxp[[mm, ww]] Tofxp[[nn, pp, qq, vv, ww]],

{ww, dim}], {vv, dim}], {qq, dim}], {pp, dim}], {nn, dim}]],

{mm, dimp}], {ll, dimp}], {kk, dimp}], {jj, dimp}], {ii, dimp}],];
Return[Tp];)];
The first line documents the usage of coordxform, whose arguments are the original covariant tensor, a list of the original coordinates, a list consisting of those coordinates as functions of the new and a list of the new coordinates. Note that the dimension of the new coordinate system need not be the same as that of the old. The function definition includes patterns in the argument list which restrict the use of coordxform to arguments which are lists. This is an effective control on the invocation of the function which can help to prevent errors. Block is used to force the scope of the variables in its first argument to be local to coordxform: another technique which can help make the code less brittle.

A series of If statements are used to ensure that coordxform is called with arguments that are appropriate to its definition and consistent with each other. Note that Mathematica stores a matrix as a list of lists; its Length is therefore the number of lists at the top level. These tests also contribute to making the code more resilient in the face of accidental misuse.

The substitution operation and the creation of the derivative matrix are essentially identical to those in the previous example. The remaining statements perform the number of matrix multiplications appropriate for the rank of the tensor being transformed. Note the index variables: the author has adopted the convention of using two letter index variables inside functions so that it is less likely that a symbol occurring in an argument will have the same name as an index variable.

Here is a Mathematica code example using the function coordxform define above. In this example, we transform the Minkowski metric in spherical coordinates to eliminate the trigonometric functions using the transformation

c = cos (theta).
The third argument to coordxform is derived from the second argument and the inverse transformation. TableForm is used to format the metrics conveniently:
?coordxform
arguments are Tl...l(x), xu, xu(x'), x'u
TableForm[ gll = {{1, 0, 0, 0}, {0, r^2, 0, 0}, {0, 0, r^2 Sin[theta]^2, 0}, {0, 0, 0, -1}}]
1000
0r200
00r2 Sin[theta]20
000-1
TableForm[ coordxform[ gll, {r, theta, phi, t}, {r, ArcCos[c], phi, t}, {r, c, phi, t}]]
1000
0- r2 / (-1+c)(1+c)00
00- (-1+c)(1+c) r20
000-1
This transformation is used extensively throughout the text in order to improve Mathematica execution times.
Tb = gb a Ta

Tb = gb a Ta

The student is encouraged to write functions "raise" and "lower" analogous to coordxform. They will each have two arguments: the original tensor and the index number affected, and can be tested by performing raise on either index of the metric, and lower on either index of the inverse metric. In either case, the result should be the Identity Matrix.
Γab c = ga e (db gc e + dc gb e - de gb c) / 2,
and its lower indices are symmetric:
Γac b = Γab c
We use this symmetry to speed up the Mathematica code to compute the connection components:
Gull = Table[ Table[ Table[ 0, {kk, dim}], {jj, dim}], {ii, dim}];

Do[ Do[ Do[

Gull[[ii, jj, kk]] = simpler[ Sum[
guu[[ii, ll]] (D[ gll[[kk, ll]], xu[[jj]]] +

D[ gll[[jj, ll]], xu[[kk]]] - D[ gll[[jj, kk]], xu[[ll]]]),

{ll, dim}] / 2],

{kk, jj, dim}], {jj, dim}], {ii, dim}];
Do[ Do[ Do[
Gull[[ii, kk, jj]] = Gull[[ii, jj, kk]],

{kk, jj + 1, dim}], {jj, dim - 1}], {ii, dim}];

Note that we must first create a zero tensor with the correct number of indices; since the remainder of the code changes individual components, Mathematica requires that they exist before being modified. This code assumes that the metric associated with this connection is stored in "gll", that the inverse metric is "guu" and that xu is again the list of coordinates. We will assume that tensor names, once adopted, are used consistently in all code examples unless otherwise stated.
d2 xa / dτ2 + Γab c (d xb / dτ) (d xc / dτ) = 0
where τ is the affine parameter. Obviously, the nonzero components control the geodesic deviation. Therefore if for some b, Γab b is only nonzero for a = b, then all curves tangent to xb are geodesics. In physical terms, their acceleration is parallel to their velocity. We will call such directions geodesic directions. If b and c are members of a subset of the coordinate directions such that Γab c is only nonzero when a is also a member of this subset, then these directions can be said to define a geodesic hypersurface: any geodesic tangent to the hypersurface is wholly contained within it.
Ra b ce = db Γea c - da Γeb c + Γfa c Γef b - Γfb c Γef a
and has the following symmetries:
Ra b c d = - Rb a c d

Ra b c d = - Ra b d c

Ra b c d = Rc d a b

R[a b c] d = 0

NB: The definition of the Riemann Curvature Tensor is independent of metric signature. However, lowering the final index will induce a relative negative sign if using negative signature metrics. This effect occurs whenever the metric is involved (ie., the cosmological constant term in Einstein's Equations, the pressure term in perfect fluid Stress-Energy Tensors, etc.).
If we consider the indices in pairs, in D dimensions each antisymmetric pair has
D (D - 1) / 2
independent components. The pairs are symmetric, so that the first three symmetries imply that there are
(D (D - 1) / 2) (D (D - 1) / 2 + 1) / 2
independent components. But the last symmetry means that we need to subtract
D ! / (4 ! (D - 4) !)
components, leaving us with
D2 (D2 - 1) / 12
independent components. This implies that all D = 1 manifolds have zero curvature.

The following Mathematica code to compute the Riemann Curvature components makes use of the first three symmetries:

Rllll = Table[ Table[ Table[ Table[ 0, {ll, dim}], {kk, dim}], {jj, dim}], {ii, dim}];

Do[ Do[ Do[ Do[

Rllll[[ii, jj, kk, mm]] = simpler[ Sum[
(D[ Gull[[ll, ii, kk]], xu[[jj]]] - D[ Gull[[ll, jj, kk]], xu[[ii]]] + Sum[
Gull[[nn, ii, kk]] Gull[[ll, nn, jj]] - Gull[[nn, jj, kk]] Gull[[ll, nn, ii]],

{nn, dim}]) gll[[ll,mm]],

{ll, dim}]],
{mm, kk + 1, dim}], {kk, dim - 1}], {jj, ii + 1, dim}], {ii, dim - 1}];
Do[ Do[ Do[ Do[
Rllll[[jj, ii, kk, ll]] = - Rllll[[ii, jj, kk, ll]],

{ll, kk + 1, dim}], {kk, dim - 1}], {jj, ii + 1, dim}], {ii, dim - 1}];

Do[ Do[ Do[ Do[
Rllll[[ii, jj, ll, kk]] = - Rllll[[ii, jj, kk, ll]],

{ll, kk + 1, dim}], {kk, dim - 1}], {jj, ii + 1, dim}], {ii, dim - 1}];

Do[ Do[ Do[ Do[
Rllll[[jj, ii, ll, kk]] = Rllll[[ii, jj, kk, ll]],

{ll, kk + 1, dim}], {kk, dim - 1}], {jj, ii + 1, dim}], {ii, dim - 1}];

computing the upper half-matrix for each pair of indicies and use the symmetries to fill in the remaining components. The final symmetry is more trouble to program than it is worth.

We note that if Ra b a b is zero, then the a-b hypersurface is flat, but none of the manifolds we will be examining are that trivial.

Da Rb c e f = da Rb c e f - Γha b Rh c e f -
Γha c Rb h e f - Γha e Rb c h f - Γha f Rb c e h

= Rb c e f ;a,

the Ricci Tensor:
Ra b = gc d Ra c b d
and the scalar curvature:
R = ga b Ra b
These tensors will be used repeatedly in the analyses below, and it makes sense for the programmer to combine all of them into a single function so that for any metric, one function call will compute them all. The inverse metric and Christoffel Symbols should of course be included in the beginning of that function. It will also be useful for the student to write a function analogous to coordxform which will return the covariant derivative of a tensor with an arbitary number of covariant indices.

The following decomposition of the Riemann Curvature [Wald] is sometimes useful:

Ra b c d = Ca b c d + Ea b c d + Ga b c d
Ea b c d = (ga c Rb d - gb c Ra d -
ga d Rb c + gb d Ra c) / (D - 2)

Ga b c d = (ga d gb c - ga c gb d) R / (D - 1)(D - 2)

If we compute E (the Einstein Curvature Tensor) and G analogously to the Riemann Tensor, the following Mathematica code will compute the Weyl Tensor in a single statement:
Cllll = Rllll - Ellll - Gllll;
We will of course be interested in Einstein's Equations:
Ra b - R ga b / 2 + Λ ga b = α Ta b
α = 2 Area (SD - 2)
where the left hand side is called the Einstein Tensor, and the stress-energy tensor T has dimensions of energy density.
Provided that ga b is dimensionless, α can be multiplied by G / c4 to normalize the dimensions of all terms to 1/length2. However, if ga b is not dimensionless, as for instance in spherically symmetric metrics, it is easier to work with Einstein's Equations with one index raised. Then gab = Iab and all is well so long as Tab has dimensions of energy density.

We have chosen to use a metric of positive signature. If one wished to work with a negative signature instead:

• ga b → - ga b
• Γab c does not change sign
• Ra b cd does not change sign
• Ra b c d → - Ra b c d
• Ra b does not change sign
• R → - R
• Ta b does not change sign

so Einstein's Equations become

Ra b - R ga b / 2 - Λ ga b = α Ta b
Taking the trace of both sides with Λ = 0, we have
(1 - D / 2) R = α T
This fact implies that for vacuum solutions (where Ta b is zero), Ra b (and all tensors and scalars derived from it) must be zero. This means that the Weyl Tensor is equal to the Riemann Tensor for such solutions and so we have the interpretation of the Weyl Tensor as that portion of the curvature which is not due to local stress-energy. It also implies that for vacuum solutions with nonzero cosmological constant, we have
R = D Λ / (D/2 - 1)

We can use the following Mathematica code to see if any given metric is a solution:

simpler[ Table[ Table[
Rll[[ii, jj]] - R gll[[ii, jj]] / 2 + L gll[[ii, jj]],

{jj, dim}], {ii, dim}] === a * Tll]

In this example, L is the cosmological constant, and a replaces α as the coefficient of the stress-energy tensor. The "===" operator ensures that Mathematica will return a simple True or False. Of course this code only tests for solutions on a specific coordinate chart; the process of solving Einstein's Equations involves identifying the manifold or submanifold covered by a given chart on which a particular metric solves the equations. Some things must still be done by human reasoning!
εa b c d ... = | g |-1/2 Sign(a,b,c,d ...)
and can use the following Mathematica code to implement it in D = 4:
epsilonuuuu[ a_, b_, c_, d_] := Signature[ {a, b, c, d}] / simpler[ Sqrt[ +-Det[ gll]]];
Since Mathematica indices must start with 1, we assume that the timelike index corresponds to 4 (or in general, D). This means that in even dimensions this function will produce an extra minus sign relative to the standard usage of 0 for the timelike index. The minus sign is selected for metrics of Lorentz signature. We can of course define the symbol analogously for any fixed D.

The next section discusses curvature invariants.