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

TL;DR

This paper introduces an analytical method using the T-matrix formalism to calculate Andreev bound state energies in long, imperfect SNS junctions modeled as two impurities, bridging the gap between numerical and analytical approaches.

Contribution

The authors develop a new analytical framework for ABS energies in long SNS junctions using two impurities and the T-matrix formalism, providing insights beyond numerical methods.

Findings

Derived a closed-form expression for ABS energies.

Showed the equivalence of an Andreev impurity to an NS junction with both Andreev and normal scattering.

Bridged the gap between perfect junction and particle-in-the-box models.

Abstract

We provide a new analytical tool to calculate the energies of Andreev bound states (ABS) in long imperfect SNS junctions, at present these can only be described by numerical tools. We model 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. The ABS energies are then calculated using the T-matrix formalism applied to a system with two Andreev impurities. Our results lie between those for a perfect long SNS junction limit described by the Andreev approximation (ABS energies depend linearly on the phase and are independent of the chemical potential) and the particle-in-the-box limit (bound state energies are independent of the phase and have a linear dependence on the chemical potential). Moreover, we recover a closed-form expression for the ABS energies by expanding around the particle-in-the-box limit.

Equations

$$|r_A{|}^2$$

$$|r{|}^2$$

$$|r_A{|}^2$$

$$|r{|}^2$$

$$-\frac{{\hslash }^2}{2m}\frac{d^2\psi }{d{x}^2}(x)+V_{\mathrm{imp}}(x)\psi (x)=E\psi (x),$$

$$V_{\mathrm{imp}}(x)=V\delta (x)\equiv \left(\begin{array}{cc}{V}_N& V_A\\ V_A^{\ast }& -V_N\end{array}\right)\delta (x).$$

$${\psi }_L(x)=A_r{e}^{ikx}+A_l{e}^{-ikx},x<0,$$

$${\psi }_R(x)=B_r{e}^{ikx}+B_l{e}^{-ikx},x>0.$$

$$-\frac{{\hslash }^2}{2m}[{\psi }_R^{\prime }(0)-{\psi }_L^{\prime }(0)]+V\psi (0)=0,$$

$$-\frac{{\hslash }^2}{2m}ik(-A_r+A_l+B_r-B_l)+V(A_r+A_l)=0.$$

$$r$$

$$r_A$$

$$G\left(k_1,k_2,i{\omega }_n\right)$$

$$+G_0\left(k_1,i{\omega }_n\right)T\left(k_1,k_2,i{\omega }_n\right)G_0\left(k_2,i{\omega }_n\right),$$

$$T\left(i{\omega }_n\right)={\left[\mathbb{I}-V\cdot \int \frac{dk}{2\pi }{G}_0\left(k,i{\omega }_n\right)\right]}^{-1}\cdot V.$$

$$\mathcal{G}(x,x^{\prime },E)={\mathcal{G}}_0(x-x^{\prime },E)+{\mathcal{G}}_0(x,E)T(E){\mathcal{G}}_0(-x^{\prime },E)$$

$$T(E)$$

$$V_{\mathrm{imp}}\left(x\right)=V_1\delta \left(x-X_1\right)+V_2\delta \left(x-X_2\right),$$

$$T_{ij}=V_i{\delta }_{ij}+V_i{\mathcal{G}}_0(X_i-X_1)T_{1j}+V_i{\mathcal{G}}_0(X_i-X_2)T_{2j},$$

$$T_{11}$$

$$T_{12}$$

$$T_{21}$$

$$T_{22}$$

$$\mathcal{G}(x,x^{\prime })={\mathcal{G}}_0(x-x^{\prime })+\sum _{ij}{\mathcal{G}}_0(x-X_i)T_{ij}{\mathcal{G}}_0(X_j-x^{\prime }).$$

