ED-04 Criticality

In this tutorial, we will look at critical spin chains and make a connection to their description in terms of conformal field theory. At a continuous quantum phase transition, a spin chain becomes scale-invariant and its long-distance physics is described by a $(1+1)$-dimensional conformal field theory (CFT). A finite chain can never be truly at the transition, but it inherits a universal fingerprint of it: the low-lying finite-size spectrum organises itself into towers labelled by the scaling dimensions $\Delta$ of the CFT’s primary operators and their descendants, with level spacings that vanish as $1/L$. Exact diagonalization is the ideal tool to see this directly, since it gives us the complete low-energy spectrum, including its momenta, for chains of a few tens of sites — enough to already resolve the universal operator content.

Ising chain

The first model we will consider is the critical Ising chain, given by the Hamiltonian

Here, the first sum runs over pairs of nearest neighbours. $\Gamma$ is referred to as transverse field; the system becomes critical for $\Gamma/J=\tfrac{1}{2}$. For $\Gamma=0$, the ground state is antiferromagnetic for $J\gt 0$ and ferromagnetic for $J \lt 0$. The system is exactly solvable (P. Pfeuty, Annals of Physics 57, 79 (1970)).

In the above equation, $\Delta$ refers to the scaling dimension of that field. The scaling fields occur in groups: the lowest one, referred to as primary field, comes with an infinite number of descendants with scaling dimension $\Delta + m$, $m \in \lbrace 1, 2, 3, ... \rbrace$.

In the exact solution of the Ising model (Eq. (3.7) in the paper by Pfeuty), the long-range correlations are found to decay as:

Additionally, we expect the scaling dimension of the identity operator to be 0.

We therefore expect scaling dimensions of 0, 1/8, 1, 1+1/8 to appear in the CFT of the Ising model. To see this, we will rescale all energies of the spectrum according to $E \rightarrow \frac{E-E_0}{(E_1-E_0)8}$. This will force the two lowest states to occur where we expect the scaling dimensions; we can then check whether the rest of the spectrum is consistent with this.

Parameters

Note that $S_z$ is not conserved once $\Gamma\neq 0$ (the transverse field flips spins), so unlike the previous tutorials there is no CONSERVED_QUANTUMNUMBERS line and sparsediag diagonalizes the full chain, only resolved by total momentum.

Lattice

The critical transverse-field Ising chain uses the same periodic chain lattice from the ALPS lattice library as the previous tutorials, now with an Ising bond $J_z$ and an on-site transverse field $\Gamma$:

    Jz    Jz    Jz          Jz
o------o------o------ ... ------o     (each o also feels an on-site field Gamma S^x)
0      1      2                L-1
|______________________________________|
                Jz   (bond L-1 -- 0, periodic)

Method

At $\Gamma=0.5$ the chain has no conserved $S_z$, so sparsediag’s Lanczos iteration is run on the full $2^L$-dimensional Hilbert space (only $1024$ and $4096$ states for $L=10,12$); requesting NUMBER_EIGENVALUES=5 returns the five lowest states needed to see the primary field and its first descendants. Because we only need a handful of low-lying states rather than the full spectrum, sparse Lanczos diagonalization remains far cheaper than full diagonalization even though no symmetry sector is used to reduce the problem size.

Using the command line

The parameter file parm_ising sets up the two system sizes:

LATTICE    = "chain lattice"
MODEL      = "spin"
local_S    = 0.5
Jxy        = 0
Jz         = -1
Gamma      = 0.5
NUMBER_EIGENVALUES = 5

{L=10}
{L=12}

parameter2xml parm_ising
sparsediag --write-xml parm_ising.in.xml

Python version

The parameter file for the Ising simulation can be found here. The simulation is run with the Python script.

We will first import some modules:

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

Then, let us set up the parameters for two system sizes. Be careful to use the transverse field $\Gamma$, not the longitudinal field $h$.

parms = []
for L in [10,12]:
    parms.append({
        'LATTICE'    : "chain lattice",
        'MODEL'      : "spin",
        'local_S'    : 0.5,
        'Jxy'        : 0,
        'Jz'         : -1,
        'Gamma'      : 0.5,
        'NUMBER_EIGENVALUES' : 5,
        'L'          : L
    })

As you can see, we will simulate two system sizes. Now let’s set up the input files and run the simulation:

prefix = 'ising'
input_file = pyalps.writeInputFiles(prefix,parms)
res = pyalps.runApplication('sparsediag', input_file)
# res = pyalps.runApplication('sparsediag', input_file, MPI=2, mpirun='mpirun')
data = pyalps.loadEigenstateMeasurements(pyalps.getResultFiles(prefix=prefix))

By uncommenting the second-last line and adapting the number of jobs and the name of your mpirun executable (which will default to mpirun), you could have ALPS use several CPUs simultaneously.

First, we will extract the lowest and first excited for each value of L and collect this into a dictionary:

E0 = {}
E1 = {}
for Lsets in data:
    L = pyalps.flatten(Lsets)[0].props['L']
    allE = []
    for q in pyalps.flatten(Lsets):
        allE += list(q.y)
    allE = np.sort(allE)
    E0[L] = allE[0]
    E1[L] = allE[1]

The above code works since we know that ALPS will load the energies in lists grouped by the simulation - data is therefore a list of lists where we have different simulations at the top level and different momenta below. Now we rescale the energies according to the equation given above and collect the data as a function of momenta:

for q in pyalps.flatten(data):
    L = q.props['L']
    q.y = (q.y-E0[L])/(E1[L]-E0[L]) * (1./8.)
