Today's just background :( I found my post getting a bit bloated, so I decided to pull most of the talking into an introduction post for the series. Keep reading to learn how to compute, diagonalize, and analyze the matrix corresponding to a many-body, strongly interacting quantum mechanical system.
If you are not a physics graduate student, still give this series a read and try not to get too intimidated. If you are a physics graduate student, still give this series a read and try not to get too intimidated. I've been working on this series for quite a while, so I hope it's worth my effort.
For further information, take a look at this set of lecture notes on the arXiv, Computational Studies of Quantum Spin Systems by Anders W. Sandvik.
Let's break down each part of that phrase.
At each site, we can have a particle pointing up $| \uparrow \rangle$, down $|\downarrow \rangle$, or some super-position of the two.
Our spin has three degrees of freedom and full $SU(2)$ symmetry. $SU(2)$ is the mathematical group that describes a spin's degrees of freedom. Once we have solved the physics of the Heisenberg case, we also have solved the XY model ($J_z=0$) and the Ising model ($J_x=J_y=0$). As we change the model, we also change the symmetry group first to $SO(2)$ for the XY model and then $\mathbb{Z}_2$ for the Ising model.
Spin only interacts with two neighbors.
Our Hamiltonian has the general form, $$ {\cal H} = \sum_i J_x S_i^x S_{i+1}^x + J_y S_i^y S_{i+1}^y + J_z S_i^z S_{i+1}^z $$ or we can restrict it to the case of $J_x=J_y$ to get the more convenient form, $$ {\cal H} = \sum_i J_{XY} \left(S^+_i S^-_{i+1}+S^-_i S^+_{i+1} \right) +J_z S_i^z S_{i+1}^z. $$ Here $S^{x,y,z}$ are our Pauli operators $$ S^x=\begin{bmatrix} 0 & 1 \\ 1 & 0 \\ \end{bmatrix} \;\;\;\;\;\;\; S^y=\begin{bmatrix} 0 & -i \\ i & 0 \\ \end{bmatrix} \;\;\;\;\;\;\; S^z=\begin{bmatrix} 1 & 0 \\ 0 & -1 \\ \end{bmatrix}, $$ and $S^{\pm}$ are the ladder operators $$ S^+=\frac{S^x+i S^y}{2} = \begin{bmatrix} 0 & 1 \\ 0 & 0\\ \end{bmatrix} \;\;\;\;\;\;\; S^- =\frac{S^x-i S^y}{2}= \begin{bmatrix} 0 & 0 \\ 1 & 0\\ \end{bmatrix}. $$
Assuming we write our basis states in the $S^z_i$ basis, we can divide the terms from the restricted Hamiltonian into on-diagonal and off-diagonal terms. The $S^z_i S^z_{i+1}$ terms compute the magnetization squared, $\vec{S} \cdot \vec{S} $, for a given state and a conserved quantity. These also lie on the diagonal of the matrix corresponding to the Hamiltonian. $$ | \Psi_{i} \rangle = H_{ii} |\Psi_i \rangle = J_{z} \sum_m S_m^z S^z_{m+1} |\Psi_i \rangle $$ The ladder terms, when applied as an operator to the state, create a new state. Thus they act as off-diagonal terms. $$ | \Psi_{j} \rangle = H_{ij} |\Psi_i \rangle = J_{XY} \sum_m \left( S^+_m S^-_{m+1} + S^-_m S^+_{m+1}\right) |\Psi_i \rangle $$
"The model is solvable, not the solution."
--- someone
The Ising, XY, and Heisenberg cases all fall into the special class of problems which have exact solutions in the infinite size limit. Very few exact solutions of quantum mechanics problems exist, so we try and get as much mileage as we can out of the couple of ones we have.
Interestingly enough, once one solution to a problem comes along, someone figures out a different way to approach the same problem.
The 1D Quantum Ising Model is equivalent to the 2D Ising model of classical statistical mechanics, exactly solved in 1944 by Lars Onsager. The solution is also equivalent to a description of free Majorana fermions. [1]
The Jordan-Wigner Transformation solves the 1D XY model by mapping spins to fermions. This transformation only works in special 1D circumstances and the Kitaev Honeycomb model. Since spins possess different anti-commutation relationships than fermions, we attach a string of operators stretching from infinity to each spin. This series of operators changes the relationship between a spin and its neighbors to fermionic. After the transformation, we get a Hamiltonian that is quadratic in the fermionic momentum operators $d_k$, $d^{\dagger}_k$, and we can see the $\cos (ka)$ dispersion relationship for the excitations, $$ {\cal \tilde{H}}_{XY}=-J \sum_k \cos (ka) d^{\dagger}_k d_k. $$ I might write a full article on this later.
Performing a Jordan-Wigner transformation on the full Heisenberg model gives a four-operator scattering term. The Bethe Ansatz, which I honestly don't know anything about, solves the full 1D Heisenberg Model, as well as some 1D Bose gas and Hubbard model problems. Come back to me in many years, or ask a Russian mathematician.
While many people look to Marie Curie, I think Emily Noether is the best female scientist/ mathematician yet. Her theorem is the single most beautiful piece of physics I have ever seen. For each symmetry, there exists a conserved quantity. Space translational symmetry gives us momentum conservation. Time translational symmetry gives us energy conservation.
Conserved quantities save us so much in classical mechanics but save us even more in condensed matter physics. By working out the symmetries of the problem, we can break the huge Hilbert space into smaller chunks of conserved quantities and only have a smaller problem with which to work.
So what are our conserved quantities?
I will only demonstrate magnetization conservation since we won't have to alter our basis to accommodate it.
In [1]:
n=4
nstates=2^n
Out[1]:
Exact Diagonalization is often memory limited. Thus we want to represent our states in the most compact format possible. Luckily, if we are dealing with spin $\frac{1}{2}$, we can just use the 0's ($|\downarrow \rangle$) and 1's ($|\uparrow \rangle$) of the machine. If you are dealing with higher spin, you can use base 3, 4, etc... Part of the reason I needed to create this separate post was to examine working with binary data.
We will keep our states stored as Int, but Julia has operations we can perform to look at the binary format and change the bits.
In [2]:
# psi is an array of all our wavefunctions
psi=convert.(Int8,collect(0:(nstates-1)) )
# Let's look at each state both in binary and base 10
println("binary form \t integer")
for p in psi
println(bitstring(p)[end-n:end],"\t\t ",p)
end
Because we will be using the powers of $2$ frequently in our calculations, we will store all them in an array. They are our placeholders, like $1,10,100,...$
In [3]:
powers2=convert.(Int8, 2.0 .^collect(0:(n-1)) )
println("binary form \t integer")
for p in powers2
println(bitstring(p)[end-n:end],"\t\t ",p)
end
Once we have magnetization for a state (a conserved quantity), we also have magnetization squared for the diagonals.
We could continue to look at the binary format of a number by calling bin, but that converts the number to an array of strings. So instead we want to perform bitwise operations to determine what the binary format looks like in terms of numbers.
Julia supports bitwise not, and, xor (exclusive or), logical shift right, arithmetic shift right, and logical/ arithmetic shift left. For our purposes, we will only be interested in and and xor .
and takes in two inputs and produces one output, given by the following logic table:
| a | b | a&b |
|---|---|---|
| 0 | 0 | 0 |
| 1 | 0 | 0 |
| 0 | 1 | 0 |
| 1 | 1 | 1 |
Julia's & is the bitwise operation and. That means if I combine two numbers, it states the overlap between the two. 1 overlaps with 1; 2 overlaps with 2; 3 overlaps with 2 and 1.
We will use this to compute magnetization.
In [4]:
println(bitstring(1)[end-2:end],"\t",
bitstring(3)[end-2:end],"\t", bitstring(1&3)[end-2:end])
In [5]:
#initializing the magnetization array
m=zeros(length(psi))
println("String \tp&powers2 \tNormed \t\tMagentization")
for i in 1:length(psi)
#Writing the magnetization
m[i]=sum(round.(Int,(psi[i].&powers2)./powers2))
println(bitstring(psi[i])[end-n:end],"\t",psi[i].&powers2,"\t"
,round.(Int,(psi[i].&powers2)./powers2),"\t",m[i])
end
The off diagonal (ladder) elements of the Hamiltonian are the permutation of two neighboring elements in the array. We can permute two indices by combining a mask number with a bitwise XOR ⊻. The symbol ⊻ is written using tab completion for the LaTeX syntax \xor. Otherwise, you can use the function xor( , ).
In [6]:
mask=Int8(3)
testp=1
println("Mask \ttest \tmasked test")
println(bitstring(mask)[end-2:end],'\t',bitstring(testp)[end-2:end]
,'\t',bitstring(testp⊻mask)[end-2:end])
But the mask 3 aka 11 only switches the spins in the first two positions. I need to switch spins in any two adjacent locations. I create this by summing together padded powers of two in order to get the 11 in the correct location.
In [7]:
mask=convert.(Int8,[0;powers2]+[powers2;0])
mask=mask[2:end-1]
println("Mask base10 \tMask Binary \tSummed from")
for i in 1:length(mask)
println(mask[i],"\t\t",bitstring(mask[i])[end-n:end],"\t\t",
bitstring(powers2[i])[end-n:end],"\t",
bitstring(powers2[i+1])[end-n:end])
end
So now let's test how the first of our three masks behaves: We know that if the mask changes a 01 for a 10 (or vice versa) that the overall magnetization will not be changed. So, we test is our mask is successful by comparing the remaining magnetization. The rows offset by two spaces have matching magnetizations.
In [8]:
println("Psi \tPsi \tMasked \t Masked\t \tmPsi \tmMasked")
for p in psi
if m[p+1]==m[p.⊻mask[1]+1]
println(" ",p,"\t ",bitstring(p)[end-n:end],"\t ",
p.⊻mask[1],"\t ",bitstring(p.⊻mask[1])[end-n:end],
"\t\t ",m[p+1],"\t ",m[p.⊻mask[1]+1])
else
println(p,'\t',bitstring(p)[end-n:end],'\t',
p.⊻mask[1],'\t',bitstring(p.⊻mask[1])[end-n:end],
"\t\t",m[p+1],"\t",m[p.⊻mask[1]+1])
end
end
Now let's try the same thing, but for the second mask. This changes the second and third spin.
In [9]:
println("Psi \tPsi \tMasked \t Masked\t \tmPsi \tmMasked")
for p in psi
if m[p+1]==m[p.⊻mask[2]+1]
println(" ",p,"\t ",bitstring(p)[end-n:end],"\t ",
p .⊻ mask[2],"\t ",bitstring(p.⊻mask[2])[end-n:end],
"\t\t ",m[p+1],"\t ",m[p.⊻mask[2]+1])
else
println(p,'\t',bitstring(p)[end-n:end],'\t',
p .⊻ mask[2],'\t',bitstring(p.⊻mask[2])[end-n:end],
"\t\t",m[p+1],"\t",m[p.⊻mask[2]+1])
end
end
That's it for background.
In the next few posts, I'll first cover breaking the matrix up according to the magnetization symmetry, and then the sorting and the searching that accompanies that change. Then we will actually diagonalize the matrix and looks at the results.
In the meantime, go back and check out my post on "Jacobi Transformation of a Symmetric Matrix" for some information on this topic.
Happy physicsing :)
[1] @book{fradkin2013field, title={Field theories of condensed matter physics}, author={Fradkin, Eduardo}, year={2013}, publisher={Cambridge University Press} }
@inproceedings{sandvik2010computational, title={Computational studies of quantum spin systems}, author={Sandvik, Anders W}, booktitle={AIP Conference Proceedings}, volume={1297}, number={1}, pages={135--338}, year={2010}, organization={AIP} }
In [ ]: