HI Ville,
I'm trying to use Kwant to investigate behavior of a lattice between two superconductors, namely trying to figure out the Josephson current. The task is basically to diagonalize the related Bogoliubov-de Gennes Hamiltonian and use the eigenvectors and eigenvalues to calculate the current in the lattice. I basically follow the example 2.6. in Kwant documentation to build up the Hamiltonian but let the both leads be superconductors and establish the system in between. I introduce the phase difference between the leads by multiplying the Delta parameter of the other lead by e^(i phi). However, I'm a bit stuck there as I don't know how to utilize the obtained eigenvectors to make wave functions which I could input to t current operator in Kwant.
To get the correct Josephson current you'll need to include the contribution from the Andreev bound states that form in the Junction. As these are true bound states you won't be able to get them by using 'kwant.wave_function', as the latter only calculates the *scattering* states. Some of the Kwant Authors recently published a paper that describes how to calculate bound states for general scattering systems [1] (i.e. systems with leads). We're working to implement this in Kwant, but it's not there yet. As far as I am aware the following are possible approaches to calculate the bound state spectrum are: + truncate the leads after some (large) length and diagonalize the resulting system. You'll have to identify "by eye" which states correspond to the true bound states of the infinite system you are approximating. The leads will need to be truncated far enough away from the scattering region that the bound state wavefunction is nearly vanishing. + Replace the leads with a self-energy term and diagonalize the resulting matrix. I can't remember what guarantees one gets about using this approach, as you make your Hamiltonian non-Hermitian by doing this. + Construct your system as 3 parts: an SN interface, an N-N interface (include any disorder or whatever here), and an NS interface. Calculate the scattering matrices of the 3 parts separately, and then the bound states can be found by combining the scattering matrices in a particular way and solving for zero determinant (I can't remember the details, but it's in section 2 of [2], and references therein). Basically everyone in our group in Delft calculates Josephson currents as a matter of course, so they're probably in a better position to tell you which way to go. I thought I'd give my 2 cents anyway. Happy Kwanting! Joe [1]: https://arxiv.org/pdf/1711.08250.pdf [2]: https://arxiv.org/pdf/cond-mat/0406127.pdf