Self-consistent theory of current injection into $d$ and $d + is$ superconductors

We present results for the steady-state nonlinear response of a $d_{x^2-y^2}$ superconducting film connected to normal-metal reservoirs under voltage bias, allowing for a subdominant $s$-wave component appearing near the interfaces. Our investigation is based on a current-conserving theory that self-consistently includes the non-equilibrium distribution functions, charge imbalance, and the voltage-dependencies of order parameters and scalar impurity self-energies. For a pure $d$-wave superconductor with [110] orientation of the interfaces to the contacts, the conductance contains a zero-bias peak reflecting the large density of zero-energy interface Andreev bound states. Including a subdominant $s$-wave pairing channel, it is in equilibrium energetically favorable for an $s$-wave order parameter component $\Delta_\mathrm{s}$ to appear near the interfaces in the time-reversal symmetry breaking combination $d+is$. The Andreev states then shift to finite energies in the density of states. Under voltage bias, we find that the non-equilibrium distribution in the contact area causes a rapid suppression of the $s$-wave component to zero as the voltage $eV\rightarrow\Delta_\mathrm{s}$. The resulting spectral rearrangements and voltage-dependent scattering amplitudes lead to a pronounced non-thermally broadened split of the zero-bias conductance peak that is not seen in a non-selfconsistent Landauer-B\"uttiker scattering approach.


