A Expressions for the Fermi-surface average

Inverse proximity effect and influence of disorder on triplet supercurrents in strongly spin-polarized ferromagnets

Abstract

We discuss the Josephson effect in strongly spin-polarized ferromagnets where triplet correlations are induced by means of spin-active interface scattering, extending our earlier work [R. Grein et al., Phys. Rev. Lett. 102, 227005 (2009)] by including impurity scattering in the ferromagnetic bulk and the inverse proximity effect in a fully self-consistent way. Our quasiclassical approach accounts for the differences of Fermi momenta and Fermi velocities between the two spin bands of the ferromagnet, and thereby overcomes an important short-coming of previous work within the framework of Usadel theory. We show that non-magnetic disorder in conjunction with spin-dependent Fermi velocities may induce a reversal of the spin-current as a function of temperature.

pacs:
72.25.Mk,74.50.+r,73.63.-b,85.25.Cp

I Introduction

Ferromagnet-superconductor hybrid systems are currently subject to intense research activity, as they were conjectured to host triplet superconductivity induced by the proximity effect(1); A. I. Buzdin (2005); (3); (4); (5). While the first successful experiments on these structures found evidence for so-called triplet correlations(6); (7); (8), whose hallmark are -oscillations A. I. Buzdin (2005) of the critical current in Josephson junction devices, the existence of equal spin triplet pairing is currently in the focus of attention(1); (5). So far, the main experimental prove of such triplet correlations is based on the Josephson effect and on phase coherent electron transport in proximity structures (9); (10); (11); (12); (13); (16); (14); (15). After two first confirmations of long-range triplet amplitudes in 2006, (9); (10) an impressive series of affirmative experimental results was published in 2010.(11); (12); (13); (14); (15) The systems under investigation varied largely in terms of materials and fabrication, the common idea being that a breaking of spin-rotation symmetry around the bulk magnetization axis must somehow be enforced at the superconductor (SC)-ferromagnet (FM) contacts - which is in line with the general theoretical concepts behind this effect(17); (1); (18); (19); (20); (21). In this article, we elaborate on our earlier work concerning long-range triplet supercurrents(22), where we used a recent extension of the quasiclassical Green function technique(23), that allows us to consider ferromagnets whose exchange splitting is of the same order of magnitude as the Fermi energy consistently within quasiclassical theory(22); (24). In the following we include additional relevant effects like the suppression of the superconducting energy gap in proximity to the interfaces, the induced magnetization in the superconductor close to the ferromagnet, and disorder in the FM bulk. These calculations are performed numerically and cover all ranges of junction length, interface transparency, and impurity scattering. In the appendix, we provide detailed discussions of the intricacies involved in solving the self-consistency equation in quasiclassical theory with spin-active boundary conditions. In particular, we point out a previously not discussed technical problem with the self consistent solution of the order parameter that arises for quasiclassical transport equations with boundary conditions where the quasiclassical propagator undergoes a jump discontinuity during reflection from the interface.

With regard to the inverse proximity effect, we show that while the suppression of the gap function in the SC-electrodes is indeed substantial for highly transparent junctions, this does not imply that higher harmonic contributions to the current-phase relation (CPR) are necessarily suppressed. As expected, impurity scattering reduces the Josephson current and in particular higher harmonic contributions to the CPR. We also find that there may be a reversal of the Josephson spin-current as a function of temperature if impurity scattering is sufficiently strong.

Ii Quasiclassical Theory of Superconductivity

Let us first summarize the theory for a system where the band structure is not spin-polarized. The quasiclassical Green function is obtained from Gor’kov’s Green function(26) by an integration over momentum(27); (28); (29); (30):

 ^g(→pF,→R,ϵn,t)=1a(→pF)∫dξp^τ3^G(→p,→R,ϵn,t), (1)

where is a quasiparticle renormalization parameter, and . The “hat” denotes a 44 matrix in combined spin- and Nambu-Gor’kov space. The propagator depends on the spatial coordinate , time , Matsubara-frequency(31) and the momentum , which lies on the Fermi-surface. The Fermi velocity corresponding to Fermi momentum is denoted by . The resulting equation of motion, which can be obtained from the Dyson-equation of Gor’kov’s Green function, is the Eilenberger-equation(27):

 iℏ→vF⋅∇→R^g+[iϵn^τ3−^Δ−^Σ,^g]=^0, (2)

which must be supplemented by a normalization condition: . denotes off-diagonal self-energies in particle-hole space and the diagonal ones.

