MC-08 Quantum Phase Transition

The goal of this tutorial is to locate and characterise quantum phase transitions in a two-dimensional spin model using finite-size scaling of quantum Monte Carlo data. Unlike classical phase transitions, quantum phase transitions occur at $T=0$ and are driven by quantum fluctuations rather than thermal ones. We will use the Binder cumulant and the spin stiffness to pin down the critical coupling, then extract the correlation-length exponent $\nu$ , the anomalous dimension $\eta$ , and the dynamical exponent $z$ — placing the transition in its universality class. The model turns out to have two quantum critical points, and we study both.

The model consists of spin- $\frac{1}{2}$ Heisenberg ladders arranged on a square lattice. Each ladder runs horizontally: sites within a ladder are connected along the legs by coupling $J_0$ and across the rungs by coupling $J_1$ . Adjacent ladders are coupled vertically by $J_2$ :

o-J0-o-J0-o-J0-o   <- leg
|         |
J1        J1        <- rung (within ladder)
|         |
o-J0-o-J0-o-J0-o   <- leg
|         |
J2        J2        <- inter-ladder coupling
|         |
o-J0-o-J0-o-J0-o   <- leg
|         |
J1        J1        <- rung (within ladder)
|         |
o-J0-o-J0-o-J0-o   <- leg

The parameter L sets the number of sites along each leg; W = L/2 sets the number of rows, giving W/2 ladders in total. In this tutorial we consider $J_0=J_1=1$ and vary the inter-ladder coupling $J_2$ (see also Fig. 1 of Wenzel and Janke, Phys. Rev. B 79, 014410 (2009) for an illustration). 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$ .

Background simulations

The finite-size scaling run across multiple system sizes takes significantly longer than the single-system scans below. Start it now so it runs in the background while you work through the rest of the tutorial.

Command line

Download parm8b and run:

parameter2xml parm8b
loop parm8b.in.xml &

Python

Run the first part of tutorial8b.py (the setup and pyalps.runApplication call) in a separate terminal or as a background process before continuing.

Identify the different phases

We begin by considering the two simple limits: decoupled ladders ( $J_2=0$ ) and 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. The square lattice, by contrast, displays long-range order with a finite staggered magnetization: this is an antiferromagnetic Néel phase.

A clear way to probe these two phases is through the magnetic susceptibility $\chi$ . Simulate an $8\times 8$ system over a range of temperatures for both cases and compare the results. For decoupled ladders the susceptibility exhibits activated behaviour at low temperature due to the spin gap; on the square lattice it tends to a finite constant as $T\to 0$ . Note that on any finite system $\chi$ will eventually tend to zero at low enough temperature due to the finite-size gap, but this is not our focus here. Run the simulation on the command line with parameter file parm8a:

parameter2xml parm8a
loop parm8a.in.xml

or with the Python script tutorial8a.py.

import pyalps
import matplotlib.pyplot as plt
import pyalps.plot
import numpy as np
import pyalps.fit_wrapper as fw

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',  # path to custom lattice definitions
              'MODEL_LIBRARY'  : 'models.xml',    # path to custom model definitions
              '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
            }
    )
    
input_file = pyalps.writeInputFiles('parm8a',parms)
pyalps.runApplication('loop',input_file)

For $J_2=0$ , the value of the spin gap can be estimated from the finite- $T$ behaviour of the magnetic susceptibility (derived in Phys. Rev. B 50, 13515 (1994)):

\chi = \frac{A}{\sqrt{T}} \exp\!\left(-\frac{\Delta}{T}\right),

where $A$ and the spin gap $\Delta$ are fitting parameters. Fit the data for $T\leq1$ and extract an estimate for the spin gap. Here is an example of how to perform this analysis in Python:

data = pyalps.loadMeasurements(pyalps.getResultFiles(pattern='parm8a.task*.out.h5'),['Staggered Susceptibility','Susceptibility'])
susc1=pyalps.collectXY(data,x='T',y='Susceptibility', foreach=['J2'])

lines = []
gap_j0 = None
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, np.array([v.mean for v in data.y]), data.x)
    prefactor = pars[0].get()
    gap = pars[1].get()
    if data.props['J2'] == 0.0:
        gap_j0 = gap
    
    lines += plt.plot(data.x, f(None, data.x, pars))
    lines[-1].set_label('$J_2=%.4s$: $\\chi = \\frac{%.4s}{\\sqrt{T}}\\exp\\left(\\frac{-%.4s}{T}\\right)$' % (data.props['J2'], prefactor,gap))

