Anticommuting Variables via NonCommutativeMultiply

I am not the first person to implement anti-commuting variables: a quick internet search lead me to the post "How can I implement anti-commuting (Grassmann) numbers (variables)?" on https://mathematica.stackexchange.com. Almost all of the code below is copied from the answers to that post, with a little editing and some changes by me. The basic idea is to implement the anticommutation rules for NonCommutativeMultiply, which is an alias for the ** operator.

Arithmetic

First, we allow the redefinition of NonCommutativeMultiply.

Unprotect[NonCommutativeMultiply];

Then, we need to implement flatness (i.e. $\psi (\chi \xi) \zeta = \psi \chi \xi \zeta$) by hand. If we don't do that first, Mathematica will get stuck in a recursive loop with some of the other rules.

(*One-way flatness:*)
ClearAttributes[NonCommutativeMultiply,Flat]
NonCommutativeMultiply[H___,NonCommutativeMultiply[M___],T___]:=NonCommutativeMultiply[H,M,T]

The first rule that we implement is the distributive property for addition and **.

(*Linearity of addition:*)
NonCommutativeMultiply[H___,Plus[A_,B__],T___]:=NonCommutativeMultiply[H,A,T]+NonCommutativeMultiply[H,Plus[B],T]

Next, we implement multiplication with (commuting) scalars.

