# On the theory of Josephson effect in a diffusive tunnel junction

## I Introduction

In recent years, considerable advances have been made in technology of preparing low-resistance tunnel junctions with a comparatively high barrier transmissivity (tunneling probability) $`\mathrm{\Gamma }`$. This primarily applies to controlled break-junctions  as well as systems based on 2D electron gas , whose conductivity undergoes a crossover from tunnel to metal type upon a change in the barrier parameters. The problem of calculation of the Josephson current through a junction with an arbitrary transmissivity in the ballistic regime (with the electron mean free path $`l`$ much greater than the coherence length $`\xi _0`$) was solved by many authors  on the basis of the model of a single-mode junction with current-carrying banks ensuring a rapid “spreading” of supercurrent and the equality of the order parameter modulus $`\mathrm{\Delta }`$ near the barrier to its bulk value (the “rigidity” condition for $`\mathrm{\Delta }`$).

In the 1D geometry (e.g., a planar junction or a superconducting channel with a tunnel barrier ), the problem is complicated considerably due to the change in the order parameter and the quasiparticle energy spectrum in the vicinity of the junction, which makes a contribution to the phase dependence of the current $`j(\mathrm{\Phi })`$. Antsygina and Svidzinskii  determined the corresponding corrections to $`j(\mathrm{\Phi })`$ of the order of $`\mathrm{\Gamma }^2`$ for a pure ($`l\xi _0`$) superconductor in the limit of low transmissivity $`\mathrm{\Gamma }1`$:

$$\delta j(\mathrm{\Phi })=\alpha (T)I(\mathrm{\Delta })\mathrm{\Gamma }(\mathrm{sin}\mathrm{\Phi }\frac{1}{2}\mathrm{sin}2\mathrm{\Phi }),\alpha (T)1,$$
(1)   
$$I(\mathrm{\Delta })=(\pi /4)e\nu _Fv_F\mathrm{\Gamma }\mathrm{\Delta }=I_c(\mathrm{\Delta })/\mathrm{tanh}(\mathrm{\Delta }/2T),$$
(2)   
where $`\nu _F`$ is the density of states, $`v_F`$ the Fermi velocity, and $`I_c(\mathrm{\Delta })`$ the critical current through the junction.

In a diffusive superconductor (the “dirty” limit $`l\xi _0=\sqrt{D/2\mathrm{\Delta }}`$, $`D=v_Fl/3`$ is the diffusion constant), the calculation of the Josephson current for an arbitrary $`\mathrm{\Gamma }`$ is hardly possible  even in a simple model disregarding the variation of the order parameter in the vicinity of the junction. As a matter of fact, the boundary conditions for isotropic Green’s functions $`\widehat{g}(𝒓,t_1t_2)`$ at the junction, obtained by Kupriyanov and Lukichev ,

$$l(\widehat{g}\widehat{g})_L=l(\widehat{g}\widehat{g})_R=\frac{3}{4}\frac{\mu d(\mu )}{r(\mu )}[\widehat{g}_L,\widehat{g}_R],$$
(3)   
where $`d(\mathrm{cos}\theta )`$ is the tunneling probability for an electron impinging the barrier at an angle $`\theta `$, and the subscripts $`R`$ and $`L`$ mark the value to the right and left of the barrier, are valid only within the first order in small angle-averaged transmissivity $`\mathrm{\Gamma }=\mu d(\mu )`$. Lambert et al.  proved that the derivation of the boundary conditions in the general case ($`d1`$) is reduced to an analysis of a system of nonlinear integral equations for the terms in the expansion of the averaged Green’s function $`\widehat{g}(𝒓,𝒑)=\widehat{g}(𝒓)+𝒑\widehat{𝒈}_1(𝒓)+\mathrm{}`$ over Legendre polynomials. This problem can be solved only for $`\mathrm{\Gamma }1`$ by expanding the right-hand side of Eq. (3) into a power series in $`\mathrm{\Gamma }`$, which was used in Ref.  for calculating the corrections to the Josephson current of the order of $`\mathrm{\Gamma }^2`$.

In this paper, we pay attention, first of all, to the fact that the problem of calculation of the current–phase relation for a diffusive junction in the 1D geometry has the sense only in the case of low transmissivity of the barrier. Indeed, simple estimates obtained on the basis of the well-known formula for $`j(\mathrm{\Phi })`$ in the first order in $`\mathrm{\Gamma }`$,

$$j_0(\mathrm{\Phi })=I(\mathrm{\Delta })\mathrm{tanh}(\mathrm{\Delta }/2T)\mathrm{sin}\mathrm{\Phi }$$
(4)   
(which coincides, according to the Anderson theorem, with the Ambegaokar–Baratoff result for a pure superconductor ), show that even for small $`\mathrm{\Gamma }l/\xi _01`$ the critical current through the junction becomes of the order of the bulk thermodynamic critical current $`n_sev_{sc}`$, where $`v_{sc}1/m\xi _0`$ is the critical velocity of the condensate, $`n_sm\nu _FD\mathrm{\Delta }`$ its density, $`m`$ the electron mass ($`\mathrm{}=1`$). Thus, for $`\mathrm{\Gamma }l/\xi _0`$ the tunnel junction does not play any longer the role of “weak link” with the jump of the order parameter phase $`\mathrm{\Phi }`$ and other features of a Josephson element. This follows even from the boundary conditions Eq. (3) if we use the estimate $`\widehat{g}\widehat{g}/\xi _0`$ in the vicinity of the junction, which leads to $`[\widehat{g}_L,\widehat{g}_R]\mathrm{sin}\mathrm{\Phi }\xi _0/l\mathrm{\Gamma }1`$ for $`\mathrm{\Gamma }l/\xi _0`$ . This criterion of weak link can be also formulated in terms of the conductance of the system in the normal state: the resistance of the junction must exceed the resistance of a metal layer of thickness $`\xi _0`$.

