# General function for the expansion of a polynomial of operators

This question is motivated by a Quantum mechanical problem - but in explaining the problem - I assume no knowledge of quantum mechanics.

I want to define a function that can expand and simplify the following operator polynomial: $(\hat{a}^\dagger + \hat{a})^n\vert0\rangle$. The usual Expand function in Mathematica assumes commutativity of the two terms in the brackets, but of course here we have $[\hat{a}, \hat{a}^\dagger] = 1$, and so care with ordering must be used.

Indeed the expansion can be done easily since the expansion can be considered binary terms, i.e. (omitting the hats on the operators A and B):

(A+B)^n = AAA + AAB + ABA + ABB + BAA + BAB + BBA + BBB

which can be identified as the numbers from 0 to 7 if we re-label $A = 0$ and $B = 1$. Indeed I have also let $A = \hat{a}^\dagger$ and $B = \hat{a}$. Now since all of the terms in this expansion then operates on the vacuum state, we have the following properties:

1. Any terms ending with $B$ is zero and so has no contribution. This is since $\hat{a}\vert0\rangle = 0$.
2. Terms that have more $B$'s than $A$'s will also evaluate to zero.

So for example the terms $AAB$, $ABB$, $BAB$, $BBA$ and $BBB$ will have no contribution to the total evaluation of $(\hat{a}^\dagger + \hat{a})^n\vert0\rangle$. I would like to apply these rules for the expansion of the polynomial for general $n$ so that the output is greatly simplified.

There will be another further rule that must be considered for higher $n$ since there will be a greater variety of terms that the above rules do not consider. For example, for $n=5$, the following term will also be zero:

1. $AABBA\vert0\rangle=0$

The easiest way to think why this is the following rules of the annihilation and creation operators: $\hat{a}\vert0\rangle = 0$ (as written above) and $\hat{a}^\dagger\vert0\rangle = \vert1\rangle$. So the operator of $B$ on $\vert0\rangle$ minuses 1 from the index, whilst $A$ adds one to the index. If the evaluation, beginning from the right, yields at any point a negative number, the term has no contribution to the expansion of the polynomial. So in the above example (Rule 3), we start from right-hand side, The A increases 0 to 1, B then decreases this to 0 and the second-B from the right decrease this to -1: hence the term has no contribution.

• (1) It often pays to define rules on NonCommutativeMultiply for this sort of thing. (2) There is a section that shows ways of handling such operators here. There is a notebook version at library.wolfram.com. Mar 30, 2015 at 14:28

You have make your own definitions to handle non-commutativity. The basic trick is to use op[x1, x2, x3, ...] to represent an operator product, and to define various rules for manipulating the op[x1, x2, x3, ...].

To solve your problem I find that the minimum set of rules is as follows (using an obvious notation):

Clear[op];
op[u___, x_ + y_, v___] := op[u, x, v] + op[u, y, v];
op[u___, w_. op[x__], v___] := w op[u, x, v];
op[u___, a, ad, v___] := op[u, ad, a, v] + op[u, v];
op[u___, a, vac] = 0;

Define the basic operator.

Evaluate the first few powers of this operator acting on the vacuum:

