Dear Chen, The problem is in the way you plot the results: syst1 has twice the number of the sites in syst2 (you can check this by len(list(syst.sites))). For syst1, you plot all the wave function but you are plotting psi_up and psi_down at the same point. (2n sites (up and down at the same position) and 2n components). Be careful, in this case, you are not plotting the sum of the density up and down but you are plotting one on top of the other. For syst2, you have n sites and 2n components but you are plotting only the n first components of the wavevector. If you want the total density you need to sum the two components for different spins: wave[::2]**2+wave[1::2]**2 I suggest for you to write two different functions for the two plots. I hope this helps Adel On Sat, Nov 17, 2018 at 11:58 AM Chen Guangze <guchen@student.ethz.ch> wrote:

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()

-- Abbout Adel