Current-phase relation in a short clean josephson junction model: application to MgB2

Motivated by recent data on high-quality MgB2 thin films implying that the smaller energy gap has l = 6 (i-wave) symmetry, we consider a simple model for an all-MgB2 symmetric Josephson Junction (JJ). The model assumes an arbitrary-strength delta-function barrier and one-dimensional current conduction. It is shown that in this context a nodal energy gap with i-wave symmetry acts as an isotropic energy gap (s-wave) with an amplitude modified by the energy-gap misalignment-angle with respect to the crystal principal axes. The corresponding exact Green’s function in momentum space is derived employing a novel approach. The ensuing current-phase relations in the strong and weak barrier-strengths limits are calculated and found to confirm known results, e.g., the Ambegaokar-Baratoff current-phase relation. Inspired by an HTS experiment that established the d-wave energy-gap symmetry, we propose a JJ-related experiment with a MgB2 bicrystal to confirm our premise that the smaller energy has i-wave symmetry.


Introduction
An important aspect of the MgB 2 superconductor is its potential technological advantages over the cuprate hightemperature superconductors (HTS). Specifically, its relatively high critical temperature ( » T 40 K C ) simplifies cryogenic requirements, its relatively long coherence length avoids the problem of grain boundaries acting as weak links, and its stoichiometry and crystal structure are simple [1][2][3][4][5]. These advantages continue to stimulate substantial work on power applications such as conductors, RF cavities, and digital applications to mention a few [6][7][8][9][10][11].
From a physics point of view MgB 2 is interesting in at least two respects. It has two energy gaps that are susceptible only to weak coupling by non-magnetic impurities and, allegedly, these two energy gaps have distinct symmetries [12][13][14][15]. By convention, the two energy gaps are labelled as the s (larger) and the p gap (smaller). Whereas the symmetry of the s energy gap is s-wave [1,3], recent low-temperature microwave-frequency measurements of the intermodulation distortion (IMD) strongly suggest for the p energy gap an i-wave (l=6) symmetry [13]. The IMD low-temperature variation provides a signature for a nodal energy gap [16,17]. This signature been observed in the cuprate YBCO which indeed has a nodal (d-wave) energy-gap. The measured low-temperature IMD of MgB 2 is strikingly similar to that of YBCO, which strongly suggests a nodal p energy gap, specifically with i-wave symmetry [15]. Other observables, such as the low-temperature penetration-depth anomaly and the temperature dependence of the microwave surface resistance provide additional support for the i-wave symmetry [13]. With its two energy gaps having different symmetries MgB 2 is so far unique.
Theoretical calculations predict and tunneling experiments confirm that the p energy gap plays a prominent role in current conduction [11,18,19]. Specifically, current conduction in both theĉ direction and in directions close but not equal to the ab plane is expected to be dominated by the p energy gap. This is consistent with the three-dimensional Fermi surface that underlies the p energy gap, in contrast to the approximately twodimensional Fermi surface underlying the s energy gap [1,20,21]. Given the π-gap role in current conduction and the evidence of its nodal symmetry, we consider here the impact of the i-wave symmetry on the current-phase relation (CPR) in all-MgB 2 Josephson junctions (JJ) [22][23][24]. To the best of our knowledge such an analysis does not exist. In view of the recent demonstration of MgB 2 /MgO/MgB 2 junctions with a high critical current, the CPR calculation of an all-MgB 2 junction seems timely [11].
The CPR is calculated in a model that assumes a delta-function barrier and one-dimensional conduction in the direction normal to the barrier's plane. One-dimensional current conduction and this type of barrier have been employed before in the context of the Blonder-Tinkham-Klapwijk (BTK) analysis of SN tunneling with conventional superconductors, and for simulating multiple Andreev reflections in the barrier domain [25][26][27]. In this sense, our work represents an extension of this research.
As shown in section 2, a key approximation underlying this work is the reduction of a nodal energy gap to a constant (isotropic) energy gap with a properly modified amplitude. In particular, for an i-wave-symmetry energy gap this reduction takes the form is the i-wave energy gap,  K denotes the relative pair-momentum in the ab plane, j is the azimuthal angle, D 0 is the energy-gap amplitude, and a is the barrier misalignment angle with respect to the crystal axes in the xy plane. When applied to a d-wave energy gap, the same reduction procedure yields In the HTS context, the reduction (1.2) yields a CPR of structure consistent to that based on group-theoretic considerations; see appendix A [28,29]. The energy-gap reductions, equations (1.1), (1.2), simplifies the analysis substantially.
In the literature there are several approaches to the calculation of the CPR for s-wave superconductors, e.g., the Bogoliubov-de Gennes equations, the Green's function approach and the Boltzmann equation [30][31][32][33]. The choice of approach is a matter of mathematical convenience for the objective at hand. In this work, we adopt the Green's function approach, originally devised to treat spatially varying perturbations in superconductors, such as a non-constant energy gap and impurities [34][35][36]. The JJ configuration fits into this category since the energy gaps on both sides of the barrier are different. There has been work to construct a real-space JJ Green's function [37,38]. It is shown here that working in momentum space yields the exact Green's function for an arbitrary barrier strength. The ensuing CPR in the weak and strong barrier-strength limits is consistent with known results, e.g., Ambegoakar-Baratoff [33]. We emphasize that while the Ambegoakar-Baratoff CPR is applicable in the limit where the coupling between the two banks is weak (or a strong barrier-strength), the present Green'sfunction approach yields, in principle, the CPR for an arbitrary barrier strengths. The premise of onedimensional current conduction normal to the barrier's plane, common to a wide body of works, is further discussed in the concluding section [39].
This paper is organized as follows. In section 2 the nodal energy-gap reduction, equation (1.1) is derived. Section 3 is devoted to the derivation of the exact Green's function in momentum space, and evaluating the ensuing CPR in the weak and strong barrier-strength limits. In section 4 we outline a generalization of the Green's function approach to accounts for weak coupling between the two energy gaps, e.g., by doping , suggestions for JJ-based experiments to confirm the i-wave p energy-gap symmetry, and a summary.

