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.
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:
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:
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.
The real space calculation is straightforward, but it does not make use of the symmetries in the Heisenberg Hamiltonian.
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:
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
Substituting the Fourier transform, and the epression for \(\mathcal{O}_j=\sqrt{2}S^z_j\) from above, setting \(r=0\), we arrive at
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:
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
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¶
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()
|