Let us consider the classical Harmonic Oscillator (a selfadjoint operator)
$$
\mathcal H=\frac12\left(-\frac{d^2}{dx^2}+x^2\right),\quad\text{with spectrum $\frac12+\mathbb N$.}
$$
This one-dimensional operator has simple eigenvalues and the eigenvector corresponding to $\frac12+n$ is the Hermite function $h_n$. The $(h_n)_{n\in \mathbb N}$ make a Hilbertian basis. Note that each $h_n$ is an entire function.
Now let's work formally and consider an operator with symbol
$$
\frac12(\xi^2+zx^2),\quad\text{$z$ is a parameter $\in \mathbb S^1, z\not=-1.$}
$$
Assuming $z=e^{2i\theta}, \vert\theta\vert<π/2$, we make the formal symplectic change of variable given by
$$
\xi=e^{i\phi/2}\eta,\ x=e^{-i\phi/2} y,\quad
\xi^2+zx^2=e^{i\phi}\bigl(\eta^2+y^2e^{-2i\phi+2i\theta}\bigr)=
e^{i\theta}\bigl(\eta^2+y^2\bigr),
$$
if $\phi=\theta$. It is now easy to verify directly that
$h_n(e^{i\theta/2} y)$ is an eigenvector of the operator
$$
\mathcal H_z=\frac12\left(-\frac{d^2}{dx^2}+zx^2\right),
$$
related to
the eigenvalue $e^{i\theta}(\frac12+n)$ and that the spectrum of $\mathcal H_z$ is thus $e^{i\theta}(\frac12+\mathbb N)$.
Of course the new eigenvectors are no longer a Hilbertian basis for lack of orthogonality but they generate $L^2(\mathbb R)$ as the $h_n$. It is somehow a piece of luck that the $h_n$ are entire functions so that it makes sense to evaluate them for complex arguments. This method is called analytic deformation and can be generalized in various geometric situations. For your precise question, you thus just have to use the classical result for $b=0$.