numerical diagonalization of tight-binding hamiltonian I would like to find the exact eigenvalues of the following tight-binding Hamiltonian, written here in second quatization:
\begin{eqnarray}
\hspace{-0.25in}{\mathcal{H}} &=& \mathcal{H}_0+ \frac{1}{2} \kappa a^2 \sum_{m, n}  (m^2+n^2) \hat{a}^\dagger_{m,n} \hat{a}_{m,n} ,  \nonumber \\ 
\hspace{-0.25in}\mathcal{H}_0 &=& - J \sum_{m, n}  \left( \hat{a}^\dagger_{m+1,n} \hat{a}_{m,n}+ e^{i \phi }  \hat{a}^\dagger_{m, n+1} \hat{a}_{m,n} \right) + \mbox{h.c.}\label{eq:hof}
\end{eqnarray}
How can I do so numerically? I know how to diagonalize the Hamiltonian matrix, but the question is how to get to that in the first place, i.e. how do I compute the matrix elements?
Edit: to clarify, I am only interested in finding the eigenvalues in the subspace of single-particle states.
 A: You can try just working out the matrix elements between the different localised states, i.e. $a_1^{\dagger}\ldots a_m^{\dagger}|0\rangle$ etc. This may be rather tedious and require you to carefully keep track of phases due to the fermionic anti-symmetry. Probably a better way is to use the Jordan-Wigner transformation, which gives you an explicit matrix representation of the fermion operators. In 2D this representation of the Hamiltonian looks ugly, but it is still perfectly useful. You will be restricted to a fairly small lattice size, i.e. 10-20 sites depending on your computer and the efficiency of your code. It will be more efficient if you use the fermion number conservation to diagonalise within each number sector separately.
In the particularly simple case of a single particle on the lattice, the problem becomes rather easy, since exchange statistics are no longer relevant. In this case, you make the replacement
$$a_{i,j}^{\dagger}a_{m,n} \to |i,j\rangle\langle m,n|,$$
where 
$$|i,j\rangle = a_{i,j}^{\dagger}|0\rangle.$$
