We begin by turning off Mathematica's spelling errors, since we will be using many variables which are spelled similarly:

Off[ General::spell, General::spell1];We now use the

We began by using the MathematicaClear[ der];der[ expr_?NumericQ, x_] := 0;

der[ expr_Plus, x_] := Sum[ der[ expr[[iind]], x], {iind, Length[expr]}];

der[ expr_Times, x_] := Sum[

Product[ expr[[jind]], {jind, iind - 1}] * der[ expr[[iind]], x] *Product[ expr[[jind]], {jind, iind + 1, Length[expr]}], {iind, Length[expr]}];

We repeat the process for the variational operator (usually denoted δ):

We further specify that the variation of a derivative is the derivative of the variation:Clear[ var];var[ expr_?NumericQ] := 0;

var[ expr_Plus] := Sum[ var[ expr[[iind]]], {iind, Length[expr]}];

var[ expr_Times] := Sum[

Product[ expr[[jind]], {jind, iind - 1}] * var[ expr[[iind]]] *Product[ expr[[jind]], {jind, iind + 1, Length[expr]}], {iind, Length[expr]}];

var[ der[ expr_, x_]] := der[ var[expr], x];Note that we dovar[ covder[ expr_, x_]] := covder[ var[expr], x];

If the variable "indi" has not been used before, it is initialzed to 1. Otherwise it is incremented and used to construct the new index symbol. The first such symbol is "aa", followed by "bb", etc. After "zz" is generated the next index will be "aaa", and so on. We now define a function whose purpose is to return the unsorted union of a list, removing redundant members:ind:= Block[ {abet, num, str, index}, (abet = {"a", "b", "c", "d", "e", "f", "g", "h", "i", "j", "k", "l","m", "n", "o", "p", "q", "r", "s", "t", "u", "v", "w", "x", "y", "z"};If[ Head[indi] === Symbol, indi = 1, indi = indi + 1];If[ indi > 26, index = Mod[indi - 1, 26] + 1, index = indi];

num = Quotient[indi + 25, 26] + 1;

str = "";

Do[ str = str <> abet[[index]], {ii, num}];

Return[ Symbol[str]];)];

