format: text

text:
# Quasi Exactly Solvable N×N-Matrix Schrödinger Operators.

## 1 Introduction

The topic of Quasi Exactly Solvable (QES) operators  has been developped in the last years. It consists of differential operators (mainly Schrödinger ones) which possess a finite dimensional invariant vector space, say $`𝒱`$, of functions. So, the restriction of the eigenvalue equation to the space $`𝒱`$ leads to an algebraic eigenvalue problem.

Although scalar QES operators have been classified in one  and several  variables, a classification of matrix QES operators is still missing. Yet, interesting examples of them have been obtained in relation with the stability analysis of soliton solutions occuring in some field theories .

This problem was first addressed in  and further developped in  and . More recently , interesting tools for the classification of N$`\times `$N matrix QES operators in one spatial dimension have been constructed and applied to the case N=2  (although potentials of the Lamé-type are not included in this classification).

Here we consider a suitable class of finite dimensional vector spaces of N polynomials in a real variable and we construct families of operators preserving sub-classes of these vector spaces. The corresponding QES equations respectively constitute “coupled channel” generalizations of the scalar QES equations.

Following the basic idea of QES operators  we consider the finite dimensional vector space of N-uples of polynomials of given degree $`n_1,n_2,\mathrm{},n_N`$ in a real variable $`x`$. We slightly generalize this vector space by setting

$$𝒱=P\left(𝒫(n_1)𝒫(n_2)\mathrm{}𝒫(n_N)\right)$$
(1)
where $`𝒫(n_i)`$, $`i=1,\mathrm{}N`$ denotes the set of real polynomials of degree at most $`n_i`$ in $`x`$ while $`P`$ is a fixed invertible $`N\times N`$ matrix operator; $`P`$ can be interpreted as a change of basis in the vector space $`𝒫(n_1)\mathrm{}𝒫(n_N)`$. With such an interpretation, and if we assume $`n_1n_2\mathrm{}n_N`$, it is reasonable to choose the matrix $`P`$ as a lower triangular matrix with $`P_{ii}=1`$.

## 2 QES Anharmonic matrix potentials

We consider an operator of the form

$$H(y)=\frac{d^2}{dy^2}1\mathrm{I}_N+M_3(y^2)$$
(2)
where $`M_3(x)`$ is a $`N\times N`$ hermitian matrix whose elements are polynomials of degree at most three in the argument $`x`$. After a standard ”gauge transformation” of $`H(y)`$ with the factor

$$\varphi (y)=\mathrm{exp}\{\frac{p_2}{2}y^4+p_1y^2\}$$
(3)
($`p_1,p_2`$ are arbitrary real parameters, $`p_2>0`$) and the change of variable $`x=y^2`$, the operator equivalent to (2) reads

$`\widehat{H}(x)`$ $`=`$ $`(4x{\displaystyle \frac{d^2}{dx^2}}+2{\displaystyle \frac{d}{dx}})1\mathrm{I}_N`$ (4)
$`+`$ $`8(p_2x+p_1)x{\displaystyle \frac{d}{dx}}1\mathrm{I}_N`$
$``$ $`(4p_2^2x^3+8p_1p_2x^2+(4p_1^26p_2)x2p_1)1\mathrm{I}_N`$
$`+`$ $`M_3(x)`$
We now determine the form of the matrix $`M_3(x)`$ such that the operator $`\widehat{H}(x)`$ possesses a finite dimensional invariant vector space of the type (1) for generic values of $`p_1,p_2`$. In this purpose, the differential operators in the second line of (4) (i.e. $`8p_2x^2d/dx`$ and $`8p_1xd/dx`$) have separately to be completed into operators which have this property. In order to make use of the results of  we conveniently rewrite $`\widehat{H}`$ according to

$`\widehat{H}(x)`$ $`=`$ $`(4x{\displaystyle \frac{d^2}{dx^2}}+2{\displaystyle \frac{d}{dx}})1\mathrm{I}_N`$ (5)
$`+`$ $`8p_2Q_++8p_1Q_0+8p_2(W_3(x)+B)`$
where the operators $`Q_+,Q_0`$ are defined by

$$Q_+=x^2\frac{d}{dx}+2xAB,Q_0=x\frac{d}{dx}+A$$
(6)
and the constant matrices $`A,B`$ are chosen to obey $`[A,B]=B`$, in such a way that

$$[Q_+,Q_0]=Q_+.$$
(7)
Without loosing generality , we can choose $`A`$ and $`B`$ in the form

$$A=\mathrm{Diag}(0,1,2,\mathrm{},N1)\frac{p}{2}1\mathrm{I}_N,B_{a,b}=c_b\delta _{a,b+1}$$
(8)
The matrix $`W_3(x)`$ is defined by identification of (4) and (5); it is symmetric and does not contain derivative. In the following, we further assume this matrix to be irreducible by $`x`$-independent changes of basis; in particular, we exclude the cases where $`W_3`$ is diagonal.

In , it was demonstrated that $`Q_+,Q_0`$ admit a finite-dimensional vector space; we now construct it explicitely. For this purpose, we define a family of vector spaces characterized by two integers $`N`$ and $`p`$

$$𝒱(N,p)=𝒫(p)𝒫(p2)\mathrm{}𝒫(p2N+2)$$
(9)
and an (invertible) $`N\times N`$ matrix operator $`P`$ with matrix elements

$`P_{ij}`$ $`=`$ $`0\mathrm{if}i<j`$ (10)
$`=`$ $`1\mathrm{if}i=j`$
$`=`$ $`\left({\displaystyle \underset{k=0}{\overset{ij1}{}}}{\displaystyle \frac{c_{j+k}}{(p+22jk)}}\right){\displaystyle \frac{1}{(ij)!}}({\displaystyle \frac{d}{dx}})^{ij}\mathrm{if}i>j`$
The vector space $`𝒱^{}P𝒱`$ can be seen as a change of basis on the space $`𝒱`$ (for brevity, we do not write the dependence on $`p,N`$ anylonger).

The following proposition, which can be checked after an algebra, provides the explicit form of the invariant vector space of $`Q_+,Q_0`$ :

Proposition 1

$`P^1BP`$ $`=`$ $`P^1(Q_++B)P(Q_++B)`$ (11)
$`P^1Q_+P`$ $`=`$ $`x^2{\displaystyle \frac{d}{dx}}+2xA`$ (12)
$`P^1Q_0P`$ $`=`$ $`Q_0`$ (13)
Because the operators on the right hand side of these equalities preserve $`𝒱`$, it follows immediately that $`Q_+`$ and $`Q_0`$ preserve the space $`𝒱^{}`$.

The requirement that the $`(W_3+B)`$-part of $`\widehat{H}`$ also preserves $`𝒱^{}`$ is guaranteed by the Proposition 2 below. For later use we will note by $`J_+,J_0,J_{}`$ the usual irreducible representation of sl(2) by $`N\times N`$ matrices, i.e.

$$[J_0,J_\pm ]=\pm J_\pm ,[J_+,J_{}]=2J_0.$$
(14)
In particular we can set $`J_0=A\frac{1}{N}\mathrm{Tr}A`$.

Proposition 2

If $`W_3(x)`$ is symmetric and irreducible, then the condition

$$P^1(W_3+B)P𝒱𝒱$$
(15)
is fulfilled if and only if

$$B=cJ_+,W_3=c(J_++J_{})$$
(16)
The values of the different parameters $`c_k`$ entering in the matrix $`B`$ are therefore fixed up to the multiplicative factor $`c`$.

So far, we proved that the condition is necessary and sufficient for $`N`$ up to seven and we are confident that it is true for arbitrary $`N`$.

The final step is to check whether the operator in the first line of (5) also preserves $`𝒱^{}`$. In this respect, it is sufficient to study the condition $`A\frac{d}{dx}𝒱^{}𝒱^{}`$. From the structure of the vector space $`𝒱^{}`$, it is quickly seen that this holds only for N=2. This shows that the generalization of the QES anharmonic oscillator to matrices is possible only for two-dimensional matrices and the 3-parameter family of operators of  is recovered.

Our results thus lead to a family of $`3`$ parameters matrix quasi exactly solvable anharmonic oscillators ($`p_1,p_2,c`$). The parameters $`p_1,p_2`$ determine the anharmonic frequencies of the $`y^6`$\- and $`y^4`$-parts of the potential $`M_3`$, while the constant $`c`$ ensures a non-trivial coupling between the different oscillators.

## 3 Clebsch-Gordan coefficients

The results presented in the previous section contains a formulation of the Clebsch-Gordan matrix of the tensor product of two representations of sl(2) in terms of differential operators. In this section we would like to further analyze some aspects of this realization.

We consider on the one hand the standard representation of the algebra sl(2) by $`N\times N`$-matrix generators noted $`J_+,J_0,J_{}`$ and on the other hand the realization expressed in terms of differential operators :

$$j_+=x^2\frac{d}{dx}+2\mu x,j_0=x\frac{d}{dx}+\mu ,j_{}=\frac{d}{dx}$$
(17)
With the choice $`\mu =n/2`$ ($`n`$ integer), this realization preserves the vector space $`𝒫(n)`$ and constitutes the corner stone of QES equations .

Then, we consider the tensor product of these representations, the generators of which are

$$\stackrel{~}{Q}_{\pm ,0}=J_{\pm ,0}+j_{\pm ,0}.$$
(18)
Again, if $`\mu =n/2`$ , it is finite dimensional and preserves, in a reducible way, the vector space

$$𝒱_0=𝒫(n)𝒫(n)\mathrm{}𝒫(n),N\mathrm{times}$$
(19)
The explicit decomposition of the representation (18), (19) into irreducible ones is achieved in two steps. We define a first transformation by means of

$$Q_{\pm ,0}=e^{xJ_{}}\stackrel{~}{Q}_{\pm ,0}e^{xJ_{}},$$
(20)
leading to

$$Q_+=x^2\frac{d}{dx}+2x(J_0+\mu 1\mathrm{I})+J_+,Q_0=x\frac{d}{dx}+J_0+\mu 1\mathrm{I},Q_{}=\frac{d}{dx}$$
(21)
which can be identified with the operators (6) if $`A=J_0+\mu 1\mathrm{I}`$, i.e. if $`\mu =\frac{Np1}{2}`$. In other words the two irreducible representations $`J_ϵ`$, $`j_ϵ`$ under investigation have dimensions $`N`$ and $`pN+2`$, respectively.

The second step leading to the desired decomposition reads

$$\overline{Q}_{\pm ,0}=P^1Q_{\pm ,0}P$$
(22)
where the operator $`P`$ is defined in eq.(10). The form of the operators $`\overline{Q}`$ is available from Proposition 1 :

$$\overline{Q}_+=x^2\frac{d}{dx}+2xA,\overline{Q}_0=x\frac{d}{dx}+A,\overline{Q}_{}=\frac{d}{dx}$$
(23)
It clearly reveals that the vector space (9) is preserved. As a consequence the matrix $`P_{cg}=\mathrm{exp}(xJ_{})P`$ constitutes the Clebsch-Gordan matrix of the decomposition.

It should be noticed that the operators $`J_{}`$, $`j_{}`$ preserve separately the vector space (19). Then, since the transformation (20) commutes with $`J_{}`$, it results that both $`J_{}`$ , $`j_{}`$ also preserve $`𝒱^{}`$. This statement in fact provides the proof that the condition (16) of the Proposition 2 above is indeed sufficient.

## 4 Lamé-type operators.

In the following the Jacobi elliptic functions

$$\mathrm{sn}(z,k),\mathrm{cn}(z,k),\mathrm{dn}(z,k).$$
(24)
of argument $`z`$ and modulus $`k`$  are abreviated respectively by $`\mathrm{sn},\mathrm{cn},\mathrm{dn}`$. These functions are periodic with period $`4K(k),4K(k),2K(k)`$ respectively ($`K(k)`$ is the complete elliptic integral of the first type).

It is well known that in order to reveal the algebraic properties of the Lamé equation

$$[\frac{d^2}{dz^2}+k^2N(N+1)\mathrm{sn}^2]\psi (z)=E\psi (z),$$
(25)
the relevant change of variable consists in posing $`x=\mathrm{sn}^2(z,k)`$. In particular the second derivative operator is transformed into

$$\frac{d^2}{dz^2}=4x(1x)(1k^2x)\frac{d^2}{d^2x}+2(3k^2x^22(1+k^2)x+1)\frac{d}{dx}$$
(26)
and Eq. (25) becomes a Fushs equation with four regular-singular points (at $`x=0,1/k^2,1,\mathrm{}`$).

To our knowledge, attempts to construct (and classify) Schrödinger QES matrix operators with this type of change of variable have not been attempted so far. Particular cases are emphasized in . The natural choice is to consider $`N\times N`$ Schrödinger matrix operators with potentials depending on the Jacobi ellipitc functions and which possess algebraic properties similar to Eq.(25). More specifically, we consider operators of the form

$$H(z)=\frac{d^2}{dz^2}1\mathrm{I}_N+V_D(z)+V_I(z)$$
(27)
with

$$V_D=\mathrm{sn}^2\mathrm{diag}(a_1,a_2,\mathrm{},a_N)+\mathrm{diag}(b_1,b_2,\mathrm{},b_N)$$
(28)
where $`a_j,b_j`$ denote real constants (without loosing generality we assume $`_{j=1}^Nb_j=0`$) and $`V_I`$ is a symmetric off-diagonal matrix of the form

