Hamiltonian of the open XXZ chain

(*Pauli matrices*)
sig[0]={{1,0},{0,1}}//N//SparseArray;
sig[1]={{0,1},{1,0}}//N//SparseArray;
sig[2]={{0,-I},{I,0}}//N//SparseArray;
sig[3]={{1,0},{0,-1}}//N//SparseArray;
(*Hamiltonia, "leng" is the chain length*)
Ham[Δ_,leng_]:=Module[{Sig},
Sig[a_,i_]:=Sig[a,i]=Apply[KroneckerProduct,Join[Table[sig[0],{l,1,i-1}],{sig[a]},Table[sig[0],{l,i+1,leng}]]];
Sum[Sum[Sig[a,i].Sig[a,i+1],{a,1,2}]+Δ Sig[3,i].Sig[3,i+1],{i,1,leng-1}
]
]

Entanglement spectrum

(*Subsystem sizes are the same (leng/2). "leng" should be even.*)
EESpec[Δ_,leng_]:=Module[{gsvec,SchMat,pre},
gsvec=Sort[Eigensystem[Ham[Δ,leng]//Normal]//Thread,#1[[1]]<#2[[1]] &][[1]][[2]];
SchMat=Module[{vec,lis,state,dim},
vec[0]={1,0}//N//SparseArray;
vec[1]={0,1}//N//SparseArray;
lis[k_]:=lis[k]=IntegerDigits[k-1,2,leng/2];
state[k_]:=state[k]=Apply[
KroneckerProduct,Table[vec[lis[k][[i]]],{i,1,leng/2}]
]//Flatten;
dim=2^(leng/2);
Table[
Flatten[KroneckerProduct[state[i],state[j]]].gsvec,{i,1,dim},{j,1,dim}
]
];
pre=Eigenvalues[(ConjugateTranspose[SchMat].SchMat)//Normal];
Table[-Log[pre[[k]]],{k,1,Length[pre]}]
]

EESpec[1.0,4]
EESpec[1.0,6]
EESpec[1.0,8]
EESpec[1.0,10]

EESpec[1.0,12]//ListPlot

Hamiltonian of the periodic XXZ chain

(*Pauli matrices*)
sig[0]={{1,0},{0,1}}//N//SparseArray;
sig[1]={{0,1},{1,0}}//N//SparseArray;
sig[2]={{0,-I},{I,0}}//N//SparseArray;
sig[3]={{1,0},{0,-1}}//N//SparseArray;
(*Hamiltonia, "leng" is the chain length*)
PHam[Δ_,leng_]:=Module[{Sig},
Sig[a_,i_]:=Sig[a,i]=Apply[KroneckerProduct,Join[Table[sig[0],{l,1,i-1}],{sig[a]},Table[sig[0],{l,i+1,leng}]]];
Sum[Sum[Sig[a,i].Sig[a,i+1],{a,1,2}]+Δ Sig[3,i].Sig[3,i+1],{i,1,leng-1}
]+(Sum[Sig[a,1].Sig[a,leng],{a,1,2}]+Δ Sig[3,1].Sig[3,leng])
]

Entanglement spectrum

(*Subsystem sizes are the same (leng/2). "leng" should be even.*)
PEESpec[Δ_,leng_]:=Module[{gsvec,SchMat,pre},
gsvec=Sort[Eigensystem[PHam[Δ,leng]//Normal]//Thread,#1[[1]]<#2[[1]] &][[1]][[2]];
SchMat=Module[{vec,lis,state,dim},
vec[0]={1,0}//N//SparseArray;
vec[1]={0,1}//N//SparseArray;
lis[k_]:=lis[k]=IntegerDigits[k-1,2,leng/2];
state[k_]:=state[k]=Apply[
KroneckerProduct,Table[vec[lis[k][[i]]],{i,1,leng/2}]
]//Flatten;
dim=2^(leng/2);
Table[
Flatten[KroneckerProduct[state[i],state[j]]].gsvec,{i,1,dim},{j,1,dim}
]
];
pre=Eigenvalues[(ConjugateTranspose[SchMat].SchMat)//Normal];
Table[-Log[pre[[k]]],{k,1,Length[pre]}]
]

PEESpec[1.0,4]
PEESpec[1.0,6]
PEESpec[1.0,8]
PEESpec[1.0,10]

PEESpec[1.0,12]//ListPlot