A basic introduction to Majorana Edge modes in a Kitaev Chain

Welcome to another blog! This blog will give a very brief introduction and background of the Majorana zero modes and way to obtain them in a simple model known as the Kitaev Chain. The aim will be to translate the Kitaev Chain Hamiltonian into a Matrix form to obtain energy spectrum and edge modes for an open chain. We will obtain these Majorana zero modes at the edges of an open chain.

Majorana Fermions

In the year 1937, a new class of particles that are its own anti-particles were hypothesized by Ettore Majorana. This new class were termed Majorana fermions or particles for obvious reasons. They can be mathematically represented by creation and anhilation operators in second quantization:

$$\gamma^\dagger_k = \gamma_k$$ where, $\gamma_k$: anhilates a particle in quantum state k; $\gamma_k^\dagger$: creates a particle in quantum state k

The equation above conveys that the operators for creating and destroying the particles are same, which is the desired property of the Majorana in contrast to Dirac fermions for which they are distinct. The above operators shall be referred as Majorana operators. There are several conventions to define these Majorana operators. We will be adopting the following($f$ is the fermionic operator):

\begin{align} f^\dagger = (\gamma_1 + i\gamma_2)/2 \nonumber \qquad f = (\gamma_1 - i\gamma_2)/2 \label{2} \end{align}

Using this definition, the following identities can be checked

1. $\gamma_k^2 = 1$
2. $\left[\gamma_i,\gamma_j\right] = 2\delta_{ij}$

There have been lot of attempts in experimental high energy physics to hunt these particles.

Majorana Zero Modes

Let’s shift gears and enter the realm of condensed matter physics. In superconducting materials, certain quasiparticle(Bogoliubov) excitations display similar properties as that of Majorana fermions. The Majorana zero mode(MZM) is a special case of Bogoliubov quasiparticle operators which are linear combinations of electron and hole operator. The MZM satisfy the properties that they square to identity and are self-Hermitian.

If we look carefully at above equations, the two Majorana modes always occurs in pair for a single fermionic mode. Can we decouple them? This was cleverly answered by Professor Alexei Kitaev in 2001, through his Kitaev chain, which is a 1D p-wave superconducting model.

Kitaev Chain

Kitaev Chain is a tight binding lattice model of spinless electrons with nearest-neighbour hopping amplitude $t$, a p-wave superconducting pairing term $\Delta$ for nearest neighbours and on-site chemical potential $\mu$. The hamiltonian for a finite and open chain in second quantization can be jotted down \begin{equation} \hat{H} = \sum_{n=1}^{N-1}\left[t\left(f_n^\dagger f_{n+1} + f_{n+1}^\dagger f_n\right) + \Delta\left(f_n f_{n+1} + f_{n+1}^\dagger f_n^\dagger\right)\right] - \sum_{n=1}^{N} \mu \left( 2f_n^\dagger f_n - 1 \right) \label{real_space_Hamiltonian} \end{equation} where $t,\Delta$ and $\mu$ are all real and $\gamma > 0$ without loss of generality. Similar to the definition above, we define the following Majorana operators \begin{align} a_{2n-1} = f_n + f_n^\dagger \quad \text{and} \quad a_{2n} = i(f_n - f_n^\dagger) \label{4} \end{align} for n = 1,2,…,N. The Kitaev Chain Hamiltonian can be rewritten in terms of these Majorana operators \begin{equation} \hat{H} = i\sum_{n=1}^{N-1} \left[J_xa_{2n}a_{2n+1} - J_ya_{2n-1} a_{2n+2}\right] + i\sum_{n=1}^{N}\mu a_{2n-1}a_{2n} \end{equation} where, $J_x = \frac{1}{2}(t - \Delta) \quad \text{and} \quad J_y = \frac{1}{2}(t + \Delta)$

Let us carefully examine the Hamiltonian in the Majorana basis. A brief look will help us solve the problem of decoupling. Let us look at two special cases in the parameter space of $t, \Delta$ and $\mu$.

Case I: $t = \Delta$ and $\mu = 0$

The conditions clearly imply $J_x=0$ and $J_y=t$ which results in \begin{equation} \hat{H} = -i\sum_{n=1}^{N-1} \left[t a_{2n-1} a_{2n+2}\right] \label{6} \end{equation}

Did we end up finding something interesting? A close attention to the equation above reveals that the even(odd) Majorana mode at the left(right) end has vanished and it is decoupled from the rest of the modes. This is exactly what we set out for! The unpaired Majorana modes.

Case II: $t = -\Delta$ and $\mu = 0$

This time $J_y = 0$ and $J_x = t$ which results in \begin{equation} \hat{H} = i\sum_{n=1}^{N-1} \left[t a_{2n}a_{2n+1}\right] \end{equation}

This time, the odd(even) Majorana mode at the left(right) end of the chain is unpaired with the rest of the chain.

Wonderful! We know a way to obtain the unpaired Majorana end modes. Let us test this out through a Python Snippet!

Let us import the necessary Python modules

import numpy as np			# Numerical computation
import numpy.linalg as la		# Linear Algebra 
import matplotlib as mp			# Generating plots
import matplotlib.pyplot as plt
plt.style.use('seaborn')		# Setting the plotting style
mp.rcParams['figure.figsize'] = (9, 7)  # Setting the size of the plots
%matplotlib notebook			# Inline plotting

Let us define the various parameters of the system

