Hubbard Model
Dmitri M Belov
belov@physics.rutgers.edu
Load packages.
> restart;
> with(combinat):
Warning, the protected name Chi has been redefined and unprotected
Prepare 4-site vectors with 2 spin up and 2 spin down
This menas (up,up,0,0)x(down,down,0,0)
> v:=[[1,1,0,0],[1,1,0,0]];
Let us make all possible permutations of this state:
>
MP[1]:=permute(v[1]);
MP[2]:=permute(v[2]);
n:=nops(MP[1]);
Here n is a total number of vectors with 2 spin up (down).
Now define the basis of the Hilbert space as
> Basis:=[seq(seq(vec(MP[1][j],MP[2][i]),i=1..n),j=1..n)];
> dimH:=n^2;
Procedures for work with Hilbert space
Inner product (assume that the Hilbert space is real)
>
int_product:=proc(x1,x2)
local i,rst;
equal:=proc(A,B)
if (A=B) then
1;
else
0;
end if;
end proc;
i_p:=proc(x,y)
local a,X,Y,j;
a:=1;
if op(0,x)=`*` then
a:=a*convert([seq(op(j,x),j=1..(nops(x)-1))],`*`);
X:=op(nops(x),x);
else
X:=x;
end if;
if op(0,y)=`*` then
a:=a*convert([seq(op(j,y),j=1..(nops(y)-1))],`*`);
Y:=op(nops(y),y);
else
Y:=y;
end if;
expand(a*convert([seq(equal(op(1,X)[i],op(1,Y)[i]),i=1..4),
seq(equal(op(2,X)[i],op(2,Y)[i]),i=1..4)],`*`));
end proc:
if (x1=0) or (x2=0) then
0;
else
if op(0,x1)=`+` then
if op(0,x2)=`+` then
rst:=0;
for i from 1 to nops(x2) do
rst:=rst+convert(map(i_p,[op(x1)],op(i,x2)),`+`);
od;
rst;
else #=`*` or nothing
convert(map(i_p,[op(x1)],x2),`+`);
end if;
else
if op(0,x2)=`+` then
convert(map(i_p,[op(x2)],x1),`+`);
else #=`*` or nothing
i_p(x1,x2);
end if;
end if;
end if;
end proc:
Warning, `equal` is implicitly declared local to procedure `int_product`
Warning, `i_p` is implicitly declared local to procedure `int_product`
Test of Inner Product
Anihilation operator for particle with spin "sigma" on site "i" acting on vector "vect"
sigma=1 --- spin up
sigma=3 --- spin down
>
ind:=proc(a)
if a=5 then
1;
else
a;
end if;
end proc:
>
Cm:=proc(i,sigma,VEct)
#sigma=1 or 2
local sgn,ss,j,m,VVV;
cmm:=proc(vect,i,sigma)
if vect=0 then
0;
else
if op(0,vect)=`*` then
m:=op(1,vect);
VVV:=op(2,vect);
else
m:=1;
VVV:=vect;
end if;
if sigma=1 then
ss:=op(1,VVV);
sgn:=(-1)^(sum(ss[j],j=1..(i-1)));
if ss[i]=1 then
ss[i]:=0;
m*sgn*vec(ss,op(2,VVV));
else
0;
end if;
else #sigma=2
ss:=op(2,VVV);
sgn:=(-1)^(sum(ss[j],j=1..(i-1))+sum(op(1,VVV)[j],j=1..4));
if ss[i]=1 then
ss[i]:=0;
m*sgn*vec(op(1,VVV),ss);
else
0;
end if;
end if;
end if;
end proc:
if op(0,VEct)=`+` then
convert(map(cmm,[op(VEct)],ind(i),sigma),`+`);
else
cmm(VEct,ind(i),sigma);
end if;
end proc:
Warning, `cmm` is implicitly declared local to procedure `Cm`
Test of Cm
Anihilation operator for particle with spin "sigma" on site "i" acting on vector "vect"
sigma=1 --- spin up
sigma=3 --- spin down
>
Cp:=proc(i,sigma,VEct)
#sigma=1 or 2
local sgn,ss,j,m,VVV;
cpp:=proc(vect,i,sigma)
if vect=0 then
0;
else
if op(0,vect)=`*` then
m:=op(1,vect);
VVV:=op(2,vect);
else
m:=1;
VVV:=vect;
end if;
if sigma=1 then
ss:=op(1,VVV);
sgn:=(-1)^(sum(ss[j],j=1..(i-1)));
if ss[i]=0 then
ss[i]:=1;
m*sgn*vec(ss,op(2,VVV));
else
0;
end if;
else #sigma=2
ss:=op(2,VVV);
sgn:=(-1)^(sum(ss[j],j=1..(i-1))+sum(op(1,VVV)[j],j=1..4));
if ss[i]=0 then
ss[i]:=1;
m*sgn*vec(op(1,VVV),ss);
else
0;
end if;
end if;
end if;
end proc;
if op(0,VEct)=`+` then
convert(map(cpp,[op(VEct)],ind(i),sigma),`+`);
else
cpp(VEct,ind(i),sigma);
end if;
end proc:
Warning, `cpp` is implicitly declared local to procedure `Cp`
Test of Cp
Test of Cp+Cm
Operator of number of particles
>
Np:=proc(i,sigma,VEct)
#sigma=1 or 2
local ss,j,m,VVV;
npp:=proc(vect,i,sigma)
if vect=0 then
0;
else
if op(0,vect)=`*` then
m:=op(1,vect);
VVV:=op(2,vect);
else
m:=1;
VVV:=vect;
end if;
if sigma=1 then
ss:=op(1,VVV);
m*ss[i]*vec(ss,op(2,VVV));
else #sigma=2
ss:=op(2,VVV);
m*ss[i]*vec(op(1,VVV),ss);
end if;
end if;
end proc:
if op(0,VEct)=`+` then
convert(map(npp,[op(VEct)],i,sigma),`+`);
else
npp(VEct,i,sigma);
end if;
end proc:
Warning, `npp` is implicitly declared local to procedure `Np`
Test of Np
>
Spin Operators
We know how to write spin operators in terms of creation/anihilation operators
"i" here means site, and "vect" means vector in the Hilbert space
>
Sz:=proc(i,vect)
1/2*(Cp(i,1,Cm(i,1,vect))-Cp(i,2,Cm(i,2,vect)));
end proc:
Sx:=proc(i,vect)
1/2*(Cp(i,1,Cm(i,2,vect))+Cp(i,2,Cm(i,1,vect)));
end proc:
Sy:=proc(i,vect)
1/2/I*(Cp(i,1,Cm(i,2,vect))-Cp(i,2,Cm(i,1,vect)));
end proc:
Test of Spin operators
Hamiltonian
Kinetic term
>
H0:=proc(vect)
local sigma,k,res;
res:=-t*convert([seq(
convert([seq(
Cp(k,sigma,Cm(k+1,sigma,vect))+Cp(k+1,sigma,Cm(k,sigma,vect))
,k=1..4)],`+`)
,sigma=1..2)],`+`);
expand(res);
end proc:
Interaction term
>
Hint:=proc(vect)
local sigma,sigmap,i,res;
res:=U*sum(Np(i,1,Np(i,2,vect)),i=1..4);
expand(res);
end proc:
Computations of the Hamiltonian as Matrix
> H:=array(1..dimH,1..dimH);
>
for i from 1 to dimH do
for j from 1 to dimH do
Hv:=H0(Basis[j])+Hint(Basis[j]);
H[i,j]:=int_product(Basis[i],Hv);
od;
od;
print(H);
>
>
>
Spin matrix
>
SiSj:=proc(i,j,vect)
Sx(i,Sx(j,vect))+Sy(i,Sy(j,vect))+Sz(i,Sz(j,vect));
end proc:
>
ind:=proc(a)
if a=5 then
1;
else
a;
end if;
end proc:
>
spin_mat:=array(1..dimH,1..dimH):
for i from 1 to dimH do
for j from 1 to dimH do
spin_mat[i,j]:=0;
for k from 1 to 4 do
spin_mat[i,j]:=spin_mat[i,j]+int_product(Basis[i],SiSj(k,ind(k+1),Basis[j]));
od:
od:
od:
print(spin_mat);
>
>
> save dimH,H,spin_mat,"hubb.m";
>
Calc. of the Ground State
> restart;
> read "hubb.m";
> with(linalg):
Warning, the protected names norm and trace have been redefined and unprotected
Let us compute Hamiltonian for t=1 and U=2,4,8
> Um:=[0,0.2,0.6,1,2,3,4,5,6,8,10,13,16,19,21,26,32];
>
for i from 1 to nops(Um) do
H||i:=map2(subs,{t=1,U=Um[i]},H):
od:
>
#[eigenvalues(H0)];
#[eigenvalues(H1)];
#[eigenvalues(H2)];
#[eigenvalues(H3)];
Let us compute eigenvalues
>
for i from 1 to nops(Um) do
E||i:=evalf(Eigenvals(H||i,vecs||i));
od:
Total energy
>
for i from 1 to nops(Um) do
Etot||i:=E||i[1]/4;
print(`Um=`,Um[i],` Etot=`,Etot||i);
od:
Ground states
>
for j from 1 to nops(Um) do
VeC||j:=array(1..dimH);
for i from 1 to dimH do
VeC||j[i]:= vecs||j[i,1];
od:
od:
Mean Kinetic and Interaction Energy
>
for i from 1 to nops(Um) do
Hmid||i:=subs({U=Um[i]},evalm(VeC||i&*H&*VeC||i))/4;
od;
>
PLOT(CURVES([seq([Um[i],diff(Hmid||i,t)],i=1..nops(Um))],THICKNESS(3),LINESTYLE(4),COLOR(RGB,1,0,0)),TEXT([10,diff(Hmid5,t)],'`Kinetic Energy`',COLOR(RGB,1,0,0)));
PLOT(CURVES([seq([Um[i],subs(t=0,Hmid||i)],i=1..nops(Um))],THICKNESS(3),LINESTYLE(3),COLOR(RGB,0,0,1)),TEXT([10,0.1],'`Interaction Energy`',COLOR(RGB,0,0,1)));
PLOT(CURVES([seq([Um[i],Etot||i],i=1..nops(Um))],THICKNESS(1),LINESTYLE(1)),TEXT([10,-.6],'`Total energy`',COLOR(RGB,0,0,0)));
Spin correlations
>
for i from 1 to nops(Um) do
spin||i:=subs({U=Um[i]},evalm(VeC||i&*spin_mat&*VeC||i))/4;
od;
>
PLOT(CURVES([seq([Um[i],spin||i],i=1..nops(Um))],THICKNESS(3),LINESTYLE(4),COLOR(RGB,1,0,0)),
TEXT([10,-0.3],'`Spin`',COLOR(RGB,1,0,0)));
>