Hi,

I still don’t understand exactly how I can implement this. From what I can understand you choose the energy not the eigenstate for the density function.

I am unclear what you mean here. In your code you do the following:

ham_mat = sys.hamiltonian_submatrix(sparse=True, params = params)

evals, evecs = sorted_eigs(sla.eigsh(ham_mat.tocsc(), k=20,
sigma=0))

kwant.plotter.map(sys, np.abs(evecs[:, 9])**2

You take the hamiltonian of the scattering region, diagonalize it, and then plot *a single eigenvector*. Your problem was that because you have >1 degree of freedom per site you cannot simply plot the eigenvector as-is. You need to first compute the density per site for this eigenvector, which you can do using the operator module:

rho = kwant.operator.Density(np.kron(-s_z, s_0)) # holes count +1 to the density, and electrons -1

kwant.plotter.map(sys, rho(evecs[:, 9))

Note that you call the density operator with a single eigenvector.

Independently of the above it's not clear to me that you should be diagonalizing just the scattering region Hamiltonian, given that your system has leads. It seems to me that you should be calculating the scattering states and the density of those.

Happy Kwanting,

Joe