xAct`xTensor`

Author

© 2002-2018, under the GNU General Public License (GPL)
     José M. Martín-García
	jose@xact.es
	http://www.xAct.es
		http://groups.google.com/group/xAct/

Helped by:
	Alfonso García-Parrado : design of vector bundles, complex structure, products and mappings.
	Guillaume Faye : design of n+m decompositions.

Intro

xTensor` extends Mathematica capabilities in abstract tensor calculus, focusing on General Relativity. It works with tensor fields with arbitrary symmetries under permutations of indices, defined on one or more vector bundles and manifolds. It computes covariant derivatives with respect to multiple connections, including gauge connections in complex vector bundles. xTensor` can work with Lie derivatives, parametric derivatives and variational derivatives. It allows the presence of a metric in each manifold and defines all the associated torsion and curvature tensors (Riemann, Ricci, Einstein, Weyl, etc.), handling automatically their definitions and relations in terms of Christoffel objects. It can work with mappings between any two manifolds and push/pull tensor fields along those maps.
xTensor` does not perform component calculations. See the twin package xCoba`.
xTensor` needs the twin package xPerm` in order to compute the canonical form of a list of indices under certain symmetry group.
xTensor` has been designed with the following priorities in mind:
	1.- Mathematical structure
	2.- Efficiency
	3.- Compliance with Mathematica style
	4.- Simplicity of input/output
In particular, concerning the mathematical structure, xTensor` imitates the usual cycle in mathematics ``let X be a (Type) with properties props; then use it to do some computations´´ with the general form of declaration
	DefType[X, props]
	computations

Load the package

This loads the package from the default directory, for example $HOME/.Mathematica/Applications/xAct/ for a single-user installation under Linux (see the installation notes for other possibilities),

<<xAct`xTensor`

0. Basics

0.1. Example session

This is a very simple example of the use of xTensor`. We start by defining a manifold and some tensor fields on its tangent bundle.

Define a 3d manifold M3:

DefManifold[M3,3,{a,b,c,d,e}]

Define a contravariant vector v and a covariant antisymmetric tensor F:

DefTensor[v[a],M3]
DefTensor[F[-a,-b],M3,Antisymmetric[{-a,-b}]]

Tensor product of both tensors, based on the abstract index notation:

F[-a,-b]v[b]

% v[a]

There is no automatic simplification. It must be explicitly asked for:

Simplification[%]

Example of dummy handling:

F[-a,-d]v[d]v[b]v[-b]+3v[b]v[c]v[-c]F[-b,d]F[-d,-a]

%//Simplification

There is not a metric and hence v cannot be covariant. Note the position of the head ERROR:

Validate[%]

Internal structure:

InputForm[%%]

Now we define a covariant derivative operator and check the first Bianchi identity (valid for any symmetric connection). Note that xTensor` uses unique (dollar-) variables as internal dummy indices to avoid index collisions. The function ScreenDollarIndices hides those indices replacing them by "normal" indices.

Covariant derivatives are defined by default with curvature but without torsion. Note also the default associated symbols: semicolon in Postfix notation and nabla in Prefix notation:

Options[DefCovD]

DefCovD[Cd[-a]]

The Leibnitz rule is automatic:

Cd[-a][ F[-b,-c]v[c] ]

% v[b]//Simplification

The tensors have been defined with their expected symmetries:

RiemannCd[-a,-b,-c,d]+RiemannCd[-b,-a,-c,d]//Simplification

but it is possible to deduce them from more basic principles:

RiemannCd[-a,-b,-c,d]+RiemannCd[-b,-a,-c,d]//RiemannToChristoffel//Simplification

RiemannCd[-a,-b,-c,d]+RiemannCd[-b,-c,-a,d]+RiemannCd[-c,-a,-b,d]//RiemannToChristoffel//Simplification

Check of the second Bianchi identity:

Antisymmetrize[Cd[-e][RiemannCd[-c,-d,-b,a]],{-c,-d,-e}]

3%//Simplification

%//CovDToChristoffel

%//RiemannToChristoffel//ScreenDollarIndices

Simplification[%]//AbsoluteTiming

Commutation of covariant derivatives. It is possible to choose the indices (use CommuteCovDs) or just let xTensor` to bring the derivative indices into canonical order (use 
SortCovDs):

Cd[-a]@Cd[-b]@v[c]

CommuteCovDs[%,Cd,{-b,-a}]//ScreenDollarIndices

SortCovDs[%%]//ScreenDollarIndices

Now we introduce a metric on the manifold. Its associated curvature tensors are defined: RiemannCD, RiciCD, EinsteinCD, etc:

Define the metric g, with negative determinant. Its Levi-Civita connection will be called CD:

DefMetric[-1,g[-a,-b],CD,SymbolOfCovD->{"|","D"}]

Standard properties of the metric and its associated objects:

CD[-a][ g[-b,-c] ]

CD[-a][ EinsteinCD[a,b] ]

g[a,-a]

g[a,b]g[-b,-c]

Let us check the antisymmetry of the second pair of indices in the Riemann tensor:

RiemannCD[-a,-b,-c,-d]+RiemannCD[-a,-b,-d,-c]

Because canonicalization is not automatic, the previous expression is not automatically zero. We expand it in derivatives of Christoffels and canonicalize. The partial derivative does not commute with the metric, so that internally the code replaces temporarily those derivatives with Levi-Civita connections:

%//RiemannToChristoffel//Simplification//ScreenDollarIndices

Metric contraction is not automatic either. The resulting expression is not obviously zero:

%//ContractMetric//Simplification

Actually it is only zero for connections deriving from a metric:

%//ChristoffelToGradMetric//Expand

%//Simplification

In this case a simple reordering of the partial derivatives (this is not automatic) would have been enough:

%%//SortCovDs

Examples of Lie derivative and bracket:

LieD[3v[a]][F[-a,-b]v[b]]

Bracket[v[a],F[a,-b]v[b]][a]

Undefine all objects on the manifold and the manifold itself:

UndefMetric[g]

UndefTensor/@{v,F};
UndefCovD[Cd]

UndefManifold[M3]

There are no symbols left in the Global` context

?Global`*

0.2. Summary

These are the objects that currently can be defined in xTensor`.
	- Manifolds
	- Mappings
	- Vector bundles
	- Indices (of different types)
	- Tensors (in particular, Metrics)
	- Covariant derivative operators (including partial derivatives)
	- Constant symbols
	- Parameters
	- Inert heads (wrappers for tensors)
	- Scalar functions
	- Charts (coordinate systems): use the twin package xCoba` (in collaboration with D. Yllanes)
	- Bases (frames): use the twin package xCoba`
We also have predefined mathematical operators:
	- Lie derivative
	- Lie bracket
	- Parametric derivative
	- Variational derivative
and useful mathematical operations:
	- Relations and rules among indexed objects
	- Indexed equations solving
	- Canonicalization and simplification of expressions
	- Automatic definition of curvature tensors, and relations among them
	- Complex conjugation
There are also specialized algorithms, frequently used in General Relativity:
	- Warped metric decomposition
	- ADM-like decomposition
	- Perturbation theory in General Relativity: use the twin package xPert` (in collaboration with D. Brizuela and G. Mena Marugán)
	- Spherical tensor harmonics: use the twin package Harmonics` (in collaboration with D. Brizuela and G. Mena Marugán)
	- Manipulation of invariants of the Riemann tensor: use the twin package Invar` (in collaboration with R. Portugal, L. Manssur and D. Yllanes)
	- Spinor calculus (in collaboration with A. García-Parrado)

0.3. Getting info

Every piece of information in Mathematica is associated to a symbol and can be obtained using the ? and ?? symbols (see the command Information in the Mathematica help). xTensor` strictly follows that rule. (In particular, we have explicitly avoided notations based on string input/output which would require a parser, much restricting the capabilities of the user.)

?DefManifold

??ERROR

0.4. Too many brackets!

This is a common complaint from most users after their first contact with xTensor`. My answer to that is the following:

1. That is how Mathematica internally works. For example, did you expect this simple derivative having internally six pairs of brackets, when it has just one in the "user-notation"?

f'[2x-y]//FullForm

2. xTensor` does not have different internal and external notations (at least not now). This helps understanding what is happening internally and allows the user to build pattern rules which always work. However, you are continually exposed to the internal notation.

3. The highest accummulation of brackets happens in derivatives. This simply follows the general idea in Mathematica (made clear in the example above) of treating a derivative as an operator on an expression, and not as a simple additional argument of the expression, as the traditional "comma-index" notation might suggest.

4. The Mathematica language is very similar to LISP, in which you have a lot of parenthesis.

5. If you are still not convinced by the previous three arguments, then I must say that it is possible to construct whatever other notation you like for a derivative, or for any other operator, and several examples of this are given in section 6.1 below.

Remove[f,x,y]

0.5. Warnings

xTensor` is a very powerful tool in the hands of an average Mathematica user, but could be difficult for an unexperienced Mathematica user. This is because I have not tried to hide Mathematica subtleties (like delayed assignments, unique symbols or upvalues), but rather I have tried to make as much use as possible of such powerful instruments. If you don't understand those three listed concepts, please spend some (not really much) time with a Mathematica tutorial. It will pay off, even if you finally don't use xTensor`.

If you think you can do your tensor computation by hand in a few lines, then xTensor` is probably too much. This system is designed to be efficient when developing big projects, not small computations: a general-purpose tensor computer algebra system always has hundreds of commands, and only big projects give you enough motivation to learn them.

xTensor` is a manipulator of abstract expressions. It does not compute (at least easily) the components of the Riemann tensor associated with a given matrix of metric components, and things like that. xTensor` is specialized in the manipulation of identities obeyed by families of objects with some properties (like being antisymmetric). Many users feel dissapointed at the beginning because they were expecting component-type computations. Use the twin package xCoba` for that.

1. Manifolds, vector bundles and parameters

1.1. Define a manifold and its tangent bundle

Just after loading,  xTensor` does not have any manifold defined. Here we see how to declare our own abstract manifolds, on which tensor fields will be defined.

Definition of a manifold.

We define a 3-dimensional manifold M3 with abstract indices {a, b, ..., h} associated to its tangent vector bundle:

DefManifold[M3,3,{a,b,c,d,e,f,g,h}]

DimOfManifold[M3]

We can ask Mathematica about this manifold using the question mark ? (command Information). We see that the data is stored in several property-functions using the concept of "upvalue" (this is why we find the  ^= sign below)

?M3

Together with the manifold we have also defined its tangent vector bundle, which stores all the information related to indices and indexed objects. By default, the tangent bundle is named by joining the symbol Tangent with the name of the manifold, but that name can be freely choosen. To conform with the standard notation, it is typeset by default as the character 𝕋 followed by the typesetting of the manifold, but the typsetting can also be freely chosen. Note that we shall often abbreviate "vector bundle" to "vbundle".

Tangent bundle of a manifold.

By default, tangent bundles are real and have no metric. The dimension of the tangent bundle is defined as that of the vector space at each point, and coincides with the dimension of the base manifold.

Tangent[M3]

%//InputForm

?TangentM3

We define a second manifold S2. Note that we can use C, D, K, N, O as indices, even though Mathematica has reserved meanings for those symbols. The capitals E and I cannot be used as indices because they are always understood as the base of natural logarithms and the square root of -1, respectively. xTensor` issues warnings to remind you of this point:

DefManifold[S2, 2, {A,B,C,D,F,G,H}]

The lists of manifolds and vbundles defined in the current session are given by the global variables $Manifolds and $VBundles, respectively:

$Manifolds

$VBundles

A manifold can be undefined (its properties are lost and the symbol is removed). In the process, the associated tangent bundle is also undefined:

UndefManifold[M3]

?M3

$Manifolds

And redefined (note the use of the command IndexRange, quite handy for long lists of indices):

DefManifold[M3, 3, IndexRange[a,h]]

$Manifolds

It is possible to define 1-dimensional manifolds, with a corresponding 1-dimensional tangent vbundle. They are treated in exactly the same way as other manifolds and vbundles. Other systems, like Ricci`, only allow a single abstract index (which can be repeated), but that is a large internal exception, and xTensor`does not follow that idea. Instead, we shall declare all indices of the same 1-dimensional space to be symmetric. See subsection 4.4 below.

The dimension of a manifold does not have to be an integer necessarily. It can be any expression. It is advisable to use constant-symbols (see 6.5 below) as manifold dimensions.

1.2. Product manifolds. V1.1.0 UPDATE

Given two manifolds, it is always possible to define another manifold as the Cartesian product of those. The tangent bundle of the product will be then constructed as the (external, or disjoint) direct sum of the tangent bundles of the submanifolds. There are two ways to approach this problem, and we will find this duality repeatedly in all of xAct`. On the one hand we can declare a new symbolic manifold using a DefType command, which will prepare everything for us in advance. On the other hand we can describe the new manifold using a CartesianProduct command, whose properties will be computed in real time whenever they are needed.

We define a product manifold as

DefManifold[M5,{M3,S2}, {μ, ν, λ, σ, η, ρ}]

?M5

whose tangent bundle is a direct sum of the corresponding subvbundles:

?TangentM5

The list of product manifolds is given by another global variable:

$ProductManifolds

The command SplitManifold also allows constructing a product manifold, or splitting a previously existing manifold.

?SplitManifold

We can also work directly with an expression representing the product manifold.

Cartesian product of manifolds.

Cartesian product of our two manifolds :

prod=CartesianProduct[M3,S2]

ManifoldQ[prod]

DimOfManifold[prod]

Tangent[prod]

%//FullForm

Clean up:

UndefManifold[M5]

Remove[prod]

1.3. Inner vector bundles

xTensor` allows the definition of vector bundles other than tangent bundles. These ``inner vbundles´´ can be real or complex, which forces the definition of the complex conjugated vbundle. Complex conjugation is controlled by the function Dagger.

Definition of a vbundle.

Define a complex inner vbundle with dimension 4 (this is the fiber dimension) and base manifold M3:

DefVBundle[InnerC,M3,4,{𝔄,𝔅,ℭ,𝔇,𝔈,𝔉,𝔊,ℌ},Dagger->Complex]

?InnerC

In the process of definition, the conjugated vbundle InnerC†, with corresponding conjugated indices, has been created. It is a servant of InnerC, and hence can only be undefined through undefinition of the latter. By default, all conjugated symbols are constructed by appending the dagger character †. (This character is stored in the goblal variable $DaggerCharacter, which can be changed.)

?InnerC†

We have already defined four different vector bundles:

$VBundles

There is also a function SplitVBundle, parallel to SplitManifold.

?SplitVBundle

1.4. Parameters

In xTensor` the indexed objects can be fields on manifolds and/or functions of real parameters (for example proper time).

Definition of a parameter

We define a real parameter that will be used below to compute parametric derivatives (see ParamD in subsection 6.6).

DefParameter[time]

?time

Manifolds and parameters are the only possible "dependencies" of the objects in xTensor`. We shall later explain this in more detail.

An important question is the relation between fields on 1d manifolds and functions of a parameter. In decompositions of ADM-type it would be nice to have a way of going from a 4d structure to a 3+1 structure and then from here to a 3d structure plus a time parameter. I still don't know how to program this idea. There are many things to change in a nontrivial way: tensor fields, covariant derivatives, metrics, etc. Please, let me know if you have suggestions.

1.5. Vector spaces: bundles over a point

It is possible to work  xTensor` with a single vector space, considered as a bundle over a manifold that happens to be a point.

Define a zero-dimensional manifold, a point. We do not need to specify indices because there is no tangent space:

DefManifold[P,0,{}]

Now we can define a vector bundle over this point. It has a single fiber, which is a vector space of the required dimension:

DefVBundle[V,P,4,{𝔞,𝔟,𝔠}]

And now we can define vectors and tensors on this space, with no dependencies:

DefTensor[X[𝔞],{}]

Undefine all objects:

UndefTensor[X]

UndefVBundle[V]

UndefManifold[P]

2. Indices

2.0. Important comments

1) Indices are the key objects in xTensor`, and the source of the generality of the system. They are carriers of contextual information, much like the decimal digits in numerical computations, or the letters in a text. An index by itself has little information (only an association to a vector space), but a group of indices together can encode a lot of information (rank, symmetries, metrics, bases, etc). In general, an index is anything which can be placed at an index-slot: that is, we shall later define four places in which we can put indices and anything which we can meaningfully put there will be called an index.

2) Starting version 0.9 of xTensor`, indices are associated to vbundles, and not to manifolds, as was the case in previous versions. This allows us to work with vbundles other than tangent bundles to manifolds.

3) We need a "language" to talk about indices, and though we usually understand what a dummy index or a covariant index is, such a language is not rigorously defined anywhere. I have invented my own language for this. Please, suggest changes if you find inconsistencies, gaps or bad choices.

4) Taken from the Internet:
	The plural of "index" for books, journals and newspapers is "indexes"; the plural of "index" for graphs and benchmarks for economics for example, is "indices"; but the plural of "appendix" is always "appendices".
I think the "index" of a tensor must follow the rule for the second group, and hence its plural will be "indices".

2.1. Properties of indices: type, character, state

1) There are five known types of index in xTensor` (other types could be added, but there is no user-defined way to do this):
	- Abstract index (AIndex). See subsection 2.2.
	- Basis index (BIndex), associated to some basis (coordinated or not) of vector fields. Described in xCoba`.
	- Component index (CIndex), also associated to some basis. Described in xCoba`.
	- Directional index (DIndex). See subsection 3.7.
	- Label index (LIndex). See subsection 3.8.
Collectively we refer to all five classes as "generalized indices" (GIndex). Additionally, we can have patterns (PIndex) for g-indices at those positions where g-indices are expected, but patterns are not considered as g-indices. In the following "indices" will mean "g-indices" unless we specify the type.

2) All indices have a "character", which can be either covariant (Down-index) or contravariant (Up-index).