The equilibrium current density in a system without spin polarization reads:

 →j(→R,t)=eNFkBT∑ϵn12Tr4⟨→vF(→pF)^τ3^g(→pF,→R,ϵn,t)⟩, (3)

where Tr denotes a trace over spin and particle-hole (Nambu-Gor’kov) degrees of freedom, and and are defined by:

 ⟨∙⟩ = 1NF∫FSd2pF(2πℏ)3|→vF(→pF)|(∙), (4) NF = ∫FSd2pF(2πℏ)3|→vF(→pF)|. (5)

is the density of states per spin at the Fermi level. In appendix A we detail our numerical Fermi-surface averaging procedure.

For a ferromagnet these expressions need to be modified in order to take into account the different Fermi surfaces for the two spin channels. As discussed in our earlier work(22); (24), we model the ferromagnet by two spin-scalar quasiclassical Green functions for each band. This is consistent with the quasiclassical approximation in the limit and naturally results in a three-channel matching problem at the SC/FM interface. We underline that an equivalent approach in the diffusive limit of quasiclassical theory (Usadel-equation (25)) is currently not available for lack of multi-channel boundary conditions. All definitions given above remain the same apart from the reduced matrix structure of the Green function and the fact that , , and now depend on the spin band in question.

We will make extensively use of a fundamental symmetry relating particle-like and hole-like quantities in quasiclassical theory, for which we introduce the -operation defined by

 ~Q(→pF,→R,ϵn,t)=[Q(−→pF,→R,ϵn,t)]∗. (6)

ii.1 Riccati parameterization

We use a particular Riccati parameterization of the quasiclassical (QC) Green function(32); (33); (34); (35); (23) which has proven very useful in the past. For Matsubara’s Green functions,(31) it consists of 2 parameters (“coherence functions”) , , describing particle-hole coherence in the superconducting state. The Green function reads in terms of these parameters(23):

 ^g=∓2πi(GF−~F−~G)±iπ^τ3, (7)

where , , and correspond to and , respectively. The transport equations are:

 (iℏ→vF⋅∇+2iϵn)γ=[γ~Δγ+Σγ−γ~Σ−Δ], (8) (iℏ→vF⋅∇−2iϵn)~γ=[~γΔ~γ+~Σ~γ−~γΣ−~Δ].

For a superconductor or normal metal the , , , , , and are 2x2 spin matrices for each Fermi surface point , as all Fermi surfaces can be regarded as doubly spin degenerate on the scale of the Fermi energy. A weak band splitting due to an external field or as in ferromagnetic alloys can be incorporated as source field in these equations, leading to spin-dependent self energies. In contrast, for a ferromagnet with strong spin splitting of the energy bands the transport equations hold separately for each spin band with corresponding Fermi surface point , and the above-mentioned quantities are all scalars. The important difference lies in the integration over in equation (1), which destroys all coherence between quasiparticle excitations living on different spin bands for the latter case of strong spin polarization, however allows for quasiparticle coherence between spin bands when the band splitting is on the low-energy scale that constitutes the phase space for quasiparticles.

The self-energies introduced here are related to those we used in Eq. (2) by

 [^Δ+^Σ]=(ΣΔ~Δ~Σ). (9)

There is a further symmetry relating the Matsubara coherence functions with positive and negative frequencies:

 γ(−ϵn)=[~γ(ϵn)]†. (10)

As a result, we only need to consider in the following. All other quantities can then be obtained from the symmetry relations. The current density in the FM spin-bands () is obtained from

 →jη(→R,t)=4eNFηkBT∑ϵn>0⟨→vFη(→pFη)Re[gη]⟩η+, (11)

with the -component of and denotes a Fermi-surface average which is taken over momenta with positive projection on the -axis only. To obtain this expression, we did not only exploit the symmetry relations but also the fact that for the spin channels an additional symmetry between the diagonal components of holds: .

ii.2 Boundary conditions

The boundary conditions for Riccati amplitudes at interfaces and surfaces are formulated in terms of the normal state scattering matrix . (35); (23) The transport equations for Riccati amplitudes can be solved by integrating them along their characteristics, which are straight lines in real space(35). Every amplitude has a unique stable direction of integration. One may hence group them into “incoming” and “outgoing” amplitudes with respect to a scattering region, which will be the interfaces here. The convention is that outgoing amplitudes are denoted by capital case letters and incoming ones by small case letters. The transport channels that participate in scattering are labeled by . In our case, these will be band indices, as the conservation of parallel momentum that we assume excludes the scattering between different momenta in the same band. The boundary conditions then relate outgoing to incoming amplitudes. For numerical calculations, the most convenient way of solving these boundary conditions is to calculate(23):

 F=(1−~γ∘γ′)−1∘γ′,G=(1−γ′∘~γ)−1, (12)

