Modeling SNS junctions and Andreev bound states using two impurities and the T -matrix formalism

We recover here the formation of Andreev bound states in SNS junctions by modeling an NS junction as a delta-function"Andreev"impurity, i.e., a localized potential which scatters an electron into a hole with opposite spin. We show using the scattering matrix formalism that, quite surprisingly, an"Andreev"impurity is equivalent to an NS junction characterized by both Andreev reflection and a finite amount of normal scattering. Our formalism cannot capture ideal NS junctions, but only imperfect ones, however, we show that this is enough to capture the generic features of Andreev bound states. Thus, we provide a new analytical tool to describe the formation of ABS using the T-matrix formalism, valid for imperfect junctions which can in general be described only by numerical tools.


I. INTRODUCTION
The complete analytical description of the formation of Andreev bound states (ABS) [1][2][3][4][5] in long SNS junctions [6][7][8][9] usually raises many technical challenges if one is interested in more than the simplest approximations. The most basic solution is using numerical tightbinding models, but if one is interested in analytical features, in general one needs to resort to many approximations and assumptions such as the Bogoliubov-de Gennes equations, 10-12 the Andreev approximation, 13 as well as other approaches. [14][15][16][17][18][19][20][21] Indeed, such limits capture well the most basic features of the ABS but remain far from describing the full behavior of these states as captured by the numerical tools (see, e.g., Ref. [22]), for instance, their bell-shaped dependence on gate-voltage. [22][23][24] Moreover, no analytical tool can capture well imperfect boundaries between the normal and the superconductor (SC), e.g., a finite amount of normal scattering at the leads, which is expected to inevitably arise in typical experimental setups.
Here we show that we can model an NS junction by considering an "Andreev" type impurity: a deltafunction localized potential that scatters an electron into a hole with opposite spin. This equivalence is demonstrated using the scattering formalism: we find that there are regimes of parameters in which the values of the reflection and Andreev reflection coefficients in the NS junction can be the same as those generated by an "Andreev" impurity. In this regime an "Andreev" impurity can thus be a good model for an NS junction characterized by both Andreev reflection and a finite amount of normal scattering.
We subsequently model an SNS junction by two "Andreev" impurities situated at a given distance, as illustrated in Fig. 1. We use the T -matrix formalism which provides an exact solution for the two-impurity problem. The resulting Green's function in the region between the two impurities corresponds to that of a normal region in an SNS junction, and thus it allows to describe the formation of ABS. We calculate the resulting dependence of the ABS energies on various parameters such as the gate voltage and the phase difference between the two SCs and we show that it is consistent with previous findings (see, for instance, the results obtained by exact diagonalization in Ref. [22]).
In Sec. II we provide the equivalence conditions between the NS junction and the "Andreev" impurity. In Sec. III we give the T -matrix formalism necessary to compute the Green's function and subsequently the local density of states (LDOS) in the presence of two impurities. We present our results in Sec. IV, leaving the conclusions to Sec. V.

