Majorana Fermions: Spinless Fermi-Hubbard Model
In this example, we show how to use the user_basis class to define Majorana fermion operators. We then show how to construct the Hamiltonian for the spinless Fermi-Hubbard model in the Majorana representation. Starting from version 0.3.5 the Majorana operator strings “x” and “y” are also available in the *fermion_basis_general classes.
Consider first a single fermion described by the creation and annihilation operators \(c^\dagger,c\). One can decompose one complex fermion into two real-valued Majorana fermion modes \(c^x, c^y\), defined by the operators
which obey the fermion commutation relations \(\{c^x, c^y\}=0\), \(\left(c^\alpha\right)^2 = 1\), for \(\alpha=x,y\).
The inverse transformation is given by
Here, we choose to denote the Majorana operators by \(c^x, c^y\), due to a formal relation with the Pauli matrices: \(c^x = \sigma^x\), \(c^y = -\sigma^y\) (NB: the extra negative sign appearing for \(c^y\) is because of the standard convention for the definition of Majorana modes), \(c^\dagger = \sigma^+\), and \(c = \sigma^-\).
One can then generalize the Majorana decomposition to every site \(j\) of the lattice: \(c^x_j, c^y_j\).
To implement Majorana modes in QuSpin, we use the versatility of the user_basis to define the on-site action of the corresponding operators (please consult this post – user_basis tutorial – for more detailed explanations on using the user_basis class). To do this, we define new fermionic operator strings x and y, corresponding to the two Majorana modes \(c^x\) and :math:c^y`, respectively. The definition of x and y follows the exact same procedure, as in the spin-1/2 basis (cf. spin-1/2 Heisenberg model in 1d), with the notable difference that one has to accomodate for the sign, arising from counting the number of fermions up to the lattice site the operator is applied.
Having access to Majorana operator strings allows us to implement Hamiltonians in the Majorana representation. To demonstrate this, we use the spinless Fermi-Hubbard model:
where \(J\) is the hopping matrix element, and \(U\) is the nearest-neighbor interaction strength.
In the Majorana representation, the same Hamiltonian reads as
The script below uses the user_basis class to define an enlarged fermion basis which allows to implement both Hamiltonians. We then compare the two matrices. Note that the user_basis readily allows one to apply symmetries to the Hamiltonian in the Majorana representation.
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 24 #
13# This example shows how to use the `user_basis` to define Majorana operators. #
14################################################################################
15from quspin.operators import hamiltonian # Hamiltonians and operators
16from quspin.basis import spinless_fermion_basis_1d # Hilbert space fermion basis_1d
17from quspin.basis.user import user_basis # Hilbert space user basis
18from quspin.basis.user import (
19 next_state_sig_32,
20 op_sig_32,
21 map_sig_32,
22 count_particles_sig_32,
23) # user basis data types signatures
24from numba import carray, cfunc, jit # numba helper functions
25from numba import uint32, int32 # numba data types
26import numpy as np
27from scipy.special import comb
28
29np.set_printoptions(suppress="True", precision=6)
30#
31N = 6 # lattice sites
32
33
34#
35############ create soinless fermion user basis object #############
36#
37@jit(
38 uint32(uint32, uint32),
39 locals=dict(
40 f_count=uint32,
41 ),
42 nopython=True,
43 nogil=True,
44)
45def _count_particles_32(state, site_ind):
46 # auxiliary function to count number of fermions, i.e. 1's in bit configuration of the state, up to site site_ind
47 # CAUTION: 32-bit integers code only!
48 f_count = state & ((0x7FFFFFFF) >> (31 - site_ind))
49 f_count = f_count - ((f_count >> 1) & 0x55555555)
50 f_count = (f_count & 0x33333333) + ((f_count >> 2) & 0x33333333)
51 return (((f_count + (f_count >> 4)) & 0x0F0F0F0F) * 0x01010101) >> 24
52
53
54#
55@cfunc(
56 op_sig_32,
57 locals=dict(s=int32, sign=int32, n=int32, b=uint32, f_count=uint32),
58)
59def op(op_struct_ptr, op_str, site_ind, N, args):
60 # using struct pointer to pass op_struct_ptr back to C++ see numba Records
61 op_struct = carray(op_struct_ptr, 1)[0]
62 err = 0
63 #
64 site_ind = N - site_ind - 1 # convention for QuSpin for mapping from bits to sites.
65 #####
66 f_count = _count_particles_32(op_struct.state, site_ind)
67 #####
68 sign = -1 if f_count & 1 else 1
69 n = (op_struct.state >> site_ind) & 1 # either 0 or 1
70 b = 1 << site_ind
71 #
72 if op_str == 120: # "x" is integer value 120 = ord("x")
73 op_struct.state ^= b
74 op_struct.matrix_ele *= sign
75
76 elif op_str == 121: # "y" is integer value 120 = ord("y")
77 op_struct.state ^= b
78 op_struct.matrix_ele *= -1.0j * sign * ((n << 1) - 1)
79
80 elif op_str == 43: # "+" is integer value 43 = ord("+")
81 op_struct.matrix_ele *= 0.0 if n else sign
82 op_struct.state ^= b # create fermion
83
84 elif op_str == 45: # "-" is integer value 45 = ord("-")
85 op_struct.matrix_ele *= sign if n else 0.0
86 op_struct.state ^= b # create fermion
87
88 elif op_str == 110: # "n" is integer value 110 = ord("n")
89 op_struct.matrix_ele *= n
90
91 elif op_str == 73: # "I" is integer value 73 = ord("I")
92 pass
93
94 else:
95 op_struct.matrix_ele = 0
96 err = -1
97 #
98 return err
99
100
101op_args = np.array([], dtype=np.uint32)
102
103
104#
105###### define symmetry maps
106#
107@cfunc(
108 map_sig_32,
109 locals=dict(
110 shift=uint32,
111 xmax=uint32,
112 x1=uint32,
113 x2=uint32,
114 period=int32,
115 l=int32,
116 f_count1=int32,
117 f_count2=int32,
118 ),
119)
120def translation(x, N, sign_ptr, args):
121 """works for all system sizes N."""
122 shift = args[0] # translate state by shift sites
123 period = N # periodicity/cyclicity of translation
124 xmax = args[1]
125 #
126 l = (shift + period) % period
127 x1 = x >> (period - l)
128 x2 = (x << l) & xmax
129 #
130 #####
131 # count number of fermions, i.e. 1's in bit configuration of x1
132 f_count1 = _count_particles_32(x1, period)
133 # count number of fermions, i.e. 1's in bit configuration of x2
134 f_count2 = _count_particles_32(x2, period)
135 #####
136 # compute fermion sign
137 sign_ptr[0] *= -1 if ((f_count1 & 1) & (f_count2 & 1) & 1) else 1
138 #
139 return x2 | x1
140
141
142T_args = np.array([1, (1 << N) - 1], dtype=np.uint32)
143
144
145#
146@cfunc(map_sig_32, locals=dict(out=uint32, s=uint32, f_count=int32))
147def parity(x, N, sign_ptr, args):
148 """works for all system sizes N."""
149 out = 0
150 s = args[0]
151 #
152 #####
153 # count number of fermions, i.e. 1's in bit configuration of the state
154 f_count = _count_particles_32(x, N)
155 #####
156 sign_ptr[0] *= -1 if ((f_count & 2) and 1) else 1
157 #
158 out ^= x & 1
159 x >>= 1
160 while x:
161 out <<= 1
162 out ^= x & 1
163 x >>= 1
164 s -= 1
165 #
166 out <<= s
167 return out
168
169
170P_args = np.array([N - 1], dtype=np.uint32)
171#
172###### construct user_basis
173# define anti-commuting bits -- fermion signs on the integer bits (not sites!) that represent a fermion degree of freedom
174noncommuting_bits = [
175 (np.arange(N), -1)
176] # fermion signs are counted w.r.t. the shift operator <<
177# define maps dict
178maps = dict(
179 T_block=(translation, N, 0, T_args),
180 P_block=(parity, 2, 0, P_args),
181)
182# maps = dict(P_block=(parity,2,0,P_args), )
183# maps = dict(T_block=(translation,N,0,T_args) )
184op_dict = dict(op=op, op_args=op_args)
185# create user basiss
186basis = user_basis(
187 np.uint32,
188 N,
189 op_dict,
190 allowed_ops=set("xy+-nI"),
191 sps=2,
192 noncommuting_bits=noncommuting_bits,
193 **maps,
194)
195#
196#
197print(basis)
198#
199############ create and compare Hamiltonians #############
200#
201##### Hamiltonian in using Majoranas
202#
203J = -np.sqrt(2.0) # hoppping
204U = +1.0 # nn interaction
205#
206hop_term_p = [[+0.5j * J, j, (j + 1) % N] for j in range(N)]
207hop_term_m = [[-0.5j * J, j, (j + 1) % N] for j in range(N)]
208density_term = [[+0.5j * U, j, j] for j in range(N)]
209int_term = [[-0.25 * U, j, j, (j + 1) % N, (j + 1) % N] for j in range(N)]
210id_term = [[0.25 * U, j] for j in range(N)]
211#
212static = [
213 ["xy", hop_term_p],
214 ["yx", hop_term_m], # kinetic energy
215 ["I", id_term],
216 ["xy", density_term],
217 ["xyxy", int_term], # nn interaction energy
218]
219dynamic = []
220#
221no_checks = dict(check_symm=False, check_pcon=False, check_herm=False)
222H_majorana = hamiltonian(static, [], basis=basis, dtype=np.float64, **no_checks)
223#
224#
225##### Hamiltonian using complex fermions
226#
227#
228hopping_pm = [[+J, j, (j + 1) % N] for j in range(N)]
229hopping_mp = [[-J, j, (j + 1) % N] for j in range(N)]
230nn_int = [[U, j, (j + 1) % N] for j in range(N)]
231#
232static = [["+-", hopping_pm], ["-+", hopping_mp], ["nn", nn_int]]
233dynamic = []
234#
235no_checks = dict(check_symm=False, check_pcon=False, check_herm=False)
236H = hamiltonian(static, [], basis=basis, dtype=np.float64, **no_checks)
237print(H.toarray())
238print()
239print(H_majorana.toarray())
240print()
241print(np.linalg.norm((H - H_majorana).toarray()))