evaluation - Simplifying expression with non-commutating entries

Say I have an expression such as

(a + X)**(b-Y)

where the variables above possess the following commutation relations [a, b] = [a, X] = [a, Y] = [b, X] = [b, Y] = 0 and [X, Y] is nonzero.

I want to use mathematica to expand out these terms. This can be done using the "Distribute" command. The input is

Distribute[(a + X)**(b-Y)]

and the output is

a ** b + a ** (-Y) + X ** b + X ** (-Y).

I now want to simplify this expression using the commutation relations given above. I'm not sure how to go about uploading these commutation relations into mathematica such that this expression can be simplified.

Answer

One way to achieve what is asked in the question, is to introduce an identity id which commutes with all quantities but can be used together with X and Y to efficiently tell NonCommutativeMultiply how to treat scalars a, b, and c:

Clear[id]
id /: NonCommutativeMultiply[id, y_] := y
id /: NonCommutativeMultiply[x_, id] := x

Unprotect[NonCommutativeMultiply];
NonCommutativeMultiply[x___, HoldPattern[Times[id, a_]], 
  y___] := a NonCommutativeMultiply[x, id, y]
NonCommutativeMultiply[x___, HoldPattern[Times[X, a_]], 
  y___] := a NonCommutativeMultiply[x, X, y]

NonCommutativeMultiply[x___, HoldPattern[Times[Y, a_]], 
  y___] := a NonCommutativeMultiply[x, Y, y]
Protect[NonCommutativeMultiply];

In defining the properties of NonCommutativeMultiply, I had to temporarily use UnProtect. The purpose of id is seen in the three last lines: I can now define the linearity of ** under regular multiplication by a scalar a in the same way for the elements id, X and Y. But by specifying that id commutes with X and Y, I designate it as the element that accompanies all scalars a: instead of writing a for a commuting element, you therefore have to write a id.

The benefit of this additional convention is that I don't need to add specific definitions for each individual variable name that is intended to be a scalar. Anything that appears in the form a id or id c etc. is a scalar by definition.

Here is what the example looks like:

Distribute[(a id + X) ** (b id - Y)]

(* ==> a b id + b X - a Y - X ** Y *)

In the result, id again appears in the one term that has only scalars in it.

To make the definitions even shorter, you could replace all three lines between Unprotect and Protect by this:

NonCommutativeMultiply[x___,HoldPattern[Times[p:(id|X|Y),a_]],y___]:=
 a NonCommutativeMultiply[x,p,y]