From this it follows that the parameter

$$W=(3\xi _0/4l)\mathrm{\Gamma }\mathrm{\Gamma }$$
(5)   
plays a fundamental role in the theory of Josephson effect for diffusive junction (the factor 3/4 is chosen for convenience of notation). We can attach to this parameter the meaning of the effective tunneling probability for Cooper pairs, which is higher than the conventional probability $`\mathrm{\Gamma }`$ of quasiparticle tunneling. Small values of $`W1`$ correspond to “weak link” conditions (Josephson effect); for $`W>1`$, the presence of a tunnel barrier virtually does not affect the supercurrent flow and the distribution of the order parameter in a diffusive superconductor. Moreover, we can expect that just $`W`$ and not $`\mathrm{\Gamma }`$ is a true parameter of the expansion of $`j(\mathrm{\Phi })`$ in the barrier transmissivity. Indeed, the dependence of the Josephson current on the mean free path is absent only within the main approximation in $`\mathrm{\Gamma }`$, Eq. (4) and, therefore, it must be manifested in higher-order terms of the expansion of $`j(\mathrm{\Phi })`$ in the emergence of additional dimensionless parameter $`\xi _0/l`$ in them, which vanishes at $`l\mathrm{}`$. An analysis of corrections to the current–phase dependence of Eq. (4), carried out in Sec. 4 of this article in the next order in $`W`$, confirms these considerations and proves that the corrections $`\mathrm{\Gamma }^2`$ to the Josephson current obtained in Ref.  and associated with the corrections to boundary conditions Eq. (3), are much smaller and insignificant in fact.

Another important result of the analysis of the current-carrying state of a diffusive Josephson junction is the conclusion concerning the emergence of localized states of electron excitations in the vicinity of the barrier. This phenomenon is well known for a ballistic tunnel junction  in which discrete energy levels

$$ϵ_n(\mathrm{\Phi })=\pm \mathrm{\Delta }(1d\mathrm{sin}^2\mathrm{\Phi }/2)^{1/2},$$
(6)   
associated with Andreev localization of electron excitations near the jump in the order parameter phase, split from the continuous spectrum in the current-carrying state. A similar phenomenon also takes place in a diffusive junction in which, however, isolated coherent energy levels cannot exist due to electron scattering at impurities and defects. In this case, the most adequate description of the variation of the energy spectrum of excitations is the deformation of their local density of states $`N(ϵ,𝒓)=\text{Re}u^R(ϵ,𝒓)`$ ($`u^R`$ is the diagonal component of the retarded Green’s function for the superconductor) which is assumed for brevity to be normalized to its value $`\nu _F`$ in the normal metal. In the absence of current, the density of states in a homogeneous superconductor has the standard form $`N_0(ϵ)=|ϵ|\mathrm{\Theta }(ϵ^2\mathrm{\Delta }^2)/\sqrt{ϵ^2\mathrm{\Delta }^2}`$ ($`\mathrm{\Theta }(x)`$ is the Heaviside function) with root singularities at the gap boundaries. In the current state, the momentum $`p_s`$ of the superfluid condensate plays the role of a depairing factor smoothing the singularities of $`N(ϵ)`$ and reducing the energy gap $`2ϵ_{}`$ by $`\mathrm{\Delta }ϵ_{}(p_s)(Dp_s^2)^{2/3}`$ . In the vicinity of a weak link, a similar (and main) factor of the energy gap suppression is the phase jump $`\mathrm{\Phi }`$ which leads to the formation of a “potential well” around the junction having a width of the order of $`\xi _0`$ and containing localized excitations with an energy $`|ϵ|<\mathrm{\Delta }`$ (see Sec. 3). In contrast to the ballistic case, the Josephson transport in a diffusive junction is performed not only by the states in the potential well, but by excitations within the whole energy region near the gap edge where the density of states differs significantly from the unperturbed value.

## II Equations for Green’s function of a low-transparent Josephson junction

In order to calculate the density of states and equilibrium supercurrent

$$j=\frac{e}{4}\nu _Fv_FD_{\mathrm{}}^+\mathrm{}𝑑ϵf_0(ϵ)\text{Tr}\sigma _z(\widehat{g}^R\widehat{g}^R\widehat{g}^A\widehat{g}^A)(ϵ)$$
(7)   
we must solve equations for the matrix retarded (advanced) Green’s functions $`\widehat{g}^{R,A}(𝒓,ϵ)`$ averaged over the ensemble of scatterers:

