Quasiclassical theory of charge transport across mesoscopic normal metal-superconducting heterostructures with current conservation

We consider the steady-state nonequilibrium behavior of mesoscopic superconducting wires connected to normal-metal reservoirs. Going beyond the diffusive limit, we utilize the quasiclassical theory and perform a self-consistent calculation that guarantees current conservation through the entire system. Going from the ballistic to the diffusive limit, we investigate several crucial phenomena such as charge imbalance, momentum-resolved nonequilbrium distributions, and the current-to-superflow conversion. Connecting to earlier results for the diffusive case, we find that superconductivity can break down at a critical bias voltage $V_\mathrm{c}$. We find that $V_\mathrm{c}$ generally increases as the interface transparency is reduced, while the dependence on the mean-free path is non-monotonous. We discussthe key differences of the ballistic and semi-ballistic regimes to the fully diffusive case.


I. INTRODUCTION
The non-equilibrium distribution of quasiparticles in hybrid nanostructures has for a long time been an important topic within the field of superconductivity [1][2][3][4][5][6][7]. The interest arises naturally due to the wide range of applications of superconducting devices in sensing, metrology, thermometry and refridgeration at the nanoscale, as well as high-speed and quantum computing [8][9][10][11][12]. The non-equilibrium distribution may be a consequence of device operation, but in many cases unwanted, or poisonous, quasiparticles may limit device performance [13,14]. Also within the field of superconducting spintronics in ferromagnetic-superconducting hybrid structures, the spin-dependent distributions have been investigated [15][16][17][18][19]. Recent advances in nano-device fabrication involving unconventional superconductors [20][21][22] also raise questions about non-equilibrium effects in such systems. One fundamental problem, from the perspective of theory, is that current conservation in hybrid structures is only guaranteed if self-consistency of the relevant self-energies is taken into account.
In many cases the neglect of current conservation is a good approximation, for instance due to a very small perturbation that enables a linear response approach. Although also in this case vertex corrections can be important to include. For the pinhole geometry between a normal metal and a superconductor, the dilution of the current flowing through the pinhole into the reservoirs leads to only small corrections if self-consistency is taken fully into account [23].
In other cases current conservation is more important. For instance in nanoscale to mesoscale devices, where the non-equilibrium distribution function and the presence of current severely affect the superconducting order parameter and other self-energies [24][25][26][27][28][29][30][31][32][33]. Then it becomes crucial to take into account the conversion between quasiparticle current and superflow through Andreev processes and the spatial variation of the order parameter on the length-scale of the superconducting coherence length.
In this paper, our aim is to establish an efficient computational strategy for calculating conserved current flow in mesoscopic hybrid structures including superconducting and normal-metal elements with mean free paths ranging from the clean limit ξ 0 to the dirty limit ξ 0 , where ξ 0 = v F /2πT c is the clean-limit superconducting coherence length defined in terms of Planck's constant , the Fermi velocity in the normal state v F , and the superconducting transition temperature T c . We utilize the quasiclassical theory of superconductivity pioneered by Eilenberger, Larkin, and Ovchinnikov [34,35]. Within this Green's function method we take self-consistency of the superconducting order parameter as well as impurity self-energies into account. We limit ourselves to nanoscale devices where the inelastic mean free path in is larger than the device size, although this is not a limitation of the theory [36]. In addition, we focus on heterostructures with only one superconducting element. In this case the non-equilibrium state is stationary [37].
In the diffusive limit, when the mean free path is much smaller than the superconducting coherence length, ξ 0 , the Usadel approximation can be made [38,39]. The Usadel equations, being diffusion-type equations, may be solved numerically for instance with a finite element method [40], or by implementing the equations of Nazarov's circuit theory [41]. Going beyond the diffusive limit, it is necessary to go back to the more general Eilenberger-Larkin-Ovchinnikov formulation.
First non-equilibrium calculations including current conservation were made in the 1990s [24][25][26][27][28][29][30]. In the work by Sols and Sánchez-Cañizares [25][26][27][28][29] a scattering approach was utilized with an approximate scheme of asymptotic self-consistency in the clean limit. The importance of spectral rearrangement in the form of Doppler shifts of quasiparticle states in the presence of superflow was pointed out. Later [30], impurity scattering was also included but only very small systems of the order of three coherence lengths were studied. This paper is organized as follows. In Section II we outline the model assumptions, give details of the quasiclassical theoretical framework, and outline our compu- A superconductor (red, grated) is connected on one end to a normal-metal reservoir (blue) via an insulating barrier (black) of transparency D. The reservoir is at a chemical potential µL, while the superconductor is assumed to be grounded on the right-hand side. (b) ISI system. The central superconducting region (red, grated) of length LS is connected to normal-metal reservoirs (blue) on both ends via insulating barriers (black). The transparencies DL,R do not have to be symmetric, in this case |µL| = |µR|. (c) NISIN system. The central superconducting region (red, grated) of length LS is connected via insulating barriers (black) to a left and right reservoir (blue) via a finite piece of normal wire (yellow) of length LN rather then directly as in case (b). In these normal-metal wires the proximity effect is taken into account. The spatial coordinate along the wire main axis is denoted z. tational strategies. In section III we report results on a range of normal metal-superconducting heterostructures, sketched in Fig. 1. Section IV summarizes the paper.

