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

\[c^x = c + c^\dagger,\qquad c^y = -i(c - c^\dagger),\]

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

\[c = \frac{1}{2}\left( c^x + i c^y \right), \qquad c^\dagger = \frac{1}{2}\left( c^x - i c^y \right).\]

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 – A tutorial on QuSpin’s user_basis – 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:

\[H = \sum_{j=0}^{L-1} J\left(c^\dagger_{j} c_{j+1} + \mathrm{h.c.}\right) + U n_{j}n_{j+1},\]

where \(J\) is the hopping matrix element, and \(U\) is the nearest-neighbor interaction strength.

In the Majorana representation, the same Hamiltonian reads as

\[H = \sum_{j=0}^{L-1} \frac{iJ}{2}\left( c^x_j c^y_{j+1} - c^y_j c^x_{j+1} \right) + \frac{U}{4}\left(1 + 2ic^x_j c^y_j - c^x_j c^y_j c^x_{j+1} c^y_{j+1} \right).\]

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

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
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
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 24                                      #
# This example shows how to use the `user_basis` to define Majorana operators. #
################################################################################
from quspin.operators import hamiltonian # Hamiltonians and operators
from quspin.basis import spinless_fermion_basis_1d # Hilbert space fermion basis_1d
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
from scipy.special import comb
np.set_printoptions(suppress='True', precision=6)
#
N=6 # lattice sites
#
############   create soinless fermion user basis object   #############
#
@jit(uint32(uint32,uint32),locals=dict(f_count=uint32,),nopython=True,nogil=True)
def _count_particles_32(state,site_ind):
	# auxiliary function to count number of fermions, i.e. 1's in bit configuration of the state, up to site site_ind
	# CAUTION: 32-bit integers code only!
	f_count = state & ((0x7FFFFFFF) >> (31 - site_ind));
	f_count = f_count - ((f_count >> 1) & 0x55555555);
	f_count = (f_count & 0x33333333) + ((f_count >> 2) & 0x33333333);
	return (((f_count + (f_count >> 4)) & 0x0F0F0F0F) * 0x01010101) >> 24
#
@cfunc(op_sig_32,
	locals=dict(s=int32,sign=int32,n=int32,b=uint32,f_count=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.
	#####
	f_count = _count_particles_32(op_struct.state,site_ind)
	#####
	sign = -1 if f_count&1 else 1
	n = (op_struct.state>>site_ind)&1 # either 0 or 1
	b = (1<<site_ind)
	#
	if op_str==120: # "x" is integer value 120 = ord("x")
		op_struct.state ^= b
		op_struct.matrix_ele *= sign

	elif op_str==121: # "y" is integer value 120 = ord("y")
		op_struct.state ^= b
		op_struct.matrix_ele *= -1.0j*sign*((n<<1)-1)
	
	elif op_str==43: # "+" is integer value 43 = ord("+")
		op_struct.matrix_ele *= (0.0 if n else sign)
		op_struct.state ^= b # create fermion

	elif op_str==45: # "-" is integer value 45 = ord("-")
		op_struct.matrix_ele *= (sign if n else 0.0)
		op_struct.state ^= b # create fermion
		
	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
#
@cfunc(map_sig_32,
	locals=dict(shift=uint32,xmax=uint32,x1=uint32,x2=uint32,period=int32,l=int32,f_count1=int32,f_count2=int32) )
def translation(x,N,sign_ptr,args):
	""" works for all system sizes N. """
	shift = args[0] # translate state by shift sites
	period = N # periodicity/cyclicity of translation
	xmax = args[1]
	#
	l = (shift+period)%period
	x1 = (x >> (period - l))
	x2 = ((x << l) & xmax)
	#
	#####
	# count number of fermions, i.e. 1's in bit configuration of x1
	f_count1 = _count_particles_32(x1,period)
	# count number of fermions, i.e. 1's in bit configuration of x2
	f_count2 = _count_particles_32(x2,period)
	#####
	# compute fermion sign
	sign_ptr[0] *= (-1 if ((f_count1&1)&(f_count2&1)&1) else 1)
	#
	return (x2 | x1)
T_args=np.array([1,(1<<N)-1],dtype=np.uint32)
#
@cfunc(map_sig_32,
	locals=dict(out=uint32,s=uint32,f_count=int32) )
def parity(x,N,sign_ptr,args):
	""" works for all system sizes N. """
	out = 0
	s = args[0]
	#
	#####
	# count number of fermions, i.e. 1's in bit configuration of the state
	f_count = _count_particles_32(x,N)
	#####
	sign_ptr[0] *= (-1 if ((f_count&2) and 1) else 1)
	#
	out ^= (x&1)
	x >>= 1
	while(x):
		out <<= 1
		out ^= (x&1)
		x >>= 1
		s -= 1
	#
	out <<= s
	return out
P_args=np.array([N-1],dtype=np.uint32)
#
######  construct user_basis
# define anti-commuting bits -- fermion signs on the integer bits (not sites!) that represent a fermion degree of freedom
noncommuting_bits = [(np.arange(N),-1)] # fermion signs are counted w.r.t. the shift operator << 
# define maps dict
maps = dict(T_block=(translation,N,0,T_args), P_block=(parity,2,0,P_args), ) 
#maps = dict(P_block=(parity,2,0,P_args), )
#maps = dict(T_block=(translation,N,0,T_args) ) 
op_dict = dict(op=op,op_args=op_args)
# create user basiss
basis = user_basis(np.uint32,N,op_dict,allowed_ops=set("xy+-nI"),sps=2,noncommuting_bits=noncommuting_bits,**maps)
#
#
print(basis)
#
############   create and compare Hamiltonians   #############
#
##### Hamiltonian in using Majoranas 
#
J=-np.sqrt(2.0) # hoppping
U=+1.0 # nn interaction
#
hop_term_p=[[+0.5j*J,j,(j+1)%N] for j in range(N)]
hop_term_m=[[-0.5j*J,j,(j+1)%N] for j in range(N)]
density_term=[[+0.5j*U,j,j] for j in range(N)]
int_term=[[-0.25*U,j,j,(j+1)%N,(j+1)%N] for j in range(N)]
id_term=[[0.25*U,j] for j in range(N)]
#
static=[['xy',hop_term_p],['yx',hop_term_m], 					# kinetic energy
		['I',id_term],['xy',density_term],['xyxy',int_term],	# nn interaction energy
		]
dynamic=[]
#
no_checks=dict(check_symm=False, check_pcon=False, check_herm=False)
H_majorana=hamiltonian(static,[],basis=basis,dtype=np.float64,**no_checks)
#
#
##### Hamiltonian using complex fermions
#
#
hopping_pm=[[+J,j,(j+1)%N] for j in range(N)]
hopping_mp=[[-J,j,(j+1)%N] for j in range(N)]
nn_int=[[U,j,(j+1)%N] for j in range(N)]
#
static=[["+-",hopping_pm],["-+",hopping_mp],["nn",nn_int]]
dynamic=[]
#
no_checks=dict(check_symm=False, check_pcon=False, check_herm=False)
H=hamiltonian(static,[],basis=basis,dtype=np.float64,**no_checks)
print(H.toarray())
print()
print(H_majorana.toarray())
print()
print(np.linalg.norm((H-H_majorana).toarray()))