Exact Diagonalization for the 1D Heisenberg Chain

9 minute read

The Heisenberg model is given by We choose nearest neighbors to get

We choose a lattice of length \(L\) and lattice spacing \(\Delta\). It is also possible to create a magnetic field (more on why later), using a new term


This will be an overly inefficient, but clear, way to generate the Hamiltonian for exact diagonalization. The steps we’ll take:

  1. Generate the basis
    • this involves finding all binary numbers with the same number of high bits, as this Hamiltonian respects \(S_z\) symmetry
  2. Generate the Hamiltonian in Matrix form
    • For each state in the basis, apply \(H|\psi\rangle\) and record the values/connections
    • This will involve a diagonal term and off diagonal term
    • put them in the appropriate spot in a (dense) matrix
  3. Solve for the matrix’s eigenvalues and eigenvectors
    • put everything together for a nice package!
  4. Analyze the output

Note: A more efficient way of doing this would be to use sparse matrices, but our \(S_z\) block is small enough to not have memory constraints


Helper Functions

These are some nify functions that help with visualization and debugging

import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as spst
def basisVisualizer(L,psi):
    '''Given psi=(#)_10, outputs the state in arrows'''
    #ex: |↓|↑↓|↑|↑|
    psi_2 = bin(psi)[2:]
    N  = len(psi_2)
    up = (L-N)*'0'+psi_2
    configStr = "|"
    uparrow   = '\u2191'
    downarrow = '\u2193'
    for i in range(L):
        blank = True
        if up[i] == '1':
            blank = False
        if up[i] == '0':
            blank = False
        if blank:
        configStr +="|"
def countBits(x):
    '''Counts number of 1s in bin(n)'''
    #From Hacker's Delight, p. 66
    x = x - ((x >> 1) & 0x55555555)
    x = (x & 0x33333333) + ((x >> 2) & 0x33333333)
    x = (x + (x >> 4)) & 0x0F0F0F0F
    x = x + (x >> 8)
    x = x + (x >> 16)
    return x & 0x0000003F 
#helper function to print binary numbers
def binp(num, length=4):
    '''print a binary number without python 0b and appropriate number of zeros'''
    return format(num, '#0{}b'.format(length + 2))[2:]

1. Generate the basis

Because there are a lot of states we’ll need to consider in a given hilbert space, we must be as concise as possible when writing a representation down. Because the computer already uses binary numbers, we’ll take advantage of binary operations and write all states as a binary number. An example of this is below


a state will have up spins as 1 and down spins as 0 since we don’t have double occupency Here is an example of |↑|↑|↓|

print("The binary number for |up up down> is: ",int('110',2))
print("In binary the state is:",binp(int('110',2),3))
print("The state",int('110',2),"can be visualized as:")
#use helper function to visualize
The binary number for |up up down> is:  6
In binary the state is: 110
The state 6 can be visualized as:

we now create a function that will return a list of lists, where each internal list contains all the appropriate binary numbers with a given $S_z$ configuration

def makeSzBasis(L):
    basisSzList = [[] for i in range(0,2*L+1,2)] #S_z can range from -L to L, index that way as well
    #this is probably a bad way to do it
    # count bits is O(log(n)) and loop is O(2**L) :(
    for i in range(2**L):
        Szi = 2*countBits(i) - L
    print("L =",L,"basis size:",2**L)
    return basisSzList

Let’s take a look at the distribution of hilbert space, specifically $S_z$

for L in [7,10]:
    basisSzList = makeSzBasis(L)
    SziVals  = []
    sizeVals = []
    for i in range(len(basisSzList)):
plt.ylabel("# of elements",fontsize=15)
L = 7 basis size: 128
L = 10 basis size: 1024


As expected, we see there are more ways to make \(S_z=0\) than anything else, and the distribution looks normal in behaivor. The problem with exact diagonalization is as we increase system size, there are exponentially more states we need to keep track of

2. Generate the Hamiltonian

To create a dense matrix of the Hamiltonian we’ll need two things. A lookup table that converts an element of the hilbert space (like 6) to its index in the Hamiltonan matrix. Then we’ll also need something that takes each state and spits out all the states its connected to, and the values of this connection. This is added into the Hamiltonian matrix. Here, I’ve hard coded it in but more general techniques for both can be used

def makeH(SzList,L,Jxy,Jz,h):
    '''Make a 1D Heisenberg chain of length L with Jxy,Jz and magnetic field h out of an SzList of states'''
    basisMap = {}
    stateID  = 0 
    #generate an ordering
    for state in SzList:
        #print(state) #basisVisualizer(L,state)
        basisMap[state] = stateID
    nH = stateID
    H = np.zeros([nH,nH])
    #now fill H
    for state in SzList:
        idxA = basisMap[state]
        H[idxA,idxA] += -h*countBits(state) # B field
        for i in range(L):
            j = (i+1)%L #nearest neighbors are hard coded here
            if (state>>i & 1) == (state>>j & 1):#matching bit check
                H[idxA,idxA] += -Jz/4
                H[idxA,idxA] -= -Jz/4
                mask = 2**(i)+2**j
                stateB= state^mask #this flips the bits at i,j
                idxB = basisMap[stateB]
                H[idxA,idxB]+= -Jxy/2
    #print(np.all(H==H.T)) #check if Hermitian and is coded properly; very slow
    return H