$$(V_I)_{ij}=\{\begin{array}{cc}\theta _{ij}\mathrm{sn}^{\alpha _{1ij}}\mathrm{cn}^{\alpha _{2ij}}\mathrm{dn}^{\alpha _{3ij}}\hfill & \text{if }ij\text{ }\hfill \\ 0\hfill & \text{ if }i=j\hfill \end{array}$$

Owing the periodicity of the Jacobi elliptic functions, the family hamiltonian above is to be considered on the Hilbert space of periodic functions on $`[0,4K(k)]`$.

Unfortunately, it has not been possible to classify the QES operators of the form (27), here we will describe the cases that we were able to treat. The difficulty of achieving the classification will appear from these few examples. The properties of the Jacobi functions that are useful to make the calculations are listed in the appendix.

### 4.1 case N=1

This case corresponds to the Lamé equation which was discussed lengthly e.g. in .

### 4.2 Case N=2

The non-diagonal potential reduces to only one component

$$V_I=\theta \mathrm{sn}^{\alpha _1}\mathrm{cn}^{\alpha _2}\mathrm{dn}^{\alpha _3}\sigma _1$$
(29)
Using a similarity transformation of the form

$$\widehat{H}(x)=U^1(z)H(z)U(z),U(z)=\mathrm{diag}(\mathrm{sn}^{\beta _1}\mathrm{cn}^{\beta _2}\mathrm{dn}^{\beta _3},\mathrm{sn}^{\gamma _1}\mathrm{cn}^{\gamma _2}\mathrm{dn}^{\gamma _3})$$
(30)
sets the operator (27) into a form with polynomial coefficients in the variable $`x=\mathrm{sn}^2`$ provided (see Appendix)

* $`\beta _j,\gamma _j=0\mathrm{or}1`$ , $`j=1,2,3`$
* $`\alpha _j\pm (\beta _j\gamma _j)=`$ non-negative even integer , $`j=1,2,3`$.

Considering first the cases where $`V_{12}`$ is a constant ($`\alpha _{1,2,3}`$ = 0) or a single linear factor in one of the Jacobi functions, we were able to show that the only possible QES operators are available for $`\theta =0`$.

For the three cases a) $`V_{12}=\theta \mathrm{sn}\mathrm{cn}`$, b) $`V_{12}=\theta \mathrm{sn}\mathrm{dn}`$ and c) $`V_{12}=\theta \mathrm{cn}\mathrm{dn}`$, the construction of non-decoupled QES operators is possible. One remarkable feature is that in each case, two sets of values of the coupling constants $`a_1,a_2,b_1,\theta `$ lead to four different algebraizations of the corresponding operator. We now discuss it in detail (posing $`bb_1=b_2`$).

#### 4.2.1 Case $`V_{12}=\theta \mathrm{sn}\mathrm{cn}`$

The two sets of values of the coupling constants leading to QES operators are noted Type 1 and Type 2.

Type 1

$`a_1=k^2(4m^2+2m+1)2b`$
$`a_2=k^2(4m^2+2m+1)+2b`$
$`\theta ^2=4b^2k^4(1+4m)^2`$
The parameters $`b`$ and $`k`$ remain free, $`m`$ is an integer. Four invariant spaces are available :

$`𝒱_1`$ $`=`$ $`\left(\begin{array}{cc}\mathrm{sn}& 0\\ 0& \mathrm{cn}\end{array}\right)\left(\begin{array}{cc}1& \kappa _1\\ 0& 1\end{array}\right)\left(\begin{array}{c}𝒫(m1)\\ 𝒫(m)\end{array}\right),\kappa _1={\displaystyle \frac{\theta }{2b+k^2(1+4m)}}`$ (37)
$`𝒱_2`$ $`=`$ $`\left(\begin{array}{cc}\mathrm{cn}& 0\\ 0& \mathrm{sn}\end{array}\right)\left(\begin{array}{cc}1& \kappa _2\\ 0& 1\end{array}\right)\left(\begin{array}{c}𝒫(m1)\\ 𝒫(m)\end{array}\right),\kappa _2={\displaystyle \frac{\theta }{2bk^2(1+4m)}}`$ (44)
$`𝒱_3`$ $`=`$ $`\left(\begin{array}{cc}\mathrm{dn}& 0\\ 0& \mathrm{sn}\mathrm{cn}\mathrm{dn}\end{array}\right)\left(\begin{array}{cc}1& \kappa _3x\\ 0& 1\end{array}\right)\left(\begin{array}{c}𝒫(m1)\\ 𝒫(m1)\end{array}\right),\kappa _3={\displaystyle \frac{1}{\kappa _2}}`$ (51)
$`𝒱_4`$ $`=`$ $`\left(\begin{array}{cc}\mathrm{sn}\mathrm{cn}\mathrm{dn}& 0\\ 0& \mathrm{dn}\end{array}\right)\left(\begin{array}{cc}1& 0\\ \kappa _4x& 1\end{array}\right)\left(\begin{array}{c}𝒫(m1)\\ 𝒫(m1)\end{array}\right),\kappa _4={\displaystyle \frac{1}{\kappa _1}}`$ (58)
Type 2

$`a_1=k^2(4m^2+6m+1)2b`$
$`a_2=k^2(4m^2+6m+1)+2b`$
$`\theta ^2=4b^2k^4(3+4m)^2`$
The corresponding invariant vector spaces read

$`𝒱_5`$ $`=`$ $`\left(\begin{array}{cc}1& 0\\ 0& \mathrm{sn}\mathrm{cn}\end{array}\right)\left(\begin{array}{cc}1& \kappa _5x\\ 0& 1\end{array}\right)\left(\begin{array}{c}𝒫(m)\\ 𝒫(m)\end{array}\right),\kappa _5={\displaystyle \frac{k^2(4m+3)2b}{\theta }}`$ (65)
$`𝒱_6`$ $`=`$ $`\left(\begin{array}{cc}\mathrm{sn}\mathrm{cn}& 0\\ 0& 1\end{array}\right)\left(\begin{array}{cc}1& 0\\ \kappa _6x& 1\end{array}\right)\left(\begin{array}{c}𝒫(m)\\ 𝒫(m)\end{array}\right),\kappa _6={\displaystyle \frac{k^2(4m+3)2b}{\theta }}`$ (72)
$`𝒱_7`$ $`=`$ $`\left(\begin{array}{cc}\mathrm{sn}\mathrm{dn}& 0\\ 0& \mathrm{cn}\mathrm{dn}\end{array}\right)\left(\begin{array}{cc}1& \kappa _7\\ 0& 1\end{array}\right)\left(\begin{array}{c}𝒫(m1)\\ 𝒫(m)\end{array}\right),\kappa _7={\displaystyle \frac{2bk^2(4m+3)}{\theta }}`$ (79)
$`𝒱_8`$ $`=`$ $`\left(\begin{array}{cc}\mathrm{cn}\mathrm{dn}& 0\\ 0& \mathrm{sn}\mathrm{dn}\end{array}\right)\left(\begin{array}{cc}1& \kappa _8\\ 0& 1\end{array}\right)\left(\begin{array}{c}𝒫(m1)\\ 𝒫(m)\end{array}\right),\kappa _8={\displaystyle \frac{2b+k^2(4m+3)}{\theta }}`$ (86)
#### 4.2.2 Case $`V_{12}=\theta \mathrm{sn}\mathrm{dn}`$

In the following cases we just write the values of the coupling constants leading to QES operators; the corresponding invariant vector spaces are similar to the ones above, their explicit form can be obtained in a straightforward way.

Type 1

$`a_1=k^2(4m^2+2m+12b)`$
$`a_2=k^2(4m^2+2m+1+2b)`$
$`\theta ^2=4k^2b^2k^2(1+4m)^2`$

Type 2

$`a_1=k^2(4m^2+6m+12b)`$
$`a_2=k^2(4m^2+6m+1+2b)`$
$`\theta ^2=4k^2b^2k^2(3+4m)^2`$

#### 4.2.3 Case $`V_{12}=\theta \mathrm{cn}\mathrm{dn}`$

Type 1

$`a_1=k^2(4m^2+2m+1)2b\frac{k^2}{1+k^2}`$
$`a_2=k^2(4m^2+2m+1)+2b\frac{k^2}{1+k^2}`$
$`\theta ^2=k^2(1+4m)^24b^2\frac{k^2}{1+k^2}`$

Type 2

$`a_1=k^2(4m^2+6m+1)2b\frac{k^2}{1+k^2}`$
$`a_2=k^2(4m^2+6m+1)+2b\frac{k^2}{1+k^2}`$
$`\theta ^2=k^2(1+4m)^24b^2\frac{k^2}{1+k^2}`$

Concerning the operators, we want to point out the following things:

* The extraction of the prefactor is done in two steps. After factorizing the appropriate products of Jacobi functions, the non-diagonal part of the potential takes the form

$$\left(\begin{array}{cc}0& x^{2ϵ}\\ x^ϵ& 0\end{array}\right),ϵ=2\mathrm{or}1\mathrm{or}0$$
(87)
which is clearly incompatible with an operator preserving a vector space of the form $`𝒫(n)𝒫(m)`$ for integers $`m,n`$. The setting of the operator in a form preserving such a vector space is realized by the second prefactor (the triangular matrix one). As an example if $`ϵ=0`$ in (87), then use is made of the relations

$$\left(\begin{array}{cc}1& A\\ 0& 1\end{array}\right)\left(\begin{array}{cc}0& 1\\ 0& 0\end{array}\right)\left(\begin{array}{cc}1& A\\ 0& 1\end{array}\right)=\left(\begin{array}{cc}0& 1\\ 0& 0\end{array}\right)$$
(88)
$$\left(\begin{array}{cc}1& A\\ 0& 1\end{array}\right)\left(\begin{array}{cc}0& 0\\ 1& 0\end{array}\right)\left(\begin{array}{cc}1& A\\ 0& 1\end{array}\right)=\left(\begin{array}{cc}A& A^2\\ 1& A\end{array}\right)$$
(89)
with $`A÷x`$ in order to eliminate the quadratic power of $`x`$ occuring in $`V_{12}`$.
* For $`V_{12}=\theta \mathrm{sn}\mathrm{cn}`$, the limit $`k0`$ is non trivial and leads to a potential of the form

$$\left(\begin{array}{cc}\mathrm{cos}^2(x)\frac{1}{2}& \mathrm{cos}(x)\mathrm{sin}(x)\\ \mathrm{cos}(x)\mathrm{sin}(x)& \mathrm{sin}^2(x)\frac{1}{2}\end{array}\right)$$
(90)
In particular, the dependence on $`m`$ disappears, so this integer can be choosen arbitrarily large and an infinite set of algebraic eigenvectors occurs in this limit. This confirms the fact that the corresponding Schrödinger equation is completely solvable and is related to the stability of some soliton solutions occuring in the Goldstone model in 1+1 dimensions  with a periodic condition for the space coordinate.
* The case $`V_{12}=\theta \mathrm{cn}\mathrm{dn}`$ with $`b=(1+k^2)/2`$ was studied lengthly in . The corresponding Schrödinger equation is related to the stability analysis of the sphaleron solution in the U(1)-Abelian Higgs model in 1+1 dimensions (again with a periodic condition for the space coordinate).

### 4.3 Case N=3

We attempted to construct 3$`\times `$3 matrix operators with some choice for the constants $`\alpha _{aij}`$ ($`a,i,j=1,2,3`$) inspired from the results above. Several trials were unsuccesfull (namely with $`V_{12}÷V_{23}÷\mathrm{sn}\mathrm{cn}`$); however the choice

$$V_I=\left(\begin{array}{ccc}0& \theta _1\mathrm{cn}\mathrm{dn}& V_{23}\\ \theta _1\mathrm{cn}\mathrm{dn}& 0& \theta _2\mathrm{cn}\mathrm{dn}\\ V_{23}& \theta _2\mathrm{cn}\mathrm{dn}& 0\end{array}\right)$$
(91)
leads to the wanted form of equations and preserves a vector space of the form

$$𝒱=\left(\begin{array}{ccc}\mathrm{cn}& 0& 0\\ 0& \mathrm{dn}& 0\\ 0& 0& \mathrm{cn}\end{array}\right).\left(\begin{array}{ccc}1& \alpha & \beta \\ 0& 1& \gamma \\ 0& 0& 1\end{array}\right).\left(\begin{array}{c}𝒫(n2)\\ 𝒫(n1)\\ 𝒫(n)\end{array}\right)$$
(92)
The condition

$$H(z)𝒱𝒱$$
(93)
leads to ten equations for the parameters $`\alpha ,\beta ,\gamma ,k,a_i,b_i`$ and $`\theta _a`$. These equations are compatible with each other only if $`k=1`$ and $`V_{23}=\theta _3\mathrm{cn}^2`$. We found it convenient to leave the parameters $`\alpha ,\beta ,\gamma `$ free and to express the coupling constants $`a_i,b_i`$ and $`\theta _a`$ in terms of them. We refrain to write these tedious expressions in general but we mention that they have the polynomial

$$(1+\gamma ^2+\beta ^2)(1+\beta ^2+\gamma ^2+\alpha ^2\gamma ^22\alpha \beta \gamma )$$
(94)
as common denominator. Also they are such that

$$a_1+a_2+a_3=2(6n^23n+4)$$
(95)
For the case with e.g. $`\alpha =\beta =\gamma =1`$ we find

$$a_1=a_2=\frac{12n^210n+11}{3},a_3=\frac{2}{3}(6n^2+n+1)$$
(96)
$$b_1=b_2=\frac{4n3}{3},b_3=2b_1$$
(97)
$$\theta _1=\frac{72n}{6},\theta _2=\frac{4n+1}{3},\theta _3=\frac{2(4n+1)}{3}$$
(98)
as another example, if $`\alpha =\beta =1,\gamma =0`$, we have

$$a_1=\frac{2(6n^27n+3)}{3},a_2=a_3=\frac{1}{3}(12n^22n+9)$$
(99)
$$b_1=\frac{2(4n+1)}{3},b_2=\frac{4n+1}{3},b_3=b_1b_2$$
(100)
$$\theta _1=\frac{34n}{3},\theta _2=\frac{20n3}{6},\theta _3=\frac{2(4n1)}{3}$$
(101)
Remarkably, replacing the first factor in (92) by $`\mathrm{diag}(\mathrm{dn},\mathrm{cn},\mathrm{dn})`$, we obtain the same solution for $`a_1,a_2,\mathrm{}`$. So the operator (91) possesses at least a double algebraization, but we have not attempted to construct the other possibly existing ones yet.

