Theoretical study of superconducting proximity–contact systems by use of the quasi–classical Green’s function

準古典的グリーン関数による超伝導近接系の研究

Yasushi Nagato Graduate School of Biosphere Sciences, Hiroshima University

January, 1994

Contents

1 Introduction 1

2 The quasi–classical Green’s function 8

2.1 The Bogoliubov–de Gennes equation and the Andreev equation . . . . . . 8 2.1.1 The Bogoliubov–de Gennes equation . . . . . . . . . . . . . . . . . 8 2.1.2 The Andreev equation . . . . . . . . . . . . . . . . . . . . . . . . . 12

2.2 The quasi–classical Green’s function for a finite triple layer system . . . . . 16 2.3 The quasi–classical Green’s function of a semi–infinite triple layer system . 22 2.4 The Green’s function and physical quantities . . . . . . . . . . . . . . . . . 25 2.5 The self–consistent pair potential . . . . . . . . . . . . . . . . . . . . . . . 29

3 Andreev reflection in the normal–normal–superconducting triple layer system 35 3.1 Andreev reflection in the normal–superconducting double infinite system . 37 3.2 Normal–normal–superconducting proximity contact system . . . . . . . . . 46

4 The local density of states of the superconducting–normal–superconducting triple layer system 59 4.1 The de Gennes–Saint-James bound state . . . . . . . . . . . . . . . . . . . 59 4.2 The superconducting–normal–superconducting semi–infinite triple layer sys-

tem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63

5 The infinite double layer system with s–wave and d–wave superconductor 73

6 Summary 95

7 Acknowledgment 98

1 Introduction

In 1911, Kameringh Onnes discovered the superconductive state in mercury at liquid helium temperature. After that, numerous superconducting materials were discovered. However, its microscopic mechanism was not realized for a long time. In 1957 Bardeen, Cooper and Schrieffer proposed their theory of superconductivity as the first successful explanation for the microscopic origin of this phenomenon[1].

It has been known that when a normal metal is connected to a superconductor, the normal metal bears superconducting properties. This phenomenon is called proximity effect because the Cooper pair in the superconductor is considered to leak into the normal metal even if the normal metal has no pairing interaction. The proximity contact has been one of the very important subjects in the study of superconductivity ( for a review, see Ref. [2]). Recent progress in technologies fabricating artificial materials has made it possible to study very clean and thin proximity-contact layers. In such a clean system, many interesting phenomena which cannot be expected in an isolated normal metal can be observed.

At the normal–superconducting ( N–S ) interfacial boundary, there occurs the so called Andreev reflection, that is, an incident electron from the N region is reflected back into a hole[5]. The N–S junction shows non–Ohmic current-voltage behavior[8][9], because the reflected hole carries away a positive charge and the current is enhanced. A theoretical attempt to explain the non-Ohmic behavior in the language of the Andreev reflection was proposed by Blonder, Tinkham and Klapwijk[10]. They have taken into account the effect of non–ideal interface but assumed a step-like pair potential. In the proximity contact system, the pair potential is in general depressed near the interface. van Son et al.[6] considered the effects by the depression on the Andreev reflection by assuming a simple analytical form for the spatial variation of the pair potential. Bruder[12] studied the Andreev scattering for a d–wave pairing superconductor, using a self-consistently solved pair potential and taking account of the interfacial reflection coefficient R. Nagato et al.[15] studied the Andreev reflection in the double infinite N–S system, and obtained the relation between the Andreev reflection coefficient and the density of states at the interface. Usually, point contact device is used in the attempt to observe the Andreev reflection through the I-V characteristics. In this system, it is expected that the voltage drop occurs at the interface between the electrode and the normal metal. It means that there is also finite reflection between the electrode and the normal metal. As a result, the multiple reflection process within the normal layer will be important and the I-V characteristics will show a geometrical resonance effect. Taking account of the geometrical resonance effect, we discuss the Andreev reflection and the differential conductance.

In a finite width normal metal connected to a superconductor, the density of states has significant structure due to the finiteness of the layer, in the energy range around the magnitude of the pair potential[28][29][30][31][32][33][34][35]. de Gennes and Saint-James first calculated the density of states in the normal metal connected to a superconductor with a spatially constant pair potential. They found that there is a state with the energy below the energy gap of the superconductor , which occurs because of the finiteness of the normal layer and the Andreev reflection at the N–S interface. This state is the so called ”de Gennes–Saint-James bound state”. The existence of the de Gennes–Saint-James bound state is the origin of the effective energy gap in the normal metal.

The density of states can be detected by a tunneling experiment, by a scanning tunneling spectroscopy(STM) and so on. In particular, the STM experiment is able to study a spatial dependence of the density of states. Recently, Inoue and Takayanagi[36] reported the STM measurement of Nb/InAs/Nb proximity contact system. Their data indicate that the local density of states in the InAs region has an effective energy gap and the gap varies spatially.

Tanaka et al.[34] studied the density of states of the S–N–S system with ideal interfaces. Hara et al.[35] studied the N–S finite system in which the N–S interface has finite reflection but the N region has no pairing interaction. We study the density of states of the normal region in the S–N–S system taking account of the effects by the finite interfacial reflection and by a pairing interaction in the normal region. The finite reflection will give rise to the geometrical resonance effect and the pairing interaction will lead to the spatial variation of the density of states.

To study the effect of the geometrical resonance on the Andreev reflection and to study the density of states in the S–N–S system, we treat the semi–infinite triple layer system as is depicted in Fig.1.

L  C  R  Z 0 L Δ(Z)

Figure 1: The semi–infinite triple layer system.

In the original BCS theory, an attractive interaction which leads to the Cooper instability is originated by the electron–phonon interaction and is an s–wave pairing interaction ( spin singlet pairing interaction ). The electron–phonon interaction may not be the only possible mechanism to obtain an attractive interaction. The possibility of unconventional pairing in superconducting metals was first discussed by Anderson and Morel[37] and Balian and Werthamer[38]. The recent discovery of ”heavy fermion superconductor” is the experimental realization of a class of superconductors with unconventional pairing. Pals et al.[39] investigated the Josephson contact between a singlet and a triplet superconductor within the tunneling Hamiltonian model. Poppe[40] observed a Josephson current between an s– wave superconductor Al and a heavy fermion superconductor CeCu_{2}Si_{2}. Ashauer et al.[41] studied a thin film of standard superconductor in proximity contact with a bulk unconventional material. It is an interesting problem to study the properties of proximity contact superconductors with different kind of pairing symmetry.

Theoretical treatments to study such proximity systems so far reported are mostly based on the ideal model. In the ideal model the pair potential in the superconductor region is assumed to be constant and the interfaces are assumed to have electron transmission coefficients of unity. The pair potential, however, cannot be constant, since it is depressed near the N–S interface due to the proximity effect. We have to treat the spatial variation on the scale of the coherence length. Moreover, the finite reflection that occurs in general at the interface of different kind metals should be taken into account. Even if the wave function can be connected smoothly at the interface, the difference of the Fermi velocity yields a finite reflection coefficient

R = |\frac{v_{F}^{N} − v_{F}^{S}}{v_{F}^{N} + v_{F}^{S}} |^{2} .  (1)

The quasi–classical Green’s function method[18][19] is able to treat the spatial variation of the pair potential and other physical observables. The conventional quasi–classical formulation[25] consists of the Eilenberger equation for the quasi–classical Green’s function ĝ, the boundary condition for ĝ and the normalization condition ĝ^{2} = −1. When one studies the proximity contact system by use of the conventional method, one has to solve the Eilenberger equation with the boundary condition and the normalization condition numerically. Since the Eilenberger equation includes a solution exploding at infinities, the numerical calculation needs sophisticated techniques. Calculations by use of the conventional quasi– classical Green’s function have been reported[12][25][26][33][41]. Kieselmann[33] obtained the self-consistent pair potential and the tunneling density of states in a N–S contact of a normal metal film with a semi–infinite superconductor. Bruder[12] studied Andreev scattering under the self–consistent pair potential in an unconventional superconductor. In the triplet superfluid ^{3}He, Kurkijärvi and Rainer [26] also studied Andreev scattering by the wall. In the finite double layer systems, however, the conventional quasi–classical technique cannot be applied[20] because the normalization condition mentioned above is valid only in bulk system.

Ashida et al.[20] have proposed a new quasi–classical formulation which can be applied to finite double layer systems. They obtained an explicit formula of the quasi–classical Green’s function which already satisfies the boundary condition. The formulation is written in a form including the evolution operators together with reflection coefficient R at the interface. This new quasi–classical Green’s function ( AAHN Green’s function ) has been applied to some systems. Ashida et al.[21] studied the transition temperature of the N–S bilayer. Hara et al.[35] studied the local density of states in the N–S bilayer. Nagato et al.[15] extended AAHN Green’s function to the system including infinite layer and discussed the Andreev reflection and the local density of states. Higashitani et al.[42] discussed the Meissner effect in the normal metal connected to the superconductor. The AAHN Green’s function is a powerful tool to study the proximity contact system. In contrast to the conventional quasi–classical formulation, the AAHN Green’s function is explicitly expressed. The explicit expression already satisfies the boundary condition. To study the proximity contact system of present interest, we extend the AAHN Green’s function to the semi–infinite proximity contact system.

This thesis is organized as follows : In Sec.2 we mainly discuss the framework of the quasi–classical Green’s function for studying the proximity contact system. We first begin with the Bogoliubov-de Gennes equation and define the Gor’kov Green’s function. Using the Andreev approximation, we derive the quasi–classical Green’s function. Following AAHN, we first obtain the quasi–classical Green’s function for the finite triple layer system. Taking the limit of the layer size to infinity, we obtain the Green’s function for the semi–infinite triple layer system. The relation between the Green’s function and the physical quantities is discussed. We also discuss on the calculation of the self–consistent pair potential. In Sec.3∼5 we apply the obtained quasi–classical Green’s function to some proximity contact systems. In Sec.3, we study the Andreev reflection for the point contact system. To explain the Andreev reflection, we first discuss it in the normal–superconducting infinite double layer. Secondly, we consider the normal–normal–superconducting model for the point contact experiment. We calculate the Andreev reflection and the differential conductance by use of the BTK formula. In Sec.4, we study the local density of states in the normal–superconducting proximity contact system. We first discuss the de Gennes–Saint-James bound state. We then obtain the local density of states in the superconducting–normal–superconducting system. In Sec.5, we discuss the supercurrent across the interface between superconductors with different pairing symmetries.

Throughout this thesis, we use the units h̄ = k_{B} = 1.

2 The quasi–classical Green’s function

In this chapter, we describe an outline of the methods to study the superconducting proximity contact system in the clean limit. We start from the BCS model Hamiltonian and derive the Bogoliubov de Gennes equation which is useful to treat an inhomogeneous superconducting system. Within the quasi–classical approximation, the Bogoliubov de Gennes equation is reduced to the Andreev equation, and a formal solution of the Andreev equation is expressed by a spatial evolution operator introduced by Ashida et al.[20]. To obtain the superconducting pair potential, the density of states and so on, we use the Green’s function method. In accordance with the quasi–classical approximation in the Bogoliubov–de Gennes equation, the same approximation to the Green’s function leads to the quasi– classical Green’s function used in this thesis. We construct the quasi–classical Green’s function for the finite triple layer system, following Ashida et al.( AAHN )[20]. The quasi– classical Green’s function can be obtained in a form including the spatial evolution operator and the reflection coefficients at the interfaces. Starting from the Green’s function of the finite width triple layer, we obtained the Green’s function in the semi–infinite triple layer system of present interest.

2.1 The Bogoliubov–de Gennes equation and the Andreev equation

2.1.1 The Bogoliubov–de Gennes equation

In the BCS theory, it is shown that a weak attractive interaction between electrons, such as that caused in second order of the electron-phonon interaction, causes an instability of the ordinary Fermi-sea ground state of the electron gas. The BCS state is characterized by an energy gap in the quasi-particle excitations of the system. The superconducting state is characterized by the energy gap ( the pair potential ). Since we treat the electrons near the Fermi surface, we choose the Fermi energy as the origin of energy. Using a creation and annihilation operator of Fermion (electron ) ψ_{α}^{†}(r), ψ_{α}(r) , the weak coupling BCS Hamiltonian is written as

H = ∑_{α} ∫ drψ_{α}^{†}(r)ξ(∇)ψ_{α}(r)
+\frac{1}{2} ∑_{α,β} ∫ drdr′v(|r − r′|)ψ_{α}^{†}(r)ψ_{β}^{†}(r′)ψ_{β}(r′)ψ_{α}(r),  (2)

ξ(∇) = −\frac{1}{2m}∇^{2} − μ,

where α, β are spin indices and μ is the chemical potential. One applies the Gor’kov approximation, i.e., mean field approximation, to Eq.(2) to obtain

H^{Gor′kov} = ∑_{α} ∫ drψ_{α}^{†}(r)ξ(∇)ψ_{α}(r)

+\frac{1}{2} ∑_{α,β} ∫ drdr′v(|r − r′|)[⟨ψ_{α}^{†}(r)ψ_{β}^{†}(r′)⟩ψ_{β}(r′)ψ_{α}(r)

+ψ_{α}^{†}(r)ψ_{β}^{†}(r′)⟨ψ_{β}(r′)ψ_{α}(r)⟩ − ⟨ψ_{α}^{†}(r)ψ_{β}^{†}(r′)⟩⟨ψ_{β}(r′)ψ_{α}(r)⟩],  (3)

The last term ⟨ ⟩⟨ ⟩ is a correction for the overcounting of the interaction. We define

Δ_{αβ}(r, r′) = v(|r − r′|)⟨ψ_{β}(r′)ψ_{α}(r)⟩,  (4)
Δ_{αβ}^{†}(r, r′) = v(|r − r′|)⟨ψ_{β}^{†}(r′)ψ_{α}^{†}(r)⟩  (5)

and obtain

H^{Gor′kov} = ∑_{α} ∫ drψ_{α}^{†}(r)ξ(∇)ψ_{α}(r)

+\frac{1}{2} ∑_{α,β} ∫ drdr′[Δ_{αβ}^{†}(r, r′)ψ_{α}(r)ψ_{β}(r′) + ψ_{α}^{†}(r)ψ_{β}^{†}(r′)Δ_{αβ}(r, r′)

−v(r − r′)⟨ψ_{α}^{†}(r)ψ_{β}^{†}(r′)⟩⟨ψ_{β}(r′)ψ_{α}(r)⟩].  (6)

The Heisenberg operators ψ_{α}(r, t) and ψ_{α}^{†}(r, t) obey equations of motion :

i∂_{t}ψ_{α}(r, t) = [ ψ_{α}(r, t), H^{Gor′kov} ]
= ∫ dr′{ξ(∇)δ(r − r′)ψ_{α}(r′, t) + ∑_{β} Δ_{αβ}(r, r′)ψ_{β}^{†}(r′, t)},  (7)
i∂_{t}ψ_{α}^{†}(r, t) = [ ψ_{α}^{†}(r, t), H^{Gor′kov} ]
= ∫ dr′{−ξ(∇)δ(r − r′)ψ_{α}^{†}(r′, t) + ∑_{β} Δ_{αβ}^{†}(r, r′)ψ_{β}(r′, t)}.  (8)

Thus we can write in a matrix form as

i∂_{t}Ψ̂(r, t) = ∫ dr′E(r, r′)Ψ̂(r′, t),  (9)

where

E(r, r′) = ( \begin{matrix}ξ(∇)δ(r − r′) & Δ(r, r′) \\ Δ^{†}(r, r′) & −ξ(∇)δ(r − r′) \end{matrix} ) ,  (10)
Ψ̂(r, t) = ( \begin{matrix}ψ_{↑}(r, t) \\ ψ_{↓}(r, t) \\ ψ_{↑}^{†}(r, t) \\ ψ_{↓}^{†}(r, t) \end{matrix} ) ,  (11)
Δ(r, r′) = ( \begin{matrix}Δ_{↑↑}(r, r′) & Δ_{↑↓}(r, r′) \\ Δ_{↓↑}(r, r′) & Δ_{↓↓}(r, r′) \end{matrix} ) , Δ^{†}(r, r′) = ( \begin{matrix}Δ_{↑↑}^{†}(r, r′) & Δ_{↑↓}^{†}(r, r′) \\ Δ_{↓↑}^{†}(r, r′) & Δ_{↓↓}^{†}(r, r′) \end{matrix} ) . (12)

In order to solve Eq.(9), we have only to consider the Bogoliubov de-Gennes equation

∫ dr′E(r, r′)Ψ(r′) = EΨ(r).  (13)

where we have regarded Ψ̂(r, t) in Eq.(9) as c-number Ψ(r, t) and have substituted Ψ(r, t) = Ψ(r) exp (−iEt) into Eq.(9) .

When one studies layered proximity contact systems where the interface and the boundaries have translational symmetry in the x, y plane, the momentum component parallel to the boundary is a conserved quantity. For a given parallel momentum p, Eq.(13) is reduced to

∫ dz′ ( \begin{matrix}[ξ_{∥} − \frac{1}{2m}∂_{z}^{2}]δ(z − z′) & Δ(p, z, z′) \\ Δ^{†}(p, z, z′) & −[ξ_{∥} − \frac{1}{2m}∂_{z}^{2}]δ(z − z′) \end{matrix} ) Ψ(z′) = EΨ(z),  (14)

where ξ_{∥} = p^{2}/2m − μ,

Ψ(z) = Ψ(p, z) = ∫ dxe^{−ip·x}Ψ(r),  (15)
Δ(p, z, z′) = ∫ d(x − x′)e^{−ip·(x−x′)}Δ(r, r′)  (16)

and

Δ(p, z, z′) = ∑_{p_{z}>0} Δ(p_{F}^{+}, z)e^{ip_{z}(z−z′)} + ∑_{p_{z}<0} Δ(p_{F}^{−}, z)e^{ip_{z}(z−z′)}  (17)

with Δ(p_{F}^{±}, z) the position dependent pair function at the Fermi momentum p_{F}^{±}. In Eq.(17), we have defined two Fermi momenta p_{F}^{+} and p_{F}^{−} associated with the parallel momentum p, as depicted in Fig.2.

0  θ p_{F}_{+} p_{F}_{−} p  p_{Fz} −p_{Fz}  z (x−y)

Figure 2: Fermi momenta p_{F}^{α} (α = ±) associated with the momentum p parallel to the interface and to the walls.

2.1.2 The Andreev equation

We first derive the Andreev equation from the Bogoliubov–de Gennes equation Eq.(14) using the quasi–classical (WKBJ) approximation[5]. Also, we discuss the boundary condition for the Andreev amplitude, following Zaǐtsev[23], Shelankov[22] and Ashida et al[20].

When the characteristic size of the system is much longer than the Fermi wave length 1/p_{F}, the quasi–classical approximation introduced by Andreev [5] is useful. Following Andreev, we write the Nambu amplitude as

Ψ(z) = Φ_{+}(z)e^{ip_{F z}z} + Φ_{−}(z)e^{−ip_{F z}z},  (18)

where p_{F z} = p_{F} cos θ is the z component of the Fermi momentum p_{F}^{+}. The slowly varying amplitude Φ_{α}(z) obeys the Andreev equation[5]

( \begin{matrix}−αiv_{F z}∂_{z} & Δ(p_{F}^{α} , z) \\ Δ^{†}(p_{F}^{α} , z) & αiv_{F z}∂_{z} \end{matrix} ) Φ_{αl} = E_{l}Φ_{αl},  (19)

where v_{F z} = v_{F} cos θ = p_{F} cos θ/m. The slowly varying amplitudes are expected to form a complete set near the given Fermi momentum in a sense that

∑_{l} Φ_{αl}(z)Φ_{αl}^{†}(z′) = δ̃(z − z′) ≡ ∫_{−p_{F z}}^{p_{F z}} dp_{z}e^{ip_{z}(z−z′)}  (20)

Since the Andreev equation is a first order differential equation, its formal solution can be written as

Φ_{α}(z) = ρ_{3}U_{α}(z, z′, E)ρ_{3}Φ_{α}(z′),  (21)

where U_{α}(z, z′, E) is a spatial evolution operator which obeys

