Graphene edge states

1. Zigzag

(*Nx: # of layers, k:momentum in the y-direction*)
Hzig[Nx_ ,k_ ]:=Module[{hoplist},
hoplist = Join[Table[{2Cos[k/2],1},{i,1,Nx-1}] // Flatten,{2Cos[k/2]}];
SparseArray[{Band[{1,2}]->hoplist,Band[{2,1}]->hoplist},{2Nx,2Nx}]
]

tab[Nx_,k_]:=Module[{pre,dim},
pre=Eigenvalues[Hzig[Nx,k]] // Sort;
dim = Length[pre];
Table[{k,m,pre[[m]]},{m,1,dim}]
]

Nx=20;
mesh = 40;

klist =  Table[2.0 π i/mesh,{i,0,mesh}];
pretab1=Flatten[Table[tab[Nx,klist[[i]]],{i,1,mesh}],1];
pretab2[m_]:=Select[pretab1,#[[2]] == m &];
newtab[m_]:=Table[{pretab2[m][[i]][[1]],pretab2[m][[i]][[3]]},{i,1,mesh}];
Engtab = Table[newtab[m],{m,1,2Nx}];
ListPlot[Engtab,Joined->True,PlotStyle->{Black}]

2. Armchair

(*Nx: # of layers, k:momentum in the y-direction*)
Harm[Nx_,k_]:=Module[{hoplist},
hoplist[1]=Join[Table[{Exp[I k/Sqrt[3]],Exp[I k/(2Sqrt[3])]},{i,1,Nx-1}] //Flatten,{Exp[I k/Sqrt[3]]}
];
hoplist[2]=Join[Table[{Exp[-I k/Sqrt[3]],Exp[-I k/(2Sqrt[3])]},{i,1,Nx-1}] //Flatten,{Exp[-I k/Sqrt[3]]}
];
hoplist[3]=Join[Table[{Exp[-I k/(2Sqrt[3])],0},{i,1,Nx-2}] //Flatten,{Exp[-I k/(2Sqrt[3])]}
];
hoplist[4]=Join[Table[{Exp[I k/(2Sqrt[3])],0},{i,1,Nx-2}] //Flatten,{Exp[I k/(2Sqrt[3])]}
];
SparseArray[{Band[{1,2}]->hoplist[1],Band[{2,1}]->hoplist[2],Band[{1,4}]->hoplist[3],Band[{4,1}]->hoplist[4]},{2Nx,2Nx}]
]

tab[Nx_,k_]:=Module[{pre,dim},
pre = Eigenvalues[Harm[Nx,k]]//Sort;
dim = Length[pre];
Table[{k,m,pre[[m]]},{m,1,dim}]
]
Nx = 20;
mesh =40;
klist =Table[-π /Sqrt[3]+ 2.0 π i/(Sqrt[3] mesh),{i,1,mesh}];
pretab1 =Flatten[ Table[tab[Nx,klist[[i]]],{i,1,mesh}],1];
pretab2[m_]:=Select[pretab1,#[[2]] == m &];
newtab[m_]:=Table[{pretab2[m][[i]][[1]],pretab2[m][[i]][[3]]},{i,1,mesh}];
Engtab = Table[newtab[m],{m,1,2Nx}];
ListPlot[Engtab,Joined->True,PlotStyle->{Black}]