Hello all,
I'm trying to implement the BHZ model already in the tight-biding model
according to the paper [Phys. Rev. B 96, 155302 (2017) or arXiv:1707.03624]. I
use the equations (5),(6) and (7), however with no magnetic field. I'm using
Kwant's documentation example 2.10 as a guideline (For what I have to do next
it's more convenient that I don't use the discretize method). After some
struggling I got the band structure of both mine and Kwant's example to be the
same, however my spin density is reversed, it's up where it's supposed to be
down. So far I haven't figured out why.
Also, the main difference between the implementation is only the make_syst
part, all the plotting in my code are snippets from Kwant's discretize.py; My
code is below:
import kwant
import matplotlib.pyplot as plt
import numpy as np
import random as rd
from math import *
import matplotlib as mpl
import tinyarray
sigma_0 = tinyarray.array([[1, 0], [0, 1]])
sigma_x = tinyarray.array([[0, 1], [1, 0]])
sigma_y = tinyarray.array([[0, -1j], [1j, 0]])
sigma_z = tinyarray.array([[1, 0], [0, -1]])
Id = np.eye(4)
G1 = np.kron(sigma_z,sigma_x)
G2 = -np.kron(sigma_0,sigma_y)
G5 = np.kron(sigma_0,sigma_z)
def make_syst(a=20,A=3.65, B=-68.6, D=-51.1, M=-0.01, C=0,L=200,W=200):
lat = kwant.lattice.square(a,norbs=4)
syst = kwant.Builder()
def square(pos):
(x, y) = pos
return (0 <= y < W and 0 <= x < L)
def onsite(site):
(x, y) = site.pos
return C*Id+M*G5-(4/a**2)*(D*Id+B*G5) #equation 5
def hoppingx(site1,site2):
return (1/a**2)*(D*Id+B*G5)+(1j/(2*a))*(A*G1) #eq 6
def hoppingy(site1,site2):
return (1/a**2)*(D*Id+B*G5)-(1j/(2*a))*(A*G2) #eq 7
syst[lat.shape(square,(0,0))] = onsite
syst[kwant.builder.HoppingKind((1, 0), lat, lat)] = hoppingx
syst[kwant.builder.HoppingKind((0,1), lat, lat)] = hoppingy
def lead_shape(pos):
(x, y) = pos
return (0 <= y < W)
#### Define the left lead. ####
lead = kwant.Builder(kwant.TranslationalSymmetry((-a, 0)))
lead[lat.shape(lead_shape,(0,0))] = onsite
# hoppings in x-direction
lead[kwant.builder.HoppingKind((1, 0), lat, lat)] =hoppingx
# hoppings in y-directions
lead[kwant.builder.HoppingKind((0, 1), lat, lat)] = hoppingy
return syst,[lead,lead.reversed()]
def main():
syst,leads = make_syst(W=1000,L=2000)
# Attach the leads to the system.
for lead in leads:
syst.attach_lead(lead)
# Finalize the system.
syst = syst.finalized()
# Check that the system looks as intended.
#kwant.plot(syst)
# Compute the band structure of lead 0.
kwant.plotter.bands(syst.leads[0],
momenta=np.linspace(-0.3, 0.3, 201), show=False)
plt.grid()
plt.xlim(-.3, 0.3)
plt.ylim(-0.05, 0.05)
plt.xlabel('momentum [1/A]')
plt.ylabel('energy [eV]')
plt.show()
wf = kwant.wave_function(syst, energy=0)
# prepare density operators
sigma_z = np.array([[1, 0], [0, -1]])
prob_density = kwant.operator.Density(syst, np.eye(4))
spin_density = kwant.operator.Density(syst, np.kron(sigma_z, np.eye(2)))
# calculate expectation values and plot them
wf_sqr = sum(prob_density(psi) for psi in wf(0)) # states from left lead
rho_sz = sum(spin_density(psi) for psi in wf(0)) # states from left lead
fig, (ax1, ax2) = plt.subplots(1, 2, sharey=True, figsize=(16, 4))
kwant.plotter.map(syst, wf_sqr, ax=ax1)
kwant.plotter.map(syst, rho_sz, ax=ax2)
ax = ax1
im = [obj for obj in ax.get_children()
if isinstance(obj, mpl.image.AxesImage)][0]
fig.colorbar(im, ax=ax)
ax = ax2
im = [obj for obj in ax.get_children()
if isinstance(obj, mpl.image.AxesImage)][0]
fig.colorbar(im, ax=ax)
ax1.set_title('Probability density')
ax2.set_title('Spin density')
plt.show()
# Call the main function if the script gets executed (as opposed to imported).
# See <http://docs.python.org/library/__main__.html>.
if __name__ == '__main__':
main()
Regards,
Victor Regis
