I wonder if there is anyone or any references that have solved the Landau level spectrum and eigenstates with respect to the following Hamiltonian:
\begin{equation} H=\frac{k_x^2-k_y^2}{m}\sigma_x+\frac{2 k_x k_y}{m}\sigma_y \end{equation}
when coupled to external magnetic field in z-direction either in Landau gauge or symmetric gauge.
I find the answer in papers that studies bilayer grapheme, e.g.
and I decided to write the answer to my own question. First we do minimal coupling to the magnetic field:
\begin{equation} H[A]=\frac{(k_x+e A_x/c)^2-(k_y+e A_y/c)^2}{m}\sigma_x+\frac{(k_x+e A_x/c)(k_y+e A_y/c)+(k_y+e A_y/c)(k_x+e A_x/c)}{m}\sigma_y \end{equation}
Notice that $k_x k_y$ should be symmetrized when replaced with canonical momentum in order to keep the Hermicity of the Hamiltonian. In Landau gauge,
\begin{equation} A_x=-B y, A_y=0 \end{equation}
then
\begin{equation} [\frac{(-i \partial_x-\frac{e B}{c} y)^2-(-i \partial_y)^2}{m}\sigma_x+\frac{(-i \partial_x-\frac{e B}{c})(-i \partial_y)+(-i \partial_y)(-i \partial_x-\frac{e B}{c})}{m}\sigma_y]\psi(\mathbf{r})=E_n \psi(\mathbf{r}) \end{equation}
Due to translational invariance in x-direction,
\begin{equation} \psi(\mathbf{r})=\frac{1}{\sqrt{L}}exp[i k x]\hat{f}_n(y) \end{equation}
one find
\begin{equation} [\frac{(k-\frac{e B}{c} y)^2-(-i \partial_y)^2}{m}\sigma_x+\frac{(k-\frac{e B}{c})(-i \partial_y)+(-i \partial_y)(k-\frac{e B}{c})}{m}\sigma_y]\hat{f}(y)=E_n \hat{f}(y) \end{equation}
Defining the creation and annihilation operator as
\begin{equation} a^-=l_B \partial_y+(l_B k-\frac{e B}{c} y/l_B), a^+=l_B \partial_y-(l_B k-\frac{e B}{c} y/l_B), \end{equation}
where $l_B=\sqrt{\frac{c}{e B}}$ is the magnetic length. We have
\begin{equation} \omega_c \begin{bmatrix} 0 & {a^+}^2 \\ {a^-}^2 & 0 \end{bmatrix} \begin{bmatrix} f_n^+(y) \\ f_n^-(y) \end{bmatrix}=E_n \begin{bmatrix} f_n^+(y) \\ f_n^-(y) \end{bmatrix} \end{equation} The spectrum and eigenstates can be solved in analogy with Harmonic oscillator problem:
\begin{equation} E_n^{\pm}=\sqrt{n(n-1)}\omega_c,n=2,3,\ldots, \hat{f}_{n,\pm}=\frac{1}{\sqrt{2}} \begin{bmatrix} \phi_n(y) \\ \pm\phi_{n-2}(y) \end{bmatrix} \end{equation}
where $\phi_n(y)$'s are eigenstates of Harmonic oscillators.