I. INTRODUCTION
Tunneling spectroscopy was one of the early experimental methods to extract information about the energy gap as predicted by the Bardeen-Cooper-Schrieffer theory of superconductivity [1][2][3]. Following the generalization to point-contact Andreev reflection spectroscopy by Blonder-Tinkham-Klapwijk (BTK) [4], where higher order Andreev processes are taken into account for more transparent point contacts, it has also served as a tool to extract the spin-polarization of ferromagnets [5][6][7] and probe the symmetry of the order parameter of unconventional superconductors [8][9][10][11][12][13][14]. A challenge is that the complicated physics at the contact may play a crucial role. This is particularly important when it comes to applying the method to unconventional superconductors, where the superconducting order parameter may be sensitive to normal reflection at the contact. The pair breaking is accompanied by the formation of interface Andreev states, that have also been taken as fingerprints of the symmetry of the order parameter [13][14][15].
For the high-temperature superconductors with d x 2 −y 2 symmetry of the order parameter, tunneling and Andreev reflection spectroscopy typically show a large zero-bias conductance peak (ZBCP) [8][9][10][11][12][13][14]. The ZBCP appears due to the large spectral weight of zero-energy Andreev bound states at surfaces oriented 45 • relative to the crystal ab-axes, also denoted [110] surfaces. These Andreev states are formed due to the different signs of the d x 2 −y 2 order parameter around the Fermi surface [9]. In a magnetic field, these Andreev states are shifted to finite energies by the screening supercurrents, leading to a fielddependent split of the ZBCP, an effect that is well established experimentally [11,12]. A spontaneous split of the ZBCP in the absence of magnetic field has been observed in several experiments, but not all [11,[16][17][18][19][20]. Such a split indicates the possibility of a subdominant component of the order parameter of either d xy or s symmetry, combined with the main d x 2 −y 2 component in a timereversal symmetry breaking d x 2 −y 2 + id xy or d x 2 −y 2 + is state, at least at the surface [12,21,22]. Other experiments either support time-reversal symmetry breaking or give severe size constraints on the subdominant order parameter [23][24][25][26][27]. One possible explanation for the absence of the split of the ZBCP that has been proposed recently is the formation of an inhomogeneous state at the edge with spontaneous circulating currents below a relatively large transition temperature T * ≈ 0.18T c0 [28,29]. This time-reversal symmetry broken state is referred to as phase crystal [30] and is robust also in presence of strong correlations and Anderson disorder [31]. If this state is formed it also suppresses the nucleation of the subdominant component if the criticial temperature for the formation of subdominant order, T s , is too low. On the other hand, if T s is sufficiently large, it will instead prevent the phase crystal to be formed [28]. We note that the large spectral weight of Andreev surface states could also lead to other instabilities, such as spontaneous magnetization [32,33]. In conclusion, the breaking of time-reversal symmetry remains a topic of high interest in this research field.
In this paper we return to the problem of tunneling and Andreev reflection spectroscopy of d-wave superconducting surfaces including the possibility of a relatively large subdominant s-wave component preventing the formation of the phase crystal. We consider the case when the normal metal probe is sufficiently invasive that the nonequilibrium distribution function in the superconductor imposed by the voltage bias has to be computed selfconsistently with the order parameter to guarantee current conservation. The voltage dependence of the order parameter, in particular the subdominant component, If there is a subdominant s-wave component present near the surface, the peak is split due to shifting of the Andreev states to finite energies (solid blue line with circles). The sharpness of the split peak and the rapid fall down back to the pure d-wave result at |eV | = 2|eVL| ≈ 0.2kBTc0 is a result of the suppression and disappearance of the s-wave component when a non-equilibrium distribution is enforced under finite voltage bias. In a non-selfconsistent Landauer-Büttiker scattering approach, the split is not visible due to thermal broadening (green dash-dotted line), here at T = 0.1Tc0. The inset shows the quasiclassical closed loop of Andreev reflections and normal reflections at the tunnel barrier that leads to the zero-energy Andreev bound state.
has a large effect on the spectral properties as well as the transport processes (normal and Andreev reflection amplitudes). Surprisingly, the resulting split of the ZBCP in the presence of the s-wave component can be greatly enhanced, as we illustrate in Fig. 1. The blue dotted line shows the split ZBCP. This non-thermally broadened split is due to the voltage dependence of scattering amplitudes in the presence of the non-equilibrium distribution function in the contact area. As we will show in this paper, the voltage dependence stems from the suppression of the subdominant order under voltage bias. In comparison, in a Landauer-Büttiker approach, where the voltage dependence of scattering amplitudes is neglected, the shift of the Andreev states to finite energy is not visible in the conductance due to thermal broadening, see the dash-dotted green line in Fig. 1. This physics becomes relevant for geometries of the normalmetal-insulator-superconductor (NIS) contact where the traditional point-contact assumption can not be made due to the largeness of the contact, with a diameter larger than the small superconducting coherence length in the high-temperature superconductor, and relatively high transparency of the tunnel barrier. In this case, the non-equilibrium distribution has to be computed selfconsistently to guarantee current conservation. The superconducting order parameters as well as impurity self-energies then render the scattering amplitudes voltagedependent, thus influencing the conductance.
The paper is organized as follows. In section II we outline the model assumptions we make for calculating the stationary conductance-voltage characteristics including current-conservation. The details of the quasiclassical theory that we use have been summarized in Appendix A. In the following section III we go, step by step, through results for the current-conserving theory for first the pure d x 2 −y 2 case for two orientations of the order parameter relative to the tunnel barrier normal. In section III A for zero misorientation (the [100] orientation) there are no zero-energy Andreev bound states. In this case we discuss the effect of charge imbalance in a d-wave superconductor as compared to the more well studied conventional s-wave case [34][35][36][37][38]. In section III B we turn to the case of 45 • misorientation (the [110] orientation), and discuss the influence of the zero-energy Andreev bound states on charge imbalance and the influence of non-equilibrium on the ZBCP. Lastly, in section III C we include the subdominant s-wave order parameter in the energetically favorable time-reversal symmetry-breaking combination d x 2 −y 2 + is. We then study in more detail the conductance-voltage characteristics shown in Fig. 1. The paper is summarized in section IV.

II. MODEL AND METHODS
A sketch of the setup that we are considering is given in Fig. 2. A d-wave superconductor (S) of length L and width W L is connected through tunnel barriers (I) to two normal metal reservoirs (N) held at voltages ±eV /2, both at a base temperature T . We assume a thin film superconductor, thin in the perpendicular z-axis direction, with thickness t and transport from left to right along the x-axis, which is parallel to the normals of the NIS interfaces. We consider injection of carriers into the superconducting film edges in the ab-plane direction, corresponding to experiments where such high-transparency contacts between YB 2 Cu 3 O 7−δ and Au have been fabricated [39,40]. The film thickness t is assumed small compared with the penetration depth λ c in the c-axis direction (along film normal), and contacts are assumed homogeneous. The current density can then be assumed to be the same in all layers of the film and we may consider a single two-dimensional layer in the following. We consider relatively low temperatures, and assume that the in-plane penetration depth λ is large compared with the coherence length ξ 0 = v F /2πk B T c0 , where is the reduced Planck constant, k B is the Boltzmann constant, v F is the normal state Fermi velocity, and T c0 is the critical temperature. In this case for a superconductor of width W fullfilling ξ 0 W λ, screening can be neglected and the charge current can be considered to flow homogeneously as function of the transverse y coordinate. The assumption of translational invariance breaks down near the lower (y = 0) and upper edges (y = W ), . Translational invariance is assumed in the y-direction and homogeneity is assumed in the perpendicular z-direction. The left reservoir is at potential µL = eV /2 while the right reservoir is at µR = −eV /2. Both are at temperature T . The angle α specifies the crystal axis misalignment with respect to the main transport axis, x. We will limit the discussion to a symmetric system with DL = DR = D. The sideview corresponds to the experimental techniques [39,40] to contact the film in the ab-plane direction with a protecting capping layer on top.
but we assume that the contributions to the total current from these edges are small compared to the large translationally invariant flow in the interior. In fact, the restriction on system size is weaker since for the thin film of thickness t < λ the pearl length λ p = λ 2 /t λ guarantees homogeneous flow in the transverse y-direction as long as ξ 0 W λ p . We utilize the non-equilibrium quasiclassical theory of superconductivity. An overview over the relevant equations is given in Appendix A, see also Refs. [38,41] for more details. In short, the Eilenberger-Larkin-Ovchinnikov equations [42][43][44] for the Keldysh, retarded, and advanced quasiclassical Green's functions, Eq. (A2), are solved self-consistently with the superconducting order parameter ∆, impurity self-energies for arbitrary mean free path , and the local electrochemical potential φ(x). The self-consistent treatment guarantees that the charge current j is conserved from source to drain contacts, i.e., it is independent of x.
The singlet d-wave order parameter is written as where the orbital basis function η d (p F ) = √ 2 cos [2 (ϕ F − α)] depends on ϕ F , which is the angle between the Fermi momentum p F and the x-axis, and the angle α giving the orientation of the d-wave clover with respect to the x-axis, see Fig. 2. The amplitude is in non-equilibrium a spatially dependent complex quantity, ∆ d (R) = |∆ d (R)| exp[iχ(R)]. Through a local gauge transformation, one may transform the order parameter to be real, and thereby obtain the superfluid momentum p s = 2 ∇χ which signals the presence of superflow.
The superconductor is connected to the normal-metal reservoirs via insulating barriers, which are assumed to be symmetric, D L = D R = D. In this paper, we include a tunnel cone such that the interface transparency D depends on the momentum angle ϕ F . Explicitly, we use . Here, D 0 is the transparency for perpendicular incidence on the surface, and the parameter β determines the narrowness of the tunnel cone. Throughout this paper, we use β = 1 corresponding to a wide tunnel cone.
The impurity self-energies for a dilute concentration n imp of impurities is computed using a self-consistent tmatrix approximation, including multiple scattering off individual impurities but neglecting crossing diagrams. As seen in Eqs. (A14)-(A15), the model contains two parameters, the impurity concentration Γ u = n imp /πN F and the scattering phase shift, δ 0 = arctan(πN F u 0 ), assuming an isotropic impurity potential u 0 .
We consider here two limits for the potential strength. In the Born limit, u 0 is small such that only the first diagram in the t-matrix series is kept. In the second unitary limit, u 0 is considered very large and the t-matrix equation is summed to infinite order for multiple scattering off a single impurity. In both cases a single parameter Γ = Γ u sin 2 δ 0 , related to the normal state mean free path = v F /2Γ = (π/Γ)ξ 0 , characterizes the impurities. In the Born limit we assume u 0 small, but n imp large, keeping Γ constant, while in the unitary limit u 0 is large such that δ 0 → π/2. This homogeneous scattering model has been widely studied in unconventional, in particular d-wave, superconductors [38,41,[45][46][47][48][49][50][51]. In the bulk superconductor, the Born limit impurities simply broaden the density of states, while the unitary limit impurities introduce a sizeable impurity band of resonance states around the Fermi energy. This low-energy band of quasiparticles influences transport. At the surface, when there are Andreev bound states at zero energy for misorientation angle α = 0, the Born limit impurities heavily broadens them, while the unitary limit impurities are less effective in increasing their lifetime.
The φ(x)-potential, see Eq. (A28), quantifies charge imbalance, i.e., the difference in chemical potential locally at x between the condensate and the quasiparticle states. This potential is determined through the assumption of local charge neutrality [52]: excess charge due to injection of quasiparticles is compensated by a local depletion of the condensate. This also means that we neglect any charging effects: the charging energy U of the central region is small compared with the energy scale /τ d (τ d dwell time) due to escape to the leads. The charge imbalance is induced at the interfaces when quasiparticles are injected from the normal metals and decays into the interior of the superconductor. The processes determining the decay away from the interfaces are Andreev reflection and impurity scattering, which is particularly important for the d-wave order parameter with gap nodes. We find that due to the gap nodes, charge imbalance does not decay to zero for any voltage |V | > 0. This is the case both for ballistic ( > L) and diffusive devices ( < L). Experimentally, for sufficiently long devices, inelastic processes will relax charge imbalance. In our treatment, we assume that the inelastic mean free path is large both compared to elastic mean free path in and system size in L. Thus, the dwell time of non-equilibrium quasiparticles in the superconductor is considered small compared with inelastic relaxation times [38].

III. RESULTS
A. Orientation α = 0 In Fig. 3 we show an overview of the spatial dependencies of all relevant physical quantities for the orientation α = 0 and a large mean free path = 100ξ 0 with impurity scattering in the Born limit for one low and one high voltage, and for the unitary limit at high voltage only. For high transparency of the interface barriers, D 0 = 1 and wide tunnel cones (β = 1), transport in the contact regions are dominated by Andreev reflection at low voltage. This leads to a conversion of quasiparticle flow [anomalous component j a (x)] at the interfaces to a dominating superflow [local equilibrium component j le (x)] far from the contacts, see Fig. 3(a). The main qualitative difference from the conventional s-wave superconducting case [38] is the presence of the d-wave gap nodes and the more strong pair-breaking effect of elastic impurity scattering. This leads to a charge imbalance throughout The energy-like mode f1, see Eq. (A21), averaged over the full Fermi surface, f1(pF, x, ε) FS , and over half the Fermi-surface with positive momenta f1(pF, x, ε) + and negative momenta f1(pF, x, ε) − with respect to the x-axis, see also Eq. (A33). In (a) at the N-S interface (x = 0), in (b) at the center of the superconductor (x = L/2), the grey shaded area in both figures indicates the voltage window of width eV = 1.5kBTc0. In (c), we show the energy-resolved spatial evolution of f1(ε, x) FS from the surface to the center of the system. In all cases, we have D0 = 1, β = 1, α = 0, Born limit impurities with = 100ξ0, and T = 0.1Tc.
the system for all voltages |V | > 0, characterized by a potential φ(x), see Fig.3(b), as well as a residual quasiparticle flow [anomalous component j a (x)] in the interior of the superconductor, see Fig. 3(a). It was shown for conventional but anisotropic s-wave superconductors that charge imbalance may be relaxed by elastic impurity scattering [35]. For the d-wave case, we find that after an initial drop of φ(x) through Andreev processes in the contact region, the gap nodes (absence of a true energy gap) and pair-breaking elastic impurity scattering result in a residual charge imbalance extending into the interior of the superconductor for all system sizes L 200ξ 0 that we have considered. Note however, that in our setup the potential profile necessarily is antisymmetric with respect to the center of the system and φ(x) is forced to pass through zero at the center. At higher bias, the superconducting phase gradient leads to substantial Doppler shifts of the continuum states. As illustrated in the inset of Fig. 3(d), rightmoving hole-like quasiparticle states are shifted down into the bias window, while right-moving electron-like quasiparticle states are shifted up. As a consequence, an electron-to hole-like quasiparticle transmission process T he becomes available in the interface region [38,53]. This effect leads to hole-like quasiparticle transport from left to right in the figure and a negative quasiparticle current j a (x) < 0 is induced in the interior of the superconductor, see Fig. 3(e). An associated sign change in right-and left-mover quasipotentials φ + (z) and φ − (z), see Eq. (A32), also appears, see Fig. 3(g). The condition of current conservation then forces the condensate to compensate for this through larger phase gradients and a resulting larger local-equilibrium current, j le (x) > j, see Fig. 3(e). Eventually, at higher voltage, but at a voltage below the maximum of the d-wave superconducting gap, superconductivity in the mesoscale device breaks down. The breakdown is due to the combined effect of Doppler shifting superflow and the highly non-equilibrium form of the distribution function f 1 (p F , R, ε) that enters the gap equation Eq. (A26), see also Fig. 7 of Ref. [38] and the discussion in Ref. [54]. As compared with the s-wave superconducting case, the break down is at a lower voltage, due to the d-wave gap nodes. We also note that the breakdown appears at a lower current compared with the usual critical current due to pure superflow [55].
In Fig. 4 we show the distribution function f 1 (p F , x, ε) for a voltage eV = 1.5k B T c0 below the maximum of the d-wave gap. As defined in Eq. (A33), the Fermi surface average has been divided into positive and negative projections of the Fermi momenta on the transport x-axis, f 1 (p F , x, ε) ± , reflecting injection of electrons from the left and right reservoirs. The distribution is modified from its equilibrium form f eq 1 (ε) = tanh(ε/2T ) into a characteristic two-step shape at the contact, where the step width is given by the local potential, here eV /2. In the interior of the superconductor, the distribution remains highly non-equilibrium and cannot be characterized by a local effective temperature, since it does not have the shape of the tanhfunction. Instead, it can be understood in simplified terms as the result of a superposition of a superconducting and a normal component due to the d-wave gap nodes. The normal component would be a constant through the superconductor, given by Fig. 4(b)] while the superconducting component has a spatial dependence such that it relaxes back to the equilibrium form tanh(ε/2T ) away from the contacts once Andreev processes are complete. The selfconsistently computed distribution in Fig. 4 is a result of an interplay between these two components. The finite quasiparticle flow in the interior of the superconductor, the finite j a (x) in Fig. 3(a), is reflected in the difference between the distributions with positive and negative projections. The superconductivity is weakened when the distribution function f 1 (p F , x, ε) is reduced substantially in a window of subgap energies, as marked by the grey shaded area in Fig. 4. Superconductivity breaks down when the window of reduced distribution reach the superconducting coherence peaks. Let us next consider strong impurity scattering. In the case of unitary-limit scattering, there is an impurity band formed around the Fermi energy of width ∼ π∆ d Γ/2 [48,56]. This enhances charge imbalance as compared with the Born limit for the same mean free path, since the superconductor is more normal metal like. The residual potential φ(x) in the interior superconductor is therefore enhanced, see Fig. 3(j).
The moving condensate of the superconductor is associated with a superconducting phase gradient. An alternative view is to use a gauge where the order parameter is real, in which case there is a spatially dependent superfluid momentum p s (x) = 2 ∂ x χ(x). As a result, a phase difference ∆χ = χ(L) − χ(0) is established across the structure, see the right-most column of Fig. 3. For larger voltages, this phase difference can grow large and surpass 2π. Since the superconductor is rather long (L = 50ξ 0 in Fig. 3), the supercurrent is still below the bulk depairing current, p s < p sc , where p sc is the critical superfluid momentum where superconductivity breaks down due to superflow [55].
To further quantify the disequilibrium throughout the superconductor we consider the local density of states (LDOS) under voltage bias, see Fig. 5. Already in equilibrium, the d-wave DOS is modified in a wide region near the contact due the inverse proximity effect [50]. Under voltage bias, the main changes in the LDOS are due to the Doppler shifts from the moving condensate. In simplified terms, for a clean system, the superconducting coherence peak at the maximum of the d-wave gap is expected to split into two peaks, for quasiparticle states co-moving and counter-moving with the condensate flow [47]. The picture changes due to impurity scattering, which mixes the co-and counter-moving quasiparticle states present in a superclean system. For more dirty systems, for instance ∼ 10ξ 0 (not shown in the figure), impurity scattering broadens the substructure into a wide peak around the gap energy. At low energies there is an enhanced density of states within a window |ε| ∼ v F p s (x), which is also due to the Doppler shifts.
At a surface of a d-wave superconductor, with the surface normal misaligned with the crystal main axes (α = 0), Andreev bound states at zero energy, the so-called midgap states, are formed [9]. For the case of α = π/4, such zero-energy states exist for all trajectory angles ϕ F . For an interface to a normal metal, the zero-energy peak in the interface LDOS acquires a width determined by the interface transparency, ∼ ∆(ϕ F )D(ϕ F ) FS . In the presence of impurity scattering, the LDOS peak is broadened further [49]. A clear signature of the zero-energy states is a peak in the differential conductance around zero voltage in a point-contact geometry, as also found experimentally [13]. Here, we go beyond the point-contact assumption and investigate the influence of zero-energy states on charge imbalance, as well as the influence of a non-equilibrium distribution on the zero-bias conductance peak.
In Fig. 6 we show the spatial dependencies of all relevant physical quantities for the case of unitary impurity scattering with = 100ξ 0 , low-transparency interfaces D 0 = 0.2, and a voltage bias of eV = 1.0k B T c0 . In this case, there are both zero-energy Andreev bound states at the interfaces to the normal metals and an impurity band around zero energy in the interior of the superconductor. The combination enhances charge imbalance and the φpotential, see Fig. 6(b). The formation of the Andreev bound states is associated with a suppression of the order parameter near the interfaces, see inset in Fig. 7(b). This weakening of superconductivity forces the condensate to produce a large phase gradient to set up the necessary superflow, see Fig. 6(d). The resulting superfluid momentum p s (x) = p s (x)x with p s (x) = 2 ∂ x χ(x) grows large near the contacts.
The effect of superflow for α = 0 considered above is best understood in terms of Doppler shifts ε → ε−v F ·p s of continuum quasiparticles propagating along trajectories characterized by a distinct v F . For surface Andreev states this is not the case when the superflow is along the surface normal p s x as in our case, since the Andreev bound states are a superposition of partial waves with both positive and negative projections of the Fermi momentum (and Fermi velocity) with respect to the x-axis, see inset in Fig. 1. This leads to different signs of the shifts v F · p s for the partial waves with opposite momen- tum projections on the x-axis, and there is no resulting Doppler shift. Instead, the finite superflow only changes the spectral weight of these states. For an illustration of this physics, we assume a clean d-wave superconductor with a perfectly reflective interface at x = 0 and neglect the suppression of the order parameter close to the surface. For a Fermi velocity with angle ϕ F ∈ (−π/2, π/2), specular scattering at the interface connects the two velocities The order parameter with α = π/4 has opposite signs along those two direciton, see Eq. (1). The resulting surface retarded Green's function reads where and z = ε + i0 + is the energy with an infinitesimal positive imaginary part. In the absence of p s , it is straightforward to show that Eq. (3) has only one first-order pole at z = 0 with a residue of π|∆(ϕ F )|, which is the spectral weight of the Andreev bound state for this trajectory angle. For a general p s = (p x s , p y s ) T , an identical analysis shows that as long as |p x s | < |∆|, the pole is shifted to z = v y F p y s , while the residue is reduced to Thus, we see that the two components of p s have very different influence on the Andreev bound state. The component normal to the interface, p x s , reduces the spectral weight of the Andreev state. This is compensated for by an enhanced spectral weight of the continuum. In contrast, the component parallel to the surface, p y s , shifts the state in energy according to ε → ε + v y F p y s . As a result, the zero-energy states shift to positive (negative) energy for trajectories with positive (negative) projection on the y-axis, and the single peak in the density of states splits into two separate peaks. This is indeed the case in the presence of an external magnetic field along the c-axis direction perpendicular to the superconducting film, where the induced screening currents Doppler shift and split the ZBCP as found experimentally [11]. In our case, though, the voltage bias only results in superflow along the x-axis and the Andreev states stay at zero energy.
In Fig. 8 we show results for the conductance as function of applied voltage for the orientation α = π/4 and two different transparencies of the interface barriers. The conductance is computed by numerical differentiation of the current defined in Eq. (A16), and is independent of the coordinate x. For an analysis of the current and conductance it is convenient to also use the Riccati parameterization and study the resulting current formula on the normal side of for instance the left NIS interface [57]. This approach is also beneficial when we compare our results with literature where a non-selfconsistent Blonder-Tinkham-Klapwijk (BTK) type of approach has been used [4,10]. Thus lifting the result of Ref. [57], the current calculated on the normal side of the left NIS interface at x = 0 can be written as where the dependencies of all quantities within the brackets on energy ε and trajectory angle ϕ F is understood, while the voltage dependence is emphasized. Here, x 1 (x 2 ) is the distribution incoming from the left (right) side of the interface, . . . + is the average over all momenta where The remaining four terms are the physical probabilities for normal reflection (R ee ) and Andreev reflection (R he ) for electrons incoming from the left normal metal reservoir, and similarly normal transmission (T ee ) and branch-conversion transmission (T he ) for electron-like quasiparticles incoming from the superconducting device side. Eq. (6) is a generalization [57] of the well-known BTK current formula for the NIS interface [4,10]. It is a generalization since all four scattering probabilities and both distribution functions, all computed within quasiclassical theory, are voltage dependent. Within traditional BTK theory [4] and its generalization to the d-wave case [10], only the incoming distribution from the normal metal side is voltage dependent according to the reservoir and point contact assumptions The incoming distribution from the superconducting side, x 2 , as well as the four scattering probabilities in Eq. (6), are given by their equilibrium values at V = 0. The conductance within this scattering approach is then simplified to In the normal state it reduces to G N = R −1 where R N is the normal state interface resistance. We call in the following such an approach the Landauer-Büttiker scattering approach, denoted with the superscript LB.
We are now ready to compare our self-consistent stationary non-equilibrium results, blue circles in Fig. 8, to the LB scattering approach. The LB calculations include both the original non-selfconsistent approach [10], here generalized to include impurity scattering, and a selfconsistent calculation of scattering amplitudes for zero voltage where the suppression of the superconducting order parameter is taken into account, as in Ref. [57]. The self-consistent LB scattering approach (solid green lines in Fig. 8) always give a lower conductance compared to the non-selfconsistent LB scattering approach (dashed orange lines in Fig. 8) because of the gap suppression near the interface. The conductance in our nonequilibrium calculation (blue filled circles) is similar in shape to the self-consistent LB result, but they do not fully agree. Corrections to the LB scattering approach are small at small voltages and vanishes at zero voltage where a linear response calculation holds. For intermediate voltages and a high-transparency interface, where the Andreev bound states are substantially broadened, corrections are generally below five percent in the given voltage range, while for low transparency, the difference between the two approaches is more substantial and of order 20 percent due to the voltage dependence of the Andreev state spectral weight, which is also reflected in the scattering amplitudes. At high voltage, where the d-wave order parameter becomes more substantially affected by the non-equilibrium distributions, the corrections grow, until relatively abruptly for these system parameters superconductivity disappears.

C. Subdominant order-parameter component
Even in the absence of an external magnetic field a splitting of the ZBCP has been observed in several experiments at very low temperatures. One mechanism proposed and extensively discussed in the literature [11,19,20] is the presence of a subdominant order parameter of d xy or s symmetry, in addition to the d x 2 −y 2 order parameter, in a combination of the form d x 2 −y 2 + id xy or d x 2 −y 2 + is. Such combinations result in a superconducting state with broken time-reversal symmetry close to the surface of the superconductor. As a result, the subdominant component becomes enhanced close to the surface and induces a splitting of the ZBCP to finite voltages [12]. Since the d xy subdominant component is very sensitive to impurity scattering, we consider in the fol- lowing a subdominant s-wave order parameter with coupling constant corresponding to a critical temperature of T s = 0.2T c0 . This subdominant component is zero in the bulk, ∆ s = 0, but reaches values of ∆ s ∼ 0.2k B T c0 close to the surfaces, see the equilibrium shape ∆ s (x) (solid blue line) in Fig. 9. Under voltage bias, the injected nonequilibrium distribution reduces both components of the order parameter. Since the subdominant component is smaller to begin with, it becomes more strongly suppressed than the dominant component already for small voltages. Fig. 9 shows the suppression of the magnitude of the subdominant component for increasing bias voltage, as computed self-consistently through Eq. (A26). At voltages above eV ≈ 0.2k B T c0 the subdominant component is fully suppressed, as seen in the inset in Fig. 9. As a result, the Andreev bound states move back towards the Fermi energy for increasing voltage, see Fig. 10. This shift and the restoration of the zero-energy peak in the surface LDOS with increasing bias voltage also affects the differential conductance.
As seen in Fig. 11, blue solid line with filled circles, the conductance peak is shifted away from zero voltage due to the equilibrium shift of the Andreev states to finite energy. Once the subdominant component is suppressed at a voltage eV ≈ 0.2k B T c0 , the conductance falls back to that of a pure d-wave superconductor without any subdominant order parameter component (dashed orange line). As comparison we have included in Fig. 11 the results of a LB calculation (green dotted line). As seen, the ZBCP is not split and the voltage dependence is smooth and similar to the case without s-wave component. This is due to thermal broadening together with the relatively large broadening of the Andreev states by impurity scattering and the coupling to the normal metal through the barrier with transparency D 0 = 0.5. A weak split of the ZBCP appears at very low temperature T = 0.01T c0 (not  11. The low-voltage differential conductance, normalized to the normal-state value, for a pure d-wave superconductor (dashed orange line) and for the case of an s-wave subdominant order parameter with Ts = 0.2Tc0 (solid blue line). Once eV = 2eVL → 0.2kBTc0 from below, the subdominant component ∆s vanishes, see also the inset of Fig. 9. The parameters are the same as in Fig. 9. As comparison we also show results from a Landauer-Büttiker scattering-approach (LB) calculation.
shown in the figure), but it is much weaker than our nonequilibrium result at T = 0.1T c0 (blue curve). Note that the self-consistent LB approach and our non-equilibrium results coincide at V = 0 where the linear response approximation holds and all quantities except the voltage perturbation of x 1 in Eq. (6) take their V = 0 equilibrium values.
The main new ingredients in our fully selfconsistent calculation are that also the right-side distribution x 2 , as well as all scattering probabilities, depend on the applied voltage. In contrast to the case without subdominant component considered above (Fig. 8), the correction to the LB approach is qualitatively important and result in a dramatic enhancement of the ZBCP split and a rapid drop in the conductance as the s-wave component is suppressed by the non-equilibrium distributions. Thermal broadening is not important near the rapid drop as is clearly seen in Fig. 11. By using the formula in Eq. (6), we find that the dominant corrections to the LB scattering approach are those due to selfconsistent changes in the amplitudes collected into the factor [1 − R ee (V ) + R he (V )] in Eq. (6). The applied bias voltage leads to a quick suppression of the subdominant order parameter component ∆ s through the highly nonequilibrium forms of the distribution functions. As a consequence, the surface Andreev bound states move back to zero energy, as seen in Fig. 10. The resonances in the scattering amplitudes mirror this spectral rearrangement. As a result the term proportional to the voltage derivative d[1 − R ee (V ) + R he (V )]/dV also gives a large con-tribution of roughly 20 percent of the total conductance. In contrast, the backflow contribution proportional to x 2 in Eq. (6) leads to small corrections of less than five percent. At a critical voltage ∼ 0.2k B T c0 , the subdominant component is fully suppressed and a rapid drop of the conductance back to the pure d-wave curve, the orange dashed line in Fig. 11), is visible. We conclude that the voltage dependence of scattering probabilities can be of great importance when the non-equilibrium distributions couple back through the order parameter self-consistency.

