Kitaev’s honeycomb Hamiltonian in vortex sectors given by the plaquette operators

This example shows how one can employ QuSpin’s user_basis class to construct the Kitaev honeycomb model Hamiltonian in the spectral sectors generated by the Wilson loop operators W. The transformation generated by W cannot be defined using QuSpin’s basis_general classes.

Consider the Kitaev honeycomb model (see reference for details of the definition):

\[H_\mathrm{Kitaev} = \sum_{\langle i,j\rangle_\alpha}J^\alpha_{ij} \sigma^\alpha_i \sigma^\alpha_j.\]

where \(\sigma\) are the Pauli matrices, and we consider 4 unit cells (8 sites) with periodic boundary conditions.

It can be shown that this Hamiltonian commutes with all Wilson loop plauqette operators, e.g., \(W=\sigma^x_0\sigma^y_1\sigma^z_2\sigma^x_3\sigma^y_4\sigma^z_5\), etc. Hence, the \(W\)-operators generate a transformation according to which \(H_\mathrm{Kitaev}\) can be block diagonalized; the corresponding sectors can be labelled by the eigenvalues of \(W\) which take the values \(\pm 1\), and determine the vortex configuration of the states in a given sector.

Here, we want to construct the Hamiltonian \(H_\mathrm{Kitaev}\) in a given vortex configuration sector using QuSpin.

It turns out, this cannot be achieved using the basis_general classes, because the operation induced by \(W\) on the computational basis states cannot be represented by a combination of lattice sites permutations and spin inversion operations. Indeed, the operators \(\sigma^z_2, \sigma^z_5\) multiply the resulting state by a state-dependent sign, while the operators \(\sigma^y_1, \sigma^y_4\) multiply the resulting state by a state-dependent factor \(\pm i\).

Notes:

  • the procedure described below is generic, and can be applied to any symmetry, provided that (in the basis chosen) the symmetry is represented by a real-valued operator.

  • the code below works only for 8 lattice sites with periodic boundary conditions; however, it should be clear how to generalize it to other system sizes.

  • this code shows an interesting example of using the sign pointer variable, inherent to the symmetry map function template of the user_basis, for spin systems [typically, sign is used to accomodate fermionic signs in fermionic models].

  • there are two more Wilson loops that wrap around the two directions of the torus; we leave it to the reader to code them up.

Please consult this tutorial post – A tutorial on QuSpin’s user_basis – for more detailed explanations on using the user_basis class.

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
166
167
168
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']='4' # set number of OpenMP threads to run in parallel
os.environ['MKL_NUM_THREADS']='4' # set number of MKL threads to run in parallel
#
quspin_path = os.path.join(os.getcwd(),"../../")
sys.path.insert(0,quspin_path)
#########################################################################
#                             example 28                                #
#  In this script we demonstrate how to use QuSpin's user_basis         #
#  to define symmetries, which are not supported by the basis_general   #
#  classes. We take an 8-site honeycomb lattice with PBC and define     #
#  the Kitaev model in the spectral sectors given by the Wilson loop /  #
#  plaquette operators W.                                               #
#########################################################################
from quspin.operators import hamiltonian # Hamiltonians and operators
from quspin.basis.user import user_basis # Hilbert space user basis
from quspin.basis.user import next_state_sig_32,op_sig_32,map_sig_32,count_particles_sig_32 # user basis data types signatures
from numba import carray,cfunc,jit # numba helper functions
from numba import uint32,int32 # numba data types
import numpy as np
#
N=8 # lattice sites
assert(N==8) # code below works only for N=8!
#
############   create spin-1/2 user basis object   #############
#
######  function to call when applying operators
@cfunc(op_sig_32,
	locals=dict(s=int32,n=int32,b=uint32), )
def op(op_struct_ptr,op_str,site_ind,N,args):
	# using struct pointer to pass op_struct_ptr back to C++ see numba Records
	op_struct = carray(op_struct_ptr,1)[0]
	err = 0
	#
	site_ind = N - site_ind - 1 # convention for QuSpin for mapping from bits to sites.
	n = (op_struct.state>>site_ind)&1 # either 0 or 1
	s = (((op_struct.state>>site_ind)&1)<<1)-1 # either -1 or 1
	b = (1<<site_ind)
	#
	if op_str==120: # "x" is integer value 120 = ord("x")
		op_struct.state ^= b

	elif op_str==121: # "y" is integer value 120 = ord("y")
		op_struct.state ^= b
		op_struct.matrix_ele *= 1.0j*s

	elif op_str==43: # "+" is integer value 43 = ord("+")
		if n: op_struct.matrix_ele = 0
		else: op_struct.state ^= b # create spin

	elif op_str==45: # "-" is integer value 45 = ord("-")
		if n: op_struct.state ^= b # destroy spin
		else: op_struct.matrix_ele = 0

	elif op_str==122: # "z" is integer value 120 = ord("z")
		op_struct.matrix_ele *= s

	elif op_str==110: # "n" is integer value 110 = ord("n")
		op_struct.matrix_ele *= n

	elif op_str==73: # "I" is integer value 73 = ord("I")
		pass

	else:
		op_struct.matrix_ele = 0
		err = -1
	#
	return err
