In the natural units, \(\hbar=c=k=1\), the Dirac equation in the presence of the vector potential reads \[\begin{equation} \tag{1} \left( \begin{array}{cc} m & \vec{\sigma}\cdot\vec{\pi} \\ \vec{\sigma}\cdot\vec{\pi} & -m \\ \end{array} \right) \left( \begin{array}{c} \chi \\ \phi \\ \end{array} \right)=E \left( \begin{array}{c} \chi \\ \phi \\ \end{array} \right), \end{equation}\] where \(\vec{\pi}= \vec{P}-q \vec{A}\) . Multiplying the matrices we have two equations: \[\begin{equation} \tag{2} m\chi+\vec{\sigma}\cdot\vec{\pi}\phi=E\chi, \quad \text{ and } \quad \vec{\sigma}\cdot\vec{\pi}\chi-m\phi=E\phi. \end{equation}\] From the second equation we get \(\phi=\frac{\vec{\sigma}\cdot\vec{\pi}\chi}{E+m}\) . Plugging this back into the first one we get the equation for \(\chi\): \[\begin{equation} \tag{3} (E^2-m^2)\chi=\left(\vec{\sigma}\cdot\vec{\pi}\right)^2\chi. \end{equation}\] We will use the identity \[\begin{equation} \tag{4} \vec{\sigma}\cdot\vec{A}\vec{\sigma}\cdot\vec{B}= I \vec{A}\cdot\vec{B}+i( \vec{A}\times\vec{B})\cdot \vec{\sigma}. \end{equation}\] For our case \(\vec{A}=\vec{B}=\vec{\pi}\) . Note that, unlike regular vectors, \(\vec{\pi}\times\vec{\pi} \neq 0.\) This is because \(\vec{\pi}\) contains a function of \(x\), \(\vec{A}(x)\) , which does not commute with \(\vec{P}\) . Let’s calculate what it is using the \(\epsilon_{i j k}\) symbol: \[\begin{eqnarray} (\vec{\pi}\times\vec{\pi})_i&=&\epsilon_{i j k}\pi_j\pi_k=\epsilon_{i j k}\left[P_j(-q A_k)+(-q A_j)P_k\right] =i q \epsilon_{i j k}\left[(\partial_jA_k)+A_k\partial_j+q A_j\partial_k\right]=i q(\vec{\nabla}\times\vec{A})_i =iq B_i \tag{5}, \end{eqnarray}\] where we dropped symmetric terms \(P_jP_k\) , \(A_j A_k\) and \(A_k\partial_j+A_j\partial_k\) since they are contracted with \(\epsilon_{i j k}\). Also note that in this notation, operators act on everything on the right side unless they are in parenthesis, that is \(\partial_jA_k=(\partial_jA_k)+A_k\partial_j\) .) Then the equation for \(\chi\) becomes \[\begin{equation} \tag{6} (E^2-m^2)\chi=\left(I \vec{\pi}^2-q\vec{B}\cdot\vec{\sigma}\right)\chi. \end{equation}\] We choose \(\vec{A}=- y B_0 \hat{i}\) for which Eq. (6) becomes \[\begin{equation} \tag{7} (E^2-m^2)\chi=\left[I ((P_x+q y B_0)^2+P_y^2+P_z^2)-q B_0\sigma_z\right]\chi. \end{equation}\] Note that Eq. (7) has no \(x\) and \(z\) dependence which suggests that the solutions must be plane waves along the \(x\) and \(z\) axis, that is \(\xi= \zeta \,e^{i p_x x+i p_z z}\) . The new equation in \(\zeta\) becomes \[\begin{equation} \tag{8} \left[I\left((P_x+q y B_0)^2-\frac{d^2}{dy^2}-E^2-m^2+P_z^2\right)+q B_0\sigma_z\right]\zeta=0. \end{equation}\] Note the changes: the operators \(P_x\) and \(P_z\) became \(p_x\) and \(p_z\) because of the plane wave expansion. Now we rewrite the equation once more, \[\begin{equation} \tag{9} \left[I\left(\frac{d^2}{dy^2}-(p_x+q y B_0)^2+E^2-m^2-p_z^2\right)+q B_0\sigma_z\right]\zeta=0. \end{equation}\] In terms of the components we have \[\begin{equation} \tag{10} \left( \begin{array}{cc} \frac{d^2}{dy^2}-(p_x+q y B_0)^2 &0 \\ +E^2-m^2-p_z^2+q B_0 & \\ & \frac{d^2}{dy^2}-(p_x+q y B_0)^2 \\ 0 & +E^2-m^2-p_z^2-q B_0 \\ \end{array} \right)\left( \begin{array}{c} \zeta^+ \\ \zeta^-\\ \end{array} \right)=0. \end{equation}\] We have two decoupled equations which differ by the sign of the \(B_0\) term only. We can actually combine them into one form introducing a parameter which simply encodes that sign, namely \[\begin{equation} \tag{11} \left( \frac{d^2}{dy^2}-(p_x+q y B_0)^2 +E^2-m^2-p_z^2+s q B_0\right)\zeta^s=0, \end{equation}\] where \(s=\pm\) . We are almost there since our equation is familiar; it is like the differential equation for Hermite polynomials. We can reveal that it is indeed the case if we define a dimensionless parameter \(\xi=\sqrt{q B_0}(y+\frac{p_x}{qB_0})\) . In this new variable we have, \[\begin{equation} \tag{12} \left( \frac{d^2}{d\xi^2}-\xi^2 +a_s\right)\zeta^s=0, \end{equation}\] where \[\begin{equation} \tag{13} a_s=\frac{E^2-m^2-p_z^2+s q B_0}{|q|B_0}. \end{equation}\] From the harmonic oscillator solution we know that Eq. (12) admits normalizable solutions only if \(a_s= 2 \nu+1\) where \(\nu=0,\,1,\,2,\cdots\) . Finally we solve Eq. (13) for \(E\) to get, \[\begin{equation} E=\sqrt{(2 \nu+1)|q|B_0+m^2+p_z^2-s q B_0}.\tag{14} \end{equation}\] Some comments on our final result: first note that there is a continuous parameter \(p_z\) in the answer, so we can change the energy continuously. But this is the trivial part of the result since nothing depends on \(z\) it is obvious that energy will change with \(p_z^2\) in the way the relativistic energy-momentum relation requires. Assume \(p_z\) is fixed. There are two parameters left, \(\nu\) and \(s\) . \(\nu\) labels the level of the harmonic oscillator, and \(s\) labels the spin. For each \(\nu\), there are two levels of energy corresponding to spin up and spin down states. Note that the energy is dominated by the mass term for non relativistic particles. Let’s expand it out to linearize it: \[\begin{equation} E=m \sqrt{1+(2 \nu+1)|q|B_0/m^2+p_z^2/m^2-s q B_0/m^2}\simeq m+ \frac{(2 \nu+1)|q|B_0}{2m} -\frac{s q}{2 m} B_0\tag{15}. \end{equation}\] \(\frac{s q}{2 m}\) term in front of \(B_0\) is the dipole moment of the fermion, and the \(1/2\) factor is correctly produced by the Dirac theory.