$$\displaystyle\Delta\rho(x)=-\frac{1}{\pi}\operatorname{Im}\operatorname{tr}\sum\limits_{ij}\mathcal{G}_{0}(x\negthickspace-\negthickspace X_{i})T_{ij}\mathcal{G}_{0}(X_{j}\negthickspace-\negthickspace x^{\prime})\Bigg{|}_{x^{\prime}=x}.$$

$${\mathcal{G}}_0(E,x)=\left(\begin{array}{cc}{\mathcal{G}}_{11}(E,x)& 0\\ 0& {\mathcal{G}}_{22}(E,x)\end{array}\right),$$

$${\mathcal{G}}_{11}(E,x)=$$

$$-\frac{{\left[\left(E_t-{\mu }_t+i\delta \right)-\mathrm{sgn}(E_t-{\mu }_t)\sqrt{{\left(E_t-{\mu }_t+i\delta \right)}^2-4}\right]}^{\frac{|x|}{a}}}{at\cdot \mathrm{sgn}(E_t-{\mu }_t)\cdot 2^{\frac{|x|}{a}}\cdot \sqrt{{\left(E_t-{\mu }_t+i\delta \right)}^2-4}}$$

$${\mathcal{G}}_{22}(E,x)=$$

$$+\frac{{\left[\left(E_t+{\mu }_t+i\delta \right)-\mathrm{sgn}(E_t+{\mu }_t)\sqrt{{\left(E_t+{\mu }_t+i\delta \right)}^2-4}\right]}^{\frac{|x|}{a}}}{at\cdot \mathrm{sgn}(E_t+{\mu }_t)\cdot 2^{\frac{|x|}{a}}\cdot \sqrt{{\left(E_t+{\mu }_t+i\delta \right)}^2-4}}$$

$${\mathcal{G}}_{11}(E,x)$$

$${\mathcal{G}}_{22}(E,x)$$

$${\xi }_k$$

$$t$$

$$V_1=\left(\begin{array}{cc}0& V_A\\ V_A& 0\end{array}\right)\text{and}{V}_2=\left(\begin{array}{cc}0& V_A{e}^{-i\phi }\\ V_A{e}^{+i\phi }& 0\end{array}\right).$$

$$E_n^{\pm A}=\pm \frac{\hslash p_F}{2mL}[\varphi -(2n+1)\pi ],n\in \mathbb{N},$$

Full text

I Introduction

The formation of Andreev bound states (ABS) Andreev (1964); Kulik (1970); Beenakker (1991); Fukuyama and Ando (1992); Beenakker (1992) in long SNS junctionsGriffin and Demers (1971); Josephson (1974); Blonder et al. (1982); de Gennes (1999) has been approached by analytical tools such as the Bogoliubov–de Gennes equations Furusaki and Tsukada (1991); Hurd and Wendin (1994); Richter et al. (2002), the Andreev approximationAffleck et al. (2000); Nazarov and Blanter (2009); Beenakker et al. (2013), as well as various other approaches Martín-Rodero et al. (1994); Levy Yeyati et al. (1995); Fazio et al. (1995); Cuevas et al. (1996); Caux et al. (2002); Bauer et al. (2007); Perfetto et al. (2009); Meng et al. (2009). However, taking into account a finite amount of backscattering at the leads in such junctions has so far only been possible using numerical tools (see, e.g., Ref. [Bena, 2012]). For instance there is no analytical description of the bell-shaped dependence of the ABS in long junctions on gate-voltage, neither of their non-linear dependence on the phaseBena (2012); Pillet et al. (2010); Dirks et al. (2011).

