Dear All, Has anyone tried to calculate bound states in a JJ? I used the package by Benoit Rossignol available at: https://gitlab.kwant-project.org/kwant/boundstate . I am sending attached the code I used to simulate it but I am having some difficulties interpreting the results. First, the algorithm could not find bound states, i.e., psi_tot returned by the solver is an empty array, for N<5. Second, it returns a total wavefunction for larger N, however the wavefunction oscillates between each site and decays very slowly in the leads. In the original proposal of this method, Istas et al (arXiv:1711.08250) solves the wavefunction of quantum billiard (Fig 5) in which this does not seem to happen. Has anyone found something similar? I do know if it a mistake in the code, an attribute of discretization or something else. Scanning over N=[1,2,3,5,10,20] and ,u=[2, Delta, 0, -2] the algorithm returns the following non-zero wavefunctions: SNS junction with N=10, S=20 and mu=2 [image: N_10_S_20_mu_2.png] SNS junction with N=20, S=20 and mu=2 [image: N_20_S_20_mu_2.png] SNS junction with N=5, S=20 and mu=-2 (in this case the amplitude of the bound states decays to zero in the lead) [image: N_5_S_20_mu_-2.png] SNSNS junction with 20/10/10/10/20 sites and mu=2. The wavefunction in the N region is similar to the SNS junction but there seems to be no decay in the middle S region. [image: N_10_S_20_mu_2.png] and of course it's mirrored version: [image: NSN_10_mu_2.png] I appreciate any comments on this subject. Best regards, Denise
