# 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/ Commented Feb 19, 2021 at 10:38
• @NiharKarve I haven't recieved meaningful answers there, that's why I'm asking it here. Commented Feb 19, 2021 at 10:40
• Crossposted to physics.stackexchange.com/q/615690/2451 Commented Feb 19, 2021 at 15:09
• Dear mathematica mods. Please merge. Commented Feb 19, 2021 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. Commented Feb 19, 2021 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]; Commented Feb 19, 2021 at 12:34
• After fixing that issue, everything turns out perfect. Commented Feb 19, 2021 at 12:37
• @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 Commented Feb 19, 2021 at 13:08
• @Kasmonium It can do p[x]^4 because that gets expanded out into p[x] ** p[x] ** p[x] ** p[x] and it knows how to deal with those. You can't expand out p[x]^n, so you need to tell it how to deal with expressions of that form. This code works purely structurally; that's why it can't deal with symbolic powers unless you give it definitions to do so. In fact, that's true for computer algebra as a whole. Commented Feb 19, 2021 at 13:25

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. Commented Feb 19, 2021 at 18:30
• Is there anyway to extend this so that commutators can interact better with other expressions? For example even a commutator like [a, e ad] where e is eulers number will not be simplified to e. Unless I unprotect numberQ and set numberQ[e]:=True but this should have been default already. It gets harder if there are integrals and conjugates of unknown functions etc. Do you have any tips for approaching this?
– Luca
Commented Nov 3, 2023 at 10:42

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

• Nice that you are using NCGBX to generate the rules. Unfortunately NCReplaceRepeated would not be able to simplify everything in one shot. The reasons is that replace is essentially pattern matching and it will not expand products that occur in an earlier expansion. By the way, this is true even for the standard Mathematica ReplaceRepeated, as in x y z /. {z -> x + y, x^2 -> 1}. Try it out! Commented Dec 11, 2021 at 22:22

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}$ Commented Feb 19, 2021 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? Commented Feb 19, 2021 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. Commented Feb 19, 2021 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. Commented Feb 19, 2021 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.... Commented Feb 19, 2021 at 13:56

This can be performed using the Quantum Mathematica add-on

http://homepage.cem.itesm.mx/lgomez/quantum/

You can implement any commutator (or even anti-commutator) algebra of your choice. See here for an example with the algebra of Pauli matrices. It works flawlessly and I have used it quite extensively while working with algebras or super-algebras with several operators.