# Calculating commutators in quantum mechanics symbolically with the help of Mathematica

Say I have a hamiltonian $$H=p_x^2+p_y^2+x^2+y^2+x^4+y^4$$, and I want to calculate the commutator $$\left[A,B\right] \equiv AB-BA$$ of arbitrary operators $$\mathcal{O}=x^ap_x^by^cp_y^d$$, where $$a,b,c,d\in\mathbb{N}$$, subject to 2 simple rules.

1. $$\left[x,y\right]=\left[p_x,p_y\right]=\left[x,p_y\right]=\left[y,p_x\right]=0\$$ and $$[x,p_x]=[y,p_y]=i$$. These are the usual commutation relations of quantum mechanics.
2. Another rule I want to impose is that all the $$p$$'s to be at the right and all the $$x$$'s to be on the left.

I want to define such a commutator operation and simply the expression $$\left[H,\mathcal{O}\right]$$ with the rules I have expressed earlier.

Are there any existing Mathematica packages that can achieve this purpose. If there is no such package, how should I go about creating my own?

• Maybe the following can help: feyncalc.github.io or math.ucsd.edu/~ncalg/DOCUMENTATION/ – Daniel Huber Feb 19 at 10:38
• @NiharKarve I haven't recieved meaningful answers there, that's why I'm asking it here. – jah Feb 19 at 10:40
• Crossposted to physics.stackexchange.com/q/615690/2451 – Qmechanic Feb 19 at 15:09
• Dear mathematica mods. Please merge. – Qmechanic Feb 19 at 15:09

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]

• That looks amazing. I will try it with other cases. Thank you for your work. – jah Feb 19 at 12:13
• So I found a problem in your code, a minor mistake. q[i_] ** p[i_] := I - p[i] ** q[i]; This line of code is not true, because it implements the anticommutator not $[x,p_x]=i$. It should be q[i_] ** p[i_] := I + p[i] ** q[i]; – jah Feb 19 at 12:34
• After fixing that issue, everything turns out perfect. – jah Feb 19 at 12:37
• But how can I also evaluate something like comm[q[x], p[x]^n], which should give i n p[x]^(n-1)? How can I work with integer exponents without entering specific numbers? – jah Feb 19 at 12:47
• @Kasmonium Ah, yes. Dumb typo there. Good thing you caught it; I just edited the answer. As for the other question: I guess you can just define that as a reduction rule? Like q[i_] ** p[i_]^n_ := I n p[i]^(n - 1) + p[i]^n ** q[i] (be sure to check my maths)? If you implement that, you might be able to leave out the last two definitions (that start with p /: and q /:) as well – Sjoerd Smit Feb 19 at 13:08

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]


• I knew I must have been re-inventing the wheel. Well, at least it was educational for my self. – Sjoerd Smit Feb 19 at 18:30

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


I don't know about a package but generally, I would consider it very easy to implement.

In order to compute the commutator relations mathematically, you would need to know how they operate on a wave function. Hence to compute a commutator relation for two operators A,B, you would calculate [A,B]psi.

comm[a_, b_, f_] := Simplify[a[b[f]] - b[a[f]]]


Here is an example to get the commutation relation for the position and momentum operator (I've set hbar to 1):

comm[x*# &, -I*D[#, x] &, f[x]]


will give you

I f[x]
`
• You misunderstood me. I don't care about how the operator acts on a wave function. I am just interested in simplifying commutators, so your introduction of a function f is totally unnecesary. Also how did you implement the canonical commutation relation in your code $[x,p_x]=i$. Thinking of the operators as differential operators acting on functions is cumbersome. I need a code which treats the operators $x,y,p_x,p_y$ as non commuting variables. Not differential operators e.g. $-i\frac{\partial}{\partial x}$ – jah Feb 19 at 10:34
• Going through your answer again, it makes more sense, because when you implement momentum as a differential operator, when you simplify the expression it automatically goes to right. But the problem is after computing the commutator, I want to see the partial derivatives as abstract momentum operators $p_x,p_y$. How can turn these derivatives into abstract operators after computing the commutator? – jah Feb 19 at 10:50
• I wouldn't say that introducing a function f is unnecessary. Let me ask this question: how do we know that [x,px]=i? It can only be established if we represent the operators in either position or momentum space and see how they act on a wave function. One has to be careful about operators as we need to always check how they act on other objects in the equation. The definition that I provide should allow you to calculate all different commutator relations. – barnold Feb 19 at 13:38
• No that is where you are wrong. You do not need momentum space or position space, e.g. the wave functions $\psi (x) =\left < x | \Psi \right >, \phi (p) = \left < p | \Psi \right >$. The momentum space and position space wave functions are just components of the abstract state ket that lives in our Hilbert space $\mathcal{H}$. Operators $x , p_x$ act on $\left | \Psi \right > \in \mathcal{H}$. You are right in saying that these operators are meaningless without anything to act upon, but you are talking about components of the state ket. You don't need wave functions in specific basis. – jah Feb 19 at 13:49
• ... (continuing) You don't need a wave function in a specific basis in order to talk about the algebra of these operators. Coming back to your question "how do we know that $[x,p_x]= i$?". You might think that the canonical commutation relation is derived from the identification $p \rightarrow -i\frac{\partial}{\partial x}$ and $[x, -i\frac{\partial}{\partial x}]f(x) = if(x)$ but that is not the case. That formula is just how the canonical commutation relation looks in the position basis. To understand how the commutation relation is derived.... – jah Feb 19 at 13:56