Dispersion of finite 2D lattice Following problem: I have the coupling matrix for an $N$-by-$N$ finite lattice of coupled masses (only nearest-neighbour coupling, periodic until terminated). I would like to numerically calculate its band structure.
Because it is 2D, I have the matrix as a pentadiagonal $N^2 \times N^2$ matrix (as opposed to tridiagonal for the 1D).
Diagonalizing it is straightforward but how do I actually assign a momentum to each of the eigenvalues I get out? Do I need to look at every single eigenmode individually or is there a way to assure that the diagonalized coupling matrix comes out in some predictably sorted way that allows me to know which eigenvalue belongs to which momentum?
 A: I assume that your $N\times N$ lattice has periodic boundary conditions and you have translation symmetry. If so, you have two translation operators $T_x$ and $T_y$ that commute with your matrix $H$ (and each other) and thus have common eigenvectors. The eigenvectors of the translation operator are labeled by momenta $k$.
In particular $T_x^N=1$, and you therefore have eigenvalues $e^{\frac{2\pi i}N k_x}$ where $k_x=0,\dots, N-1$. Similarly for $T_y$. If you have an eigenvector of $H$
$$ Hv_\lambda = \lambda v_\lambda$$
with no degeneracy, then you can measure the momenta as follows
$$ T_xv_\lambda = e^{\frac{2\pi i}N k_x}v_\lambda \qquad T_yv_\lambda = e^{\frac{2\pi i}N k_y}v_\lambda.$$
This eigenvector $v_\lambda$ has momenta $(k_x, k_y)$. If you have a degeneracy, ie several eigenvectors with eigenvalue $\lambda$: $$v^{(1)}_\lambda, \dots, v^{(r)}_\lambda$$
then these $H$ eigenvectors are not necessarily eigenvectors of $T_x$ and $T_y$. In order to momentum-resolve these, you need to find the correct superposition of these. A straightforward, although probably not the most numerically efficient, way to do this is to construct projection operators.
For example something like
$$ P_x(k_x) = \frac 1N\sum_{n=1}^N e^{-\frac{2\pi i}N k_x n} T_x^n. $$
If I wrote it correctly, this should satisfy $P_x(k_x)^2=P_x(k_x)$ and be Hermitian (thus is a projection). Furthermore you have something like
$$ T_xP_x(k_x) = e^{\frac{2\pi i}N k_x}P_x(k_x).$$
This means that for any vector $v$, $P(k_x)v$ is either zero or an eigenstate of $T_x$ with eigenvalue $^{\frac{2\pi i}N k_x}$.
You have a similar projector for $T_y$. You can thus momentum-resolve a degenerate set of vectors $v^{(1)}_\lambda, \dots, v^{(r)}_\lambda$, by acting with $P_x(k_x)$ and $P_y(k_y)$ on these vectors until you find $r$ linearly independent vectors. These will all have well-defined momentum.
An easier and numerically more efficient way to do this might be the following. Instead of diagonalizing $H$, you should instead diagonalize
$$ H' = H + \delta_x T_x + \delta_y T_y,$$
for some very small numbers $\delta_x$ and $\delta_y$. This will make it more likely that the computer spits out momentum resolved eigenvectors, and you can thus avoid projectors on degenerate state.
