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 – user_basis tutorial – for more detailed explanations on using the user_basis class.

Script

download 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"] = "4"  # set number of OpenMP threads to run in parallel
  8os.environ["MKL_NUM_THREADS"] = "4"  # set number of MKL threads to run in parallel
  9#
 10
 11#########################################################################
 12#                             example 28                                #
 13#  In this script we demonstrate how to use QuSpin's user_basis         #
 14#  to define symmetries, which are not supported by the basis_general   #
 15#  classes. We take an 8-site honeycomb lattice with PBC and define     #
 16#  the Kitaev model in the spectral sectors given by the Wilson loop /  #
 17#  plaquette operators W.                                               #
 18#########################################################################
 19from quspin.operators import hamiltonian  # Hamiltonians and operators
 20from quspin.basis.user import user_basis  # Hilbert space user basis
 21from quspin.basis.user import (
 22    next_state_sig_32,
 23    op_sig_32,
 24    map_sig_32,
 25    count_particles_sig_32,
 26)  # user basis data types signatures
 27from numba import carray, cfunc, jit  # numba helper functions
 28from numba import uint32, int32  # numba data types
 29import numpy as np
 30
 31#
 32N = 8  # lattice sites
 33assert N == 8  # code below works only for N=8!
 34
 35
 36#
 37############   create spin-1/2 user basis object   #############
 38#
 39######  function to call when applying operators
 40@cfunc(
 41    op_sig_32,
 42    locals=dict(s=int32, n=int32, b=uint32),
 43)
 44def op(op_struct_ptr, op_str, site_ind, N, args):
 45    # using struct pointer to pass op_struct_ptr back to C++ see numba Records
 46    op_struct = carray(op_struct_ptr, 1)[0]
 47    err = 0
 48    #
 49    site_ind = N - site_ind - 1  # convention for QuSpin for mapping from bits to sites.
 50    n = (op_struct.state >> site_ind) & 1  # either 0 or 1
 51    s = (((op_struct.state >> site_ind) & 1) << 1) - 1  # either -1 or 1
 52    b = 1 << site_ind
 53    #
 54    if op_str == 120:  # "x" is integer value 120 = ord("x")
 55        op_struct.state ^= b
 56
 57    elif op_str == 121:  # "y" is integer value 120 = ord("y")
 58        op_struct.state ^= b
 59        op_struct.matrix_ele *= 1.0j * s
 60
 61    elif op_str == 43:  # "+" is integer value 43 = ord("+")
 62        if n:
 63            op_struct.matrix_ele = 0
 64        else:
 65            op_struct.state ^= b  # create spin
 66
 67    elif op_str == 45:  # "-" is integer value 45 = ord("-")
 68        if n:
 69            op_struct.state ^= b  # destroy spin
 70        else:
 71            op_struct.matrix_ele = 0
 72
 73    elif op_str == 122:  # "z" is integer value 120 = ord("z")
 74        op_struct.matrix_ele *= s
 75
 76    elif op_str == 110:  # "n" is integer value 110 = ord("n")
 77        op_struct.matrix_ele *= n
 78
 79    elif op_str == 73:  # "I" is integer value 73 = ord("I")
 80        pass
 81
 82    else:
 83        op_struct.matrix_ele = 0
 84        err = -1
 85    #
 86    return err
 87
 88
 89op_args = np.array([], dtype=np.uint32)
 90
 91
 92#
 93######  define symmetry maps
 94#
 95@jit(
 96    uint32(uint32, uint32, int32),
 97)
 98def _compute_occupation(state, site_ind, N):
 99    # auxiliary function to check occupation of state at site_ind
