Why harmonic oscillator levels for noninteracting electrons in 2D in an applied magnetic field?

Why should harmonic oscillator levels arise out of the blue?

Edit after reading new comments: if you write down the Hamiltonian of the system using the standard minimal coupling of the momentum $vec{p} to vec{p} + q vec{A}$ and you write down $vec{A}$ for a uniform magnetic field using the gauge that you prefer, you get exactly the Hamiltonian of a harmonic oscilator. You should not be surprised, cause circular motion is a two-dimensional harmonic motion!

Is there a relation between classical circular trajectories and quantum harmonic oscillator levels?

Yes! The classical motion of a particle of mass $m$ and electric charge $q$ in an external magnetic field $B$ is given by Newton's second law: $$ m frac{v^2}{r} = q v B ;;; rightarrow ;;; r = frac{mv}{qB}. $$ Let me write this radius in terms of kinetic energy $E=frac{1}{2}mv^2$, which is more convenient for comparison with quantum motion and let me write $r^2$ rather than $r$: $$ r^2 = frac{2mE}{(qB)^2}. $$

Now in quantum mechanics you can prove that the energy eigenvalues are $$ E_n = frac{hbar q B}{m} left( n + frac{1}{2} right), $$ where $n$ is a non negative integer labeling the levels; and the corresponding wave function is $$ psi_n(x,y) = N_n (x + iy)^n e^{-qB(x^2+y^2)/4hbar}, $$ $N_n$ being a normalization constant.

Now if you compute the expectation value of the operator $r^2 = x^2+y^2$ you get the following result $$ leftlangle r^2 rightrangle = frac{2hbar}{qB} ( n + 1). $$ You can compare the classical result with the quantum result at high energies, so you can assume $n gg 1$ and simplify $ leftlangle r^2 rightrangle approx frac{2hbar n}{qB}$, and similarly $E = frac{hbar q B n}{m}$. If you combine these two results to get rid of $n$, you immediately get $$ leftlangle r^2 rightrangle = frac{2 m E}{(qB)^2}, $$ which is exactly the classical result!