Hi Joe, Thank you, now it works!
But I am a bit confused about the physical meaning of the current, created by modes in lead n. How to plot the real current between the two leads? Imagine, as the simplest case, that my geometry is inversion-symmetric, so that lead 0 goes into lead 1 under inversion. Then I expect the current to be inversion-symmetric as well.
Summing the current contributions for *all* the scattering states originating from lead L at an energy E gives you the current you would measure if the fermi levels in all the leads are at energy E and you apply an infinitesimal voltage to lead L.
For example if we want to see the current profile if we add an infinitesimal voltage dV to lead 0:
wfs = kwant.wave_function(syst, energy=E, params=...) J_0 = dV * sum(J(psi) for psi in wfs(0))
This is analogous to the transmission obtained by kwant.smatrix:
smatrix = kwant.smatrix(syst, energy=E, params=...) I_10 = dV * smatrix.transmission(1, 0)
I_10 is the current we would measure in lead 1 after applying dV to lead 0.
In the above I have elided the e^2/h factor for brevity, but hopefully I have been clear enough to get the point across.
Hope that helps,