Anisotropic three spin problem

©  Prof. Dr. Jürgen Schnack, Universität Bielefeld, 24. 09. 2013

Initialization

NS=3;
spin={2,2,2}; 
MMax =spin[[1]]+spin[[2]]+spin[[3]];
DimH=(2*spin[[1]]+1)*(2*spin[[2]]+1)*(2*spin[[3]]+1)

Konstanten

J=-1*(({{0, -1.0, -1.0}, {-1.0, 0, -1.0}, {-1.0, -1.0, 0}}))  (* Kelvin (Subscript[k, B]) *);
d=-0.2; (* easy axis *)
d1=d(* Kelvin (Subscript[k, B]) *);
d2=d(* Kelvin (Subscript[k, B]) *);
d3=d(* Kelvin (Subscript[k, B]) *);
gmuBDurchkB=1.3434 (* K/T *);

Create basis

Basis=Table[{0,0,0},{n,1,DimH}];

For[i=2,i<DimH+1,i++,
For[j=1,j<NS+1,j++,
Basis[[i,j]]=Basis[[i-1,j]]];
Basis[[i,1]]=Basis[[i,1]]+1;
j=1;
While[j<NS,
If[Basis[[i,j]]==(2*spin[[j]]+1),Basis[[i,j]]=0;Basis[[i,j+1]]=Basis[[i,j+1]]+1;j++,j=NS]]
];

MatrixForm[Basis];

Output table

EValues=Table[{0,0},{n,1,1}];

Work in subspaces H(M)

For[M=MMax,M>=0,M--,
MBasis=Table[{0,0,0,0,0,0,0},{n,1,1}];
For[i=1,i<DimH+1,i++,
ThisA=(∑_(j=1)^NS)Basis[[i,j]];
If[ThisA==MMax-M,AppendTo[MBasis,Basis[[i]]]]];
MBasis=Delete[MBasis,1];
DimHM=Length[MBasis];
(* Print[TableForm[MBasis]] *)
HMatrix=Table[0.0,{i,1,DimHM},{j,1,DimHM}];
For[i=1,i<DimHM+1,i++,
For[k=1,k<NS+1,k++,
For[k2=k+1,k2<NS+1,k2++,
HMatrix[[i,i]]=HMatrix[[i,i]]-2*J[[k,k2]]*
          (spin[[k]]-MBasis[[i,k]])*(spin[[k2]]-MBasis[[i,k2]])]]];
For[i=1,i<DimHM+1,i++,
For[j=1,j<DimHM+1,j++,
For[k=1,k<NS+1,k++,
For[k2=k+1,k2<NS+1,k2++,
HilfsElement=Table[MBasis[[j,n]],{n,1,NS}];
HilfsElement[[k]]=HilfsElement[[k]]-1;
HilfsElement[[k2]]=HilfsElement[[k2]]+1;
HMatrix[[i,j]]=HMatrix[[i,j]]-2*J[[k,k2]]*0.5*
Product[KroneckerDelta[MBasis[[i,m]]-HilfsElement[[m]]],{m,1,NS}]*
√((MBasis[[j,k]])*(2*spin[[k]]-MBasis[[j,k]]+1))*
√((2*spin[[k2]]-MBasis[[j,k2]])*(MBasis[[j,k2]]+1))
]]]];
For[i=1,i<DimHM+1,i++,
For[j=1,j<DimHM+1,j++,
For[k=1,k<NS,k++,
For[k2=k+1,k2<NS+1,k2++,
HilfsElement=Table[MBasis[[j,n]],{n,1,NS}];
HilfsElement[[k]]=HilfsElement[[k]]+1;
HilfsElement[[k2]]=HilfsElement[[k2]]-1;
HMatrix[[i,j]]=HMatrix[[i,j]]-2*J[[k,k2]]*0.5*
Product[KroneckerDelta[MBasis[[i,m]]-HilfsElement[[m]]],{m,1,NS}]*
√((MBasis[[j,k2]])*(2*spin[[k2]]-MBasis[[j,k2]]+1))*
√((2*spin[[k]]-MBasis[[j,k]])*(MBasis[[j,k]]+1))
]]]];
EV0=Sort[Chop[Eigenvalues[HMatrix]]];
(* Print[TableForm[EV0]] *)
For[j=1,j<DimHM+1,j++,
		HMatrix[[j,j]]=HMatrix[[j,j]]+
d1*(spin[[1]]-MBasis[[j,1]])^2+
d2*(spin[[2]]-MBasis[[j,2]])^2+
d3*(spin[[3]]-MBasis[[j,3]])^2
];
EV=Sort[Chop[Eigenvalues[HMatrix]]];
(* Print[TableForm[EV]] *)
For[i=1,i<Length[EV]+1,i++,
AppendTo[EValues,{M,EV[[i]]}]];
If[M>0,For[i=1,i<Length[EV]+1,i++,
AppendTo[EValues,{-M,EV[[i]]}]]]
];

Eigenvalues and Graphics

EValues=Delete[EValues,1];
Print[TableForm[EValues]]

ListPlot[EValues,BaseStyle->{FontFamily->"Helvetica",FontSize->18},PlotStyle->{{PointSize[Large],Red}},
AxesLabel->{"M","E/((Subscript[k, B])K)"},PlotLabel->"Energieeigenwerte"]

Tasks

- af coupling, easy axis
- af coupling, hard axis
- f coupling, easy axis
- f coupling, hard axis
- strong/weak exchange limit