$$[\sigma _zϵ+\mathrm{\Delta }\mathrm{exp}(i\sigma _z\chi )i\sigma _y,\widehat{g}]=iD(\widehat{g}\widehat{g}),$$
(8)   
Here $`\mathrm{\Delta }`$ and $`\chi `$ are the modulus and phase of the order parameter and $`f_0(ϵ)=(1/2)(1+\mathrm{tanh}(ϵ/2T))`$ is the equilibrium distribution function.

According to the normalization condition $`\widehat{g}^2=1`$ for the Green’s function, the matrix $`\widehat{g}`$ can be presented as $`\widehat{g}=𝝈𝒖`$, where $`𝝈`$ is the vector formed by Pauli matrices. Using the well-known relations $`(𝝈𝒂)(𝝈𝒃)=𝒂𝒃+i𝝈[𝒂\times 𝒃]`$, $`[\sigma _z,𝝈]=2i[𝝈\times 𝒔]`$, where $`𝒔`$ is the unit vector of “isotopic spin” directed along the $`z`$-axis in the space of Pauli matrices, we can obtain from Eqs. (3) and (8) the following equations and the boundary conditions for the vector Green’s function $`𝒖`$:

$$ϵ[𝒔\times 𝒖]+i\mathrm{\Delta }[𝝌\times 𝒖]=(D/2)[𝒖\times 𝒖],𝒖^2=1,$$
(9)   
$$\xi _0[𝒖\times 𝒖]_{L,R}=2W[𝒖_L\times 𝒖_R],$$
(10)   
where $`𝝌=(\mathrm{sin}\chi ,\mathrm{cos}\chi ,0)`$ is the symbolic vector of the order parameter phase.

Singling out the component of the vector $`𝒖`$ along the direction $`𝒔`$: $`𝒖=𝒔u+i𝒗`$ ($`𝒗𝒔=0`$), we project Eq. (9) onto the $`(x,y)`$-plane in the space of Pauli matrices:

$$ϵ𝒗\mathrm{\Delta }u𝝌=(iD/2)(u𝒗𝒗u),u^2𝒗^2=1$$
(11)   
and introduce the unit vector $`𝝍=(\mathrm{sin}\psi ,\mathrm{cos}\psi ,0)`$ directed along $`𝒗`$: $`𝒗=𝝍v`$, where $`\psi (𝒓,ϵ)`$ is the phase of “anomalous” Green’s function $`v`$ ($`𝝍=[𝝍\times 𝒔]\psi `$). The obtained system of scalar equations is a possible representation of Usadell equations:

$$ϵv\mathrm{\Delta }u\mathrm{cos}(\psi \chi )=\frac{i}{2}D[(uvvu)uv(\psi )^2],$$
(12)   
$$\mathrm{\Delta }v\mathrm{sin}(\psi \chi )=(iD/2)(v^2\psi ),$$
(13)   
$$u^2v^2=1,$$
(14)   
and its solutions determine the supercurrent

$$j(\mathrm{\Phi })=e\nu _FD_{\mathrm{}}^+\mathrm{}𝑑ϵf_0\text{Im}(v^R)^2\psi ^R.$$
(15)   
Choosing the coordinate axis $`x`$ orthogonally to the contact plane $`x=0`$ ($`\chi (+0)=\chi (0)=\mathrm{\Phi }/2`$) and taking into account the continuity of Green’s function and antisymmetry of their derivatives, we can easily obtain from Eq. (10) the boundary conditions to Eqs. (12), (13) for $`x+0`$:

$$\xi _0(uvvu)(0)=4Wu(0)v(0)\mathrm{sin}^2\psi (0),$$
(16)   
$$\xi _0\psi (0)=2W\mathrm{sin}2\psi (0).$$
(17)   
Far away from the junction, the behavior of the order parameter and Green’s function phases is described by linear asymptotic form corresponding to the given value of current

$$\chi (+\mathrm{})=\psi (+\mathrm{})=\chi _{\mathrm{}}+2p_sx,p_s=(W/\xi _0)\mathrm{sin}\mathrm{\Phi },$$
(18)   
i.e., of the superfluid momentum $`p_s`$ whose magnitude is determined in the main approximation by the condition of equality of the current Eq. (4) through the junction to its value $`j=\pi e\nu _FDp_s\mathrm{\Delta }\mathrm{tanh}(\mathrm{\Delta }/2T)`$ in the bulk of the metal. The Green’s functions tend to their asymptotic values satisfying Eqs. (12)–(14) for $`\psi =\chi `$ and $`u=v=0`$.

Using the parametrization $`u=\mathrm{cosh}\theta ,v=\mathrm{sinh}\theta `$, which takes into account the normalization condition Eq. (14), we can put in correspondence to the vector Green’s function $`𝒖`$ the following geometrical image . The unit vector $`𝒖`$ in a normal metal is directed along the isospin axis $`z`$ (which corresponds to a purely electron or hole state of excitation of a Fermi gas), while in a superconductor this vector is deflected from the axis through an imaginary angle $`i\theta `$ and turned around it through the azimuthal angle $`\psi `$. In the spatially homogeneous case, this angle obviously coincides with the phase of the order parameter ($`\psi =\chi `$), and the scalar Green’s functions $`u`$ and $`v`$ are described by the formulas

