Boson commutation relations

Question

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 *)

Answer 1

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]

Now your function:

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

Answer 2

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.