I had recently been studying the Dirac equation, and as an example of how the equation is used, I was given a problem about the Landau levels of graphene (but I personally have no knowledge about graphene).
I was given the Hamiltonian
$$H = \begin{pmatrix} 0 & v_F(p_x - i(p_y + \frac{e}{c}B x)) \\ v_F(p_x - i(p_y + \frac{e}{c}B x)) & 0 \end{pmatrix}$$
using the Landau gauge $\textbf{A} = (0,Bx,0)$ ($B$ has only a z-component) and $v_F$ is the Fermi velocity. And I would like to solve for the energy levels and the 2-component wave functions.
I personally solved this problem by explicitly multiplying out the matrix and found that this Hamiltonian gives two decoupled harmonic oscillator equations. But the official solution solved it by introducing a pair of ladder operator right from the beginning
$$a = \sqrt{\frac{c}{2e\hbar B}}(p_x - i(p_y + \frac{e}{c}B x)) $$ $$a^\dagger = \sqrt{\frac{c}{2e\hbar B}}(p_x + i(p_y + \frac{e}{c}B x)) $$
and rewriting the Hamiltonian as
$$ H = \sqrt{\hbar \omega}\begin{pmatrix} 0 & a \\ a^\dagger & 0 \end{pmatrix} , \ \omega = \frac{2eBv_F^2}{c}$$
Then the energy levels are found to be $E = \pm \sqrt{n \hbar \omega}$ by exploiting the ladder operators.
While I fully understood how to use the algebra of $a$ and $a^{\dagger}$ to obtain the wave function and the energy levels, but I just have no clue how one motivates them in the first place. More precisely, I would like to ask
How can I be sure that the $a$ and $a^{\dagger}$ I introduced here will actually behave like in the harmonic oscillator problem, i.e they are indeed the same raising and lowering operators? (except from the fact that the Hamiltonian is actually a pair of harmonic oscillator when multiplied out)
How do I decide the constant in front of my $a$ and $a^{\dagger}$? The choice seems really arbitrary to me, and the official solution did not explain how it come up with this constant factor.
I had also been given a second problem for a massive Hamiltonian ("effective mass" $m^*$)
$$H = -\frac{1}{2m^*}\begin{pmatrix} 0 & (p_x - i(p_y + \frac{e}{c}B x))^2 \\ (p_x - i(p_y + \frac{e}{c}B x))^2 & 0 \end{pmatrix}$$
and I would like to know do the same method applies directly to this Hamiltonian as well (using the same $a$ and $a^{\dagger}$ as above?) ? In this case when the Hamiltonian is multiplied out, it will no longer gives a harmonic oscillator, but something of power 4, and I don't see how $a$ and $a^{\dagger}$ can still be useful.
I had nonetheless tried it out, and got the energy levels to be $E = n\hbar \omega$, where $\omega = \frac{eB}{m^*c}$ and the wave function $\psi = (|n>, -(n-1)|n-2>)$, where $|n>$ are the n-th harmonic oscillator wave function. It will also be very helpful if I can know whether my results are correct.
This is the first time I encountered this kind of method, on top of the already difficult enough topic of Dirac equation, and I am really confused. Thanks for any help.