Quasiclassical circuit-theory of contiguous disordered multiband superconductors

We consider a general problem of a Josephson contact between two multiband superconductors with coexisting superconducting and magnetic phases. As a particular example, we use the quasiclassical theory of superconductivity to study the properties of a Josephson contact between two disordered $s^{\pm}$-wave superconductors allowing for the coexistence between superconductivity and spin-density-wave orders. The intra- and inter-band scattering effects of disorder are treated within the self-consistent Born approximation. We calculate the spatial profile of the corresponding order parameters on both sides of the interface assuming that the interface has finite reflection coefficient and use our results to evaluate the local density of states at the interface as well as critical supercurrent through the junction as a function of phase or applied voltage. Our methods are particularly well suited for describing spatially inhomogeneous states of iron-based superconductors where controlled structural disorder can be created by an electron irradiation. We reveal the connection between our theory and the circuit-theory of Andreev reflection and extend it to superconducting junctions of arbitrary nature. Lastly, we outline directions for further developments in the context of proximity circuits of correlated electron systems.


I. INTRODUCTION
In many practical cases superconductivity occurs in the form of a spatially inhomogeneous state [1,2]. This can be triggered intrinsically due to thermodynamic reasons or created extrinsically by forming contacts of superconductors with other materials. The fundamental example of the first kind of inhomogeneity is given by the Abrikosov vortex state, which brings about the spatial modulation of the order parameter [1,3]. Josephson junction is the example of the other kind, where inhomogeneity is created near the contact area when two superconductors are brought into proximity via a tunnel barrier or other type of the weak link [4,5]. In both of these cases, and many other physical situations, the spatial inhomogeneity extends over the length scale of superconducting coherence length that is large as compared to electron Fermi wavelength. Under this condition, the semiclassical theory of superconductivity based on Eilenberger [6] and Usadel [7] equations become applicable. These two methods were developed to treat relatively clean and strongly disordered superconductors respectively. The solutions of the Eilenberger and Usadel equations relate the observable properties of a superconducting structure, for example critical current, to microscopic characteristics of the materials forming the junction and its geometry.
An alternative method to describe mesoscopic superconducting structures is based on the random matrix and scattering matrix theories [8]. In this approach, all microscopic details are condensed into symmetry properties of the scattering matrix representing a disordered region of the junction which is typically parametrized by a set of transmission eigenvalues. An observable of interest is then expressed in terms of these transmissions similar to the Landauer-Buttiker transport theory. This phenomenology is more straightforward and intuitive than semiclassical kinetic theory, but it is more restrictive in terms of conditions when it applies. Yet there is a parameter range when both methods work, however, the connection between them is not immediately obvious.
This link has been provided by the circuit-theory of Andreev reflection developed originally by Nazarov [9], later reviewed and extended by several authors [8,10,11]. Circuit theory can be formulated as the finite set of rules for connectors and nodes of a given superconducting devices, analogous in spirit to Kirchhoff's rules. It also gives a prescription to deal with boundary conditions and in particular average over the transmission eigenvalues [12,13], which are in general random for a disordered or chaotic junction between superconductors. In recent years, we witnessed the emergence of novel classes of multiband unconventional superconductors, primarily the large family of iron-pnictides (see reviews [14,15] and references therein). Semiclassical methods of superconductivity were successfully applied to describe their properties including the proximity and Josephson effects [16][17][18][19][20][21][22][23][24][25][26][27][28][29][30] but circuit-theory has not been derived for these systems. The motivation for this work is to put forward a detailed theory of superconducting contacts, where the material constituents forming the junction harbor complex superconducting phases and competing order parameters.
This paper is organized as follows. In Section II, we formulate the simplest two-band model that allows for the coexistence of superconducting (SC) and spindensity-wave (SDW) orders, and derive the Eilenberger equations, which form the technical basis for our work. arXiv:1908.07953v2 [cond-mat.str-el] 28 Dec 2019 In the Section III, we employ the method developed by Yip [31] to solve the quasiclassical equations and, at the same time, satisfy full nonlinear boundary conditions derived by Zaitsev [32]. In addition, we have arbitrary transparencies and shapes of potential barriers forming the constriction. We demonstrate, that the special auxiliary decomposition of nonlinear constraints naturally leads to the circuit-theory boundary conditions as elaborated by Nazarov [10]. In Section IV, we present the results for the numerical solution for the spatially dependence of the superconducting order parameter, local density of states at the interface and Josephson current. In Section V, we briefly review several universal examples of the Josephson effect in mesoscopic superconductornormal-superconductor (SNS) structures and related devices with insulating barriers and micro-constrictions. We discuss how circuit-theory captures in a unified fashion multiple results for the Josephson current-phase relationships that were previously known from the separate semiclassical calculations and extend that to the case of proximity junctions of correlated electrons. Section VI is devoted to the discussion of the results and outlook for further developments. Lastly, in Appendices A, B & C, we provide the details on the derivations of the expressions that we used to obtain the solution of the Eilenberger equations.

