Supplementary Note for the lecture “Non-Hermitian Topology: Symmetry, Gap conditions, and classification”
Aug. 31, 2023 @National Yang Ming Chiao Tung University

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*)=50;
spectra[hop,Lx]

N = 0 (n∈Z)

(1) No symmetry

{n(*matrix size*),range(*hopping range*),ξ(*Scale of exp decaying*)}={1,2,2};
makehop[];
Lx(*Site length*)=50;
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*)}={4,2,2};
makehop[];
{ϕ,η,c}={1,1,-1};U=KroneckerProduct[IdentityMatrix[n/2],PauliMatrix[3]];
symmetrizehop[ϕ,η,c];
Lx(*Site length*)=100;
spectra[hop,Lx]

(4) Chiral Symmetry

{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]

(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]

(Back up) d=2, Spectra of open boundary condition with x: PBC and y: OBC

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

makehop2D[]:=Module[{},
ClearAll[hop];
Do[hop[nx,ny]=rmat[n]*E^(-Sqrt[nx^2+ny^2]/ξ),{nx,-range,range},{ny,-range,range}];
];

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

(0) Hermitian system

ClearAll[kx,ky,x,y,ihop,hop];
$Assumptions={kx,ky,x,y}∈Reals;
{n(*matrix size*),range(*hopping range*),ξ(*Scale of exp decaying*)}={4,2,2};
makehop2D[];
{ϕ,η,c}={1,-1,1};U=KroneckerProduct[IdentityMatrix[n/4],PauliMatrix[0],PauliMatrix[0]];
symmetrizehop2D[ϕ,η,c];
Ly(*Site length of y*)=20;
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];
Eval=Flatten[Table[Table[{ikx,Re[e],Im[e]},{e,Eigenvalues[hOBC/.{kx->ikx}]}],{ikx,-Pi,Pi,2 Pi/200}],1];
ListPlot[Eval[[All,{2,3}]],PlotRange->All,AxesLabel->{"Re[E]","Im[E]"}]
ListPointPlot3D[Eval,PlotRange->All,AxesLabel->{"(Subscript[k, x])","Re[E]","Im[E]"}]

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

ClearAll[kx,ky,x,y,ihop,hop];
$Assumptions={kx,ky,x,y}∈Reals;
{n(*matrix size*),range(*hopping range*),ξ(*Scale of exp decaying*)}={4,2,2};
makehop2D[];
(*AIII+S- = AIII+η-*)
{ϕ,η,c}={1,-1,1};U=KroneckerProduct[IdentityMatrix[n/4],PauliMatrix[1],PauliMatrix[0]];
symmetrizehop2D[ϕ,η,c];
{ϕ,η,c}={1,1,-1};U=KroneckerProduct[IdentityMatrix[n/4],PauliMatrix[3],PauliMatrix[0]];
symmetrizehop2D[ϕ,η,c];
Ly(*Site length of y*)=20;
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];
Eval=Flatten[Table[Table[{ikx,Re[e],Im[e]},{e,Eigenvalues[hOBC/.{kx->ikx}]}],{ikx,-Pi,Pi,2 Pi/200}],1];
ListPlot[Eval[[All,{2,3}]],PlotRange->All,AxesLabel->{"Re[E]","Im[E]"}]
ListPointPlot3D[Eval,PlotRange->All,AxesLabel->{"(Subscript[k, x])","Re[E]","Im[E]"}]