Chern Insulator

(Press Alt+A and then Shift+Enter to start)

Pauli matrices

sigma1={{0,1},{1,0}};

sigma2={{0,-I},{I,0}};

sigma3={{1,0},{0,-1}};

Bloch Hamiltonian

ChernBloch[kx_,ky_,m_]:=Sin[kx]sigma1+Sin[ky]sigma2+(m-Cos[kx]-Cos[ky])sigma3

Dispersion relation (gap closing at m=-2,0,2)

Manipulate[Plot3D[Sort[Round[Eigenvalues[ChernBloch[kx,ky,m]],10^-3],Greater],{kx,-π,π},{ky,-π,π},PlotRange->{-5,5},ViewPoint->{2,1,0.2},ImageSize->600],{m,-3,3,0.1}]

Bloch Hamiltonian as map to (ℝ^3)

v[kx_,ky_,m_]:={Sin[kx],Sin[ky],m-Cos[kx]-Cos[ky]}

Visualising the image of the momentum torus

Module[{m=1.3},
Manipulate[
Show[
ParametricPlot3D[v[kx,ky,m],
{kx,-π,π},{ky,-π,max},
PlotRange->{{-2,2},{-2,2},{-5,5}},
PlotPoints->20,
Axes->False,
Boxed->False,
Mesh->False,
PlotStyle->{If[Round[m,0.01]==-2||Round[m,0.01]==0||Round[m,0.01]==2,Red,Gray],Opacity[1],FaceForm[Blue,Green]},
Lighting->"Neutral",
ImageSize->500,
ViewPoint->{1,0.5,0.5}
],
Graphics3D[{Black,Sphere[{0,0,0},.1]}],
Graphics3D[Line[{{-2,0,0},{2,0,0}}]],
Graphics3D[Line[{{0,-2,0},{0,2,0}}]],
Graphics3D[Line[{{0,0,-5},{0,0,5}}]]
],
{{max,-π+0.1},-π+0.1,π,0.02}
]
]

Varying the parameter m

Manipulate[
Show[
ParametricPlot3D[v[kx,ky,m],
{kx,-π,π},{ky,-π,π},
PlotRange->{{-5,5},{-5,5},{-5,5}},
PlotPoints->20,
Axes->False,
Boxed->False,
Mesh->False,
PlotStyle->{If[Round[m,0.01]==-2||Round[m,0.01]==0||Round[m,0.01]==2,Red,FaceForm[Blue,Green]],Opacity[0.8]},
Lighting->"Neutral",
ImageSize->700
],
Graphics3D[{Black,Sphere[{0,0,0},.1]}],
Graphics3D[Line[{{-5,0,0},{5,0,0}}]],
Graphics3D[Line[{{0,-5,0},{0,5,0}}]],
Graphics3D[Line[{{0,0,-5},{0,0,5}}]]
],
{{m,0},-2.5,2.5,0.1}]

Homotopy invariant

“Flattened” Hamiltonian as a map to (S^2) (equivalent to ground state map)

vNorm[kx_,ky_,m_]:=v[kx,ky,m]/(v[kx,ky,m].v[kx,ky,m])^0.5

Mapping degree (= Chern number)

ListPlot[
Table[
{m,1/(4π)*NIntegrate[vNorm[kx,ky,m].Cross[D[vNorm[kx,ky,m],kx],D[vNorm[kx,ky,m],ky]],{kx,-π,π},{ky,-π,π},AccuracyGoal->3]},
{m,-3,3,0.123}],
PlotStyle->PointSize[0.02]
]//Quiet

Boundary states

Defining the Hamiltonian on a strip with translation invariance in the y-direction

hopping[L_]:=SparseArray[Table[{i+1,i}->1,{i,1,L-1}],L,0]

ChernStrip[kx_,m_,L_]:=KroneckerProduct[hopping[L],1/2(I*sigma2-sigma3)]+
Conjugate[Transpose[KroneckerProduct[hopping[L],1/2(I*sigma2-sigma3)]]]+
KroneckerProduct[IdentityMatrix[L],(m-Cos[kx])sigma3+Sin[kx]sigma1]

Varying m

Manipulate[Module[{x=Table[{k,Sort[Round[Eigenvalues[ChernStrip[k,m,20]],0.001],Greater]},{k,-π,π,π/48.}]},ListLinePlot[Table[Table[{x[[i]][[1]],x[[i]][[2]][[j]]},{i,1,Length[x]}],{j,1,Length[x[[1]][[2]]]}],PlotStyle->{RGBColor[{0,0,0}]},ImageSize->500,PlotRange->{{-π,π},{-3,3}}]],
{m,-2.5,2.5,0.1}]

