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
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
from __future__ import print_function, division
#
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
#
quspin_path = os.path.join(os.getcwd(),"../../")
sys.path.insert(0,quspin_path)
#######################################################################
#                            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
	for i,omega in enumerate(omegas):
		lhs = LHS(Hq,omega,eta,E0)
		x,*_ = sp.linalg.bicg(lhs,psiA)
		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()