Lecture 3 - Quantum mechanics of fermions

by Leonid Pryadko

Initialization

Pauli Matrices

(Subscript[σ, x])={{0,1},{1,0}};(Subscript[σ, y])={{0,-I},{I,0}};(Subscript[σ, z])={{1,0},{0,-1}}; e=IdentityMatrix[2];

Direct product

Unprotect[NonCommutativeMultiply];

NonCommutativeMultiply[A_?MatrixQ,B_?MatrixQ]:=KroneckerProduct[A,B];

Spin-1/2 spin matrices

sx=(Subscript[σ, x])/2;sy=(Subscript[σ, y])/2; sz=(Subscript[σ, z])/2;

Boson operators truncated to NN

NN=7;

bp:=Block[{i,j},Table[If[j==i-1,Sqrt[j],0],{i,1,NN},{j,1,NN}]];

bm:=Transpose[bp];

bn:=bp.bm;

be:=IdentityMatrix[NN];

Identical particles: spinless fermions

Non-interacting fermions on a chain

EE[k_]:=-2t Cos[k]/.t->1;

Plot[EE[x],{x,-Pi,Pi}]

EE[2 Pi/4]//N

Here is one way to generate the spectrum of up to 4 non-interacting spinless fermions on an open chain of length nn:

DoFerm[nn_Integer]:=Block[{k,k1,k2,k3,k4,vals={}},
Print["N=0, k=",0," energ=",0//N];
For[k=1,k<=nn,k++,Print["N=1, k=",k," energ=",(EE[ Pi k/(nn+1)])//N]];
For[k=1,k<=nn,k++,vals={};
For[k1=1,k1<=nn,k1++,
k2=k-k1;
If[k2<=0,k2+=nn];
If[k2>k1,vals=Append[vals,(EE[Pi k1/(nn+1)]+EE[ Pi k2/(nn+1)])//N];
]];
Print["N=2, k=",k," energ=",vals];
];
For[k=1,k<=nn,k++,vals={};
For[k1=1,k1<=nn,k1++,
For[k2=1,k2<=nn,k2++,
k3=k-k1-k2;
If[k3<=0,k3+=nn];
If[k3<=0,k3+=nn];
If[k3>k2 && k2>k1,vals=Append[vals,(EE[ Pi k1/(nn+1)]+EE[ Pi k2/(nn+1)]+EE[ Pi k3/(nn+1)])//N];
]]];
Print["N=3, k=",k," energ=",vals];
];
For[k=1,k<=nn,k++,vals={};
For[k1=1,k1<=nn,k1++,
For[k2=1,k2<=nn,k2++,
For[k3=1,k3<=nn,k3++,
k4=k-k1-k2-k3;
If[k4<=0,k4+=nn];
If[k4<=0,k4+=nn];
If[k4<=0,k4+=nn];
If[k4>k3&& k3>k2 && k2>k1,vals=Append[vals,(EE[Pi k1/(nn+1)]+EE[ Pi k2/(nn+1)]+EE[ Pi k3/(nn+1)]+EE[Pi k4/(nn+1)])//N];
]]]];
Print["N=4, k=",k," energ=",vals];
];
];

DoFerm[3]

Energy mapping with hard-core Bosons

2-site

NN=4;

Hb=-bp**bm-bm**bp+(U/2) ((bn.(bn-be))**be+be**(bn.(bn-be))); Hb//MatrixForm

Length[Hb]

Eigenvalues[Hb/.U->100.]//Sort//Chop

Take[Eigenvalues[Hb/.U->1000000000.]//Sort//Chop,5]

3-site

NN=4;

Hb=-(bp**bm+bm**bp)**be-be**(bp**bm+bm**bp)+(U/2) ((bn.(bn-be))**be**be+be**(bn.(bn-be))**be+be**be**(bn.(bn-be)));

Take[Eigenvalues[Hb/.U->10000000.]//Sort//Chop,9]

Hb=-(bp**bm+bm**bp)**be-be**(bp**bm+bm**bp)+(U/2) ((bn.(bn-be))**be**be+be**(bn.(bn-be))**be+be**be**(bn.(bn-be)))+0.01 (bn**be**be+be**bn**be+be**be**bn);

Take[Eigenvalues[Hb/.U->10000000.]//Sort//Chop,9]

DoFerm[3]

4-site

NN=4;

Hb=-(bp**bm+bm**bp)**be**be-be**(bp**bm+bm**bp)**be-be**be**(bp**bm+bm**bp)+(U/2) ((bn.(bn-be))**be**be**be+be**(bn.(bn-be))**be**be+be**be**(bn.(bn-be))**be+be**be**be**(bn.(bn-be)));

Take[Chop[Eigenvalues[Hb/.U->1000000.],0.00001]//Sort,17]

DoFerm[4]

Now, consider a periodic chain (length 3).

This is for fermions (up to 4 particles)

DoFermPer[nn_Integer]:=Block[{k,k1,k2,k3,k4,vals={}},
Print["N=0, k=",0," energ=",0//N];
For[k=0,k<nn,k++,Print["N=1, k=",k," energ=",(EE[ 2 Pi k/(nn)])//N]];
For[k=0,k<nn,k++,vals={};
For[k1=0,k1<nn,k1++,
k2=k-k1;
If[k2<0,k2+=nn];
If[k2>k1,vals=Append[vals,(EE[2 Pi k1/nn]+EE[ 2 Pi k2/nn])//N];
]];
Print["N=2, k=",k," energ=",vals];
];
For[k=0,k<nn,k++,vals={};
For[k1=0,k1<nn,k1++,
For[k2=0,k2<nn,k2++,
k3=k-k1-k2;
If[k3<0,k3+=nn];
If[k3<0,k3+=nn];
If[k3>k2 && k2>k1,vals=Append[vals,(EE[ 2Pi k1/nn]+EE[2 Pi k2/nn]+EE[ 2 Pi k3/nn])//N];
]]];
Print["N=3, k=",k," energ=",vals];
];
For[k=0,k<nn,k++,vals={};
For[k1=0,k1<nn,k1++,
For[k2=0,k2<nn,k2++,
For[k3=0,k3<nn,k3++,
k4=k-k1-k2-k3;
If[k4<0,k4+=nn];
If[k4<0,k4+=nn];
If[k4<0,k4+=nn];
If[k4>k3&& k3>k2 && k2>k1,vals=Append[vals,(EE[2 Pi k1/nn]+EE[ 2 Pi k2/nn]+EE[ 2 Pi k3/nn]+EE[2 Pi k4/nn])//N];
]]]];
Print["N=4, k=",k," energ=",vals];
];
];

For periodic chain length 3.

This is for bosons

Hb=-(bp**bm+bm**bp)**be-be**(bp**bm+bm**bp)-bp**be**bm-bm**be**bp+(U/2) ((bn.(bn-be))**be**be+be**(bn.(bn-be))**be+be**be**(bn.(bn-be)));

Take[Chop[Eigenvalues[Hb/.U->50000000.],0.0001]//Sort,9]

Same with fermions, k=2 π m/3, m=-1,0,1

DoFermPer[3]

Clearly,  the level degeneracy is different.

Try to fix the degeneracy, adding π/nn to the lattice momentum when N is even

DoBosPer[nn_Integer]:=Block[{k,k1,k2,k3,k4,vals={}},
Print["N=0, k=",0," energ=",0//N];
For[k=0,k<nn,k++,Print["N=1, k=",k," energ=",(EE[ 2 Pi k/(nn)])//N]];
For[k=0,k<nn,k++,vals={};
For[k1=0,k1<nn,k1++,
k2=k-k1;
If[k2<0,k2+=nn];
If[k2>k1,vals=Append[vals,(EE[Pi/nn+2 Pi k1/nn]+EE[Pi/nn+ 2 Pi k2/nn])//N];
]];
Print["N=2, k=",k," energ=",vals];
];
For[k=0,k<nn,k++,vals={};
For[k1=0,k1<nn,k1++,
For[k2=0,k2<nn,k2++,
k3=k-k1-k2;
If[k3<0,k3+=nn];
If[k3<0,k3+=nn];
If[k3>k2 && k2>k1,vals=Append[vals,(EE[ 2Pi k1/nn]+EE[2 Pi k2/nn]+EE[ 2 Pi k3/nn])//N];
]]];
Print["N=3, k=",k," energ=",vals];
];
For[k=0,k<nn,k++,vals={};
For[k1=0,k1<nn,k1++,
For[k2=0,k2<nn,k2++,
For[k3=0,k3<nn,k3++,
k4=k-k1-k2-k3;
If[k4<0,k4+=nn];
If[k4<0,k4+=nn];
If[k4<0,k4+=nn];
If[k4>k3&& k3>k2 && k2>k1,vals=Append[vals,(EE[Pi/nn+2 Pi k1/nn]+EE[Pi/nn+ 2 Pi k2/nn]+EE[Pi/nn+ 2 Pi k3/nn]+EE[Pi/nn+2 Pi k4/nn])//N];
]]]];
Print["N=4, k=",k," energ=",vals];
];
];

DoBosPer[3]

For chain n = 4 particles

Hb=-(bp**bm+bm**bp)**be**be-be**(bp**bm+bm**bp)**be-be**be**(bp**bm+bm**bp)-bp**be**be**bm-bm**be**be**bp+(U/2) ((bn.(bn-be))**be**be**be+be**(bn.(bn-be))**be**be+be**be**(bn.(bn-be))**be+be**be**be**(bn.(bn-be)));

Chop[Take[Sort[Eigenvalues[Hb/.U->5000000.]],17],10^-5]

Add a bit of chemical potential to distinguish energies by the particle number

Hb=-(bp**bm+bm**bp)**be**be-be**(bp**bm+bm**bp)**be-be**be**(bp**bm+bm**bp)-bp**be**be**bm-bm**be**be**bp+(U/2) ((bn.(bn-be))**be**be**be+be**(bn.(bn-be))**be**be+be**be**(bn.(bn-be))**be+be**be**be**(bn.(bn-be)))+0.01 (bn**be**be**be+be**bn**be**be+be**be**bn**be+be**be**be**bn);

Chop[Take[Sort[Eigenvalues[Hb/.U->5000000.]],17],10^-5]

DoBosPer[4]

DoFermPer[4]