Inspection of zero energy states

MinVecs2 gives the two last eigenvectors (which have the smallest absolute value)

MinVecs2[matrix_]:=Module[{x=Eigensystem[matrix]},Take[x[[2]],{Length[matrix]-1,Length[matrix]}]]

Probability density of zero energy states in y-direction

Manipulate[ListLinePlot[Table[Abs[MinVecs2[ChernStrip[If[m≤0,π,0],m,20]]][[i]]^2,{i,1,2}],PlotRange->{0,1}],{{m,1.7},-2.,2.,0.1}]

Quantum Spin Hall Effect

PauliOld={sigma1,sigma2,sigma3};

sigma1new:=KroneckerProduct[sigma3,sigma1]

sigma2new:=KroneckerProduct[IdentityMatrix[2],sigma2]

sigma3new:=KroneckerProduct[IdentityMatrix[2],sigma3]

PauliNew={sigma1new,sigma2new,sigma3new};

4-by-4 Hamiltonian by replacing Pauli matrices:

QSH[kx_,ky_,m_]:=vNorm[kx,ky,m].PauliNew

Chern number for m=1.3:

Module[{m=1.3},
Round[
1/(16π I)*NIntegrate[Tr[QSH[kx,ky,m].(D[QSH[kx,ky,m],ky].D[QSH[kx,ky,m],kx]-D[QSH[kx,ky,m],kx].D[QSH[kx,ky,m],ky])],
{kx,-π,π},{ky,-π,π},
AccuracyGoal->3],
10^-5]]

Boundary states

Defining the Hamiltonian on a strip with translation invariance in the y-direction

hopping[L_]:=SparseArray[Table[{i+1,i}->1,{i,1,L-1}],L,0]

QSHStrip[kx_,m_,L_]:=KroneckerProduct[hopping[L],1/2(I*sigma2new-sigma3new)]+
Conjugate[Transpose[KroneckerProduct[hopping[L],1/2(I*sigma2new-sigma3new)]]]+
KroneckerProduct[IdentityMatrix[L],(m-Cos[kx])sigma3new+Sin[kx]sigma1new]

Varying m

Manipulate[Module[{x=Table[{k,Sort[Round[Eigenvalues[QSHStrip[k,m,20]],0.001],Greater]},{k,-π,π,π/48.}]},ListLinePlot[Table[Table[{x[[i]][[1]],x[[i]][[2]][[j]]},{i,1,Length[x]}],{j,1,Length[x[[1]][[2]]]}],PlotStyle->{RGBColor[{0,0,0}]},ImageSize->500,PlotRange->{{-π,π},{-3,3}}]],
{m,-2.5,2.5,0.1}]

Inspection of zero energy states

MinVecs4 gives the four last eigenvectors (which have the smallest absolute value)

MinVecs4[matrix_]:=Module[{x=Eigensystem[matrix]},
Take[x[[2]],{Length[matrix]-3,Length[matrix]}]]

Probability density of zero energy states in y-direction

Manipulate[ListLinePlot[Table[Abs[MinVecs4[QSHStrip[If[m≤0,π,0],m,20]]][[i]]^2,{i,1,4}],PlotRange->{0,1}],{{m,1.7},-2.,2.,0.1}]

Breaking time reversal with strength B

sigmaB:=KroneckerProduct[sigma1,{{1,0},{0,0}}]

(Only last summand B·sigmaB added)

QSHStripB[kx_,m_,L_,B_]:=KroneckerProduct[hopping[L],1/2(I*sigma2new-sigma3new)]+
Conjugate[Transpose[KroneckerProduct[hopping[L],1/2(I*sigma2new-sigma3new)]]]+
KroneckerProduct[IdentityMatrix[L],(m-Cos[kx])sigma3new+Sin[kx]sigma1new+B sigmaB]

Varying B

Manipulate[Module[{x=Table[{k,Sort[Round[Eigenvalues[QSHStripB[k,1.3,20,B]],0.001],Greater]},{k,-π,π,π/48.}]},ListLinePlot[Table[Table[{x[[i]][[1]],x[[i]][[2]][[j]]},{i,1,Length[x]}],{j,1,Length[x[[1]][[2]]]}],PlotStyle->{RGBColor[{0,0,0}]},ImageSize->500,PlotRange->{{-π,π},{-3,3}}]],
{B,0,1,0.05}]