Direct Calculation of spectral functions using symmetries

This example demonstrates how to use the general basis function method Op_shift_sector() to compute spectral functions directly without Fourier -transforming auto-correlation functions. The idea here is to use SciPy’s iterative solvers for sparse linear systems to calculate the action of an inverted operator on a state.

The spectral function in ground state $| 0 ⟩$ can be written as:

C (ω) = - \frac{1}{π} ⟨ 0 | A^{†} \frac{1}{ω + i η + E_{0} - H} A | 0 ⟩,

where $η$ is the boradening factor for the spectral peaks. For an operator $A$ that has quantum numbers that differ from the ground state $| 0 ⟩$ we can use Op_shift_sector() to calculate the a new ket: $| A ⟩ = A | 0 ⟩$ . With this new ket, we can construct the Hamiltonian $H$ in the new sector, and solve the following equation:

(ω + i η + E_{0} - H) | x (ω) ⟩ = | A ⟩ .

This equation can be solved using the bi-conjugate gradient (cf. scipy.sparse.linalg.bicg) method that is implemented in SciPy. In order to solve that equation with SciPy we define a custom class (see code below) that produces the action $(ω \pm i η + E_{0} - H) | u ⟩$ for an arbitrary state $| u ⟩$ . Now we can use $| x (ω) ⟩$ to calculate the spectral function a given value of $ω$ as follows:

C (ω) = - \frac{1}{π} ⟨ A | x (ω) ⟩ .

Note that this method can also be used to calculate spectral functions for any state by replacing $| 0 ⟩ \to | ψ ⟩$ . This may be useful for calculating spectral functions away from equilibrium by using $| ψ (t) ⟩$ that has undergone some kind of time evolution (in the Schrödinger picture).

In the example below, we look at the Heisenberg chain ( $J = 1$ ) with periodic boundary conditions. The periodic boundary condition allows us to use translation symmetry to reduce the total Hilbert space into blocks labeled by the many-body momentum quantum number. We limit the chain lengths to $L = 4 n$ for $n = 1, 2, 3. . .$ , since this is required to have the ground state in the zero-momentum sector (for other system sizes, the ground state of the isotropic Heisenberg chain has finite momentum).

Below, we calculate two different spectral functions:

G_{z z} (ω, q) = ⟨ 0 | S_{- q}^{z} \frac{1}{ω + i η + E_{0} - H} S_{q}^{z} | 0 ⟩

G_{+ -} (ω, q) = ⟨ 0 | S_{- q}^{-} \frac{1}{ω + i η + E_{0} - H} S_{q}^{+} | 0 ⟩

where we have defined:

S_{q}^{z} = \frac{1}{\sqrt{L}} \sum_{r = 0}^{L - 1} \exp (- i \frac{2 π q r}{L}) S_{r}^{z}

S_{q}^{\pm} = \frac{1}{\sqrt{2 L}} \sum_{r = 0}^{L - 1} \exp (- i \frac{2 π q r}{L}) S_{r}^{\pm}

Because the model has full SU(2) symmetry, we expect that the two spectral functions should give the same result. We also exclude the spectral function for $q = L / 2$ ( $π$ -momentum) as the spectral peak is very large due to the quasi long-range antiferromagnetic order in the ground state. The variable on_the_fly can be used to switch between using a hamiltonian object to a quantum_LinearOperator object to calculate the matrix vector product. One can also change the total spin of the local spin operators by changing the variable S to compute the spectral function for higher-spin Heisenberg models.

Script

download script

#
import sys, os

os.environ["KMP_DUPLICATE_LIB_OK"] = (
    "True"  # uncomment this line if omp error occurs on OSX for python 3
)
os.environ["OMP_NUM_THREADS"] = "1"  # set number of OpenMP threads to run in parallel
os.environ["MKL_NUM_THREADS"] = "1"  # set number of MKL threads to run in parallel
#

#######################################################################
#                            example 26                               #
# This example shows how to use the `Op-shit_sector` method of the    #
# general basis class to compute spectral functions using symmetries. #
#######################################################################
from quspin.basis import spin_basis_general
from quspin.operators import hamiltonian, quantum_LinearOperator
import scipy.sparse as sp
import numexpr, cProfile
import numpy as np
import matplotlib.pyplot as plt


#
#
# define custom LinearOperator object that generates the left hand side of the equation.
#
class LHS(sp.linalg.LinearOperator):
    #
    def __init__(self, H, omega, eta, E0, kwargs={}):
        self._H = H  # Hamiltonian
        self._z = omega + 1j * eta + E0  # complex energy
        self._kwargs = kwargs  # arguments

    #
    @property
    def shape(self):
        return (self._H.Ns, self._H.Ns)

    #
    @property
    def dtype(self):
        return np.dtype(self._H.dtype)

    #
    def _matvec(self, v):
        # left multiplication
        return self._z * v - self._H.dot(v, **self._kwargs)

    #
    def _rmatvec(self, v):
        # right multiplication
        return self._z.conj() * v - self._H.dot(v, **self._kwargs)


