Lecture 1 - Quantum mechanics of few-spin systems: translational symmetry

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;

Periodic Heisenberg spin chain (4 spins)

H=J(sx**sx**e**e+sx**e**e**sx+e**sx**sx**e+e**e**sx**sx+
sy**sy**e**e+sy**e**e**sy+e**sy**sy**e+e**e**sy**sy+
sz**e**e**sz+sz**sz**e**e+e**sz**sz**e+e**e**sz**sz); H//MatrixForm

Eigenvalues[H/.J->1]//N//Sort

Translation Operator

Convenient functions in Mathematica

IntegerDigits[100435,2,20]

RotateRight[%]

FromDigits[%,2]

For[j=0,j<2^4,j++,
Print[j,"   ",IntegerDigits[j,2,4],"  ",
RotateRight[IntegerDigits[j,2,4]],"   ",
FromDigits[ RotateRight[IntegerDigits[j,2,4]],2]]]

This function constructs translation operator for a chain of n spins

Transl[n_Integer]:=Block[{matr,NN,i,j,i1},
NN=(2^n);
matr=ConstantArray[0,{NN,NN}];
For[i=1,i≤NN,i++,j=i-1;
i1=1+FromDigits[RotateRight[IntegerDigits[j,2,n]],2];
matr[[i1,i]]=1];matr]

Try it out for n=4 spins

Transl[3]//MatrixForm

TT=Transl[4]; TT//MatrixForm

Translation repeated 4 times gives identity

(TT.TT.TT.TT)//MatrixForm

Eigenvalues have magnitude 1

Eigenvalues[TT]

Eigenvalues[Transl[5]]

Translation operator commutes with the Hamiltonian

(H.TT-TT.H)//MatrixForm

Take the ground state of the AF Hamiltonian

{vals,vecs}=Eigensystem[H/.J->1.];

vals

vecs

v0=vecs[[1]]

v1=vecs[[2]]

(TT.v0)-v0

Chop[(TT.v1)+v1]

Reduced Hamiltonian for k=0

One=IdentityMatrix[Length[TT]];

Projector to zero-momentum state

P=(One+TT+TT.TT+TT.TT.TT)/4;

Consider only the symmetric states with respect to translations.

H0=(P.H.P)/.J->1.;H0//MatrixForm

r=MatrixRank[P]

