# Contracting indices in lists

Using an answer to this question (Working with Dirac Gamma Matrices using FeynCalc - A simple problem), I defined

Clear[γ]
SetAttributes[γ, Listable]
γ[μ_] := If[μ == 0, KroneckerProduct @@ PauliMatrix[{1, μ}],
I KroneckerProduct @@ PauliMatrix[{2, μ}]
]

So γ is an object labelled by μ, which is an index whose value runs from 0 to 3 (0, 1, 2, 3). I would like to contract γ with γ, i.e.

,

where the summation over repeated indices is understood. (For physicists: I also defined γ with lower indices, which has the spatial components with a different sign, but this is not the point now).

The question is: how can I do this? I don't think I can use TensorProduct and/or TensorContract, since γ is not a tensor.

And what if I want to contract one or more indices of a γ-like object which is labelled by two indices or more?

P.S.: I do not want to use FeynCalc.

• Well, what a coincidence that the "answer" you mentioned was just given by me years ago :). May 18, 2020 at 5:01

index[_[x__]] := x (* <-- Extracts arguments from an expression. *)

muteIndexSum[list_] := Module[
{indices, repeatedIndices, result},
(*ALL indices: *)
indices = Table[index[gamma], {gamma, list}];
(*REPEATED indices: *)
repeatedIndices = {};
Do[If[Count[indices, i] == 2 && ContainsNone[repeatedIndices, {i}],
AppendTo[repeatedIndices, i]], {i, indices}];
(*Dot product of gamma matrices: *)
result = Apply[Dot, list];
Do[result = Sum[result /. i -> k, {k, 0, 3}], {i,
repeatedIndices}];
(* If there are no repeated indices,
this will just return the dot product of what you passed. *)

Return[result];
]

Since you have already used γ, I am using γm as a notation (using γ would mess things up):

muteIndexSum[{γm[μ], γm[ν], γm[μ]}]

Output:

γm[0].γm[ν].γm[0] + γm[1].γm[ν].γm[1] + γm[2].γm[ν].γm[2] + γm[3].γm[ν].γm[3]

If you contract all indices,

muteIndexSum[{γm[μ], γm[ν], γm[μ], γm[ν]}]

The output will be:

γm[0].γm[0].γm[0].γm[0] + γm[0].γm[1].γm[0].γm[1] + γm[0].γm[2].γm[0].γm[2] + γm[0].γm[3].γm[0].γm[3] + γm[1].γm[0].γm[1].γm[0] + γm[1].γm[1].γm[1].γm[1] + γm[1].γm[2].γm[1].γm[2] + γm[1].γm[3].γm[1].γm[3] + γm[2].γm[0].γm[2].γm[0] + γm[2].γm[1].γm[2].γm[1] + γm[2].γm[2].γm[2].γm[2] + γm[2].γm[3].γm[2].γm[3] + γm[3].γm[0].γm[3].γm[0] + γm[3].γm[1].γm[3].γm[1] + γm[3].γm[2].γm[3].γm[2] + γm[3].γm[3].γm[3].γm[3]

If you want to use your pre-defined γ matrices, just use the rule to substitute:

muteIndexSum[{γm[μ], γm[ν], γm[μ], γm[ν]}]/.{γm -> γ}

Output:

{{4, 0, 0, 0}, {0, 4, 0, 0}, {0, 0, 4, 0}, {0, 0, 0, 4}}
• Thank you, it works perfectly! I have a silly question although: if I define Clear[[Gamma]alti] SetAttributes[[Gamma]alti, Listable] [Gamma]alti[[Mu]_] := If[[Mu] == 0, KroneckerProduct @@ PauliMatrix[{1, [Mu]}], I KroneckerProduct @@ PauliMatrix[{2, [Mu]}]] and use for example muteIndexSum[{γm[μ], γm[ν], γm[μ], γm[ν]}]/.{γm -> γalti}, it does not work... How can the name be a problem?
– FB20
May 17, 2020 at 15:02
• Renaming should not cause problems. But sometimes Mathematica misbehaves. Try saving the .nb file, restarting Mathematica and evaluating the notebook. If changing the name is all you did, it should work just fine. May 17, 2020 at 15:10
• P.s.: If this solved your problem, you can mark it as the correct answer. May 17, 2020 at 15:31
• You were right, I restarted and it's all fine now. But I just noticed that your code doesn't work if I want to contract two different objects, like γalti and γbassi. I became aware of that since γ_mu γ^mu should be equal to 4I, when I is the 4x4 identity matrix. Any solution for that? In physics, when the metric is not euclidean, tensors with upper indices are not equal to tensors with lower indices. I'm not able to evaluate if for this other problem a new question is necessary.
– FB20
May 17, 2020 at 17:17
• *where I is the 4x4 identity matrix.
– FB20
May 17, 2020 at 17:28