II. FORMULATION OF THE PROBLEM
In what follows, we introduce the model Hamiltonian and write down the quasiclassical equations for the correlation functions which are used to determine the spatial profile of the superconducting and magnetic order parameters across the interface. We will consider short junctions at an arbitrary transparency between two multiband superconductors. Generally, disorder in these systems induces both intra-band and inter-band scattering. In addition, the symmetry of the order parameter can be unconventional. As a guiding example, we will study a superconductor with the s ± symmetry of the order parameter relevant for some classes of iron-pnictides. We will also consider a more complicated case, when the superconducting state coexists with another order such as spin-density wave. Under such circumstances, it will be impossible to write the Josephson current of such junctions just in terms of transmission eigenvalues. However, it is still possible to derive a closely related circuit-theory expression written in the form of semiclassical Green functions.

A. Model Hamiltonian
Following the discussion in Refs. [33,34], we consider a model with two cylindrical Fermi surfaces. One Fermi surface has an electron-type (c) and the other one has a hole-type (f) excitations. We introduce the following

eight-component spinor in momentum representation
Hereσ y is a Pauli matrix,ψ † pa = (a † p↑ , a † p↓ ), (a = c, f ) andψ T denotes the transpose of the operator. The form of (1) ensures the correction definition of the spin density operator at point r: In this paper, however, we will limit our discussion to the z-component of the spin operator (2), and as it turns out it will be more convenient to work with the following spinorΨ † The Hamiltonian for our problem can be written down using the mean-field approximation, The noninteracting partĤ 0 has the standard form pertinent to our choice of the basis spinor (3): where µ is a chemical potential and the mass anisotropy between hole-and electron-like bands was ignored. The remaining mean-field part contains two terms which account for the superconducting pairing in s ± symmetry channel with the amplitude ∆ and spin-density wave order parameter, M = M e z : In these expressions, we use the Pauli matrices τ i , ρ i and σ i (i = 1, 2, 3) defined in the subspace of band, Nambu and spin degrees of freedom correspondingly. Lastly, the third term in (4) describes the effects of disorder induced by chemical substitution at lattice sites The first term in the brackets accounts for the intraband scattering, while the second term describes the scattering between the bands. Having defined the Hamiltonian, next we outline the steps which lead to the equations for the quasiclassical Green's functions.

B. Eilenberger equation
For simplicity, let us first assume that disorder is the only source of spatial inhomogeneities. To derive the equations for the quasi-classical correlation function, one starts with the Dyson equation for the single-particle Green's function in the imaginary time representation averaged over various disorder realizations: with ω n = πT (2n+1) being the fermionic Matsubara frequency. Within the self-consistent Born approximation, the corresponding expression for the self-energy readŝ where Γ 0,π ∝ ν F |u 0,π | 2 are the corresponding disorder scattering rates and ν F is the density of states at the Fermi level.
The quasiclassical Eilenberger function is defined according tô The equation for the function can be obtained from the Dyson equation above (9) by eliminating the single particle spectrum, ξ p . In the spatially inhomogeneous case, which naturally arises in nonzero external magnetic field or in the presence of a contact between two superconductors, functionŝ G(p, iω n ),Ĝ(iω n ) as well as self-energy Σ(iω n ) and the order parameters ∆, M will also depend on the 'centerof-mass' coordinate R = (r + r )/2. Thus, functionĜ (R, iω n , v) satisfies the following equation: where [X;Ŷ ] stands for the commutator of matrices and v F = v F n. In equations above, we have omitted writing the dependence on R and ω n in relevant functions for brevity. Importantly, since the quasiclassical equations are linear inĜ one needs to specify the constraint condition for this function to avoid an ambiguity. Simple calculation shows that the quasiclassical function must satisfy the nonlinear normalization condition In addition, given the problem at hand, the Eilenberger equation above needs to be supplemented with the boundary conditions.

C. Boundary conditions
To determine the spacial variation of the order parameters ∆(x) and M (x) through the interface, we need to solve (12) on each side of the interface and then match the quasiclassical functions at the interface (x = 0) with the use of the boundary conditions [32]: Here the dependence of the quasiclassical functions on Matsubara frequency has been suppressed,Ĝ s(a) (0) = (Ĝ(v x , 0) ±Ĝ(−v x , 0))/2 andĜ s± (0) = (Ĝ s (+0) ± G s (−0))/2. New parameter D is the transparency coefficient for the interface, while R = 1 − D. Note that the boundary conditions are non-linear. This happens because the interference between the quasiparticle paths on both sides of an interfaces has been completely ignored. The effects of the interference between the trajectories go beyond the scope of this work and will be reported elsewhere [35].

III. ANALYSIS OF THE QUASICLASSICAL EQUATIONS
In this Section, we first discuss the approach to solving the Eilenberger equations and then show the results of our numerical solution for the spatial variation of ∆(x) and M (x) across the junction.

A. Ansatz for the quasiclassical functions
The task of solving matrix Eilenberger equation (12) presents a major challenge. In addition to being supplemented by the nonlinear boundary conditions (14) for a given ∆(x) and M (x), these functions must be determined self-consistently via relations [36]: where Ĝ = Ĝ (x, ω n , v x ) and averaging is performed over the directions of the quasiparticle trajectories. Here, Ω Λ is the energy scale of an ultraviolet cutoff, while λ sc,sdw are the corresponding coupling constants and we employed the standard notation a + = a x + ia y . Clearly, to make further progress we need to specify the matrix structure of the functionĜ that will respect the nonlinear normalization constraint. We now present the ansatz for the functionĜ, that follows the constraint: (16) where ζ = (x, ω n , v F n x ). The first term in (16) is the quasiclassical function for the normal component, so that in the absence of an interface and when ∆ = M = 0 it obtains g z (iω n ) = sign(ω n ). The second term accounts for the superconducting correlations: Here the anomalous f z component must be constant in the bulk, while f x is only nonzero in a close proximity to an interface and is an odd function of unit vector n x . SimilarlyĜ Finally, the last term in (16), as will be shown below, appears only when both ∆ = 0 and M = 0: and it vanishes in the bulk on the both sides of an interface.
After substituting these expressions into (12) and equating the terms proportional to the same combination of the direct matrix productsτ iρjσk we find the following equations: To keep concise notations we have introduced an additional self-energy functions with Γ m = Γ 0 − Γ π and Γ t = Γ 0 + Γ π , and with implicit averaging that is performed over all possible values of unit vector n x .

B. Quasiclassical function components in the bulk
Eilenberger equation acquires the simplest form in the bulk when the gradient term can be discarded. For definiteness we consider the bulk of a superconductor at x > 0. According to our discussion above, only three functions g z , f z and s z are non-zero. Simple calculation yields where superscript b in all the functions implies value of that function taken in the bulk of a sample, namely g b z = g z (x → ∞) etc. Furthermore, it is easy to show that the third equation is redundant. However, as we will see below, in the vicinity to the interface an analogue of this equation will determine the spatial variation of the function g x , Eq. (19).
Numerical solution of the first two equations together with the self-consistency equations (15), produces the well known phase diagram of SC-SDW coexistence shown in Fig. 1 for a certain choice of parameters (compare that to Refs. [33,34,36]). This model reveals the dome-like structure of superconductivity overlapping with SDW state. Bending of the superconducting dome in the nonmagnetic phase occurs due to finite Γ π that serves as an effective pair-breaking factor for s ± superconductivity. Suppression of magnetic order already occurs at the level of intra-band scattering that is governed by Γ 0 . The width of the coexistence region max[T N (∆) − T c (M )] can be controlled by the ratio between the scattering rates Γ π /Γ 0 and can change substantially, however within this model it always remains rather narrow.