#
##### calculate action without constructing the Hamiltonian matrix
#
on_the_fly = (
    False  # toggles between using a `hamiltnian` or a `quantum_LinearOperator` object
)
# chain length
L = 12
# on-site spin size
S = "1/2"
# translation transformation on the lattice sites [0,...,L-1]
T = (np.arange(L) + 1) % L
# this example does not work under these conditions because ground-state sector is not q=0
if (L // 2) % 2 != 0:
    raise ValueError(
        "Example requires modifications for Heisenberg chains with L=4*n+2."
    )
if L % 2 != 0:
    raise ValueError(
        "Example requires modifications for Heisenberg chains with odd number of sites."
    )
# construct basis
basis0 = spin_basis_general(L, S=S, m=0, pauli=False, kblock=(T, 0))
# construct static list for Heisenberg chain
Jzz_list = [[1.0, i, (i + 1) % L] for i in range(L)]
Jxy_list = [[0.5, i, (i + 1) % L] for i in range(L)]
static = [["zz", Jzz_list], ["+-", Jxy_list], ["-+", Jxy_list]]
# construct operator for Hamiltonian in the ground state sector
if on_the_fly:
    H0 = quantum_LinearOperator(static, basis=basis0, dtype=np.float64)
else:
    H0 = hamiltonian(static, [], basis=basis0, dtype=np.float64)
# calculate ground state
[E0], psi0 = H0.eigsh(k=1, which="SA")
psi0 = psi0.ravel()
# define all possible momentum sectors excluding q=L/2 (pi-momentum) where the peak is abnormally large
qs = np.arange(-L // 2 + 1, L // 2, 1)
# define frequencies to calculate spectral function for
omegas = np.arange(0, 4, 0.05)
# spectral peaks broadening factor
eta = 0.1
# allocate arrays to store data
Gzz = np.zeros(omegas.shape + qs.shape, dtype=np.complex128)
Gpm = np.zeros(omegas.shape + qs.shape, dtype=np.complex128)
# loop over momentum sectors
for j, q in enumerate(qs):
    print("computing momentum block q={}".format(q))
    # define block
    block = dict(qblock=(T, q))
    #
    ####################################################################
    #
    # ---------------------------------------------------- #
    #                 calculation but for SzSz             #
    # ---------------------------------------------------- #
    #
    # define operator list for Op_shift_sector
    f = lambda i: np.exp(-2j * np.pi * q * i / L) / np.sqrt(L)
    Op_list = [["z", [i], f(i)] for i in range(L)]
    # define basis
    basisq = spin_basis_general(L, S=S, m=0, pauli=False, **block)
    # define operators in the q-momentum sector
    if on_the_fly:
        Hq = quantum_LinearOperator(
            static,
            basis=basisq,
            dtype=np.complex128,
            check_symm=False,
            check_pcon=False,
            check_herm=False,
        )
    else:
        Hq = hamiltonian(
            static,
            [],
            basis=basisq,
            dtype=np.complex128,
            check_symm=False,
            check_pcon=False,
            check_herm=False,
        )
    # shift sectors
    psiA = basisq.Op_shift_sector(basis0, Op_list, psi0)
    #
    ### apply vector correction method
    #
    # solve (z-H)|x> = |A> solve for |x>  using iterative solver for each omega
    if np.linalg.norm(psiA, ord=np.inf) < 1e-10:
        # since the vector is zero, we can skip the calculation as the result will be zero
        continue
        
    for i, omega in enumerate(omegas):
        lhs = LHS(Hq, omega, eta, E0)
        
        x, *_ = sp.linalg.bicgstab(lhs, psiA, maxiter=1000)
        Gzz[i, j] = -np.vdot(psiA, x) / np.pi
    #
    #####################################################################
    #
    # ---------------------------------------------------- #
    #            same calculation but for S-S+             #
    # ---------------------------------------------------- #
    #
    # divide by extra sqrt(2) to get extra factor of 1/2 when taking sandwich: needed since H = 1/2 (S^+_i S^-_j + h.c.) + S^z_j S^z_j
    f = lambda i: np.exp(-2j * np.pi * q * i / L) * np.sqrt(1.0 / (2 * L))
    Op_list = [["+", [i], f(i)] for i in range(L)]
    # change S_z projection up by S for action of S+ operator
    S_z_tot = 0 + eval(S)
    # calculate magnetization density from S_z_tot as defined in the documentation
    # m = S_z_tot / (S * N), S: local spin (as number), N: total spins
    m = S_z_tot / (eval(S) * L)
    basisq = spin_basis_general(L, S=S, m=m, pauli=False, **block)
    # define operators in the q-momentum sector
    if on_the_fly:
        Hq = quantum_LinearOperator(
            static,
            basis=basisq,
            dtype=np.complex128,
            check_symm=False,
            check_pcon=False,
            check_herm=False,
        )
    else:
        Hq = hamiltonian(
            static,
            [],
            basis=basisq,
            dtype=np.complex128,
            check_symm=False,
            check_pcon=False,
            check_herm=False,
        )
    # shift sectors
    psiA = basisq.Op_shift_sector(basis0, Op_list, psi0)
    #
    ### apply vector correction method
    #
    # solve (z-H)|x> = |A> solve for |x>  using iterative solver for each omega
    for i, omega in enumerate(omegas):
        lhs = LHS(Hq, omega, eta, E0)
        x, *_ = sp.linalg.bicg(lhs, psiA)
        Gpm[i, j] = -np.vdot(psiA, x) / np.pi
#
##### plot results
#
ks = 2 * np.pi * qs / L  # compute physical momentum values
f, (ax1, ax2) = plt.subplots(2, 1, figsize=(3, 4), sharex=True)
ax1.pcolormesh(ks, omegas, Gzz.imag, shading="nearest")
ax2.pcolormesh(ks, omegas, Gpm.imag, shading="nearest")
ax2.set_xlabel("$k$")
ax2.set_ylabel("$\\omega$")
ax1.set_ylabel("$\\omega$")
ax1.set_title("$G_{zz}(\\omega,k)$")
ax2.set_title("$G_{+-}(\\omega,k)$")
# plt.show()
plt.close()