The fit returns the gap $\Delta$ and the prefactor $A$ ; the resulting curves are shown below.

plt.figure()
pyalps.plot.plot(susc1)
plt.xlabel(r'$T$')
plt.ylabel(r'$\chi$')
plt.title('spin gap $\\Delta \\approx$ %.4s' % gap_j0)
plt.legend()
plt.show()

For $J_2=0$ the susceptibility should drop steeply at low $T$ , and the fit curve should follow the data closely, yielding a spin gap $\Delta \approx 0.5$ . For $J_2=1$ the susceptibility should approach a finite constant as $T \to 0$ , reflecting the absence of a gap on the square lattice.

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. The background run started above scans the coupling range $J_2 \in [0.2,0.4]$ for system sizes $L=8,10,12,16$ at an inverse temperature $\beta=2L$ (see parm8b / tutorial8b.py for the full parameter set). Once it has finished, load and analyse the results as described below.

Choice of inverse temperature

The simulations are run at a finite inverse temperature $\beta = 2L$ rather than strictly at $T=0$ . The choice $\beta \propto L$ is not arbitrary: at this quantum phase transition the dynamical critical exponent is $z=1$ , which means that time and space scale the same way. Keeping $\beta = 2L$ ensures that the imaginary-time extent of the path integral grows with system size in the same proportion, so the simulation samples the ground-state physics rather than thermal excitations. You will verify this choice implicitly when you observe that the $\rho_s L$ curves cross at a single point — this crossing only works cleanly if $z=1$ and $\beta \propto L^z$ .

Note also that the computational time of the loop algorithm scales linearly with $\beta$ : since no finite- $T$ QMC algorithm can scale better than linearly with the space-time volume $\beta L^d$ , this scaling is optimal even at a quantum critical point.

To confirm that results are not affected by the finite temperature, change $\beta=2L$ to $\beta=4L$ and repeat the simulation — the stiffness and Binder cumulant should be unchanged. Try also $\beta=L/4$ to observe how results degrade when the system is far from the ground state.

Staggered magnetization, Binder cumulant, spin stiffness

As in the classical Monte Carlo tutorial, we pinpoint the phase transition using two observables.

The first is the Binder cumulant of the staggered magnetization $m_s$ , the order parameter of the antiferromagnetic phase:

U_4 = \frac{\langle m_s^4\rangle}{\langle m_s^2\rangle^2}.

The second is the spin stiffness (Wenzel and Janke, Phys. Rev. B 79, 014410 (2009)), which has smaller finite-size corrections than the Binder cumulant crossing for this model:

\rho_s = \frac{3}{4\beta} \langle w_x^2 + w_y^2\rangle,

where $w_x$ , $w_y$ are the winding numbers of worldlines along the $x$ - and $y$ -directions. At the quantum critical point $\rho_s \propto L^{d-2-z}$ , where $d$ is the spatial dimension and $z$ the dynamical critical exponent. With $z=1$ the combination $\rho_s L$ is dimensionless at criticality, so curves for different $L$ all cross at $J_2^c$ . The fact that both observables cross at a single point — rather than one diverging — indicates a continuous transition rather than a first-order one.

Load and plot the observables with:

