Operators

Pauli matrices

id=(({{1, 0}, {0, 1}}));σx=(({{0, 1}, {1, 0}})); σy=(({{0, -I}, {I, 0}}));σz=(({{1, 0}, {0, -1}}));σp=(({{0, 1}, {0, 0}}));σm=(({{0, 0}, {1, 0}}));

Fermion operators

(* a: site index,
   b: spin index,
   c: orbital index,
   L: # possible single-electron states,
  Mc: annihilation op.,
  Mcdag: creation op.
*)
Mc[a_,b_,c_,L_]:= Module[{int},
int=FromDigits[{a,b,c},2]+1;
Apply[KroneckerProduct,Join[Table[-σz,{l,1,int-1}],{σp},Table[id,{l,int+1,L}]]]
]//SparseArray;
Mcdag[a_,b_,c_,L_]:=Transpose[Mc[a,b,c,L]]

2-site Hubbard model

dim=8;
(*number op.*)
Nmat=Sum[Mcdag[a,b,c,dim].Mc[a,b,c,dim],{a,0,1},{b,0,1},{c,0,1}];
(*Hamiltonian, 
  t: hopping amplitude,
  U: on-site Coulomb int.,
 n0: el. number at site 0,
 n1: el. number at site 1
*)
Ham=Module[{n0,n1,Umat,pre,Tmat},
n0=Sum[Mcdag[0,b,c,dim].Mc[0,b,c,dim],{b,0,1},{c,0,1}];
n1=Sum[Mcdag[1,b,c,dim].Mc[1,b,c,dim],{b,0,1},{c,0,1}];
Umat=U (n0.(n0-IdentityMatrix[2^dim])/2+n1.(n1-IdentityMatrix[2^dim])/2)//SparseArray;
pre=-t Mcdag[1,0,0,dim].Mc[0,0,0,dim]-t Mcdag[1,0,1,dim].Mc[0,0,1,dim]-t Mcdag[1,1,0,dim].Mc[0,1,0,dim]-t Mcdag[1,1,1,dim].Mc[0,1,1,dim];
Tmat=Simplify[pre+ConjugateTranspose[pre],Element[t,Reals]]//SparseArray;
Umat+Tmat
];

(*In the following, we consider the subspace where the total number of electrons is 2.
The lowest-energy manifold in this subspace is labeled by*)
loweng=Module[{el2,low},
el2=Position[Table[Nmat[[i,i]],{i,1,2^dim}],2];
low=Position[Table[Ham[[i,i]],{i,1,2^dim}],0];
Intersection[el2,low]//Flatten
];
(*The second-lowest manifold in the subspace is labeled by*)
seceng=Module[{el2,sec},
el2=Position[Table[Nmat[[i,i]],{i,1,2^dim}],2];
sec=Position[Table[Ham[[i,i]],{i,1,2^dim}],U];
Intersection[el2,sec]//Flatten
];
(*Then the effective hamiltonian for the low-energy manifold is obtained as*)
effHam=Module[{tm},
tm=Part[Ham,loweng,seceng];
Simplify[-tm.ConjugateTranspose[tm]/U,Element[t,Reals]]
];

effHam//MatrixForm

SU(4) Kugel-Khomskii Hamiltonian

(*Kugel-Khomskii Hamiltonian PLUS some constant*)
sig[0]:=id;
sig[1]:={{0,1},{1,0}};
sig[2]:={{0,-I},{I,0}};
sig[3]:={{1,0},{0,-1}};
kro[A__]:=KroneckerProduct[A]
HamKK=(t^2/(2 U))Sum[kro[sig[a],sig[0],sig[a],sig[0]],{a,0,3}].Sum[kro[sig[0],sig[b],sig[0],sig[b]],{b,0,3}]-(2 t^2/U) IdentityMatrix[16];

HamKK-effHam//MatrixForm

How to read off the coefficients using Hilbert-Schmidt inner product

coe1[a_,b_,c_,d_]:=Tr[effHam.ConjugateTranspose[kro[sig[a],sig[b],sig[c],sig[d]]]]/2^4//FullSimplify
ind[0]:=0; ind[1]:=x; ind[2]:=y; ind[3]:=z;
Sum[coe1[a,b,c,d] ((σ1)^(ind[a])) ((τ1)^(ind[b]))((σ2)^(ind[c])) ((τ2)^(ind[d])),{a,0,3},{b,0,3},{c,0,3},{d,0,3}]

coe2[a_,b_,c_,d_]:=Tr[HamKK.ConjugateTranspose[kro[sig[a],sig[b],sig[c],sig[d]]]]/2^4//FullSimplify
ind[0]:=0; ind[1]:=x; ind[2]:=y; ind[3]:=z;
Sum[coe2[a,b,c,d] ((σ1)^(ind[a])) ((τ1)^(ind[b]))((σ2)^(ind[c])) ((τ2)^(ind[d])),{a,0,3},{b,0,3},{c,0,3},{d,0,3}]