How to manipulate gauge theory in Mathematica?

I want to know if there is a way of typing into Mathematica an expression like the following,

$$\epsilon^{\mu \nu \lambda} f^{abc} A^a_\mu A^b_\nu A^c_\lambda + g\epsilon^{\mu \nu \lambda} A^a_\mu \partial_\lambda A^a_\nu + \bar{\psi}(\gamma^\mu(\partial_\mu + gA^a_\mu T^a))\psi$$

(..repeated indices are understood to be summed over..)

where $g$ is a number, $A^a_\mu$ can be thought of as matrices with $a,b,c = {1,2,3,..,N}$ for some $N$ and $\mu, \nu, \lambda = \{ 0,1,2\}$. So the $\partial_\mu$ are partial derivatives as $\partial_\mu = \frac{\partial}{\partial x^\mu}$. $f^{abc}$ is a set of numbers depending on the values of a,b and c and it is completely cyclic and anti-symmetric in it. $\epsilon^{\mu \nu \lambda}$ evaluates to $0$ if any two or more of its indices are equal and evaluates to 1 or -1 depending on whether the three distinct entries are in cyclic or anti-cyclic order.

$\psi$ should also be thought of as a matrix $\psi^a_i$ where $i,j = \{0,1,2\}$. $\gamma^\mu$ are a chosen set of $3\times 3$ matrices. Each of $T^a$ is a $N \times N$ matrix. Then the terms involving $\psi$ when expanded out look like, $$\bar{\psi}\gamma^\mu \partial _\mu \psi = (\psi^\dagger)^a_i (\gamma^0\gamma^\mu \partial_\mu )_{ij}\psi_j^a \quad{ \rm and }\quad\bar{\psi}\gamma^\mu A^a_\mu T^a \psi = (\psi^\dagger)^a_i(\gamma^\mu)_{ij}(A^c_\mu T^c)^{ab} \psi^b_j$$

• I would like to be able to input the above expression into Mathematica without having to explicitly specify the numbers $f^{abc}$ and the matrices $T^a$. I would like to be able to manipulate the expression with the matrices $T,A,\psi$ and the numbers $f^{abc}$,$g$ as being variables.

If the above is possible then I would eventually like to do something like write $A^a_\mu = B^a_\mu + C^a_\mu$ and $\psi^a_i = \eta ^a _i + \xi ^a _i$ and expand the expression in terms of B,C,$\eta$ and $\xi$.

• What you displayed is a gauge theory Langrangian. The real question is, what operators (or expectation values in the path integral formulation) you need to work with. In any case, I am not at all sure that the generic tensors would help you much. What I would do is to write a custom package which would encode the properties of gauge fields (group properties, Bianchi identities for structure constants), spinors (grassman variables), etc. You may want to search for packages for grassman variables and symbolic manipulations with SU(N), and also look at how FeynCalc was implemented. Commented Nov 28, 2012 at 8:02
• I did something similar for a much simpler Schwinger model (which is a 1+1 U(1) gauge theory), to simplify some diagrammatic computations, and it worked quite well for me, but of course that is a way simpler task. Commented Nov 28, 2012 at 8:04
• As Leonid mentioned, have a look at feyncalc.org/FeynCalcBook/QuantumField, feyncalc.org/FeynCalcBook/ExpandPartialD, etc. The only thing I never bothered to implement is a Levi-Civita tensor with three instead of four indices. But you can easily do that yourself Commented Nov 28, 2012 at 23:31

Please take a look at the new Symbolic Tensors functionality introduced in version 9. Specifically:

• Interesting preview of version 9 - too bad it isn't available yet...
– Jens
Commented Nov 28, 2012 at 4:38
• yes, but this is in version 9, which is not yet available, at the moment. Besides, it seems that the new tensor capabilities do not use index notation, so most physicists will not find them very practical (IMHO) Commented Nov 28, 2012 at 4:38
• Congrats with your 20k milestone! Commented Nov 28, 2012 at 10:32
• Thanks for the pointer to the new documentation Vitaliy - looks very useful, indeed! Commented Nov 28, 2012 at 13:38
• @Jens it is now ;) Commented Nov 28, 2012 at 16:41

What you need is a tensor package for Mathematica. The most advanced one - freely available -is the xAct suite.

Without a tensor package, you will not be able to input/manipulate indexed objects easily.

I have started writing a function to implement (a very simple version of) Einstein's notation either formally or within a manifold of a given dimension.