NestList[op[s, #] &, op[s, vac], 4]

(*
{
}
*)

The question deals with a particular operator characterized by the commutation relation $[a, a^\dagger] = 1$. The polynomial to be expanded is also of a special form, $(a+a^\dagger)^n$. Also, the whole polynomial is supposed to be applied to the unique state vector $\vert 0\rangle$ that has the property $a \vert 0\rangle = 0$ (which is the null vector).

So this is a rather specialized function. The only thing that needs to be general here is that it should work for any power $n$. The rules in the question are derived from these very special circumstances. What I'm doing below is an approach that is in principle much more general because it would work with small modifications also for other polynomials and other state vectors to which the operators are applied.

You can solve this problem without using replacement rules. I'll just base the calculation purely on linear algebra, knowing what the matrix of the operator $a$ in the basis spanned by $\vert 0\rangle$, $a^{\dagger}\vert 0\rangle$, etc. looks like:

reduceHO[n_] :=
Module[{a = SparseArray[{Band[{1, 2}] -> Sqrt[Range[n]]}, {n, n} + 1]},
Total[(SuperDagger["a"])^Range[0,n]
MatrixPower[a + Transpose[a], n].UnitVector[n + 1, 1]/Sqrt[Range[0, n]!]]]

The argument of reduceHO is the power $n$ in the polynomial $(a+a^\dagger)^n$ to be expanded. Here is a test:

reduceHO[2]

$$\left(\text{a}^{\dagger }\right)^2+1$$

reduceHO[5]

$$\left(\text{a}^{\dagger }\right)^5+10 \left(\text{a}^{\dagger }\right)^3+15 \text{a}^{\dagger }$$

This is the expanded operator polynomial which is equivalent to the original one, when applied to the state $\vert 0\rangle$. It is just an attempt at a pretty output format - the detailed formatting could be customized further.

This approach is based on converting the operators to matrices, and then using MatrixPower to do the polynomial. The dimension of the matrix just has to be equal to the largest power of the raising operator $a^{\dagger}$ that can appear in the polynomial.

The basis states have normalization factors in them, which give rise to the Sqrt factors in the calculation of the matrix elements. These are divided out again at the end, leaving only the amplitude factors with which the different powers of the (unnormalized) vectors $(a^\dagger)^\nu\vert 0\rangle$ appear. The MatrixPower first yields a vector when applied to the UnitVector representing the ground state $\vert 0\rangle$, which is converted back to a symbolic expression inside the Total which sums up the list created from the amplitude factors and the symbolic powers $(a^\dagger)^\nu$, written as strings so that they remain inert. The last step could be modified depending on what you intend to do with the result.

Explanations

The matrix elements on which the power is based are taken from here. To show what the matrix has to do with the rules of the original question, here is its form for the case $n=6$:

With[{n = 6},
a = Normal@
SparseArray[{Band[{1, 2}] -> Sqrt[Range[n]]}, {n, n} + 1]];
MatrixForm[a]

$$\left( \begin{array}{ccccccc} 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & \sqrt{2} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & \sqrt{3} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 2 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & \sqrt{5} & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & \sqrt{6} \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \end{array} \right)$$

The state $\vert 0\rangle$ is now represented by the canonical unit vector $\{1, 0, 0, 0, 0, 0, 0\}$. This makes it clear that whenever $a$ is directly applied to $\vert 0\rangle$, the result is the null vector $\{0, 0, 0, 0, 0, 0, 0\}$, because the off-diagonal form of the above matrix shifts the entries of a vector one step to the left. It also multiplies the result by scale factors, but they are needed in order to insure the commutation relation:

MatrixForm[a.Transpose[a] - Transpose[a].a]

$$\left( \begin{array}{ccccccc} 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & -6 \\ \end{array} \right)$$

The last entry is not 1 because of the truncation to a finite dimension; the actual commutation relation can only hold in an infinite-dimensional space (in any finite dimension, taking the trace of $[a, a^\dagger] = 1$ on both sides would cause a contradiction). As long as one avoids such truncation effects by choosing the dimension of the matrix large enough, all the desired algebraic properties of $a$ follow from its matrix form above.

The other matrix appearing in the problem is $a^\dagger = a^T$, the Transpose of $a$. Its effect is to shift a vector in the opposite direction (multiplying it by scale factors again). This is why it's possible to predict that a product with more factors $a$ than $a^\dagger$ will map the state $\vert 0\rangle$ onto the null vector: the net shift is to the left.

Nice idea. I give this example to my 3rd year students—essentially Heisenberg's matrix mechanics.

Just a minor comment; you can use Dot instead of Total:

reduceHO[n_]:=
Module[
{a=SparseArray[{Band[{1,2}]->Sqrt[Range[n]]},{n,n}+1]},
((SuperDagger["a"])^Range[0,n]/Sqrt[Range[0,n]!]) .
MatrixPower[a+Transpose[a],n].UnitVector[n+1,1]]
• You probably meant this as a comment to my answer, but accidentally posted it in a new answer box. Thanks for the suggestion, anyway...
– Jens
Feb 5, 2016 at 7:36
• Yes—but in a comment the code does not format ... Feb 9, 2016 at 8:47