# 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. –  Leonid Shifrin Nov 28 '12 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. –  Leonid Shifrin Nov 28 '12 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 –  Rolf Mertig Nov 28 '12 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 Nov 28 '12 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) –  magma Nov 28 '12 at 4:38
Congrats with your 20k milestone! –  Sjoerd C. de Vries Nov 28 '12 at 10:32
Thanks for the pointer to the new documentation Vitaliy - looks very useful, indeed! –  Mark McClure Nov 28 '12 at 13:38
@Jens it is now ;) –  Vitaliy Kaurov Nov 28 '12 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]


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(?)

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