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]];

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]);

MP[1] := [[1, 1, 0, 0], [1, 0, 1, 0], [1, 0, 0, 1],...

MP[2] := [[1, 1, 0, 0], [1, 0, 1, 0], [1, 0, 0, 1],...

n := 6

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)];

Basis := [vec([1, 1, 0, 0],[1, 1, 0, 0]), vec([1, 1...
Basis := [vec([1, 1, 0, 0],[1, 1, 0, 0]), vec([1, 1...
Basis := [vec([1, 1, 0, 0],[1, 1, 0, 0]), vec([1, 1...
Basis := [vec([1, 1, 0, 0],[1, 1, 0, 0]), vec([1, 1...
Basis := [vec([1, 1, 0, 0],[1, 1, 0, 0]), vec([1, 1...
Basis := [vec([1, 1, 0, 0],[1, 1, 0, 0]), vec([1, 1...
Basis := [vec([1, 1, 0, 0],[1, 1, 0, 0]), vec([1, 1...
Basis := [vec([1, 1, 0, 0],[1, 1, 0, 0]), vec([1, 1...

> dimH:=n^2;

dimH := 36

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);

H := array(1 .. 36,1 .. 36,[])

> 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);

matrix([[2*U, -t, 0, 0, t, 0, -t, 0, 0, 0, 0, 0, 0,...

>

>

>

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);

matrix([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...

>

>

> 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];

Um := [0, .2, .6, 1, 2, 3, 4, 5, 6, 8, 10, 13, 16, ...

> 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:

`Um=`, 0, `  Etot=`, -1.000000002

`Um=`, .2, `  Etot=`, -.9634925615

`Um=`, .6, `  Etot=`, -.8960396740

`Um=`, 1, `  Etot=`, -.8352119050

`Um=`, 2, `  Etot=`, -.7071067868

`Um=`, 3, `  Etot=`, -.6061072348

`Um=`, 4, `  Etot=`, -.5256871202

`Um=`, 5, `  Etot=`, -.4610721975

`Um=`, 6, `  Etot=`, -.4086507600

`Um=`, 8, `  Etot=`, -.3300587250

`Um=`, 10, `  Etot=`, -.2749694425

`Um=`, 13, `  Etot=`, -.2185721200

`Um=`, 16, `  Etot=`, -.1807160750

`Um=`, 19, `  Etot=`, -.1537559648

`Um=`, 21, `  Etot=`, -.1397620191

`Um=`, 26, `  Etot=`, -.1137282884

`Um=`, 32, `  Etot=`, -.9285262000e-1

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;

Hmid1 := -1.000000002*t

Hmid2 := .3553754645e-1-.9990301082*t

Hmid3 := .9598825382e-1-.9920279280*t

Hmid4 := .1444967726-.9797086782*t

Hmid5 := .2265409204-.9336477038*t

Hmid6 := .2693875392-.8754947668*t

Hmid7 := .2873253722-.8130124922*t

Hmid8 := .2899793245-.7510515182*t

Hmid9 := .2838105058-.6924612685*t

Hmid10 := .2599837825-.5900425205*t

Hmid11 := .2328693620-.5078388052*t

Hmid12 := .1967130871-.4152851990*t

Hmid13 := .1681304288-.3488465098*t

Hmid14 := .1459171634-.2996731352*t

Hmid15 := .1338442420-.2736062488*t

Hmid16 := .1105130644-.2242413703*t

Hmid17 := .9109318095e-1-.1839457866*t

> 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)));

[Maple Plot]

[Maple Plot]

[Maple Plot]

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;

spin1 := -.5675967218e-1

spin2 := -.2915082015

spin3 := -.3103831545

spin4 := -.3272913218

spin5 := -.3625589445

spin6 := -.3900293412

spin7 := -.4115160680

spin8 := -.4282873800

spin9 := -.4413607025

spin10 := -.4595885645

spin11 := -.4709940952

spin12 := -.4811267228

spin13 := -.4868842922

spin14 := -.4904104325

spin15 := -.4920400870

spin16 := -.4946899528

spin17 := -.4964430442

> 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)));

[Maple Plot]

>