where denotes a matrix product in channel space and:

 ~γkk′=δkk′~γk,γkk′=δkk′γk,γ′=S∘γ∘~S. (13)

is the normal state scattering matrix, which – thanks to the conservation of parallel momentum – is block-diagonal here, and the blocks can be labeled by . The outgoing Riccati amplitudes are then obtained from:

 Γk=G−1kkFkk. (14)

As explained at length in our earlier publications,(23); (24) the -matrices will be either or , depending on whether is larger or smaller than the Fermi wave vector of the minority band in the FM. Thus, all matrices defined above will have these dimensions, as channels with different do not mix.

Iii Solution of the Eilenberger equation

iii.1 Impurities

The simplest way to consider impurity scattering is to look at elastic, non-magnetic scattering in the self-consistent Born-approximation. This gives the following self-energies in the ferromagnet(21):

 (ΣΔ~Δ~Σ)η=12πτη⟨^gη⟩η. (15)

Note that here , denote all off-diagonal self-energies in Nambu-space, rather than a superconducting gap. The Riccati transport equations then need to be solved numerically, and in each step the impurity self-energy is updated.

The scattering time is related to the scattering length by . Since corresponds to the mean free path of the quasiparticles, it is directly related to the impurity concentration and should be the same for both spin-bands (we consider here spin-inactive scattering). The Fermi-velocities, however, are different and hence . This argument can also be made quantitative by noting that the impurity self energy in Born approximation is proportional to , with impurity concentration , impurity scattering potential , and density of states at the Fermi level . As this defines , one obtains for the two spin directions . As shown in appendix A, , and thus (assuming equal effective masses for the spin bands) , i.e. the mean free path for both spin bands is equal.

iii.2 Gap-equation

The mean-field gap-equation in the SC is

 Δ(z)⋅ln(TTc)=∑n(⟨f(ϵn,→kF,z)⟩−πΔ|ϵn|), (16)

where the coupling constant is eliminated in favor of . In appendix B we discuss in detail some sophisticated problems when iterating this equation at an interface. We show that whenever the quasiclassical Green function undergoes a jump-discontinuity under reflection from the interface, the only self-consistent solution of the above equation at the interface is . The discontinuity of the quasiclassical propagator in turn results from the fact that the microscopic reflection process at the interface falls outside the range of applicability of QC theory. To remedy this issue, we introduce a length cut-off and calculate the gap-equation only for distances larger than from the interface. We show that if the cut-off length is chosen small enough, neither its precise value nor the profile of the SC-gap for distances smaller than influences the results. We choose for our numerical calculations .

iii.3 Numerical solution

We found that a very stable way of solving the Riccati differential equations for spin-polarized systems is to resort to analytical solutions for constant self-energies. I.e., instead of using a standard solver for ordinary differential equations along each trajectory, we approximate the self-energies by a constant in each interval that corresponds to the grid-spacing of our discretization in the -direction, and use the analytical solution for homogeneous self-energies (23):

 γ(ρ)=γh+eiρΩ1[γ0−γh]{eiρΩ2+C(ρ)⋅[γ0−γh]}−1, (17)

where is the homogeneous bulk solution, is the initial value and parameterizes the respective trajectory as . Moreover, , and

 C(ρ)=C0eiρΩ1−eiρΩ2C0, (18)

where is given by the solution of:

 ~Δ=C0Ω1−Ω2C0. (19)

The value of at grid point is then calculated from the analytical solution with the value at grid point as initial value. These general formulas simplify inside the FM, because all quantities are scalar and commute. In the SC they simplify as well, since we only consider the order parameter self-energy there and thus always have and . The homogeneous bulk solution is for these cases

 γh=−ΔE±i√−Δ~Δ−E2, (20)

where .

As we include the self-consistent equations for the impurity potential and the order parameters, we have to deal with a non-linear problem, in which the solutions for different trajectories and Matsubara frequencies are coupled. We solve this problem iteratively as as illustrated in Fig. 2 (see also appendix B.2).