Nsites = 25              # Number of lattice sites
Nprime = 2*Nsites
e_threshold = 1E-6	 # Threshold for finding zero eigenstates
params = {
't' : 2.0,               # Nearest neighbor hopping
'Delta' : 2.0,           # Superconducting pairing term
'mu' : 0.0               # Chemical potential
}

We define the Kitaev Chain Hamiltonian in the Majorana basis

def kitaev_ham(Nsites,params):
    Hmat = np.zeros([Nprime,Nprime])		# Declare a 2Nx2N matrix
    Jx = 0.5*(params['t'] - params['Delta'])	
    Jy = 0.5*(params['t'] + params['Delta'])
    
    for n in range(Nsites-1):
        Hmat[2*n,2*n+1]   = Jx
        Hmat[2*n+1,2*n]   = -Jx
        Hmat[2*n-1,2*n+2] = -Jy
        Hmat[2*n+2,2*n-1] = Jy
        Hmat[2*n-1,2*n]   = params['mu']
        Hmat[2*n,2*n-1]   = -params['mu']

    Hmat[2*(Nsites-1)-1,2*(Nsites-1)] = params['mu']
    Hmat[2*(Nsites-1),2*(Nsites-1)-1] = -params['mu']
    Hmat = 1j*Hmat

    return Hmat

Now that we have defined our Hamiltonian in a Matrix representation, let us visualise it.

Visualize the Hamiltonian Matrix

def visualise_mat(Hmat):
    plt.imshow(Hmat.imag) 	# The real part of the matrix is a zero matrix
    plt.colorbar()
    plt.show()

We have fixed the system parameters such that we realise case I. Let us see if we get any zero energy modes.

Diagonalise the Hamiltonian matrix and plot the spectrum

def plot_spectrum(Hmat):
    evals,evecs = la.eigh(Hmat)
    evals = evals.real
    plt.scatter(np.arange(len(evals)),evals)
    plt.title('Energy Spectrum of Chain with {} Sites'.format(Nsites))
    plt.show()
    
plot_spectrum(kitaev_ham(Nsites,params))

Hurray! We do see two energy modes at zero energy. Rest of the spectrum lies above and below the zero energy and also well separated. We can take it further and verify if the eigenvectors corresponding to these zero energy modes are localised at the edges of the chain.

Check for zero energy modes

# Extract the indices of energy modes close to zero
def check_zeromodes(evals):
    nzmodes = 0
    zmodes_ind = np.where(abs(evals) <= e_threshold)[0]
    return zmodes_ind,len(zmodes_ind)

Plot the probability distribution of the Majorana zero modes(edge modes)

def plot_zeromodes(evals,evecs,params):

    param_info = '\n'.join((
    r'$\mu=%.2f$' % (params['mu']),
    r'$t=%.2f$' % (params['t']),
    r'$\Delta=%.2f$' % (params['Delta'])))

    zmodes_ind,cnt_zmodes = check_zeromodes(evals)
    if cnt_zmodes > 0:
        fig,ax = plt.subplots(1,cnt_zmodes,figsize=(20, 10))
        fig.suptitle('Probability distribution of Zero modes',fontsize=20, fontweight='bold')
        for cnt in range(cnt_zmodes):
            ax1 = ax[cnt]
            ax1.plot(np.abs(evecs[:,zmodes_ind[cnt]])**2)
            ax1.set_title('Edge mode {}'.format(cnt+1),fontsize=20)
            ax1.set_xlabel('Site Number',fontsize=20)
            ax1.set_ylabel('$|\psi|^2$',fontsize=20)
            ax1.text(0.43, 0.95, param_info, transform=ax1.transAxes, fontsize=16,
        verticalalignment='top', bbox=dict(boxstyle="square",facecolor="white"))
            ax1.tick_params(axis='both', which='major', labelsize=16)
        #plt.savefig('Edge_modes_Kitaev.pdf')
        plt.show()
        
evals,evecs = la.eigh((kitaev_ham(Nsites,params)))
plot_zeromodes(evals,evecs,params)

Indeed, we see that wavefunction($|\psi|^2$ to be precise) corresponding to these zero energy modes are localised at the edges of the chain. We can also repeat the same exercise for Case II.

Wonderful! We know a way to obtain the unpaired Majorana end modes. The Majorana modes in the each case have zero energy and are localised at the end of the chain. Thus, it makes sense to coin the term, Majorana zero modes(MZMs). It is fascinating that the MZMs appear when tweaking the parameters! However, there might be some skepticism here as we obtained MZMs at a very special points in the parameter space. Let us vary the parameters of our system and check the fate of these zero energy modes.

var_mu = np.linspace(0,4,101)
var_energy = np.zeros([len(var_mu),Nprime])

for i in range(len(var_mu)):
    var_energy[i] = la.eigh(kitaev_ham(Nsites,params = {'t' : 2.0,'Delta' : 2.0, 'mu' : var_mu[i]}))[0]

plt.title("Energy Spectrum as a function of $\mu/t$ (N = 15)")
for i in range(Nprime):
    plt.plot(var_mu,var_energy[:,i])
plt.ylabel('Energy')
plt.xlabel('$\mu/t$')
plt.show()

It is a very interesting result. The zero energy modes exist even when we tweak the chemical potential. In other words, the MZMs don’t just appear by just fine tuning to a special parameter values but isolated MZMs stay protected as long as the bulk gap is finite. The reason behind this exciting observation is due to a concept known as Symmetry Protected Topological Phases(SPTP) which is a story for some other blog. For more details on this, refer to the following online course - Topology in Condensed Matter.

Thank you for your time! Happy learning!