Skip to content
Tutorials
Monte Carlo Simulations
MC-07 Phase Transition

MC-07 Phase Transition

The goal of this tutorial is to learn how to detect a second-order phase transition from finite-size simulations, using the classic example of the two-dimensional Ising model. No true phase transition can occur on a finite system, but finite-size simulations show clear precursor signs that — combined with finite-size scaling — allow a precise determination of the critical temperature and the universality class of the transition.

Almost everything is known about the phase transition in the 2D Ising model since it is exactly solvable. Here we recover the location of the critical point and the critical exponents as if they were unknown, to illustrate the methods. Extracting precise estimates requires long simulations, so we start the fine-grid run in the background at the outset and analyse the coarse-grid run first.

Background simulations

Start the fine-grid simulations now so they run while you work through the rest of the tutorial.

Command line

Download parm7b and run:

parameter2xml parm7b
spinmc --Tmin 10 parm7b.in.xml &

The --Tmin 10 flag sets a checkpoint interval of 10 seconds, allowing the simulation to be safely interrupted and resumed.

Python

The first part of tutorial7b.py sets up and launches the same simulation:

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

parms = []
for l in [32, 48, 64]:
    for t in [2.24, 2.25, 2.26, 2.27, 2.28, 2.29, 2.30, 2.31, 2.32, 2.33, 2.34, 2.35]:
        parms.append(
            {
                'LATTICE'        : "square lattice",
                'T'              : t,
                'J'              : 1,
                'THERMALIZATION' : 5000,
                'SWEEPS'         : 150000,
                'UPDATE'         : "cluster",
                'MODEL'          : "Ising",
                'L'              : l
            }
        )

input_file = pyalps.writeInputFiles('parm7b', parms)
pyalps.runApplication('spinmc', input_file, Tmin=5)

Rough scan: locating the phase transition

We first make a coarse temperature scan on small systems to locate the critical region.

Command line

Download parm7a and run:

parameter2xml parm7a
spinmc --Tmin 5 parm7a.in.xml

Python

Using tutorial7a.py:

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

parms = []
for l in [4, 8, 16]:
    for t in [5.0, 4.5, 4.0, 3.5, 3.0, 2.9, 2.8, 2.7]:
        parms.append(
            {
                'LATTICE'        : "square lattice",
                'T'              : t,
                'J'              : 1,
                'THERMALIZATION' : 1000,
                'SWEEPS'         : 400000,
                'UPDATE'         : "cluster",
                'MODEL'          : "Ising",
                'L'              : l
            }
        )
    for t in [2.6, 2.5, 2.4, 2.3, 2.2, 2.1, 2.0, 1.9, 1.8, 1.7, 1.6, 1.5, 1.2]:
        parms.append(
            {
                'LATTICE'        : "square lattice",
                'T'              : t,
                'J'              : 1,
                'THERMALIZATION' : 1000,
                'SWEEPS'         : 40000,
                'UPDATE'         : "cluster",
                'MODEL'          : "Ising",
                'L'              : l
            }
        )

input_file = pyalps.writeInputFiles('parm7a', parms)
pyalps.runApplication('spinmc', input_file, Tmin=5)  # Tmin: minimum wall time in seconds between checkpoints

Magnetization, susceptibility, and specific heat

The order parameter of the Ising transition is the magnetization per site mm. On a finite system ⟨m⟩=0\langle m \rangle = 0 by symmetry, so we plot the average absolute value ⟨∣m∣⟩\langle |m| \rangle instead. Load the results and plot them:

pyalps.evaluateSpinMC(pyalps.getResultFiles(prefix='parm7a'))  # computes derived quantities (susceptibility, specific heat, Binder cumulant) from raw MC output
data = pyalps.loadMeasurements(
    pyalps.getResultFiles(prefix='parm7a'),
    ['|Magnetization|', 'Connected Susceptibility', 'Specific Heat', 'Binder Cumulant', 'Binder Cumulant U2']
)
magnetization_abs = pyalps.collectXY(data, x='T', y='|Magnetization|',          foreach=['L'])
connected_susc    = pyalps.collectXY(data, x='T', y='Connected Susceptibility',  foreach=['L'])
spec_heat         = pyalps.collectXY(data, x='T', y='Specific Heat',             foreach=['L'])
binder_u4         = pyalps.collectXY(data, x='T', y='Binder Cumulant',           foreach=['L'])
binder_u2         = pyalps.collectXY(data, x='T', y='Binder Cumulant U2',        foreach=['L'])