What the shape of $$\gamma_\mu \gamma^\mu$$ do you expect? I am not sure whether the code below works to your satisfaction or not:

Total @ MapThread[Dot, {γ[Range[0, 3]], γ[Range[0, 3]]}]

Update

OK, let me make it in a more formal and natural way. First let some points made clear:

• Contraction of a pair of two same Greek indexes is accompanied with Minkowski metric signatured $$\mathrm{diag} g_{\mu\nu} = \{+, -, -, -\}$$ (g = DiagonalMatrix[SparseArray @ {1, -1, -1, -1}];);
• Contraction of a pair of two same Latin indexes is the same as that of a matrix dot product, or say the metric is just the identity matrix;
• Einstein convention is employed.

Then the quantity of interest in fact is

$$\gamma^\mu \gamma_\mu = (\gamma^\mu)^{mn} (\gamma_\mu)_{nl} = (\gamma^{^{\substack{\color{red}{1} \\ \mu}}})^{^{\substack{2\; \color{blue}{3} \\ mn}}} g_{_{\substack{\mu\nu \\ \color{red}{4}\,\color{cyan}{5}}}} (\gamma^{^{\substack{\color{cyan}{6} \\ \nu}}})_{_{\substack{nl \\ \color{blue}{7}8}}}$$

So it should work by using TensorProduct with TensorContract (pay attention to which two indexes are paired to contract):

γμ = γ[Range[0, 3]];
TensorContract[TensorProduct[γμ, g, γμ], {{1, 4}, {5, 6}, {3, 7}}]

In the same spirit $$\sigma^{\mu\nu}\sigma_{\nu\rho} = (\sigma^{\mu\nu})^{mn}\ g_{\nu\alpha}\ (\sigma^{\alpha\beta})_{nl}\ g_{\beta\rho} \\ \sigma^{\mu\nu}\sigma_{\mu\nu} = (\sigma^{\mu\nu})^{mn}\ g_{\mu\alpha}\ (\sigma^{\alpha\beta})_{nl}\ g_{\beta\nu}$$ should correspond to, respectively,

Clear[σ]
σ[μ_, ν_] := I (γ[μ].γ[ν] - γ[ν].γ[μ]) / 4
σμν = Outer[σ, Range[0, 3], Range[0, 3]] // SparseArray;
TensorContract[TensorProduct[σμν, g, σμν, g], {{2, 5}, {4, 9}, {6, 7}, {8, 11}}]
TensorContract[TensorProduct[σμν, g, σμν, g], {{1, 5}, {2, 12}, {4, 9}, {6, 7}, {8, 11}}]

Other contractions can be done similarly.

• Thank you for your answer! It works, but I have problem with two-indices objects. In my code, I also use: Clear[[Sigma]] SetAttributes[[Sigma], Listable] [Sigma][[Mu]_, [Nu]_] := I ([Gamma][[Mu]].[Gamma][[Nu]] - [Gamma][[Nu]].\ [Gamma][[Mu]])/4 Now what if I want to contract two indices of four, or every indices? It would be fantastic to use such a short and quick command like yours.
– FB20
May 17, 2020 at 14:55
• @FB20 Plz see my update. May 18, 2020 at 4:44
• Thanks a lot! I'm sorry to bother you again, but I have another doubt; take as an example your TensorContract[TensorProduct[σμν, g, σμν, g], {{2, 5}, {4, 9}, {6, 7}, {8, 11}}]. This is an object with two free (spacetime) indices: mu and rho. I would like to define a two-indices tensor which is precisely this quantity. I tried with the standard definition ":=" but I failed...
– FB20
May 18, 2020 at 20:39
• @FB20 Maybe you can use Part to extract from what is obtained there. May 19, 2020 at 6:30
• Ok I solved with a little trick but I have another problem ahaha. I made another question.
– FB20
May 19, 2020 at 10:22