# How to do Einstein Summation in Mathematica?

I am working with the de Sitter metric which takes the form $$g=\tau^{-2}\left(-\left(\Lambda / 3-\tau^{2}\right)^{-1} d \tau^{2}+\left(\Lambda / 3-\tau^{2}\right) d t^{2}+g_{\mathbb{S}^{2}}\right)$$ in the co-ordinates $$\{\tau, t, \theta, \phi\}.$$ Define $$F=-Q d\tau \wedge dt.$$ Then I want to define the following function which takes a two tensor (metric) and returns a two tensor as well.

\begin{align} L(h)_{\mu \kappa}&=h^{\nu \lambda}\left(F_{\mu \nu} F_{\kappa \lambda}-\frac{1}{2} g_{\mu \kappa} F_{\rho v} F^{\rho}{ }_{\lambda}\right)+\frac{1}{2} h_{\mu \kappa}|F|_{g}^{2} \\ &= h^{\nu \lambda}\left(F_{\mu \nu} F_{\kappa \lambda}-\frac{1}{2} g_{\mu \kappa} F_{\rho v} F^{\rho}{ }_{\lambda}\right)-\frac{Q^2 \tau^4}{2} h_{\mu \kappa}\\ &= g^{\nu a} g^{\lambda b}h_{ab}\left(F_{\mu \nu} F_{\kappa \lambda}-\frac{1}{2} g_{\mu \kappa} g^{\rho c} F_{\rho v} F_{c \lambda}\right)-\frac{Q^2 \tau^4}{2} h_{\mu \kappa}. \end{align}

Here we identify two forms with a two tensor in the following sense that $$F=-Qd\tau\otimes dt+Qdt \otimes d\tau.$$

My goal is to compute $$L$$ for $$d\tau^2/\tau^2$$. Here is the code that I wrote to compute $$L(\tau^{-2}d\tau^2).$$ Note that $$f=(\Lambda/3-\tau^2).$$

g = 1/\[Tau]^2 {{-1/f, 0, 0,
0}, {0, f, 0, 0}, {0, 0, 1, 0}, {0,
0, 0, Sin[\[Theta]]^2}}

invg = Simplify[Inverse[g]]

h = 1/\[Tau]^2 {{1, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}}

F = {{0, -Q, 0, 0}, {Q, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}};

n = 4

expr = Table[Simplify[Sum[invg[[v, a]] invg[[l, b]] h[[a,b]] (F[[u, v]] F[[k, l]] -(1/2)g[[u, k]] invg[[r,c]] F[[r,v]] F[[c,l]])-(1/2)Q^2 \[Tau]^4 h[[u,k]]
, {v, 1, n}, {l, 1, n}, {r, 1, n}, {a, 1, n}, {b, 1, n}, {c, 1, n}]], {u, 4}, {k, 4}]


When I run the above, I get the following output

{{-(4095/2) Q^2 \[Tau]^2, 0, 0, 0}, {0, 31/2 f^2 Q^2 \[Tau]^2, 0,
0}, {0, 0, -(1/2) f Q^2 \[Tau]^2, 0}, {0, 0,
0, -(1/2) f Q^2 \[Tau]^2 Sin[\[Theta]]^2}}


which is a bit weird with especially large constants. I computed by hand the first entry of the matrix, i.e. $$L(h)_{11} = -\frac{1}{2}Q^2 \tau^2 f^2$$

whereas Mathematica gives me a factor of $$-4095$$. Not sure why this is happening?

• Welcome to Mathematica.SE! I hope you will become a regular contributor. To get started, 1) take the introductory tour now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. Feb 21 at 13:24
• @bbgodfrey thanks for the introduction! Feb 21 at 13:25
• If you search this site, you will find several questions relating to relativity, some referencing specialized packages relevant to your question. Feb 21 at 15:31
• You have to put it as mathematica code. In the meantime have a look at TensoriaCalc Feb 22 at 15:09
• @Sumit thanks. I will make the edit soon. Feb 22 at 16:17

Disclaimer$$_1$$: The following is a full answer to the OP, but before that I provided some useful -in my opinion- guidelines on how the issue should be tackled.

Disclaimer$$_2$$: I think I have understood the OP correctly, however I might be wrong. In such case, please, leave a comment and I will remove my answer.

Disclaimer$$_3$$: As you understand I will try to explain this as generally as possible and do a practical application in the end for illustrational purposes. The title that the author chose motivated me, as I do believe that with a general understanding of some basic steps we can build some tensor related quantities relatively easily.

