# Help with diagonalizing a BdG Hamiltonian

I am working on a simulation which diagonalizes a BdG Hamiltonian; let's start with one of the form: $$\hat{H}=-t\sum_{,\sigma}\hat{c}^\dagger_{r\sigma}\hat{c}_{r'\sigma}+\sum_{r}\Big(\Delta_o c^\dagger_{r+}c^\dagger_{r-}+h.c.\Big)$$ So a uniform s-wave superconductor with a nearest neighbor hopping term.

This is very easy to diagonalize analytically with a Fourier transform. If we say we are working with a 2D NxN square lattice with PBCs, we can easily see that we have $$\vec{k}$$ restricted to the first Brillouin zone, and the energies are of the form (I set lattice constant to 1 below): $$E(\vec{k})=\pm\sqrt{4t^2(cos(k_x)+cos(k_y))^2+|\Delta_o|^2}$$

I wrote a computer program which looks at the real space Hamiltonian; we think of $$\hat{H}$$ in the following form (dropping hats on creation/annihilation operators): $$\hat{H}=\begin{pmatrix}\vec{c}_+^\dagger\vec{c}_-^\dagger \vec{c}_+\vec{c}_-\end{pmatrix}H\begin{pmatrix}\vec{c}_+\\\vec{c}_-\\\vec{c}_+^\dagger\\\vec{c}_-^\dagger \end{pmatrix}$$ The vectors stand for $$\vec{c}_+$$ a column vector with the $$N^2$$ spin up annihilation operators for each lattice site etc.

We just choose a H such that the product equals the sum above, then use Matlab, in my case, to diagonalize the Hamiltonian. The eigenvalues should, and for the most part do, correspond to the analytic solution, but there are a few main issues.

1.)There seems to be a 4-fold degeneracy in my computer program, which is too much. Are there any reasons that my solutions are not independent?

2.)The solutions of the analytic solution are always in the list of eigen values I find from my code, but the reverse is not true: I get eigen values from my code which are not present in my analytic solution. I sum over all combination of $$k_x$$'s and $$k_y$$'s and still do not get all the same answers. To the very least does my process sound ok? If so, any ideas why I get extra values in one case but not the other.

*I set $$k_i\in[-\pi,\pi-\Delta k]$$ where $$\Delta k=\frac{2\pi}{N}$$. Think I have this right too...

Any help will be greatly appreciated!

No answers a few years later but I did solve this a while back. The issue was I artificially doubled my solution set. This problem is much simpiler to solve with the following Nambu spinor: $$\psi=\begin{bmatrix} \vec{c}_\uparrow \\\vec{c}\;^\dagger_\downarrow \end{bmatrix}$$