3. Solve for the eigenvalues and eigenvectors

After making \(H\) for each block in \(S_z\), we’ll record all the energies and where the ground state is

def getSpectrum(L,Jxy,Jz,h):
    '''Returns lowestEnergy, 
               Sz sector of the GS, 
               GS eigenvector, 
               and all energies'''
    basisSzList  = makeSzBasis(L)
    energies     = []
    lowestEnergy = 1e10
    for Szi,SzList in enumerate(basisSzList):
        #print("Sz sector:",-L+2*Szi)
        H     = makeH(SzList,L,Jxy,Jz,h)
        lam,v = np.linalg.eigh(H)
        #keep track of GS
        if min(lam) < lowestEnergy:
            lowestEnergy  = min(lam)
            GSSector      = -L+2*Szi
            GSEigenvector = v[:,lam.argmin()]
    print("Energies assembled!")
    print("Lowest energy:",lowestEnergy)
    print("The ground state occured in Sz=",GSSector)
    return (lowestEnergy,GSSector,GSEigenvector,energies)
for i in energies:
L = 4 basis size: 16
Energies assembled!
Lowest energy: -1.9999999999999996
The ground state occured in Sz= 0
[-1.00000000e+00 -2.25514052e-17  0.00000000e+00  1.00000000e+00]
[-2.00000000e+00 -1.00000000e+00 -6.19908411e-17  0.00000000e+00
  6.84227766e-49  1.00000000e+00]
[-1.00000000e+00 -2.25514052e-17  0.00000000e+00  1.00000000e+00]

Plotting routines:

