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.