I would like to use sage to solve problems for the 1-dimensional quantum harmonic oscillator. The algebraic solution to the problem is to define operators $a$ and $a^{+}$ that satisfy the commutator $[a,a^{+}]$=1. The Hamiltonian can then be written as $$ H = \hbar \omega (a^+ a +1/2). $$ What I would like to do with sage is simplify expressions like $a(a^{+})^n$ by using the commutator to write $aa^{+}=1+a^{+}a$, and use this repeatedly to 'move' $a$ to the right, through all of the $a^{+}$ operators. In other words, I would like to do the following calculation automatically.
$$ a(a^{+})^n = aa^+(a^+)^{n-1} $$ $$ = (1+a^{+}a)(a^+)^{n-1} $$ $$ = a^+(1+aa^{+})(a^+)^{n-2} $$ $$ ... $$ $$ = n(a^+)^{n-1}+(a^+)^na$$
This type of manipulation is necessary to calculate normalization coefficients for the wave functions. I have searched extensively but still don't know how to set up a simple algebra with the $a$'s like this with sage. (Note that sympy has some functionality to do this, but I want to use the algebraic apparatus within sage, because I want to explore more sophisticated algebras once I understand things better.) Any help is greatly appreciated!