IV. SUMMARY
In summary we have presented a comprehensive study of the stationary non-equilibrium response of a d-wave superconductor coupled to two normal-metal reservoirs under a voltage bias, taking into account different orientations of the order parameter relative to the interfaces to the reservoirs, formation of interface zero-energy Andreev bound states, scalar impurity scattering, and a possible s-wave subdominant component of the order parameter. In all cases we ensure current conservation by computing all self-energies selfconsistently taking into account the non-equilibrium distribution functions. For the case of a pure d-wave order parameter, we have found that charge imbalance extends into the bulk of the superconductor due to the presence of the nodes of the d-wave order parameter and the pair breaking effects of scalar impurities. Charge imbalance is enhanced for orientations where zero-energy interface Andreev bound states are formed and it is also enhanced for unitary limit impurity scattering which induce a low-energy band of quasiparticles states. The non-equilibrium distribution function induced in the superconductor suppresses the dwave order parameter and superconductivity disappears for voltages approaching the gap voltage. Nevertheless, despite the induced non-equilibrium state, the resulting conductance-voltage dependencies for the case of pure d x 2 −y 2 superconducting order resembles the results of a non-selfconsistent Landauer-Büttiker approach, if the original approach [10] is corrected for the zero-voltage (V = 0) equilibrium suppression of the order parameter near the interfaces and the broadening effect of impurities. This holds until superconductivity is suppressed at voltages approaching the gap voltage.
When we allow for a subdominant component with s-wave symmetry, forming the time-reversal symmetry breaking combination d x 2 −y 2 + is near the interfaces, the effects of the non-equilibrium distribution is much more severe. In this case, the zero-bias conductance peak split is dramatically enhanced compared with the nonselfconsistent Landauer-Büttiker approach. This effect is due to the suppression of the subdominant s-wave order parameter when the voltage is applied. This leads to spectral rearrangements and corrections to the scattering amplitudes. As a result, a non-thermally broadened split of the zero-bias conductance peak appears.
The Landauer-Büttiker approach holds for pointcontact geometries, where the contact radius is smaller than the coherence length ξ 0 . Since ξ 0 is small in hightemperature superconductors such contacts are challenging to fabricate. Our approach is applicable to the complementary geometry with very wide and transparent contacts where the induced non-equilibrium distribution is important. Such contacts to superconducting films are experimentally feasible [39,40] and could serve as an interesting probe of the rich surface physics of unconventional superconductors, including the possibility of subdominant order parameters and time-reversal symmetry breaking.