Here we provide the first fully analytical tool that can take into account the normal scattering at the leads in long junctions (in the limit of energies much smaller than the superconducting gap). When there is no normal scattering at the leads, the physics of the ABS in this limit is described by the Andreev approximation Affleck et al. (2000); Nazarov and Blanter (2009); Beenakker et al. (2013), yielding a linear dependence of the ABS energies on the phase; in this limit their energies are independent of the chemical potential. In the opposite limit, when the normal scattering at the leads is very large, we recover the particle-in-the-box limit, i.e., the normal region is completely isolated and the energies of the bound states correspond to the quantized levels of a finite-size system; the energies of these states depend linearly on the gate voltage and are independent of the SC phase difference. For the intermediate regime, in which the normal scattering is finite, we show that the energies of the ABS lie between the two limits and we obtain a closed-form expression for these energies by expanding around the particle-in-the box limit. This expression allows to extract information for example about the periodicity and the amplitude of the ABS oscillations with the chemical potential, the amplitude of the ABS oscillations with the phase, and the width of the ABS levels as a function of the scattering strength; such information will be useful in order to analyze experimental observations of the ABS in long imperfect SNS junctions, once such observations become available.

We will first show that we can model an imperfect NS junction by considering an “Andreev” type impurity: a delta-function 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. Our formalism cannot describe the ideal NS junctions, a finite amount of normal scattering at the leads is automatically included.

However, the point of this work is not to model a sole NS junction, for which many alternative analytical works exist, but to model an imperfect SNS junction, or equivalently, two imperfect SN junctions situated at a finite and large distance from each other, for which it is to our knowledge impossible to obtain fully analytically the form of the local density of states or the energy dependence of the ABS. We thus model the SNS junction by considering 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 us to describe the formation of ABS. We calculate the resulting dependence of the ABS energies on various parameters such as the gate voltage, the phase difference between the two SCs and the normal scattering at the leads, and we show that it is consistent with previous findings (see, for instance, the results obtained by numerical diagonalization in Ref. [Bena, 2012]).

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 theoryBlonder et al. (1982) indicates that the reflection and Andreev reflection coefficients, resulting from injecting an electron from the normal side of the junction (see Fig. 2 for the schematics of the junction), are given by

[TABLE]

where $\Delta$ 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

[TABLE]

For $Z=0$ we have a perfect Andreev junction, i.e., $r_{A}=1$, while for $Z\gg 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

[TABLE]

