Return to Answer

Final note to answer Question 3: This functionality cannot be used for an infinite number of generators and/or relations, because it relies on the explicit computation of a Groebner basis, and then the explicit computation of the reduction of the polynomial relative to this basis. This can only be done, of course, if the numbers of generators and relations are both finite.

One last variation. Let's show how to make a three-oscillator algebra. The generators are

gens = Join[a /@ Range[3], aD /@ Range[3]]
(* {a[1], a[2], a[3], aD[1], aD[2], aD[3]} *)

We need the relations $\hat{a}_n\hat{a}_m^{\dagger} - \hat{a}_m^{\dagger}\hat{a}_n =\delta_{nm}$ and $\hat{a}_n\hat{a}_m - \hat{a}_m\hat{a}_n =0$, and we implement them as

rels = Join[
         Table[{a[n] ** aD[m] - aD[m] ** a[n] - KroneckerDelta[m, n]}, {n, 1, 3}, {m, 1, 3}],
         Table[{a[n] ** a[m] - a[m] ** a[n]}, {n, 1, 3}, {m, n + 1, 3}]
        ] // Flatten;

We can generate a Groebner basis from this as

gbrel = NonCommutativeGroebnerBasis[rels, gens];

but in this case, the relations form a basis already. Then, finally, define

toNormalOrder[expr_] :=  NonCommutativePolynomialReduce[expr, rels, gens(*, alg*)]

as usual (with alg in there in case you need to add some scalars, as above).

Final note to answer Question 3: This functionality cannot be used for an infinite number of generators and/or relations, because it relies on the explicit computation of a Groebner basis, and then the explicit computation of the reduction of the polynomial relative to this basis. This can only be done, of course, if the numbers of generators and relations are both finite.

I'm working on implementing a function that will compute general matrix elements in the oscillator eigenstates, but that can be fiddly, so perhaps I'll add it later.

To compute arbitrary matrix elements in the oscillator eigenbasis, we exploit the definitions of the eigenstates in terms of the raising operator to write $$ \langle m \lvert f(\hat{a},\hat{a}^{\dagger}) \rvert n \rangle = \langle 0 \lvert \left( \frac{1}{\sqrt{m!}}\hat{a}^m\right) f(\hat{a},\hat{a}^{\dagger}) \lvert \left( \frac{1}{\sqrt{n!}}\left(\hat{a}^{\dagger}\right)^n\right) \rvert 0 \rangle\,, $$ and so take the vacuum expectation of the modified operator that include the extra powers of the raising and lowering operators.

To do this, we have to sort-of kluge something together, and this is mainly because the default for GeneralizedPower is to leave

GeneralizedPower[NonCommutativeMultiply, a, 0]

unevaluated. I assume this is because we can't assume that the identity is 1 in the algebra we're working with. Here is the function:

Clear@matrixElement
matrixElement[n_ /; n \[Element] PositiveIntegers, expr_, m_ /; m \[Element] PositiveIntegers] := vacuumExpectation[(1/Sqrt[n!] GeneralizedPower[NonCommutativeMultiply, a, n]) ** expr ** (1/Sqrt[m!] GeneralizedPower[NonCommutativeMultiply, aD, m])]
matrixElement[0, expr_, m_ /; m \[Element] PositiveIntegers] := vacuumExpectation[expr ** (1/Sqrt[m!] GeneralizedPower[NonCommutativeMultiply, aD, m])]
matrixElement[m_ /; m \[Element] PositiveIntegers, expr_, 0] := vacuumExpectation[(1/Sqrt[m!] GeneralizedPower[NonCommutativeMultiply, a, m]) ** expr]
matrixElement[0, expr_, 0] := vacuumExpectation[expr]

and, as an example,

matrixElement[0, xx, 2]
(* ℏ/(Sqrt[2] m ω) *)

NonCommutativePolynomialReduce[xx, rels, gens, alg]
(* {{-((ℏ #1)/(2 m ω)) &}, -(ℏ/(2 m ω)) + (ℏ GeneralizedPower[NonCommutativeMultiply, a, 2])/(2 m ω) + (ℏ GeneralizedPower[NonCommutativeMultiply, aD, 2])/(2 m ω) + (ℏ a ** aD)/(m ω)} *)

toNormalOrder[expr_] := NonCommutativePolynomialReduce[expr_, rels, gens, alg]

NonCommutativePolynomialReduce[xx, rels, Reverse@gens, alg]
(* {{-((ℏ #1)/(2 m ω)) &}, -(ℏ/(2 m ω)) + (ℏ GeneralizedPower[NonCommutativeMultiply, a, 2])/(2 m ω) + (ℏ GeneralizedPower[NonCommutativeMultiply, aD, 2])/(2 m ω) + (ℏ a ** aD)/(m ω)} *)

toNormalOrder[expr_] := NonCommutativePolynomialReduce[expr, rels, gens, alg]