3) Some indices can be used to represent contractions, following the Einstein convention: two repeated indices in the same tensor or tensor product, but each having a different character.  We call them "contractible" or "Einstein" indices (EIndex) and currently only abstract and basis indices are allowed to be contractible. (Note that having contractible indices with the same character is a syntactic error in xTensor`.) E-indices have a "state": they can be Free or contracted (aka Dummy). For completeness, the state of components, directions and labels is said to be always "Blocked", what in particular means that they can be repeated (i.e. truly repeated, not simply staggered).

4) Basis and component indices are known by part of the routines of xTensor` (for example the formatting routines, the index selectors, the canonicalizer or the constructors of rules), but used in the twin package xCoba`. From now on we shall not consider those two types of indices. See xCobaDoc.nb for further information.
5) There is a further property of an index, only used internally: its metric-state, saying whether an index can be raised or lowered using a given metric (not whether it has been or not actually raised already). This will be important in section 7.3.

2.2. Covariant and contravariant abstract indices

Abstract indices in xTensor` are valid atomic symbols without numeric value (excluding therefore E and I) and with no special output (excluding symbols like Space, Tab, etc.) One of the most important decisions when working with tensors is the notation for covariant and contravariant indices. Given the binary character of the problem there is only one natural simple choice: the use of signs + and - . We shall denote contravariant indices as +a or a, and covariant indices as -a. This choice is very simple, but not particularly convenient when doing pattern matching, due to the asymmetry of a being a symbol but -a being a composite expression.

Validation of indices.

We can check whether a symbol has been registered as an abstract index using the function AIndexQ

AIndexQ[A]

AIndexQ[-b]

AIndexQ[q]

Remove[q]

In particular, we can specify the vbundle it must belong to:

AIndexQ[A,TangentS2]

AIndexQ[𝔄,TangentM3]

AIndexQ[𝔄,InnerC]

AIndexQ[𝔄†,InnerC]

AIndexQ[𝔄†,InnerC†]

Each vbundle has two special functions that select abstract indices on that manifold. For example for the vbundle TangentM3 they are TangentM3`Q and TangentM3`pmQ. The former only accepts contravariant indices; the latter accepts both characters. The use of context notation (the ` ) is purely historical : before using upvalues to store information I tried to use contexts, but it is less convenient.

TangentM3`Q[a]

TangentM3`Q[-a]

TangentM3`Q[A]

TangentM3`pmQ[a]

TangentM3`pmQ[-a]

2.3. Get (new) indices

Sometimes we need to use many indices. There are two ways to do it. We can register more symbols using AddIndices or we can ask xTensor` to generate new indices for a vbundle using NewIndexIn (it uses a combination of the last registered symbol and integer numbers). User-defined indices and computer-defined indices are kept in two different lists.

Functions that return or change the available indices.

Currently there are only user-defined indices. We add some more:

IndicesOfVBundle[TangentM3]

AddIndices[TangentM3,{i,j,k,l}]

IndicesOfVBundle[TangentM3]

We can also remove indices, but this is a very dangerous operation because some previous expressions could get corrupted if they contain removed indices. There is no built-in command to check whether this is happening or not.

RemoveIndices[TangentM3,{j,i,l,k}]

IndicesOfVBundle[TangentM3]

xTensor` generates new indices using NewIndexIn

NewIndexIn[TangentM3]

IndicesOfVBundle[TangentM3]

The user-function to get any number of indices is GetIndicesOfVBundle. Note that the last argument is a list of symbols that we do not want to get:

GetIndicesOfVBundle[TangentM3,14,{d,e,h3}]

IndicesOfVBundle[TangentM3]

The indices of complex conjugated vbundles are paired, and therefore if we generate an index in one of them, automatically the complex conjugated index is generated in the other one:

NewIndexIn[InnerC]

IndicesOfVBundle[InnerC]

IndicesOfVBundle[InnerC†]

GetIndicesOfVBundle[InnerC†,10]

IndicesOfVBundle[InnerC]

Note the information associated to an abstract index:

?𝔄

3. Tensors and tensor slots

3.1. Define a tensor

Now we can define tensors. In xTensor` we always work with tensor fields on zero, one or several manifolds, having indices on zero, one or several of their vbundles.

Definition of tensors.

Define a tensor field T on M3 with two contravariant indices and one covariant index, and two other tensors. The actual indices used in the definition are irrelevant, and are only meant to specify the vbundle associated to each slot

DefTensor[T[a,b,-c],M3]

DefTensor[S[a,b],M3]

DefTensor[v[a],M3]

Define a scalar field r (do not forget the empty pair of brackets!)

DefTensor[r[],M3]

Define a tensor U on M3 and S2 with two contravariant indices and one covariant index, The tensor is also a function of the parameter time:

DefTensor[U[a,b,-C],{M3,S2,time}]

The tensor is identified by a symbol and we shall associate to that symbol all definitions related to the tensor. This is internally done using UpSet (^=) rather than Set (=). The definitions can then be collected using ?

?U

A tensor can be undefined. Its properties are lost, and the symbol is removed.

UndefTensor[U]

?U

An object with indices on a given vbundle is considered to be necessarily a field on the base manifold of that vbundle. xTensor` checks that point. See the discussion below on Constant objects.

The tensor U with indices on vbundles TangentM3 and TangentS2 is defined to live only on M3. xTensor` adds the manifold S2 to the list DependenciesOfTensor[U].

DefTensor[U[a,A],M3]

?U

UndefTensor[U]

Note that, from version 1.0 of xTensor`, it is now possible to work with tensors whose head is not a symbol, but a more complicated expression. See the update in section 3.7.

3.2. Standard output

It is important to have nice-looking output expressions, mainly when we deal with very large expressions. xTensor` has built-in rules to convert tensors and derivatives into 2-dimensional boxes for StandardForm output. Those boxes can be cut and pasted, using InterpretationBox. Additionally, translation to LATEX will be shown below.

The output of tensors in StandardForm has been defined with 2-dimensional boxes, but the internal structure is kept in 1-dimensional form (as usual in Mathematica):

{T[a,b,-c],S[a,b],v[-a],r[]}

InputForm[%]

No rules have been defined for OutputForm because the output is simple enough:

OutputForm[%%]

The ouput representation of the heads can be given at definition time, or changed later. Again, this is only possible in StandardForm. When using this option, it is very important to remember that the input and output names of the same tensor are different!

Output of defined heads.

The tensor TT will be output as τ. PrintAs must give a string or a function that applied on the defined symbol gives a string

DefTensor[TT[a,-b], M3, PrintAs->"τ"]

TT[a,-c]

InputForm[%]

As an example of a printas function we can use:

FirstCharacter[symbol_]:=First[Characters[ToString[symbol]]]

DefTensor[Force[a], M3, PrintAs->FirstCharacter]

Force[a]

The formatting of a tensor can be manually changed at any time:

PrintAs[Force]^="ℱ"

Force[a]

The 2-dimensional boxes are based on a combination of three built-in's: 
	SubsuperscriptBox/GridBox, which construct the sub/super-indices structures,
	StyleBox, which arranges the indices with the proper interspaces, and
	InterpretationBox, which stores both the box-form and the input-form of the expression. Note the option Editable->False, which prevents desynchronization. 
Before version 1.0.0 we always used a formatting based on SubsuperscriptBox. Now we use a formatting based on GridBox, which allows vertical alignment of indices. This can be controlled with a global variable:

$TensorBoxes

Force[a]

ToBoxes[%]

$TensorBoxes="ScriptBox"

Force[a]

ToBoxes[%]

$TensorBoxes="GridBox";

We can cut-and-paste the output, but not edit it, nor construct a tensor directly in box-form (do this manually; the "help" format of this notebook does not allow cut and paste through different sessions):

(Subscript[ℱ,  ]^a)

InputForm[%]

Clean up:

UndefTensor[Force]

UndefTensor[TT]

Remove[FirstCharacter,symbol]

3.3. Expressions and validation

Once we use abstract indices, there is no need to have a special TensorTimes. The usual Times product is enough and has one important advantage: we can use many built-in rules in Mathematica to simplify Times expressions. As usual, contractions are given using Einstein rule of staggered up/down repeated indices.

The order of factors is irrelevant in Times expressions:

v[-a]T[a,b]

v[-a]T[a,b]+T[a,b]v[-a]

The syntax of any expression can be checked using the function Validate. This function works recursively from the most external structures to the most internal ones, calling private functions ValidateTensor, ValidateCovD, etc. That gives you an idea of which problems will be detected first.

Functions that deal with validation.

If we get no errors, the input is syntactically correct:

Validate[T[a,b,-c]S[c,d]v[-d]+r[]^3S[a,b]]

Validate[T[a,b,c]S[-c,d]v[-d]+r[]^3S[a,b]]

Some errors cannot be always localized in a single object. We do not return the expression:

Validate[T[a,b,-c]v[-c]]

Validate[S[a,b]v[-b]+v[-a]]

Because validation takes up some time, it is not automatic. xTensor` has the command Validate to be used by the user as required. However, if validation is needed in all computations it is possible to switch on automatic validation, by acting on the $Pre variable:

Here, there is not automatic validation:

S[a,b]v[b]

Validate[%]

This switches on automatic validation:

$Pre

$Pre=Validate;

S[a,b]v[b]

This switches off automatic validation:

$Pre=.

S[a,b]v[b]

Note the HoldFirst  attribute of Validate, to make sure that we validate the input expression, and not whatever it evaluates to. That means that Evaluate must be used when the expression is not explicit:

Attributes[Validate]

expr=S[a,b]v[b];

Validate[expr]

Validate[Evaluate[expr]]

Remove[expr]

UndefTensor[S]

3.4. Symmetries

There are two kinds of symmetries for a tensor. We can have monoterm (or permutation) symmetries, such that a tensor remains equal or changes sign under a permutation of its indices. Or we can have multiterm (for example cyclic) symmetries, such that a linear combination of several of those permuted tensors is equal to zero. In the former case the symmetry of a tensor of n indices can be described as a subgroup of the group {1,-1}x S_n. In the latter, the symmetry can be described as an algebra. Currently xTensor` can deal with an arbitrary case (n not too large: say 200 indices at most) of monoterm symmetries. However, it does not implement any general algorithm for multiterm symmetries, because no efficient algorithm is known for that case.
The twin package xPerm` is used to work with permutations. See xPermDoc.nb for further information.

Functions related to symmetries

Any permutation group on a set of n points can be given using a set of generators (head GenSet) or, better, a strong set of generators (head StrongGenSet). The functions Symmetric and Antisymmetric are example: From version 1.0.5, following Mathematica 9, the expressions Symmetric[inds] and Antisymmetric[inds] are inert. The use of a second argument is needed to produce actually a strong generating set. The second argument will be in red, but that is a harmless issue.

Symmetric[{a,b},Cycles]

Antisymmetric[{a,b},Cycles]

Define a symmetric tensor S on M3

DefTensor[S[a,b], M3,Symmetric[{a,b}]]

?S

Define a totally antisymmetric tensor U on M3. We can use the names of the indices or their positions in the tensor

DefTensor[U[-a,-b,-c], M3,Antisymmetric[{1,2,3}]]

?U

For example, for the tensor U the group is generated by the two signed permutations -Cycles[{1,2}] (representing a change of sign under the exchange of slots 1 and 2) and -Cycles[{2,3}] (representing a change of sign under the exchange of slots 2 and 3). Again, Cycles is colored in red because there is a Mathematica symbol with the same name.

Note however that not every generating set is valid:

For example, in a group of three indices, antisymmetry in the first two indices is not consistent with symmetry in the second and third indices, unless the whole tensor is zero:

DefTensor[W[a,b,c],M3,GenSet[-Cycles[{1,2}],Cycles[{2,3}]]]

The symmetries must be consistent with the character of the indices (irrespectively of whether there is a metric or not). This is checked.

We cannot define a symmetric tensor with two indices of different character:

DefTensor[Q[a,-b],M3,Symmetric[{1,2}]]

unless we force it (so that the check is not performed):

DefTensor[Q[a,-b],M3,Symmetric[{1,2}],ForceSymmetries->True]

For symmetry groups of order larger than 10000 the check would take too long and is not performed:

DefTensor[supereta[a,b,c,d,e,-f,-g,-h,-h1],M3,Antisymmetric[Range[8]]]

UndefTensor/@{Q,supereta,U};

Remove[W]

3.5. delta, Gdelta and Sdelta. V1.1.0 UPDATE

The identity tensor on vbundles is called delta and is the same for all vbundles. This has changed from old versions of xTensor`, in which each vbundle had its own identity tensor.
The delta is a very special object. For example Penrose & Rindler use a delta object with one index on top of the other, rather than staggered, not following the "sorted slots" convention. That means that indices of such an object cannot be moved up and down with a metric. Such a construction is not natural in xTensor and hence we must look for an alternative. Up to version 0.9.6 delta was considered to be always a symmetric tensor, though this didn't make much sense in the absence of a metric, because an object with mixed indices cannot have symmetry. The problem becomes even worse in the presence of an antisymmetric metric (spinor and symplectic calculus), and hence this has now changed. The main requirement is that delta is a tensor and has a preferred index configuration, which will be taken as delta[-a, b], which is more natural in the usual conventions of spinor calculus. In the presence of a metric other configurations are possible, and actually delta becomes equal to the metric. That means that in the presence of a symmetric metric delta is symmetric, and in the presence of an antisymmetric metric delta is antisymmetric.

Identity tensors.

The delta tensor on the the vector bundle TangentM3:

delta[-a,b]

The opposite configuration is forbidden without a metric:

delta[a,-b]//Validate

There is automatic simplification with most of the expected rules for the delta tensors:

delta[-a,a]

delta[-a,b]delta[-b,c]

delta[-b,a]v[-a]

When both indices can be contracted, the second index of delta is contracted first:

delta[-a,b]S[a,-b]

The indices of delta can never have the same character, unless automatic conversion to a metric is possible:

delta[a,b]//Catch

In version 0.9.0 xTensor` added a generalized delta tensor, with two sets of antisymmetric slots. In version 1.1.0 its contraction properties have been improved. In particular we now try to avoid expanding generalized deltas into single deltas whenever possible.

This is the generalized tensor with 2 lower indices and 2 upper indices:

Gdelta[-a,-b,c,d]

It is antisymmetric in the first group of indices and separately antisymmetric in the second group of indices. Note that there is no group-exchange symmetry because there is no metric on this vbundle yet.

{Gdelta[-b,-a,c,d],Gdelta[-a,-b,d,c],Gdelta[-d,-c,b,a]}//ToCanonical

Staggered indices are automatically contracted (there is no need of a metric for this):

Gdelta[-a,-b,-c,a,b,c]

Gdelta[-a,-b,-c,d,e,c]

Contractions in the same group give zero. A metric is actually needed here, and Gdelta does not check its existence.

Gdelta[-a,-b,a,-c,d,e]

It is always possible to convert it into a determinant of normal delta tensors:

ExpandGdelta[Gdelta[-a,-b,d,e]]

Conversion into delta is automatic:

Gdelta[-a,b]//InputForm

We can contract two generalized deltas into another generalized delta:

Gdelta[-a,-b,c,d]Gdelta[-c,-d,e,f]

Gdelta[-a,-b,-c,d,e,f]Gdelta[-d,-e,-f,g,h,h1]

Gdelta[-a,-b,-c,d,e,f]Gdelta[-d,-e,-g,h,h1,h2]

Let us check this last result.

ExpandGdelta[%]

%-ExpandGdelta[Gdelta[-a,-b,-c,d,e,f]]ExpandGdelta[Gdelta[-d,-e,-g,h,h1,h2]]

%//Expand

We can also contract generalized deltas with antisymmetric tensors, without expanding into standard deltas:

DefTensor[U[a,b,c],M3,Antisymmetric[{1,2,3}]]

U[a,b,c]Gdelta[-a,-b,-c,e,f,g]

Again, let us check the result:

U[a,b,c]ExpandGdelta[Gdelta[-a,-b,-c,e,f,g]]

%//Expand

%//ToCanonical

UndefTensor[U]

In version 1.0.0, xTensor` has added yet another delta tensor, the "superdelta", which accepts any symmetry group. This object is still under development.

This tensor projects any other tensor onto an expression with the symmetries of the Riemann tensor :

Sdelta[RiemannSymmetric[{1,2,3,4}]][-a,-b,-c,-d,e,f,g,h]

%//ExpandSdelta

3.6. Directional indices

From the abstract point of view, covariant (contravariant) indices of tensors represent slots where vectors (covectors) can be "contracted". It is interesting sometimes to write that contraction explicitly: xTensor` uses the head Dir to denote that kind of "directional index". The introduction of a special head is a technical decision: we expect this feature to be seldom used and therefore it makes no sense to force xTensor` to check continuously whether a given slot contains an abstract index or a vector. All definitions concerning directional-indices are associated to Dir, and hence they are not even considered unless the head Dir is found.

Directional indices.

This is the usual representation for a contraction:

v[-b]S[a,b]

But it can also be given as (note the different results of the validation routine below)

S[a,Dir[v[-c]] ]

Validate[%]

Directional slots containing vectors are represented downstairs, and slots containing covectors are represented upstairs:

S[a,Dir[v[c]] ]

Validate[%]

S[a,Dir[v[-c]] ]

Validate[%]

It has the expected properties for a tensorial slot:

DefTensor[w[a],M3]

S[a,Dir[3v[-c]+w[-c]] ]

UndefTensor[w]

In particular, xTensor` understands the whole expression as a vector, and not as a 2-tensor:

UpVectorQ[S[a,Dir[v[-c]]]]

It is important to realize that the free index c in Dir[v[c]] is completely irrelevant except for its sign and the tangent bundle it belongs to. It is there for consistency reasons: it identifies the argument of Dir as an upvector on TangentM3. I call those indices "ultraindices".

If you find annoying to write Dir[v[c]] many times, you can name this expression:

dirv=Dir[v[c]];
S[a,dirv]

Finally, there is a pair of functions to convert the Dir indices into normal contractions, and viceversa :

SeparateDir[%,v]//ScreenDollarIndices

ContractDir[%,v]

Remove[dirv]

3.7. Label indices

Sometimes we want to have indices with no tensorial meaning. For instance the labels l, m of the spherical harmonics can be defined as label indices. This is done using a new head LI, following similar ideas to those bringing us to the head Dir. There are no definitions associated to LI; the user can use LI for any purpose he likes.

Label indices.

We can define the spherical harmonics as

DefTensor[Y[LI@l,LI@m],S2]

?Y

The up/down character of the labels only indicates the output position of the indices. Note that the signs inside LI are part of the label, and do not indicate the character.

Y[LI@L,LI@M]

{Y[-LI[L],-LI[2]],Y[LI[-L],LI[-2]]}

UndefTensor[Y]

In version 1.0.0 xTensor` has added the possibility of having tensors with nonatomic heads. However, this is not accepted in DefTensor, which means that we need to declare by hand all properties of the tensors.

For example, we could have (symmetric) tensor spherical harmonics, with a variable number of indices, expressed as (note the use of upvalues)

xTensorQ[Y[l_,m_,s_]]^=True;
SlotsOfTensor[Y[l_,m_,s_]]^:=Table[-TangentS2,{s}];
DependenciesOfTensor[Y[l_,m_,s_]]^={S2};
SymmetryGroupOfTensor[Y[l_,m_,s_][inds___]]^:=Symmetric[Range[Length[{inds}]]];
PrintAs[Y[l_,m_,s_]]^:=SubscriptBox["Y",ToBoxes[{l,m,s}]];

Y[3,-2,1][a,b,c]

Sometimes we do not need to define values for all functions. The only essential functions are xTensorQ, SlotsOfTensor and DependenciesOfTensor, that identify the tensor field character of the object. Then SymmetryGroupOfTensor is needed if the tensor has a symmetry, and PrintAs if special formatting is needed. The only other functions that might be needed are Dagger and WeightOfTensor, the latter having a default of zero. The remaining function are used by xTensor` for bookkeping purposes, namely DefInfo, HostsOf and TensorID. These non-atomic-head tensors cannot be defined via DefTensor. They must necessarily be defined by hand.

Remove[Y,l,m,s,inds,L,M]

3.8. Complex conjugation

Complex conjugation is performed by the function Dagger, except on indices, where we shall use the variant DaggerIndex, to avoid overloading Dagger too much.

Complex conjugation

A tensor with indices on a complex vbundle always has a dagger pair:

DefTensor[V[a,𝔅],TangentM3,Dagger->Complex]

?V

?V†

Complex conjugation acts on the name of the tensor and on its indices as well:

Dagger[ V[a,ℭ] ]

UndefTensor[V]

Conjugation of indices (this should be used only for programming):

DaggerIndex/@{a,-b,𝔄,-𝔅†}

Tensors can be defined with different properties against complex conjugation.

On real vbundles tensors are real by default:

Options[DefTensor]

DefTensor[V[a],M3]

Dagger[3I V[a]]

UndefTensor[V]

However, it is possible to "complexify" a real vbundle, in such a way that we duplicate its real dimension. Then vectors can have nontrivial complex conjugates:

DefTensor[V[a],M3,Dagger->Complex]

Dagger[3I V[a]]

UndefTensor[V]

If the vbundle is complex then all tensors on it must be complex:

DefTensor[V[𝔄],M3]

DefTensor[V[𝔄],M3,Dagger->Complex]

Dagger[3I V[𝔄]]

UndefTensor[V]

On complex vbundles it is possible to define Hermitian tensors if they have equal number of indices on both conjugated vbundles (though there is no built-in operation of Hermitian conjugation):

DefTensor[V[𝔄,𝔅,ℭ†,𝔇†],M3,Dagger->Hermitian]

Dagger[V[𝔄,𝔅,ℭ†,𝔇†]]

Dagger[V†[𝔄†,𝔅†,ℭ,𝔇]]

DefTensor[W[𝔄,𝔅,ℭ†],M3,Dagger->Hermitian]

UndefTensor[V];

Note that the transposition involved for Hermitian objects does not require any particular ordering of indices of the tensors. However it follows the convention of exchanging the first slot of a given complex vbundle with the first slot of its complex conjugate:

DefTensor[V[𝔄,𝔅,ℭ†,𝔇†,𝔈,𝔉†],M3,Dagger->Hermitian]

Dagger[V[𝔄,𝔅,ℭ†,𝔇†,𝔈,𝔉†]]

UndefTensor[V]

The symbol for conjugation can be changed using the global variable $DaggerCharacter, but of course the command for complex conjugation will still be called Dagger.

$DaggerCharacter

Remove[W]

3.9. Tensors with a variable number of slots

Sometimes it is desirable to work with tensors with a variable number of slots, and xTensor` supports it. For example, this is important when dealing with tensor harmonics, and actually the special package xAct`Harmonics` makes use of it.
In general, once a tensor has been defined with a number of slots and the symmetry properties of those slots, it is possible to use the tensor with more slots, and the system will assume that there are no additional symmetry properties on the added slots. The problem, however, comes when we want all slots to participate in the symmetry of the tensor. The solution is very simple:
The symmetry of a tensor is given by the function SymmetryGroupOfTensor, which stores the symmetry as an upvalue for the name of the tensor. However, this function tries first to get the symmetry from the whole expression of the tensor (indices included), and only if this is not defined it will use the stored definition. Therefore we can give additional symmetry definitions for the tensor which will overwrite the initial symmetry assignment. From version 0.9.4 we include a special notation to indicate from the beginning the presence of a variable number of slots in vbundle.

Define a tensor with a variable number of indices on a given vbundle:

DefTensor[Z[AnyIndices@TangentM3],M3]

?Z

The tensor is given no symmetry. By default we return the identity group:

SymmetryGroupOfTensor[Z]

SymmetryGroupOfTensor[Z[]]

Now make it symmetric for any number of indices:

SymmetryGroupOfTensor[Z[inds___]]^:=Symmetric[Range[Length[{inds}]]]

SymmetryGroupOfTensor[Z[]]

SymmetryGroupOfTensor[Z[a,b]]

SymmetryGroupOfTensor[Z[a,b,c,d]]

Clean up:

UndefTensor[Z]

Remove[inds]

In case you need symmetries other  than total symmetry or total antisymmetry, you will need to use the general notation for permutation groups (see xPermDoc.nb for a description).

3.10. Clean up

It is good practice in xTensor, and in all of xAct, to remove those symbols that are no longer needed. This makes the symbols available for other uses, and also frees memory that may be used to store properties of those symbols.
This documentation notebook is quite large, and hence at the end of some sections we will have “Clean up” sections removing symbols used in examples in the previous sections.

Undefine tensors:

$Tensors

UndefTensor/@%;

We cannot yet undefine the manifold M3 because other geometric structures (in this case an inner vector bundle) are defined on it. We call them “visitors” of the “host” structure, the manifold M3 in this case.

Catch[UndefManifold[M3]]

Undefine the inner vbundle explicitly:

UndefVBundle[InnerC]

Now we can undefine the manifolds:

UndefManifold[S2]

UndefManifold[M3]

And the parameter:

UndefParameter[time]

Our session is empty again:

?Global`*

TODO: Where are those two symbols coming from? They should not be there...

4. Canonicalization

4.1. Basics

In any calculation it is essential to bring all tensors and products of tensors down to a canonical form, in order to be compared. This is done by the function ToCanonical, probably the most important ingredient of xTensor`, and certainly the most complicated function.
Do not confuse the canonicalization proces, by which an expression is brought to a canonical form (not necessarily simple) from the simplification process, by which an expression is rewritten as an equivalent, but simpler form. The former is uniquely defined once the canonical form of every sintactycally correct expression is decided. The latter is largely subjective, and generally more difficult to work with. Fortunately Mathematica has a very good algorithm for simplification (Simplify). This is one of the main reasons to avoid introducing TensorTimes or TensorPlus, which would force us to build a new TensorSimplify, and another strong reason to use abstract index notation. In xTensor` we have defined a command Simplification, which simply applies ToCanonical first, and then Simplify.
The hardest part of the process (the canonicalization of indices) is performed by the companion package xPerm`. The canonicalization code is duplicated in that package: there is first a pure Mathematica version of the code, and then there is an external C-executable, which is much faster, but requires a MathLink connection, which does not work for some operating systems. Which code is used is controlled with several switches, as we will see later on.

Simplification functions

Define a manifold and some tensors fields:

DefManifold[M3,3,{a,b,c,d,e,f,g,h,i,j}]

DefTensor[U[-a,-b,-c],M3,Antisymmetric[{1,2,3}]]

DefTensor[S[a,b],M3]

DefTensor[v[a],M3]

Canonical form of the totally antisymmetric tensor U: sort indices in alphabetical order:

ToCanonical[ U[-a,-d,-b] ]

U[-a,-b,-c]+U[-c,-b,-a]//ToCanonical

Dummy indices are replaced so that the total number of different dummies is minimal:

S[a,b]v[-b]+S[a,c]v[-c]//ToCanonical

The process of canonicalization is far from trivial (see next section). We can monitor it using Verbose options.

Verbosing options of ToCanonical.

These three cells show the output of the three Verbose options of ToCanonical in a simple case. The output of xPermVerbose can be very large even for small expressions. We do not use the external executable because there is not verbose information from it:

SetOptions[CanonicalPerm,MathLink->False]

ToCanonical[U[b,c,a]v[-b],TimeVerbose->True]

ToCanonical[U[b,c,a]v[-b],Verbose->True]

ToCanonical[U[b,c,a]v[-b],xPermVerbose->True]

The canonicalization of expressions involving a large number of indices is a slow process, mainly when most of those indices are dummies. In those cases the Mathematica code for xPerm`, which is interpreted, becomes too slow. A C executable called xperm is used to compute the hardest part of the calculations; xTensor` communicates with it using a MathLink connection. Currently this executable exists for Linux and Mac for any version of Mathematica, and for MS Windows starting Mathematica 6.0.

In order to be able to connect to the external executable, the following global variable must return True. If you get False, contact JMM to see what is the problem, and whether it can be solved.

$xpermQ

MathLink connection.

Let us show an example of the difference in timings between the interpreted Mathematica code and the external executable.

We define an antisymmetric tensor Anti:

DefTensor[Anti[a,b],M3,Antisymmetric[{a,b}],PrintAs->"A"]

Now we fake that there is a metric on TangentM3. We shall later explain how to define a metric properly.

MetricsOfVBundle[TangentM3]^={metricg};
SymmetryGroupOfTensor[metricg]^=Symmetric[{1,2}];

By default the external executable is used, if the link was possible,

Options[CanonicalPerm]

Let us temporarily use the internal interpreted Mathematica code to see the difference in timings :

SetOptions[CanonicalPerm,MathLink->False]

Powers of an odd number of antisymmetric tensors are 0 :

Options[ToCanonical]

AbsoluteTiming[ToCanonical[Anti[a,-a],TimeVerbose->True]]

AbsoluteTiming[ToCanonical[Anti[a,b]Anti[-b,-a],TimeVerbose->True]]

AbsoluteTiming[ToCanonical[Anti[a,b]Anti[-b,-c]Anti[c,-a],TimeVerbose->True]]

AbsoluteTiming[ToCanonical[Anti[a,b]Anti[-b,-c]Anti[c,d]Anti[-d,-a],TimeVerbose->True]]

AbsoluteTiming[ToCanonical[Anti[a,b]Anti[-b,-c]Anti[c,d]Anti[-d,-e]Anti[e,-a],TimeVerbose->True]]

AbsoluteTiming[ToCanonical[Anti[a,b]Anti[-b,-c]Anti[c,d]Anti[-d,-e]Anti[e,f]Anti[-f,-a],TimeVerbose->True]]

AbsoluteTiming[ToCanonical[Anti[a,b]Anti[-b,-c]Anti[c,d]Anti[-d,-e]Anti[e,f]Anti[-f,-g]Anti[g,-a],TimeVerbose->True]]

These are the corresponding timings with the C executable:

SetOptions[CanonicalPerm,MathLink->True]

AbsoluteTiming[ToCanonical[Anti[a,b]Anti[-b,-c]Anti[c,d]Anti[-d,-e]Anti[e,-a]]]

AbsoluteTiming[ToCanonical[Anti[a,b]Anti[-b,-c]Anti[c,d]Anti[-d,-e]Anti[e,f]Anti[-f,-a]]]

AbsoluteTiming[ToCanonical[Anti[a,b]Anti[-b,-c]Anti[c,d]Anti[-d,-e]Anti[e,f]Anti[-f,-g]Anti[g,-a]]]

AbsoluteTiming[ToCanonical[Anti[a,b]Anti[-b,-c]Anti[c,d]Anti[-d,-e]Anti[e,f]Anti[-f,-g]Anti[g,h]Anti[-h,-a]]]

AbsoluteTiming[ToCanonical[Anti[a,b]Anti[-b,-c]Anti[c,d]Anti[-d,-e]Anti[e,f]Anti[-f,-g]Anti[g,h]Anti[-h,-i]Anti[i,-a]]]

4.2. Technical details

Canonicalization:
The canonicalization process has essentially six steps, of which only steps 4 and 5 are really hard. Assume we start from a generic syntactically correct input in xTensor` on which we act with ToCanonical:
	1) Expand the expression to convert it into a sum (Plus) of terms, each being a product (Times) of indexed objects, sharing dummy indices. This can be considered a polynomial in tensor variables, and as with any normal polynomial, it cannot be canonicalized without first expanding it. We use the Mathematica command Expand. It might happen that during the expansion some parts of the expression are evaluated and give new products to expand. That means that Expand must be used repeatedly until the expression does no longer change (we use the Mathematica command FixedPoint for that).
	2) Minimize the number of different dummies in the whole expression. We use the xTensor` command SameDummies.
	3) ToCanonical is threaded on each term. Steps 4 and 5 cannot compare different terms in the expression. From now on we only deal with products of objects.
	4) The objects in each term must be sorted according to certain priorities depending on the types and properties of the objects, but not on the actual indices they have. The term is then considered to be a single tensor, whose properties can be deduced from the properties of the elementary objects composing it. This idea is taken from R. Portugal (put journal ref. here).
	5) We extract the list of indices of that composite tensor and canonicalize that list according to the symmetries of the composite object and comparing with the "ideal" ordering of those indices (given in xTensor` by the function IndexSort). This "ideal" ordering is also decided by a number of priorities depending on the properties of the indices (type, character and state), as explained below. In general, the indices can change slot, but there is no creation of new indices, the exception being the treatment of derivatives which generate new terms. A very important issue at this step is the interaction between metrics and derivatives, discussed below. The problem will be even harder if frames are involved. The canonicalization algorithms for tensors are basically those of R. Portugal et al (put journal refs here).
	6) By now each term is already canonicalized. The function Plus will add up equal terms and the result is returned, in canonical form. No simplification (e.g. factorization) is attempted by ToCanonical.