plt.figure()
pyalps.plot.plot(magnetization_abs)
plt.xlabel('Temperature $T$')
plt.ylabel('Magnetization $|m|$')
plt.title('2D Ising model')

The magnetization rises from zero at high temperature to its saturation value at low TT. The upturn sharpens with increasing system size, but the precise transition temperature is hard to read off directly.

To locate the transition more precisely, examine fluctuations. The connected susceptibility χ=βN(⟨m2⟩−⟨∣m∣⟩2)\chi = \beta N (\langle m^2 \rangle - \langle |m| \rangle^2) measures magnetization fluctuations and diverges at the critical point:

plt.figure()
pyalps.plot.plot(connected_susc)
plt.xlabel('Temperature $T$')
plt.ylabel('Connected Susceptibility $\chi_c$')
plt.title('2D Ising model')

A marked peak appears around T=2.2T = 2.2–2.42.4, and the peak grows with system size. The specific heat cv=β2N(⟨e2⟩−⟨e⟩2)c_v = \beta^2 N (\langle e^2 \rangle - \langle e \rangle^2), where ee is the internal energy per site, shows a similar but less pronounced peak in the same range:

plt.figure()
pyalps.plot.plot(spec_heat)
plt.xlabel('Temperature $T$')
plt.ylabel('Specific Heat $c_v$')
plt.title('2D Ising model')
plt.show()

The specific heat peak near T≈2.3T \approx 2.3 grows only weakly with system size — much more slowly than the susceptibility — indicating that the exponent α\alpha is close to zero.

Binder cumulant

The Binder cumulant provides a sharper tool for locating the critical temperature. Define the ratio

U4=⟨m4⟩⟨m2⟩2,U_4 = \frac{\langle m^4 \rangle}{\langle m^2 \rangle^2},

which equals 1 in the fully ordered low-temperature phase and 3 for a Gaussian (high-temperature) order-parameter distribution. At the critical point its value is universal: it is the same for all system sizes, so curves for different LL all cross at TcT_c.

plt.figure()
pyalps.plot.plot(binder_u4)
plt.xlabel('Temperature $T$')
plt.ylabel('Binder Cumulant $U_4$')
plt.title('2D Ising model')
plt.show()

From the crossing of the curves we can read off Tc∈[2.2,2.3]T_c \in [2.2, 2.3]. The second cumulant U2=⟨m2⟩/⟨∣m∣⟩2U_2 = \langle m^2 \rangle / \langle |m| \rangle^2 has the same crossing property; plotting it from the already-loaded binder_u2 data is left as an exercise.

Precise determination of TcT_c and critical exponents

Using the finer temperature grid and larger system sizes from parm7b, we can extract TcT_c and the critical exponents more accurately.

Load the results:

pyalps.evaluateSpinMC(pyalps.getResultFiles(prefix='parm7b'))  # computes derived quantities from raw MC output
data = pyalps.loadMeasurements(
    pyalps.getResultFiles(prefix='parm7b'),
    ['|Magnetization|', 'Connected Susceptibility', 'Specific Heat', 'Binder Cumulant', 'Binder Cumulant U2']
)
magnetization_abs = pyalps.collectXY(data, x='T', y='|Magnetization|',          foreach=['L'])
connected_susc    = pyalps.collectXY(data, x='T', y='Connected Susceptibility',  foreach=['L'])
spec_heat         = pyalps.collectXY(data, x='T', y='Specific Heat',             foreach=['L'])
binder_u4         = pyalps.collectXY(data, x='T', y='Binder Cumulant',           foreach=['L'])
binder_u2         = pyalps.collectXY(data, x='T', y='Binder Cumulant U2',        foreach=['L'])

Binder cumulant crossing 

plt.figure()
pyalps.plot.plot(binder_u4)
plt.xlabel('Temperature $T$')
plt.ylabel('Binder Cumulant $U_4$')
plt.title('2D Ising model')
plt.show()

The crossing of the L=32L = 32, 4848, and 6464 curves gives a refined estimate of TcT_c.

Data collapse and the correlation-length exponent ν\nu

Finite-size scaling predicts that the Binder cumulant obeys the scaling form

U4=f ⁣(L1/ν T−TcTc),U_4 = f\!\left(L^{1/\nu}\,\frac{T - T_c}{T_c}\right),

