Linear operator algebra: How to distinguish scalars and operators?

Question

In the book "A Physicist' s Guide to Mathematica" Patrick Tam introduces symbolic rules for an algebra of linear operators. Mathematically the rules are as follows, where upper (lower) case letters denote operators (scalars).

identity operator U
$AU = U A = A$

Distributive property
$A (B + C) = AB + AC$
$(A + B) C = AC + BC$

Associative property
$A(BC) = (AB)C$

"Scalar property"
$A(aB) = a AB$
$(aA)B = a AB$

He implements these rules based on the built-in function NonCommutativeMultiply as it is already associative and in general noncommutative as required. In Mathematica this might be written as:

Unprotect[NonCommutativeMultiply];

A_ ** U := A;
U ** A_ := A;

A_ ** (B_ + C_) := A ** B + A ** C;
(A_ + B_) ** C_ := A ** C + B ** C;

A_ ** ((x_. y_^n_. /; number3Q[x, y, n]) B_) := ((x y^n) A ** B);
((x_. y_^n_. /; number3Q[x, y, n]) A_) ** B_ := ((x y^n) A ** B);

Protect[NonCommutativeMultiply];

number3Q[x_, y_, n_] := NumberQ[x] && NumberQ[y] && NumberQ[n];

To tell Mathematica that a is a scalar we can write NumberQ[a]^=True. Thanks to the pattern x_. y_^n_. it will recognise some composite expressions involving a as scalars, but for example not Sqrt[a+1] or Sin[a]. Hence, I have the following questions.

1. Is there a more robust way to identify scalar expressions, that only include combinations of mathematical functions and scalars? A scalar is a symbol with NumberQ[symbol]==True.

2. Is there a manual way to tell that a full expression is a "number", that it will give NumberQ[expr]==True?

3. Is it more advisable to label operators with a special head and treat everything else as scalars, instead of trying to identify scalars? How would an implementation look like?

Answer 1

Using NumericQ instead of NumberQ:

Unprotect[NonCommutativeMultiply];

A_ ** U := A;
U ** A_ := A;

A_ ** (B_ + C_) := A ** B + A ** C;
(A_ + B_) ** C_ := A ** C + B ** C;

A_ ** (a_ B_) /; NumericQ[a] := a (A ** B)
(a_ A_) ** B_ /; NumericQ[a] := a (A ** B)
A_ ** (B_ a_) /; NumericQ[a] := a (A ** B)
(A_ a_) ** B_ /; NumericQ[a] := a (A ** B)

NumericQ[a] ^= True

Protect[NonCommutativeMultiply];

Now A ** (B Sin[a]) gives,

(* A ** B Sin[a] *)

which if you look at the FullForm is actually Times[NonCommutativeMultiply[A,B],Sin[a]], i.e. the scalar rule has been applied.