# 1D Spin Chain Values and Vectors

Check out Part 1 for the theoretical background. Today is all programming.

Remember to keep a reasonable number of spins, $n$. The number of states goes as $2^n$, and the size of the Hamiltonian will go as $2^n \times 2^n$. A 10-spin chain will have over a million entries, before taking into account any memory reducing tricks.

Here’s the details of the largest spin-chains that fit on my machine which has 16GB of RAM and a 3.2 GHz Intel i5 processor:

n | m_z | t (min) | Memory (Gb) |
---|---|---|---|

16 | 8 | 1.5 | 2 |

18 | 9 | ? | >16 |

18 | 8 | 32.1 | 15.5 |

I have included a file in this directory, ED.jl, that is just the necessary executable parts of this Jupyter notebook. For large $n$, I recommend running ED.jl.

```
using LinearAlgebra
```

Here, we input one parameter `n`

, the number of spins in our chain.

The program automatically calculates the parameter `nstates`

.

```
n=4
nstates=2^n
```

Now, let’s write out all of our possible states in the $S^z$ basis.

```
psi=convert.(Int8, collect(0:(nstates-1)) )
for p in psi
println(bitstring(p)[end-n:end],' ',p)
end
```

```
00000 0
00001 1
00010 2
00011 3
00100 4
00101 5
00110 6
00111 7
01000 8
01001 9
01010 10
01011 11
01100 12
01101 13
01110 14
01111 15
```

As in Part 1, we will be using the powers of 2 to compute magnetization, and masks to flip spins. To not have to calculate them each time, we just store them in memory.

```
powers2=collect(0:(n-1));
powers2=convert.(Int8, 2.0 .^powers2 );
mask=[0;powers2]+[powers2;0];
mask=convert.(Int8, [mask[2:end-1];[1+2^(n-1)]])
for m in mask
println(bitstring(m)[end-n:end])
end
```

```
00011
00110
01100
01001
```

In Part 1, I used the number of up-spins as a proxy for magnetization. Here, we need the *actual* magnetization, not a proxy. An up-spin is $+1/2$ and a down-spin is $-1/2$. We modify our magnetization by

\begin{equation} m=\frac{1}{2}\left(n_{\uparrow}-n_{\downarrow} \right) = \frac{1}{2} \left( n_{\uparrow}- \left( n - n_{\uparrow} \right) \right), \end{equation}

\begin{equation} m = n_{\uparrow} - \frac{n}{2}. \end{equation}

```
m=zeros(Int8,length(psi))
for i in 1:nstates
m[i]=sum((psi[i].&powers2)./(powers2))
end
m=m.-n/2
```

```
16-element Array{Float64,1}:
-2.0
-1.0
-1.0
0.0
-1.0
0.0
0.0
1.0
-1.0
0.0
0.0
1.0
0.0
1.0
1.0
2.0
```

## Grouped by Magnetization

Now that we have the magnetizations corresponding to each state, we perform some trickery that allows us to reduce the difficulty of our problem dramatically.

Magnetization is a conserved quantity. By Noether’s theorem, we know that the Hamiltonian is not going to mix states of different magnetizations. We only deal with one magnetization at a time, which is a much smaller problem.

```
# The possible values for magnetization
ma=collect(0:1:n).-n/2
```

```
5-element Array{Float64,1}:
-2.0
-1.0
0.0
1.0
2.0
```

Now let’s just pick out a single magnetization quantum number $mz$ and only work with that matrix for the rest of the post.

```
# The magnetic quantum number
mz=3
# An array of states with the correct magnetization
psi_mz=psi[m.==ma[mz]]
[psi_mz bitstring.(psi_mz) m[psi_mz.+1]]
```

```
6×3 Array{Any,2}:
3 "00000011" 0.0
5 "00000101" 0.0
6 "00000110" 0.0
9 "00001001" 0.0
10 "00001010" 0.0
12 "00001100" 0.0
```

And now creating the matrix.

Stuff goes here