Disclaimer$$_4$$: From some point of view this could be considered -almost- a dupe of the following two:

First related thread; it is closed

The reasoning behind me suggesting that it is -almost- a dupe is precisely that the philosophy in the approach is the same. I will make an honest effort to make this answer useful to the author of the OP and future users that might want to deal with tensor structures -at least some basic stuff thereof.

Disclaimer$$_5$$: We are not showing the explicit computation of the desired quantities in the OP, but rather of the affine connection; the Christoffel symbols with three indices. The structures and code lines are similar -I think I am not missing anything- and hopefully said computation is enough to demonstrate the general approach. I am also pretty certain that the Christoffel symbols for de Sitter are known explicitly, and the author can check for potential mistakes easily.

Disclaimer$$_6$$: You can find in the linked threads that I did not create the whole code from scratch. I found something online, I tweaked, re-wrote, and goes on and on.

Discussion and understanding of the aspects of the task at hand:

In the OP a four-dimensional de Sitter metric is discussed and we are sticking to that, though the approach is easily generalized to other cases. Since there is no clear mention of the the notation $$\dot{g}$$ I will assume that the dot stands for a derivative w.r.t time; that is $$\tau$$ in the conventions used.

Finally, I will also be using an explicit parameterization of the two-dimensional metric given by the following:

$$\begin{equation} g_{S^2} = d\theta^2 + \sin^2 \theta d\phi^2 \, . \end{equation}$$

So, the first thing that we need to address is how to go about coding rank-2 index in four-dimensions; the metric. That's not very hard to do. Your metric looks exactly like this

$$\begin{pmatrix} g_{\tau \tau} & g_{\tau t} & g_{\tau \theta} & g_{\tau \phi} \\ g_{t \tau} & g_{t t} & g_{t \theta} & g_{t \phi} \\ g_{\theta \tau} & g_{\theta t} & g_{\theta \theta} & g_{\theta \phi} \\ g_{\phi \tau} & g_{\phi t} & g_{\phi \theta} & g_{\phi \phi} \end{pmatrix}$$

Ok, so in order to implement this into Mathematica all we need to do is to consider a list of lists. Not too difficult to write down. It looks like the following:

ds4 = {{gττ, gτt, gτθ,
gτφ}, {gtτ, gtt, gtθ,
gtφ}, {gθτ, gθt, gθθ,
gθφ}, {gφτ, gφτ, gφθ,
gφφ}};


and you can display it if you want to have a look as follows:

ds4 // MatrixForm


Note: keep in mind that MatrixForm is for display purposes ONLY and not to be used in the definitions of the quantities. It will cause issues. As a next step we need to understand how we are able to take derivatives of the metric w.r.t a coordinate. But since we are dealing with lists, Part is our friend. Something like the following should do the trick:

D[metric[[s, j]], coord[[k]]]


where in the above metric is the defined list of lists and coord is the coordinate system of your choice which can be passed into Mathematica as a list of the desired symbols.

Finally, how do we make contractions? That's not too bad either to implement. We are dealing with lists and we want to consider the summation over all allowed values.

As far as I can tell that's all you need from the perspective of coding. Maybe I am wrong, but the expression for $$F$$ you have provided is just a rank-2 tensor in four-dimensions -like the metric- that is antisymmetric -unlike the metric. But we have already discussed how to implement a general tensor of said form without any symmetry properties, so you can just write it as you wish.

As a side comment for the interested reader who is wondering how I worked out that $$F$$ has to be antisymmetric there is a fundamental property in differential geometry that for a p-form $$A^{(p)}$$ and a q-form $$B^{(q)}$$ it holds that:

$$\begin{equation} A^{(p)} \wedge B^{(q)} = (-1)^{p \cdot q} B^{(q)} \wedge A^{(p)} \, , \end{equation}$$ and it is obvious that for $$p=1=q$$ we have antisymmetry.

A practical explicit example:

Step 1: We define our coordinates and the dimensionality of our spacetime.

(*Define a list of the coordinates*)
coord = {τ, t, θ, φ};
(*The dimension n of the spacetime*)
n = Length[coord];


Step 2: We define explicitly the metric with which we want to work. Note that we define is the following:

$$\begin{equation} ds^2 = G_{AB} dx^A dx^B \end{equation}$$

in other words the metric elements we define as metric are the ones with the indices down; $$G_{MN}$$ and you can take the Inverse to raise the indices. We define our metric like so:

