Dear Kwant developers and users, I write concerning a doubt that I encountered trying to study by Kwant a model for the quantum Hall/quantum spin Hall effect. In particular, I am focusing on the two-dimensionakl BHZ model, first proposed in the famous paper https://science.sciencemag.org/content/314/5806/1757 Notice that first I keep only one of the two time-reversal-connected blocks forming the quantum spin Hall system (when the mass parameter M is positive); it is known that each block describes a quantum Hall system (again for M >0, otherwise a trivial insulator occurs). I attach below my defining code (included in the % lines) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% import kwant from kwant.digest import uniform # a (deterministic) pseudorandom number generator import math from matplotlib import pyplot as plt import cmath import scipy.sparse.linalg as sla import scipy.linalg as la import numpy from numpy import linalg as LA # For plotting from matplotlib import pyplot # For matrix support import tinyarray #import progressbar #----------------------------------------------------- # We define the variables and the constants a = 1 #lattice step #t = 1 #along x and y #----------------------------------------------------- nullmatrix = tinyarray.array([[0, 0], [0, 0]]) matid = tinyarray.array([[1, 0], [0, 1]]) sigma_0 = tinyarray.array([[1, 0], [0, -1]]) sigma_xp = tinyarray.array([[1, -1j], [-1j, -1]]) sigma_yp = tinyarray.array([[1, 1/2], [-1/2, -1]]) ----------------------------------------------------------------------------------------------------------------- def sistema_BHZ_2e(Lx, Ly, M) : syst = kwant.Builder() lat = kwant.lattice.square(norbs=2) for i in range(0,Lx): #syst[lat(i, 0), lat(i, L - 1)] = sigma_yp for j in range(0,Ly): # On-site Hamiltonian syst[lat(i, j)] = (M - 4)*sigma_0 # Hopping in y-direction if i > 0: syst[lat(i, j), lat(i - 1, j)] = sigma_xp # Hopping in z-direction if j > 0: syst[lat(i, j), lat(i, j - 1)] = sigma_yp return syst %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Correspondingly to what described above, I would expect to see, when M

0, edge modes. Indeed, I obtain them calculating the local densities via the KPM (approximate) method working in Kwant:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% # Plot several local density of states maps in different subplots def plot_ldos(syst, densities, file_name=None): step = 0.1 fig, axes = plt.subplots(1, len(densities), figsize=(7*len(densities), 7)) for ax, (title, rho) in zip(axes, densities): kwant.plotter.density(syst,rho.real, colorbar=True, ax=ax, relwidth = step) ax.set_title(title) ax.set(adjustable='box', aspect='equal') plt.show() plt.clf() -------------------------------------------------------------------------------------------------------------- M = 2 L = 16 Lx = L Ly = L syst = sistema_BHZ_2e(Lx, Ly, M) syst = syst.finalized() hamat = syst.hamiltonian_submatrix(sparse = False) evals, evecs = la.eigh(hamat) kwant_op = kwant.operator.Density(syst, sum=False) local_dos = kwant.kpm.SpectralDensity(syst, operator=kwant_op) print("first negative energy =", evals[int(Lx*Ly-1)]) print("first_positive_energy =", evals[int(Lx*Ly)]) first_negative_energy_ldos = local_dos(energy=evals[int(Lx*Ly-1)]) first_positive_energy_ldos = local_dos(energy=evals[int(Lx*Ly)]) plot_ldos(syst, [('first negative energy', first_negative_energy_ldos),('first positive energy', first_positive_energy_ldos)]) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% I attach below a typical output (KPM.pdf). However, when I calculate the same densities not using the KPM approach (then via the exact function kwant.operator.Density(syst) only ) , that means by the code: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% M = 2 L = 16 Lx = L Ly = L syst = sistema_BHZ_2e(Lx, Ly, M) syst = syst.finalized() hamat = syst.hamiltonian_submatrix(sparse = False) evals, evecs = la.eigh(hamat) kwant_op = kwant.operator.Density(syst, sum=False) print("first negative energy =", evals[int(Lx*Ly-1)]) print("first_positive_energy =", evals[int(Lx*Ly)]) first_negative_energy_ldos = kwant_op(evecs[int(Lx*Ly-1)]) first_positive_energy_ldos = kwant_op(evecs[int(Lx*Ly)]) plot_ldos(syst, [('first negative energy', first_negative_energy_ldos),('first positive energy', first_positive_energy_ldos)]) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% I obtain a totally different density profile, not displaying edge localization (see No_KPM.pdf attached below). I do not understand the reason of such a mismatch. Can you help me ? Thank you very much and best, L.