C. Normalization condition
Simple algebraic manipulations with equations (20) show that components ofĜ satisfy for any value of x. By sending x → ±∞, it immediately follows that the constant must be equal to one. On the other hand we can use Eq. (13) directly with (16) to find where we already took into account (22). Again, as it can be checked by the direct calculation the second term here is actually a constant Constant appearing in this equation must be zero due to the vanishing of functions g x , s x and f x in the bulk. Thus, we have derived the matrix form of the quasiclassical function and have demonstrated that normalization condition for the functionĜ holds.

IV. APPLICATIONS AND RESULTS
In this Section, we present the results of our analysis of the quasiclassical equations and use these results to compute the observables: local density of states and critical current.

A. Order parameters and local density of states
The numerical solution of the quasiclassical equations (20) for the geometry of a junction illustrated in Fig. 2, is plotted in Fig. 3. These plots are one of the main results of this paper concerning the nature of the proximity effect in a complex superconducting phase. One important observation that we can make in regards to the spacial changes of the superconducting order parameter is that it varies substantially only in an immediate vicinity of the boundary between the two superconductors. On the contrary, the spin-density-wave order parameter changes on a somewhat larger length scale. Therefore, our results formally justify the often used approximation of constant order parameter on both sides of an interface.
It is important to point out here that precisely this aspect of the problem leads to practically universal predictions for the current-phase relations of S SDW -N-S SDW Josephson junctions. It should be noted, however, that disorder model considered here is not the only one that gives coexistence scenario. In the band models [37,38], coexistence region can be significantly broader in parameter space so that proximity problem in principle may also have qualitatively different behavior, in particular displaying longer coherence lengths.
With the solution of the Eilenberger equations, we can easily determine the local density of states at the interface (x = +0), using the well known expression [39,40] ρ LDOS (ω, x) = Tr Re τ 3ρ3σ0Ĝ (ω + i0, x, n x ) and the averaging is performed over all directions of unit vector n. For a fixed position from an interface, order parameter ∆(x = +0) and components of G(iω n , x = +0, n x ) are known from the numerical solution. Upon averaging over n x , we obtain Ĝ (iω n , x, n x ) and by performing an analytic continuation to real frequencies, iω n → ω + i0, we then, are able to compute the

B. Josephson effects
In analogy with the local density of states, Josephson current through the junctions also admits representation in terms of the Eilenberger function: In order to computeĜ a (0), one generally speaking needs to consider equations (12) with the complex ∆(x) on each side of the interface. However, in the case when there is no magnetic field the problem can be significantly simplified by using the unitary transformation, which eliminates the phase from the order parameter. For example, assume that the order parameter on the left side of the interface is ∆ 1 (x) = |∆ 1 (x)|e iχ(x) . Then, we introduce a unitary transformation S(χ) = cos(χ/2)τ 3ρ0σ0 − i sin(χ/2)τ 3ρ3σ0 .
It is easy to verify thatŜ∆ 1Ŝ † = |∆ 1 |. If one now implements this unitary transformation for the quasiclassical functions, it follows that the Eilenberger equation acquires essentially the same form as the one with purely real order parameter, except for the extra term proportional to ∂χ/∂x. This term, however, can be ig-nored for one does expect the phase to vary substantially across the junction. Furthermore, by performing the inverse unitary transformation, one can determine function G a (iω n , 0, v F sin φ).
After somewhat lengthy calculation, we found Here χ is the global phase difference between the order parameters on the both sides of the interface and the superscripts l and r mean that the functions should be evaluated either on the left (x = −δ) or the right (x = +δ) sides of the interface. Our results for the Josephson current-phase relation are shown in Fig. 6. Quite naturally, we find that for small D, the Josephson current will be proportional to sin χ. We note that the quasiclassical functions which account for the magnetic order do not explicitly enter into the expression for the Josephson current. Motivated by ideas and practical implementation of Josephson Scanning Tunneling Spectroscopy (JSTS) as a diagnostic of unconventional superconductivity [41][42][43][44][45][46], we briefly consider this effect in our model. To determine the dependence of the critic current on external voltage V we use the usual expression [41,47]: Here F 1(2) (ω) are the quasiclassical anomalous Green's functions for the first (second) superconductor and R N is the contact resistance. In the JSTS setup, one usually uses SC tip with known properties as a reference point and another SC as a study system. For this reason, we choose F 1 (ω) in the form which describes a clean BCS superconductor with the pairing amplitude ∆ BCS : As for the function F 2 (ω 2 ), it can be directly obtained from f b z (iω n ) by performing an analytic continuation. For simplicity we consider f b z (iω n ) calculated away from the contact that creates spatial inhomogeneity, but obviously the calculation can be done for any spatial location of the tunneling tip with respect to the junction.
The results of our numerical calculations based on Eq. (28), are presented on Fig. 7. To make a contrast with the textbook example, we also plotted the critical current for a contact between two BCS superconductors. It is worth reminding the reader that in the case of two identical BCS superconductors, the real part of the critical current has a logarithmic divergence at the threshold voltage eV = 2∆ BCS , while the imaginary part has a square-root singularity at the same value of external voltage. Now, if one of the BCS superconductors is replaced with an unconventional disordered superconductor, the presence of disorder and nonzero magnetization produces the smearing of the above mentioned singularities and lead to the substantial broadening of the dependence I c (V ).

V. OVERVIEW OF UNIVERSAL JOSEPHSON CURRENT-PHASE RELATIONS
The physics of the dc-Josephson effect is ultimately related to the sub-gap states that carry the supercurrent.  These states form as a result of Andreev reflections that electrons undergo when impinging on superconducting interfaces. Location of these states inside the superconducting gap depends on the superconducting phase difference across the junction. Kulik solved the first microscopic model of superconductor-normal-superconductor (SNS) model and derived the spectrum of Andreev states in various limits (see e.g. Refs. [5,47]). In a way, this work marked the beginning of intensive studies of various kinds of Josephson junctions that spanned over multiple decades. The interest to this problem has been continuously sustained to the present day not only due to the fundamental physics involved and applications of this effect, but also emergence of the new classes of unconventional superconductors. The most elegant way to derive the spectrum of Andreev bound states is by using scattering matrix approach. Beenakker derived the general determinant for-mula [48], which has very transparent and intuitive meaning. In the limit of the short junction, when length of the link separating superconducting leads is small compared to the coherence length (L ξ), this formula simplifies to a famous expression for a pair of Andreev levels per-channel of the junction E n (χ) = ±∆ 1 − D n sin 2 (χ/2). (30) In this theory, the junction is modeled as a multi-mode conductor where each conduction channel labeled by an index n has certain transmission coefficient D n . The Josephson current J(χ) carried by these states is given by With this formula at hand, one can recover multiple special cases. Indeed, at temperatures close to the critical, T T c , the superconducting gap is small, ∆ T , so that one can expand the thermal factor of hyperbolic tangent at small argument, which gives as a result where we introduced the total normal state resistance of the junction by means of the Landauer formula R −1 N = (2e 2 /h) n D n . This sinusoidal current-phase relationship for a superconductor-constriction-superconductor (ScS) junction was originally derived by Aslamazov and Larkin from the Ginzburg-Landau theory [49]. This result is universal in the sense that it applies to any kind of constriction. It is also a generic property that Josephson current is harmonic (sinusoidal) near T c . Alternatively, one can consider a tunnel barrier, D n 1, which corresponds to a class of superconductor-insulatorsuperconductor (SIS) type junctions. This yields the following expression for the current in the form that was derived first by Ambegaokar and Baratoff from the tunneling Hamiltonian [50]. In the opposite limit of fully transparent channels, D n = 1 for n = 1, . . . , N , one recovers the model of quantum point contact, namely S-QPC-S junction. In this case, the spectrum of Andreev levels simplifies to E n = ∆ cos(χ/2) for any channel and the corresponding current is This formula was derived first by Kulik and Omelyanchuk from the Eilenberger equations [51].
In realistic contacts of actual devices, transmissions are neither fully ballistic nor of tunneling type, rather there is a continuous distribution ρ(D) of transmission eigenvalues. There are several generic contact types that have been discussed in the literature. Their distributions are described by the function of the form The case with the power exponent p = 1/2 corresponds to two ballistic connectors with equal conductances in series. The case with p = 1 corresponds to the Dorokhov function of a diffusive connector [12]. The case with p = 3/2 was considered by Schep and Bauer [13] and corresponds to an interface with a high density of randomly distributed scatterers. The normalization factors N p are chosen in such a way as to ensure the total resistance of the junction to be R −1 Consider p = 1 first: the normalization is N 1 = 1/2, and the ensemble averaged Josephson current is (taken at zero temperature for simplicity) The integral can be found in elementary functions with the final result that corresponds to Kulik-Omelyanchuk computation carried out for the disordered SNS junction based on the Usadel equations [52]. For the case with p = 3/2, an analogous averaging yields which after the substitution D = sin 2 (x) reduces to the complete elliptic integral of the second kind This current-phase relationship was obtained first by Kupriyanov and Lukichev from Usadel equations in the context of SINIS junction (see Ref. 53 for review). Its connection to Eq. (31) with subsequent averaging over ρ(D) was pointed out by Brinkman and Golubov [54].
Lastly, the case with p = 1/2 corresponds to a chaotic cavity/quantum dot that supports current as was studied by Brouwer and Beenakker [55]. The corresponding current-phase relationship can be written as a combination of elliptic functions of the first and second kind All these examples give different functional form of the Josephson current, yet all of the them support parametrically the same critical current, J c ∆/eR N , which is governed by the total conductance of the junction in the normal state and size of the gap in the leads. In that regard these results are universal. In the extended junctions, when the size of the weak link is large as compared to the coherence length L ξ, the situation is different as critical current will decay as a power law or even exponentially with L depending on temperature. The decay of the current is related to the large dwell time needed for quasiparticles to travel across the junction. Additional features may appear due to the complexities of the proximity effect related to induced spectral gaps, including secondary gaps near ∆, that ultimately modify current amplitude and its dependence on the phase [56].
Interestingly, the family of such (almost universal) results for mesoscopic systems can be extended to include more complex proximity junctions of correlated electrons, such as S SDW NS SDW or S SDW INIS SDW . Indeed, thanks to the exact numerical results we have for the spatial profile of the order parameters, one can take the step function model to the leading approximation. Assuming the symmetric case one then finds for the trace in Eq.
For a multi-mode junction without inter-mode scattering, one can directly average this expression over ρ(D) which gives for the Josephson current a compact formula where we introduced dimensionless function At zero temperature, the Matsubara sum can be converted into an integral over the real axis of continuous frequencies 2πT ωn → dω and remaining calculations can be carried out for any of the discussed models of transmission distributions (35). For instance, for the S SDW INIS SDW junction (model with p = 3/2) one finds where c p is the numerical factor of the order of one. We notice that in the part of the phase diagram where SDW competes with SC, M ∆, supercurrent is suppressed J c ∼ ∆ 2 /(eR N M ). Other models of contacts can be analyzed in the similar way.

