Disorder-robust $p$-wave pairing with odd frequency dependence in normal metal-conventional superconductor junctions

We investigate the induced superconducting pair correlations in junctions between a conventional spin-singlet $s$-wave superconductor and a disordered normal metal. Decomposing the pair amplitude based on its symmetries in the time domain, we demonstrate that the odd-time, or equivalently odd-frequency, spin-singlet $p$-wave correlations are both significant in size and entirely robust against random non-magnetic disorder. We find that these odd-frequency correlations can even be generated by disorder. Our results show that anisotropic odd-frequency pairing represent an important fraction of the proximity-induced correlations in disordered superconducting hybrid structures.


I. INTRODUCTION
Superconductivity arises from the formation and condensation of a macroscopic number of electron pairs, or Cooper pairs. The superconducting state is characterized by the Cooper pair amplitude which, due to the fermionic nature of electrons, is antisymmetric under the exchange of all degrees of freedom describing the paired electron states. These symmetries include the positions, spin, and relative time separation of the two electrons.
In BCS theory of superconductivity the interparticle interaction is instantaneous and, hence, only static, equal-time, pair amplitudes contribute to the order parameter. Such equal-time pair symmetries are then constrained to be either even in space and odd in spin, as in the conventional spin-singlet s-wave state, or odd in space and even in spin (e.g. spin-triplet p-wave). However, even for BCS superconductors, a growing body of evidence suggests that dynamical pair correlations that are odd in the relative time [1][2][3][4][5] can play a role in the physics of superconducting heterostructures [6][7][8]. Due to their odd time-dependence, these dynamical pair correlations must be either even in space and even in spin or odd in space and odd in spin. Such exotic odd-time pairs are equivalently also referred to as odd-frequency (odd-ω) pairs [1].
Odd-ω pairing can also emerge in normal metalconventional superconductor (NS) junctions. Here the * Corresponding author: jorge.cayao@physics.uu.se interface allows for both the original even-ω spin-singlet s-wave (isotropic) and an induced odd-ω spin-singlet pwave (anisotropic) symmetry in the N region [8,[25][26][27][28]. However, unlike FS junctions, in disordered NS systems odd-ω pairing has so far been largely ignored. One oftencited reason for this neglect is the argument that p-wave amplitudes are killed in the presence of disorder [29][30][31]. For example, in the quasiclassical Usadel formulation for dirty systems, only isotropic s-wave solutions are explicitly considered, while any anisotropic part is assumed to be negligible [32]. Additionally, Anderson's theorem [33] has been invoked to defend the assertion of fragility of p-wave pairs in disordered NS junctions, even though it does not technically concern proximity-induced effects.
In this work we explicitly show that in NS junctions odd-ω spin-singlet anisotropic p-wave pair correlations are not only robust against disorder but as disorder strength increases, the odd-ω amplitudes constitute a growing fraction of the proximity-induced pair correlations in the N region. Remarkably, we even find odd-ω pairing generated by disorder. We do this by considering a NS junction between a disordered N region and a conventional spin-singlet s-wave superconductor, see Fig. 1, explicitly calculating the fully quantum mechanical time evolution of the proximity-induced superconducting correlations. By showing that odd-ω p-wave pair correlations are ubiquitous and robust, our results provide a completely different view of superconducting pair correlations in disordered heterostructures, with potentially far reaching consequences for the physics, such as a paramagnetic Meissner effect.
The remainder of this paper is organized as follows. In Sec. II we present the model and outline the method to obtain the superconducting pair amplitudes. In Sec. III we present the results for the clean and disordered regimes in 1D and also show some consequences in 2D. We conclude in Sec. IV. For completeness, in the Appendices we analytically obtain the induced pair amplitudes in 1D for the clean regime analytically using perturbation theory. Furthermore, here we also present further numerical calculations for 2D. where filled blue and red circles represent sites in N (x < 0) and S (x > 0) regions. The relevant spatial orbital symmetries, s-wave, extended s-wave, and px,y-wave, are depicted in blue (red) for positive (negative) amplitudes.

II. MODEL AND METHODS
To demonstrate the robustness of odd-time/odd-ω pair correlations we study the proximity-induced pair amplitudes in both one-dimensional (1D) and 2D NS junctions between a disordered N region and a conventional spinsinglet s-wave S region, with the interface at x = 0, see where c jσ is the destruction operator for a particle with spin σ at site j = (j x , j y ) in 2D or j = j x in 1D, . . . denotes nearest neighbor sites, t the nearest neighbor hopping amplitude, and µ the chemical potential that controls the band filling. The conventional spinsinglet s-wave superconductor is captured by H ∆ = i θ(x)∆c † i↑ c † i↓ + h.c., with ∆ being the order parameter and θ the Heaviside step function. We further introduce random scalar, i.e. non-magnetic, potential disorder on the normal side of the junction using the Anderson disorder model: where the uncorrelated potentials V i are drawn from the uniform box distribution on the interval [−δ, δ] and added to every site of our discrete lattice model [34]. The disorder strength is controlled by the width of the distribution δ/t which allows us to tune the disorder from the clean limit or weakly disordered regimes (all sites have none or weak scatterers) to extremely strong disorder (some sites host strong scatterers, while other stays weak), where the disorder dominates over the kinetic energy. This type of disorder can arise from non-magnetic charge impurities or fluctuations in the chemical potential, both inevitable in real experimental samples.

A. Pair amplitudes
We are interested in the superconducting correlations in NS junctions. These are characterized by the pair amplitudes which can be obtained from the time-ordered anomalous Green's function F ijαβ (t 2 , t 1 ) = T c iα (t 2 )c jβ (t 1 ) , where T is the time-ordering operator. The pair amplitude F ijαβ (t 2 , t 1 ) accounts for the propagation of quasiparticles from the (iα)-state at time t 1 to the (jβ)-state at time t 2 . Since F ijαβ (t 2 , t 1 ) only depends on the relative time difference τ = t 2 − t 1 , we focus on the forward time propagation, where the even-(e) and odd-time (o), or equivalently even-and odd-ω, components are extracted as F (0) )/2. Moreover, with no spin-active fields in NS junctions, the pair amplitudes exhibit the same spin symmetry as the parent superconductor. Hence, the pair amplitudes F e(o) ijαβ (τ ) can be decomposed into the pair wave amplitudes [35], where we focus our study on the lowest rotation symmetric wave representation projections with |r| ≤ 1: w S (r) = δ |r|=0 , w Sext (r) = (1/z)δ |r|=1 , w Pi (r) = [(î·r)/2]δ |r|=1 are the s-, extended s-wave, and p i -wave amplitudes. Here z is the coordination number, σ y the y-Pauli matrix, and the center of mass R = (i+j)/2, and relative r = (i−j)/2 spatial coordinates. While we expect that higher order representations also have a finite magnitude [25], to simplify the presentation we focus on |r| ≤ 1 which allows for a targeted and controlled investigation of the interplay between disorder and anisotropic odd-ω pair correlations. Evaluation of the pair wave amplitudes in each symmetry channel in Eq. (2) only requires computation of timedependent expectation values of the form c iα (τ )c jβ (0) . We perform these calculations numerically using the EPOCH (Equilibrium Propagator by Orthogonal polynomial CHain) method for computing the equilibrium propagators, i.e. Green's functions, directly in the timedomain that we have developed in Ref. [36]. In the EPOCH method the equilibrium propagator is expanded in Legendre polynomial mode transients that are recursively computed, allowing us to study both the even-and odd-ω pair wave amplitudes directly in real space as well as the time domain. As such, this analysis notably goes beyond previous studies of odd-ω pairing [8,25,26,28], since the treatment is explicitly in the more natural time domain and also fully quantum mechanical.

III. PAIR AMPLITUDES IN CLEAN AND DISORDERED NS JUNCTIONS
To ultimately understand the pair wave amplitude behavior in the presence of disorder in NS junctions, we first study its time evolution without disorder, δ/t = 0. Then, we analyze the disordered regime by varying the disorder strength δ/t.

A. Clean regime
In Fig. 2 we show space and time evolution of the magnitude of the even-ω s-wave and odd-ω p-wave pair wave amplitudes, |d e,o (τ, R)| in a clean 1D NS junction. We explicitly find no other finite amplitudes with |r| ≤ 1, apart from even-ω extended s-wave correlations, which show similar behavior as the s-wave pairing. Notably, the p-wave term is not only non-local and odd in space, but also non-local and odd in the relative time coordinate τ , fully consistent with previous results studying the timedependence of odd-ω pair amplitudes [37]. At τ = 0 the even-ω s-wave magnitude exhibits a large and constant value in S, due to the finite order parameter ∆ in S, as see in Fig. 2(a). The odd-ω pair magnitude in (b), on the other hand, is necessarily zero at τ = 0. For finite τ both pair wave amplitudes develop a fan-shaped profile in both S and N regions, albeit with smaller odd-ω magnitudes. The fan profile is due to wave packet propagation in time and space, with a decay as we move away from the NS interface, since the superconducting proximity effect diminishes with distance [38]. There is also an oscillatory pattern, seen as resonances in Fig. 2(c,d), due to interference between incident and reflected quasiparticle states at the NS interface, as noted before [28,37,39,40]. In particular, the resonance peaks in the N region are a result of the following time evolution process: an electron at x 1 in N has to first propagate towards the NS interface, where it is reflected back as a hole through Andreev reflection, and then propagate back to x 2 in N in order to contribute to the pair wave amplitude [41]. A simple analytic calculation using a 1D continuum model yields the expected time of this process to τ 0 = d/v F , where d = |x 1 |+|x 2 | is the round trip distance and v F the Fermi velocity, see Supplemental Material (SM) for a complete derivation [42]. This value of τ 0 is fully consistent with the results in Fig. 2(c,d). The overall result is that large pair magnitudes disperse linearly away from the interface in both the N and S regions, in agreement with the role of the NS interface as the generator of even-and odd-ω correlations [8,25,26,28].

B. Disorder regime
Next we turn to our main topic on how the induced pair correlations are affected by disorder in the N region. Our goal is to ascertain if the odd-ω p-wave correlations are very fragile to disorder as previously assumed [29][30][31].
To this end, we investigate how the odd-ω p-wave correlations are affected by random scalar potential disorder. As a control, we also investigate how both the on-site and extended s-wave pair correlations are affected, since both of these isotropic wave states are assumed to be disorder robust. The extended s-wave also provides additional control in that it has the same spatial extent as the odd-ω p-wave, both being nearest-neighbor correlations. Our main results are presented in Fig. 3, where we plot the induced pair correlations in a 1D NS system for a wide range of disorder strength δ/t, from the clean limit δ/t = 0 to the strongly disordered limit δ/t = 2. For any finite δ/t, the disorder produces a finite elastic mean free path l e ≈ 4.3a(t/δ) 2 , see SM [42]. Our range of disorder therefore spans both the regimes where the mean free path is significantly shorter and longer than the superconducting coherence length ξ = v F /∆ ∼ 20a.
We first present the time evolution of the disorderaveraged pair magnitudes with the on-site even-ω s-wave in Fig. 3(a) and the odd-ω p-wave in Fig. 3(b). Both magnitudes are taken at the fixed spatial position indicated by the solid white line in Fig. 2. At the vanishing disorder strength, we recover the same situation as discussed in Fig. 2, where both the even and odd-ω pair magnitudes are finite in the N region and exhibit strong oscillations in time. For both magnitudes, the oscillatory patterns in time gradually diminish with increased disorder, but the changes in the magnitudes are clearly different between the even and odd-ω correlations. Extracting the pair magnitudes at their first resonance peak, indicated by the circle markers in Fig. 3(a,b) and plotted in Fig. 3(c) as a function of disorder strength, we show that, while the on-site and extended s-wave correlations are both monotonically decreasing with the disorder strength, the odd-ω p-wave correlations are non-monotonic with an initial enhancement with increasing disorder. Thus, as we present in Fig. 3(d), the ratio of odd-ω p-wave to even-ω s-wave correlations actually increases with the disorder strength. In particular, the odd to even ratio undergoes a rapid growth around the crossover l e ∼ ξ, appearing

25.
10 3 Pair Magnitude (c)  near δ/t ≈ 0.5, and is monotonically increasing for all l e /ξ, all the way to the dirty limit l e ξ. This is in complete contrast to the expected disorder fragility of all anisotropic non-s-wave states [29][30][31] and clearly shows that the p-wave pairs cannot be ignored even in the dirty limit. The relative increase of the odd-ω p-wave correlations is sustained throughout the N region, even beyond the superconducting coherence length ξ ∼ 20a, as is clearly evident from the ratios calculated at multiple distances away from the NS junction in Fig. 3(d). This rules out a mere interface effect.
To examine the stability of these disorder effects for different parameter values, we present in Fig. 3 (e,f) the disordered averaged pair magnitudes as a function of the chemical potential µ, again at the first resonance peak (circles in Fig. 3(a,b)). In the clean regime the even-ω swave pair magnitude is almost constant for all values of µ, while the odd-ω p-wave amplitude identically vanishes at half-filling, µ = 0. The latter is due to particle-hole symmetry in the normal state at µ = 0, resulting in a cancellation of odd-time correlations [43]. While this cancellation of the odd-ω amplitude at µ = 0 in the clean limit occurs in a fine-tuned situation, what is highly relevant for the general understanding of p-wave pair correlations, is that they still become finite at µ = 0 when introducing disorder. Thus disorder can even generate p-wave pairing. Furthermore, as seen in Fig. 3(f), where we plot the p-wave to s-wave ratio, the p-wave pair correlations are large, at around 40%, and essentially constant for disorder ranging from moderately weak to very strong. These results highlight the crucial role of non-magnetic disor-der, with disorder not suppressing odd-ω p-wave pairing, but even generating it. Additionally, we have also verified that qualitatively both the clean and disordered behaviors are independent of the size of the superconducting order parameter for all ∆ < t The exceptional disorder robustness of the odd-ω pwave correlations is a consequence of the specific mechanism by which they appear. In the clean limit, p-wave pair correlations in NS junctions are generated due to the spatial translation symmetry breaking at the interface. Fermion statistics dictates that these p-wave correlations must have an odd-ω symmetry. Introducing nonmagnetic disorder in the N region cause two diametrically opposing effects. First, disorder reduces the propagation of Cooper pairs, especially anisotropic pairs. Secondly, disorder further disrupts the translation symmetry and thus generates more odd-ω p-wave correlations. The exact balance between these two effects can only be handled numerically and our results show that the generation even outweighs the suppression in propagation of p-wave pairs. Thus odd-ω spin-singlet p-wave pairing is remarkably robust to disorder in dirty NS junctions. This result is in stark contrast to the dirty limit, or Usadel, formulation of quasiclassical theory, where only isotropic s-wave pairing states are explicitly considered, while any anisotropic pairing terms are assumed to vanish in the dirty limit l e /ξ 1 [29][30][31]. It hence follows that any proper quasiclassical treatment has to include anisotropic pair channels using e.g. the Eilenberger equations [44,45]. We also stress that our results are not in contradiction to Anderson's theorem [32], as the pair amplitudes in N are proximity-induced, while Anderson's theorem only formally applies to bulk superconductors.

C. Beyond 1D
In the above we focused on the results for 1D NS junctions, but our main findings also hold for NS junctions in higher dimensions. In the SM [42] we show plots that are equivalent to Fig. 2-3 but for a 2D NS junction and that are analogously obtained by calculating the time evolution of Eq. (1). We have likewise verified similar results for 3D junctions. The fact that we find similar results across multiple spatial dimensions means that our main results are independent of any localization effects that can otherwise be pronounced in low dimensional systems [46], since the localization lengths are of different orders of magnitude in different dimensions [47].
While the results in higher dimensions are strikingly similar to those in 1D, and also more stable due to increased self-averaging, there are some importance differences. In addition to the on-site and extended even-ω s-waves and the odd-ω p x -wave correlations previously seen in 1D, there exist in 2D also a odd-ω p y -wave correlations that is oriented along the junction interface, see Fig.1. In Fig. 4 we plot the induced disorder-averaged pair magnitudes versus disorder strengths in a 2D SN junction, now also with the p y -wave correlations. The plot can be directly compared to the analogous plot in Fig.3 for 1D, given the very similar elastic mean free path in the 2D, l e ≈ 4.7a(t/δ) 2 , see SM [42]. In the clean limit δ/t = 0, the p y -wave magnitude is completely absent, because the interface only breaks the translation invariance along the x-direction and therefore no correlations with a p y -wave symmetry can be generated. Still, at finite disorder, we see that finite odd-ω p y -wave pair amplitudes are generated alongside the initial enhancement of the already finite odd-ω p x -wave amplitude. In contrast, the presumed disorder robust isotropic s-wave states instead show a pronounced monotonic suppression with increasing disorder. This demonstrates again that the odd-ω p-wave pair amplitudes are highly disorder robust and always constitute a significant portion of the induced pairing in disordered NS junctions. The emergence of odd-ω correlations along the y-axis with disorder can generally be understood as arising from the additional breaking of translational invariance along the y-axis by the disorder and its resulting scattering of pair wave amplitudes into the p y -wave channel, although we point out that this alone is not enough. In fact, random charge potential disorder cannot by itself generate odd-ω pair correlations, as shown in Ref. 37. What is needed in addition is the NS junction and its associated spatial gradient in the superconducting order parameter, which therefore plays an important role in the disorder generation of the odd-ω p y -wave pair amplitudes.

IV. CONCLUSIONS
In conclusion, we have demonstrated that odd-ω spinsinglet p-wave pair correlations are robust against random scalar, non-magnetic, disorder in conventional and dirty NS junctions. In fact, the p-wave pair correlations display a non-monotonic dependence on the disorder strength, with an initial enhancement and also constitute a growing fraction of the induced pair correlations in the N region, notably even in the highly disordered regime where the elastic mean free path is much smaller than the superconducting coherence length. The p-wave pair correlations therefore display an even greater disorder robustness than the isotropic s-wave components that are instead monotonically decreasing with an increased disorder. As a consequence, odd-ω p-wave pair correlations always constitute a significant portion of the proximity-induced superconductivity in dirty NS junctions. With properties, such as the Josephson and Meissner effects, being directly dependent on the superconducting pair amplitudes, we expect our results to have interesting consequences for various NS systems and superconducting hybrid structures.

V. ACKNOWLEDGEMENTS
We thank A. Balatsky, D. Chakraborty, and T. Löfwander for useful discussions. We acknowledge financial support from the Swedish Research Council Supplemental Material for "Disorder-robust p-wave pairing with odd frequency dependence in normal metal-conventional superconductor junctions" In this supplementary material we provide details to further support the results and conclusions of the main text.

I. ANALYTIC TREATMENT OF PAIR AMPLITUDES IN CLEAN NORMAL METAL-SUPERCONDUCTOR JUNCTIONS
In the main text we presented numerical results for the proximity-induced superconducting pair amplitudes in the normal region of normal metal-superconductor (NS) junctions, modeled using a tight-binding Hamiltonian. An important advantage of our numerical approach is that it allows us to easily consider both the clean and disordered regimes, with no approximations. However, since numerical analyses can only be performed for a finite number of parameters, insight can often be gained from a complementary approximate analytic calculation. In this Section we study the induced pair amplitudes in a one-dimensional (1D) NS junction in the clean regime analytically using perturbation theory.

A. Bare Green's functions
To facilitate our analytic treatment we assume that both sides of the NS junction are infinite and well described by continuum models with dispersions given by ξ k = k 2 /2m − µ and E k = ξ 2 k + ∆ 2 , with = 1, for the N and S regions, respectively. Here, m denotes the effective mass of the electrons in both the N and S regions, k is the wave vector, µ is the chemical potential, and ∆ is the superconducting order parameter in the S region. Given these dispersions, and first assuming no tunneling between the two regions, it is straightforward to show that the Nambu space Green's functions in the N and S regions are given by [1,2]: whereτ i (τ 0 ) is the i th Pauli matrix (identity matrix) in the 2×2 particle-hole space and z is a complex frequency allowing us to keep track of Matsubara (z = iω n ), retarded (z = ω + i0 + ), and advanced (z = ω − i0 − ) Green's functions, simultaneously. Note that the pair amplitude, F , is simply given by the off-diagonal matrix element of the Nambu space Green's function, and, therefore, it is trivially zero in the N region and given by the term proportional toτ 1 in the S region in the limit of no tunneling between the two regions.

B. NS junction
To study proximity-induced pairing in the N region, we assume the NS interface allows local quasiparticle tunneling between the two systems at x = 0. In this case, the Dyson equation for the full Nambu-space Green's function in the N region may be written in real space as: This equation can be iterated to obtain a power series expression forĜ N (x 1 , x 2 ; z) to arbitrary order in the tunneling, t 0 . Since we are only interested in the qualitative behavior of these expressions we truncate this expansion to leading order in t 0 :Ĝ The bare Green's functions on the right-hand side of Eq. (3) may be obtained by Fourier transforming the k-space expressions in Eqs. (1) to real space giving: where These expressions can now be inserted into Eq. (3) to find analytic expressions for the Nambu space Green's function in the N region. Note that the bare Green's functions in real space given by Eqs. (4) depend solely on one coordinate as they are the Fourier transform of Eqs. (1) which describe continuum (infinite) systems where translational invariance is not broken.

C. Proximity-induced pair amplitudes
As we are primarily interested in the induced pair amplitudes in the N region, we focus on the anomalous component of the Green's function in Eq. (3), which is given by: where v F = k F /m. This expression for the proximity-induced pairing in the N region is exact to order t 2 0 . However, we would like to simplify this expression to gain a better understanding of the physical properties of these pair amplitudes. To proceed, we note that the pair amplitudes are strongly peaked for frequencies z ∼ ∆. Moreover, we recall that, in typical superconductors, ∆ << µ. Therefore, we assume that z/µ, ∆/µ << 1 and expand all terms to leading order in these quantities. After a fair amount of algebra we obtain the approximate expression: which is accurate to leading order in t 0 , z/µ, and ∆/µ. Now that we have a relatively compact expression for the proximity-induced pairing in the N region, we wish to make contact with the numerical results in the main text by Fourier-transforming to the time domain. First, we recall that, in general, a time-ordered propagator is related to the retarded and advanced propagators by: Using this relation, together with the expression in Eq. (7) the time-ordered pair amplitude in the N region may be written as: where we define τ 0 = |x1|+|x2| v F , which corresponds to the total time for a quasiparticle at position x 1 in N to travel towards the NS interface and then back to position x 2 in N, as discussed in the main text. By inspection of the terms within the curly braces in Eq. (9), we see that the proximity-induced pair amplitude possesses both even and odd components in the time variable, t. These even-and odd-time components are given by: Focusing on the odd-time/odd-frequency component, and changing integration variables to u = ω/∆, we find: where in the second line we have written the exact integral on the interval [0, ∞) as two approximate integrals: (i) the interval where we assume ω/∆ ≤ 1, and (ii) the interval where we assume ω/∆ >> 1. Simplifying these expressions we can write the odd-time/odd-frequency pair amplitude in the N region as: where E 1 (z) = ∞ 1 e −tz /tdt is the exponential integral. One of the most important features to observe from Eq. (12) is that the last term gives rise to a sharp peak in time, at t = τ 0 = |x1|+|x2| v F , in full agreement with the numerical results presented in the main text. In fact, this peak has an intuitive origin in the physical processes giving rise to all pair amplitudes in the N region: τ 0 is the amount of time it takes a quasiparticle to complete a trip from x 1 to the interface, reflect as a hole at the interface, and travel back to the point x 2 . Additionally, we note that the expressions given by Eqs. (10) for the even-and odd-time pair amplitudes reveal an oscillatory behavior in both time and space, consistent with the numerical results presented in the main text for the clean regime in 1D. Furthermore, we see here that the periodicity of the spatial oscillations is determined by the Fermi wave vector k F , while the periodicity of the temporal oscillations is determined by ∆. Lastly, we point out that a similar analysis can be carried out for higher dimensions, where the main features discussed in this part are expected to hold.

S4
E. Comparison to numerical results in 1D In Fig. 2(a,b) in the main text we showed that in the clean regime both the even-and odd-time pair amplitudes develop resonance peaks for a finite time separation, which then propagate away from the NS interface, resulting in a fan-like profile. In this Section we show that the time separation at which the first resonance peak develops, defining the outer structure of the fan, is the travel time for making the round trip distance to the junction at the Fermi velocity of the normal state. We first show in Fig. S1 (a) that with an increasing time separation, the resonance peaks of Fig. 2 in the main text move linearly away from the junction, implying that the resonance peaks have a definite speed, here called v R . According the the analytical calculations presented in the previous section, this speed is the Fermi velocity v F . To test numerically if v F is indeed the speed in the numerical calculations, it does not suffice to just compute the time separation and position of the resonance peaks as in Fig. S1(a), but rather we need to also manipulate the Fermi velocity v F , because only then can we observe if the resonance propagation of the resonances is dependent on the Fermi velocity.
The Fermi velocity v F is defined as the slope of the band dispersion, v F = |d (k)/dk| at the Fermi energy. By changing the chemical potential µ the filling of the band is changed and thus the Fermi energy intersects the energy band (k) at new point, with a different slope. Thus, the Fermi velocity v F is manipulated by the chemical potential µ. If the N region had been infinite, it would have had the band dispersion (k) = 2 cos(k) with Fermi velocity v F = |d (k)/dk| = 2| sin(k)|. Even though the N region is finite in the numerical model, we still approach our problem assuming that these ideal expression approximates the band dispersion. With this assumption the Fermi velocity is related to the chemical potential µ through v F = 4 − µ 2 . We can then numerical calculate the speeds of the resonance peaks v R (µ) as the slope of the line in Fig. S1(a) by changing µ. Recognizing that R/a on the y-axis in Fig. S1(a) is half the round-trip distance, we in Fig. S1(b) show how the ideal v F (µ) relates to 2v R (µ). We see that the agreement is almost perfectly v F (µ) = 2v R (µ), exactly as suggested by our analytical results. The minor deviation is easily caused by the fact that changing µ not only changes v F but also the low-energy density of states and thus the overall size of the pair amplitudes. To conclude, it is the Fermi velocity that sets the speed of the resonance peaks at the NS junction, and the time-separation for the peaks is given by the round-trip distance 2R. In closing, we also show in Fig. S1(c) that height of the resonance peaks, i.e. the pair magnitude, slowly decays with distance away from the junction for both odd-ω and even-ω due to the inevitable decay of the proximity effect away from the interface.

II. ADDITIONAL RESULTS FOR 2D NS JUNCTIONS
In the main text we mainly presented the paring amplitudes for 1D NS junctions. Here we show that the main conclusions that we draw from the 1D system also hold in higher dimensions. Starting with the clean regime of a 2D junction, we show, in complete analogy with Fig. 2 in the main text, the propagation of even-ω s-wave and the odd-ω p x -wave pair magnitudes |d e,o (τ, R)| in Fig. S1. As for the 1D case, the pair magnitudes are obtained from Eq. (2) in the main text. As seen, most aspects of the propagation remain unchanged between 1D and 2D junctions. At equal times τ = 0, the even-ω magnitudes in Fig. S2(a) have a largely constant amplitude inside the S region, which spreads into the N region by proximity effect. In contrast, the odd-ω magnitude at τ = 0 is identically zero by symmetry. At finite time separations τ > 0, we find a ridge of large pair magnitudes, a resonance, linearly emanating away from the NS interface, also in agreement with the 1D results. Thus the broad features remain unchanged between 1D and 2D, although there are also some differences. Most notable from Fig. S2 is the absence of the fan-like oscillatory pattern seen in 1D, and, in addition, the pair resonances are much broader, as clearly seen in the position cuts presented in Fig. S2(c,d). Both of these qualitative differences can be directly interpreted as resulting from the change in dimension. To understand this we return to the physical picture suggested by the pair amplitude time evolution: a quasiparticle making a round-trip from the position R in the N region to the NS interface and then back to R. In the 1D junction, there is only one trajectory for this trip. However, in 2D there are a plethora trajectories besides the shortest path, which is the straight line joining the R point and the point on the NS interface closest to it. These 2D trajectories include loops and braids that take full advantage of the extra dimension. However, the total length of these trajectories are always longer than the shortest path and thus their associated travel times will be longer. Still, among the likely paths, most will arrive shortly after the round-trip time for the shortest path. As a consequence, similarly as in the 1D junction, the first resonance peak reflects the travel time of the shortest path. The added variations in travel times in 2D however does two things; first, it considerably broadens the first resonance peak; second, it washes out the sharp fan-like oscillatory pattern seen in 1D. This explains all qualitative differences seen in Fig. S2 compared to the 1D junction. As this picture suggests, we also see in Fig. S2 that the width of the first resonance peak widens with R, reflecting that longer paths have more opportunities to deviate. Less evident from Fig. S2, but still worth mentioning, is that in 2D the height of the pair amplitude at the resonance peaks decay faster than in 1D, partly reflecting that the peak is broadened, but also that the pair amplitude wave disperses geometrically in 2D unlike in 1D.
Finally we also compare the disordered regimes for 1D and 2D NS junctions. In Fig. S3 we show how the disorderaveraged pair magnitudes | d e(o) s(p) | evolve with disorder strength δ, similarly as in Fig. 3 in the main text. The main difference is that in Fig. S3 we show all wave projections, including the even-ω s-wave (a) and extended s-wave (b) states, as well as the odd-ω p x -wave (c) and also the odd-ω p y -wave (d) states. We see in agreement with the 1D results in the main text that the even-ω waves are monotonically decreasing (apart from statistical fluctuation) with increasing disorder, while the odd-ω p x -wave instead shows an initial increase. Additionally, we find that the odd-ω p y -wave, which by symmetry is identically zero in the clean region, is quickly generated by disorder. Thus, considering the relative abundance of odd-ω p-wave to the even-ω s-wave pair magnitudes in Fig. S3(e), we see that the ratio of odd-to even-ω pairing rapidly increases with disorder, with the p y ratio quickly approaching that of the original odd-ω p x -wave wave.

III. ELASTIC MEAN FREE PATH
In this Section we discuss and numerically calculate the elastic mean free path in the disordered normal metal side (N) of the SN junction. Without the superconductor but with a finite disorder strength δ/t, the zero temperature conductivity of the N side of the junction, is given within the Kubo-Greenwood formalism by where l e (E F ) is the elastic mean free path, v F (E F ) is the Fermi velocity, and ρ(E F ) is the density of states, all measured at the Fermi energy E F . Following Refs. [3,4], the transport properties and l e (E F ) are efficiently computed in a fully quantum mechanical manner from the temporal dynamics of a wave-packet by its mean square displacement ∆X 2 (E F , τ ) = δ(E F − H)[X(τ ) −X(0)] 2 /ρ(E F ). This is the case because of the way the wave packet spreads in space is directly related to the diffusion coefficient, which is in turn related to the conductivity through the time-dependent Einstein relation, σ(E F , τ ) = e 2 ρ(E F )D(E F , τ ) .
Therefore, by calculating the diffusion coefficient D, it is possible to obtain the elastic mean free path l e (E F ). In Fig. S4(a) we present the diffusion coefficient D(E F , τ ) computed from the time evolution of a stochastic wavepacket on a L = 10 5 a long one dimensional (1D) lattice according to the methods and algorithms of Refs. [3,4] using 2500 polynomial moments. In a disorder free system, the conduction is ballistic and D(E F , τ ) = v 2 F (E F )τ increases linearly with the elapsed time τ . In contrast, the disorder present in Fig. S4(a) triggers a diffusive regime where, after a short initial time ballistic motion, the diffusion coefficient D(E F , τ ) saturates due to scattering. From the initial short-time ballistic motion, we calculate the Fermi velocity v F (E F ) by a least squares fit, and from the saturated maximum value of the diffusion coefficient we extract the elastic mean free path as To get further insights into the behavior of the mean free path, we plot in Fig. S4(b) the elastic mean free path obtained from Eq. (16) as a function of different disorder strength δ/t. As we can see, l e (E F ) exhibits a monotonically decreasing behavior. The dashed line shows a fit to the inverse square of the disorder strength: l e (E F ) = α(t/δ) 2 . This functional dependence is based on the scattering rate from the Fermi golden rule [5][6][7] and yields a very good fit for α ≈ 4.3a . The same relationship also holds for a 2D lattice but there we find α ≈ 4.7a using E F = 0.5t. A significant difference between the 1D and 2D cases is, however, the duration over which D(E, τ ) is diffusive before entering a localized regime with its telltale decay from D max (E F ). This occurs when the mean square displacement saturates, which defines the localization length λ(E F ) [4,8]. We find an order of magnitude longer localization lengths λ(E F ) in 2D compared with 1D, consistent with earlier results [9], but despite the very similar elastic mean free paths. We use the mean free paths obtained and discussed here to compare with the superconducting coherence length in the main text, Section III B.