Nodal energy-gap reduction
In momentum space and cylindrical coordinates the nodal i-wave energy gap, ( ) ,the corresponding real-space functional form is , sin , and denoting by g the azimuthal angle of the vector | |( ( ) ( ))    r r r g g -¢ = D cos , sin , it follows that · ( Employing the expansion [40] 3) x 0 and k F denote the = T 0 K coherence length and the Fermi momentum, respectively. Consider first the azimuthal angle g. The vector   r r -¢ may point in any direction in the xy plane.
However, under the assumption of one-dimensional current flow normal to the barrier's plane (at = x 0) it follows that only vectors » sin 6 sin 6 . Next consider the divergent dimensionless integral in equation (2.3), plotted in figure 1 for values pertinent to MgB 2 , i.e., x = k 40 F 0 and a chosen large cutoff [3]. As figure 1 shows, this integral diverges near the origin, however not at the origin. This real space divergence resembles the divergence of s-wave energy-gap in real space, as embodied by a delta function [34]. Here, however, and as expected for a nodal energy gap in real space, it has a finite range on atomic scale [36]. In the present context of a global property, i.e., the CPR, we adopt the approximation of ignoring the finite-range aspect of the nodal energy gap. Specifically, the divergent integral in equation (2.3) is approximated as˜( where the k F 2 factor in (2.4) has been added to keep the correct dimensionality of the integral (ℓ -2 where ℓ denotes the length dimension).
The overall minus sign in equation (2.3) is of no consequence since, as shown in section 3 and appendix A only products of two energy gaps enter the CPR expression. Consequently the minus sign in equation (2.3) is henceforth ignored. All these considerations yield for the i-wave (l=6) energy gap, equation (2.3), the reduced form Equation (2.5) implies that in the present context, a nodal i-wave energy gap acts as a isotropic (s-wave) energy gap, however, with an amplitude modified according to the misalignment angle a. When applying the same line of arguments for a d-wave energy gap, as in YBCO, we obtain equation (1.2). As shown in next section and appendix A, the CPR is proportional to the product of the two energy-gaps on both JJ banks. Hence equation (1.2) leads to a CPR expression identical to that obtained from group-theoretic considerations [29].  and upper integrationlimit cutoff at˜= Q 10. Increasing the cutoff sharpens and increase the peak near (however, never at)˜= Q 0. [25,26]. To the best of our knowledge, our approach to the CPR calculation differs from the large body of literature on s-wave JJ [33,39]. Specifically, starting from the integral form of the Gorkov equation we solve for the exact JJ Green's function for an arbitrary barrier strength. Subsequently we derive the ensuing CPR in the common limits of strong and weak barrier strength [35]. A bonus of the present calculation is an explicit expression for the JJ critical current in terms of the barrier's strength. In this section we outline the approach and quote results. Details are deferred to appendix A.

