Square-lattice topo. ins. & edge states

(See arXiv:0802.3537 and eqs. (14, 15) therein.)

(*Pauli matrices*)
sig[0]={{1,0},{0,1}};
sig[1]={{0,1},{1,0}};
sig[2]={{0,-I},{I,0}};
sig[3]={{1,0},{0,-1}};
(*diagonal element*)
Hamdiag[mass_,ky_]:=Sin[ky] sig[2]+(mass+Cos[ky]) sig[3]
(*off-diagonal element*)
Hamoff[ky_]:=(sig[3]-I sig[1])/2
(*Hamiltonian for an Lx-layer system with a mass parameter "m"*)
OpenHam[mass_,ky_,Lx_]:=KroneckerProduct[SparseArray[{{i_,i_}->1},{Lx,Lx}],Hamdiag[mass,ky]]+KroneckerProduct[SparseArray[{{i_,j_}/;i-j==-1->1},{Lx,Lx}],Hamoff[ky]]+KroneckerProduct[SparseArray[{{i_,j_}/;i-j==1->1},{Lx,Lx}],ConjugateTranspose[Hamoff[ky]]]

mesh=30;
evtable[mass_,Lx_]:=Module[{tab,klist,pretab,fin,newtab},
tab[ky_]:=Module[{pre,dim},
pre=Eigenvalues[OpenHam[mass,ky,Lx]]//Sort;
dim=Length[pre];
Table[{ky,m,pre[[m]]},{m,1,dim}]
];
klist=Table[-π+2.0 π i/mesh,{i,1,mesh}];
pretab[1]:=tab[klist[[1]]];
pretab[k_]:=pretab[k]=Join[pretab[k-1],tab[klist[[k]]]];
fin=pretab[Length[klist]];
newtab[k_]:=Module[{pre},
pre=Select[fin,#[[2]] ==k &];
Table[{pre[[m]][[1]],pre[[m]][[3]]}//Flatten,{m,1,Length[pre]}]
];
Table[newtab[k],{k,1,2 Lx}]
]

ListPlot[evtable[-1.0,25],Joined-> True,PlotStyle->{Black}]

ListPlot[evtable[-2.0,25],Joined-> True,PlotStyle->{Black}]

ListPlot[evtable[-3.0,25],Joined-> True,PlotStyle->{Black}]