Exact diagonalization : lecture 1 ¶
Introduction : Quantum spins systems¶
Calculations and graphs on the blackboard
Reminders on spin quantization¶
Spin algebra $\hbar=1$¶
Commutation relations : $$ \vec{S}\wedge\vec{S} = i\vec{S} $$ Quantum numbers and states $|{m,S}\rangle$: $$ \vec{S}^2|{m,S}\rangle = S(S+1)|{m,S}\rangle \qquad S^z|{m,S}\rangle = m|{m,S}\rangle\qquad -S\leq m\leq S$$ Ladder operators : $$ S^{\pm} = S^x \pm iS^y\qquad S^{\pm}|{m,S}\rangle = \sqrt{S(S+1)-m(m\pm 1)} |{m\pm 1,S}\rangle $$
Example : spin $1/2$¶
In the $\{|{\uparrow}\rangle,|{\downarrow}\rangle\}$ basis : $$ S^z = \frac12 \begin{pmatrix} {1} & {0} \\ {0} & {-1} \\ \end{pmatrix} \quad S^x = \frac12 \begin{pmatrix} {0} & {1} \\ {1} & {0} \\ \end{pmatrix} \quad S^y = \frac12 \begin{pmatrix} {0} & {-i} \\ {i} & {0} \\ \end{pmatrix}$$ The ladder operators read $$ S^+ = \begin{pmatrix} {0} & {1} \\ {0} & {0} \\ \end{pmatrix} \quad S^- = \begin{pmatrix} {0} & {0} \\ {1} & {0} \\ \end{pmatrix}$$
Diagonal terms : "classical" energy
Off-diagonal terms : "quantum fluctuations" or "tunneling" elements, generate superposition of states and entanglement between the classical states.
Two-body interaction between spins¶
Ising part : $J_z S^z_1S^z_2$, diagonal term, gives $\pm 1/4$, $J_z>0$ is antiferromagnetic, $J_z<0$ is ferromagnetic.
Transverse part : $J_{xy}[S^x_1S^x_2+S^y_1S^y_2] = \frac{J_{xy}}{2}[S^+_1S^-_2+S^-_1S^+_2]$ induces spin-flip and quantum fluctuations.
XXZ model : combines both the Ising and spin-flip parts, usually written $$ \Delta S^z_1S^z_2 + \frac{1}{2}[S^+_1S^-_2+S^-_1S^+_2] $$
Heisenberg model : it corresponds to the $\Delta=1$ limit of the XXZ model, so that the coupling reads $$ J \vec{S}_1\cdot\vec{S}_2 $$ with $J>0$ for antiferro and $J<0$ for ferromagnetic. The model is clearly invariant under rotation symmetry in the spin space (scalar product). This corresponds to a global SU(2) symmetry so that the total spin is conserved. For two spins, the total spin is $\vec{S}_{tot} = \vec{S}_1 + \vec{S}_2$. As $ 2\vec{S}_1\cdot\vec{S}_2 = \vec{S}_{tot}^2 - \vec{S}_1^2 - \vec{S}_2^2$, which eigenvalues are $S_{tot}(S_{tot}+1) - 2S(S+1)$.
Example with spin-$1/2$¶
In this case, either $S_{tot} = 0$ (singlet state) or $S_{tot}=1$ (triplet states). Thus, the eigenvalues of the Heisenberg coupling are $-3J/4$ (singlet) and $J/4$ (triplets).
We can derive the eigenstates by simply writting the coupling in the basis $\{|{\uparrow\uparrow}\rangle,|{\uparrow\downarrow}\rangle,|{\downarrow\uparrow}\rangle,|{\downarrow\downarrow}\rangle\}$ $$ J \vec{S}_1\cdot\vec{S}_2 = J \begin{pmatrix} {1/4} & {0} & {0} & {0} \\ {0} & {-1/4}& {1/2} & {0} \\ {0} & {1/2}& {-1/4} & {0} \\ {0} & {0}& {0} & {1/4} \end{pmatrix} $$ The eigenstates are clearly $$ \frac{1}{\sqrt{2}}[|{\uparrow\downarrow}\rangle - |{\downarrow\uparrow}\rangle] \quad\text{(singlet)}$$ and $$ \frac{1}{\sqrt{2}}[|{\uparrow\downarrow}\rangle + |{\downarrow\uparrow}\rangle],\quad|{\uparrow\uparrow}\rangle ,\quad|{\downarrow\downarrow}\rangle \quad\text{(triplets)}$$ Notice that the Néel state $|{\uparrow\downarrow\uparrow\downarrow\uparrow\cdots}\rangle$ is not an eigenstate of the Heisenberg coupling while the ferromagnetic state $|{\uparrow\uparrow\uparrow\uparrow\cdots}\rangle$ is.
External magnetic field¶
An external magnetic field couples via Zeeman effect to the spins. If all physical constants (gyromagnetic factor, Bohr magneton) are included in the definition of the field $\vec{h}$, this coupling reads $$ -\vec{h}\cdot \vec{S} = -h_xS^x - h_yS^y - h_zS^z$$ $h_z$ is a longitudinal field (diagonal) while $h_{x,y}$ are transverse fields (off-diagonal).
The Ising in a transverse field (ITF) model¶
It combines a competition between the Ising term (here taken ferromagnetic with $J>0$ in the notations), and a transverse field along $x$, written $h$ : $$ H = - J \sum_{\langle{i,j}\rangle} S_i^zS_j^z -h\sum_i S_i^x $$
Quantum phases transitions¶
Qualitative temperature versus transverse external field diagram for a ferromagnet.
Behavior of the order parameter. Limiting ground-states at zero temperature.
Mean-field (or classical) analysis of the 1D ITF model, existence of a critical point. One obtains $h^c = 2JS$ ($h^c = J$ if extrapolated to a spin-$1/2$) and $$ \langle{S^x}\rangle = Sh/h^c \qquad \langle{S^z}\rangle = S\sqrt{1-(h/h^c)^2}$$
Qualitative picture for the correlation functions and the gap across the transition.
A quick start¶
Binary operators¶
bit representation of spin configurations in the $S^z$ basis
conf = int("100110",base=2)
print(conf)
conf = 0b100110
print(conf,bin(conf))
print("XOR operator")
for i in [0,1]:
for j in [0,1]:
print(str(i)+"^"+str(j),":",i^j)
print("OR operator")
for i in [0,1]:
for j in [0,1]:
print(str(i)+"|"+str(j),":",i|j)
print("AND operator")
for i in [0,1]:
for j in [0,1]:
print(str(i)+"&"+str(j),":",i&j)
import numpy as np
print("NOT operator")
for i in [0,1]:
print("~"+str(i),":",np.binary_repr(~i,64))
print()
L=20
site9 = 2**9 # tenth from the right, count starts at 0
block9 = 2**9-1
print(np.binary_repr(site9,L))
print(np.binary_repr(~site9,L))
print(np.binary_repr(block9,L))
print(np.binary_repr(1<<9,L))
print(np.binary_repr(site9>>9,L))
print(np.binary_repr(block9<<6,L))
A first implementation¶
We want to find the ground-state of the Ising in a transverse field Hamiltonian $$ H = -J\sum_{i=1}^L S_i^zS_{i+1}^z - h\sum_i S_i^x $$ we represent each conf by a binary number from 0 to $2^L-1$, each conf can be its own index in the Hilbert space.
def binconf(c): return np.binary_repr(c,L)
def readsite(conf,i): return (conf&(1<<i))>>i
L = 3
hilbertsize = 2**L
for conf in range(hilbertsize):
print(binconf(conf))
def SzSz(conf,i,j):
si = readsite(conf,i)-1/2
sj = readsite(conf,j)-1/2
return si*sj
def IsingHam(conf):
return sum([SzSz(conf,i,(i+1)%L) for i in range(L) ])
for conf in range(hilbertsize):
print(IsingHam(conf))
def Sx(conf,i):
return 1/2, conf^(1<<i)
for conf in range(hilbertsize):
print(binconf(conf),binconf(Sx(conf,0)[1]))
def Spinflip(conf,i,j):
"SxSx + SySy = 0.5(S+S- + S-S+) term"
if readsite(conf,i) != readsite(conf,j):
return 0.5, conf^((1<<i)^(1<<j))
else: return 0.0, conf
for conf in range(hilbertsize):
print(binconf(conf),binconf(Spinflip(conf,0,1)[1]))
hx = 1.0
# diagonal part
Ham = np.diag([IsingHam(conf) for conf in range(hilbertsize)])
# off-diagonal part
for conf in range(hilbertsize):
for i in range(L):
value, newconf = Sx(conf,i)
Ham[newconf,conf] -= hx*value
print(Ham)
np.set_printoptions(precision=2,suppress=True)
En, Phin = np.linalg.eigh(Ham)
for i in range(len(En)):
print(round(En[i],5),"\t",Phin[:,i])
Object oriented programming implementation¶
So far, we only used functions and builtins or library objects and tools. In order to make life easier, we would like to build quantum mechanics objects in order to manipulate them as we do on the paper. Therefore, we use object oriented programming which is quite easy in Python.
Basics of classes of overloading of operators¶
classes, instance, overloading of operators, inheritance.
The goal¶
is to be able to write user friendly scripts that are easy to understand and allows us to write algorithms that directly manipulate the objects. Here is an example of what we would like to write :
from quantumising import *
L, hx = 8, 0.3
HamH = GeneralOp([ [ -hx, Sx(i)] for i in range(L) ])
HamIsing = GeneralOp([ [ 1.0, SzSz(i,(i+1)%L)] for i in range(L) ])
Ham = HamIsing + HamH
print(Ham)
hilbert = IsingHilbert(L=L)
psi0 = Ham.Groundstate(hilbert)
print(psi0.measure(Ham), psi0.measure(Sx(0)), psi0.measure(SzSz(0,L//2)))
This may look like a waste of time for the Ising model but the same structure can be reused for more complicated Hilbert spaces encountered later.
Hilbert space¶
First an helper class which only provides a namespace for functions manipulating configurations. Using namespace prevents one from overwriting functions. The methods here are class methods and not instance method. Python automatically builds help from functions documentation !
import sys
class Conf:
"A namespace for functions on configurations"
def binconf(c,L):
"binary representation of a conf integer"
return np.binary_repr(c,L)
def ket(conf,L):
"quantum ket formatting with first site on the left"
sconf = Conf.binconf(conf,L)
return "|"+"".join([s for s in sconf[::-1]])+">"
def readsite(conf,i):
"reads quantum number on site i"
return (conf&(1<<i))>>i
help(Conf)
Now a class for Ising hilbertspace. Hilbert spaces are meant to be instances.
class IsingHilbert(object):
"Hilbert space for Ising in a transverse field model"
def __init__(self, L=1):
"description is a list containing the local hilbert space at site i"
self.L = L
self.hilbertsize = 2**self.L
# no need for that here, because each conf is its own index for Ising
# but we introduce it for latter purpose
self.Confs = { i:i for i in range(self.hilbertsize) }
def __repr__(self):
"prints the first 30 confs of hilbert space"
from itertools import islice
num = min(30,self.hilbertsize)
s = "Printing first "+str(num)+" configurations among "+str(self.hilbertsize)+"\n"
s += ", ".join([Conf.ket(c,self.L) for c in list(islice(self.Confs.keys(),num))])
return s
def index(self,conf):
"returns the index of a configuration"
return self.Confs[conf]
hilbert = IsingHilbert(L=6)
print(hilbert)
Wave function¶
A class describing a wave-function. It's an array of complex (or real) numbers, but one needs the meaning of indices, ie. must keep a link to the Hilbert space. In order to benefit from all builtins numpy arrays functions, we use inheritence from the np.ndarray object.
import numpy as np
def cpp_convert(x,p):
if type(x) is complex or type(x) is np.complex64 or type(x) is np.complex128:
return "("+str(round(x.real,p))+","+str(round(x.imag,p))+")"
else:
return str(round(x,p))
class Wavefunction(np.ndarray):
def __new__(cls,hilbert,dtype=float):
obj = np.ndarray.__new__(cls,(hilbert.hilbertsize,),dtype)
obj.hilbert = hilbert
return obj
def __array_finalize__(self, obj):
if obj is None: return
self.hilbert = getattr(obj,'hilbert',None)
def setvalue(self, conf, value):
"Sets the coefficient in front of conf"
self[self.hilbert.index(conf)] = value
def getvalue(self, conf):
"gets the coefficient in front of conf"
return self[self.hilbert.index(conf)]
def readable(self,precision=4,eps=1e-4):
from itertools import islice
num = min(30,self.hilbert.hilbertsize)
L = self.hilbert.L
s = ""
for conf,i in islice(self.hilbert.Confs.items(),num):
coef = self[i]
if abs(coef)>eps:
symb = " + "
if type(coef) is np.float64 and np.sign(coef) == -1.0:
symb = " - "
coef = abs(coef)
s += symb + cpp_convert(coef,precision) + Conf.ket(conf,L)
if s[:3] == " + ": s = s[3:]
return "|psi> = "+s
def norm(self):
return np.linalg.norm(self)
def normalize(self):
norm = self.norm()
try: self /= norm
except: print("can't normalize a vector of norm", norm)
def randomize(self,a=-1.0,b=1.0):
self[:] = (b-a)*np.random.random_sample(self.size).reshape(self.shape)+a
if self.dtype == np.complex128:
self.imag = (b-a)*np.random.random_sample(self.size).reshape(self.shape)+a
self.normalize()
def measure(self,op,hermitic=True):
"computes mean-value of an hermitic operator"
res = scalar(self,op*self)
if abs(res.imag) > 1e-15:
print("Warning, issue with hermiticity")
return res.real
def scalar(psi1,psi2):
"computes overlaps of psi1 and psi2 with complex conjugation on psi1"
return np.vdot(psi1,psi2)
def zeroWf(hilbert,dtype=float):
psi = Wavefunction(hilbert,dtype=dtype)
psi[:] = 0.0
return psi
def basisvectorWf(hilbert,index,dtype=float):
psi = zeroWf(hilbert,dtype=dtype)
psi[index] = 1.0
return psi
def randomWf(hilbert,dtype=float):
psi = Wavefunction(hilbert,dtype=dtype)
psi.randomize()
return psi
def copy(psi):
res = Wavefunction(psi.hilbert,dtype=psi.dtype)
res[:] = psi[:]
return res
hilbert = IsingHilbert(L=2)
psi = randomWf(hilbert)
print(psi)
print(psi.readable(precision=2))
psi1 = randomWf(hilbert)
psi2 = randomWf(hilbert)
print(psi1.readable(),psi2.readable(),sep="\n")
psi3 = Wavefunction(hilbert)
psi3[:] = 0.3*(psi1 + 0.1*psi2)
print(psi3.readable())
Operators¶
Now, we would like to build operators and be able to apply them to wave-functions.
class Operator(object):
def __mul__(self,psi):
res = zeroWf(psi.hilbert,dtype=psi.dtype)
for conf, i in psi.hilbert.Confs.items():
coef, newconf = self.Apply(conf)
j = psi.hilbert.Confs[newconf]
res[j] += coef * psi[i]
return res
def __str__(self):
return self.name
Specialized to spin $1/2$
class Sz(Operator):
def __init__(self,site):
self.site = site
self.name = "Sz_"+str(site)
def Apply(self,conf):
return Conf.readsite(conf,self.site)-1/2, conf
class SzSz(Operator):
def __init__(self,s1,s2):
self.s1 = s1
self.s2 = s2
self.name = "Sz_"+str(s1)+"Sz_"+str(s2)
def Apply(self,conf):
si = Conf.readsite(conf,self.s1)-1/2
sj = Conf.readsite(conf,self.s2)-1/2
return si*sj, conf
class Sx(Operator):
def __init__(self,site):
self.site = site
self.name = "Sx_"+str(site)
def Apply(self,conf):
return 1/2, conf^(1<<self.site)
class Spinflip(Operator):
def __init__(self,i,j):
self.i = i
self.j = j
self.name = "Sx_"+str(i)+"Sx_"+str(j)+"Sy_"+str(i)+"Sy_"+str(j)
def Apply(self,conf):
if Conf.readsite(conf,self.i) != Conf.readsite(conf,self.j):
return 0.5, conf^((1<<self.i)^(1<<self.j))
else: return 0.0, conf
hilbert = IsingHilbert(L=2)
singlet = zeroWf(hilbert)
singlet.setvalue(int("10",2), 1/np.sqrt(2))
singlet.setvalue(int("01",2),-1/np.sqrt(2))
print(singlet.readable())
print(singlet.measure(Sz(0)), singlet.measure(Sz(1)))
print(singlet.measure(SzSz(0,1)), singlet.measure(Spinflip(0,1)))
triplet = zeroWf(hilbert)
triplet.setvalue(int("10",2), 1/np.sqrt(2))
triplet.setvalue(int("01",2), 1/np.sqrt(2))
print(triplet.measure(SzSz(0,1)), triplet.measure(Spinflip(0,1)))
class GeneralOp(list):
"""
A list of 2-element in the format [ [coef, Operator] ]
corresponding to \sum_<i,j> coef_ij * Op_ij
warning: assumes hermitic operators
"""
def __init__(self,iterable=[]):
list.__init__(self,iterable)
self.having_matrix = False
def __str__(self):
if not len(self): return ""
s = "("+str(self[0][0])+")*"+str(self[0][1])
for coef, op in self[1:]:
s += "+("+str(coef)+")*"+str(op)
return s
def __iadd__(self,other):
return GeneralOp(self[:]+other[:])
def __add__(self,other):
return GeneralOp(self[:]+other[:])
def __mul__(self,psi):
res = copy(psi)
if self.having_matrix:
res[:] = np.dot(self.matrix,psi)
else:
res[:] = 0.0
for coef, op in self:
res += coef * (op*psi)
return res
def createMatrix(self,hilbert,dtype=float):
if self.having_matrix: return
if hilbert.hilbertsize < 2: return
self.createFullMatrix(hilbert,dtype)
self.having_matrix = True
def createFullMatrix(self,hilbert,dtype=float):
size = hilbert.hilbertsize
self.matrix = np.zeros(shape=(size,size), dtype = dtype)
for index in np.arange(size):
self.matrix[:,index] = self*basisvectorWf(hilbert,index)
def diagonalize(self,with_eigvec=True):
"Assumes an hermitic operator by default"
if not self.having_matrix:
print("Error, don't have the matrix")
return None
else:
if with_eigvec:
self.En, self.Phin = np.linalg.eigh(self.matrix)
else:
self.En = np.linalg.eigvalsh(self.matrix)
def Groundstate(self,hilbert,dtype=float):
"computes the groundstate and returns it"
if not self.having_matrix: self.createMatrix(hilbert,dtype)
if not hasattr(self,"Phin"): self.diagonalize()
gsIndex = np.argmin(self.En)
gs = Wavefunction(hilbert,dtype=dtype)
gs[:] = self.Phin[:, gsIndex]
return gs
Finding the ground-state and doing basic measurements¶
L, h = 3, 0.3
HamH = GeneralOp([ [ -h, Sx(i)] for i in range(L) ])
HamIsing = GeneralOp([ [ -1.0, SzSz(i,(i+1)%L)] for i in range(L) ])
Ham = HamIsing + HamH
print(Ham)
hilbert = IsingHilbert(L=L)
psi = randomWf(hilbert)
print((Ham*psi).readable())
psi0 = Ham.Groundstate(hilbert)
print(psi0.measure(Ham), psi0.measure(Sx(0)), psi0.measure(SzSz(0,L//2)))
Use a module¶
Now we can also put all the class definitions in a single quantumising.py file that we import at the begining of the script.
Important : if you are using the IPython notebook and have imported a module, and if you modify externally the module, the new version will not be reloaded in the next call even though you have called the from module import * command. A inelegant solution is to restart the Kernel. Other solutions can be found on the IPython documentation. Yet, prefer testing using online python command and include your code in the notebook when ready.
Downloading the module¶
from quantumising import *
L, h = 8, 3
HamH = GeneralOp([ [ -h, Sx(i)] for i in range(L) ])
HamIsing = GeneralOp([ [ -1.0, SzSz(i,(i+1)%L)] for i in range(L) ])
Ham = HamIsing + HamH
print(Ham)
hilbert = IsingHilbert(L=L)
psi0 = Ham.Groundstate(hilbert)
print(psi0.measure(Ham), psi0.measure(Sx(0)))
from pylab import *
plot(list(range(L)),[psi0.measure(SzSz(0,r)) for r in range(L)],'-o')
show()
