Hi Zhan,
Thanks very much for your reply. Looking forward to the details in your review article soon after. I wonder how kwant.physics.modes obtain propagating modes, which can be used to calculate transport observables.
As known, let \lambda=exp^{ik}, propagating modes have |\lambda|=1 while the others have |\lambda|\ne 1. One can solve a (generalized) eigenproblem to obtain all \lambda and idenity propagating modes as those having |\lambda| equal to 1 within a small tolerance. If so,for an infinitely long wire with a square cross section possessing N*N lattice site, one has to solve all the eigenvalues of a matrix whose size is no smaller than N^2 by N^2, which is a very difficult numerical task for a large N (say N=1000). Does Kwant solve all \lambda as I describe or in a fast way that directly single out the propogating modes?
Kwant uses dense linear algebra for the lead modes calculation, so if you are trying to solve the modes problem for a system with a 1M orbital cross-section in the unit cell, you will indeed not even be able to store the Hamiltonian in memory (it would require ~16TB) let alone do any calculations with it.
Kwant uses a (generalized) Schur decomposition to solve the modes problem, and then extracts the propagating eigenvectors, and a basis of Schur vectors (guaranteed to be orthogonal) for the decaying subspace.
I am not aware of any algorithms for sparse (generalized) Schur decomposition, so it is not clear how we could generalize this part of Kwant to use sparse matrices.
My understanding is that it is done this way for stability reasons (which I will try and outline here).
In order to solve the scattering problem we need both the
propagating *and* decaying modes. Even though we need the
propagating modes (eigenvectors) directly (to separate incoming
and outgoing modes, and interpret the scattering amplitudes) we do
*not* need the decaying eigenvectors themselves (because we do not
typically use the "scattering" amplitudes into the decaying modes
for transport calculations); we need only an orthogonal basis for
the decaying subspace. In addition the decaying modes are often
almost collinear (e.g. in graphene), so we would run into
stability problems when solving the scattering problem if we used
the decaying modes directly.
Even if you didn't use the Schur decomposition, and calculated
the eigenvectors directly, you'd still have to compute all the
propagating and decaying eigenvectors, and sparse iterative
methods are typically best at finding a few eigenvectors only. I
suspect that the best way to proceed would be highly tailored to
the specific system that you're studying.
Happy Kwanting,
Joe