II. MODEL AND METHODS
In this paper we shall study a number of normal metalsuperconducting devices with steady-state current flow between source and drain contacts, see Fig. 1. The superconductor is considered to be at ground potential in equilibrium, while in non-equilibrium we also compute the local electrochemical potential φ(z) self-consistently. This potential describes the local difference in potential between quasiparticles and the condensate. Our study is thus directly related to the charge imbalance problem studied for a long time and recently reviewed in Ref. [7]. Here we revisit this problem and compute all quantities self-consistently to guarantee current conservation. Moreover, we go beyond the diffusive approximation and cover the whole range of mean free paths from the clean to the dirty limits.
We shall consider a number of set-ups, starting in Section III A with a normal metal reservoir tunnel coupled to a superconductor which we assume continues to a grounded reservoir on the other side, see Fig. 1(a). This is in principle the usual set-up in the scattering approach to tunneling into a superconductor, pioneered for the ballistic case by Blonder-Tinkham-Klapwijk (BTK) in Ref. [23]. It has been extended and widely used also for unconventional superconductors, see, e.g., the reviews in Refs. [42,43]. The approach typically neglects selfconsistency and current conservation, which is a good approximation when there is a Sharvin, or pinhole-like, contact to the superconductor that guarantees current dilution into the contacts. We are here interested in the corrections when this assumption is not strictly valid, similar to the situation in the experiments in Refs. [31,33].
We limit our study to small mesoscopic devices at a low temperature T = 0.01T c with good contacts to the reservoirs. More specifically, the system size L is assumed sufficiently small and the contacts to the reservoirs are assumed to be sufficiently good that the inelastic scattering time is the longest time-scale in the problem and for instance electron-phonon relaxation and recombination processes can be neglected [44][45][46], see discussion below. As the main consequence, we cannot properly describe charge imbalance relaxation when the applied voltage is sufficiently high that a considerable amount of quasiparticles are injected into continuum states in the superconductor in the setup shown in Fig. 1(a). The potential φ(z) becomes finite far into the bulk of the superconductor and charge imbalance does not decay. In this case, the setup in Fig. 1(a) becomes ill-defined in the sense that we can not calculate the distribution function far in the superconductor. We therefore consider the well-defined setup in Fig. 1(b), where there is a normal metal reservoir also on the right side. In this case we study charge imbalance self-consistently on distances smaller than the inelastic scattering length L in . Finally, we study the case when parts of the normal metal leads to the superconductor are influenced by the proximity effect, see Fig. 1(c). This leads to a series resistance that depends on both system size and elastic mean free path.
To justify that inelastic scattering can be neglected at low temperature and small device sizes, we should compare the electron-phonon scattering times, in particular the recombination time τ r , and the time for diffusion from one contact to the other, or dwell time in the device, τ d . Kaplan et al. [44] showed that in a bulk superconductor with thermal distributions the electron-phonon recombination time is exponentially long at low temperature, τ r ∝ τ 0 e ∆/kBT , where ∆ is the superconducting gap, while the electron-phonon scattering time follows a power-law, τ s ∝ τ 0 (∆/T ) 7/2 . The typical time-scale τ 0 ∼ 0.1 − 1 ns for strong-coupling superconductors, but is as long as τ 0 ∼ 500 ns in weak-coupling superconductors. When the quasiparticle distribution is not thermal, experiments show a wide range of low-temperature recombination time scales ranging form 10 µs up to ∼ 1 ms [13,14,47]. Catelani et al. [46] showed that in the diffusive limit and at low voltage, the non-equilibrium population in a small superconducting device induced by quasiparticle injection from normal metal reservoirs decays even weaker through electron-phonon scattering processes if the quasiparticle life-time due to escape to reservoirs is small τ d τ 0 ∆/T e ∆/T . In our studies below we shall see that the non-equilibrium occupation of continuum states is indeed very small, as in Ref. [46], and the moving condensate (finite superfluid momentum below) carries the current in the interior of the devices as a result of Andreev processes near the interfaces. We note in passing, that such processes do not lead to heat flowing into the superconductor [48]. Let us compare the recombination time τ r to the escape time to reservoirs τ d . In the normal state we use the Thouless energy E T = D/L 2 with the diffusion constant D = 1 3 v F to estimate τ d = /E T ∼ 20 ns, assuming a typical Fermi velocity v F ∼ 10 6 m/s, system size L ∼ 100ξ 0 , with the clean limit ξ 0 ∼ 1 µm, and mean free path ∼ ξ 0 . Since we assume high transparencies of the barriers to the leads, the escape time is dominated by the diffusion time. Thus, in the normal state the escape time to reservoirs is much shorter than the electron-phonon scattering time, τ d τ 0 . In the superconducting state the low group velocity of quasiparticle states leads to a reduction of the diffusion constant as D → D √ ε 2 − ∆ 2 /|ε| and the escape time τ d is enhanced at the gap edge. Experimentally it was measured to give approximately a factor of 10 reduction of the effective diffusion constant at low temperatures [49]. At the same time it should be compared with τ r typically enhanced at low temperatures in the superconducting state by several orders of magnitude.
In conclusion, we assume good contacts to reservoirs allowing for efficient Andreev processes, a weak-coupling superconductor not in the extreme dirty limit (we let ≥ 0.1ξ 0 ), and system size L small enough that τ d τ r , τ s holds and neglect electron-phonon processes.

