Dear all, 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!
Bests, Guangze The code is as follows: # -*- 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()