Creation & Annihilation Operators with Non-Commutative Algebra Mathematica 14.3

Here's how to use the NonCommutativeAlgebra functionality to implement an algebra with a finite number of generators and relations, illustrated by implementing the oscillator algebra. This will address Questions 1 and 2 in the OP (I'll wrap up at the end explaining why this won't work for a countably infinite number of generators or relations).

First define a set of generators, which are the raising and lowering operators $\hat{a}$ and $\hat{a}^{\dagger}$ for the algebra:

gens = {a, aD};

The defining relation for the algebra is that the commutator is the identity, i.e., $\hat{a}\hat{a}^{\dagger} - \hat{a}^{\dagger}\hat{a} -1 = 0,$ so we define a set of relations as

reps = {Commutator[a,aD] - 1};

This is a noncommutative Groebner basis:

NonCommutativeGroebnerBasis[rels, gens]
(* {-1 + a**aD - aD**a} *)

so it can be used to reduce noncommutative polynomials in the generators to a canonical form. Since I have specified the generators in the order {a, aD}, this canonical form will always be a noncommutative polynomial in which each monomial is of the form $(\hat{a}^{\dagger})^m\hat{a}^n$. Reversing the order of the generators reverses the order of these two factors in the canonical form.

Before I show these details, I need to define the algebra. The default algebra is

NonCommutativeAlgebra[]
(* NonCommutativeAlgebra[<|"Multiplication" -> NonCommutativeMultiply, "Addition" -> Plus, "Unity" -> 1, "Zero" -> 0, "CommutativeVariables" -> {}, "ScalarVariables" -> {}|>] *)

most of which is good. However, the only assumed scalars are explicit complex numbers, so I'll add some scalar variables that I know I'll need for later calculations for the harmonic oscillator algebra:

alg = NonCommutativeAlgebra[<|"ScalarVariables" -> {ℏ, m, ω}|>];

This is the minimum needed to do calculations with the oscillator algebra. As an example, let's look at the operator $\hat{x}^2$, which is

xx = (Sqrt[ℏ/(2 m ω)] (a + aD)) ** (Sqrt[ℏ/(2 m ω)] (a + aD));

To write this in canonical form, we reduce the polynomial relative to the generators, relations, and algebra using NonCommutativePolynomialReduce. Note that this function uses NonCommutativeExpand behind the scenes, so we don't need to call that function explicitly.

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

The second term is the remainder, and this is exactly the canonical form of $\hat{x}^2$: notice that the lowering operators are all on the right in each term. Note that if we reverse the order of the generators, we get the reverse with all the lowering operators on the left:

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

It might be worth packaging this up into a function:

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

so that we can instead call

toNormalOrder[xx]

to get the same things as before.

This does everything we need it to! But, we can add functionality with our rules. For instance, let's suppose we wanted to implement a vacuum expectation value. We can "abuse" NonCommutativePolynomialReduce to make this happen.

Clear@vacuumExpectation
vacuumExpectation[{{__Function}, expr_}] := NonCommutativePolynomialReduce[expr, gens, gens, alg]
vacuumExpectation[expr_] := vacuumExpectation[NonCommutativePolynomialReduce[expr, gbrel, gens, alg]][[2]]

The second part of this definition takes an expression which can be expanded into a noncommutative polynomial in our generators and reduces it relative to the generators and relations, leaving a sum of terms that are in normal order (lowering operators all to the right). Since the vacuum expectation value will render all non-constant terms zero in this canonical form, we can call NonCommutativePolynomialReduce again, but this time replace the relations with just the generators, which renders any non-constant term in the canonical form equal to zero, leaving just the constant term. This is what the first part of the definition does.

As an example, we can do

vacuumExpectation[x]
(* ℏ/(2 m ω) *)

which should be the right answer!

Note that to ensure things are more robust, it probably makes more sense to define toNormalOrder and vacuumExpectation in such a way that you feed in the generators and such explicitly. However, the trick of calling NonCommutativePolynomialReduce again with the generators in place of the relations already assumes a set of generators, relations, and canonical form. So perhaps we just need to be careful.

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

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.