Dear Prof Akhmerov, Dear Dr Istas, Many thanks for your email. We have since then found an older implementation of the bound state algorithm hosted on github<https://gitlab.kwant-project.org/jbweston/kwant/-/blob/boundstate/kwant/phys...>, which follows the original paper more faithfully and can calculate wavefunctions in the leads. The older version has some problems with the root finder, and is limited in that it can only find at most one bound state in a given energy range. However, re-importing the energies calculated from the featured boundstate algorithm show that the two implementations are mostly consistent with each other. However, for our system of interest, which has a BdG Hamiltonian, we have been getting unphysical lead wavefunctions in 3 senses: i) they are discontinuous at the lead-scattering region interface; ii) their amplitude do not vary smoothly as a function of model parameters for the same branch of solution; and iii) their norm do not sum up to 1 despite having been passed through a normalisation routine in the boundstate code. This makes our calculations of dipole matrix elements impossible. We have reduced the model to the simplest form that still retains the errors we have been getting. We have also put in another example where the boundstate code works as a baseline. We wonder if you would be happy to quickly have a look at the problem, described in the attached notebook, and point to where we might find potential bugs or fixes? Many thanks for your help and assistance. Yours, Sincerely Chi & Ryan Imperial College London ________________________________ From: Anton Akhmerov <anton.akhmerov+kd@gmail.com> Sent: 04 February 2022 15:12 To: Zhang, Chi <chi.zhang118@imperial.ac.uk> Cc: kwant-discuss@kwant-project.org <kwant-discuss@kwant-project.org>; Tiew, Hoe R <hoe.tiew18@imperial.ac.uk>; Moors, Kristof <k.moors@fz-juelich.de> Subject: Re: [Kwant] Inquiry about extracting boundstate eigenvectors/eigenvalues from kwant ******************* This email originates from outside Imperial. Do not click on links and attachments unless you recognise the sender. If you trust the sender, add them to your safe senders list https://spam.ic.ac.uk/SpamConsole/Senders.aspx to disable email stamping for this address. ******************* Hello Chi and Ryan, The vecs and vecslmbainv arrays are indeed linear superpositions of the lead modes, which means that the evanescent parts are all mixed. vecs and vecslmbdainv are therefore related by multiplication of an upper-triangular part of the Schur decomposition. Or rather it's a bit more complicated—the propagating modes are diagonalized, the evanescent modes are kept orthogonal (and are therefore superpositions of different translation eigenvectors). Therefore to extract the wave functions in further unit cells you'll need to dig into the algorithm. In these lines you see the Schur form being reordered to isolate the evanescent modes. You'll want to also keep the Schur form of the eigenproblem. The wave function is given by the Schur vectors times the Schur matrix to the power of the number of the unit cell. Regarding the part with normalization by hopping—that's another tricky bit if your hopping isn't full rank. In general, we project the wave function on the nonzero singular values of the hopping, and define the extract function to obtain the full wave function (see https://gitlab.kwant-project.org/kwant/kwant/-/blob/v1.4.2/kwant/physics/lea...). Does this give you enough pointers in the right direction? Best, Anton On Fri, 4 Feb 2022 at 11:42, Zhang, Chi <chi.zhang118@imperial.ac.uk> wrote:

Dear Devs,

I hope this email finds you well.

We are two Master's students trying to calculate Andreev/Majorana transition elements using the kwant boundstate algorithm (doi: 10.21468/SciPostPhys.4.5.026). For our purposes we need to extract the evanescent tails of the boundstates in the leads as well, which we had hoped to obtain from the StabilisedModes object in kwant.InfiniteSystem.modes(). However, in general, the vecs and vecslmbdainv object are not related by a constant eigenvalue, as advertised in the documentation. We wonder if you could tell us how to transform from the vecs and vecslmbdainv objects back to the eigenvectors*eigenvalues and eigenvectors object for our calculations?

Many thanks in advance for your response!

Sincerely, Chi Zhang & Ryan Tiew MSci Physics with Theoretical Physics Imperial College London

----------------------------------------------------- Some technical details of what we have tried here. Please ignore this if you don't have enough time to read.

We are aware that kwant is transforming vecslmbdainv by sqrt(norm(hop))*eigenvectors, and vecs by (hop/sqrt(norm(hop)))@eigenvectors@eigenvalues. However, this does not seem to be the full story and we are not able to reproduce what is stored in vecs and vecslmbdainv when there are more than 1 evanescent modes. In particular, we have tried to do a real generalised Schur decomposition of the two matrices in the lead schroedinger's equation, grabbed the vectors stored in the right unitary, and put them through the basis transformation described above. Still the there is no sign of better agreements. We therefore we believe it would be more productive to inquire directly from you how to transform from vecs and vecslmbdainv to actual eigenvectors and eigenvalues. Many thanks in advance for your reply.