VI. DISCUSSION AND OUTLOOK
By using the quasiclassical theory of superconductivity, we have performed the fully self-consistent treatment of the Josephson junctions between two two-band superconductors in which nodeless order parameter with s ±symmetry, coexists with an itinerant SDW order. By solving the corresponding quasiclassical equations for the Eilenberger functions, we have found the variation of the superconducting and magnetic order parameters across the interface with arbitrary transparency. Using the results of the numerical solution, we have computed the local density of states ρ LDOS , Josephson current-phase relations J(χ), as well as dependence of the nonequilibrium critical current on external voltage, I c (V ). The features pertaining to the presence of the magnetic order, are clearly pronounced in the local density of states, on distances of the order of the coherence length from the interface provided that pairing amplitude in at least one of the superconductors exceeds the magnetic order parameter. For the critical current, we find that (i) suppression in the parts of the phase diagram where SDW dominates superconductivity, and (ii) disorder leads to smearing of the various sharp features in I c (V ) found for the contact between BCS superconductors.
Our work can be further extended in multiple directions. It is of practical importance to consider effects of disorder for more realistic Fermi surfaces including ellipticity for example. It is of special interest to consider three-band models that is the minimal setting for the appearance of nematic order that has to be included in the Eilenberger semiclassical scheme. There is also clear motivation to extend this semiclassical theory to real time axis to address dynamical responses of superconductors with competing orders.  Auxiliary quasiclassical functionsĜ >,< (u) are considered to be the functions of parameter u = x/vxτ∆, where τ∆ is the relaxation time of the order parameter. Top panel: the trial order parameter is chosen to correspond to the physical order parameter of a superconductor to the left of the interface and the particle's velocity is assumed to be negative, vx < 0. Therefore, the diverging solution of the quasiclassical equations at x → ∞ (u → −∞) is denoted byĜ > 1 (u), while the diverging solution at x → −∞ (u → ∞) is denoted byĜ < 1 (u). Bottom panel:Ĝ > (u) andĜ < (u) now diverge on opposite sides of the interface and trial form of the order parameter is chosen to match the bulk value of the order parameter on the right hand side of the interface.

