Calculating commutators in quantum mechanics symbolically with the help of Mathematica
I think I've figured it out, but I haven't checked it super carefully so buyer beware:
Unprotect[NonCommutativeMultiply];
ClearAll[NonCommutativeMultiply];
NonCommutativeMultiply[] := 1;
NonCommutativeMultiply[a_] := a;
NonCommutativeMultiply[first___, const_?NumericQ*b_, rest___] :=
const*NonCommutativeMultiply[first, b, rest];
NonCommutativeMultiply[first___, const_?NumericQ, rest___] :=
const*NonCommutativeMultiply[first, rest];
MakeBoxes[
NonCommutativeMultiply[first___,
args : Longest@Repeated[x_, {2, \[Infinity]}], rest___], form_] :=
RowBox[Flatten@{
If[Length[{first}] > 0,
{MakeBoxes[NonCommutativeMultiply[first], form], "**"},
Nothing
],
SuperscriptBox[MakeBoxes[x, form], ToBoxes[Length[{args}], form]],
If[Length[{rest}] > 0,
{"**", MakeBoxes[NonCommutativeMultiply[rest], form]},
Nothing
]
}];
MakeBoxes[NonCommutativeMultiply[arg_], form_] := MakeBoxes[arg, form]
SetAttributes[NonCommutativeMultiply, Flat];
q[i_] ** p[i_] := I + p[i] ** q[i];
q[i_] ** p[j_] := p[j] ** q[i];
q[i_] ** q[j_] /; ! OrderedQ[{i, j}] := q[j] ** q[i]
p[i_] ** p[j_] /; ! OrderedQ[{i, j}] := p[j] ** p[i]
a_ ** (b_ + c_) := a ** b + a ** c;
(b_ + c_) ** a_ := b ** a + c ** a;
(* Allowing for powers in input and output *)
p /: p[i_]^n_Integer :=
NonCommutativeMultiply @@ ConstantArray[p[i], n];
q /: q[i_]^n_Integer :=
NonCommutativeMultiply @@ ConstantArray[q[i], n];
As you can see, I went with p
and q
to denote the conjugate variables to distinguish them from xyz as the directions. You will need to use **
to multiply them together, but you can use powers as well:
comm[a_, b_] := a ** b - b ** a
comm[q[x], p[x]]
comm[q[x], p[y]]
comm[p[x] ** q[x], p[x]]
comm[q[x] ** p[x]^3, p[x]]
I - 2 p[x] ** q[x]
0
-2 p[x]^2 ** q[x] + I p[x]
-I p[x]^3 + 2 p[x]^4 ** q[x]
Edit
You can add the following definition if you want to be able to raise to symbolic powers:
q[i_] ** p[i_]^n_ := I n p[i]^(n - 1) + p[i]^n ** q[i]
James F. Feagin's Quantum Methods with Mathematica book has an elegant implementation of this in chapter 15.1 Commutator Algebra.
It's along the lines of @Sjoerd's answer (but figured I'd provide the reference to the book above), first defining typical identities for the NonCommutativeMultiply
symbol:
Unprotect[NonCommutativeMultiply];
A_ ** (B_ + C_) := A ** B + A ** C
(B_ + C_) ** A_ := B ** A + C ** A
A_ ** c_?NumberQ := c A
c_?NumberQ ** A_ := c A
A_ ** (B_ c_?NumberQ) := c A ** B
(A_ c_?NumberQ) ** B_ := c A ** B
A_ ** (B_ c_Rational) := c A ** B
(A_ c_Rational) ** B_ := c A ** B
A_ ** (B_ c_Power) := c A ** B
(A_ c_Power) ** B_ := c A ** B
and then defining a fundamental commutation expression. e.g. for OP's case:
commutator[A_, B_] := A ** B - B ** A
fundamentalCommutation[expr_] := ExpandAll[expr //. p[i_] ** q[i_] :> q[i] ** p[i] - I h]
which indeed recovers the derivative action of the momentum operator:
h /: NumberQ[h] = True;
{commutator[p[x]/(-I h), q[x]],
commutator[p[x]/(-I h), q[x] ** q[x]],
commutator[p[x]/(-I h), q[x] ** q[x] ** q[x]]} //fundamentalCommutation
{1, 2 q[x], 3 q[x] ** q[x]}
It is then easy to use a different fundamental commutation expression, e.g. for working with raising and lowering operators:
fundamentalComm[expr_] := ExpandAll[expr //. a ** ad :> ad ** a + 1]
{commutator[a, ad], commutator[ad, a], commutator[a ** ad, ad]} // fundamentalComm
{1, -1, ad}
That can be done with the NCAlgebra (non-commutative algebra) package, see the documentation.
Example:
(* Import package *)
<< NC`
<< NCAlgebra`
<< NCGBX`
SetNonCommutative[x, y, px, py]
SetMonomialOrder[x, y, px, py] (* x to the left, p to the right *)
NCSetOutput[NonCommutativeMultiply -> True] (* pretty output *)
(* commutation relations *)
gb = NCMakeGB[{
x ** y - y ** x,
px ** py - py ** px,
x ** py - py ** x,
y ** px - px ** y,
x ** px - px ** x - I,
y ** py - py ** y - I}, 20];
(* Define expression that should be simplified *)
Comm[A_, B_] := A ** B - B ** A;
H = px ** px + py ** py + x ** x + y ** y + x ** x ** x ** x + y ** y ** y ** y;
expression = Comm[H, px] // NCExpand;
NCReplaceRepeated[expression, gb] // NCExpand
--> I x - x ** px ** x + x ** x ** px + I x ** x ** x - x ** px ** x ** x ** x + x ** x ** x ** x ** px
According to the docs, I think the above should have worked, but apparently it doesn't. It seems NCReplaceRepeated and NCExpand have to be applied multiple times until the result "converges":
foo[X_] := NCExpand@NCReplaceRepeated[X, gb]
expression // foo // foo // foo // foo
--> 2 I x + 4 I x ** x ** x