where ff is a universal function. Shifting the horizontal axis to (T−Tc)/Tc(T - T_c)/T_c should make the curves cross near zero. Multiplying by LaL^a should then collapse all curves onto a single master curve when a=1/νa = 1/\nu. Try values of aa close to 1:

Tc = ...   # your estimate from the Binder crossing
a  = ...   # try values near 1.0

for d in binder_u4:
    d.x -= Tc
    d.x  = d.x / Tc
    l    = d.props['L']
    d.x  = d.x * pow(float(l), a)

plt.figure()
pyalps.plot.plot(binder_u4)
plt.xlabel('$L^{1/\\nu}(T-T_c)/T_c$')
plt.ylabel('Binder Cumulant $U_4$')
plt.title('2D Ising model')
plt.show()

When the collapse is good, the correlation-length exponent is ν=1/a\nu = 1/a. This determination of ν\nu is independent of all other critical exponents.

Peak scaling of susceptibility and specific heat

Finite-size scaling also predicts how the peak values of the susceptibility and specific heat grow with system size at TcT_c:

χ(Tc)∼Lγ/ν,Cv(Tc)∼Lα/ν.\chi(T_c) \sim L^{\gamma/\nu}, \qquad C_v(T_c) \sim L^{\alpha/\nu}.

Note that the peak positions drift slightly with LL as Tc(L)=Tc+A L−1/νT_c(L) = T_c + A\,L^{-1/\nu}, so you can evaluate the peaks either at the true TcT_c from the Binder crossing or at the numerically determined peak temperature for each LL.

plt.figure()
pyalps.plot.plot(connected_susc)
plt.xlabel('Temperature $T$')
plt.ylabel('Connected Susceptibility $\chi_c$')
plt.title('2D Ising model')

plt.figure()
pyalps.plot.plot(spec_heat)
plt.xlabel('Temperature $T$')
plt.ylabel('Specific Heat $c_v$')
plt.title('2D Ising model')
plt.show()

Susceptibility peak fit

Extract the peak values and fit a power law in LL:

cs_mean = []
for q in connected_susc:
    cs_mean.append(np.array([d.mean for d in q.y]))

peak_cs       = pyalps.DataSet()
peak_cs.props = pyalps.dict_intersect([q.props for q in connected_susc])
peak_cs.y     = np.array([np.max(q) for q in cs_mean])
peak_cs.x     = np.array([q.props['L'] for q in connected_susc])

sel       = np.argsort(peak_cs.x)
peak_cs.y = peak_cs.y[sel]
peak_cs.x = peak_cs.x[sel]

pars = [fw.Parameter(1), fw.Parameter(1)]
f    = lambda self, x, pars: pars[0]() * np.power(x, pars[1]())
fw.fit(None, f, pars, peak_cs.y, peak_cs.x)
gamma_nu = pars[1].get()

plt.figure()
plt.plot(peak_cs.x, f(None, peak_cs.x, pars))
pyalps.plot.plot(peak_cs)
plt.xlabel('System Size $L$')
plt.ylabel('Connected Susceptibility $\\chi_c(T_c)$')
plt.title('2D Ising model, $\\gamma/\\nu$ = %.4s' % gamma_nu)
plt.show()

The fitted exponent gives γ/ν\gamma/\nu. Using the hyperscaling relation γ/ν=2−η\gamma/\nu = 2 - \eta, you can read off the anomalous dimension η\eta.

Specific heat peak fit

Repeat the same analysis for the specific heat:

sh_mean = []
for q in spec_heat:
    sh_mean.append(np.array([d.mean for d in q.y]))

peak_sh       = pyalps.DataSet()
peak_sh.props = pyalps.dict_intersect([q.props for q in spec_heat])
peak_sh.y     = np.array([np.max(q) for q in sh_mean])
peak_sh.x     = np.array([q.props['L'] for q in spec_heat])

sel       = np.argsort(peak_sh.x)
peak_sh.y = peak_sh.y[sel]
peak_sh.x = peak_sh.x[sel]

pars = [fw.Parameter(1), fw.Parameter(1)]
f    = lambda self, x, pars: pars[0]() * np.power(x, pars[1]())
fw.fit(None, f, pars, peak_sh.y, peak_sh.x)
alpha_nu = pars[1].get()

