Fermion operators

Basic operators

σ0=(({{1, 0}, {0, 1}}))//SparseArray;
σp=(({{0, 1}, {0, 0}}))//SparseArray;
σm=(({{0, 0}, {1, 0}}))//SparseArray;
σ3=(({{1, 0}, {0, -1}}))//SparseArray;

Spinless fermion and number operators

(*annihilation operators*)
cop[i_,M_]:=Apply[KroneckerProduct,Join[Table[σ3,{l,1,i-1}],{σp},Table[σ0,{l,i+1,M}]]]//SparseArray;
(*creation operators*)
cdg[i_,M_]:=Transpose[cop[i,M]];
(*number operators*)
nop[i_,M_]:=cdg[i,M].cop[i,M];

(*Define anticommutator between a and b*)
AC[a_,b_]:=a.b+b.a;
(*Check the anticommutation relations*)
M=4;
Table[
{{i,j},If[
AC[cop[i,M],cdg[j,M]]==0 IdentityMatrix[2^M],
0,TrueQ[AC[cop[i,M],cdg[j,M]]==IdentityMatrix[2^M]]
]},{i,1,M},{j,1,M}
]

Fock space

v[0]={1,0};
v[1]={0,1};
(*vacuum state*)
vac=Apply[KroneckerProduct,Table[v[0],{i,1,M}]]//Flatten
cop[1,M].vac
(*1-fermion state*)
cdg[1,M].vac
(*normalization*)
vac.cop[1,M].cdg[1,M].vac
(*2-fermion states*)
cdg[1,M].cdg[3,M].vac
(*normalization*)
vac.cop[3,M].cop[1,M].cdg[1,M].cdg[3,M].vac

Majorana fermions

γ1[i_,M_]:=If[Mod[i,2]==0,cop[i,M]+cdg[i,M],-I (cop[i,M]-cdg[i,M])];
γ2[i_,M_]:=If[Mod[i,2]==0,-I (cop[i,M]-cdg[i,M]),cop[i,M]+cdg[i,M]];

(*Check Hermiticity*)
M=4;
Table[
If[
ConjugateTranspose[γ1[i,M]]-γ1[i,M]==0 IdentityMatrix[2^M],
True,False
],
{i,1,M}
]
Table[
If[
ConjugateTranspose[γ2[i,M]]-γ2[i,M]==0 IdentityMatrix[2^M],
True,False
],
{i,1,M}
]

(*Check anticommutation relations*)
Table[
{{i,j},If[
AC[γ1[i,M],γ1[j,M]]==0 IdentityMatrix[2^M],
0,TrueQ[AC[γ1[i,M],γ1[j,M]]==2 IdentityMatrix[2^M]]
]},{i,1,M},{j,1,M}
]
Table[
{{i,j},If[
AC[γ2[i,M],γ2[j,M]]==0 IdentityMatrix[2^M],
0,TrueQ[AC[γ1[i,M],γ1[j,M]]==2 IdentityMatrix[2^M]]
]},{i,1,M},{j,1,M}
]
Table[
{{i,j},If[
AC[γ1[i,M],γ2[j,M]]==0 IdentityMatrix[2^M],
0,TrueQ[AC[γ1[i,M],γ1[j,M]]==2 IdentityMatrix[2^M]]
]},{i,1,M},{j,1,M}
]

Tight-binding Hamiltonian

(*2-site problem*)
M=2;
Hhop=t (cdg[1,M].cop[2,M]+cdg[2,M].cop[1,M]);
Hhop//Eigenvalues
Eigensystem[Hhop]//Thread//Sort

(*10-site chain with random hopping*)
M=10;
Hhop=Sum[RandomReal[{-1,1}] (cdg[i,M].cop[i+1,M]+cdg[i+1,M].cop[i,M]),{i,1,M-1}];
Eigenvalues[Hhop//Normal]//Sort;
(*Distribution of the energy eigenvalues*)
Histogram[%]