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 the last part of the ED-04 tutorial, where the Heisenberg chain was considered. We will add a next-nearest neighbour coupling term $J_2 \sum_{\langle \langle i,j \rangle \rangle} S_i \cdot S_j$ to the Hamiltonian.

This is of course indication of a phase transition at some intermediate $J_1/J_2$ \in $(0,1/2)$. (insert references from http://pre.aps.org/pdf/PRE/v76/i6/e061108 once APS website works again)

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.

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.

Python

We start with the usual imports:

import pyalps
import pyalps.plot
from pyalps.dict_intersect import dict_intersect
import numpy as np
import matplotlib.plot 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_QUANTUMNUMBER' : '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 paramters - 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 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()

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

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.

Python

The new parameters are:

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

The remainder of the script is analogous to the ED-04 cases and can be found in [3].

ED-04 Criticality ED-06 Full Diagonalization