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

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–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–7] and probe the symmetry of the order parameter of unconventional superconductors [8–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–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–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 field-dependent 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–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 time-reversal 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–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\,^{\ast}\approx 0.18T_\mathrm{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_\mathrm{s}$, is too low. On the other hand, if $T_\mathrm{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 non-equilibrium distribution function in the superconductor imposed by the voltage bias has to be computed self-consistently with the order parameter to guarantee current conservation. The voltage dependence of the order parameter, in particular the subdominant component, 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 figure 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 figure 1. This physics becomes relevant for geometries of the normal-metal–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 self-consistently to guarantee current conservation. The superconducting order parameters as well as impurity self-energies then render the scattering amplitudes voltage-dependent, thus influencing the conductance.

Zoom In Zoom Out Reset image size

The paper is organized as follows. In section 2 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. In the following section 3 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 3.1 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–38]. In section 3.2 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 3.3 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 figure 1. The paper is summarized in section 4.

A sketch of the setup that we are considering is given in figure 2. A d-wave superconductor (S) of length L and width $W\gg L$ is connected through tunnel barriers (I) to two normal metal reservoirs (N) held at voltages $\pm 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₂Cu₃O$_{7-\delta}$ and Au have been fabricated [39, 40]. The film thickness t is assumed small compared with the penetration depth $\lambda_\mathrm{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 $\xi_0 = \hbar v_\mathrm{F}/2{\pi}k_\mathrm{B}T_\mathrm{c0}$, where $\hbar$ is the reduced Planck constant, $k_\mathrm{B}$ is the Boltzmann constant, $v_\mathrm F$ is the normal state Fermi velocity, and $T_\mathrm{c0}$ is the critical temperature. In this case for a superconductor of width W fullfilling $\xi_0\ll W\ll\lambda$, 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), 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\lt\lambda$ the pearl length $\lambda_\mathrm{p} = \lambda^2/t\gg\lambda$ guarantees homogeneous flow in the transverse y-direction as long as $\xi_0\ll W\ll\lambda_\mathrm{p}$.

Zoom In Zoom Out Reset image size

We utilize the non-equilibrium quasiclassical theory of superconductivity. An overview over the relevant equations is given in appendix, see also [38, 41] for more details. In short, the Eilenberger–Larkin–Ovchinnikov equations [42–44] for the Keldysh, retarded, and advanced quasiclassical Green's functions, equation (A2), are solved self-consistently with the superconducting order parameter Δ, impurity self-energies for arbitrary mean free path $\ell$, and the local electrochemical potential $\phi(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

$$\begin{equation} \Delta(\mathbf{p}_\mathrm{F},\mathbf{R}) = \Delta_\mathrm{d}(\mathbf{R})\eta_\mathrm{d}(\mathbf{p}_\mathrm{F}), \end{equation} \tag{ 1 }$$

where the orbital basis function $\eta_\mathrm{d}(\mathbf{p}_\mathrm{F}) = \sqrt{2}\cos\left[2\left(\varphi_\mathrm{F}-\alpha\right)\right]$ depends on $\varphi_\mathrm{F}$, which is the angle between the Fermi momentum $\mathbf{p}_\mathrm{F}$ and the x-axis, and the angle α giving the orientation of the d-wave clover with respect to the x-axis, see figure 2. The amplitude is in non-equilibrium a spatially dependent complex quantity, $\Delta_\mathrm{d}(\mathbf{R}) = |\Delta_\mathrm{d}(\mathbf{R})|\exp[i\chi(\mathbf{R})]$. Through a local gauge transformation, one may transform the order parameter to be real, and thereby obtain the superfluid momentum $\mathbf{p}_\mathrm{s} = \tfrac{\hbar}{2}\nabla\chi$ 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_\mathrm{L} = D_\mathrm{R} = D$. In this paper, we include a tunnel cone such that the interface transparency D depends on the momentum angle $\varphi_\mathrm{F}$. Explicitly, we use $D(\varphi_\mathrm{F}) = D_0(e^{-\beta \sin^2~\varphi_\mathrm{F}} - e^{-\beta})/(1 - e^{-\beta})$. Here, D₀ 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_\mathrm{imp}$ of impurities is computed using a self-consistent t-matrix approximation, including multiple scattering off individual impurities but neglecting crossing diagrams. As seen in equations (A14) and (A15), the model contains two parameters, the impurity concentration $\Gamma_u = n_\mathrm{imp}/\pi\mathcal{N}_\mathrm{F}$ and the scattering phase shift, $\delta_0 = \arctan(\pi \mathcal{N}_\mathrm{F} u_0)$, assuming an isotropic impurity potential u₀.

We consider here two limits for the potential strength. In the Born limit, u₀ is small such that only the first diagram in the t-matrix series is kept. In the second unitary limit, u₀ 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 $\Gamma = \Gamma_u \sin^2\delta_0$, related to the normal state mean free path $\ell = \hbar v_\mathrm{F}/2\Gamma = (\pi/\Gamma)\xi_0$, characterizes the impurities. In the Born limit we assume u₀ small, but $n_\mathrm{imp}$ large, keeping Γ constant, while in the unitary limit u₀ is large such that $\delta_0\rightarrow\pi/2$. This homogeneous scattering model has been widely studied in unconventional, in particular d-wave, superconductors [38, 41, 45–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 $\phi(x)$-potential, see equation (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 $\hbar/\tau_\mathrm{d}$ ($\tau_\mathrm{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|\gt0$. This is the case both for ballistic ($\ell\gt L$) and diffusive devices ($\ell\lt 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 $\ell_\mathrm{in}\gg\ell$ and system size $\ell_\mathrm{in}\gg L$. Thus, the dwell time of non-equilibrium quasiparticles in the superconductor is considered small compared with inelastic relaxation times [38].

3.1. Orientation α = 0

In figure 3 we show an overview of the spatial dependencies of all relevant physical quantities for the orientation α = 0 and a large mean free path $\ell = 100\xi_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^\mathrm{a}(x)$] at the interfaces to a dominating superflow [local equilibrium component $j^\mathrm{le}(x)$] far from the contacts, see figure 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 system for all voltages $|V|\gt0$, characterized by a potential $\phi(x)$, see figure 3(b), as well as a residual quasiparticle flow [anomalous component $j^\mathrm{a}(x)$] in the interior of the superconductor, see figure 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 $\phi(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\lesssim 200\xi_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 $\phi(x)$ is forced to pass through zero at the center.

Zoom In Zoom Out Reset image size

At higher bias, the superconducting phase gradient leads to substantial Doppler shifts of the continuum states. As illustrated in the inset of figure 3(d), right-moving 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_\mathrm{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^\mathrm{a}(x)\lt0$ is induced in the interior of the superconductor, see figure 3(e). An associated sign change in right- and left-mover quasipotentials $\phi_{+}(z)$ and $\phi_{-}(z)$, see equation (A32), also appears, see figure 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^\mathrm{le}(x)\gt j$, see figure 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(\mathbf{p}_\mathrm{F},\mathbf{R},\varepsilon)$ that enters the gap equation (A26), see also figure 7 of [38] and the discussion in [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 figure 4 we show the distribution function $f_1(\mathbf{p}_\mathrm{F},x,\varepsilon)$ for a voltage $eV = 1.5k_\mathrm{B}T_\mathrm{c0}$ below the maximum of the d-wave gap. As defined in equation (A33), the Fermi surface average has been divided into positive and negative projections of the Fermi momenta on the transport x-axis, $\left\langle\, f_1(\mathbf{p}_\mathrm{F},x,\varepsilon) \right\rangle_{\pm}$, reflecting injection of electrons from the left and right reservoirs. The distribution is modified from its equilibrium form $f\;_1^\mathrm{eq}(\varepsilon) = \tanh(\varepsilon/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 tanh-function. 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 $f\,_1^\mathrm{N}(\varepsilon) = \tfrac{1}{2}\left( \tanh[(\varepsilon-eV/2)/(2T)] +\tanh[(\varepsilon+eV/2)/2T)] \right)$, (red dotted line in figure 4(b)) while the superconducting component has a spatial dependence such that it relaxes back to the equilibrium form $\tanh(\varepsilon/2T)$ away from the contacts once Andreev processes are complete. The self-consistently computed distribution in figure 4 is a result of an interplay between these two components. The finite quasiparticle flow in the interior of the superconductor, the finite $j^\mathrm{a}(x)$ in figure 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(\mathbf{p}_\mathrm{F},x,\varepsilon)$ is reduced substantially in a window of subgap energies, as marked by the grey shaded area in figure 4. Superconductivity breaks down when the window of reduced distribution reach the superconducting coherence peaks.

Zoom In Zoom Out Reset image size

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 ${\sim}\sqrt{\pi\Delta_\mathrm{d}\Gamma/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 $\phi(x)$ in the interior superconductor is therefore enhanced, see figure 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_\mathrm{s}(x) = \tfrac{\hbar}{2}\partial_x\chi(x)$. As a result, a phase difference $\Delta\chi = \chi(L)-\chi(0)$ is established across the structure, see the right-most column of figure 3. For larger voltages, this phase difference can grow large and surpass 2π. Since the superconductor is rather long ($L = 50\xi_0$ in figure 3), the supercurrent is still below the bulk depairing current, $p_\mathrm{s}\lt p_\mathrm{sc}$, where $p_\mathrm{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 figure 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 $\ell\sim 10\xi_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 $|\varepsilon|\sim v_\mathrm{F}p_\mathrm{s}(x)$, which is also due to the Doppler shifts.

Zoom In Zoom Out Reset image size

3.2. Orientation $\alpha = \pi/$4

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 $\alpha = \pi/4$, such zero-energy states exist for all trajectory angles $\varphi_\mathrm{F}$. For an interface to a normal metal, the zero-energy peak in the interface LDOS acquires a width determined by the interface transparency, ${\sim}\left\langle \Delta(\varphi_\mathrm{F})D(\varphi_\mathrm{F}) \right\rangle_{\scriptscriptstyle\mathrm{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 ZBCP.

In figure 6 we show the spatial dependencies of all relevant physical quantities for the case of unitary impurity scattering with $\ell = 100\xi_0$, low-transparency interfaces $D_0 = 0.2$, and a voltage bias of $eV = 1.0k_\mathrm{B}T_\mathrm{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 figure 6(b). The formation of the Andreev bound states is associated with a suppression of the order parameter near the interfaces, see inset in figure 7(b). This weakening of superconductivity forces the condensate to produce a large phase gradient to set up the necessary superflow, see figure 6(d). The resulting superfluid momentum $\mathbf{p}_\mathrm{s}(x) = p_\mathrm{s}(x){\hat x}$ with $p_\mathrm{s}(x) = \tfrac{\hbar}{2}\partial_x\chi(x)$ grows large near the contacts.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The effect of superflow for α = 0 considered above is best understood in terms of Doppler shifts $\varepsilon\rightarrow\varepsilon-\mathbf{v}_\mathrm{F}\cdot\mathbf{p}_\mathrm{s}$ of continuum quasiparticles propagating along trajectories characterized by a distinct $\mathbf{v}_\mathrm{F}$. For surface Andreev states this is not the case when the superflow is along the surface normal $\mathbf{p}_\mathrm{s}\parallel{\hat 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 figure 1. This leads to different signs of the shifts $\mathbf{v}_\mathrm{F}\cdot\mathbf{p}_\mathrm{s}$ for the partial waves with opposite momentum 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 $\varphi_\mathrm{F} \in (-\pi/2, \pi/2)$, specular scattering at the interface connects the two velocities

$$\begin{align} \mathbf{v}_\mathrm{F}^\mathrm{out}\! = \! \begin{pmatrix} v\,_\mathrm{F}^x\\[4pt] v_\mathrm{F}^y \end{pmatrix} \equiv v_\mathrm{F}\begin{pmatrix} \cos \varphi_\mathrm{F}\\[4pt] \sin \varphi_\mathrm{F} \end{pmatrix} , \mathbf{v}_\mathrm{F}^\mathrm{in} = v_\mathrm{F} \begin{pmatrix} -\cos \varphi_\mathrm{F}\\[4pt] \sin \varphi_\mathrm{F}\\ \end{pmatrix}. \end{align} \tag{ 2 }$$

The order parameter with $\alpha = \pi/4$ has opposite signs along those two direciton, see equation (1). The resulting surface retarded Green's function reads

$$\begin{align} \frac{g^\mathrm{R}(z,\varphi_\mathrm{F})}{-i\pi} = \frac{(\mathbf{v}_\mathrm{F}^{\mathrm{in}} - \mathbf{v}_\mathrm{F}^{\mathrm{out}}) \cdot \mathbf{p}_\mathrm{s} + i \left[\Omega(\mathbf{v}_\mathrm{F}^{\mathrm{out}}) + \Omega(\mathbf{v}_\mathrm{F}^{\mathrm{in}}) \right]}{2z - (\mathbf{v}_\mathrm{F}^{\mathrm{in}} + \mathbf{v}_\mathrm{F}^{\mathrm{out}}) \cdot \mathbf{p}_\mathrm{s} + i \left[\Omega(\mathbf{v}_\mathrm{F}^{\mathrm{out}}) - \Omega(\mathbf{v}_\mathrm{F}^{\mathrm{in}}) \right]}, \end{align} \tag{ 3 }$$

where

$$\begin{align} \Omega(\mathbf{v}_\mathrm{F}) \equiv \sqrt{\Delta^2(\mathbf{p}_\mathrm{F}) - (z - \mathbf{v}_\mathrm{F} \cdot \mathbf{p}_\mathrm{s} )^2 }, \end{align} \tag{ 4 }$$

and $z = \varepsilon+i0^+$ is the energy with an infinitesimal positive imaginary part. In the absence of $\mathbf{p}_\mathrm{s}$, it is straightforward to show that equation (3) has only one first-order pole at z = 0 with a residue of $\pi |\Delta(\varphi_\mathrm{F})|$, which is the spectral weight of the Andreev bound state for this trajectory angle. For a general $\mathbf{p}_\mathrm{s} = (p_\mathrm{s}^x,p_\mathrm{s}^y)^\mathrm{T}$, an identical analysis shows that as long as $|p_\mathrm{s}^x| \lt |\Delta|$, the pole is shifted to $z = v_\mathrm{F}^y p_\mathrm{s}^y$, while the residue is reduced to

$$\begin{align} \pi\sqrt{\Delta^2(\mathbf{p}_\mathrm{F}) - (v\,_\mathrm{F}^x p_\mathrm{s}^x )^2 }. \end{align} \tag{ 5 }$$

Thus, we see that the two components of $\mathbf{p}_\mathrm{s}$ have very different influence on the Andreev bound state. The component normal to the interface, $p_\mathrm{s}^x$, 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_\mathrm{s}^y$, shifts the state in energy according to $\varepsilon \rightarrow \varepsilon + v_\mathrm{F}^y p_\mathrm{s}^y$. 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 figure 8 we show results for the conductance as function of applied voltage for the orientation $\alpha = \pi/4$ and two different transparencies of the interface barriers. The conductance is computed by numerical differentiation of the current defined in equation (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 BTK type of approach has been used [4, 10]. Thus lifting the result of [57], the current calculated on the normal side of the left NIS interface at x = 0 can be written as

$$\begin{align} j_x(V) & = 2e\mathcal{N}_\mathrm{F} \int\limits_{-\varepsilon_\mathrm{c}}^{\varepsilon_\mathrm{c}} \mathrm{d}\varepsilon \bigl\langle v\,_\mathrm{F}^x \bigl[\textsf{x}_1( V) \left\{1 - R_\mathrm{ee} (V) + R_\mathrm{he} (V) \right\}\nonumber\\ & \quad + \textsf{x}_2( V)\left\{\bar{T}_\mathrm{he} (V) - \bar{T}_\mathrm{ee} (V) \right\}\bigr]\bigr\rangle_{+}, \end{align} \tag{ 6 }$$

where the dependencies of all quantities within the brackets on energy ε and trajectory angle $\varphi_\mathrm{F}$ is understood, while the voltage dependence is emphasized. Here, $\textsf{x}_1\,\,(\textsf{x}_2)$ is the distribution incoming from the left (right) side of the interface, $\langle \dots\rangle_{+}$ is the average over all momenta where $v\,_\mathrm{F}^x \equiv v_\mathrm{F} \cos \varphi_\mathrm{F} \gt 0$. The remaining four terms are the physical probabilities for normal reflection ($R_\mathrm{ee}$) and Andreev reflection ($R_\mathrm{he}$) for electrons incoming from the left normal metal reservoir, and similarly normal transmission ($\bar{T}_\mathrm{ee}$) and branch-conversion transmission ($\bar{T}_\mathrm{he}$) for electron-like quasiparticles incoming from the superconducting device side. Equation (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 $\textsf{x}_1 = \tanh[(\varepsilon-eV)/2k_\mathrm{B}T_\mathrm{c0}]$. The incoming distribution from the superconducting side, $\textsf{x}_2$, as well as the four scattering probabilities in equation (6), are given by their equilibrium values at V = 0. The conductance within this scattering approach is then simplified to

$$\begin{align} G^\mathrm{LB}(V) = 2e\mathcal{N}_\mathrm{F} \int\limits_{-\varepsilon_\mathrm{c}}^{\varepsilon_\mathrm{c}} \!&\mathrm{d}\varepsilon \left\langle v\,_\mathrm{F}^x \frac{\mathrm{d}\textsf{x}_1(V)}{\mathrm{d}V} \left[ 1 - R_\mathrm{ee}(0) + R_\mathrm{he}(0) \right] \right\rangle_+. \end{align} \tag{ 7 }$$

In the normal state it reduces to $G_\mathrm{N} = R_\mathrm{N}^{-1} = 2e^2\mathcal{N}_\mathrm{F} v_\mathrm{F}$ $\left\langle D(\varphi_\mathrm{F})\cos\varphi_\mathrm{F} \right\rangle_{+}$, where $R_\mathrm{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.

Zoom In Zoom Out Reset image size

We are now ready to compare our self-consistent stationary non-equilibrium results, blue circles in figure 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 self-consistent calculation of scattering amplitudes for zero voltage where the suppression of the superconducting order parameter is taken into account, as in [57]. The self-consistent LB scattering approach (solid green lines in figure 8) always give a lower conductance compared to the non-selfconsistent LB scattering approach (dashed orange lines in figure 8) because of the gap suppression near the interface. The conductance in our non-equilibrium 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.

3.3. 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 following a subdominant s-wave order parameter with coupling constant corresponding to a critical temperature of $T_\mathrm{s} = 0.2T_\mathrm{c0}$. This subdominant component is zero in the bulk, $\Delta_\mathrm{s} = 0$, but reaches values of $\Delta_\mathrm{s} \sim 0.2k_\mathrm{B} T_\mathrm{c0}$ close to the surfaces, see the equilibrium shape $\Delta_\mathrm s (x)$ (solid blue line) in figure 9.

Zoom In Zoom Out Reset image size

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. Figure 9 shows the suppression of the magnitude of the subdominant component for increasing bias voltage, as computed self-consistently through equation (A26). At voltages above $eV \approx 0.2k_\mathrm{B} T_\mathrm{c0}$ the subdominant component is fully suppressed, as seen in the inset in figure 9. As a result, the Andreev bound states move back towards the Fermi energy for increasing voltage, see figure 10. This shift and the restoration of the zero-energy peak in the surface LDOS with increasing bias voltage also affects the differential conductance.

Zoom In Zoom Out Reset image size

As seen in figure 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 \approx 0.2k_\mathrm{B} T_\mathrm{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 figure 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_\mathrm{c0}$ (not shown in the figure), but it is much weaker than our non-equilibrium result at $T = 0.1T_\mathrm{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 $\textsf{x}_1$ in equation (6) take their V = 0 equilibrium values.

Zoom In Zoom Out Reset image size

The main new ingredients in our fully selfconsistent calculation are that also the right-side distribution $\textsf{x}_2$, as well as all scattering probabilities, depend on the applied voltage. In contrast to the case without subdominant component considered above (figure 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 figure 11. By using the formula in equation (6), we find that the dominant corrections to the LB scattering approach are those due to self-consistent changes in the amplitudes collected into the factor $[1-R_\mathrm{ee}(V) + R_\mathrm{he}(V)]$ in equation (6). The applied bias voltage leads to a quick suppression of the subdominant order parameter component $\Delta_\mathrm{s}$ through the highly non-equilibrium forms of the distribution functions. As a consequence, the surface Andreev bound states move back to zero energy, as seen in figure 10. The resonances in the scattering amplitudes mirror this spectral rearrangement. As a result the term proportional to the voltage derivative $\mathrm{d}[1 - R_\mathrm{ee}(V) + R_\mathrm{he}(V)]/\mathrm{d}V$ also gives a large contribution of roughly 20 percent of the total conductance. In contrast, the backflow contribution proportional to x₂ in equation (6) leads to small corrections of less than five percent. At a critical voltage ${\sim}0.2k_\mathrm{B}T_\mathrm{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 figure 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.

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 d-wave 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 ZBCP split is dramatically enhanced compared with the non-selfconsistent 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 ZBCP appears.

The Landauer–Büttiker approach holds for point-contact geometries, where the contact radius is smaller than the coherence length ξ₀. Since ξ₀ is small in high-temperature 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.

The data that support the findings of this study are available upon reasonable request from the authors.

We thank F Lombardi for valuable discussions. 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.

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 [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 as

$$\begin{equation} \check{g}(\mathbf{p}_\mathrm{F},\mathbf{R},\varepsilon) = \begin{pmatrix} \hat{g}^\mathrm{R}(\mathbf{p}_\mathrm{F},\mathbf{R},\varepsilon) & \hat{g}^\mathrm{K}(\mathbf{p}_\mathrm{F},\mathbf{R},\varepsilon) \\ 0 & \hat{g}^\mathrm{A}(\mathbf{p}_\mathrm{F},\mathbf{R},\varepsilon) \end{pmatrix}. \end{equation} \tag{ A1 }$$

For brevity we often drop the explicit references to the dependencies on momentum direction $\mathbf{p}_\mathrm{F}$, coordinate R, and energy ɛ. All three elements in equation (A1) are matrices in Nambu space, as indicated by the $\hat{~}$. In the steady-state, $\check{g}$ is time-independent and obeys the so-called Eilenberger equation,

$$\begin{equation} i \hbar \mathbf{v}_\mathrm{F} \cdot \nabla \check{g} + \left[ \varepsilon \hat{\tau}_3\check{1} - \check{h}, \check{g} \right] = 0, \end{equation} \tag{ A2 }$$

where $[\check A,\check B]$ denotes a commutator between matrices $\check A$ and $\check B$, as well as the normalization condition $\check{g}^2 = - \pi^2~\check{1}$. 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 $\check{h}$ has the same Keldysh structure as the propagator in equation (A1). The components in Nambu space are

$$\begin{align} \hat{h}_\mathrm{}^\mathrm{R,A}(\mathbf{R},\varepsilon) = \begin{pmatrix} \Sigma^\mathrm{R,A}(\mathbf{R},\varepsilon) & \Delta^\mathrm{R,A}(\mathbf{R},\varepsilon) \\ \tilde\Delta^\mathrm{R,A}(\mathbf{R},\varepsilon) & \tilde\Sigma^\mathrm{R,A}(\mathbf{R},\varepsilon) \end{pmatrix}, \end{align} \tag{ A3 }$$

and

$$\begin{align} \hat{h}^\mathrm{K}(\mathbf{R},\varepsilon) \equiv \begin{pmatrix} \Sigma^\mathrm{K}(\mathbf{R},\varepsilon) & \Delta^\mathrm{K}(\mathbf{R},\varepsilon)\\ - \tilde{\Delta}^\mathrm{K}(\mathbf{R},\varepsilon) & - \tilde{\Sigma}^\mathrm{K}(\mathbf{R},\varepsilon)\end{pmatrix} . \end{align} \tag{ A4 }$$

These self-energies are defined below.

The elements of $\check{g}$ can be parametrized using coherence amplitudes $\gamma, \tilde{\gamma}$ [58–60] and generalized distribution functions $\textsf{x}, \tilde{\textsf{x}}$ [61, 62]. Using the definitions of $\mathcal{G}^\mathrm R \equiv ( 1 - \gamma^\mathrm{R} \tilde{\gamma}^\mathrm{R})^{-1}$ and $\mathcal{F}^\mathrm{R} \equiv \mathcal{G}^\mathrm{R}\gamma^\mathrm{R}$, we have for the retarded element

$$\begin{align} \hat{g}^\mathrm{R} = -2\pi i \! \begin{pmatrix} \mathcal{G}& \mathcal{F} \\ -\tilde{\mathcal{F}} & -\tilde{\mathcal{G}} \end{pmatrix}^\mathrm{R} + i \pi \hat{\tau}_3, \end{align} \tag{ A5 }$$

and an analogous expression for the advanced function $\hat{g}^\mathrm A$. The Keldysh component reads

$$\begin{align} \hat{g}^\mathrm{K} & = \begin{pmatrix} g& f \\ - \tilde{f} & - \tilde{g} \end{pmatrix}^\mathrm{K} = - 2~\pi i \begin{pmatrix} \mathcal{X} & \mathcal{Y} \\ \tilde{\mathcal{Y}} & \tilde{\mathcal{X}} \end{pmatrix}^\mathrm{K}, \nonumber \\ & \equiv -2~\pi i \begin{pmatrix} \mathcal{G} & \mathcal{F} \\ -\tilde{\mathcal{F}} & -\tilde{\mathcal{G}} \end{pmatrix}^\mathrm{R} \begin{pmatrix} \textsf{x} & 0 \\ 0 & \tilde{\textsf{x}} \end{pmatrix} \begin{pmatrix} \mathcal{G} & \mathcal{F} \\ -\tilde{\mathcal{F}} & -\tilde{\mathcal{G}} \end{pmatrix}^\mathrm{A}. \end{align} \tag{ A6 }$$

An important symmetry in the theory is particle-hole conjugation, which can be expressed as

$$\begin{align} \tilde{A}(\mathbf{p}_\mathrm{F},\mathbf{R},\varepsilon) = A^*(-\mathbf{p}_\mathrm{F},\mathbf{R},-\varepsilon^*). \end{align} \tag{ A7 }$$

A set of coupled transport equations for the parametrizing functions can be derived from equation (A2). The Riccati equation for the coherence amplitude reads

$$\begin{align} \left( i\hbar\mathbf{v}_\mathrm{F}\cdot\nabla + 2~\varepsilon \right) \gamma^{\mathrm{R},\mathrm{A}} \! = \!\bigl( \gamma \tilde{\Delta} \gamma + \Sigma \gamma - \gamma \tilde{\Sigma} - \Delta \bigr)^{\mathrm{R},\mathrm{A}}\!, \end{align} \tag{ A8 }$$

while the equation for the distribution function $\textsf{x}$ is

$$\begin{align} &i \hbar \mathbf{v}_\mathrm{F} \cdot \nabla \textsf{x} - \left[ \gamma \tilde{\Delta} + \Sigma \right]^\mathrm{R} \textsf{x} - \textsf{x} \left[ \Delta \tilde{\gamma} - \Sigma \right]^\mathrm{A}\nonumber\\ & \quad = - \gamma^\mathrm{R} \tilde{\Sigma}^\mathrm{K} \tilde{\gamma}^\mathrm{A} + \Delta^\mathrm{K} \tilde{\gamma}^\mathrm{A} + \gamma^\mathrm{R} \tilde{\Delta}^\mathrm{K} - \Sigma^\mathrm{K}. \end{align} \tag{ A9 }$$

Application of equation (A7) to these two equations gives the corresponding equations for $\tilde{\gamma}^\mathrm{R,A}$ and $\tilde{\textsf{x}}$. 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 $\mathbf{v}_\mathrm{F}$, in our case fully determined by the momentum orientation angle $\varphi_\mathrm{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 [61–63].

Having obtained solutions to equations (A8) and (A9), we can update all self-energies. They consist of the mean field order parameter and scalar impurity self energy as

$$\begin{equation} \check{h}(\mathbf{p}_\mathrm{F},\mathbf{R},\varepsilon) = \check{h}_\mathrm{mf}(\mathbf{p}_\mathrm{F},\mathbf{R}) + \check{h}_\mathrm{s}(\mathbf{R},\varepsilon). \end{equation} \tag{ A10 }$$

The Keldysh matrix structure of the mean field order parameter is simple, $\check{h}_\mathrm{mf} = \hat\Delta\check{1}$, while the Nambu structure is $\hat\Delta = \Re(\Delta)\hat\tau_1-\Im(\Delta)\hat\tau_2$. The order parameter consists of a sum over the d-wave and subdominant s-wave components. Each component satisfies a gap equation

$$\begin{equation} \Delta_\Gamma(\mathbf{R})\! = \!\mathcal{N}_\mathrm{F}\lambda_\Gamma\!\int\limits_{-\varepsilon_\mathrm{c}}^{\varepsilon_{\text{c}}}\!\frac{d\varepsilon}{8\pi i} \left\langle {Tr}\left[i\sigma_2\eta_\Gamma(\mathbf{p}_\mathrm{F}) f\,^\mathrm{K}(\mathbf{p}_\mathrm{F}\, ,\mathbf{R},\varepsilon)\right] \right\rangle_{\scriptscriptstyle\mathrm{FS}}, \end{equation} \tag{ A11 }$$

where $\Gamma = \mathrm{d},\,\mathrm{s}$. The basis functions are $\eta_\mathrm{d}(\mathbf{p}_\mathrm{F}) = \sqrt{2}\cos[2(\varphi_\mathrm{F}-\alpha)]$ and $\eta_\mathrm{s}(\mathbf{p}_\mathrm{F}) = 1$. The pairing interactions $\lambda_\Gamma$ are eliminated in favor of the clean-limit critical temperatures for each channel by the replacement

$$\begin{equation} \frac{1}{\mathcal{N}_\mathrm{F}\lambda_\Gamma}\rightarrow \ln\left(\frac{T}{T_{\mathrm{c}\Gamma}}\right) +\int^{\varepsilon_\mathrm{c}}_{\varepsilon_\mathrm{c}} \frac{d\varepsilon}{2} \frac{1}{{\varepsilon}} \tanh\left(\frac{\varepsilon}{2k_\mathrm{B} T}\right). \end{equation} \tag{ A12 }$$

In the main text, since the d-wave is the dominant component, we let $T_\mathrm{c0} = T_\mathrm{cd}$ and use $k_\mathrm{B}T_\mathrm{c0}$ as the natural energy scale. The critical temperature of the subdominant s-wave is for short denoted $T_\mathrm{s}\lt T_\mathrm{c0}$.

The impurity self-energies are calculated in the non-crossing approximation for the t-matrix equation [45], discussed in detail in [41]. The matrix structure is

$$\begin{align} \check{h}_\mathrm s = n_\mathrm{i} \check{t} \equiv n_\mathrm{i} \begin{pmatrix} \hat{t}^\mathrm{\,R} & \hat{t}^\mathrm{\,K} \\ 0 & \hat{t}^\mathrm{\,A} \end{pmatrix}. \end{align} \tag{ A13 }$$

For scattering that is isotropic in momentum space with an s-wave scattering potential u₀ the elements of $\check{t}$ satisfy the equations

$$\begin{align} \hat{t}^\mathrm{\,R,A} & = \frac{u_0~\hat{1} + u_0^2~\mathcal{N}_\mathrm{F} \left\langle \hat{g}^\mathrm{R,A} \right\rangle_{\scriptscriptstyle\mathrm{FS}}}{\hat{1} - \left[ u_0~\mathcal{N}_\mathrm{F} \left\langle \hat{g}^\mathrm{R,A} \right\rangle_{\scriptscriptstyle\mathrm{FS}}\right]^2 }, \end{align} \tag{ A14 }$$

$$\begin{align} \hat{t}^\mathrm{\,K} & = \mathcal{N}_\mathrm{F} \hat{t}^\mathrm{\,R} \left\langle \hat{g}^\mathrm{K} \right\rangle_{\scriptscriptstyle\mathrm{FS}} \hat{t}^\mathrm{\,A}.\end{align} \tag{ A15 }$$

This procedure of solving equations (A8) and (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_\mathrm{I}|/|j_\mathrm{I}| \lt 5 \times 10^{-3}$, so the current everywhere deviates less than half a percent from $j_\mathrm{I}$, the current at the interface to the normal-metal reservoirs. The charge current is determined by the Keldysh component $\hat{g}^\mathrm K$,

$$\begin{align} \mathbf{j}(\mathbf{R}) & = e\mathcal{N}_\mathrm{F}\!\int\limits_{-\infty}^{\infty}\! \frac{\mathrm{d}\varepsilon}{8\pi i} \left\langle \mathrm{Tr}\left[ \mathbf{v}_\mathrm{F} \hat{\tau}_3~\hat{g}^\mathrm{K}(\mathbf{p}_\mathrm{F},\mathbf{R},\varepsilon) \right] \right\rangle_{\scriptscriptstyle\mathrm{FS}}, \end{align} \tag{ A16 }$$

as is the local electrochemical potential

$$\begin{equation} \phi(\mathbf{R}) = \frac{1}{2e} \int\limits_{-\infty}^{\infty} \frac{\mathrm{d}\varepsilon}{8~\pi i} \left\langle \mathrm{Tr}~\hat{g}^\mathrm{K}(\mathbf{p}_\mathrm{F},\mathbf{R},\varepsilon) \right\rangle_{\scriptscriptstyle\mathrm{FS}}. \end{equation} \tag{ A17 }$$

The retarded component $\hat{g}^\mathrm R$ determines the normalized, momentum-resolved local density of states per spin

$$\begin{equation} \mathcal{N}(\mathbf{p}_\mathrm{F}, \mathbf{R}, \varepsilon) = -\frac{1}{4\pi}\mathrm{Im}~\mathrm{Tr}\left[\hat{\tau}_3~\hat{g}^\mathrm R(\mathbf{p}_\mathrm{F}, \mathbf{R}, \varepsilon)\right], \end{equation} \tag{ A18 }$$

and thus the full density of states

$$\begin{equation} N(\mathbf{R},\varepsilon) = 2\mathcal{N}_\mathrm{F}\left\langle \mathcal{N}(\mathbf{p}_\mathrm{F}, \mathbf{R}, \varepsilon) \right\rangle_{\scriptscriptstyle\mathrm{FS}}. \end{equation} \tag{ A19 }$$

The Keldysh component $\hat{g}^\mathrm{K}$ can, in addition to equation (A6), also be parametrized as

$$\begin{align} \hat{g}^\mathrm{K} = \hat{g}^\mathrm{R} \hat{f} - \hat{f}\, \hat{g}^\mathrm{A}, \end{align} \tag{ A20 }$$

where the distribution matrix $\hat{f}$ reads

$$\begin{align} \hat{f} = \begin{pmatrix} h & 0 \\ 0 & -\tilde{h} \end{pmatrix} = f_1~\hat{1} + f_3~\hat{\tau}_3. \end{align} \tag{ A21 }$$

The non-equilibrium electron distribution function is then given by

$$\begin{align} f_\mathrm{e} = \tfrac{1}{2}(1 - h) = \tfrac{1}{2}\left(1 - f_1 - f_3\right), \end{align} \tag{ A22 }$$

which becomes a Fermi–Dirac distribution function in equilibrium. We refer to f₁ (f₃) as the energy-like (charge-like) mode, connecting to the nomenclature established for diffusive superconductors [64]. After Fermi-surface averaging, they have different symmetries as function of energy, namely

$$\begin{align} \left\langle\, f_1~\right\rangle_{\scriptscriptstyle\mathrm{FS}}(\varepsilon) = - \left\langle\, f_1~\right\rangle_{\scriptscriptstyle\mathrm{FS}}(-\varepsilon), \end{align} \tag{ A23 }$$

$$\begin{align} \left\langle\, f_3~\right\rangle_{\scriptscriptstyle\mathrm{FS}}(\varepsilon) = \left\langle\, f_3~\right\rangle_{\scriptscriptstyle\mathrm{FS}}(-\varepsilon), \end{align} \tag{ A24 }$$

which corresponds to the symmetries obtained in the diffusive case. The function h can be related to $\textsf{x}$ by

$$\begin{align} \textsf{x} = h + \gamma^\mathrm{R} \tilde{h} \tilde{\gamma}^\mathrm{A}. \end{align} \tag{ A25 }$$

The equation of motion for $\textsf{x}$ and $\tilde{\textsf{x}}$ is easier to solve, while the final results can be easier to interpret in terms of h, or $f_{\mathrm{e}/1/3}$. As an example, we can write the gap equations in (A11) as

$$\begin{equation} \frac{\Delta_\Gamma(\mathbf{R})}{\mathcal{N}_\mathrm{F}\lambda_\Gamma/2}\! = \! \!\int\limits_{-\varepsilon_\mathrm{c}}^{\varepsilon_\mathrm{c}}\!\mathrm{d}\varepsilon \left\langle \eta_\Gamma(\mathbf{p}_\mathrm{F}) \left[\, f_1( \mathcal{F}^\mathrm R + \mathcal{F}^\mathrm A) - f_3( \mathcal{F}^\mathrm R - \mathcal{F}^\mathrm A) \right] \right\rangle_{\scriptscriptstyle\mathrm{FS}}, \end{equation} \tag{ A26 }$$

and the charge current, equation (A16), as

$$\begin{align} \mathbf{j}(\mathbf{R})\! = \!-e\mathcal{N}_\mathrm{F} \int\limits_{-\infty}^{\infty} \mathrm{d}\varepsilon \left\langle \mathbf{v}_\mathrm{F} f_1(\mathbf{p}_\mathrm{F},\mathbf{R},\varepsilon) \mathcal{N}(\mathbf{p}_\mathrm{F},\mathbf{R},\varepsilon) \right\rangle_{\scriptscriptstyle\mathrm{FS}}. \end{align} \tag{ A27 }$$

The natural unit for charge current is then $j_0 = e v_\mathrm{F}\mathcal{N}_\mathrm{F}k_\mathrm{B}T_\mathrm{c0}$. The electrochemical potential, equation (A17), can similarly be written as

$$\begin{equation} \phi(\mathbf{R}) = -\frac{1}{2e} \int\limits_{-\infty}^{\infty} \mathrm{d}\varepsilon \left\langle\, f_3(\mathbf{p}_\mathrm{F},\mathbf{R},\varepsilon) \mathcal{N}(\mathbf{p}_\mathrm{F},\mathbf{R},\varepsilon) \right\rangle_{\scriptscriptstyle\mathrm{FS}}.\end{equation} \tag{ A28 }$$

The distribution function h can be split into a local-equilibrium function $h^\mathrm{le}$, and an anomalous function $h^\mathrm{a}$,

$$\begin{align} h = h^\mathrm{le} + ( h - h^\mathrm{le}) \equiv h^\mathrm{le} + h^\mathrm a, \end{align} \tag{ A29 }$$

where $h^\mathrm{le}$ is given by

$$\begin{equation} h^\mathrm{le}(\mathbf{R},\varepsilon) = \tanh\frac{\varepsilon-e\phi(\mathbf{R})}{2T}. \end{equation} \tag{ A30 }$$

The splitting introduced in equation (A29) directly translate to an analogous splitting of the function $\textsf{x}$ by using equation (A25), and for the two modes f₁ and f₃ by using equation (A21). As a consequence observables such as the charge current can be similarly split,

$$\begin{equation} \mathbf{j}(\mathbf{R}) = \mathbf{j}^\mathrm{le}(\mathbf{R}) + \mathbf{j}^\mathrm{a}(\mathbf{R}). \end{equation} \tag{ A31 }$$

Since $h^\mathrm{le}$, and hence $f\,_\mathrm{1,3}^\mathrm{le}$, are independent of $\mathbf{p}_\mathrm{F}$, equation (A27) shows that the local-equilibrium current $\mathbf{j}^\mathrm{le}$ can only be non-zero if there is an asymmetry in the angle-resolved spectrum $\mathcal{N}(\mathbf{p}_\mathrm{F},\mathbf{R},\varepsilon)$. In our case, this asymmetry is created by the self-consistently determined superflow $\mathbf{p}_\mathrm{s}$ in the superconductor. The anomalous current $\mathbf{j}^\mathrm{a}$, due to the distribution $h^\mathrm{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 $\phi_+$ ($\phi_-$), defined as

$$\begin{equation} \phi_{\pm} :\, = \phi -\frac{1}{2e}\int\limits_{-\infty}^{\infty}\! \frac{d\varepsilon}{2}~ \left( \langle\mathcal{X}^\mathrm{a} \rangle_{\pm} + \langle \tilde{\mathcal{X} }^\mathrm{a} \rangle_\mp \right), \end{equation} \tag{ A32 }$$

where $\mathcal{X}^\mathrm a$ is obtained by using $\textsf{x}^\mathrm{a}$ in equation (A6). In equation (A32), the $\langle \dots \rangle_\pm$ denotes a partial Fermi-surface average

$$\begin{align} \langle A \rangle_\pm &\equiv \int\limits_{0}^{2\pi}\! \frac{d\varphi_\mathrm{F}}{\pi} \!~ A(\varphi_\mathrm{F}) \Theta(\pm \cos \varphi_\mathrm{F}), \end{align} \tag{ A33 }$$

where the Heaviside step function $\Theta(\pm\cos\varphi_\mathrm{F})$ gives unity if $v\,_\mathrm{F}^x$, as defined in equation (2), is positive (+) or negative (−), and zero otherwise. In the normal state, the left-mover and right-mover quasipotentials satisfy $\phi = (\phi_+ + \phi_-)/2$, which connects to the concepts used in normal-state ballistic systems [65].