How can I identify momentum eigenstates in a tight binding model with degenerate energy eigenstates?

Question

Summary:

I numerically diagonalize a tight binding Hamiltonian to get energy eigenvectors, some of which are degenerate. However, the numerically diagonalized degenerate eigenvectors are not necessarily eigenvectors of momentum as well. How can I numerically get eigenvectors of both momentum and the Hamiltonian?

More Details:

I am using a single orbital tight binding model in graphene, and applying a magnetic field so that Landau Levels appear. Specifically, I am using the tight binding model

$$H = t \sum_{r, r' \in NN}a^{\dagger}(r)b(r')e^{i\theta_{r, r'}} + h.c.$$

where $\theta_{r,r'} = q \int_r^{r'}A(r) \cdot dr$. A is the vector potential, and q is the particle charge. $a^\dagger, a, b^\dagger, b$ are the creation and annihilation operators on the A and B sublattices of graphene. Picking the Landau Gauge, $A = -B y \hat x$, and $\theta$ becomes $ \theta_{r, r'} = \frac{1}{2} q B (x'-x)(y'+y)$

I want to find exact solutions to this Hamiltonian, so I numerically diagonalize the tight binding model on a lattice with periodic boundary conditions. As expected, I get Landau levels. Of course, each Landau level has high degeneracy. Therefore when I numerically diagonalize the tight binding Hamiltonian, I get many eigenvectors with the same eigenvalue.

Translational symmetry in the x-direction is maintained, so $k_x$ should be a good quantum number. Additionally, translation by $\Delta_y = \frac{4\pi}{qB}$ leave the Hamiltonian invariant as well, so $k_y$ should be a good quantum number as well, although the Bruillon zone for valid $k_y$ is reduced from the case where there is no magnetic field.

When I numerically diagonalize the tight binding Hamiltonian, I get eigenvectors of the Hamiltonian, but not necessarily eigenvectors with quantum numbers $k_x$ and $k_y$. Is there a way for me to numerically find simultaneous eigenvectors of $p_x$, $p_y$, and $H$?

Alternatively, is there a way to express the momentum operator $p_x$ and $p_y$ in the tight binding basis? If I can express the operators $p_x$ and $p_y$ in the tight binding basis, I can easily find simultaneous eigenvectors of $p_x$, $p_y$, and $H$ numerically.

I am aware that the low-energy graphene Hamiltonian can be expanded around the K and K' points, and exact solutions for the landau levels as a function of $k_x$ can be obtained. However, I need full, exact solutions to the complete tight binding Hamiltonian so this is not helpful for me.

Answer 1

Recall that when you diagonalize your Hamiltonian, that corresponds to changing basis into that of $a_E, a^\dagger_E$, which are the annihilation and creation operators of particles with energy $E$, meaning your Hamiltonian becomes

$$H = \sum_E E\ (a^\dagger_E a_E + b^\dagger_E b_E).$$

Similarly you can find your momentum operator as $p_a = \sum_{p} p\ a^\dagger_p a_p$, and similarly for your $b$ particles. If you then Fourier transform these operators from momentum to position you will have your momentum operator in the tight binding basis.