$$u^{R,A}=\mathrm{cosh}\theta _s=\frac{ϵ}{\sqrt{(ϵ\pm i0)^2\mathrm{\Delta }^2}},v^{R,A}=\mathrm{sinh}\theta _s,$$
(19)   
where $`\pm i0`$ defines the position of singularities of the retarded (advanced) Green’s function in the complex plane $`ϵ`$, and the square root in Eq. (20) is defined so that $`u^{R,A}\pm 1`$ for $`ϵ+\mathrm{}`$.

Eqs. (12)–(14) for Green’s functions should be supplemented by the self-consistency conditions for the modulus and phase of the order parameter:

$$\mathrm{\Delta }=\lambda _{\mathrm{}}^+\mathrm{}𝑑ϵf_0\text{Re}v^R,$$
(20)   
$$_{\mathrm{}}^+\mathrm{}𝑑ϵf_0\text{Re}v^R\mathrm{sin}(\psi ^R\chi )=0,$$
(21)   
where $`\lambda `$ is the constant of superconducting interaction. Taking into account the current conservation law, Eqs. (13) and (21), it is convenient to calculate the value of current at the barrier ($`x+0`$) by expressing $`\psi (0)`$ in Eq. (15) with the help of Eq. (17) through the phase jump $`2\psi (0)`$:

$$j(\mathrm{\Phi })=\frac{e}{2}\nu _Fv_F\mathrm{\Gamma }_{\mathrm{}}^+\mathrm{}𝑑ϵf_0\text{Im}(v^R(0))^2\mathrm{sin}2\psi ^R(0),$$
(22)   
which allows us to single out explicitly the small parameter of the theory, i.e., the barrier transmissivity $`\mathrm{\Gamma }`$. It can easily be verified that in the main approximation using the unperturbed values of Green’s function of Eq. (19) and phase $`\psi (0)\chi (0)=\mathrm{\Phi }/2`$, Eq. (22) leads to the result of Eq. (4).

A simplifying factor in the case of a low transmissivity of the barrier is that the quantities $`\psi \chi `$ and $`\psi `$ proportional to the current through the junction are small (see Eqs. (17) and (13)), and hence we can omit in Eq. (12) the terms quadratic in $`W`$ and containing the phase gradients. Replacing $`\psi (0)\chi (0)=\mathrm{\Phi }/2`$ in the boundary conditions Eqs. (16), (17), to the same degree of accuracy, we obtain the equation and the boundary conditions for the parameter $`\theta `$:

$$ϵ\mathrm{sinh}\theta \mathrm{\Delta }(x)\mathrm{cosh}\theta =(iD/2)^2\theta ,$$
(23)   
$$\xi _0\theta (0)=2W\mathrm{sinh}2\theta (0)\mathrm{sin}^2\mathrm{\Phi }/2,\theta (+\mathrm{})=\theta _s.$$
(24)   
Direct application of the perturbation theory to the solution of Eq. (23) ($`\theta (x)=\theta _0+\theta _1(x)`$, $`\mathrm{\Delta }(x)=\mathrm{\Delta }_0+\mathrm{\Delta }_1(x)`$) leads to an expression for the correction $`\theta _1(x)`$ containing nonintegrable singularities at the gap boundaries, and as a consequence, to the divergence of the corresponding correction to the Josephson current Eq. (4). This is associated with the emergence of localized states of quasiparticles at a tunnel junction in the current-carrying state mentioned in Introduction and considered in the next section.

## III Localized states at a tunnel barrier

It will be proved below that the depth of the “potential well” in the vicinity of the barrier is much larger than the scale of variation of the order parameter. Consequently, it is sufficient to confine an analysis of the behavior of the density of states to the model with a constant $`\mathrm{\Delta }`$, in which Eq. (23) has a simple solution describing the attenuation of perturbations of Green’s functions at a distance $`\xi _0`$ from the barrier:

$$\mathrm{tanh}\frac{\theta (x)\theta _s}{4}=\mathrm{tanh}\frac{\theta (0)\theta _s}{4}\mathrm{exp}(k_ϵ|x|),$$
(26)   
$$k_ϵ^2=i\xi _0^2\mathrm{sinh}\theta _s,\text{Re}k_ϵ>0.$$
(27)   
The quantity $`\theta (0)`$ satisfies the boundary condition following from Eqs. (24) and (25):

$$k_ϵ\xi _0\mathrm{sinh}\frac{\theta _s\theta (0)}{2}=\gamma \mathrm{sinh}2\theta (0),\gamma =W\mathrm{sin}^2\frac{\mathrm{\Phi }}{2}1,$$
(28)   
which can be reduced to the eighth-power algebraic equation in $`z=\mathrm{exp}\theta (0)`$:

$$2z^3(zz_s)^2=i\gamma ^2(z_s^21)(z^41)^2,z_s=\mathrm{exp}\theta _s.$$
(29)   
In the general case (for an arbitrary $`ϵ`$), the solution of Eq. (27) can be obtained only numerically, but the presence of the small parameter $`\gamma `$ in (26) and (27) makes it possible to apply the perturbation theory. Far away from the spectrum boundary, we can put $`\theta (0)=\theta _s`$ on right-hand side of (26), which leads to the following expression for the correction to the density of states at the barrier:

$$N(ϵ,0)N_0(ϵ)=2\gamma \text{Re}\left(\sqrt{i\mathrm{sinh}^3\theta _s}\mathrm{sinh}2\theta _s\right),$$
(30)   
that becomes obviously inapplicable for $`|ϵ|\mathrm{\Delta }`$, where $`|\theta _s|\mathrm{}`$. In this region, we must apply the improved perturbation theory (IPT) by putting $`|z|,|z_s|1`$ for an arbitrary (not necessarily small) value of $`zz_s`$. This not only reduces the power of the general Eq. (27), but also allows us to write it in a universal form which does not contain the depairing parameter $`\gamma `$:

$$(y\sqrt{E}1)^2=iy^5,$$
(32)   
$$y=z/\beta \sqrt{2},E=\beta ^2(ϵ\mathrm{\Delta })/\mathrm{\Delta },\beta =(2/\gamma )^{1/5}1,$$
(33)   
Relations Eq. (29) show that the increase in the density of states is bounded by a quantity of the order of $`\beta W^{2/5}`$ as we approach the spectrum boundary. Thus, the range of applicability of the conventional perturbation theory, Eq. (28), is determined by the condition $`(ϵ\mathrm{\Delta })/\mathrm{\Delta }\beta ^2`$ and overlaps with the region of applicability $`(ϵ\mathrm{\Delta })/\mathrm{\Delta }1`$ of the IPT. The boundary $`ϵ_{}`$ of the spectrum (the position of the bottom of the potential well), below which the density of states vanishes, corresponds to the emergence of a purely imaginary root of Eq. (29a) at the point $`E_{}=(25/6)(2/3)^{1/5}3.842`$:

$$ϵ_{}(\mathrm{\Phi })=\mathrm{\Delta }[1C(W\mathrm{sin}^2\frac{\mathrm{\Phi }}{2})^{4/5}],C=\frac{25}{36^{1/5}}5.824.$$
(34)   
The dependence of the position of the spectrum boundary on the phase jump at the junction is illustrated by Fig. 1 in which a similar dependence of the position of the Andreev level Eq. (6) in a junction between pure superconductors is shown for comparison. It should be noted that the scale of variation of $`ϵ_{}(\mathrm{\Phi })`$ is much larger than the splitting of the Andreev level from the boundary of the continuous spectrum for the same barrier transmissivity. This is associated with the large value of the depairing parameter $`\gamma `$ in the diffusive junction as compared to the splitting parameter $`\mathrm{\Gamma }`$ of the Andreev level as well as with the large numerical value of the constant $`C`$ defining the shift of the spectrum boundary Eq. (30). Fig. 2 shows the results of numerical calculation of the density of states at the junction on the basis of the general formula Eq. (27) for different values of the depairing parameter, which show that in addition of the root singularity ($`\sqrt{ϵϵ_{}}`$) at the spectrum boundary, the quantity $`N(ϵ)`$ has a “beak-type” root singularity for $`ϵ=\mathrm{\Delta }`$. Its physical nature is associated with an infinite increase in the attenuation length $`k_ϵ^1`$ of the perturbation of Green’s function in the bulk of the metal, Eq. (25), within the vicinity of the gap boundary.

For $`ϵ_{}<ϵ<\mathrm{\Delta }`$, the density of states decreases exponentially with increasing distance from the junction (Fig. 3), which corresponds qualitatively to the image of the potential well of depth $`\mathrm{\Delta }ϵ_{}`$ and of width $`\xi _0`$ with excitations localized in it.

It is well known that the Josephson current is carried through a ballistic junction by localized excitations only and can be presented in the following form:

$$j(\mathrm{\Phi })=2e\underset{n}{}\frac{ϵ_n(\mathrm{\Phi })}{\mathrm{\Phi }}\mathrm{tanh}\frac{ϵ_n(\mathrm{\Phi })}{2T},$$
(35)   
where the index $`n`$ labels Andreev levels. At the same time, Eq. (22) for current expressed in the IPT approximation in terms of the reduced variables of Eq. (29),

$$j(\mathrm{\Phi })I(\mathrm{\Delta })\mathrm{tanh}\frac{\mathrm{\Delta }}{2T}\mathrm{sin}\mathrm{\Phi }\underset{E_{}(\mathrm{\Phi })}{\overset{\mathrm{}}{}}\frac{dE}{\pi }\text{Im}\left(y^R\right)^2=j_0(\mathrm{\Phi }),$$
(36)   
shows that the charge transfer in a diffusive junction is performed not only by the states within the potential well ($`E<0,ϵ<\mathrm{\Delta }`$), but also by the excitations with energy $`ϵ>\mathrm{\Delta }`$ in the region $`ϵ\mathrm{\Delta }\mathrm{\Delta }\beta ^2`$, where the density of states differs significantly from the unperturbed value $`N_0(ϵ)`$. It should be noted in this connection that Argaman  proposed an analog of Eq. (31) for a diffusive system, which can be obtained by the replacement of the energy $`ϵ_n(\mathrm{\Phi })`$ of Andreev levels by the local value $`ϵ(\xi ,\mathrm{\Phi },x)`$ of the excitation energy for $`x=0`$, which is adiabatically deformed by supercurrent, using instead of the discrete number $`n`$ the continuous variable