This function is an example of "programming by unintended consequences". We eschew this programming technique in general but borrow it here from the Mathematica Help Browser because it was too arcane to pass up (see the Further Examples under the Reap built-in function).unsun[ l_List] := Reap[ Sow[ 1, Flatten[l]], _, #1&][[2]];

We will be analyzing expressions which contain a number of both free and contracted indices. In many
cases, simplification of these expressions will be very hard for Mathematica because it does not know
to rename indices in order to facilitate factoring. We therefore define a function **iclean**, which
will *attempt* to clean up contracted indices and do preliminary simplification. *Note that there
is no easy way to do this in a foolproof manner*; some manual simplification will inevitably be necessary.
The arguments to **iclean** will be an expression to "clean up" and a function to be Collected on.

We begin by defining an identity function **symind** whose eventual purpose will be to apply
symmetrization rules to an expression. We define it here so that **iclean** will not have to know about
the specific tensors we will be working with; it will be redefined below as we teach Mathematica about
General Relativity:

"ex" contains either the expansion of the argument expression, or if it was a simple product, a list whose only member is the argument. The loop processes each term in the expansion. "modex" is a modified version of "ex", where all occurences of arithmetic, derivative or determinant functions have been eliminated; this leaves only those function invocations which we wish to simplify. "ftab" is a table of the functions appearing in this term, with the collected function first. "funcs" is a list of the function names, and "inds" is a list of the indices which are arguments to them. "uinds" is used to select from inds those indices which are repeated; in this way we do not modify free indices. Setting "indi" to zero so that the indices generated for each term bysymind[ expr_] := expr;

iclean[ expr_, coll_] := Block[ {ex, modex, len, res, ftab, funcs, inds, uinds},(ex = Expand[ expr];If[ Head[ex] === Plus, len = Length[ex], (ex = {ex}; len = 1;)];

res = 0;

Do[(

modex = ex[[iind]] /. {Times -> List, Plus -> List, Power -> List} /.Return[{der -> List, covder -> List} /. det[g_] -> 1;ftab = Select[unsun[ Join[funcs =Select[ modex, (Head[#1] === coll)&],(!(NumericQ[#1] || (Head[#1] === Symbol)))&];Select[ modex, (Head[#1] =!= coll)&]]

/. var -> Sequence],

unsun[ Table[ Head[ ftab[[jind]]], {jind, Length[ftab]}]];inds = Cases[ ftab //. Thread[ funcs -> Sequence], _Symbol];

uinds =

unsun[ inds];inds = Select[ uinds, (Count[inds, #1] > 1)&];

indi = 0;

res = res + (ex[[iind]] /.

ReleaseHold[ Thread[ inds ->Hold[ind]]]);),{iind, len}];

symind[ Collect[res, coll[al__]]]];)];

Now it is time to teach Mathematica something about General Relativity. We begin by defining the
variation of the square root of the determinant of the metric, using the **ind** function to
generate summation indices:

var[ Sqrt[- det[g]]] := - Sqrt[- det[g]] g[ aind =This is of course based on the definition of the determinant as the exponentiation of the the trace of the log of the metric tensor:ind, bind =ind]var[ gi[ aind, bind]] / 2;

δ det(g) = δ exp( tr( ln(g)))We now redefine= exp( tr( ln(g))) δ tr( ln(g))= det(g) g

^{a b}δ g_{a b}

Again note the use ofsymind[ expr_] := Map[ Collect[#1, Sqrt[- det[g]]]&,ReleaseHold[ expr /.g[a_, b_] -> g[Hold[ Sort[{a, b}]]] /.gi[a_, b_] -> gi[

Hold[ Sort[{a, b}]]] /.G[a_, b_, c_] -> G[

Hold[ Join[ {a}, Sort[{b, c}]]]] /.R[a_, b_] -> R[

Hold[ Sort[{a, b}]]]] /. List -> Sequence];

We can now begin to derive Einstein's Equations. We define the **Einstein-Hilbert Lagrangian**:

L := - Sqrt[- det[g]] Rs;where "Rs" is the scalar curvature. We next find the variation of that Lagrangian in terms of the Ricci Tensor and the inverse metric:

varL =For arbitrary variation of the (inverse) metric, the quantity in parentheses must vanish for the Lagrangian to be stationary: this gives us the vacuum Einstein Equations. However, the variation of the Ricci Tensor must also vanish, and so we express it in terms of the Christoffel Symbols (by way of the Riemann Tensor):iclean[ var[L] /. Rs -> gi[a, b] R[a, b], var] /. gi[a_, b_] R[a_, b_] -> RsSqrt[- det[g]] (1/2 Rs g[aa, bb] - R[aa, bb]) var[ gi[aa, bb]]- Sqrt[- det[g]] gi[aa, bb] var[ R[aa, bb]]

varR =We next re-express the derivatives in terms of covariant derivatives (since the variation of the Christoffel Symbols is a tensor):iclean[ Sqrt[- det[g]] gi[aa, bb] var[ R[aa, bb]] /.R[a_, b_] -> R[a, cind =ind, b, cind] /.R[a_, b_, c_, d_] -> der[ G[d, a, c], x[b]] - der[ G[d, b, c], x[a]] +

G[find =

ind, a, c] G[d, find, b] - G[find, b, c] G[d, find, a], gi]Sqrt[- det[g]] gi[ aa, bb] (der[ var[ G[cc, aa, bb]], x[cc]] -der[ var[ G[cc, bb, cc]], x[aa]] +G[cc, cc, dd] var[ G[ dd, aa, bb]] - G[cc, bb, dd] var[ G[dd, aa, cc]] -

G[cc, aa, dd] var[ G[dd, bb, cc]] + G[cc, aa, bb] var[ G[dd, cc, dd]])

varRa =Since the covariant derivative of the (inverse) metric is zero (we are assuming a symmetric metric), we can bring it inside the covder function:iclean[ varR /. der[ var[ G[a_, b_, c_]], x[d_]] -> covder[ var[ G[a, b, c]], x[d]] -G[a, eind =ind, d] var[ G[eind, b, c]] + G[eind, c, d] var[ G[a, eind, b]] +G[eind, b, d] var[ G[a, eind, c]], gi]

(covder[ var[ G[cc, aa, bb]], x[cc]] -covder[ var[ G[cc, bb, cc]], x[aa]]) Sqrt[- det[g]] gi[aa, bb]

varRb = Expand[ varRa] /. covder[a_, b_] gi[c__] -> covder[ gi[c]a, b]These are covariant divergences of vectors which can converted to ordinary divergences by noting thatcovder[ gi[aa, bb] var[ G[cc, aa, bb]], x[cc]] Sqrt[- det[g]] -covder[ gi[aa, bb] var[ G[cc, bb, cc]], x[aa]] Sqrt[- det[g]]

(- det(g))where the comma denotes an ordinary derivative. We examine the left hand side of this equation, expanding the Christoffel Symbol in terms of the metric and its inverse:^{1/2}V^{a}_{; a}= ((- det(g))^{1/2}V^{a})_{, a}

and see that it is the same as the right hand side by using the variation of the determinant as an ordinary derivative:iclean[ Sqrt[- det[g]] (der[ vu[a], x[a]] + G[a, bind =ind, a] vu[bind]) /.G[a_, b_, c_] -> gi[a, find =ind] * (der[ g[c, find], x[b]] +der[ g[b, find], x[c]] - der[ g[b, c], x[find]]) / 2, Sqrt]

Sqrt[- det[g]] (der[ vu[aa], x[aa]] +1/2 der[ g[aa, bb], x[cc]] gi[aa, bb] vu[cc])

Applying this equivalence, we see that the variation of the Ricci Tensor is an ordinary divergence, which vanishes on a manifold without boundary (or if the variations of the metric vanish on the boundary):iclean[ Simplify[ der[ Sqrt[- det[g]] vu[a], x[a]] /.der[ Sqrt[- det[g]], x[a]] -> - der[ det[g], x[a]] / (2 Sqrt[- det[g]]) /.der[ det[g], x[a]] -> det[ g] der[ g[aa, bb], x[a]] gi[aa, bb]], Sqrt]

Sqrt[- det[g]] (der[ vu[aa], x[aa]] +1/2 der[ g[aa, bb], x[cc]] gi[aa, bb] vu[cc])

varRb /. Sqrt[- det[g]] covder[a_, b_] -> d[ Sqrt[- det[g]]a, b]While these techniques are powerful, it is important to examine every result carefully: it is still easy to accidentally make a mess of the indices, which can result in corruption or loss of terms!d[ Sqrt[- det[g]] gi[aa, bb] var[ G[cc, aa, bb]], x[cc]] -d[ Sqrt[- det[g]] gi[aa, bb] var[ G[cc, bb, cc]], x[aa]]

ClearAttributes[ Times,If this is not done, or if you accidentally reset the attribute, you will soon have operators commuting which should not be:Orderless];Expand[ - z y x /. (y x) -> x y - Comm[ x, y]]

- z x y + z Comm[ x, y]

SetAttributes[ Times,Orderless];Expand[ - z y x /. (y x) -> x y - Comm[ x, y]]

- x y z + z Comm[ x, y]

The next appendix is about parallel algebraic computation.

- Table of Contents
- Index:

©2005, Kenneth R. Koehler. All Rights Reserved. This document may be freely reproduced provided that this copyright notice is included.

Please send comments or suggestions to the author.