Nonlinear Resonance of Superconductor/Normal Metal Structures to Microwaves
aa r X i v : . [ c ond - m a t . s up r- c on ] A p r Nonlinear Resonance of Superconductor/Normal Metal Structures to Microwaves
E. Kandelaki, A. F. Volkov,
1, 2
K. B. Efetov, and V. T. Petrashov Theoretische Physik III, Ruhr-Universit¨at Bochum, D-44780 Bochum, Germany Institute of Radioengineering and Electronics of the Russian Academyof Sciences, 103907 Moscow, Russia Department of Physics, Royal Holloway, University of London,Egham, Surrey TW20 0EX, United Kingdom (Dated: October 30, 2018)We study the variation of the differential conductance G = dj/dV of a normal metal wire ina Superconductor/Normal metal heterostructure with a cross geometry under external microwaveradiation applied to the superconducting parts. Our theoretical treatment is based on the quasiclas-sical Green’s functions technique in the diffusive limit. Two limiting cases are considered: first, thelimit of a weak proximity effect and low microwave frequency, second, the limit of a short dimension(short normal wire) and small irradiation amplitude. PACS numbers: 74.45.+c, 74.50.+r, 85.25.Dq, 03.67.Lx, 85.25.Cp
I. INTRODUCTION
Superconductor/normal metal (S/N) nanostructures,where the proximity effect (PE) plays an important role,have been studied very actively during last two decades.Interesting phenomena have been discovered in thecourse of these studies. Perhaps, the most remarkableone is an oscillatory dependence of the conductance ofa normal wire attached to two superconductors whichare incorporated into a superconducting loop . Thisphenomenon was observed in the so-called “Andreevinterferometers”, i.e. in multi-terminal SNS junctions(see as well as reviews and references therein).The reason for this oscillatory behavior of the differ-ential conductance G = dj/dV is a modification of thetransport properties of the n wire due to the PE, i.e. dueto the condensate induced in the n wire. The density ofthe induced condensate is very sensitive to an appliedmagnetic field H and oscillates with increasing H .Theory was successful in explaining the exper-iments and predicting new phenomena, including there-entrance of the conductance to the normal state inmesoscopic proximity conductors and transitions tothe π -state in the voltage-biased Andreev interferometersdue to non-equilibrium effects . The non-monotonicbehavior of the conductance in SN point contactsand controllable nanostructrures has been observed inRefs. , and the change of the sign of the criticalJosephson current in multiterminal SNS junctions hasbeen found in Refs. . Many important results of thestudy of the SN mesoscopic structures are reviewed inRefs. .The so-called π -states have also been realized in equi-librium Josephson SFS junctions with a ferromagnetic(F) layer between superconductors or in SIS Joseph-son junctions of high- T c , d-wave superconductors .A number of new phenomena have been discovered in thin one-dimensional N and S wires (see also fora recent review and references therein).Mesoscopic SNS structures proved to be a promisingalternative to Superconducting Quantum InterferenceDevices (SQUIDs) for certain applications, includingmagnetic flux measurements and read-out of quantumbits (qubits) with a potential to achieve higher thanstate-of-the-art fidelity, sensitivity and read-out speed.To achieve such challenging aims extensive investigationsof high frequency properties of S/N nanostructures on ascale similar to that of SQUIDs are in order.Studies undertaken to date concerned mainly the sta-tionary properties of S/N structures. Experimental dataon S/N structures under microwave radiation appearedonly recently . As to theoretical studies, one canmention two papers where the ac impedance of aS/N structure was calculated. However, measuring thefrequency dependence of the ac conductance is not aneasy task. It is more convenient to measure a nonlineardc response (dc conductance) to a microwave radiation.Recently, a numerical calculation of the dependenceof the critical Josephson current I c in SNS junctionon the amplitude of an external ac radiation has beenperformed .In this paper, using a simple model we calculate thedc conductance of a normal ( n ) wire in an S/N structure(cross geometry) as a function of the frequency Ω andthe amplitude of the external microwave radiation. Weconsider the limiting cases of a long and a short n wireand show that the response has a resonance peak at afrequency Ω close to ε s / ¯ h , where ε s is the energy of asubgap in the n wire induced by the PE. Our theorypredicts novel resonances and can help to optimize quan-tum devices based on hybrid SNS nanostructures .We employ the quasiclassical Green’s function tech-nique in the diffusive limit. This means that we will y x N N − V + V SS + ϕ ( t ) − ϕ ( t ) n n - L x L x L y - L y FIG. 1: Structure under consideration solve the Usadel equation for the retarded (advanced)Green’s function ˆ g R ( A ) and the corresponding equationfor the Keldysh matrix function ˆ g K (section 2). First,a weak PE will be considered when the Usadel equationcan be linearized (section 3). We calculate the dc con-ductance of the n wire in this limit, assuming that thefrequency of the ac radiation is low (Ω ≪ T ). In sec-tion 4, the opposite limiting case of a short n wire willbe analyzed for arbitrary frequencies Ω. We present thefrequency dependence of the correction to the dc conduc-tance caused by ac radiation. In section 5, we discuss theobtained results. II. MODEL AND BASIC EQUATIONS
We consider an S/N structure shown in Fig. 1. It con-sists of a n wire or n film which connects two N and S reservoirs ( n and N stand for a normal metal, S means asuperconductor). The superconducting reservoirs may beconnected by a superconducting contour. The transversedimensions of the n wire are supposed to be smaller thancharacteristic dimensions of the problem, but larger thanthe Fermi wave length and the mean free path l (diffu-sive case). This implies that all quantities depend onlyon coordinates along the wire (the x − coordinate in thehorizontal direction and the y − coordinate in the verti-cal direction). The dc voltage 2 V is applied between thenormal N reservoirs, and the phase difference 2 ϕ existsbetween the superconducting reservoirs. The phase ϕ isassumed to be time-dependent ϕ ( t ) = ϕ + ϕ Ω cos(Ω t ) (1)and related to the magnetic flux Φ inside the supercon-ducting contour: ϕ ( t ) = π Φ( t ) / Φ with Φ( t ) = H ( t ) S ,where H ( t ) is an applied magnetic field and S is the areaof the superconducting contour; that is, the magneticfield contains not only a constant component, but alsoan oscillating one. For simplicity, we assume the structure to be symmet-ric both in the horizontal and vertical directions. Thisimplies, in particular, that the interface resistances R nN at x = ± L x are equal to each other (correspondingly, R nS ( L y ) = R nS ( − L y )). Our aim is to calculate the dif-ferential dc conductance G between the N reservoirs G = djdV (2)as a function of the amplitude of the ac signal ϕ Ω andthe frequency Ω.The calculations will be carried out on the basis of thewell developed quasiclassical Green’s functions technique(see the reviews ) which successfully was applied tothe study of S/N structures (see for example ).In this technique, all types of Green’s functions (the ”nor-mal” and Gor’kov’s functions as well as the retarded, ad-vanced and Keldysh functions) are matrix elements of a4 × g = (cid:18) ˆ g R ˆ g K g A (cid:19) (3)where ˆ g R ( A ) are matrices of the retarded (advanced)Green’s functions, and ˆ g K is a matrix of the Keldyshfunctions. The first matrices describe thermodynamicalproperties of the system (the density of states, super-current etc), whereas the matrix ˆ g K is used to describedissipative transport and nonequilibrium properties.The matrix ˇ g satisfies the normalization condition (ˇ g ◦ ˇ g ) ( t, t ′ ) = δ ( t − t ′ ) (4)where ” ◦ ” denotes the integral product (ˇ g ◦ ˇ g ) ( t, t ′ ) = R dt ˇ g ( t, t ) · ˇ g ( t , t ′ ) and ” · ” is the conventional matrixproduct. The Fourier transform performed as ˇ g ( ε, ε ′ ) = R dt dt ′ e iεt − iε ′ t ′ ˇ g ( t, t ′ ) yields (ˇ g ◦ ˇ g ) ( ε, ε ′ ) = 2 πδ ( ε − ε ′ )where now (ˇ g ◦ ˇ g ) ( ε, ε ′ ) = R dε π ˇ g ( ε, ε ) · ˇ g ( ε , ε ′ ).The matrix of Keldysh functions ˆ g K can be expressedin terms of the matrices ˆ g R ( A ) and a matrix of distribu-tion functions ˆ F : ˆ g K = ˆ g R ◦ ˆ F − ˆ F ◦ ˆ g A (5)where the matrix ˆ F can be assumed to be diagonal :ˆ F = ˆ τ F + + ˆ τ F − . (6)Here ˆ τ is the identity matrix and ˆ τ the third Paulimatrix. The function F − describes the charge imbalance(premultiplied with the DOS and integrated over all en-ergies it gives the local voltage), while F + characterizesthe energy distribution of quasiparticles.Due to the general relation ˆ g A ( ε, ε ′ ) = − ˆ τ ˆ g R † ( ε ′ , ε )ˆ τ (7)one can immediately calculate ˆ g A after finding thematrix ˆ g R . That means that knowing the matrices ˆ g R and ˆ F we can determine all entries of ˇ g .The Green’s functions in N and S reservoirs are as-sumed to have an equilibrium form corresponding to thevoltages ± V and phases ± ϕ ( t ). For example, the re-tarded (advanced) Green’s functions in the upper S reser-voir are ˆ g R ( A ) S ( t, t ′ ) = ˆ S ( t ) · ˆ g R ( A ) S ( t − t ′ ) · ˆ S † ( t ′ ) (8)where ˆ S ( t ) = exp[ i ˆ τ ϕ ( t ) /
2] is a unitary transformationmatrix and the Fourier transform of ˆ g R ( A ) S ( t − t ′ ) equalsˆ g R ( A ) S ( ε ) = 1 ξ R ( A ) ε (cid:18) ε ∆ − ∆ − ε (cid:19) (9)with ξ R ( A ) ε = ± p ( ε ± i − ∆ , i.e. the matrix ˆ g R ( A ) S describes the BCS superconductor in the absence ofphase. The retarded (advanced) Green’s functions in thelower S reservoir are determined in the same way withthe replacement ϕ ( t ) → − ϕ ( t ). The matrix ˆ g R ( A ) N in theright (left) N reservoirs is equal to ˆ g R ( A ) S with ∆ = 0,i.e. ˆ g R ( A ) N = ± ˆ τ .In the reservoirs the matrix ˆ F ( t, t ′ ) can be repre-sented through the equilibrium distribution F eq ( ε ) =tanh( ε/ T ) via Eq. (8)ˆ F ( t, t ′ ) = ˆ S ( t ) · F eq ( t − t ′ ) ˆ S † ( t ′ ) . (10)The phase ϕ ( t ) in the upper S reservoir is given byEq. (1), and for ϕ N ( t ) in the right N reservoir, we have: ϕ N ( t ) = 2 eV t . Therefore in the normal reservoir (atthe right) the matrix distribution function has diagonalelements ˆ F N ( ε ) , = tanh (cid:0) T ( ε ± eV ) (cid:1) , and can bewritten as ˆ F N ( ε ) = ˆ τ F N + ( ε ) + ˆ τ F N − ( ε ) F N ± ( ε ) = 12 (cid:20) tanh ε + eV T ± tanh ε − eV T (cid:21) . (11)Thus, all Green’s functions in the reservoirs are definedabove.Our task is to find the matrix ˇ g in the n wire. In theconsidered diffusive limit it obeys the equation ˇ τ · ∂ ˇ g∂t + ∂ ˇ g∂t ′ · ˇ τ + i ( eV n ( t )ˇ g − ˇ geV n ( t ′ )) − D ∇ (ˇ g ◦∇ ˇ g ) = 0 , (12)where ˇ τ is a diagonal matrix with equal elements(ˇ τ ) , = ˆ τ , V n is a local electrical potential in the n wire. We dropped the inelastic collision term suppos-ing that E T h = D/L max ≫ τ − inel , where D is the diffusioncoefficient, L max = max { L x,y } and τ inel is an inelasticscattering time. This equation is complemented by theboundary condition ˇ g ◦ ∂ x,y ˇ g | x,y = ± L x,y = ± κ N,S [ˇ g, ˇ g N,S ] ◦ (13) where κ N,S = 1 / (2 σR nN,nS ) , R nN,nS are the nN and nS interface resistances per unit area and σ is the conduc-tivity of the n wire. Here we introduced the commutator[ˇ g, ˇ g N,S ] ◦ = ˇ g ◦ ˇ g N,S − ˇ g N,S ◦ ˇ g . The current in the n wireis determined by the formula j = σ e Tr { ˆ τ · π (ˇ g ◦ ∂ x ˇ g ) ( t, t ) } (14)The matrix element (ˇ g ◦ ∂ x ˇ g ) is the Keldysh compo-nent that equals (ˇ g ◦ ∂ x ˇ g ) = ˆ g R ◦ ∂ x ˆ g K + ˆ g K ◦ ∂ x ˆ g A .Even in a time-independent case, an analytical solu-tion of the problem can by found only under certainassumptions . In the considered case of a time-dependent phase difference, the problem becomes evenmore complicated. In order to solve the problem analyt-ically, we consider two limiting cases: a) weak proximityeffect and slow phase variation in time; b) strong prox-imity effect in a short n wire and arbitrary frequency Ωof the phase oscillations. III. WEAK PROXIMITY EFFECT; SLOWPHASE VARIATION
In this section we will assume that the proximity ef-fect is weak and the phase difference ϕ ( t ) is almost con-stant in time. The latter assumption means that thefrequency of phase variation satisfies the condition Ω ≪ T / ¯ h . The weakness of the PE means that the anoma-lous (Gor’kov’s) part ˆ f R ( A ) of the retarded and advancedGreen’s functions in the n wire ˆ g R ( A ) = g R ( A ) ˆ τ + ˆ f R ( A ) can be assumed to be small | ˆ f R ( A ) | ≪ . (15)The matrix ˆ f R ( A ) contains only off-diagonal elements.The diagonal part obtained from the normalization is g R ( A ) ˆ τ ≈ ± ˆ τ (cid:18) − (cid:16) ˆ f R ( A ) (cid:17) (cid:19) . (16)Now we can linearize Eq. (12) for the component 11(22),that is, for the retarded (advanced) Green’s functions.Then we obtain a simple linear equation ∇ ˆ f R ( A ) − κ ε ˆ f R ( A ) = 0 (17)where κ R ( A ) ε = p ∓ iε/D . The boundary conditions (13)for the matrices ˆ f R ( A ) acquire the form h ∂ x ˆ f R ( A ) + 2 κ N ˆ f R ( A ) i (cid:12)(cid:12)(cid:12) x =+ L x = 0; (18) h ∂ x ˆ f R ( A ) − κ N ˆ f R ( A ) i (cid:12)(cid:12)(cid:12) x = − L x = 0; (19) h ∂ y ˆ f R ( A ) − κ S (cid:16) ˆ f R ( A ) S, + ϕ ∓ g R ( A ) S · ˆ f R ( A ) (cid:17)i (cid:12)(cid:12)(cid:12) y =+ L y = 0; (20) h ∂ y ˆ f R ( A ) + 2 κ S (cid:16) ˆ f R ( A ) S, − ϕ ∓ g R ( A ) S · ˆ f R ( A ) (cid:17)i (cid:12)(cid:12)(cid:12) y = − L y = 0 . (21)As follows from Eq. (8) the functions g S , ˆ f S,ϕ are g R ( A ) S = ε/ξ R ( A ) ε , (22)ˆ f R ( A ) S,ϕ = ( i ˆ τ cos ϕ + i ˆ τ sin ϕ )∆ /ξ R ( A ) ε . (23)We took into account that ϕ ( t ) is almost constant in time.One can show that the solution for ˆ f R in the horizontalpart of the n wire is ˆ f R = i ˆ τ f ( x ) , (24) f ( x ) = C cosh( θ x x/L x ) + sgn( x ) S sinh( θ x x/L x ) (25)where sgn( x ) is the sign function. Dropping the index R of the quantities κ ε , g S ( ε ) , ξ ε the integration constants C and S can be written as C ( ε, ϕ ) = ( θ x cosh θ x + r N sinh θ x ) · r S ∆ cos ϕ/ξ ε D ( ε ) , (26a) S ( ε, ϕ ) = − ( r N cosh θ x + θ x sinh θ x ) · r S ∆ cos ϕ/ξ ε D ( ε ) (26b)where D ( ε ) = ( r S g S ( ε ) θ x + r N θ y ) cosh( θ x + θ y ) +( r S g S ( ε ) r N + θ x θ y ) sinh( θ x + θ y ), r N,S = 2 κ N,S L x,y , θ x,y = κ ε L x,y .Knowing the function ˆ f R ( x ), one can find the correc-tion to the conductance of the n wire due to the PE.In order to obtain the current, we take the (12) com-ponent (the Keldysh component) of Eq. (12), multiplythis component by ˆ τ and take the trace. In the Fourierrepresentation we get (compare with Eq. (2) of Ref. ) M − ( ε, ϕ, x ) ∂ x F − ( x ) = c ( ε, ϕ ) . (27)where the function M − ( ε, ϕ, x ) = 1 + ( f ( x ) + f ∗ ( x )) determines the correction to the conductivity caused bythe PE and c ( ε, ϕ ) is an integration constant that is re-lated to the current: j = σ e Z ∞−∞ dε c ( ε ) . (28)It is determined from the boundary condition that canbe obtained from Eq. (13) M − ( ε, ϕ, x ) ∂ x F − ( x ) = c ( ε, ϕ ) = ν [ F N − − F − ( L x )] . (29)where ν ( ε, ϕ ) = ℜ{ f ( L x ) } is the density of statesin the n wire near the nN interface. Finding F − ( L x )and c ( ε ) from Eq. (29), we obtain for the current density(compare with Eq. (13) of Ref. ) j ( ϕ ) = 12 e Z ∞−∞ dε F N − R nN /ν + R n h M ( ε, ϕ ) − i (30)Here F N − is defined according to Eq. (11), R n = L x /σ isthe resistance of the n wire of the length L x in the normal ∆ R nN R nN ∆ R n R n D- - - ∆ R * R * ´ FIG. 2: Bias voltage dependence of the normalized variationsof the resistance contributions δR nN /R nN and δR n /R n . Pa-rameter values: ϕ = π/ L y /L x = 1, ε N / ∆ = 2 . · − , ε S / ∆ = 5 · − , R n /R nN = 1. state, and h M ( ε, ϕ ) − i = (1 /L x ) R L x dx M − ( ε, ϕ, x ) − .The first term in the denominator is the nN interfaceresistance and the second term is the resistance ofthe (0 , L x ) section of the n wire modified by the PE.The expressions for the DOS ν ( ε, ϕ ) and the function h M ( ε, ϕ ) − i are given in appendix.For the differential conductance G = dj/dV at zerotemperature we obtain G ( V, ϕ ( t )) = { R nN /ν + R n h M ( eV, ϕ ) − i} − . (31)In Fig. 2 we show the dependence of the nN interfaceresistance variation δR nN = R nN /ν − R nN and the resis-tance variation of the n wire δR n = R n h M ( eV, ϕ ) − i− R n on the bias voltage V for a fixed phase difference. Itcan be seen that the δR nN is either positive or negativedepending on V , while δR n is always negative, i.e. thePE leads to voltage dependent changes of the interfaceresistance (caused by the changes of the DOS in the n wire) and to a decrease of the resistance of the n wire.The conductance variation δG = G ( V, ϕ ) − G n , isshown in Fig. 3 for various values of R nN /R N , where G n = { R nN + R N } − is the conductance of the n -wirein the normal state. These results have been obtainedearlier .We are interested in the dc conductance variation aver-aged in time: δG av = (Ω / π ) R π/ Ω0 dt δG ( V, ϕ ( t )). First,from Eqs. (25-26) we can extract the dependence of thefunction f on the phase ϕ : f ( x, ϕ ) = f ( x,
0) cos ϕ . Hencewe obtain M − ( ε, ϕ, x ) = 1 + δM − ( ε, , x ) cos ϕ where δM − ( ε, ϕ, x ) = M − ( ε, ϕ, x ) −
1. At the same time, ν ( ε, ϕ ) = 1 + δν ( ε,
0) cos ϕ with δν ( ε, ϕ ) = ν ( ε, ϕ ) − δG ( V, ϕ ( t )) = δG ( V,
0) cos ϕ ( t ) (32) abc 0.02 0.04 0.06 0.08 0.10eV D- - ∆ GG n ´ FIG. 3: Bias voltage dependence of the normalized con-ductance variation δG/G n . Parameter values: ϕ = π/ L y /L x = 1, ε N / ∆ = 2 . · − , ε S / ∆ = 5 · − . Differentcases: a) R n /R nN = 0 .
5, b) R n /R nN = 1, c) R n /R nN = 2. j W ∆ G av G n ´ FIG. 4: Bias voltage dependence of the normalized time-averaged conductance variation δG av /G n . Parameter values: ϕ = π/ L y /L x = 1, ε N / ∆ = 2 . · − , ε S / ∆ = 5 · − , eV / ∆ = 5 · − , R n /R nN = 1. which by averaging over time yields δG av = δG ( V, ·
12 (1 + J (2 ϕ Ω ) cos(2 ϕ )) (33)where J is the Bessel function of the first kind andzeroth order. This oscillatory behavior of the time-averaged (dc) conductance variation δG av as a functionof the ac amplitude can be seen in Fig. 4.Thus, the calculations carried out in this section un-der assumption of adiabatic phase variations allow us toobtain the dependence of the conductance change δG av on the amplitude ϕ Ω , but provide no information aboutthe frequency dependence of δG av . This dependence willbe found in the next section. IV. STRONG PE; SHORT NORMAL WIRE
In this section we analyze the limiting case of a short n wire when the Thouless energy E T h = D/L x is muchlarger than characteristic energies: E T h ≫ Dκ N,S , T, eV .In this case all the functions in Eq. (12) are almost con-stant in space and we can integrate this equation from { x, y } = ± { x, y } = ± L x,y over x and y coordinates.The term ˆ τ · ∂ t ˇ g + ∂ t ′ ˇ g · ˆ τ (in the Fourier representation − iε ˆ τ · ˇ g + iε ′ ˇ g · ˆ τ ) is considered as a constant and theterm with the voltage V is omitted because we neglectthe voltage drop over the n wire; the voltage drops acrossthe nN, nS interfaces and is set to be zero in the n wire.Performing this integration and summing up the results,we obtain2( L x + L y ) ˇ A = ˇ J x ( L x ) − ˇ J x ( − L x ) − ˇ J x (+0) + ˇ J x ( − J y ( L y ) − ˇ J y ( − L y ) − ˇ J y (+0) + ˇ J y ( −
0) (34)where ˇ J x ( x ) = D ˇ g ◦ ∂ x ˇ g | y =0 , ˇ J y ( y ) = D ˇ g ◦ ∂ y ˇ g | x =0 , andˇ A = − i ( ε ˇ τ · ˇ g − ˇ g · ˇ τ ε ′ ). Integration around the point( x, y ) = (0 ,
0) yields the conservation of the ”currents”(using terminology of the circuit theory )ˇ J x (+0) + ˇ J y (+0) = ˇ J x ( −
0) + ˇ J y ( −
0) (35)Combining Eqs. (34-35) and the boundary condi-tions (13), we arrive at the equation ε ˇ τ · ˇ g − ˇ g · ˇ τ ε ′ = iε N [ˇ g, ˇ g N + ] ◦ + iε S [ˇ g, ˇ g S + ] ◦ (36)Here ε N,S = D/ (2 R nN,nS σL ) is a characteristic energyrelated to the interface transparencies, L = L x + L y .The energy ε N determines the damping in the spectrumof the n wire and the energy ε S is related to a subgapinduced in the n wire due to the PE. The matrices ˇ g N,S + are equal to: ˇ g N,S + = [ˇ g N,S ( L x,y ) + ˇ g N,S ( − L x,y )].In the limit of the short n wire considered in thissection, we need to find only the retarded (advanced)Green’s functions. Indeed, let us rewrite the expressionfor the current (14) using the boundary condition (13)at the right nN interface and concentrating on the dccomponent of the current: j = 116 eR nN Tr { ˆ τ · ∞ Z −∞ dε (ˆ g R · ˆ g KN +ˆ g K · ˆ g AN − ˆ g RN · ˆ g K − ˆ g KN · ˆ g A ) } (37)where ˆ g R ( A ) N = ± ˆ τ and Tr { ˆ τ · (ˆ g R · ˆ g KN ) } = 4 g R F N − . Thedistribution function F N − in the N reservoir is definedin Eq. (11). The integral over energies from the secondand third term is zero because it is proportional to thevoltage in the n wire which is set to be zero. Thereforethe current can be written as j = 12 eR nN ∞ Z −∞ dε ν ( ε ) F N − ( ε ) (38)where ν ( ε ) = ( g R − g A ) = ℜ{ g R ( ε ) } . This formulahas an obvious physical meaning - the current throughthe nN interface is determined by the product of theDOS in the n wire and N reservoir ( ν N = 1) and thedistribution function in the N reservoir (the distributionfunction F − in the n wire is zero).Using Eqs. (2), (11), (38) we arrive at the followingexpression for the differential conductance: G = 12 R nN ∞ Z −∞ dε T ν ( ε ) " ε + eV T + 1cosh ε − eV T (39)In order to find the matrix ˆ g R , we can write the (11)component of Eq. (36) in the form˜ ε ˆ τ · ˆ g R − ˆ g R · ˆ τ ˜ ε ′ = iε S [ˆ g R , ˆ g RS + ] ◦ (40)where ˜ ε = ε + iε N , ˜ ε ′ = ε ′ + iε N .According to Eqs. (1), (8) the matrix ˆ g RS + is a functionof two times, ˆ g RS + ( t, t ′ ), that is, in the Fourier representa-tion it is function of two energies: ε, ε ′ . Therefore, to findthe matrix ˆ g R ( ε, ε ′ ) in a general case is a formidable task.However, we can assume that the amplitude of the accomponent of the phase ϕ Ω is small and obtain the solu-tion making an expansion in powers of ϕ Ω :ˆ g R = ˆ g R + δ ˆ g R + δ ˆ g R + . . . (41)Here and later all matrix Green’s functions writtenwithout arguments are functions of two energies ( ε, ε ′ ).Those of them which are diagonal in energy may be also(explicitly) written with a single energy argument, e.g.ˆ g RS = ˆ g RS ( ε )2 πδ ( ε − ε ′ ).Similar to Eq. (41) we represent the matrix ˆ g RS + (up tothe second order in ϕ Ω ) as ˆ g RS + = ˆ g RS + δ ˆ g RS + + δ ˆ g RS + and find from Eq. (8) for the stationary part ˆ g RS and thecorrections δ ˆ g RS + (first order in ϕ Ω ) and δ ˆ g RS + (secondorder in ϕ Ω ):ˆ g RS = 2 πδ ( ε ˆ τ + i ∆ cos ϕ ˆ τ ) ξ − ε (42) δ ˆ g RS + = − i ˆ τ π ϕ Ω ∆ sin ϕ ( ξ − ε + ξ − ε ′ ) ( δ Ω + δ − Ω ) (43) δ ˆ g RS + = − i ˆ τ π ϕ ∆ cos ϕ ( P + P ) (44)where we used the notation δ ω ≡ δ ( ε − ε ′ + ω ), ξ ε ≡ ξ Rε and defined the functions P = δ (2 ξ − ε ′ + ξ − ε ′ +Ω + ξ − ε ′ − Ω ) ,P = ( δ + δ − ) ( ξ − ε ′ + ξ − ε + 2 ξ − ( ε + ε ′ ) ) . (45)Using the expressions for δ ˆ g RS + and δ ˆ g RS + givenabove we can calculate the corrections to ˆ g R up to the second order in ϕ Ω and the corresponding modificationof the DOS δν in the n wire.In the zeroth-order approximation, i.e. for ϕ Ω = 0 weobtain from Eq. (40) ˆ g R ( ε, ε ′ ) = ˆ g R ( ε )2 πδ ( ε − ε ′ ) wherethe matrix ˆ g R ( ε ) obeys the equation[ˆ τ E Rε + i ˆ τ E Rsg , ˆ g R ] = 0 (46)containing E Rε = ˜ ε + iε S g RS ( ε ) = ε + iε N + iε S g RS ( ε ), E Rsg = iε S cos ϕ f RS ( ε ), g RS ( ε ) = ε/ξ Rε , f RS ( ε ) = ∆ /ξ Rε .The solution of this equation is ˆ g R ( ε ) = ˆ τ g R ( ε ) + i ˆ τ f R ( ε ); g R ( ε ) = E Rε /ζ Rε , f R ( ε ) = E Rsg,ϕ /ζ Rε (47)where ζ Rε = q ( E Rε ) − (cid:0) E Rsg (cid:1) . The quantity E Rsg determines a subgap induced in the n wire due to thePE. Indeed, consider the most interesting case of smallenergies assuming that { ε, ε S } ≪ ∆; then, ξ Rε ≈ i ∆, f RS ( ε ) ≈ − i and ζ Rε ≈ p ( ε + iε N ) − ( ε S cos ϕ ) . Thismeans that the spectrum of the n wire has the sameform as in the BCS superconductor with a damping ε N and a subgap ε S | cos ϕ | , which depends on the nS interface transparency and phase difference.Note that the formula for the subgap induced in theN metal due to the PE in a tunnel SIN junction wasobtained by McMillan .The obtained results for the functions g R ( ε ) and f R ( ε )can be easily generalized for the case of asymmetric nS interfaces with different interface resistances R nS , (cor-respondingly, ε S , ). In the limit ε S , ≪ ∆, we obtainfor the subgap ε sg ε sg ( ϕ ) = q ε S + ε S + 2 ε S ε S cos 2 ϕ (48)This formula shows that that the subgap as a functionof the phase difference ϕ varies from | ε S − ε S | for ϕ = π/ ( ε S + ε S ) for ϕ = 0.We proceed finding the corrections of the first ( δ ˆ g R )and second ( δ ˆ g R ) order in ϕ Ω for ˆ g R in a way similar tothe one used in . The correction of the first order δ ˆ g (for brevity, we drop the index R ) obeys the equation ζ ε ˆ g ( ε ) · δ ˆ g − δ ˆ g · ˆ g ( ε ′ ) ζ ε ′ = iε S [ˆ g , δ ˆ g S + ] ◦ (49)which contains all terms of the first order in ϕ Ω from Eq. (40). Note, that we are making use of therelation ˆ g ( ε ) = ζ − ε (˜ ε ˆ τ + iε S ˆ g S ( ε )) evident fromEqs. (40), (46-47).In order to solve Eq. (49), it is useful to employ thenormalization condition (4) for ˆ g ≡ ˆ g R which for thefirst-order term of ˆ g ◦ ˆ g yieldsˆ g ( ε ) · δ ˆ g + δ ˆ g · ˆ g ( ε ′ ) = 0 (50)From Eqs. (49-50), we find δ ˆ g = iε S δ ˇ g S + − ˆ g ( ε ) · δ ˆ g S + · ˆ g ( ε ′ ) ζ ε + ζ ε ′ (51)We determine the correction δ ˆ g in the same manner.Reading off the second-order terms in Eq. (40) gives[ ζ ε ˆ g , δ ˆ g ] ◦ = iε S (cid:16) [ˆ g , δ ˆ g S + ] ◦ + [ δ ˆ g, δ ˆ g S + ] ◦ (cid:17) (52)The second-order part of the normalization condition isˆ g ( ε ) · δ ˆ g + δ ˆ g · ˆ g ( ε ′ ) = − δ ˆ g ◦ δ ˆ g (53)Thus, we obtain the second-order correction δ ˆ g = iε S δ ˇ g S + − ˆ g ( ε ) · δ ˆ g S + · ˆ g ( ε ′ ) ζ ε + ζ ε ′ + (cid:26) iε S [ δ ˆ g S + , δ ˆ g ] ◦ ζ ε + ζ ε ′ − ζ ε ( δ ˆ g ◦ δ ˆ g ) ζ ε + ζ ε ′ (cid:27) · ˆ g ( ε ′ ) . (54)In order to calculate the correction to the dc conduc-tance caused by the ac radiation, δG , we need to findTr { ˆ τ · δ ˆ g } and Tr { ˆ τ · δ ˆ g } and take their parts propor-tional to 2 πδ ( ε − ε ′ ). By inspection of Eqs. (43), (51) onerecognizes that the first-order correction contains onlyterms proportional to δ ( ε − ε ′ ± Ω) and therefore onlycontributes to the ac current. This is the fundamental reason why the second-order analysis is needed todetermine the variation of the dc conductance.As a result we just have to find the multiple of2 πδ ( ε − ε ′ ) contained in Tr { ˆ τ · δ ˆ g } which we denote as2 δ dc g ( ε ), that is δ dc g ( ε )2 πδ ( ε − ε ′ ) := Tr { ˆ τ · δ ˆ g } dc .We represent the function δ dc g ( ε ) as a sum δ dc g ( ε ) = δ (0) dc g ( ε ) + δ (Ω) dc g ( ε ) (55)where the function δ (0) dc g ( ε ) originates from the first termin Eq. (54) and the function δ (Ω) dc g ( ε ) from the secondand third terms. If we consider the case when the subgap ε S | cos ϕ | is much less than ∆, i.e. ε S | cos ϕ | ≪ ∆ (56)then, at low energies ε < ∼ ε S , the function δ (0) dc g ( ε ) isalmost independent of Ω, whereas the function δ (Ω) dc g ( ε )depends strongly on Ω at ε ≈ ε S | cos ϕ | . Assuming thevalidity of Eq. (56) we obtain δ (0) dc g ( ε ) = − ε S ϕ cos ϕ g ( ε ) ζ ε , (57) δ (Ω) dc g ( ε, Ω) = 14 ε S ϕ sin ϕ X ± Ω g ( ε )[1 + f ( ε ) f ( ε +Ω) + g ( ε ) g ( ε +Ω)][ ζ ε + ζ ε +Ω ] + f ( ε )[ g ( ε ) f ( ε +Ω) + f ( ε ) g ( ε +Ω)] ζ ε [ ζ ε + ζ ε +Ω ] (58)where the functions g ( ε ), f ( ε ), ζ ε are defined inEq. (47). The sum sign index ” ± Ω” in Eq. (58) meansthat the given expression is added to the same one withthe negative frequency ( − Ω).Using the function δ dc g ( ε, Ω) we can calculate a cor-rection to the DOS δν ( ε, Ω) due to the PE and with theaid of Eq. (39) find the correction δG ( V, Ω) to the differ-ential dc conductance. As follows from Eq. (39), at zerotemperature the normalized differential dc conductance˜ G ( V, Ω) = G ( V, Ω) R nN is equal to˜ G ( V, Ω) ≡ ˜ G ( V ) + δ ˜ G ( V, Ω) = ν ( eV ) + δν ( eV, Ω) (59)with the definitions ν ( eV ) = ℜ ( g ( eV )) and δν ( eV, Ω) = ℜ ( δ dc g ( eV, Ω)).Using the obtained formula for g ( ε ) and δ dc g ( ε ) wecan calculate the conductance G and its correction δG due to microwave radiation for different values ofparameters (damping ε N , phase difference 2 ϕ , etc.).The dependence of the conductance in the absence of radiation G versus the applied voltage V is presentedin Fig. 5. We see that this dependence follows theenergy dependence of a SIN junction. In our case the n wire with an induced subgap plays a role of ”S” with adamping ε N in the ”superconductor”. Since the valueof the induced subgap, ε sg = ε S | cos ϕ | , depends on thephase difference 2 ϕ , the position of the peak in thedependence G ( V ) is shifted downwards with increasing ϕ . Note that in an asymmetrical system ( ε S = ε S )the lowest value of the subgap is not zero (cf. Eq. (48)).In Fig. 6 we show the bias voltage dependence ofthe conductance correction due to ac radiation δ dc G (coefficient in front of ϕ ) for different values of ϕ . Themagnitude and the position of the arising peaks dependstrongly on the values of the parameters, e.g. ϕ .By varying the stationary phase difference ϕ or thedamping ε N one can change the frequency dependenceof the correction δ dc G considerably. This is shown inFig. 7 and Fig. 8 respectively. One can see that if ε N ≪ ε sg ( ϕ ), then the dependence δ dc G (Ω) has a peak abc0.00 0.05 0.10 0.15 0.20 0.25 0.30eV D Ž FIG. 5: Normalized stationary differential conductance ˜ G versus bias voltage V . Parameter values: T / ∆ = 10 − , ε S / ∆ = 0 . ε N / ∆ = 10 − . Different cases: a) ϕ = π/ ϕ = π/
4, c) ϕ = 3 π/ abc0.05 0.10 0.15 0.20 0.25 0.30eV D- - ∆ G ~ j W FIG. 6: Normalized second-order correction of differentialconductance δ ˜ G versus bias voltage V . Parameter values: T / ∆ = 10 − , ε S / ∆ = 0 . ε N / ∆ = 10 − , Ω / ∆ = 5 · − .Different cases: a) ϕ = π/
8, b) ϕ = π/
4, c) ϕ = 3 π/ located at ≈ ε sg ( ϕ ) and split into two subpeaks. Thesplitting becomes more and more distinct with increasingbias voltage V. With decreasing ϕ and increasing ε N ,the form of this dependence changes significantly. Forexample, the resonance curve becomes broader withincreasing damping. Increasing temperature leads to asimilar effect as one can see in Fig. 9.In Fig. 10 we plot the normalized conductance correc-tion δ dc ˜ G ( ϕ ) as a function of ϕ for different values ofthe bias voltage V . At large V this dependence is closeto sinusoidal, but at smaller voltages the form of the pe-riodic function δ dc G ( ϕ ) becomes more complicated. V. CONCLUSION
We have calculated the change of the conductancein an S/N structure of the cross geometry under the abc0.00 0.05 0.10 0.15 0.20 WD ∆ G ~ j W FIG. 7: Normalized second-order correction of differentialconductance δ ˜ G versus ac frequency Ω. Parameter values: T / ∆ = 2 · − , ε S / ∆ = 0 . ε N / ∆ = 10 − , eV / ∆ = 10 − .Different cases: a) ϕ = π/
6, b) ϕ = π/
4, c) ϕ = π/ abc0.00 0.05 0.10 0.15 0.20 WD ∆ G ~ j W FIG. 8: Normalized second-order correction of differentialconductance δ ˜ G versus ac frequency Ω. Parameter values: T / ∆ = 2 · − , ε S / ∆ = 0 . eV / ∆ = 10 − , ϕ = π/ ε N / ∆ = 5 · − , b) ε N / ∆ = 10 − ,c) ε N / ∆ = 2 · − . influence of microwave radiation. The calculations havebeen carried out on the basis of quasiclassical Green’sfunctions in the diffusive limit. Two different limitingcases have been considered: a) a weak proximity effectand low frequency Ω of radiation; b) a strong proximityeffect and small amplitude of radiation.In the case a), the conductance change δG consistsof two parts. One is related to a change of the nN interface resistance due to a modification of the DOSof the n wire. At small applied voltages V N , it isnegative. Another part is caused by a modification ofthe conductance of the n wire due to the PE. This partis positive and consists of two competing contributions.One contribution, which is negative, stems from thea modification of the DOS of the n wire. Another abcd0.00 0.05 0.10 0.15 0.20 WD ∆ G ~ j W FIG. 9: Normalized second-order correction of differentialconductance δ ˜ G versus ac frequency Ω, Parameter values: ε S / ∆ = 0 . ε N / ∆ = 10 − , eV / ∆ = 10 − , ϕ = π/
3. Differ-ent cases: a) T = 0, b) T / ∆ = 2 · − , c) T / ∆ = 6 · − ,d) T / ∆ = 10 − . abc0.5 1.0 1.5 2.0 2.5 3.0 j - - ∆ G ~ j W FIG. 10: Normalized second-order correction of differentialconductance δ ˜ G versus stationary phase difference ϕ . Pa-rameter values: T / ∆ = 2 · − , ε S / ∆ = 0 . ε N / ∆ = 10 − ,Ω / ∆ = 10 − . Different cases: a) V = 0, b) eV / ∆ = 2 · − ,c) eV / ∆ = 4 · − . contribution is positive and caused by a term, which issimilar to the Maki-Thompson term . The conduc-tance change δG oscillates and decays with increasingamplitude of radiation.In the case b) a short n wire was considered so that theresistance of the n wire is negligible in comparison withthe resistance of the nN interface. The correction δG has been found under assumption of a small amplitudeof the radiation. We found that at small applied voltages V , the dependence δG (Ω) has a resonance form. It hasa maximum when the frequency Ω is of the order of ε sg = ε S | cos ϕ | where ε sg is a subgap in the spectrumof the n wire induced by the PE. With increasing V , thepeak in the dependence δG (Ω) splits into two peaks andoverall form of this dependence becomes complicated. We assumed that the nS interface resistance islarger than the resistance of the n wire, that is: ρL ≪ R nS . This inequality can be written in the form ε S ≪ ε T h ≡ D/L , that is, the subgap energy in the n wire is much smaller than the Thouless energy ε T h .In the opposite limit, ε S ≫ ε T h , a gap of the order of ε T h is induced in the n wire. This limit can be studiednumerically. However, one can expect that in this limitthe resonance should take place at ω res ≈ ε T h / ¯ h . Exper-iment performed in Ref. corresponds to this limit. Thefrequency corresponding to the Thouless energy in ex-periment is equal to ε T h /h ≈ π (400 / − ) s − ≈ GHz ,whereas the resonance frequency is ν res ≈ GHz .Note that we considered a simplified model. For exam-ple, we did not account for the change of the distributionfunction in the n film (heating effects). One can giveestimations when the ”heating” can be neglected. Thechange of an ”effective” electron temperature δT in the n wire is approximately given by: δT ≈ τ e − ph σE /c e ,where τ e − ph is the electron-phonon inelastic scatteringtime, E < ∼ δV S R L /R b L = ¯ h Ω( R L /eLR b ) ϕ Ω is theac electric field in the n wire and c e ≈ T · n/ε F isthe heat capacity of electron gas with concentration n . Taking into account that ( R b σ ) − ≈ Z /l , wefind that δT /T < ∼ ( ε /T ) ( τ e − ph /τ ) Z ϕ , where Z is the dimensionless coefficient of electron penetrationthrough the SN interface, which is assumed to be small, l = vτ is the mean free path in the n wire. There-fore, the heating would be very small if the condition ϕ Ω ≪ ( ε /T ) p ( τ /τ e − ph ) Z − is fulfilled.The obtained results are useful for understanding theresponse of the considered and analogous SN systems tomicrowave radiation which can be used, for example, inQ-bits.We would like to thank SFB 491 for financial support.One of us (VTP) was supported by the British EPSRC(Grant EP/F016891/1). VI. APPENDIX
The DOS in Eqs. (29-30) is given by the formula ν ( ε, ϕ ) = ℜ [1 + f ( L x ) ] with f ( L x ) defined in Eq. (25).Making use of the weak proximity effect approximationthe function we rewrite h M ( ε, ϕ ) − i in Eq. (30) as h M ( ε, ϕ ) − i = 1 − hℜ{ f + f f ∗ }i (60)0Using Eq. (25) one can easily calculate h f i = C + S θ x θ x + C − S CS sinh θ x θ x (61) h f f ∗ i = | C | + | S | θ ′ x θ ′ x + | C | − | S | θ ′′ x θ ′′ x + ℜ (cid:26) C ∗ S (cid:18) sinh θ ′ x θ ′ x + i sin θ ′′ x θ ′′ x (cid:19)(cid:27) (62) where θ ′ x and θ ′′ x are the real and imaginary parts of θ x respectively, i.e. θ x = θ ′ x + iθ ′′ x .We use these expressions for calculating the function h M ( ε, ϕ ) − i and conductance G . V. T. Petrashov, V. N. Antonov, and M. Persson, PhysicaScripta
T42 , 136 (1992), see also ”Low Dimensional Prop-erties of Solids” (Proceedings of the Nobel Jubilee Sympo-sium, G¨oteborg, Sweden, December 4-7,1991) ed. M. Jon-son and T. Claeson, World Scientific Publishing,1992. V. T. Petrashov, V. N. Antonov, P. Delsing, and R. Clae-son, Phys. Rev. Lett. , 347 (1993), ibid , 5268 (1995). H. Pothier, S. Gu´eron, D. Esteve, and M. H. Devoret, Phys.Rev. Lett. , 2488 (1994). P. G. N. de Vegvar, T. A. Fulton, W. H. Mallison, andR. E. Miller, Phys. Rev. Lett. , 1416 (1994). A. Dimoulas, J. P. Heida, B. J. v. Wees, T. M. Klapwijk,W. v. d. Graaf, and G. Borghs, Phys. Rev. Lett. , 602(1995). J. Eom, C.-J. Chien, and V. Chandrasekhar, Phys. Rev.Lett. , 437 (1998). C. W. J. Beenakker, in
Mesoscopic Quantum Physics (North-Holland, Amsterdam, 1995), ed. by E. Akkermans,G. Montambaux, J.-L. Pichard and J. Zinn-Justin. C. J. Lambert and R. Raimondi, J. Phys.: Condens. Mat-ter , 901 (1998). Y. V. Nazarov, Superlattices and Microstructures , 1221(1999). W. Belzig, F. K. Wilhelm, C. Bruder, G. Sch¨on, and A. D.Zaikin, Superlattices and Microstructures , 1251 (1999). A. Volkov, A. Zaitsev, and T. Klapwijk, Physica C: Super-conductivity , 21 (1993). F. W. J. Hekking and Y. V. Nazarov, Phys. Rev. Lett. ,1625 (1993). A. Zaitsev, Physica B: Condensed Matter , 274 (1994). Y. V. Nazarov and T. H. Stoof, Phys. Rev. Lett. , 823(1996). S. N. Artemenko, A. F. Volkov, and A. V. Zaitsev, SolidState Communications , 771 (1979). A. Volkov, N. Allsopp, and C. J. Lambert, Journal ofPhysics: Condensed Matter , L45 (1996). A. F. Volkov, Phys. Rev. Lett. , 4730 (1995). F. K. Wilhelm, G. Sch¨on, and A. D. Zaikin, Phys. Rev.Lett. , 1682 (1998). S.-K. Yip, Phys. Rev. B , 5803 (1998). V. N. Gubankov and N. M. Margolin, JETP Lett. , 673(1979). P. Charlat, H. Courtois, P. Gandit, D. Mailly, A. F. Volkov,and B. Pannetier, Phys. Rev. Lett. , 4950 (1996). V. Petrashov, R. Shaikhaidarov, and I. Sosnin, JETP Let-ters , 839 (1996). J. J. A. Baselmans, A. F. Morpurgo, B. J. van Wees, andT. M. Klapwijk, Nature , 43 (1999). R. Shaikhaidarov, A. F. Volkov, H. Takayanagi, V. T. Pe- trashov, and P. Delsing, Phys. Rev. B , R14649 (2000). A. A. Golubov, M. Y. Kupriyanov, and E. Il’ichev, Rev.Mod. Phys. , 411 (2004). A. I. Buzdin, Rev. Mod. Phys. , 935 (2005). F. S. Bergeret, A. F. Volkov, and K. B. Efetov, Rev. Mod.Phys. , 1321 (2005). C. C. Tsuei and J. R. Kirtley, Rev. Mod. Phys. , 969(2000). D. J. Van Harlingen, Rev. Mod. Phys. , 515 (1995). A. Bezryadin, C. N. Lau, and M. Tinkham, Nature ,971 (2000). M. Zgirski, K.-P. Riikonen, V. Touboltsev, and K. Aru-tyunov, Nano Letters , 1029 (2005). M. Tian, N. Kumar, S. Xu, J. Wang, J. S. Kurtz, andM. H. W. Chan, Phys. Rev. Lett. , 076802 (2005). K. Arutyunov, D. Golubev, and A. Zaikin, Physics Reports , 1 (2008). V. T. Petrashov, K. G. Chua, K. M. Marshall, R. S.Shaikhaidarov, and J. T. Nicholls, Phys. Rev. Lett. ,147001 (2005). F. Chiodi, M. Aprili, and B. Reulet (2009),arXiv:0908.1070. C. Checkley, A. Iagallo, R. Shaikhaidarov, J. T. Nicholls,and V. T. Petrashov (2010), arXiv:1003.2176. D. Golubev and A. Zaikin, Europhys. Lett. , 37009(2009). P. Virtanen, T. T. Heikkil¨a, F. S. Bergeret, and J. C.Cuevas (2010), arXiv:1001.5149. F. Giazotto, T. T. Heikkila, A. Luukanen, A. M. Savin,and J. P. Pekola, Rev. Mod. Phys. , 217 (2006). K. D. Usadel, Phys. Rev. Lett. , 507 (1970). A. I. Larkin and Y. N. Ovchinnikov, Sov. Phys. JETP ,1200 (1969), [Zh. Eksp. Teor. Fiz. , 2262 (1968)]. J. Rammer and H. Smith, Rev. Mod. Phys. , 323 (1986). N. B. Kopnin,
Theory of Nonequilibrium Superconductivity (Clarendon Press, Oxford, 2001). V. R. Kogan, V. V. Pavlovskii, and A. F. Volkov, Euro-phys. Lett. , 875 (2002). A. F. Volkov and V. V. Pavlovskii, in
Proceedings of the21st Rencontres de Moriond (Les Arcs, France, 1996), p.267. A. A. Golubov, F. K. Wilhelm, and A. D. Zaikin, Phys.Rev. B , 1123 (1997). S. N. Artemenko and A. F. Volkov, Soviet Physics Uspekhi , 295 (1979). A. I. Larkin and Y. N. Ovchinnikov, Sov. Phys. JETP ,960 (1975), [Zh. Eksp. Teor. Fiz. , 1915 (1975)]. M. Y. Kupriyanov and V. F. Lukichev, Sov. Phys. JETP , 1163 (1988), [Zh. Eksp. Teor. Fiz. , 139 (1988)]. W. L. McMillan, Phys. Rev. , 537 (1968). A. V. Zaitsev, A. F. Volkov, S. W. D. Bailey, and C. J. Lambert, Phys. Rev. B60