A. Quasiclassical theory
To study stationary current flow in normal metalsuperconducting heterostructures under voltage bias, we utilize the nonequilibrium quasiclassical Green's function method [34,35,39,50,51]. All observables are then expressed in terms of three quasiclassical propagators. The retarded and advanced Green's functions contain information about the spectral properties, while the Keldysh Green's function also contains information about the nonequilibrium distribution function. In superconducting systems, the self-consistency condition on the order parameter leads to a coupling of spectrum and occupation that complicates the physical interpretation. Nevertheless, within the quasiclassical formulation, the two may be disentangled at the end of the calculation. The retarded (R), advanced (A), and Keldysh (K) propagators depend on momentum direction p F on the Fermi surface, coordinate R, and energy ε. We arrange them into a matrix in Keldysh space, denoted by a check, Each propagator is in turn a matrix in combined particlehole (Nambu) and spin spaces. We mark Nambu matrices by a hat (ˆ) and denote the three Pauli matrices in Nambu space asτ i and in spin space as σ i , where in both cases i = 1, 2 and 3. The charge-current density is obtained from the Keldysh propagator: Tr where e = −|e| is the electron charge, N F is the density of states at the Fermi level in the normal state, and v F is the Fermi velocity. The trace is over both Nambuand spin spaces. In the second line we also introduced the spectral current density j(R, ε). In three dimensions with a spherical Fermi surface, as we will assume in this paper, the average is given by integration over polar and azimuthal angles, For a more complicated Fermi surface, see, e.g., Ref. [52]. The local electro-chemical potential is determined by a requirement of local charge neutrality [53]. It has the form In equilibrium this potential is zero, while under steadystate non-equilibrium conditions it signals charge imbalance induced by injection of quasiparticles from normalmetal reservoirs. To achieve current conservation, it is required to compute φ(R) self-consistently with the selfenergies. The Green's functions are obtained by solving the quasiclassical transport equation, in combination with the normalization condition, Above and in the following, when there is no cause of confusion, we drop the explicit reference to the dependences on p F , R, and ε.
The matrixȟ with self-energies appearing in Eq. (6) has the follwing elements in Nambu space: Objects with a tilde are related to non-tilde objects through particle-hole conjugation: In this paper, we include self energies for the mean-field superconducting order parameter as well as scalar and spin-flip impurities. They add up according tǒ h =ȟ mf +ȟ s +ȟ sf .
The superconducting order parameter is obtained through the self-consistency equation where λ < 0 is the coupling strength, and ε c is an energy cut-off. The spin-singlet symmetry of the order parameter is guaranteed by including the factor iσ 2 in the trace. Using the standard trick, where the linearized gap equation is added and subtracted, we eliminate the coupling strength λ and the high-energy cutoff c from the gap equation in favor of the clean-limit superconducting transition temperature T c . In Nambu space we now havê while the Keldysh part is zero,ĥ K mf = 0. The complexvalued order parameter can be written in terms of its magnitude and phase, ∆ 0 (R) = |∆ 0 (R)| exp(iχ(R)). In practice we choose to work in a gauge where the order parameter is real and the phase enters through the superfluid momentum, defined as We assume that the lateral dimensions of the heterostructure are small compared with the penetration depth. In this case the back-coupling of the electromagnetic gauge field can be neglected. Scattering against scalar impurities are taken into account within the self-consistent Born approximation, where the scattering time τ is related to the mean free path via = v F τ . Additional spin-flip impurities are included in normal-metal regions through the self-energy The mean free path for spin-flip scattering sf = v F τ sf will always be large compared with the mean free path due to scalar impurities, sf , and is included in order for the proximity effect to decay. In this way we restrict ourselves to setups with normal-metal leads in Fig. 1(c) where superconducting coherence has decayed at the places we connect our device to reservoirs.
In the numerical implementation we use the Riccati parametrization of the Green's function [50]. The normalization condition in Eq. (7) is then automatically fulfilled. The components of the retarded and advanced Green's functions are written aŝ where The coherence functions are local amplitudes for conversion between particlelike and holelike branches (e ↔ h) and therefore express superconducting coherence. These functions satisfy the following Riccati equations i The components ofĝ K in particle-hole space are labeled asĝ We use a parametrization of these elements in terms of distribution functions x(p F , R, ε) andx(p F , R, ε), The distribution function x satisfies the following transport equation: The corresponding equation forx can be obtained through Eq. (9). We note that the coherence functions are in general matrices in spin space, but for the singlet superconducting case we study here they are simply proportional to iσ 2 . The distribution functions, on the other hand, are trivially proportional to a unit matrix in spin space. Note that the spatial coordinate along the system main axis is denoted z, as in Fig. 1, and should not be confused with the distribution function x(p F , R, ε).
In Eqs. (19)- (20) and Eq. (24) each velocity direction v F (p F ) defines a trajectory in real space. We obtain the coherence functions and the distribution functions by integrating the first-order differential equations starting from initial conditions at the reservoirs. The equations for γ R and x are stable to integrate in the direction of v F , while the equations forγ R andx are stable to integrate in the opposite direction. The opposite holds for the advanced coherence functions. We use the operator method presented in Ref. [51] for stepping on a spatial grid with piece-wise constant self-energies [54].
Internal boundaries, such as tunnel barriers, have to be included through boundary conditions. These boundary conditions are given as Eqs. (32)- (35) in Ref. [50].
The key input is the interface scattering matrix, which is expressed in terms of transmission and reflection amplitudes for incoming electrons in the normal state: where R(p F ) + D(p F ) = 1. The angle between the incoming momentum directionp F and the interface normal n define a directional cosine asp F ·n = cos θ F . The angular dependence of the transmission function is taken as a tunnel cone with the form The parameter β tunes the width of the cone, while D 0 sets the transparency for perpendicular incidence.

