ED-05 ED Phase Transition

Critical point of the Heisenberg chain with next-nearest-neighbour interaction

In this tutorial, we will follow up on the last part of the ED-04 tutorial, where the Heisenberg chain was considered. We will add a next-nearest-neighbour coupling term to the Hamiltonian,

where the first sum runs over nearest-neighbour bonds and the second over next-nearest-neighbour bonds on the same periodic chain. This is the celebrated $J_1$–$J_2$ (or Majumdar-Ghosh) chain.

In the limit of $J_2 = 0$, this model reduces to the critical Heisenberg chain, which is solvable by Bethe ansatz. At $J_2/J_1=0.5$, the model is also exactly solvable: it sits at the Majumdar-Ghosh point, where the ground state is an exact product of nearest-neighbour singlets (C.K. Majumdar and D.K. Ghosh, J. Math. Phys. 10, 1388 (1969)),

This dimerized, doubly-degenerate ground state is qualitatively different from the gapless, translationally-invariant ground state at $J_2=0$, which is of course indication of a phase transition at some intermediate $J_2/J_1 \in (0,1/2)$ — numerically known to occur at $J_2/J_1 \approx 0.2411$.

In the first part of this tutorial, we will locate the position of the critical point by looking at how the spectrum, in particular the gap in different symmetry sectors, changes as we tune the couplings. In the second part, we will revisit the CFT content of the critical chain. Analytically, it can be shown that the model at criticality is described by the same CFT as the Heisenberg chain, but the weight of the marginal operator which lead to the logarithmically vanishing finite-size corrections is zero and therefore the scaling dimensions can be found much more accurately.

Parameters

Lattice

The nnn chain lattice extends the periodic chain with a second set of bonds connecting each site to its next-nearest neighbour:

        J          J          J
   o---------o---------o---------o---(periodic)
   0         1         2         3
    \___________________________/
        J1 (=J2)   next-nearest-neighbour bonds

This built-in lattice is documented, along with the rest of the ALPS lattice library, as an extension of the plain chain lattice used in the previous tutorials.

Method

As in ED-02, we diagonalize the $S_z=0$ and $S_z=1$ sectors separately with sparsediag to get the singlet and triplet gaps; the largest sector here ($L=12$, $S_z=0$, dimension 924) is trivial for the Lanczos algorithm. We sweep $J_2$ densely enough to see the gaps cross, which locates the finite-size estimate of the critical point, and then diagonalize once more, at fixed $J_2/J_1=0.25$ close to that crossing, requesting more eigenstates to resolve the CFT tower.

So first, let us plot the energy of the ground state and the first excited state as well as the gap in the singlet ($S_z = 0$) and triplet ($S_z=1$) sector.

Using the command line

parm5a sweeps $J_2/J_1\in[0,0.5]$ for both $L=6$ and $L=8$, in both the $S_z=0$ and $S_z=1$ sectors:

MODEL="spin"
LATTICE="nnn chain lattice"
CONSERVED_QUANTUMNUMBERS="Sz"
local_S=1/2
J=1
NUMBER_EIGENVALUES=2
Sz_total=0
{ L=6; J1=0.0 }
{ L=6; J1=0.1 }
{ L=6; J1=0.2 }
{ L=6; J1=0.3 }
{ L=6; J1=0.4 }
{ L=6; J1=0.5 }
{ L=8; J1=0.0 }
{ L=8; J1=0.1 }
{ L=8; J1=0.2 }
{ L=8; J1=0.3 }
{ L=8; J1=0.4 }
{ L=8; J1=0.5 }
Sz_total=1
{ L=6; J1=0.0 }
{ L=6; J1=0.1 }
{ L=6; J1=0.2 }
{ L=6; J1=0.3 }
{ L=6; J1=0.4 }
{ L=6; J1=0.5 }
{ L=8; J1=0.0 }
{ L=8; J1=0.1 }
{ L=8; J1=0.2 }
{ L=8; J1=0.3 }
{ L=8; J1=0.4 }
{ L=8; J1=0.5 }

parameter2xml parm5a
sparsediag --write-xml parm5a.in.xml

Python

To automate the sweep and the data analysis, we use the script tutorial5a.py. We start with the usual imports:

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

Again, we use the $S_z$ quantum number, but now we will run simulations in different sectors $(S_z=0,1)$. We run for the system sizes $L=6,8$, since the effect we’re looking for occurs already at very small system sizs.

prefix = 'alps-nnn-heisenberg'
parms = []
for L in [6,8]:
    for Szt in [0,1]:
        for J1 in np.linspace(0,0.5,6):
            parms.append({
                'LATTICE'              : "nnn chain lattice",
                'MODEL'                : "spin",
                'local_S'              : 0.5,
                'J'                    : 1,
                'NUMBER_EIGENVALUES'   : 2,
                'CONSERVED_QUANTUMNUMBERS' : 'Sz',
                'Sz_total'             : Szt,
                'J1'                   : J1,
                'L'                    : L
            })

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