100    # CAUTION: 32-bit integers code only!
101    site_ind = (
102        N - site_ind - 1
103    )  # compute lattice site in reversed bit configuration (cf QuSpin convention for mapping from bits to sites)
104    return (state >> site_ind) & 1  # occupation: either 0 or 1
105
106
107#
108@jit(
109    uint32(uint32, uint32, int32),
110    locals=dict(
111        b=uint32,
112    ),
113)
114def _flip_occupation(state, site_ind, N):
115    # auxiliary function to flip occupation of state at site_ind
116    # CAUTION: 32-bit integers code only!
117    site_ind = (
118        N - site_ind - 1
119    )  # compute lattice site in reversed bit configuration (cf QuSpin convention for mapping from bits to sites)
120    b = 1
121    b <<= site_ind  # compute a "mask" integer b which is 1 on site site_ind and zero elsewhere
122    return state ^ b  # flip occupation on site site_ind
123
124
125#
126@cfunc(
127    map_sig_32,
128    locals=dict(
129        out=uint32,
130    ),
131)
132def W_symm(state, N, sign_ptr, args):
133    """
134    applies plaquette operator W = "xyzxyz" on sites defined in args.
135    # CAUTION: 32-bit integers code only!
136
137    - the two 'z'-operators in W are easiest to apply, since just return the occupation number in the spin basis.
138    - the two 'x'-operators in W are alo easy -- they only flip the occupation
139    - the two 'y'-operators in W do both: they flip the occupation and mutiply the state by +/-i (each)
140    """
141    assert N == 8  # works only for N=8!
142    out = 0
143    #
144    # compute sign by applying operators from W on state
145    sign_ptr[0] *= -1 + 2 * _compute_occupation(
146        state, args[1], N
147    )  # 'y'; factor i taken into account below
148    sign_ptr[0] *= -1 + 2 * _compute_occupation(state, args[2], N)  # 'z'
149    sign_ptr[0] *= -1 + 2 * _compute_occupation(
150        state, args[4], N
151    )  # 'y'; factor i taken into account below
152    sign_ptr[0] *= -1 + 2 * _compute_occupation(state, args[5], N)  # 'z'
153    sign_ptr[0] *= -1  # -1=i**2, yy
154    #
155    # flip occupation of state on sites 0,1, 3,4
156    out = _flip_occupation(state, args[0], N)  # 'x'
157    out = _flip_occupation(out, args[1], N)  # 'y'
158    out = _flip_occupation(out, args[3], N)  # 'x'
159    out = _flip_occupation(out, args[4], N)  # 'y'
160    #
161    return out
162
163
164# the argument list stores the sites for the four W operators, defined clockwise
165W_args_1 = np.array([0, 1, 2, 3, 4, 5], dtype=np.uint32)  # plaquette 1
166W_args_2 = np.array([4, 3, 6, 1, 0, 7], dtype=np.uint32)  # plaquette 2
167W_args_3 = np.array([6, 7, 0, 5, 2, 1], dtype=np.uint32)  # plaquette 3
168W_args_4 = np.array([2, 5, 4, 7, 6, 3], dtype=np.uint32)  # plaquette 4
169#
170######  construct user_basis
171#
172# define maps dict
173q_nums = [0, 0, 0, 0]  # not all sectors contain a finite-number of states
174maps = dict(
175    W1_block=(W_symm, 2, q_nums[0], W_args_1),
176    W2_block=(W_symm, 2, q_nums[1], W_args_2),
177    W3_block=(W_symm, 2, q_nums[2], W_args_3),
178    W4_block=(W_symm, 2, q_nums[3], W_args_4),
179)
180op_dict = dict(op=op, op_args=op_args)
181# create user basis
182basis = user_basis(np.uint32, N, op_dict, allowed_ops=set("+-xyznI"), sps=2, **maps)
183# print(basis)
184#
185##### define Kitaev Hamiltonian
186#
187### coupling strengths
188Jx, Jy, Jz = 1.0, 1.0, 1.0
189#
190### site-coupling lists
191zz_int = [[Jz, 0, 1], [Jz, 4, 3], [Jz, 2, 5], [Jz, 6, 7]]
192yy_int = [[Jy, 2, 3], [Jy, 0, 5], [Jy, 6, 1], [Jy, 4, 7]]
193xx_int = [[Jx, 2, 1], [Jx, 4, 5], [Jx, 6, 3], [Jx, 0, 7]]
194#
195### construct Hamitonian
196no_checks = dict(check_pcon=False, check_herm=False, check_symm=False)
197#
198Hzz = hamiltonian(
199    [
200        ["zz", zz_int],
201    ],
202    [],
203    basis=basis,
204    dtype=np.complex128,
205    **no_checks,
206)
207Hyy = hamiltonian(
208    [
209        ["yy", yy_int],
210    ],
211    [],
212    basis=basis,
213    dtype=np.complex128,
214    **no_checks,
215)
216Hxx = hamiltonian(
217    [
218        ["xx", xx_int],
219    ],
220    [],
221    basis=basis,
222    dtype=np.complex128,
223    **no_checks,
224)
225#
226H_Kitaev = Hzz + Hyy + Hxx
227#
228### diagonalize H_Kitaev
229#
230# E, V = H_Kitaev.eigsh(k=4, which='SA')
231E, V = H_Kitaev.eigh()
232print("energy spectrum:", E[:4])
233#
234### construct plaquette operator and check its e'state expectation
235#
236W_int = [[1.0, *list(W_args_1)]]
237W = hamiltonian(
238    [
239        ["xyzxyz", W_int],
240    ],
241    [],
242    basis=basis,
243    dtype=np.float64,
244    **no_checks,
245)
246print(
247    "check expectation of W-operator:",
248    W.expt_value(
249        V[:, :4],
250    ).real,
251)