第68回物性若手夏の学校
「結晶対称性とトポロジカル絶縁体」 
補足資料

質問等は
ken.shiozaki@yukawa.kyoto - u.ac.jp
にご連絡ください。

各種関数の定義

Pfaffianの実装については以下をコピーした。
https : // en.wikipedia.org/wiki/Talk : Pfaffian

(*ランダムにn×n複素行列を生成する*)
rmat[n_]:=RandomReal[{-1,1},{n,n}]+I*RandomReal[{-1,1},{n,n}];
(*Pfaffian*)
Pf[A_] := 
  Switch[Length[A], 0, 1, _?OddQ, 0, _?EvenQ, xPf[A, 1]];
xPf[A_, p0_] := Module[{A0, n, pivot, sign = 1, A1, p1},
   n = Length[A]/2;
   If[n != 1, A0 = A;
    pivot = First[Ordering[Normal[Abs[A0[[2 n - 1, All]]]], -1]];
    If[pivot != 2 n, A0[[{pivot, 2 n}, All]] = A0[[{2 n, pivot}, All]];
     A0[[All, {pivot, 2 n}]] = A0[[All, {2 n, pivot}]];
     sign = -1;];
    p1 = A0[[2 n - 1, 2 n]];
    A1 = p1 A0[[1 ;; 2 n - 2, 1 ;; 2 n - 2]];
    A1 += (# - Transpose[#]) &@
     Outer[Times, A0[[1 ;; 2 n - 2, 2 n]], 
      A0[[1 ;; 2 n - 2, 2 n - 1]]];
    A1 /= p0;
    sign xPf[A1, p1], A[[1, 2]]]];

空間１次元、Wanner関数

ハミルトニアンをランダムに生成

(*Hamiltonian*)
Clear[kx,x,hop];
$Assumptions={kx,x}∈Reals;
n(*行列サイズ*)=4;
range(*ホッピングのレンジ*)=3;
ξ(*ホッピングの減衰スケール*)=2;
Do[hop[nx]=rmat[n]*E^(-Abs[nx]/ξ),{nx,0,range}];
hop[0]=1/2*(hop[0]+hop[0]\[ConjugateTranspose]);
Do[hop[-nx]=hop[nx]\[ConjugateTranspose],{nx,1,range}];
(*Do[hop[nx]=hop[nx]-1/n*Tr[hop[nx]],{nx,-range,range}];*)
h=Sum[hop[nx]*E^(-I*nx*kx),{nx,-range,range}]//Chop;
h=h+100*IdentityMatrix[n](*定数シフト*);
(*h=1/2*(h+h\[ConjugateTranspose])//Simplify//Chop(*エルミート化*);*)
Plot[Eigenvalues[h]-100,{kx,-Pi,Pi},Ticks->{{-Pi,0,Pi},Automatic},AxesLabel->{"(Subscript[k, x])","E((Subscript[k, x]))"}]
(*ListPlot[Table[{ikx,e},{e,Eigenvalues[h/.{kx->ikx}]},{ikx,-Pi,Pi,2 Pi/100}],Ticks->{Automatic,Automatic},AxesLabel->{"(Subscript[k, x])","E((Subscript[k, x]))"}]*)

レベル反発により縮退点は存在しないことに注意。

Wannier関数の計算

x = 0 に局在するWannier関数を構成する。

数値的に得られたセル周期関数をそのままフーリエ変換する場合

Eigensystemの仕様に注意：絶対値の大きな固有値から順番に並べられる。

Mx(*波数空間の格子点の数*)=1000;
U(*最低エネルギーのバンドのみ取り出す*)=Table[Normalize/@Eigensystem[N[h/.{kx->ikx}]][[2,-1;;]]\[Transpose],{ikx,0,2*Pi,2*Pi/Mx}];
U[[-1]]=U[[1]];
L=50;
W=Table[1/Mx*Sum[U[[j]]*E^(2*Pi*I*(j-1)/Mx*x),{j,1,Mx+1}],{x,-L,L}];
ListPlot[Table[{x,(Norm/@W)[[x+L+1]]},{x,-L,L}],AxesLabel->{"x","|(Subscript[W, 0])(x)|"},PlotRange->All,Filling->Axis,FillingStyle->Thickness[0.005]]
ListPlot[Table[{x,(Log/@Norm/@W)[[x+L+1]]},{x,-L,L}],AxesLabel->{"x","log|(Subscript[W, 0])(x)|"},PlotRange->All]

隣り合った波数における内積 <u(k+δk)|u(k)>  のU(1)位相を一定にする場合。

(*U(1)位相を<u|u+δu> =E^(i γ/δk)なるように修正*)
γ(*Berry位相を計算*)=Arg[Product[Det[U[[j+1]]\[ConjugateTranspose].U[[j]]],{j,1,Mx}]];
Print["Berry位相 = ",γ];
Do[
θ=Arg[Det[U[[j+1]]\[ConjugateTranspose].U[[j]]]];
U[[j+1]]=U[[j+1]]*E^(I*θ)*E^(-I*γ/Mx);
,{j,1,Mx}];
U[[-1]]=U[[1]];
W=Table[1/Mx*Sum[U[[j]]*E^(2*Pi*I*(j-1)/Mx*x),{j,1,Mx+1}],{x,-L,L}];
(*BarChart[Norm/@W,ChartLabels->Range[-L,L],AxesLabel->{"x","|(Subscript[W, 0])(x)|"}]*)
ListPlot[Table[{x,(Norm/@W)[[x+L+1]]},{x,-L,L}],AxesLabel->{"x","|(Subscript[W, 0])(x)|"},PlotRange->All,Filling->Axis,FillingStyle->Thickness[0.005]]
ListPlot[Table[{x,(Log/@Norm/@W)[[x+L+1]]},{x,-L,L}],AxesLabel->{"x","log|(Subscript[W, 0])(x)|"},PlotRange->All]

指数関数的局在が確認できる。
20 サイト程度で数値誤差に達していると思われる。

開放端条件におけるエネルギー固有値

-> ingap stateが出現することがある。

Lx(*サイト数*)=100;
hOBC=Table[0,{x1,Lx},{x2,Lx}];
Do[hOBC[[x+nx,x]]=hop[nx],{nx,0,range},{x,Lx-nx}];
Do[hOBC[[x,x+nx]]=hop[-nx],{nx,1,range},{x,Lx-nx}];
hOBC=ArrayFlatten[hOBC];
ListPlot[Eigenvalues[hOBC+100*IdentityMatrix[n*Lx]]-100]

(T^2)=-1 なる時間反転対称性を課した場合

-> 時間反転対称点 (Subscript[k, x])=0,π  において縮退（Kramers縮退）する。

(*Hamiltonian*)
Clear[kx,x,hop];
$Assumptions={kx,x}∈Reals;
matsize(*行列サイズ*)=4;
range(*ホッピングのレンジ*)=3;
ξ(*ホッピングの減衰スケール*)=2;
Do[hop[nx]=rmat[matsize]*E^(-Abs[nx]/ξ),{nx,0,range}];
hop[0]=1/2*(hop[0]+hop[0]\[ConjugateTranspose]);
Do[hop[-nx]=hop[nx]\[ConjugateTranspose],{nx,1,range}];
(*Do[hop[nx]=hop[nx]-1/n*Tr[hop[nx]],{nx,-range,range}];*)
h=Sum[hop[nx]*E^(-I*nx*kx),{nx,-range,range}]//Chop;
UT(*TRS*)=KroneckerProduct[IdentityMatrix[matsize/2],{{0,1},{-1,0}}];
h=1/2*(h+UT.(h/.{kx->-kx})\[Conjugate].UT\[ConjugateTranspose])//Simplify;
h=h+100*IdentityMatrix[matsize](*定数シフト*);
(*h=1/2*(h+h\[ConjugateTranspose])//Simplify//Chop(*エルミート化*);*)
Plot[Eigenvalues[h]-100,{kx,-Pi,Pi},Ticks->{{-Pi,0,Pi},Automatic},AxesLabel->{"(Subscript[k, x])","E((Subscript[k, x]))"}]
(*ListPlot[Table[{ikx,e},{e,Eigenvalues[h/.{kx->ikx}]},{ikx,-Pi,Pi,2 Pi/100}],Ticks->{Automatic,Automatic},AxesLabel->{"(Subscript[k, x])","E((Subscript[k, x]))"}]*)

Wannier関数

- x = 0 に局在するWannier関数。
 - 最低エネルギー固有状態のみを用いる。 -> 縮退点を含むことに注意する。

数値的に得られたセル周期関数をそのままフーリエ変換する場合

Mx(*波数空間の格子点の数*)=1000;
U(*最低エネルギーのバンドのみ取り出す*)=Table[Normalize/@Eigensystem[N[h/.{kx->ikx}]][[2,-1;;]]\[Transpose],{ikx,0,2*Pi,2*Pi/Mx}];
U[[-1]]=U[[1]];
L=50;
W=Table[1/Mx*Sum[U[[j]]*E^(2*Pi*I*(j-1)/Mx*x),{j,1,Mx+1}],{x,-L,L}];
ListPlot[Table[{x,(Norm/@W)[[x+L+1]]},{x,-L,L}],AxesLabel->{"x","|(Subscript[W, 0])(x)|"},PlotRange->All,Filling->Axis,FillingStyle->Thickness[0.005]]
ListPlot[Table[{x,(Log/@Norm/@W)[[x+L+1]]},{x,-L,L}],AxesLabel->{"x","log|(Subscript[W, 0])(x)|"},PlotRange->All]

隣り合った波数における内積 <u(k+δk)|u(k)>  のU(1)位相を一定にする場合。

(*U(1)位相を<u|u+δu> =E^(i γ/δk)なるように修正*)
γ(*Berry位相を計算*)=Arg[Product[Det[U[[j+1]]\[ConjugateTranspose].U[[j]]],{j,1,Mx}]];
Print["Berry位相 = ",γ];
Do[
θ=Arg[Det[U[[j+1]]\[ConjugateTranspose].U[[j]]]];
U[[j+1]]=U[[j+1]]*E^(I*θ)*E^(-I*γ/Mx);
,{j,1,Mx}];
U[[-1]]=U[[1]];
W=Table[1/Mx*Sum[U[[j]]*E^(2*Pi*I*(j-1)/Mx*x),{j,1,Mx+1}],{x,-L,L}];
ListPlot[Table[{x,(Norm/@W)[[x+L+1]]},{x,-L,L}],AxesLabel->{"x","|(Subscript[W, 0])(x)|"},PlotRange->All,Filling->Axis,FillingStyle->Thickness[0.005]]
ListPlot[Table[{x,(Log/@Norm/@W)[[x+L+1]]},{x,-L,L}],AxesLabel->{"x","log|(Subscript[W, 0])(x)|"},PlotRange->All]

-> バンド縮退のため、指数関数的に局在するWanner状態が得られない。

空間２次元

ランダムにハミルトニアンを生成する。

E = 0 において一様なエネルギーギャップが出現しやすいように、ハミルトニアンを手でトレースレスに再定義する。

ClearAll[kx,ky,x,y,ihop,hop];
$Assumptions={kx,ky,x,y}∈Reals;
matsize(*行列サイズ*)=2;
range(*ホッピングのレンジ*)=1;
ξ(*ホッピングの減衰スケール*)=3;
Do[ihop[nx,ny]=rmat[matsize]*E^(-Sqrt[nx^2+ny^2]/ξ),{nx,-range,range},{ny,-range,range}];
Do[hop[nx,ny]=1/2*(ihop[nx,ny]+ihop[-nx,-ny]\[ConjugateTranspose]),{nx,-range,range},{ny,-range,range}];
Do[hop[nx,ny]=hop[nx,ny]-1/matsize*Tr[hop[nx,ny]],{nx,-range,range},{ny,-range,range}];
h=Sum[hop[nx,ny]*E^(-I*(nx*kx+ny*ky)),{nx,-range,range},{ny,-range,range}]//Chop;
(*h=h-1/n*Tr[h]//Simplify;*)
h=h+100*IdentityMatrix[matsize](*定数シフト*);
(*h=1/2*(h+h\[ConjugateTranspose])//Simplify//Chop(*エルミート化*);*)
Plot3D[Eigenvalues[h]-100,{kx,-Pi,Pi},{ky,-Pi,Pi},Ticks->{Automatic,Automatic,Automatic},AxesLabel->{"(Subscript[k, x])","(Subscript[k, y])","E((Subscript[k, x]),(Subscript[k, y]))"}]

Chern数

(*Chern数*)
M(*波数空間の刻み数*)=30;
l(*下から何番目までのバンドのChern数を計算するか*)=1;
(*Berry flux over BZ*)
U=Table[Normalize/@Eigensystem[h/.{kx->ikx,ky->iky}][[2,-l;;]]\[Transpose],{ikx,-Pi,Pi,2*Pi/M},{iky,-Pi,Pi,2*Pi/M}];
SumFlux=Total[Table[Arg[Det[U[[i,j]]\[ConjugateTranspose].U[[i,j+1]].U[[i,j+1]]\[ConjugateTranspose].U[[i+1,j+1]].U[[i+1,j+1]]\[ConjugateTranspose].U[[i+1,j]].U[[i+1,j]]\[ConjugateTranspose].U[[i,j]]]],{i,1,M},{j,1,M}]
//Flatten];
SumFlux/(2*Pi)

y方向の分極 (Subscript[P, y])(kx) ∈ [0,1] の計算

(*Flow*)
Mx=100;
My=20;
(*Berry flux over BZ*)
U=Table[Normalize/@Eigensystem[h/.{kx->ikx,ky->iky}][[2,-l;;]]\[Transpose],{ikx,-Pi,Pi,2*Pi/Mx},{iky,-Pi,Pi,2*Pi/My}];
Do[U[[i,-1]]=U[[i,1]],{i,1,Mx+1}];
γy=Table[Arg[Product[Det[U[[i,j+1]]\[ConjugateTranspose].U[[i,j]]],{j,1,My}]],{i,1,Mx}];
γy=Table[{2*Pi*(nx-1)/Mx-Pi,γy[[nx]]/(2*Pi)},{nx,1,Mx}];
γy=AppendTo[γy,γy[[1]]];
ListPlot[γy,PlotRange->{{-Pi,Pi},{-1/2,1/2}},AxesLabel->{"(Subscript[k, x])","(Subscript[P, y])((Subscript[k, x]))"}]

y方向開放端条件、x方向周期境界条件における、エネルギー固有状態の計算

y ∈ [1, Ly/4] なる領域に波動関数の95 % が局在している固有状態を左端に局在する状態、
 y ∈ [Ly - Ly/4, Ly] なる領域に波動関数の95 % が局在している固有状態を右端に局在する状態とする。

Ly(*y方向のサイト数*)=60;
hOBC=Table[0,{x1,Lx},{x2,Lx}];
Do[hOBC[[y+ny,y]]=Sum[hop[nx,ny]*E^(-I*nx*kx),{nx,-range,range}],{ny,0,range},{y,Ly-ny}];
Do[hOBC[[y,y+ny]]=Sum[hop[nx,-ny]*E^(-I*nx*kx),{nx,-range,range}],{ny,1,range},{y,Ly-ny}];
hOBC=ArrayFlatten[hOBC];
Esys(*wave number, {Eigenvalues},{Eigenvectors}*)=Table[Join[{ikx},Eigensystem[hOBC/.{kx->ikx}]],{ikx,-Pi,Pi,2 Pi/100}];
(*局在状態の計算*)
{LeftEdgeModes,RightEdgeModes}={{},{}};
Do[
Do[
vec=Esys[[k,3,n]];
If[Norm[vec[[;;matsize*Round[Lx/4]]]]>0.95,AppendTo[LeftEdgeModes,{Esys[[k,1]],Esys[[k,2,n]]}];];
If[Norm[vec[[-matsize*Round[Lx/4];;]]]>0.95,AppendTo[RightEdgeModes,{Esys[[k,1]],Esys[[k,2,n]]}];];
,{n,Length[Esys[[k,2]]]}];
,{k,Length[Esys]}];
ListPlot[
{Flatten[Table[{Esys[[i,1]],Esys[[i,2,j]]},{i,Length[Esys]},{j,Length[Esys[[i,2]]]}],1],
LeftEdgeModes,
RightEdgeModes},PlotStyle->{{PointSize[0.005]},{PointSize[0.005],Red},{PointSize[0.005],Green}},PlotLegends->{"Bulk","Left edge","Right edge"}]

x方向の分極 (Subscript[P, x])(ky) ∈ [0,1] の計算

(*Flow*)
Mx=20;
My=100;
(*Berry flux over BZ*)
U=Table[Normalize/@Eigensystem[h/.{kx->ikx,ky->iky}][[2,-l;;]]\[Transpose],{ikx,-Pi,Pi,2*Pi/Mx},{iky,-Pi,Pi,2*Pi/My}];
Do[U[[-1,i]]=U[[1,i]],{i,1,My+1}];
γx=Table[Arg[Product[Det[U[[i+1,j]]\[ConjugateTranspose].U[[i,j]]],{i,1,Mx}]],{j,1,My}];
γx=Table[{2*Pi*(ny-1)/My-Pi,γx[[ny]]/(2*Pi)},{ny,1,My}];
γx=AppendTo[γx,γx[[1]]];
ListPlot[γx,PlotRange->{{-Pi,Pi},{-1/2,1/2}},AxesLabel->{"(Subscript[k, y])","(Subscript[P, x])((Subscript[k, y]))"}]

x方向開放端条件、y方向周期境界条件における、エネルギー固有状態の計算

x ∈ [1, Lx/4] なる領域に波動関数の95 % が局在している固有状態を左端に局在する状態、
 x ∈ [Lx - Lx/4, Lx] なる領域に波動関数の95 % が局在している固有状態を右端に局在する状態とする。

Lx(*x方向のサイト数*)=60;
hOBC=Table[0,{x1,Lx},{x2,Lx}];
Do[hOBC[[x+nx,x]]=Sum[hop[nx,ny]*E^(-I*ny*ky),{ny,-range,range}],{nx,0,range},{x,Lx-nx}];
Do[hOBC[[x,x+nx]]=Sum[hop[-nx,ny]*E^(-I*ny*ky),{ny,-range,range}],{nx,1,range},{x,Lx-nx}];
hOBC=ArrayFlatten[hOBC];
Esys(*wave number, {Eigenvalues},{Eigenvectors}*)=Table[Join[{iky},Eigensystem[hOBC/.{ky->iky}]],{iky,-Pi,Pi,2 Pi/100}];
(*局在状態の計算*)
{LeftEdgeModes,RightEdgeModes}={{},{}};
Do[
Do[
vec=Esys[[k,3,n]];
If[Norm[vec[[;;matsize*Round[Lx/4]]]]>0.95,AppendTo[LeftEdgeModes,{Esys[[k,1]],Esys[[k,2,n]]}];];
If[Norm[vec[[-matsize*Round[Lx/4];;]]]>0.95,AppendTo[RightEdgeModes,{Esys[[k,1]],Esys[[k,2,n]]}];];
,{n,Length[Esys[[k,2]]]}];
,{k,Length[Esys]}];
ListPlot[
{Flatten[Table[{Esys[[i,1]],Esys[[i,2,j]]},{i,Length[Esys]},{j,Length[Esys[[i,2]]]}],1],
LeftEdgeModes,
RightEdgeModes},PlotStyle->{{PointSize[0.005]},{PointSize[0.005],Red},{PointSize[0.005],Green}},PlotLegends->{"Bulk","Left edge","Right edge"}]

空間２次元、時間反転対称性あり

時間反転対称性は (T^2)=-1 （クラスAII）とする。波数空間表示のハミルトニアンに対して、
(Subscript[U, T]) H(SuperStar[((Subscript[k, x]),(Subscript[k, y]))]) (Subscript[U, T])†=H((Subscript[k, x]),(Subscript[k, y])), (Subscript[U, T])(SuperStar[(Subscript[U, T])])=-1
なる条件が課される。

ランダムにハミルトニアンを生成する。

E = 0 において一様なエネルギーギャップが出現しやすいように、ハミルトニアンをトレースレスに再定義する。

ClearAll[kx,ky,x,y,ihop,hop];
$Assumptions={kx,ky,x,y}∈Reals;
matsize(*行列サイズ. 偶数*)=4;
σy={{0,-I},{I,0}};
UT(*時間反転対称性の行列*)=KroneckerProduct[IdentityMatrix[matsize/2],I*σy];
range(*ホッピングのレンジ*)=1;
ξ(*ホッピングの減衰スケール*)=3;
Do[ihop[nx,ny]=rmat[matsize]*E^(-Sqrt[nx^2+ny^2]/ξ),{nx,-range,range},{ny,-range,range}];
Do[hop[nx,ny]=1/2*(ihop[nx,ny]+ihop[-nx,-ny]\[ConjugateTranspose]),{nx,-range,range},{ny,-range,range}];
Do[hop[nx,ny]=hop[nx,ny]-1/matsize*Tr[hop[nx,ny]],{nx,-range,range},{ny,-range,range}];
(*時間反転対称性を課す*)
Do[hop[nx,ny]=1/2*(hop[nx,ny]+UT.hop[nx,ny]\[Conjugate].UT\[ConjugateTranspose]),{nx,-range,range},{ny,-range,range}];
h=Sum[hop[nx,ny]*E^(-I*(nx*kx+ny*ky)),{nx,-range,range},{ny,-range,range}]//Chop;
h=h+100*IdentityMatrix[matsize](*定数シフト*);
Plot3D[Eigenvalues[h]-100,{kx,-Pi,Pi},{ky,-Pi,Pi},Ticks->{Automatic,Automatic,Automatic},AxesLabel->{"(Subscript[k, x])","(Subscript[k, y])","E((Subscript[k, x]),(Subscript[k, y]))"}]

Dirac模型

h = (Sin[kx]*σx + Sin[ky]*σy).τz + (Cos[kx] + Cos[ky] - m)*τx;
UT = I σy;

(*Dirac Model*)
{τ0,τx,τy,τz}=Table[KroneckerProduct[PauliMatrix[0],PauliMatrix[i]],{i,0,3}];
{σ0,σx,σy,σz}=Table[KroneckerProduct[PauliMatrix[i],PauliMatrix[0]],{i,0,3}];
UT=I*σy;
matsize=4;
range=1;
m=1;
hop[0,0]=-m τx//N;
hop[1,0]=-σx.τz/(2*I)+τx/2//N;
hop[-1,0]=σx.τz/(2*I)+τx/2//N;
hop[0,1]=-σy.τz/(2*I)+τx/2//N;
hop[0,-1]=σy.τz/(2*I)+τx/2//N;
hop[1,1]=0*IdentityMatrix[4]//N;
hop[1,-1]=0*IdentityMatrix[4]//N;
hop[-1,1]=0*IdentityMatrix[4]//N;
hop[-1,-1]=0*IdentityMatrix[4]//N;
h=Sum[hop[nx,ny]*E^(-I*(nx*kx+ny*ky)),{nx,-range,range},{ny,-range,range}]//Chop;
h=h+100*IdentityMatrix[matsize](*定数シフト*);
(*h=1/2*(h+h\[ConjugateTranspose])//Simplify//Chop(*エルミート化*)*)
Plot3D[Eigenvalues[h]-100,{kx,-Pi,Pi},{ky,-Pi,Pi},Ticks->{Automatic,Automatic,Automatic},AxesLabel->{"(Subscript[k, x])","(Subscript[k, y])","E((Subscript[k, x]),(Subscript[k, y]))"}]

Kane = Male Z2不変量の計算

Mx=20(*even*);
My=10(*even*);
(*
(*Polarization*)
U=Table[Normalize/@Eigensystem[h/.{kx->ikx,ky->0}][[2,n+1;;]]\[Transpose],{ikx,-Pi,Pi,2*Pi/(2*M)}];
U[[-1]]=U[[1]];
γ0=Arg[Det[Product[U[[j]]\[ConjugateTranspose].U[[j+1]],{j,1,2*M}]]]
U=Table[Normalize/@Eigensystem[h/.{kx->ikx,ky->Pi}][[2,n+1;;]]\[Transpose],{ikx,-Pi,Pi,2*Pi/(2*M)}];
U[[-1]]=U[[1]];
γpi=Arg[Det[Product[U[[j]]\[ConjugateTranspose].U[[j+1]],{j,1,2*M}]]]
*)
U=Table[Normalize/@Eigensystem[N[h/.{kx->ikx,ky->iky}]][[2,-2;;]]\[Transpose],{ikx,-Pi,Pi,2*Pi/Mx},{iky,0,Pi,Pi/My}];
Do[U[[Mx+1,j]]=U[[1,j]],{j,1,My+1}];
(*時間反転分極*)
γT0=Arg[Product[Det[U[[i+1,1]]\[ConjugateTranspose].U[[i,1]]],{i,Mx/2+1,Mx}]*Pf[U[[Mx/2+1,1]]\[ConjugateTranspose].UT.U[[Mx/2+1,1]]\[Conjugate]]/Pf[U[[Mx+1,1]]\[ConjugateTranspose].UT.U[[Mx+1,1]]\[Conjugate]]];
γTpi=Arg[Product[Det[U[[i+1,My+1]]\[ConjugateTranspose].U[[i,My+1]]],{i,Mx/2+1,Mx}]*Pf[U[[Mx/2+1,My+1]]\[ConjugateTranspose].UT.U[[Mx/2+1,My+1]]\[Conjugate]]/Pf[U[[Mx+1,My+1]]\[ConjugateTranspose].UT.U[[Mx+1,My+1]]\[Conjugate]]];
(*Berry flux over BZ/2*)
SumFlux=Total[Table[Arg[Det[U[[i,j]]\[ConjugateTranspose].U[[i,j+1]].U[[i,j+1]]\[ConjugateTranspose].U[[i+1,j+1]].U[[i+1,j+1]]\[ConjugateTranspose].U[[i+1,j]].U[[i+1,j]]\[ConjugateTranspose].U[[i,j]]]],{i,1,Mx},{j,1,My}]//Flatten];
Print["Z2 invariant = ",(SumFlux-2*γT0+2*γTpi)/(2*Pi)];

Wilson線の固有値の流れ

Mx=20(*even*);
My=100(*even*);
U=Table[Normalize/@Eigensystem[h/.{kx->ikx,ky->iky}][[2,-2;;]]\[Transpose],{ikx,-Pi,Pi,2*Pi/Mx},{iky,0,Pi,Pi/My}];
Do[U[[-1,j]]=U[[1,j]],{j,1,My+1}];
W=Table[
U0=IdentityMatrix[2];
Do[U0=U[[i+1,j]]\[ConjugateTranspose].U[[i,j]].U0,{i,1,Mx}];
U0,{j,1,My+1}];
γT=Arg/@Eigenvalues/@W;
PT=Reap[Do[Sow[{γT[[ny,j]]/(2*Pi),Pi*(ny-1)/My}],{ny,1,My+1},{j,2}]][[2,1]];
ListPlot[PT,PlotRange->{{-1/2,1/2},{0,Pi}},AxesLabel->{"(Subscript[P, x])((Subscript[k, y]))","(Subscript[k, y])"}]

y方向開放端条件、x方向周期境界条件における、エネルギー固有状態の計算

y ∈ [1, Ly/4] なる領域に波動関数の95 % が局在している固有状態を左端に局在する状態、
 y ∈ [Ly - Ly/4, Ly] なる領域に波動関数の95 % が局在している固有状態を右端に局在する状態とする。

Ly(*y方向のサイト数*)=60;
hOBC=Table[0,{x1,Lx},{x2,Lx}];
Do[hOBC[[y+ny,y]]=Sum[hop[nx,ny]*E^(-I*nx*kx),{nx,-range,range}],{ny,0,range},{y,Ly-ny}];
Do[hOBC[[y,y+ny]]=Sum[hop[nx,-ny]*E^(-I*nx*kx),{nx,-range,range}],{ny,1,range},{y,Ly-ny}];
hOBC=ArrayFlatten[hOBC];
Esys(*wave number, {Eigenvalues},{Eigenvectors}*)=Table[Join[{ikx},Eigensystem[hOBC/.{kx->ikx}]],{ikx,-Pi,Pi,2 Pi/100}];
(*局在状態の計算*)
{LeftEdgeModes,RightEdgeModes}={{},{}};
Do[
Do[
vec=Esys[[k,3,n]];
If[Norm[vec[[;;matsize*Round[Lx/4]]]]>0.95,AppendTo[LeftEdgeModes,{Esys[[k,1]],Esys[[k,2,n]]}];];
If[Norm[vec[[-matsize*Round[Lx/4];;]]]>0.95,AppendTo[RightEdgeModes,{Esys[[k,1]],Esys[[k,2,n]]}];];
,{n,Length[Esys[[k,2]]]}];
,{k,Length[Esys]}];
ListPlot[
{Flatten[Table[{Esys[[i,1]],Esys[[i,2,j]]},{i,Length[Esys]},{j,Length[Esys[[i,2]]]}],1],
LeftEdgeModes,
RightEdgeModes},PlotStyle->{{PointSize[0.005]},{PointSize[0.005],Red},{PointSize[0.005],Green}},PlotLegends->{"Bulk","Left edge","Right edge"}]

x方向開放端条件、y方向周期境界条件における、エネルギー固有状態の計算

x ∈ [1, Lx/4] なる領域に波動関数の95 % が局在している固有状態を左端に局在する状態、
 x ∈ [Lx - Lx/4, Lx] なる領域に波動関数の95 % が局在している固有状態を右端に局在する状態とする。

Lx(*x方向のサイト数*)=60;
hOBC=Table[0,{x1,Lx},{x2,Lx}];
Do[hOBC[[x+nx,x]]=Sum[hop[nx,ny]*E^(-I*ny*ky),{ny,-range,range}],{nx,0,range},{x,Lx-nx}];
Do[hOBC[[x,x+nx]]=Sum[hop[-nx,ny]*E^(-I*ny*ky),{ny,-range,range}],{nx,1,range},{x,Lx-nx}];
hOBC=ArrayFlatten[hOBC];
Esys(*wave number, {Eigenvalues},{Eigenvectors}*)=Table[Join[{iky},Eigensystem[hOBC/.{ky->iky}]],{iky,-Pi,Pi,2 Pi/100}];
(*局在状態の計算*)
{LeftEdgeModes,RightEdgeModes}={{},{}};
Do[
Do[
vec=Esys[[k,3,n]];
If[Norm[vec[[;;matsize*Round[Lx/4]]]]>0.95,AppendTo[LeftEdgeModes,{Esys[[k,1]],Esys[[k,2,n]]}];];
If[Norm[vec[[-matsize*Round[Lx/4];;]]]>0.95,AppendTo[RightEdgeModes,{Esys[[k,1]],Esys[[k,2,n]]}];];
,{n,Length[Esys[[k,2]]]}];
,{k,Length[Esys]}];
ListPlot[
{Flatten[Table[{Esys[[i,1]],Esys[[i,2,j]]},{i,Length[Esys]},{j,Length[Esys[[i,2]]]}],1],
LeftEdgeModes,
RightEdgeModes},PlotStyle->{{PointSize[0.005]},{PointSize[0.005],Red},{PointSize[0.005],Green}},PlotLegends->{"Bulk","Left edge","Right edge"}]

以下Back up

2D inversion

(*Hamiltonian*)
$Assumptions={kx,ky}∈Reals;
n=2;
σz={{1,0},{0,-1}};
UI=KroneckerProduct[IdentityMatrix[n],σz];
r[]:=RandomReal[{-1,1},{2*n,2*n}]+I*RandomReal[{-1,1},{2*n,2*n}];
h=r[]+r[]*E^(-I kx)+r[]*E^(I kx)+r[]*E^(-I ky)+r[]*E^(I ky);
(*h=r[]+r[]*E^(-I kx)+r[]*E^(I kx)+r[]*E^(-I ky)+r[]*E^(I ky)+0.5*r[]*E^(-2*I kx)+0.5*r[]*E^(2*I kx)+0.5*r[]*E^(-2*I ky)+0.5*r[]*E^(2*I ky)+0.5*r[]*E^(-2*I (kx+ky))+0.5*r[]*E^(2*I (kx+ky))+0.5*r[]*E^(-2*I (kx-ky))+0.5*r[]*E^(2*I (kx-ky));*)
h=(h+h\[ConjugateTranspose]//Simplify);
h=(h+UI.(h/.{kx->-kx,ky->-ky}).UI\[ConjugateTranspose])//Simplify;
h=h+100*IdentityMatrix[2*n];
Plot3D[Eigenvalues[h],{kx,-Pi,Pi},{ky,-Pi,Pi},Ticks->{Automatic,Automatic,None},AxesLabel->{"(Subscript[k, x])","(Subscript[k, y])","E((Subscript[k, x]),(Subscript[k, y]))"}]
(*Flow*)
Mx=30(*even*);
My=100(*even*);
U=Table[Normalize/@Eigensystem[h/.{kx->ikx,ky->iky}][[2,-2;;]]\[Transpose],{ikx,-Pi,Pi,2*Pi/Mx},{iky,-Pi,Pi,2*Pi/My}];
Do[U[[-1,j]]=U[[1,j]],{j,1,My+1}];
W=Table[
U0=IdentityMatrix[2];
Do[U0=U[[i+1,j]]\[ConjugateTranspose].U[[i,j]].U0,{i,1,Mx}];
U0,{j,1,My+1}];
γ=Arg/@Eigenvalues/@W;
P=Reap[Do[Sow[{γ[[ny,j]]/(2*Pi),2*Pi*(ny-1)/My-Pi}],{ny,1,My+1},{j,2}]][[2,1]];
ListPlot[P,PlotRange->{{-1/2,1/2},{-Pi,Pi}},AxesLabel->{"(Subscript[P, x])((Subscript[k, y]))","(Subscript[k, y])"}]
(*Flow*)
Mx=100(*even*);
My=30(*even*);
U=Table[Normalize/@Eigensystem[h/.{kx->ikx,ky->iky}][[2,-2;;]]\[Transpose],{ikx,-Pi,Pi,2*Pi/Mx},{iky,-Pi,Pi,2*Pi/My}];
Do[U[[i,-1]]=U[[i,1]],{i,1,My+1}];
W=Table[
U0=IdentityMatrix[2];
Do[U0=U[[i,j+1]]\[ConjugateTranspose].U[[i,j]].U0,{j,1,My}];
U0,{i,1,Mx+1}];
γ=Arg/@Eigenvalues/@W;
P=Reap[Do[Sow[{γ[[nx,j]]/(2*Pi),2*Pi*(nx-1)/Mx-Pi}],{nx,1,Mx+1},{j,2}]][[2,1]];
ListPlot[P,PlotRange->{{-1/2,1/2},{-Pi,Pi}},AxesLabel->{"(Subscript[P, y])((Subscript[k, x]))","(Subscript[k, x])"}]
ϕ=U[[Mx/2+1,My/2+1]];
Eigenvalues[ϕ\[ConjugateTranspose].UI.ϕ]
ϕ=U[[-1,My/2+1]];
Eigenvalues[ϕ\[ConjugateTranspose].UI.ϕ]
ϕ=U[[Mx/2+1,-1]];
Eigenvalues[ϕ\[ConjugateTranspose].UI.ϕ]
ϕ=U[[-1,-1]];
Eigenvalues[ϕ\[ConjugateTranspose].UI.ϕ]

(*偶然得られたfragile TIっぽくてそうでない例*)
h={{101.90491440507972` -0.3285930629328826` (E^(-I kx))-0.3285930629328826` (E^(I kx))-0.12495342481728589` (E^(-I ky))-0.12495342481728589` (E^(I ky)),(-1.2157182523586476`-1.5678214195509295` I) (E^(-I kx))+(1.2157182523586476` +1.5678214195509295` I) (E^(I kx))+(0.9416274073439248` -0.8588163367063357` I) (E^(-I ky))-(0.9416274073439248` -0.8588163367063357` I) (E^(I ky)),(-0.05608449156392403`-3.3649434566941956` I)+(0.6988158630450041` +0.21948907817640428` I) (E^(-I kx))+(0.6988158630450041` +0.21948907817640428` I) (E^(I kx))+(0.6000162610724269` -1.1331640437345638` I) (E^(-I ky))+(0.6000162610724269` -1.1331640437345638` I) (E^(I ky)),(-0.5152083728739338`+0.1026163806458209` I) (E^(-I kx))+(0.5152083728739338` -0.1026163806458209` I) (E^(I kx))-(2.1776750386138093` +0.6121284391886888` I) (E^(-I ky))+(2.1776750386138093` +0.6121284391886888` I) (E^(I ky))},{(1.2157182523586476` -1.5678214195509295` I) (E^(-I kx))-(1.2157182523586476` -1.5678214195509295` I) (E^(I kx))-(0.9416274073439248` +0.8588163367063357` I) (E^(-I ky))+(0.9416274073439248` +0.8588163367063357` I) (E^(I ky)),98.46093662867843` -0.6761753002252666` (E^(-I kx))-0.6761753002252666` (E^(I kx))-0.4985054137339917` (E^(-I ky))-0.4985054137339917` (E^(I ky)),(1.5665649108993969` +0.008662194950003421` I) (E^(-I kx))-(1.5665649108993969` +0.008662194950003421` I) (E^(I kx))-(2.1436676360125193` +1.6285731529800245` I) (E^(-I ky))+(2.1436676360125193` +1.6285731529800245` I) (E^(I ky)),(-1.0560525375580214`-0.2150644294212931` I)-(0.15884978738464905` -0.8288833106205993` I) (E^(-I kx))-(0.15884978738464905` -0.8288833106205993` I) (E^(I kx))-(0.4578733396769632` -0.3506204045229082` I) (E^(-I ky))-(0.4578733396769632` -0.3506204045229082` I) (E^(I ky))},{(-0.05608449156392403`+3.3649434566941956` I)+(0.6988158630450041` -0.21948907817640428` I) (E^(-I kx))+(0.6988158630450041` -0.21948907817640428` I) (E^(I kx))+(0.6000162610724269` +1.1331640437345638` I) (E^(-I ky))+(0.6000162610724269` +1.1331640437345638` I) (E^(I ky)),(-1.5665649108993969`+0.008662194950003421` I) (E^(-I kx))+(1.5665649108993969` -0.008662194950003421` I) (E^(I kx))+(2.1436676360125193` -1.6285731529800245` I) (E^(-I ky))-(2.1436676360125193` -1.6285731529800245` I) (E^(I ky)),97.68182319754189` -0.21353521730101654` (E^(-I kx))-0.21353521730101654` (E^(I kx))-1.9403470680962815` (E^(-I ky))-1.9403470680962815` (E^(I ky)),(-0.5366388109331233`-1.4726237499143071` I) (E^(-I kx))+(0.5366388109331233` +1.4726237499143071` I) (E^(I kx))+(2.5112486154735056` -0.7925231211396229` I) (E^(-I ky))-(2.5112486154735056` -0.7925231211396229` I) (E^(I ky))},{(0.5152083728739338` +0.1026163806458209` I) (E^(-I kx))-(0.5152083728739338` +0.1026163806458209` I) (E^(I kx))+(2.1776750386138093` -0.6121284391886888` I) (E^(-I ky))-(2.1776750386138093` -0.6121284391886888` I) (E^(I ky)),(-1.0560525375580214`+0.2150644294212931` I)-(0.15884978738464905` +0.8288833106205993` I) (E^(-I kx))-(0.15884978738464905` +0.8288833106205993` I) (E^(I kx))-(0.4578733396769632` +0.3506204045229082` I) (E^(-I ky))-(0.4578733396769632` +0.3506204045229082` I) (E^(I ky)),(0.5366388109331233` -1.4726237499143071` I) (E^(-I kx))-(0.5366388109331233` -1.4726237499143071` I) (E^(I kx))-(2.5112486154735056` +0.7925231211396229` I) (E^(-I ky))+(2.5112486154735056` +0.7925231211396229` I) (E^(I ky)),102.63883528690909` -0.3132212143108566` (E^(-I kx))-0.3132212143108566` (E^(I kx))+0.5117915491477429` (E^(-I ky))+0.5117915491477429` (E^(I ky))}};

(*Dirac Model plus perturbation*)
{τ0,τx,τy,τz}=Table[KroneckerProduct[PauliMatrix[0],PauliMatrix[i]],{i,0,3}];
{σ0,σx,σy,σz}=Table[KroneckerProduct[PauliMatrix[i],PauliMatrix[0]],{i,0,3}];
UI=τz;
n=2;
h=(Sin[kx]*σx+Sin[ky]*σy).τx+(Cos[kx]+Cos[ky]-1)*τz;
r[]:=RandomReal[{-1,1},{2*n,2*n}]+I*RandomReal[{-1,1},{2*n,2*n}];
V=r[]+r[]*E^(-I kx)+r[]*E^(I kx)+r[]*E^(-I ky)+r[]*E^(I ky);
V=(V+UI.(V/.{kx->-kx,ky->-ky}).UI\[ConjugateTranspose])//Simplify;
V=V+V\[ConjugateTranspose];
h=h+0.1*V;
h=h+100*IdentityMatrix[2*n];
Plot3D[Eigenvalues[h],{kx,-Pi,Pi},{ky,-Pi,Pi}]
(*Flow*)
Mx=30(*even*);
My=100(*even*);
U=Table[Normalize/@Eigensystem[h/.{kx->ikx,ky->iky}][[2,-2;;]]\[Transpose],{ikx,-Pi,Pi,2*Pi/Mx},{iky,-Pi,Pi,2*Pi/My}];
Do[U[[-1,j]]=U[[1,j]],{j,1,My+1}];
W=Table[
U0=IdentityMatrix[2];
Do[U0=U[[i+1,j]]\[ConjugateTranspose].U[[i,j]].U0,{i,1,Mx}];
U0,{j,1,My+1}];
γ=Arg/@Eigenvalues/@W;
P=Reap[Do[Sow[{γ[[ny,j]]/(2*Pi),2*Pi*(ny-1)/My-Pi}],{ny,1,My+1},{j,2}]][[2,1]];
ListPlot[P,PlotRange->{{-1/2,1/2},{-Pi,Pi}},AxesLabel->{"(Subscript[P, x])((Subscript[k, y]))","(Subscript[k, y])"}]
(*Flow*)
Mx=100(*even*);
My=30(*even*);
U=Table[Normalize/@Eigensystem[h/.{kx->ikx,ky->iky}][[2,-2;;]]\[Transpose],{ikx,-Pi,Pi,2*Pi/Mx},{iky,-Pi,Pi,2*Pi/My}];
Do[U[[i,-1]]=U[[i,1]],{i,1,My+1}];
W=Table[
U0=IdentityMatrix[2];
Do[U0=U[[i,j+1]]\[ConjugateTranspose].U[[i,j]].U0,{j,1,My}];
U0,{i,1,Mx+1}];
γ=Arg/@Eigenvalues/@W;
P=Reap[Do[Sow[{γ[[nx,j]]/(2*Pi),2*Pi*(nx-1)/Mx-Pi}],{nx,1,Mx+1},{j,2}]][[2,1]];
ListPlot[P,PlotRange->{{-1/2,1/2},{-Pi,Pi}},AxesLabel->{"(Subscript[P, y])((Subscript[k, x]))","(Subscript[k, x])"}]

ϕ=U[[Mx/2+1,My/2+1]];
Eigenvalues[ϕ\[ConjugateTranspose].UI.ϕ]
ϕ=U[[-1,My/2+1]];
Eigenvalues[ϕ\[ConjugateTranspose].UI.ϕ]
ϕ=U[[Mx/2+1,-1]];
Eigenvalues[ϕ\[ConjugateTranspose].UI.ϕ]
ϕ=U[[-1,-1]];
Eigenvalues[ϕ\[ConjugateTranspose].UI.ϕ]

Back Up

(*U(2) Wilson loop of BZ/2*)
Mx=40(*even*);
My=20(*even*);
U=Table[Normalize/@Eigensystem[N[h/.{kx->ikx,ky->iky}]][[2,2n-1;;]]\[Transpose],{ikx,-Pi,Pi,2*Pi/Mx},{iky,0,Pi,Pi/My}];
Do[U[[Mx+1,j]]=U[[1,j]],{j,1,My+1}];
U0=IdentityMatrix[2];
Do[U0=U0.U[[i,1]]\[ConjugateTranspose].U[[i+1,1]],{i,1,Mx}];
Do[U0=U0.U[[Mx+1,j]]\[ConjugateTranspose].U[[Mx+1,j+1]],{j,1,My}];
Do[U0=U0.U[[i+1,My+1]]\[ConjugateTranspose].U[[i,My+1]],{i,Mx,1,-1}];
Do[U0=U0.U[[1,j+1]]\[ConjugateTranspose].U[[1,j]],{j,My,1,-1}];
U0//MatrixForm
WUn=Arg[Tr[U0]]
WU1=Arg[Det[U0]]
2*WUn-WU1

(*U(n) Wilson loop of BZ/2. With Polar decomposition*)
Mx=40(*even*);
My=20(*even*);
U=Table[Normalize/@Eigensystem[N[h/.{kx->ikx,ky->iky}]][[2,n+1;;]]\[Transpose],{ikx,-Pi,Pi,2*Pi/Mx},{iky,0,Pi,Pi/My}];
Do[U[[Mx+1,j]]=U[[1,j]],{j,1,My+1}];
U0=IdentityMatrix[n];
Do[SVD=SingularValueDecomposition[U[[i,1]]\[ConjugateTranspose].U[[i+1,1]]];Q=SVD[[1]].SVD[[3]]\[ConjugateTranspose];U0=U0.Q;,{i,1,Mx}];
Do[SVD=SingularValueDecomposition[U[[Mx+1,j]]\[ConjugateTranspose].U[[Mx+1,j+1]]];Q=SVD[[1]].SVD[[3]]\[ConjugateTranspose];U0=U0.Q;,{j,1,My}];
Do[SVD=SingularValueDecomposition[U[[i+1,My+1]]\[ConjugateTranspose].U[[i,My+1]]];Q=SVD[[1]].SVD[[3]]\[ConjugateTranspose];U0=U0.Q;,{i,Mx,1,-1}];
Do[SVD=SingularValueDecomposition[U[[1,j+1]]\[ConjugateTranspose].U[[1,j]]];Q=SVD[[1]].SVD[[3]]\[ConjugateTranspose];U0=U0.Q,{j,My,1,-1}];
U0//MatrixForm
WUn=Arg[1/2*Tr[U0]]
WU1=Arg[Det[U0]]
2*WUn-WU1

(*Wilson loop of BZ*)
Mx=40(*even*);
My=40(*even*);
U=Table[Normalize/@Eigensystem[N[h/.{kx->ikx,ky->iky}]][[2,n+1;;]]\[Transpose],{ikx,-Pi,Pi,2*Pi/Mx},{iky,-Pi,Pi,2*Pi/My}];
Do[U[[Mx+1,j]]=U[[1,j]],{j,1,My+1}];
Do[U[[i,My+1]]=U[[i,1]],{i,1,Mx+1}];
U0=IdentityMatrix[n];
Do[U0=U0.U[[i,1]]\[ConjugateTranspose].U[[i+1,1]],{i,1,Mx}];
Do[U0=U0.U[[Mx+1,j]]\[ConjugateTranspose].U[[Mx+1,j+1]],{j,1,My}];
Do[U0=U0.U[[i+1,My+1]]\[ConjugateTranspose].U[[i,My+1]],{i,Mx,1,-1}];
Do[U0=U0.U[[1,j+1]]\[ConjugateTranspose].U[[1,j]],{j,My,1,-1}];
U0
Tr[U0]
Det[U0]