The data analysis is slightly more involved in this case than in the previous ones. In particular, we will rely heavily on the feature of hierarchical datasets. To understand the physics, it is actually sufficient to look only at the ground and first excited states - so if you feel that the calculation of the gaps is too confusing, don’t worry about it too much.

In a first step, we join all energies for a given set of J1, L, Sz_total and sort them. First, we group by the parameters - J1, L, Sz_total. Each item in the loop over grouped will therefore contain a list of datasets for different momenta. We will join those and then use the dict_intersect function to find the properties for the resulting dataset; this function just takes a list of dictionaries and returns the part that is equal for all of them. We use numpy’s argsort function to obtain the list of indices that will sort y; this allows us to sort x accordingly, although that will probably not be needed.

grouped = pyalps.groupSets(pyalps.flatten(data), ['J1', 'L', 'Sz_total'])
nd = []
for group in grouped:
    ally = []
    allx = []
    for q in group:
        ally += list(q.y)
        allx += list(q.x)
    r = pyalps.DataSet()
    sel = np.argsort(ally)
    r.y = np.array(ally)[sel]
    r.x = np.array(allx)[sel]
    r.props = dict_intersect([q.props for q in group])
    nd.append( r )
data = nd

Next, we have to remove states that occur in the $S_z=1$ sector from the $S_z=0$ sector. We group them by J1, L, such that each group contains the spectra for the two different Sz_total sectors. We then use the function subtract_spectrum, which removes elements from the dataset passed as first argument which are also contained in the second argument. As an optional argument, this function accepts a maximum relative difference.

grouped = pyalps.groupSets(pyalps.flatten(data), ['J1', 'L'])
nd = []
for group in grouped:
    if group[0].props['Sz_total'] == 0:
        s0 = group[0]
        s1 = group[1]
    else:
        s0 = group[1]
        s1 = group[0]
    s0 = pyalps.subtract_spectrum(s0, s1)
    nd.append(s0)
    nd.append(s1)
data = nd

Now, we create a new list of datasets (sector_E) that will contain only the energy of the ground state (‘gs’) or first excited state (‘fe’). We will store this onto the property which. That will subsequently allow us to use the collectXY function to create a plot of the gs and fe energy vs the coupling for each L.

sector_E = []
grouped = pyalps.groupSets(pyalps.flatten(data), ['Sz_total', 'J1', 'L'])
for group in grouped:
    allE = []
    for q in group:
        allE += list(q.y)
    allE = np.sort(allE)
    
    d = pyalps.DataSet()
    d.props = dict_intersect([q.props for q in group])
    d.x = np.array([0])
    d.y = np.array([allE[0]])
    d.props['which'] = 'gs'
    sector_E.append(d)
    
    d2 = copy.deepcopy(d)
    d2.y = np.array([allE[1]])
    d2.props['which'] = 'fe'
    sector_E.append(d2)

sector_energies = pyalps.collectXY(sector_E, 'J1', 'Energy', ['Sz_total', 'which', 'L'])
plt.figure()
pyalps.plot.plot(sector_energies)
plt.xlabel('$J_1/J$')
plt.ylabel('$E_0$')
plt.legend(prop={'size':8})

In the last step, we calculate the singlet and triplet gap.These are defined as the energy difference between the lowest state of the system and a) the first excited state in the singlet ($S_z=0$ sector, b) the lowest state in the triplet ($S_z=1$) sector.

grouped = pyalps.groupSets( pyalps.groupSets(pyalps.flatten(data), ['J1', 'L']), ['Sz_total'])

gaps = []
for J1g in grouped:
    totalmin = 1000
    for q in pyalps.flatten(J1g):
        totalmin = min(totalmin, np.min(q.y))
    
    for Szg in J1g:
        allE = []
        for q in Szg:
            allE += list(q.y)
        allE = np.sort(allE)
        d = pyalps.DataSet()
        d.props = pyalps.dict_intersect([q.props for q in Szg])
        d.props['observable'] = 'gap'
        print totalmin,d.props['Sz_total']
        if d.props['Sz_total'] == 0:
            d.y = np.array([allE[1]-totalmin])
        else:
            d.y = np.array([allE[0]-totalmin])
        d.x = np.array([0])
        d.props['line'] = '.-'
        gaps.append(d)

gaps = pyalps.collectXY(gaps, 'J1', 'gap', ['Sz_total', 'L'])

plt.figure()
pyalps.plot.plot(gaps)
plt.xlabel('$J_1/J$')
plt.ylabel('$\Delta$')
plt.legend(prop={'size':8})

