Functions of ladder operators - Identities for or methods of solution I would like to find the eigenstates of a potential with terms like
$$
\left(\frac{1}{4}A^{4}x^{4}-A^{2}x^{2}+1\right)^{2}
$$
and am planning to use ladder operators to find the solution in the basis of harmonic oscillator eigenstates. In that case (absorbing new constants into $B$) I must work with terms like
$$
\left(\frac{1}{4}B^{4}\left(a^{+}+a\right)^{4}-B^{2}\left(a^{+}+a\right)^{2}+1\right)^{2}
$$
I can just expand this out and work from there, but that carries with it the risk of computational errors which may be difficult to track down. It would be safer to use an identity if one is available.
Is there a better method to approach this than just brute-forcing it? I thought there was an identity for the action of a function of ladder operators upon a given state $|j\rangle$, but if there is, I can't seem to remember it.
EDIT: just to clarify, I do not expect a closed-form solution. I just want to expand in the basis of states $|j\rangle$ and diagonalize in a truncated (e.g. $j=0,1,2,...100$) basis. If I can figure out the matrix elements of terms like this, i.e. $\langle k|V|j\rangle$ where $V$ is the potential term above, that's all I need.
 A: You need to think about the wavefunction $\psi(x) = Ae^{-x^2}$, and the momentum operator definition $p=-i\frac{d}{dx}$. If you apply the momentum operator on the wavefunction, you get Hermite polynomials (which are in the form of the potential). The question wants the expectation value of the function in x basis:
$$<|V|> =A^2 \int_{-\infty }^{\infty}e^{-\alpha x^2}V(x)e^{-\alpha x^2}dx$$ Where $\alpha = \frac{m\omega}{2 \hbar}$. You will see that multiple applications of the momentum operator on the wavefunction will produce the potential in question. Off the top of my head, $H_2(x) = 4x^2-2$, which results from the momentum operator applied to the wave function twice. if raised to the 4th power will produce something similar to your problem. If you put the momentum operator in terms of $p=\frac{i\hbar \omega}{2}\sqrt{a^\dagger - a }$ then:
$$A^2 \int_{-\infty }^{\infty}e^{-\alpha x^2}(a^\dagger-a )^4 e^{-\alpha x^2}dx$$
When you expand this, only the $a^\dagger a$ terms will survive for i=j, and you will get some coupling terms too. 