plt.figure()
plt.plot(peak_sh.x, f(None, peak_sh.x, pars))
pyalps.plot.plot(peak_sh)
plt.xlabel('System Size $L$')
plt.ylabel('Specific Heat $c_v(T_c)$')
plt.title('2D Ising model, $\\alpha/\\nu$ = %.4s' % alpha_nu)
plt.show()

Note that α\alpha can be positive or negative: the specific heat need not diverge as a power law at a continuous transition.

Data collapse for susceptibility and magnetization

The full scaling forms are

χ=L2−η g ⁣(L1/ν T−TcTc),∣m∣=L−β/ν h ⁣(L1/ν T−TcTc),\chi = L^{2-\eta}\, g\!\left(L^{1/\nu}\,\frac{T-T_c}{T_c}\right), \qquad |m| = L^{-\beta/\nu}\, h\!\left(L^{1/\nu}\,\frac{T-T_c}{T_c}\right),

where gg and hh are universal functions. Using your estimates of ν\nu and TcT_c from the Binder collapse, tune 2−η2-\eta to collapse the susceptibility curves:

for d in connected_susc:
    d.x -= Tc
    d.x  = d.x / Tc
    l    = d.props['L']
    d.x  = d.x * pow(float(l), a)

two_minus_eta = ...   # your estimate
for d in connected_susc:
    l   = d.props['L']
    d.y = d.y / pow(float(l), two_minus_eta)

plt.figure()
pyalps.plot.plot(connected_susc)
plt.xlabel('$L^{1/\\nu}(T-T_c)/T_c$')
plt.ylabel('$L^{-(2-\\eta)}\\chi_c$')
plt.title('2D Ising model')
plt.show()

And tune β/ν\beta/\nu to collapse the magnetization curves:

for d in magnetization_abs:
    d.x -= Tc
    d.x  = d.x / Tc
    l    = d.props['L']
    d.x  = d.x * pow(float(l), a)

beta_over_nu = ...   # your estimate
for d in magnetization_abs:
    l   = d.props['L']
    d.y = d.y / pow(float(l), -beta_over_nu)  # dividing by L^(-β/ν) = multiplying by L^(β/ν)

plt.figure()
pyalps.plot.plot(magnetization_abs)
plt.xlabel('$L^{1/\\nu}(T-T_c)/T_c$')
plt.ylabel('$L^{\\beta/\\nu}|m|$')
plt.title('2D Ising model')
plt.show()

Exact values and comparison

The exact critical exponents of the 2D Ising model are

ν=1,η=14,β=18,α=0.\nu = 1, \quad \eta = \tfrac{1}{4}, \quad \beta = \tfrac{1}{8}, \quad \alpha = 0.

For α=0\alpha = 0: the specific heat diverges logarithmically rather than as a power law, so a power-law fit returns an exponent close to zero but not exactly zero.

For more precise estimates, fix the critical temperature to the exact value Tc=2/ln⁡(1+2)≈2.2692T_c = 2/\ln(1+\sqrt{2}) \approx 2.2692 and simulate larger system sizes. With TcT_c fixed, a useful way to extract ν\nu is via the derivative of the Binder cumulant. At TcT_c, finite-size scaling implies dU4/dT∼L1/νdU_4/dT \sim L^{1/\nu}, so a power-law fit of this derivative as a function of LL gives ν\nu directly. The derivative can be obtained by numerical differentiation of the Monte Carlo data, or measured as a thermodynamic average during the simulation (which requires good statistics).

Questions

  • At what temperature do the Binder cumulant curves cross? How does your estimate compare with the exact value Tc≈2.2692T_c \approx 2.2692?
  • How do your numerical estimates of ν\nu, η\eta, β\beta, and α\alpha compare with the exact values?
  • What value of ν\nu do you extract from the Binder cumulant data collapse? Compare with the exact value ν=1\nu = 1.
  • The susceptibility peak grows strongly with system size while the specific heat peak grows much more slowly. What does this imply about α\alpha?
  • Using γ/ν=2−η\gamma/\nu = 2 - \eta, what value of η\eta do you obtain? Compare with the exact value η=1/4\eta = 1/4.
  • Try increasing SWEEPS for the largest system sizes. How does this improve the quality of the data collapse?
  • (Advanced) Compute dU4/dTdU_4/dT numerically and verify the L1/νL^{1/\nu} scaling.
MC-06 QWLMC-08 Quantum Phase Transition