The 3$`\times `$3-matrix QES potential (91) strongly contrasts with the 2$`\times `$2 ones obtained above. While the former depends on three free parameters but exists only on the full line (since $`k=1`$), the latter can have an arbitrary period (equal to $`4K(k)`$) but has only one free parameter (for instance noted by $`b`$) for a fixed $`k`$. The fact that the operator $`H(z)`$ for a 3$`\times `$3-matrix QES potential is quasi exactly solvable only for $`k^2=1`$ was unexpected to us. Accordingly, this operator corresponds to a $`3\times 3`$ matrix version of the Pöschl-Teller operator.

## 5 Concluding remarks

In the first part of this paper we obtained the matrix generalization of the celebrated sextic QES anharmonic oscillator. The examples of operators presented in the second part reflect the difficulty to classify the coupled-channel (or matrix) QES Schrödinger equations when the change of variable involves several singular points like in (26). We hope that this note will motivate further investigation of this problem.

## 6 Appendix A

For $`k=0`$ we have $`K(0)=\frac{\pi }{2}`$ and the Jacobi functions reduce to standard trigonometric functions : $`\mathrm{sn}(z,0)=\mathrm{sin}(z)`$, $`\mathrm{cn}(z,0)=\mathrm{cos}(z)`$, $`\mathrm{dn}(z,0)=1`$. For $`k=1`$, we have $`K(1)=\mathrm{}`$ and the Jacobi functions reduce to elementary functions :

$$\mathrm{sn}(z,1)=\mathrm{tanh}(z)$$
(102)
$$\mathrm{cn}(z,1)=\mathrm{dn}(z,1)=\frac{1}{\mathrm{cosh}(z)}$$
(103)
In this limit, the Lamé equation becomes a Poschl-Teller equation. For generic values of $`k`$, the Jacobi functions obey the following relations :

$$\mathrm{cn}^2+\mathrm{sn}^2=1,\mathrm{dn}^2+k^2\mathrm{sn}^2=1$$
(104)
$$\frac{d}{dz}\mathrm{sn}=\mathrm{cn}\mathrm{dn},\frac{d}{dz}\mathrm{cn}=\mathrm{sn}\mathrm{dn},\frac{d}{dz}\mathrm{dn}=k^2\mathrm{sn}\mathrm{cn}$$
(105)
These identities as well as the following ones are useful to establish the equations in the variable $`x=\mathrm{sn}^2`$ after the prefactor including the Jacobi functions has been extracted :

$`f`$ $`f^{\prime \prime }/f`$ $`(\mathrm{sn}\mathrm{cn}\mathrm{dn})f^{}/f`$
$`1`$ $`0`$ $`0`$
$`\mathrm{sn}`$ $`2k^2x(1+k^2)`$ $`k^2x^2(1+k^2)x+1`$
$`\mathrm{cn}`$ $`2k^2x1`$ $`k^2x^2x`$
$`\mathrm{dn}`$ $`2k^2xk^2`$ $`k^2x^2k^2x`$
$`\mathrm{cn}\mathrm{dn}`$ $`6k^2x(1+k^2)`$ $`2k^2x^2(1+k^2)x`$
$`\mathrm{sn}\mathrm{dn}`$ $`6k^2x(1+4k^2)`$ $`2k^2x^2(1+2k^2)x+1`$
$`\mathrm{sn}\mathrm{cn}`$ $`6k^2x(4+k^2)`$ $`2k^2x^2(2+k^2)x+1`$
$`\mathrm{sn}\mathrm{cn}\mathrm{dn}`$ $`12k^2x4(1+k^2)`$ $`3k^2x^22(1+k^2)x+1`$
No similar identities are available (to our knowledge) with different choices of the function $`f`$.


text:
# Coherent multiple Andreev reflections and current resonances in SNS junctions

## I Introduction

Transport properties of small conducting structures are strongly influenced by size effects. Oscillation of magnetoresistance in thin metallic films, and quantization of conductance in narrow wires and point contacts are examples of such effects. Size effects in superconducting tunneling have attracted attention since early experiments by Tomasch . In these experiments, oscillations of the tunnel conductance as a function of applied voltage were found for tunneling from a superconductor to a thin superconducting film of an NS proximity bilayer. The geometric resonance nature of the effect was clearly indicated by the dependence of the period of oscillations on the thickness of the superconducting film. Similar conductance oscillations for tunneling into a normal metal film of NS bilayers were reported by Rowell and McMillan . Later on an even more pronounced effect – steps on the current-voltage characteristics of SINS junctions at applied subgap voltages, $`eV<2\mathrm{\Delta }`$ – was observed by Rowell (for a review see Ref. ). In addition to the dependence on the thickness of the N-film, the period of the current steps also shows temperature dependence which scales with the temperature dependence of the superconducting gap $`\mathrm{\Delta }(T)`$. The current steps occur at applied subgap voltages, $`eV<2\mathrm{\Delta }`$, and they are understood as resonant features due to quasiparticle tunneling through superconducting bound states existing in INS wells at energies lying within the superconducting gap, $`|E|<\mathrm{\Delta }`$, de Gennes—Saint-James levels.

Recently, properties of superconducting bound states have attracted new attention in connection with observations of conductance anomalies in mesoscopic NS structures. Resonant oscillations of the subgap conductance in mesoscopic quasi-ballistic NS junctions have been reported by Morpurgo et al. The zero bias conductance peak, a related resonant phenomenon in diffusive NS junctions, was discovered even earlier and then intensively studied both theoretically and experimentally (for a review see Ref. ).

In these recent studies of mesoscopic junctions, attention was switched from the single-particle current through superconducting bound states to the two-particle (Andreev) current. The traditional view of subgap current transport in proximity NINS and SINS structures considers single-particle tunneling into bound states in the normal region of the INS well, implicitly assuming that the normal region of the INS well plays the role of equilibrium reservoir. Such a model is appropriate for low-transmission tunnel junctions with low tunneling rate compared to the inelastic relaxation rate. However, transparent mesoscopic structures are in a different transport regime where the bound levels are well decoupled from the superconducting reservoirs, and where injected quasiparticles escape from the INS well via Andreev reflection. Resonant two-particle current in quantum NINS junctions has been theoretically studied in Refs.. In SNS junctions, the quasiparticles may undergo multiple Andreev reflections (MAR) before they escape into the reservoirs. In a number of recent experiments with ballistic SNS devices fabricated with high mobility 2D electron gas (2DEG)  the highly coherent MAR transport regime has been realized.

The purpose of this paper is to develop a theory of coherent multiple Andreev reflections which will be applicable to 2DEG-based SNS junctions. In 2DEG devices the separation of the superconducting electrodes $`L`$ is typically larger than 200 nm, which is of the same order of magnitude as the superconducting coherence length, $`\xi _0=\mathrm{}v_F/\mathrm{\Delta }`$ ($`v_F`$ is the Fermi velocity of the 2D electrons), and superconducting bound states are formed well inside the energy gap. The presence of bound states in the junctions of finite length gives rise to resonances in the MAR transport, which dramatically affects the subgap current. Furthermore, it is possible in 2DEG devices to use electrostatic gates to reach the quantum transport regime with a small number of electron modes and variable transmissivity. A theory of coherent MAR has earlier been developed for short superconducting junctions, $`L\xi _0`$, where superconducting bound states do not play any significant role. Such a theory is consistent with the physical situation in atomic-size superconducting point contacts,  where quantization of conduction modes has turned out to be very helpful for detailed comparison between theory and experiment. The purpose of our study is to investigate the interplay between superconducting bound state resonances and coherent MAR in long SNS junctions, $`L\xi _0`$, in the quantum transport regime.

In a number of publications, the coherent MAR approach has been applied to long SNS junctions. However, these studies were restricted to fully transparent junctions where the bound states are strongly washed out and the resonances are not pronounced (in fact, as we will show, at zero temperature the current in such junctions does not show any structures). We will study junctions with arbitrary transmissivity, $`0<D<1`$ and pay special attention to the low-transparency limit, $`D1`$, where the resonance effects are most pronounced.

The paper is organized as follows. In Sec. II we derive a 1D model for a gated ballistic 2DEG device with one transport mode. In Sec. III we construct a scheme for calculating MAR amplitudes in terms of wave propagation in energy space. In Sec. IV, single current resonances are studied, and Sec. V is devoted to a discussion of the interplay between resonances in multi-particle currents. The properties of the total subgap current is discussed in Sec. VI.

## II 1D model for quantum SNS junctions

We consider an SNS junction similar to the one discussed by Takayanagi et al. schematically shown in Fig. 1. The junction consists of a normal conducting channel fabricated with a high mobility 2DEG, which is confined between superconducting electrodes. The distance between the electrodes is comparable to the superconducting coherence length and small compared to the elastic and inelastic mean free paths and to the normal electron dephasing length. The superconductor-2DEG interfaces are highly transmissive, the transmission coefficient typically exceeding a value 0.75, and the number of conducting modes in the 2DEG channel is controlled by a split gate.

Under these conditions, electrons ballistically move from one electrode to the other while occasionally being scattered by rare impurities or junction interfaces. Under a voltage bias applied to the junction, the transport regime corresponds to fully coherent multiple Andreev reflections (MAR). To calculate the dc current we will apply the scattering theory approach  generalized for superconducting junctions; see Refs.  and references therein.

The normal electron propagation through the junction is generally described by the $`N`$-channel scattering matrix. By assuming the split gate to select only one transport mode, we will characterize the transport through this mode by energy-dependent transmission amplitude, $`d(E)`$, and reflection amplitudes, $`r(E),r^{}(E)`$ (the energy $`E`$ is counted from the Fermi energy). The scattering amplitudes satisfying the unitarity relations, $`dr^{}+d^{}r^{}=0`$, $`|d|^2+|r|^2=1`$. The energy dispersion of the scattering amplitudes will introduce the normal electron (Breit-Wigner) and superconducting (Andreev) resonances in the scattering problem. The effect of narrow Breit-Wigner resonances on coherent MAR was earlier studied by Johansson et al. and Levi Yeyati et al. Here we will focus on the effect of Andreev resonances and only consider Breit-Wigner resonances which are wide on the scale of the energy gap. This will allow us to neglect the energy dispersion of the junction transparency, $`D=|d|^2`$const. However, the scattering phases may depend on the energy, which yields the Andreev resonances. In the simplest case, this dependence is a linear function within the energy interval $`|E|\mathrm{\Delta }`$, and we will write it on the form

$$d(E)=d_0e^{iaE},r(E)=r_0e^{ibE},$$
(1)
where $`a,b`$ are constant. In this case, the scattering properties of the normal channel are similar to those of a 1D NIN junction. Indeed, the corresponding 1D transfer matrix,

$$\widehat{T}(E)=\left(\begin{array}{cc}e^{iaE}/d_0\hfill & r_0^{}e^{i(ab)E}/d_0^{}\hfill \\ r_0e^{i(ab)E}/d_0\hfill & e^{iaE}/d_0^{}\hfill \end{array}\right),$$
(2)
can be decomposed into a product of three transfer matrices,

$$\widehat{T}(E)=e^{i\sigma _zL_lE/\mathrm{\Delta }\xi _0}\widehat{T}(0)e^{i\sigma _zL_rE/\mathrm{\Delta }\xi _0}e^{i\sigma _zk(E)L_l}\widehat{T}e^{i\sigma _zk(E)L_r},$$
(3)
where $`\sigma _z`$ is a Pauli matrix. The first and the last matrices describe ballistic propagation of an electron, with wave vector $`k(E)=\sqrt{2m(E_F+E)}/\mathrm{}k_F+E/\mathrm{\Delta }\xi _0`$, through the right and left N-regions of an effective junction with lengths $`L_r=b\mathrm{\Delta }\xi _0/2`$ and $`L_l=(ab/2)\mathrm{\Delta }\xi _0`$ respectively (from right to left), and the matrix $`\widehat{T}=e^{i\sigma _zk_FL_l}\widehat{T}(0)e^{i\sigma _zk_FL_r}`$ describes an effective barrier (I).

Quasiparticle propagation through the effective 1D SNINS junction is described by means of the time-dependent Bogoliubov-de Gennes (BdG) equation,

$$\left(\begin{array}{cc}H_0(x)& \mathrm{\Delta }(x)e^{\text{sgn}(x)ieVt/\mathrm{}}\\ \mathrm{\Delta }(x)e^{\text{sgn}(x)ieVt/\mathrm{}}& H_0(x)\end{array}\right)\left(\begin{array}{c}u(x,t)\\ v(x,t)\end{array}\right)=i\mathrm{}_t\left(\begin{array}{c}u(x,t)\\ v(x,t)\end{array}\right),$$
(4)
where $`H_0=\widehat{p}^2/2m\mu +U(x)\text{sgn}(x)eV/2`$ is the normal electron Hamiltonian, $`U(x)`$ is the impurity potential and $`V`$ is the applied voltage. The superconducting order parameter $`\mathrm{\Delta }(x)`$ is constant within the superconducting electrodes and zero within the normal region $`L_l<x<L_r`$(see Fig. 2). In the further calculations, the impurity potential is described by the transfer matrix $`\widehat{T}`$ in Eq. (3). The spatial distribution of the applied potential along the channel is modeled with a step-like function $`\pm eV/2`$. In fact, the actual spatial distribution of the potential does not play any role in this system: it can be included in the transfer matrix in Eq. (3), leading to an additional energy-independent shift in the scattering phases in the matrix $`\widehat{T}`$. As we will see later \[comment after Eq. (17)\], the energy-independent phases in the $`\widehat{T}`$-matrix do not affect the current, and can therefore be excluded.

