Stationary states of n-TASEP (PBC)

leng=4(*chain length*);
n=3(*# species*);
v[α_]:=Join[Table[0,{i,0,α-1}],{1},Table[0,{i,α+1,n}]](*fundamental vectors*)
kro[A__]:=KroneckerProduct[A]
Em[α_,β_]:=kro[v[α],v[β]]//SparseArray
id:=IdentityMatrix[n+1]//SparseArray;

(*Sector spanned by permutations of {a1, a2, ...}*)
Sector[A_]:=(FromDigits[#,n+1]+1)&/@Permutations[A]

Sector[{0,1,2,3}]

(*Markov Matrix*)
Mark=Module[{ham,bterm},
ham[1]=id;
ham[2]=Sum[Sum[-kro[Em[β,β],Em[α,α]]+kro[Em[α,β],Em[β,α]],{α,0,β-1}],{β,0,n}]//SparseArray;
ham[k_]:=ham[k]=kro[ham[1],ham[k-1]]+kro[ham[2],SparseArray[{{i_, i_}->1}, {(n+1)^(k-2),(n+1)^(k-2)}]];
bterm=Sum[Sum[-kro[Em[α,α],SparseArray[{{i_, i_}->1}, {(n+1)^(leng-2),(n+1)^(leng-2)}],Em[β,β]]+kro[Em[β,α],SparseArray[{{i_, i_}->1}, {(n+1)^(leng-2),(n+1)^(leng-2)}],Em[α,β]],{α,0,β-1}],{β,0,n}];
ham[leng]+bterm
];

(*Specify the sector.*)
sec=Sector[{0,1,2,3}];
(*Get the coefficient of each configuration.*)
SecMark=Part[Mark,sec,sec];
clist=Module[{dim,coef,zero,pre},
dim=Length[sec];
coef=Join[{1},Table[c[k],{k,2,dim}]];
zero=Table[0,{k,1,dim}];
pre=coef/.Solve[SecMark.coef==zero];
pre/Min[pre]//Flatten
];
Table[{clist[[k]],IntegerDigits[sec[[k]]-1,n+1,leng]},{k,1,Length[sec]}]

Table[{clist[[k]],IntegerDigits[sec[[k]]-1,n+1,leng]/.{0->Style[○,FontSize->18] ,1-> Style[●,FontColor->Red,FontSize->18],2->Style[●,FontColor->Green,FontSize->18] ,3->Style[●,FontColor->Blue,FontSize->18] }},{k,1,Length[sec]}]