I am trying to plot the density and current density of a chiral edge state in a 2D quantum Hall system. I tried to build up the system in two ways:
1. define two sublattices,
2. define one lattice, use matrix form for the on-site potential and hoppings
I expected that these two ways should yield the same result. However, the result for the chiral edge state is not the same. In terms of density, method 1 gives the correct picture whereas method 2 doesn't. In terms of current density, the resulting picture
is the same, whereas len(current) differs by a factor of 2.
Is there anything wrong in one of these methods of building a system or in my method of plotting density and current density? Thank you!
# -*- coding: utf-8 -*-
import kwant
import tinyarray as ta
lat_u = kwant.lattice.square(a=1,name='u',norbs=1)
lat_d = kwant.lattice.square(a=1,name='d',norbs=1)
def make_system():
syst = kwant.Builder()
W=30
x,y=W,W
M=0.5
k_0=1
for i in range(x):
for j in range(y):
syst[lat_u(i,j)] = M*(k_0**2-4)
syst[lat_d(i,j)] = (M*(k_0**2-4))*(-1)
if i>0:
syst[lat_u(i,j),lat_u(i-1,j)] = M
syst[lat_d(i,j),lat_d(i-1,j)] = -M
syst[lat_u(i,j),lat_d(i-1,j)] = 1j/2
syst[lat_d(i,j),lat_u(i-1,j)] = 1j/2
if j>0:
syst[lat_u(i,j),lat_u(i,j-1)] = M
syst[lat_d(i,j),lat_d(i,j-1)] = -M
syst[lat_u(i,j),lat_d(i,j-1)] = 1/2
syst[lat_d(i,j),lat_u(i,j-1)] = -1/2
syst1=syst.finalized()
return syst1
def make_system2():
syst = kwant.Builder()
lat = kwant.lattice.square(a=1,name='u',norbs=2)
sigma_x = ta.array([[0, 1], [1, 0]])
sigma_y = ta.array([[0, -1j], [1j, 0]])
sigma_z = ta.array([[1, 0], [0, -1]])
W=30
x,y=W,W
M=0.5
k0=1
for i in range(x):
for j in range(y):
syst[lat(i,j)] = M*(k0**2-4)*sigma_z
if i>0:
syst[lat(i,j),lat(i-1,j)] = 1j/2*sigma_x+M*sigma_z
if j>0:
syst[lat(i,j),lat(i,j-1)] = 1j/2*sigma_y+M*sigma_z
syst1=syst.finalized()
return syst1
def plot_data(syst,n):
import scipy.linalg as la
ham = syst.hamiltonian_submatrix()
evecs = la.eigh(ham)[1]
psi=evecs[:, n]
J_0 = kwant.operator.Current(syst)
current=J_0(psi)
print(len(current))
kwant.plotter.current(syst, current, colorbar=False)
wf = abs(psi)**2
def site_size(i):
return wf[i] / wf.max()
kwant.plot(syst, site_size=site_size, site_color=(0, 0, 1, 0.3),
hop_lw=0.1)
def main():
syst=make_system()
plot_data(syst,900)
syst2=make_system2()
plot_data(syst2,900)
if __name__ == '__main__':
main()