The priorities for sorting objects in step 4 are fixed, even though they could be easily changed if a user needs it, but the priorities for sorting indices in step 5 can be configured by the user. In this section we shall describe what those priorities are and how some of them affect the efficiency of the whole process.

Sorting objects:
Objects are sorted according to seven considerations, most based in lexicographic ordering (this is also what Mathematica does). Objects are sorted with respect to the first of these criteria applying:
	1) Type: use alphabetical order in the types:
		Constant
		ConstantSymbol
		Parameter
		Scalar (see subsection 8.7 for a description of this head)
		Tensor
	     In subexpressions sometimes we need the type of a composite object: the type of a product is the product of types, and the type of a sum is the sum of types. Both Times and Plus sort their arguments lexicographically.
	2) Order (in the sense of number of elements) of the group of symmetry of the objects. More symmetric (higher order) objects come last (to improve efficiency for tensors with lots of indices), though I have found that the opposite choice is aesthetically better in general.
	3) Names of objects. Every object is given a name which depends on the names and the structure of its components, but not on its indices.
	4) Number of indices of objects. Sort first the objects with less indices. This is also more efficient.
	5) Number of free indices (of the whole term) of objects. Sort first the objects with more free indices.
	6) If the global variable $CommuteFreeIndices is False then sort first the object with the least free index according to IndexSort. If it is True (the default) then do not take into account the actual free indices and jump to step 6.
	7) If several objects cannot be sorted according to the previous considerations, then they are considered equivalent (I call them "commuting") and left then marked as CommutingObjects[o1, o2, ...], what will be used by the permutations algorithms to add new symmetries among the indices of those equivalent tensors.

Sorting indices:
These are all possible things we can consider when sorting indices, in no special order.
	- Relative position of free and dummy indices.
	- Relative position of the members of the same dummy pair.
	- Relative position of up and down indices.
	- Alphabetic order of indices.
	- The order that the indices have in IndicesOfManifold (the order at definition time).
The default is the following, chosen "experimentally" as the most efficient, in average (I don't really have an explanation of why this seems to be more efficient than other choices):
	1) free indices come before dummy indices.
	2) sort the indices according to lexicographic order.
	3) up-indices come before down-indices.
Those priorities can be changed. We need three priorities, that must be chosen from four pairs of strings, not choosing two of the same pair: "free" / "dummy", "up" / "down", "lexicographic" / "antilexicographic" and "positional" / "antipositional".

Functions that sort indices.

?IndexSort

Indices of manifold TangentM3:

IndicesOfVBundle[TangentM3]

Study this case, taking into account the default priorities given above:

IndexSort[{a,b,i,j,-e,d,c,-c,-j,e,f}]

Now we change priorities. The list is sorted differently

SetIndexSortPriorities["down","free","antilexicographic"]

?IndexSort

IndexSort[{a,b,i,j,-e,d,c,-c,-j,e,f}]

Note that once "positional" or "lexicographic" is chosen (or their opposite), then the priority "free"/"dummy" becomes irrelevant:

SetIndexSortPriorities["positional","down","dummy"]

IndexSort[{a,b,i,j,-e,d,c,-c,-j,e,f}]

The choice of priorities can make a huge difference in the efficiency of the dummy-canonicalization process. In principle, it is not possible to predict which set of priorities will be better in a given situation, but experimentally we have found that trying to match pairs of dummies with symmetry pairs of S gives better results:

expr=Anti[a,b]Anti[-b,-c]Anti[c,d]Anti[-d,-e]Anti[e,f]Anti[-f,-a]

Having one of "up"/"down" before "free"/"up" seems terrible in the Dummy algorithm:

SetIndexSortPriorities["up","free","lexicographic"]

AbsoluteTiming[ToCanonical[expr,TimeVerbose->True,MathLink->False]]

This is the default choice, instead. Note that the final result is different. This is not a problem: what we want is a uniquely defined process of canonicalization, but which one we choose is "only" a matter of efficiency and aesthetics.

SetIndexSortPriorities["free","lexicographic","up"]

AbsoluteTiming[ToCanonical[expr,TimeVerbose->True,MathLink->False]]

Remove[expr]

4.3. Benchmarking

This section contains a series of examples performed in a Linux box with a 1.7 GHz Dual core processor. At least 256 Mb RAM are required to repeat some of these examples.
We use the external executable xperm, compiled from the code xperm.c and linked from xPerm` through the MathLink protocol (xperm.tm template). The code xperm.c contains some C99 extensions, supported by GNU gcc, and therefore Linux and Mac can always use this executable. However, all other Windows C-compilers that MathLink can communicate to (Borland, MS Visual 2003, etc) are still C90 complier only. The gcc compiler under Windows (under the cygwin system) is compatible with MathLink only from version 6.0 of Mathematica, but not before and therefore pre-6.0 versions of Mathematica cannot use xperm.

This variable says whether the connection to the external executable has been performed correctly. If you get False, don't try to repeat these examples; they will take too long.

$xpermQ

Example 1:

Let us come back to the problem of products of antisymmetric tensors.

This simple code constructs the products we want to canonicalize:

productAnti[n_]:=Times@@Apply[Anti,Partition[RotateLeft@Flatten@Transpose[{-#,#}&@GetIndicesOfVBundle[TangentM3,n]],2],{1}]

productAnti[1]

productAnti[7]

We shall need up to 59 different indices:

GetIndicesOfVBundle[TangentM3,59,{}]

And this canonicalizes it, giving the AbsoluteTiming spent. Note that we can get products of up to 25 tensors in less than half a second (in that case the dummy group D has more than 5 10^32 elements).

canAnti[n_]:=Flatten[{n,With[{product=productAnti[n]},AbsoluteTiming[ToCanonical[product]]]}]

canAnti[6]

canAnti[7]

canAnti[25]

Let us see how the timings scale with the number of antisymmetric tensors. Evaluating this cell takes about a minute in a 1.8GHz Dual core processor. The executable xperm.linux requires here up to 30 Mbytes. Each computation is done only once, and that invariably introduces small deviations depending on which other tasks the computer is busy with.

data=canAnti/@Range[1,59,2]

TableForm[data,TableHeadings->{None,{"n","time","result"}}]

This is a log-log plot. A straight line represents a power-law. The fit on the last half of points seems to favour a fourth power-law scaling.

fit[x_]=Fit[Log[data[[Range[15,30],{1,2}]]],{1,x},x]

Show[ListLogLogPlot[data[[All,{1,2}]],PlotStyle->{PointSize[0.02],Hue[0.7]},DisplayFunction->Identity],Plot[fit[x],{x,0,5},DisplayFunction->Identity],DisplayFunction->$DisplayFunction,Frame->True,Axes->False]

This is a linear-log plot. A straight line represents an exponential:

ListLogPlot[data[[All,{1,2}]],PlotStyle->{PointSize[0.02],Hue[0.7]},Frame->True,Axes->False]

It is difficult to decide from the plots whether the behaviour is polynomic or exponential, but seems polynomic. This is misleading. The Dummy algorithm is based on the Intersection algorithm, which is known to be exponential in the worst case, though highly efficient. The meaning of this will become clear in the second example below.
This linear-log plot (not produced in this notebook) compares the timings of several packages on the same problem. On the x-axis, the number of antisymmetric tensors to canonicalize (both when the result is zero and when it is not zero). On the y-axis, the timings in Seconds. The color code is the following:
	Red: MathTensor
	Yello: dhPark
	Green: Tools of Tensor Calculus
	Blue (upper): xTensor with pure Mathematica code in xPerm
	Magenta: Canon (R. Portugal's own implementation in Maple of his canonicalization algorithms)
	Blue (lower): xTensor with the external C-executable xperm
The efficiency of the permutation-based algorithms with respect to more classic algorithms is clear.

UndefTensor[Anti];
Remove[productAnti,canAnti]

Example 2:

As a second example, let us work with the canonicalization of products of Riemann tensors.

Define the Riemann tensor:

DefTensor[R[a,b,c,d],M3,RiemannSymmetric[{1,2,3,4}]]

with the expected permutation symmetries:

R[a,-a,b,c]//ToCanonical

R[-b,-c,-d,-a]//ToCanonical

The function productR constructs random products of Riemann tensors. The second argument gives the list of free indices of the expression (always with even length). There is no defined conversion to Ricci or RicciScalar tensors.

canonindices[npairs_,frees_]:=Join[frees,Flatten[{#,-#}&/@GetIndicesOfVBundle[TangentM3,npairs,Join[frees,-frees]]]];
productR[n_,frees_:{}]:=Apply[Times,Apply[R,Partition[canonindices[2n-Length[frees]/2,frees][[First@RandomPerm[4n]]],4],1]]/;EvenQ[Length[frees]];

productR[1]

productR[10]

productR[3,{a,b,c,-d}]

productR[3,{a,b,c}]

This second function canonicalizes the products of Riemanns, timing the process:

canR[n_,frees_:{}]:=Flatten[{n,With[{product=productR[n,frees]},AbsoluteTiming[ToCanonical[product]]]}]

Examples:

canR[1]

canR[2]

canR[3]

canR[4]

canR[10]

canR[20]

canR[25]

Scaling of timings. Evaluating next cell could take between half a minute and a minute:

data1=canR/@Range[1,25];
data2=canR/@Range[1,25];
data3=canR/@Range[1,25];

plot0=ListLogPlot[Join[data1,data2,data3][[All,{1,2}]],PlotRange->All,Frame->True,PlotStyle->{PointSize[0.02],Hue[0.7]}]

We repeat the computation, but now leaving two and four free indices in the expressions:

data1=canR[#,{a,b}]&/@Range[1,25];
data2=canR[#,{a,b}]&/@Range[1,25];
data3=canR[#,{a,b}]&/@Range[1,25];

plot1=ListLogPlot[Join[data1,data2,data3][[All,{1,2}]],PlotRange->All,Frame->True,PlotStyle->{PointSize[0.02],Hue[0]}]

data1=canR[#,{a,b,c,d}]&/@Range[1,25];
data2=canR[#,{a,b,c,d}]&/@Range[1,25];
data3=canR[#,{a,b,c,d}]&/@Range[1,25];

plot2=ListLogPlot[Join[data1,data2,data3][[All,{1,2}]],PlotRange->All,Frame->True,PlotStyle->{PointSize[0.02],Hue[0.3]}]

Combining the three plots, we do not see any important difference in the timings:

Show[plot0,plot1,plot2]

Example 3:

Now let us concentrate on the products of seven Riemanns.

This figure (not constructed in this notebook) shows an histogram of the timings of canonicalization of a hundred thousand ramdom products of seven Riemanns. Note in the smaller plot that there are some cases which take much more time (and much more memory) to canonicalize. This is a consequence of the algorithm not being polynomial, but exponential in the worst case. Those are cases in which the symmetry is much higher than usual, for example because the product of Riemanns can be separated in products of decoupled monomials. There are cases in which, however, the reason is not at all clear, so that in general it is not possible to prepare the canonicalization algorithms to avoid all possible "worst cases".

Clean up

UndefTensor[R]
Remove[canonindices,productR,canR,data1,data2,data3,plot0,plot1,plot2]

MetricsOfVBundle[TangentM3]^={};
Clear[metricg];

Summarizing, the general canonicalization algorithms of xTensor` are highly efficient in most cases, much faster than the algorithms used by other tensor computer algebra systems. However, from time to time you will find special cases which take several orders of magnitude longer to canonicalize than other similar cases (for example a simple reordering of the indices keeping the same tensors), usually also swallowing up a large chunk of memory in your computer. This is normal behaviour, and it is due to the fact that the Group Intersection algorithm unavoidably behaves like that. Fortunately the special cases are always that: special; if you understand the meaning of "special" in your problem, you can prepare the system in advance to switch to "special" alternative algorithms.
In hindsight, this is the price we pay for using a general algorithm, but for the case of permutation symmetries it is worth paying it, because only in a very small fraction of cases the algorithm behaves exponentially, instead of polynomially. The case with multiterm symmetries is precisely the opposite: all known algorithms behave exponentially in most cases.

4.4. Canonicalization in a 1-dimensional space

xTensor` deals with 1-dimensional spaces in exactly the same way as with any other space. In particular, indices cannot be repeated, even though there is just one indiependent direction. This is solved by imposing full symmetry at canonicalization time.

Define a 1-dimensional space:

DefManifold[M1,1,{a1,b1,c1,d1}]

and some tensors in this space :

DefTensor[W[a1,b1],M1]

DefTensor[R[a1,b1],M1]

Now something like this is zero :

W[a1,b1]R[c1,d1]-W[d1,a1]R[c1,b1]//ToCanonical

Clean up :

UndefTensor/@{W,R};
UndefManifold[M1]

5. Rules and patterns

5.1. Unique dummies and screening

One of the most important issues in a tensor packages is how to deal with index contractions, i.e. with dummy indices. Until now this was not a real issue for us because all input expressions were simply canonicalized, a process which does not introduce new indices, but rather, eliminates indices. Now we face the opposite problem: expanding expressions into new expressions with more dummy indices, for example the conversion of a covariant derivative of a tensor in terms of a different covariant derivative. There are three problems with conversions that we must address.
This subsection explains in detail what those three problems are and how they will be solved. Though the explanations might seem a bit technical, it is important to understand them, because they are based on core concepts of Mathematica.
As an aside, there are two ways to define conversions of expressions into new expressions in Mathematica: assignments (Set and its variants SetDelayed, TagSet, UpSet, etc.) and rules (Rule and its variant RuleDelayed). Asignments can be considered as global or automatic rules. Rules can be considered as local assignments. See the Mathematica Help.
First problem. Let us see a simple example using Set:

Define the vector w in terms of S and v. Introduce a simple dummy:

DefTensor[w[a],M3]

w[a_]=S[a,b]v[-b]

[Note that we have used a Pattern on the right hand side. That is normal Mathematica, and a helpful feature which allows us to indicate the generality of the rule (which indices or tensors are accepted in the rule)].

Then these expressions would be totally wrong:

w[b]

w[a]w[-a]

Mathematica recommends the following solution to this problem:

Use a Module construct, such that the dummy index b is replaced by a unique symbol denoted with a $ sign, and a unique integer number. Note that we must now use SetDelayed (:=), and not Set (=), to avoid inmediate evaluation of Module:

w[a_]:=Module[{b},S[a,b]v[-b]]

w[b]

w[a]w[-a]

This solves the first problem but generates two more, which are solved by "screening" the dollar-indices:

Screening of indices.

First, each time we use the assignment we get a different object, what complicates comparisons. The following is not automatically recognized as zero!

w[a]-w[a]

%//InputForm

The second problem is purely aesthetical: the dollar-indices are terribly ugly. xTensor` has the function ScreenDollarIndices to convert the dollar-indices into normal-looking indices:

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

ScreenDollarIndices[%]

InputForm[%]

The Module construct is not the solution for all problems in tensor conversions. There are two more problems, as seen in the following examples. Second problem:

Objects and indices cannot share the same symbol names. Suppose there were a tensor also called a:

w[a_]:=a[a]

w[b]

In xTensor`, this is solved by enforcing type declarations. Every object must be defined supplying a symbol, and that symbol is given a type. A symbol cannot have two different types. This can be seen as a restriction, but it is actually a useful safety feature.

The tensor a[b] cannot be defined:

DefTensor[a[b],M3]

The final (third) problem involves dummies again, and happens in the few cases in which the same symbol represents a dummy on both sides and there are additional dummies on the right. The only possible solution is changing the name of the conflicting dummies either on the left or on the right hand side.

Define the following relation with dummies both on the left and right hand sides. Note that the b dummy is present in both sides and is a pattern on the left hand side (Mathematica already points out the error by reddening the dummy):

w[a_,b_,-b_]:=Module[{b,c},S[a,b,-b,c,-c]]

w[a,b,-b]

w[a,c,-c]

For example, we can change the name of the dummy index b on the right hand side:

w[a_,b_,-b_]:=Module[{d,c},S[a,d,-d,c,-c]]

w[a,b,-b]

w[a,c,-c]

As we said, the second problem is solved in xTensor` simply by enforcing type declarations. The first and third problems will be solved introducing a new collection of Index* assignment and rule functions which properly prepare the required Module constructs. They are IndexSet, IndexSetDelayed, IndexRule and IndexRuleDelayed, imitating their respective Mathematica counterparts, and they are described in detail in the next subsections.

We finish this section coming back to the function ScreenDollarIndices. We recommend to automate its use by assigning one of the global variables $Post or $PrePrint to it. There is a small difference between these two global variables (mentioned below), and it is not clear that one of them is always better than the other for screening, but I generally prefer $PrePrint:

Assign $PrePrint:

$PrePrint=ScreenDollarIndices;

Now the screening is automatic:

w[a,c,-c]