B. Distribution functions
In order to simplify the numerics, as well as getting more insight into the physics, it is convenient to split the distribution function into a local-equilibrium part x le (p F , R, ε) and a so-called anomalous part x a (p F , R, ε) by writing The local-equilibrium part is defined according to where the familiar equilibrium form of the distribution function is The remaining anomalous part of the distribution function, x a , captures the essential non-equilibrium effects.
The splitting of x according to Eq. (27) carries over to a corresponding splitting of the Keldysh part of the Green's function:ĝ This naturally leads to a separation of all observables that involveĝ K into a local equilibrium part and an anomalous part. When this splitting is used, only the anomalous part of the distribution has to be stepped along trajectories since x le is entirely determined by local quantities. For the stationary state, the equation of motion for x a is given by Eq. (24) with the replacements [51] x → x a , The distribution function x(p F , R, ε) can be related to the distribution function h(p F , R, ε) commonly used in literature [3,51]. The two are related through which has the inverse In a normal metal with γ = 0, the two distributions are identical, while in the presence of superconducting coherence this is no longer the case. The distribution h can be divided into two parts with different structure in electron-hole space: One obtains an electron distribution function through which in equilibrium reduces to the Fermi-Dirac distribution function. The distribution x can not be split in a similar fashion because it mixes particle and hole distributions. It is convenient to focus on x andx in the numerical implementation, while for interpreting the final results it may be beneficial to transform to the distribution h, or f e/1/3 . In terms of h, the Keldysh Green's functionĝ K has the form Correspondingly, all observables can be expressed either in terms of x or h types of distributions. This also makes the role of the two components f 1 and f 3 more transparent. The identities relevant for this paper are for Eq. (11), for Eq. (5), and for Eq. (2). For the current we use as natural unit is a normalized momentum-resolved local density of states per spin: The total local density of states is obtained as Finally, it is important to remember that all nonequilibrium distribution functions discussed above depend on momentum direction p F . When impurity scattering is introduced, they become increasingly isotropic as the mean free path is reduced. It is then beneficial to study the Fermi surface averaged distributions: They also have a certain symmetry with respect to energy. Following the nomenclature in Ref. [3], f L is the longitudinal mode that is an odd function of energy, while f T is the transverse mode that is an even function of energy. In the diffusive limit, as in the Usadel formulation, these isotropic distributions also appear in observables and are the natural objects to study. In our formulation valid for arbitrary mean free path, only momentum resolved distributions, such as x, h, or f e/1/3 , are used in the calculations. Their properties are different than f L and f T appearing in the diffusive Usadel theory. For instance, their symmetries are governed by Eq. (9).

C. Local chemical potentials
The local chemical potential φ(R) is given by Eq. (5). The local value φ(R) determines the local-equilibrium component x le (p F , R, ε) in Eq. (28), and the driving terms for the anomalous component x a (p F , R, ε) in Eqs. (31)- (33). We further define local right-mover φ + (R) and left-mover φ − (R) quasi-potentials as where the label + (−) indicates positive (negative) momentum projection along the z-axis and X a indicates that only the anomalous distribution x a is used in Eq. (22). The averages in Eq. (46) are defined as where ξ F = cos θ F . The two quasipotentials φ ± (R) in Eq. (46) are measures of occupation of right-and leftmoving quasiparticle states relative to the local chemical potential φ(R). For systems not in the diffusive (Usadel) limit, these two quantities are more useful to develop an understanding of the transport behavior than the averaged chemical potential φ, see also the discussion in Ref. [55]. These potentials do not imply that the non-equilibrium distribution for left-and right-movers are simply Fermi distributions centered around these quasipotentials. Note also that the averages over half the Fermi surface in Eqs. (47) are missing a factor of 1/2 in comparison to Eq. (4). With this definition, in a normal metal meaning that the measurable chemical potential φ is the average of left-and right-mover quasipotentials [55]. Note however, that Eq. (48) is not satisfied in the superconducting state because of the proximity and inverse proximity effects.

D. Voltage boundary condition
In the case of voltage bias, we treat the normal regions as reservoirs at a chemical potential µ L = eV /2 and µ R = −eV /2, respectively. They are assumed to be unaffected by the contact to the central superconductor of length L. Thus, we are neglecting both the proximity effect as well as the effect of current flow from or into the normal leads. Assuming such reservoirs, the distribution functions have the boundary conditions The coherence functions incoming from the normal metal are assumed to vanish. Within the quantities introduced above, these distribution functions enforce a certain chemical potential of right-movers (left-movers) at the left (right) edge of the system.