(*Scalars come out. Scalar variables are numeric or contain only commuting variables*)
$CommutingVariables={};
ScalarQ[f_]:=Or[NumericQ[f],MemberQ[$CommutingVariables,f],AllTrue[Variables[f],MemberQ[$CommutingVariables,#]&]]
NonCommutativeMultiply[H___,Times[c_,A__],T___]:=c NonCommutativeMultiply[H,Times[A],T]/;ScalarQ[c]
NonCommutativeMultiply[H___,c_,T___]:=c NonCommutativeMultiply[H,T]/;ScalarQ[c]

Before performing computations, all symbols that represent commuting variable must be added to the list $CommutingVariables. Then ScalarQ[c] checks if c is a number, a commuting variable, or an expression that contains only commuting variables. If the check returns true, then $\psi (c \chi)$ becomes $c (\psi \chi)$. Some examples:

In:= 2**3
Out= 6 NonCommutativeMultiply[]

In:= 2**3**x
Out= 6 NonCommutativeMultiply[x]

In:= $CommutingVariables={a,b,c};
In:= a**x**(b+y)
Out= a b NonCommutativeMultiply[x] + a x**y

You can also use Block to localize the value of $CommutingVariables.

If a term contains at most one non-commutative variable, we get rid of NonCommutativeMultiply.

(*Simplify NonCommutativeMultiply[\[Psi]] to \[Psi] and NonCommutativeMultiply[Subscript[\[Psi], i]] to Subscript[\[Psi], i]:*)
NonCommutativeMultiply[H_]:=H/;(Head[H]===Symbol||Head[H]===Subscript)

(*Empty product is 1:*)
NonCommutativeMultiply[]:=1

Lastly, we implement the anti-commutation relations by ordering the terms in a product alphabetically (using OrderedQ) and adding a sign every time we switch two neighbouring symbols. So $\beta \alpha$ becomes $-\alpha \beta$. If an expression contains the same symbol twice, we use a separate rule for $\alpha^2 = 0$.

(*Canonical ordering:*)
NonCommutativeMultiply[H___,B_,A_,T___]:=-NonCommutativeMultiply[H,A,B,T]/;Not[OrderedQ[{B,A}]]

(*Squares vanish:*)
NonCommutativeMultiply[H___,A_,M___,A_,T___]:=0

Note that the "One-way flatness" rule MUST come before "Canonical ordering" and "Squares vanish", otherwise Mathematica would do the following (wrong) simplification: $$ (\beta \gamma) \alpha = -\alpha (\beta \gamma) = -\alpha \beta \gamma $$ when the correct result is $$ (\beta \gamma) \alpha = \beta \gamma \alpha = -\beta \alpha \gamma = \alpha \beta \gamma. $$ This example also illustrates that products of two (or any even number of) anti-commuting variables behaves like a commuting variable (but note that $(\alpha \beta)^2 = 0$).

With $CommutingVariables={}, the following tests should all evaluate to True:

u**(v**x)**y===u**v**x**y
u**(v+x)**y===u**v**y+u**x**y
v**u===-u**v
u**u===0
u**v**u==0
(u-v)**(u+v)===2u**v
(u-v)**(u+v)**(u+v)===0
(u+v)**(v+x)**(x+u)===2u**v**x
(u**v+x**y)**(u**v+x**y)===2 u**v**x**y
2**u===2u
1**2**3===6

I used === to compare the full form of the expressions. They look, for example, like this:

In:= FullForm[2 u**v**x**y]
Out//FullForm= Times[2,NonCommutativeMultiply[u,v,x,y]]

Anticommuting derivatives/Grassmann integrals

The derivatives are linear operators that satisfy $$ \partial_{\psi} c = c \partial_{\psi} 1 = 0, \quad \partial_{\psi} \psi = 1 $$ where the term $c$ does not contain $\psi$. To implement the rules, I define a new function DerivAntiCom that takes a symbol (the variable for $\partial$) and an expression. The following code is taken from this answer.

(*Linearity of addition:*)
DerivAntiCom[\[Chi]_Symbol,Plus[A_,B__]]:=DerivAntiCom[\[Chi],A]+DerivAntiCom[\[Chi],Plus[B]]

(*Constants come out:*)
DerivAntiCom[\[Chi]_Symbol,Times[c_,A__]]:=c DerivAntiCom[\[Chi],A]/;(FreeQ[c,\[Chi]]&&FreeQ[c,NonCommutativeMultiply])

(*Derivative of a constant is zero:*)
DerivAntiCom[\[Chi]_Symbol,c_]:=0/;FreeQ[c,\[Chi]]

(*Derivative of var with respect to var is one:*)
DerivAntiCom[\[Chi]_,\[Chi]_]:=1

Note that we also want $\partial_{\chi} \psi \chi = -\partial_{\chi} \chi \psi = -\psi$. Additionally, for products of derivatives I choose the convention that they are evaluated right-to-left, for example: $$ \partial_{\chi} \partial_{\psi} \chi \psi = -\partial_{\chi} (\partial_{\psi} \psi) \chi = -1 \partial_{\chi} \chi = -1. $$ This is implemented by

(*Derivative of A**var**B with respect to var is (-1)^k A**B, where k is the number of anticommuting variables in A:*)
DerivAntiCom[\[Chi]_Symbol,expr_NonCommutativeMultiply]:=
  Block[{NonCommutativeMultiply},
    Replace[expr,NonCommutativeMultiply[A___,\[Chi],B___]:>(-1)^Length[{A}]NonCommutativeMultiply[A,B]]
  ]

(*Derivative with respect to a list of variables is done starting with the last variable in the list.*)
DerivAntiCom[\[Chi]s_List,expr_]:=If[Length[\[Chi]s]==1,
  DerivAntiCom[\[Chi]s[[1]],expr],
  DerivAntiCom[Most[\[Chi]s],DerivAntiCom[Last[\[Chi]s],expr]]
]

With $CommutingVariables={a,b}, these tests should all evaluate to True:

DerivAntiCom[x,2x+3y+5z]===2
DerivAntiCom[y,2x+3y+5z]===3
DerivAntiCom[z,2x+3y+5z]===5
DerivAntiCom[y,x**y]===-x
DerivAntiCom[{x,y},x**y]===-DerivAntiCom[{y,x},x**y]
DerivAntiCom[x,(x-y)**(a+b x)]===a+b y
DerivAntiCom[y,(x-y)**(a+b x)]===-a-b x
DerivAntiCom[{x,y},(x-y)**(a+b  x)]===-b
DerivAntiCom[{y,x},(x-y)**(a+b  x)]===b

Some convenience functions

The non-commutative analogs of x^p (where p is an integer) and Series[expr, {x, x0, n}] are

PowerNCM[expr_,p_Integer]:=If[p==0,1,Nest[#**expr&,expr,p-1]]

SeriesNCM[expr_,{x_,x0_,n_},\[Chi]_]:=Module[{powers,series},
  powers=NestList[#**\[Chi]&,1,n];
  series=Series[expr,{x,x0,n}];
  Normal[series]/.{x:>\[Chi],x^p_:>powers[[p+1]]}
]

For example (with $CommutingVariables={}):

In:= Table[PowerNCM[a+b+c+a**b,p],{p,0,5}]
Out= {1,a+b+c+a**b,2a**b**c,0,0,0}

In:= SeriesNCM[Exp[x],{x,0,5},a+b+c+a**b]
Out= 1+a+b+c+a**b+a**b**c

Note that SeriesNCM must use a temporary commuting variable x to obtain the Taylor-expansion before x is replaced by an expression containing anti-commuting variables.

The Dot and Det functions cannot be used for vectors and matrices containing anti-commuting variables, because it relies on Times instead of NonCommutativeMultiply. The following functions should cover most cases:

DotNCM[vec1_,vec2_]/;If[
  VectorQ[vec1]&&VectorQ[vec2]&&Length[vec1]==Length[vec2],
  True,
  Message[DotNCM::shape,vec1, vec2];False
]:=Inner[NonCommutativeMultiply,vec1,vec2,Plus]
DotNCM::shape="Vectors `1` and `2` have incompatible shapes.";

MatVecNCM[mat_,vec_]/;If[
  MatrixQ[mat]&&VectorQ[vec]&&Dimensions[mat][[2]]==Length[vec],
  True,
  Message[MatVecNCM::shape,mat, vec];False
]:=Inner[NonCommutativeMultiply,mat,vec,Plus]
MatVecNCM::shape="Matrix `1` and vector `2` have incompatible shapes.";

MatMatNCM[mat1_,mat2_]/;If[
  MatrixQ[mat1]&&MatrixQ[mat2]&&Dimensions[mat1][[2]]==Dimensions[mat2][[1]],
  True,
  Message[MatMatNCM::shape,mat1, mat2];False
]:=Inner[NonCommutativeMultiply,mat1,mat2,Plus]
MatMatNCM::shape="Matrix `1` and `2` have incompatible shapes.";

VecMatVecNCM[vec1_,mat_,vec2_]/;If[
  VectorQ[vec1]&&MatrixQ[mat]&&VectorQ[vec2]&&Dimensions[mat]=={Length[vec1],Length[vec2]},
  True,
  Message[VecMatVecNCM::shape,vec1,mat, vec2];False
]:=Inner[NonCommutativeMultiply,Inner[NonCommutativeMultiply,vec1,mat,Plus],vec2,Plus]
(*Alternative: Module[{i,j},Sum[vec1[[i]]**mat[[i,j]]**vec2[[j]],{i,Dimensions[mat][[1]]},{j,Dimensions[mat][[2]]}]]*)
VecMatVecNCM::shape="Vector `1` and matrix `2` and vector `3` have incompatible shapes.";

DetNCM[mat_]/;If[SquareMatrixQ[mat],True,Message[DetNCM::shape,mat];False]:=
  Module[{i,perm,n=Length[mat]},
  Sum[
    Signature[perm]Apply[NonCommutativeMultiply,Table[mat[[i,perm[[i]]]],{i,n}]],
    {perm,Permutations[Range[n]]}
  ]
]
DetNCM::shape="Matrix `1` is not a square matrix.";

Use

Examples:

In:= DotNCM[{a,d},{c,b}]
Out= a**c-b**d

In:= MatVecNCM[{{a,b},{c,d}},{e,f}]
Out= {a**e+b**f,c**e+d**f}

In:= DetNCM[{{d,c},{b,a}}]
Out= d**a-c**b