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.

Ising chain

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

$$ H=J_{z} \sum_{\langle i,j \rangle} S^i_z S^j_z + \Gamma \sum_i S^i_x $$

Here, the first sum runs over pairs of nearest neighbours. $\Gamma$ is referred to as transverse field; the system becomes critical for $\Gamma/J=\frac{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-90 (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 P. Pfeuty), the long-range correlations are found to decay as: $$ \langle S^i_z S^{i+n}_z \rangle \sim n^{-2\times 1/8} $$ $$ \langle S^i_y S^{i+n}_y \rangle \sim n^{-2\times(1+1/8)} $$ $$ \langle S^i_x S^{i+n}_x \rangle \sim n^{-2\times 1} $$ 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.

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.plot 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()

Heisenberg chain

Now let us consider a more complicated example: the antiferromagnetic Heisenberg chain for spin-1/2 degrees of freedom, described by $$ H = \sum_{\langle i,j \rangle} \mathbf{S}^i \cdot \mathbf{S}^j $$

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. The parameter file can be found here. 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.

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 et al 1989 J. Phys. A: Math. Gen. 22 511 and the next tutorial.