Autocorrelation functions using symmetries

This example demonstrates how to use the general basis function method Op_shift_sector() to compute autocorrelation functions of operators in the Heisenberg model on a 1d chain.

\[H = J\sum_{j=0}^{L-1} S^+_{j+1}S^-_{j} + \mathrm{h.c.} + S^z_{j+1}S^z_j,\]

where \(S_j\) is the spin-1/2 operator on lattice site \(j\); we use periodic boundary conditions. We are interested in the following autocorrelation function:

\[C(t) = \langle\psi_\mathrm{GS}|\mathcal{O}^\dagger(t)\mathcal{O}(0)|\psi_\mathrm{GS}\rangle = \langle\psi_\mathrm{GS}|e^{+i t H}\mathcal{O}^\dagger e^{-i t H}\mathcal{O}|\psi_\mathrm{GS}\rangle\]

where \(\mathcal{O} = \sqrt{2}S^z_{j=0}\).

For the purpose of computing the autocorrelation function, it is advantageous to split the calculation as follows:

\[C(t) = \langle\psi_\mathrm{GS}(t)|\mathcal{O}^\dagger|(\mathcal{O}\psi_\mathrm{GS})(t)\rangle,\]

where \(|\psi_\mathrm{GS}(t)\rangle = e^{-i t H}|\psi_\mathrm{GS}\rangle\) (this is a trivial phase factor, but we keep it here for generality), and \(|(\mathcal{O}\psi_\mathrm{GS})(t)\rangle = e^{-i t H}\mathcal{O}|\psi_\mathrm{GS}\rangle\).

In the example below, we compute \(C(t)\) (i) in real space, and (ii) in momentum space.

  1. The real space calculation is straightforward, but it does not make use of the symmetries in the Heisenberg Hamiltonian.

  2. The momentum-space calculation is interesting, because the operator \(\mathcal{O}\) carries momentum itself; thus, when it acts on the (time-evolved) ground state, it changes the momentum of the state. In QuSpin, this operation can be done using the general basis class method Op_shift_sector(), provided we know ahead of time exactly by how much the momentum of the state will be changed. To understand how the symmetry calculation works, consider the more general nonequal-space, nonequal-time correlation function:

\[C_r(t) = \langle\psi_\mathrm{GS}|\mathcal{O}_r^\dagger(t)\mathcal{O}_0(0)|\psi_\mathrm{GS}\rangle = \frac{1}{L}\sum_{j=0}^{L-1}\langle\psi_\mathrm{GS}|\mathcal{O}_{j+r}^\dagger(t)\mathcal{O}_{j}(0)|\psi_\mathrm{GS}\rangle,\]

where in the second equality we explicitly used translation invariance. Using the Fourier transform \(\mathcal{O}_q(t) = 1/\sqrt{L}\sum_{j=0}^{L-1}\mathrm{e}^{-i \frac{2\pi q}{L} j}\mathcal{O}_j(t)\), we arrive at

\[C_r(t) = \frac{1}{L}\sum_{q=0}^{L-1} \mathrm{e}^{+i\frac{2\pi q}{L} r}\mathcal{O}^\dagger_q(t)\mathcal{O}_q(0).\]

Substituting the Fourier transform, and the epression for \(\mathcal{O}_j=\sqrt{2}S^z_j\) from above, setting \(r=0\), we arrive at

