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\rangle\) can be written as:

\[C(\omega) = -\frac{1}{\pi}\langle 0|A^\dagger\frac{1}{\omega+i\eta + E_0 - H} A|0\rangle,\]

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

\[(\omega+i\eta + E_0 - H)|x(\omega)\rangle = |A\rangle.\]

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 \((\omega \pm i\eta + E_0 - H)|u\rangle\) for an arbitrary state \(|u\rangle\). Now we can use \(|x(\omega)\rangle\) to calculate the spectral function a given value of \(\omega\) as follows:

\[C(\omega) = -\frac{1}{\pi}\langle A|x(\omega)\rangle.\]

Note that this method can also be used to calculate spectral functions for any state by replacing \(|0\rangle\rightarrow|\psi\rangle\). This may be useful for calculating spectral functions away from equilibrium by using \(|\psi(t)\rangle\) 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=4n\) 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_{zz}(\omega,q) = \langle 0|S^{z}_{-q}\frac{1}{\omega+i\eta + E_0 - H}S^{z}_q|0\rangle\]
\[G_{+-}(\omega,q) = \langle 0|S^{-}_{-q}\frac{1}{\omega+i\eta + E_0 - H}S^{+}_q|0\rangle\]

where we have defined:

\[S^{z}_q = \frac{1}{\sqrt{L}}\sum_{r=0}^{L-1} \exp\left(-i\frac{2\pi q r}{L}\right) S^z_r\]
\[S^{\pm}_q = \frac{1}{\sqrt{2L}}\sum_{r=0}^{L-1} \exp\left(-i\frac{2\pi q r}{L}\right) S^{\pm}_r\]

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\) (\(\pi\)-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

  1#
  2import sys, os
  3
  4os.environ["KMP_DUPLICATE_LIB_OK"] = (
  5    "True"  # uncomment this line if omp error occurs on OSX for python 3
  6)
  7os.environ["OMP_NUM_THREADS"] = "1"  # set number of OpenMP threads to run in parallel
  8os.environ["MKL_NUM_THREADS"] = "1"  # set number of MKL threads to run in parallel
  9#
 10
 11#######################################################################
 12#                            example 26                               #
 13# This example shows how to use the `Op-shit_sector` method of the    #
 14# general basis class to compute spectral functions using symmetries. #
 15#######################################################################
 16from quspin.basis import spin_basis_general
 17from quspin.operators import hamiltonian, quantum_LinearOperator
 18import scipy.sparse as sp
 19import numexpr, cProfile
 20import numpy as np
 21import matplotlib.pyplot as plt
 22
 23
 24#
 25#
 26# define custom LinearOperator object that generates the left hand side of the equation.
 27#
 28class LHS(sp.linalg.LinearOperator):
 29    #
 30    def __init__(self, H, omega, eta, E0, kwargs={}):
 31        self._H = H  # Hamiltonian
 32        self._z = omega + 1j * eta + E0  # complex energy
 33        self._kwargs = kwargs  # arguments
 34
 35    #
 36    @property
 37    def shape(self):
 38        return (self._H.Ns, self._H.Ns)
 39
 40    #
 41    @property
 42    def dtype(self):
 43        return np.dtype(self._H.dtype)
 44
 45    #
 46    def _matvec(self, v):
 47        # left multiplication
 48        return self._z * v - self._H.dot(v, **self._kwargs)
 49
 50    #
 51    def _rmatvec(self, v):
 52        # right multiplication
 53        return self._z.conj() * v - self._H.dot(v, **self._kwargs)
 54
 55
 56#
 57##### calculate action without constructing the Hamiltonian matrix
 58#
 59on_the_fly = (
 60    False  # toggles between using a `hamiltnian` or a `quantum_LinearOperator` object
 61)
 62# chain length
 63L = 12
 64# on-site spin size
 65S = "1/2"
 66# translation transformation on the lattice sites [0,...,L-1]
 67T = (np.arange(L) + 1) % L
 68# this example does not work under these conditions because ground-state sector is not q=0
 69if (L // 2) % 2 != 0:
 70    raise ValueError(
 71        "Example requires modifications for Heisenberg chains with L=4*n+2."
 72    )
 73if L % 2 != 0:
 74    raise ValueError(
 75        "Example requires modifications for Heisenberg chains with odd number of sites."
 76    )
 77# construct basis
 78basis0 = spin_basis_general(L, S=S, m=0, pauli=False, kblock=(T, 0))
 79# construct static list for Heisenberg chain
 80Jzz_list = [[1.0, i, (i + 1) % L] for i in range(L)]
 81Jxy_list = [[0.5, i, (i + 1) % L] for i in range(L)]
 82static = [["zz", Jzz_list], ["+-", Jxy_list], ["-+", Jxy_list]]
 83# construct operator for Hamiltonian in the ground state sector
 84if on_the_fly:
 85    H0 = quantum_LinearOperator(static, basis=basis0, dtype=np.float64)
 86else:
 87    H0 = hamiltonian(static, [], basis=basis0, dtype=np.float64)
 88# calculate ground state
 89[E0], psi0 = H0.eigsh(k=1, which="SA")
 90psi0 = psi0.ravel()
 91# define all possible momentum sectors excluding q=L/2 (pi-momentum) where the peak is abnormally large
 92qs = np.arange(-L // 2 + 1, L // 2, 1)
 93# define frequencies to calculate spectral function for
 94omegas = np.arange(0, 4, 0.05)
 95# spectral peaks broadening factor
 96eta = 0.1
 97# allocate arrays to store data
 98Gzz = np.zeros(omegas.shape + qs.shape, dtype=np.complex128)
 99Gpm = np.zeros(omegas.shape + qs.shape, dtype=np.complex128)
100# loop over momentum sectors
101for j, q in enumerate(qs):
102    print("computing momentum block q={}".format(q))
103    # define block
104    block = dict(qblock=(T, q))
105    #
106    ####################################################################
107    #
108    # ---------------------------------------------------- #
109    #                 calculation but for SzSz             #
110    # ---------------------------------------------------- #
111    #
112    # define operator list for Op_shift_sector
113    f = lambda i: np.exp(-2j * np.pi * q * i / L) / np.sqrt(L)
114    Op_list = [["z", [i], f(i)] for i in range(L)]
115    # define basis
116    basisq = spin_basis_general(L, S=S, m=0, pauli=False, **block)
117    # define operators in the q-momentum sector
118    if on_the_fly:
119        Hq = quantum_LinearOperator(
120            static,
121            basis=basisq,
122            dtype=np.complex128,
123            check_symm=False,
124            check_pcon=False,
125            check_herm=False,
126        )
127    else:
128        Hq = hamiltonian(
129            static,
130            [],
131            basis=basisq,
132            dtype=np.complex128,
133            check_symm=False,
134            check_pcon=False,
135            check_herm=False,
136        )
137    # shift sectors
138    psiA = basisq.Op_shift_sector(basis0, Op_list, psi0)
139    #
140    ### apply vector correction method
141    #
142    # solve (z-H)|x> = |A> solve for |x>  using iterative solver for each omega
143    if np.linalg.norm(psiA, ord=np.inf) < 1e-10:
144        # since the vector is zero, we can skip the calculation as the result will be zero
145        continue
146        
147    for i, omega in enumerate(omegas):
148        lhs = LHS(Hq, omega, eta, E0)
149        
150        x, *_ = sp.linalg.bicgstab(lhs, psiA, maxiter=1000)
151        Gzz[i, j] = -np.vdot(psiA, x) / np.pi
152    #
153    #####################################################################
154    #
155    # ---------------------------------------------------- #
156    #            same calculation but for S-S+             #
157    # ---------------------------------------------------- #
158    #
159    # 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
160    f = lambda i: np.exp(-2j * np.pi * q * i / L) * np.sqrt(1.0 / (2 * L))
161    Op_list = [["+", [i], f(i)] for i in range(L)]
162    # change S_z projection up by S for action of S+ operator
163    S_z_tot = 0 + eval(S)
164    # calculate magnetization density from S_z_tot as defined in the documentation
165    # m = S_z_tot / (S * N), S: local spin (as number), N: total spins
166    m = S_z_tot / (eval(S) * L)
167    basisq = spin_basis_general(L, S=S, m=m, pauli=False, **block)
168    # define operators in the q-momentum sector
169    if on_the_fly:
170        Hq = quantum_LinearOperator(
171            static,
172            basis=basisq,
173            dtype=np.complex128,
174            check_symm=False,
175            check_pcon=False,
176            check_herm=False,
177        )
178    else:
179        Hq = hamiltonian(
180            static,
181            [],
182            basis=basisq,
183            dtype=np.complex128,
184            check_symm=False,
185            check_pcon=False,
186            check_herm=False,
187        )
188    # shift sectors
189    psiA = basisq.Op_shift_sector(basis0, Op_list, psi0)
190    #
191    ### apply vector correction method
192    #
193    # solve (z-H)|x> = |A> solve for |x>  using iterative solver for each omega
194    for i, omega in enumerate(omegas):
195        lhs = LHS(Hq, omega, eta, E0)
196        x, *_ = sp.linalg.bicg(lhs, psiA)
197        Gpm[i, j] = -np.vdot(psiA, x) / np.pi
198#
199##### plot results
200#
201ks = 2 * np.pi * qs / L  # compute physical momentum values
202f, (ax1, ax2) = plt.subplots(2, 1, figsize=(3, 4), sharex=True)
203ax1.pcolormesh(ks, omegas, Gzz.imag, shading="nearest")
204ax2.pcolormesh(ks, omegas, Gpm.imag, shading="nearest")
205ax2.set_xlabel("$k$")
206ax2.set_ylabel("$\\omega$")
207ax1.set_ylabel("$\\omega$")
208ax1.set_title("$G_{zz}(\\omega,k)$")
209ax2.set_title("$G_{+-}(\\omega,k)$")
210# plt.show()
211plt.close()