The screening happened after the result was assigned to Out[...] and hence if we temporarily remove $PrePrint, we still see the dollar-indices. (That wouldn't happen using $Post.)

$PrePrint=.

%%

Clean up

w[a_,b_,-b_]=.

From now on, in this notebook, we shall use automatic screening:

$PrePrint=ScreenDollarIndices;

5.2. IndexSet and IndexSetDelayed

The xTensor` functions IndexSet and IndexSetDelayed imitate the behaviour of Set and SetDelayed (which delays the evaluation of the rhs to the moment when the lhs is converted into the rhs).

Set functions for indexed expressions.

Example of use of IndexSet:

IndexSet[w[a_],S[a,b]v[-b]]

?w

w[b]

In this case the code finds the conflicting dummies, and then replaces the offending dummy on the right hand side:

DefTensor[T[a,b,c],M3]

IndexSet[T[a_,b_,-b_],S[a,b,-b,c,-c]]

?T

T[a,c,-c]

Example of use of IndexSetDelayed. Note that we only use a pattern in the first index of U and that w is still not evaluated:

IndexSetDelayed[T[-a_,-b,-c],w[-a]w[-b]w[-c]]

?T

T[-a,-b,-c]

We can change the definition of w and the conversion from S to w will be still valid:

IndexSet[w[a_],S[a,b]S[-b,-c]v[c]]

T[-a,-b,-c]

Clean up:

T[a_,b_,-b_]=.

T[-a_,-b,-c]=.

w[a_]=.

5.3. Index patterns

In the previous definitions of the vector w[a] we always used the pattern a_ which stands for any possible input, including up-indices and down-indices of any manifold, or even directional or label (or basis or component) indices. Sometimes this is what we want, but often it is not. In the following definitions there are no dummies or conflicts between indices and tensor names, so that we can safely use the usual SetDelayed (:=) definitions to simplify the examples.
The notation for tensors and indices in xTensor` is very simple, but that complicates working with patterns. On the other hand, there are always several ways to represent the same pattern.

Define, for the sake of simplicity:

dirvup=Dir[v[d]];
dirvdown=Dir[v[-d]];

This pattern stands for everything. Every w is converted into T:

w[a_]:=S[a]

{w[a],w[-a],w[dirvdown],w[dirvup],w[A],w[-A]}

w[a_]=.

Let us first work with the character of the indices: use the functions

Index character.

This pattern matches only up-indices:

w[a_?UpIndexQ]:=S[a]

{w[a],w[-a],w[dirvdown],w[dirvup],w[A],w[-A]}

w[a_?UpIndexQ]=.

This pattern matches only down-indices:

w[a_?DownIndexQ]:=S[a]

{w[a],w[-a],w[dirvdown],w[dirvup],w[A],w[-A]}

w[a_?DownIndexQ]=.

And now we work with the type of index: use the functions

Detect types of indices.

Introduce another manifold, with a new set of indices:

DefManifold[S2,2,{A,B,C,D,F,G}]

This pattern matches only abstract indices, from any vector bundle:

w[a_?AIndexQ]:=S[a]

{w[a],w[-a],w[dirvup],w[dirvdown],w[A],w[-A]}

w[a_?AIndexQ]=.

There is a simple and efficient pattern for abstract up-indices or abstract down-indices:

w[a_Symbol]:=S[a]

{w[a],w[-a],w[dirvup],w[dirvdown],w[A],w[-A]}

w[a_Symbol]=.

w[-a_Symbol]:=S[-a]

{w[a],w[-a],w[dirvup],w[dirvdown],w[A],w[-A]}

w[-a_Symbol]=.

This pattern matches only directional indices. Recall that directions do not have a sign in front. A message is sent complaining about the character of the vector v.

w[a_?DIndexQ]:=S[a]

{w[a],w[-a],w[dirvup],w[dirvdown],w[A],w[-A]}

w[a_?DIndexQ]=.

There is also a simpler version (which does not find the problem with the position of the ultraindex):

w[a_Dir]:=S[a]

{w[a],w[-a],w[dirvup],w[dirvdown],w[A],w[-A]}

w[a_Dir]=.

The functions AIndexQ, BIndexQ, CIndexQ, DIndexQ and GIndexQ admit a second argument restricting the manifold to which the indices must belong.

This pattern only matches abstract indices on TangentM3:

w[a_?(AIndexQ[#,TangentM3]&)]:=S[a]

{w[a],w[-a],w[dirvup],w[dirvdown],w[A],w[-A]}

w[a_?(AIndexQ[#,TangentM3]&)]=.

Again there are simpler and cleaner versions, also discriminating between upindices and downindices (Q-functions) or not (pmQ-functions). The use of Symbol is not needed now.

w[a_?TangentM3`Q]:=S[a]

{w[a],w[-a],w[dirvup],w[dirvdown],w[A],w[-A]}

w[a_?TangentM3`Q]=.

w[-a_?TangentM3`Q]:=S[-a]

{w[a],w[-a],w[dirvup],w[dirvdown],w[A],w[-A]}

w[-a_?TangentM3`Q]=.

w[a_?TangentM3`pmQ]:=S[a]

{w[a],w[-a],w[dirvup],w[dirvdown],w[A],w[-A]}

w[a_?TangentM3`pmQ]=.

In case of doubt use the function PatternIndex, which constructs patterns of the required form.

These are examples of different patterns for abstract indices named a. Note the minus sign in front of a for down-indices. The pattern is for the symbol of the index, and not for the whole index.

PatternIndex[a,AIndex]

PatternIndex[a,AIndex,Up]

PatternIndex[a,AIndex,Down,TangentM3]

PatternIndex[a,AIndex,Null,TangentM3]

These are patterns for directional indices named d:

PatternIndex[d,DIndex]

PatternIndex[d,DIndex,Up]

PatternIndex[d,DIndex,Down,TangentM3]

PatternIndex[d,DIndex,Null,TangentM3]

5.4. IndexRule and IndexRuleDelayed

The function IndexRule relates to Rule exactly in the same way that IndexSet relates to Set. It just introduces Module constructs. The symbol ↦ has been introduced. It is input as ∖[RightTeeArrow], but this input sequence doesn't work with the style used in this notebook (I don't know why).

Rule  functions for indexed expressions.

Here we have again the problem of repeated indices b:

{w[a],w[b]}/.w[a_]->S[a,b]v[-b]

Using IndexRule the problem disappears:

{w[a],w[b]}/.IndexRule[w[a_],S[a,b]v[-b]]

{w[a],w[b]}/.w[a_]↦S[a,b]v[-b]

5.5. MakeRule

Apart from the coding aspects solved by the Index* functions, there are other, more mathematical issues which must be adressed. For instance, given a rule LHS->RHS how can we apply the rule on expressions which are equivalent to LHS but indices are permuted due to symmetries or metric shifting? Can we apply the rule on all vbundles or just on some of them? These kind of options are controlled by MakeRule, the second most complicated function in xTensor` after ToCanonical. This function uses the properties of the metric tensors, and hence we delay its description to subsection 9.1 below.

6. Derivatives

6.1. Covariant derivatives

We understand a covariant derivative as a covariant derivative operator (see Wald). In particular, ordinary ("partial") derivative operators are obtained as covariant derivatives with the options Curvature->False and Torsion->False. Covariant derivatives are real operators acting on (possibly complex) vector bundles. The treatment of derivatives on xTensor` is fully abstract, and hence we do not explain here how to take components of a derivative of a tensor. This is dealt with in the companion package xCoba`, and explained in its documentation.

Definition of a covariant derivative.

Default options for DefCovD:

Options[DefCovD]

Define a covariant derivative Cd on TangentM3 (the default vbundle obtained from the index given in the derivative). Automatically the Christoffel symbols and the Riemann and Ricci tensors are defined. By default xTensor` does not assume that this derivative is associated to a metric and therefore the symmetries of the Riemann tensor are not as expected: there is just antisymmetry in the first pair of indices. The Ricci tensor has no symmetry at all. There is no Ricci scalar.

DefCovD[Cd[-a]]

SymmetryGroupOfTensor[RiemannCd]

RiemannCd[-a,-b,-c,b]

SymmetryGroupOfTensor[RicciCd]

The names of the associated tensors are determined by the function GiveSymbol. Their output symbols are determined by the function GiveOutputString.If you want to change any of those, modify the corresponding function before defining the object.

GiveSymbol[Riemann,Cd]

GiveOutputString[Riemann,Cd]

Derivatives can be represented in two ways, encoded as the strings "Prefix" or "Postfix" in the global variable $CovDFormat. The default value is "Prefix". Directional derivatives are always given in "Prefix" notation.

$CovDFormat

Cd[-a][T[a,b,-c]]

$CovDFormat="Postfix";

Cd[-a][T[a,b,-c]]

Cd[Dir[v[d]]][T[a,b,-c]]

Derivatives are linear operators and implement automatically the Leibnitz rule. The final ordering of terms is decided by Mathematica.

DefTensor[r[],M3]

Cd[-a][7r[]v[a]+S[a,b]v[-b]]

Dummy indices are, by default, pushed to the right by the canonicalization process

Cd[-a][T[f,c,-b]U[-c,-e,-d]]//ToCanonical

Commutation of derivatives.

By default, a special ordinary derivative called PD is defined. It is one of the many covariant derivative operators without curvature and without torsion, but no further assumptions are made about it. It can be used as a generic "partial derivative" in most computations, but its meaning is different. Note that they are not automatically sorted:

PD[-a][PD[-b][T[c,d,-e]]]

$CovDFormat="Prefix";

PD[-a][PD[-b][T[c,d,-e]]]

PD[-a][PD[-b][T[c,d,-e]]]-PD[-b][PD[-a][T[c,d,-e]]]

%//SortCovDs

We can automate the commutation of derivatives:

SortCovDsStart[PD]

PD[-a][PD[-b][T[c,d,-e]]]-PD[-b][PD[-a][T[c,d,-e]]]

SortCovDsStop[PD]

All symmetric (torsionless) connections commute (with themselves) on scalar expressions. That can be avoided switching off the variable $CommuteCovDsOnScalars.

The input of high-order derivatives is certainly slow due to the sharp syntactic restrictions of xTensor`. There are several options to avoid this problem. Of course, you cannot use Validate unless you modify it to recognize the new syntactic extensions.

The standard way to input a third order derivative is this:

PD[-a][PD[-b][PD[-c][T[c,d,-e]]]]

We can use the prefix notation of Mathematica:

PD[-a]@PD[-b]@PD[-c]@T[c,d,-e]

or we can define our own rules to input derivatives, breaking the standard syntax. The simplest possibility, based on the fact that the first bracket of PD always contains a single index, is:

Unprotect[PD];

PD[expr_,a__,b_]:=PD[b][PD[expr,a]]
PD[expr_,a_]:=PD[a][expr]

PD[T[c,d,-e],-c,-b,-a]

InputForm[%]

PD[expr_,a__,b_]=.
PD[expr_,a_]=.

or an intermediate situation, based on the same fact, is:

PD[a_,b__][expr_]:=PD[a][PD[b][expr]]

PD[-a,-b,-c][T[c,d,e]]

InputForm[%]

PD[a_,b__][expr_]=.

A totally different notation could be

PD/:tensor_?xTensorQ[indices___,PD[index_]]:=PD[index][tensor[indices]]

T[c,d,-e,PD[-c],PD[-b],PD[-a]]

InputForm[%]

PD/:tensor_?xTensorQ[indices___,PD[index_]]=.

Protect[PD];

6.2A. Implosion and explosion (I)

It is important to understand the following notational problem with derivatives. From a mathematical point of view, the derivative of a function or a tensor is another function or tensor, and hence should share its same notation. However, we must input the derivative using an operator, and hence a notational change should happen at some point. For instance, in Mathematica the operator D[f[x], x] is inmediately converted into Derivative[1][f][x]. The former doesn't have the "function form" head[x] but the latter does, with the head being Derivative[1][f]. There are many obvious advantages in this change, because now we can use a consistent treatment for all functions, but there are also disadvantages.

The main disadvantage of the change is, a mystery for all new users of Mathematica, that this does not work:

D[f[x],x]/.f[x]->Sin[x]

but this does work:

D[f[x],x]/.f->Sin

In xTensor` we work essentially with the operator form only, and hence rules are simpler to use, but we have the problem that tensors and their covariant derivatives are treated differently. To alleviate that problem, a new device has been introduced in version 0.9.6: we now have the concepts of "implosion" and "explosion":

Suppose we start from the covariant derivative

Cd[-a][v[b]]

InputForm[%]

Now we "implode" this expression, resulting in a new tensor being defined (there is no message reporting the creation of this new tensor, but that can be activated setting the global variable $ImplodeInfo = Info). Compare with the previous expression, and please pay attention to the positions of indices!

Implode[%]

InputForm[%]

The expression can be "exploded", returning to the original:

Explode[%]

Implosion involves the construction of new tensors, and in this process their symmetries must also be constructed.

This object is symmetric:

PD[-a]@PD[-b]@v[c]

SymmetryOf[%]

If we implode it, that symmetry must be stored into the new tensor:

Implode[%%]

SymmetryGroupOfTensor[PDPDv]

Other torsionless derivatives are commuted on scalars :

Cd[-a]@Cd[-b]@r[]

Implode[%]

SymmetryGroupOfTensor[CdCdr]

but not on other tensors :

Cd[-a]@Cd[-b]@v[a]

Implode[%]

SymmetryGroupOfTensor[CdCdv]

The problem of implosion becomes more complicated when metrics (and bases) are present, but we defer the corresponding analysis to the next section, after metrics are introduced.

6.2B. TensorDerivative

Now that we can have tensors with non-atomic heads, it is natural to absorb the imploded information in the head of a new tensor. Following Derivative in Mathematica we introduce the head TensorDerivative.

Take a tensor expression containing derivatives:

Cd[-a][v[b]]

Convert it into the new notation (note that the derivative index is placed last):

%//ToTensorDerivative

%//InputForm

and back:

%%//FromTensorDerivative

%//InputForm

6.3. Change of covariant derivative

In xTensor` we can work simultaneously with any number of covariant derivatives.

The list of all covariant derivatives currently defined is

$CovDs

We define another one

DefCovD[CD[-a],SymbolOfCovD->{"|","D"}]

The difference between two covariant derivatives of the same tensor can be expressed in terms of Christoffel tensors (the C tensors of Wald), and in xTensor` we shall always use this point of view of Christoffels: they are always tensors, but associated to two covariant derivatives. The Christoffel tensor defined together with each derivative is the tensor associated to that derivative and the derivative PD associated to a generic chart.

Christoffel constructs the Christoffel tensor relating two derivatives. It is antisymmetric in its derivatives arguments.

Christoffel[CD,PD][a,-b,-c]

%//InputForm

Christoffel[PD,CD][a,-b,-c]

%//InputForm

The Christoffel tensor relating two non-PD derivatives is constructed automatically whenever needed, with (lexicographically) sorted derivatives, to avoid duplicity:

Christoffel[CD,Cd][a,-b,-c]

%//InputForm

The origin of this curious conversion from Christoffel[Cd,CD] to ChristoffelCdCD is again the fact that we cannot associate information to the former.

Using this structure of Christoffel tensors we can relate any two derivatives of any tensor.

Suppose this derivative:

expr=CD[-d][U[-a,b,-c]]

ChangeCovD (aka CovDToChristoffel in xTensor` version 0.7) changes by default to PD:

ChangeCovD[expr]

That is equivalent to

ChangeCovD[expr,CD]

but we can also change to any desired derivative:

ChangeCovD[expr,CD,Cd]

ChangeCovD[%,Cd,CD]

If you do not like the double-derivative Christoffels, you can always break them to the CovD-PD Christoffels:

%%//BreakChristoffel

The expansion is recursive:

ChangeCovD[Cd[-e][Cd[-a][T[f,c,-b]]],Cd]

Change of covariant derivative.

If the relation between two covariant derivatives is fully described by a Christoffel tensor, then the curvature and torsion tensors associated to them must be also related by those Christoffel tensors.

This is the curvature tensor of the derivative CD:

RiemannCD[-a,-b,-c,d]

ChangeCurvature (aka RiemannToChristoffel in xTensor` version 0.7) changes any curvature tensor of a derivative to the curvature tensor of other derivative (by default PD, with zero curvature).

ChangeCurvature[%]

but we can convert to any derivative:

ChangeCurvature[%%,CD,Cd]

ChangeCurvature[RiemannCD[-a,-b,-c,a],CD,Cd]

Define a covariant derivative with torsion, but not metric-compatible. Now the Christoffel tensor is non-symmetric:

DefCovD[CDT[-a],SymbolOfCovD->{"#","DT"},Torsion->True]

Now the formulas involving CDT will contain torsion terms:

ChangeCurvature[RiemannCd[-a,-b,-c,a],Cd,CDT]

%//InputForm

The difference between the torsion tensors of two derivatives is given by the antisymmetric part of the Christoffel relating them. In this case the torsion of Cd is zero. The change is performed by ChangeTorsion (aka TorsionToChristoffel in xTensor` version 0.7):

ChangeTorsion[TorsionCDT[a,-b,-c],CDT,Cd]

Induced changes in curvature and torsion tensors.

Undefine some derivatives:

UndefCovD/@{CD,CDT};

Note that xTensor` also allows metric-compatible connections with torsion. The symmetry properties of the associated curvature tensors are not complete. We shall later illustrate this type of derivatives, after defining metric fields.

6.4. General Bianchi identities

Let us now check the Bianchi identities for general covariant derivatives with torsion. Note that these identities are not directly encoded in xTensor`, but can be easily computed.

Define a covariant derivative with torsion. Note the special symmetry properties of the defined tensors:

DefCovD[CD[-a],Torsion->True]

This is the general form of the first Bianchi identity:

RiemannCD[-a,-b,-c,d]+CD[-a][ TorsionCD[d,-b,-c] ]-TorsionCD[e,-a,-b]TorsionCD[d,-c,-e]

6Antisymmetrize[%,{-a,-b,-c}]

The computation can be performed by transforming the derivative CD and its associated tensors into PD and Christoffel tensors:

%//ChangeCovD

%//RiemannToChristoffel

%//TorsionToChristoffel

%//Expand

This is the general form of the second Bianchi identity:

CD[-a][ RiemannCD[-b,-c,-d,e] ]-TorsionCD[f,-a,-b]RiemannCD[-c,-f,-d,e]

6Antisymmetrize[%,{-a,-b,-c}]

The computation proceeds along the same lines:

%//ChangeCovD

%//RiemannToChristoffel

%//TorsionToChristoffel

%//Expand

%//ToCanonical

Clean up:

UndefCovD[CD]

6.5. Constant symbols

It is important to discuss the concept of a constant in xTensor`. A scalar field c[] is a constant if it is defined on no manifold. However a nonscalar tensor field cannot be constant on the base manifold of its indices (actually xTensor` automatically defines the tensor to be a field on the base manifolds on its indices). On the other hand we have to use symbols denoting constants (as Newton's constant). We shall use the command DefConstantSymbol to define those constants.

Functions that operate with constants.

This defines a constant scalar field:

DefTensor[const[],{}]

?const

By default, it is not automatic to check that a derivative on an object living on a different manifold is zero because that is a slow process that could be internally happening too often.

Cd[-a][const[]]

CheckZeroDerivative[%]

We can automate it for a particular derivative (or for all):

CheckZeroDerivativeStart[PD]

PD[-a][const[]]

CheckZeroDerivativeStop[PD]

An arbitrary symbol is not understood by default as a constant. We must define it as a constant:

Cd[-a][GNewton]

DefConstantSymbol[GNewton]

Cd[-a][GNewton]

UndefConstantSymbol[GNewton]

PD[-a][GNewton]

Remove[GNewton]

6.6. Parameter derivatives

Given a parameter defined with DefParameter we can take parametric derivatives of expressions which depend on that parameter, using ParamD. For historical reasons xTensor` has a generic parametric derivative called OverDot such that every field is considered to depend on the corresponding (undefined) parameter.

Parametric derivatives.

The generic parametric derivative OverDot (Dot is used by Mathematica as inner scalar product) is represented by a (very) small over dot:

OverDot[r[]^2]

OverDot has been defined to commute with partial derivatives, but not with other covariant derivatives or with Lie derivatives:

OverDot[PD[-a][r[]]]//InputForm

OverDot[Cd[-a][r[]]]//InputForm

All fields currently defined do not depend on external parameters, so that we define new objects:

DefParameter[time,PrintAs->"t"]

ParamD[time][r[]]

DefTensor[Q[-a],{M3,time}]

?Q

V1.0 UPDATE: From version 1.0.2 there are several posible representations of a ParamD derivative, controlled by a global variable called $ParamDFormat, with default value "Prefix". Let us first simplify the printing of the parameter:

$ParamDFormat

ParamD[time][3Q[-a]v[a]]

ParamD[time,time][3Q[-a]v[a]]

$ParamDFormat="SinglePrefix";

ParamD[time][3Q[-a]v[a]]

ParamD[time,time][3Q[-a]v[a]]

$ParamDFormat="Postfix";

ParamD[time][3Q[-a]v[a]]

ParamD[time,time][3Q[-a]v[a]]

Clean up:

UndefTensor[Q]

UndefParameter[time]

6.7. Lie derivatives

xTensor` does not define independent Lie derivatives, but a single command LieD which depends both on a contravariant vector field and on the expression that we want to derive.

Lie derivatives.

Lie derivatives of expressions with respect to a contravariant vector field are not automatically expanded in terms of derivatives.

LieD[v[a]][S[-b,c]]

However, we can expand them using any derivative (PD is the default).

LieDToCovD[LieD[w[a]][S[-b,c]],PD]

LieDToCovD[LieD[w[a]][S[-b,c]],Cd]

The derivative can also be expanded directly:

LieD[w[a],Cd][S[-b,c]]

LieD knows the Leibniz rule:

LieD[w[a]][S[a]T[-a] ]

We can have any expression in the "directional argument". Recall however that it is not a tensorial slot (it hence cannot be described with an abstract index):

LieD[7r[]w[a]+v[a]][T[a]]

6.8. Directional derivatives

From the mathematical point of view, the simplest concept of a derivative is that of a "derivation", a mapping from scalar fields to scalar fields. In fact, the concept of tanget space at a point is constructed as the space of derivations at that point. We can work with operators acting as derivations (or "directional derivatives") using the Dir head.

The object PD[Dir[v[a]]] can be interpreted as the directional derivative along the vector v[a]:

PD[Dir[w[a]]][r[]]

and it has all the expected properties:

PD[Dir[ 3S[a,Dir[v[-b]]]+w[a] ]][ 1/r[]^2 ]

Actually, it is possible to work directly with the derivations acting on scalars, saving 6 brackets (!) in the notation. For example we define the derivation wder associated to the vector w[a]:

wder:=PD[Dir[w[a]]]

wder[ r[] ]

Or even better, in prefix notation:

wder@r[]

ScalarQ[%]

Remove[wder]

6.9. Lie brackets

The Lie bracket of two vector fields is a third vector field.

Using abstract indices a bracket can be expressed as

Bracket[v[b],w[b]][a]

Note that the indices of the vector fields are ultraindices. The index of the resulting vector field is represented outside the bracket.

It has almost all expected properties: bilinearity, antisymmetry, derivation, ..., but not yet the Jacobi rule.

Bracket[v[b]+w[b],3w[b]+r[]v[b]][a]

UndefTensor[w];

6.10. Derivatives on inner vbundles

It is possible to define a derivative which acts on an inner vector bundle (as used in Yang-Mills gauge theories).

Define a complex vector bundle:

DefVBundle[InnerC,M3,4,{𝔄,𝔅,ℭ,𝔇,𝔈,𝔉,𝔊,ℌ},Dagger->Complex]

Define a connection on a complex inner vbundle:

DefCovD[ICD[-a],InnerC,SymbolOfCovD->{"*","D"}]

We see that a single torsion tensor has been defined, but two Christoffel and two Riemann tensors. The inner Christoffel is called AChristoffelICD and the inner Riemann tensor is called FRiemannICD. They are both complex tensors and therefore have their complex conjugates. Note the different vbundles of their indices.

?AChristoffelICD

?FRiemannICD

?AChristoffelICD†

?FRiemannICD†

The derivative ICD has curvature in both the tangent bundle TangentM3 and in the inner bundle InnerC. That becomes most apparent when changing to a different derivative:

DefTensor[X[a,𝔅],M3,Dagger->Complex]

ICD[-c][X[a,𝔅]]

ChangeCovD[%,ICD]

ChangeCovD[%,PD,ICD]

or when commuting two derivatives:

ICD[-c]@ICD[-d]@X[a,𝔅]

SortCovDs[%]

For their complex conjugates:

ICD[-c][X[a,𝔅]]//Dagger

ChangeCovD[%,ICD]

ICD[-c]@ICD[-d]@X[a,𝔅]//Dagger

SortCovDs[%]

UndefTensor[X]

An important concept, introduced in version 0.9.1 of xTensor`, is that of extension of a derivative.  Given a derivative Cd it is possible to define a second one CD in such a way that the associated tensors of CD on its tangent bundle will be those of Cd. We say then that CD extends Cd, or that CD has been constructed by extension of Cd.

Let us take the Cd derivative as base covd:

$CovDs

We now extend that derivative to CD:

DefCovD[CD[-a],InnerC,SymbolOfCovD->{"#","D"},ExtendedFrom->Cd]

Now the tensors associated to the extended derivative will be those of the base derivative:

{ChristoffelCD[a,-b,-c],RiemannCD[-a,-b,-c,d],TorsionCD[a,-b,-c]}

And then we have the inner curvature of the extended derivative:

{AChristoffelCD[𝔄,-b,-ℭ],FRiemannCD[-a,-b,-ℭ,𝔇]}

We can undefine the extended derivative, without removing any of the base derivative objects:

UndefCovD[CD]

{ChristoffelCd[a,-b,-c],TorsionCd[a,-b,-c],RiemannCd[-a,-b,-c,d],RicciCd[-a,-b]}

{ChristoffelCD[a,-b,-c],TorsionCD[a,-b,-c],RiemannCD[-a,-b,-c,d],RicciCD[-a,-b]}

7. Metrics

7.1. Definition

Now we introduce metrics on the vbundles.

Definition of a metric.

We define a metric metricg with negative determinant. The associated covariant derivative will be called CD and will be represented with a semicolon under IndexForm. All the needed associated tensors are defined (note that epsilon has a suffix metricg and the other tensors have a suffix CD). Because we use indices of 𝕋M it is understood that metricg will be a metric on M3.

DefMetric[-1,metricg[-a,-b],CD,PrintAs->"g"]

Now we have

RiemannCD[-a,b,-b,-c]

%//InputForm

RiemannCD[-a,-b,b,a]

%//InputForm

CD[-c][metricg[-a,-b]]

Cd[-c][metricg[-a,-b]]

CD[-c][epsilonmetricg[a,b,c]]

CD[-a][EinsteinCD[b,c]]

CD[-a][EinsteinCD[a,b]]

The determinant of the basis in the default basis AIndex (the double tilde denotes weight +2 in that basis; this is Astekhar notation):

Detmetricg[]

Its Weyl tensor is zero because the manifold is 3d:

WeylCD[-a,-b,-c,-d]

There is also the traceless part of the Ricci tensor:

TFRicciCD[-a,-b]

TFRicciCD[a,-a]

Now xTensor` accepts index structures that do not correspond to the original definition, as long as they only involve defined metrics:

SlotsOfTensor[U]

Validate[U[a,b,c]]

DefTensor[X[a,b,A],{M3,S2}]

SlotsOfTensor[X]

Validate[X[a,b,-A]]

UndefTensor[X]

Expression of the Christoffel symbols (with respect to PD) in terms of the metric:

CD[-a][v[-b]]//ChangeCovD

%//ChristoffelToGradMetric

7.2. Contraction

Metrics are not contracted by default because in general that obscures the manipulation of indices and the canonicalization process becomes slower. However, xTensor` has a number of commands that perform metric contraction. When the expression is a pure tensor contraction is a simple process, but when there are derivatives which are not associated to the (first) metric of the vbundle then problematic situations could develop. There are two situations we need to worry about: one is the process of contraction itself, shown in this subsection; the other is the influence of the metric in the process of canonicalization, studied in the next subsection.

Metric contraction.

Metric factors are automatically contracted with themselve:

metricg[a,b]metricg[-b,-c]

metricg[a,-a]

However they are only contracted with other objects if required. This is to have more control on expressions formed by metric factors, and mostly to avoid unexpected index-shiftings.

metricg[a,b]S[-a,c]

ContractMetric[%,metricg]

Derivatives of a metric are automatically converted so that the indices of the metric are covariant, introducing new dummies. Note that derivatives of mixed-indices metrics are zero. That means that all non-zero derivatives of metrics will have covariant indices:

PD[-a][metricg[b,c]]

PD[-a][metricg[b,-c]]

Contraction of derivative indices is not allowed by default, but can be forced using the option AllowUpperDerivatives:

ContractMetric[metricg[a,b]PD[-b][S[c,d]],metricg]

ContractMetric[metricg[a,b]PD[-b][S[c,d]],metricg,AllowUpperDerivatives->True]

With ContractMetric a metric can be contracted through its Levi-Civita derivative, but not through other derivatives:

ContractMetric[metricg[a,b]CD[-d][S[-a,c]],metricg]

ContractMetric[metricg[a,b]PD[-d][S[-a,c]],metricg]

In those cases we need to use a different function which forces contraction over derivatives, but at the price of producing new terms with derivatives of the metric:

ContractMetric[metricg[a,b]PD[-d][S[-a,c]],metricg,OverDerivatives->True]

7.3A. Canonicalization with a metric. IMPORTANT BUG

In the presence of a metric new symmetries can be added to the process of canonicalization. The command ToCanonical has an option UseMetricOnVBundle that specifies on which vbundles the metric (if there is one) must be used during canonicalization. Possible values are All (default), None or a list of (maybe one) vbundles.

Tensors S and T have no symmetries and therefore their canonicalization without metric is trivial:

ToCanonical[S[c,a]T[-a,-b,b],UseMetricOnVBundle->None]

However with metric the canonical form is different because pairs of dummies are sorted such that the up-index comes first:

ToCanonical[S[c,a]T[-a,-b,b]]

When derivatives are present canonicalization with a metric is harder, unless the only derivative is its Levi-Civita associated operator. In other cases, ToCanonical shows a warning message and changes to a different algorithm rather than using the default algorithm which just moves indices around and that could get easily confused. There are two possible internal algorithms, controlled with the option Method of ToCanonical. Before version 0.9.8 of xTensor` the only possibility was (and still is the default) Method->ChangeCovD, which changes the non-compatible derivatives into Levi-Civita connections. Now there is a second option, Method->Implode, which uses the new idea of converting all derivatives into tensors. The latter is faster and more general, but hasn't been tested much. Use it with care. Both methods considerably slow down the process of canonicalization because they require several internal calls to ToCanonical.

We use the same tensors S and T, with metric-compatible derivatives:

ToCanonical[S[c,a]CD[-d][T[-a,-b,b]],UseMetricOnVBundle->None]//AbsoluteTiming

ToCanonical[S[c,a]CD[-d][T[-a,-b,b]]]//AbsoluteTiming

But now if we use PD instead:

ToCanonical[S[c,a]PD[-d][T[-a,-b,b]],UseMetricOnVBundle->None]//AbsoluteTiming

ToCanonical[S[c,a]PD[-d][T[-a,-b,b]]]//AbsoluteTiming

This form of canonicalization leaves lots of Christoffel hanging around. You can get rid of them in three steps: first use ChristoffelToGradMetric; then use ContractMetric,

Last[%]//ChristoffelToGradMetric//ContractMetric

Finally, canonicalize again, but this time without moving indices, to avoid bringing back the Christoffels:

ToCanonical[%,UseMetricOnVBundle->None]

The whole process is automatic with the option ExpandChristoffel:

ToCanonical[S[c,a]PD[-d][T[-a,-b,b]],ExpandChristoffel->True]//AbsoluteTiming

The alternative algorithm based on implosion is similar (note that it would require a further use of ContractMetric to arrive at the same result):

ToCanonical[S[c,a]PD[-d][T[-a,-b,b]],Method->Implode]//AbsoluteTiming

This is another, somehow harder, example. Partial derivatives do not commute when having upper indices:

PD[-b]@PD[a]@v[-a]-PD[a]@PD[-b]@v[-a]//Simplification

Now we can drop the two second derivatives, and remove the Christoffels:

%//SortCovDs

result1=Simplification[%//ChristoffelToGradMetric//ContractMetric,UseMetricOnVBundle->None]//AbsoluteTiming

In this case the newer algorithm is faster, and again differs from the previous result just in the use of ContractMetric:

ToCanonical[PD[-b]@PD[a]@v[-a]-PD[a]@PD[-b]@v[-a],Method->Implode]//AbsoluteTiming

General recommendations:
	1. The simplest situation is having just one tangent bundle, with one metric and its Levi-Civita connection. In that case the metric can be always contracted, and we recommend to do it, just to decrease the number of indices.
	2. If there are several tangent bundles, each one having one metric with its Levi-Civita connection, then the metrics can be freely contracted, and the only additional problem is that you decide whether to use CheckZeroDerivative at some intermediate steps of the calculation, or else automate that point.
	3. If there is one vbundle, one metric with its Levi-Civita connection and an additional derivative (ordinary or not, that is irrelevant) then the canonicalizer must change to the slower algorithm. I recommend not to change the character of the slots of the tensors. To do that set UseMetricOnVBundle->False in ToCanonical and never use ContractMetric. The metric factors will be always explicit, and the number of indices in the expressions will be larger, but overall the canonicalizations will be faster.
	4. If there are several vbundles with several derivatives, combine comments 2 and 3.
	5. If furthermore there are several metrics on the vbundles, each with its own Levi-Civita connection, then before start crying, see next subsection.

The canonicalizer also takes into account whether indices can be moved up and down or not. This is done internally using the concept of metric-state for indices. Study next example carefully:

This expression is clearly antisymmetric under the exchange of derivatives:

U[a,b,c]PD[-a][v[d]]PD[-b][v[e]]metricg[-d,-e]

and it is thus zero:

ToCanonical[%]

We contract the metric through a derivative. The expression however is now different and not zero:

U[a,b,c]PD[-a][v[-d]]PD[-b][v[e]]metricg[d,-e]

but the system detects it and change to the alternative algorithm (THIS IS BROKEN NOW):

ToCanonical[%,Method->ChangeCovD,ExpandChristoffel->True]

ToCanonical[%%,Method->Implode]//ContractMetric

Another example:

T[a,-b,c]PD[-d][v[-a]]PD[-e][v[b]]metricg[d,e]

ToCanonical[%]

ToCanonical[%%,ExpandChristoffel->True]

ToCanonical[%%%,Method->Implode]

7.3B. Canonicalization with a metric using TensorDerivative

The head TensorDerivative allows an easier handling of this problem.

Take the same example we had before. We arrive at a result much faster:

PD[-b]@PD[a]@v[-a]-PD[a]@PD[-b]@v[-a]//ToTensorDerivative

%//ToCanonical

Other examples shown before:

U[a,b,c]PD[-a][v[d]]PD[-b][v[e]]metricg[-d,-e]

%//ToTensorDerivative

ToCanonical[%]

Another example:

U[a,b,c]PD[-a][v[-d]]PD[-b][v[e]]metricg[d,-e]

%//ToTensorDerivative

ToCanonical[%]

Another example:

T[a,-b,c]PD[-d][v[-a]]PD[-e][v[b]]metricg[d,e]

%//ToTensorDerivative

ToCanonical[%]

7.4. Several metrics

It is possible to work with several metrics on the same vbundle. Of course, not all of them can raise and lower the indices of other tensors because we would inmediately loose track of the metric that was used to move a given index. The situation in xTensor` is as follows. The list of metrics on a vbundle is stored as an upvalue of the vbundle for the function MetricsOfVBundle. The metrics are placed in that list in order of definition. The first metric is considered special and it is the only one which can move indices. All other metrics (called "frozen" metrics) are simply symmetric two-tensors with a Levi-Civita associated connection and a number of associated tensors. For a frozen metric, many of the expected rules do not work.
It is possible to remove the first metric, suddenly converting the second metric into the first of the list. I cannot think of any situation where that could be valid or just safe.
I have never seriously tested this part of the code. Please use it with maximum care!

Currently there is only one metric on TangentM3:

MetricsOfVBundle[TangentM3]

We define a second (frozen) metric:

DefMetric[1,frozen[-a,-b],CD2,SymbolOfCovD->{"|","D"},PrintAs->"f"]

MetricsOfVBundle[TangentM3]

Compare these two results:

{metricg[a,-b],frozen[a,-b]}

That is because these three expressions behave very differently:

{frozen[a,c]frozen[-c,-b],metricg[a,c]frozen[-c,-b],metricg[a,c]metricg[-c,-b]}

ContractMetric[%]

In particular the inverse of frozen[-a,-b] is not frozen[a,b], but Invfrozen[a,b]:

Invfrozen[a,b]frozen[-b,-c]

It is problematic to generalize typical expressions. For example this is not zero now:

CD2[-a][EinsteinCD2[a,b]]

but this is zero:

Invfrozen[a,b]CD2[-a][EinsteinCD2[-b,-c]]

V1.0 UPDATE: An important problem is that of the symmetries of the Riemann tensor. Only if we lower the upper index of the Riemann tensor with the appropriate metric we get the correct answer. But when that appropriate metric is not the first metric then the all-lower-indices Riemann tensor cannot be identified with the original tensor directly. Hence, we introduce the RiemannDown tensor.

Note the different symmetries of these two tensors:

SymmetryGroupOfTensor[RiemannCD2[-a,-b,-c,d]]

SymmetryGroupOfTensor[RiemannDownCD2[-a,-b,-c,-d]]

It is possible to convert between those two tensors (note the different formatted representation):

RiemannCD2[-a,-b,-c,d]

%//RiemannToRiemannDown

%//RiemannDownToRiemann

UndefMetric[frozen]

Other commands and algorithms have been updated in V1 .0 to work properly with frozen metrics.

7.5. Curvature tensors of the metric

When we define a metric, a number of curvature tensors are defined. Currently they are: Riemann, Ricci, RicciScalar, Einstein, TFRici, Weyl. Note that they are all associated to the Levi-Civita connection of the metric, and not directly to the metric. There are a number of functions to change among them.

First we have some automatic contractions. These can be prevented with the option CurvatureRelations of DefCovD. By default, contractions of the Riemann tensor are automatically replaced by the Ricci tensor:

RiemannCD[-a,-b,-c,a]

%//InputForm

And contractions of the Ricci tensor are replaced by the Ricci scalar:

metricg[c,d]RicciCD[-d,-c]

%//InputForm

RicciCD[a,-a]

Then we can change among the tensors:

EinsteinCD[-a,-b]//EinsteinToRicci

%//RicciToEinstein

In a 3d manifold the Weyl tensor is zero:

WeylCD[-a,-b,-c,-d]

RiemannCD[-a,-b,-c,-d]//RiemannToWeyl

We can also change between Ricci and TFRicci:

RicciCD[-a,-b]//RicciToTFRicci

%//TFRicciToRicci//ToCanonical

Relations among curvature objects.

As we said, we can work with connections which are compatible with a given metric field but have torsion too.

We define a new connection associated to our metric field, but this time with torsion. The Riemann tensor is antisymmetric in both pairs, but now those pairs cannot be exchanged. Hence the Ricci tensor is not symmetric. Currently Weyl is not even defined:

DefCovD[CDT[-a],SymbolOfCovD->{"#","D"},Torsion->True,FromMetric->metricg]

The difference of covariant derivatives

CDT[-a][v[b]]-CD[-a][v[b]]

is given by a Christoffel tensor

ChangeCovD[%,CDT,CD]

which can be expressed uniquely in terms of torsion and metric fields:

%//ChristoffelToGradMetric

UndefCovD[CDT]

7.6. Product metrics

Spacetimes with a high degree of symmetry can be sometimes locally decomposed as products of reduced manifolds and the orbits of symmetry. The metric of the whole manifold can then be given in terms of metrics on simpler vbundles. Currently there is a stupid limitation to two subvbundles. It will be dropped soon.

Define a product manifold:

DefManifold[M5,{M3,S2}, {μ, ν, λ, σ, η, ρ}]

Define a product metric. We already have a scalar r[] on M3. Now we define another scalar w[] on S2:

DefTensor[w[],S2]

DefProductMetric[g5[-μ,-ν],{{TangentM3,w[]},{TangentS2,r[]}},Cd5,SymbolOfCovD->{"#","𝔇"}]

DefMetric[1,gamma[-A,-B],CdS,SymbolOfCovD->{":","D"},PrintAs->"γ"]

We define the warped metric      w[]^2 MetricOfVBundle[TangentM3]+r[]^2 MetricOfVBundle[TangentS2] :

DefProductMetric[g5[-μ,-ν],{{TangentM3,w[]},{TangentS2,r[]}},Cd5,SymbolOfCovD->{"#","𝔇"}]

Computations with product metrics.

The delta tensors are automatically expanded:

{g5[a,-b],g5[A,-B]}

InputForm/@%

Now we can compute any object on M5 in terms of objects of M3 and S2:

ExpandProductMetric[{g5[-a,-b],g5[-a,-B],g5[-A,-B]},g5]

ExpandProductMetric[{ChristoffelCd5[A,-B,-C],ChristoffelCd5[a,-B,-C],ChristoffelCd5[A,b,-C],ChristoffelCd5[a,-b,-C],ChristoffelCd5[A,-b,-c],ChristoffelCd5[a,-b,-c]},g5]//ContractMetric

RiemannCd5[-a,-b,-c,-d]//ExpandProductMetric//ContractMetric//Simplify

RicciScalarCd5[]//ExpandProductMetric//ContractMetric//Simplification

By default derivative indices cannot be contracted with metric tensors. This behaviour can be changed using the option AllowUpperDerivatives:

ContractMetric[%,AllowUpperDerivatives->True]//Simplification

Clean up

UndefMetric[g5]

UndefMetric[metricg]

UndefMetric[gamma]

UndefManifold[M5]

7.7. Flat and Cartesian metrics

A flat metric is one without curvature. We define its Christoffel and epsilon objects, and then zero curvature tensors.

We can define a flat metric using the option FlatMetric:

DefMetric[-1,gflat[-a,-b],Cdflat,SymbolOfCovD->{"%","∂"},FlatMetric->True]

UndefMetric[gflat]

A delicate issue is that of a Cartesian metric. The concept of a flat metric is covariantly defined and therefore perfectly fits in the current structure of xTensor`. However we say that a metric in a given basis of vectors is Cartesian if all metric components in that basis are constant (not necessarily of unit modulus). For a flat metric coordinate systems of that kind always exist. We introduce the possibility of associating PD as the covariant derivative of a flat metric. That means that PD is the partial derivative of one of those coordinate systems with respect to that metric and hence PD is not general anymore. This is an ugly trick for a problem which is only properly solved in the twin package xCoba`.

We redefine the flat metric. In this case not even the Christoffel symbol is defined. By definition it is zero.

DefMetric[-1,gflat[-a,-b],PD,SymbolOfCovD->{",","∂"},FlatMetric->True]

PD[-a][gflat[-b,-c]]

UndefMetric[gflat]

7.8. Induced metrics

xTensor` has a simple way to define the induced metric obtained from projection of another metric along hypersurface-orthogonal vector fields, as it is typically done in the ADM 3+1 decomposition of the spacetime metric:

We first define the global ambient metric:

DefMetric[-1,metricg[-a,-b],CD,PrintAs->"g"];

Define the timelike normal to the slices. The use of MakeRule and AutomaticRules will be explained later.

DefTensor[n[a],M3]

AutomaticRules[n,MakeRule[{n[a]n[-a],-1}]]

AutomaticRules[n,MakeRule[{metricg[-a,-b]n[a]n[b],-1}]]

We choose the standard conventions in Numerical Relativity:

$ExtrinsicKSign=-1;
$AccelerationSign=-1;

Now we can define the induced metric associated to the pair {ambient metric, orthogonal vector field}. Note that even though there is another metric on the tangent bundle, an induced metric is not considered to be a frozen metric.

DefMetric[1,metrich[-a,-b],cd,SymbolOfCovD->{"|","D"},InducedFrom->{metricg,n},PrintAs->"h"]

For a hypersurface orthogonal vector we can express the acceleration vector in terms of a lapse function:

DefTensor[lapse[],M3,PrintAs->"α"]

LapseRule=Accelerationn[a_]->cd[a][lapse[]]/lapse[]

We define a generic vector field and a spatial vector field:

DefTensor[V[-a],M3]

DefTensor[W[-a],M3,OrthogonalTo->{n[a]},ProjectedWith->{metrich[a,-b]}]

The first operation we need is the decomposition of any expression along the normal vector. This is given by the command InducedDecomposition, which requires as second argument the pair {projector, normal}. The decomposition process does not change the input; it simply rewrites it in a different way:

Any expression can be decomposed in projected and orthogonal parts with respect to the vector field. Note the special output of the head Projectormetrich:

InducedDecomposition[W[-a],{metrich,n}]

InducedDecomposition[V[-a],{metrich,n}]

InducedDecomposition[7V[-a]V[-b]+W[-a]W[-b],{metrich,n}]

InducedDecomposition[RiemannCD[-a,-b,-c,-d],{metrich,n}]

ToCanonical[%]

Projection along the normal vector is automatic:

n[d]%//Expand

n[c]%//Expand

n[b]%%//Expand

The head Projectormetrich is the formal projector, but we now need a command which actually performs the projection. This is done with ProjectWith[metrich]. Do not confuse them. They have been constructed as a matched pair, and it is very useful to go from the former to the latter.

As long as we need to work only with the projected part of RicciCD, we can use this object:

Projectormetrich[RicciCD[-a,-b]]

However, if we need to perform the projection, we can do this:

%/.Projectormetrich->ProjectWith[metrich]

Finally, we can expand the projectors as

%//ProjectorToMetric

%//Expand

The projection process can be easily undone using InducedDecomposition:

InducedDecomposition[%,{metrich,n}]//ToCanonical

Projection operators.

The second issue is the manipulation of derivatives. The projected metric has its own Levi-Civita connection, which is a proper covariant derivative only when acting on objects on the projected manifold. The relation between the original and the new covariant derivatives is given by  a kind of "Christoffel" encoded in two objects: the extrinsic curvature ExtrinsicK tensor and the Acceleration vector:

This is the derivative of the normal (note the sign convention):

CD[-a][n[b]]

%//GradNormalToExtrinsicK

It is useful to remove derivatives in favour of these two tensors whenever possible. For example for "spatial" vectors:

n[-a]CD[-b][W[a]]

%//GradNormalToExtrinsicK//ContractMetric

The opposite command also exists:

ExtrinsicKmetrich[-a,-b]

%//ExtrinsicKToGradNormal//ContractMetric

Let us now project the derivative on a "spatial" vector:

CD[-c][W[-b]]

%//ProjectWith[metrich]

%//ProjectorToMetric

%//Expand

%//GradNormalToExtrinsicK

%//Simplification

A different way to do it is this:

cd[-c][W[-b]]

ProjectDerivative[%]

%/.Projectormetrich->ProjectWith[metrich]

Induced derivatives.

A very important detail is that of the Leibnitz rule for induced derivatives. It is not valid on expressions not on the projected vbundle!

The vector W lives on the projected vbundle and therefore the Leibnitz rule is valid:

cd[-a][W[b]W[c]]

If there are free non-spatial indices we send an error message, and the input expression is returned unevaluated:

cd[-a][V[b]V[c]]

With dummy indices we can expand both spatial indices and non-spatial indices:

cd[-a][W[a]W[b]]

cd[-a][V[b]V[-b]]//ToCanonical

In other cases, we prefer to show a warning and leave the expression untouched:

cd[-a][V[a]V[b]]

The most important application of this in GR is the ADM-type decomposition of a spacetime. In the rest of this section we reproduce some equations from the well-organized notes on ADM by M. Choptuik.

1) Foliations and normals. We use his same notations for the normal vector and the acceleration vector:

n[a]CD[-a][n[b]]//GradNormalToExtrinsicK//Expand

2) The projection tensor and the spatial metric. We do not follow the sign switch of York. The orthogonal projector is our pair Projectormetrich / ProjectWith[metrich] as explained above.

3) The spatial derivative operator and curvature tensor.

cd[-a][metrich[-b,-c]]

(cd[-a]@cd[-b][#]-cd[-b]@cd[-a][#])&@W[-c]

%//SortCovDs//ToCanonical

Riemanncd[-a,-b,-c,d]{n[a],n[b],n[c],n[-d]}

Riemanncd[-a,-c,-b,c]

The contraction of Ricci into RicciScalar does not work:

{Riccicd[-a,a],metrich[a,b]Riccicd[-a,-b],metricg[a,b]Riccicd[-a,-b]}

4) The extrinsic curvature tensor

CD[-a][n[-b]]//GradNormalToExtrinsicK//ContractMetric

Antisymmetrize[ProjectWith[metrich][%],{-a,-b}]//ContractMetric//ToCanonical

ProjectWith[metrich][-1/2LieD[n[c],CD][metricg[-a,-b]]//GradNormalToExtrinsicK]//ContractMetric//ToCanonical

-1/2LieD[n[c],CD][metrich[-a,-b]]//ProjectorToMetric//GradNormalToExtrinsicK//ContractMetric//ToCanonical

5) The Gauss-Codazzi equations. Single derivative:

UndefTensor[v]

DefTensor[v[a],M3,OrthogonalTo->{n[-a]},ProjectedWith->{metrich[-a,b]}]

CD[-c][v[-b]]

%//ProjectWith[metrich]

%//ProjectorToMetric

%//Expand

%//GradNormalToExtrinsicK

%//Simplification

We compare with expression (47) of Choptuik's notes:

-%+CD[-c][v[-b]]-n[-b]v[-f]CD[-c][n[f]]+n[-c]n[e]CD[-e][v[-b]]-n[-c]n[-b]v[-f]Accelerationn[f]

%//GradNormalToExtrinsicK

%//Simplification

Second derivatives:

ProjectWith[metrich][CD[-d][CD[-c][v[-b]]]]-cd[-d][cd[-c][v[-b]]]

%//ProjectDerivative//ProjectDerivative

%/.Projectormetrich->ProjectWith[metrich]

%//ProjectorToMetric//Expand

%//GradNormalToExtrinsicK//Expand

%//ContractMetric

%//Simplification

First Gauss-Codazzi equation:

expr1=Projectormetrich[CD[-d][CD[-c][v[-b]]]-CD[-c][CD[-d][v[-b]]]]

expr1A=expr1//SortCovDs

expr1/.Projectormetrich->ProjectWith[metrich]//ProjectorToMetric//Expand

expr1B=%//GradNormalToExtrinsicK//Simplification

expr2=cd[-d][cd[-c][v[-b]]]-cd[-c][cd[-d][v[-b]]]

expr2A=expr2//SortCovDs

expr2//ProjectDerivative//ProjectDerivative

%/.Projectormetrich->ProjectWith[metrich]//ProjectorToMetric//Expand

%//GradNormalToExtrinsicK//ContractMetric

expr2B=%//Simplification

(expr1A-expr1B)-(expr2A-expr2B)//Simplification

IndexSet[GCRiemann1[-a_,-b_,-c_,-d_],Simplification[(%/.v[_]->1//ScreenDollarIndices)]]

If we decompose the Riemann tensor along the vector n

inducedRiemann[-a_,-b_,-c_,-d_]=InducedDecomposition[RiemannCD[-a,-b,-c,-d],{metrich,n}]//ToCanonical

then we can express it as

inducedRiemann[-a_,-b_,-c_,-d_]=inducedRiemann[-a,-b,-c,-d]-GCRiemann1[-a,-b,-c,-d]//ToCanonical

Second Gauss-Codazzi equation:

expr1=Projectormetrich[CD[-d][CD[-c][n[-b]]]-CD[-c][CD[-d][n[-b]]]]

expr1A=expr1//SortCovDs

expr1/.Projectormetrich->ProjectWith[metrich]//ProjectorToMetric//Expand

expr1B=%//GradNormalToExtrinsicK//ContractMetric//GradNormalToExtrinsicK//Simplification

expr2A=cd[-c][ExtrinsicKmetrich[-b,-d]]-cd[-d][ExtrinsicKmetrich[-b,-c]]

expr2A//ProjectDerivative

%/.Projectormetrich->ProjectWith[metrich]//ProjectorToMetric//Expand

expr2B=%//GradNormalToExtrinsicK//ContractMetric//Simplification

IndexSet[GCRiemann2[-b_,-c_,-d_],(expr1A-expr1B)-(expr2A-expr2B)//Simplification]

We can further decompose the Riemann tensor now:

inducedRiemann[-a_,-b_,-c_,-d_]=inducedRiemann[-a,-b,-c,-d]-n[-d]GCRiemann2[-c,-a,-b]+n[-c]GCRiemann2[-d,-a,-b]-n[-b]GCRiemann2[-a,-c,-d]+n[-a]GCRiemann2[-b,-c,-d]//ToCanonical

"Third Gauss-Codazzi equation":

expr1=Projectormetrich[n[d](CD[-d][CD[-c][n[-b]]]-CD[-c][CD[-d][n[-b]]])]

expr1A=expr1//SortCovDs

expr1/.Projectormetrich->ProjectWith[metrich]//ProjectorToMetric//Expand

expr1B=%//GradNormalToExtrinsicK//ContractMetric//GradNormalToExtrinsicK//Simplification

expr2A=-cd[-c][Accelerationn[-b]]-Projectormetrich[n[d]CD[-d][ExtrinsicKmetrich[-b,-c]]]

expr2A//ProjectDerivative

%/.Projectormetrich->ProjectWith[metrich]//ProjectorToMetric//Expand

expr2B=%//GradNormalToExtrinsicK//ContractMetric//Simplification

Now we have (note that this pair of equations is equivalent to (105) of Choptuik):

IndexSet[GCRiemann3[-b_,-c_],(expr1A-expr1B)-(expr2A-expr2B)//Simplification]

ProjectWith[metrich][LieD[n[a],CD][ExtrinsicKmetrich[-b,-c]]-n[a]CD[-a][ExtrinsicKmetrich[-b,-c]]//GradNormalToExtrinsicK]//Simplification//ContractMetric

Finally we can combine the three Gauss-Codazzi equations into a single expression:

inducedRiemann[-a_,-b_,-c_,-d_]=inducedRiemann[-a,-b,-c,-d]+n[-b]n[-d]GCRiemann3[-a,-c]-n[-b]n[-c]GCRiemann3[-a,-d]-n[-a]n[-d]GCRiemann3[-b,-c]+n[-a]n[-c]GCRiemann3[-b,-d]//ToCanonical

GaussCodazziRiemannRule=IndexRule[RiemannCD[-a,-b,-c,-d],inducedRiemann[-a,-b,-c,-d]]

This expression is stored in xTensor` as

GaussCodazzi[RiemannCD[-a,-b,-c,-d],metrich]

%-Last[GaussCodazziRiemannRule]//ToCanonical

Most of the nontrivial computations in this area can now be easily performed using the Gauss-Codazzi formula. For example, let us compute the projections of the Ricci tensor and the Ricci scalar:

Decompose:

inducedRicci[-a_,-b_]=InducedDecomposition[RicciCD[-a,-b],{metrich,n}]//NoScalar

First formula:

zero=RicciCD[-a,-b]-inducedRiemann[-a,-c,-b,-d] metricg[c,d]//ContractMetric//NoScalar

zeroA=ProjectWith[metrich][%]//ContractMetric//ToCanonical

zeroB=Projectormetrich[RicciCD[-a,-b]]-ProjectWith[metrich][RicciCD[-a,-b]]

IndexSet[GCRicci1[-a_,-b_],zeroA+zeroB//ToCanonical]

inducedRicci[-a_,-b_]=inducedRicci[-a,-b]-GCRicci1[-a,-b]//ToCanonical

Second formula:

zero=(RicciCD[-a,-c]-inducedRiemann[-a,-b,-c,-d] metricg[b,d])n[c]//ContractMetric//NoScalar

zeroA=ProjectWith[metrich][%]//MetricToProjector//ToCanonical

zeroB=Projectormetrich[RicciCD[-a,-b]n[b]]-ProjectWith[metrich][RicciCD[-a,-b]n[b]]

IndexSet[GCRicci2[-a_],zeroA+zeroB//ToCanonical]

inducedRicci[-a_,-b_]=inducedRicci[-a,-b]+GCRicci2[-a]n[-b]+GCRicci2[-b]n[-a]//ToCanonical

Third formula:

zero=(RicciCD[-a,-c]-inducedRiemann[-a,-b,-c,-d] metricg[b,d])n[a]n[c]//Expand//ReplaceDummies

inducedRicci[-a_,-b_]=inducedRicci[-a,-b]-n[-a]n[-b]zero//Expand//SameDummies

Finally, the Ricci scalar:

inducedRiemann[-a,-b,-c,-d]metricg[b,d]metricg[a,c]//ContractMetric//NoScalar

%/.Projectormetrich->ProjectWith[metrich]//ContractMetric//ToCanonical

%//ProjectorToMetric//ContractMetric

%//GradNormalToExtrinsicK

With a lapse:

inducedRicciScalar[]=%/.LapseRule//Expand

UndefMetric[metrich]

UndefTensor/@{V,W};

Remove[zero,zeroA,zeroB,expr1,expr1A,expr1B,expr2,expr2A,expr2B,GCRiemann1,GCRiemann2,GCRiemann3,GCRicci1,GCRicci2,GaussCodazziRiemannRule,LapseRule,inducedRiemann,inducedRicci,inducedRicciScalar]

7.9. Products of epsilon tensors

In many applications we need to manipulate products of two epsilon tensors. They are automatically converted into products of metrics:

Products without contracted indices:

epsilonmetricg[a,-b,-c]epsilonmetricg[-d,-e,f]

Or with contracted indices:

epsilonmetricg[a,b,c]epsilonmetricg[-d,-e,-c]

epsilonmetricg[a,b,c]epsilonmetricg[-b,-e,-c]

epsilonmetricg[a,b,c]epsilonmetricg[-a,-b,-c]

7.10. Implosion and explosion (II)

Implosion of a tensor is only possible directly when its indices have the correct characters (those specified at definition time). In all other cases we need to pull out a number of metric factors first, to correct them.

The tensor v was defined with an upper index. Hence here we get an additional derivative of the metric, which is also imploded:

Cd[-a][v[-b]]

Implode[%]

InputForm[%]

Explosion now does not return to the original expression:

Explode[%%]

but this expression is equivalent to the original one:

ContractMetric[%,OverDerivatives->True]

7.11. Variational derivatives

There is not yet a concept of integration in xTensor`. Instead of working with a variational derivative of a functional, we shall assume that such a functional is the integral of some integrand (or "Lagrangian") and that the integral is computed with respect to a volume form which vanishes under the action of some derivative (this is the derivative which will be "integrated by parts").

For instance, let us compute the RicciScalar field of our metric:

rs=RicciScalarCD[]//RiemannToChristoffel//NoScalar//ChristoffelToGradMetric//Expand

Working intensively with partial derivatives, it is interesting to set

SetOptions[ToCanonical,UseMetricOnVBundle->None]

Now we perform a direct variational derivative of the Einstein-Hilbert Lagrangian with respect to the inverse metric. It is a long computation:

AbsoluteTiming[Length[result=ToCanonical@VarD[metricg[a,b],PD][Sqrt[-Detmetricg[]]rs]]]

The expected result is the Einstein tensor:

EinsteinCD[-a,-b]//EinsteinToRicci//RiemannToChristoffel//NoScalar//ChristoffelToGradMetric//ToCanonical

result-Sqrt[-Detmetricg[]]%//ToCanonical

Let us now define the electromagnetic fields:

DefTensor[MaxwellA[-a],M3,PrintAs->"A"]

DefTensor[MaxwellF[-a,-b],M3,Antisymmetric[{1,2}],PrintAs->"F"]

MaxwellF[a_,b_]:=PD[a][MaxwellA[b]]-PD[b][MaxwellA[a]]

This is the electromagnetic Lagrangian:

metricg[a,b]metricg[c,d]MaxwellF[-a,-c]MaxwellF[-b,-d]/4

and these are the Maxwell equations on a curved background:

result=VarD[MaxwellA[a],PD][Sqrt[-Detmetricg[]]%]

metricg[b,c]CD[-b][MaxwellF[-a,-c]]//ChangeCovD//ChristoffelToGradMetric//ToCanonical

result-Sqrt[-Detmetricg[]]%//ToCanonical

UndefTensor/@{MaxwellA,MaxwellF};

Restore standard options:

SetOptions[ToCanonical,UseMetricOnVBundle->All]

7.12. Clean up

Again, we remove all objects in our session. This must be done in the correct order, to respect the relations among the symbols. Otherwise, xTensor will refuse to remove the corresponding symbols, to avoid breaking internal relations.

Start by undefining the metric fields:

UndefMetric[metricg]

Undefine tensors with imploded variants before removing the derivative:

UndefTensor/@{r,v,T};

Undefine the derivatives:

UndefCovD[ICD]

UndefCovD[Cd]

Undefine the rest of the tensors:

UndefTensor/@$Tensors;

Undefine inner vector bundles:

UndefVBundle[InnerC]

Finally undefine the manifolds:

UndefManifold[S2]

UndefManifold[M3]

There are still symbols in the Global` , not under control of xTensor. We can remove them all now:

?Global`*

Remove["Global`*"]

Our session is clean again:

?Global`*

7. Mappings and tensors along mappings. V1.1.0 NEW

8.1. Diffeomorphic mapping

Since version 1.1.0 of xTensor we can define mappings between manifolds and use them to push/pull tensor fields from one vbundle to another. As in the rest of xAct, objects must be considered only locally. All manifolds in xTensor are real, and so are the mappings between any two of them.

Definition of mappings.

Define two manifolds of the same dimension, with their respective tangent bundles:

DefManifold[M,4,{a,b,c,d}]

DefManifold[P,4,{A,B,C,D}]

Define a mapping from the manifold M to the manifold P. By default mappings between manifolds of the same dimension are assumed to be diffeomorphisms (we will see later how to change this).

DefMapping[ϕ,M->P]

Together with the mapping ϕ we have defined its inverse mapping, and also their differentials, i.e., the linear operators relating the tangent bundles. Let us first have a look at the mappings themselves, and their properties.

Properties of mappings.

A mapping has a domain manifold and an image (or codomain) manifold.

MappingDomain[ϕ]->MappingImage[ϕ]

This is the inverse mapping of ϕ. By default it is named by prepending i:

InverseMapping[ϕ]

The inverse mapping is defined from the image to the domain of the original mapping:

MappingDomain[iϕ]->MappingImage[iϕ]

Both mappings have all good properties in the (local) diffeomorphic case:

{ImmersionQ[ϕ],SubmersionQ[ϕ],DiffeomorphismQ[ϕ]}

{ImmersionQ[iϕ],SubmersionQ[iϕ],DiffeomorphismQ[iϕ]}

Because ϕ is a submersion, the image manifold coincides with the range:

ϕ[M]

The same can be said of its inverse:

iϕ[P]

8.2. Non-diffeomorphic mapping

Mappings do not have to relate manifolds of the same dimension. Or even if they do, they might fail to be immersions or submersions.

Define a manifold of higher dimension:

DefManifold[L,5,{α,β,γ,δ}]

Now we can define another mapping, this time being only an immersion:

DefMapping[ψ,P->L,SubmersionQ->False]

These are the properties of the mapping ψ:

{ImmersionQ[ψ],SubmersionQ[ψ],DiffeomorphismQ[ψ]}

MappingDomain[ψ]->MappingImage[ψ]

The range does not coincide with all of the image manifold. It is just a submanifold of the image:

ψ[P]

SubmanifoldQ[L,ψ[P]]

The inverse mapping is defined on all of the image of ψ, because we assume it has been extended in some smooth arbitrary way. This is just a simplification to avoid working with subvbundles of TangentL, and might be improved in some future version of xTensor.

MappingDomain[iψ]->MappingImage[iψ]

Properties of the inverse mapping iψ:

{ImmersionQ[iψ],SubmersionQ[iψ],DiffeomorphismQ[iψ]}

It is only a left-inverse mapping, but not a right-inverse mapping of ψ:

iψ[ψ[P]]

ψ[iψ[L]]

8.3. Mapping composition

Mappings can be composed, resulting in other mappings. To avoid interferences with Composition, and to have nicer typesetting and an alias, we use SmallCircle to represent mapping composition.

Composition of mappings.

Compose the mappings ϕ and ψ, following the standard right-to-left convention. Note the alias esc-sc-esc for SmallCircle:

map=ψ∘ϕ

MappingDomain[map]->MappingImage[map]

Inverse of the composition. This requires all individual mappings to be invertible (at least left- or right- invertible). Note the reversal of ordering.

imap=InverseMapping[map]

Properties of the composition:

{ImmersionQ[map],SubmersionQ[map],DiffeomorphismQ[map]}

Because the composed map is an immersion but not a submersion, its inverse is a left-inverse but not a right-inverse:

imap∘map

%//InputForm

map∘imap

More explicitly:

imap[map[M]]

map[imap[L]]

The identity mapping is denoted by IdentityMapping[manifold], explicitly indicating the manifold where it acts. It has all the expected properties.

The identity mapping on the manifold M:

idM=IdentityMapping[M]

MappingDomain[idM]->MappingImage[idM]

It has all good properties, and it is its own inverse:

{ImmersionQ[idM],SubmersionQ[idM],DiffeomorphismQ[idM]}

InverseMapping[idM]

idM[M]

It is the identity of composition:

ϕ∘idM

IdentityMapping[P]∘ϕ

8.4. The pullback bundle and its tensor fields

Given a vector bundle VB with base manifold P and a mapping ϕ: M->P we can construct another bundle with base M and fibers those of VB. It is called the pullback bundle of VB by ϕ. This is effectively an abstraction of the concept of reparametrization.

Pull back vector bundle.

This is the tangent bundle of the manifold P, with 4-dimensional fibers:

Tangent[P]

{VBundleQ[%],BaseOfVBundle[%],DimOfVBundle[%]}

This is its pullback bundle by the mapping ϕ, now with base M and the same 4-dimensional fibers:

PullBackVBundle[TangentP,ϕ]

{VBundleQ[%],BaseOfVBundle[%],DimOfVBundle[%]}

Note its internal non-symbolic representation

%%//InputForm

Even though the fibers are isomorphic, the associated indices must be different because indices must identify vbundles uniquely in xAct. To avoid requesting new indices from the user, we prepend the name of the mapping to the indices:

IndicesOfVBundle[TangentP]

IndicesOfVBundle[PullBackVBundle[TangentP,ϕ]]

Tensor fields defined on the pullback bundles are called “tensors along a map” or “pullback sections”. They will carry distinctive “starred indices”, as a reminder that dependencies have been changed (this is a bit like the DeWitt index notation). In this way, for example, we will not be able to contract tensors belonging to different vbundles. The most important example is the differential of the mapping, described in the next section. We will then see two other examples: precomposition and linear push of tensor fields.

Since version 1.1.0 of xTensor, tensors can be natively tensors along a map. This is a section of the pullback bundle ϕ*𝕋P, because we use an index of that vbundle:

DefTensor[W[ϕA],M]

SlotsOfTensor[W]

DependenciesOfTensor[W]

Indices of a pullback bundle are represented with a star. We need to look at their input form to see the actual mapping.

W[ϕA]

Algebraic operations proceed normally:

W[ϕA]delta[-ϕA,ϕB]-7W[ϕB]

Undefine the field:

UndefTensor[W]

Note: in our scheme of things, “pull back vbundle” is not the best name. It should have been “precomposed vbundle”, but the former is standard, so that is what we use.

8.5. Differential of a mapping

Associated to each smooth mapping we have its differential, also known as “tangent map”. This is a linear operator mapping vectors on a vbundle to vectors on a pullback bundle of the same base. It will hence be represented in xTensor as a 2-tensor with a normal index and a starred index.

Mapping differentials.

This is the tangent map of ϕ, by default constructed by prepending d to the name of the mapping:

TangentTensor[ϕ]

It is a tensor field on the manifold M:

xTensorQ[dϕ]

DependenciesOfTensor[dϕ]

Its two slots belong to different vbundles, the second one being a pullback bundle:

SlotsOfTensor[dϕ]

Both bundles have the same base, which coincides with the dependency of the tangent map (that’s why it would be incorrect to have 𝕋P as a slot of dϕ: its dependency would be P):

BaseOfVBundle/@%

With indices we have the following expression, where we need a starred index.

dϕ[-a,ϕB]

The inverse mapping iϕ has its own tangent map:

TangentTensor[iϕ]

This is now a field on P, and has slots in -𝕋P (this is xAct notation for the dual (SuperStar[𝕋])P of 𝕋P) and the pullback by iϕ of 𝕋M:

(Jose: The -1 of the inverse is inconvenient. Use iϕ instead)

DependenciesOfTensor[diϕ]

SlotsOfTensor[diϕ]

BaseOfVBundle/@%

Using the corresponding indices we have:

diϕ[-A,iϕb]

dϕ and diϕ are NOT mutual inverses. That’s clear in two facts: 1) they are fields on different manifolds, and 2) they have indices on different vbundles. However there is a close relation. We use Inv to compute the inverse of a 2-tensor:

Inv[dϕ][-ϕA,b]

%//InputForm

If the mapping is not a diffeomorphism then its differential may be only left-invertible or right-invertible or none.

(Jose: The second argument of InvertibleQ is ugly. Use separate functions LeftInvertibleQ and RightInvertibleQ)

Because ϕ is a diffeomorphism its tangent map is invertible from both sides:

{InvertibleQ[dϕ,Left],InvertibleQ[dϕ,Right],InvertibleQ[dϕ]}

And so is the differential of its inverse:

{InvertibleQ[diϕ,Left],InvertibleQ[diϕ,Right],InvertibleQ[diϕ]}

However, the mapping ψ is only an immersion (note the unfortunate fact that ψ is left-invertible but its differential is right invertible, due to our choice down-up for the delta, instead of the standard up-down):

{InvertibleQ[dψ,Left],InvertibleQ[dψ,Right],InvertibleQ[dψ]}

dψ[-A,ψα]Inv[dψ][-ψα,B]

Inv[dψ][-ψα,A]dψ[-A,ψβ]

That’s a rather complicated expression, but it is what it is.

The properties of the differential of iψ are reversed with respect to those of ψ:

{InvertibleQ[diψ,Left],InvertibleQ[diψ,Right],InvertibleQ[diψ]}

diψ[-ψα,A]Inv[diψ][-A,ψβ]

Inv[diψ][-A,ψα]diψ[-ψα,B]

8.6. Precomposition

Using precomposition we modify the dependencies of a tensor, but do not change essentially its slot structure (we only “add stars”). In particular, this operation is still nontrivial on scalar fields.
Precomposition is very similar to composition (SmallCircle). We need to distinguish them because mappings and tensor fields are two different types of objects in xTensor. That is, a mapping is composed with another mapping, and a tensor field is precomposed with a mapping. This is a fundamental design decision.

Precomposition.

By precomposition we change the properties of a tensor. For example:

DefTensor[T[A,-B],P]

DependenciesOfTensor[T]

SlotsOfTensor[T]

If we precompose the tensor we get:

DependenciesOfTensor[Precompose[T,ϕ]]

SlotsOfTensor[Precompose[T,ϕ]]

That explains why the properties of the precomposition of diϕ are compatible with those of dϕ:

DependenciesOfTensor[Precompose[diϕ,ϕ]]

SlotsOfTensor[Precompose[diϕ,ϕ]]

DependenciesOfTensor[dϕ]

SlotsOfTensor[dϕ]

And hence now we can contract in four different ways:

dϕ[-a,ϕB]Inv[dϕ][-ϕB,c]

Inv[dϕ][-ϕA,b]dϕ[-b,ϕC]

diϕ[-A,iϕb]Inv[diϕ][-iϕb,C]

Inv[diϕ][-iϕa,B]diϕ[-B,iϕc]

Another possible example of the need to be careful with indices is in the computation of the differential of a composition. The differential of ψ∘ϕ is the product of their differentials, but these are fields in different vbundles, so there must be precomposition.

We will work with a double pullback vbundle,

PullBackVBundle[TangentL,ψ]

PullBackVBundle[%,ϕ]

%//InputForm

whose indices are “double mapped”:

IndicesOfVBundle[%%]

The order of mappings is reversed in the index names (as a result of pull back). It would have been more natural to write the mappings on the right or the indices, to keep the order, but we already have the dagger, dollar and number postfixes for indices, and I prefer not to mess up with those other more important marks.

$PrePrint=ScreenDollarIndices;

This is the differential of the composition ψ∘ϕ (as a contraction of linear operators, the order is reversed unfortunately, as we mentioned before):

TangentTensor[ψ∘ϕ][-a,ϕψα]

As we had anticipated, we need precomposition of one of the differentials in order to have consistent dependencies:

DependenciesOf[%]

8.7. Pullback

Here we find the first fundamental operation of working with tensors and mappings: the pullback of a tensor field on the image manifold to a tensor field on the domain manifold. This operation is always defined for covariant tensor fields, and can be extended to contravariant tensor fields if the mapping is an immersion (QUESTION: is the latter the right statement?).

TODO:
	- Document PullBackTensor.
	- Document PullBackCovD and complete its encoding.
	- Encode how standard operations (like Lie derivatives and brackets) behave under pull back and precomposition.

Pull back.

Declare a covector field on -𝕋P:

DefTensor[Ω[-A],P]

DependenciesOfTensor[Ω]

SlotsOfTensor[Ω]

We can pull it back to -𝕋M. We need to specify how the free indices must be renamed:

pbexpr=PullBack[Ω[-A],ϕ,{-A->-a}]

%//InputForm

DependenciesOfTensor[PullBackTensor[Ω,ϕ]]

SlotsOfTensor[PullBackTensor[Ω,ϕ]]

We can expand the pulled-back expression in terms of precomposition and contraction with the tangent maps:

pbexpr//ExpandPullBack

%//InputForm

A higher-rank covariant field:

DefTensor[Υ[-A,-B,-C],P]

PullBack[Υ[-A,-B,-C],ϕ,{-A->-a,-B->-b,-C->-c}]

%//ExpandPullBack

DependenciesOf[%]

Now with contravariant slots and a contraction. We only need to specify the change of free indices. Note the need of the differential inverse, but there is no need to use the inverse mapping:

DefTensor[Λ[A,B,-C,-D],P]

PullBack[Λ[A,B,-B,-C],ϕ,{A->a,-C->-c}]

%//ExpandPullBack

A further contraction gives a precomposed scalar:

% delta[-a,c]

8.8. Linear push with the mapping differential

This operation is dual with respect to precomposition. Here we do not change the dependencies of a tensor field, but we do act linearly on its indices. Effectively, this is nothing but a condensed notation to denote contractions with the mapping differentials. This operation is trivial on scalar fields.
The name is far from ideal, but I don’t see what else to choose.

TODO:
	- Analyze what happens to the covariant derivatives of linear push fields.
	- Encode how the standard operations (like Lie derivatives and brackets) behave under linear push and pushforward.

Linear push.

Take a tensor field in 𝕋M:

DefTensor[v[a],M]

Push it to the pullback bundle (SuperStar[ϕ])𝕋P:

LinearPush[v,dϕ][ϕA]

It is a vector field on the pullback bundle:

SlotsOfTensor[Head[%]]

It is still a field on M:

DependenciesOf[%%]

Expand a LinearPush expression:

LinearPush[v,dϕ][ϕA]//ExpandPushForward

%//InputForm

Now take a covector field in -𝕋M:

DefTensor[ω[-a],M]

It is possible to push covector fields in this case because the mapping differential is left-invertible:

LinearPush[ω,dϕ][-ϕA]

SlotsOfTensor[Head[%]]

That is clearer if we expand the expression into precomposition and contraction:

LinearPush[ω,dϕ][-ϕA]//ExpandPushForward

%//InputForm

8.9. PushForward. MINOR BUG

And this is our second fundamental operation of working with tensors and mappings: the pushforward of a tensor field on the domain manifold to a tensor field on the image manifold. This operation always requires an invertible mapping, because we need to use precomposition along the inverse mapping to change the dependencies. QUESTION: Fully invertible or is it enough with being a submersion?

TODO:
	- Document PushForwardTensor.

Push forward.

Push forward the vector field v on 𝕋M to a vector field on 𝕋P. Again, we need to specify the renaming of the free indices:

PushForward[3v[a],ϕ,{a->A}]

%//InputForm

It is now a vector field on P: (TODO: THIS IS WRONG!!!)

DependenciesOf[%]

We see that the actual operation requires the use of the inverse mapping:

%%//ExpandPushForward

In fact, push forward is equivalent to a pull back along the inverse mapping:

PullBack[3v[a],iϕ,{a->A}]

%//ExpandPullBack

Push forward a covector field. Now we need both the inverse mapping and its differential:

PushForwardTensor[ω,ϕ][-A]

%//ExpandPushForward

Push forward a mixed tensor, and with contractions:

DefTensor[S[a,b,-c],M]

PushForward[S[a,b,-b],ϕ,{a->A}]

SlotsOfTensor[Head[%]]

8.10. Clean up. TRIVIAL BUG

Remove the various objects we have declared in this section:

$Mappings

$Tensors

$Manifolds

It must be done in this order: first the mappings (to remove their differentials too), then the tensors and finally the manifolds:

UndefMapping/@{ϕ,ψ};

UndefTensor/@$Tensors;

UndefManifold/@$Manifolds;

Remove other symbols used:

Remove[idM,map,imap,pbexpr]

TODO: WE DID A TERRIBLE JOB AT REMOVING INDICES... WHAT ARE THE mapTangentM symbols?

?Global`*

8.11. Graphical summary

In the following diagram we have the connections between the manifolds M and P and their respective tangent bundles 𝕋M and 𝕋P. Fields in this diagram are lines connecting a manifold to a bundle. Vertical lines have “normal” indices, while diagonal lines have “starred” indices. We see how starting with v we can move forward its (upper) head using contraction dϕ, while starting with W we can move backward its tail using precomposition with ϕ. Adding the lines for iϕ is trivial. Moving the vertical line of v to W would be a full “push forward” and involves both moving the head and the tail. Moving W to v would be a full “pull back”, also involving both operations.

LayeredGraphPlot[{
{Style["M",Large]->Style["P",Large],Style["ϕ",Large]},
{Style["M",Large]->Style["𝕋M",Large],Style["v",Large]},
{Style["𝕋M",Large]->Style["𝕋P",Large],Style["dϕ",Large]},
{Style["P",Large]->Style["𝕋P",Large],Style["W",Large]},
{Style["M",Large]->Style["𝕋P",Large],Style["dϕ·v or W∘ϕ",Large]}
},
VertexLabeling->True,FormatType->StandardForm,
VertexCoordinateRules->{{0,0},{1,0},{0,1},{1,1}}
]/.Network`GraphPlot`wrapper->Identity

9. More on rules

9.1. MakeRule

Defining rules is one of the most important components of a tensor package. The function IndexRule takes care of dummy generation, but there are other points to consider when defining a rule:
	1. Checking consistency between left and right hand sides.
	2. Defining rules that work only for a given vbundle.
	3. Extending the rule to pattern expressions with equivalent index-structure under symmetries.
	4. Extendig the rule to patterns with a different up/down character if metrics are present.
We define the function MakeRule, which takes care of all these things. Note that there is not a parallel function for Set, but below we present a way to automatize a collection of rules.

Construction of rules.

Define again the environment for examples. This time we quite the info messages:

$DefInfoQ=False;

DefManifold[M3,3,{a,b,c,d,e,f,g,h}]
DefManifold[S2,2,{A,B,C,D}]//Quiet
DefMetric[-1,metricg[-a,-b],CD,PrintAs->"g"]
DefTensor[v[a],M3]
DefTensor[T[a,b],M3]
DefTensor[U[-a,-b,-c],M3,Antisymmetric[{1,2,3}]]

$PrePrint=ScreenDollarIndices;

The simplest way to define a rule would easily give corrupted answers:

rule=v[a_]->T[a,-b]v[b]

v[b]/.rule

We have seen that changing to IndexRule we can solve those problems (underlined indices denote patterns):

rule=v[a_]↦T[a,b]v[-b]

v[b]/.rule

%/.rule

We can also define the same rule using MakeRule. Note that there is no pattern on the lhs; by default all indices on the LHS are converted into patterns:

rule=MakeRule[{v[a],T[a,-b]v[b]},MetricOn->None]

T[-a,-b]v[a]v[b]/.rule

InputForm[rule]

Due to the option TestIndices->True, membership of the indices of the LHS is checked. The function M3`Q checks that an up-index belongs to M3, but it does not accept down-indices (this is due to the option MetricOn->None):

v[a]/.rule

v[-a]/.rule

v[A]/.rule

In order to allow both up/down-indices of a given vbundle we use MetricOn. We need to have a metric first. (If the vbundle of the index did not have a metric MakeRule would refuse to change the up/down character):

rule=MakeRule[{v[a],T[a,-b]v[b]},MetricOn->{a},ContractMetrics->True]

InputForm[rule]

v[a]/.rule

v[-a]/.rule

If we have symmetries the situation gets even more complicated:

(rule=MakeRule[{U[a,b,c]v[-c],v[a]T[b]-v[b]T[a]},MetricOn->All,ContractMetrics->True])

U[-a,b,c]v[-b]/.rule

U[-a,b,c]v[a]/.rule

U[-a,c,-c]v[a]/.rule//ToCanonical

It is possible to automate rules, declaring them as definitions for a symbol, in this case the tensor U:

AutomaticRules[U,rule]

U[-a,b,c]v[-b]

U[-a,b,c]v[a]

?U

The simplest way to remove those definitions for U is removing the tensor and defining it again:

UndefTensor[U]

DefTensor[U[-a,-b,-c],M3,Antisymmetric[{1,2,3}]]

9.2. Riemann expressions

As we said, there is no algorithm at present in xTensor` to canonicalize expressions with multiterm symmetries. That is an obstacle in GR when one wants to work with the Riemann tensor, which has a cyclic symmetry. A simple, brute-force, solution is to prepare a number of rules which are only valid for the Riemann tensor. We follow MathTensor's discussion of this point (see MathTensor book).
There are 48 rules in total, numbered from 1 to 40, some of them having subcases a, b, c, ...
14 of those rules are simple consequences of the monoterm symmetries of the Riemann tensor, and hence can be considered already implemented in ToCanonical:

RiemannRule1:

RiemannCD[-a,a,-c,-d]//ToCanonical

RiemannCD[-a,-b,-c,c]//ToCanonical

RiemannRule2:

RiemannCD[-a,-b,a,-d]//ToCanonical

RiemannRule3:

RiemannCD[-a,-b,-c,a]//ToCanonical

RiemannRule4:

RiemannCD[-a,-b,b,-d]//ToCanonical

RiemannRule5:

RiemannCD[-a,-b,-c,b]//ToCanonical

RiemannRule6:

RiemannCD[-a,-b,-c,-d]RicciCD[a,b]//ToCanonical

RiemannCD[-a,-b,-c,-d]RicciCD[c,d]//ToCanonical

RiemannRule31:

CD[-e]@RiemannCD[-a,-b,-c,-d]RicciCD[a,b]//ToCanonical

CD[-e]@RiemannCD[-a,-b,-c,-d]RicciCD[c,d]//ToCanonical

RiemannRule32:

CD[-e]@RiemannCD[-a,-b,-c,-d]CD[-f]@RicciCD[a,b]//ToCanonical

CD[-e]@RiemannCD[-a,-b,-c,-d]CD[-f]@RicciCD[c,d]//ToCanonical

RiemannRule33:

CD[-c]@RicciCD[-a,-b]RiemannCD[a,b,d,e]//ToCanonical

CD[-c]@RicciCD[-a,-b]RiemannCD[d,e,a,b]//ToCanonical

There are 2 more rules which can be obtained from the fact that the Einstein tensor is divergence free:

RiemannRule16:

CD[b][RicciCD[-a,-b]]//RicciToEinstein//ContractMetric

RiemannRule17:

CD[a][RicciCD[-a,-b]]//RicciToEinstein//ContractMetric

The other 32 rules are actual consequences of the cyclic symmetry of the Riemann tensor, and therefore must be programmed independently. Using MakeRule we can avoid typing equivalent rules (up to permutations of indices), reducing the number to 17 independent rules.
We first study 13 rules having derivatives of a single curvature tensor on the LHS. They are in fact only 6 independent rules:
	RuleDivRiemann
	RuleDivGradRicci
	RuleBoxRiemann
	RuleDivGradRiemann
	RuleD4RicciScalar
	RuleD4Ricci

The four rules 20--23 are all the same:

RuleDivRiemann=MakeRule[{CD[a]@RiemannCD[-a,-b,-c,-d],-CD[-d]@RicciCD[-b,-c]+CD[-c]@RicciCD[-b,-d]},UseSymmetries->True,MetricOn->All];

RiemannRule20:

CD[a]@RiemannCD[-a,-b,-c,-d]/.RuleDivRiemann

RiemannRule21:

CD[b]@RiemannCD[-a,-b,-c,-d]/.RuleDivRiemann

RiemannRule22:

CD[c]@RiemannCD[-a,-b,-c,-d]/.RuleDivRiemann

RiemannRule23:

CD[d]@RiemannCD[-a,-b,-c,-d]/.RuleDivRiemann

Create more indices:

GetIndicesOfVBundle[TangentM3,14]

Two rules with second derivatives of Ricci can be combined in a single rule:

RuleDivGradRicci=MakeRule[{CD[b]@CD[-c]@RicciCD[-a,-b],1/2CD[-c]@CD[-a]@RicciScalarCD[]+RicciCD[-h1,-a]RicciCD[-c,h1]-RicciCD[-h1,-h2]RiemannCD[-a,h2,-c,h1]},UseSymmetries->True,MetricOn->All];

RiemannRule18:

CD[b]@CD[-c]@RicciCD[-a,-b]/.RuleDivGradRicci

RiemannRule19:

CD[a]@CD[-c]@RicciCD[-a,-b]/.RuleDivGradRicci

Rules 26 to 30 are now two rules in xTensor`:

RuleBoxRiemann=MakeRule[{CD[e]@CD[-e]@RiemannCD[-a,-b,-c,-d],
CD[-b]@CD[-d]@RicciCD[-a,-c]-CD[-b]@CD[-c]@RicciCD[-a,-d]-
CD[-a]@CD[-d]@RicciCD[-b,-c]+CD[-a]@CD[-c]@RicciCD[-b,-d]+
RicciCD[-h1,-b]RiemannCD[-a,h1,-c,-d]-
RicciCD[-h1,-a]RiemannCD[-b,h1,-c,-d]-
RiemannCD[-h1,-h2,-c,-d]RiemannCD[-a,h1,-b,h2]-
RiemannCD[-h1,-b,-h2,-d]RiemannCD[-a,h1,-c,h2]+
RiemannCD[-h1,-b,-h2,-c]RiemannCD[-a,h1,-d,h2]+RiemannCD[-h1,-a,-h2,-d]RiemannCD[-b,h1,-c,h2]-RiemannCD[-h1,-a,-h2,-c]RiemannCD[-b,h1,-d,h2]-
RiemannCD[-h1,-a,-h2,-b]RiemannCD[-c,-d,h1,h2]},UseSymmetries->True,MetricOn->All];

RuleDivGradRiemann=MakeRule[{CD[a]@CD[-e]@RiemannCD[-a,-b,-c,-d],-CD[-e]@CD[-d]@RicciCD[-b,-c]+CD[-e]@CD[-c]@RicciCD[-b,-d]-RicciCD[-h1,-e]RiemannCD[-b,h1,-c,-d]+RiemannCD[-h1,-h2,-c,-d]RiemannCD[-b,h1,-e,h2]-RiemannCD[-h1,-b,-h2,-d]RiemannCD[-c,h2,-e,h1]+RiemannCD[-h1,-b,-h2,-c]RiemannCD[-d,h2,-e,h1]},UseSymmetries->True,MetricOn->All];

RiemannRule26:

CD[e]@CD[-e]@RiemannCD[-a,-b,-c,-d]/.RuleBoxRiemann

RiemannRule27:

CD[a]@CD[-e]@RiemannCD[-a,-b,-c,-d]/.RuleDivGradRiemann

RiemannRule28:

CD[b]@CD[-e]@RiemannCD[-a,-b,-c,-d]/.RuleDivGradRiemann

RiemannRule29:

CD[c]@CD[-e]@RiemannCD[-a,-b,-c,-d]/.RuleDivGradRiemann

RiemannRule30:

CD[d]@CD[-e]@RiemannCD[-a,-b,-c,-d]/.RuleDivGradRiemann

Fourth derivative of the Ricci scalar: conversion to Box^2:

RuleD4RicciScalar=MakeRule[{CD[b]@CD[a]@CD[-b]@CD[-a]@RicciScalarCD[],1/2CD[-h1]@RicciScalarCD[]CD[h1]@RicciScalarCD[]+CD[h2]@CD[-h2]@CD[h1]@CD[-h1]@RicciScalarCD[]+CD[-h2]@CD[-h1]@RicciScalarCD[]RicciCD[h1,h2]},UseSymmetries->True,MetricOn->All];

RiemannRule24:

CD[b]@CD[a]@CD[-b]@CD[-a]@RicciScalarCD[]/.RuleD4RicciScalar

Fourth derivative of the Ricci tensor:

RuleD4Ricci=MakeRule[{CD[b]@CD[a]@CD[c]@CD[-c]@RicciCD[-a,-b],1/2CD[-h1]@RicciScalarCD[]CD[h1]@RicciScalarCD[]-3CD[-h3]@RicciCD[-h1,-h2]CD[h3]@RicciCD[h1,h2]+4CD[-h3]@RicciCD[-h1,-h2]CD[h2]@RicciCD[h1,h3]+1/2CD[h2]@CD[-h2]@CD[h1]@CD[-h1]@RicciScalarCD[]+2CD[-h2]@CD[-h1]@RicciScalarCD[]RicciCD[h1,h2]-CD[h3]@CD[-h3]@RicciCD[-h1,-h2]RicciCD[h1,h2]+2RicciCD[-h1,-h2]RicciCD[-h3,h1]RicciCD[h2,h3]-2CD[-h4]@CD[-h3]@RicciCD[-h1,-h2]RiemannCD[h1,h3,h2,h4]-2RicciCD[-h1,-h2]RicciCD[-h3,-h4]RiemannCD[h1,h3,h2,h4]},UseSymmetries->True,MetricOn->All];

RiemannRule25:

CD[b]@CD[a]@CD[c]@CD[-c]@RicciCD[-a,-b]/.RuleD4Ricci

Now we study six rules with products of two Riemann tensors on the LHS. There are only two different rules, depending on whether there are two or three contracted indices:
	RuleTwoRiemann3
	RuleTwoRiemann2

These are the independent rules:

RuleTwoRiemann3=MakeRule[{RiemannCD[-a,-b,-c,-d]RiemannCD[-e,a,b,c],-1/2RiemannCD[-f,-d,-g,-h]RiemannCD[-e,f,g,h]},UseSymmetries->True,MetricOn->All];

RuleTwoRiemann2=MakeRule[{RiemannCD[-a,-b,-c,-d]RiemannCD[-e,a,-g,b],1/2RiemannCD[-f,-h,-c,-d]RiemannCD[-e,-g,f,h]},UseSymmetries->True,MetricOn->All];

RiemannRule7:

RiemannCD[-a,-b,-c,-d]RiemannCD[-e,a,b,c]/.RuleTwoRiemann3

RiemannRule8:

RiemannCD[-a,-b,-c,-d]RiemannCD[a,c,b,h]/.RuleTwoRiemann3

RiemannRule9a:

RiemannCD[-a,-b,-c,-d]RiemannCD[-e,a,-g,b]/.RuleTwoRiemann2

RiemannRule9b:

RiemannCD[-a,-b,-c,-d]RiemannCD[e,a,b,h]/.RuleTwoRiemann2

RiemannRule9c:

RiemannCD[-a,-b,-c,-d]RiemannCD[a,f,b,h]/.RuleTwoRiemann2

RiemannRule9d:

RiemannCD[-a,-b,-c,-d]RiemannCD[a,f,g,b]/.RuleTwoRiemann2

There are seven more rules involving two curvature tensors and derivatives, corresponding to five independent rules:
	RuleD01TwoRiemann
	RuleD11TwoRiemann
	RuleD2RicciRiemann
	RuleD2RicciDRiemann
	RuleD11bisTwoRiemann

RuleD01TwoRiemann=MakeRule[{CD[-e]@RiemannCD[a,f,b,c]RiemannCD[-a,-b,-c,-d],-1/2CD[-e]@RiemannCD[h1,f,h2,h3]RiemannCD[-h1,-d,-h2,-h3]},UseSymmetries->True,MetricOn->All];

RuleD11TwoRiemann=MakeRule[{CD[-e]@RiemannCD[-a,-b,-c,-d]CD[-f]@RiemannCD[a,g,b,c],-1/2CD[-e]@RiemannCD[-h1,-d,-h2,-h3]CD[-f]@RiemannCD[h1,g,h2,h3]},UseSymmetries->True,MetricOn->All];

RuleD2RicciRiemann=MakeRule[{CD[-d]@CD[-c]@RicciCD[-a,-b]RiemannCD[a,d,b,c],CD[-h4]@CD[-h3]@RicciCD[-h1,-h2]RiemannCD[h1,h3,h2,h4]},UseSymmetries->True,MetricOn->All];

RuleD2RicciDRiemann=MakeRule[{CD[-e]@RiemannCD[a,d,b,c]CD[-d]@CD[-c]@RicciCD[-a,-b],CD[-e]@RiemannCD[h1,h2,h3,h4]CD[-h4]@CD[-h2]@RicciCD[-h1,-h3]},UseSymmetries->True,MetricOn->All];

This rule has a typo in MathTensor's book. I assume the index c on the RHS is actually f:

RuleD11bisTwoRiemann=MakeRule[{CD[-g]@RiemannCD[-a,-b,-c,-d]CD[g]@RiemannCD[c,e,d,f],1/2CD[-h3]@RiemannCD[-h1,-h2,-a,-b]CD[h3]@RiemannCD[h1,h2,f,e]},UseSymmetries->True,MetricOn->All];

RiemannRule34:

CD[-e]@RiemannCD[a,f,b,c]RiemannCD[-a,-b,-c,-d]/.RuleD01TwoRiemann

RiemannRule35:

CD[-e]@RiemannCD[-a,-b,-c,-d]CD[-f]@RiemannCD[a,g,b,c]/.RuleD11TwoRiemann

RiemannRule36:

CD[-e]@RiemannCD[a,c,b,f]RiemannCD[-a,-b,-c,-d]/.RuleD01TwoRiemann

RiemannRule37:

CD[-e]@RiemannCD[-a,-b,-c,-d]CD[-f]@RiemannCD[a,c,b,g]/.RuleD11TwoRiemann

RiemannRule38:

CD[-g]@RiemannCD[-a,-b,-c,-d]CD[g]@RiemannCD[c,e,d,f]/.RuleD11bisTwoRiemann

RiemannRule39:

CD[-e]@RiemannCD[a,d,b,c]CD[-d]@CD[-c]@RicciCD[-a,-b]/.RuleD2RicciDRiemann

RiemannRule40:

CD[-d]@CD[-c]@RicciCD[-a,-b]RiemannCD[a,d,b,c]/.RuleD2RicciRiemann

Finally we study six rules with products of three curvature tensors, only four of which are independent:
	RuleRicciTwoRiemann
	RuleThreeRiemann
	RuleThreeRiemannB
	RuleThreeRiemannC

Constructing the independent rules takes rather long because there are lots of equivalent cases:

AbsoluteTiming[Length[RuleRicciTwoRiemann=MakeRule[{RicciCD[-a,-b]RiemannCD[-c,-d,-e,a]RiemannCD[b,c,d,e],-1/2RicciCD[-h1,-h2]RiemannCD[-h3,-h4,-h5,h1]RiemannCD[h2,h5,h3,h4]},UseSymmetries->True,MetricOn->All]]]

AbsoluteTiming[Length[RuleThreeRiemann=MakeRule[{RiemannCD[-a,-b,-c,-d]RiemannCD[-e,-f,a,b]RiemannCD[c,e,d,f],1/2RiemannCD[-h1,-h2,-h3,-h4]RiemannCD[-h5,-h6,h1,h2]RiemannCD[h3,h4,h5,h6]},UseSymmetries->True,MetricOn->All]]]

AbsoluteTiming[Length[RuleThreeRiemannB=MakeRule[{RiemannCD[-a,-b,-c,-d]RiemannCD[-e,a,-f,b]RiemannCD[c,e,d,f],1/4RiemannCD[-h1,-h2,-h3,-h4]RiemannCD[-h5,-h6,h3,h4]RiemannCD[h1,h2,h5,h6]},UseSymmetries->True,MetricOn->All]]]

AbsoluteTiming[Length[RuleThreeRiemannC=MakeRule[{RiemannCD[-a,-b,-c,-d]RiemannCD[-e,a,-f,c]RiemannCD[b,f,d,e],RiemannCD[-h1,-h2,-h3,-h4]RiemannCD[-h5,h1,-h6,h3]RiemannCD[h2,h5,h4,h6]-1/4RiemannCD[-h1,-h2,-h3,-h4]RiemannCD[-h5,-h6,h1,h2]RiemannCD[h3,h4,h5,h6]},UseSymmetries->True,MetricOn->All]]]

RiemannRule10:

RicciCD[-a,-b]RiemannCD[-c,-d,-e,a]RiemannCD[b,c,d,e]/.RuleRicciTwoRiemann

RiemannRule11:

RiemannCD[-a,-b,-c,-d]RiemannCD[-e,-f,a,b]RiemannCD[c,e,d,f]/.RuleThreeRiemann

RiemannRule12:

RiemannCD[-a,-b,-c,-d]RiemannCD[-e,a,-f,b]RiemannCD[c,e,d,f]/.RuleThreeRiemannB

RiemannRule13:

RiemannCD[-a,-b,-c,-d]RiemannCD[-e,a,-f,b]RiemannCD[c,d,e,f]/.RuleThreeRiemann

RiemannRule14:

RiemannCD[-a,-b,-c,-d]RiemannCD[-e,a,-f,c]RiemannCD[b,d,e,f]/.RuleThreeRiemannB

RiemannRule15:

RiemannCD[-a,-b,-c,-d]RiemannCD[-e,a,-f,c]RiemannCD[b,f,d,e]/.RuleThreeRiemannC

Remove rules:

Remove[RuleBoxRiemann,RuleD01TwoRiemann,RuleD11bisTwoRiemann,RuleD11TwoRiemann,RuleD2RicciDRiemann,RuleD2RicciRiemann,RuleD4Ricci,RuleD4RicciScalar,RuleDivGradRicci,RuleDivGradRiemann,RuleDivRiemann,RuleRicciTwoRiemann,RuleThreeRiemann,RuleThreeRiemannB,RuleThreeRiemannC,RuleTwoRiemann2,RuleTwoRiemann3]

10. Manipulation of expressions

10.1. Find indices

Given an expression, very often we want to extract all indices of the expression, or those of certain type, or character or state. We have already seen how to select indices of any type and any character. In this subsection we shall define functions to extract all indices according to a number of ` `selectors´ ´.
There are two sets of functions. There is first the internal function FindIndices, which returns all indices of an expression, without first evaluating it, and its relatives FindFreeIndices, FindDummyIndices and FindBlockedIndices. There is then the user driver for finding indices, called IndicesOf, which admits a number of criteria for index selection, and evaluates its input (as the user typically expects).

Define some more objects:

DefTensor[r[],M3]
DefTensor[S[LI[_],a,b],M3,Symmetric[{2,3}]]
DefTensor[n[a],M3]

$DefInfoQ=True;

Let us start with this expression

expr=4T[a,b,-c,Dir[v[d]]]v[-a]r[]^2CD[-e][S[LI[1],-b,d]]

The three Find*Indices functions have attribute HoldFirst, and hence require the use of Evaluate when we do not feed the expression itself:

FindIndices[Evaluate[expr]]

FindFreeIndices[Evaluate[expr]]

FindDummyIndices[Evaluate[expr]]

FindBlockedIndices[Evaluate[expr]]

This attribute is given to be able to extract indices of expressions before they get evaluated:

n/:n[a_]n[-a_]:=-1

n[a]n[-a]

FindIndices[n[a]n[-a]]

Note that the returned lists of indices have head IndexList, in order to avoid confusion with basis and component indices:

InputForm[%]

Finding indices. Internal functions

Then we have the user driver, with two pairs of brackets:

With no selectors, we get all indices in the expression:

IndicesOf[][expr]

We can select indices by type:

IndicesOf[AIndex][expr]

IndicesOf[DIndex][expr]

IndicesOf[LIndex][expr]

or by character:

IndicesOf[Up][expr]

IndicesOf[Down][expr]

or by state:

IndicesOf[Free][expr]

IndicesOf[Dummy][expr]

IndicesOf[Blocked][expr]

We can also specify a vbundle (or a basis):

IndicesOf[TangentM3][expr]

or a tensor:

IndicesOf[S][expr]

We can combine several selectors, in two possible ways: a list of selectors represents the union, and a sequence represents the intersection.

Directions or labels:

IndicesOf[{DIndex,LIndex}][expr]

Free indices on the tensor S:

IndicesOf[Free,S][expr]

Members of a selectors list are conmutative, but not members of a sequence: the construction of the final list is performed from left to right. Compare with the previous result:

IndicesOf[S,Free][expr]

This powerful function is very flexible and has more possibilities once we allow the presence of bases (see xCobaDoc.nb).

Finding indices. User driver

10.2. Replace indices

A second, commonly required operation is the change of indices. In principle, the user should never need to replace indices directly in an expression. But sometimes it is needed, either because two expressions with different free indices must be compared, or to avoid problems with the dummy indices. This subsection introduces the basic function ReplaceIndex, and five more functions based on it.
Replacing indices is a basic critical operation: done without care would inmediately lead to wrong results. Whenever possible, we recommend the use of delta contractions to change indices, which is a safer and mathematically clearer operation.

Above all, avoid the temptation of changing indices using normal rules, unless you are sure that index-collisions will not appear.

Index replacements.

Let us go back to our expression:

expr

For similar reasons, ReplaceIndex has attribute HoldFirst, and hence we need Evaluate. We can replace a single index:

ReplaceIndex[Evaluate[expr],a->f]

Note that the paired index -a has not been replaced. We need separate rules for that (this gives greater flexibility):

ReplaceIndex[Evaluate[expr],{a->f,-a->-f}]

We can change the character of an index, or even the type:

ReplaceIndex[Evaluate[expr],{a->-f,-a->f,Dir[v[d]]->LI[2]}]

If the rules do not apply, nothing happens:

ReplaceIndex[Evaluate[expr],{f->a}]

Pattern indices can be changed, but this should never be required:

ReplaceIndex[T[a_],a->b]

The replacement is performed with great care. In the following list T is a tensor, but X is not a tensor:

ReplaceIndex[{a,T[a],T[-a],X[a]},a->b]

As I said, the user should try to refrain from using ReplaceIndex. If a free index must be changed, use contractions with delta tensors: if we want to change from T[a] to T[b] we can do this:

T[a]delta[-a,b]

In some cases we need to change the dummy indices of an expression. For example when we want to multiply a previous expression by itself. ReplaceDummies is the function which calls ReplaceIndex to manipulate all possible cases with dummies. A particular case, always required in the canonicalization process, is implemented in SameDummies.

Suppose now this scalar expression:

expr=T[a,b]T[-a,-b]

This is incorrect:

expr expr

The following is correct. ReplaceDummies changes the dummy by the given indices

expr ReplaceDummies[expr,IndexList[e,f]]

If not enough indices (or no index at all) are given, the dummies are replaced by unique  dollar-indices:

expr ReplaceDummies[expr]

IndicesOf[][%]

A particular case of dummy manipulation is the minimization of the number of dummies in a sum of terms. This is achieved using SameDummies, which is simply a combined call to FindDummyIndices and ReplaceDummies.

SameDummies[T[a,b,d,-d]v[-b]+T[a,c,e,-e]v[-c]]

In parallel, sometimes we need to change the free indices of an expression. This happens for example when we need to compare the results of several independent computations. If we do not know the index translation table which brings one expression to the other, then we need to search among all permutations, what could take quite some time...

This changes the free indices a,c to b,e. The index b was a dummy in the original expression, and hence it must be replaced (by a unique dollar-dummy) before actually changing the free indices:

ChangeFreeIndices[T[a,b,c,d,-b,-d],{b,e}]

It is not obvious to see whether these two expressions are the same or not under a renaming of free indices:

expr1=4U[-a,-b,-c]U[b,-d,-e]T[c,d]
expr2=2U[-d,-b,-a]U[a,-f,g]T[d,-g]

The function EqualExpressionsQ gives True if one is simply a constant times the other, and False otherwise:

EqualExpressionsQ[expr1,expr2]

EqualExpressionsQ[expr1,T[]expr2]

Finally, we introduce the operation of index splitting. The key idea is the concept of a splitrule of the form   a -> IndexList[ a1, a2, ... ].  The function SplitIndex is rather efficient because it uses ReplaceIndex only once, replacing the index by a temporary variable, which is later used to do all other replacements with ReplaceAll, which is much faster.

We can construct a list of replacements as follows:

SplitIndex[expr1,-a->IndexList[-A,-B,-C]]

or an array:

SplitIndex[expr1,{-a->IndexList[-A,-B],-b->IndexList[C,D]}]

Dimensions[%]

There is also a combined notation for multiple expansions over the same range. Again, this is performed through a single call to ReplaceIndex.

SplitIndex[expr1,IndexList[-a,-b]->IndexList[LI[1],LI[2],LI[3]]]

Remove[expr1,expr2]

10.3. Scalar and scalar functions

Another important problem is the following. In the same way that we are using Times for a product of tensors, we shall use Power, Sqrt, Exp and other functions acting on scalar tensor expressions. There is a problem: there are built-in definitions for those functions that could be incorrect for tensors. For example we have (a b)^n = a^n b^n. However (v[a] v[-a])^n is not equal to v[a]^n v[-a]^n. To prevent those problems we have two options: we could define new functions IndexPower, IndexSqrt, etc, or we can introduce a new head Scalar such that the arguments of those functions must be always wrapped with that head.

Scalar expressions.

The direct use of Power could lead to unconsistent expressions:

(v[a]v[-a])^2

Instead, the right way to write it would be

Scalar[v[a]v[-a]]^2

xTensor` can work with those objects as usual:

Scalar[v[a]v[-a]]^2 Scalar[v[b]v[-b]]^3/Scalar[v[c]v[-c]]^4

Simplification[%]

Note that Scalar effectively shields its argument from the rest of an expression and therefore there can be repeated indices:

v[a]Scalar[v[a]v[-a]]

Validate[%]

By default the function Simplification does not add new Scalar heads, because that takes some time for large expressions:

v[b]v[-b]/Scalar[v[a]v[-a]]//Simplification

You can force it using PutScalar:

PutScalar[%]

Simplification[%]

Two other functions to manipulate Scalar expressions are BreakScalars and NoScalar:

Scalar[v[a]v[-a]v[b]v[-b]]

%//BreakScalars

%//Simplification

%//NoScalar

Definition of a scalar function.

A scalar function is one that only accepts scalar arguments and returns scalars.

$ScalarFunctions

We define a scalar function SF:

DefScalarFunction[SF]

ScalarQ[SF[v[a]v[-a]]]

$CovDs

CD[-a][1/SF[Scalar[v[a]v[-a]]]]

10.4. Inert heads

Sometimes we need to wrap tensors with some head. This is an extremely general concept and it is prepared to play the role of an arbitary head with the only property of being "transparent" to the symmetry and canonicalization routines. Dagger and ERROR are both inert heads.

Definition of an inert head.

A prototypical example could be traceless, transverse tensors: We define an inert head called TTPart:

DefInertHead[TTPart]

TTPart[U[-b,-c,-a]]//ToCanonical

UpVectorQ[TTPart[v[a]]]

UndefInertHead[TTPart]

With help from Guillaume Faye, this concept has been improved in version 0.9.2 of xTensor` in three different ways:

1) There is an option LinearQ to make the inert head linear with respect to constants. The default value is False.

DefInertHead[lhead,LinearQ->True]

LinearQ[lhead]

lhead[3T[a,b]+S[a,b]]

2) There is an option ContractThrough to specify a list of metrics, projectors and/or the delta tensor which can be contracted through the inert head. In this example only the delta tensor can be contracted:

DefInertHead[chead,ContractThrough->{delta}]

DefInertHead[mhead,ContractThrough->{delta,metricg}]

?mhead

Compare now these cases:

delta[a,-b]{lhead[T[b]],chead[T[b]],mhead[T[b]]}

ContractMetric[metricg[-a,-b]{lhead[T[b]],chead[T[b]],mhead[T[b]]}]

3) It is possible to have additional arguments, including arguments with synchronized lists of indices. This does not even have to be specified as an argument. The indices must be surrounded with the head IndexList. The system keeps a one to one relation between the indices in the tensorial expression and the indices in the lists of arguments. This correspondence is based on the slots: for example in the following expression the IndexList marks the {3, 2} slots of the tensor U, and this is kept through the use of the xTensor` functions:

expr=U[-a,-b,d]mhead[U[b,c,a],{a,c},IndexList[a,c]];

ReplaceIndex[Evaluate[expr],{c->f}]//InputForm

ToCanonical[expr]//InputForm

metricg[-f,-c]expr//ContractMetric//InputForm

Perhaps there should be something more general than slot-syncronization, but we do not know of a complete set of possibilities, and no other case has been asked for until now.

UndefInertHead[lhead]

UndefInertHead[chead]

UndefInertHead[mhead]

10.5. Monomials

A monomial is defined as a product of tensors such that it cannot be further reduced into products with different dummies. The function BreakInMonomials decomposes an arbitrary expression in monomials with (inert) head Monomial.

We break this expressions in monomials. Note that scalars are left outside:

7r[]^2T[a,b,c]U[-b,-c,-e]v[e]T[f,g,h]v[-g]v[-h]

BreakInMonomials[%]

Now for example each monomial could be canonicalized independently:

%/.expr_Monomial:>ToCanonical[expr]

Monomials.

10.6. Symmetrization of an expression. V1.1.0 UPDATE

Symmetrization.

Symmetrization of an expression under a given group of permutations. Numbers on the third argument refer to positions in the second argument of ImposeSymmetry. Note that ToCanonical is always applied after the symmetrization functions in order to take into account the antisymmetry of the tensor U. It is not automatically implemented into ImposeSymmetry because for large numbers of indices it would take ages to canonicalize all terms!

ImposeSymmetry[U[-a,-b,-c],{-a,-b},Symmetric[{1,2}]]//ToCanonical

ImposeSymmetry[U[-a,-b,-c],{-a,-b},Antisymmetric[{1,2}]]//ToCanonical

ImposeSymmetry[U[-a,-b,-c]U[-d,-e,-f],{-a,-d},Antisymmetric[{1,2}]]//ToCanonical

xAct`xTensor`Symmetrize[U[-a,-b,-c]U[-d,-e,-f],{-a,-d}]//ToCanonical

xAct`xTensor`Symmetrize[U[-a,-b,-c]U[-d,-e,-f],{-e,-d}]//ToCanonical

Antisymmetrize[U[-a,-b,-c]U[-d,-e,-f],{-a,-d}]//ToCanonical

Antisymmetrize[U[-a,-b,-c]U[-d,-e,-f],{-e,-d}]//ToCanonical

ps=PairSymmetrize[U[-a,-b,-c]U[-d,-e,-f],{{-a,-b},{-c,-d},{-e,-f}}]//ToCanonical

pa=PairAntisymmetrize[U[-a,-b,-c]U[-d,-e,-f],{{-a,-b},{-c,-d},{-e,-f}}]//ToCanonical

ps-ReplaceIndex[Evaluate[ps],{-a->-c,-b->-d,-c->-a,-d->-b}]//Simplification

pa+ReplaceIndex[Evaluate[pa],{-a->-c,-b->-d,-c->-a,-d->-b}]//Simplification

Working with Gdelta implements certain identities involving many terms, which are even unknown to ToCanonical. For example, for a Gdelta with antisymmetric blocks of n indices, antisymmetrizing over more than n indices gives zero. ToCanonical does not know it, but the contraction rules encode that information. We develop this example in arbitrary dimension.

Define a manifold with metric in arbitrary dimension:

$DefInfoQ=$UndefInfoQ=False;
DefConstantSymbol[dim]
DefManifold[Mdim,dim,IndexRange[\[FormalA],\[FormalZ]]]
DefMetric[1,metric[-\[FormalA],-\[FormalB]],cd,PrintAs->"g"]

Take a Gdelta, with antisymmetric blocks of 3 indices:

Gdelta[-\[FormalA],-\[FormalB],-\[FormalC],\[FormalD],\[FormalE],\[FormalF]]

Antisymmetrize in four indices, two in each block:

Antisymmetrize[%,{-\[FormalA],-\[FormalB],\[FormalE],\[FormalF]}]

ToCanonical cannot see that this is zero:

%//ToCanonical

We prove it by expansion into elementary deltas:

%//ExpandGdelta

%//ToCanonical

However, a different way of antisymmetrizing is contracting with a Gdelta, and this gives zero 
immediately, independently of the dimension:

Gdelta[-\[FormalA],-\[FormalB],-\[FormalC],\[FormalD],\[FormalE],\[FormalF]]Gdelta[\[FormalA],\[FormalB],-\[FormalE],-\[FormalF],\[FormalG],\[FormalH],\[FormalI],\[FormalJ]]

Clean up:

UndefMetric[metric]
UndefManifold[Mdim]
UndefConstantSymbol[dim]
$DefInfoQ=$UndefInfoQ=True;

10.7. Symmetric Trace-Free part

The command STFPart returns the symmetric, trace-free part of an expression with only free indices. In order for a symmetric object to have a trace we need to use a metric.

This is the STF part of a 3-index tensor without symmetry, and with respect to the metricg:

STFPart[T[a,b,c],metricg]

The result is not simplified at all. We must ask for explicit metric contraction and simplification:

%//ContractMetric//Simplification

10.8. Mapping at specified positions

We use the package ExpressionManipulation` (written by David Park & Ted Ersek) to have full control on the positions of an expression

For versions of xTensor` before 0.9.8 this package was automatically loaded, but now we need to do it explicitly, to avoid confusing it with the xAct` packages. This package has no copyright messages:

<<xAct`ExpressionManipulation`

Suppose we have an expression like this, which is clearly 0 because U is antisymmetric and S is symmetric:

DefTensor[G[a,b],M3,Symmetric[{1,2}]]

expr=U[a,b,c]v[-a]v[-b]v[-c]+U[a,b,c]v[-a]G[-b,-c]+T[-b,-c]v[b]v[c]-T[-d,-e]v[d]v[e]

Simplification[expr]

Each object in the expression can be identified in a tree-like form by its position. In this case there are four positions at level 1:

expr//ColorTerms

We can evaluate a given function at one or several positions:

MapAt[Simplification,expr,{2}]

MapAt[Simplification,expr,{{1},{2}}]

In order to map a function over several terms simultaneously we need to make use of the concept of ExtendedPosition (also from the ExpressionManipulation` package) and the new function NewMapAt. See notes for ExtendedPosition and eP.

MapAt[f,expr,{{3},{4}}]

NewMapAt[f,expr,eP[{},{3,4}]]

We can also localize a given pattern:

expr//ColorPositionsOfPattern[U[__]]

or map a function onto the occurrences of that pattern:

Position[expr,U[__]]

MapAt[f,expr,%]

In case the expression is held, we can use EvaluateAt (also from the package ExpressionManipulation`) instead of MapAt:

Hold[U[a,b,c]v[-a]v[-b]v[-c]+U[a,b,c]v[-a]S[-b,-c]+T[-b,-c]v[b]v[c]-T[-d,-e]v[d]v[e]+1+2]

MapAt[Simplification,%,{1,1}]

EvaluateAt[{1,1},Simplification][%%]

10.9. Coefficients

Given an expression, sometimes we need to get the (possibly indexed) coefficient multiplying a given tensor. When the objects have symmetries the coefficient might be more complicated than expected.

Coefficients of indexed expressions.

The basic relation is this:

IndexCoefficient[T[a,b],T[c,d]]

Metrics and symmetries must be taken into account:

IndexCoefficient[T[-a,-b],T[c,d]]

IndexCoefficient[metricg[-a,-b],metricg[-c,-d]]

IndexCoefficient[T[a,b]v[-a]v[-b],v[-c]]//Simplification

IndexCoefficient[T[a,b]v[-a]v[-b],v[-c]v[-d]]

Imitating Coefficient, dependencies must be explicit:

IndexCoefficient[T[a,b],metricg[a,b]]

Frequently we need further manipulation of the result:

IndexCoefficient[T[a,b]v[-b]v[c]+T[d,-d]v[a]v[c],v[e]]

%//ContractMetric//Simplification

A closely related function is IndexCollect, which acts as a simple recursive driver for IndexCoefficient. This pair has been designed to follow closely the Mathematica pair Collect / Coefficient.

Suppose an expression like this:

expr=T[a,b,-c]v[c]+U[a,b,c,d]v[-c]v[-d]+metricg[a,b]

IndexCollect[expr,{v[d],v[e]},Simplification]

10.10. Tensor equations

It is not simple to deal with tensor equations because there are many things to worry about simultaneously. Currently xTensor` only works with linear equations where the x tensor is not contracted. The concept of equality (==) has been generalized to include tensor equalities.

From the point of view of Mathematica these two expressions are not related:

T[a,b]v[-b]==T[a,c]v[-c]

However xTensor` recognizes them as equal:

Simplification[%]

These two are clearly different:

Simplification[T[a,b]v[-b]==T[-a,c]v[-c]]

Solving equations. The answer is a tensor rule. Options to IndexSolve are actually options to MakeRule:

rule=IndexSolve[v[a]v[-a]T[b,c,d]==v[b]v[c]v[a]S[-a,d],T[b,c,d], MetricOn->{b},ContractMetrics->True]

{T[a,b,c],T[-a,b,c],T[a,-b,c],T[-a,-b,c]}/.rule

10.11. Clean up

UndefMetric[metricg]

UndefTensor/@$Tensors;

UndefManifold[S2]

UndefManifold[M3]

UndefScalarFunction[SF]

11. Final comments

11.1. ToDo list. Version

There are many things that xTensor` cannot do yet. Here there is a list with some of them:
	- Differential forms
	- Grasmann variables

Note: For further information about xTensor`, and to be kept informed about new releases, you may contact the author electronically at jose@xact.es. This is xTensorDoc.nb, the docfile of xTensor`, currently in version 1.1.3.

The latest version of the package can always be downloaded from   http://www.xAct.es/

11.2. Frequently Asked Questions

1. Why is delta contraction automatic? I want to have full control on every single detail of the process.
    Why is metric contraction not automatic? It is boring to ask continuously for trivial tasks.

There are two opposite approaches to the issue of automatic simplification: Too little automatic simplification makes the computations hard and slow, while too much automatic simplification makes the process too restrictive and the result seems just magic. Following Mathematica, we have tried to automatize all those things which are commonly required, like contraction of any tensor with delta, or the expansion of derivatives of products using the Leibnitz rule. However, most of the processes are not automatic: simplification, metric contraction, commutation of derivatives, etc.
There are some processes for which the user can choose whether they are automatic or not: look for the pairs commandStart / commandStop.

2. I want to define several tensors at the same time. Why not allowing this syntax? 
			DefTensor[{T[a], S[b,c], R[-a,-b,-c,-d}, M]

We do not include new definitions which do not save time or space with respect to their Mathematica counterparts. This helps learning Mathematica, and save space in xTensor`. In that particular example we recommend to use 
		Map[ DefTensor[#, M]&, {T[a], S[b,c], R[-a,-b,-c,-d]} ]
or simply threading
		Thread[ DefTensor[{T[a], S[b,c], R[-a,-b,-c,-d}, M] ]

3. In MathTensor many commands have an additional argument to map them only to a certain part of an expression. Why don't you have the same thing?

The answer is identical to that of question 2. Use the syntax MapAt[f, expr, n] to map the function f on the n-th term of expr.

4. xTensor is great for abstract tensor calculus, but I want to do component computations. Can I do that with xTensor?

No, there is no support in xTensor for component tensor computer algebra. The twin package xCoba is being developed for that purpose, but it does not have the degree of maturity of xTensor.

5. I want to do 3+1 decompositions (ADM-type). xTensor has a four-dimensional approach to that. Can I do the same with the traditional {i, 0} index splitting formalism?

Yes, this is possible, but there are no predefined tools to help in the process. You need to code them up using the building blocks: define three manifolds: 4d, 3d and 1d. Then use SplitIndex and TraceProductDummy for the computations.

11.3. Statistics

Collection of public symbols in xTensor`:

Names["xAct`xTensor`*"]

Length[%]

TimeUsed[]

MemoryInUse[]/1000/1024//N

MaxMemoryUsed[]/1000/1024//N