At each step, we assume a given set of initial values for the coherence functions at the outer boundaries of the system and at the interfaces and given values for the impurity potential and order parameters. We then integrate the transport equations and solve the boundary conditions at the interfaces to update all coherence functions and calculate the initial values for the next step of the iteration. Subsequently, all self-energies are updated by evaluating the respective self-consistent equation with the new coherence functions. The initial values at the outer boundaries are always assumed to be the respective SC bulk-solution. The order-parameter phase difference enters via these boundary values. At the initial step, we chose the coherence functions and the impurity potential inside the FM to be zero and the SC gap to be constant up to the interface.

Iv Scattering Matrix

The normal-state scattering matrices entering the quasiclassical boundary conditions cannot be obtained within QC theory itself. We calculate them from plane wave matching in the normal-state using a simple model for the interface scattering potential that captures all relevant effects.

 HI=∑→kμνc†→kμ(ξI→kδμν+→h⋅→σ)c→kν, (21)

where is the interface exchange field, and:

 ξI→k=ℏ2→k22m+VI−EF. (22)

plays the role of a spin-independent interface potential. Accordingly, the normal-state Hamiltonians of the ferromagnet and the superconductor read:

 HFM =∑→kμνc†→kμ(ξFM→kδμν+J⋅→σ)c→kν, (23) HSC =∑→kμνc†→kμ(ξSC→kδμν)c→kν. (24)

Obviously, is invariant under rotations in spin-space, while both and break spin-rotation symmetry. For we assume for simplicity a parabolic dispersion with the same effective mass in all cases. The dispersions and and are shifted by a constant energy and , respectively, compared to . The unitary matrix , which maps the spin-eigenbasis of the interface to that of the ferromagnet, reads:

 U=(cosα2−sinα2e−iφsinα2eiφcosα2). (25)

The scattering matrix is obtained by diagonalizing the Hamiltonians and inferring the respective values of at the Fermi-level for a given . The corresponding wave-functions in the three layers are then matched at the SC/I and the I/FM interface respectively, choosing a spin-quantization axis in the SC aligned with that of the interface, while at the I/FM interface, the rotation must be taken into account.

The requirement that must be conserved across the interface implies that even if the ferromagnet is not fully spin-polarized, some trajectories will have an evanescent solution in the minority band of the FM, corresponding to the half-metallic case. In the case were this minority band solution is propagating we obtain a -scattering matrix, if it is evanescent, the matrix is as discussed in appendix C.

V Results

Our scattering geometry is shown in in Fig. 3. In what follows, the Fermi-momentum of the majority (spin-up) band is assumed to be identical to the Fermi-momentum of the SC. This is simply a matter of convenience, as we must only distinguish two and not three types of trajectories in this case. The third type of trajectories would describe total reflection from the interface, which does not contribute to transport. Our methods can easily applied to the general case as well. The value of the minority band Fermi-momentum is determined by the exchange field , the magnitude of which is assumed to be the same in the interface and in the FM bulk (however not the direction).

In what follows we consider two parameter sets for the interface. A “high”-transparency interface with thickness and and a “low”-transparency interface with , and . If and are the sub-matrices of related to transmission from the SC to the spin-up band of the FM and vice versa, then the relevant transmission quantity for inducing triplet correlations in the spin-up band is:

 T↑=|S↑,SCiσySSC,↑|, (26)

and analogously we have for the spin-down band. These quantities are plotted in Fig. 4 for the two interfaces discussed here. Their functional dependence on was investigated extensively in Ref. (24).

In order to discuss various symmetry components of the induced triplet pair correlations, We use the following decomposition of the quasiclassical Green functions:

 g=g0+→σ⋅→g,f=(fs+→σ⋅→ft)iσy, (27)

and obtain the induced magnetization from

 →M=2μBkBTNF∑n⟨→g(ϵn,→pF)⟩. (28)

As measures of odd-frequency and odd-momentum triplet amplitudes we define

 →fϵ =kBT2∑n>0⟨→ft(ϵn,→pF)−→ft(−ϵn,→pF)⟩, (29) →fp =kBT2∑n[⟨→ft(ϵn,→pF)⟩+−⟨→ft(ϵn,→pF)⟩−]. (30)

Here, is the s-wave component of the odd-frequency correlations. does not correspond to a particular spherical harmonic. Our model system only breaks translational invariance in the -direction and is rotationally invariant in the --plane. We therefore project out the correlations which are odd under .