E. Current boundary condition
In left-right asymmetric systems, for instance an ISI structure with different interface barrier transparencies, or an NISIN structure with different lengths of the normal metals, we can not a priori know the potential profile. In this case we enforce a certain current I b through the system boundaries and compute the resulting chemical potentials in left (µ L ) and right (µ R ) reservoirs, as well as the potential drop through the system self-consistently by requiring current conservation. The resulting voltage drop is eV = µ L − µ R , and would experimetally correspond to the externally applied voltage that result in the current I b . We illustrate the procedure for the left system edge, where we need a boundary condition for the incoming distribution function x + .
To simplify the explanation, we assume a singletrajectory model so that only the two trajectories with ξ F = ±1 enter the calculation. The generalization to a full Fermi-surface average only adds trivial prefactors. Assuming that the normal-metal region is large enough for the proximity effect to decay, the current at the system edge can be written as where we removedx by symmetry. A possible choice for x + is where φ b + is the boundary value that will be specified in the following. Using this choice and the requirement I ! = I b , we rewrite Eq. (51) as This can be used to determine φ b + and thus the incoming distribution function x + (0) iteratively starting from the initial guess x − (0) = 0. In each iteration we obtain a new value of x − (0) that adjusts φ b + so that the total current at the edge is equal to I b . At the right system edge, the roles of x + and x − are swapped.
Once a self-consistent solution is obtained, we have charge-current conservation I(z) = I b everywhere in the structure. The procedure determines the chemical potential at the left edge µ L = φ b + and right edge µ R = φ b − requiring a specified current I b throughout the structure. In a completely normal-metal system the difference between the edge potentials will, for a given current, be dependent on the total resistance between the two ends.
In the case of a central superconducting region, the respective edge potential will similarly be dependent on the resistance of the normal metal regions and interface resistances before the current becomes transformed into supercurrent in the superconductor, i.e., a finite p s (z). This approach is particularly useful for asymmetric structures, where the bias drop over the structure does not have to be symmetric. This scheme is then a theoretical description of a four-point measurement: for a certain current flowing through the structure from source to drain leads, we determine the potential drop that would be measured by an additional high-resistance voltage probe. This scheme is thus highly relevant to experiments where a superconductor is dc-current biased and the potential drop is measured.

F. Details on numerical procedure
Performing a fully self-consistent nonequilibrium calculation is numerically challenging since Green's functions and self-energies have to be determined on the real axis. Here, we wish to outline some details on our numerical approach.
Firstly, all energies obtain a small imaginary part ε → ε ± iη, when calculating retarded (+) and advanced (−) quantities. We chose η = 10 −3 k B T c . In principle, this determines the required energy resolution to be ∆E ≈ η, and energies up to a cutoff energy E c ∆ 0 , Γ, eφ have to be included. For the main results of our paper, we chose E c = 500k B T c . In order to reduce the resulting large number of energy points, we use a non-uniform grid that is more dense in an energy interval ε ∈ (−5∆, 5∆). Most of the physically relevant information is contained in this energy window. For example, the distribution function x a is only nonzero in the bias window |ε| ≤ µ L,R .
The number of iterations necessary to achieve selfconsistency is greatly reduced by iterating p s (z), keeping the order parameter real, instead of iterating a complexvalued order parameter. The real order parameter is updated through the gap equation in Eq. (11). For the quasi-1D setup, we update the superfluid momentum p s (z) using the condition of current conservation: In the case of the current boundary condition discussed in Sect. II E, the constant on the right-hand-side is the enforced boundary current j b . On the other hand, if a potential boundary condition is used then the right-hand side will be given by the current j int at the interface between a normal and a superconducting region. Exactly at the interface grid points, both ∆ 0 and p s have to be kept zero [61]. On the normal side we specify both distribution and coherence functions while on the superconducting side the incoming functions are the result of propagating both functions through the self-energy landscape. An equivalent formulation of Eq. (54) for the potential boundary condition is then with an analogous expression for the current boundary condition. The local current deviation before selfconsistency, δj(z) = 0, can then be used to update p s (z) throughout the superconductor. Starting from an initial guess p s = 0 everywhere in the superconductor, we use where n is the iteration index and q is the update step size of order unity. The interface current j int will typically change during iterations but will eventually reach a fixpoint. In the results we present here, current is conserved up to a local relative error of δj(z)/j int < 5 · 10 −3 but higher accuracy can be achieved.
In a majority of our calculations the Fermi surface average does not alter our results. We therefore use a onetrajectory model for simplicity, i.e., we keep θ F = 0 and π only. In relation to the tunnel cone in Eq. (26), this implies D = D 0 . In the few cases where the full angular dependence affects our results, we have chosen a wide tunnel cone with β = 1.

G. Example: Normal metal
Let us start with text book examples [55,56] in order to set our formalism into perspective. In a device as in Fig. 1(b) but consisting of normal-metal regions only, all superconducting coherence functions γ R vanish and transport is described by the distribution functions only. The transport equation for x simplifies to the Boltzmann equation for a normal metal, Assuming homogeneous current flow in the transverse direction, the current flowing along the wire main axis, quantified by a coordinate z, simplifies to where A is the wire cross sectional area. When only elastic impurity scattering is considered, combining Eqs. (57)-(59) implies conservation of spectral current, j(ε) = −eN F v F f 1 FS . From these equations it is straight forward to perform the standard linear response calculation to an applied electric field along the z-axis, E = −∂ z φ(z), with the local equilibrium distribution in Eq. (29). The resulting charge conductivity is σ N = 2e 2 N F D for a dirty normal metal.
For the device geometries, we go beyond linear response and perform calculations using assumptions of scattering theory with the device coupled to reservoirs [56]. In Fig. 2 we show how the chemical potential, as well as left-and right-mover quasi-potentials, drop across a normal metal piece between two reservoirs placed at z = 0 and z = 5ξ 0 , for different mean free paths. We note that the natural length scale is the mean free path, but although there is no superconductor here we use the coherence length to be consistent with later sections of the paper.
In the clean case (l → ∞), see Fig. 2(a), the incoming left-and right-moving populations simply shoot through the structure and do not mix, hence φ + = µ L , φ − = µ R , and φ = 0 everywhere. This corresponds to a Landauer-Büttiker wavefunction based scattering approach for the trivial case of no scattering in the device. The applied bias is eV /Tc = 1, and the mean free path l/ξ0 is (a) ∞, i.e., ballistic limit, (b) 10, (c) 1, and (d) 0.1, diffusive limit.
In the diffusive limit, in contrast, the chemical potential φ at the system edges equals the reservoir chemical potentials and in between interpolates linearly between the two. There is no difference between left-and right-movers quasi-potentials in the diffusive limit, see Fig. 2(d), since isotropization is to lowest order approximation local (l L). The fully diffusive limit can also be obtained from the Usadel equation [32]. The solution for the isotropic distribution function reads where x L and x R are the left-edge and right-edge boundary values, respectively. This leads to so that in the diffusive limit, the boundary condition for φ ± at the edges becomes effectively boundary conditions for φ. The solution is then indeed in agreement with Fig. 2(d).
For mean free paths in between these two limiting cases, as in Fig.2(b) and Fig. 2(c), the behaviour of φ as well as φ ± interpolate between the ballistic and diffusive cases. We note that since we compute the chemical potential φ everywhere in the device, its edge values will disagree with the reservoir chemical potentials. It is well known [55,56] that this differences correspond to the contact potentials, or spreading resistance [31], via the adiabatic and ballistic leads to the reservoirs. Since we go beyond linear response in our treatment, we need to have well-defined reservoirs with well-defined incoming distribution functions (φ b ± ), in agreement with the Landauer-Büttiker approach [56].