iv_{F z}∂_{z}U_{α}(z, z′, E) = −α(E − Δ̂_{α}(z))ρ_{3}U_{α}(z, z′, E),  (22)
Δ̂_{α}(z) = ( \begin{matrix}0 & Δ(p_{F}^{α} , z) \\ Δ^{†}(p_{F}^{α} , z) & 0 \end{matrix} )  (23)

and

U_{α}(z, z) = 1.  (24)

In the above, ρ_{1}, ρ_{3} are the Pauli matrices in the particle-hole space. The properties of the evolution operator is discussed in detail in the appendix.

L  C  R  Z 0  L −L_{L}  L + L_{R} Δ(Z)

Figure 3: The finite triple layer system.

To discuss the semi–infinite triple layer system, we start from the finite width triple layer system as is depicted in Fig.3. Let us discuss the boundary condition at the L − C, C − R interfaces and the layer ends at z = −L_{L}, z = L + L_{R}. We use superscripts or subscripts L, C and R to denote the quantity in the L, C and R layer, respectively. Since the parallel momentum p is conserved,

p_{F}^{L} sin θ_{L} = p_{F}^{C} sin θ_{C} = p_{F}^{R} sin θ_{R}.  (25)

The boundary conditions for Φ_{α}’s at the interfaces[20][22][23] can be obtained by the following arguments.

For the L–C interface boundary located at z = 0, the proximity interface is assumed to be confined within a narrow range −δ < z < δ, where δ is of order of 1/p_{F}. It follows that the Nambu amplitude Ψ should be the asymptotic solution of the interface problem. Since the interface process is a high-energy ( ∼ E_{F} ) and short range ( ∼ 1/p_{F} ) process, it is governed by the rapidly varying part e^{ip_{F} z} of the Nambu amplitudes and is consequently common to both the superfluid and the normal phases. The interface process is, therefore, characterized by the two independent asymptotic solutions Ψ_{1} and Ψ_{2} for the electrons at the Fermi level :