II. SCATTERING FORMALISM AND CONDITIONS OF EQUIVALENCE BETWEEN THE NS JUNCTION AND THE "ANDREEV" IMPURITY
For an NS junction the BTK theory 8 indicates that the reflection and Andreev reflection coefficients, resulting from injecting an electron from the normal side of the arXiv:1906.10139v1 [cond-mat.mes-hall] 24 Jun 2019 junction (see Fig. 2 for the schematics of the junction), are given by Figure 2. Schematics of an NS junction where we consider E to be inferior to the SC gap, i.e., the transmission coefficients t, tA = 0.
where ∆ is the SC gap, Z is the dimensionless barrier strength introducing a finite amount of normal reflection at the interface. We will focus only on energies much smaller than the gap, thus we have For Z = 0 we have a perfect Andreev junction, i.e., r A = 1, while for Z 1 we have a bad junction with a lot of normal reflection. In order to verify the equivalence between an "Andreev" impurity and an NS junction we will calculate the reflection and transmission coefficients, as well as the Andreev reflection and transmission for the "Andreev" impurity, and check in which regime they correspond to the traditional BTK values mentioned above. Thus, for a delta-function Andreev impurity (see Fig. 3) we write down the Schrödinger equation where ψ(x) = (ψ e (x), ψ h (x)) T is a two-component wave function with the upper and lower components standing for the electron and hole wave functions, respectively. We assume a quadratic dispersion with m denoting the quasiparticle mass. Hence, the impurity potential is a 2×2 matrix in the electron-hole space, and can be written as In the most general case a normal reflection component V N = 0 should be present, however, we have shown that it hinders the equivalence, and therefore, hereinafter we will set V N = 0. For x = 0 the solution of Eq. (5) is just a linear combination of two-component vectors multiplied by rightmoving and left-moving plane waves, e ikx and e −ikx , correspondingly, where the wavevector k ≡ √ 2mE/ . The most general solution in the presence of the impurity is We focus on the simple case of an injected electron incoming on the barrier from the left, thus A r = (1, 0) T , and no injected electron or hole from the right, i.e., B l = (0, 0) T . Hence A l = (r, r A ) T , where r is the regular reflection coefficient, and r A is the Andreev reflection coefficient. Moreover, B r = (t, t A ) T with t being regular transmission and t A Andreev transmission. By assuming that the wave function ψ(x) = ψ L (x)Θ(−x) + ψ R (x)Θ(x) is continuous at x = 0, we obtain 1 + r = t, r A = t A . We then write down the continuity equation for the derivative of the wave function. The delta-function potential gives rise to a discontinuity at the origin, proportional to the value of the impurity potential, in occurence here: which yields We note that the solution to this problem is energydependent, but since we are interested in energies very close to the Fermi level, we can take k to be constant k = k F , and we denote α ≡ yield In Fig. 4 we plot |r| 2 + |r A | 2 and |r A |/|r| as a function of V A , while setting α = 1. This shows that for large enough Andreev potentials the impurity models the NS junction asymptotically well, i.e., |r| 2 +|r A | 2 ≈ 1. On the other hand, we see that with increasing V A the ratio between the Andreev and the regular reflection is decreasing, e.g., for V A = 3.5 for which |r| 2 + |r A | 2 ≈ 0.92, we barely have |r A |/|r| ≈ 0.3. This corresponds to a value of Z ≈ 1 for the NS junction (see Eq. (3)). While the corresponding junction is definitely far from perfect Andreev reflection it can still support Andreev bound states, and we describe their behavior in Sec. IV. This situation describes probably quite accurately the realistic parameters for many experimental NS interfaces.