where $\psi(x)=\left(\psi_{e}(x),\psi_{h}(x)\right)^{\operatorname{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\times 2$ matrix in the electron-hole space, and can be written as

[TABLE]

In the most general case a normal reflection component $V_{N}\neq 0$ should be present, however, we have shown that it hinders the equivalence, and therefore, hereinafter we will set $V_{N}=0$. Note that this does not imply that in the equivalent SN junction there is no normal scattering, in fact $V_{A}$ alone is generating both an Andreev and a normal component in the equivalence.

For $x\neq 0$ the solution of Eq. (5) is just a linear combination of two-component vectors multiplied by right-moving and left-moving plane waves, $e^{ikx}$ and $e^{-ikx}$, correspondingly, where the wavevector $k\equiv\sqrt{2mE}/\hbar$. The most general solution in the presence of the impurity is

[TABLE]

We focus on the simple case of an injected electron incoming on the barrier from the left, thus $A_{r}=(1,0)^{\operatorname{T}}$, and no injected electron or hole from the right, i.e., $B_{l}=(0,0)^{\operatorname{T}}$. Hence $A_{l}=(r,r_{A})^{\operatorname{T}}$, where $r$ is the regular reflection coefficient, and $r_{A}$ is the Andreev reflection coefficient. Moreover, $B_{r}=(t,t_{A})^{\operatorname{T}}$ with $t$ being regular transmission and $t_{A}$ Andreev transmission.

By assuming that the wave function $\psi(x)=\psi_{L}(x)\Theta(-x)+\psi_{R}(x)\Theta(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:

[TABLE]

which yields

[TABLE]

We note that the solution to this problem is energy-dependent, 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 $\alpha\equiv\frac{\hbar^{2}k_{F}^{2}}{m}$. The above equations yield

[TABLE]

In Fig. 4 we plot $|r|^{2}+|r_{A}|^{2}$ and $|r_{A}|/|r|$ as a function of $V_{A}$, while setting $\alpha=1$. This shows that for large enough Andreev potentials the impurity models the NS junction asymptotically well, i.e., $|r|^{2}+|r_{A}|^{2}\approx 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}\approx 0.92$, we barely have $|r_{A}|/|r|\approx 0.3$. This corresponds to a value of $Z\approx 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 $T$-matrix formalism. We start with the case of a single impurity. Denoting the momentum-space Hamiltonian of a given system $\mathcal{H}_{k}$, we define the unperturbed Matsubara Green’s function (GF) as: $G_{0}\left(k,i\omega_{n}\right)=\left[i\omega_{n}-\mathcal{H}_{k}\right]^{-1}$, where $\omega_{n}$ denote the Matsubara frequencies. In the presence of an impurity, the Green’s function is modified to:

[TABLE]

where the $T$-matrix $T\left(k_{1},k_{2},i\omega_{n}\right)$ embodies impurity scattering processes Balatsky et al. (2006); Mahan (2000). For a delta-function impurity $V_{\mathrm{imp}}\left(x\right)=V\delta\left(x\right)$, the form of the $T$-matrix in 1D is momentum independent and is given by Balatsky et al. (2006); Ziegler et al. (1996); Salkola et al. (1996); Bena (2008):

[TABLE]

To compute physical quantities such as, e.g., the local density of states, in what follows we use the retarded GF $\mathcal{G}(k_{1},k_{2},E)$ obtained by the analytical continuation of the Matsubara GF $G(k_{1},k_{2},i\omega_{n})$ (i.e., by setting $i\omega_{n}\rightarrow E+i\delta$, with $\delta\rightarrow 0^{+}$). The real space equivalents of Eqs. (13) and  (14) can be written as

[TABLE]

and

[TABLE]

respectively, where $\mathcal{G}_{0}(x=0,E)$ is equivalent to integrating the Green’s function over all momenta (cf. Eq. (14)).

In the presence of two delta-function impurities with amplitudes $V_{i}$, localized at $x=X_{i}$, $i\in\{1,2\}$, i.e.,

[TABLE]

the full $T$-matrix can be found from the following equation using the unperturbed retarded Green’s function real-space form Choi (2004, 2005):

[TABLE]

where $i,j\in\{1,2\}$. For brevity we omit energy-dependence in all functions. Solving the system of four equations in Eq. (16), we get:

[TABLE]

where we defined $T^{(0)}_{i}\equiv\left[\mathbb{I}-V_{i}\mathcal{G}_{0}(0)\right]^{-1}V_{i}$. 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 $T$-matrix elements:

[TABLE]

The correction to the local density of states due to the impurities can be expressed as

[TABLE]

IV Results

The normal part of the junction is described using a simple lattice Hamiltonian $\mathcal{H}_{k}=\xi_{k}\tau_{z},$ where $\tau_{z}$ is the Pauli matrix acting in the particle-hole subspace, $\xi_{k}\equiv\mu-2t\cos ka$, where $a$ is the lattice constant, $t$ denotes the hopping parameter and $\mu$ 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

[TABLE]

where

[TABLE]

if $E^{2}+\mu^{2}\neq 0$, and

[TABLE]

for $E_{t}^{2}+\mu_{t}^{2}=0$. The expressions in Eqs. (22) and (23) are obtained directly from the integral defining the Fourier transform of the momentum-space Green’s function. In other words, we first set in the integrand $E_{t}=\mu_{t}=0$ and then we perform the integral over momenta to obtain the real-space formulae. We expressed the energy and the chemical potential in terms of the hopping amplitude, i.e., $E_{t}\equiv E/t$ and $\mu_{t}=\mu/t$. The positive infinitesimal shift of energy, $+i\delta$, $\delta\to+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\in\mathbb{Z}$.

The results that we 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 $\xi_{k}$ in a quadratic form:

[TABLE]

This allows us to extract:

[TABLE]

In what follows we use values of the chemical potential around $\mu=1.5$ so that $\alpha\approx 1$. We also set by default the hopping parameter $t=1$, therefore, making $\mu_{t}$ and $E_{t}$ equivalent to $\mu$ and $E$, respectively. The broadening is hereinafter set to $\delta=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}\to V_{A}e^{\pm i\varphi}$ 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:

[TABLE]

Note that since we assume that the energy is much smaller than the SC gap our results apply to the long SNS junction limit. As mentioned in Sec. II we focus on some intermediate value of $V_{A}\approx 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))111Extra insight and a more detailed analysis of this exact result in the limit of the continuum model is provided in the Appendix. 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\in\left[0,\,200\right]$.