The phase difference between the two superconductors follows from the Josephson relation ($`\dot{\varphi }=2eV/\mathrm{}`$) and introduces time dependence into the problem. The superconducting electrodes are considered to be equilibrium reservoirs where the quasiparticle wave function is a superposition of electron- and hole-like plane waves,

$$e^{\pm i\stackrel{~}{k}^exi(E\pm \sigma _zeV/2)t/\mathrm{}}\begin{array}{cc}\left(\begin{array}{c}u\\ v\end{array}\right),& e^{\pm i\stackrel{~}{k}^hxi(E\pm \sigma _zeV/2)t/\mathrm{}}\left(\begin{array}{c}v\\ u\end{array}\right)\end{array}.$$
(5)
In this equation, the $`\pm `$ signs in the time-dependent factors refer to the left/right electrode, $`\stackrel{~}{k}^{e,h}(E)`$ is the wave vector of electron/hole-like quasiparticles, and $`u(E)`$, $`v(E)`$ are the Bogoliubov amplitudes. The ratio of the Bogoliubov amplitudes equals the amplitude of Andreev reflection for particles incoming from the neighboring normal region,

$$\frac{v}{u}=a(E)=\{\begin{array}{cc}\left(Esgn(E)\sqrt{E^2\mathrm{\Delta }^2}\right)/\mathrm{\Delta },\hfill & \hfill |E|>\mathrm{\Delta }\\ \left(Ei\sqrt{\mathrm{\Delta }^2E^2}\right)/\mathrm{\Delta },\hfill & \hfill |E|<\mathrm{\Delta }\end{array}.$$
(6)
Since the time dependencies of the wave functions in the two reservoirs are different, the quasiparticle scattering by the junction is inelastic, and one has to consider a superposition of plane waves with different energies in order to construct scattering states.

## III Calculation of current using scattering states

### A Recursion relations for MAR amplitudes

We will now proceed with the construction of recurrences for the scattering amplitudes following the method suggested by Johansson et al. . To this end we introduce the wave functions in the left/right normal region ($`l/r`$) of the junction with respect to the position of the impurity. A particular scattering state, labeled with the energy $`E`$ of the incoming quasiparticle, will consist of a superposition of plane waves with energies $`E_n=E+neV`$, where $`n`$ is an integer, $`\mathrm{}<n<\mathrm{}`$,

$$\mathrm{\Psi }_l(E)=\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}e^{i\left(E_n+\sigma _zeV/2\right)t/\mathrm{}}\left[\left(\begin{array}{c}c_{n+}^{,l}e^{ik_n^ex}+c_{n+}^{,l}e^{ik_n^ex}\\ 0\end{array}\right)+\left(\begin{array}{c}0\\ c_n^{,l}e^{ik_n^hx}+c_n^{,l}e^{ik_n^hx}\end{array}\right)\right]$$
(7)
$$\mathrm{\Psi }_r(E)=\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}e^{i\left(E_n\sigma _zeV/2\right)t/\mathrm{}}\left[\left(\begin{array}{c}c_n^{,r}e^{ik_n^ex}+c_n^{,r}e^{ik_n^ex}\\ 0\end{array}\right)+\left(\begin{array}{c}0\\ c_{n+}^{,r}e^{ik_n^hx}+c_{n+}^{,r}e^{ik_n^hx}\end{array}\right)\right].$$

The normal electron/hole wave vector $`k_n^{e,h}`$ is here defined as $`k_n^{e,h}=k(\pm E_n)`$, $`k(E)=\sqrt{2m(E_F+E)}/\mathrm{}k_F+E/\mathrm{}v_F`$. The meaning of the labels for the scattering (MAR) amplitudes $`c_n`$ will be explained below.

Continuity of the scattering state wave function across the left and right NS interfaces determines the relation between the electron and hole amplitudes in the vicinity of each interface,

$$c_{n+}^{}=a_nc_n^{},c_{n+}^{}=a_n^1c_n^{},n0,a_n=a(E_n),$$
(8)
which describes elastic Andreev reflection (l/r indices are omitted). It is convenient to consider scattering amplitudes near the impurity (at $`x=\pm 0`$) rather than at the NS interfaces, and to rewrite Eq. (8) for such amplitudes, combining the amplitudes of the ballistic propagation through the normal regions with the Andreev reflection amplitude. Then, in vector notation,

$$\widehat{c}_{n\pm }=\left(\begin{array}{c}c_{n\pm }^{}\\ c_{n\pm }^{}\end{array}\right),$$
(9)
the modified relation (8) takes the form

$$\widehat{c}_{n+}=\widehat{U}_n\widehat{c}_n,n0,$$
(10)
where

$$\widehat{U}_n=e^{i\sigma _zE_nL_{l,r}/\mathrm{\Delta }\xi _0}\left(\begin{array}{cc}a_n& 0\\ 0& a_n^1\end{array}\right)e^{i\sigma _zE_nL_{l,r}/\mathrm{\Delta }\xi _0}e^{i\sigma _z\phi _n}.$$
(11)
The phase $`\phi _n=2E_nL_{l,r}/\mathrm{\Delta }\xi _0\mathrm{arccos}(E_n/\mathrm{\Delta })`$, characterizing $`\widehat{U}_n`$, is real inside the energy gap, $`|E_n|\mathrm{\Delta }`$, where it describes the total energy dependent phase shift due to ballistic propagation and Andreev reflection. Outside the gap, $`|E_n|\mathrm{\Delta }`$, the phase $`\phi _n`$ has an imaginary part which describes leakage into the superconducting reservoirs due to incomplete Andreev reflections.

By matching harmonics with the same time dependence in Eq. (7), we derive a relation between scattering amplitudes at the left and the right side of the barrier:

$$\widehat{c}_{(n+1)}^l=\widehat{T}\widehat{c}_{n+}^r,\widehat{c}_{(n+1)}^r=\widehat{T}^1\widehat{c}_{n+}^l,$$
(12)
where the effective barrier transfer matrix $`\widehat{T}`$ is defined in Eqs. (2), (3).

The recursion relations in Eqs. (10) and (12) couple the scattering amplitudes $`\widehat{c}_{n\pm }`$ into an infinitely large equation system. This equation system describing coherent MAR is illustrated by the MAR diagram in Fig. 3.

The electron part of the quasiparticle injected at the left NS interface propagates upwards along the energy axis, the amplitudes for this propagation being labeled with $`c^{}`$. At the injection energy $`E=E_0`$ (amplitude $`c_{0+}^{}`$), the quasiparticle is accelerated across the barrier (I), where the potential drops. Thus, it enters the right normal part of the junction with energy $`E_1`$ ($`c_1^{}`$), undergoes Andreev reflection and goes back as a hole ($`c_{1+}^{}`$), entering the left normal part of the junction having been accelerated to energy $`E_2`$ ($`c_2^{}`$), and is then again converted into an electron ($`c_{2+}^{}`$). The $`\pm `$ indices label the amplitudes after ($`+`$) and before ($``$) the Andreev reflection for propagation upwards along the energy axis. There is a similar trajectory of injected holes, which descends in energy, with the MAR amplitudes labeled with $`c^{}`$. Due to electron back scattering at the barrier, the upward and downward propagating waves are mixed, e.g. $`c_{0+}^{}`$ being not only forward scattered into $`c_1^{}`$, but also back-scattered into $`c_{0+}^{}`$, which opens up the possibility of interference.

Injection from the left reservoir, shown in Fig. 3, generates a MAR path which only connects even side bands at the left side of the junction with odd side bands at the right side. Injection from the right reservoir will generate a different MAR path, with even side bands at the right side of the junction, i.e. the diagram in Fig. 3 will effectively be mirrored around the barrier (I). Thus, there are two independent equation systems for the MAR amplitudes: injection from the left and from right. The $`l,r`$ labels in the MAR amplitudes can then be omitted since they are uniquely defined by the source term and the side band index.

The transport along the energy axis generated by MAR, from energy $`E`$ to $`E_n`$, is conveniently described by the effective transfer matrix $`\widehat{M}_{n0}`$,

$$\widehat{c}_n=\widehat{M}_{n0}\widehat{c}_{0+},n>0\widehat{c}_{n+}=[\widehat{M}_{0n}]^1\widehat{c}_0,n<0$$
(13)
$$\widehat{M}_{nm}=\widehat{T}_{n1}\widehat{U}_{n1}\mathrm{}\widehat{U}_{m+1}\widehat{T}_m,n>m,$$
(14)
where $`\widehat{T}_{2k}=\widehat{T}^1`$ and $`\widehat{T}_{2k+1}=\widehat{T}`$ for the injection from the left (for the injection from the right, the even and odd side band indices are interchanged). For paths within the superconducting gap, $`|E_n,E_m|<\mathrm{\Delta }`$, the matrix $`\widehat{M}_{n0}`$ satisfies the standard transfer matrix equation, $`\widehat{M}_{nm}^{}\sigma _z\widehat{M}_{nm}^{}=\sigma _z`$, which provides conservation of probability current along the energy axis,

$$j_{n\pm }^p=\widehat{c}_{n\pm }^{}\sigma _z\widehat{c}_{n\pm }^{},j_{n\pm }^p=j_{m\pm }^p,|E_{n,m}|<\mathrm{\Delta }.$$
(15)
An important consequence of the coherence of MAR is the possibility of transmission resonances in energy space. From the form of the $`\widehat{M}`$-matrix,

$$\widehat{M}_{n0}=\mathrm{}\widehat{T}^1e^{i\sigma _z\phi _k}\widehat{T}\mathrm{},$$
(16)
it is evident that when $`\phi _k=m\pi `$, the two matrices $`\widehat{T}^1`$ and $`\widehat{T}`$ will cancel each other, and the probability of transmission through this part will be unity, which leads to resonant enhancement of MAR. The solutions $`E^{(m)}`$ of the resonance equation,

$$\phi _k^{(m)}=\phi _km\pi =\frac{2E_kL_{l,r}}{\mathrm{\Delta }\xi _0}\mathrm{arccos}\frac{E_k}{\mathrm{\Delta }}m\pi =0,$$
(17)
coincide with the spectrum of the de Gennes–Saint-James levels localized in INS quantum wells . The corresponding bound states are located either on the left or the right side of the junction.

Without loss of generality, the calculations can be performed with a real matrix $`\widehat{T}`$. The transformation to such a real matrix is given by $`\widehat{T}\widehat{V}_1\widehat{T}\widehat{V}_2`$, with diagonal unitary matrices $`\widehat{V}_{1,2}`$ whose elements are constructed with the scattering phases, which are energy-independent. It is clear from Eq. (16) that since these energy-independent matrices commute with the matrices $`\widehat{U}_n`$, they cancel each other, and the matrix $`\widehat{M}_{n0}`$ undergoes a similar transformation. This will lead to an overall phase shift of the scattering state which does not affect the current.

It is interesting to consider the special case of fully transparent junctions, $`D=1`$, which has been studied in the literature . In this case, all matrices $`\widehat{T}_n`$ in Eq. (14) are equal to the unity matrix, and the $`\widehat{M}`$-matrix takes the simple form $`\widehat{M}_{n0}=\mathrm{exp}(i\sigma _z_{m=1}^{n1}\phi _m)`$. The length of the junction then enters only through the phase of the MAR amplitudes, which drops out of the side band current. Thus the dc current of fully transparent SNS junctions is independent of length and equal to the current in quantum constrictions . In particular, at zero temperature this current does not show any structures in the subgap current. It is also worth mentioning that in this particular case of fully transparent SNS junctions, the $`\widehat{M}`$-matrix is diagonal and therefore a closed set of recursive relations can be derived for the MAR probabilities (not just for the MAR amplitudes, as in the general case), equivalent to the equations for distribution functions derived in the original paper by Klapwijk, Blonder and Tinkham .

Equation (13) describes “source-free” propagation along the MAR ladder. To complete the set of equations for the MAR amplitudes we need to take into account quasi-particle injection, which introduces a source term in Eq. (13). To this end, let us consider a quasiparticle incoming from the left superconducting electrode with energy $`E`$, having a wave function of the form

$$\mathrm{\Psi }_l^S(E)=e^{iEt/\mathrm{}i\sigma _zeVt/2\mathrm{}}\left[\delta _{\nu e}e^{i\stackrel{~}{k}^ex}\left(\begin{array}{c}u\\ v\end{array}\right)+\delta _{\nu h}e^{i\stackrel{~}{k}^hx}\left(\begin{array}{c}v\\ u\end{array}\right)\right].$$
(18)
The two terms in this equation refer to electron- ($`\nu =e`$) and hole-like ($`\nu =h`$) injected quasiparticles. We now include this wave function into the continuity condition at the NS interface at energy $`E`$, which gives us the following relation between the MAR amplitudes $`\widehat{c}_{0+}`$ and $`\widehat{c}_0`$,

$$\widehat{c}_{0+}=\widehat{U}_0\widehat{c}_0+\widehat{Y}$$
(19)
$$\widehat{Y}(E)=(u^2v^2)\left(\begin{array}{c}\delta _{\nu e}/u\\ \delta _{\nu h}e^{2iEL_l/\mathrm{\Delta }\xi _0}/v\end{array}\right).$$

For quasiparticles injected from the right, a similar equation holds with the substitutions $`eh`$, and $`L_lL_r`$. Equations (13) and (19) give a complete set of equations for the MAR amplitudes with the boundary conditions $`\widehat{c}_\pm \mathrm{}=0`$ at infinity.

### B Calculation of MAR amplitudes

A formal solution of Eqs. (13) and (19), which is useful both for numerical calculations and analytical investigations, can be constructed by reducing this infinite set of recursion relations to a finite set by representing the MAR process above $`E_n`$ and below $`E_0`$ by boundary conditions involving reflection amplitudes $`r_{n+}`$ and $`r_0`$, defined as $`c_{n+}^{}=c_{n+}^{}r_{n+}^{}`$ and $`c_0^{}=c_0^{}r_0^{}`$. This gives the following representation for the vectors in Eq. (9),

