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
