Hi Adel,


Thanks for the very astute observations! I admit this had me worried for a little bit, but I'm now confident that I understand what is going on. While this is not a bug in Kwant, it is confusing and the tutorial is not showing what it is meant to, so thank you for pointing this out



I was consulting the example of the current in a closed system and I am a little bit confused:


The figure showing the current flowing inside a circle [1] is obtained for a magnetic field B=0 !

But we know that for B=0, the Hamiltonian is a real symmetric matrix and thus diagonalized by a real orthogonal matrix, so the eigenvectors are real (any phase is just artificial.)

By using the formula giving the current we should find zero !
moreover we can ask why we have a current flowing clockwise and not anticlockwise (I remind that B=0) !

If I look at the spectrum at B=0 I see that some eigenvalues are doubly degenerate. So if ψ_a and ψ_b are real energy eigenvectors with energy E, then any linear combination is also an eigenvector with the same energy. Even if ψ_a and ψ_b are real (up to a global phase), and hence carry no current, there is no such restriction on ψ_c. Unfortunately when there is a degeneracy the eigenvalue solver will in general pick an arbitrary basis in this 2D subspace, which will mean that the eigenvectors you have will not necessarily respect time reversal symmetry, even if the Hamiltonian does. I verified that for eigenvectors with no degenerate partner @ B=0 the current is indeed 0 everywhere (as it should be!). This is *not*, however, what we wanted to show in the tutorial!



When we put B=0.05,  we obtain  the figures below for the current and the wave function.
the wave function has circular symmetry but not the current. Moreover, x and y play exactly the same role so the current should be the same on x and on y (independently from the used gauge)

Am I missing something ?

This is unfortunately a pitfall of using the operators as they are currently: when evaluating an operator on a wavefunction, you *must* pass the same parameters as you did when getting the Hamiltonian that produced that wavefunction:

    J = kwant.operator.Current(syst)

    H = syst.hamiltonian_submatrix(sparse=True, args=[B])
    _, evecs = scipy.sparse.linalg.eigsh(H, ...)
    current = J(evecs[:, 9], args=[B])  # note the same 'args' are passed here too!

If you do the following:

    J = kwant.operator.Current(syst)

    H = syst.hamiltonian_submatrix(sparse=True, args=[B])
    _, evecs = scipy.sparse.linalg.eigsh(H, ...)
    current = J(evecs[:, 9])  # forgot to pass the 'args'

then you will get the incorrect result. This is because the operator 'J' needs to evaluate the Hamiltonian in order to evaluate the current, and there is no way for the operator to "know" that 'evecs' was created from a hamiltonian evaluated with magnetic field 'B'!

There is a warning about this in the documentation [1], but it's not obvious from the documentation that the same problem can arise in other contexts.

In the future we want to make it easier to work with these different quantities without having to put the burden of bookkeeping on the user, but it's difficult to do this while trying to keep things simple. If anybody has any suggestions we are very open to them!

Happy Kwanting,

Joe


[1]: https://kwant-project.org/doc/1/tutorial/operators#using-bind-for-speed