# 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.

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$.

• Thanks @MengCheng, that's exactly what I was looking for. – Aegon Nov 4 '15 at 19:12
• It will work and this is completely equivalent to the method described in your answer. The reason is that $h$ is not a random Hermitian matrix: it has particle-hole symmetry. – Meng Cheng Nov 4 '15 at 20:11
• Actually, once I have periodic boundary conditions in one direction, it's probably easiest to use Nambu spinors i.e., $\Psi = (....,c_{k,n},c_{-k,n}^{\dagger},....)$ and then numerically diagonalize the resulting matrix. Of course, I could then do another unitary transformation to bring it to the form @MengCheng suggested in terms of the quasi-particle creation-annihilation operators. – Aegon Nov 5 '15 at 2:04
• @Aegon With periodic boundary condition you will never have edge modes. – Meng Cheng Nov 5 '15 at 2:07
• @MengCheng It's a 2D system and I would only have periodic boundary conditions in one-direction and open-boundary conditions in the other. See arxiv.org/pdf/1105.4700.pdf for instance, where they explicitly demonstrate this in Fig. 27. I'm just trying to redo their calculation. – Aegon Nov 5 '15 at 2:10

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

• It seems the two methods you and MengCheng proposed are almost the same. Nevertheless, I'm wondering about the computational efficiency of yours. Is it simpler to find the Youla form than to diagonalise the Hamiltonian directly ? See en.wikipedia.org/wiki/Skew-symmetric_matrix#Spectral_theory for Youla's form (this is the matrix you gave in your answer) I mean in term of computational ressource. – FraSchelle Nov 5 '15 at 10:18
• I believe (but am not sure) the computational complexity is about the same. – Jahan Claes Nov 5 '15 at 15:26