In order to compare the behavior of these bound states with that of previously studied ABS Bena (2012); Pillet et al. (2010); Dirks et al. (2011), 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. In the plots we choose to average the LDOS over a few sites (i.e., 6 sites) to avoid any fluctuation effects due to the short wavelength oscillations visible in Fig. 5. 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. [Bena, 2012]).

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.

In order to understand this observed dependence of the ABS energies on the phase and chemical potential we remind the reader that the focus of our work is the study of long SNS junctions in the imperfect long junction limit. Also we focus on the limit of small energies with respect to the gap, $E\ll\Delta$. In what follows we write down the solutions for the two extreme limits characterizing this junction, that have both been studied previously in the literature. First, for the long perfect SNS junctions (no normal backscattering at the leads), the corresponding results are described by the Andreev approximation. In this limit the energies of the ABS are linear with the phase and independent of chemical potential, and they are given by

[TABLE]

To the other extreme, in the limit of zero Andreev scattering and infinite normal scattering at the leads our setup is equivalent to the particle-in-a-box problem, which for quadratically-dispersed quasiparticles has the following well-known solution:

[TABLE]

independent of phase, and linear with the chemical potential.

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 $\mu$, using the previously introduced definition of $\alpha\equiv p_{F}^{2}/m=2(2t-\mu)$. Since we have taken $t=1$ we need to set $p_{F}=\sqrt{2-\mu}$. Thus we can rewrite

[TABLE]

where $t=1$ and $L$ is expressed in units of the lattice constant $a$.

The corresponding limits for the bound states in the junction are depicted in Figs. 6 and 7 by the dashed red lines (Andreev limit) and the dotted blue lines (particle-in-a-box limit). Note that our results interpolate between the two regimes; our approach is unique in the sense that it is the only analytical approach to capture this transition between a long SNS junction with perfect contacts and a perfectly isolated system. Moreover, as we will show in what follows, by expanding about the particle-in-the-box limit, we can obtain a fully closed form for the ABS dependence of energy depicted in Figs. 6 and 7, valid for all the values of the chemical potential with the exception of a countable set of points.

We also plot the dependence of the ABS on chemical potential and phase difference for parameter values that take us closer to the two ideal limits, for example we take a smaller $|r_{A}|/|r|$ (for example, take $V_{A}=10$) such that we recover a junction very close to the particle-in-the-box scenario (see Fig. 8).

Note that in this limit the Andreev reflection goes to zero and the bell-shaped oscillations transform into almost fully-linear crossings (reflecting the linear energy dependence of the bound states in a quantum dot on gate voltage, corresponding to Eq. (28) and denoted in blue). Moreover, the dependence of the ABS energy on the SC phase difference becomes almost insignificant, as expected.

On the other hand, we 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}\ll 1$ (see Fig. 9). Nevertheless, the main features of the ABS are generally preserved, and we see that we get closer to the Andreev limit denoted in red.

Moreover, we can obtain a fully closed form for the ABS energy levels if we employ the continuum model introduced in Sec. II, namely we consider a system with quadratic dispersion (see the Appendix) and expand around the particle-in-the-box-limit. The zeroth order contribution to the ABS level is thus given by Eq. (28). Subsequently we calculate the first-order correction to this solution in terms of the small dimensionless parameter $(p_{F}/mV_{A})^{2}$:

[TABLE]

We are leaving the technical details of the derivation 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.