We first investigate the inverse proximity effect, as shown in Fig. 5. Spin-active scattering at the SC/FM interface induces triplet-pairing correlations and a magnetization inside the SC-electrodes close to the interface. Moreover, the order parameter is suppressed at the interface. While this is also caused by spin-active scattering, here the predominant cause for this suppression is the transparency of the interface. This can be seen in the topmost plots of Fig. 5 (left), where we compare a highly transparent interface to an interface in the tunneling limit. The reason why the spin-mixing effect plays a minor role here, is that under the assumption of a box shaped scattering potential the magnitude of the interface spin-mixing parameters that are responsible for creating triplet-pair correlations are typically small (see our extensive discussion of this effect in Ref. (24)). This can be seen in the top right panel, where the absolute value of the induced magnetization is shown. The transparency is largely irrelevant here, as the inverse proximity effect is due to spin-active back-scattering from the superconducting side, i.e. a property of the spin-dependent scattering phases in the reflection parts of the scattering matrix. The induced magnetization decays on the scale of the superconducting coherence length, . Note that this induced magnetization is a result of the superconducting inverse proximity effect (production of triplet pair correlations by the magnetic interface) and absent in the normal state.

In the lower panels of Fig. 5 we project our various symmetry components of the triplet pair correlations. Since the SC itself is spin-rotation invariant, it makes sense to plot the invariant quantity rather than the individual components of , which are not invariant. Since we consider in Fig. 5 a rather clean superconductor (), odd-momentum and odd-frequency correlations are comparable in magnitude and decay on the same length scale . In the case of the high transparency interface, the induced triplet-correlations and the magnetization show a slight dependence on the length of the junction (not shown), i.e. superconducting correlations that propagated through the FM-layer also contribute to the inverse proximity effect in that case.

We now study the influence of impurity scattering on the triplet pair correlations. In Fig. 6, we plot both the spatial profiles of the absolute value of the total momentum averaged equal-spin triplet correlations, (s-wave, odd frequency) and of the correlations which are odd under , (odd momentum), for each spin direction inside the FM-slab. As the impurity concentration increases, these correlations decay faster into the FM and the relative magnitude of decreases. However, the magnitude of the latter only becomes negligible for very high impurity concentrations. This shows that the Usadel approximation does require very short scattering lengths in order to be able to neglect odd-parity correlations in the ferromagnet. Right next to the interfaces, to suppress to about 1/10 of requires a scattering length as short as . One needs to carefully estimate if this is then already on the inter-atomic length scale, a length scale outside the scope of quasiclassical theory (both Eilenberger-Larkin-Ovchinnikov and Usadel theory), in which case the boundary layer must be treated microscopically. One option is to incorporate this region as an ’isotropization’ zone in effective boundary conditions for the Usadel equation.

We next discuss the Josephson effect in the structure shown in Fig. 1. With regard to the CPR, we first observe that higher harmonics, i.e. a non-sinusoidal CPR, can still be present in the high-transparency case, even if the suppression of pairing-correlations in the SC is properly taken into account. To show this, we consider a junction with high transparency interfaces in Fig. 7. The CPR is calculated as a function of the impurity concentration. Decreasing the scattering length leads to a suppression of the critical current, but also of higher harmonics contributions, as can be seen in the normalized plots on the right-hand side.

In Fig. 8, we consider a junction in the tunneling limit (, ) and calculate the critical current as a function of temperature. For a clean junction, , we find a strong spin-polarization of the current. For higher impurity concentrations, the current is not only suppressed globally, but the spin current, , also decreases relative to the total current. At high impurity concentrations, there may even be a crossover to a negative spin-current at small temperatures. This can be understood as follows. From the clean-limit solution, we see that the loss of particle-hole coherence due to temperature penalizes the spin-down current more than the spin-up current. The mean free path due to impurity scattering is, however, the same for both bands, and the cumulative effect of the scattering only depends on the length of the trajectory. This can easily be understood by rewriting the transport equation as:

 (iℏ^→vFη∇+i2ϵn/vFη)γ=[γ~Δγ+Σγ−γ~Σ−Δ]/vFη. (31)

Since all self-energies are proportional to (the superconducting gap is zero in the FM), the factor cancels on the right hand side of this equation, which describes impurity scattering, and only remains in – which is the term that can be interpreted as the suppression due to finite temperatures, and which clearly is seen to be spin-dependent.

Vi Conclusions

