Dear Luca,
I forgot to add the kwant-discuss email for my first answer.
First, we shall shortly describe the main characteristics of KPM and some specifics of its implementation in general (not in kwant so I might get wrong). In fact, solving a tight-binding problem implies the diagonalization of the Hamiltonian matrix. In this approach we mainly avoid diagonalization and instead use approximative methods to calculate the properties of interest.
Within KPM approach is often expressed over an energy broadening which lead to obtain spatially dependent properties ( let us say local density of states or Green’s function) which are computed separately for each spatial position. In this purpose the KPM is seen as orthogonal to conventional eigenvalue solvers or simply say sparse diagonalization problem.
Comming back to you question:
I guess that the kwant.kpm.SpectralDensity need to have the energy broadening arguments of your system. So the KPM computes the entire spectrum simultaneously (compute the spatialLDOS at multiple energy values in one calculation ) for this purpose we have to set like you did
kwant_op = kwant.operator.Density(syst, sum=False)
local_dos = kwant.kpm.SpectralDensity(syst, operator=kwant_op)
It was stated in kwant. Pdf that " kpm.SpectralDensity calculates the spectral density of an operator This class makes use of the kernel polynomial method (KPM), presented in [R1], to obtain the spectral density ρ A (e), as a function of the energy e, of some operator A that acts on a kwant system or a Hamiltonian.”
If not may be in your case,in kwant.operator.Density you have to normalize over the sum of all states.
For all these statements and based on my understanding, I think from physical point of view the kwant.operator.Density and kwant.kpm.SpectralDensity are not the same and the your result would confirm my thought. So, I would say that your question relies on the key difference between the notion of Density and spectral density. You can check the class of kwant.operator to figure out the difference.
Dear Adel,
thank you very much for your mail and sorry for the
length of my previous one.
I found very clear and helpful the link that you sent to me,
allowing me to perform various checks.
In particular, I checked the correctness of my Hamiltonian.
Beyond the doubts about the quantum spin Hall effect and its
Hamiltonian toy model,
I remain with the trouble to compute edge states.
Actually, I do not get the difference between doing it via the KPM algorithm
(now let me skip the details about the definition of the Hamiltonian)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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)
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)])
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
yielding nice edge states in the topological phase at the mass parameter M < 0
(I checked that the energie are the correct ones for the edge states,
around zero),
and the (exact) calculation of the spatial density
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
..................
kwant_op = kwant.operator.Density(syst, sum=False)
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)])
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
yielding states (with the same energies of those above)
not localized on the edges, as they should.
Do you have some hint about the mismatch ?
Thank you very much again
L.
Il giorno sab 29 feb 2020 alle ore 19:48 Adel Belayadi
<adelphys@gmail.com> ha scritto:
>
> Dear Luca,
> Good day.
> I guess the link below would help you out.
> https://gitlab.kwant-project.org/qt/topocm/blob/7a2b68a616ef2d5c75222162f77abd793a291ae3/src/w5_qshe/qshe_experiments.md
> best, Adel
>
> Le ven. 28 févr. 2020 à 09:02, Luca Lepori <llepori81@gmail.com> a écrit :
>>
>> 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.