Asymmetric system: current boundary condition
Next let us study an asymmetric device of total length L = 15ξ 0 , with two internal barriers at z = 5ξ 0 and z = 10ξ 0 of transparencies D L = 0.3 and D R = 0.8, see Fig. 3. In this case, since we a priori do not know how the voltage drops across the system, we need to enforce a current and let the chemical potentials at left and right edges float up and down to reach self-consistent values while enforcing current conservation across the system. An example of the resulting total voltage drop eV , chemical potential profile φ(z), and the corresponding electronic distribution function f e ( , z) are shown in Fig. 3 for an intermediate mean free path l/ξ 0 = 1 and an enforced current I b = 0.1I 0 , where I 0 = j 0 A. This result for the distribution function agrees well with the experiment by Pothier et al. [57].
Having established how our formalism works in the normal state, we shall from now on focus on the case when the central region (in this example z/ξ 0 ∈ [5, 10]) is superconducting.

A. I-S system
The first setup in Fig. 1(a)   chemical potential in the normal metal to µ L = eV relative to the superconducting reservoir µ R = 0. The contact specifies a boundary condition for the incoming distribution function x + (z = 0 − , ε) according to Eq. (49). The boundary condition for the spectral part is γ R (z = 0 − , ε) = 0, i.e., the proximity effect is neglected in the normal-metal side in agreement with the reservoir assumption. If self-consistency and current conservation is neglected in the superconductor we obtain the BTK result for the interface conductance [23]. We go beyond this approximation and examine the effect of current conservation as well as impurity scattering on transport.
The first thing to note is that Anderson's theorem [58], which states that the superconducting order parameter is unaffected by scalar impurities, relies on time-reversal symmetry. Therefore it will not hold in the presence of current flow. Even far from the contact, the current flow induces Doppler shifts of right-and left-moving quasiparticles according to where the superfluid momentum p s is set by the gradient of the order parameter phase χ as defined in Eq. (13). These Doppler shifts lead to violation of Anderson's theorem and the order parameter depends on the mean free path. In Fig. 4 we show results of a self-consistent calculation for one particular set of parameters and for one voltage. The boundary condition for the distribution function enforces a certain anomalous current I a at the interface that slowly decays into the superconductor, see Fig. 4

(b).
Since current is conserved, this decay is compensated for by a corresponding increase in the local-equilibrium current I le . The conversion is associated with the appearance of the chemical potential φ(z) near the interface, see Fig. 4(a). The local-equilibrium current is carried by a finite, position-dependent superflow in the superconductor, meaning that the order parameter phase varies with coordinate and the superfluid momentum is non-zero, see Fig. 4(d). The suppression of the order parameter close to the interface, shown in Fig. 4(c), leads to a peak in the superfluid momentum p s (z) before both reach their bulk values roughly 10ξ 0 away from the interface. The absolute value of the total current depends on interface barrier transparency, impurity concentration in the superconductor, and the applied voltage.
The spectral current and the distribution functions are shown in Fig. 5. The anomalous current in Fig. 5(a) is due to the quasiparticle injection in the subgap region within a voltage window 2eV . Through Andreev reflection this current is converted to superflow, which is carried by continuum states, see Fig. 5(b). In the inter-face region, where the potential φ is non-zero, also the transverse distribution function f T is non-vanishing as shown in Fig. 5(d). Well inside the superconductor it has decayed to zero. At the same time the longitudinal distribution decays back to the equilibrium form and the supercurrent is only due to the Doppler shifts in Eq. (62) of continuum states.
In Fig. 6 we display the total current as function of applied bias voltage for three different interface transparencies, and three different mean free paths. For the clean case (solid lines, = 10ξ 0 ) the current-voltage characteristics agree well with the BTK results (dotted lines), in particular at low voltage. Corrections appear at higher voltage and are mainly due to the finite Doppler shifts. The solid lines end where superconductivity near the interface starts to break down, an effect that we will study in more detail below.
In the case of a highly transmissive barrier, D = 0.8 in Fig. 6(a), the current gets reduced with shorter mean free path. In the presence of impurities, quasiparticles that are transmitted through the interface can be scattered back before they get Andreev-reflected. For a given injection energy eV , Andreev reflection happens on a lengthscale ξ AR = v F / √ ∆ 2 − eV 2 , while scattering happens on the scale of . The reduction of Andreev reflection and current is thus stronger for higher impurity concen- tration and at higher voltages.
In contrast, the current is increasing with reduced mean free path for barriers with lower transparency, see the cases of D = 0.5 and D = 0.3 in Fig. 6(b) and Fig. 6(c). This behavior can be understood from two points of view. First, by including impurity scattering the subgap local density of states becomes finite near the interface, even in equilibrium, due to the so-called inverse proximity effect that occurs on the coherence length scale [31,59]. From a transport perspective, Andreev reflection is already reduced by the low barrier transparency. On the other hand, normal reflection gets reduced by impurity scattering. The impurity scattering filling the subgap energies with states is associated with a finite normal electron transmission amplitude locally at the interface, which is not present in the ballistic case. This enables a larger current for the same applied voltage. Note that in the bulk superconducting region there is only superflow, meaning that the local finite amplitude of normal transmission at the interface is not associated with quasiparticle flow in the asymptotic region, which also remains fully gapped.
In the experiment presented in Ref. [31], the resistance of an Al wire was measured as function of temperature. The fit to Usadel theory was good, but at low temperature the resistance was lower than what was predicted by theory. According to Ref. [31] the transparencies of the interfaces to the normal-metal reservoirs were very high, corresponding to our Fig. 6(a). We predict that the correction to the Usadel result from ballistic effects is a decrease of the IS resistance [enhanced conductance in Fig. 6(a)]. This follows the trend seen in the experiment at low temperature.