metric = 1/τ^2 {{-(1/(Λ/3 - τ^2)), 0, 0,
0}, {0, Λ/3 - τ^2, 0, 0}, {0, 0, 1, 0}, {0,
0, 0, Sin[θ]^2}};


with its inverse being given by:

inversemetric = Simplify[Inverse[metric]];


and at this point we can perform a sanity check:

metric.inversemetric // FullSimplify;
% // MatrixForm Step 3: Code and display the Christoffel symbols. Before doing so, for the reader's convenience we remind their mathematical definition

$$\begin{equation} \Gamma^{\lambda}_{\mu \nu} = \frac{1}{2} g^{\lambda \alpha} \left( \partial_{\mu} g_{\alpha \nu} + \partial_{\nu} g_{\mu \alpha} - \partial_{\alpha} g_{\mu \nu} \right) \end{equation}$$

This can be coded as:

affine = Simplify[
Table[(1/2)*
Sum[(inversemetric[[i, s]])*(D[metric[[s, j]], coord[[k]]] +
D[metric[[s, k]], coord[[j]]] -
D[metric[[j, k]], coord[[s]]]), {s, 1, n}], {i, 1, n}, {j, 1,
n}, {k, 1, n}]];


and the following provides a neat display; note that Γ[a,b,c] stands for $$\Gamma^{a}_{b,c}$$ and we are only displaying the independent components -there is a symmetry in the lower indices. Final comment, as it is customary I did the display such that $$0$$ is the temporal component, which is reflected as a minus 1 in the ToString command. If you don't like having indices $$0,1,2,3$$ and you prefer $$1,2,3,4$$ just remove that -1.

listaffine :=
Table[If[UnsameQ[affine[[i, j, k]],
0], {ToString[Γ[i - 1, j - 1, k - 1]],
affine[[i, j, k]]}], {i, 1, n}, {j, 1, n}, {k, 1, j}]
TableForm[Partition[DeleteCases[Flatten[listaffine], Null], 2],
TableSpacing -> {2, 2}] Edit: Let's do the contraction of the rank-2 antisymmetric tensor. In the OP it reads $$|F|^2_g$$. Not quite sure what the index $$g$$ denotes here, so I am computing what I know to be the square, that is:

$$\begin{equation} F^2 = F^{ab}F_{ab} = g^{ac}g^{bd}F_{ab}F_{cd} \end{equation}$$

and to do that I write in Mathematica

Ftnsr = {{0, -QQ, 0, 0}, {QQ, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}};


and then:

1/2 Sum[inversemetric[[a, c]] inversemetric[[b, d]] Ftnsr[[a,
b]] Ftnsr[[c, d]], {a, 1, n}, {b, 1, n}, {c, 1, n}, {d, 1,
n}] // FullSimplify


which yields

-QQ^2 τ^4

The factor $$\tfrac{1}{2}$$ comes from counting the indices and is not to be confused with the one written in the OP. The full thing should be coded as:

1/2 1/2 Sum[
inversemetric[[a, c]] inversemetric[[b, d]] Ftnsr[[a, b]] Ftnsr[[c,
d]], {a, 1, n}, {b, 1, n}, {c, 1, n}, {d, 1,
n}] // FullSimplify


which returns

-(1/2) QQ^2 τ^4

as it should.

Edit 1: link to the original Christoffel related stuff.

Edit 2: The code I provided for the Christoffel symbols -I mean matching of the indices explicitly- is as follows:

$$\begin{equation} \Gamma^{i}_{kj} = \frac{1}{2}g^{is}\left(\partial_k g_{sj} + \partial_j g_{sk} - \partial_s g_{jk} \right) \end{equation}$$

and the related comment that I would like to add at this point is that non-contracted indices are not INSIDE the Sum command, but ONLY in the Table.

I proceed now to analyze one term of your expression. We show the term below

 $$\begin{equation} \frac{1}{2} F^2 \dot{g}_{mk} = \frac{1}{4} g^{ac}g^{bd}F_{ab}F_{cd} \dot{g}_{mk} \end{equation}$$ 

Useful comments: $$g^{xx}$$ is actually inversemetric[[x, x]], $$F_{xx}$$ is Ftnsr[[x, x]] and $$g_{xx}$$ is metric[[x, x]]. Hence, we need to code the following -note I am providing the code from scratch

Working assumption: I am assuming here that $$\dot{g}_{xx}$$ is just the metric components and this is a slightly confusing notation.