The exact JJ Green's function
In the one-dimensional coordinate space along the x axis the integral form of the Gorkov equation for the Green's function ( , ; n and current-density expression ( ) j x are , for all n values.
is the unperturbed Green's function that pertains to all interactions except pairing. In particular it includes the barrier. The symbol ( ) D x denotes the spatially dependent energy gap, q S is the singlecarrier charge (positive or negative), denote the Boltzmann constant and temperature, respectively,  is the normalized Planck's constant, and w n denotes the Matsubara frequencies [41]. The phase-discontinuity variable in (3.1) is suppressed to simplify notation. The barrier ( ) The barrier transparency is determined by the barrier's strength parameter v , B i.e., a 'large' v B value corresponds to weak transparency and vice-versa for 'small' v B value. Motivated by the expressions in appendix A and previous work these limits are quantified according to the dimensionless barrier-strength parameter Z where m M , denote the chemical potential and the band-structure carrier's mass, respectively [e.g., [25,42]. Correspondingly,  Z 1 corresponds to the weak barrier transparency while  » Z 1 0corresponds to the strong barrier transparency limits.
The two inputs to equation and ( ) D x to which we turn now. For a clean metal and the delta-function barrier (3.2), the explicit expression of the Green's function is [43] ( ) For the sake of conciseness, the w n dependencies ofl l , in equation (3.4) have been suppressed. The arbitrary square-root sign ofl l , has been chosen to guarantee a negative imaginary part for all w n frequencies. The other input to equation (3.1) is the energy gap ( ) D x . In view of the approximation (2.5) we adopt the where ( ) q x denotes Heaviside step-function, y is the phase-discontinuity across the barrier, and D D , L R are the energy-gap amplitudes on the left and right sides of the barrier, respectively. While the approach below applies for arbitrary D D , L R values (appendix A) here we restrict ourselves to the particularly interesting cases when The h = -1 case in (3.6) corresponds to two misalignment angles such that a positive energy gap on one JJ bank is aligned opposite to an equal and negative energy gap on the other bank. Such p junctions have been demonstrated in for YBCO, with a d-wave nodal energy-gap [44]. Thus at the = x 0 plane there are three discontinuities: That of the delta-function barrier, the phase y and the energy-gap amplitude.
Solving Gorkov's equation (3.1) in real space is difficult [37,38]. However, in momentum space an exactsolution of ( ) w y G k q , ; ; n is possible, see appendix A. Recalling the definition of the double Fouriertransform of a function ( ) i k x q x it follows from (3.1) that the CPR is given by In equation (3.8) the current density is identified with its value at the barrier's location, = x 0. In principle the current density is location independent provided self-consistent energy gaps are employed [32,38]. However, when a phenomenological energy gap is employed, as is commonly the case, the charge-current density continuity equation has a source-term that involves the product ( ) ( ) To avoid this term and yet satisfy the continuity equation the current density should be evaluated at a point where ( ) D = x 0, i.e. at the barrier's location = x 0. Following appendix A and equation (3.6), the explicit expression for (   for an arbitrary barrier strength Z. The Green's function, equation (3.9), is a central result of this work. The corresponding CPR is given by equation (3.8) . Here we consider the CPR in the common limits when  Z 1 and  » Z 1 0,corresponding to short SIS and SNS junctions respectively, Details of the CPR calculation are deferred to appendix A. The final expression turns out simple, see below, since most terms in (3.8) drop out due to one of the following three reasons: 1. All terms involving odd k q , powers drop in the k q , integration. 2. When transforming the integration in the remaining terms to relative- drop in the xintegration.
3. All terms odd in the Matzubara frequency w n drop upon the frequency summation. As expected, the only surviving terms involve the phase discontinuity y.
Consider first the weak barrier-transparency limit when  Z 1. As shown in appendix A, in this limit, to the lowest order in  m D 1 0 and Z 1 parameters it follows that  » R 0 and   » V b (equation (3.11). The ensuing calculation yields the result: are defined in equations (3.6) and (1.1), respectively. The temperature dependence of ( ) a D 0 is suppressed to keep notation concise. The CPR (3.12) is the text-book Ambegaokar-Baratoff result with the added bonus of an explicit expression for the critical current in terms of the barrier's strength [45]. Note from (3.12) that for h = -1, i.e. for p junctions, the current flows in the 'wrong' direction [46].
In the opposite limit, when  Z 0, i.e., the case of strong barrier transparency as in a SNS junctions with a thin, clean, metallic barrier the calculation is quite different. In this limit we consider only the h = 1 case. As shown in equation (A.10) and thereafter, only two diagonal  R matrix-elements equation n and the Matsubara-frequencies summation in (3.13) yields The CPR (3.14) is identical with that derived in previous work [30]. The fact that the exact Green's function (3.9) yield the correct results in the two common barrier-strength limits lends support to its veracity.

Generalization to include intra-band coupling
The Green's function derived above pertains to a JJ with single energy gap superconductors. This is justified for MgB 2 in the presence of non-magnetic impurities since the two energy gaps have distinct underlying band structure [12]. In instances where that coupling cannot be ignored the approach based on the Gorkov equation (3.1) requires generalization. Details are deferred to the appendix B. Labelling the two energy gaps by 1 and 2, this generalization takes the form of two integral equations, formally written as: . For 'small' | | I , 2 as is the case in MgB 2 , the correction terms in equation (4.1) can be approximated by using the uncoupled-system solutions.

Proposed experiments
Whereas in our view the data supporting an i-wave symmetry for the p energy gap is convincing, this may be the place to suggest further experiments to test this symmetry [13]. Based on experience gained with the exploration of the energy-gap symmetry in HTS, the suggestions are: 1. to make a corner SQUID made of a single MgB 2 crystal and measure the current variation as function of the magnetic flux through the SQUID's loop and 2. to make a p-JJ with crystallographically properly aligned banks and measure the spontaneous current in a loop that contains one such JJ [46][47][48][49]. The corner-SQUID approach requires a single crystal sufficiently large and thick to allow attachment of two electrodes on points on itsãb facets where the energy gap is either positive or negative. The presumed i-wave energy-gap symmetry, equation (1.1), implies energy-gap maxima at misalignment angles    15 , 75 , 135 ,... off the principal crystalx -axis, and corresponding energy-gap minima at - - - 15 , 75 , 135 ... The spread of either maxima or minima is  30 , hence the precision tolerance on the electrodes' attachment-points seems achievable. A schematic p-JJ is depicted in figure 2. In analogy to the tricrystal experiments with YBCO, the two JJ banks are single crystals where the two banks are aligned such that an energy-gap maximum on one bank faces an energy-gap minimum on the other bank [44,50]. This experiment requires adequately aligned substrates on which to grow epitaxial MgB 2 .

Summary
The simplicity of the present JJ model, i.e., the assumption of one-dimensional conduction and a delta-function barrier, comes at a price. The one-dimensional conduction ignores the angular-spread of tunneling carriers around the normal to the barrier plane. Denoting that angular-spread by ( ) d q , an ad-hoc way to account for this off-normal tunneling is to multiply the one-dimensional current-density expression by the factor ( ) ( ) d q p k 2 F 2 since only the forward Fermi-sphere hemisphere contributes [51]. Regarding the assumed deltafunction barrier, this assumption excludes processes within the barrier such as dissipation, multiple Andreev reflections and magnetic-field effects such as vortex motion [11,52,53].
In summary, in the context of a JJ and under the approximation of one-dimensional current conduction normal to the barrier's plane, we show that a nodal energy gap with an i-wave symmetry is equivalent to an isotropic energy gap with an amplitude modified according to the misalignment angle with respect to the crystal's axes. This result and mathematical convenience suggest a focus on a JJ model with a delta function barrier of arbitrary strength. In this respect, the present work generalizes the BTK model, commonly employed for the analysis of SN tunnel junctions [25]. An explicit expression for the exact JJ Green's function is derived. The ensuing CPR is calculated in the two common limits of weak and strong barrier-strengths where we recover known CPRs, e.g., the Ambegaokar-Baratoff relation, in the strong barrier-strength limit, with an explicit expression for the critical current in terms of the barrier's strength. While this work presupposes a single energy gap, as suggested by weak coupling between the energy gaps induced by non-magnetic defects, we outline an approach to account for such weak coupling. Following HTS experiments that established its energy-gap symmetry, we propose corresponding experiments with MgB 2 single crystals.
The authors acknowledge the support of the US Office of Naval Research and the Carderock Division of the Naval Surface Warfare Center's In-house Laboratory Independent Research Program sponsored by the Office of Naval Research administered under Program Element 0601152N.

Solution of Gorkov's equation (3.1)
To evaluate the complicated kernel ( ) w y K x x , ; ; , n 1 equation (3.1), it is first decomposed into spatial sectors by writing it in the form . A.1 breaks into the sum of a kernel pertaining to the underlying superconductor, K , S C and a term pertaining to the JJ, K , J J that includes the phase discontinuity.
Explicitly  Decomposition (A.2) guarantees that in the absence of the phase-discontinuity y the underlying superconductor kernel is automatically recovered.
In momentum space, equation (3.7), the Gorkov's equation (3.1) takes the form where dependencies on some variables in have been suppressed to simplify notation. Employing (3.4), the integrations in (A.2) yield the following explicit expressions:    A key advantage of the above Fourier transformation is that all components in equations (A.6), (A.7) are separable in the k, q variables. This property allows for the solution of the unknowns (3.10) in closed form, see equation (3.9).
In the 'Born' approximation we take , implies that the delta-function term dominates for low and high Z values. Correspondingly, the Born approximation (3.11) in this temperature regime is Next we evaluate the matrix  R in equation (3.11). Consider first the weak transparency regime,  ¥ Z . In that limit ( Consequently, for the matrix  R, equation (3.11), = -= - Likewise in the opposite limit, when estimates of the remaining matrix elements = -R R P S Q S , , implies that they can be neglected as well. All in all it follows that for ( )   ¥ » R Z 0 and the Born expressions for ( ( ) ( ) ( ))  = V S q P q Q q , , provide an excellent approximation. The corresponding CPR is the Ambegaokar-Baratoff result.
Consider now the opposite limit, when  Z 0, e.g., SNS junctions with a thin N barrier. We consider only the h = 1 case. In this limit, the Born approximation, eqaution (A.8), remains unchanged. However, the unknown ( ) S q drops out since to lowest order in Z we have ( ) = F v 0 B and the matrix  R, equation (3.11), The terms that contribute to the matrix  R are those that involve the phase y. Explicitly, where the relations have been employed. This yields for I D the estimate which leads to the CPR equation (3.14).
In this part of the appendix we outline generalization of the above approach to the general case where. Equation (3.5) holds yet not equation (3.6). Define the following additional (divergent) quantities As the K JJ term in equation (A.14) indicates, the CPR is proportional to the product to D D , L R consistent with previous results based on group-theoretic considerations [29]. After some rearrangement the kernel K SC can be rewritten as follows  The first term in the RHS of (A.15) has an analytical Fourier transform, equation (A.6). Consider now the particular energy-gap relation, equation (3.6). When h = 1, which is the symmetrical JJ, the second and third terms on the RHS of (A.15) vanish. When h = -1, which is a symmetrical p-JJ, the term proportional to ( ) d D 2 vanishes while the third term on the RHS of (A.15) is proportional to -D = -D 2 2 , L R 2 2 as in equation (3.9). In the corresponding  ¥ Z , the structure of the CPR, equation (3.12) remains unchanged. When the energy gaps of both banks are not related as in equation (3.6), e.g., in a JJ where one bank is nodal while the other is a proper s-wave superconductor the above-generalized approach should be applied [47].
where t = it denotes the imaginary time, and t T is the time-ordering operator [35]. where pairing between carriers belonging to different bands is neglected [54]. The symbol (ˆ)  H p B in equation (B.2) denotes the band-structure Hamiltonian.