This might suit your purpose as a starting point?

Clear[einstein];
einstein::usage =
"converts double indices to sums using einstein's convention";

einstein[p_: n][x_] := x /. {
Subscript[a__] Subscript[b__] Subscript[c__] Subscript[d__] :> (
idx =
Tally[Flatten[{Rest[List[a]], Rest[List[b]], Rest[List[c]],
Rest[List[d]]}]];
idx = Select[idx, #[[2]] == 2 &];
If[idx != {},
idx = idx // Transpose // First;
isum = Map[{#, p} &, idx];
Sum[Subscript[a] Subscript[b] Subscript[c] Subscript[d],
Sequence @@ isum // Release],
Subscript[a] Subscript[b] Subscript[c] Subscript[d]] ),
f_[Subscript[a__], d___] g_[Subscript[b__], e___] Subscript[
c__] :> (
idx =
Tally[Flatten[{Rest[List[a]], Rest[List[b]], Rest[List[c]]}]];
idx = Select[idx, #[[2]] == 2 &];
If[idx != {},
idx = idx // Transpose // First;
isum = Map[{#, p} &, idx];
Sum[f[Subscript[a], d] g[Subscript[b], e] Subscript[c],
Sequence @@ isum // Release],
Subscript[a] Subscript[b] Subscript[c]] ),
f_[Subscript[a__], d___] Subscript[b__] Subscript[c__] :> (
idx =
Tally[Flatten[{Rest[List[a]], Rest[List[b]], Rest[List[c]]}]];
idx = Select[idx, #[[2]] == 2 &];
If[idx != {},
idx = idx // Transpose // First;
isum = Map[{#, p} &, idx];
Sum[f[Subscript[a], d] Subscript[b] Subscript[c],
Sequence @@ isum // Release],
Subscript[a] Subscript[b] Subscript[c]] ),
Subscript[a__] Subscript[b__] Subscript[c__] :> (
idx =
Tally[Flatten[{Rest[List[a]], Rest[List[b]], Rest[List[c]]}]];
idx = Select[idx, #[[2]] == 2 &];
If[idx != {},
idx = idx // Transpose // First;
isum = Map[{#, p} &, idx];
Sum[Subscript[a] Subscript[b] Subscript[c],
Sequence @@ isum // Release],
Subscript[a] Subscript[b] Subscript[c]] ),
f_[ Subscript[a__], c___] Subscript[b__] :> (
idx = Tally[Flatten[{Rest[List[a]], Rest[List[b]]}]];
idx = Select[idx, #[[2]] == 2 &];
If[idx != {},
idx = idx // Transpose // First;
isum = Map[{#, p} &, idx];
Sum[f[Subscript[a], c] Subscript[b],
Sequence @@ isum // Release], Subscript[a] Subscript[b]] ),
Subscript[a__] Subscript[b__] :> (
idx = Tally[Flatten[{Rest[List[a]], Rest[List[b]]}]];
idx = Select[idx, #[[2]] == 2 &];
If[idx != {},
idx = idx // Transpose // First;
isum = Map[{#, p} &, idx];
Sum[Subscript[a] Subscript[b], Sequence @@ isum // Release],
Subscript[a] Subscript[b]] ),
Subscript[a__] :> (
idx = Tally[Flatten[{Rest[List[a]]}]];
idx = Select[idx, #[[2]] == 2 &];
If[idx != {},
idx = idx // Transpose // First;
isum = Map[{#, p} &, idx];
Sum[Subscript[a] , Sequence @@ isum // Release],
Subscript[a] ] )}


So that it behaves the following way

Subscript[a, i, i] Subscript[b, k, k] // einstein[]


Subscript[a, i, j] Subscript[b, j] // einstein[]


Subscript[a, i, j] Subscript[b, j] // einstein[3]


Subscript[ϵ, i, j, k] Subscript[a, k, l] Subscript[b, j, l] //einstein[]


This function might be useful for simple minded calculations?

Note however that

 Subscript[a, i] Subscript[a, i] // einstein[]


currently fails; to do with a HoldForm problem which is probably easy to fix(?)

Update

The package VEST may be of interest too.

• May you shed some light here: mathematica.stackexchange.com/questions/39068 Commented Dec 19, 2013 at 23:49
• @Idear It seems to me it would require an Anti-Einstein rule to re convert sums into Einstein's notation? Commented Dec 20, 2013 at 8:20