フェルミループが存在する場合のHall係数の格子近似計算

(*Berry phase*)
Clear[m]
$Assumptions={kx,ky}∈Reals;
σ[i_]:=PauliMatrix[i];
L=100;
Eref=1;
h=Sin[kx] σ[1]+Sin[ky] σ[2]+(m-Cos[kx]-Cos[ky]) σ[3];
m=1.5;

n=3;
r[n_]:=RandomReal[{-1,1},{n,n}]+I*RandomReal[{-1,1},{n,n}];
h=Sum[r[n]*E^(I (nx kx+ny ky)),{nx,-1,1},{ny,-1,1}];
h=1/2*(h+h\[ConjugateTranspose])//Simplify;
h=h-1/n*Tr[h]*IdentityMatrix[n](*Trを差っ引く*);

Plot3D[{Eref+100,Eigenvalues[h+100*IdentityMatrix[n]]},{kx,-Pi,Pi},{ky,-Pi,Pi}]
hlis=Table[h/.{kx->nx*2 Pi/L,ky->ny*2 Pi/L},{nx,0,L},{ny,0,L}]//N;
(*Eval=Map[Re[Eigenvalues[#]]&,hlis,{2}];
nocc=Map[Length[Select[#,#<Eref&]]&,Eval,{2}]
*)


(*Eref以下の固有状態．プラケットで計算*)
(*
Us=Table[
U=Table[
{Eval,Evec}=Eigensystem[hlis[[nx,ny]]];
pos=Position[Eval,_?(Re[#]<Eref&)]//Flatten;
If[Length[pos]==0,{},Normalize/@Evec[[pos]]],
{nx,1,L+1},{ny,1,L+1}];
Do[U[[nx,L+1]]=U[[nx,1]],{nx,1,L+1}];
Do[U[[L+1,ny]]=U[[1,ny]],{ny,1,L+1}];
U,{Eref,-5,5,0.5}];
*)

γ=Table[
U=Table[
{Eval,Evec}=Eigensystem[hlis[[nx,ny]]];
pos=Position[Eval,_?(Re[#]<Eref&)]//Flatten;
If[Length[pos]==0,{},(Normalize/@Evec[[pos]])\[Transpose]]
,{nx,1,L+1},{ny,1,L+1}];
Do[U[[nx,L+1]]=U[[nx,1]],{nx,1,L+1}];
Do[U[[L+1,ny]]=U[[1,ny]],{ny,1,L+1}];
θs=Table[
u={U[[nx,ny]],U[[nx+1,ny]],U[[nx+1,ny+1]],U[[nx,ny+1]]};
If[FreeQ[u,{}],
W=u[[1]]\[ConjugateTranspose].u[[4]].u[[4]]\[ConjugateTranspose].u[[3]].u[[3]]\[ConjugateTranspose].u[[2]].u[[2]]\[ConjugateTranspose].u[[1]];
eval=Eigenvalues[W]//Chop;
eval=Select[eval,Abs[#]>10^-5&];
Arg[Product[z,{z,eval}]],
0]
,{nx,1,L},{ny,1,L}];
Total[θs//Flatten]/(2 Pi)
,{Eref,-5,5,0.01}]//Chop
ListPlot[γ,PlotRange->All]
(*
Do[Print[Map[Length,Us[[n]],{2}]//MatrixForm],{n,Length[Us]}]
*)

ListLinePlot[γ,PlotRange->All]

(*Berry phase. w/ Class AII TRS*)
Clear[m]
$Assumptions={kx,ky}∈Reals;
σ[i_]:=KroneckerProduct[PauliMatrix[0],PauliMatrix[i]];
L=20;
Eref=1;
h=Sin[kx] σ[1]+Sin[ky] σ[2]+(m-Cos[kx]-Cos[ky]) σ[3];
m=1.;

r[n_]:=RandomReal[{-1,1},{n,n}]+I*RandomReal[{-1,1},{n,n}];
n=2;
h=Sum[r[2*n]*E^(I (nx kx+ny ky)),{nx,-1,1},{ny,-1,1}];
h=h+h\[ConjugateTranspose]//Simplify;
h=h-1/(2*n)*Tr[h]*σ[0];
h=1/2*(h+σ[2].(h/.{kx->-kx,ky->-ky})\[Conjugate].σ[2])//Simplify;

Plot3D[{Eref+100,Eigenvalues[h+100*IdentityMatrix[2*n]]},{kx,-Pi,Pi},{ky,-Pi,Pi}]
hlis=Table[h/.{kx->nx*2 Pi/L-Pi,ky->ny*2 Pi/L-Pi},{nx,1,L+1},{ny,1,L+1}];
(*Eval=Map[Re[Eigenvalues[#]]&,hlis,{2}];
nocc=Map[Length[Select[#,#<Eref&]]&,Eval,{2}]
*)

(*Eref以下の固有状態．プラケットで計算*)
γ=Table[
U=Table[
{Eval,Evec}=Eigensystem[hlis[[nx,ny]]];
pos=Position[Eval,_?(Re[#]<Eref&)]//Flatten;
If[Length[pos]==0,{},(Normalize/@Evec[[pos]])\[Transpose]]
,{nx,1,L+1},{ny,1,L+1}];
Do[U[[nx,L+1]]=U[[nx,1]],{nx,1,L+1}];
Do[U[[L+1,ny]]=U[[1,ny]],{ny,1,L+1}];
θs=Table[
u={U[[nx,ny]],U[[nx+1,ny]],U[[nx+1,ny+1]],U[[nx,ny+1]]};
If[FreeQ[u,{}],
W=u[[1]]\[ConjugateTranspose].u[[4]].u[[4]]\[ConjugateTranspose].u[[3]].u[[3]]\[ConjugateTranspose].u[[2]].u[[2]]\[ConjugateTranspose].u[[1]];
eval=Eigenvalues[W]//Chop;
eval=Select[eval,Abs[#]>10^-5&];
Arg[Product[z,{z,eval}]],
0]
,{nx,1,L},{ny,1,L}];
Total[θs//Flatten]/(2 Pi)
,{Eref,-5,5,0.1}]
ListPlot[γ]