How to do Einstein Summation in Mathematica?

Question

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?

Answer 1

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

Second related thread

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[[1]]

τ

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

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

The full term is coded as

lastterm = 
  1/2 1/2 Table[
     D[metric[[m, k]], coord[[1]] ] 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