I am studying lattice field theory and would like to understand the momentum operator for free Dirac fermions on a square lattice.
In this case one needs to discretize the momentum operator (which otherwise would have a continuous spectrum):
$$P_\mu^\text{(free)} = - i\hbar \int \text{d}x^3 \psi(x)^\dagger\partial_\mu \psi(x) \to -i\hbar\sum_x \psi_x^\dagger \left( \dfrac{\psi_{x+\hat{\mu}} - \psi_{x-\hat{\mu}}}{2 a} \right)$$
Where $\psi$ is a Dirac spinor, $x$ denotes the lattice site, $a$ the lattice constant, $\hat{\mu}$ the unit vector in direction $\mu$ and I used central differences to discretize the derivative.
My question now is: What is the spectrum of this $P_\mu$ on the lattice (assuming periodic boundary conditions)?
My take on it: I would expect that since the phase of the Dirac Fermion is ~$e^{ip \cdot x}$ ($\hbar = 1$), then the periodic boundary conditions on the lattice require for each dimension: $e^{ip_\mu L} = 1$. Here $L = N_\mu a$ denotes the size of my lattice.
From this I would deduce $p_\mu \in \left( - \frac{\pi}{a}, -\frac{\pi}{Na} (n-1), ... , \frac{\pi}{a} \right)$ with integers $n$ as is typical in condensed matter physics. However, when I try to calculate the spectrum of the discretized $P_\mu$ numerically for small lattices, I get different results (sometimes the spectrum is not even linear). Where did I go wrong? Is my expected spectrum correct?