ACKNOWLEDGMENTS
We thank F. Lombardi for valuable discussion. We acknowledge the Swedish research council for financial support. The computations were enabled by resources provided by the Swedish National Infrastructure for Computing (SNIC) at NSC partially funded by the Swedish Research Council through grant agreement no. 2018-05973.

Appendix A: Quasiclassical theory
Here, we review the key set of equations of the quasiclassical theory of superconductivity that we base our calculations on. For further details we refer to Refs. [38,41].
In the stationary non-equilibrium case under consideration we need to calculate the retarded and advanced Green's functions as well as the Keldysh component. They are organized into a matrix in Keldysh space aš For brevity we often drop the explicit references to the dependencies on momentum direction p F , coordinate R, and energy ε. All three elements in Eq. (A1) are matrices in Nambu space, as indicated by theˆ. In the steady-state,ǧ is time-independent and obeys the socalled Eilenberger equation, where [Ǎ,B] denotes a commutator between matricesǍ andB, as well as the normalization conditionǧ 2 = −π 21 . These two equations were first derived by Eilenberger [42], and separately Larkin and Ovchinnikov [43], and the generalization to nonequilibrium is due to Eliashberg [44]. The quasiclassical Hamiltonianȟ has the same Keldysh structure as the propagator in Eq. (A1). The components in Nambu space arê These self-energies are defined below. The elements ofǧ can be parametrized using coherence amplitudes γ,γ [58][59][60] and generalized distribution functions x,x [61,62]. Using the definitions of G R ≡ (1 − γ RγR ) −1 and F R ≡ G R γ R , we have for the retarded element and an analogous expression for the advanced function g A . The Keldysh component readŝ An important symmetry in the theory is particle-hole conjugation, which can be expressed as A set of coupled transport equations for the parametrizing functions can be derived from Eq. (A2). The Riccati equation for the coherence amplitude reads while the equation for the distribution function x is Application of Eq. (A7) to these two equations gives the corresponding equations forγ R,A andx. Solutions to both transport equations are found by propagating from a start point to an end point on a trajectory determined by the Fermi momentum v F , in our case fully determined by the momentum orientation angle ϕ F . At interfaces between the superconductor and the normal-metal reservoirs, both the coherence and distribution functions have to be connected by boundary conditions. The latter are derived from the scattering matrices for an insulating barrier between a normal metal and a superconductor, see Refs. 61-63.
Having obtained solutions to Eqs. (A8)-(A9), we can update all self-energies. They consist of the mean field order parameter and scalar impurity self energy aš The Keldysh matrix structure of the mean field order parameter is simple,ȟ mf =∆1, while the Nambu structure is∆ = (∆)τ 1 − (∆)τ 2 . The order parameter consists of a sum over the d-wave and subdominant s-wave components. Each component satisfies a gap equation and η s (p F ) = 1. The pairing interactions λ Γ are eliminated in favor of the clean-limit critical temperatures for each channel by the replacement In the main text, since the d-wave is the dominant component, we let T c0 = T cd and use k B T c0 as the natural energy scale. The critical temperature of the subdominant s-wave is for short denoted T s < T c0 . The impurity self-energies are calculated in the noncrossing approximation for the t-matrix equation [45], discussed in detail in Ref. 41. The matrix structure iš For scattering that is isotropic in momentum space with an s-wave scattering potential u 0 the elements ofť satisfy the equationst This procedure of solving Eqs. (A8)-(A9) and updating the selfenergies is iterated until a fully self-consistent solution is found and current is conserved throughout the structure. In the present paper, we allow for a maximum relative error |j(x) − j I |/|j I | < 5 · 10 −3 , so the current everywhere deviates less than half a percent from j I , the current at the interface to the normal-metal reservoirs. The charge current is determined by the Keldysh componentĝ K , Tr as is the local electrochemical potential The retarded componentĝ R determines the normalized, momentum-resolved local density of states per spin N (p F , R, ε) = − 1 4π Im Tr τ 3ĝ R (p F , R, ε) , and thus the full density of states The Keldysh componentĝ K can, in addition to Eq. (A6), also be parametrized aŝ where the distribution matrixf readŝ The non-equilibrium electron distribution function is then given by which becomes a Fermi-Dirac distribution function in equilibrium. We refer to f 1 (f 3 ) as the energy-like (charge-like) mode, connecting to the nomenclature established for diffusive superconductors [64]. After Fermisurface averaging, they have different symmetries as function of energy, namely which corresponds to the symmetries obtained in the diffusive case. The function h can be related to x by The equation of motion for x andx is easier to solve, while the final results can be easier to interpret in terms of h, or f e/1/3 . As an example, we can write the gap equations in Eq. (A11) as and the charge current, Eq. (A16), as The natural unit for charge current is then j 0 = ev F N F k B T c0 . The electrochemical potential, Eq. (A17), can similarly be written as The distribution function h can be split into a localequilibrium function h le , and an anomalous function h a , where h le is given by The splitting introduced in Eq. (A29) directly translate to an analogous splitting of the function x by using Eq. (A25), and for the two modes f 1 and f 3 by using Eq. (A21). As a consequence observables such as the charge current can be similarly split, j(R) = j le (R) + j a (R).
Since h le , and hence f le 1,3 , are independent of p F , Eq. (A27) shows that the local-equilibrium current j le can only be non-zero if there is an asymmetry in the angleresolved spectrum N (p F , R, ε). In our case, this asymmetry is created by the self-consistently determined superflow p s in the superconductor. The anomalous current j a , due to the distribution h a , then contains all current flow due to differences in occupation of left-moving and right-moving quasiparticles. A measure of this difference in occupation are the right-mover (left-mover) quasipotential φ + (φ − ), defined as where X a is obtained by using x a in Eq. (A6). In Eq. (A32), the . . . ± denotes a partial Fermi-surface average where the Heaviside step function Θ(± cos ϕ F ) gives unity if v x F , as defined in Eq. (2), is positive (+) or negative (-), and zero otherwise. In the normal state, the left-mover and right-mover quasipotentials satisfy φ = (φ + + φ − )/2, which connects to the concepts used in normal-state ballistic systems [65].