```
dim=length(psi_mz)
M=ma[mz]*(ma[mz]+1)*Matrix{Float64}(I,dim,dim)
#M=zeros(Float64,dim,dim); use this for XY model
```

Though we have significantly reduced the size of the matrix by restricting to one magnetization, we no longer have our states in `1,2,3,4...`

order. Their position in an array no longer corresponds to their value. Therefore, we need a function to determine their index once we know their value.

We can find the index of the flipped state multiple different ways, but the simplest is by the Midpoint method. We split the interval in half, and see if the value we are looking for is higher or lower than the middle point. Then we get a new interval.

```
function findstate(state,set)
#Lower bound of interval
imin=1
#Upper bound of interval
imax=length(set)
# checking if the lower bound is what we are looking for
if set[imin] == state
return imin
end
# checking if the upper bound is what we are looking for
if set[imax] == state
return imax
end
# Initializing variables
# looking to see if we've found it yet
found=false
# how many times we've gone around the while loop
count=0
while found==false && count < length(set)
count+=1
tester=floor(Int,imin+(imax-imin)/2)
if state < set[tester]
imax=tester-1
elseif state > set[tester]
imin=tester+1
else
found=true
return tester
end
end
if found == false
println("findstate never found a match")
println("Are you sure the state is in that Array?")
end
return 0
end
```

Now time to generate the matrix.

For each state and for each pair of adjacent spins within that state, we apply the operator that flips adjacent spins, `$mask`

. Sometimes the adjacent spins will take on the same value, `00`

or `11`

. In this circumstance, the off-diagonal part of the Hamiltonian would not act on those spins. The state generated by the operator would have a different magnetization, and we can neglect that pair.

If the new state produced by this process has the same magnetization, we know that the flip exists in the Hamiltonian, and add the entry to the matrix.

In this algorithm, we do end up going over each pair twice, but I have yet to figure out how to take advantage of the degeneracy to cut the calculation in half. Let me know if you have a better way to write this.

```
mp=sum(psi_mz[1].&powers2./powers2)
for ii in 1:length(psi_mz)
p=psi_mz[ii]
for jj in 1:n
flipped=p⊻mask[jj]
if sum((flipped.&powers2)./powers2) == mp
tester=findstate(flipped,psi_mz)
M[ii,tester]=.5
M[tester,ii]=.5
println(bitstring(p)[end-n:end],'\t',bitstring(flipped)[end-n:end])
end
end
end
```

```
00011 00101
00011 01010
00101 00110
00101 00011
00101 01001
00101 01100
00110 00101
00110 01010
01001 01010
01001 00101
01010 01001
01010 01100
01010 00110
01010 00011
01100 01010
01100 00101
```

```
M
```

```
6×6 Array{Float64,2}:
0.0 0.5 0.0 0.0 0.5 0.0
0.5 0.0 0.5 0.5 0.0 0.5
0.0 0.5 0.0 0.0 0.5 0.0
0.0 0.5 0.0 0.0 0.5 0.0
0.5 0.0 0.5 0.5 0.0 0.5
0.0 0.5 0.0 0.0 0.5 0.0
```

```
F=eigen(M)
display(F.values)
display(F.vectors)
```

```
6-element Array{Float64,1}:
-1.4142135623730943
0.0
0.0
0.0
1.9984014443252818e-15
1.414213562373095
6×6 Array{Float64,2}:
0.353553 0.0 -0.211325 0.788675 -0.288675 -0.353553
-0.5 0.707107 0.0 0.0 -5.1279e-16 -0.5
0.353553 0.0 0.788675 -0.211325 -0.288675 -0.353553
0.353553 0.0 -0.57735 -0.57735 -0.288675 -0.353553
-0.5 -0.707107 0.0 0.0 -5.1279e-16 -0.5
0.353553 0.0 0.0 0.0 0.866025 -0.353553
```

Now we have eigenvalues and eigenvectors! You just solved the Heisenburg Spin Chain!

In my next post, I will analyze what this tells us about the system and what we can do with the information.