# Boson commutation relations

I have tried to adapt this answer to my problem of calculating some bosonic commutation relations, but there are still some issues.

The way I'm implementing the commutator is straightforward:

commutator[x_, y_] :=x**y-y**x


Example: if I want to compute $[b^\dagger b,ab^\dagger]$ I write

commutator[BosonC["b"]**BosonA["b"], BosonA["a"]**BosonC["b"]]


and the output is $ab^\dagger$ as it should be. However this fails when I compute $[a^\dagger a,ab^\dagger]$ (which should be $-ab^\dagger)$:

commutator[BosonC["a"]**BosonA["a"], BosonA["a"]**BosonC["b"]]
Out: a†** a^2 ** b†- a ** b†** a†** a


How can I modify the code in this answer to have it work properly?

EDIT Building on the answers of @QuantumDot and @evanb, I came up with this solution. First I implement the commutator, with Distribute.

NCM[x___] = NonCommutativeMultiply[x];
SetAttributes[commutator, HoldAll]

NCM[] := 1
commutator[NCM[x___], NCM[y___]] := Distribute[NCM[x, y] - NCM[y, x]]
commutator[x_, y_] := Distribute[NCM[x, y] - NCM[y, x]]


Then I implement two tools, one for teaching Mathematica how to swap creation and annihilation operators and one is for operator ordering:

dag[x_] := ToExpression[ToString[x] ~~ "†"]

mode[x_] := With[{x†= dag[x]},
NCM[left___, x, x†, right___] := NCM[left, x†, x, right] + NCM[left, right]]

leftright[L_, R_] := With[{R† = dag[R], L† = dag[L]},
NCM[left___, pr : R | R†, pl : L | L†, right___] := NCM[left, pl, pr, right]]


Now I can use it like this: after evaluating the definitions I input (for instance)

mode[a]
mode[b]
leftright[a,b]


And finally I can evaluate commutators, for instance

commutator[NCM[a†,a] + NCM[b†,b], NCM[a,b†]]
(* 0 *)

• Can you write out the reasoning behind $[a^{\dagger} a, a b^{\dagger}] = -a b^{\dagger}$? It's not clear what rules you apply to get that, from this question or the linked one. Sep 27, 2015 at 14:46
• The dagger is \[Dagger] or <escape> dg <escape> Sep 27, 2015 at 16:38
• Can you post the full code that actually produces a ** b\[Dagger] from commutator[BosonC["b"]**BosonA["b"], BosonA["a"]**BosonC["b"]]? Sep 28, 2015 at 11:32
• @PatrickStevens. See the relations here: bosons (not fermions). Sep 28, 2015 at 17:39
• perhaps you would like to have a look at SNEG package :) nrgljubljana.ijs.si/sneg Sep 28, 2015 at 17:55

The function NonCommutativeMultiply has too long of a name, so I make a short-hand version of it (NCP stands for non-commutative product):

NCP[x___] := NonCommutativeMultiply[x];


Now, here's the code

NCP[] := 1
NCP[left___, a, a†, right___] := NCP[left, a†, a, right] + NCP[left, right]
NCP[left___, b, b†, right___] := NCP[left, b†, b, right] + NCP[left, right]
NCP[left___, pl : a | a†, pr : b | b†, right___] := NCP[left, pr, pl, right]


SetAttributes[commutator, HoldAll]
commutator[NCP[x___], NCP[y___]] := NCP[x, y] - NCP[y, x]


Let's give it a try:

commutator[NCP[b†, b], NCP[a, b†]]


b† ** a

• Great answer too, I'm testing them out. Here I would need linearity to make it perfect. Maybe I can fiddle around with Distribute Sep 30, 2015 at 15:32

You get get far by just implementing the linearity rules for the commutator (and that 0 and 1 "escape" the matrix product. The first two aren't strictly needed, but I include them for completeness:

commutator[Plus[a_,A__],B_]:=commutator[a,B]+commutator[Plus[A],B]
commutator[A_,Plus[b_,B__]]:=commutator[A,b]+commutator[A,Plus[B]]
commutator[A_**B_,C_]:=A**commutator[B,C]+commutator[A,C]**B
commutator[A_,B_**C_]:=B**commutator[A,C]+commutator[A,B]**C
commutator[A_,A_]:=0
commutator[___,0,___]:=0

Unprotect[NonCommutativeMultiply];
NonCommutativeMultiply[___,0,___]:=0
NonCommutativeMultiply[H___, 1, T___] := NonCommutativeMultiply[H, T]


Now evaluating commutator[adag**a,a**bdag] gives

(a**commutator[adag,bdag]+commutator[adag,a]**bdag)**a+adag**a**commutator[a,bdag]


You can get further by implementing the fact that different oscillators commute. For example:

commutator[adag, bdag] = 0
commutator[a, bdag] = 0


Reevaluating commutator[adag**a,a**bdag] gives

commutator[adag, a] ** bdag ** a


Finally, you can implement the canonical relation:

commutator[adag, a] = 1;


To get what you want.

You can of course decorate / generalize what's above so that you don't have to write every relation out by hand.

• This is great! How can I tell MMA that the commutator is antisymmetric (so I can avoid writing the rules for, say, both [a,a†] and [a†,a], etc..?) Sep 30, 2015 at 15:26
• Ah, that's a little bit tricky. You have to make sure you don't create a pattern that always binds, giving an infinite loop. One simple solution is commutator[A_, B_] := -commutator[B, A] /; Not[OrderedQ[{A, B}]]. You can come up with more complicated predicates if you want particular orderings, of course. The above works, for example, if you specify [a,adag]=-1, but it won't deduce that relation from [adag, a]=1. Oct 5, 2015 at 3:49
• Another approach is to use Simplify with a custom TransformationFunction. Assuming the definitions in my answer: tryOtherOrder[e_] := (e /. {commutator[A_, B_] :> -commutator[B, A]}) commutator[a, adag] (* doesn't change *) Simplify[commutator[a, adag], TransformationFunctions -> tryOtherOrder] (* gives -1*) Oct 5, 2015 at 4:00