Ψ_{1} = \{ \begin{matrix}e^{ip_{F z}^{L}z} + r_{0}e^{−ip_{F z}^{L}z} & for z < −δ \\ d_{0}e^{ip_{F z}^{C}z} & for z > δ, \end{matrix}  (26)
Ψ_{2} = \{ \begin{matrix}d̃_{0}e^{−ip_{F z}^{L}z} & for z < −δ \\ e^{−ip_{F z}^{C}z} + r̃_{0}e^{ip_{F z}^{C}z} & for z > δ, \end{matrix}  (27)

where r_{0} and d_{0} are the reflection and the transmission amplitude at the L–C interface, respectively, and d̃_{0} = d_{0}v_{F z}^{C}/v_{F z}^{L}, r̃_{0} = −r_{0}^{∗}d_{0}/d_{0}^{∗} which are obtained by use of the constantness of the Wronskian. Noting that the Nambu amplitudes Ψ^{L} and Ψ^{C} are given by appropriate linear combinations of Ψ_{1} and Ψ_{2},

Ψ = AΨ_{1} + BΨ_{2}

we obtain the boundary conditions. For the C–R interface boundary located at z = L, the boundary condition can be obtained in the same way. Thus, the boundary conditions for the Andreev amplitudes at z = 0 and z = L are written as

Φ_{+}^{L}(0) = \frac{1}{d_{0}} Φ_{+}^{C}(0) + \frac{r_{0}^{∗}}{d_{0}^{∗}} Φ_{−}^{C}(0),  (28)
Φ_{−}^{L}(0) = \frac{r_{0}}{d_{0}} Φ_{+}^{C}(0) + \frac{1}{d_{0}^{∗}} Φ_{−}^{C}(0),  (29)
Φ_{+}^{C}(L) = \frac{1}{d_{L}} Φ_{+}^{R}(L) + \frac{r_{L}^{∗}}{d_{L}^{∗}} Φ_{−}^{R}(L),  (30)

Φ_{−}^{C}(L) = \frac{r_{L}}{d_{L}} Φ_{+}^{R}(L) + \frac{1}{d_{L}^{∗}} Φ_{−}^{R}(L),  (31)

where r and d are the reflection and the transmission amplitude for the electronic state at the Fermi level and they can be expressed in terms of the reflection coefficients R_{0}, R_{L} of the interface as

r_{0} = \sqrt{R_{0}} e^{iθ_{r_{0}}},  (32)
d_{0} = \sqrt{\frac{v_{F z}^{L}}{v_{F z}^{C}} (1 − R_{0})} e^{iθ_{d_{0}}},  (33)
r_{L} = \sqrt{R_{L}} e^{iθ_{r_{L}}+2ip_{F z}^{C}L},  (34)
d_{L} = \sqrt{\frac{v_{F z}^{C}}{v_{F z}^{R}} (1 − R_{L})} e^{iθ_{d_{L}}+i(p_{F z}^{C}−p_{F z}^{R})L},  (35)

where θ_{r}, θ_{d} are the phase of r, d, respectively. As we shall eventually show, the final results depend only on the reflection coefficients R_{0} = |r_{0}|^{2}, R_{L} = |r_{L}|^{2} and not on the phases θ_{r} and θ_{d}. The boundary conditions at the ends of the layers are written as

Φ_{+}^{L}(−L_{L})e^{−ip_{F z}^{N}L_{L}} + Φ_{−}^{L}(−L_{L})e^{ip_{F z}^{N}L_{L}+iη_{L}} = 0,  (36)
Φ_{+}^{R}(L + L_{R})e^{ip_{F z}^{S}(L+L_{R})+iη_{R}} + Φ_{−}^{R}(L + L_{R})e^{−ip_{F z}^{R}(L+L_{R})} = 0,  (37)

where η_{L}, η_{R} are possible phase shift at the boundaries. When exp iη = 1, one has a fixed end condition and when exp iη = −1 one has a free end condition. As we shall see, the physical quantities does not depend on these phases. One can study an inhomogeneous superconducting system by use of the Bogoliubov de Gennes equation Eq.(14) or by use of the Andreev equation Eq.(19) derived here. In this thesis, to discuss the proximity system, we have mainly treat it by use of the quasi–classical Green’s function discussed in later sections.

2.2 The quasi–classical Green’s function for a finite triple layer system

In this section we closely follow the prescription proposed by Ashida et al.[20].

We consider a geometry as shown in Fig. 3 . The widths of the layers L_{L}, L and L_{R} are assumed to be still much longer than the Fermi wavelength 1/p_{F}, but not necessarily longer than the coherence length ξ. The interface and boundaries are assumed to have translational symmetry in the x, y plane. The momentum component parallel to be the boundary is, therefore, a conserved quantity. From now on, we use superscripts or subscripts L, C and R to denote the quantity in the L, C and R layer, respectively.

The Green’s function method[3][4] is a useful method to study the superconducting proximity contact system. The Gor’kov Green’s function is written in terms of the Nambu amplitude Ψ(r) as

G(r, r′, ε) = ∑_{l} \frac{Ψ_{l}(r)Ψ_{l}^{†}(r′)}{ε − E_{l}}.  (38)

Since the momentum component p parallel to the interface is a conserved quantity in the present system, we have only to consider the spatial variation in z–direction. Hence

G(r, r′, ε) = \frac{1}{(2π)^{2}} ∫ ∫ d^{2}p G(p, z, z′, ε)e^{ip·(x−x′)}.  (39)

Then, the Green’s function, which depends on the z direction, is written as

G(z, z′, ε) = ∑_{l} \frac{Ψ_{l}(z)Ψ_{l}^{†}(z′)}{ε − E_{l}},  (40)

which satisfies the Gor’kov equation

∫ dz′′ [ε − ( \begin{matrix}[ξ_{∥} − \frac{1}{2m}∂_{z}^{2}]δ(z − z′′) & Δ(p, z, z′′) \\ Δ^{†}(p, z, z′′) & −[ξ_{∥} − \frac{1}{2m}∂_{z}^{2}]δ(z − z′′) \end{matrix} )] G(z′′, z′, ε) = δ(z − z′). (41)

Substituting Eq.(18) into Eq.(40), one finds that Gor’kov function is decomposed as

G(p, z, z′, ε) = G_{++}(z, z′, ε) e^{ip_{F z}(z−z′)}

+ G_{+−}(z, z′, ε) e^{ip_{F z}(z+z′)}
+ G_{−+}(z, z′, ε) e^{−ip_{F z}(z+z′)}
+ G_{−−}(z, z′, ε) e^{−ip_{F z}(z−z′)},  (42)

where

G_{αβ}(z, z′, ε) = ∑_{l} \frac{Φ_{αl}(z)Φ_{βl}^{†}(z′)}{ε − E_{l}},  (43)

which satisfies the equation

[ε − ( \begin{matrix}−αiv_{F z}∂_{z} & Δ(p_{F}^{α} , z) \\ Δ^{†}(p_{F}^{α} , z) & αiv_{F z}∂_{z} \end{matrix} )] G_{αβ}(z, z′, ε) = δ_{αβ}δ̃(z − z′).  (44)

It is useful for later use to introduce the notion of ”directional space” which is a two dimensional space spanned by α = ± . From Eq.(44), one finds that the diagonal ( in the directional space ) elements G_{++} and G_{−−} have a jump at z = z′, i.e.,

G_{αα}(z + 0, z) − G_{αα}(z − 0, z) = −i\frac{α}{v_{F z}} ρ_{3}  (45)

but off–diagonal elements have no jump. Here ρ_{3} is a Pauli matrix in particle–hole space. Noting the jump given by Eq.(45), we define the quasi–classical Green’s function ĝ_{αβ} by

ĝ_{αβ}(z) ± i(γ_{3})_{αβ} = −2πv_{F z}ρ_{3}G_{αβ}(z ± 0, z),  (46)

where γ_{3} is a Pauli matrix in the directional space. Then, the position diagonal element of the Gor’kov Green’s function is written in terms of the quasi–classical Green’s function as

G(p, z, z, ε) = −\frac{1}{2v_{F z}} ρ_{3} [ ĝ_{++}(z)
+ ĝ_{+−}(z) e^{2ip_{F z}z}
+ ĝ_{−+}(z) e^{−2ip_{F z}z}
+ ĝ_{−−}(z) ]. (47)

The quasi–classical Green’s function ĝ_{+−} and ĝ_{−+} do not contribute to observed quantities because they are accompanied by a rapidly changing phase factor e^{±2ip_{F z}z} that vanishes under the p summation. In the conventional quasi-classical Green’s function theory,[25] therefore, only ĝ_{++} and ĝ_{−−} are considered. To solve the boundary problem, however, it is more convenient to treat all the components of ĝ_{αβ}(z) on equal footing, because the boundary problem then can be solved within linear algebra.

The generalized set of the quasi-classical Green’s functions obeys the Eilenberger equation[18],[19]

iv_{F z}∂_{z}ĝ_{αβ} = −α(ε − Δ̂_{α}(z))ρ_{3}ĝ_{αβ} + ĝ_{αβ}β(ε − Δ̂_{β}(z))ρ_{3}.  (48)

Since the Eilenberger equation is a first order differential equation, a formal solution can be expressed in terms of the evolution operator given by Eq.(23) in the previous section as follows :

ĝ_{αβ}(z) = U_{α}(z, z′)ĝ_{αβ}(z′)U_{β}(z′, z).  (49)

In treating the Green’s function, we use the evolution operator with complex ε.

The boundary conditions for the quasi-classical Green’s function at the interfaces have been given by Zaǐtsev[23], Millis et al.[24], Shelankov[22] and Ashida et al.[20]. Using the boundary condition Eq.(31) for the amplitude of the Andreev equation and the definition of the quasi–classical Green’s function, one can obtain the boundary condition of ĝ at the interface boundaries as follows :

ĝ_{αβ}^{L}(0) = ∑_{γδ} M̂_{0}_{αγ} ĝ_{γδ}^{C}(0) M̂_{0}_{δβ}^{†},
ĝ_{αβ}^{C}(0) = ∑_{γδ} M̂_{L}_{αγ} ĝ_{γδ}^{R}(0) M̂_{L}_{δβ}^{†},  (50)

where M̂_{0} and M̂_{L} are matrices

M̂_{0} = \sqrt{\frac{v_{F z}^{L}}{v_{F z}^{C}}} ( \begin{matrix}\frac{1}{d_{0}} & \frac{r_{0}^{∗}}{d_{0}^{∗}} \\ \frac{r_{0}}{d_{0}} & \frac{1}{d_{0}^{∗}} \end{matrix} ) ,  (51)
M̂_{L} = \sqrt{\frac{v_{F z}^{C}}{v_{F z}^{R}}} ( \begin{matrix}\frac{1}{d_{L}} & \frac{r_{L}^{∗}}{d_{L}^{∗}} \\ \frac{r_{L}}{d_{L}} & \frac{1}{d_{L}^{∗}} \end{matrix} ) .  (52)

In the same way, using Eq.(37), the boundary conditions at the ends of the layers are

ĝ_{++}^{L}(−L_{L}) ± i = ĝ_{−−}^{L}(−L_{L}) ± i = −e^{∓iη̂_{L}}ĝ_{±∓}^{L}(−L_{L}),  (53)
ĝ_{++}^{R}(L + L_{R}) ± i = ĝ_{−−}^{R}(L + L_{R}) ± i = −e^{∓iη̂_{R}}ĝ_{±∓}^{R}(L + L_{R}),  (54)

where η̂_{L} = 2p_{F z}^{L}L_{L} + η_{L}, η̂_{R} = 2p_{F z}^{R}(L + L_{R}) + η_{R} and η_{L}, η_{R} are possible phase shift at the boundaries.

Using the boundary conditions Eqs.(50)–(54) and the spatial evolution operator of the ĝ given by Eq.(49), one can obtain the quasi–classical Green’s function at arbitrary position. The explicit expression of the quasi–classical Green’s function can be written as follows :

ĝ_{αα}^{L}(z) = U_{α}^{L}(z, −L_{L}) h^{L} U_{α}^{L}(−L_{L}, z),  (55)
ĝ_{αα}^{R}(z) = U_{α}^{R}(z, L + L_{R}) h^{R} U_{α}^{R}(L + L_{R}, z),  (56)
ĝ_{αα}^{C}(z) = U_{α}^{C}(z, L) h_{α}^{C} U_{α}^{C}(L, z),  (57)

where

h^{L} = (−i)\frac{e^{iϕ_{L}} + A_{L}}{e^{iϕ_{L}} − A_{L}} ,  (58)
A_{L} = (U_{+}^{L})^{−1} \frac{A_{C} − \sqrt{R_{0}}e^{iϕ_{C}}}{e^{iϕ_{C}} − \sqrt{R_{0}}A_{C}} U_{−}^{L},  (59)
A_{C} = U_{+}^{C} \frac{A_{R} − \sqrt{R_{L}}e^{iϕ_{R}}}{e^{iϕ_{R}} − \sqrt{R_{L}}A_{R}} (U_{−}^{C})^{−1},  (60)
A_{R} = U_{+}^{R}(U_{−}^{R})^{−1}  (61)

and

h^{R} = (−i)\frac{e^{iϕ_{R}} + B_{R}}{e^{iϕ_{R}} − B_{R}} ,  (62)
B_{R} = (U_{−}^{R})^{−1} \frac{B_{C} + \sqrt{R_{L}}e^{iϕ_{C}}}{e^{iϕ_{C}} + \sqrt{R_{L}}B_{C}} U_{+}^{R},  (63)
B_{C} = (U_{−}^{C})^{−1} \frac{B_{L} + \sqrt{R_{0}}e^{iϕ_{L}}}{e^{iϕ_{L}} + \sqrt{R_{0}}B_{L}} U_{+}^{C},  (64)
B_{L} = U_{−}^{L}(U_{+}^{L})^{−1}  (65)

and

h_{±}^{C} = (1 − R_{L})^{−1}[U_{±}^{R}h^{R}(U_{±}^{R})^{−1}R_{L}U_{∓}^{R}h^{R}(U_{−}^{R})^{−1}  (66)
−\sqrt{R_{L}}(e^{iϕ_{R}}U_{−}^{R}(h^{R} + i)(U_{+}^{R})^{−1} + e^{−iϕ_{R}}U_{+}^{R}(h^{R} − i)(U_{−}^{R})^{−1})]  (67)

where

U_{±}^{L} = U_{±}^{L}(0, −L_{L}),  (68)
U_{±}^{C} = U_{±}^{C}(0, L),  (69)
U_{±}^{R} = U_{±}^{R}(L, L + L_{R}),  (70)
ϕ_{L} = 2p_{F z}^{L}L_{L} + η_{L} + θ_{r_{0}},  (71)
ϕ_{C} = 2p_{F z}^{C}L + θ_{r_{L}} − θ_{r_{0}} + 2θ_{d_{0}},  (72)
ϕ_{R} = 2p_{F z}^{R}(L + L_{L}) + η_{R} − θ_{r_{L}} + 2θ_{d_{L}}.  (73)

The phase factors e^{iϕ_{L}}, e^{iϕ_{C}} and e^{iϕ_{R}} are rapidly varying functions of layer sizes L_{L}, L and L_{R} and also of the polar angles of the Fermi momentum, because p_{F} L ≫ 1. We are not interested, however, in the size accuracy of order 1/p_{F} nor in the accuracy of polar angle of order 1/p_{F} L. Apart from the correction of order 1/p_{F} L, therefore, the physical quantities of interest are obtained from the quasi–classical Green’s function averaged over the phases ϕ_{L}, ϕ_{C} and ϕ_{R}. Moreover, in actual systems the phases θ_{r_{0}(r_{L})}, θ_{d_{0}(d_{L})} and η_{L(R)} will be random variables reflecting the microscopic irregularities at the interfaces and at the walls. It follows that the averages over the phases ϕ_{L}, ϕ_{C} and ϕ_{R} can be performed independently in spite of the fact that the polar angles of the Fermi momenta in the L, C and R layers are connected by Eq.(25). As a result, the averaged Green’s function do not depend on the phase of the interface reflection amplitudes but are determined only by the reflection coefficients R_{0} and R_{L}.

Then, to calculate the physical quantities from the quasi–classical Green’s function the finite triple layers system, one has to average it over the phases. AAHN analytically averaged the Green’s function of the finite double layers over the two phases. In this thesis, however, we do not consider the finite triple layers system, but are interested in semi– infinite triple layers systems, such as one point contact system and so on. In the next section, to treat the semi–infinite triple systems, therefore, we take the limit of L_{L} and L_{R} to infinity in the Green’s function obtained in this section. As a result, the Green’s functions taken the limit do not depend on the phases ϕ_{L} and ϕ_{R}. We have only to average those over the phase ϕ_{C}.

From now on, we confine ourselves to the singlet superconducting proximity system, although the quasi–classical Green’s function obtained in this section can also be used to investigate the triplet superconducting system. For a singlet superconductor, the pair potential is written in the spin space as

Δ(p_{F}^{α} , z) = ( \begin{matrix}0 & Δ_{↑↓} \\ Δ_{↓↑} & 0 \end{matrix} ) ,  (74)

therefore, we can decouple the spin space. For a singlet superconductor, we have only to treat a 2 × 2 matrix in the particle–hole space. Properties of the evolution operator in the 2 × 2 matrix space are discussed in the appendix.

In equilibrium states without supercurrent, we can take the pair potential real. Then, it can be shown that the evolution operator satisfies

U_{α}(z, z′) = ρ_{1}U_{−α}(z, z′)ρ_{1}.  (75)

( See appendix. ) Using this relation Eq.(75), we find

ĝ_{++} = ρ_{3}^{T}ĝ_{−−}ρ_{3}.  (76)

Therefore one has only to treat ĝ_{++}.

2.3 The quasi–classical Green’s function of a semi–infinite triple layer system

In this section we derive the quasi–classical Green’s function for the semi–infinite geometry as shown in Fig. 1 [15]. As we have noted at the end of the previous section, we consider a singlet superconducting proximity contact system in this thesis. We assume that the pair potential is real. The case when the pair potential is not real will be discussed in Sec.5. To obtain the quasi–classical Green’s function of the semi–infinite triple layer system, we take the limit of L_{L}, L_{R} to infinity in the Green’s function of the finite triple layer system obtained in the previous section[15]. In this case, the pair potential Δ(z) will tend to the bulk value Δ_{bulk} at sufficiently large |z|. Hence, the evolution operator U can be divided into a growing part and a damping part :

U_{+}(z, z′, ε) = Λ_{+}(z, z′)e^{−iκ(z−z′)} + Λ_{−}(z, z′)e^{iκ(z−z′)},  (77)
U_{−}(z, z′, ε) = ρ_{1}Λ_{+}(z, z′)ρ_{1}e^{−iκ(z−z′)} + ρ_{1}Λ_{−}(z, z′)ρ_{1}e^{iκ(z−z′)},  (78)

where κ = Ω/v_{F z} ≡ \sqrt{ε^{2} − Δ_{bulk}^{2}}/v_{F z} and the square root is defined to have a positive imaginary part in the complex ε plane. The operator Λ_{±} has a projection operator like properties as is shown in detail in the appendix. When L_{R} → ∞, since Im κ > 0, the evolution operator U_{+}(L, L + L_{R}) becomes

U_{+}(L, L + L_{R}) = Λ_{+}(L, L + L_{R})e^{iκL_{R}} + Λ_{−}(L, L + L_{R})e^{−iκL_{R}},

⟶ Λ_{−}(L, L + L_{R})e^{−iκL_{R}}|_{L_{R}→∞} a divergent term.

Retaining the most divergent terms when L_{L} → ∞ and L_{R} → ∞, we can find the quasi–classical Green’s function of the semi–infinite triple layer system.

L side ( z < 0 )

ĝ_{++}^{L}(z) = i\frac{2A_{L}^{∞}(z) − trA_{L}^{∞}(z)}{trA_{L}^{∞}(z)},  (79)
A_{L}^{∞}(z) = Ǔ_{+}^{L}(z, 0)(A_{C}^{∞} + R_{0}ρ_{2}^{T}A_{C}^{∞}ρ_{2} + 2\sqrt{R_{0}R_{L}}trA_{R}^{∞} cos ϕ_{C})
×ρ_{1}Λ_{+}^{L}(0, −∞)ρ_{3}^{T}Λ_{+}^{L}(z, −∞)ρ_{3}ρ_{1},  (80)
A_{C}^{∞} = Ǔ_{+}^{C}(0, L)(A_{R}^{∞} + R_{L}ρ_{2}^{T}A_{R}^{∞}ρ_{2})ρ_{3}^{T} Ǔ_{+}^{C}(0, L)ρ_{3},  (81)
A_{R}^{∞} = ϕ̌_{−}^{R}(L)^{T}ϕ̌_{−}^{R}(L)ρ_{3},  (82)
R side ( z > L )
ĝ_{++}^{R}(z) = i\frac{2B_{R}^{∞}(z) − trB_{R}^{∞}(z)}{trB_{R}^{∞}(z)},  (83)
B_{R}^{∞}(z) = Λ_{−}^{R}(z, ∞)ρ_{3}^{T}Λ_{−}^{R}(L, ∞)ρ_{3}
×(B_{C}^{∞} + R_{L}ρ_{2}^{T}B_{C}^{∞}ρ_{2} + 2\sqrt{R_{0}R_{L}}trB_{L}^{∞} cos ϕ_{C})ρ_{2}^{T} Ǔ_{+}^{R}(z, L)ρ_{2},  (84)
B_{C}^{∞} = ρ_{3}^{T} Ǔ_{+}^{C}(0, L)ρ_{3}(B_{L}^{∞} + R_{0}ρ_{2}^{T}B_{L}^{∞}ρ_{2})Ǔ_{+}^{C}(0, L),  (85)
B_{L}^{∞} = ρ_{1}ϕ̌_{+}^{L}(0)^{T}ϕ̌_{+}^{L}(0)ρ_{3}ρ_{1}  (86)

and C side ( 0 < z < L )

ĝ_{++}^{C}(z) = i\frac{2C_{C}^{∞}(z) − trC_{C}^{∞}(z) + 2i\sqrt{R_{0}R_{L}}trA_{R}^{∞}trB_{L}^{∞} sin ϕ_{C}}{trC_{C}^{∞}(z) + 2\sqrt{R_{0}R_{L}}trA_{R}^{∞}trB_{L}^{∞} cos ϕ_{C}} ,  (87)

C_{C}^{∞}(z) = Ǔ_{+}^{C}(z, L)(A_{R}^{∞} + R_{L}ρ_{2}^{T}A_{R}^{∞}ρ_{2})
×ρ_{3}^{T} Ǔ_{+}^{C}(0, L)ρ_{3}(B_{L}^{∞} + R_{0}ρ_{2}^{T}B_{L}^{∞}ρ_{2})
×Ǔ_{+}^{C}(0, L)ρ_{2}^{T} Ǔ_{+}^{C}(z, L)ρ_{2}.  (88)

where

Ǔ_{+}^{C}(z, L) = U_{+}^{C}(z, L)e^{−iκ^{C}(z−L)},  (89)
Ǔ_{+}^{R}(z, L) = U_{+}^{R}(z, L)e^{iκ^{R}(z−L)},  (90)
Ǔ_{+}^{L}(z, 0) = U_{+}^{L}(z, 0)e^{−iκ^{L}z}  (91)

and ϕ̌_{+} and ϕ̌_{−} are decomposed elements of the evolution operator defined in the appendix.

In the course of taking the limits, the phase factors, e^{iϕ_{L}} and e^{iϕ_{R}}, disappear. Since the C layer has finite width, the phase factor e^{iϕ_{C}} remains in the Green’s function.

e^{iϕ_{C}} = e^{i(2p_{F z}^{C}L+θ_{r_{L}}−θ_{r_{0}}+2θ_{d_{0}})}.  (92)

As we have noted in the previous section, this phase factor is a rapidly varying function of layer size L and also of the polar angles of the Fermi momentum. We are not interested, however, in the size accuracy of order 1/p_{F} and in the accuracy of polar angle of order 1/p_{F} L. Apart from the correction of order 1/p_{F} L, therefore, the physical quantities of interest are obtained from the quasi–classical Green’s function averaged over the phase ϕ_{C}. As a result, the averaged Green’s function do not depend on the phase of the interface reflection amplitudes but are determined only by the reflection coefficients R_{0} and R_{L}. The averaged Green’s function is defined by

⟨ĝ_{++}(z)⟩ ≡ \frac{1}{2π} ∫_{0}^{2π} dϕ_{C} ĝ_{++}(z)  (93)

For example, the averaged Green’s function of the C layer is written as

⟨ĝ_{++}(z)⟩ = i\frac{2C_{C}^{∞}(z) − trC_{C}^{∞}(z)}{\sqrt{(trC_{C}^{∞}(z))^{2} − 4R_{0}R_{L}(trA_{R}^{∞}trB_{L}^{∞})^{2}}}  (94)

We study some proximity contact systems by use of the averaged Green’s function.

2.4 The Green’s function and physical quantities

We discuss how the physical quantities can be calculated from the Green’s function. We start from the Bogoliubov–de Gennes equation Eq.(13) ,

∫ dr′E(r, r′)Ψ_{l}(r′) = E_{l}Ψ_{l}(r).  (95)

Using positive energy solution of this equation,

Ψ_{l}(r) = ( \begin{matrix}u_{l↑}(r) \\ u_{l↓}(r) \\ v_{l↑}(r) \\ v_{l↓}(r) \end{matrix} ) ,  (96)

we define a matrix U_{l}, which diagonalize the Hamiltonian Eq.(6),

U_{l}(r) = ( \begin{matrix}u_{l↑}(r) & v_{l↑}(r)^{∗} \\ u_{l↓}(r) & v_{l↓}(r)^{∗} \\ v_{l↑}(r) & u_{l↑}(r)^{∗} \\ v_{l↓}(r) & u_{l↓}(r)^{∗} \end{matrix} ) .  (97)

This matrix satisfies

∑_{l}′ U_{l}(r)U_{l}^{†}(r′) = δ(r − r′),

where the symbol _{∑l}′ indicate the sum over positive energy states. Also, matrix E(r, r′) can be reduced to

E(r, r′) = ∑_{l}′ U_{l}(r) ( \begin{matrix}E_{l} & 0 \\ 0 & −E_{l} \end{matrix} ) U_{l}^{†}(r′).  (98)

Using Eq.(10), we reduce Eq.(6) to

H^{Gor′kov} = ∫ dr ∫ dr′Ψ̂^{†}(r)E(r, r′)Ψ̂(r′) + const. .  (99)

The Hamiltonian can be rewritten by Eq.(98), i.e.,

H^{Gor′kov} = ∑_{l}′ ∫ dr(U_{l}^{†}(r)Ψ̂(r))^{†} ( \begin{matrix}E_{l} & 0 \\ 0 & −E_{l} \end{matrix} ) ∫ dr′U_{l}^{†}(r′)Ψ̂(r′) + const.
= ∑_{l}′ γ̂_{l}^{†} ( \begin{matrix}E_{l} & 0 \\ 0 & −E_{l} \end{matrix} ) γ̂_{l} ,  (100)
γ̂_{l} = ( \begin{matrix}γ_{l} \\ γ_{l}^{†} \end{matrix} ) ,  (101)

where γ_{l} and γ_{l}^{†} are an annihilation and a creation operators of the Bogolon ( a quasi– particle ), respectively. In terms of the operator γ̂_{l}, the electron operator Ψ̂ is written as

Ψ̂(r) = ∑_{l}′ U_{l}(r)γ̂_{l},  (102)

more explicitly

ψ_{α}(r) = ∑_{l}′(u_{l,α}(r)γ_{l} + v_{l,α}(r)^{∗}γ_{l}^{†}),  (103)
ψ_{α}^{†}(r) = ∑_{l}′(v_{l,α}(r)γ_{l} + u_{l,α}(r)^{∗}γ_{l}^{†}).  (104)

One calls these equations the Bogoliubov transformation.

(i) The gap equation The pair potential Δ_{αβ}(r, r′) was defined by

Δ_{αβ}(r, r′) = v(|r − r′|)⟨ψ_{β}(r′)ψ_{α}(r)⟩.  (105)

The gap equation can be rewritten in terms of the new operators γ_{l} and γ_{l}^{†},

Δ_{αβ}(r, r′) = v(|r − r′|)⟨∑_{l,l′}′(u_{l,β}(r′)γ_{l} + v_{l,β}(r′)^{∗}γ_{l}^{†})
×(u_{l′,α}(r)γ_{l′} + v_{l′α}(r)^{∗}γ_{l′}^{†})⟩,

= v(|r − r′|) ∑_{l}′ {u_{l,β}(r′)v_{l,α}(r)^{∗}⟨γ_{l}γ_{l}^{†}⟩
+v_{l,α}(r′)^{∗}u_{lβ}(r)⟨γ_{l}^{†}γ_{l}⟩}
= v(|r − r′|) ∑_{l}′ {u_{l,β}(r′)v_{l,α}(r)^{∗}(1 − f(E_{l})) + v_{l,α}(r′)^{∗}u_{l,β}(r)f(E_{l})}
= T v(|r − r′|) ∑_{ω_{n}} ∑_{l}′ \{\frac{u_{l,α}(r)v_{l,β}(r′)^{∗}}{iω_{n} − E_{l}} + \frac{v_{l,α}(r)^{∗}u_{l,β}(r′)}{iω_{n} + E_{l}}\} ,  (106)

where T is the temperature, and we have used

⟨γ_{l}^{†}γ_{l}⟩ = f(E_{l}) = \frac{1}{1 + e^{βE_{l}}}
= T ∑_{ω_{n}} \frac{e^{iω_{n}δ}}{iω_{n} − E_{l}} ,  (107)

ω_{n} = πT (2n + 1) is the Matsubara frequency ( n = 0, ±1, ±2, ... ), and e^{iω_{n}δ} is a convergent factor omitted above. Rewriting the sum of the second term of Eq.(106) to the sum over the negative energy state, one can obtain

Δ_{αβ}(r, r′) = T v(|r − r′|) ∑_{l} ∑_{ω_{n}} \frac{u_{l,α}(r)v_{l,β}(r′)^{∗}}{iω_{n} − E_{l}} ,  (108)

where ∑_{l} means the sum over the all states.

Comparing Eq.(108) with the definition of the Gor’kov Green’s function,

G(r, r′, ε) = ∑_{l} \frac{Ψ_{l}(r)Ψ_{l}^{†}(r′)}{ε − E_{l}},  (109)

one can obtain the gap equation in terms of the Gor’kov Green’s function G, i.e.,

Δ(r, r′) = T v(|r − r′|) ∑_{ω_{n}} G(r, r′, iω_{n})|_{1,2},  (110)

where G|_{1,2} is an off–diagonal (1, 2) element of G in the particle–hole space.

When the system has the translational symmetry in the x–y plane, using Eqs.(16)(17), Eq.(39) and Eq.(47), one can obtain the gap equation in terms of the quasi–classical Green’s function

Δ(p_{F} , z) = −\frac{T N(0)}{4} ∑_{ω_{n}} ∑_{α=±} ∫_{0}^{2π} dϕ^{α} ∫0^{\frac{π}{2}} sin θ^{α}dθ^{α}v_{p_{F} ·p_{F}^{α}} ĝ_{αα}(iω_{n}, z)|_{1,2},  (111)

where θ^{α} and ϕ^{α} are the polar angles of the Fermi momentum p_{F}^{α}, v_{p_{F} ·p_{F}^{α}} is the pairing interaction, and N(0) is the density of states at the Fermi surface.

(ii) The electric current We derive an expression of the electric current in terms of the Green’s function. As is well known, the electric current J(r) is defined as

J(r) = ∑_{α}⟨\{\frac{ie}{2m} [(∇_{r}ψ_{α}^{†}(r))ψ_{α}(r) − ψ_{α}^{†}(r)∇_{r}ψ_{α}(r)] − \frac{e^{2}}{mc}Aψ_{α}^{†}(r)ψ_{α}(r)\}⟩,  (112)

where A is a vector potential, ψ_{α} and ψ_{α}^{†} are operators of an electron with spin α and the symbol ⟨ ⟩ denotes the thermal average. From now on, we ignore the last term, which is proportional to a vector potential A, and discuss on

J(r) = \frac{ie}{2m} ∑_{α}(∇_{r′} − ∇_{r})⟨ψ_{α}^{†}(r′)ψ_{α}(r)⟩|_{r′→r}.  (113)

The term ⟨ψ^{†}(r′)ψ(r)⟩ is rewritten using the quasi–particle operators as follows :

∑_{α}⟨ψ_{α}^{†}(r′)ψ_{α}(r)⟩ = ∑_{α} ∑_{l,l′}′⟨(v_{l,α}(r′)γ_{l} + u_{l,α}(r′)^{∗}γ_{l}^{†})(u_{l′,α}(r)γ_{l′} + v_{l′,α}(r)^{∗}γ_{l′}^{†})⟩
= T ∑_{ω_{n}} ∑_{α} ∑_{l}′ [\frac{u_{l,α}(r′)^{∗}u_{l,α}(r)}{iω_{n} − E_{l}} + \frac{v_{l,α}(r′)v_{l,α}(r)^{∗}}{iω_{n} + E_{l}} ]
= T ∑_{ω_{n}} tr_{S} G(r, r′, iω_{n})|_{1,1},  (114)

where we have used Eq.(107), G|_{1,1} is a diagonal (1, 1) element of G in the particle–hole space, and tr_{S} G|_{1,1} denotes a trace of G|_{1,1} in the spin–space. Then, one finally obtains the electric current in terms of the Gor’kov Green’s function as

J(r) = \frac{ieT}{2m} [(∇_{r′} − ∇_{r}) ∑_{ω_{n}} tr_{S} G(r, r′, iω_{n})|_{1,1}]_{r′→r} .  (115)

In terms of the quasi–classical Green’s function one can write the z–component of the electric current as

J_{z}(z) = −ev_{F} \frac{T N(0)}{2} ∫_{0}^{2π} dϕ ∫_{0}^{π/2} dθ sin θ cos θ ∑_{ω_{n}} tr_{S} [ĝ_{++}(ω_{n}, z) − ĝ_{−−}(ω_{n}, z)]_{11} . (116)

2.5 The self–consistent pair potential

In the superconducting proximity contact system, the pair potential near the interface varies on the scale of the superconducting coherence length ξ. To study the proximity effect, one has to take into account the self–consistency for the pair potential.

To obtain the self–consistent pair potential, we have to solve the gap equation. From Eq.(111), the gap equation can be written as

Δ(p_{F} , z) = \frac{T N(0)}{4} ∑_{ω_{n}} ∑_{α=±} ∫_{0}^{2π} dϕ^{α} ∫0^{\frac{π}{2}} sin θ^{α}dθ^{α}(2l + 1)g_{l}P_{l}(p_{F} · p_{F}^{α})ĝ_{αα}(iω_{n}, z)|_{1,2}, (117)

where N(0) is the density of states at the Fermi surface, and g_{l} denotes the strength of the pairing interaction for the l–th partial wave. We have defined g_{l} such that is positive when the interaction is attractive. Here, l determines the symmetry type of the order parameter of the superconductor, i.e., l = 0 for s-wave, l = 2 for d-wave superconductivity. P_{l} are Legendre polynomials and are expanded in terms of spherical harmonic functions Y_{lm}, i.e.,

P_{l}(p_{F} · p_{F}^{α}) = \frac{4π}{2l + 1} ∑_{m=−l}^{m=l} Y_{lm}(θ, ϕ)Y_{lm}^{∗}(θ^{α}, ϕ^{α}),  (118)

where θ and ϕ are the polar angles of the Fermi momentum. We confine ourselves to treat the singlet superconductors.

For the Matsubara frequency ε = iω_{n}, if only Δ(p_{F}^{+}) = Δ(p_{F}^{−}), the evolution operator satisfies

U_{+}(ω_{n}) = ρ_{2}U_{−}(ω_{n})^{∗}ρ_{2} and (119)
U_{+}(ω_{n}) = ρ_{1}U_{+}(−ω_{n})^{∗}ρ_{1},  (120)

even when the pair potential is not real. From these relations, we can easily find

ĝ_{++}(ω_{n}) = −^{T} ĝ_{−−}(ω_{n})^{∗}  and  (121)
ĝ_{++}(ω_{n}) = −ρ_{1}ĝ_{++}(−ω_{n})^{∗}ρ_{1}.  (122)

To eliminate the coupling constant g_{l} from the gap equation in favor of the transition temperature T_{C}, we introduce T_{C} by letting Δ → 0 in the gap equation :

N(0)g_{l} = \frac{1}{2πT_{C} ∑_{ω_{n}_{c}}^{ω_{C}} \frac{1}{ω_{n}_{c}}}  ω_{n}_{c} = πT_{C}(2n + 1)
= \frac{1}{∑_{n=0}^{ω_{C}/2πT_{C}} \frac{1}{n + 1/2}}
= \frac{1}{log \frac{T}{T_{C}} + ∑_{n=0}^{ω_{C}/2πT} \frac{1}{n + 1/2}}  (123)

Then, we obtain

Δ(p_{F} , z) = \frac{\frac{T}{4} ∑_{ω_{n}}^{ω_{C}} ∑_{α=±} ∫_{0}^{2π} dϕ^{α} ∫_{0}^{π/2} sin θ^{α}dθ^{α}(2l + 1)P_{l}(p_{F} · p_{F}^{α})ĝ_{αα}(iω_{n}, z)⌋_{1,2}}{log \frac{T}{T_{C}} + ∑_{n=0}^{ω_{C}/2πT} \frac{1}{n + 1/2}} , (124)

We obtain the self–consistent pair potential by solving this equation iteratively. The right hand side of this equation depends on the cutoff energy ω_{C} but converges for large ω_{C}. So far the ω_{C} used in the conventional quasi–classical theory, was at most ω_{C} ∼ 10T_{C}. In this thesis, we chose it 30T_{C} except for section 5 in which we used 10T_{C}. In contrast to the conventional quasi–classical method in which one should solve the Eilenberger equation Eq.(48) numerically under the boundary condition and normalization condition, we have only to solve the equation for the evolution operator in the present formulation since this formulation already satisfies the boundary condition. Numerical work have been considerably reduced in the present formulation. The iteration procedure is repeated until the difference of Δ(z)/Δ_{bulk} at subsequent iterations becomes smaller than 10^{−5} at every point.

Appendix A. Properties of the evolution operator

In this appendix, we consider the properties of the evolution operator, which is a 2 × 2 matrix in the particle–hole space. Here, we assume the pair potential to be real. The case of complex pair potential will be discussed in the appendix of Sec.5.

The evolution operator obeys an equation

∂_{z}U_{α}(z, z′) = αAU_{α}(z, z′) (α = ±1)  (125)
A = \frac{i}{v_{F z}} ( \begin{matrix}ε & Δ(z) \\ −Δ(z) & −ε \end{matrix} )  (126)

and the boundary condition,

U_{α}(z, z) = 1.  (127)

From the structure of the equations, it can be proved that

det U_{α}(z, z′) = 1.  (128)

Since ρ_{1}Aρ_{1} = −A, therefore

U_{−α}(z, z′) = ρ_{1}U_{α}(z, z′)ρ_{1},  (129)

we have only to consider U_{+}(z, z′). The matrix A also has an important property

Ā ≡ ρ_{2}^{T}Aρ_{2} = −A,  (130)

where ^{T}A is the transpose of A and Ā is just the cofactor matrix of A.

Let us consider a differential equation

∂_{z} ( \begin{matrix}u \\ v \end{matrix} ) = A ( \begin{matrix}u \\ v \end{matrix} ) .  (131)

It has two independent solutions

ϕ_{+}(z) = ( \begin{matrix}u_{+}(z) \\ v_{+}(z) \end{matrix} ) and ϕ_{−}(z) = ( \begin{matrix}u_{−}(z) \\ v_{−}(z) \end{matrix} )  (132)

with the nonvanishing Wronskian

W = ^{T}ϕ_{+}ρ_{2}ϕ_{−} = −^{T}ϕ_{−}ρ_{2}ϕ_{+} = −i(u_{+}(z)v_{−}(z) − u_{−}(z)v_{+}(z)) = const. ≠ 0.  (133)

The evolution operator U_{+}(z, z′) can be written in a linear combination of the above linearly independent solutions as

U_{+}(z, z′) = Q_{+}(z, z′) + Q_{−}(z, z′),  (134)
Q_{+}(z, z′) = \frac{−1}{W}ϕ_{+}(z)^{T}ϕ_{−}(z′)ρ_{2},  (135)
Q_{−}(z, z′) = \frac{1}{W} ϕ_{−}(z)^{T}ϕ_{+}(z′)ρ_{2}.  (136)

It follows that

Q_{+}(z, z”)Q_{+}(z”, z′) = Q_{+}(z, z′),  (137)
Q_{−}(z, z”)Q_{−}(z”, z′) = Q_{−}(z, z′),  (138)
Q_{+}(z, z”)Q_{−}(z”, z′) = Q_{−}(z, z”)Q_{+}(z”, z′) = 0  (139)

and also that

det Q_{±}(z, z′) = 0.  (140)

From the above results together with Eq.(127), it can be readily shown that Q_{+}(z, z) and Q_{−}(z, z) are the projection operators that project out ϕ_{+} and ϕ_{−}, respectively.

Moreover, we can show that

Q_{±}(z′, z) = Q̄_{∓}(z, z′) ≡ ρ_{2}^{T}Q_{∓}(z, z′)ρ_{2},  (141)

from which the inverse operator of U_{+}(z, z′) is given by

U_{+}^{−1}(z, z′) = U_{+}(z′, z) = Q_{+}(z′, z) + Q_{−}(z′, z)  (142)
= Q̄_{−}(z, z′) + Q̄_{+}(z, z′)  (143)
= Ū_{+}(z, z′) .  (144)

This is a natural result when we recall the fact that det U_{+}(z, z′) = 1.

Now we consider a system in which the superconductor has semi-infinite width. The pair potential Δ(z) tends to the bulk value Δ_{bulk} at z → ∞. In that case, we can choose the independent solutions so that

ϕ_{±} = ϕ̌_{±}e^{∓iκz} ≡ ( ǔ_{±} ; v̌_{±} ) e^{∓iκz},  (145)

where κ = Ω/v_{F z} ≡ \sqrt{ε^{2} − Δ_{bulk}^{2}}/v_{F z} and when z → ∞

( ǔ_{±} ; v̌_{±} ) ⟶ const. × ( \begin{matrix}ε ∓ Ω \\ −Δ_{bulk} \end{matrix} ) .  (146)

In the above, we have defined the square root Ω = \sqrt{ε^{2} − Δ_{bulk}^{2}} such that has positive imaginary part in the complex ε-plane. Thus, ϕ_{+} and ϕ_{−} are growing and damping solutions of Eq. (131), respectively. In this case,

U_{+}(z, z′) = Λ_{+}(z, z′)e^{−iκ(z−z′)} + Λ_{−}(z, z′)e^{iκ(z−z′)},  (147)
Λ_{+}(z, z′) = \frac{−1}{W}ϕ̌_{+}(z)^{T}ϕ̌_{−}(z′)ρ_{2},  (148)
Λ_{−}(z, z′) = \frac{1}{W} ϕ̌_{−}(z)^{T}ϕ̌_{+}(z′)ρ_{2}.  (149)

All the properties of Q_{±} shown above are preserved by Λ_{±}.

From Eqs. (125) and (145), it can be shown that

∂_{z} ( ǔ_{±} ; v̌_{±} ) = \frac{i}{v_{F z}} ( \begin{matrix}ε ± Ω & Δ(z) \\ −Δ(z) & −ε ± Ω \end{matrix} ) ( ǔ_{±} ; v̌_{±} ) .  (150)

Since in the theory of the Green’s function in the semi-infinite geometry only the ratio v̌_{±}/ǔ_{±} is important, we define

D_{±}(z) = −i\frac{v̌_{±}(z)}{ǔ_{±}(z)},  (151)

which satisfies

v_{F z}∂_{z}D_{±} = −2iεD_{±} + Δ(z) (D_{±}^{2} − 1)  (152)

and at z → ∞

D_{±}(∞) = \frac{iΔ_{bulk}}{ε ∓ Ω} .  (153)

3 Andreev reflection in the normal–normal–superconducting triple layer system

The current–voltage I–V characteristics of a normal metal–superconducting contact show a non–Ohmic behavior. The non–Ohmic properties of the N–S contact are considered to be due to the Andreev reflection in which an incident electron ( hole ) from the N region is retroreflected into a hole ( electron ). ( See Fig.4. ) Most of non–Ohmic properties have

Figure 4: The Andreev reflection ( right–hand side ) and the normal reflection ( left–hand side ). The arrows give the directions of the group velocity of incoming and outgoing particle.

been analyzed using the theory by Blonder, Tinkham and Klapwijk ( BTK )[10]. We briefly review the BTK formula for a differential conductance of the normal–superconducting junction.

On the basis of the semiconductor model, the electrical current I^{nn} through the normal metal–normal metal ( N–N ) junction applied the potential difference eV between the two metal is written as

I^{nn} = AN(0)ev_{F} ∫ dE{ f(E − eV) − R_{N}^{nn}f(E − eV) − (1 − R_{N}^{nn})f(E) } (154)

in the left side,

where f(E) is the Fermi function, R_{N}^{nn} is the (normal) reflection coefficient and A is a constant to be determined by the geometry of the junction. Here, it is assumed that the distribution functions of all the incoming particles are given by equilibrium Fermi function determined by the reservoirs. In Eq.(154), three contributions to the electric current mean the distribution of right going electrons, of reflected left going electrons, and of electrons injected into the left side from the right side, respectively. The current is rewritten as

I^{nn} = AN(0)ev_{F} ∫ dE(1 − R_{N}^{nn}){f(E − eV) − f(E)} .  (155)

An intrinsic difference between the normal–superconducting junction and the N–N junction is the existence of the Andreev reflection in which an incident electron is reflected into a hole at the N–S interface. Since the reflected hole has positive charge, the electric current is enhanced by the existence of the Andreev reflection. Taking account of the Andreev reflection, Blonder et al.[10] obtained the electrical current I^{ns} through the normal– superconducting junction

I^{ns} = AN(0)ev_{F} ∫ dE(1 − R_{N}(E) + R_{A}(E))(f(E − eV) − f(E)) ,  (156)

where R_{A}(E) and R_{N}(E) are the Andreev reflection and the normal reflection coefficient, respectively. It is noted that the existence of the Andreev reflection increases the electric current. This expression was also obtained by Furusaki[11]. From the above equation, the differential conductance can be written as

\frac{dI}{dV} = \frac{1}{R(0)(1 − R_{N}(∞))}(1 − R_{N}(eV) + R_{A}(eV)) ,  (157)

where R(0) is the normal state resistance. For instance, in a completely transmissive contact, since R_{A} = 1 and R_{N} = 0 below the gap voltage (Δ_{bulk}/e), the differential conductance becomes twice the normal conductance.

We first discuss the Andreev reflection in the double infinite proximity contact system in Sec.3.1. And, taking account of a geometrical resonance effect in point contact device, we discuss the Andreev reflection coefficient and the differential conductance in Sec.3.2.

3.1 Andreev reflection in the normal–superconducting double infinite system

N  S 0  Z Δ(Z)

Figure 5: The normal–superconducting infinite double layer system.

First we discuss the Andreev reflection in the normal–superconducting ( N–S ) double infinite layer system[14][15]. The left ( right ) side region is assumed to be a normal– metal ( a superconductor ) as depicted in Fig.5. We consider a scattering problem for an incident electron with an incident energy E from the N region toward the N–S interface. Solving the Andreev equation(19) with appropriate boundary conditions, we can discuss the Andreev scattering. In the N region ( z < 0 ), since the pair potential will tend to zero at sufficiently large |z| even if the pairing interaction is non-zero, the asymptotic behavior of the N region wave function can easily be solved. ( If the pairing interaction g is zero, the pair potential Δ is also zero, since Δ = gF , where F is pair amplitude. If the pair potential is zero, the Andreev equation can easily be solved. ) The asymptotic behavior of the wave function Φ_{α}^{N} with the energy E are written as

Φ_{+}^{N}(z) ⟶^{z→∞} ( \begin{matrix}1 \\ 0 \end{matrix} ) e^{iκ^{N} z} + ( \begin{matrix}0 \\ r_{A} \end{matrix} ) e^{−iκ^{N} z},
Φ_{−}^{N}(z) ⟶^{z→∞} ( \begin{matrix}r_{N} \\ 0 \end{matrix} ) e^{−iκ^{N} z},  (158)

where κ^{N} = E/v_{F z}^{N}, r_{A} and r_{N} are the Andreev and the normal reflection amplitude, respectively. In the S–side ( z > 0 ) , the pair potential will tend to be the bulk value at sufficiently large distance. The asymptotic solutions of the wave function are written as

Φ_{+}^{S}(z) ⟶^{z→∞} c_{+} ( \begin{matrix}α_{∞} \\ β_{∞} \end{matrix} ) = c_{+} ( \begin{matrix}E + Ω \\ Δ_{bulk}^{S} \end{matrix} ) e^{iκ^{S}z},
Φ_{−}^{S}(z) ⟶^{z→∞} c_{−}ρ_{1} ( \begin{matrix}α_{∞} \\ β_{∞} \end{matrix} ) ,  (159)

where κ^{S} = Ω/v_{F z}^{S} ≡ \sqrt{E^{2} − Δ^{S}_{bulk}^{2}}/v_{F z}^{S} , ρ_{1} is the Pauli matrix and we have defined the square root to have positive imaginary part.

Since the Andreev equation is a first order differential equation, its formal solution can be written as

Φ_{α}(z) = ρ_{3}U_{α}(z, z′, E)ρ_{3}Φ_{α}(z′)  (160)

where U is the spatial evolution operator which we have defined in the previous chapter. Since the pair potential have been assumed to be real in this section, the evolution operator satisfies

U_{+}(z, z′, E) = ρ_{1}U_{−}(z, z′, E)ρ_{1}.  (161)

Using the asymptotic solutions Eqs.(158)(159), the interfacial boundary condition Eqs.(31) and the evolution operator, we can find the formal solution of the wave functions. One has only to solve a following equation :

U(0, −∞) ( \begin{matrix}Φ_{+}^{N}(−∞) \\ Φ_{−}^{N}(−∞) \end{matrix} ) = M̂_{0} ( \begin{matrix}Φ_{+}^{S}(0) \\ Φ_{−}^{S}(0) \end{matrix} )

= M̂_{0}U(0, ∞) ( \begin{matrix}Φ_{+}^{S}(∞) \\ Φ_{−}^{S}(∞) \end{matrix} ) ,  (162)

where

U(z, z′) ≡ ( \begin{matrix}ρ_{3}U_{+}(z, z′)ρ_{3} & 0 \\ 0 & ρ_{3}U_{−}(z, z′)ρ_{3} \end{matrix} ) = ( \begin{matrix}ρ_{3}U_{+}(z, z′)ρ_{3} & 0 \\ 0 & ρ_{2}U_{+}(z, z′)ρ_{2} \end{matrix} ) , (163)

where we have used Eq.(161). Solving this equation, it can be found that each amplitude coefficient of the wave function are generally

r_{A} = \frac{(αγ^{∗} + βδ^{∗})(αδ + βγ) − R_{0}(αγ + βδ)(αδ^{∗} + βγ^{∗})}{(αγ^{∗} + βδ^{∗})^{2} − R_{0}(αδ^{∗} + βγ^{∗})^{2}},  (164)
r_{N} = r_{0} \frac{(α^{2} − β^{2})}{(αγ^{∗} + βδ^{∗})^{2} − R_{0}(αδ^{∗} + βγ^{∗})^{2}} ,  (165)
c_{+} = d_{0} \frac{(αγ^{∗} − βδ^{∗})}{(αγ^{∗} + βδ^{∗})^{2} − R_{0}(αδ^{∗} + βγ^{∗})^{2}} ,  (166)
c_{−} = −r_{0}d_{0}^{∗} \frac{(α^{2} − β^{2})}{(αδ^{∗} + βγ^{∗})^{2} − R_{0}(αδ^{∗} + βγ^{∗})^{2}} ,  (167)

where

U_{+}^{N}(0, −∞) = ( \begin{matrix}γ & δ^{∗} \\ δ & γ^{∗} \end{matrix} ) ,  (168)
( \begin{matrix}α \\ β \end{matrix} ) = U_{+}^{S}(0, ∞) ( \begin{matrix}α_{∞} \\ β_{∞} \end{matrix} ) .  (169)

Using the time dependent Andreev equation,

i∂_{t}Φ_{α}(z) = ( \begin{matrix}−αiv_{F z}∂_{z} & Δ(z) \\ Δ(z) & αiv_{F z}∂_{z} \end{matrix} ) Φ_{α}(z),  (170)

the flux of the quasi–particles are given by

∑_{α=±} αv_{F z}Im (Φ_{α}^{†}ρ_{3}∂_{z}Φ_{α}).  (171)

Using this expression for the flux and noting that |α|^{2} − |β|^{2} is a conserved quantity in the Andreev equation, we obtain the reflection and transmission coefficients

R_{A} = |\frac{(αγ^{∗} + βδ^{∗})(αδ + βγ) − R_{0}(αγ + βδ)(αδ^{∗} + βγ^{∗})}{(αγ^{∗} + βδ^{∗})^{2} − R_{0}(αδ^{∗} + βγ^{∗})^{2}} |^{2}  (172)

Andreev reflection coefficient,

R_{N} = R_{0} |\frac{(α^{2} − β^{2})}{(αγ^{∗} + βδ^{∗})^{2} − R_{0}(αδ^{∗} + βγ^{∗})^{2}} |^{2}  (173)

Normal reflection coefficient,

T_{E} = (1 − R_{0})\frac{(|α|^{2} − |β|^{2})|αγ^{∗} − βδ^{∗}|^{2}}{|(αγ^{∗} + βδ^{∗})^{2} − R_{0}(αδ^{∗} + βγ^{∗})^{2}|^{2}}  (174)

transmission coefficient for a quasi–electron,

T_{H} = R_{0}(1 − R_{0})\frac{(|α|^{2} − |β|^{2})|αδ^{∗} + βγ^{∗}|^{2}}{|(αδ^{∗} + βγ^{∗})^{2} − R_{0}(αδ^{∗} + βγ^{∗})^{2}|^{2}}  (175)

transmission coefficient for a quasi–hole,

These coefficients satisfy the flux conservation law, i.e.,

R_{A} + R_{N} + T_{E} + T_{H} = 1.  (176)

Using this expression, the Andreev reflection coefficients can be calculated numerically. Typical results are shown in Figs.6 and 7. In Fig.6, we show the results when the N region has no pairing interaction. The self–consistent pair potential is plotted in Fig.6(a). The energy dependence of the Andreev reflection shows no qualitative difference from the BTK results, except for a small shift of the peak energy. The origin of the peak shift has been discussed in detail by Nagato et al.[15]. In Fig.7, the results when the N region has finite pairing interaction ( attractive and repulsive ). Due to the proximity effect, the normal region has a pair potential when the N region has a pairing interaction. The pair potential (a)

−5  0  5 0

0.2 0.4 0.6

R_{0} = 0.0

R_{0} = 0.2

R_{0} = 0.5

Z/ξ

Δ(Z)/πT_{C}^{S}

(b)

0 1 1

R_{0} = 0.0

R_{0} = 0.2

R_{0} = 0.5

E/Δ_{bulk}

R_{A}

Figure 6: The R_{0} dependence of (a) the pair potential and of (b) the Andreev reflection for the N–S infinite double layer. Temperature T = 0.2T_{C}^{S}, the critical temperature of the normal region t_{C}^{N} ≡ T_{C}^{N}/T_{C}^{S} = 0.0. The superconducting coherence length ξ ≡ v_{F} /πT_{C}^{S}.

(a)

−10  −5  0  5

0 0.2 0.4 0.6

t_{C}^{N} = 0.01
t_{C}^{N} = 0.0
g^{N}/g^{S} =−0.5

Z/ξ

Δ(Z)/πT_{C}^{S}

(b)

0 1 1

t_{C}^{N} = 0.0

t_{C}^{N} = 0.01

g^{N}/g^{S} =−0.5

E/Δ_{bulk}

R_{A}

Figure 7: The pairing interaction dependence of (a) the pair potential and of (b) the Andreev reflection for the N–S infinite double layer. Temperature T = 0.2T_{C}^{S}, the interfacial reflection coefficient R_{0} = 0.2. of the normal region increases as the interfacial reflection R_{0} decreases and as the critical temperature of the N increases. When the pairing interaction of the N is repulsive, the N metal has the pair potential of opposite sign to that of the S. In the case R_{0} = 0, the Andreev reflection coefficient R_{A} at the lower energies than Δ_{bulk}, is unity. The R_{A} decreases as R_{0} increases at lower energies. When the pair interaction is repulsive, the Andreev reflection at low energy is enhanced. When the pair interaction is attractive, it is suppressed inversely. In particular, at the zero energy limit, the Andreev reflection coefficient R_{A} can be estimated to be

R_{A}(0) = |\frac{1 − R_{0} − (1 + R_{0}) tanh 2δ}{1 + R_{0} − (1 − R_{0}) tanh 2δ} |^{2},  (177)
δ = \frac{1}{v_{F z}^{N}} ∫_{0}^{−∞} dzΔ^{N}(z).  (178)

The Andreev reflection R_{A} at lower energies is given in terms of R_{0} and of the integral of the pair potential over the entire region of the N side. By use of Eq.(157), the differential conductances corresponding to the cases in Figs.6 and 7 are shown in Fig.8. We have plotted the normalized differential conductance K such that tends to unity at the high energy limit, i.e.,

K(eV) ≡ \frac{dI/dV}{dI/dV |_{eV →∞}} = \frac{1 − R_{N}(eV) + R_{A}(eV)}{1 − R_{0}}  for the N–S case.  (179)

Before concluding this subsection, we make some comments on the wave packet picture of the Andreev reflection proposed in Refs. [11], [13] and [7] to explain the stair-like structure observed in the differential conductance. According to this picture, an electron wave packet with energy E changes the sign of its group velocity at the position that satisfies E = Δ(z), thus is converted to a hole wave packet. It follows that, if some structure in the energy dependence of the Andreev reflection is found, the structure can be assigned to singular (a)

0 1 2 4

R_{0} = 0.0
R_{0} = 0.2
R_{0} = 0.5
eV/Δ_{bulk}

K(eV) (b)

0 1 1 2

t_{C}^{N} = 0.0

t_{C}^{N} = 0.01
g^{N}/g^{S} =−0.5
eV/Δ_{bulk}

K(eV)

Figure 8: The differential conductance (a) and (b) for the cases shown in Figs. 6 and 7 , respectively. points in Δ(z), say Δ(0). The peak around Δ(0) may be interpreted in this way. They suggested that the conductance enhancement at lower energies is due to the wave packet reflection by the induced pair potential in the N region. Such a picture, however, cannot be applied to the N-S contact system, in particular at low temperatures. In order to form a wave packet of the Bogoliubov quasi-particle which has a group velocity with definite sign, the range of the momenta involved should be sufficiently small, i.e., less than T_{C}/v_{F}. As a result, the width of the wave packet becomes longer than the coherence length ξ = v_{F} /πT_{C}. On the other hand, the depression range of the pair potential Δ(z) is of order ξ as can be seen from Fig. 7. It is difficult, therefore, to say at which point the wave packet is reflected. In fact, the present analysis shows that the induced pair potential in the N region does not lead to a significant enhancement of the Andreev reflection at lower energies.

3.2 Normal–normal–superconducting proximity contact system

The experiments to observe the Andreev reflection are usually performed using a point contact system. A needle-like electrode is attached to the normal layer deposited on an semi-infinite superconductor. A typical geometry of the system is shown in Fig.9. Some of the experimental results of the differential conductance shows a stair-like structure in the voltage dependence, which can not be explained by the BTK framework discussed in the last section. There are some attempts to explain the stair-like structure by taking into account the possible finite pair potential in the N-side. But we have shown in the last section that the effects by the finite pair potential is not so significant as to change the structure in the voltage dependence.

We instead are interested in the effect by the finite reflection of electrons at the point contact. Since the voltage drop is expected to occur at the point contact, electrons and holes which traverse the normal layer will be reflected not only by the N-S interface but also by the point contact. The multiple reflection within the normal layer is expected to lead to some geometrical resonance effects in the differential conductance.

We consider a model normal-normal-superconducting triple layer system as shown in Fig.10, to study the geometrical resonance effects. For simplicity, we assume the electrode to be a normal metal. The pairing interaction g in the normal metal is not necessarily g = 0 but can be g ≠ 0.

Let us consider the Andreev scattering in the normal–normal–superconducting triple layer ( N–N′–S ). Using the same methods as in the N–S case discussed in Sec.3.1, we can find the formal solution of the wave functions. We have only to solve a following equation,

electrode  _{N}  _{S}

Figure 9: The point contact device.

N  _{N′}  _{S} ( electrode )  _{0} _{L}  _{Z}_{Δ(Z)}_{R_{0}}  _{R_{L}}

Figure 10: The normal–normal–superconducting ( N–N′–S ) proximity contact system. instead of Eq.(162),

U(0, −∞) ( \begin{matrix}Φ_{+}^{N}(−∞) \\ Φ_{−}^{N}(−∞) \end{matrix} ) = M̂_{0} ( \begin{matrix}Φ_{+}^{N′}(0) \\ Φ_{−}^{N′}(0) \end{matrix} )
= M̂_{0}U(0, L) ( \begin{matrix}Φ_{+}^{N′}(L) \\ Φ_{−}^{N′}(L) \end{matrix} )
= M̂_{0}U(0, L)M̂_{L} ( \begin{matrix}Φ_{+}^{S}(L) \\ Φ_{−}^{S}(L) \end{matrix} )
= M̂_{0}U(0, L)M̂_{L}U(L, ∞) ( \begin{matrix}Φ_{+}^{S}(∞) \\ Φ_{−}^{S}(∞) \end{matrix} ) ,  (180)

where Φ_{±}^{N}(−∞) and Φ_{±}^{S}(∞) are defined by Eqs.(158) and (159). From this equation, we can obtain the amplitudes of the wave function.

r_{A} = {(ν^{∗}β + μα)(μ^{∗}β + να) − R_{0}(τ^{∗}β + σα)(σ^{∗}β + τα)
−R_{L}(ν^{∗}α + μβ)(μ^{∗}α + νβ) + R_{0}R_{L}(τ^{∗}α + σβ)(σ^{∗}α + τβ)
+2\sqrt{R_{0}R_{L}}(η^{∗}ζ − ηζ^{∗})(α^{2} − β^{2}) cos ϕ_{C}\}/D,  (181)
r_{N} = e^{iθ_{r_{0}}}\{\sqrt{R_{L}}(e^{iϕ_{C}} + R_{e}^{−iϕ_{C}})
+\sqrt{R_{0}}[α^{2}(γ^{2} − δ^{2}) − β^{2}(γ^{∗}^{2} − δ^{∗}^{2}) + 2αβ(γδ^{∗} − γ^{∗}δ)]
+\sqrt{R_{0}}R_{L}[α^{2}(γ^{∗}^{2} − δ^{∗}^{2}) − β^{2}(γ^{2} − δ^{2}) + 2αβ(γ^{∗}δ − γδ^{∗})]\}/D, (182)
c_{+} = d_{0}d_{L}(να + μ^{∗}β + \sqrt{R_{0}R_{L}}e^{−iϕ_{C}}(σβ + τ^{∗}α))/D,  (183)
c_{−} = −r_{L}d_{0}d_{L}^{∗}(ν^{∗}α + μβ + \sqrt{\frac{R_{0}}{R_{L}}} e^{−iϕ_{C}}(σα + τ^{∗}β))/D,  (184)
D = {(ν^{∗}β + μα)^{2} − R_{0}(τ^{∗}β + σα)^{2}
−R_{L}(ν^{∗}α + μβ)^{2} + R_{0}R_{L}(τ^{∗}α + σβ)^{2}
+2\sqrt{R_{0}R_{L}}(μτ^{∗} − ν^{∗}σ)(α^{2} − β^{2}) cos ϕ_{C}\},  (185)

where

μ = ζ^{∗}γ − η^{∗}δ, ν = ζδ − ηγ,

τ = ζγ − ηδ, σ = ζ^{∗}δ − η^{∗}γ,
U_{+}^{N}(0, −∞) = ρ_{3} ( \begin{matrix}ζ & η^{∗} \\ η & ζ^{∗} \end{matrix} ) ρ_{3},  (186)
U_{+}^{N′}(0, L) = ρ_{3} ( \begin{matrix}γ & δ^{∗} \\ δ & γ^{∗} \end{matrix} ) ρ_{3},  (187)
( \begin{matrix}α \\ β \end{matrix} ) = ρ_{3}U_{+}(L, ∞)ρ_{3} ( \begin{matrix}α_{∞} \\ β_{∞} \end{matrix} ) .  (188)

The above equations contain now the phase factor e^{iϕ_{C}}. This describes the multiple reflection processes in the finite layer. Apart from a correction of order 1/p_{F} L the Andreev scattering coefficients are obtained by averaging over ϕ_{C}. All the coefficients of the Andreev scattering are obtained in the forms

R_{A} = \frac{1}{2π} ∫_{0}^{2π} dϕ_{C}|r_{A}|^{2},  (189)
R_{N} = \frac{1}{2π} ∫_{0}^{2π} dϕ_{C}|r_{N}|^{2},  (190)
T_{E} = \frac{1}{2π} ∫_{0}^{2π} dϕ_{C} \frac{v_{F z}^{S}}{v_{F z}^{N}} (|α|^{2} − |β|^{2})|c_{+}|^{2},  (191)
T_{H} = \frac{1}{2π} ∫_{0}^{2π} dϕ_{C} \frac{v_{F z}^{S}}{v_{F z}^{N}} (|α|^{2} − |β|^{2})|c_{−}|^{2},  (192)

(193)

where we note that |α_{z}|^{2} − |β_{z}|^{2} is a quasi–particle flux conserved in the Andreev equation. These coefficients satisfy the flux conservation law, i.e.,

R_{A} + R_{N} + T_{E} + T_{H} = 1.  (194)

Typical results are shown in Figs.11. The differential conductance K normalized by its high energy limit

K(eV) = \frac{1 − R_{0}R_{L}}{(1 − R_{0})(1 − R_{L})}\{1 − R_{N}(eV) + R_{A}(eV)\}.  (195)

(a)

0 0.2 0.4

R_{0} = 0.2
R_{0} = 0.5
R_{0} = 0.8

R_{L} = 0.5

−4  −2  0  2  4
Z/ξ

Δ(Z)/πT_{C}^{S}

(b)

0  0.4  0.8  1.2

0.2 0.4 0.6

E/Δ_{bulk}

R_{A}

Figure 11: The typical result of (a) the pair potential, of (b) the Andreev reflection and of (c) the differential conductance for the N–N′–S system. The Andreev reflection as a function of the incident energy E has double peak structure. T = 0.2T_{C}^{S}, the scaled N′ layer size l ≡ L/ξ = 1.0 and t_{C}^{N} = t_{C}^{N′} = 0.01.

0  0.4  0.8  1.2

1 2

3 R_{L} = 0.5

R_{0} = 0.2
R_{0} = 0.5

R_{0} = 0.8

eV/Δ_{bulk}^{S}

K(eV)

is plotted. The Andreev reflection coefficient and the differential conductance in the energy range below Δ_{bulk} has double peak structure when both the interfaces have finite reflection. When either reflection coefficient ( R_{0} or R_{L} ) vanishes, the system is equivalent to that of the infinite N-S system discussed in the last section.

The double peak structure can be interpreted as follows. In the N′ region, a right going electron ( hole ) is partially reflected into a hole ( electron ) at the N′-S interface by the Andreev process. The reflected hole ( electron ) is reflected into a hole ( electron ) at the N-N′ interface by the normal reflection process. The hole ( electron ) is again converted into an electron ( hole ) at the N′-S interface and so on. This sequence leads to the de Gennes–Saint-James bound state which will be discussed in detail in the next section. The bound state is broadened by the finiteness of the reflection at both the interfaces. The double structure found in the Andreev reflection can be interpreted as resonance states formed in a valley between the partially reflecting interfaces.

Kieselmann[33] has calculated the density of states of the finite width normal layer on the semi–infinite superconductor by use of the conventional quasi–classical Green’s function method. He showed that the density of states also has double peak structure at the lower energies than Δ_{bulk}. This double peak is also related to the broadened de Gennes–Saint James bound state. Ashida et al.[20] have shown the existence of the similar bound states in the finite N–S double layer. van Son et al.[13] have studied the case when R_{0} ≠ 0 and R_{L} = 0 using a model pair potential.

The peak structure shifts to the lower energy as the N′ layer size L increases as is shown in Fig.12. When the layer size L is sufficiently large, there appear several double peaks.

We consider a role by the pairing interaction in the normal region for the probability of Andreev reflection. Figure 13 is for the case that t_{C}^{N} = 0 and t_{C}^{N′} ≠ 0. Energy level of the virtual bound state in the N′ is lifted by the pair potential in the N′. The two peaks of R_{A} are shifted to higher energy as t_{C}^{N′} increases. Figure 14 is for the case that t_{C}^{N} ≠ 0 and t_{C}^{N′} = 0. In a similar manner to the N–S infinite double layer system, the probability of the Andreev reflection in low energy region is slightly enhanced as t_{C}^{N} increases. But the effect by finite t_{C}^{N} or t_{C}^{N′} is not so large as to change the structures.

In the one point contact experiment, the differential conductance has a stair-like structure with respect to an incident energy[6][7]. In the N–N′–S system with finite reflection coefficients, we obtained the probability of Andreev reflection having two peak as a function of an incident energy. However, the enhancement of the differential conductivity at zero bias in the above experiments can not be explained. We can not reproduce the stair-like structure found in the point contact system. Nevertheless, we expect that the double peak structure will be observed in the coplanar geometry. We have also shown that the finite pairing interaction in the N′ layer does not explain the stair-like structure in contrast to the claim by Ref.[11], [13] and [7].

(a) 0 1 0.2 0.4 0.6

l = 1

l = 2

l = 5
E/Δ_{bulk}^{S}

R_{A} (b)

0 1 1 2 3

l = 1

l = 2

l = 5
eV/Δ_{bulk}^{S}

K(eV)

Figure 12: The N′ layer size L dependence of (a) the Andreev reflection, and of (b) the differential conductance. T = 0.2T_{C}^{S}, R_{0} = R_{L} = 0.5 and t_{C}^{N} = t_{C}^{N′} = 0.01.

We did not consider a possible supercurrent flowing across the N′-S interface. When the normal metal has a pairing interaction, there exists a pair potential in the normal region as was shown above. It is possible that Josephson like supercurrent may flow across the N-S interface in addition to the excess current carried by the Andreev reflected hole. The excess low voltage conductance observed in some experiments[16] may be interpreted as being due to the supercurrent. In order to discuss this possibility in a unified way, it is necessary to develop a non-equilibrium quasi-classical formulation. This is an interesting problem to be further studied.

(a)

0 0.2 0.4

t_{C}^{N′} = 0.0 t_{C}^{N′} = 0.01

g^{N′}/g^{S} =−0.5

−4  −2  0  2  4
Z/ξ

Δ(Z)/πT_{C}^{S}

(b)

0  0.4  0.8  1.2

0.2 0.4

E/Δ_{bulk}^{S}

R_{A}

Figure 13: The N′ layer has a pairing interaction. (a) The pair potential, (b) the Andreev reflection, and (c) the differential conductance. T = 0.2T_{C}, t_{C}^{N} = 0.0, l = 1.0 and R_{0} = R_{L} = 0.5.

(c)

0  0.4  0.8  1.2

1 2

t_{C}^{N′} = 0.01

t_{C}^{N′} = 0.0

g^{N′}/g^{S} = −0.5 eV/Δ_{bulk}^{S} K(eV)

(a)

0 0.2 0.4 0.6

t_{C}^{N} = 0.0

t_{C}^{N} = 0.01

g^{N}/g^{S} = −0.5

−4  −2  0  2  4
Z/ξ

Δ(Z)/πT_{C}^{S}

(b)

0  0.4  0.8  1.2

0.2 0.4

E/Δ_{bulk}^{S}

R_{A}

Figure 14: The N layer has a pairing interaction. (a) The pair potential, (b) the Andreev reflection, and (c) the differential conductance. T = 0.2T_{C}, t_{C}^{N′} = 0.0, l = 1.0 and R_{0} = R_{L} = 0.5.

(c)

0  0.4  0.8  1.2

1 2

t_{C}^{N} = 0.01 t_{C}^{N} = 0.0 g^{N}/g^{S} = −0.5

eV/Δ_{bulk}^{S}

K(eV)

4 The local density of states of the superconducting– normal–superconducting triple layer system

Recently, in the superconducting–normal metal–superconducting ( S–N–S ) proximity contact system the local density of states was studied using a scanning tunneling microscope by Inoue and Takayanagi[36]. Their data indicates that the local density of states in the normal region has an effective energy gap and depends on a distance from the S–N interface.

Tanaka et al.[34] studied the density of states of the S–N–S semi–infinite system with the self–consistently solved pair potential but with ideal interfaces by numerically solving the Bogoliubov–de Gennes equation. They showed that the density of states spatially varies when there exists the pairing interaction in the N region. Also, Hara et al.[35] investigated the density of states in the N–S finite double layer system in which the N–S interface has finite reflection but the N region has no pairing interaction. They showed that the density of states has an energy gap–like structure, although their model does not have a pair potential because having no pairing interaction in the N region. We study the density of states of the S-N-S system having realistic boundaries based on the self–consistent pair potential.

4.1 The de Gennes–Saint-James bound state

We first consider the bound state, which was discussed by de Gennes and Saint-James[28], in a finite width normal layer connected to a bulk superconductor as depicted in Fig.15. It was assumed that the pair potential is spatially constant and the N–S interface is perfectly transmissive. They considered the states below the superconducting pair potential energy Δ.

0  L Δ  S N  Z

Figure 15: A normal metal connected to a bulk superconductor with a constant pair potential. The width of the normal metal is L. The normal–superconducting interface is assumed to be perfectly transmissive.

Within the Andreev approximation, general solutions Ψ of the Bogoliubov de Gennes

equation Eq.(13) with the energy E can be easily solved as

Ψ^{N}(z) = Φ_{+}^{N}(z)e^{ip_{F z}z} + Φ_{−}^{N}(z)e^{−ip_{F z}z}
= (( \begin{matrix}a \\ 0 \end{matrix} ) e^{iκ^{N} z} + ( \begin{matrix}0 \\ b \end{matrix} ) e^{−iκ^{N} z}) e^{ip_{F z}z} + (( \begin{matrix}c \\ 0 \end{matrix} ) e^{−iκ^{N} z} + ( \begin{matrix}0 \\ d \end{matrix} ) e^{iκ^{N} z}) e^{−ip_{F z}z},(196)
Ψ^{S}(z) = Φ_{+}^{S}(z)e^{ip_{F z}z} + Φ_{−}^{S}(z)e^{−ip_{F z}z}
= (e ( \begin{matrix}Δ \\ E − iS \end{matrix} ) e^{−κ^{S}z}) e^{ip_{F z}z} + (f ( \begin{matrix}Δ \\ E + iS \end{matrix} ) e^{−κ^{S}z}) e^{−ip_{F z}z},  (197)
κ^{N} = E/v_{F z},  (198)
κ^{S} = S/v_{F z},  (199)
S = \sqrt{Δ^{2} − E^{2}},  (200)

where a, b, c, d, e and f are constants and Φ_{+}^{S} and Φ_{−}^{S} have been chosen to tend to zero at

sufficiently large z. On the other hand, the boundary conditions at the normal layer end z = 0 and at the N–S interface z = L are written as

Φ_{+}^{N}(0) + e^{iη_{0}}Φ_{−}^{N}(0) = 0,  (201)
Φ_{±}^{N}(L) = Φ_{±}^{S}(L),  (202)

where η_{0} is a phase shift at the N layer end. From the general solutions and the boundary conditions, one can get an equation to determine the bound states.

(nπ + φ) cos θ = \frac{2L}{ξ_{0}} cos φ n : integer
cos φ = \frac{E}{Δ},  (203)

where ξ_{0} = v_{F} /Δ. Consider first the case θ = 0 (propagation along z). For fixed L, equation (203) has a finite number of solutions φ_{n}(E_{n}(0)) below the energy gap Δ. The number of solutions increases with L/ξ_{0}. There is always at least one solution below Δ ( even when L/ξ_{0} ≪ 1 ).

The total density of states including all the polar angles θ of the Fermi momenta is defined as

n(ω) = \frac{A}{(2π)^{2}} ∫ dk_{x}dk_{y} ∑_{n} δ(ω − E_{n}(θ))  (204)

where A is the area of the interfacial surface. Performing integration of k_{x}, k_{y} , one can obtain

n(ω) = 2N(0)π \frac{2L}{ξ_{0}} ∑_{n=n_{c}} cos φ\frac{1}{(nπ + φ)^{2}} (1 + \frac{cot φ}{nπ + φ}) ,  (205)
cos φ = \frac{ω}{Δ},

where n_{c} is a lower limit of n–summation, i.e., n_{c} is the smallest integer which satisfies

nπ + φ ≥ \frac{2L}{ξ_{0}} cos φ.

de Gennes and Saint-James obtained the density of states with discontinuities at ε = E_{n}(0) ( Fig.16 ).

0.5  1 048  y = 2L/ξ_{0} y = 2.5 y = 8.0 y = 20.0 ω/Δ n(ω)/N(0)

Figure 16: The total density of states below the gap energy Δ in the normal region of the system shown in Fig. 15 .

When \frac{ω}{Δ} ≪ 1 and \frac{ω}{Δ} \frac{L}{ξ}_{0} ≪ 1, the density of states n(ω) is

n(ω) ∼ 2N(0)π \frac{L}{ξ_{0}} \frac{ω}{Δ}.

The density of states is linear with respect to the energy and to the normal layer size around the low energy. Similar results was obtained by Hara et al.[35] in the finite normal– superconducting double layer solved self–consistently with a finite interfacial reflection.

When the layer size L is large, there occur several discontinuities in the density of states as is shown in Fig.16. When \frac{ω}{Δ} ≪ 1 and \frac{ω}{Δ} \frac{L}{ξ}_{0} ≫ 1, the density of states is

n(ω) ∼ 2N(0) .

In this case, the density of states becomes the normal state density of states for both projections.

It is found that the density of states in the normal region with a superconductor has energy spectrum below the gap energy Δ because of the finiteness of the normal layer.

4.2 The superconducting–normal–superconducting semi–infinite triple layer system

In this section, using a self–consistently solved pair potential, we discuss the density of states in the S–N–S system with finite interfacial reflection as depicted in Fig.17. We assume R ≡ R_{0} = R_{L}. We also take account of the pairing interaction in the N region. In terms of the diagonal elements of the Green’s function for the real frequency ε = ω, the local density of states associated with the Fermi momentum p_{F}^{α} are defined by

n(ω, p_{F}^{α} , z) ≡ ∑_{l} δ(ω − E_{l})Φ_{αl}(z)Φ_{αl}^{†}(z)
= \frac{1}{4πiv_{F z}} ρ_{3} [ĝ_{αα}(ω + i0, z) − ĝ_{αα}(ω − i0, z)] .  (206)

From the expression for the quasi–classical Green’s function Eqs.(79)–(88), the Green’s function can be in general written in such a form as

ĝ = i(2X − trX) = i ( \begin{matrix}a & b \\ c & −a \end{matrix} ) .

where X is a matrix. For the real frequency ε = ω, the evolution operator satisfies

U_{α}(ω + i0) = ρ_{1}U_{α}(ω − i0)^{∗}ρ_{1} .  (207)

From this equation, it can be found that the real frequency Green’s function satisfies

ĝ_{αα}(ω + i0) = −ρ_{1}ĝ_{αα}(ω − i0)^{∗}ρ_{1} .  (208)

−2  0  2  4 0 0.2 0.4 0.6  R =0.5 0.2 0.0 S  S N Z/ξ Δ(Z)/πT_{C}^{S}

Figure 17: The pair potential for the superconducting–normal–superconducting semi– infinite triple layer. T = 0.2T_{C}^{S}, t_{C}^{N} = 0.0, l ≡ L/ξ = 1.0 and ξ = v_{F} /πT_{C}^{S}.

1  2 024  R = 0.0 0.2 0.5 ω/Δ_{bulk}^{S} 2πv_{F} n^{N}(ω)

Figure 18: The density of states as a function of the energy ω for the case shown in Fig. 17 . This density of states is for the Fermi momentum normal to the interface, so called the tunneling density of states.

Then, ĝ(ω + i0) − ĝ(ω − i0) is written in a form

ĝ_{αα}(ω + i0) − ĝ_{αα}(ω − i0) = i ( \begin{matrix}a & b \\ c & −a \end{matrix} ) − iρ_{1} ( \begin{matrix}a & b \\ c & −a \end{matrix} )^{∗} ρ_{1}
= i ( \begin{matrix}a + a^{∗} & b − c^{∗} \\ c − b^{∗} & −a − a^{∗} \end{matrix} ) .  (209)

To obtain the density of states, therefore, we have only to treat the Green’s function for the real frequency ε = ω + i0 and we obtain

n(ω, p_{F}^{α} , z) = \frac{1}{2πv_{F z}} ρ_{3}Imĝ_{αα}(ω + i0).  (210)

Using the self–consistently solved pair potential, we numerically calculated the density of states when the polar angle θ = 0. When the interfacial reflection coefficients are finite, the density of states has double peak as a function of the energy ω as is shown in Fig.18. In some systems this double peak structure of the density of states has been reported theoretically. Kieselmann[33] obtained the density of states in a N–S contact of a finite normal layer with a semi–infinite superconductor. Ashida et al.[20] and Hara et al.[35] obtained the density of states in the finite N–S contact system.

The density of states at the lower energies than Δ_{bulk} are related to the de Gennes– Saint-James bound states[28] discussed in the previous section. The de Gennes–Saint James bound states for the polar angle θ = 0 becomes considerably broad when the interfacial reflection coefficients are finite as is shown in Fig.18. The width of the broadened bound states increases with the reflection coefficient R.

When the pairing interaction of the normal region is zero, it can be easily shown that the density of states of the N region is spatially constant as follows : the spatial variation of the density of states is expressed by the evolution operator U, i.e.,

n(z) = \frac{1}{2πv_{F z}} ρ_{3}Im U_{±}(z, z′)ĝ_{±±}(z′)U_{±}(z′, z) .

The U_{±} of the normal region without pair potential is written as

U_{±}^{N}(z, z′) = ( \begin{matrix}e^{±iω(z−z′)/v_{F z}} & 0 \\ 0 & e^{∓iω(z−z′)/v_{F z}} \end{matrix} ) .  (211)

Then, in the normal region the density of states does not spatially change. On the other hand, if the pairing interaction is finite in the normal region, the density of states spatially changes, but its peak position does not change spatially. Figure 19 shows the density of states at the center of the N region when the N metal has a pairing interaction. When the pairing interaction is attractive (repulsive), the peak level of the density of states shifts to higher (lower) energy side.

The bound state level shifts down as increasing the normal layer size L in Fig.20. It corresponds to the usual bound state problem of an electron in a potential well. When the layer size is sufficiently large, several broad bound states appear, in a similar manner to the de Gennes–Saint-James’ result shown in the previous subsection.

It has been predicted that the density of states of the normal region have energy gap at the lower energy than Δ_{bulk}. This state can be detected by use of a tunneling experiment, a scanning tunneling microscope ( STM ) measurement and so on. In the above, we calculated the density of states for the polar angle θ = 0 of the Fermi momentum, so called ”the tunneling density of states” which is detected by the tunneling experiment[17]. The STM experiment, however, can detect rather ”the total density of states” which is averaged over all direction of the Fermi momentum, than the tunneling density of states. Next, we calculate the total density of states in the normal region.

The total density of states is defined by

n_{total}(z) ≡ N(0) ∫_{0}^{π/2} dθ sin θ Imĝ(z, ω + i0).  (212)

The total density of states was first obtained by de Gennes and Saint-James in the ideal (a)

−2  0  2  4

0 0.2 0.4 0.6

t_{C}^{N} = 0.01

t_{C}^{N} = 0.0 g^{N}/g^{S} =−0.5

Z/ξ

Δ(Z)/πT_{C}^{S} (b)

0  0.5  1  1.5

5 10

t_{C}^{N} = 0.0 t_{C}^{N} = 0.01 g^{N}/g^{S} =−0.5

ω/Δ_{bulk}^{S}

2πv_{F} n^{N}(ω)

Figure 19: (a) The pair potential and (b) the N region tunneling density of states when the N metal has a pairing interaction. T = 0.2T_{C}^{S}, l = L/ξ = 1.0 and R = 0.2.

0.5  1  1.5 05 10  l = 1 l = 2  ω/Δ_{bulk}^{S} 2πv_{F} n^{N}(ω)

Figure 20: N layer size dependence of the density of states. l = L/ξ, T = 0.2T_{C}^{S}, t_{C}^{N} = 0.0 and R = 0.2.

N–S contact system ( see the previous section. ) Recently, the total density of states was calculated by Hara et al.[35] in the finite N–S contact with finite interfacial reflection coefficient and by Tanaka et al.[34] in the semi–infinite S–N–S contact without the interfacial reflection coefficients.

Typical result is shown in Fig.21. In the normal region the particle with the polar angle θ ≠ 0 moves effectively the more long distance than the normal layer width L. Since the total density of states includes contributions from all the polar angles of the Fermi momentum, it includes the states with lower energy than the lowest energy level of the tunneling density of states. The total density of states has a cusp and a shoulder. The cusp and the shoulder are located at the two peak positions in the tunneling density of states.

The finite density of states in the energy range below Δ_{bulk} is originated from the de Gennes–Saint-James bound state in the N–S system discussed in the previous subsection.

0.5  1 0 2.55  ω/Δ_{bulk}^{S} n_{total}(ω)/N(0)

Figure 21: The total density of states including all the Fermi momentum components in the N region. A solid line shows the total density of states and a dashed line shows the tunneling density of states. T = 0.2T_{C}^{S}, t_{C}^{N} = 0.0, l = L/ξ = 1.0 and R = 0.2.

0.5  1 01234  t_{C}^{N} = 0.0 t_{C}^{N} = 0.1 g^{N}/g^{S} = −0.5 ω/Δ_{bulk}^{S} n_{total}(ω, Z =0.5ξ)/N(0)

Figure 22: The total density of states when the N metal has a pairing interaction. Z = 0.5ξ, T = 0.2T_{C}^{S}, l = L/ξ = 1.0 and R = 0.2.

In the N region, a particle is reflected into a hole at the N–S interface by the Andreev reflection. The backward hole is reflected into an electron at the S–N interface. The de Gennes–Saint James bound state is nothing but this closed sequence of particles and holes. When the interfacial reflection coefficient is zero, the phase of the wave function earned along the closed path is rather small because of a cancellation between the particle path and the hole path. When the interfacial reflection coefficient is not zero, however, a closed path including only particle and/or only hole is mixed. The latter path earns a phase an amount of 2p_{F} L cos θ. When R is finite, therefore, the de Gennes–Saint James bound state energy becomes very sensitive to the polar angle of the Fermi momentum. It forms a band with width of order v_{F} π/L when contributions from very small polar angle range | cos θ − 1| ≤ π/p_{F} L are taken into account.

When the normal region has the attractive pairing interaction, the density of states in the normal region has an energy gap as is shown in Fig.22. When the normal region has the repulsive pairing interaction, the energy level shifts to the lower energy side. For any pairing interaction, it can be analytically shown that the density of states at ω = 0 is zero. When ω = 0, the evolution operator can be written as

U_{+}(z, z′) = ( \begin{matrix}cosh δ_{z,z′} & −i sinh δ_{z,z′} \\ i sinh δ_{z,z′} & cosh δ_{z,z′} \end{matrix} ) ,  (213)
δ_{z,z′} = −\frac{1}{v_{F z}} ∫_{z′}^{z} dz′′Δ(z′′).  (214)

Using this expression, we can obtain the quasi–classical Green’s function in the C layer

ĝ_{++}^{C}(z, ω = 0) = i\frac{2C(z) − trC(z)}{\sqrt{(trC(z))^{2} − 4R_{0}R_{L}}} ,  (215)
C(z) = ((1 + R_{0}R_{L}) cosh 2δ_{0,L} + (1 − R_{0}R_{L}) sinh 2δ_{0,L}) ρ_{0}
+ ((1 − R_{0}R_{L}) cosh 2δ_{0,L} + (1 + R_{0}R_{L}) sinh 2δ_{0,L}) ρ_{2}.  (216)

The Green’s function with ω = 0, therefore, does not depend on the position. It is also shown that the density of states, which is obtained from the above Green’s function, at ω = 0 is zero for any pair potential. This result is different from numerical results obtained by Tanaka et al.[34]. They calculated the total density of states in the S–N–S system with ideal interfaces by numerically solving the Bogoliubov–de Gennes equation. Their results show that when the N metal has a repulsive interaction, the density of states at ω = 0 is not zero and depends on the position. In contrast to their numerical calculation, we can explicitly show that the density of states at ω = 0 is zero and does not depend on the position.

As we have noted, when the pairing interaction of the normal metal is zero, the density of states of the normal metal is spatially constant. When the pairing interaction is not zero, however, the density of states spatially varies as is shown in Fig.23. The density of states is larger at the point of which the pair potential is smaller. Also, when the layer size L of the central normal metal is large, the density of states has several peaks at the lower energy than Δ_{bulk} as is shown in Fig.24.

Recently, Inoue and Takayanagi studied the local density of states of Nb/InAs/Nb system by use of STM. It was reported that the local density of states in the InAs region has an energy gap. Also, It was reported that the local density of states depends on a position in the InAs region. In the present model, when the normal metal has a pairing interaction, the density of states depends on the position in the normal region. The InAs, therefore, may have a finite pairing interaction. Since a realistic system is not in the clean limit, however, it seems necessary to take account of the impurity effect for more quantitative comparison with the experiments.

0.5  1 01234  Z = 0ξ Z = 0.16ξ Z = 0.5ξ ω/Δ_{bulk}^{S} n_{total}(ω)/N(0)

Figure 23: Position dependence of the N region total density of states when the N metal has a pairing interaction t_{C}^{N} = 0.1. T = 0.2T_{C}^{S}, l = L/ξ = 1.0 and R = 0.2.

0.5  1 024  ω/Δ_{bulk}^{S} n_{total}(ω, Z =5ξ)/N(0)

Figure 24: The density of states when the N layer size L is long, L/ξ = 10.0. A solid line shows the total density of states and a dashed line shows the tunneling density of states. Z = 5ξ, T = 0.2T_{C}^{S}, t_{C}^{N} = 0.1 and R = 0.2.

5 The infinite double layer system with s–wave and d–wave superconductor

So far we discussed s–wave superconducting proximity contact systems. In this last section, we study the proximity contact system which consists of superconductors with different symmetry.

It is believed that ”heavy fermion superconductor” and ”high Tc superconductor”, which were recently discovered, are classified as superconductors with unconventional pairing. Pals et al.[39] investigated the Josephson contact between a singlet and a triplet superconductor using the tunneling Hamiltonian. Poppe[40] observed a Josephson current between an s–wave superconductor Al and a heavy fermion superconductor CeCu_{2}Si_{2}. Ashauer et al.[41] studied theoretically a thin film of standard superconductor in proximity contact with a bulk unconventional material by use of the conventional quasi–classical method.

We study the double infinite proximity contact system between superconductors with different pairing interactions, that is, with s–wave pairing and d–wave pairing interaction as depicted in Fig.25. Since both the superconductors are singlet superconductor, we have only to treat 2 × 2 matrices. The d–wave pairing interaction is expanded in terms of the spherical harmonic functions Y_{2m},

v_{d}(p_{F} · p′_{F}) = −5g_{2}P_{2}(p_{F} · p′_{F})
= −4πg_{2} ∑_{m=−2}^{m=2} Y_{2m}(θ, ϕ)Y_{2m}^{∗}(θ′, ϕ′),  (217)

where θ, ϕ are the polar angles of the Fermi momentum. Therefore, the pair potential of

s-wave  _{R} _{0}  d-wave  _{Z}

Figure 25: The s–wave superconducting and the d–wave superconducting double infinite proximity contact system.

the d–wave superconductor can be written in the form

Δ(p_{F} , z) = ∑_{m=−2}^{m=2} Δ_{m}(z)Y_{2m}(θ, ϕ).  (218)

The d–wave superconducting pair potential in general can not be chosen real. The evolution operator U_{α}, therefore, does not satisfy the relation U_{α} = ρ_{1}U_{−α}ρ_{1}. Instead, it can be found that it satisfies for the Matsubara frequency ε = iω_{n}

U_{α}(ω_{n}) = ρ_{2}U_{−α}(ω_{n})^{∗}ρ_{2},  (219)
U_{α}(ω_{n}) = ρ_{1}U_{α}(−ω_{n})^{∗}ρ_{1}.  (220)

Therefore, the quasi–classical Matsubara Green’s function satisfies

ĝ_{++}(ω_{n}) = −^{T} ĝ_{−−}(ω_{n})^{∗},  (221)
ĝ_{++}(ω_{n}) = −ρ_{1}ĝ_{++}(−ω_{n})^{∗}ρ_{1},  (222)

even when the pair potential is not real.

Setting R_{L} = 0 in the quasi–classical Green’s function of the finite triple layer system and taking the limits of L_{L}, L_{R} to infinity, the quasi–classical Green’s function for the double infinite system can be obtained. Here, we rewrite R_{0} to R. Noting the above relations valid for a complex pair potential, one can obtain the quasi–classical Green’s functions, i.e.,

ĝ_{++}^{L}(z) = i\frac{2A_{L}(z) − trA_{L}(z)}{trA_{L}(z)},  (223)
ĝ_{++}^{R}(z) = i\frac{2B_{R}(z) − trB_{R}(z)}{trB_{R}(z)},  (224)
A_{L}(z) = Ǔ_{+}^{L}(z, 0) [V_{−}^{R}(0) + Rρ_{2}^{T}V_{−}^{R}(0)ρ_{2}] ρ_{2}^{T}V_{+}^{L}(z)ρ_{2},  (225)
B_{R}(z) = V_{−}^{R}(z) [ρ_{2}^{T}V_{−}^{L}(0)ρ_{2} + RV_{+}^{L}(0)] ρ_{2}^{T}Ǔ_{+}^{R}(z, 0)ρ_{2},  (226)
V_{+}^{L}(z) = ϕ̌_{+}^{L}(z)ϕ̌_{+}^{L}(0)^{†},  (227)
V_{−}^{R}(z) = ϕ̌_{−}^{R}(z)ϕ̌_{−}^{R}(0)^{†},  (228)

where

Ǔ_{+}^{L}(z, 0) = U_{+}^{L}(z, 0)e^{E_{n}^{L}z/v_{F z}^{L}},  (229)
Ǔ_{+}^{R}(z, 0) = U_{+}^{R}(z, 0)e^{−E_{n}^{R}z/v_{F z}^{R}}  (230)

and ϕ̌_{±}(z) is a solution of the equation

∂_{z}ϕ̌_{±}(z) = \frac{1}{v_{F z}} ( \begin{matrix}−ω_{n} ∓ E_{n} & iΔ(z) \\ −iΔ(z)^{∗} & ω_{n} ∓ E_{n} \end{matrix} ) ϕ̌_{±}(z),  (231)
E_{n} = \sqrt{ω_{n}^{2} + |Δ(z)|^{2}},

with the boundary condition

∂_{z}ϕ̌_{±}(z) ⟶^{|z|→∞} 0 .

We consider a simple case where the superconducting critical temperature in both the sides are same. Anderson and Morel[37] showed that in an isotropic system the most stable d–wave superconductor pair potential is

Δ(p_{F} , z) = Δ^{AM}(z)\sqrt{\frac{4π}{5}} (Y_{2,0} + i\frac{1}{\sqrt{2}}(Y_{2,2} + Y_{2,−2})).  (232)

We consider the d–wave superconductor with the pair potential defined by Eq.(232). Figure.26 shows a profile of the pair potential near the interface. The solid line is Δ^{S}(z) and Δ^{AM}(z) in the left side and the right side, respectively. The dashed line is Δ^{S} and Δ^{AM} in the case when the other side metal is a normal metal without the pairing interaction. The pair potential is suppressed around the interface when R < 1 due to the proximity effect. When the other side is rather a superconducting metal than a normal metal, the suppression of the pair potential is small. On the other hand, we considered another d–wave superconductor defined by

Δ(p_{F} , z) = Δ_{xy} sin^{2} θ cos 2ϕ  (233)
= Δ_{xy} \frac{\sqrt{96π}}{6\sqrt{5}} (Y_{2,2} + Y_{2,−2}).

In Fig.27, the pair potential is suppressed near the interface, but the suppression in the case where the both side metals are superconductor is almost equal to that in the case where the other side metal is a normal metal. The d–wave superconductor defined by Eq.(232) and defined by Eq.(233) are different in that the pair potential of Eq.232 contains s–wave component. The s–wave pair amplitude in these d–wave superconductor, F^{S}, can be obtained by

F^{S} = \frac{T N(0)}{4} ∑_{ω_{n}} ∑_{α=±} ∫_{0}^{2π} dϕ^{α} ∫_{0}^{π/2} sin θ^{α}dθ^{α}P_{0}(p_{F} · p_{F}^{α})ĝ_{αα}(iω_{n}, Δ^{d})⌋_{1,2},  (234)

−4  0  4 0 0.4 0.8  Z/ξ Δ(Z)/πT_{C}^{S}

Figure 26: The s–wave superconductor vs the d–wave superconductor defined by Eq.(232). The dashed line is Δ^{S} and Δ^{AM} in the case when the other side metal is a normal metal without the pairing interaction. T = 0.5T_{C} and R = 0.0.

−2  0  2 0 0.4 0.8  Z/ξ Δ(Z)/πT_{C}^{S}

Figure 27: The s–wave superconductor vs the d–wave superconductor defined by Eq.(233). The dashed line is Δ^{S} and Δ_{xy} in the case when the other side metal is a normal metal without the pairing interaction. T = 0.5T_{C} and R = 0.0. where P_{0} = 1. In the d–wave superconductor defined by Eq.(233) the s–wave pair amplitude, F^{S} vanishes after the polar angle integral. On the other hand, the superconductor defined by Eq.(232) has finite s–wave amplitudes. If the d-wave superconductor has the s–wave pair amplitude, the d–wave pair potential near the interface is more enhanced than in the case where the neighbor metal is a normal metal. Due to this enhancement of the d–wave superconductor, the neighbor s–wave superconductor near the interface is also enhanced.

We next consider the case when a supercurrent flows in this proximity contact system. For simplicity, we only treat the d–wave superconductor defined by

Δ^{d}(p_{F} , z) = Δ_{20}(z)(3 cos^{2} θ − 1) ∝ Y_{2,0},  (235)

which is obtained by a linear combination of Eq.(232) with its complex conjugate. When the supercurrent flows, the pair potential has a spatially varying phase factor.

Δ^{d}(z) = |Δ_{20}(z)|e^{iφ_{z}^{d}}(3 cos^{2} θ − 1).  (236)
Δ^{S}(z) = |Δ^{S}(z)|e^{iφ_{z}^{S}}.  (237)

As is well known, when the supercurrent flows, a spatial derivative of the phase is a finite constant value in the bulk region.

\frac{d}{dz} φ_{z} → q  in the bulk region .  (238)

From Eq.(116), the electric current along the z direction is written in terms of the quasi– classical Matsubara Green’s function as

J_{z}(z) = −ev_{F} \frac{T N(0)}{2} ∫_{0}^{2π} dϕ ∫_{0}^{π/2} dθ sin θ cos θ ∑_{ω_{n}} [ĝ_{++}(ω_{n}, z) − ĝ_{−−}(ω_{n}, z)]_{11} . (239)

Noting the conservation of the electric current and setting the spatial derivative q of the phase at sufficiently large z to an appropriate value, we numerically calculated a self– consistent pair potential containing a spatially varying phase factor.

When the supercurrent flows, the quasi–classical Green’s functions can be written as

ĝ_{++}^{L}(z) = i\frac{2A_{L}(z) − trA_{L}(z)}{trA_{L}(z)},  (240)
ĝ_{++}^{R}(z) = i\frac{2B_{R}(z) − trB_{R}(z)}{trB_{R}(z)},  (241)
A_{L}(z) = Ǔ_{+}^{L}(z, 0) [V_{−}^{R}(0) + Rρ_{2}^{T} V_{−}^{R}(0)ρ_{2}] ρ_{2}^{T} V_{+}^{L}(z)ρ_{2},  (242)
B_{R}(z) = V_{−}^{R}(z) [ρ_{2}^{T} V_{−}^{L}(0)ρ_{2} + RV_{+}^{L}(0)] ρ_{2}^{T}Ǔ_{+}^{R}(z, 0)ρ_{2},  (243)
V_{+}^{L}(z) = φ̂_{z}^{L}ϕ̃_{+}^{L}(z)(φ̂_{0}^{L}ϕ̃_{+}^{L}(0))^{†},  (244)
V_{−}^{R}(z) = φ̂_{z}^{R}ϕ̃_{−}^{R}(z)(φ̂_{0}^{R}ϕ̃_{−}^{R}(0))^{†},  (245)
φ̂ = ( \begin{matrix}e^{−iφ_{z}/2} & 0 \\ 0 & e^{iφ_{z}/2} \end{matrix} ) ,  (246)

where

Ǔ_{+}^{L}(z, 0) = U_{+}^{L}(z, 0)e^{Ẽ_{n}^{L}z/v_{F z}^{L}},  (247)
Ǔ_{+}^{R}(z, 0) = U_{+}^{R}(z, 0)e^{−Ẽ_{n}^{R}z/v_{F z}^{R}}  (248)

and ϕ̃_{±}(z) is a solution of the equation

∂_{z}ϕ̃_{±}(z) = \frac{1}{v_{F z}} ( \begin{matrix}−ω̃_{n} ∓ Ẽ_{n} & i|Δ_{z}| \\ −i|Δ_{z}| & ω̃_{n} ∓ Ẽ_{n} \end{matrix} ) ϕ̃_{±}(z),  (249)
ω̃_{n} = ω_{n} + i\frac{v_{F z}∂_{z}φ_{z}}{2} ,
Ẽ_{n} = \sqrt{ω̃_{n}^{2} + |Δ(z)|^{2}},

with the boundary condition

∂_{z}ϕ̃_{±}(z) ⟶^{|z|→∞} 0.

−5  0  5 −1 01  Z/ξ |Δ|/πT_{C} and φ_{Z}

Figure 28: The pair potential when the supercurrent flows. The solid line is an absolute value of the pair potential and the dashed line is the phase of the pair potential. T = 0.1T_{C}, R = 0.0 and q_{S}ξ = 0.1.

−4  −2  0  2  4 −1 01  Z/ξ |Δ|/πT_{C} and φ_{Z}

Figure 29: The pair potential when the supercurrent flows in the s–wave vs s–wave Josephson junction. T = 0.8T_{C}, R = 0.5 and q_{S}ξ = 0.3.

−5  _{0}  _{5} Z/ξ J(Z) (a.u.)

Figure 30: The position dependence of the supercurrent numerically obtained by use of the self–consistent pair potential. T = 0.2T_{C}, R = 0.0 and q_{S}ξ = 0.1.

The gap equation can be written in terms of the off–diagonal terms of ĝ_{++}.

Δ(p_{F} , z) = \frac{πT ∑_{ω_{n}>0}^{ω_{C}} ∫_{0}^{π/2} sin θ′dθ′ṽ_{S(d)} [ĝ_{++}(iω_{n}, z)⌋_{1,2} − ĝ_{++}(iω_{n}, z)^{∗}⌋_{2,1}]}{log \frac{T}{T_{C}} + ∑_{n=0}^{ω_{C}/2πT} \frac{1}{n + 1/2}}, (250)
ṽ_{S} = 1 for s–wave,
ṽ_{d} = \frac{5}{4}(3 cos^{2} θ − 1)(3 cos^{2} θ′ − 1) for d–wave defined by Eq.(235). (251)

We obtained the self–consistent pair potential numerically by use of the Green’s function. Typical result is shown in Fig.28. We have expressed the magnitude of the supercurrent J by the spatial derivative of the phase q_{bulk}^{S}ξ at the bulk region of the s–wave superconductor. The obtained self–consistent solution shows an unexpected behavior. The phase gradient changes its sign near the interface. This result is quite different from the usual s–wave vs s–wave Josephson junction. In the s–wave vs s–wave superconductor, the spatial derivative of the phase q_{z} was larger than the bulk q_{bulk}. This can be understood from the usual G–L equation. The electric current is written as

J ∝ |Δ(z)|^{2} q_{z}.

Since |Δ(z)| near the interface are suppressed, q_{z} near the interface is larger than the bulk q_{bulk} as a result of current conservation. Numerical calculations using the quasi–classical Green’s function shown in Fig.29 for the s–wave vs s–wave junction show the same results as the G–L arguments. To check the conservation of the electric current in the s–wave vs d–wave junction, we calculated the current by use of the self–consistently solved pair potential. In Fig.30, the electric current is found to be conserved everywhere. The junction between the different pairing superconductors cannot be expressed in a form of the usual Josephson junction. This phase anomaly near the interface is remarkable at lower temperatures.( Fig.31. ) Also, under a fixed current, the phase difference at the interface does not monotonically varies as is shown in Fig.32, when the interfacial reflection coefficient R increases.

The physical origin of the above results are as yet not clarified. We try to discuss whether the pair potential profile obtained above can be obtained from the Ginzburg–Landau ( G–L ) equation. By expanding the gap equation with respect to the pair potential, we discuss the system in the G–L region. Higashitani[42] studied the Josephson junction of s–wave vs s–wave superconducting junction in the G–L expansion of the quasi–classical Green’s function. We follow the prescription proposed by Higashitani. As is well known, in the G–L region 1 − T/T_{C} ≪ 1, the zeroth, the first and the second derivatives of the pair potential are of the order of

Δ ∼ O(t),
Δ′ ∼ O(t^{2}),

Δ′′ ∼ O(t^{3}),  (252)
where t ≡ (1 − T/T_{C})^{\frac{1}{2}}

where a prime of Δ′ means the spatial derivative of Δ. Since t ≪ 1, we retain the terms up to the order of O(t^{3}) in the G–L expansion of the gap equation.

−2  0  2 0 0.4 0.8  T =0.1T_{C} T =0.2T_{C} T =0.5T_{C} Z/ξ |Δ(Z)|/πT_{C}

−2  0  2 −0.8 −0.4 0 0.4 0.8  T =0.1T_{C} T =0.2T_{C} T =0.5T_{C} Z/ξ φ_{Z}

Figure 31: The pair potential in the presence of the supercurrent at various temperatures . T = 0.1, 0.2, 0.5T_{C}, R = 0.0 and q_{S}ξ = 0.1.

−5  0  5 0 0.2 0.4 0.6  R=0.6 R=0.2 R=0.0 Z/ξ |Δ(Z)|/πT_{C}

−5  0  5 −0.5 0 0.5  R = 0.0 R = 0.2 R = 0.6 Z/ξ φ_{Z} − q · Z

Figure 32: Pair potential profile of the s–wave vs d–wave junction for various interfacial reflection coefficients under a fixed current. The phase difference of the pair potential at the interface does not monotonically varies as a function of the interfacial reflection coefficient. T = 0.1T_{C} and q_{S}ξ = 0.1.

We rewrite the Green’s function in terms of ϕ_{±} given in the appendix. Since ĝ_{++}^{R} can be shown to have the same form as ĝ_{−−}^{L} except for the superscripts R and L and for the sign of v_{F z}, we only treat ĝ_{++}^{R} . Before expanding the Green’s function ĝ_{++}^{R}, we write it in terms of ϕ_{±}

ĝ_{++}(z) = i (2A_{R}(z) − 1) ,  (253)
A_{R}(z) = φ̂_{z}^{†} \frac{ϕ̃_{−}^{R}(z)^{T}ϕ̃_{+}^{R}(z)ρ_{2}}{^{T}ϕ̃_{+}^{R}ρ_{2}ϕ̃_{−}^{R}} φ̂_{z}
+φ̂_{z}^{†} \frac{[^{T}ϕ̃_{−}^{R^{∗}}(0)φ̂_{0}^{R}R̂_{L}φ̂_{0}^{†}ϕ̃_{+}^{R}(0)]ϕ̃_{−}^{R}(z)^{T}ϕ̃_{−}^{R}(z)ρ_{2}}{[^{T}ϕ̃_{−}^{R^{∗}}(0)φ̂_{0}^{R}R̂_{L}φ̂_{0}^{†}ϕ̃_{−}^{R}(0)]^{T}ϕ̃_{−}^{R}ρ_{2}ϕ̃_{+}^{R}} φ̂_{z}e^{−2Ẽ_{n}z},  (254)
R̂_{L} = φ̂_{0}^{L^{†}} [ρ_{2}ϕ̃_{+}^{L^{∗}}(0)^{T}ϕ̃_{+}^{L}(0)ρ_{2} + Rϕ̃_{+}^{L}(0)^{T}ϕ̃_{+}^{L^{∗}}(0)] φ̂_{0}^{L},  (255)

where we have used the properties of ϕ̃ discussed in the appendix. It is useful to define ϕ̃_{±} as follows :

ϕ̃_{−}(z) = F_{−}(z) ( \begin{matrix}1 \\ iD_{−}(z) \end{matrix} ) ,  (256)
ϕ̃_{+}(z) = F_{+}(z)^{∗} ( \begin{matrix}iD_{+}(z)^{∗} \\ 1 \end{matrix} ) ,  (257)

where D_{±} and F_{±} satisfy

∂_{z}D_{−} = \frac{1}{v_{F z}} (2ω̃_{n}D_{−} − |Δ(z)|(1 − D_{−}^{2})),  (258)
∂_{z}F_{−} = −\frac{1}{v_{F z}} (ω̃_{n} − Ẽ_{n} + |Δ(z)|D_{−})F_{−},  (259)
∂_{z}D_{+}^{∗} = −\frac{1}{v_{F z}} (2ω̃_{n}D_{+}^{∗} − |Δ(z)|(1 − D_{+}^{∗^{2}})),  (260)
∂_{z}F_{+}^{∗} = \frac{1}{v_{F z}} (ω̃_{n} − Ẽ_{n} + |Δ(z)|D_{+}^{∗})F_{+}^{∗}  (261)

with the boundary condition

D_{−}|_{bulk} = D_{+}^{∗}|_{bulk} = \frac{|Δ_{bulk}|}{ω̃_{n}^{bulk} + Ẽ_{n}} ω̃_{n}^{bulk} = ω + i\frac{v_{F z}q}{2} .

Using the equation(257), one can write A_{R}(z) in terms of D and F ,

A_{R}(z) = (φ̂_{z}^{R})^{†} [A_{R}^{1}(z) + A_{R}^{2}(z) (\frac{F_{−}^{R}(z)}{F_{−}^{R}(0)})^{2} e^{−2Ẽ_{n}z}] φ̂_{z}^{R},  (262)
A_{R}^{1}(z) = \frac{( \begin{matrix}1 & −iD_{+}^{R}(z)^{∗} \\ iD_{−}^{R}(z) & D_{−}^{R}(z)D_{+}^{R}(z)^{∗} \end{matrix} )}{1 + D_{−}^{R}(z)D_{+}^{R}(z)^{∗}},  (263)
A_{R}^{2}(z) = \frac{( \begin{matrix}D_{−}^{R}(z) & i \\ iD_{−}^{R}(z)^{2} & −D_{−}^{R}(z) \end{matrix} )}{1 + D_{−}^{R}(0)D_{+}^{R}(0)^{∗}}K_{0}^{L},  (264)
K_{0}^{L} = [(D_{+}^{R}(0)^{∗} − D_{+}^{L}(0)^{∗}e^{iφ_{LR}})(1 + D_{−}^{R}(0)^{∗}D_{+}^{L}(0)e^{−iφ_{LR}})
−R(D_{−}^{R}(0)^{∗} − D_{+}^{L}(0)^{∗}e^{iφ_{LR}})(1 + D_{+}^{R}(0)^{∗}D_{+}^{L}(0)e^{−iφ_{LR}})]
/ [|1 + D_{−}^{R}(0)^{∗}D_{+}^{L}(0)e^{−iφ_{LR}}|^{2} + R|1 + D_{+}^{R}(0)^{∗}D_{+}^{L}(0)e^{−iφ_{LR}}|^{2}] , (265)

where φ_{LR} = φ_{L} − φ_{R} . Following Higashitani, from the differential equation(258) one can expand D as follows :

D_{−}(z) = \frac{|Δ(z)|}{2ω̃_{n}} + \frac{v_{F z}∂_{z}|Δ(z)|}{(2ω̃_{n})^{2}} + \frac{v_{F z}^{2}∂_{z}^{2}|Δ(z)|}{(2ω̃_{n})^{3}} − i\frac{|Δ(z)|v_{F z}^{2}∂_{z}q_{z}}{(2ω̃_{n})^{3}} − (\frac{|Δ(z)|}{2ω̃_{n}})^{3} + ... (266)

D_{+} can be obtained by changing the sign of v_{F z} in the expansion of D_{−}.

The Green’s function consists of two terms. ( See Eqs.(253)∼(255). ) One of them spatially varies on the scale of the order of the coherence length ξ_{T} (t) ∼ O(t^{−1}) since it consists of ϕ_{±} only. The other term varies on the scale of the order of ξ_{0} = v_{F} /πT_{C}, since it has an exponential damping factor. In the G–L region in which ξ_{T} is sufficiently larger than ξ_{0}, the second term is finite only in a narrow region near the interface. The gap equation can be written in such a form as

Δ(z) = ⟨⟨g_{1}(z) + g_{2}(z)e^{Ω_{z}}⟩⟩  (267)

where g_{1} and g_{2} are functions varying on the scale of the order ξ_{T} which are related to A_{R}^{1} and A_{R}^{2}, respectively. The short range behavior is described by the exponential factor e^{Ω_{z}} damping within the range of ξ_{0}. In the above, ⟨⟨...⟩⟩ means the average over the polar angles and the sum over ω_{n}. In Eq.(267), the first term g_{1} describes the long range profile of the pair potential and the second term g_{2}e^{Ω_{z}} is considered to give an information on the boundary condition. Following Higashitani’s method, we consider the equations :

Δ(z) = ⟨⟨g_{1}(z)⟩⟩,  (268)
0 = ∫_{0}^{∞} dz⟨⟨g_{2}(z)e^{Ω_{z}}⟩⟩.  (269)

Equation(268) gives the gap equation and Eq.(269) gives the boundary condition. Higashitani obtained the same boundary condition from Eq.(269) as that which de Gennes obtained in the S–S′ contact system by different approach.

Substituting Eq.(266) into the Green’s function, we obtain the G–L equation after some manipulation :

−Δ^{r}(z) = ξ_{T}^{2} [⟨Θ_{R} cos^{2} θ⟩_{R}(Δ^{r}(z)q^{r}(z)^{2} − Δ^{r}(z)′′) + ⟨Θ_{R}^{3}⟩_{R}(Δ^{r}(z))^{3}/v_{F}^{2} ] , (270)
−Δ^{l}(z) = ξ_{T}^{2} [⟨Θ_{L} cos^{2} θ⟩_{L}(Δ^{l}(z)q^{l}(z)^{2} − Δ^{l}(z)′′) + ⟨Θ_{L}^{3}⟩_{L}(Δ^{l}(z))^{3}/v_{F}^{2} ] ,  (271)
0 = 2Δ^{r}(z)′q^{r}(z) + Δ^{r}(z)q^{r}(z)′,  (272)
0 = 2Δ^{l}(z)′q^{l}(z) + Δ^{l}(z)q^{l}(z)′,  (273)

where we have defined

|Δ^{R(L)}(z)| = Δ^{r(l)}(z)Θ_{R(L)},  (274)
Θ_{R(L)} = \{ \begin{matrix}1, & for S region \\ 3 cos^{2} θ − 1, & for d region \end{matrix}  (275)
ξ_{T}^{2} = \frac{7ζ(3)}{16 log T/T_{C}} (\frac{v_{F}}{πT} )^{2} ,  (276)
⟨...⟩ ≡ ∫_{0}^{π/2} dθ sin θ × v_{S(d)}.... ,  (277)
v_{S} = 1 and v_{d} = \frac{5}{4}Θ_{d} .  (278)

The boundary condition are

0 = R̃⟨Θ_{R} cos θ⟩_{R}Δ^{r}(0) − a^{R}(Δ^{r}(0)′ + iΔ^{r}(0)q^{r}(0))⟨Θ_{R} cos^{2} θ⟩_{R}
−e^{iφ_{LR}} [R̃⟨Θ_{L} cos θ⟩_{R}Δ^{l}(0) − a^{L}R̃(Δ^{l}(0)′ + iΔ^{l}(0)q^{l}(0))⟨Θ_{L} cos^{2} θ⟩_{R}] , (279)
0 = R̃⟨Θ_{R} cos θ⟩_{L}Δ^{r}(0) + a^{R}R̃(Δ^{r}(0)′ + iΔ^{r}(0)q^{r}(0))⟨Θ_{R} cos^{2} θ⟩_{L}
−e^{iφ_{LR}} [R̃⟨Θ_{L} cos θ⟩_{L}Δ^{l}(0) + a^{L}(Δ^{l}(0)′ + iΔ^{l}(0)q^{l}(0))⟨Θ_{L} cos^{2} θ⟩_{L}] , (280)
a^{R(L)} = \frac{7ζ(3)}{2π^{2}} \frac{v_{F}^{R(L)}}{πT}  (281)

and

R̃ = \frac{1 − R}{1 + R}.

We have taken account of terms up to order of O(t^{3}) for the G–L equation. On the other hand, the boundary condition includes terms up to order of O(t^{2}), in accordance with the boundary condition obtained by de Gennes or by Higashitani. The obtained boundary conditions is linear. They can be rewritten as

( \begin{matrix}R̃D_{r} \\ a^{R}D′_{r} \end{matrix} ) = ( \begin{matrix}N_{11} & N_{12} \\ N_{21} & N_{22} \end{matrix} ) ( \begin{matrix}R̃D_{l} \\ a^{L}D′_{l} \end{matrix} ) ,  (282)
D_{r} = α_{r}Δ̃_{r} = [⟨Θ_{R} cos θ⟩_{L}⟨Θ_{R} cos^{2} θ⟩_{R} + R̃⟨Θ_{R} cos θ⟩_{R}⟨Θ_{R} cos^{2} θ⟩_{L}]^{\frac{1}{2}} Δ̃_{r}, (283)
D_{l} = α_{l}Δ̃_{l} = [⟨Θ_{L} cos θ⟩_{R}⟨Θ_{L} cos^{2} θ⟩_{L} + R̃⟨Θ_{L} cos θ⟩_{L}⟨Θ_{L} cos^{2} θ⟩_{R}]^{\frac{1}{2}} Δ̃_{l}, (284)
D′_{r} = α_{r}Δ̃′_{r} , D′_{l} = α_{l}Δ̃′_{l} , (285)
N_{11} = \frac{1}{α_{r}α_{l}} [⟨Θ_{L} cos θ⟩_{L}⟨Θ_{R} cos^{2} θ⟩_{R} + R̃⟨Θ_{L} cos θ⟩_{R}⟨Θ_{R} cos^{2} θ⟩_{L}] ,  (286)
N_{12} = \frac{1}{α_{r}α_{l}} [⟨Θ_{L} cos^{2} θ⟩_{L}⟨Θ_{R} cos^{2} θ⟩_{R} − R̃^{2}⟨Θ_{L} cos^{2} θ⟩_{R}⟨Θ_{R} cos^{2} θ⟩_{L}] ,  (287)
N_{21} = \frac{1}{α_{r}α_{l}} [⟨Θ_{L} cos θ⟩_{L}⟨Θ_{R} cos θ⟩_{R} − ⟨Θ_{L} cos θ⟩_{R}⟨Θ_{R} cos θ⟩_{L}] ,  (288)
N_{22} = \frac{1}{α_{r}α_{l}} [⟨Θ_{L} cos^{2} θ⟩_{L}⟨Θ_{R} cos θ⟩_{R} + R̃⟨Θ_{L} cos^{2} θ⟩_{R}⟨Θ_{R} cos θ⟩_{L}] ,  (289)

where N satisfies det N = 1 and

Δ̃_{r} ≡ \frac{Δ^{R}(z)}{Θ_{R}} |_{z=0} ,  (290)
Δ̃_{l} ≡ \frac{Δ^{L}(z)}{Θ_{L}} |_{z=0} ,  (291)
Δ̃′_{r} ≡ ∂_{z} \frac{Δ^{R}(z)}{Θ_{R}} |_{z=0} ,  (292)
Δ̃′_{l} ≡ ∂_{z} \frac{Δ^{L}(z)}{Θ_{L}} |_{z=0} .  (293)

We obtained the linear boundary condition similar to de Gennes’ one. It is, however, found that when the superconductors have different pairing symmetries, this discussion is not consistent. From Eq.(239), the current in the G–L temperature range in the bulk region can be written as

J_{z}(z) = ev_{F} \frac{7ζ(3)}{(2πT)^{2}} N(0)⟨Θ⟩_{j}Δ(z)^{2}v_{F} q(z),  (294)

where

⟨...⟩_{j} = ∫_{0}^{π/2} dθ sin θ cos θ ... .

From the requirement that the current should be the same in both the bulk region, we obtain

⟨Θ_{R}^{2}⟩_{j}v_{F} Δ^{r}(0)^{2}q^{r}(0) = ⟨Θ_{L}^{2}⟩_{j}v_{F} Δ^{l}(0)^{2}q^{l}(0).  (295)

On the other hand, using the linear boundary condition Eq.(282), we obtain the condition

α_{r}^{2}v_{F} Δ^{r}(0)^{2}q_{r}(0) = α_{l}^{2}v_{F} Δ^{l}(0)^{2}q_{l}(0).  (296)

When the pairing symmetry of the both side superconductors is different, the current conservation condition Eq.(295) is not equivalent to the condition Eq.(296) derived from the linear boundary condition. When the pairing of both the sides is of the same type, as the S–S′ system discussed by Higashitani, there occurs no such problem. When the supercurrent flows across the junction between different symmetry superconductors, the above discussed G–L argument seems not to be consistent. We point out some problems at issue.

1. We have derived the boundary condition Eq.(269) by dividing Eq.(267) into two parts. It is not clear whether the way employed to obtain the boundary condition Eq.(269) is unique in the G–L region, although we never find out a more reasonable way than that used above which was successful to obtain the boundary condition in the s-wave vs s-wave system.

2. The source of the anomaly should be the exponential damping term g_{2}e^{Ω_{z}} in Eq.(267).

exp(Ω_{z}) ∼ exp z/ξ_{0}.

The spatial derivative of this term is of order of ξ_{0}^{−1}O(t) although a range of the anomaly is negligible in the G–L region ( ξ_{0} ≪ ξ_{T} ). The spatial derivative of the damping term may not be negligible in the present case and we may have to take into account the higher order terms in the G–L expansion of the damping term.

These problems are to be examined in the future study.

Appendix B. The properties of the evolution operator ( in the presence of the supercurrent flows. )

In this appendix, we consider the properties of the evolution operator when the pair potential is not real. When the supercurrent flows, the phase of the pair potential spatially varies

Δ(z) = |Δ(z)|e^{iφ_{z}}  (297)

→ |Δ_{bulk}|e^{i(q_{bulk}z+α)} at the bulk region  (298)

where q_{bulk} = ∂_{z}φ_{z}⌋_{bulk} and α is a constant. We introduce a modified evolution operator defined as

Ũ_{α}(z, z′) = φ̂_{z} U_{α}(z, z′) φ̂_{z}^{†}  (299)
φ̂_{z} = ( \begin{matrix}e^{−iφ_{z}/2} & 0 \\ 0 & e^{iφ_{z}/2} \end{matrix} )  (300)
Δ(z) = |Δ(z)|e^{iφ_{z}}.  (301)

Then, the evolution operator satisfies the differential equation

∂_{z}Ũ_{α}(z, z′) = αA_{α}Ũ_{α}(z, z′) (α = ±1)  (302)
A_{α} = \frac{i}{v_{F z}} ( \begin{matrix}ε̃_{α} & |Δ(z)| \\ −|Δ(z)| & −ε̃_{α} \end{matrix} )  (303)
ε̃_{α} = ε − α\frac{v_{F z}∂_{z}φ_{z}}{2}  (304)

The matrix A has also an important property

ρ_{2}^{T}A_{α}ρ_{2} = −A_{α}.  (305)

It is to be noted that when the supercurrent flows, the evolution operator U does not satisfy

U_{α} = ρ_{1}U_{−α}ρ_{1}.

Except for this, the above defined evolution operator Ũ has the same properties as U discussed in the appendix of Sec.2 . Thus, Ũ can be decomposed as follows :

Ũ_{α}(z, z′) = Λ_{+}^{α}(z, z′)e^{−iκ^{α}(z−z′)} + Λ_{−}^{α}(z, z′)e^{iκ^{α}(z−z′)},  (306)
Λ_{+}^{α}(z, z′) = \frac{−1}{W_{α}} ϕ̃_{+}^{α}(z)^{T}ϕ̃_{−}^{α}(z′)ρ_{2},  (307)
Λ_{−}^{α}(z, z′) = \frac{1}{W_{α}} ϕ̃_{−}^{α}(z)^{T}ϕ̃_{+}^{α}(z′)ρ_{2}.  (308)
W_{α} = ^{T}ϕ̃_{+}^{α}ρ_{2}ϕ̃_{−}^{α} = −^{T}ϕ̃_{−}^{α}ρ_{2}ϕ̃_{+}^{α} = const. ≠ 0  (309)

and ϕ̃_{±}^{α} are the solution of the equation

∂_{z}ϕ̃_{±}^{α} = \frac{αi}{v_{F z}} ( \begin{matrix}ε_{α} ± Ω_{α} & |Δ(z)| \\ −|Δ(z)| & −ε_{α} ± Ω_{α} \end{matrix} ) ϕ̃_{±}^{α},  (310)
κ^{α} = \sqrt{ε_{α}^{2} − |Δ(z)|^{2}} = Ω_{α}/v_{F z}  (311)

with the boundary condition

∂_{z}ϕ̃_{±}^{α} ⟶ 0 at sufficient large z,

that is,

ϕ̃_{±}^{α} ⟶ const. × ( \begin{matrix}ε_{α} ∓ Ω_{α} \\ −|Δ_{bulk}| \end{matrix} ) .  (312)

In general, one can not obtain the relation between U_{+} and U_{−} if the pair potential is not real. In the case of the Matsubara frequency ε = iω_{n}, however, one can obtain the relation of U_{±}.

Evolution operator for the Matsubara frequency ε = iω_{n}

For the Matsubara frequency ε = iω_{n}, the evolution operator satisfies

U_{α}(ω_{n}) = ρ_{2}U_{−α}(ω_{n})^{∗}ρ_{2}  (313)

and

U_{α}(ω_{n}) = ρ_{1}U_{α}(−ω_{n})^{∗}ρ_{1},  (314)

when Δ(p_{F}^{+}) = Δ(p_{F}^{−}) is satisfied. Therefore, we have only to treat U_{+} for the Matsubara frequency. We write down some useful properties of the evolution operator for the Matsubara frequency E = iω_{n} :

U_{+}(z, z′) = φ̂_{z}^{†}Ũ_{+}(z, z′)φ̂_{z′},  (315)

U_{−}(z, z′) = φ̂_{z}^{†}ρ_{2}Ũ_{+}(z, z′)^{∗}ρ_{2}φ̂_{z′},  (316)
Ũ_{+}(z, z′) = Λ_{+}^{ω_{n}}(z, z′)e^{Ẽ_{n}(z−z′)} + Λ_{−}^{ω_{n}}(z, z′)e^{−Ẽ_{n}(z−z′)},  (317)
Λ_{±}^{ω_{n}}(z, z′) = \frac{∓1}{W}ϕ̃_{±}^{ω_{n}}(z)^{T}ϕ̃_{∓}^{ω_{n}}(z′)ρ_{2},  (318)
W = ^{T}ϕ̃_{+}^{ω_{n}}ρ_{2}ϕ̃_{−}^{ω_{n}} = −^{T}ϕ̃_{−}^{ω_{n}}ρ_{2}ϕ̃_{+}^{ω_{n}} = const. ≠ 0,  (319)

and ϕ̃_{±}^{ω_{n}} satisfies an equation

∂_{z}ϕ̃_{±}^{ω_{n}} = \frac{1}{v_{F z}} ( \begin{matrix}−ω̃_{n} ∓ Ẽ_{n} & i|Δ(z)| \\ −i|Δ(z)| & ω̃_{n} ∓ Ẽ_{n} \end{matrix} ) ϕ̃_{±}^{ω_{n}},  (320)
ω̃_{n} = ω_{n} + i\frac{v_{F z}∂_{z}φ_{z}}{2} ,  (321)
Ẽ_{n} = \sqrt{ω̃_{n}^{2} + |Δ_{bulk}|^{2}},  (322)

with the boundary condition

ϕ̃_{±}^{ω_{n}} ⟶ const. × ( \begin{matrix}ω̃_{n} ∓ Ẽ_{n} \\ −i|Δ_{bulk}| \end{matrix} ) .  (323)

6 Summary

We have derived the quasi–classical Green’s function in the triple layer system including semi–infinite superconductor by extending the AAHN formulation for double layer system in the clean limit and have studied some superconducting proximity contact systems.

Following AAHN, we first constructed a solution of the Green’s function for the superconducting finite triple layer system in a form including the spatial evolution operator within the quasi–classical approximation. Taking the limit of both the layer sizes L_{L} and L_{R} to infinity ( keeping the center layer size L finite ), we have obtained the solution of the quasi–classical Green’s function for the semi–infinite triple layer system. The present formulation has a great advantage in computing the self–consistent pair potential. In the conventional quasi–classical Green’s function method, one has to solve the Eilenberger equation under the restriction of the normalization condition as well as of the boundary condition[23]. In numerical calculations according to that program, one needs sophisticated techniques to find converging solutions at infinities. In the present formulation, we have obtained an explicit form of the Green’s function which already satisfies the boundary condition and is written by quantities converging at infinities. This reduces the numerical efforts very much.

One of the applications of the present formulation is a study of the point contact experiment. Taking account of the reflection coefficients at the point contact and at the normal–superconducting interface, we have calculated the Andreev reflection. It is found that the Andreev reflection, which is closely related to the differential conductance, can have double peak structure as a function of the incident energy below the bulk energy gap Δ_{bulk}^{S}. This happens because of a finiteness of the normal layer and finite interfacial reflection coefficients. It was found that the differential conductance also can has double peak structure as a function of the bias voltage.

Next, we have obtained the density of states of the N region for the superconducting– normal–superconducting proximity contact system. It is found that the density of states has a structure below the bulk energy gap Δ_{bulk}^{S}, in a similar manner to the Andreev reflection in the N–N′–S system. This structure comes from the de Gennes–Saint-James bound state which is originated by the Andreev reflection at the N–S interface. Also we have calculated the total density of states. This total density of states can be detected by scanning tunneling spectroscopy. The spatial variation of the density of states, which was detected in a STM experiment, has been realized by taking account of the pairing interaction of the N region.

We have also studied the s–wave superconductor and the d–wave superconductor junction. The self-consistent calculation of the pair potential can be achieved by use of the present quasi-classical formulation. The pair potential near the interface are suppressed due to the proximity effect. When the supercurrent flows through the junction, it has been found from the numerical results, that the spatial variation of the phase φ_{z} of the pair potential in the s–wave vs d–wave junction is different from that in the s–wave vs s–wave junction near the interface. A sign of the spatial derivative of the phase near the interface is opposite to that of the bulk region. Also, even if the reflection coefficient is zero, the phase shows a jump at the interface. To interpret these unordinary behavior, we have tried to analyze this junction by use of the G–L expansion devised for the S–S junction. It was, however, found that an obtained boundary condition is not consistent. The anomalous behavior of the pair potential is remarkable at lower temperature. It indicates that such behaviors may not be reproduced by the G–L expansion. This problem is still to be examined in future study.

We have studied some systems in the clean limit and in equilibrium. All actual systems, however, are not in the clean limit and not in equilibrium. It is important to take account of the impurity effect and extend the formulation to the non–equilibrium system.

7 Acknowledgment

I would like to thank Professor Katsuhiko Nagai for various instruction, suggestion, stimulating discussion and critical reading of the manuscript. I would like to thank Doctor Hiroshi Shimahara for stimulating discussions and helpful advices. I would like to thank Professor Jun’ichiro Hara of Yamaguchi University for helpful discussions. I would like to thank Professor Mitsuo Watabe, Professor Hironobu Fujii and Professor Kozo Hoshino for various advices. I would like to thank Seiji Higashitani for various discussions. I am very grateful to my colleagues in Professor Nagai’s group for discussions, instruction and encouragement.

References

[1] J.Bardeen, L.N.Cooper, and J.R.Schrieffer , Phys.Rev.108, 1175(1957).

[2] G.Deutcher and P.G. de Gennes, in Superconductivity ed. R.D. Parks ( Marcel Dekker, Inc., New York, 1969), p.1005.

[3] A.A.Abrikosov, L.P.Gorkov, and I.E.Dzyaloshinski Methods of Quantum Field Theory in Statistical Physics

[4] 中嶋貞雄 ”超伝導入門” 培風館

[5] A.F.Andreev , Zh.Eksp.Teor.Fiz.46,1823(1964)[Sov.Phys.-JETP 19,1228(1964)]

[6] P.C.van Son, H.van Kempen and P.Wyder , Phys.Rev.Lett.59, 2226(1987)

[7] T.Nishino, M.Hatano, H.Hasegawa, T.Kure and F.Murai , Phys.Rev.B41,7274(1990)

[8] P.A.M.Benistant, H.van Kempen, and P.Wyder , Phys.Rev.Lett.51, 817(1983)

[9] H.F.C.Hoevers, P.J.M.van Bentum, L.E.C.van de Leemput , H.van Kempen , A.J.G.Schellingerhout and D.van der Marel , Physica C152, 105(1988)

[10] G.E.Blonder, M.Tinkham, and T.M.Klapwijk , Phys.Rev.B25, 4515(1982)

[11] A.Furusaki , Master thesis, University of Tokyo, 1990

[12] Chr.Bruder , Phys.Rev.B41, 4017(1990)

[13] P.C.van Son, H.van Kempen and P.Wyder , Phys.Rev.B37, 5015(1988)

[14] K.Nagai, Y.Nagato, S.Higashitani, M.Ashida and J.Hara , Physica C 185-189, 1451(1991)

[15] Y.Nagato, K.Nagai and J.Hara , J.Low Temp.Phys.93, 33(1993)

[16] A.Kastalsky, A.W.Kleinsasser, L.H.Greene, R.Bhat , F.P.Milliken and J.P.Harbison, Phys.Rev.Lett.67, 3026(1991)

[17] A.Kastalsky, L.H.Greene, J.B.Barner and R.Bhat , Phys.Rev.Lett.64, 958(1990)

[18] G.Eilenberger , Z.Phys.214, 195(1968)

[19] A.I.Larkin and Yu.N.Ovchinnikov , Zh.Eksp.Teor.Fiz.55, 2262(1968) [Soviet.Phys.- JETP28, 1200(1969)

[20] M.Ashida, S.Aoyama, J.Hara and K.Nagai , Phys.Rev.B40, 8673(1989)

[21] M.Ashida, J.Hara, and K.Nagai , Phys.Rev.B45, 828(1992)

[22] A.L.Shelankov , Fiz.Tverd.Tela26, 1615(1984) [Soviet Phys.-Solid State26, 981(1984)]

[23] A.V.Zaǐtsev , Zh.Eksp,Teor.Fiz.86, 1742(1984) [Soviet.Phys.-JETP59, 1015(1984)]

[24] A.Millis, D.Rainer, and J.A.Sauls , Phys.Rev.B38, 4504(1988)

[25] J.W.Serene and D.Rainer , Phys.Rep.101, 221(1983)

[26] J.Kurkijärvi and D.Rainer , in Modern P roblems in Condensed Matter Sciences , edited by W.P.Harperin and L.P.Pitaevskii (Elsevier,Amsterdam,1989)

[27] S.Higashitani and K.Nagai , to be published in the Proceedings of LT20 (University of Oregon,1993).

[28] P.G.deGennes and D.Saint-James , Phys.Lett.4, 151(1963)

[29] W.L.McMillan , Phys.Rev.175, 537(1968)

[30] O.Entin-Wohlman and J.Bar-Sagi , Phys.Rev.B18, 3174(1978)

[31] G.B.Arnold , Phys.Rev.B18, 1076(1978)

[32] W.J.Gallagher , Phys.Rev.B22, 1233(1980)

[33] G.Kieselmann , Phys.Rev.B35, 6762(1987)

[34] Y.Tanaka, H.Yamagami, and M.Tsukada , Solid state Commun.79, 349(1991)

[35] J.Hara, M.Ashida, and K.Nagai , Phys.Rev.B47, 11263(1993)

[36] K.Inoue and H.Takayanagi , Phys.Rev.B43, 6214(1991)

[37] P.W.Anderson and P.Morel , Phys.Rev.123, 1911(1961)

[38] R.Balian and N.R.Werthamer , Phys.Rev.131, 1553(1963)

[39] J.A.Pals, W. van Haeringen, and M.H.van Maaren , Phys.Rev.B15, 2592(1977)

[40] U.Poppe , J.Mag.Mag.Mater.52, 157(1985)

[41] B.Ashauer, G.Kieselmann, and D.Rainer , J.Low Temp.Phys.63, 349(1986)

[42] S.Higashitani , Master thesis, Yamaguchi Univ. (1991)