How should I interpret the eigenvalues of a momentum operator discretised with a central-difference scheme? I have a simple 1D momentum operator $P = -i\frac{\partial}{\partial x}$. I discretise it on a grid with spacing $h$ with a central-difference scheme like so: $\left[PF\right]_{x=jh} = \left[\mathbf{P}\cdot\mathbf{F}\right]$ where
$$\mathbf{P} =
\begin{bmatrix}0 & \frac{-i}{2h} && \cdots& \frac{i}{2h}\\
\frac{i}{2h}& 0 &\frac{-i}{2h}&&\\
&\frac{i}{2h}&0&\frac{-i}{2h}&\\
&&\frac{i}{2h}&0&\frac{-i}{2h}\\
\frac{-i}{2h}&\cdots&&\frac{i}{2h}&0\end{bmatrix}$$
When I diagonalise this Hamiltonian, I get pairs of degenerate eigenvalues. However, I expected eigenvalues approximating $\lambda_n = \frac{2\pi n}{L}$
What is the meaning of the degenerate eigenvalues?
[edit] - I must have made a mistake somewhere. When I wrote the script from scratch the degenerate eigenvalues went away
import numpy as np                                                                           
                                                                                              
def main():                                                                                  
                                                                                              
    n = 50                                                                                   
    L = 1.0e0                                                                                
    h = L/n                                                                                  
                                                                                              
    # Building momentum operator                                                             
    #                                                                                        
    P = np.zeros((n, n), dtype=np.cdouble)                                                   
    for i in range(n-1):                                                                     
        P[i,i+1] = -0.5e0j/h                                                                 
        P[i+1,i] = 0.5e0j/h                                                                  
                                                                                              
    P[0,n-1] = 0.5e0j/h                                                                      
    P[n-1,0] = -0.5e0j/h                                                                     
                                                                                              
    # Solve                                                                                  
    #                                                                                        
    L, V = np.linalg.eig(P)                                                                  
                                                                                              
    # Sort by eigenvalues
    #                                                                    
    _map = np.argsort(L)                                                                     
    L = L[_map]                                                                              
    V = V[:,_map]                                                                            
                                                                                              
    # Print eigenvalues                                                                      
    for i, l in enumerate(L):                                                                
        print(f"{i} {l.real}")                                                               
                                                                                              
                                                                                              
if __name__ == '__main__':                                                                   
    main()

 A: This is a routine application of Sylvester's shift matrix (1882), applied by Weyl (1927) to his celebrated quantum mechanics around the clock.
Without loss of generality, take n=5 and see how the clockwork ticks. (You should be able to generalize to arbitrary n trivially).
Define $\omega=\exp(i2\pi/5)$, so $1+\omega+ \omega^2+ \omega^3+ \omega^4=0$.
The shift matrix is
$$
S= 
\begin{bmatrix}0 &  1&0& 0&  0\\
 0& 0 & 1&0&0\\
0& 0&0&0& 1\\
 1& 0&0& 0&0\end{bmatrix},$$
with an evident cycle of 5, $S^5=I$.
Its eigenvalues are, of course, $ \omega, \omega^2,  \omega^3, \omega^4, 1$ summing to zero, evident from above, as well as the trace of the matrix. The corresponding eigenvectors are $v_1=(1,\omega,  \omega^2, \omega^3, \omega^4)$,  $v_2=(1,\omega^2,  \omega^4, \omega,\omega^3)$, $v_3=(1,\omega^3,  \omega,\omega^4, \omega^2)$, $v_4=(1,\omega^4,  \omega^3, \omega^2, \omega)$, $v_5=(1, 1,1,  1,1)$.  Observe how complex conjugation pairs them up, $v_1^*= v_4$, $v_2^*=  v_3$, as well as their eigenvalues, (with the minus sign to follow below further flipping them).
You then have
$$hP={1\over 2i} ( S-S^\dagger),$$
with eigenvectors $v_5;v_1,v_4;v_2,v_3,$ for eigenvalues 0; $\pm\sin (2\pi/5)$; $\pm\sin (4\pi/5)$. Plug in and check.
The pattern $\pm\sin (2k\pi/n)$ for generic n should be evident. The pairing of the eigenvalues is a straightforward consequence of the hermiticity with the antisymmetry/complex conjugation symmetry you noted.
QM around the clock is a highly recommended area you might wish to master.
A: You have discovered the notorious fermion doubling problem.  There is no way to discretize the first order momentum operator and at the same time preserve hermiticity and a low energy spectrum that looks
like $p$.
