Tensor calculus with open-source software: the SageManifolds project

The SageManifolds project aims at extending the mathematics software system Sage towards differential geometry and tensor calculus. Like Sage, SageManifolds is free, open-source and is based on the Python programming language. We discuss here some details of the implementation, which relies on Sage's parent/element framework, and present a concrete example of use.


Introduction
Computer algebra for general relativity (GR) has a long history, which started almost as soon as computer algebra itself in the 1960s.The first GR program was GEOM, written by J.G. Fletcher in 1965 [1].Its main capability was to compute the Riemann tensor of a given metric.In 1969, R.A. d'Inverno developed ALAM (for Atlas Lisp Algebraic Manipulator ) and used it to compute the Riemann and Ricci tensors of the Bondi metric.According to [2], the original calculations took Bondi and collaborators 6 months to finish, while the computation with ALAM took 4 minutes and yielded the discovery of 6 errors in the original paper.Since then, numerous packages have been developed: the reader is referred to [3] for a review of computer algebra systems for GR prior to 2002, and to [4] for a more recent review focused on tensor calculus.It is also worth to point out the extensive list of tensor calculus packages maintained by J. M. Martin-Garcia at [5].
As far as tensor calculus is concerned, the above packages can be distinguished by the type of computation that they perform: abstract calculus (xAct/xTensor, Ricci, Cadabra, Redberry), or component calculus (xAct/xCoba, DifferentialGeometry, GRTensorII, Atlas 2).In the first category, tensor operations such as contraction or covariant differentiation are performed by manipulating the indices themselves rather than the components to which they correspond.In the second category, vector frames are explicitly introduced on the manifold and tensor operations are carried out on the components in a given frame.
3. An overview of Sage Sage [14] is a free, open-source mathematics software system, which is based on the Python programming language.It makes use of over 90 open-source packages, among which are Maxima and Pynac (symbolic calculations), GAP (group theory), PARI/GP (number theory), Singular (polynomial computations), and matplotlib (high quality 2D figures).Sage provides a uniform Python interface to all these packages; however, Sage is much more than a mere interface: it contains a large and increasing part of original code (more than 750,000 lines of Python and Cython, involving 5344 classes).Sage was created in 2005 by W. Stein [15] and since then its development has been sustained by more than a hundred researchers (mostly mathematicians).Very good introductory textbooks about Sage are [16,17,18].
Apart from the syntax, which is based on a popular programming language and not a custom script language, a difference between Sage and, e.g., Maple or Mathematica is the usage of the parent/element pattern.This framework more closely reflects actual mathematics.For instance, in Mathematica, all objects are trees of symbols and the program is essentially a set of sophisticated rules to manipulate symbols.On the contrary, in Sage each object has a given type (i.e. is an instance of a given Python class1 ), and one distinguishes parent types, which model mathematical sets with some structure (e.g.algebraic structure), from element types, which model set elements.Moreover, each parent belongs to some dynamically generated class that encodes informations about its category, in the mathematical sense of the word (see [19] for a discussion of Sage's category framework).Automatic conversion rules, called coercions, prior to a binary operation, e.g.x + y with x and y having different parents, are implemented.

Aim and scope
Sage is well developed in many areas of mathematics but very little exists for differential geometry and tensor calculus.One may mention differential forms defined on a fixed coordinate patch, implemented by J. Vankerschaver [20], and the 2-dimensional parametrized surfaces of the 3dimensional Euclidean space recently added by M. Malakhaltsev and J. Vankerschaver [21].
The aim of SageManifolds [22] is to introduce smooth manifolds and tensor fields in Sage, with the following requirements: (i) one should be able to introduce various coordinate charts on a manifold, with the relevant transition maps; (ii) tensor fields must be manipulated as such and not through their components with respect to a specific (possibly coordinate) vector frame.
Concretely, the project amounts to creating new Python classes, such as Manifold, Chart, TensorField or Metric, to implement them within the parent/element pattern and to code mathematical operations as class methods.For instance the class Manifold, devoted to real smooth manifolds, is a parent class, i.e. it inherits from Sage's class Parent.On the other hand, the class devoted to manifold points, ManifoldPoint, is an element class and therefore inherits from Sage's class Element.This is illustrated by the inheritance diagram of diagram, each class at the base of some arrow is a subclass (also called derived class) of the class at the arrowhead.Note however that the actual type of a parent is a dynamically generated class taking into account the mathematical category to which it belongs.For instance, the actual type of a smooth manifold is not Manifold, but a subclass of it named Manifold with category, reflecting the fact that Manifold is declared in the category of Sets2 .Note also that the class Manifold inherits from Sage's class UniqueRepresentation, which ensures that there is a unique manifold instance for a given dimension and given name.

Implementation of charts
Given a smooth manifold M of dimension n, a coordinate chart on some open subset U ⊂ M is implemented in SageManifolds via the class Chart, whose main data is a n-tuple of Sage symbolic variables (x 1 , . . ., x n ), each of them representing a coordinate.In general, more than one (regular) chart is required to cover the entire manifold.For instance, at least 2 charts are necessary for the n-dimensional sphere S n (n ≥ 1) and the torus T 2 and 3 charts for the real projective plane RP 2 (see Fig. 6 below).Accordingly, SageManifolds allows for an arbitrary number of charts.To fully specify the manifold, one shall also provide the transition maps (changes of coordinates) on overlapping chart domains (SageManifolds class CoordChange).

Implementation of scalar fields
A scalar field on manifold M is a smooth mapping  where U is some open subset of M. A scalar field has different coordinate representations F , F , etc. in different charts X, X, etc. defined on U : These representations are stored in some attribute of the class ScalarField, namely a Python dictionary3 whose keys are the various charts defined on U : f. express = X : F, X : F , . . . .
Each representation F is an instance of the class FunctionChart, which resembles Sage native symbolic functions, but involves automatic simplifications in all arithmetic operations.Given an open subset U ⊂ M, the set C ∞ (U ) of scalar fields defined on U has naturally the structure of a commutative algebra over R: it is clearly a vector space over R and it is endowed with a commutative ring structure by pointwise multiplication: The algebra C ∞ (U ) is implemented in SageManifolds via the parent class ScalarFieldAlgebra, in the category CommutativeAlgebras.
The corresponding element class is of course ScalarField (cf.Fig. 2).

Modules and free modules
Given an open subset U ⊂ M, the set X (U ) of all smooth vector fields defined on U has naturally the structure of a module over the algebra C ∞ (U ).Let us recall that a module is similar to a vector space, except that it is based on a ring (here C ∞ (U )) instead of a field (usually R or C in physical applications).Of course, every vector space is a module, since every field is a ring.There is an important difference though: every vector space has a basis (as a consequence of the axiom of choice), while a module does not necessarily have any.When it possesses one, it is called a free module.Moreover, if the module's base ring is commutative, it can be shown that all bases have the same cardinality, which is called the rank of the module (for vector spaces, which are free modules, the word dimension is used instead of rank ).
If X (U ) is a free module (examples are provided in Sec.4.5 below), a basis of it is nothing but a vector frame (e a ) 1≤a≤n on U (often called a tetrad in the context of 4-dimensional GR): (5) The rank of X (U ) is thus n, i.e. the manifold's dimension 4 .At any point p ∈ U , Eq. ( 5) gives birth to an identity in the tangent vector space T p M: which means that the set (e a (p)) 1≤a≤n is a basis of T p M. Note that if U is covered by a chart (x a ) 1≤a≤n , then (∂/∂x a ) 1≤a≤n is a vector frame on U , usually called coordinate frame or natural basis.Note also that, being a vector space over R, the tangent space T p M represents another kind of free module which occurs naturally in the current context.It turns out that so far only free modules with a distinguished basis were implemented in Sage.This means that, given a free module M of rank n, all calculations refer to a single basis of M .This amounts to identifying M with R n , where R is the ring over which M is based.This is unfortunately not sufficient for dealing with smooth manifolds in a coordinate-independent way.For instance, there is no canonical isomorphism between T p M and R n when no coordinate system is privileged in the neighborhood of p. Therefore we have started a pure algebraic part of SageManifolds to implement generic free modules, with an arbitrary number of bases, none of them being distinguished.This resulted in (i) the parent class FiniteRankFreeModule, within Sage's category Modules, and (ii) the element class FiniteRankFreeModuleElement.Then both classes VectorFieldFreeModule (for X (U ), when it is a free module) and TangentSpace (for T p M) inherit from FiniteRankFreeModule (see Fig. 3).

Implementation of vector fields
Ultimately, in SageManifolds, vector fields are described by their components with respect to various vector frames, according to Eq. ( 5), but without any vector frame being privileged, leaving the freedom to select one to the user, as well as to change coordinates.A key point is that not every manifold admits a global vector frame.A manifold M, or more generally an open subset U ⊂ M, that admits a global vector frame is called parallelizable.Equivalently, M is parallelizable if, and only if, X (M) is a free module.In terms of tangent bundles, parallelizable manifolds are those for which the tangent bundle is trivial: Examples of parallelizable manifolds are [23] • the Cartesian space R n for n = 1, 2, . .., • the sphere S 3 SU(2), as any Lie group, • the sphere S 7 , • any orientable 3-manifold (Steenrod theorem [24]).• the sphere S 2 (as a consequence of the hairy ball theorem), as well as any sphere S n with n ∈ {1, 3, 7}, • the real projective plane RP 2 .
Actually, "most" manifolds are non-parallelizable.As noticed above, if a manifold is covered by a single chart, it is parallelizable (the prototype being R n ).But the reverse is not true: S 1 and T 2 are parallelizable and require at least two charts to cover them.
If the manifold M is not parallelizable, we assume that it can be covered by a finite number N of parallelizable open subsets U i (1 ≤ i ≤ N ).In particular, this holds if M is compact, for any compact manifold admits a finite atlas.We then consider the restrictions of vector fields to the U i 's.For each i, X (U i ) is a free module of rank n = dim M and is implemented in SageManifolds as an instance of VectorFieldFreeModule (cf.Sec.4.4 and Figs. 3 and 4).Each vector field v ∈ X (U i ) has different sets of components (v a ) 1≤a≤n in different vector frames (e a ) 1≤a≤n introduced on U i [cf.Eq. ( 5)].They are stored as a Python dictionary whose keys are the vector frames: v. components = {(e) : (v a ), (ê) : (v a ), . ..} . (7)

Implementation of tensor fields
The implementation of tensor fields in SageManifolds follows the strategy adopted for vector fields.Consider for instance a tensor field T of type (1,1) on the manifold M. It can which defines T a b , is meaningful only when a vector frame (e a ) exists5 .Therefore, one first decomposes the tensor field T into its restrictions T | U i on parallelizable open subsets of M, U i (1 ≤ i ≤ N ) and then considers the components on various vector frames on each subset U i .For each vector frame (e a ), the set of components (T a b ) is stored in a devoted class (named Components), which takes into account all the tensor monoterm symmetries: only non-redundant components are stored, the other ones being deduced by (anti)symmetry.This is illustrated in Fig. 5, which depicts the internal storage of tensor fields in SageManifolds.Note that each component T a b is a scalar field on U i , according to the formula Accordingly, the penultimate level of Fig. 5 corresponds to the scalar field storage, as described by (3).The last level is constituted by Sage's symbolic expressions (class Expression).

Parallelization
To improve the reactivity of SageManifolds and take advantage of multicore processors, some tensorial operations are performed by parallel processes.The parallelization is implemented by means of the Python library multiprocessing, via the built-in Sage decorator @parallel.Using it permits to define a function that is run on different sub-processes.If n processes are used, given a function and a list of arguments for it, any process will call the function with an element of the list, one at time, spanning all the list.Currently6 , the parallelized operations are tensor algebra, tensor contractions, computation of the connection coefficient and computation of Riemann tensor.
The parallelization of an operation is achieved by first creating a function which computes the required operation on a subset of the components of a tensor; second, by creating a list of 2n (twice the number of used processes) arguments for this function.Then applying this function to the input list, the calculation is performed in parallel.At the end of the computation a fourth phase is needed to retrieve the results.The choice to divide the work in 2n is a compromise between the load balancing and the cost of creating multiple processes.The number of processors to be used in the parallelization can be controlled by the user.

SageManifolds at work: Kerr spacetime and Simon-Mars tensor
We give hereafter a short illustration of SageManifolds focused on tensor calculus in 4-dimensional GR.Another example, to be found at [25], is based on the manifold S 2 and focuses more on the use of multiple charts and on the treatment of non-parallelizable manifolds.Yet another example illustrates some graphical capabilities of SageManifolds: Figure 6 shows the famous immersion of the real projective plane RP 2 into the Euclidean space R 3 known as the Boy surface.This figure has been obtained by means of the method plot() applied to three coordinate charts covering RP 2 , the definition of which is related to the interpretation of RP 2 as the set of straight lines ∆ through the origin of R 3 : (i) in red, the chart X 1 covering the open subset of RP 2 defined by all lines ∆ that are not parallel to the plane z = 0, the coordinates of X 1 being the coordinates (x, y) of the intersection of the considered line ∆ with the plane z = 1; (ii) in green, the chart X 2 covering the open subset defined by all lines ∆ that are not parallel to the plane x = 0, the coordinates of X 2 being the coordinates (y, z) of intersection with the plane x = 1; (iii) in blue, the chart X 3 covering the open subset defined by all lines ∆ that are not parallel to the plane y = 0, the coordinates of X 2 being the coordinates (z, x) of intersection with the plane y = 1. Figure 6 actually shows the coordinate grids of these three charts through the Apéry map [26], which realizes an immersion of RP 2 into R 3 .This example, as many others, can be found at [25].
Let us consider a 4-dimensional spacetime, i.e. a smooth 4-manifold M endowed with a Lorentzian metric g.We assume that (M, g) is stationary and denote by ξ the corresponding Killing vector field.The Simon-Mars tensor w.r.t.ξ is then the type-(0,3) tensor field S defined by [27] S where The Simon-Mars tensor provides a nice characterization of Kerr spacetime, according the following theorem proved by Mars [27]: if g satisfies the vacuum Einstein equation and (M, g) contains a stationary asymptotically flat end M ∞ such that ξ tends to a time translation at infinity in M ∞ and the Komar mass of ξ in M ∞ is non-zero, then S = 0 if, and only if, (M, g) is locally isometric to a Kerr spacetime.In what follows, we use SageManifolds to compute the Simon-Mars tensor according to formula (10) for the Kerr metric and check that we get zero (the "if" part of the above theorem).The corresponding worksheet can be downloaded from http://sagemanifolds.obspm.fr/examples/html/SM_Simon-Mars_Kerr.html.For the sake of clarity, let us recall that, as an object-oriented language, Python (and hence Sage) makes use of the following postfix notation:

function(arguments)
In a functional language, this would correspond to result = function(object,arguments).For instance, the Riemann tensor of a metric g is obtained as riem = g.riemann()(in this case, there is no extra argument, hence the empty parentheses).With this in mind, let us proceed with the computation by means of SageManifolds.In the text below, the blue color denotes the outputs as they appear in the Sage notebook (note that all outputs are automatically L A T E Xformatted by Sage).
The first step is to declare the Kerr spacetime (or more precisely the part of the Kerr spacetime covered by Boyer-Lindquist coordinates) as a 4-dimensional manifold: The standard Boyer-Lindquist coordinates (t, r, θ, φ) are introduced by declaring a chart X on M, via the method chart(), the argument of which is a string expressing the coordinates names, their ranges (the default is (−∞, +∞)) and their L A T E X symbols: X.<t,r,th,ph> = M.chart('t r:(0,+oo) th:(0,pi):\\theta ph:(0,2*pi):\\phi') print X ; X chart (M, (t, r, th, ph)) (M, (t, r, θ, φ)) We define next the Kerr metric g by setting its components in the coordinate frame associated with Boyer-Lindquist coordinates.Since the latter is the current manifold's default frame (being the only one defined at this stage), we do not need to specify it when referring to the components by their indices: The Levi-Civita connection ∇ associated with g is obtained by the method connection(): nab = g.connection() ; print nab Levi-Civita connection 'nabla g' associated with the Lorentzian metric 'g' on the 4-dimensional manifold 'M' As a check, we verify that the covariant derivative of g with respect to ∇ vanishes identically: nab(g).view() As mentionned above, the default vector frame on the spacetime manifold is the coordinate basis associated with Boyer-Lindquist coordinates:

Figure 2 .
Figure 2. Python classes for scalar fields on a manifold.

Figure 3 .
Figure 3. Python classes for modules.For each of them, the class of the base ring is indicated, as well as the class for the elements.

Figure 4 .
Figure 4. Python classes implementing tensors and tensor fields.

Figure 5 .
Figure 5.Storage of tensor fields in SageManifolds.Each red box represents a Python dictionary; the dictionary values are depicted by yellow boxes, the keys being indicated at the left of each box.

Figure 6 .
Figure 6.Boy surface depicted via the grids of 3 coordinate charts covering RP 2 (see the text for the color code).