Quantum spin-chains I. Extensive energy

The quantum antiferromagnetic Heisenberg chain is a chain of sites with dipole-dipole interactions between nearest neighbors

$H=J\sum_{i=1} \bm{S}_i\cdot\bm{S}_{i+1}=J\sum_{i}(S_{x,i}S_{x,i+1}+S_{y,i}S_{y,i+1}+S_{z,i}S_{z,i+1})$

We will work in a tensor basis in which the state of each site is

$v_1=|1/2,1/2\rangle=\left(\begin{array}{c}1\\ 0\end{array}\right)$

$v_2=|1/2,-1/2\rangle=\left(\begin{array}{c}0\\ 1\end{array}\right)$

Let $\bm{I}$ be the $2\times 2$ unit matrix,

$S_x={1\over 2}\left(\begin{array}{cc}0&1\\ 1&0\end{array}\right), \quad S_y={1\over 2}\left(\begin{array}{cc}0&-i\\ i&0\end{array}\right)$

$S_z={1\over 2}\left(\begin{array}{cc}1&0\\ 0&-1\end{array}\right)$

be the single particle spin operators in this basis. Then the spin operators for the $i^{th}$ site (acting only on the spin in the $i^{th}$ site) are

$S_{x,i}=\bm{I}\otimes \cdots\otimes \bm{I}\otimes S_x\otimes \bm{I}\otimes\cdots\otimes \bm{I}$

with the $S_x$ in the $i^{th}$ place in the tensor, $i=1$ copies of $\bm{I}$ to the left of it, $\ell-i-1$ copies of $\bm{I}$ to the right of it.

Julia provides us with some nice tools for manipulating tensors (these work the same way as their octave, python and lisp counterparts, so learn on, learn them all)

using LinearAlgebra, Kronecker
psi1=[1;0]
2-element Vector{Int64}:
 1
 0

psi2=[2 3]
1×2 Matrix{Int64}:
 2  3

psi1⊗psi1
4×1 Kronecker.KroneckerProduct{Int64, Matrix{Int64}, Matrix{Int64}}:
 1
 0
 0
 0

psi1⊗psi2
2×2 Kronecker.KroneckerProduct{Int64, Matrix{Int64}, Matrix{Int64}}:
 2  3
 0  0

psi2⊗psi2
1×4 Kronecker.KroneckerProduct{Int64, Matrix{Int64}, Matrix{Int64}}:
 4  6  6  9

psi3=[3;4]
2-element Vector{Int64}:
 3
 4

julia> psi1⊗psi3
4×1 Kronecker.KroneckerProduct{Int64, Matrix{Int64}, Matrix{Int64}}:
 3
 4
 0
 0

psi3⊗psi1
4×1 Kronecker.KroneckerProduct{Int64, Matrix{Int64}, Matrix{Int64}}:
 3
 0
 4
 0

OK, you get the drift. Let’s write the Hamiltonian in a nice form

$(S_{x,1}S_{x,2}+S_{y,1}S_{y,2}+S_{z,1}S_{z,2})=P={1\over 2}\left(\begin{array}{cccc} 1 & 0 &0 &0\\ 0 &-1 & 2 &0\\ 0 & 2 & -1 & 0\\ 0 & 0 & 0 & 1\end{array}\right)$

and deduce that

$(S_{x,i}S_{x,i+1}+S_{y,i}S_{y,i+1}+S_{z,i}S_{z,i+1})$

$=\bm{I}\otimes \cdots\otimes \bm{I}\otimes P\otimes \bm{I}\otimes\cdots\otimes \bm{I}=P_{i,i+1}$

with the $P$ in the $i^{th}$ and $i^{th}+1$ places in the tensor, $i=1$ copies of $\bm{I}$ to the left of it,
$\ell-i-2$ copies of $\bm{I}$ to the right of it.
so we have established that

$H=J\sum_{i} P_{i,i+1}$

Next we write a code that will create the Hamiltonian for a chain of length $N$

Hamiltonian=function(N::Int64)
sx=[0 1;1 0]
sy=[0 -im; im 0]
sz=[1 0; 0 -1]
u=[1 0; 0 1]
P=real.(sx⊗sx +sy⊗sy+sz⊗sz);
H=P;
Lend=P
for k=1:(N-2)
H=H⊗u+u⊗Lend;
Lend=u⊗Lend;
end
H/4.0
end

Time for some fun: we have previously talked about the Lanczos method for finding the largest magnitude eigenvalue of a matrix, so look under “Algorithms” and re-read that post. If you don’t want to, no problem, we can use the eigenvalue routine in LinearAlgebra. We write a code that finds the largest magnitude eigenvalue for chains of length $N=2,3,\cdots, N_{max}$ where $N_{max}$ is the largest array that does not kill the julia REPL on your machine. On the RockPro64 SBC that I am running this on (4 Gb RAM) this is $N_{max}=12$ , a $4096 \times 4096$ Hamiltonian. What we will do is show that the growth of this eigenvalue is linear in $N$

using LinearAlgebra, Kronecker, DataFrames, Gaston
Heisenberg=function(N::Int64,prec::Int64=8)
sx=[0 1;1 0]
sy=[0 -im; im 0]
sz=[1 0; 0 -1]
u=[1 0; 0 1]
P=real.(sx⊗sx +sy⊗sy+sz⊗sz);
H=P;
Lend=P
for k=1:(N-2)
H=H⊗u+u⊗Lend;
Lend=u⊗Lend;
end
H=-H/4.0;
maximum(round.(eigvals(H),digits=prec))
end


Dims=[2 3 4 5 6 7 8 9 10 11 12]
maxeigs=Heisenberg.(Dims,8)

We load the output and input values into a dataframe, run a linear regression to get the best fit parameters, and plot it. The energy of a many-body system should scale with the size of the system, this is the property of extensiveness. The maximal energy magnitude has this property:

dat=DataFrame(L=vec(float.(Dims)),E=vec(maxeigs))
11×2 DataFrame
 Row │ L        E       
     │ Float64  Float64 
─────┼──────────────────
   1 │     2.0  0.75
   2 │     3.0  1.0
   3 │     4.0  1.61603
   4 │     5.0  1.92789
   5 │     6.0  2.49358
   6 │     7.0  2.83624
   7 │     8.0  3.37493
   8 │     9.0  3.73632
   9 │    10.0  4.25804
  10 │    11.0  4.63209
  11 │    12.0  5.14209

LG=function(x,y)
       N=length(x)
       Den=N*sum(x .* x)-sum(x)^2
       a=N*sum(x .* y)-sum(x)*sum(y)
       b=sum(x .* x)*sum(y)-sum(x)*sum(x .* y)
       a/Den, b/Den
       end

a,b=LG(Dims[5:11],maxeigs[5:11])
(0.4435839410714287, -0.2103568610714412)

line=b .+ a .* Dims
dat.Line=vec(line)
plot(dat.L,dat.E)
plot!(dat.L,dat.Line, plotstyle="linespoints")
save(term="png",output="Heisenberg.png")

If you take time to check, you will see that the largest magnitude energy level is negative, so it is the ground state. If the number of spins in the chain is even and $J\ge 0$ (we set $J=1$ ) the model is antiferromagnetic and the ground state is unique. You have therefor found the $T=0$ Helmholtz free energy,since