Exact Diagonalization of a BdG Hamiltonian on a Finite Lattice I would like to numerically find the edge modes of a $p_x$ + $i p_y$ BdG Hamiltonian. The lattice version is given by
H = $\sum\left[-t \left(c_{m+1,n}^{\dagger} c_{m,n} + \text{h.c} \right) - t\left(c_{m,n+1}^{\dagger} c_{m,n} + \text{h.c} \right) - \mu\,c_{m,n}^{\dagger}c_{m,n} + \left(\Delta c_{m+1,n}^{\dagger} c_{m,n}^{\dagger} + \Delta^* c_{m,n} c_{m+1,n}\right) + \left(i\Delta c_{m,n+1}^{\dagger} c_{m,n}^{\dagger} -i \Delta^* c_{m,n} c_{m,n+1}\right)\right]$
where $c_{m,n}$ is the annihilation operator for a spin polarised fermion on site (m,n).
While I understand how to take this system and put it on a finite lattice if there are only hopping terms ($c^{\dagger}c$ terms), how would I do that for terms such as $c^{\dagger}c^{\dagger}$? 
Specifically, I want to find the spectrum of this system if it has periodic boundary conditions in one direction and open boundary conditions in the other.
 A: So: I assume you want to diagonalize this problem by rewriting the Hamiltonian as $H=\sum E_id_i^\dagger d_i$, where $d_i$ are quasiparticle operators which obey the Fermionic commutation relations.
If we only had $c^\dagger c$ terms, we would be able to write H as
$$
H=H_{ij}c_i^\dagger c_j
$$
We could then prove that if $\{c_i\}$ obey the Fermion commutation relations, then so too do $\{U_{ij}c_j\}$, where $U$ is any unitary matrix. We would then just diagonalize $H_{ij}$ with unitary matrices, and get our quasiparticle operator.
With the $\Delta$ term, we can no longer do this, because our quasiparticle operators will in general have to combine $c_i$ and $c_i^\dagger$. However, there is a clever way to get around this.
Define majorana operators $\gamma_{2j-1}=c_j+c_j^\dagger$, $\gamma_{2j}=\frac{c_j-c_j^\dagger}{i}$. A simple calculation shows that $\{\gamma_a, \gamma_b\}=2\delta_{ab}$, and $\gamma^\dagger=\gamma$. You can then rewrite your Hamiltonian as $H=\sum H_{ij} \gamma_i\gamma_j$.
A simple calculation shows that if $O$ is any orthogonal matrix, then $\{O_{ij}\gamma_j\}$ also obey the same commutation relations. So, you can freely diagonalize $H_{ij}$ using orthogonal matrices, get a set of $\tilde{\gamma}_i$s, then transform these back into fermionic operators. 
As a side note, you don't quite want to get $H_{ij}$ into a diagonal form using the orthogonal operators. You actually want to get it in the form
$$
\left(\begin{array}{cccc}
0&e_i&...&0 \\
-e_i&0&...&0 \\
...&...&...&... \\
0&0&...&0 \\
\end{array}\right)
$$
as this will ensure that H is diagonal when you transform back into the fermion operators.
I've just given a rough sketch, it's a bit involved but it will work. For more, see here: http://arxiv.org/pdf/cond-mat/0010440v2.pdf
A: First you need to bring it into the following form:
$H=\Psi^\dagger h \Psi$
Here $\Psi$ is a big column vector:
$\Psi=(\dots, c_{m,n}, \dots, c_{m,n}^\dagger, \dots)^T$
Basically, the first half of $\Psi$ are all annihilation operators, and the second half are all creation ones. If the number of sites is $N$, the size of $\Psi$ is $2N$. So $h$ is a $2N\times 2N$ matrix. To bring it into this form, one has to do a little bit of work, to rewrite all $c_i^\dagger c_j$ term as $-c_j c_i^\dagger$, etc. But this is not too difficult.
If you do everything correctly, $h$ is a Hermitian matrix and you can now go and diagonalize it. The results are of course the energies of the Bogoliubov quasiparticles and their forms are given by the unitary transformation that diagonalize $h$.