plt.show()

Output data

Diagonalizing the sectors independently gives the following singlet and triplet gaps as a function of $J_2/J_1$ (here labelled J1 in the ALPS parameter, as noted above):

For both system sizes the singlet gap (initially the larger of the two, since $J_2=0$ is gapless with the lowest excitation approaching zero only as $L\to\infty$) drops below the triplet gap between $J_2/J_1=0.2$ and $0.3$, crossing almost exactly at $J_2/J_1\approx0.25$ — already close to the accepted thermodynamic-limit value $J_2/J_1\approx0.2411$ even at these small sizes. This level crossing is the finite-size signature of the transition into the dimerized phase; note also that the singlet gap vanishes exactly at $J_2/J_1=0.5$, the Majumdar-Ghosh point, consistent with its exactly two-fold degenerate ground state.

Heisenberg chain with next-nearest neighbour coupling: CFT assignments

The finite-size corrections can be significantly reduced by tuning to a critical point in a frustrated J1-J2 chain. Despite the different couplings, this model can be shown to have the same continuum critical field theory and we can therefore extract the scaling dimensions in this limit. For a detailed discussion of this point, take a look at the reference I. Affleck, D. Gepner, H.J. Schulz and T. Ziman, J. Phys. A: Math. Gen. 22, 511 (1989).

Compare to the spectrum you have obtained above: you will see that the correspondence to the expected scaling dimensions is much easier to see and that they converge much faster for increasing system size.

Using the command line

parm5b fixes $J_2/J_1=0.25$, close to the crossing located above, and requests 5 eigenstates in the $S_z=0$ sector for two system sizes:

MODEL="spin"
LATTICE="nnn chain lattice"
CONSERVED_QUANTUMNUMBERS="Sz"
local_S=1/2
J=1
J1=0.25
NUMBER_EIGENVALUES=5
Sz_total=0
{ L=10 }
{ L=12 }

parameter2xml parm5b
sparsediag --write-xml parm5b.in.xml

Python

The new parameters, used in the script tutorial5b.py, are:

parms_ = {
    'LATTICE'              : "nnn chain lattice",
    'MODEL'                : "spin",
    'local_S'              : 0.5,
    'J'                    : 1,
    'J1'                   : 0.25,
    'NUMBER_EIGENVALUES'   : 5,
    'CONSERVED_QUANTUMNUMBERS' : 'Sz',
    'Sz_total' : 0
}
prefix = 'nnn-heisenberg'
parms = []
for L in [10,12]:
    parms_.update({'L':L})
    parms.append(copy.deepcopy(parms_))

The rest of the script proceeds exactly as in the Heisenberg-chain CFT analysis of ED-04: write and run the input files, load the eigenstate measurements, rescale by the gap between the two lowest states in the $S_z=0$ sector, and overlay the expected scaling dimensions:

input_file = pyalps.writeInputFiles(prefix,parms)
res = pyalps.runApplication('sparsediag', input_file)
data = pyalps.loadEigenstateMeasurements(pyalps.getResultFiles(prefix=prefix))

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]

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

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

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

Rescaling the independently-diagonalized $S_z=0$-sector spectrum at $J_2/J_1=0.25$ the same way as above gives:

Compare this to the plain Heisenberg chain ($J_2=0$) values from ED-04, reproduced here for the same sizes: 0.880 ($L=10$) and 0.857 ($L=12$) for the level that should converge to 1. At $J_2/J_1=0.25$ the same level is already at 0.947 and 0.960 — both closer to the target and moving in the right direction as $L$ increases, exactly as expected once the marginal operator responsible for the log-corrections in ED-04 has been tuned away.

Summary

Comparing the singlet and triplet gaps of the frustrated $J_1$–$J_2$ chain locates its dimerization transition, by finite-size level crossing, within a few percent of the accepted value $J_2/J_1\approx0.2411$ using only $L=6,8$ chains; diagonalizing near that point additionally shows that the same $c=1$ CFT operator content identified in ED-04 emerges far more cleanly here, since tuning $J_2$ removes the marginal operator responsible for the slow logarithmic finite-size corrections of the unfrustrated chain.

Questions

• From the table of singlet and triplet gaps, estimate the crossing point $J_2/J_1$ separately for $L=6$ and $L=8$ by linear interpolation between the bracketing points. Does the estimate move towards $0.2411$ as $L$ increases?
• Why does the singlet gap, and not the triplet gap, vanish exactly at the Majumdar-Ghosh point $J_2/J_1=0.5$?
• Compare the $L=10,12$ rescaled level-3 values at $J_2/J_1=0.25$ to the corresponding ED-04 values at $J_2=0$. Which is closer to the CFT prediction of 1, and which converges faster with $L$?