MC-06 QWL

The goal of this tutorial is to compute thermodynamic properties of quantum spin systems using the quantum Wang-Landau (QWL) algorithm. Unlike path-integral QMC methods that simulate at a fixed temperature, QWL performs a random walk in the space of expansion orders of the partition function to estimate the density of states directly. From a single simulation one can then evaluate thermodynamic observables — energy, specific heat, entropy, susceptibility — over the full temperature range without re-running.

Thermodynamics of quantum Heisenberg spin chains

The ferromagnetic Heisenberg chain

The parameter file parm6a sets up a QWL simulation of the $S=\frac{1}{2}$ Heisenberg ferromagnet on a chain of 40 sites:

LATTICE="chain lattice"
MODEL="spin"
local_S = 1/2
L       = 40
CUTOFF  = 500
{J = -1}

CUTOFF sets the maximum expansion order of the partition function sampled by the algorithm; 500 is sufficient for a 40-site chain over the temperature range studied here.

Command line

After preparing the input files and running the qwl code:

parameter2xml parm6a
qwl parm6a.in.xml

produce XML plot files for all thermodynamic observables using qwl_evaluate:

qwl_evaluate --T_MIN 0.1 --T_MAX 10 --DELTA_T 0.1 parm6a.task1.out.xml

This generates one XML file per observable:

parm6a.task1.plot.energy.xml
parm6a.task1.plot.free_energy.xml
parm6a.task1.plot.entropy.xml
parm6a.task1.plot.specific_heat.xml
parm6a.task1.plot.uniform_structure_factor.xml
parm6a.task1.plot.staggered_structure_factor.xml
parm6a.task1.plot.uniform_susceptibility.xml

Extract the data to plain text with plot2text:

plot2text parm6a.task1.plot.energy.xml

The tools plot2xmgr (Grace) and plot2gp (Gnuplot) produce equivalent output in those formats. The preferred method for analysis is Python, described next.

Python

The script tutorial6a.py sets up and runs the simulation, then evaluates all observables in one call:

import pyalps
import matplotlib.pyplot as plt

data = pyalps.evaluateQWL(pyalps.getResultFiles(prefix='parm6a'),
                          DELTA_T=0.1, T_MIN=0.1, T_MAX=10.0)

for s in pyalps.flatten(data):
    plt.figure()
    plt.title("Ferromagnetic Heisenberg chain")
    pyalps.plot.plot(s)
    plt.show()

For the ferromagnet ($J=-1$) you should see a broad specific-heat peak at low temperature (a Schottky-like anomaly) and a uniform susceptibility that diverges as $T\to 0$, consistent with ferromagnetic order developing at zero temperature in 1D.

The antiferromagnetic Heisenberg chain

To simulate the antiferromagnetic chain set $J=1$ instead of $J=-1$. The parameters are in parm6b and the Python script in tutorial6b.py; everything else is identical to the ferromagnetic case.

For the antiferromagnet the uniform susceptibility remains finite as $T\to 0$ (a signature of the spin-liquid ground state of the 1D antiferromagnet), while the specific-heat peak shifts and broadens differently.

The three-dimensional Heisenberg antiferromagnet

Simulating the 3D quantum Heisenberg antiferromagnet

The parameter file parm6c sets up a QWL simulation of the $S=\frac{1}{2}$ Heisenberg antiferromagnet on a simple cubic lattice with $4^3=64$ sites. The Python script is tutorial6c.py. Run and evaluate exactly as for the chain above.

The staggered structure factor $S(\pi,\pi,\pi)$ should start rising steeply near $T\approx 1$, signalling the onset of antiferromagnetic correlations. The specific heat shows a corresponding peak and the entropy decreases rapidly in the same temperature range.

Finite-size scaling analysis to determine the critical point

Finite-size scaling predicts that the staggered structure factor at the critical point scales as $S(L) \propto L^{2-\eta}$ with $\eta \approx 0.034$ (3D classical Heisenberg universality class). A plot of $S(L)/L^{2-\eta}$ vs.\ temperature should show curves for different $L$ crossing at the critical temperature $T_c$.

The parameter file parm6d (or tutorial6d.py) runs the cubic antiferromagnet for two system sizes ($L=4$ and $L=6$) with a larger CUTOFF=1000 to maintain accuracy at lower temperatures. After running, load the results:

import pyalps
import matplotlib.pyplot as plt
import copy

results = pyalps.evaluateQWL(pyalps.getResultFiles(prefix='parm6d'),
                             DELTA_T=0.05, T_MIN=0.5, T_MAX=1.5)

Extract the staggered structure factor for each system size, rescale by $L^{-(2-\eta)}$, and label by $L$:

eta = 0.034
data = []
for s in pyalps.flatten(results):
    if s.props['ylabel'] == 'Staggered Structure Factor per Site':
        d = copy.deepcopy(s)
        l = s.props['L']
        d.props['label'] = 'L=' + str(l)
        d.y = d.y * pow(float(l), -(2 - eta))  # rescale by L^{-(2-eta)}
        data.append(d)

Then plot the scaling curves:

plt.figure()
plt.title("Scaling plot for cubic lattice Heisenberg antiferromagnet")
pyalps.plot.plot(data)
plt.legend()
plt.xlabel('Temperature $T/J$')
plt.ylabel('$S(\pi,\pi,\pi)\, L^{-(2-\eta)}$')
plt.show()

The curves for different $L$ should cross near $T_c \approx 0.946$.

Questions

• Compare the specific heat and uniform susceptibility of the ferromagnetic ($J=-1$) and antiferromagnetic ($J=1$) chains. Where are the differences most pronounced, and why are the two cases nearly identical at high temperature?
• What is the entropy of the chain at $T=0$ and as $T\to\infty$? (If the $L=40$ data are noisy at low $T$, repeat with $L=8$.) Is this consistent with the third law of thermodynamics?
• Why does the uniform susceptibility diverge as $T\to 0$ for the ferromagnet but remain finite for the antiferromagnet?
• Why does the staggered structure factor of the 3D antiferromagnet start increasing near $T\approx 1$? What other thermodynamic signatures accompany this?
• Do the rescaled structure-factor curves $S(L)\,L^{-(2-\eta)}$ cross at a single temperature? What is your estimate of $T_c$? Compare with the literature value $T_c = 0.946$.
• How could you obtain a more precise estimate of $T_c$? What would you need to change in the simulation?
• Would you expect the critical temperature of the quantum Heisenberg ferromagnet on a cubic lattice to be the same as for the antiferromagnet? How would you check?