III. T -MATRIX FORMALISM
In order to take into account the cumulated effect of all-order impurity scattering processes, we employ the Tmatrix formalism. We start with the case of a single impurity. Denoting the momentum-space Hamiltonian of a given system H k , we define the unperturbed Matsubara Green's function (GF) as: where ω n denote the Matsubara frequencies. In the presence of an impurity, the Green's function is modified to: where the T -matrix T (k 1 , k 2 , iω n ) embodies impurity scattering processes. 25,26 For a delta-function impurity V imp (x) = V δ (x), the form of the T -matrix in 1D is momentum independent and is given by To compute physical quantities such as, e.g., the local density of states, in what follows we use the retarded GF G(k 1 , k 2 , E) obtained by the analytical continuation of the Matsubara GF G(k 1 , k 2 , iω n ) (i.e., by setting iω n → E + iδ, with δ → 0 + ). The real space equivalents of Eqs. (13) and (14) can be written as and respectively.
In the presence of two delta-function impurities with the full T -matrix can be found from the following equation using the unperturbed retarded Green's function real-space form: 30,31 where i, j ∈ {1, 2}. For brevity we omit energydependence in all functions. Solving the system of four equations in Eq. (16), we get: The poles of the T -matrix ||T ij || yield the energies of the bound states in the system. Furthermore, we can also express the full perturbed Green's function in terms of the Tmatrix elements: The correction to the local density of states due to the impurities can be expressed as The normal part of the junction is described using a simple lattice Hamiltonian H k = ξ k τ z , where τ z is the Pauli matrix acting in the particle-hole subspace, ξ k ≡ µ−2t cos ka, where a is the lattice constant, t denotes the hopping parameter and µ is the chemical potential. The retarded Green's function in real space is computed as the Fourier transform of its momentum-space representation, and is given by where for E 2 t + µ 2 t = 0. We expressed the energy and the chemical potential in terms of the hopping amplitude, i.e., E t ≡ E/t and µ t = µ/t. The positive infinitesimal shift of energy, +iδ, δ → +0, corresponds to an inverse quasiparticle lifetime and is generally related to the width of the energy levels. Note that since the expressions above are obtained within the lattice model, the results are valid only for x = na, where n ∈ Z.
The results that we will present in this section are evaluated using this full tight-binding model. However, in order to test the validity of our approximations, we establish also a correspondence between the continuum model (used in Sec. II to make the connection between an NS junction and an "Andreev impurity") and the lattice model. We thus expand ξ k in a quadratic form: This allows us to extract: In what follows we use values of the chemical potential around µ = 1.5 so that α ≈ 1. We also set by default the hopping parameter t = 1, therefore, making µ t and E t equivalent to µ and E, respectively. The broadening is hereinafter set to δ = 0.001. To model the SNS junction we introduce two Andreev impurities, as in Eq. (6). In order to take into account the phase difference between the superconductors we replace V A → V A e ±iϕ for one of the impurities, choosing the signs differently for the 12 and 21 components, to preserve the Hermitian character of the Hamiltonian. Using the notations from Eq. (15), we can thus write: As mentioned in Sec. II we focus on some intermediate value of V A ≈ 3.5 corresponding to a ratio of Andreev and regular reflection of 0.3, and we study the formation of Andreev bound states using the formalism described in Sec. III, the T -matrix for two delta-function impurities. This yields an exact formula for the perturbed Green's function of the system (see Eq. (17)), which allows to extract the local density of states (see Eq. (18)) 32 . We first plot the latter in Fig. (5) as a function of energy and position, for a fixed value of the chemical potential. Note the formation of bound states in the region between the two impurities, i.e., for x ∈ [0, 200].
In order to compare the behavior of these bound states with that of previously studied ABS, [22][23][24] we focus on a given position and we plot in Fig. (6) the dependence of the LDOS as of function of energy and chemical potential. Indeed, we recover the oscillatory bell-shaped dependence of the ABS energy on the chemical potential (see, e.g., Figs. 4a and 4d in Ref. [22]). Furthermore, in Fig. 7 we plot the dependence of the LDOS as a function of energy and the phase difference between the two SCs, for a fixed chemical potential. We note that the amplitude of the phase oscillations, while not very large, is still significant, marking the presence of nonzero Andreev reflection and of the corresponding phase coherence.
Above we have reproduced the results for good SNS junctions, in what follows we consider decreasing |r A |/|r| (for example, take V A = 10) such that we recover a non-Andreev junction (see Fig. 8). Note that in this limit the Andreev reflection goes to zero and the bell-shaped oscillations transform into fully-linear crossings (reflecting a typical linear energy dependence of the bound states in a quantum dot on gate voltage). Moreover, the dependence of the ABS energy on the SC phase difference becomes almost insignificant, as expected.
On the other hand we could try to focus on smaller V A values (for instance, V A = 1), artificially increasing |r A |/|r|, however the results we obtain are no longer physical and consistent with the NS junction since |r A | 2 + |r| 2 1 (see Fig. 9). Nevertheless, the main features of the ABS are generally preserved.
Note that the broadening of the energy levels is governed by two parameters, δ and V A . The former enters as an artificial energy broadening in Eqs. (20)(21), originating from the definition of the retarded Green's function (the source of this term is usually a finite quasiparticle lifetime, etc.). However, in our calculation this value is very small (δ = 0.001) and cannot explain  the observed width of the ABS levels. We rather find that the ABS energy broadening increases when the value of the Andreev impurity V A is reduced. This is qualitatively similar to what happens for an SNS junction, i.e., the larger the value of the SC gap and of the normal-SC coupling, the broader the ABS levels (see, e.g., Figs. 4a and 4d in Ref. [22]). The resemblance between our results and Ref. [22] is not accidental since the impurity potential amplitude V A controls the amplitude of the Andreev reflection; we thus expect that when V A is smaller the Andreev reflection increased and correspondingly the ABS levels widen. In the extreme limit in which the Andreev reflection goes to zero (V A → ∞) the width of the bound state levels is governed solely by δ. In order to obtain more insight into the behavior of the ABS energy levels, we employ the continuum model introduced in Sec. II, namely we consider a system with quadratic dispersion (see the Appendix). Since the Andreev impurity potentials required to satisfy the equivalence to a SN junction are quite large, we start with a zeroth order approximation corresponding to infinitely strong Andreev impurities. As expected, in this limit we recover the particle-in-a-box solutions for the boundstate energies, i.e., where we have denoted p F ≡ k F . Subsequently we calculate the first-order correction to this solution in terms of the small dimensionless parameter (p F /mV A ) 2 : We are leaving the technical details, as well as the form of the next higher order correction to the Appendix. This correction to the particle-in-a-box energies has both real and imaginary parts. The real part is responsible for the energy shifts, whereas the imaginary part embodies the energy level broadening.
To make the correspondence with the lattice-model results we choose m = 0.5, and we express p F in terms of the chemical potential µ, using the previously introduced definition of α ≡ p 2 F /m = 2(2t − µ). Since we have taken t = 1 we need to set p F = √ 2 − µ. This allows us to plot the energies E ± 0 + δE ± as a function of the chemical potential µ (see Fig. 10). There is a very good agreement between these closed-form analytical results for the continuum model and those presented in Fig. 6 corresponding to a direct evaluation of the T-matrix and the corresponding LDOS for a full tight-binding model. Note that our approximation does not work for values of the chemical potential close to the singularities in Eq. (26). A subsequent simplification of Eq. (26) allows us to extract some other important information about the ABS level behavior. First, we recall that our approach is valid at small energies, Hence we can take π 2 n 2 p 2 F L 2 ≈ 1, and 2 − π 2 n 2 p 2 F L 2 ≈ 1 in Eq. (26). Thus, we can extract the periodicity of the Andreev levels to be δp F · L ≈ 2π and consequently |δµ| = 4π For the values considered here (L = 200, µ ≈ 1.5) we see that this corresponds to δµ = 0.045, which is observed in both Fig. 10 and Fig. 6. On the other hand, the singularities in Eq. (26) are governed by 2δp F L ≈ 2π, yielding |δµ| = 0.022 corresponding to the spacing between the red dashed lines in Fig. 10.
Moreover, Eq. (26) allows us to see that the amplitude of the phase oscillations, i.e., the coefficient of the cos(ϕ) term, is proportional to (p F /mV A ) 2 . Indeed, this goes to zero when there is no Andreev reflection V A → ∞, and increases with reducing V A . Eqs. (3) and (12) allow to make a direct correspondence between the amplitude of these oscillations and the dimensionless barrier strength Z, namely, it is decreasing roughly as 1/Z 4 .
Furthermore, the width of the ABS levels can also be extracted from the imaginary part of Eq. (26), and is thus equal to 1 ≈ 0.001, consistent with Figs. 6 and 7. Note that with decreasing V A these estimates become inaccurate and the higher order terms become important (see Appendix A), for example, for V A = 3.5 the next order correction which we neglect here is of the order of ≈ 16%.
We should note that no alternative simple expression for our Eq. (26) providing a closed form for the energies of the ABS levels as a function of phase difference and chemical potential, has ever been derived for long SNS junctions with imperfect contacts. This demonstrates the strength of our approach to obtain analytical insight into problems for which there is no other analytical alternative.