data = pyalps.loadMeasurements(pyalps.getResultFiles(pattern='parm8b.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']

plt.figure()
pyalps.plot.plot(stiffness)
plt.xlabel(r'$J_2$')
plt.ylabel(r'Stiffness $\rho_s L$')
plt.title('coupled ladders')

plt.figure()
pyalps.plot.plot(binder)
plt.xlabel(r'$J_2$')
plt.ylabel(r'$g(m_s)$')
plt.title('coupled ladders')
plt.show()

The $\rho_s L$ curves for different $L$ should rise from zero on the left and cross near a single point, giving a first estimate of $J_2^c \approx 0.30$ – $0.31$ . The Binder cumulant curves cross in the same region, though with larger finite-size corrections that make the crossing point less sharp at these system sizes. Both observables crossing (rather than one diverging) confirms that the transition is continuous.

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 requires a more precise determination of $J_2^c$ .

This is done by running larger system sizes on a finer grid of $J_2$ values, as set up in parm8d and tutorial8d.py. Note that these simulations are CPU-intensive and are left as an exercise. Plot the Binder cumulant $U_4$ and the rescaled stiffness $\rho_s L$ for different system sizes; the crossing point gives a refined estimate of $J_2^c$ . To extract $\nu$ , consider how the derivatives of these quantities with respect to $J_2$ , evaluated at $J_2^c$ , scale with system size. These derivatives can in principle be measured directly in the Monte Carlo, but for this tutorial it is sufficient to compute them by numerical differentiation using the fine $J_2$ grid.

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:

\left.\frac{dU_4}{dJ_2}\right|_{J_2^c} \propto L\left.\frac{d\rho_s}{dJ_2}\right|_{J_2^c} \propto L^{1/\nu}.

Extract $\nu$ from a power-law fit. Both quantities should yield consistent estimates close to $\nu \approx 0.71$ , the 3D classical Heisenberg value.

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 from the finite-size scaling of the staggered susceptibility at the critical point:

\chi_s(J_2^c) \sim L^{2-\eta}.

Note that one must use the staggered susceptibility $\chi_s$ — the susceptibility associated with fluctuations of the order parameter — rather than the uniform susceptibility $\chi$ . Plot $\chi_s$ at $J_2^c$ as a function of system size and extract $\eta$ from the slope. The expected value is $\eta \approx 0.034$ , close to zero, meaning the staggered susceptibility grows nearly as $L^2$ at the critical point.

Both quantum phase transitions in this tutorial belong to the universality class of the finite-temperature transition of the 3D classical Heisenberg model.

Locate the second quantum critical point

When $J_2$ is very large, the inter-ladder coupling dominates over the intra-ladder couplings $J_0 = J_1 = 1$ , and the system rearranges into effective singlet dimers along the $J_2$ bonds. This is again a gapped spin liquid phase, so the phase diagram has the structure: spin liquid → Néel AFM → spin liquid as $J_2$ increases from 0 to large values. There must therefore be a second quantum critical point $J_2^{c_2}$ at which the Néel order is destroyed.

We repeat the finite-size scaling analysis in the parameter regime $J_2 \in [1.8, 2.1]$ using the parameter file parm8c or the script tutorial8c.py, which use the same system sizes and $\beta=2L$ as parm8b but scan the higher- $J_2$ range.

Command line

parameter2xml parm8c
loop parm8c.in.xml

Python

import pyalps
import matplotlib.pyplot as plt
import pyalps.plot
import numpy as np

parms = []
for l in [8,10,12,16]:
    for j2 in [1.8,1.85,1.9,1.95,2.,2.05,2.1]:
        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
            }
        )

input_file = pyalps.writeInputFiles('parm8c', parms)
pyalps.runApplication('loop', input_file)

Load and plot the stiffness $\rho_s L$ and the Binder cumulant using the same analysis commands as before, replacing parm8b with parm8c in the file pattern.

The $\rho_s L$ curves should again cross at a single point, now near $J_2^{c_2} \approx 1.91$ . The Binder cumulant crossing should appear in the same region. This second transition belongs to the same universality class as the first — the 3D classical Heisenberg model — so the same critical exponents $\nu$ and $\eta$ apply.

Questions

Plot the magnetic susceptibility for $J_2=0$ and $J_2=1$ . How do the low-temperature behaviours differ, and what does each tell you about the ground state?
For $J_2=0$ , fit $\chi$ to $A/\sqrt{T}\exp(-\Delta/T)$ and extract the spin gap $\Delta$ . How does your estimate compare with values from the literature (Phys. Rev. Lett. 73, 886 (1994); Phys. Rev. Lett. 77, 1865 (1996))?
From the crossings of the $\rho_s L$ and Binder cumulant curves, what is your estimate of the first quantum critical point $J_2^c$ ?
Change $\beta=2L$ to $\beta=4L$ and repeat. Are the stiffness and Binder cumulant affected? Now try $\beta=L/4$ : do the results change noticeably? What does this tell you about the role of temperature in these simulations?
From the power-law scaling $dU_4/dJ_2|_{J_2^c} \propto L^{1/\nu}$ , what value of $\nu$ do you obtain?
From the finite-size scaling of the staggered susceptibility $\chi_s(J_2^c) \sim L^{2-\eta}$ , what value of $\eta$ do you obtain?
Compare your estimates of $\nu$ and $\eta$ with the 3D classical Heisenberg universality class values reported in Phys. Rev. B 65, 144520 (2002).
(Advanced) Perform data collapses for $U_4$ and $\rho_s L$ using your estimates of $J_2^c$ and $\nu$ . How sensitive is the quality of the collapse to the precise value of $J_2^c$ ?
Scan $J_2 \in [1.8, 2.1]$ to locate the second quantum critical point. How does your estimate compare with the high-precision results of Wenzel and Janke, Phys. Rev. B 79, 014410 (2009)?

MC-07 Phase Transition