PT symmetry breaking

Without PT symmetry

rmat[n_]:=RandomReal[{-1,1},{n,n}]+I*RandomReal[{-1,1},{n,n}];
Clear[h];
n=20;
h=(1-t)*rmat[n]+t*rmat[n];
Manipulate[ListPlot[{Re[#],Im[#]}&/@Eigenvalues[h/.{t->it}],PlotRange->{{-3,3},{-3,3}},AxesLabel->{"Re","Im"}],{it,0,1}]

With PT symmetry

rmat[n_]:=RandomReal[{-1,1},{n,n}]+I*RandomReal[{-1,1},{n,n}];
Clear[h];
n=20;
h=(1-t)*rmat[n]+t*rmat[n];
h=1/2*(h+h\[Conjugate]);
Manipulate[ListPlot[{Re[#],Im[#]}&/@Eigenvalues[h/.{t->it}],PlotRange->{{-3,3},{-3,3}},AxesLabel->{"Re","Im"}],{it,0,1}]

Absence of Non-Hermitian skin effect with inversion symmetry
U H(k) U†=H(-k), (U^2)=1

(*Random hoping model*)
rmat[n_]:=RandomReal[{-1,1},{n,n}]+I*RandomReal[{-1,1},{n,n}];
Clear[kx,x,hop];
$Assumptions={kx,x}∈Reals;
n(*matrix size*)=2(*even*);
U=KroneckerProduct[IdentityMatrix[n/2],PauliMatrix[3]];
range(*hopping range*)=3;
ξ(*Scale of exp decaying*)=2;
Do[hop[nx]=rmat[n]*E^(-Abs[nx]/ξ),{nx,-range,range}];
hop[0]=1/2*(hop[0]+U.hop[0].U\[ConjugateTranspose]);
Do[hop[-nx]=U.hop[nx].U\[ConjugateTranspose],{nx,1,range}];
Lx(*Site length*)=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];
{Eval,Evec}=Eigensystem[hOBC];
lis=RandomInteger[{1,n*Lx},5];
hPBC=Sum[hop[nx]*E^(-I*nx*kx),{nx,-range,range}];
combinedData={
{Re[#],Im[#]}&/@Eval,
Flatten[Table[{Re[#],Im[#]}&/@Eigenvalues[hPBC/.{kx->2*Pi*px/(10*Lx)}],{px,1,10*Lx}],1]
};
(*Detection of edge states*)
edge={};
Do[
v=Evec[[i]];
If[Norm[v[[;;n*Round[Lx/10]]]]^2+Norm[v[[-n*Round[Lx/10];;]]]^2>0.9,AppendTo[edge,i]];
,{i,Evec//Length}];
ListPlot[combinedData,PlotLabel->StringJoin["Lx = ",ToString[Lx],", n = ",ToString[n]],PlotStyle->{Blue,Green},PlotLegends->PointLegend[{"OBC","PBC"},LegendMarkers->{"●",15}],AxesLabel->{"Re","Im"},PlotRange->All,Epilog->(Text[Style[ToString[Position[lis,#][[1,1]]],Red,Bold,20],{Re[Eval[[#]]],Im[Eval[[#]]]}]&/@lis)~Join~(Text[Style["e",Red,Bold,20],{Re[Eval[[#]]],Im[Eval[[#]]]+0.2}]&/@edge)]
Table[ListPlot[#^2&/@Norm/@ArrayReshape[Evec[[lis[[i]]]],{Lx,n}],PlotRange->All,AxesLabel->{"x","||(Subscript[v, x])((||)^2)"},PlotLabel->Text[Style[ToString[i],Red,Bold,20]]],{i,lis//Length}]//Print;
Table[ListPlot[#^2&/@Norm/@ArrayReshape[Evec[[edge[[i]]]],{Lx,n}],PlotRange->All,AxesLabel->{"x","||(Subscript[v, x])((||)^2)"},PlotLabel->Text[Style["edge state",Red,20]]],{i,edge//Length}]//Print;

Non-Hermitian skin effect for 38 symmetry classes

rmat[n_]:=RandomReal[{-1,1},{n,n}]+I*RandomReal[{-1,1},{n,n}];

makehop[]:=Module[{},
ClearAll[hop];
Do[hop[nx]=rmat[n]*E^(-Abs[nx]/ξ),{nx,-range,range}];
];

symmetrizehop[ϕ_,η_,c_]:=Switch[{ϕ,η},
{1,1},
Do[hop[nx]=1/2*(hop[nx]+c*U.hop[nx].U\[ConjugateTranspose]),{nx,-range,range}];,
{-1,1},
Do[hop[nx]=1/2*(hop[nx]+c*U.hop[nx]\[Conjugate].U\[ConjugateTranspose]),{nx,-range,range}];,
{1,-1},
hop[0]=1/2*(hop[0]+c*U.hop[0]\[ConjugateTranspose].U\[ConjugateTranspose]);Do[hop[-nx]=c*U.hop[nx]\[ConjugateTranspose].U\[ConjugateTranspose],{nx,1,range}];,
{-1,-1},
hop[0]=1/2*(hop[0]+c*U.hop[0]\[Transpose].U\[ConjugateTranspose]);Do[hop[-nx]=c*U.hop[nx]\[Transpose].U\[ConjugateTranspose],{nx,1,range}];
];

spectra[hop_,Lx_]:=Module[{kx},
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];
{Eval,Evec}=Eigensystem[hOBC];
lis=RandomInteger[{1,n*Lx},5];
hPBC=Sum[hop[nx]*E^(-I*nx*kx),{nx,-range,range}];
combinedData={
{Re[#],Im[#]}&/@Eval,
Flatten[Table[{Re[#],Im[#]}&/@Eigenvalues[hPBC/.{kx->2*Pi*px/(10*Lx)}],{px,1,10*Lx}],1]
};
ListPlot[combinedData,PlotLabel->StringJoin["Lx = ",ToString[Lx],", n = ",ToString[n]],PlotStyle->{Blue,Green},PlotLegends->PointLegend[{"OBC","PBC"},LegendMarkers->{"●",15}],AxesLabel->{"Re","Im"},PlotRange->All,Epilog->(Text[Style[ToString[Position[lis,#][[1,1]]],Red,Bold,20],{Re[Eval[[#]]],Im[Eval[[#]]]}]&/@lis)]//Print;
Table[ListPlot[#^2&/@Norm/@ArrayReshape[Evec[[lis[[i]]]],{Lx,n}],PlotRange->All,AxesLabel->{"x","||(Subscript[v, x])((||)^2)"},PlotLabel->Text[Style[ToString[i],Red,Bold,20]]],{i,lis//Length}]//Print;
];

spectrawithEdge[hop_,Lx_]:=Module[{kx},
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];
{Eval,Evec}=Eigensystem[hOBC];
lis=RandomInteger[{1,n*Lx},5];
hPBC=Sum[hop[nx]*E^(-I*nx*kx),{nx,-range,range}];
combinedData={
{Re[#],Im[#]}&/@Eval,
Flatten[Table[{Re[#],Im[#]}&/@Eigenvalues[hPBC/.{kx->2*Pi*px/(10*Lx)}],{px,1,10*Lx}],1]
};
(*Detection of edge states*)
edge={};
Do[
v=Evec[[i]];
If[Norm[v[[;;n*Round[Lx/10]]]]^2+Norm[v[[-n*Round[Lx/10];;]]]^2>0.9,AppendTo[edge,i]];
,{i,Evec//Length}];
ListPlot[combinedData,PlotLabel->StringJoin["Lx = ",ToString[Lx],", n = ",ToString[n]],PlotStyle->{Blue,Green},PlotLegends->PointLegend[{"OBC","PBC"},LegendMarkers->{"●",15}],AxesLabel->{"Re","Im"},PlotRange->All,Epilog->(Text[Style[ToString[Position[lis,#][[1,1]]],Red,Bold,20],{Re[Eval[[#]]],Im[Eval[[#]]]}]&/@lis)~Join~(Text[Style["e",Red,Bold,20],{Re[Eval[[#]]],Im[Eval[[#]]]+0.2}]&/@edge)]//Print;
Table[ListPlot[#^2&/@Norm/@ArrayReshape[Evec[[lis[[i]]]],{Lx,n}],PlotRange->All,AxesLabel->{"x","||(Subscript[v, x])((||)^2)"},PlotLabel->Text[Style[ToString[i],Red,Bold,20]]],{i,lis//Length}]//Print;
Table[ListPlot[#^2&/@Norm/@ArrayReshape[Evec[[edge[[i]]]],{Lx,n}],PlotRange->All,AxesLabel->{"x","||(Subscript[v, x])((||)^2)"},PlotLabel->Text[Style["edge state",Red,20]]],{i,edge//Length}]//Print;
];

(0) Hermitian systems

{n(*matrix size*),range(*hopping range*),ξ(*Scale of exp decaying*)}={2,2,2};
makehop[];
{ϕ,η,c}={1,-1,1};U=IdentityMatrix[n];
symmetrizehop[ϕ,η,c];
Lx(*Site length*)=100;
spectra[hop,Lx]

N = 0 (n∈Z)

(1) No symmetry

{n(*matrix size*),range(*hopping range*),ξ(*Scale of exp decaying*)}={3,2,1};
makehop[];
Lx(*Site length*)=50;
spectra[hop,Lx]
Lx(*Site length*)=500;
spectra[hop,Lx]

N = 1 (n∈2Z)

(2) Pseudo Hermiticity

{n(*matrix size, should be even*),range(*hopping range*),ξ(*Scale of exp decaying*)}={4,2,2};
makehop[];
{ϕ,η,c}={1,-1,1};U=KroneckerProduct[IdentityMatrix[n/2],PauliMatrix[3]];
symmetrizehop[ϕ,η,c];
Lx(*Site length*)=100;
spectra[hop,Lx]

(3) Sublattice Symmetry

{n(*matrix size, should be even*),range(*hopping range*),ξ(*Scale of exp decaying*)}={2,1,2};
makehop[];
{ϕ,η,c}={1,1,-1};U=KroneckerProduct[IdentityMatrix[n/2],PauliMatrix[3]];
symmetrizehop[ϕ,η,c];
(*hop[0]=0*hop[0];*)
Lx(*Site length*)=50;
spectra[hop,Lx]

(4) Chiral Symmetry

{n(*matrix size, should be even*),range(*hopping range*),ξ(*Scale of exp decaying*)}={2,2,2};
makehop[];
{ϕ,η,c}={1,-1,-1};U=KroneckerProduct[IdentityMatrix[n/2],PauliMatrix[3]];
symmetrizehop[ϕ,η,c];
Lx(*Site length*)=100;
spectra[hop,Lx]

(5) AI

{n(*matrix size, should be even*),range(*hopping range*),ξ(*Scale of exp decaying*)}={4,2,2};
makehop[];
{ϕ,η,c}={-1,1,1};U=KroneckerProduct[IdentityMatrix[n/2],PauliMatrix[3]];
symmetrizehop[ϕ,η,c];
Lx(*Site length*)=100;
spectra[hop,Lx]

(6) AII

{n(*matrix size, should be even*),range(*hopping range*),ξ(*Scale of exp decaying*)}={4,2,2};
makehop[];
{ϕ,η,c}={-1,1,1};U=KroneckerProduct[IdentityMatrix[n/2],PauliMatrix[2]];
symmetrizehop[ϕ,η,c];
Lx(*Site length*)=100;
spectra[hop,Lx]

(7) AI†

{n(*matrix size, should be even*),range(*hopping range*),ξ(*Scale of exp decaying*)}={4,2,2};
makehop[];
{ϕ,η,c}={-1,-1,1};U=KroneckerProduct[IdentityMatrix[n/2],PauliMatrix[3]];
symmetrizehop[ϕ,η,c];
Lx(*Site length*)=100;
spectrawithEdge[hop,Lx]

(8) AII†

{n(*matrix size, should be even*),range(*hopping range*),ξ(*Scale of exp decaying*)}={4,2,2};
makehop[];
{ϕ,η,c}={-1,-1,1};U=KroneckerProduct[IdentityMatrix[n/2],PauliMatrix[2]];
symmetrizehop[ϕ,η,c];
Lx(*Site length*)=100;
spectra[hop,Lx]

(9) D

{n(*matrix size, should be even*),range(*hopping range*),ξ(*Scale of exp decaying*)}={4,2,2};
makehop[];
{ϕ,η,c}={-1,-1,-1};U=KroneckerProduct[IdentityMatrix[n/2],PauliMatrix[3]];
symmetrizehop[ϕ,η,c];
Lx(*Site length*)=100;
spectra[hop,Lx]

(10) C

{n(*matrix size, should be even*),range(*hopping range*),ξ(*Scale of exp decaying*)}={4,2,2};
makehop[];
{ϕ,η,c}={-1,-1,-1};U=KroneckerProduct[IdentityMatrix[n/2],PauliMatrix[2]];
symmetrizehop[ϕ,η,c];
Lx(*Site length*)=100;
spectra[hop,Lx]

(11) D†

{n(*matrix size, should be even*),range(*hopping range*),ξ(*Scale of exp decaying*)}={4,2,2};
makehop[];
{ϕ,η,c}={-1,1,-1};U=KroneckerProduct[IdentityMatrix[n/2],PauliMatrix[3]];
symmetrizehop[ϕ,η,c];
Lx(*Site length*)=100;
spectra[hop,Lx]

(12) C†

{n(*matrix size, should be even*),range(*hopping range*),ξ(*Scale of exp decaying*)}={4,2,2};
makehop[];
{ϕ,η,c}={-1,1,-1};U=KroneckerProduct[IdentityMatrix[n/2],PauliMatrix[2]];
symmetrizehop[ϕ,η,c];
Lx(*Site length*)=100;
spectra[hop,Lx]

N = 2 (n∈4Z)

(13) BDI

{n(*matrix size, should be even*),range(*hopping range*),ξ(*Scale of exp decaying*)}={4,2,2};
makehop[];
{ϕ,η,c}={-1,1,1};U=KroneckerProduct[IdentityMatrix[n/4],PauliMatrix[0],PauliMatrix[3]];
symmetrizehop[ϕ,η,c];
{ϕ,η,c}={-1,-1,-1};U=KroneckerProduct[IdentityMatrix[n/4],PauliMatrix[3],PauliMatrix[0]];
symmetrizehop[ϕ,η,c];
Lx(*Site length*)=100;
spectra[hop,Lx]

(14) DIII

{n(*matrix size, should be 4Z*),range(*hopping range*),ξ(*Scale of exp decaying*)}={4,2,2};
makehop[];
{ϕ,η,c}={-1,1,1};U=KroneckerProduct[IdentityMatrix[n/4],PauliMatrix[0],PauliMatrix[2]];
symmetrizehop[ϕ,η,c];
{ϕ,η,c}={-1,-1,-1};U=KroneckerProduct[IdentityMatrix[n/4],PauliMatrix[3],PauliMatrix[0]];
symmetrizehop[ϕ,η,c];
Lx(*Site length*)=100;
spectra[hop,Lx]

(15) CI

{n(*matrix size, should be even*),range(*hopping range*),ξ(*Scale of exp decaying*)}={4,2,2};
makehop[];
{ϕ,η,c}={-1,1,1};U=KroneckerProduct[IdentityMatrix[n/4],PauliMatrix[0],PauliMatrix[3]];
symmetrizehop[ϕ,η,c];
{ϕ,η,c}={-1,-1,-1};U=KroneckerProduct[IdentityMatrix[n/4],PauliMatrix[2],PauliMatrix[0]];
symmetrizehop[ϕ,η,c];
Lx(*Site length*)=100;
spectra[hop,Lx]

(16) CII

{n(*matrix size, should be 4Z*),range(*hopping range*),ξ(*Scale of exp decaying*)}={4,2,2};
makehop[];
{ϕ,η,c}={-1,1,1};U=KroneckerProduct[IdentityMatrix[n/4],PauliMatrix[0],PauliMatrix[2]];
symmetrizehop[ϕ,η,c];
{ϕ,η,c}={-1,-1,-1};U=KroneckerProduct[IdentityMatrix[n/4],PauliMatrix[2],PauliMatrix[0]];
symmetrizehop[ϕ,η,c];
Lx(*Site length*)=100;
spectra[hop,Lx]

(17) BDI†

{n(*matrix size, should be even*),range(*hopping range*),ξ(*Scale of exp decaying*)}={4,2,2};
makehop[];
{ϕ,η,c}={-1,-1,1};U=KroneckerProduct[IdentityMatrix[n/4],PauliMatrix[0],PauliMatrix[3]];
symmetrizehop[ϕ,η,c];
{ϕ,η,c}={-1,1,-1};U=KroneckerProduct[IdentityMatrix[n/4],PauliMatrix[3],PauliMatrix[0]];
symmetrizehop[ϕ,η,c];
Lx(*Site length*)=100;
spectrawithEdge[hop,Lx]

(18) DIII†

{n(*matrix size, should be 4Z*),range(*hopping range*),ξ(*Scale of exp decaying*)}={4,2,2};
makehop[];
{ϕ,η,c}={-1,-1,1};U=KroneckerProduct[IdentityMatrix[n/4],PauliMatrix[0],PauliMatrix[2]];
symmetrizehop[ϕ,η,c];
{ϕ,η,c}={-1,1,-1};U=KroneckerProduct[IdentityMatrix[n/4],PauliMatrix[3],PauliMatrix[0]];
symmetrizehop[ϕ,η,c];
Lx(*Site length*)=100;
spectra[hop,Lx]

(19) CI†

{n(*matrix size, should be even*),range(*hopping range*),ξ(*Scale of exp decaying*)}={4,2,2};
makehop[];
{ϕ,η,c}={-1,-1,1};U=KroneckerProduct[IdentityMatrix[n/4],PauliMatrix[0],PauliMatrix[3]];
symmetrizehop[ϕ,η,c];
{ϕ,η,c}={-1,1,-1};U=KroneckerProduct[IdentityMatrix[n/4],PauliMatrix[2],PauliMatrix[0]];
symmetrizehop[ϕ,η,c];
Lx(*Site length*)=100;
spectrawithEdge[hop,Lx]

(20) CII†

{n(*matrix size, should be 4Z*),range(*hopping range*),ξ(*Scale of exp decaying*)}={4,2,2};
makehop[];
{ϕ,η,c}={-1,-1,1};U=KroneckerProduct[IdentityMatrix[n/4],PauliMatrix[0],PauliMatrix[2]];
symmetrizehop[ϕ,η,c];
{ϕ,η,c}={-1,1,-1};U=KroneckerProduct[IdentityMatrix[n/4],PauliMatrix[2],PauliMatrix[0]];
symmetrizehop[ϕ,η,c];
Lx(*Site length*)=100;
spectra[hop,Lx]

(21) TRS+TRS† ((Subscript[u, 1])(SuperStar[(Subscript[u, 1])])=1, (Subscript[u, 2])(SuperStar[(Subscript[u, 2])])=1)

{n(*matrix size, should be even*),range(*hopping range*),ξ(*Scale of exp decaying*)}={4,2,2};
makehop[];
{ϕ,η,c}={-1,1,1};U=KroneckerProduct[IdentityMatrix[n/4],PauliMatrix[0],PauliMatrix[3]];
symmetrizehop[ϕ,η,c];
{ϕ,η,c}={-1,-1,1};U=KroneckerProduct[IdentityMatrix[n/4],PauliMatrix[3],PauliMatrix[0]];
symmetrizehop[ϕ,η,c];
Lx(*Site length*)=100;
spectrawithEdge[hop,Lx]

(22) TRS+TRS† ((Subscript[u, 1])(SuperStar[(Subscript[u, 1])])=-1, (Subscript[u, 2])(SuperStar[(Subscript[u, 2])])=1)

{n(*matrix size, should be even*),range(*hopping range*),ξ(*Scale of exp decaying*)}={4,2,2};
makehop[];
{ϕ,η,c}={-1,1,1};U=KroneckerProduct[IdentityMatrix[n/4],PauliMatrix[0],PauliMatrix[2]];
symmetrizehop[ϕ,η,c];
{ϕ,η,c}={-1,-1,1};U=KroneckerProduct[IdentityMatrix[n/4],PauliMatrix[3],PauliMatrix[0]];
symmetrizehop[ϕ,η,c];
Lx(*Site length*)=100;
spectrawithEdge[hop,Lx]

(23) TRS+TRS† ((Subscript[u, 1])(SuperStar[(Subscript[u, 1])])=1, (Subscript[u, 2])(SuperStar[(Subscript[u, 2])])=-1)

{n(*matrix size, should be even*),range(*hopping range*),ξ(*Scale of exp decaying*)}={4,2,2};
makehop[];
{ϕ,η,c}={-1,1,1};U=KroneckerProduct[IdentityMatrix[n/4],PauliMatrix[0],PauliMatrix[3]];
symmetrizehop[ϕ,η,c];
{ϕ,η,c}={-1,-1,1};U=KroneckerProduct[IdentityMatrix[n/4],PauliMatrix[2],PauliMatrix[0]];
symmetrizehop[ϕ,η,c];
Lx(*Site length*)=100;
spectra[hop,Lx]

(24) TRS+TRS† ((Subscript[u, 1])(SuperStar[(Subscript[u, 1])])=-1, (Subscript[u, 2])(SuperStar[(Subscript[u, 2])])=-1)

{n(*matrix size, should be even*),range(*hopping range*),ξ(*Scale of exp decaying*)}={4,2,2};
makehop[];
{ϕ,η,c}={-1,1,1};U=KroneckerProduct[IdentityMatrix[n/4],PauliMatrix[0],PauliMatrix[2]];
symmetrizehop[ϕ,η,c];
{ϕ,η,c}={-1,-1,1};U=KroneckerProduct[IdentityMatrix[n/4],PauliMatrix[2],PauliMatrix[0]];
symmetrizehop[ϕ,η,c];
Lx(*Site length*)=100;
spectra[hop,Lx]

(25) TRS+PHS† ((Subscript[u, 1])(SuperStar[(Subscript[u, 1])])=1, (Subscript[u, 2])(SuperStar[(Subscript[u, 2])])=1)

{n(*matrix size, should be even*),range(*hopping range*),ξ(*Scale of exp decaying*)}={4,2,2};
makehop[];
{ϕ,η,c}={-1,1,1};U=KroneckerProduct[IdentityMatrix[n/4],PauliMatrix[0],PauliMatrix[3]];
symmetrizehop[ϕ,η,c];
{ϕ,η,c}={-1,1,-1};U=KroneckerProduct[IdentityMatrix[n/4],PauliMatrix[3],PauliMatrix[0]];
symmetrizehop[ϕ,η,c];
Lx(*Site length*)=100;
spectra[hop,Lx]

(26) TRS+PHS† ((Subscript[u, 1])(SuperStar[(Subscript[u, 1])])=-1, (Subscript[u, 2])(SuperStar[(Subscript[u, 2])])=1)

{n(*matrix size, should be even*),range(*hopping range*),ξ(*Scale of exp decaying*)}={4,2,2};
makehop[];
{ϕ,η,c}={-1,1,1};U=KroneckerProduct[IdentityMatrix[n/4],PauliMatrix[0],PauliMatrix[2]];
symmetrizehop[ϕ,η,c];
{ϕ,η,c}={-1,1,-1};U=KroneckerProduct[IdentityMatrix[n/4],PauliMatrix[3],PauliMatrix[0]];
symmetrizehop[ϕ,η,c];
Lx(*Site length*)=100;
spectra[hop,Lx]

(27) TRS+PHS† ((Subscript[u, 1])(SuperStar[(Subscript[u, 1])])=1, (Subscript[u, 2])(SuperStar[(Subscript[u, 2])])=-1)

{n(*matrix size, should be even*),range(*hopping range*),ξ(*Scale of exp decaying*)}={4,2,2};
makehop[];
{ϕ,η,c}={-1,1,1};U=KroneckerProduct[IdentityMatrix[n/4],PauliMatrix[0],PauliMatrix[3]];
symmetrizehop[ϕ,η,c];
{ϕ,η,c}={-1,1,-1};U=KroneckerProduct[IdentityMatrix[n/4],PauliMatrix[2],PauliMatrix[0]];
symmetrizehop[ϕ,η,c];
Lx(*Site length*)=100;
spectra[hop,Lx]

(28) TRS+PHS† ((Subscript[u, 1])(SuperStar[(Subscript[u, 1])])=-1, (Subscript[u, 2])(SuperStar[(Subscript[u, 2])])=-1)

{n(*matrix size, should be even*),range(*hopping range*),ξ(*Scale of exp decaying*)}={4,2,2};
makehop[];
{ϕ,η,c}={-1,1,1};U=KroneckerProduct[IdentityMatrix[n/4],PauliMatrix[0],PauliMatrix[2]];
symmetrizehop[ϕ,η,c];
{ϕ,η,c}={-1,1,-1};U=KroneckerProduct[IdentityMatrix[n/4],PauliMatrix[2],PauliMatrix[0]];
symmetrizehop[ϕ,η,c];
Lx(*Site length*)=100;
spectra[hop,Lx]

(29) Pseudo Hermiticity + Sublattice Symmetry ((Subscript[u, 1])(Subscript[u, 2])=(Subscript[u, 2])(Subscript[u, 1]))

{n(*matrix size, should be even*),range(*hopping range*),ξ(*Scale of exp decaying*)}={4,2,2};
makehop[];
{ϕ,η,c}={1,-1,1};U=KroneckerProduct[IdentityMatrix[n/4],PauliMatrix[0],PauliMatrix[3]];
symmetrizehop[ϕ,η,c];
{ϕ,η,c}={1,1,-1};U=KroneckerProduct[IdentityMatrix[n/4],PauliMatrix[3],PauliMatrix[0]];
symmetrizehop[ϕ,η,c];
Lx(*Site length*)=100;
spectra[hop,Lx]

(30) Pseudo Hermiticity + Sublattice Symmetry ((Subscript[u, 1])(Subscript[u, 2])=-(Subscript[u, 2])(Subscript[u, 1]))

{n(*matrix size, should be even*),range(*hopping range*),ξ(*Scale of exp decaying*)}={4,2,2};
makehop[];
{ϕ,η,c}={1,-1,1};U=KroneckerProduct[IdentityMatrix[n/4],PauliMatrix[1],PauliMatrix[0]];
symmetrizehop[ϕ,η,c];
{ϕ,η,c}={1,1,-1};U=KroneckerProduct[IdentityMatrix[n/4],PauliMatrix[3],PauliMatrix[0]];
symmetrizehop[ϕ,η,c];
Lx(*Site length*)=100;
spectra[hop,Lx]

N = 3 (n∈8Z)

(31) TRS+TRS†+PHS† ((Subscript[u, 1])(SuperStar[(Subscript[u, 1])])=1, (Subscript[u, 2])(SuperStar[(Subscript[u, 2])])=1, (Subscript[u, 3])(SuperStar[(Subscript[u, 3])])=1)

{n(*matrix size, should be even*),range(*hopping range*),ξ(*Scale of exp decaying*)}={8,2,2};
makehop[];
{ϕ,η,c}={-1,1,1};U=KroneckerProduct[IdentityMatrix[n/8],PauliMatrix[3],PauliMatrix[0],PauliMatrix[0]];
symmetrizehop[ϕ,η,c];
{ϕ,η,c}={-1,-1,1};U=KroneckerProduct[IdentityMatrix[n/8],PauliMatrix[0],PauliMatrix[3],PauliMatrix[0]];
symmetrizehop[ϕ,η,c];
{ϕ,η,c}={-1,1,-1};U=KroneckerProduct[IdentityMatrix[n/8],PauliMatrix[0],PauliMatrix[0],PauliMatrix[3]];
symmetrizehop[ϕ,η,c];
Lx(*Site length*)=100;
spectrawithEdge[hop,Lx]

(32) TRS+TRS†+PHS† ((Subscript[u, 1])(SuperStar[(Subscript[u, 1])])=-1, (Subscript[u, 2])(SuperStar[(Subscript[u, 2])])=1, (Subscript[u, 3])(SuperStar[(Subscript[u, 3])])=1)

{n(*matrix size, should be even*),range(*hopping range*),ξ(*Scale of exp decaying*)}={8,2,2};
makehop[];
{ϕ,η,c}={-1,1,1};U=KroneckerProduct[IdentityMatrix[n/8],PauliMatrix[2],PauliMatrix[0],PauliMatrix[0]];
symmetrizehop[ϕ,η,c];
{ϕ,η,c}={-1,-1,1};U=KroneckerProduct[IdentityMatrix[n/8],PauliMatrix[0],PauliMatrix[3],PauliMatrix[0]];
symmetrizehop[ϕ,η,c];
{ϕ,η,c}={-1,1,-1};U=KroneckerProduct[IdentityMatrix[n/8],PauliMatrix[0],PauliMatrix[0],PauliMatrix[3]];
symmetrizehop[ϕ,η,c];
Lx(*Site length*)=100;
spectrawithEdge[hop,Lx]

(33) TRS+TRS†+PHS† ((Subscript[u, 1])(SuperStar[(Subscript[u, 1])])=1, (Subscript[u, 2])(SuperStar[(Subscript[u, 2])])=-1, (Subscript[u, 3])(SuperStar[(Subscript[u, 3])])=1)

{n(*matrix size, should be even*),range(*hopping range*),ξ(*Scale of exp decaying*)}={8,2,2};
makehop[];
{ϕ,η,c}={-1,1,1};U=KroneckerProduct[IdentityMatrix[n/8],PauliMatrix[3],PauliMatrix[0],PauliMatrix[0]];
symmetrizehop[ϕ,η,c];
{ϕ,η,c}={-1,-1,1};U=KroneckerProduct[IdentityMatrix[n/8],PauliMatrix[0],PauliMatrix[2],PauliMatrix[0]];
symmetrizehop[ϕ,η,c];
{ϕ,η,c}={-1,1,-1};U=KroneckerProduct[IdentityMatrix[n/8],PauliMatrix[0],PauliMatrix[0],PauliMatrix[3]];
symmetrizehop[ϕ,η,c];
Lx(*Site length*)=100;
spectra[hop,Lx]

(34) TRS+TRS†+PHS† ((Subscript[u, 1])(SuperStar[(Subscript[u, 1])])=1, (Subscript[u, 2])(SuperStar[(Subscript[u, 2])])=1, (Subscript[u, 3])(SuperStar[(Subscript[u, 3])])=-1)

{n(*matrix size, should be even*),range(*hopping range*),ξ(*Scale of exp decaying*)}={8,2,2};
makehop[];
{ϕ,η,c}={-1,1,1};U=KroneckerProduct[IdentityMatrix[n/8],PauliMatrix[3],PauliMatrix[0],PauliMatrix[0]];
symmetrizehop[ϕ,η,c];
{ϕ,η,c}={-1,-1,1};U=KroneckerProduct[IdentityMatrix[n/8],PauliMatrix[0],PauliMatrix[3],PauliMatrix[0]];
symmetrizehop[ϕ,η,c];
{ϕ,η,c}={-1,1,-1};U=KroneckerProduct[IdentityMatrix[n/8],PauliMatrix[0],PauliMatrix[0],PauliMatrix[2]];
symmetrizehop[ϕ,η,c];
Lx(*Site length*)=100;
spectrawithEdge[hop,Lx]

(35) TRS+TRS†+PHS† ((Subscript[u, 1])(SuperStar[(Subscript[u, 1])])=-1, (Subscript[u, 2])(SuperStar[(Subscript[u, 2])])=-1, (Subscript[u, 3])(SuperStar[(Subscript[u, 3])])=1)

{n(*matrix size, should be even*),range(*hopping range*),ξ(*Scale of exp decaying*)}={8,2,2};
makehop[];
{ϕ,η,c}={-1,1,1};U=KroneckerProduct[IdentityMatrix[n/8],PauliMatrix[2],PauliMatrix[0],PauliMatrix[0]];
symmetrizehop[ϕ,η,c];
{ϕ,η,c}={-1,-1,1};U=KroneckerProduct[IdentityMatrix[n/8],PauliMatrix[0],PauliMatrix[2],PauliMatrix[0]];
symmetrizehop[ϕ,η,c];
{ϕ,η,c}={-1,1,-1};U=KroneckerProduct[IdentityMatrix[n/8],PauliMatrix[0],PauliMatrix[0],PauliMatrix[3]];
symmetrizehop[ϕ,η,c];
Lx(*Site length*)=100;
spectra[hop,Lx]

(36) TRS+TRS†+PHS† ((Subscript[u, 1])(SuperStar[(Subscript[u, 1])])=-1, (Subscript[u, 2])(SuperStar[(Subscript[u, 2])])=1, (Subscript[u, 3])(SuperStar[(Subscript[u, 3])])=-1)

{n(*matrix size, should be even*),range(*hopping range*),ξ(*Scale of exp decaying*)}={8,2,2};
makehop[];
{ϕ,η,c}={-1,1,1};U=KroneckerProduct[IdentityMatrix[n/8],PauliMatrix[2],PauliMatrix[0],PauliMatrix[0]];
symmetrizehop[ϕ,η,c];
{ϕ,η,c}={-1,-1,1};U=KroneckerProduct[IdentityMatrix[n/8],PauliMatrix[0],PauliMatrix[3],PauliMatrix[0]];
symmetrizehop[ϕ,η,c];
{ϕ,η,c}={-1,1,-1};U=KroneckerProduct[IdentityMatrix[n/8],PauliMatrix[0],PauliMatrix[0],PauliMatrix[2]];
symmetrizehop[ϕ,η,c];
Lx(*Site length*)=100;
spectrawithEdge[hop,Lx]

(37) TRS+TRS†+PHS† ((Subscript[u, 1])(SuperStar[(Subscript[u, 1])])=1, (Subscript[u, 2])(SuperStar[(Subscript[u, 2])])=-1, (Subscript[u, 3])(SuperStar[(Subscript[u, 3])])=-1)

{n(*matrix size, should be even*),range(*hopping range*),ξ(*Scale of exp decaying*)}={8,2,2};
makehop[];
{ϕ,η,c}={-1,1,1};U=KroneckerProduct[IdentityMatrix[n/8],PauliMatrix[3],PauliMatrix[0],PauliMatrix[0]];
symmetrizehop[ϕ,η,c];
{ϕ,η,c}={-1,-1,1};U=KroneckerProduct[IdentityMatrix[n/8],PauliMatrix[0],PauliMatrix[2],PauliMatrix[0]];
symmetrizehop[ϕ,η,c];
{ϕ,η,c}={-1,1,-1};U=KroneckerProduct[IdentityMatrix[n/8],PauliMatrix[0],PauliMatrix[0],PauliMatrix[2]];
symmetrizehop[ϕ,η,c];
Lx(*Site length*)=100;
spectra[hop,Lx]

(38) TRS+TRS†+PHS† ((Subscript[u, 1])(SuperStar[(Subscript[u, 1])])=-1, (Subscript[u, 2])(SuperStar[(Subscript[u, 2])])=-1, (Subscript[u, 3])(SuperStar[(Subscript[u, 3])])=-1)

{n(*matrix size, should be even*),range(*hopping range*),ξ(*Scale of exp decaying*)}={8,2,2};
makehop[];
{ϕ,η,c}={-1,1,1};U=KroneckerProduct[IdentityMatrix[n/8],PauliMatrix[2],PauliMatrix[0],PauliMatrix[0]];
symmetrizehop[ϕ,η,c];
{ϕ,η,c}={-1,-1,1};U=KroneckerProduct[IdentityMatrix[n/8],PauliMatrix[0],PauliMatrix[2],PauliMatrix[0]];
symmetrizehop[ϕ,η,c];
{ϕ,η,c}={-1,1,-1};U=KroneckerProduct[IdentityMatrix[n/8],PauliMatrix[0],PauliMatrix[0],PauliMatrix[2]];
symmetrizehop[ϕ,η,c];
Lx(*Site length*)=100;
spectra[hop,Lx]

ハミルトニアンのグループ分け

Ref : Long = Zhang, PRL 130, 036601 (2023)
（摂動に対する安定性は実装していない）

r[n_]:=RandomReal[{-1,1},{n,n}]+I*RandomReal[{-1,1},{n,n}];
makeQ[h_]:=Module[{Λ,U},
{Λ,U}=Eigensystem[h];
Return[U\[Transpose].DiagonalMatrix[Sign[Λ]].U\[Conjugate]];];
δ=10^(-8);
connectQ[h0_,h1_]:=Module[{Q0,Q1},
{Q0,Q1}=makeQ/@{h0,h1};
If[Min[Abs[Eigenvalues[(Q0+Q1)/2]]]>δ,Return[1];,Return[0];]
];
{σx,σy,σz}=Table[PauliMatrix[i],{i,3}];

グラフ・ラプラシアンを用いる方法 （効率悪い）

n=6;
Nsample=100;
hs=Table[h=r[n];h=1/2*(h+h\[ConjugateTranspose]),Nsample];
Mcon(*断熱的につながれば1,繋がらなければ0*)=Table[connectQ[hs[[i]],hs[[j]]],{i,Nsample},{j,Nsample}];
deg=DiagonalMatrix[Total/@Mcon-1];
adj=Mcon-IdentityMatrix[Nsample];
L(*グラフラプラシアン*) =deg-adj;
NullSpace[L]
AdjacencyGraph[adj]

順番にグループ分けする方法
断熱変形可能性は同値関係なので，代表ハミルトニアンの対で判定すれば良い．
断熱的に繋がらないハミルトニアンが見つかればラベルを増やす.

Class A

n=4;
Nsample=20000;
hs=Table[h=r[n];h=1/2*(h+h\[ConjugateTranspose]),Nsample];
class={};
Do[
h0=hs[[1]];
tmp=Reap[Do[If[connectQ[h0,h1]==0(*繋がらない*),Sow[h1]],{h1,hs[[2;;]]}];];
If[Length[tmp[[2]]]==0,AppendTo[class,Length[hs]];Break[];,AppendTo[class,Length[hs]-Length[tmp[[2,1]]]];];
hs=tmp[[2,1]];
,{re,100}];
Print[class,Total[class]]
BarChart[class]

Class AIII

n=4;
Nsample=10000;
Γ=KroneckerProduct[IdentityMatrix[n/2],σz];
hs=Table[h=r[n];h=1/2*(h+h\[ConjugateTranspose]);h=1/2*(h-Γ.h.Γ),Nsample];
class={};
Do[
h0=hs[[1]];
tmp=Reap[Do[If[connectQ[h0,h1]==0(*繋がらない*),Sow[h1]],{h1,hs[[2;;]]}];];
If[Length[tmp[[2]]]==0,AppendTo[class,Length[hs]];Break[];,AppendTo[class,Length[hs]-Length[tmp[[2,1]]]];];
hs=tmp[[2,1]];
,{re,100}];
Print[class,Total[class]]
BarChart[class]

Class AI

n=4;
Nsample=10000;
hs=Table[h=r[n];h=1/2*(h+h\[ConjugateTranspose]);h=1/2*(h+h\[Conjugate]),Nsample];
class={};
Do[
h0=hs[[1]];
tmp=Reap[Do[If[connectQ[h0,h1]==0(*繋がらない*),Sow[h1]],{h1,hs[[2;;]]}];];
If[Length[tmp[[2]]]==0,AppendTo[class,Length[hs]];Break[];,AppendTo[class,Length[hs]-Length[tmp[[2,1]]]];];
hs=tmp[[2,1]];
,{re,100}];
Print[class,Total[class]]
BarChart[class]

Class BDI

n=4;
Nsample=10000;
uc=KroneckerProduct[IdentityMatrix[n/2],σz];
hs=Table[h=r[n];h=1/2*(h+h\[ConjugateTranspose]);h=1/2*(h+h\[Conjugate]);h=1/2*(h-uc.h\[Conjugate].uc\[ConjugateTranspose]),Nsample];
class={};
Do[
h0=hs[[1]];
tmp=Reap[Do[If[connectQ[h0,h1]==0(*繋がらない*),Sow[h1]],{h1,hs[[2;;]]}];];
If[Length[tmp[[2]]]==0,AppendTo[class,Length[hs]];Break[];,AppendTo[class,Length[hs]-Length[tmp[[2,1]]]];];
hs=tmp[[2,1]];
,{re,100}];
Print[class,Total[class]]
BarChart[class]

Class D

n=4;
Nsample=10000;
hs=Table[h=r[n];h=1/2*(h+h\[ConjugateTranspose]);h=1/2*(h-h\[Transpose]),Nsample];
class={};
Do[
h0=hs[[1]];
tmp=Reap[Do[If[connectQ[h0,h1]==0(*繋がらない*),Sow[h1]],{h1,hs[[2;;]]}];];
If[Length[tmp[[2]]]==0,AppendTo[class,Length[hs]];Break[];,AppendTo[class,Length[hs]-Length[tmp[[2,1]]]];];
hs=tmp[[2,1]];
,{re,100}];
Print[class,Total[class]]
BarChart[class]

Class DIII

n=4;
Nsample=10000;
ut=KroneckerProduct[IdentityMatrix[n/2],σy];
hs=Table[h=r[n];h=1/2*(h+h\[ConjugateTranspose]);h=1/2*(h+ut.h\[Conjugate].ut\[ConjugateTranspose]);h=1/2*(h-h\[Conjugate]),Nsample];
class={};
Do[
h0=hs[[1]];
tmp=Reap[Do[If[connectQ[h0,h1]==0(*繋がらない*),Sow[h1]],{h1,hs[[2;;]]}];];
If[Length[tmp[[2]]]==0,AppendTo[class,Length[hs]];Break[];,AppendTo[class,Length[hs]-Length[tmp[[2,1]]]];];
hs=tmp[[2,1]];
,{re,100}];
Print[class,Total[class]]
BarChart[class]

Class AII TRS

n=4;
Nsample=10000;
u=KroneckerProduct[IdentityMatrix[n/2],PauliMatrix[2]];
hs=Table[h=r[n];h=1/2*(h+h\[ConjugateTranspose]);h=1/2*(h+u.h\[Conjugate].u\[ConjugateTranspose]),Nsample];
class={};
Do[
h0=hs[[1]];
tmp=Reap[Do[If[connectQ[h0,h1]==0(*繋がらない*),Sow[h1]],{h1,hs[[2;;]]}];];
If[Length[tmp[[2]]]==0,AppendTo[class,Length[hs]];Break[];,AppendTo[class,Length[hs]-Length[tmp[[2,1]]]];];
hs=tmp[[2,1]];
,{re,100}];
Print[class,Total[class]]
BarChart[class]

Class CII

n=8;
Nsample=10000;
ut=KroneckerProduct[IdentityMatrix[n/4],σy,IdentityMatrix[2]];
uc=KroneckerProduct[IdentityMatrix[n/4],IdentityMatrix[2],σy];
hs=Table[h=r[n];h=1/2*(h+h\[ConjugateTranspose]);h=1/2*(h+ut.h\[Conjugate].ut\[ConjugateTranspose]);h=1/2*(h-uc.h\[Conjugate].uc\[ConjugateTranspose]),Nsample];
class={};
Do[
h0=hs[[1]];
tmp=Reap[Do[If[connectQ[h0,h1]==0(*繋がらない*),Sow[h1]],{h1,hs[[2;;]]}];];
If[Length[tmp[[2]]]==0,AppendTo[class,Length[hs]];Break[];,AppendTo[class,Length[hs]-Length[tmp[[2,1]]]];];
hs=tmp[[2,1]];
,{re,100}];
Print[class,Total[class]]
BarChart[class]

Class C PHS

n=4;
Nsample=10000;
u=KroneckerProduct[IdentityMatrix[n/2],PauliMatrix[2]];
hs=Table[h=r[n];h=1/2*(h+h\[ConjugateTranspose]);h=1/2*(h-u.h\[Transpose].u\[ConjugateTranspose]),Nsample];
class={};
Do[
h0=hs[[1]];
tmp=Reap[Do[If[connectQ[h0,h1]==0(*繋がらない*),Sow[h1]],{h1,hs[[2;;]]}];];
If[Length[tmp[[2]]]==0,AppendTo[class,Length[hs]];Break[];,AppendTo[class,Length[hs]-Length[tmp[[2,1]]]];];
hs=tmp[[2,1]];
,{re,100}];
Print[class,Total[class]]
BarChart[class]

Class CI

n=4;
Nsample=10000;
uc=KroneckerProduct[IdentityMatrix[n/2],σy];
hs=Table[h=r[n];h=1/2*(h+h\[ConjugateTranspose]);h=1/2*(h+h\[Conjugate]);h=1/2*(h-uc.h\[Conjugate].uc\[ConjugateTranspose]),Nsample];
class={};
Do[
h0=hs[[1]];
tmp=Reap[Do[If[connectQ[h0,h1]==0(*繋がらない*),Sow[h1]],{h1,hs[[2;;]]}];];
If[Length[tmp[[2]]]==0,AppendTo[class,Length[hs]];Break[];,AppendTo[class,Length[hs]-Length[tmp[[2,1]]]];];
hs=tmp[[2,1]];
,{re,100}];
Print[class,Total[class]]
BarChart[class]

2 バンド模型

r[n_]:=RandomReal[{-1,1},{n,n}]+I*RandomReal[{-1,1},{n,n}];
makeQ[h_]:=Module[{Λ,U},
{Λ,U}=Eigensystem[h];
Return[U\[Transpose].DiagonalMatrix[Sign[Λ]].U\[Conjugate]];];

0 次元

h=r[2];
h=1/2*(h+h\[ConjugateTranspose]);
h=h-1/2*Tr[h];
sgnh=makeQ[h];
n=Table[Tr[PauliMatrix[i].sgnh]/2,{i,3}]//Chop;
Show[Graphics3D[
{{Opacity[0.2],Sphere[{0,0,0}]},(*透明度と球の指定*)
{Arrow[{{0,0,0},{1.2,0,0}}]},Text["x",{1.3,0,0}],
{Arrow[{{0,0,0},{0,1.2,0}}]},Text["y",{0,1.3,0}],
{Arrow[{{0,0,0},{0,0,1.2}}]} ,Text["z",{0,0,1.3}]}
],
ListPointPlot3D[{n},PlotStyle->PointSize[0.03]],Boxed->False (*立方体の枠を消す*)]

１次元

(*Hamiltonian*)
Clear[kx,x,hop];
$Assumptions={kx,x}∈Reals;
n(*行列サイズ*)=2;
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=1/2*(h+h\[ConjugateTranspose])//Simplify//Chop(*エルミート化*);*)
Plot[Eigenvalues[h+100*IdentityMatrix[n]]-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]))"}]*)
num=1000;
ns=Table[
sgnh=makeQ[h];
Table[Tr[PauliMatrix[i].sgnh]/2,{i,3}]//Chop,{kx,-Pi,Pi,2 Pi/num}];
tmp=Table[
Show[Graphics3D[
{{Opacity[0.2],Sphere[{0,0,0}]},(*透明度と球の指定*)
{Arrow[{{0,0,0},{1.2,0,0}}]},Text["x",{1.3,0,0}],
{Arrow[{{0,0,0},{0,1.2,0}}]},Text["y",{0,1.3,0}],
{Arrow[{{0,0,0},{0,0,1.2}}]} ,Text["z",{0,0,1.3}]}
],
ListPointPlot3D[ns[[;;i]],PlotStyle->PointSize[0.02]],Boxed->False (*立方体の枠を消す*)]
,{i,1,num+1,50}];
ListAnimate[tmp]

対称性を課す
ex: class D PHS: 
(Subscript[σ, x])H(k)* (Subscript[σ, x])=-H(-k)
->k=0,piは北極 or 南極に

(*Hamiltonian*)
Clear[kx,x,hop];
$Assumptions={kx,x}∈Reals;
n(*行列サイズ*)=2;
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=1/2*(h+h\[ConjugateTranspose])//Simplify//Chop(*エルミート化*);*)
(*class D PHS*)
σx=PauliMatrix[1];
h=1/2*(h-σx.(h/.{kx->-kx})\[Conjugate].σx);
Plot[Eigenvalues[h+100*IdentityMatrix[n]]-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]))"}]*)
num=1000;
ns=Table[
sgnh=makeQ[h];
Table[Tr[PauliMatrix[i].sgnh]/2,{i,3}]//Chop,{kx,-Pi,Pi,2 Pi/num}];
tmp=Table[
Show[Graphics3D[
{{Opacity[0.2],Sphere[{0,0,0}]},(*透明度と球の指定*)
{Arrow[{{0,0,0},{1.2,0,0}}]},Text["x",{1.3,0,0}],
{Arrow[{{0,0,0},{0,1.2,0}}]},Text["y",{0,1.3,0}],
{Arrow[{{0,0,0},{0,0,1.2}}]} ,Text["z",{0,0,1.3}]}
],
ListPointPlot3D[ns[[;;i]],PlotStyle->PointSize[0.02]],Boxed->False (*立方体の枠を消す*)]
,{i,1,num/2+1,50}];
ListAnimate[tmp]

２ 次元

ClearAll[kx,ky,x,y,ihop,hop];
$Assumptions={kx,ky,x,y}∈Reals;
matsize(*行列サイズ*)=2;
range(*ホッピングのレンジ*)=2;
ξ(*ホッピングの減衰スケール*)=1;
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+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]))"}]

(*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];
Print["Chern数=",SumFlux/(2*Pi)//Chop];

h=h-100*IdentityMatrix[matsize](*定数シフトを戻す*);
num=50;
ns=Table[
sgnh=makeQ[h];
Table[Tr[PauliMatrix[i].sgnh]/2,{i,3}]//Chop,{kx,-Pi,Pi,2 Pi/num},{ky,-Pi,Pi,2 Pi/(num*5)}];
tmp=Table[
Show[Graphics3D[
{{Opacity[0.2],Sphere[{0,0,0}]},(*透明度と球の指定*)
{Arrow[{{0,0,0},{1.2,0,0}}]},Text["x",{1.3,0,0}],
{Arrow[{{0,0,0},{0,1.2,0}}]},Text["y",{0,1.3,0}],
{Arrow[{{0,0,0},{0,0,1.2}}]} ,Text["z",{0,0,1.3}]}
],
ListPointPlot3D[ns[[;;i]],PlotStyle->PointSize[0.02]],Boxed->False (*立方体の枠を消す*)]
,{i,1,num}];
ListAnimate[tmp]

ハミルトニアンのブロック対角化

gtab=FiniteGroupData[{"SymmetricGroup",4},"MultiplicationTable"];
gtab//MatrixForm
n=Length[gtab];
(*正則表現*)
reg=Table[KroneckerDelta[h,gtab[[g,k]]],{g,n},{h,n},{k,n}];
h=RandomReal[{-1,1},{n,n}]+I*RandomReal[{-1,1},{n,n}];
(*対称化*)
h=1/n*Sum[reg[[g]].h.reg[[g]]\[ConjugateTranspose],{g,n}];
{Λ,U}=Eigensystem[h];
U=U\[Transpose];
Print[Λ];

data={Re[#],Im[#]}&/@Λ;

(*固有値の縮退度を計算*)
degeneracies=Tally[Round[data,0.0001]];

(*プロット用のデータとテキストラベルを生成*)
plotData={};
textLabels={};

For[i=1,i≤Length[degeneracies],i++,AppendTo[plotData,First[degeneracies[[i]]]];
AppendTo[textLabels,Text[ToString[Last[degeneracies[[i]]]],First[degeneracies[[i]]],{0,-2}]];]

(*点とテキストを同じグラフ上に表示*)
Show[ListPlot[plotData,PlotStyle->PointSize[0.02],PlotRange->All],Graphics[textLabels],PlotRangePadding->{0.5,0.5} (*これを追加*)]