$$\widehat{c}_{n+}=c_{n+}^{}\left(\begin{array}{c}1\\ r_{n+}\end{array}\right),\widehat{c}_0=c_0^{}\left(\begin{array}{c}r_0\\ 1\end{array}\right).$$
(20)
The reflection amplitudes $`r_{n+}`$ and $`r_0`$ are independent of the injection, in contrast to the coefficients $`c_{n+}^{},c_0^{}`$. Furthermore, they are determined by the boundary conditions $`\widehat{c}_\pm \mathrm{}=0`$ and can be expressed in terms of the matrix elements of $`\widehat{M}_{Nn}`$ and $`\widehat{M}_{0(N)}`$, where $`N\mathrm{}`$,

$$\underset{N\mathrm{}}{lim}\widehat{M}_{Nn}\left(\begin{array}{c}1\\ r_{n+}\end{array}\right)=0,\underset{N\mathrm{}}{lim}[\widehat{M}_{0(N)}]^1\left(\begin{array}{c}r_0\\ 1\end{array}\right)=0.$$
(21)
In other words, the vectors in Eq. (21) are equal to the asymptotical values of the eigenvectors of $`\widehat{M}`$-matrices corresponding to the eigenvalues which decrease when $`N`$ goes to infinity. The advantage of introducing the reflection amplitudes $`r_{n+}`$ and $`r_0`$ is that although they have to be calculated numerically, the recurrences which they obey do not contain resonances and converge rather quickly. This is in contrast to the matrix $`\widehat{M}_{n0}`$, which does possess resonances, but which can be calculated analytically in a straight forward way for any given $`n`$.

The solutions of the recursion equations Eqs. (10) and (12) can now be explicitly written down. For any given energy $`E`$ we get four different sets of solutions for four scattering states including electron/hole injection from the left and the right. Using the formal expression in Eq. (20) and the matrix elements of $`M_{n0}=\left(\begin{array}{cc}m_{11}& m_{12}\\ m_{21}& m_{22}\end{array}\right)`$, the solutions for injection from the left ($`n>0`$) have the form,

$$\widehat{c}_{n+}=\frac{u(1a_0^2)e^{i\phi _n}\left[\delta _{\nu e}+a_0r_0\delta _{\nu h}\right]}{m_{22}+m_{21}r_0e^{2i\phi _0}m_{12}r_{n+}e^{2i\phi _n}m_{11}r_0r_{n+}e^{2i\phi _0+2i\phi _n}}\left(\begin{array}{c}1\\ r_{n+}\end{array}\right).$$
(22)
The solutions for injection from the right can be found through interchanging $`eh`$ and calculating all quantities with respect to injection from the right. The solutions for $`n<0`$ are calculated in a similar manner.

### C Calculation of current

Now turning our attention to the current, we calculate it in the normal region next to the barrier, using the wave function in this region, $`\mathrm{\Psi }`$, and assuming quasiparticle equilibrium within the electrodes. The current then takes the form,

$$I(t)=\frac{e}{hk_F}_{\mathrm{}}^\mathrm{\Delta }𝑑E(u^2v^2)^1\underset{e/h,l/r}{}\text{Im}\left\{\mathrm{\Psi }^{}\frac{}{x}\mathrm{\Psi }^{}\right\}\mathrm{tanh}\frac{|E|}{2k_BT},$$
(23)
where $`(u^2v^2)^1=|E|/\sqrt{E^2\mathrm{\Delta }^2}=|E|/\xi `$ is the superconducting density of states, and the sum is over the four scattering states at a given energy $`E`$ associated with the electron- and hole-like quasiparticles (e/h) injected from the left and right (l/r). The current can be divided into parts with different time dependence and expressed as a sum over harmonics:

$$I(t)=\underset{N}{}I_Ne^{2iNeVt/\mathrm{}}$$
(24)
Focusing on the dc ($`N=0`$) component, and calculating the contribution of each scattering state at the injection side of the junction, we express the current spectral density $`J(E)`$ through the probability currents of electrons and holes at energies $`E_{2n}`$ (Fig. 3),

$$I_{dc}=\frac{e}{h}_{\mathrm{}}^\mathrm{\Delta }𝑑EJ(E),$$

$$J(E)=\underset{e/h,l/r}{}\frac{|E|}{\xi }\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}\left(\widehat{c}_{2n}^{}\sigma _z\widehat{c}_{2n}^{}+\widehat{c}_{2n+}^{}\sigma _z\widehat{c}_{2n+}^{}\right).$$
(25)
These currents coinside with the probability currents $`j_{n\pm }^p`$, Eq.(15), flowing along the energy axis.

It is convenient to introduce a leakage current $`J_n`$, defined as the difference of the probability currents before and after Andreev reflection,