Eq. (29) where we set $p_{F}=\sqrt{2-\mu}$ allows us to plot the energies $E_{n}^{\pm B}+\delta E_{n}^{\pm B}$ as a function of the chemical potential $\mu$ (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 perturbative treatment within the continuum model does not work for values of the chemical potential close to the singularities in Eq. (29), which constitute a countable set of points.

A subsequent simplification of Eq. (29) allows us to extract some other important information about the ABS level behavior. First, we recall that our approach is valid at small energies, $\delta E^{\pm B}_{n}/\mu\ll 1$. Hence we can take $\frac{\pi^{2}n^{2}}{p_{F}^{2}L^{2}}\approx 1$, and $\sqrt{2-\frac{\pi^{2}n^{2}}{p_{F}^{2}L^{2}}}\approx 1$ in Eq. (29). Thus, we can extract the periodicity of the Andreev levels to be $\delta p_{F}\cdot L\approx 2\pi$ and consequently $|\delta\mu|=\frac{4\pi}{L}\sqrt{2-\mu}$. For the values considered here ($L=200$, $\mu\approx 1.5$) we see that this corresponds to $\delta\mu=0.045$, which is observed in both Fig. 6 and Fig. 10. On the other hand, the singularities in Eq. (29) are governed by $2\delta p_{F}L\approx 2\pi$, yielding $|\delta\mu|=0.022$ corresponding to the spacing between the green vertical lines in Fig. 10.

Moreover, Eq. (29) allows us to see that the amplitude of the phase oscillations, i.e., the coefficient of the $\cos(\varphi)$ term, is proportional to $(p_{F}/mV_{A})^{2}$. Indeed, this goes to zero when there is no Andreev reflection $V_{A}\rightarrow\infty$, 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. (29), and is thus equal to $\frac{1}{p_{F}L}\left(\frac{p_{F}}{mV_{A}}\right)^{2}\approx 0.001$ for $V_{A}=3.5$, 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 $\approx 16\%$.

Note that the broadening of the ABS levels is governed by two physical phenomena. The first is the finite quasiparticle lifetime, which is quite large in regular systems making the levels very sharp. We introduce by hand an artificial broadening, described by the imaginary contribution $i\delta$, with $\delta=0.001$, into the energy in the Green’s functions. It is proportional to the inverse quasiparticle lifetime and it is the main factor of level broadening in regular systems. However, in a setup consisting of a central region coupled to the leads, there is a second physical factor responsible for the level broadening, originating from the coupling to the leads. When the coupling to the leads is purely Andreev-like, the width of the ABS is not affected by its size, see Fig. 7 in Ref. [Bena, 2012]. However, when there is also a normal coupling we expect that the ABS level width is going to be affected by this normal coupling, same as the quantum dot levels becoming wider when the dot is better coupled to the leads; in the perfect coupling limit such levels can become very wide, while in the bad coupling limit the levels are very sharp (see for example Figs. 4a and 4b in Ref. [Bena, 2012] for details). In our model, it is the parameter $V_{A}$ that determines the barrier height and thus the coupling to the SC leads. Therefore, we expect that if the coupling is very small, i.e., $V_{A}\to\infty$, then the energy level broadening is controlled solely by $\delta$ and we obtain the ‘particle in a box solutions’. Indeed, this is exactly what happens (see Eq. (A13-A15) in Appendix A), and in this limit the levels are very sharp. On the other hand, when $V_{A}$ is small we have an additional broadening due to the normal coupling to the leads and this is indeed what we observe.

We should note that no alternative simple expression for our Eq. (29) 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 $T$-matrix formalism, we can model a long imperfect SNS junction and the formation of Andreev bound states. Most importantly, we have obtained analytical expressions of the local density of states, as well as of the energies of the ABS on different parameters of the system, including the chemical potential and the phase difference between the SCs. In our model the Andreev to normal scattering ratio can be controlled by varying the Andreev impurity potential amplitude. We show that our method interpolates between the Andreev limit corresponding to a perfect long SNS junction for which the energies of the ABS are independent of chemical potential and linear with phase, and the particle-in-the-box limit in which the energy levels are independent of phase and linear with chemical potential. One of the most spectacular results recovered by this technique is a closed-form expression of the energies of the ABS levels, as well as of their widths, obtained by a first-order expansion around the particle-in-the-box limit, as a function of the chemical potential, the SC phase difference, and the interface barrier strength, till now only possible by non-analytical tools such as numerical tight-binding exact diagonalization. This makes our method a unique analytical tool to understand and characterize ABS. Moreover, our technique can be easily generalized to other more complex setups, for example, systems with asymmetric leads, or in which the spin of the electrons also plays a role. We have checked for example that in our formalism the Zeeman field only splits the energies of the ABS levels, consistent with an intuitive understanding. However, more complex factors such as Rashba spin-orbit coupling would require a modification of the equations of motion, and of all the conditions for the equivalence with the Andreev impurity; if successful, such an analysis will be addressed in a separate work.

Our approach thus provides an easily accessible non-numerical tool to describe the formation of ABS in long SNS junctions. While at present ABS have been directly observed only in short junctions—quantum-dot-type setups for which only a few states can be observed inside the SC gap—we can expect that longer and cleaner wires will also be realized in the future, for which models such as the one presented here will be necessary in order to investigate the underlying physics. Already some progress has been made in this directionBeenakker et al. (2013), and some experimental setups in which a short but non-zero length plays a role have been examinedTosi et al. (2019). In order to make contact with such experiments our approach may eventually need to be extended to include other factors, such as the presence of interactions and the spin-splitting of the levels.

Acknowledgements.

V.K. would like to acknowledge the ERC Starting Grant No. 679722 and the Roland Gustafsson foundation for theoretical physics.

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:

[TABLE]

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

[TABLE]

where $\varphi$ 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:

[TABLE]

To solve the problem we need to Fourier-transform the Green’s function to real space as follows

[TABLE]

Two integrals are sufficient to define the real-space form of the Green’s function:

[TABLE]

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

[TABLE]

Now we can proceed to solving the Schrödinger equation for the Hamiltonian in Eq. (30):

[TABLE]

In the Fourier space we get:

[TABLE]

Going back to the real space we have:

[TABLE]

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$):

