Right now crossing resolution in modes accounts for numerical precision in an inconsistent manner, which may result in false identification of propagating modes as evanescent or vice versa. This is in particular important when the user works in physical units, and absolute values of matrix elements are large.

Example reproducing the bug (in a notebook, comment out the first line for running from a script):

%matplotlib inline
from matplotlib import pyplot
import kwant
import tinyarray as ta

W = 10
mu = 1
scale = 100000
sz = ta.array([[1, 0], [0, -1]])
syst = kwant.Builder(kwant.TranslationalSymmetry((1, 0)))
lat = kwant.lattice.square()
syst[(lat(0, j) for j in range(10))] = scale * mu * sz
syst[lat.neighbors()] = scale * sz
syst = syst.finalized()

kwant.plotter.bands(syst, momenta=401)
len(syst.modes()[0].momenta)  # prints "0"

The system here is designed to have a lot of crossings at zero energy. Check that setting scale = 1 results in modes finding 28 propagating modes, while scale shouldn't change the number of propagating modes.

@michaelwimmer any ideas for workaround? Rescale h and t inside the modes function?

PS credit for discovering goes to @basnijholt