$$\xi =_{ϵ_{}\left(\mathrm{\Phi }\right)}^{ϵ(\xi ,\mathrm{\Phi },x)}𝑑ϵ^{}N(ϵ^{},\mathrm{\Phi },x)$$
(37)   
viz., the number of states with an energy smaller than $`ϵ`$ ($`\xi =\mathrm{\Theta }(ϵ^2\mathrm{\Delta }^2)\sqrt{ϵ^2\mathrm{\Delta }^2}`$ for a homogeneous superconductor) . One can assume that the contributions from the bound and delocalized states to the Josephson current are taken into account simultaneously by the formula

$$j(\mathrm{\Phi })=2e\nu _F_0^{\mathrm{}}𝑑\xi \frac{ϵ(\xi ,\mathrm{\Phi },0)}{\mathrm{\Phi }}\mathrm{tanh}\frac{ϵ(\xi ,\mathrm{\Phi },0)}{2T},$$
(38)   
which, however, leads to correct results only in the case of a homogeneous current-carrying state (where $`\chi `$ plays the role of $`\mathrm{\Phi }`$) or a wide $`SNS`$-junction (with a width $`L\xi _0`$ of the normal layer) and is inapplicable for a narrow bridge and tunnel junction. Nevertheless, the consideration of the function $`ϵ(\xi ,\mathrm{\Phi },x)`$ is useful in these cases also since this allows us to visualize the variation of the energy distribution of quasiparticle states in the vicinity of the junction (Fig. 4).

## IV Current–phase dependence for a junction in the second order in $`W`$

Although the modified perturbation theory for Green’s function in the energy representation described in the preceding section is the most physically obvious method operating with actual excitation energies, it leads to considerable formal difficulties in the calculation of corrections to the Josephson current Eq. (4). Indeed, it was shown in the previous section that the expression for $`j(\mathrm{\Phi })`$ calculated on the basis of the IPT for Green’s functions, Eq. (32), coincides with Eq. (4) since the small IPT parameter $`\beta ^2`$ cancels out as we go over to the reduced variables of Eq. (29). Thus, in order to calculate the corrections to Eq. (4) we are interested in, we must leave the approximation of Eq. (29) that describes the behavior of Green’s functions correctly only in a narrow vicinity of singularity in the density of states. For this purpose, it is convenient to use the formalism of temperature Green’s functions by going over from integration over energy in Eqs. (20)–(22) to summation over the Matsubara frequencies $`\omega _n=\pi T(2n+1),n=0,\pm 1,\pm 2,\mathrm{}`$:

$$j(\mathrm{\Phi })=\pi e\nu _Fv_F\mathrm{\Gamma }T\underset{\omega _n>0}{}\text{Re}v^2(0)\mathrm{sin}2\psi (0),$$
(39)   
$$\mathrm{\Delta }(x)=2\pi \lambda T\underset{\omega _n>0}{}\text{Im}v(x)$$
(40)   
and making the substitution $`ϵi\omega _n`$ in Eq. (23). This allows us to avoid divergences of the type of Eq. (28) in the perturbation theory which, unlike the IPT, makes it possible to take into account the coordinate dependence $`\mathrm{\Delta }(x)`$.

It is expedient to use as the main approximation in the asymptotic expansion $`\theta =\theta _0+\theta _1+\mathrm{}`$ the “adiabatic” value of Green’s function corresponding to the local value of $`\mathrm{\Delta }(x)=\mathrm{\Delta }+\mathrm{\Delta }_1(x)`$, ($`\mathrm{\Delta }_1(\mathrm{})=0`$):

$$u_0(x)=\mathrm{cosh}\theta _0(x)=\frac{\omega _n}{\stackrel{~}{\omega }_n(x)},v_0(x)=\mathrm{sinh}\theta _0(x)=\frac{\mathrm{\Delta }}{\stackrel{~}{\omega }_n(x)},$$
(41)   
where $`\stackrel{~}{\omega }_n(x)=\sqrt{\omega _{n}^{}{}_{}{}^{2}+\mathrm{\Delta }^2(x)}`$. In this case, the correction $`\theta _1(x)`$ satisfies the nonhomogeneous equation

$$^2\theta _1k_\omega ^2\theta _1=^2\theta _0,k_\omega ^2=2\stackrel{~}{\omega }_n/D$$
(42)   
with the boundary conditions $`\theta _1(+0)=2W\mathrm{sinh}2\theta _s\times \mathrm{sin}^2\mathrm{\Phi }/2`$, $`\theta _1(\mathrm{})=0`$, where $`\mathrm{cosh}\theta _s=\omega _n/\stackrel{~}{\omega }_n`$ is the value of the Green’s function far away from the junction with the unperturbed value of $`\mathrm{\Delta }`$, and $`\stackrel{~}{\omega }_n=\sqrt{\omega _{n}^{}{}_{}{}^{2}+\mathrm{\Delta }^2}`$.

The self-consistency condition for $`\mathrm{\Delta }_1(x)`$ following from Eq. (20),

$$\mathrm{\Delta }_1(q)T\underset{\omega _n>0}{}\frac{\mathrm{\Delta }^2}{\stackrel{~}{\omega }_{n}^{}{}_{}{}^{3}}=T\underset{\omega _n>0}{}\frac{\omega _n}{\stackrel{~}{\omega }_n}\text{Im}\theta _1(i\omega _n,q)$$
(43)   
completes the system of equations for determining the corrections $`\theta _1`$ and $`\mathrm{\Delta }_1`$, whose solution in the Fourier representation has the form