$$J_n=\underset{e/h,l/r}{}\frac{|E|}{\xi }\left(j_n^pj_{n+}^p\right).$$
(26)
$`J_n`$ represents the amount of probability current from all the scattering states injected at energy $`E`$ and leaking out of the junction at energy $`E_n`$ (Fig. 3). The leakage current is zero inside the energy gap due to complete Andreev reflection, $`J_n=0`$, $`|E_n|<\mathrm{\Delta }`$ \[cf. Eq. (15\].

The explicit expression for the leakage current for $`n0`$ follows from Eq. (26) after insertion of Eqs. (22) and (10),

$$J_n=\underset{l/r}{}\frac{(1|a_0|^2)(1|a_n|^2)(1+|r_0a_0|^2)(1+|r_{n+}a_n|^2)}{\left|m_{22}+e^{2i\phi _0}r_0m_{21}e^{2i\phi _n}r_{n+}m_{12}e^{2i\phi _0}e^{2i\phi _n}r_0r_{n+}m_{11}\right|^2}.$$
(27)
It follows from Eq. (27) that the leakage currents are positive for all $`n0`$, $`J_n0`$. One can also show that they satisfy the inequality $`_{n0}J_n4`$, which is a consequence of the conservation of probability current: the leakage current of all side bands except of the side band $`n=0`$ does not exceed the probability current injected into four scattering states. Furthermore, the leakage current satisfies the important detailed balance equation ,

$$J_n(E)=J_n(E_n)$$
(28)
i.e. the leakage at energy $`E_n`$ due to the injection at energy $`E`$ is the same as the leakage at energy $`E`$ due to injection at energy $`E_n`$.

Using the continuity of current across the barrier, $`j_{n+}^p=j_{(n+1)}^p`$, guaranteed by the transfer matrix $`\widehat{T}`$, we can express the probability currents in Eq. (25) through the leakage current, Eq. (26),

$$\underset{e/h,l/r}{}\frac{|E|}{\xi }j_n^p=\underset{k=n}{\overset{\mathrm{}}{}}J_k,n>0$$
(29)
$$\underset{e/h,l/r}{}\frac{|E|}{\xi }j_{n+}^p=\underset{k=n}{\overset{\mathrm{}}{}}J_k,n<0,$$

by adding and subtracting consecutive terms in the sum. The spectral density of the dc charge current Eq. (25) can then be written on the form

$$J(E)=\underset{n}{}nJ_n(E),$$
(30)
since $`J_n`$ appears in $`n`$ probability currents. This formula has a clear physical meaning: the contribution to the charge current of the $`n`$-th side band is proportional to the leakage current of the side band times the effective transferred charge $`ne`$.

The detailed balance of the leakage currents, Eq. (28), allows us explicitly to prove that at zero temperature the scattering processes between (occupied) states with negative energies, $`E,E_n\mathrm{\Delta }`$ do not contribute to the current, in agreement with the Pauli exclusion principle. Indeed, by separating the contributions from side bands with $`n<0`$ and remembering that the leakage current is zero within the gap, we get for zero temperature,

$$I_{dc}=\frac{e}{h}_{\mathrm{}}^\mathrm{\Delta }𝑑E\underset{n0}{}nJ_n(E)=\underset{n>0}{}\frac{ne}{h}\left(_{\mathrm{}}^{\mathrm{\Delta }eV}𝑑EJ_n(E)+_{\mathrm{\Delta }eV}^\mathrm{\Delta }𝑑EJ_n(E)_{\mathrm{}}^\mathrm{\Delta }𝑑EJ_n(E)\right),$$
(31)
where the first and the third terms cancel each other by virtue of Eq. (28). At finite temperature, these two terms produce current of thermal excitations while the second term gives the current of real excitations created by the voltage source. Keeping only this term which dominates at low temperature, we finally get,

$$I_{dc}=\underset{n>0}{}I_n,I_n=\frac{ne}{h}\theta (neV2\mathrm{\Delta })_{\mathrm{\Delta }neV}^\mathrm{\Delta }𝑑EJ_n(E)\mathrm{tanh}(|E|/2k_BT).$$
(32)
We end this section by noting a technically useful symmetry in the current density, namely $`J_n(E)=J_n(EneV)`$, seen from the explicit form of the $`\widehat{M}_{n0}`$-matrix. This allows us to reduce the integration interval in Eq. (32) to $`neV/2<E<\mathrm{\Delta }`$.

## IV Current in terms of $`n`$-particle processes

The approach formulated above provides necessary foundations for numerical calculation of the current for arbitrary transparency and length. However, to get a full understanding of the rich subgap structure in the current-voltage characteristics, which may seem quite random, especially for intermediate transparencies and lengths (see Figs. 14-16), we will conduct a detailed analytical study of the limit of low transparency $`D1`$. The separation of currents into $`n`$-particle currents, Eq. (32), is our basis for analysis, and we will study each current $`I_n`$ separately.

As explained in the previous section, the de Gennes–Saint-James levels, Eqs. (16) and (17), are important for the current transport through the junction, leading to resonant enhancement of the current. Our main attention in this and the next section is on the calculation of the position, height and width of the main current peaks and oscillations which have the magnitude of order $`D`$. To simplify notations, left/right injection indices are omitted in most cases.

### A Single-particle current

The single-particle current, which dominates at large applied voltages, has, according to Eq. (32), an onset at $`eV=2\mathrm{\Delta }`$. The full numerical solution for the single-particle current is plotted in Fig. 4. The current shows pronounced oscillations, and the magnitude of the slope at the current onset strongly depends on the junction length.

To understand this behavior, we analyze Eq. (27) in the limit of small transparency $`D1`$, i.e. in the tunnel limit. First we note (see appendix A) that the reflection amplitudes $`r_{n+}`$ and $`r_0`$ may be expanded as

$$r_{n+}=(1)^n\sqrt{R}+O(a_{n+1}^2D)$$
(33)
$$r_0=\sqrt{R}+O(a_1^2D),$$

After inserting the explicit form of $`\widehat{M}_{10}=\widehat{T}`$ together with the expansion (33) into Eq. (27) and putting $`R=1`$, we can write the single particle current on the form

$$I_1=\frac{eD}{h}\theta (eV2\mathrm{\Delta })_{\mathrm{\Delta }eV}^\mathrm{\Delta }𝑑E\left[N^l(E)N^r(E_1)+N^r(E)N^l(E_1)\right]\mathrm{tanh}(|E|/2k_BT),$$
(34)
where

$$N^{l,r}(E)=\frac{|E|\sqrt{E^2\mathrm{\Delta }^2}}{E^2\mathrm{\Delta }^2+\mathrm{\Delta }^2\mathrm{sin}^2\left(2EL_{l,r}/\mathrm{\Delta }\xi _0\right)}.$$
(35)
In analogy with the tunnel formula for the current,  $`N^{l,r}`$ is identified as the tunneling density of states (DOS) on the left/right side of the junction. In Fig. 5 the energy dependence of the DOS is presented. The deviation of this DOS from the normal metal density of states is a manifestation of the proximity effect. The expression (35) for the DOS has earlier been derived for proximity NS sandwiches . Note that the DOS in our case is constant throughout the N-regions. In junctions with arbitrary length, the DOS usually approaches zero at the gap edge $`|E|=\mathrm{\Delta }`$ (Ref. ). Exceptions are junctions with lengths $`L_{l,r}=m\pi \xi _0/2`$ where a bound state splits off from the gap edge. In this case, the DOS diverges at the gap edge. The quantum well structure of the SNS junctions also give rise to quasi-bound states in the continuum spectrum, $`|E|>\mathrm{\Delta }`$, seen as oscillations in the DOS.

The single particle current in Eq. (34) is written as the integral over the product of the DOS at the entrance energy $`E`$ and the exit energy $`E+eV`$. The latter depends on the applied voltage, as well as the integration interval, and therefore the DOS oscillations produce oscillations of the current $`I_1`$ as a function of voltage (Rowell-McMillan oscillations ). The oscillations become more pronounced when the junction is sufficiently long, and the differential conductance may even become negative. It is also clear that the DOS oscillates as a function of the length of the junction, which give rise to oscillations also in $`I_1`$.

In short junctions, $`L\xi _0`$, the current onset at $`eV=2\mathrm{\Delta }`$ is very steep, see Fig. 4. In junctions with finite length, the current onset is smeared and replaced with a smooth oscillating behavior. This can be directly related to the smearing of the singularity in the DOS at the gap edge. The length where the crossover between these two behaviors occurs can be taken as a measure of when finite length effects become important. To estimate this length, we write equation (27) for small lengths $`L\xi _0`$, near the threshold, $`eV=2\mathrm{\Delta }+\mathrm{\Omega }`$, $`\mathrm{\Omega }\mathrm{\Delta }`$, keeping the first order terms in $`D`$ in the denominator. For a symmetric junction, $`L_l=L_r=L/2`$, we get

$$I_1=\frac{e\mathrm{\Delta }\mathrm{tanh}(\mathrm{\Delta }/k_BT)}{h}_0^\pi \frac{D\mathrm{sin}^2zdz}{\left(\mathrm{sin}z+{\displaystyle \frac{D\mathrm{\Delta }}{4\mathrm{\Omega }}}\right)^2+{\displaystyle \frac{L^2\mathrm{\Delta }}{\xi _0^2\mathrm{\Omega }}}\left(1+{\displaystyle \frac{D\mathrm{\Delta }}{4\mathrm{\Omega }}}\right)}.$$
(36)
From this formula it is clear that for short junctions ($`L=0`$), the current onset has the width $`\mathrm{\Omega }\mathrm{\Delta }D/4`$. If $`L`$ is of the order of $`\xi _0\sqrt{D}/2`$, the size of the onset has substantially diminished, and there is no visible onset at $`eV=2\mathrm{\Delta }`$ when $`L\xi _0\sqrt{D}/2`$. This crossover between steep onset and smooth behavior, which happens already for quite short lengths if $`D`$ is small, can be interpreted in terms of a bound state which is situated exactly at the gap edge in short junctions ($`L=0`$), and which moves down into the gap when $`L>0`$, the effect becoming fully pronounced when the distance from the gap edge, $`\mathrm{}v_F/L`$, exceeds the dispersion of the Andreev state, $`\sqrt{D}\mathrm{\Delta }`$, in symmetric junctions.

When $`L_{l,r}`$ approaches $`\pi \xi _0/2`$, the lowest quasi-bound state in the continuum spectrum approaches the gap edge. This leads to an accumulation of the spectral weight at the gap edge and reappearance of the singularity in the DOS, which results in the reappearance of a sharp current onset at $`eV=2\mathrm{\Delta }`$, but with smaller magnitude; see Fig. 4 ($`L_{l,r}=\pi \xi _0`$).

It is of interest to note that in our calculations, based on the scattering theory approach, the bound states are not directly involved in the single-particle transport, which therefore is non-resonant and shows no subgap resonance peaks. Within the tunnel model approach the situation is qualitatively different: the DOS in Eq. (34) usually includes the contribution of the broadened bound states, and therefore the single-particle current exists and has pronounced resonant features at subgap voltages $`eV<2\mathrm{\Delta }`$. This difference results from the fact that, within the tunnel model approach, the superconducting bound states are implicitly assumed to be connected to the reservoirs (broadening due to inelastic interaction), which allows a stationary current to flow through the bound states. In contrast, within the scattering approach, the bound states are disconnected from the reservoirs and have zero intrinsic width. In this case the bound states obtain their width only due to higher order tunneling processes involving Andreev reflections, which are manifested by the resonant multi-particle currents. In practice, the relevance of the multi-particle versus single-particle mechanism of the subgap current transport is determined by physics and depends on the ratio of the corresponding dwelling and relaxation times. In this paper, the inelastic relaxation time $`\tau _i`$ which determines the width of the single-particle resonances is assumed to be much larger than the dwelling time of the most important two-particle current, $`\tau _i\mathrm{}v_F/LD`$.

### B Two-particle current

The two-particle current $`I_2`$ in quantum point contacts, ($`L\xi _0`$) is of order $`D^2`$ when $`eV<2\mathrm{\Delta }`$ and of order $`D^2\mathrm{ln}D`$ when $`eV>2\mathrm{\Delta }`$ (Ref. ). For finite length junctions, the situation is different. For the MAR paths where the energy of the Andreev reflection coincides with a bound state, the current spectral density $`j_2^p`$ is of order unity, due to resonant transmission through this state. For low transparency $`D1`$, this gives a sharp concentration of the current density around the resonant energies. In this limit, the two-particle current is well described by the sum of contributions from these resonances, and to evaluate them we examine the energy dependence of $`J_2`$ close to the resonant energies, $`E_1=E^{(m)}+\delta E`$. Let us consider the contribution to the leakage current $`\left[J_2\right]^l`$ from quasi-particles injected from the left. As shown in appendix B, in this case Eq. (27) reduces to the standard Breit-Wigner resonance form

$$\left[J_2\right]^l=\frac{\mathrm{\Gamma }_0^{(m)}\mathrm{\Gamma }_2^{(m)}}{\left({\displaystyle \frac{\delta E\delta E^{(m)}}{\mathrm{\Delta }}}\right)^2+\left({\displaystyle \frac{\mathrm{\Gamma }_0^{(m)}+\mathrm{\Gamma }_2^{(m)}}{2}}\right)^2},$$
(37)
where the tunneling rates $`\mathrm{\Gamma }_n^{(m)}`$ are given by $`\mathrm{\Gamma }_n^{(m)}=N^l(E_n)D/2\eta ^{(m)}`$, $`n=0,2`$, and

$$\eta ^{(m)}=\mathrm{\Delta }\frac{\phi }{E}|_{E=E^{(m)}}=\frac{2L_r}{\xi _0}+\frac{\mathrm{\Delta }}{\sqrt{\mathrm{\Delta }^2(E^{(m)})^2}},$$
(38)
and the position of the resonance is shifted by

$$\delta E^{(m)}=\frac{D\mathrm{\Delta }}{4\eta ^{(m)}}\text{Im}\left\{\frac{1+e^{2i\phi _0}}{1e^{2i\phi _0}}+\frac{1+e^{2i\phi _2}}{1e^{2i\phi _2}}\right\}.$$
(39)
An analogous result is valid for quasiparticles injected from the right.

After integrating over energy, the two-particle current in the resonance approximation may be written on the form

$$I_2(eV)=\underset{i=l,r}{}\underset{m0}{}\frac{2e}{h}\theta \left(eVE^{(m)}\mathrm{\Delta }\right)\frac{2\pi D\mathrm{\Delta }}{\eta ^{(m)}}\frac{N^i\left(E^{(m)}eV\right)N^i\left(E^{(m)}+eV\right)}{N^i\left(E^{(m)}eV\right)+N^i\left(E^{(m)}+eV\right)}f^{(m)}(T,V),$$
(40)
where the summation is over the positive bound level energies, $`0<E^{(m)}<\mathrm{\Delta }`$, and the DOS $`N^i`$ should be calculated at the injection side of the junction and $`f^{(m)}(T,V)=(1/2)\left[\mathrm{tanh}\left((eVE^{(m)})/2k_BT\right)+\mathrm{tanh}\left((eV+E^{(m)})/2k_BT\right)\right]`$. According to Eq. (40), the two-particle current $`I_2(eV)`$ increases in a step-like manner in the voltage region $`\mathrm{\Delta }<eV<2\mathrm{\Delta }`$. The steps occur at every voltage where a new resonant channel through a bound state opens up, at $`eV^{(m)}=\mathrm{\Delta }+E^{(m)}`$. We note that the step positions depend on temperature and approximately scale with $`\mathrm{\Delta }(T)`$. Each current step has the height of order $`D`$. As seen from Eq. (40), the contribution to the current of a particular bound state $`E^{(m)}`$ is modulated, as a function of voltage, by the oscillations of the density of states at the entrance and exit energies, $`N(E^{(m)}\pm eV)`$. In other words, the pronounced oscillations of the two-particle current seen in Fig. 6 reflect how close the entrance and exit energies $`E^{(m)}\pm eV`$ are to a quasi-bound state in the continuum. For $`eV>2\mathrm{\Delta }`$, the two-particle current $`I_2`$ oscillates around a constant value with an amplitude of oscillation decreasing as $`\mathrm{\Delta }^2/eV^2`$ for large voltages.

It is interesting to compare the resonant structures of the two-particle current with the resonant structures in NINS junctions. In NINS junctions, the resonant current steps occur at $`eV=E^{(m)}`$, and they do not have any modulation because the DOS on the normal side of the junction is constant.

The distance between the resonances and the resonance widths are proportional to the bound level spacing, and they decrease in long junctions. For sufficiently long junctions, the two-particle current may thus give the appearance to include a series of peaks, as shown on Fig. 7. In symmetric junctions, the bound state energies at both sides of the barrier will coincide, reducing the number of steps with a factor of two, and giving current steps of double height.

We will conclude this subsection by noting that the difference between the full numerical calculation of the two-particle current and the resonant approximation given in Eq. (40) is rather small already when $`D=0.1`$, as can be seen in Fig. 7.

### C Excess current

Excess current in SNS junctions, i.e. the difference between the current in superconducting junction and in the normal junction at large voltage,

$$I^{exc}=IG_NV+O(\mathrm{\Delta }/eV),$$
(41)
is commonly considered as a measure of the intensity of Andreev reflection. In tunnel SIS junctions and low-transmissive point contacts the excess current is small, $`I^{exc}D^2e\mathrm{\Delta }/\pi \mathrm{}`$, $`D1`$, while in fully transparent contacts the excess current is large, $`I^{exc}=8e\mathrm{\Delta }/3\pi \mathrm{}`$, $`D=1`$. Accordingly, one would expect large excess current in long SNS junctions due to the resonant enhancement of the two-particle current. However, the excess current is small because of a large deficiency, of order $`D`$, of the single-particle current caused by the broadening of the current onset at the threshold. As we will show, the single-particle and two-particle currents undergo a fine cancellation, yielding small net excess current of order $`D^2`$ when $`D1`$.

The excess current has contributions only from the single- and two-particle currents, since all higher order currents include at least one Andreev reflection outside the gap whose probability is of order $`(\mathrm{\Delta }/eV)^2`$. In the limit of large voltage, $`eV\mathrm{\Delta }`$, the relevant part of the current in Eq. (32) then takes the form

$$I_1=\frac{4De}{h}_{eV/2}^\mathrm{\Delta }𝑑E\frac{(1a_0^2)(1+Ra_0^2)}{1+R^2a_0^42R\text{Re}\left\{e^{2i\phi _0}\right\}}$$
(42)
$$I_2=\frac{8D^2e}{h}_{eV}^\mathrm{\Delta }𝑑E\frac{|a_1|^2}{1+R^2|a_1|^42R\text{Re}\left\{e^{2i\phi _1}\right\}}.$$

These equations are written for symmetric junctions, $`L_r=L_l=L/2`$, and for zero temperature; small Andreev reflection amplitudes $`|a(eV/2)|1`$ have been neglected in Eq. (27). The behavior of the current in Eq. (42) as a function of voltage is presented in Fig. 8 for different lengths. It is clearly seen that the limiting value of the excess current is approached much faster in finite length SNS junctions compared to point contacts ($`L=0`$). In Fig. 9 the excess current behavior with respect to the junction length is presented for different transparencies.

To analytically examine the excess current in the limit of small transparency, $`D1`$, it is convenient to start with equations (34) and (40). To first order of $`D`$ the excess current assumes the form ($`T=0`$),

$$I^{exc}=I_1^{exc}+I_2^{exc},$$
(43)
$$I_1^{exc}=\frac{4eD\mathrm{\Delta }}{h}+\frac{2eD}{h}_\mathrm{\Delta }^{\mathrm{}}\left[N^l(E)+N^r(E)2\right]𝑑E,$$

$$I_2^{exc}=\underset{l/r,m0}{}\frac{2\pi De\mathrm{\Delta }}{h\eta ^{(m)}}.$$

Let us consider the contributions to the single-particle current from the left electrode,

$$\left[I_1^{exc}\right]^l=\frac{2eD\mathrm{\Delta }}{h}+\frac{2eD}{h}_\mathrm{\Delta }^{\mathrm{}}\left[N^l(E)1\right]𝑑E.$$
(44)
Inserting $`N^l(E)`$ from Eq. (35), this equation can be transformed to the form,

$$\left[I_1^{exc}\right]^l=\frac{2eD}{h}_\mathrm{\Delta }^{\mathrm{}}\left(\frac{E\xi }{\xi ^2+\mathrm{\Delta }^2\mathrm{sin}^2(2EL_l/\mathrm{\Delta }\xi _0)}\frac{E}{\xi }\right)𝑑E=\frac{eD}{h}_{\mathrm{}}^{\mathrm{}}𝑑\xi \frac{\mathrm{sin}^2(2EL_l/\mathrm{\Delta }\xi _0)}{\xi ^2+\mathrm{sin}^2(2EL_l/\mathrm{\Delta }\xi _0)},$$
(45)
where $`\xi =\sqrt{E^2\mathrm{\Delta }^2}`$. It is now possible to analytically continue the integral in the upper half plane which will reduce the integral to a sum over the residues of the poles given by the equation $`\xi ^2+\mathrm{sin}^2(2EL_l/\mathrm{\Delta }\xi _0)`$. Comparing this equation with Eq. (17) we find that the poles coincide with the energies of the bound states in the gap. The excess current contribution from the left-injected single particle current is thus

$$\left[I_1^{exc}\right]^l=\frac{2D\pi e\mathrm{\Delta }}{h}\underset{m0}{}\frac{1}{\eta ^{(m)}}=\left[I_2^{exc}\right]^l,$$
(46)
where $`\left[I_2^{exc}\right]^l`$ is the contribution to the two-particle current from the bound state resonances at the left electrode. A similar relation is derived for current from the right electrode. Thus, there is exact cancellation of the excess single-particle and two-particle currents to first order in $`D`$.

It is interesting to note that the cancellation effect is related to the conservation of the number of states in a proximity normal metal compared to the conventional normal metal. It follows from Eq. (44) that $`I_1h/2eD\mathrm{\Delta }`$ is equal to the difference between the number of continuum states in the proximity metal and the total number of states in a conventional metal, while, on the other hand, the number of the bound states, is equal to

$$_0^\mathrm{\Delta }𝑑E\underset{m0}{}\delta (\phi (E)m\pi )=_0^\mathrm{\Delta }𝑑E\underset{m0}{}\delta (EE^{(m)})/\eta ^{(m)}=I_2h/2eD\mathrm{\Delta },$$
(47)
according to Eq. (43).

## V Interplay between resonances

For processes with several Andreev reflections ($`n3`$), the possibilities for resonances increase. Every Andreev reflection energy may coincide with a bound state energy and thus be resonant. For some specific voltages, more than one resonance is important, creating a situation of overlapping resonances, which can enhance the current giving peaks in the current-voltage characteristics at these voltages.

### A Three-particle current

The three-particle current $`I_3`$ has a non-resonant value of order $`D^3`$. However, $`I_3`$ is enhanced to order $`D^2`$ when the energy of one of the two Andreev reflections coincides with a bound state energy.

For the applied voltage equal to the difference between two bound state energies, $`eV^{(km)}=E^{(m)}E^{(k)}`$, two resonances occur simultaneously, i.e. form a resonance consisting of two overlapping single resonances; see the inset in Fig. 10. This will enhance the current to order $`D`$ close to this voltage, giving a peak in the IVC. The number of peaks is equal to the number of bound state pairs. The peaks are located in the voltage interval $`2\mathrm{\Delta }/3<eV<2\mathrm{\Delta }`$; we note that the peak positions are weakly dependent on temperature.

To evaluate the height and the width of these peaks, we study the contribution from overlapping resonances at $`E_1E^{(k)}<0`$ and at $`E_2E^{(m)}>0`$. Close to these energies, the phases $`\phi _1^{(k)}`$ and $`\phi _2^{(m)}`$, defined in Eq. (17), are close to zero, and we find the current spectral density for injection of a quasiparticle from the left (see Appendix C),

$$\left[J_3\right]^{(km),l}=\frac{D^3N^l(E)N^r(E_3)}{\left|D4\phi _1^{(k)}\phi _2^{(m)}+iD\left(\phi _1^{(k)}N^r(E_3)+\phi _2^{(m)}N^l(E)\right)\right|^2}.$$
(48)
We now expand $`\phi _1^{(k)}`$, $`\phi _2^{(m)}`$ in the deviation from perfect overlap in energy, $`\delta E=(E_1E^{(k)}+E_2E^{(m)})/2`$, and in voltage, $`\delta V=VV^{(km)}`$, and find, using $`D1`$, from Eq. (48)

$$\left[J_3(E)\right]^{(km),l}=\frac{D\mathrm{\Gamma }_0^{(k)}\mathrm{\Gamma }_3^{(m)}}{\left(\delta E_+\delta E_{}/\mathrm{\Delta }^2D/4\eta ^{(k)}\eta ^{(m)}\right)^2+\mathrm{\Lambda }^2},$$
(49)
where $`\mathrm{\Lambda }=(\mathrm{\Gamma }_3^{(m)}\delta E_++\mathrm{\Gamma }_0^{(k)}\delta E_{})/\mathrm{\Delta }`$, $`\delta E_\pm =\delta E\pm e\delta V/2`$. The energy dependence of the current in Eq. (49) has the form of two resonant peaks with width $`\mathrm{\Gamma }\mathrm{\Delta }D\mathrm{\Delta }/\eta `$ split by the energy interval $`\sqrt{D}\mathrm{\Delta }/\eta `$ at $`\delta V=0`$, the peak splitting increasing with increasing $`\delta V`$. After integration over energy, the overlapping resonances give a current contribution in the form of a current peak ($`k_BT\mathrm{\Delta }`$),

$$I_3^{(km)}(\delta V)=\frac{3De}{h}\frac{\pi \mathrm{\Delta }}{1+\eta ^{(k)}\eta ^{(m)}\left({\displaystyle \frac{e\delta V}{\sqrt{D}\mathrm{\Delta }}}\right)^2}\frac{2N^lN^r}{\eta ^{(k)}N^r+\eta ^{(m)}N^l}.$$
(50)
In this equation, the densities of states $`N^{l,r}`$ are taken at the entrance and exit energies, $`N^l\left(2E^{(k)}E^{^{(m)}}\right)`$ and $`N^r\left(2E^{(m)}E^{(k)}\right)`$, and the temperature is taken to be zero. A factor of $`2`$ has been included in Eq. (50) to take into account the similar resonant process for injection from the right, where $`E_1=E^{(m)}`$ and $`E_2=E^{(k)}`$.

The curve for the three-particle current versus voltage thus consists of peaks with heights of order $`D`$ and half-width $`\mathrm{\Gamma }_V=\sqrt{D}\mathrm{\Delta }/\eta `$ on top of a background of order $`D^2`$. The background current increases with voltage in the interval $`2\mathrm{\Delta }/3<eV<\mathrm{\Delta }`$ as more single resonances come into the integration region. In the interval $`\mathrm{\Delta }<eV`$, the background current decreases due to broadening of the resonances because of leakage associated with incomplete Andreev reflection outside the gap.

In long symmetric junctions the current peaks form an interesting triangular pattern. To see this, we first note that if the bound state spectrum were perfectly linear, several of the peaks described by Eq. (48)-(50) will be situated at the same voltage since $`eV^{(km)}=E^{(m)}E^{(k)}=E^{(m+1)}E^{(k+1)}`$ (see also the inset in Fig. 11), and thus the total number of peaks will be reduced while their respective height will be increased. Since the bound state spectrum is not linear, the peaks show splitting. However, the deviation from linearity is small and in practice the peaks form clusters, giving combined peaks with height roughly equal to the number of clustering peaks. In the interval $`2\mathrm{\Delta }/3<eV<\mathrm{\Delta }`$, the number of peaks in a cluster increases in steps of three from $`1`$ to $`4`$, etc., up to the number of bound states. In the interval $`\mathrm{\Delta }<eV<2\mathrm{\Delta }`$ the number of peaks in a cluster decreases in steps of one. This gives an appearance of a “peak triangle” for very long junctions, shown in Fig. 11.

This “peak triangle” is further enhanced by the background current, which has a similar triangular form, as explained above.

### B Four-particle current

The four-particle current has a non-resonant value of order $`D^4`$, which is enhanced to order $`D^3`$ when the energy of one of the three Andreev reflections coincides with a bound state energy. Similar to the three-particle current, overlapping resonances can enhance the magnitude of the current $`I_4`$ to the order $`D`$ for those voltages where both the first and the third Andreev reflections coincide with the bound states, as shown in the upper inset in Fig. 12. Indeed, it is clear from the explicit form of $`\widehat{M}_{40}=\widehat{T}e^{i\sigma _z\phi _3}\widehat{T}^1e^{i\sigma _z\phi _2}\widehat{T}e^{i\sigma _z\phi _1}\widehat{T}^1`$ that when $`\phi _1=k\pi `$ and $`\phi _3=m\pi `$, then $`\widehat{M}_{40}=(1)^{k+m}e^{i\sigma _z\phi _2}`$, i.e. the transparency of the MAR trajectory is enhanced to unity. Other combinations of the resonances, e.g. when the first and the second Andreev reflection occur at bound state energies, will produce peaks of order $`D^2`$ or smaller, as described in Appendix D.

Focusing on the double resonances which produce large ($`D`$) current peaks, we find that in short junctions with just one pair of bound states, $`\pm E^{(0)}`$, the double resonance will occur at voltage $`eV=E^{(0)}`$, provided the energy of the bound state is within the interval $`\mathrm{\Delta }/2E^{(0)}\mathrm{\Delta }`$. The spectral density of the current has a form similar to the one in Eq. (49), the major difference being the small peak splitting, $`D\mathrm{\Delta }/\eta `$. The height of the resulting current peak ($`k_BT\mathrm{\Delta }`$) is

$$\left[I_4\right]_{max}=\frac{\pi De\mathrm{\Delta }}{h}\frac{1\left[a(2E^{(0)})\right]^4}{1+\left[a(2E^{(0)})\right]^4},$$
(51)
where $`a(2E^{(0)})`$ is the Andreev reflection amplitude at energy $`2E^{(0)}`$.

For longer junctions, there are many possibilities to have overlapping resonances. Two bound states at one side of the junction with energies $`E^{(k)}<0`$ and $`E^{(m)}>0`$ can give a peak in $`I_4`$ if $`(E^{(m)}E^{(k)})/2=eV\mathrm{\Delta }/2`$. Although the height of all peaks is roughly proportional to $`D`$, numerically the heights (and widths) of the peaks may vary considerably depending on the position of the second Andreev reflection. If the second Andreev reflection does not occur at the energy of a bound state, the situation is similar to the one described above; see lower inset in Fig. 12. However, if a bound state is close to the energy of the second Andreev reflection, then the current spectral density $`I_4^p(E)`$ consists of the three full-transmission peaks with widths $`D\mathrm{\Delta }/\eta `$ which are split within the interval $`\sqrt{D}\mathrm{\Delta }/\eta `$ (triple resonance). The triple resonance has larger spectral weight compared to the double resonance, which results in the larger height and width of the current peak.

Rigorously speaking, a triple resonance can only occur in asymmetric junctions because it requires equal distance between neighboring resonances, while the bound state spectrum in symmetric junctions is not equidistant. However, in long junctions, the deviation from the equidistant spectrum is small, and quasi-triple resonances may therefore occur also in long symmetric junctions.

This effect can be observed in Fig. 12, where the four-particle current for a symmetric junction with length $`L=7\xi _0>2\pi \xi _0`$ consists of three peaks with different heights: the central peak corresponding to the quasi-triple resonance while the two side peaks corresponding to the double resonances with the heights given by Eq. (51).

Finally, it is worth noting that, similar to the situation for the three-particle current, the peaks will form clusters, giving a smaller number of current peaks than the number of pairs of bound states in long junctions.

### C High order currents

The studied properties of multiple resonances in three- and four-particle currents allow us to make some general conclusions about resonant behavior of the high order multi-particle currents which determine the total current at small voltage. The non-resonant magnitude of an $`n`$-particle current is of order $`D^n`$ at the threshold voltage, $`eV_n=2\mathrm{\Delta }/n`$, and therefore the total non-resonant current exponentially decreases with the applied voltage (in transparent junctions, $`D1`$, the current is exponentially small at $`eV<\mathrm{\Delta }\sqrt{1D}`$). However, multiple resonances may enhance the magnitude of the current by several orders of $`D`$. The major question of interest here concerns the maximum value of the resonant current, in particular whether it can be of order $`D`$ at arbitrary small voltage.

To obtain such large current at small voltage, it is necessary to achieve a transmission probability through a high order MAR path equal to unity, which implies that the energy of at least every other Andreev reflection must coincide with a bound state (cf. the discussion in the previous subsection). For $`n>4`$, this means that three or more bound states must be approximately equidistant in energy. Since the bound state spectrum is non-equidistant, Eq. (17), this is generally not possible if the resonances are narrow; therefore, in junctions with arbitrary geometry and small transmissivity there are no large current peaks below the voltage $`eV=\mathrm{\Delta }/2`$.

However, the possibility of a large resonant current exists for junctions with sufficiently large transparency. To find the relevant transparency, let us consider a very long symmetric junction and assume for the moment that the bound state spectrum is equidistant, $`E^{(m+1)}E^{(m)}=const`$. Then, from mapping of the $`n`$-th order MAR process on a 1D multi-barrier structure (see Fig.13), it is clear that if the applied voltage is commensurate with the level spacing, e.g., $`eV=E^{(m+1)}E^{(m)}`$, the multi-barrier structure is periodic, and full transmission is achieved leading to a current peak. This conclusion is valid also for a non-equidistant spectrum if the variation of the interlevel distance does not exceed the width of the full-transmission band. The deviation of the bound state spectrum from the best linear fit does not exceed the value 0.33$`\mathrm{\Delta }\xi _0/L`$, Fig.13. On the other hand, the width of the full-transmission energy band is $`\sqrt{D}\mathrm{\Delta }/\eta `$ for equidistant spectrum and for $`n\mathrm{}`$. Thus one should expect large current structures in long symmetric junctions with transparency $`D>0.1`$ to occur at voltages $`eV>\mathrm{\Delta }\xi _0/L`$. In junctions with smaller transparency, large current structures may appear only at $`eV>\mathrm{\Delta }/2`$, as explained before; see Fig. 16. It is also easy to see that in asymmetric junctions, where the width of the full transmission band for an equidistant spectrum is $`D\mathrm{\Delta }/\eta `$ (since the relevant resonances at one side of the junction are weakly coupled to each other through the MAR process), large resonant current at small voltage may exist if $`D>0.33`$. Our numerical investigations confirm that in symmetric junctions when $`D`$ is of the order $`10^3`$, the multiple resonances are completely blocked and current peaks are exponentially suppressed at $`eV<\mathrm{\Delta }/2`$.

## VI Summary

Adding up the contributions to the current calculated in this paper, we arrive at a rather complex form of current-voltage characteristics (IVC) at subgap voltages, as shown in Figs. 14-16. Nevertheless, the analysis of the tunnel limit allows us to classify various subgap current structures. Here we will summarize the results of this classification. As a reference system we will take a short ($`L=0`$) junction where the form of the IVC is well studied . The current structures in short junctions can be interpreted as resonant features due to quasi-bound states situated at the edges of the energy gap , the resonant conditions selecting voltages equal to the gap subharmonics, $`eV=2\mathrm{\Delta }/n`$. This subharmonic gap structure of the short junction gradually changes with increasing junction length as bound states move down into the gap, giving rise to IVC structures with steps, oscillations and peaks. Major points are:

(i) The current in the subgap region is considerably enhanced, compared to the short junction case. This effect is present as soon as the effective length $`L/\xi _0`$ is comparable to, or larger than, the square root of transparency of the junction, $`L/\xi _0\sqrt{D}`$.

(ii) The main onset of the current in short junctions at $`eV=2\mathrm{\Delta }`$ shifts downwards in voltage to the value $`eV=\mathrm{\Delta }+E^{(0)}`$ where $`E^{(0)}`$ is the energy of the bound state. This shift is caused by the resonant two-particle current giving a contribution to the total current of the order of the single-particle current.

(iii) For longer junctions, the current onset transforms into a staircase within the voltage interval $`\mathrm{\Delta }<eV<2\mathrm{\Delta }`$ with the number of steps corresponding to the number of bound states, the step positions being given by $`eV=\mathrm{\Delta }+E^{(m)}`$. This is due to the resonances in the two-particle current transported through bound states. Resonant channels open up, one by one, as the voltage increases and bound states enter the “energy window” available for two-particle processes. The current plateaus are not flat but modulated because of oscillations of the density of continuum states. The period of the modulation is roughly equal to the interlevel distance, and it decreases with the junction length. The amplitude of the modulation, on the other hand, increases with the junction length. Thus, in long junctions, the current structures take the form of a series of peaks (see Fig. 14) within the voltage interval $`\mathrm{\Delta }<eV<2\mathrm{\Delta }`$. The position of the peaks has pronounced temperature dependence, scaling with the temperature dependence of the order parameter, while the distance between peaks has a weak temperature dependence.

(iv) There is another series of the current peaks whose positions only weakly depend on temperature and are entirely determined by the bound state spectrum: $`eV=E^{(m)}E^{(k)}`$ and $`eV=(E^{(m)}E^{(k)})/2`$. These peaks are caused by the overlap of two resonances in the three- and four-particle currents and they exist in the intervals of applied voltage $`2\mathrm{\Delta }/3<eV<2\mathrm{\Delta }`$ and $`\mathrm{\Delta }/2<eV<\mathrm{\Delta }`$ respectively. The heights of these peaks are comparable with the heights of the two-particle current structures ($`D`$).

(v) At voltage smaller than $`eV=\mathrm{\Delta }/2`$ the resonant current structures generally become smaller in magnitude (at least by one order in $`D`$) if the junction transparency is sufficiently small ($`D0.1`$), and the current decays exponentially when $`eV`$ approaches zero (although for some particular junction lengths there could be huge ($`D`$) current peaks caused by multiple resonances). This qualitative difference of the IVC below and above $`eV=\mathrm{\Delta }/2`$ allows one to expect a cross over from power to exponential dependence of IVC in multichannel junctions.

(vi) In transparent junctions, all current structures will persist but become smooth; appreciable current will appear below $`eV=\mathrm{\Delta }/2`$ as soon as $`D1/3`$. The current structures completely disappear in fully transparent junctions, $`D=1`$, where the IVC does not depend on the junction length; see Fig. 16.

(vii) At voltage larger than $`2\mathrm{\Delta }`$, the current undergoes oscillations, similar to Rowell-McMillan oscillations , and the excess current is approached much faster than in short junctions. In low transparency junctions the excess current is small, $`I^{exc}D^2,D1`$, despite strong Andreev reflection and large pair current $`I_2D`$.

## Acknowledgment

We thank E. N. Bratus’, J. Lantz and T. Löfwander for discussions. Support from NFR and NUTEK (Sweden) and from NEDO (Japan) is gratefully acknowledged.

## A Approximation for $`r_0`$ and $`r_{n+}`$

In this appendix, the expansion is derived for the reflection amplitudes in Eq. (33) for a quasiparticle injected from the left. From the definition of $`r_0`$ and $`r_{(1)}`$, Eq. (20), we know

$$\widehat{c}_0=c_0^{}\left(\begin{array}{c}r_0\\ 1\end{array}\right)$$
(A1)
$$\widehat{c}_{(1)}=c_{(1)}^{}\left(\begin{array}{c}r_{(1)}\\ 1\end{array}\right).$$
(A2)
They are related as

$$\widehat{c}_0=\widehat{M}_{01}\widehat{U}_1\widehat{c}_{(1)},$$
(A3)
where $`\widehat{M}_{01}=\widehat{T}`$ and $`U_1=e^{i\sigma _z\phi _1}`$. From this relation, we find $`r_0`$ in terms of $`r_{(1)}`$ as

$$r_0=\frac{\sqrt{R}+r_{(1)}e^{2i\phi _1}}{1+\sqrt{R}r_{(1)}e^{2i\phi _1}}=\sqrt{R}\left(\frac{1+x/\sqrt{R}}{1x}\right),$$
(A4)
where $`x=(1\sqrt{R})r_{(1)}e^{2i\phi _1}/(1+r_{(1)}e^{2i\phi _1})`$. When $`|x|1`$, we can make an expansion in this parameter to get to the form

$$r_0=\sqrt{R}+D\frac{r_{(1)}e^{2i\phi _1}}{1+r_{(1)}e^{2i\phi _1}}=\sqrt{R}+O(a_1^2D).$$
(A5)
Similarly we also get

$$r_{n+}=(1)^n\sqrt{R}+D\frac{r_{(n+1)+}e^{2i\phi _{n+1}}}{1+(1)^nr_{(n+1)+}e^{2i\phi _{n+1}}}=(1)^n\sqrt{R}+O(a_{n+1}^2D).$$
(A6)
## B Resonance in two-particle current

In this appendix, we derive the resonant form of the two-particle current, Eq. (37), for a quasiparticle injected from the left. The definition of $`\widehat{M}_{20}`$ is $`\widehat{M}_{20}=\widehat{T}e^{\sigma _z\phi _1}\widehat{T}^1`$, which using the pseudo-unitarity of the transfer matrices $`\sigma _z\widehat{T}^{}\sigma _z=\widehat{T}^1`$ can be written on the form

$$\widehat{M}_{20}=\frac{2i}{\sqrt{D}}\mathrm{sin}\phi _1\widehat{T}\sigma _z+e^{i\sigma _z\phi _1}.$$
(B1)
It simplifies in the limit $`D1`$, $`|\phi _1^{(m)}|=|\phi _1m\pi |1`$ to

$$\widehat{M}_{20}=\frac{(1)^k}{D}\left(2i\phi _1^{(m)}(1+\sigma _x)+D\right).$$
(B2)
Inserting the simplified expansion of $`\widehat{M}_{20}`$ and the expansion of $`r_{2+}`$ and $`r_0`$ from Eq. (33) into Eq. (27), as well as putting $`R=1`$, the leakage current density takes the form

$$\left[J_2(E)\right]^l=\frac{{\displaystyle \frac{1|a_0|^4}{\left|1e^{2i\phi _0}\right|^2}}{\displaystyle \frac{1|a_2|^4}{\left|1e^{2i\phi _2}\right|^2}}D^2}{\left|2i\phi _1^{(m)}+{\displaystyle \frac{D}{2}}\left({\displaystyle \frac{1+e^{2i\phi _0}}{1e^{2i\phi _0}}}+{\displaystyle \frac{1+e^{2i\phi _2}}{1e^{2i\phi _2}}}\right)\right|^2}.$$
(B3)
We make an expansion of the phase $`\phi _1^{(m)}=\eta ^{(m)}(EE^{(m)})/\mathrm{\Delta }=\eta ^{(m)}\delta E/\mathrm{\Delta }`$, where

$$\eta ^{(m)}=\mathrm{\Delta }\frac{\phi }{E}|_{E=E^{(m)}}=\frac{2L_r}{\xi _0}+\frac{\mathrm{\Delta }}{\sqrt{\mathrm{\Delta }^2(E^{(m)})^2}}.$$
(B4)
The two-particle current density now takes a Breit-Wigner form

$$\left[J_2(E)\right]^l=\frac{\mathrm{\Gamma }_0^{(m)}\mathrm{\Gamma }_2^{(m)}}{\left({\displaystyle \frac{\delta E\delta E^{(m)}}{\mathrm{\Delta }}}\right)^2+\left({\displaystyle \frac{\mathrm{\Gamma }_0^{(m)}+\mathrm{\Gamma }_2^{(m)}}{2}}\right)^2}$$
(B5)
where the tunneling rates are given by $`\mathrm{\Gamma }_{0,2}^{(m)}=N^l\left(E_{0,2}\right)D/2\eta ^{(m)}`$ where

$$N^l\left(E_{0,2}\right)=\text{Re}\left\{\frac{1+e^{2i\phi _{0,2}}}{1e^{2i\phi _{0,2}}}\right\}=\frac{1|a_{0,2}|^4}{\left|1e^{2i\phi _{0,2}}\right|^2},$$
(B6)
are equal to the DOS, Eq. (35) at energy $`E_{0,2}`$. The resonance is slightly shifted from $`E^{(m)}`$ with

$$\delta E^{(m)}=\frac{D\mathrm{\Delta }}{4\eta ^{(m)}}\text{Im}\left\{\frac{1+e^{2i\phi _0}}{1e^{2i\phi _0}}+\frac{1+e^{2i\phi _2}}{1e^{2i\phi _2}}\right\}.$$
(B7)
## C Resonance in three-particle current

In this appendix, the resonant form of the three-particle current, Eq. (48), is derived. The $`\widehat{M}_{30}`$-matrix, which by definition is

$$\widehat{M}_{30}=\widehat{T}^1e^{i\sigma _z\phi _2}\widehat{T}e^{i\sigma _z\phi _1}\widehat{T}^1,$$
(C1)
can be transformed using Eq. (B1) to

$$\widehat{M}_{30}=\widehat{T}^1e^{i\sigma _z\phi _2}\widehat{T}\widehat{T}^1\widehat{T}e^{i\sigma _z\phi _1}\widehat{T}^1=\left(\frac{2i}{\sqrt{D}}\mathrm{sin}\phi _2\widehat{T}^1\sigma _z+e^{i\sigma _z\phi _2}\right)\widehat{T}^1\left(\frac{2i}{\sqrt{D}}\mathrm{sin}\phi _1\widehat{T}\sigma _z+e^{i\sigma _z\phi _1}\right),$$
(C2)
which can be written in the form

$$\widehat{M}_{30}=4\mathrm{sin}\phi _1\mathrm{sin}\phi _2\left(\widehat{1}/\sqrt{D}+\widehat{T}^1/D\right)+2i\sigma _z\mathrm{sin}\left(\phi _1+\phi _2\right)+e^{i\sigma _z\phi _2}\widehat{T}^1e^{i\sigma _z\phi _1}.$$
(C3)
It simplifies in the limit of $`D1`$, $`|\phi _1^{(k)}|=|\phi _1k\pi |1`$ and $`|\phi _2^{(m)}|=|\phi _2m\pi |1`$ to

$$\widehat{M}_{30}=\frac{(1)^{k+m}}{D^{3/2}}\left[\left(D4\phi _1^{(k)}\phi _2^{(m)}\right)(1\sigma _x)+Di\sigma _z\left(\phi _1^{(k)}(1\sigma _x)+\phi _2^{(m)}(1+\sigma _x)\right)\right].$$
(C4)
Inserting this form of the $`\widehat{M}_{30}`$-matrix and the expansion (33) for $`r_0`$ and $`r_{3+}`$ into Eq. (27), as well as putting $`R=1`$, the probability current density for injection of a quasiparticle from the left takes the form

$$\left[J_3(E)\right]^l=\frac{(1|a_0|^4)(1|a_3|^4)D^3}{\left|1e^{2i\phi _0}\right|^2\left|1e^{2i\phi _3}\right|^2\left|Q\right|^2}$$
(C5)
$$Q=\left(D4\phi _1^{(k)}\phi _2^{(m)}\right)+iD\left(\phi _1^{(k)}\frac{1+e^{2i\phi _3}}{1e^{2i\phi _3}}+\phi _2^{(m)}\frac{1+e^{2i\phi _0}}{1e^{2i\phi _0}}\right)$$

where $`D1`$ is once again used.

Since $`|\phi _1^{(k)}|1`$ and $`|\phi _2^{(m)}|1`$ and the DOS at energies $`E_{0,3}`$, Eq. (35), are equal to

$$N^l(E)=\text{Re}\left\{\frac{1+e^{2i\phi _0}}{1e^{2i\phi _0}}\right\}=\frac{1|a_0|^4}{\left|1e^{2i\phi _0}\right|^2},$$
(C6)
$$N^r(E_3)=\text{Re}\left\{\frac{1+e^{2i\phi _3}}{1e^{2i\phi _3}}\right\}=\frac{1|a_3|^4}{\left|1e^{2i\phi _3}\right|^2},$$
(C7)
we arrive at the form

$$\left[J_3(E)\right]^l=\frac{N^l(E)N^r(E_3)D^3}{\left|D4\phi _1^{(k)}\phi _2^{(m)}+iD\left(\phi _1^{(k)}N^r(E_3)+\phi _2^{(m)}N^l(E)\right)\right|^2}.$$
(C8)
## D Resonance in four-particle current

In this appendix, we discuss the structure of the resonance in the four-particle current. The matrix

$$\widehat{M}_{40}=\widehat{T}e^{i\sigma _z\phi _3}\widehat{T}^1e^{i\sigma _z\phi _2}\widehat{T}e^{i\sigma _z\phi _1}\widehat{T}^1,$$
(D1)
can be written as

$$\widehat{M}_{40}=\frac{i\sigma _z}{D^2}[8\mathrm{sin}\phi _1\mathrm{sin}\phi _2\mathrm{sin}\phi _3\sqrt{D}\widehat{T}^1+D\mathrm{sin}\phi _1\mathrm{sin}\phi _2\mathrm{sin}\phi _3+D^2\mathrm{sin}(\phi _1+\phi _3\phi _2)+$$
(D2)
$$+2D\mathrm{sin}\phi _1\mathrm{cos}(\phi _3\phi _2)\sqrt{D}\widehat{T}^1+2D\mathrm{sin}\phi _3\mathrm{cos}(\phi _1\phi _2)\sqrt{D}\widehat{T}^1]+$$

$$+\frac{1}{D^2}[4D\mathrm{sin}\phi _1\mathrm{sin}\phi _3\mathrm{cos}\phi _2+2D\mathrm{sin}\phi _3\mathrm{sin}(\phi _1\phi _2)\sqrt{D}\widehat{T}+$$

$$+2D\mathrm{sin}\phi _1\mathrm{sin}(\phi _3\phi _2)\sqrt{D}\widehat{T}^1+D^2\mathrm{cos}(\phi _1+\phi _3\phi _2)].$$

From Eq. (D2) it is clear that in general $`\widehat{M}_{40}1/D^2`$. When both $`\phi _1`$ and $`\phi _3`$ are close to a multiple of $`\pi `$, $`\widehat{M}_{40}1`$, while close to other double resonances, e.g. when $`\phi _1`$ and $`\phi _2`$ are close to a multiple of $`\pi `$, $`\widehat{M}_{40}1/D`$.