def generatePlots(L,Jxy,Jz,h):
     GSEigenvector,energies) = getSpectrum(L,Jxy,Jz,h)
    total_energies = [en for szlist in energies for en in szlist]
    maxE           = np.max(total_energies)
    offset = 0
    for i in range(len(energies)):
        if len(energies)-4>i>2:
            if i%2==0:

    plt.xlabel("Arbitrary Order",fontsize=15)
    plt.title(r"XXZ model with $L="+str(L)+",\,\,\, J_z="+str(Jz)+",\,\,\, J_{xy}="+str(Jxy)+"$",fontsize=15)
    plt.plot([0,offset],[lowestEnergy,lowestEnergy],'--',label="Ground State")
    plt.legend(loc='lower right')
    plt.xlabel("state order",fontsize=15)
    plt.title("Ground State Eigenvector \nwith $L="+str(L)+",\,\,\, J_z="+str(Jz)+",\,\,\, J_{xy}="+str(Jxy)+"$",fontsize=15)
    basisSzList   = makeSzBasis(L)
    GSEigenvector = np.abs(np.round(GSEigenvector,10)) #rounding errors
    bigStatesID   = np.argwhere(np.abs(GSEigenvector) == np.max((GSEigenvector))).reshape((1,-1))
    #Get the states
    print("The biggest states are:")
    for state in bigStatesID[0]:
        bigStates = basisSzList[(GSSector+L)//2][state]

4. Physics

We’ll look at two systems, an antiferromagnet and a ferromagnet


#Constants of the problem
L   = 12
Jxy = 1     # when J<0 it is antiferromagnetic 
Jz  = 1     # J_z = J_xy = J is the normal Heisenberg model instead of XXZ model
h   = -1e-5 # small h field to break degeneracy
L = 12 basis size: 4096
Energies assembled!
Lowest energy: -3.0
The ground state occured in Sz= -12




L = 12 basis size: 4096
The biggest states are:

Ground State Eigenvector

This is exactly the maximal ferromagnetic state! We added the tiny magnetic field so that we could easily choose the up one over the down one, something that would naturally happen in nature but is hard to get on the computer. There are also degenerate product states in every \(S_z\) sector that are now gapped by the magnetic field.


All we have to do is swap the sign on \(J_z,J_{xy}\)

#Constants of the problem
L   = 12
Jxy = -1
Jz  = -1
h   = 0
L = 12 basis size: 4096
Energies assembled!
Lowest energy: -5.387390917445205
The ground state occured in Sz= 0




L = 12 basis size: 4096
The biggest states are:
CPU times: user 2.73 s, sys: 191 ms, total: 2.93 s
Wall time: 2.24 s

Notice how the maximal states of the ground state are the most antiferromagnetic states, but there is no difference to the model over the first and second one.

There’s also lower and lower energy, the lower total \(|S_z|\) is, which is expected as we want the least ferromagnetic configuration we can create.

Anderson Tower of States

Secretly, in the thermodynamic limit there is sponanteous symmetry breaking with the antiferromagnetic ground state, along with gapless goldstone boson modes that come with it. Exact Diagonalization cannot see these as the system size is finite, however we can see the “hidden” effects of it.

Given an eigenstate from our computed spectrum, we’ll need to calculate the total spin expectation value \(\langle S^2\rangle\) or

def calcTotalS(L,stateList,basisMap,eigvector):
    '''Given an eigenvector, calculate <S^2>=S*(S+1)'''
    totalS = 0
    for i in range(len(eigvector)):
        currState = stateList[i]
        amp       = eigvector[i]
        for ix in range(L):
            for jx in range(L):
                if (currState>>ix & 1) == (currState>>jx & 1):
                    totalS += 0.25*amp**2 #Sz, 1/4 = 1/2*1/2
                    if ix == jx: totalS+= 0.5*amp**2
                    totalS -=0.25*amp**2 #Sz -1/4 = +1/2 * -1/2
                    mask      = 2**(ix)+2**jx
                    flipState = currState^mask #this flips the bits at i,j
                    idxFlip   = basisMap[flipState]
                    totalS   += 0.5*amp*eigvector[idxFlip]
    return totalS
#unit tests
if np.abs(calcTotalS(2,[3],{3:0},[1]) - 1*(1+1)) > 1e-6:
    raise ValueError("Bad <S^2>")
if np.abs(calcTotalS(2,[0],{0:0},[1]) - 1*(1+1)) > 1e-6:
    raise ValueError("Bad <S^2>")
if np.abs(calcTotalS(2,[1,2],{1:0,2:1},[1/np.sqrt(2),-1/np.sqrt(2)]) - 0*(0+1)) > 1e-6:
    raise ValueError("Bad <S^2>")
if np.abs(calcTotalS(2,[1,2],{1:0,2:1},[1/np.sqrt(2),1/np.sqrt(2)]) - 1*(1+1)) > 1e-6:
    raise ValueError("Bad <S^2>")
Jxy = -1
Jz  = -1
h   = 0
minEList = []
minEList_goldstone = []
minSList = []
LList = [8,10,12,14] 
for L in LList:
    basisSzList = makeSzBasis(L)
    SList = [] #this will hold all the S^2 values and what their energies are
    EList = []
    for sZ in range(len(basisSzList)//2+1):#exploit the Sz symmetry
        SzList   = basisSzList[sZ]
        basisMap = {}
        stateID  = 0 
        #generate the ordering
        for state in SzList:
            basisMap[state] = stateID
        H     = makeH(SzList,L,Jxy,Jz,h)
        lam,v = np.linalg.eigh(H)
        numStates = min(len(lam),10) #eigenvalues are sorted smallest to largest
        EList     = EList+list(lam[:numStates])
        for i in range(0,numStates):
    SList,EList = zip(*sorted(zip(SList,EList))) #sort them by S^2
    #now find all the S^2 values and determine the minimum
    idx   = 1
    start = 0
    end   = 1
    uniqueS           = [SList[0]]
    uniqueE           = []
    uniqueE_goldstone = []
    while idx!=len(SList):
        if uniqueS[-1] != SList[idx]:
            end = idx+1
            ESort = sorted(EList[start:end])
            if len(ESort)>1:
    ESort = sorted(EList[start:])
    if len(ESort)>1:
for i in range(len(minEList)):
plt.ylabel("Lowest Energy",fontsize=15)
plt.title("Anderson Tower of States",fontsize=15)
for i in range(len(minEList)):
    plt.plot(minSList[i],[minEList[i][0]-en for en in minEList[i]],'-o',label="L="+str(LList[i]))
plt.ylabel("Difference\nfrom GS",fontsize=15)

#plot the slopes
slopes = []
for i in range(len(minEList)):
    slope, intercept, r_value, p_value, std_err = spst.linregress(minSList[i][1:],minEList[i][1:])
plt.plot([1/L for L in LList],slopes,'o-')
L = 8 basis size: 256
L = 10 basis size: 1024
L = 12 basis size: 4096
L = 14 basis size: 16384



CPU times: user 1min 18s, sys: 2.15 s, total: 1min 20s
Wall time: 1min 14s

What is going on? We see that the degenerate states that occur in the thermodynamic limit are slowly coming down. How fast are they coming down? \(1/L\), which is shown in the second plot.

%load_ext watermark

CPython 3.7.0
IPython 6.4.0

compiler   : Clang 9.0.0 (clang-900.0.39.2)
system     : Darwin
release    : 18.2.0
machine    : x86_64
processor  : i386
CPU cores  : 4
interpreter: 64bit