[TABLE]

Or rewritten in a block-matrix form:

[TABLE]

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:

[TABLE]

where we have introduced a dimensionless parameter $\gamma=2mE/p^{2}_{F}\equiv 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}\to\infty$, in which Eq. (41) is easily solvable:

[TABLE]

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}\to\infty$. This can be done straightforwardly using Eqs. (37) while plugging in the energies obtained in Eq. (42): For electrons we have:

[TABLE]

whereas for holes we obtain:

[TABLE]

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. (42). Hence, we rewrite Eq. (41) in the following form:

[TABLE]

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 [math]-th approximation we have the same solutions as those for a particle in a box (see Eq. (42)):

[TABLE]

In what follows we find the first correction to these two series of solutions. We demonstrate it using $\gamma=\gamma_{0}^{+}$. First, we divide Eq. (45) by $\left(1-e^{2ip_{F}L\sqrt{1-\gamma}}\right)$, since by our choice of root this factor is never equal to zero:

[TABLE]

We substitute the [math]-th approximation solution into the right-hand side of the Eq. (47). The left-hand side we represent as Taylor series around the point $\gamma=\gamma_{0}^{+}$ up to the first-order correction, i.e.,

[TABLE]

where $\delta\gamma^{+}$ is the sought-for correction to $\gamma_{0}^{+}$. Computing the derivative and using Eq. (47) we get:

[TABLE]

Similarly we have:

[TABLE]

Finally, the full solution can be written as

[TABLE]

Note, that $\gamma^{\pm}$ 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\in\mathbb{N}$ such that

[TABLE]

Applying this condition to the corrections obtained above, we get:

[TABLE]

Note that at values $p_{F}L=\pi q$, with $q\in\mathbb{Z}$ this correction diverges, setting the limits of validity of this perturbative approach.