Statistical detection of Josephson, Andreev, and single quasiparticle currents in scanning tunneling microscopy

We present a method to identify distinct tunneling modes in tunable superconducting tunnel junction composed of superconducting tip and sample in scanning tunneling microscope. Combining the measurement of the relative decay constant of tunneling current extracted from I-V-z spectroscopy with its statistical analysis over the atomic disorders in the sample surface, we identified the crossover of tunneling modes between single quasiparticle tunneling, multiple Andreev reflection, and Josephson tunneling with respect to the bias voltage. The method enables one to determine the particular tunneling regime independently of the spectral shapes, and to reveal the intrinsic modulation of Andreev reflection and Josephson current that would be crucial for quantum device application of superconductors.

The tunnel junction between superconductors is the heart of modern quantum information devices.
Superconductor qubits, which are one of the most promising scalable qubits right now, heavily rely on the Josephson effect, which is quantum tunneling of Cooper pairs between the two superconductors [1,2].
The superconducting tunnel junction can also probe the nature of superconductivity, because tunneling behavior is highly dependent on the type of superconductors, such as either s-wave or d-wave [3,4], and trivial or topological [5].
Inside the superconducting tunnel junction, the interplay between electrons and Cooper pairs results in various modes of tunneling, i.e., single quasiparticle tunneling, multiple Andreev reflection (MAR) and Josephson tunneling [6]. One important factor that defines the mode of tunneling is the tunneling barrier, like its height [7,8] and capacitance [9,10]. However, in superconducting tunnel devices, tunneling barrier is not controllable, and so the systematic study of the crossover between these tunneling modes for the change in barrier properties is difficult. Scanning tunneling microscopy (STM), in contrast, can precisely control the tunneling distance down to picometer resolution by controlling its tip height, and determine the geometry of the tunneling junction in atomic scale by imaging the surface of superconductor. Thus, STM provides a fitting platform to study the behavior of various tunneling modes with respect to the tunneling barrier [11][12][13][14][15][16].
In this letter, we used STM with a superconducting tip to make a superconducting tunnel junction with vacuum tunneling barrier that is controllable to picometer precision. The spectroscopy based on I-V-z characteristics allows us to extract the relative decaying constant of the tunneling current with respect to the tip height and display the crossover between the single quasiparticle tunneling, MAR, and Josephson tunneling. Moreover, the relative decay constant of different tunneling modes displays distinct spatial distribution relative to the atomic disorders, whose statistical analysis gives complementary information for the transitions between the tunneling modes. Theoretical calculation with tight-binding simulation shows the broad variability of tunneling probability of Andreev reflection as a function of quasiparticle energy, whose resonant tunneling at the superconducting gap edge ultimately translates into intermediate decay rate-constant between single quasiparticle and Josephson tunneling. Our results indicate that the superconducting tunnel junction can be viewed as a parallel resistor circuit of various tunneling regimes, which can be turned on/off or amplified by the junction geometry and tunneling energy.
The experiment was performed in SPECS JT-STM operated at UHV condition (< 10 -10 mbar) and base operating temperature of 1.2 K. The surface of Pb(110) single crystal was cleaned by repeating sputtering-annealing cycles several times. The condition for sputtering was ion energy of 1 kV, Ne gas pressure of about 1 × 10 -5 mbar, emission current of 10 mA, ion current of 12 μA, and duration of 15 min, and the annealing condition was 150 °C for 40 min. The STM tip was coated with Pb by applying bias voltage of 100 V to the tungsten tip and then dip it into the Pb(110) single crystal while limiting the maximum current to be 100 μA with a 10 MΩ resistor in series [17]. Conversely, Pb on the tip can be removed by field emission on Au(111) with bias voltage of 100 V. The I-V curves were acquired by sweeping DC bias without any additional electrical signals. dI/dV curves were obtained by numerical differentiating the I-V curves and Gaussian smoothing.  tip height z is controllable in a picometer precision, which enables precise measurement of the decay of tunneling current I with respect to z. Theoretical models based on Bardeen's formalism [18] predicted that single electron tunneling probability T N decays exponentially with z in the asymptotic case of large z [19].
In case of superconducting tunnel junction, single quasiparticle tunneling dominates when the bias is larger than the sum of the sample superconducting gap Δ s and the tip superconducting gap Δ t [ Fig. 1 [11,20]. The decay constant of single quasiparticle tunneling is same as the one in normal state [18], which we denote as κ N . However, when the bias is lower than Δ s + Δ t , single quasiparticle tunneling is  [9,14,15]. In this regime, inelastic tunneling process adds another factor of E J /E C to the tunneling probability, where E J is Josephson energy and E C is capacitive charging energy [23].
The Ambegaokar-Baratoff formula shows that E J is proportional to the tunneling probability T N [7], while E C stays almost constant with sub-nanometer shift of z because it mostly depends on the macroscopic tip geometry [14,24], so the decay constant becomes 2κ N .    Fig. 2(b)], as expected from the Δ s ≈ 1.35 meV in Pb and Δ t slightly less than Δ s [11,15]. Inside the superconducting gap, the structure of dI/dV first rapidly decreases with approaching the Fermi level, but also reveals an increase of tunneling conductance at the Fermi level particularly for spectra acquired with closest proximity between the tip and the surface. The bias lower than the superconducting gap indicates that the origin of the intra-gap features is MAR or Josephson tunneling, or perhaps both. We note that the feature at the Fermi level is broader compared to the previous reports of Josephson currents in STM junction [11,16], possibly due to the effective junction temperature of 3.2 K that is slightly higher than the base temperature (as estimated by fitting the superconducting gap of Pb(110) [25]).
The decay constant κ(V b ) was extracted by fitting I-z or dI/dV-z curves to the exponential function [25].
We further defined the relative decay constant as κ/κ N , where κ N ≡ κ(V max ). In Fig. 2(c), we plot κ/κ N extracted from I-z curves of Pb tip -Pb(110) junction. It is clearly shown that κ/κ N = 1 for e|V b | > Δ s + Δ t as expected for single quasiparticle tunneling, but as e|V b | becomes smaller than Δ s + Δ t , it increases from 1 to almost exactly 2 around V b = 0 mV, as expected for Josephson tunneling [14,23]. The transition between κ/κ N = 1 and κ/κ N = 2 proceeds through two nearly flat plateaus at around |V b | = 1.0 ~ 1.5 mV with κ/κ N = 1.7 ± 0.05. The spike feature right at 0 mV is an artifact due to very small currents and the resulting uncertainty in the fitting of the I-z curves.
The value of κ/κ N larger than 1 demonstrates that the tunneling process involves Cooper pairs, like the case of MAR and Josephson tunneling. The expected biases for single Andreev reflection and Josephson tunneling are e|V b | ~ Δ s ≈ 1.35 meV and e|V b | ~ 0 meV, respectively [15], so we can tentatively assign the plateaus in κ/κ N accordingly. However, κ/κ N = 1.7 for single Andreev reflection deviates from the expected value of 2 [21,22], and so we need further confirmation that this value is indeed from Andreev reflection and explore the origin of the deviation. In a control measurement of κ/κ N with Pb tip -Au(111) junction, κ/κ N = 1 outside the superconducting gap but rose to 1.2 inside the gap. This observation qualitatively supports our assignment of 1 < κ/κ N < 2 for single Andreev reflection, and that the double step feature in Pb tip -Pb(110) junction is not an artifact of the Pb coated tip. To gain complementary insight into the crossover of tunneling mechanisms for different κ/κ N regimes, we observe the distinct behavior of Andreev reflection and Josephson tunneling with respect to disorder. It is known that MAR strongly depends on the detailed atomic structure of the junction [8,26], while Josephson tunneling does not because Cooper pairs are delocalized over the coherence length (~ 80 nm in Pb) [27]. Therefore, they should be discernable via spatial distribution in heterogeneous samples.
To probe the feasibility of differentiating the mechanisms via disorder, we acquired the distance dependent spectroscopy in a wide area of the Pb(110) surface with natural population of defects [ Fig.   3(a)]. The STM topograph displays two types of surface defects, labeled d1 and d2, and one subsurface defect that generates scattering pattern around the center. Figure 3 highly depends on the atomic structure of the defects, as demonstrated in the κ/κ N of defect d1 whose value is lower than 2 for all bias. However, as V b approaches 0 mV, κ/κ N in all positions converges to 2, whose spatial uniformity supports the transition from MAR to Josephson tunneling.
For statistical analysis of the effect of disorder, we then acquired I-V-z curves on the regularly spaced grid.
The maps of κ/κ N at specific bias in Fig. 3(c-f) and their corresponding histograms in Fig. 3g  with still relatively broad distribution (see supplemental movie in [25] for the full dataset of κ/κ N maps and histograms). Note that there was a tip change at about 2/5 point in y axis of the κ/κ N maps, but the qualitative behavior of κ/κ N that we described above does not change.
The crossover of mechanisms can be effectively captured by calculating the Kullback-Leibler (KL) divergence between the probability density of distributions of the κ/κ N values at varying bias and the maximum bias, such as the ones in Fig. (g) from the numerical approximation of κ/κ N histograms, after centering them at the mean value [ Fig. 3(h)] [28]. The remarkably sharp transitions as a function of bias are evident. Most notably, the transition between MAR and Josephson regime occurs at |V b | ~ 0.5 mV in the KL divergence as a steep rise, which is also consistent with the bias where κ/κ N maps become uniform and loses atomic features [25]. We note that in the calculation of the KL divergence we intentionally subtracted the mean value of the distributions, so as to focus the comparison on the variation rather than averaged value of the κ/κ N . Thereby, the statistical analysis provides independent confirmation of our hypothesis of the distinct behavior of MAR and Josephson current over disorder and displays the transition in κ/κ N distribution more clearly. To simulate Andreev reflection as a function of basic properties of a tunnel junction, we utilized Kwant code for tight-binding transport calculations [29]. Following a basic algorithm, the Bogoliubov-de Gennes particle-hole symmetric Hamiltonian was implemented in a spinless system without magnetic field on a square lattice. The system is schematically shown in Fig. 4(a). The scattering region is localized between columns 3 and 4, with the superconducting gap being non-zero in columns 4-10 (blue dots), and the normal lead in columns 1-3 (yellow dots). The shaded dots on each side denote the beginning of infinite leads attached to each region, with the same corresponding Hamiltonians as in the regions. The conductance G calculated for such geometry clearly exhibits a superconducting gap and resonant Andreev reflection on the gap-edge [ Fig. 4(b)]. Subsequently, we have systematically varied the properties of the tunneling barrier to calculate conductance as a function of barrier properties. The analysis of the apparent decay constants was then carried out similarly to the experiment.
As seen in Fig. 4(c), the κ/κ N exhibits a broad range, from near 0 to about 3, across the superconducting gap. The values near 0 originate from resonant Andreev tunneling at the gap-edge, where the transmission probability is enhanced as is common for generic resonant tunneling. The effect of the resonance decays as the energy decreases toward the middle of the gap. Eventually κ/κ N exceeds 2 and reaches values as high as 3. The results for Andreev reflection are similar for the cases where we tune the barrier height and barrier width, although the resonance effect is most pronounced in the barrier-height dependent calculation. It is also important to consider the effect of broadening on these results, which were calculated at 0K. A simplest broadening of ~ 2kT applied to the calculated κ/κ N reveals that the resonance-mediated decrease below 1 is quickly smeared [ Fig. 4(d)]. However, the effect of the resonance is still present, producing intermediate κ/κ N between 1 and 2, over the width of the about half-the gap.
Overall, these results are qualitatively similar to the experimental observations in Fig. 2(c).
Both larger and smaller than 2 values for κ/κ N in the Andreev regime allows us to propose a simple conceptual picture of the crossover of tunneling mechanisms as three parallel channels whose resistance depends on tip height z and the tunneling energy e|V b |. When e|V b | > Δ s + Δ t and the density of states is considerable, single quasiparticle tunneling is dominant even though the MAR and Josephson tunneling is possible. Meanwhile, in the region of the superconducting gaps, the preponderance of Josephson vs Andreev tunneling will be determined by their decay constant, with the smaller decay constant determining the dominant current. Just below the gap, the resonance effect favors single Andreev reflection. Toward zero bias, higher order MARs appear successively, but also Josephson tunneling starts to compete with MAR, and eventually crossover to Josephson tunneling happens because MAR has higher decay constant than Josephson tunneling. The crossover will be sensitive to many specifics of the tunneling junction and the measurement conditions. For example, the role of prefactors in exponential decay of tunneling current remains to be understood in future analysis. Conversely, detailed understanding of the crossover is likely to provide a new window into the properties of impurities in superconductors.
In summary, we combined the I-V-z spectroscopy of the STM-based superconducting tunnel junction with the statistical analysis over atomic disorders to provide the clear illustration of crossovers between single quasiparticle tunneling, MAR, and Josephson tunneling. The analysis enabled one to determine the particular tunneling regime independently of the spectral shapes, which provides a valuable complement to other methods of analysis of superconducting junctions. In particular, rigorously comparing and correlating the effects of disorder in different regimes will deepen our understanding on the nature of superconductivity. For example, revealing intrinsic modulations of the MAR and Josephson currents would provide a pathway to identify pairing symmetry [12], inhomogeneity of superfluid density [16,30,31], and possible existence of exotic quasiparticles [32,33], at the heart of modern quest for quantum and topological computing.