Take[Eigenvalues[H0]//N,r]

Quantum Ising model in a field

Same Hilbert space.   H = (∑_i) (Subscript[(S^z), i]) (Subscript[(S^z), i+1])+h (Subscript[(S^x), i])

This is how it looks for n = 3

HH3=4 (sz**sz**e+e**sz**sz+sz**e**sz)+2(sx**e**e+e**sx**e+e**e**sx);HH3//MatrixForm

Here is a program to construct the Hamiltonian directly

e

e**{{1}}

DoFullH[h_,n_]:=Block[{i,ee,lis,ans},
 ee={{1}};lis={ee};
For[i=1,i<=n,i++,
ee=e**ee;
lis=Append[lis,ee]];
ans=4 sz**lis[[n-1]]**sz;
For[i=1,i<=n-1,i++,
ans+=4 lis[[i]]**sz**sz**lis[[n-i]]];
For[i=1,i<=n,i++,
ans+=2 h*lis[[i]]**sx**lis[[n+1-i]]];
ans]

DoFullH[1.,3]//MatrixForm

Eigenvalues[DoFullH[1.,3]]//Sort//Chop

Eigenvalues[HH3]//N//Sort

Now, write a program to construct subspaces with different k

Find the optimal translation

DoSearch[i_Integer, n_Integer]:=Block[{m,ans,val},
ans={0,i};
For[m=1,m<n,m++,
val=FromDigits[ RotateRight[IntegerDigits[i,2,n],m],2];
 (* Print[ans,"  ",m,"  ",val];  *)
If[val<ans[[2]],ans={m,val}]];
ans];

DoSearch[16,5]

For[i=0,i<2^5,i++,Print[i,"  ",DoSearch[i,5]," ",IntegerDigits[i,2,5]]]

This flips a spin at position m of n

FlipSpin[i_,m_Integer,n_Integer]:=Block[{lis},
lis=IntegerDigits[i,2,n];
(* Print["was:   ",lis];  *)
lis[[m]]=1-lis[[m]];
 (* Print["became:",lis]; *) 
FromDigits[lis,2]]

FlipSpin[11,2,7]

IntegerDigits[11,2,7]

IntegerDigits[43,2,7]

This computes the PE for a chain of length n

DoPE[i_,n_Integer]:=Block[{ans=0,lis,m},
lis=IntegerDigits[i,2,n];
(*  Print[lis]; *)
If[lis[[1]]==lis[[n]],ans++,ans--]
For[m=1,m<n,m++,
If[lis[[m]]==lis[[m+1]],ans++,ans--]];
ans]

DoPE[0,5]

This function returns the degeneracy of the corresponding state.

Degen[i_,n_Integer]:=Block[{ans,lis,m},
lis=IntegerDigits[i,2,n];
ans=1;
For[m=1,m<n,m++,
If[lis==RotateRight[lis,m],ans++]];
n/ans]

This function checks whether the symmetrization returns zero

CheckProj[i_,k_Integer,n_Integer]:=Block[{ans,m},
(* For a canonical state, check that it does not correspond to zero *)
lis=IntegerDigits[i,2,n];
ans=1;
For[m=1,m<n,m++,
{shift,val}=DoSearch[FromDigits[RotateRight[lis,m],2],n];
(*Print[m,"  ",shift];*)
If[val≠i,Print["error!"]];
ans+=Exp[I 2 Pi  (shift+m) k/n];
];
ans=FullSimplify[ans]/n
]

CheckProj[0,4,5]

This does the actual Hamiltonian in the subspace with the momentum k

DoHamQI[h_,k_Integer,n_Integer]:=Block[{i,m,energ,j,j1,in,mn,jn,jn1,val,ans,pos,dim,transl,fromlab={}},
(* First, calculate the dimension of the Hilbert space needed & construct the forward translation table *)
For[i=0;dim=0,i<(2^n),i++,If[DoSearch[i,n][[1]]==0&&CheckProj[i,k,n]≠0,dim++;fromlab=Append[fromlab,i]]];
Print["fromlab=",fromlab];
(* now construct the backward translation table *)
transl=Array[0&,(2^n)];For[j=1,j≤dim,j++,transl[[1+fromlab[[j]]]]=j];
Print["dim=",dim,",  transl=",transl];
(* initialize an empty array for the Hamiltonian *)
ans=ConstantArray[0,{dim,dim}];
For[i=1,i≤dim,i++,j=fromlab[[i]]; (* j is the binary encoding of the corresponding vector *)
Print["i=",i," state=",IntegerDigits[j,2,n]];
ans[[i,i]]+=DoPE[j,n]; (* add the potential energy *)
For[pos=1,pos≤n,pos++, (* go over positions for (σ^x) operators *)
in=FlipSpin[j,pos,n]; (* this is the result of flip-spin operation *)
{mn,jn}=DoSearch[in,n]; (* displacement and the canonic representation *)
Print["pos=",pos," stateX=",IntegerDigits[jn,2,n]," mn=",mn];
jn1=transl[[jn+1]]; (* the corresponding column in the Hilbert state *)
If[jn1≠0,Print["jn1=",jn1," i=",i];ans[[jn1,i]]+=((h (E^((I 2 π k mn)/n)) (Sqrt[Degen[j,n]]))/(Sqrt[Degen[jn,n]]));]]];
ans]

ham=DoHamQI[1, 0,3];ham//MatrixForm

This is the same function without the printouts

DoHamQI[h_,k_Integer,n_Integer]:=Block[{i,m,energ,j,j1,in,mn,jn,jn1,val,ans,pos,dim,transl,fromlab={}},For[i=0;dim=0,i<(2^n),i++,
If[DoSearch[i,n][[1]]==0&&CheckProj[i,k,n]≠0,
dim++;
fromlab=Append[fromlab,i]]];
transl=Array[0&,(2^n)];
For[j=1,j≤dim,j++,transl[[1+fromlab[[j]]]]=j];
ans=ConstantArray[0,{dim,dim}];
For[i=1,i≤dim,i++,j=fromlab[[i]];
ans[[i,i]]+=DoPE[j,n];
For[pos=1,pos≤n,pos++,in=FlipSpin[j,pos,n];{mn,jn}=DoSearch[in,n];jn1=transl[[jn+1]];If[jn1≠0,ans[[jn1,i]]+=((h (E^((I 2 π k mn)/n)) (Sqrt[Degen[j,n]]))/(Sqrt[Degen[jn,n]]));]]];ans]

For[k=0,k<5,k++,
Print["k=",k,"  ",Eigenvalues[DoHamQI[1.,k,5]]//Chop//Sort]]

What remains is to compare with the spectrum constructed explicitly.  Do it for n=5:

HH5=DoFullH[1,5];HH5//MatrixForm

Eigenvalues[HH5]//N//Sort

Also, for n = 4

HH4=4 (sz**sz**e**e+e**sz**sz**e+e**e**sz**sz+sz**e**e**sz)+
2(sx**e**e**e+e**sx**e**e+e**e**sx**e+e**e**e**sx); HH4//MatrixForm

Eigenvalues[HH4//N]//Chop//Sort

For[k=0,k<4,k++,
Print["k=",k,"  ",DoHamQI[1,k,4]//MatrixForm,"  ",Eigenvalues[DoHamQI[1.,k,4]]//Sort//Chop]]

Also, for n = 3

HH=4 (sz**sz**e+e**sz**sz+sz**e**sz)+2(sx**e**e+e**sx**e+e**e**sx);HH//MatrixForm

Eigenvalues[ HH]//N//Sort

For[k=0,k<3,k++,
Print["k=",k,"  ",DoHamQI[1,k,3]//MatrixForm,"  ",Eigenvalues[DoHamQI[1.,k,3]]//Sort//Chop]]

Also, for n = 2 (here the potential energy has to be doubled to account for two ways to connect)

For[k=0,k<2,k++,
Print["k=",k,"  ",DoHamQI[1,k,2]//MatrixForm,"  ",Eigenvalues[DoHamQI[1.,k,2]]//Sort]]

HH=8(sz**sz)+2(sx**e+e**sx);HH//MatrixForm

Eigenvalues[HH]//N//Sort

Also, for n = 10 (check how long it takes the usual way)

Timing[HH10=DoFullH[1,10];Eigenvalues[HH10//N]//Sort]

And also by diagonalizing each sector separately

Timing[For[k=0,k<10,k++,
Print["k=",k,"  ",Eigenvalues[DoHamQI[1.,k,10]//N]//Chop//Sort]]]

This is why we went to all this trouble - compare the calculation time for 12 spins

Timing[For[k=0,k<12,k++,
Print["k=",k,"  ",Eigenvalues[DoHamQI[1.,k,12]//N]//Chop//Sort]]]

Timing[HH12=DoFullH[1,12];Eigenvalues[HH12//N]//Sort]

42 seconds vs. 237 for full - matrix calculation