op_args=np.array([],dtype=np.uint32)
#
######  define symmetry maps
#
@jit(uint32(uint32,uint32,int32),)
def _compute_occupation(state,site_ind,N):
	# auxiliary function to check occupation of state at site_ind
	# CAUTION: 32-bit integers code only!
	site_ind = N - site_ind - 1 # compute lattice site in reversed bit configuration (cf QuSpin convention for mapping from bits to sites)
	return (state>>site_ind)&1 # occupation: either 0 or 1
#
@jit(uint32(uint32,uint32,int32),locals=dict(b=uint32,))
def _flip_occupation(state,site_ind,N):
	# auxiliary function to flip occupation of state at site_ind
	# CAUTION: 32-bit integers code only!
	site_ind = N - site_ind - 1     # compute lattice site in reversed bit configuration (cf QuSpin convention for mapping from bits to sites)
	b = 1; b <<= site_ind   	    # compute a "mask" integer b which is 1 on site site_ind and zero elsewhere      			
	return state^b # flip occupation on site site_ind
#
@cfunc(map_sig_32,
	locals=dict(out=uint32,) )
def W_symm(state,N,sign_ptr,args):
	""" 
	applies plaquette operator W = "xyzxyz" on sites defined in args. 
	# CAUTION: 32-bit integers code only!

	- the two 'z'-operators in W are easiest to apply, since just return the occupation number in the spin basis.
	- the two 'x'-operators in W are alo easy -- they only flip the occupation
	- the two 'y'-operators in W do both: they flip the occupation and mutiply the state by +/-i (each) 
	"""
	assert(N==8) # works only for N=8!
	out = 0
	#
	# compute sign by applying operators from W on state
	sign_ptr[0] *= -1 + 2*_compute_occupation(state,args[1],N) # 'y'; factor i taken into account below
	sign_ptr[0] *= -1 + 2*_compute_occupation(state,args[2],N) # 'z'
	sign_ptr[0] *= -1 + 2*_compute_occupation(state,args[4],N) # 'y'; factor i taken into account below
	sign_ptr[0] *= -1 + 2*_compute_occupation(state,args[5],N) # 'z'
	sign_ptr[0] *= -1 # -1=i**2, yy 
	#
	# flip occupation of state on sites 0,1, 3,4
	out = _flip_occupation(state,args[0],N) # 'x'
	out = _flip_occupation(out,  args[1],N) # 'y'
	out = _flip_occupation(out,  args[3],N) # 'x'
	out = _flip_occupation(out,  args[4],N) # 'y'
	#
	return out
# the argument list stores the sites for the four W operators, defined clockwise
W_args_1=np.array([0,1,2,3,4,5],dtype=np.uint32) # plaquette 1
W_args_2=np.array([4,3,6,1,0,7],dtype=np.uint32) # plaquette 2
W_args_3=np.array([6,7,0,5,2,1],dtype=np.uint32) # plaquette 3
W_args_4=np.array([2,5,4,7,6,3],dtype=np.uint32) # plaquette 4
#
######  construct user_basis 
#
# define maps dict
q_nums=[0,0,0,0] # not all sectors contain a finite-number of states
maps = dict( W1_block=(W_symm,2,q_nums[0],W_args_1), 
			 W2_block=(W_symm,2,q_nums[1],W_args_2), 
			 W3_block=(W_symm,2,q_nums[2],W_args_3), 
			 W4_block=(W_symm,2,q_nums[3],W_args_4),
			) 
op_dict = dict(op=op,op_args=op_args)
# create user basis
basis = user_basis(np.uint32,N,op_dict,allowed_ops=set("+-xyznI"),sps=2,**maps)
#print(basis)
#
##### define Kitaev Hamiltonian
#
### coupling strengths
Jx,Jy,Jz=1.0,1.0,1.0 
#
### site-coupling lists
zz_int = [[Jz,0,1], [Jz,4,3], [Jz,2,5], [Jz,6,7] ] 
yy_int = [[Jy,2,3], [Jy,0,5], [Jy,6,1], [Jy,4,7] ]
xx_int = [[Jx,2,1], [Jx,4,5], [Jx,6,3], [Jx,0,7] ] 
#
### construct Hamitonian
no_checks=dict(check_pcon=False, check_herm=False, check_symm=False)
#
Hzz = hamiltonian([['zz',zz_int],], [], basis=basis, dtype=np.complex128, **no_checks) 
Hyy = hamiltonian([['yy',yy_int],], [], basis=basis, dtype=np.complex128, **no_checks) 
Hxx = hamiltonian([['xx',xx_int],], [], basis=basis, dtype=np.complex128, **no_checks) 
#
H_Kitaev = Hzz + Hyy  + Hxx
#
### diagonalize H_Kitaev
#
#E, V = H_Kitaev.eigsh(k=4, which='SA')
E, V = H_Kitaev.eigh()
print('energy spectrum:', E[:4])
#
### construct plaquette operator and check its e'state expectation
#
W_int = [[1.0, *list(W_args_1) ]]
W = hamiltonian([['xyzxyz',W_int],], [], basis=basis, dtype=np.float64, **no_checks) 
print('check expectation of W-operator:', W.expt_value(V[:,:4],).real)