Critical voltage
All of the above current-voltage curves end at certain voltage points that we refer to as critical voltages V c following Ref. [32]. At the critical voltage the order parameter either vanishes throughout the system, i.e. superconductivity breaks down, or self-consistency can not be achieved. The latter case happens in the clean limit for high barrier transparency, where the current is sufficiently high that the Doppler shifts lead to injection into continuum states. In this case, the potential φ(z) also extends out into the bulk of the superconductor. We will return to this case in the next section, when we have two normal-metal reservoirs and a more well defined set-up.
The breakdown behavior can be understood by a reexamination of the self-consistency equation for the order parameter ∆ 0 . Using the parametrization ofĝ K in terms of distribution functions f 1 and f 3 in Eq. (36), the selfconsistency equation for ∆ 0 can be written as in Eq. (39).
Both terms consist of a spectral contribution, F R ±F A , and either the energy-like mode f 1 or the charge-like mode f 3 . Note that both distributions are real functions. We use a gauge where ∆ 0 is real and obtain p s , which gives the spectral rearrangements. For the real part of the kernel in the right hand side, we find that the first term gives the dominating contribution to ∆ 0 . Figure 7 shows the real part of the first term of the kernel in Eq. (39) for a voltage close to V c . For comparison, the equilibrium forms are included as dotted lines. In Fig. 7(a) we show the factor Re F R + F A , while in Fig. 7(b) we display the distribution function f 1 for right-and left-moving states. The main effect of the current injection is that the distribution function is reduced in an energy window of width 2eV around zero energy. In addition, the coherence peaks in Re F R + F A are Doppler shifted towards lower energies. As a result, the total kernel gets reduced with increasing voltage, see Fig. 7(c). For voltages approaching V c , the coherence peaks gets shifted inside the energy window where f 1 is strongly reduced. This leads to a suppression of the order parameter and eventually turns the system normal.
We summarize in Fig. 8 the critical voltage as function of barrier transparency for three mean free paths from the clean to the dirty limit. The general trend is that lower transparency barriers lead to a higher critical voltage. This is due to the reduced current and reduced Doppler shifts for the same applied voltage at low transparency. For transparencies D 0.75, we observe an increase in the critical voltage with decreasing mean free path. The reduction in current with decreasing mean free path, shown in Fig. 6(a) above for D = 0.8, means that p s goes down and the coherence peaks get shifted less, allowing for a larger critical voltage. In the ballistic case, a high-transparency interface can have steady-state solutions with non-decaying chemical potential in the bulk of the superconductor. In this case we inject into continuum states and the non-equilibrium population does not decay deep in the superconductor since we neglect inelastic relaxation processes. This situation for the IS set-up is beyond the scope of this paper, but we return to this case in the next section where we have two reservoirs and a well defined set-up. The corresponding points where injection into continuum state appear are marked by empty symbols rather than filled ones in Fig. 8.
In contrast, the critical voltage decreases with impurity concentration for transparencies D 0.75. We understand this result in the light of Fig. 6. For smaller transparencies (D = 0.5 and D = 0.3) the current for a particular applied voltage increases for smaller mean free paths. This leads to more pair breaking and a lower critical voltage.
Finally, we note that in the diffusive limit, the critical voltage is maximal in the tunneling limit and decreases towards a minimum at full transparency. For the case of D = 1, we recover the critical voltage for the diffusive NSN system in the zero-temperature limit [32], indicated by the red, dashed line.