We carried out a comprehensive study of the triplet-Josephson effect and inverse proximity effect in SC/FM/SC junctions with highly spin-polarized FM-layers taking into account impurity scattering inside the FM in the self-consistent Born-approximation. We showed that odd-momentum triplet correlations are generically present in these junctions and are only suppressed at very high impurity concentrations. We found an induced magnetization in the SC-electrodes resulting from the inverse proximity effect and a suppression of the SC-gap which is mainly related to the transparency of the interface for the interface model we considered here. Larger spin-mixing phases would, however, also result in a substantial suppression of the gap close to the interface (21). These can be realized by considering a smooth scattering potential at the interface. Even if this suppression is taken into account, the CPR of high-transparency junctions may have higher-harmonic contributions. Regarding the critical current, we see that impurity scattering suppresses the Josephson spin-current and may even lead to an inversion of it at small temperatures and short scattering lengths.

Acknowledgements.
RG acknowledges financial support from the Karlsruhe House of Young Scientists. TL acknowledges support from the Swedish Research Council (VR). ME acknowledges support from EPSRC under grant reference EP/J010618/1.

Appendix A Expressions for the Fermi-surface average

For the spherical bands that we assume, the density of states at the Fermi-level is:

 NF,η =∫FS,ηd2pF,η(2πℏ)3|→vF,η(→pF,η)|=2m(2π)2ℏ3⋅pF,η. (32)

The FS-average for the self-energies is calculated from:

 ⟨x(→pF)⟩η=12NF,η∫pF,η0d[arcsinp||pF,η] p||(x+η+x−η), (33)

with and . And the FS-average for the current:

 ⟨vF,ηg⟩η,+=2π⋅^ez(2πℏ)3NF,η∫pF,η0dp|| p||⋅Re[g+]. (34)

The number of -points we use in our calculation (one hundred) is to small to make these expressions converge numerically. Thus, we use a numerical normalization factor which is calculated from Eq. (33) as . This ensures that all our FS-averages are properly normalized. To be precise, the numerical integration is implemented as a trapezoidal rule over one hundred -points which are equidistant in .

Appendix B Numerical iteration of gap and boundary conditions

b.1 Discretization of the gap-profile

The numerical solver we use is not a standard Runge-Kutta scheme but relies on the analytical solution for constant self-energies. The only numerical approximation made is assuming a step-wise variation of the self-energies between grid-points. Naturally, the question arises how this stepwise profile is to be interpolated when knowing the self-energies at a given number of grid-points. A subtlety here is that the coupled system of transport and self-consistency equations becomes fundamentally inconsistent if the condition , where denotes the th grid-point, is violated. This is because the asymptotic solution of the transport equation approaches the local homogeneous bulk-solution as the Matsubara frequency goes to infinity. As we will see in appendix B.3, this is crucial for formulating the gap-equation in a cut-off independent way. The interpolated profile we use is illustrated in Fig. 9 (left).

b.2 Iterative scheme

We assume that at the beginning of the iteration step, the incoming solutions at all interfaces and the self-energies are known from the previous step. The procedure works as follows:

1. Solve the boundary conditions at all interfaces.

2. Integrate all outgoing coherence functions up to the next interface or outer boundary.

3. Update all self-energies. We use the secant method to accelerate convergence.

4. Update all incoming coherence functions.

We use two abortion criteria, since the iteration in the SCs may converge a lot faster than in the central FM layer if the transparency of the interfaces is small. If the gap-profile is well-converged, the iteration in the SC is stopped and continues only inside the FM.

b.3 Elimination of the high-energy cut-off

Generally, technical energy cut-offs that appear in expressions for self-energies or observables in quasiclassical theory must be eliminated in favor of phenomenological parameters. This is usually achieved by adding and subtracting terms in Matsubara sums or energy integrals that cancel the divergent behavior of these expressions and allow to formally extend the sum or integral to infinity (29). The other half of such an “added zero” can typically be included in a phenomenological (measurable) parameter. While this procedure is standard for the weak-coupling gap equation in quasiclassical theory, it turns out to be problematic at interfaces. We are not aware of any discussion of this subtle issue in the literature before. We here discuss a method for dealing with this issue. We start from the self-consistency equation:

 Δ(→pF,→R)=kBT∑|ϵn|<ϵc⟨V(→pF,→p′F)f(→p′F,→R,ϵn)⟩. (35)

In what follows, we limit the discussion to singlet, s-wave superconductivity, i.e. , but the argumentation can easily be generalized. The interaction constant can be eliminated by linearizing the above equation close to (). Since the interaction constant is a bulk-property, it is justified to consider the homogeneous bulk case, i.e.:

 f(ϵn)=πΔ√ϵ2n+Δ2\lx@stackrelϵn≫Δ≈πΔ|ϵn|. (36)

