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