# Entanglement spectrum of a superconductor

Consider the Hamiltonian $$H = \sum_{k} \epsilon_k c_{k}^{\dagger} c_{k} + (\Delta_k c_{k}c_{-k} + h.c.)$$ where $$\epsilon_k$$ and $$\Delta_k$$ are that of p-wave superconductor. So we have essentially a p-wave superconductor.

I want to compute the real space entanglement spectrum. So I divide the system into two parts A and B, and decompose all the operators in two parts. Then using these decomposed operators I write the Hamiltonian in the following basis $$\begin{pmatrix} c_{k,A}^{\dagger} & c_{k,B}^{\dagger} & c_{-k,A} & c_{-k,B} \end{pmatrix}$$. Then I use the correlation-matrix method described in arXiv:cond-mat/0212631 to find the entanglement spectrum. I compute the following correlation matrix $$\begin{pmatrix} \langle c_{k,A}^{\dagger} c_{k,A} \rangle & \langle c_{-k,A} c_{k,A} \rangle \\ \langle c_{k,A}^{\dagger} c_{-k,A}^{\dagger} \rangle & \langle c_{-k,A} c_{-k,A}^{\dagger} \rangle \end{pmatrix}$$ The python function for computing the correlation matrix is

import numpy as np
from mpmath import *
from Hamiltonian import Hamiltonian

def correlation_matrix(ky,n_cutoff,CA,CB):
H = Hamiltonian(ky,n_cutoff,CA,CB)
eigvals, eigvecs = eigh(H)
U = eigvecs.transpose()
U = U**-1
N = int(len(U)/2)
C = matrix(N,N)
for i in range(N):
for j in range(N):
for k in range(2*N):
C[i,j] += conj(U[2*i,k]) * U[2*j,k] * np.heaviside(-np.float(eigvals[k]),0)
return C


I do this entire calculation on the surface of a sphere so I expect to get a chiral Majorana mode in my entanglement spectrum. However I end up getting two counter-propagating Majorana modes. Can anybody help me spot the problem? Or guide me in computing the entanglement spectrum of superconductors!

• Isn't there a contradiction in your use of the wavenumber $k$, which is non-local in space, with a decomposition into two spatially distinct parts? Aug 27 at 13:54
• You can define an operator $c_{k,A}^{\dagger}$ which creates an electron in the subspace A with momentum k. So the decomposition of $c_k^{\dagger}$ would be $c_k^{\dagger} = \alpha_k c_{k,A}^{\dagger} + \beta_k c_{k,B}^{\dagger}$ for some coefficients $\alpha_k$ and $\beta_k$. It is like writing a wavefunction as a sum of two wavefunctions, each of which is non-zero only in one of the subspaces. Aug 27 at 14:27