spectrum = pyalps.collectXY(data, 'TOTAL_MOMENTUM', 'Energy', foreach=['L'])

For comparison, let us also show the primary fields and their first few descendants:

for SD in [0.125, 1, 1+0.125, 2]:
    d = pyalps.DataSet()
    d.x = np.array([0,4])
    d.y = SD+0*d.x
    spectrum += [d]

Finally we create the plot:

pyalps.plot.plot(spectrum)
plt.legend(prop={'size':8})
plt.xlabel("$k$")
plt.ylabel("E_0")
plt.xlim(-0.02, math.pi+0.02)
plt.show()

Output data

Independently diagonalizing the same Hamiltonian gives the following rescaled low-energy levels (momentum-summed, sorted) for the two chain lengths:

By construction the first two rescaled levels sit exactly at 0 and 0.125; the third and fourth levels are genuine predictions of the diagonalization and are already within a few percent of the expected CFT values 1 and 1.125 at these small sizes, converging further as $L$ increases — direct numerical evidence for the operator content $\{0, \tfrac18, 1, 1+\tfrac18, \ldots\}$ of the Ising CFT.

Heisenberg chain

Now let us consider a more complicated example: the antiferromagnetic Heisenberg chain for spin-1/2 degrees of freedom, described by

The critical theory for this model has a central charge $c=1$ with primary fields 0, 0.5 and 1. As opposed to the Ising model, finite-size corrections vanish only logarithmically and are therefore quite pronounced at the system sizes we reach.

Parameters

Lattice

Same periodic chain lattice as above (see ALPS lattice library), now with isotropic Heisenberg coupling and no transverse field:

    J     J     J           J
o-------o-------o--- ... ---o
0       1       2          L-1
|_________________________________|
                J   (bond L-1 -- 0, periodic)

Method

Unlike the Ising point, the Heisenberg chain conserves $S_z$, so we can restrict sparsediag to the $S_z=0$ sector — this already contains the ground state and enough of the low-lying tower to read off the CFT content, while reducing the sector dimension from $2^{12}=4096$ to $\binom{12}{6}=924$ for $L=12$.

Using the command line

The parameter file sets up the two system sizes:

LATTICE    = "chain lattice"
MODEL      = "spin"
local_S    = 0.5
Jxy        = 1
Jz         = 1
NUMBER_EIGENVALUES = 5
CONSERVED_QUANTUMNUMBERS = 'Sz'
Sz_total   = 0

{L=10}
{L=12}

parameter2xml parm_heisenberg
sparsediag --write-xml parm_heisenberg.in.xml

Python version

The Python file is analogous to the above description and will therefore not be described here in detail. The primary difference is that we can make use of the U(1) symmetry of the system; we also run the simulations only for the $S_z = 0$ sector, since this contains all the relevant states. The rescaling used in heisenberg.py fixes the first excited level at 0.5 (the spin-operator scaling dimension) rather than 0.125, since the leading primary field here is $\Delta_1=1/2$.

Looking at the spectrum, try to identify the different scaling fields and how they approach the correct values as the system size is increased. You will notice that this is quite difficult. For a detailed discussion of the system, refer to I. Affleck, D. Gepner, H.J. Schulz and T. Ziman, J. Phys. A: Math. Gen. 22, 511 (1989) and the next tutorial.

Output data

Diagonalizing independently, the $S_z=0$-sector levels rescaled by $(E-E_0)/\Delta_{S_z=1}\times\tfrac12$ (where $\Delta_{S_z=1}=E_0(S_z{=}1)-E_0(S_z{=}0)$ is the triplet gap) come out as:

The third level, expected to converge to the descendant at scaling dimension 1, instead drifts slowly away from it, from 0.880 ($L=10$) to 0.825 ($L=16$), rather than approaching it monotonically the way the Ising levels did. This sluggish, non-monotonic behaviour is the practical symptom of the marginally irrelevant operator mentioned above, whose log-suppressed corrections dominate the size dependence at accessible $L$ and only fade at astronomically larger sizes. ED-05 shows how tuning to a nearby frustrated point removes this marginal operator and restores fast, textbook $1/L$ convergence.

Summary

Both the exactly solvable transverse-field Ising chain and the harder-to-converge Heisenberg chain confirm, from finite-size exact-diagonalization spectra alone, the universal CFT operator content ($c=1/2$ with dimensions $\{0,\tfrac18,1,\ldots\}$ for Ising; $c=1$ with dimensions $\{0,\tfrac12,1,\ldots\}$ for Heisenberg) predicted by their respective field theories — with the Heisenberg case also illustrating how a marginally irrelevant operator can make that convergence numerically deceptive.

Questions

• For the Ising chain, why do the rescaled energies come in degenerate pairs starting from the third level (e.g. two states near 1.107)? (Hint: think about the momentum quantum number of the corresponding descendant operators.)
• Repeat the Ising CFT analysis for a value of $\Gamma$ away from the critical point (e.g. $\Gamma=0.2$ or $\Gamma=1.0$). Does the rescaled spectrum still look CFT-like as $L$ increases?
• For the Heisenberg chain, the third rescaled level moves away from its CFT target as $L$ grows from 10 to 16, instead of settling towards it. Compute a couple of further sizes (e.g. $L=18,20$) if you can — does the drift eventually turn around, and how would you distinguish “very slow monotonic convergence” from “genuinely non-monotonic” behaviour with only a handful of accessible system sizes?