Cubic perturbation to coupled quantum harmonic oscillators

I recently came across this two-dimensional problem of a particle in a potential of the form $$V = \displaystyle{\frac{1}{2}m \omega^2} \big(y^2 + x^2y \big) - \alpha y,$$ where $x$ and $y$ are known to be small, and I was trying to solve the time-independent Schroedinger equation,
$$-\frac{\hbar^2}{2m} \left( \frac{\partial^2 \psi}{\partial x^2}+\frac{\partial^2 \psi}{\partial y^2} \right) + (y^2 + x^2y ) \psi - \alpha y\, \psi = E \psi,$$
to obtain the wavefunction $\psi$ and the corresponding allowed energy levels.

My initial attempt at a solution was to treat the $x^2y$ term as a perturbation. Retaining just the $y^2 - \alpha y$ term, I obtain the zeroth-order wavefunction and energy levels as

$$E^{(0)}_n = \left( n + \frac{1}{2}\hbar \omega \right) + \frac{p_x^2}{2m},$$

$$\psi^{(0)}_{n} = \mathrm{e}^{i\, p_x / \hbar} \;\mathrm{H}_n \bigg(\sqrt{\frac{m \omega}{\hbar}} \rho \bigg) \mathrm{e}^{\displaystyle{\frac{-m \omega \rho^2}{2 \hbar}}},$$

where $\rho$ equals $y$ minus some easily calculable constant. However, the main problem arises when I try to calculate the first-order energy-levels $$E^{(1)}_n = \displaystyle{\frac{1}{2}m \omega^2}\langle \psi^{(0)}_n \,\lvert \, x^2y \,\rvert\, \psi^{(0)}_n \rangle.$$ This integral obviously cannot be evaluated since the $x^2$ term, when integrated from $-\infty$ to $\infty$, diverges. Moreover, since the additional term is cubic, I also cannot employ a change of coordinates to decouple the system.

Hence, I was wondering if there exists some clean method by which one can solve this problem. I would really appreciate any help in this regard. Thanks!

• If I see this correctly, the third order term is not bound from below and there is no ground state. One could probably still calculate the decay time of a wave packet, as it "leaks" out of the potential? – CuriousOne Dec 15 '14 at 8:40
• It is indeed true that for the third-order term alone, there is no ground state. However calculation of the decay time, as done herein (arxiv.org/pdf/quant-ph/0703234.pdf) would involve the initial length of the box in which the particle is confined. What if this is on an unbounded domain? – ACuriousMind Dec 15 '14 at 8:45
• I think it would make sense to perform time dependent perturbation theory and to turn the third order term on at $t>0$ using the unperturbed solutions as initial wave function. – CuriousOne Dec 15 '14 at 8:55
• @CuriousOne, I am not very familiar with time-dependent perturbation theory so could you please expand upon your comment in an answer? It would be really very helpful if you could explain how to proceed with the problem. – ACuriousMind Dec 15 '14 at 9:11
• You can evaluate the action of the perturbation terms on the eigenstates of the unperturbed solutions. Since the ground state decays very quickly in both x and y directions, the integrals should all be convergent (just guessing here). That will give you the time scale on which those states decay into free (in this case accelerated along the y-axis?) solutions. – CuriousOne Dec 15 '14 at 9:26