\[C_{r=0}(t) = \sum_{q=0}^{L-1} \left(\frac{1}{L}\sum_{j=0}^{L-1}\mathrm{e}^{-i \frac{2\pi q}{L} j} \sqrt{2}S^z_{j}(t) \right) \times \left(\frac{1}{L}\sum_{j'=0}^{L-1}\mathrm{e}^{-i \frac{2\pi q}{L} j' } \sqrt{2}S^z_{j'}(0) \right)\]

which is the expression we use in the code snippet below (note that since \(S^z\) is hermitian, it does not matter whether we use \(\mathrm{e}^{-i\dots}\) or \(\mathrm{e}^{+i\dots}\) here).

More generally, the operator Fourier decomposition of an operator \(\mathcal{O}_l\) with respect to any discrete symmetry transformation \(Q\) of periodicity/cyclicity \(m_Q\) (\(Q^{m_Q}=1\)), and eigenvalues \(\{\exp(-2\pi i q/m_Q)\}_{q=0}^{m_Q}\), is given by:

\[\mathcal{O}_{q} = \frac{1}{\sqrt m_Q}\sum_{j=0}^{m_Q-1} \mathrm{e}^{-i \frac{2\pi q}{m_Q} j} (Q^j)^\dagger \mathcal{O}_{l} Q^j.\]

For instance, if \(Q\) is the translation operator then \((Q^j)^\dagger \mathcal{O}_{l} Q^j = \mathcal{O}_{l+j}\); if \(Q\) is the reflection about the middle of the chain: \((Q)^\dagger \mathcal{O}_{l} Q = \mathcal{O}_{L-1-l}\), etc. The most general symmetry expression for the correlator then reads

\[\begin{split}C_{r}(t) &= \langle\psi_\mathrm{GS}|\mathcal{O}_{l+r}^\dagger(t)\mathcal{O}_l(0)|\psi_\mathrm{GS}\rangle = \frac{1}{m_Q}\sum_{j=0}^{m_Q-1}\langle\psi_\mathrm{GS}|(Q^j)^\dagger\mathcal{O}_{r+l}^\dagger(t)Q^j(Q^j)^\dagger\mathcal{O}_{l}(0)Q^j|\psi_\mathrm{GS}\rangle \\ &= \sum_{q=0}^{m_Q-1} \mathrm{e}^{+i\frac{2\pi q}{L} r} \left(\frac{1}{m_Q}\sum_{j=0}^{m_Q-1}\mathrm{e}^{+i \frac{2\pi q}{L} j} (Q^j)^\dagger \mathcal{O}_{l}^\dagger(t) Q^j \right) \times \left(\frac{1}{m_Q}\sum_{j'=0}^{m_Q-1}\mathrm{e}^{-i \frac{2\pi q}{m_Q} j' } (Q^{j'})^\dagger \mathcal{O}_{0}(l) Q^{j'} \right)\end{split}\]

This allows to exploit more symmetries of the Heisenberg model, if needed. An example of how this works for Op_shift_sector(), for reflection symmetry, is shown here.

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
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 19                                #	
# This exampled shows how to use the Op_shift_sector() method of the   #
# general basis classes to compute autocorrelation functions in the    #
# Heisenberg model.                                                    #
########################################################################
from quspin.basis import spin_basis_general
from quspin.operators import hamiltonian
from quspin.tools.evolution import expm_multiply_parallel
import numpy as np
import matplotlib.pyplot as plt # plotting library
#
##### Heisenberg model parameers
#
L=10
#
# hamiltonian parameters
J_list = [[1.0,i,(i+1)%L] for i in range(L)]
static = [[op,J_list] for op in ["-+","+-","zz"]]
# time vector
times = np.linspace(0.0,5.0,101)
#
##### straightforward autocorrelation function without using symmetries
#	
def auto_correlator(L,times,S="1/2"):
	# construct basis in zero magnetization sector: no lattice symmetries
	basis = spin_basis_general(L,S=S,m=0,pauli=False)
	# define Heisenberg Hamiltonian
	no_checks = dict(check_symm=False,check_herm=False,check_pcon=False)
	H = hamiltonian(static,[],basis=basis,dtype=np.float64,**no_checks)
	# compute GS
	E,V = H.eigsh(k=1,which="SA")
	psi_GS = V[:,0]
	# evolve GS under H (gives a trivial phase factor)
	psi_GS_t = H.evolve(psi_GS,0.0,times)
	#
	###### define operator O to compute the autocorrelation function of
	#
	op_list = [["z",[0],np.sqrt(2.0)]]
	# use inplace_Op to apply operator O on psi_GS
	Opsi_GS = basis.inplace_Op(psi_GS,op_list,np.float64)
	# time evolve Opsi_GS under H
	Opsi_GS_t = H.evolve(Opsi_GS,0.0,times)
	# apply operator O on time-evolved psi_t
	O_psi_GS_t = basis.inplace_Op(psi_GS_t,op_list,np.float64)
	# compute autocorrelator
	C_t = np.einsum("ij,ij->j",O_psi_GS_t.conj(),Opsi_GS_t)
	#
	return C_t
#
def auto_correlator_symm(L,times,S="1/2"):
	# define momentum p sector of the GS of the Heisenberg Hamiltonian
	if (L//2)%2:
		p = L//2 # corresponds to momentum pi
		dtype = np.complex128
	else:
		p = 0
		dtype = np.float64
	#
	# define translation operator
	T = (np.arange(L)+1)%L
	# compute the basis in the momentum sector of the GS of the Heisenberg model
	basis_p = spin_basis_general(L,S=S,m=0,kblock=(T,p),pauli=False)
	# define Heisenberg Hamiltonian 
	no_checks = dict(check_symm=False,check_herm=False,check_pcon=False)
	H = hamiltonian(static,[],basis=basis_p,dtype=dtype,**no_checks)
	# compute GS 
	E,V = H.eigsh(k=1,which="SA")
	psi_GS = V[:,0]
	# evolve GS under symmetry-reduced H (gives a trivial phase factor)
	psi_GS_t = H.evolve(psi_GS,0,times)
	#
	##### compute autocorrelation function foe every momentum sector
	Cq_t = np.zeros((times.shape[0],L),dtype=np.complex128)
	#
	for q in range(L): # sum over symmetry sectors
		#
		###### define operator O_q, sum over lattice sites
		op_list = [["z",[j],(np.sqrt(2.0)/L)*np.exp(-1j*2.0*np.pi*q*j/L)] for j in range(L)]
		# compute basis in the (q+p)-momentum sector (the total momentum of O_q|psi_GS> is q+p)
		basis_q = spin_basis_general(L,S=S,m=0,kblock=(T,p+q),pauli=False)
		# define Hamiltonian in the q-momentum sector
		Hq = hamiltonian(static,[],basis=basis_q,dtype=np.complex128,**no_checks)
		# use Op_shift_sector apply operator O_q to GS; the momentum of the new state is p+q
		Opsi_GS = basis_q.Op_shift_sector(basis_p,op_list,psi_GS)
		# time evolve Opsi_GS under H_q
		Opsi_GS_t = Hq.evolve(Opsi_GS,0.0,times)
		# apply operator O on time-evolved psi_t
		O_psi_GS_t = basis_q.Op_shift_sector(basis_p,op_list,psi_GS_t)
		# compute autocorrelator for every momentum sector
		Cq_t[...,q] = np.einsum("ij,ij->j",O_psi_GS_t.conj(),Opsi_GS_t)
	#
	return np.sum(Cq_t,axis=1) # sum over momentum sectors
#
##### compute autocorrelation function
C_t = auto_correlator(L,times)
C_t_symm = auto_correlator_symm(L,times)
#
##### plot result
#
plt.plot(times,C_t.real,'-b',label='no symm.: $\\mathrm{Re}\\;C(t)$')
plt.plot(times,C_t.imag,'-r',label='no symm.: $\\mathrm{Im}\\;C(t)$')
#
plt.plot(times,C_t_symm.real,'ob',label='symm.: $\\mathrm{Re}\\;C(t)$')
plt.plot(times,C_t_symm.imag,'or',label='symm.: $\\mathrm{Im}\\;C(t)$')
#
plt.xlabel('time $t$')
plt.legend()
#
#plt.show()