B. I-S-I system
In the ISI case, we choose a symmetric system with D L = D R = D, and a symmetric bias of µ L = −µ R = eV /2. Similar to the IS system, we find a critical voltage where superconductivity breaks down which, in turn, depends on interface transparency and mean free path, see Fig. 9. The main new result compared with the IS case is that for relatively clean systems with hightransparency barriers D 0.8, we can study solutions with non-vanishing φ everywhere in the superconductor. Such solutions typically have oscillations in both the order parameter, the superfluid momentum, and the division between anomalous and local-equilibrium currents, see into the bias window. Under voltage bias these states become occupied through an electron-hole transmission process. Within the 1D model, these states interfere with returning quasielectron continuum states which results in oscillations analogous to Tomasch oscillations, with a voltage-dependent wavelength that does not depend on system size. For shorter mean free paths of the order of the oscillation period, the interference is suppressed. In the diffusive limit they are absent.
Going beyond the 1D model by performing a full circular Fermi-surface average, the oscillations survive for barriers with a wide tunnel cone, while they are suppressed for narrow ones. The narrow tunnel cone leads effectively to a lower barrier transparency and less dramatic Doppler shifts. The transition from non-oscillating to oscillating solutions becomes less sharp as function of voltage. The quasihole states are not available for all trajectories at once due to the angular dependency of the Doppler shift in Eq. (62). Typically, small-scale oscillations set in first, and increase in magnitude as more trajectories enter the bias window with increasing applied voltage. In summary, only the quasi-1D model with a wide tunnel cone can display pronounced oscillations.

C. N-I-S-I-N system
In this section we include the effect of the proximity effect in the normal metal sides. This means that we examine what happens when the central superconducting region is not immediately connected to perfect reservoirs or perfectly ballistic leads, but those reservoirs are located at some distance away from the barriers. We study first a symmetric system, with equal transmissions of the two insulating barriers connecting the superconductor to the normal-metal regions, D L = D R = D, and consider a symmetric bias µ L = −µ R = eV /2. In Fig. 11 the potential drops in the normal-metal regions are shown for different mean free paths. In the ballistic case, Fig. 11(a), the lead is only weakly proximitized and the potential drop is almost linear. For intermediate to small mean free paths, the linear potential drops are increasingly altered by the presence of the proximity effect. The interior of the superconductor is at this applied voltage at ground potential. The enhanced proximity effect in the normal metals reduces the local potential φ(z) near the N-I-S interface as compared to the normal case, i.e., the N-I-N system studied in Section II G. As a consequence, the proximity effect leads to a larger voltage drop in the normal metal. Fig. 12 shows a comparison of the critical voltages for the NISIN system and the corresponding ISI system. For transparencies D 0.75, we see a slight, almost constant increase of the critical voltage in the NISIN system. The normal "lead" gives rise to an additional resistance compared to the ISI case where only the interface resistance determines the current flow.
At higher transparencies, D 0.8, the critical voltage starts to increase again in the NISIN system and plateaus for D 0.9, while it monotonically decreases in the ISI case. In this limit, the resistance of the normal lead dominates over the small interface resistance. For this set of parameters, the behavior changes at D ≈ 0.75. In general, the turning point will depend on the lead resistance compared to the interface resistance and will thus depend on the length of the normal lead and the mean free path.
Lastly, we assume that the NISIN structure is not symmetric but rather has different interface transparencies on the two sides. In this case, the potential profile is not symmetric and we must use the current-bias scheme introduce in Sect. II E. We specify a boundary current I b and find the potentials φ b ± that have to be enforced at the system edges. Fig. 13 shows an example result of one such calculation. As can be seen, the potential applied on the two sides will be different depending on the interface transparencies. This leads to spectral currents of different widths 2φ b ± in the normal-region leads, as shown in Fig. 14(a). The asymmetry between left and right interfaces are also reflected in the distribution functions, see Fig. 14(b)-(c). Still, in the bulk interior, the distributions have decayed to equilibrium forms and the current is only due to spectral rearrangements due to p s .

IV. SUMMARY
In this paper, we studied charge transport through mesoscopic, normal metal-superconductor wires using quasiclassical theory. Such hybrid systems with arbitrary mean free paths were considered, extending the fully ballistic or fully diffusive limits studied in literature. We performed charge-current conserving calcula-tions and studied phenomena such as charge imbalance, the conversion of quasiparticle current to superflow, and the critical voltage of the superconductor. For normal metal-superconductor interfaces that are not of pinhole-type, we observed a critical voltage V c at which the superconductor turns normal. This connects with a similar observation for the fully-diffusive system with a fully-transparent interface [32]. The critical voltage was found to be the result of an interplay between Doppler shifts and the injected nonequilibrium energy mode. We investigated the influence of both the mean free path and interface transparency on the critical voltage. In general, V c increased with smaller interface transparency as this reduced the current through the structure. The effect of impurities on the critical voltage was found to depend on the interface transparency. Higher impurity concentration increased V c for very transmissive interfaces but reduced it for all transparencies below D = 0.75. We obtained similar results for superconductors connected on both ends to normal-metal leads. However, we find additional oscillating solutions in the ballistic regime of > ξ 0 , which led to an increased critical voltage in such systems. Lastly, short pieces of normalmetal leads were included to study the influence of the proximity effect on transport. We found only small corrections for weakly-proximitized ballistic systems. For intermediate to diffusive systems, the additional resistance of the normal lead increased the critical voltage compared to when the proximity effect is neglected.
It should be noted, that the estimate of V c from our calculations are obtained neglecting many phenomena that may become of importance near breakdown, such as for instance the electron-phonon interaction. It would be of interest to improve the model and consider the highly non-equilibrium distribution function in Fig. 7 and its influence on inelastic processes, but this is beyond the scope of the present paper.
The method introduced in this paper is applicable to any kind of mesoscopic superconducting systems for arbitrary mean free path. It can thus be used to study other phenomena, such as spin-charge density separation [17,60], heat flow [62], or unconventional superconductors with other order-parameter symmetries [42,43], where the surface and interface physics may be nontrivial even in the equilibrium state [63][64][65].