V. CONCLUSIONS
We have shown that, by introducing two "Andreev" impurities into a normal metal, and by employing the Tmatrix formalism, we can model an SNS junction and the formation of Andreev bound states. Most importantly, we have obtained analytical expressions allowing to study the LDOS, as well as the dependence of the energies of the ABS on different parameters of the system, including the chemical potential and the phase difference between the SCs. We have found that our approach recovers the same properties of ABS as previously calculated by numerical methods. Moreover it provides the advantage of being able to model imperfect junctions whose Andreev to normal scattering ratio can be controlled by varying the Andreev impurity potential amplitude. The most spectacular result recovered by this technique is a closedform expression of the dependence of the energies of the ABS levels, as well as of their widths, on the chemical potential, the SC phase difference, and the interface barrier strength.
We thus conclude that our technique allows to get a unique analytical insight into the physics of SNS junctions and Andreev bound states, and it is a versatile tool that can be also generalized to other setups in order to obtain exact expressions for complex situations that cannot be explored alternatively except by numerical methods. Our approach should thus provide an easily accessible non-numerical tool useful to describe ongoing experiments which often involve complex setups and imperfect junctions.
result in the limit of the continuum model is provided in the Appendix.

Appendix A: Andreev bound states in the continuum model
We consider a 1D spinless normal metal with two Andreev impurities localised at x = 0 and x = L. We can write down a simple model in the basis of electrons and holes as follows: Above p F denotes the Fermi momentum, m is the effective mass of electrons, V 1,2 are the amplitudes of the impurity potentials. The Andreev impurities in this model can be written as where ϕ corresponds to the phase difference between the SCs. Such impurity potentials reflect electrons into holes, and vice versa. In what follows we study the resulting impurity-induced bound states. We start by computing the unperturbed retarded Green's function of the superconductor. In momentum space it is given by: To solve the problem we need to Fourier-transform the Green's function to real space as follows Two integrals are sufficient to define the real-space form of the Green's function: The infinitesimal shift in the denominator originates from the definition of the retarded Green's function. In terms of the functions defined above the Green's function in coordinate space becomes Now we can proceed to solving the Schrödinger equation for the Hamiltonian in Eq. (A1): In the Fourier space we get: Going back to the real space we have: First, we need to find the wave function values at x = 0 and x = L, and then define the full coordinate dependence using those. Thus, we have a system of equations to solve (in what follows we omit x = 0, L in the argument of the wave function and just write 0, L): Or rewritten in a block-matrix form: In order to ensure that this equation has non-trivial solutions we set the determinant of the block matrix to zero, which in turn yields the equation for the energies of bound states: where we have introduced a dimensionless parameter γ = 2mE/p 2 F ≡ E/E F . The equation above is transcendental, and in the most general case does not have an analytical solution.
First, we turn to the limit of V A → ∞, in which Eq. (A12) is easily solvable: There is nothing surprising about this result: we have just obtained the energy levels of a particle in a box (since the limit of infinite potential corresponds to that). What is left is to find the wave functions at x = 0 and x = L in order to obtain the final expression for the wave function in the limit of V A → ∞. This can be done straightforwardly using Eqs. (A8) while plugging in the energies obtained in Eq. (A13): For electrons we have: whereas for holes we obtain: It is easy to verify that these wave functions correspond to those of a particle in a box.
In order to extract the Fermi energy dependence of Andreev bound state energies, we need to find the first-order correction to the particle-in-a-box solutions in Eq. (A13). Hence, we rewrite Eq. (A12) in the following form: For the parameter range we are interested in this equation can be considered perturbatively with p F /mV A being the small parameter. In the 0-th approximation we have the same solutions as those for a particle in a box (see Eq. (A13)): In what follows we find the first correction to these two series of solutions. We demonstrate it using γ = γ + 0 . First, we divide Eq. (A16) by 1 − e 2ip F L √ 1−γ , since by our choice of root this factor is never equal to zero: We substitute the 0-th approximation solution into the right-hand side of the Eq. (A18). The left-hand side we represent as Taylor series around the point γ = γ + 0 up to the first-order correction, i.e., where δγ + is the sought-for correction to γ + 0 . Computing the derivative and using Eq. (A18) we get: Similarly we have: Finally, the full solution can be written as Note, that γ ± are complex, and in order to obtain the energies of the bound states we need to take the real parts of those expressions, whereas the imaginary parts yield the broadening. It is also worth mentioning that we are interested in bound states forming close to zero energy, in other words, we consider values of n ∈ N such that π 2 n 2 p 2 F L 2 ∼ 1.
Applying this condition to the corrections obtained above, we get: Note that at values p F L = πq, with q ∈ Z this correction diverges, setting the limits of validity of this perturbative approach.