Constant magnetic field applied to a quantum harmonic oscillator

You should be using cylindrical coordinates because the symmetry of the problem demands it. Doing anything else is just playing hard-headed.

You should split the hamiltonian into a $z$ component and a cylindrical radial component: $$ \begin{align*} \hat {\mathcal H}_{\rho} &= \frac{\hat p_x ^2+\hat p_y^2}{2m} + \frac{1}{2}m \left ( \omega_\text{c}^2 + \omega_0 ^2 \right)\left(\hat x^2+\hat y^2\right)-\omega_\text{c} \hat L_z,\\ \hat {\mathcal H}_{z} &= \frac{\hat p_z ^2}{2m} + \frac{1}{2}m \omega_0 ^2 \hat z^2, \end{align*} $$ where $\omega_\text{c}=\frac{qB}{2m}$ is the relevant cyclotron frequency. To solve this without any calculation:

• The $z$ component is immediate.

• The radial component has a 2D central potential $\mathcal H_{\rho,0}=\frac{\hat p_x ^2+\hat p_y^2}{2m} + \frac{1}{2}m \left ( \omega_\text{c}^2 + \omega_0 ^2 \right)\left(\hat x^2+\hat y^2\right)$ which commutes with $\hat L_z$, so you can simply add the energies: that is, $\mathcal H_{\rho,0}$ has a common eigenbasis with $\omega_\text{c} \hat L_z$, and you can diagonalize them simultaneously; this 2D harmonic oscillator is relatively standard but you can still solve it explicitly.

  To do that, you can find the eigenenergies of $\mathcal H_{\rho,0}$ by splitting it into two 1D harmonic oscillators, but the product basis of those oscillators does not give you an angular-momentum eigenbasis. Instead, the thing to do is to note that the set of eigenfunctions with a fixed total excitation number $n$, i.e. the subspace $$\mathcal L_n=\mathrm{span} \left\{|n_x\rangle \otimes |n_y\rangle : n_x+n_y=n\right\}$$ is invariant under rotations (because the hamiltonian is) so you can diagonalize $\hat L_z$ in this subspace.

  How do you do that, then? Well, you notice that your basis functions, of the form $H_{n_x}(x)H_{n_y}(y)$ once you remove the gaussian, are polynomials of degree $n$ in $x,y$ and with definite parity, and you're looking to decompose them as linear combinations of $(x\pm iy)^m$ for $m=\ldots,n-4,n-2,n$, which are the eigenfunctions of $\hat L_z$.

  That then tells you what eigenvalues $m$ of $\hat L_z$ are allowed for each eigenspace of $\mathcal H_{\rho,0}$ with excitation number $n$, i.e. $m\leq n$ with both of the same parity. That's sufficient to fix the spectrum.