This equation also shows why a technical cut-off must be included in (35), as the Matsubara sum diverges. This problem does not arise in the microscopic theory, where the asymptotic behavior of is given by .(36) The linearized equation reads:

 1V0=kBTcπ∑|ϵn|<ϵc1|ϵn|. (37)

Given that , we then obtain the well-known standard result in terms of the digamma-function ,

 1V0=ψ(12+ϵc2πkBTc)−ψ(12). (38)

At this point, one can see that keeping a technical cut-off is, at least in the Matsubara case, not tenable. should of course be an integer number. Obviously one can choose so that is integer here, but must be the same for all temperatures and since the spacing of Matsubara frequencies is temperature dependent, it will generally not be commensurate with . One can verify numerically that if such a cut-off is kept in , the gap-profile one finds is the BCS-relation with a jigsaw-pattern superimposed to it, which vanishes as . The standard way of dealing with this problem is to regularize the Matsubara sum. This is achieved by subtracting the asymptotic form of on both sides of (35):

 Δ(1V0−kBTπ∑|ϵn|<ϵc1|ϵn|) (39) =kBT∑|ϵn|<ϵc(⟨f(→p′F,→R,ϵn)⟩−πΔ|ϵn|).

Both sides of this equation can now be regularized, if the expression (38) for is taken into account. The cut-off can be taken to infinity, which yields

 Δln(TTc)=kBT∑ϵn(⟨f(→p′F,→R,ϵn)⟩−πΔ|ϵn|). (40)

For the boundary value problem one may object that the counter term was adapted to the homogeneous bulk-case, but expanding the solution of the Riccati-equations in orders of , one finds (33):

 γ(→pF,→R,ϵn) =−Δ(→R)2iϵn−iℏ(→vF∇)Δ(→R)(2ϵn)2+O(1/ϵ3n), (41) ~γ(→pF,→R,ϵn) =~Δ(→R)2iϵn−iℏ(→vF∇)~Δ(→R)(2ϵn)2+O(1/ϵ3n),

and therefore:

 f =−2iπ∑n(γ~γ)nγ (42) =πΔ(→R)|ϵn|∓2πℏ(→vF∇)Δ(→R)(2ϵn)2+O(1/ϵ3n), ϵn≷0.

Thus, even if the gap-profile is inhomogeneous, the asymptotic behavior remains the same as long as the gap is a smoothly varying function.

Up to this point, the procedure is well-known and established. A problem does, however, arise at interfaces where the boundary conditions leads to a discontinuous behavior of the coherence functions under reflection. To see this, we consider the coherence function to leading order in . The incoming coherence function at the interface is then given by . To illustrate the point we here assume that the solution of the boundary conditions at the interface is:

 Γ=RSCγ(z=0)~RSC, (43)

where denotes the reflection matrix on the SC-side of the interface. This is not fully generic, since we assume that the transmission from the other side of the interface can be neglected. However, for the case we consider here this is reasonable as the coherence functions on the FM side are exponentially suppressed by propagation through the FM layer and can thus be neglected at high energies. Now it is easy to see that the asymptotic behavior of at the interface is given by:

 ⟨f(→pF,ϵn)⟩ =⟨f(→pF,ϵn)⟩++f(→pF,ϵn)⟩− =π2Δ|ϵn|(1−iσy⟨RSCiσy~RSC⟩). (44)

Here denotes an FS-average over trajectories pointing towards/away from the interface. It is obvious that this -function will only have the required asymptotic behavior if or . This implies that the right-hand side of equation (40) will diverge if that is not the case. Adapting the cancellation term will not help either, since the coupling constant is fixed by the bulk-behavior, which implies that the left-hand side of (40) would diverge in that case. So eventually, one arrives at the conclusion that the only consistent solution for is . The problem here is the discontinuous jump of the -function upon reflection at the interface, which is directly related to the fact that the interface region cannot be described within quasiclassical theory. While the use of effective boundary conditions is sufficient to circumvent this problem if the self-consistency relation is neglected, one must go a step further if it is retained. Our solution to this problem is to assume that the self-consistency equation must not be calculated if the distance to the interface is smaller than a length cut-off . To see how this resolves our issue, we have to calculate and . Here, is the incoming coherence function and obtained by integrating the Riccati equation from some initial value inside the bulk. Since the gap-profile is smooth for , this implies that it does have the correct asymptotic behavior at high energies. originates from the interface. If we assume that the gap has a constant value in some environment of , where can be arbitrarily small, the solution for will be given by (23)

 Γc=γh,c+eiηΩ1[γ0−γh,c]{eiηΩ2+C(ρ)⋅[γ0−γh,c]}−1, (45)