coord = {τ, t, θ, φ};
n = Length[coord];
metric = 1/τ^2 {{-(1/(Λ/3 - τ^2)), 0, 0,
0}, {0, Λ/3 - τ^2, 0, 0}, {0, 0, 1, 0}, {0,
0, 0, Sin[θ]^2}};
inversemetric = Simplify[Inverse[metric]];
Ftnsr = {{0, -QQ, 0, 0}, {QQ, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}};


and then we run

lastterm =
1/2 1/2 Table[
metric[[m, k]] Sum[
inversemetric[[a, c]] inversemetric[[b, d]] Ftnsr[[a,
b]] Ftnsr[[c, d]], {a, 1, n}, {b, 1, n}, {c, 1, n}, {d, 1,
n}], {m, 1, n}, {k, 1, n}] // FullSimplify;


We can check the result

lastterm // TableForm Edit 3: In case that $$\dot{g}_{xx}$$ is the derivative w.r.t $$\tau$$.

The only change here is how to code one small bit. From all the coordinates, $$\tau$$ is the first one, namely

coord[]


τ

In order to take the derivative of the metric w.r.t the above coordinate you need to consider

D[metric[[m, k]], coord[] ]


The full term is coded as

lastterm =
1/2 1/2 Table[
D[metric[[m, k]], coord[] ] Sum[
inversemetric[[a, c]] inversemetric[[b, d]] Ftnsr[[a,
b]] Ftnsr[[c, d]], {a, 1, n}, {b, 1, n}, {c, 1, n}, {d, 1,
n}], {m, 1, n}, {k, 1, n}] // FullSimplify;


in this case and yields

lastterm // TableForm Edit 4: Going for the bounty

Below I am providing he full code that computes what you want. I also think you missed a part of the $$L(h)_{11}$$ expression; see below.

Input is

Quit[]


and now

coord = {τ, t, θ, φ};
dim = Length[coord];
metric = 1/τ^2 {{-(1/(Λ/3 - τ^2)), 0, 0,
0}, {0, Λ/3 - τ^2, 0, 0}, {0, 0, 1, 0}, {0,
0, 0, Sin[θ]^2}};
inversemetric = Simplify[Inverse[metric]];
Ftnsr = {{0, -QQ, 0, 0}, {QQ, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}};
htnsr = 1/τ^2 {{1, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0,
0, 0}};


Coding and checking term by term the expression given for $$L(h)$$. We have

term1 = Table[
Sum[inversemetric[[n, a]] inversemetric[[l, b]] htnsr[[a,
b]] Ftnsr[[m, n]] Ftnsr[[k, l]], {n, 1, dim}, {a, 1, dim}, {l,
1, dim}, {b, 1, dim}], {m, 1, dim}, {k, 1, dim}] //
FullSimplify;
term1 // TableForm term2 = -(1/2) Table[
Sum[inversemetric[[n, a]] inversemetric[[l, b]] metric[[m,
k]] inversemetric[[r, c]] htnsr[[a, b]] Ftnsr[[r,
n]] Ftnsr[[c, l]], {n, 1, dim}, {a, 1, dim}, {l, 1, dim}, {b,
1, dim}, {c, 1, dim}, {r, 1, dim}], {m, 1, dim}, {k, 1,
dim}] // FullSimplify;
term2 // TableForm term3 = -(1/2) (QQ τ^2)/
4 Table[htnsr[[m, k]], {m, 1, dim}, {k, 1, dim}] // FullSimplify;
term3 // TableForm And finally the sum of the above contributions

term1 + term2 + term3 // Expand // TableForm • Thanks for your detailed response, I tried to follow your strategy and wrote the code above, but I get some weird constants in the summation. Could you perhaps suggest why that could happen. Thanks for your time! Feb 28 at 10:25
• @ShreyAryan is this the code that has appeared in the OP, now? If yes, I can try to have a look a bit later, just let me know.
– kcr
Feb 28 at 17:39
• yeah, I just added it, thanks for looking at it! Feb 28 at 17:39
• thanks a lot for taking the time to write this answer. This is very helpful for my work. Mar 4 at 20:24
• yeah the 1/2 factor is wrong. However, for the third term is simply $Q\tau^4/2$ multiplied by the input 2-tensor. So if you make that change then the $Q^2/8$ term will become $Q^2 \tau^2/2$ which makes the $L(h)_{11}=0$ Mar 4 at 20:44