Simplifying expression involving Boson operators

Question

I want to do calculations involving Boson operators $a(k), a^{\dagger}(k)$, e.g $(x + a(k))(y+a^{\dagger}(k))$. My code is:

NCM[x___] := NonCommutativeMultiply[x];
NCM[] := 1
NCM[___, 0, ___] := 0
NCM[H___, 1, T___] := NCM[H, T]
NCM[left___, a[k1_], SuperDagger[a][k2_], right___] := 
NCM[left, SuperDagger[a][k2], a[k1], right] + 
NCM[left, right]*KroneckerDelta[k1 - k2]
Unprotect[NCM];
NCM[left___, HoldPattern[Times[id, x_]], right___] := 
x  NCM[left, right]
NCM[left___, HoldPattern[Times[a[k_], x_]], right___] := 
x   NCM[left, a[k], right]
NCM[left___, HoldPattern[Times[SuperDagger[a][k_], x_]], right___] := 
x NCM[left, SuperDagger[a][k], right]
Protect[NCM];

where based on Simplifying expression with non-commutating entries I introduced the unit id and based on Boson Commutation relations's first reply I constructed the $[a(k),a^{\dagger}(k')]=\delta_{kk'}$ relation.

Now if I try:

Distribute[(x id + a[k]) ** (y id - SuperDagger[a][k])]
(* I get *)
(id x) ** (id y) + (id x) ** (-SuperDagger[a][k]) + a[k] ** (id y) + 
a[k] ** (-SuperDagger[a][k])

and not x y -x SuperDagger[a][k] + y a[k] - a[k] **SuperDagger[a][k]

I can't figure out why...