where and . This solution has the correct asymptotic behavior, irrespective of how the gap-profile varies for , as it is given by the homogeneous solution at sufficiently high energies. We found that the sum in (40) must be, however, calculated up to very high energies, since should be very short. Thus, we need to discuss high-energy contributions next.

b.4 High-energy contribution

To numerically compute (40), we consider two numerical energy cut-offs on the right-hand side:

 ln(TTc)Δ =∑|ϵn|<ϵc,low(⟨f(→p′F,→R,ϵn)⟩−πΔ|ϵn|) (46) +∑ϵc,high>|ϵn|≥ϵc,low(⟨f(→p′F,→R,ϵn)⟩−πΔ|ϵn|).

Up to we calculate the full coherence functions in the whole structure without any approximations. The high-energy contribution is only calculated in the SC-electrodes. To efficiently calculate this contribution, we consider the linearized Riccati-equation:

 i∂ργ(ρ)+2iϵnγ=−Δ(ρ), (47)

which is easily solved exactly by a variation of constants:

 γ(ρ) =[C(ρ)+γ0]β(ρ), β=exp[−2ϵnρ], C(ρ) =∫ρ0dρ′iΔ(ρ′)β(ρ′). (48)

Again assuming that the gap-profile is constant, this yields:

 γ(ρ)=γh+[γ0−γh]β(ρ). (49)

We also assume that the boundary conditions reduce to (43). It is easy to convince oneself that the linearized Riccati-equation is exact up to second order in . Up to second order, we also have:

 f=−2iπγ, (50)

from which we calculate the contribution to the gap function.

As discussed above, the need for taking into account this high-energy correction factor arises from interface scattering. It is therefore not surprising that it vanishes rapidly inside the SC-bulk. This fortunately implies that we must only calculate it very close to the interface.

b.5 Extension to zero temperature

For , we have

 ϵn+1−ϵn=2πkBT→0, (51) kBT∑ϵn→12π∫∞−∞dϵ. (52)

The zero-temperature Green function is thus obtained by replacing and all Matsubara sums must be replaced with the above integral. Care must be taken in regard to the normalization factor in the self-consistency equation:

 ln(TTc)+nc∑n=11n−1/2, (53)

with , since both terms in this sum become infinite. For , the second term approaches , and we have:

 ln(TTc)+nc∑n=11n−1/2=ln(ϵc2πkBTc)−ψ(12). (54)

Appendix C Calculation of the S-matrix

c.1 HM-case

First, we need the wave-vectors which are obtained from solving the one-dimensional Schrödinger-equation for given , . These are and for , for , and and for , where is associated to the evanescent mode. These wave-vectors may be real or imaginary, corresponding to either a conducting or an insulating interface. Moreover, we have the Fermi-velocities and in the SC and the majority band of the FM, respectively, which are proportional to and respectively. After defining the auxiliary quantities

 a± =cosk±a+ipk±sink±a, (55) b± =cosk±a−ipk±sink±a, (56) c± =ik±k↑sink±a+pk↑cosk±a, (57) d± =ik±k↑sink±a−pk↑cosk±a, (58)

where is the thickness of the interface layer, and the matrices:

 P(x) =U†(x+00x−)U, (59)

where is the spin-rotation from the spin-basis of the interface to that of the FM-bulk, we obtain:

 SR=U−1V (60)

with

 U =⎛⎜ ⎜ ⎜ ⎜ ⎜⎝−√v↑vP11(b)−√v↑vP12(b)1−√v↑vP11(d)−√v↑vP12(d)1P21(d)−iκk↑P21(d)P22(d)−iκk↑P22(b)0⎞⎟ ⎟ ⎟ ⎟ ⎟⎠, (61) V =⎛⎜ ⎜ ⎜ ⎜ ⎜⎝√v↑vP11(a)√v↑vP12(a)−1√v↑vP11(c)√v↑vP12(c)1iκk↑P21(a)−P21(c)iκk↑P22(a)−P22(c)0⎞⎟ ⎟ ⎟ ⎟ ⎟⎠. (62)

c.2 FM case

In the FM case, the minority band of the FM has a propagating mode. We write and for the wave-vectors of the plane-waves in the FM-bands and , for the corresponding group velocities. We use the same definitions for , and as above, but redefine:

 c±=ik±sink±a+pcosk±a, (63) d±=ik±sink±a