MC-08 Quantum Phase Transition
In this tutorial we will learn how to detect quantum critical points in a quantum spin model. The model we are going to look at is the square lattice quantum Heisenberg model with dimerization in the form of ladder arrangements. Ladders with coupling $J_0$ on the legs and $J_1$ on the rungs are coupled together with coupling strength $J_2$. The model is depicted in Fig. 1 of Wenzel and Janke, Phys. Rev. B 79, 014410 (2009), albeit with different notations. In this tutorial we will consider the case $J_0=J_1=1$, and vary the interladder coupling $J_2$. Even though there are no phase transitions at finite temperatures in 2D Heisenberg models (Mermin-Wagner Theorem), transitions between different ground-states can occur at $T=0$.
(missing picture)
Identify the different phases
First of all, we consider the two simple limits of decoupled ladders ($J_2=0$) and of the isotropic square lattice ($J_2=1$). The decoupled ladders have a ground-state with short-range correlations and exhibit a finite spin gap : this is a spin liquid phase. On the other hand, the square lattice displays long-range order with a finite staggered magnetization : this is an antiferromagnetic Néel phase.
A simple and illustrative way of probing these two different physics is by looking at the magnetic susceptibility $\chi$. Let us simulate an 8x8 system using the following set of temperatures in the two different cases. Plot and compare the magnetic susceptibility in both the decoupled ($J_2=0$) and isotropic ($J_2=1$) situations. For decoupled ladders, the susceptibility exhibits an activated behaviour at low temperature due to the presence of the spin gap, whereas on the square lattice the susceptibility tends to a constant at low T. Please note that on a finite system, $\chi$ will always eventually tend to zero at small enough temperature due to the presence of a finite-size gap - this is however not our topic of interest here. You can run the simulation on the command line using the parameter file parm8a:
parameter2xml parm8a
loop parm8a.in.xml
or by using the python script tutorial8a.py
import pyalps
import matplotlib.pyplot as plt
import pyalps.pyplot
import numpy as np
import pyalps.fit_wrapper as fw
from math import sqrt
prepare the input parameters
parms = []
for j2 in [0.,1.]:
for t in [0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0]:
parms.append(
{
'LATTICE' : "coupled ladders",
'LATTICE_LIBRARY': 'lattices.xml',
'MODEL_LIBRARY' : 'models.xml',
'local_S' : 0.5,
'ALGORITHM' : 'loop',
'SEED' : 0,
'T' : t,
'J0' : 1 ,
'J1' : 1,
'J2' : j2,
'THERMALIZATION' : 5000,
'SWEEPS' : 50000,
'MODEL' : "spin",
'L' : 8,
'W' : 4
}
)
write the input file and run the simulation
input_file = pyalps.writeInputFiles('parm8a',parms)
pyalps.runApplication('loop',input_file)
For $J_2=0$, the value of the spin gap can be even further estimated by the finite-T behaviour of the magnetic susceptibility, using the following formula (derived in Phys. Rev. B 50, 13515 (1994)) $\chi=A/\sqrt{T}\exp(-\Delta/T)$ where $A$ and the spin gap $\Delta$ are fitting parameters. Fit the data for $T\leq1$ and extract an estimate for the spin gap. How does it compare to estimates available in the literature (such as in Phys. Rev. Lett. 73, 886 (1994) or Phys. Rev. Lett. 77, 1865 (1996)) ? Here is an example how you can load perform this analysis in python:
lines = []
for data in susc1:
pars = [fw.Parameter(1), fw.Parameter(1)]
data.y= data.y[data.x < 1]
data.x= data.x[data.x < 1]
f = lambda self, x, pars: (pars[0]()/np.sqrt(x))*np.exp(-pars[1]()/x)
fw.fit(None, f, pars, [v.mean for v in data.y], data.x)
prefactor = pars[0].get()
gap = pars[1].get()
print prefactor,gap
lines += plt.plot(data.x, f(None, data.x, pars))
lines[-1].set_label('$J_2=%.4s$: $\chi = \\frac{%.4s}{T}\exp(\\frac{-%.4s}{T})$' % (data.props['J2'], prefactor,gap))
Locate the phase transition
Having identified two different phases at $J_2=0$ and $J_2=1$, there must be (at least) one quantum phase transition separating them. We scan the coupling range $J_2 \in [0.2,0.4]$ for system sizes $L=8,10,12,16$ and simulate the model at an inverse temperate $\beta=2L$ using the parameter-file parm8b or the python script tutorial8b.py:
import pyalps
import matplotlib.pyplot as plt
import pyalps.pyplot
import numpy as np
prepare the input parameters
parms = []
for l in [8,10,12,16]:
for j2 in [0.2,0.25,0.3,0.35,0.4]:
parms.append(
{
'LATTICE' : "coupled ladders",
'LATTICE_LIBRARY': 'lattices.xml',
'MODEL_LIBRARY' : 'models.xml',
'local_S' : 0.5,
'ALGORITHM' : 'loop',
'SEED' : 0,
'BETA' : 2*l,
'J0' : 1 ,
'J1' : 1,
'J2' : j2,
'THERMALIZATION' : 5000,
'SWEEPS' : 50000,
'MODEL' : "spin",
'L' : l,
'W' : l/2
}
)
write the input file and run the simulation
input_file = pyalps.writeInputFiles('parm8a',parms)
pyalps.runApplication('loop',input_file)
Staggered magnetization, Binder Cumulant, spin stiffness
As in the classical Monte Carlo Tutorial we pinpoint the phase transition by an analysis of the Binder cumulant $U_4=\langle m_s^4\rangle /\langle m_s^2\rangle^2$ of the staggered magnetization $m_s$, which is the order parameter of the antiferromagnetic phase. Since the crossing point of the Binder cumulant shows large deviations at small system sizes for the model studied in this tutorial, we will also consider the spin stiffness, which has smaller finite-size corrections (Wenzel and Janke, Phys. Rev. B 79, 014410 (2009)). This observable is given by $\rho_s = \frac{3}{4\beta} \langle w_x^2 + w_y^2\rangle$ with winding numbers $w_x$,$w_y$ of worldlines along the $x$- and $y$-direction and it scales as $\rho_s \propto L^{d-2-z}$ at the quantum critical point, where $d$ is the dimension of the system and $z$ is the dynamical critical exponent. With $z=1$ the quantity $\rho_sL$ crosses at the critical point for different system system sizes $L$. Note that the fact that Binder cumulant and spin stiffness cross actually indicate that the phase transition is continuous, and not first order.
You can load and plot the observables by using the following lines:
data = pyalps.loadMeasurements(pyalps.getResultFiles(pattern='parm8a.task*.out.h5'),['Binder Ratio of Staggered Magnetization','Stiffness'])
binder=pyalps.collectXY(data,x='J2',y='Binder Ratio of Staggered Magnetization', foreach=['L'])
stiffness =pyalps.collectXY(data,x='J2',y='Stiffness', foreach=['L'])
for q in stiffness:
q.y = q.y*q.props['L']
make plot
plt.figure()
pyalps.pyplot.plot(stiffness)
plt.xlabel(r'$J2$')
plt.ylabel(r'Stiffness $\rho_s L$')
plt.title('coupled ladders')
plt.figure()
pyalps.pyplot.plot(binder)
plt.xlabel(r'$J_2$')
plt.ylabel(r'$g(m_s)$')
plt.title('coupled ladders')
plt.show()
What is your estimate of the quantum critical point?
Finite temperature effects
You may have noticed that the simulations are not performed at zero temperature, but at a finite value of inverse temperature $\beta=T^{-1}=2L$. One needs to be sure that the physics (and in particular the estimate of the quantum critical point) is not affected by finite-temperature effects. A brute-force yet simple way of checking this is to decrease the temperature and perform the same simulations: if results are not affected, then the simulations are indeed converged. Change $\beta=2L$ to $\beta=4L$, and check whether the stiffness and Binder cumulants are affected. Try now $\beta=L/4$: are results different?
Two remarks are important here: First, you may have noticed that the computational time scales approximatively linearly with $\beta$. In fact, given the path integral representation used, this means that the scaling of the loop algorithm used here is optimal (indeed, a finite-T QMC algorithm cannot scale faster than with the space-time volume $\beta L^d$), even at a quantum phase transition.
Second, why did we choose inverse temperature $\beta$ to be proportionnal to $L$? In fact, this linear relation between time and space scales originates from the value $z=1$ of the dynamical critical exponent for this quantum phase transition. You have implicitely checked it by looking at the crossings of $\rho_s L$. In general, $z$ does not need to be equal to unity, and the correct scaling of temperature with system size (in order to ensure ground-state sampling) has to be checked.
Estimates of critical exponents
You have obtained a rough estimate of the quantum critical point $J_2^c$. As in the classical case, extracting the critical exponents require more work and in particular a more precise determination of $J_2^c$.
We will obtain one by considering larger system sizes on a finer grid of $J_2^c$, the parameters are speficied in parm8d. Please note that these simulations will take quite some CPU time and we therefore we leave it to you as an exercise. Plot again the Binder cumulant of the staggered magnetization $U_4$ as well as the stiffness multiplied by system size $\rho_s L$ for different system sizes. The crossings of these curves should allow a more precise estimate of $J_2^c$. To obtain the critical exponent $\nu$ related to the divergence of the correlation length, it is useful to consider the scaling with system size of the derivative (with respect to $J_2^c$) of these quantities, when taken precisely at $J_2^c$. These derivatives $\frac{dU_4}{d J_2}$ and $L \frac{d\rho_s}{d J_2}$ can be obtained in principle as a Monte Carlo measurement, however for this tutorial, it is sufficient to perform a numerical differentiation which is possible thanks to the fine grid in $J_2$.
Perform the numerical differentiations for the different system sizes for both quantities, and plot their values at $J_2^c$ as a function of system size. Data should scale as a power law : $\frac{dU_4}{d J_2}(J_2^c) \propto L \frac{d\rho_s}{d J_2}(J_2^c) \propto L^{1/\nu}$. Which value of $\nu$ do you obtain?
** Exercise:** You can visually determine the quality of your estimate by performing data collapses, as in the classical case. The scaling forms for $U_4$ and $\rho_s L$ are identical to the ones for the Binder cumulant in the classical phase transition example.
Besides $z$ and $\nu$, the last independent exponent $\eta$ can be obtained in the following way. As in the classical case, the susceptibility at the critical point should scales as $\chi_s (J_2^c) \sim L^{2-\eta}$. Note that in this last expression one needs to consider the susceptibility associated to fluctuations of the order parameter: therefore the scaling with system size of the staggered susceptibility $\chi_s$, and not of the uniform susceptibility $\chi$, has to be performed here. Plot $\chi_s$ at $J_2^c$ as a function of system size, which estimate of $\eta$ do you obtain?
The quantum phase transition studied in this tutorial belongs to the universality class of the phase transition at finite temperature of the 3d classical Heisenberg model. Compare your estimates of critical exponents $\nu$ and $\eta$ with the ones reported in Phys. Rev. B 65, 144520 (2002) for this universality class.
Locate the second quantum critical point
This model not only undergoes one quantum phase transition - at higher values of $J_2$ you will find another quantum phase transition and we will repeat the same analysis for this parameter regime. Replace the line
for j2 in [0.2,0.25,0.3,0.35,0.4]:
by
for j2 in [1.8,1.85,1.9,1.95,2.,2.05,2.1]:
and run the simulation and the plotting commands again, see tutorial8c.py. Compare your results to the high-precision results presented in Wenzel and Janke, Phys. Rev. B 79, 014410 (2009).
Running the simulations in Vistrails
As usual you will find the Vistrails workflow for this tutorial in the same directory as the parameter files and the python scripts or you can download it here mc-08-quantum-phase-transition.vt