Appendix A: Method of an auxiliary solution
In what follows, we will briefly review a theoretical approach, first proposed by Yip [31], that allows one to circumvent the issues associated with the nonlinearity of the boundary conditions. This approach makes it possible to write down the expressions for the quasiclassical functions, which match each other on the interface with the finite reflection coefficient (see also Refs. [57][58][59][60][61]). The only assumption which goes into making this procedure work is that superconductors extended to distances on which the correlations functions and the corresponding order parameters recover their bulk values. The trick we will use, consists of expressing the physical solution in terms of the unphysical (i.e. divergent) ones and implementing the resulting relations to simplify the boundary conditions.
Let us make some general observations. First, it is clear that the Eilenberger equation (4) can be formally written as v x ∂ xĜ = L ;Ĝ . (A1) IfĜ sol is a solution of (A1), thanĜ 2 sol is also a solution. Furthermore, due to the normalization condition, it follows that once the boundary conditions are taken into account, our quasiclassical function matrix G at the both sides of the interface, followsĜ 2 =Ĝ 2 0 , where G 0 is a quasiparticle correlator in the bulk,Ĝ 2 0 =Î. Instead of solving a problem on both sides of the interface and then trying to match the corresponding solutions, one ignores the interface and solves independently two problems with the parameters matching those at each side of the interface. Solution of each of these two problems, requires a profile of the pairing amplitude as an input, so for the both problems one can use the bulk value of the order parameters at the each side. Within each of these two problems, we need to solve separately for incoming (n x > 0) and outgoing (n x < 0) trajectories, so generally speaking we are solving four problems in total, Fig. 8.
The logic behind simplifying the boundary condition is as follows. Let us consider two unphysical solutions of the Eilenberger equation without an interface, see Fig. 8. Specifically, for v x > 0 we introduce the physical pairing field in the region x > 0 and an auxiliary field∆ aux in the region x < 0 and consider the solution denoted bŷ G > ∼ e −xλ/vx (λ > 0) which diverges (it must diverge since this is an unphysical solution) at x → −∞ and it vanishes as x → +∞, so that this solution vanishes as u = x/v x → ∞. Similarly,Ĝ < ∼ e xλ/vx vanishes at x → −∞ but it diverges as x → +∞ with the physical value of the pairing field for x < 0 and auxiliary order parameter for x > 0. It is easy to show that the product of these two solutions is also a solution of (A1). From the two diverging solutions we can construct the bounded solution:Ĝ b = a Ĝ <Ĝ> −Ĝ >Ĝ< . (A2) Here the normalization constant a is determined by the normalization condition for the physical solution of the Eilenberger equation: a = (Ĝ <Ĝ> +Ĝ >Ĝ< ) −1 , where we took into account the matrix structure of the quasiclassical functions, i.e. the anticommutator ofĜ < andĜ > must be proportional to the unit matrix. Thus, the physical solution (i.e. the one which remains finite in the bulk) in terms of the two unphysical ones readŝ An important property of these auxiliary matrices is that After somewhat lengthy but otherwise straightforward calculation (see below), we obtain the following expressions for the values of the quasiclassical functions at the interface as:Ĝ (A5) where coordinate dependence on the r.h.s. is suppressed.
Remarkably, these relations have the form of the circuittheory boundary conditions of Andreev refection [9,10]. Thus we were able to express the values of the quasiclassical functions, which determine the physical properties of the junctions, in terms of the correlations functions found using the auxiliary solutions. In principle, one can use these values to setup the boundary value problem and solve the Eilenberger equations anew. Indeed, the general solution of the quasiclassical equations with given ∆(x) can always be written as a linear combination of the bulk solution and the solution of the auxiliary problem, (see Eq. (D1) of Appendix D).