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#
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 19 #
13# This exampled shows how to use the Op_shift_sector() method of the #
14# general basis classes to compute autocorrelation functions in the #
15# Heisenberg model. #
16########################################################################
17from quspin.basis import spin_basis_general
18from quspin.operators import hamiltonian
19from quspin.tools.evolution import expm_multiply_parallel
20import numpy as np
21import matplotlib.pyplot as plt # plotting library
22
23#
24##### Heisenberg model parameers
25#
26L = 10
27#
28# hamiltonian parameters
29J_list = [[1.0, i, (i + 1) % L] for i in range(L)]
30static = [[op, J_list] for op in ["-+", "+-", "zz"]]
31# time vector
32times = np.linspace(0.0, 5.0, 101)
33
34
35#
36##### straightforward autocorrelation function without using symmetries
37#
38def auto_correlator(L, times, S="1/2"):
39 # construct basis in zero magnetization sector: no lattice symmetries
40 basis = spin_basis_general(L, S=S, m=0, pauli=False)
41 # define Heisenberg Hamiltonian
42 no_checks = dict(check_symm=False, check_herm=False, check_pcon=False)
43 H = hamiltonian(static, [], basis=basis, dtype=np.float64, **no_checks)
44 # compute GS
45 E, V = H.eigsh(k=1, which="SA")
46 psi_GS = V[:, 0]
47 # evolve GS under H (gives a trivial phase factor)
48 psi_GS_t = H.evolve(psi_GS, 0.0, times)
49 #
50 ###### define operator O to compute the autocorrelation function of
51 #
52 op_list = [["z", [0], np.sqrt(2.0)]]
53 # use inplace_Op to apply operator O on psi_GS
54 Opsi_GS = basis.inplace_Op(psi_GS, op_list, np.float64)
55 # time evolve Opsi_GS under H
56 Opsi_GS_t = H.evolve(Opsi_GS, 0.0, times)
57 # apply operator O on time-evolved psi_t
58 O_psi_GS_t = basis.inplace_Op(psi_GS_t, op_list, np.float64)
59 # compute autocorrelator
60 C_t = np.einsum("ij,ij->j", O_psi_GS_t.conj(), Opsi_GS_t)
61 #
62 return C_t
63
64
65#
66def auto_correlator_symm(L, times, S="1/2"):
67 # define momentum p sector of the GS of the Heisenberg Hamiltonian
68 if (L // 2) % 2:
69 p = L // 2 # corresponds to momentum pi
70 dtype = np.complex128
71 else:
72 p = 0
73 dtype = np.float64
74 #
75 # define translation operator
76 T = (np.arange(L) + 1) % L
77 # compute the basis in the momentum sector of the GS of the Heisenberg model
78 basis_p = spin_basis_general(L, S=S, m=0, kblock=(T, p), pauli=False)
79 # define Heisenberg Hamiltonian
80 no_checks = dict(check_symm=False, check_herm=False, check_pcon=False)
81 H = hamiltonian(static, [], basis=basis_p, dtype=dtype, **no_checks)
82 # compute GS
83 E, V = H.eigsh(k=1, which="SA")
84 psi_GS = V[:, 0]
85 # evolve GS under symmetry-reduced H (gives a trivial phase factor)
86 psi_GS_t = H.evolve(psi_GS, 0, times)
87 #
88 ##### compute autocorrelation function foe every momentum sector
89 Cq_t = np.zeros((times.shape[0], L), dtype=np.complex128)
90 #
91 for q in range(L): # sum over symmetry sectors
92 #
93 ###### define operator O_q, sum over lattice sites
94 op_list = [
95 ["z", [j], (np.sqrt(2.0) / L) * np.exp(-1j * 2.0 * np.pi * q * j / L)]
96 for j in range(L)
97 ]
98 # compute basis in the (q+p)-momentum sector (the total momentum of O_q|psi_GS> is q+p)
99 basis_q = spin_basis_general(L, S=S, m=0, kblock=(T, p + q), pauli=False)
100 # define Hamiltonian in the q-momentum sector
101 Hq = hamiltonian(static, [], basis=basis_q, dtype=np.complex128, **no_checks)
102 # use Op_shift_sector apply operator O_q to GS; the momentum of the new state is p+q
103 Opsi_GS = basis_q.Op_shift_sector(basis_p, op_list, psi_GS)
104 # time evolve Opsi_GS under H_q
105 Opsi_GS_t = Hq.evolve(Opsi_GS, 0.0, times)
106 # apply operator O on time-evolved psi_t
107 O_psi_GS_t = basis_q.Op_shift_sector(basis_p, op_list, psi_GS_t)
108 # compute autocorrelator for every momentum sector
109 Cq_t[..., q] = np.einsum("ij,ij->j", O_psi_GS_t.conj(), Opsi_GS_t)
110 #
111 return np.sum(Cq_t, axis=1) # sum over momentum sectors
112
113
114#
115##### compute autocorrelation function
116C_t = auto_correlator(L, times)
117C_t_symm = auto_correlator_symm(L, times)
118#
119##### plot result
120#
121plt.plot(times, C_t.real, "-b", label="no symm.: $\\mathrm{Re}\\;C(t)$")
122plt.plot(times, C_t.imag, "-r", label="no symm.: $\\mathrm{Im}\\;C(t)$")
123#
124plt.plot(times, C_t_symm.real, "ob", label="symm.: $\\mathrm{Re}\\;C(t)$")
125plt.plot(times, C_t_symm.imag, "or", label="symm.: $\\mathrm{Im}\\;C(t)$")
126#
127plt.xlabel("time $t$")
128plt.legend()
129#
130# plt.show()