$$\mathrm{\Delta }_1(q)=8W\mathrm{\Delta }\frac{B(q)}{\xi _0A(q)}\mathrm{sin}^2\mathrm{\Phi }/2,$$
(45)   
$$\theta _1(i\omega _n,q)=8W\mathrm{\Delta }\frac{i\omega _n}{\stackrel{~}{\omega }_n}\frac{1}{q^2+k_\omega ^2}\frac{A(0)}{\xi _0A(q)}\mathrm{sin}^2\mathrm{\Phi }/2,$$
(46)   
$$A(q)=A(0)+q^2B(q),A(0)=2\pi T\underset{\omega _n>0}{}\frac{\mathrm{\Delta }^2}{\stackrel{~}{\omega }_{n}^{}{}_{}{}^{3}},$$
(48)   
$$B(q)=2\pi T\underset{\omega _n>0}{}\frac{\omega _n^2}{\stackrel{~}{\omega }_{n}^{}{}_{}{}^{3}}\frac{1}{q^2+k_\omega ^2}.$$
(49)   
$`(\theta _1(i\omega _n,x),\mathrm{\Delta }_1(x))={\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{dq}{2\pi }}e^{iqx}(\theta _1(i\omega _n,q),\mathrm{\Delta }_1(q)).`$

As regards the correction to the asymptotic value Eq. (18) of the phase $`\psi (x)`$ of the Green’s function, it is equal to zero in this approximation. In order to prove this, we introduce the quantity $`\phi =\psi \chi 1`$, which, according to Eq. (13), obeys the equation

$$^2\phi k_\omega ^2\phi =^2\chi _1,$$
(50)   
where $`\chi _1=\chi (x)\chi (\mathrm{})`$ is a correction to Eq. (18) localized near the junction. Taking into account the boundary condition $`\phi (0)=\chi _1(0)`$ following from Eqs. (17) and (18), we find that this equation has the simple solution $`\phi (i\omega _n,q)=q^2\chi _1(q)/(q^2+k_\omega ^2)`$ which leads, after the substitution into the self-consistency condition Eq. (21), to the homogeneous integral equation for $`\chi _1(q)`$:

$$T\underset{\omega _n>0}{}\frac{\mathrm{\Delta }}{\stackrel{~}{\omega }_n}_{\mathrm{}}^+\mathrm{}𝑑q\frac{q^2\mathrm{cos}qx}{q^2+k_\omega ^2}\chi _1(q)=0.$$
(51)   
The only nonsingular solution of Eq. (43) is $`\chi _1(q)0`$, which proves the absence of a correction to the Josephson current due to the deviation of the behavior of the phases of the order parameter and Green’s functions from the linear law Eq. (18). This result can be explained as follows. The correction $`\chi _1(x)`$ is obviously of the order of the small correction $`p_{s1}(x)`$ to the constant value $`p_s`$ of Eq. (18) in the vicinity of the junction, that ensures the conservation of the current upon a change in $`N(ϵ)`$ and $`\mathrm{\Delta }`$. Since the value of $`p_sW`$, the correction to this quantity, and hence $`\chi _1(x)`$ and $`\phi `$ have a higher order of smallness ($`W^2`$) than the corrections of the order of $`W`$ we are interested in.

Substituting Eqs. (40), (41) into Eq. (22), we obtain the required correction to the Josephson current:

$`\delta j=j(\mathrm{\Phi })j_0(\mathrm{\Phi })={\displaystyle \frac{4T}{\mathrm{\Delta }}}I(\mathrm{\Delta })\mathrm{sin}\mathrm{\Phi }{\displaystyle \underset{\omega _n>0}{}}\text{Re}\left(v^2+{\displaystyle \frac{\mathrm{\Delta }^2}{\stackrel{~}{\omega }_{n}^{}{}_{}{}^{2}}}\right)=`$

$$=I(\mathrm{\Delta })W_0Z(T)\left(\mathrm{sin}\mathrm{\Phi }\frac{1}{2}\mathrm{sin}2\mathrm{\Phi }\right),$$
(52)   
$$Z(T)=\frac{16}{\pi }\sqrt{\mathrm{\Delta }\mathrm{\Delta }_0}T\underset{\omega _n>0}{}\frac{\omega _n^2}{\stackrel{~}{\omega }_{n}^{}{}_{}{}^{4}}\underset{\mathrm{}}{\overset{+\mathrm{}}{}}\frac{dk}{k^2+\stackrel{~}{k}_\omega ^2}[1+\frac{\stackrel{~}{k}_\omega ^2B(k)}{A(k)}],$$
(53)   
where $`\stackrel{~}{k}_\omega =\stackrel{~}{\omega }_n/\mathrm{\Delta }`$, $`A(k)`$ and $`B(k)`$ are defined by Eqs. (41) upon the substitution $`k_\omega \stackrel{~}{k}_\omega `$, and $`W_0`$ and $`\mathrm{\Delta }_0`$ are the values of $`W`$ and $`\mathrm{\Delta }`$ at $`T=0`$.

At low temperatures ($`T\mathrm{\Delta }`$), the summation over $`\omega _n`$ in Eqs. (41) and (45) can be replaced by integration with respect to the continuous variable $`\omega `$:

$`A(0)=1,B(k)={\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{\mathrm{tanh}^2vdv}{k^2+\mathrm{cosh}v}}=`$

$`={\displaystyle \frac{1}{k^4}}\left({\displaystyle \frac{\pi }{2}}2\sqrt{1k^2}\mathrm{arctan}\sqrt{{\displaystyle \frac{1k^2}{1+k^2}}}k^2\right),`$

which leads to the following asymptotic value of the function $`Z(T)`$ for $`T0`$:

$$Z(T)=\frac{8}{\pi ^2}_0^{\mathrm{}}𝑑k\left[\frac{\pi k^2}{(1+k^2)^{9/4}}+\frac{2B^2(k)}{1+k^2B(k)}\right]2.178.$$
(54)   
In the vicinity of critical temperature ($`\mathrm{\Delta }T`$), the quantity $`A(0)7\zeta (3)\mathrm{\Delta }^2/4\pi ^2T^2`$ is small, and the main contribution to integral of Eq. (45) comes from the region of small wave vectors $`k\mathrm{\Delta }/T`$ corresponding to damping of perturbations at large distances of the order of $`\xi (T)(T_cT)^{1/2}`$. This allows us to replace the function $`B(k)`$ by its value $`\pi \mathrm{\Delta }/4T`$ for $`k=0`$:

$`Z(T)={\displaystyle \frac{32\sqrt{\mathrm{\Delta }\mathrm{\Delta }_0}}{\pi ^3T}}{\displaystyle \underset{n0}{}}{\displaystyle \frac{1}{(2n+1)^2}}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{B(0)dk}{A(0)+k^2B(0)}}=`$

$$=2\pi \sqrt{\frac{\pi \mathrm{\Delta }_0}{7\zeta (3)T_c}}5.099.$$
(55)   
The results of numerical calculations of the $`Z(T)`$ dependence within the entire temperature range $`0<T<T_c`$ are presented in Fig. 5.

Similarly, we can calculate by using Eqs. (40) and (41) the asymptotic values of the correction $`\mathrm{\Delta }_1(0)`$ to the unperturbed value of the order parameter at the junction:

$$\frac{\mathrm{\Delta }_1}{\mathrm{\Delta }_0}=\alpha (T)W_0\mathrm{sin}^2\frac{\mathrm{\Phi }}{2},\alpha (0)=3.037,\alpha (T_c)=5.782.$$
(56)   
The dependence of the order parameter $`\mathrm{\Delta }(0)`$ on the phase jump at the junction at $`T=0`$ presented in Fig. 1 shows that the main contribution to the energy gap suppression comes from the depairing mechanism considered in Sec. 3, and the change in the order parameter is smaller than the variation of $`ϵ_{}(\mathrm{\Phi })`$.

The structure of the phase and temperature dependences of the correction to the Josephson current of Eq. (44) in a diffusive superconductor virtually coincide with expression Eq. (1) for a junction between pure metals except the following circumstance noted in Introduction: the parameter of the expansion of $`j(\mathrm{\Phi })`$ in the transmissivity of the junction for $`l\xi _0`$ is not the tunneling probability $`\mathrm{\Gamma }`$, but a considerably larger parameter $`W`$, Eq. (5). This allows one to observe higher harmonics of the current–phase dependence in diffusive tunnel junction with a comparatively high resistance. Koops et al.  apparently reported on the first experimental results in this field.

The theory discussed above describes the current–phase dependence for a diffusive Josephson junction in the whole temperature range $`0T<T_c`$ except a narrow neighborhood of $`T_c`$ in which $`\mathrm{\Delta }/T_cW_0`$ ($`\mathrm{\Delta }/T_c\mathrm{\Gamma }`$ in a pure superconductor), and the magnitude of corrections Eqs. (44) and (1) becomes comparable with $`j_0(\mathrm{\Phi })`$, while the correction Eq. (48) to $`\mathrm{\Delta }`$ becomes of the order of its unperturbed value. This means that in the definition Eq. (5) of the parameter $`W`$ near $`T_c`$ the coherence length $`\xi _0(T)`$ describing the characteristic scale of spatial variations of Green’s function and density of states should be replaced by the characteristic length $`\xi (T)`$ of variation of the order parameter (healing length) in the Ginzburg–Landau theory, whose order of magnitude is the same as $`\xi _0`$ far away from $`T_c`$. Taking into account the results of calculations of $`j(\mathrm{\Phi })`$ for a pure superconductor in the vicinity of $`T_c`$ , we can obtain the following interpolation estimate of the effective transmissivity $`W`$ suitable for any temperatures and mean free paths:

$$W\mathrm{\Gamma }\xi (T)\left(\frac{1}{l}+\frac{1}{\xi (0)}\right).$$
(57)   
As we approach $`T_c`$, the value of $`W`$ increases infinitely, which is accompanied with a decrease in the phase jump for a given external current bounded by its critical value. Thus, in the 1D geometry for an arbitrarily large normal resistance of the junction, there exists a narrow region near $`T_c`$ in which the phase difference of the order parameter at the junction is small up to values of current of the order of the bulk critical current.

The authors are grateful to T.N. Antsygina and V.S. Shumeiko for fruitful discussions.

This research was supported by the Foundation for Fundamental Studies at the National Academy of Sciences of the Ukraine (grant No. 2.4/136).