Precursors of Majorana modes and their length-dependent energy oscillations probed at both ends of atomic Shiba chains

Isolated Majorana modes (MMs) are highly non-local quantum states with non-Abelian exchange statistics, which localize at the two ends of finite-size 1D topological superconductors of sufficient length. Experimental evidence for MMs is so far based on the detection of several key signatures: for example, a conductance peak pinned to the Fermi energy or an oscillatory peak splitting in short 1D systems when the MMs overlap. However, most of these key signatures were probed only on one of the ends of the 1D system, and firm evidence for an MM requires the simultaneous detection of all the key signatures on both ends. Here we construct short atomic spin chains on a superconductor—also known as Shiba chains—up to a chain length of 45 atoms using tip-assisted atom manipulation in scanning tunnelling microscopy experiments. We observe zero-energy conductance peaks localized at both ends of the chain that simultaneously split off from the Fermi energy in an oscillatory fashion after altering the chain length. By fitting the parameters of a low-energy model to the data, we find that the peaks are consistent with precursors of MMs that evolve into isolated MMs protected by an estimated topological gap of 50 μeV in chains of at least 35 nm length, corresponding to 70 atoms.

We can apply this method to the spectroscopic data presented in the main text. In Extended Data Fig. 1a, we provide additional dI/dV maps similar to the ones of Fig.2c of the main text but acquired at larger energies much closer to the gap edge Nb of the substrate. They reveal additional states with an integer number of maxima . These states are related to other multiorbital YSR states, most probably the -YSR state 1,4 , which is different from the δ-YSR state investigated in the main text, as we will substantiate in the following. In order to extract the dispersions of all sub-gap bands, we calculate the line-wise absolute squares of the Fourier transformations of dI/dV line-profiles similar to Fig. 3d-f of the main text. We sum over MnN chains of various lengths from 10 ≤ ≤ 32 to collect as many points as possible. The result can be seen in Extended Data Fig. 1b. Most notably, close to F , the FFT is dominated by a single faint band crossing F at F ≈ ± 2 ⁄ , which we can assign to the δ-YSR band. This again supports the single-band low-energy model described in the main text and the choice of the Fermi wavevector for this band. All other dispersing structures remain at higher energies. In particular, the states shown in Extended Data Fig. 1a belong to a relatively flat band that stays close to the gap edge Nb as it can be seen at energies 0.5 meV < | | < 1.5 meV in Extended Data Fig. 1b. Most probably it derives from hybridization of the -YSR states as can be seen in Fig. 3g of the main text, where the features within this energy range gradually emerge from the -YSR state of the individual atoms (cf. Fig. 1a of the main text) with increasing chain length. Moreover, these features have a similar particle-hole asymmetric spectral weight as the single atom -YSR state. These states are also responsible for the clear band formation in Mn chains on Nb(110) along the [001] direction 1 where this band crosses the Fermi level and opens a topological gap in contrast to the chains along the [11 ̅ 0] direction investigated here. It is plausible that the bandwidth of this band is reduced in chains along [11 ̅ 0] with respect to the [001] chains, since the inter-atomic distance is enhanced with respect to that in the chains along [001]. Since this band is confined to an energy window of 0.5 meV < | | < 1.5 meV (cf. Extended Data Fig. 1b and Fig. 3g of the main text) and thus well separated from F it does not contribute to the topological classification of the chains along [11 ̅ 0] within a low-energy effective model. Thus, by choosing a different chain geometry, we were able to controllably remove the influence of the -YSR band on the low-energy physics of our atomic chains.
In addition to the δand -YSR states discussed in the main manuscript text and here, there are two other additional less intense YSR states of the single Mn atom 1,4 , see e.g. the weak resonances at ±0.8 meV in Fig.1a. Similar to the -YSR states, these additional YSR states can be disregarded for the low-energy effective model here for two reasons. Firstly, they are much farther away from F as compared to the δ-YSR state. Secondly, unlike the δ-YSR state, they do not have lobes along [11 ̅ 0]. Therefore, in chains built along [11 ̅ 0], these additional states aren't expected to hybridize sufficiently strong to form Shiba bands which cross F .

Supplementary Note 2 | Magnetic structure of the chains
In Extended Data Fig. 2a, we show dI/dV maps of a Mn11 chain along [11 ̅ 0] obtained with a spin-polarized Cr tip. While the tip magnetization has been shown to be stable using a magnetic reference chain presented in Ref. 5, the dI/dV signal on the chain changes with an applied out-of-plane field z . This indicates that the chain's magnetization aligns with z . We find a homogeneous contrast along the chain, indicating a ferromagnetic alignment of the magnetic moments. Furthermore, we repeated this measurement using a superconducting Nb tip functionalized with Fe atoms, inducing YSR states on the tip. We can exploit the spin-polarization of the YSR states to map the magnetic structure of our sample 5 . The result is shown in Extended Data Fig. 2b, measured on a Mn15 chain along [11 ̅ 0]. The YSR-asymmetry parameter YSR,norm ( ⃗ t ) at the tip position ⃗ t is defined in Ref. 5 and resembles the sample's spin-polarization. Again, we find a finite and spatially homogeneous spin-contrast, confirming the ferromagnetic structure of the chains. The observation that the magnetization can be aligned by a relatively small out-of-plane magnetic field further indicates that the easy axis of the chain is perpendicular to the surface.

Supplementary Note 3 | Additional model calculations for YSR chains
Studying the model introduced in the Methods section of the main text, we find that with the parameters used in Fig. 4 of the main text, the chain is deep in a topological phase, although with a comparably small topological gap of p = 80 µeV (see Fig. 1e of the main text). Thus, we expect isolated MMs to form at the chains' ends for → ∞. We study this limit theoretically in Supplementary Fig. 1.
In Supplementary Fig. 1a, we plot the energy resolved LDOS along a chain of N = 200 using the parameters from Fig. 4 of the main text. This shows that the exponential decay length of the Majorana wavefunction is long, which explains their significant residual interaction in experimentally accessible chains of N < 50. In Supplementary Figs Supplementary Fig. 1c. The temperature broadening was reduced to an equivalent of = 200 mK in this plot for better visibility. In Supplementary Fig. 1d, we depict the lowest energy wavefunction for N = 200, clearly illustrating the exponential decay of the MM for long chains. Although the MM has significant weight especially on the chains' end sites, there is strong spectral weight spread over the first ~70 sites. Furthermore, we can compute the MM energy splitting vs. chain length, as plotted in Supplementary Fig. 1e. As visible in the plot, the MMs converge to energies below the topological gap , which is defined by the energy of the energetically lowest bulk state in infinite chains, for N > 70.
Note, however, that the decay constant strongly depends on the strength of the Rashba spin-orbit coupling (SOC) term h , which cannot be directly measured in the experiment. The effect of the Rashba SOC strength h is presented in Supplementary Fig. 2 for chains up to N = 180 sites. For the exact value h = 0.00, only conventional finite-size quantized bulk states are found to be filling the gap. Here, the topological index is not well defined since the system is entirely gapless. However, for any arbitrarily small h , the chain is in the topologically nontrivial phase. As it is visible in Supplementary Fig. 2, the lowest energy states separate more and more from the higher energy states with increasing h . This separation is even slightly visible for very small values (e.g. h = 0.01 π/d in Supplementary Fig. 2a), in agreement with the bulk-boundary correspondence. For larger h , the topological gap increases and the convergence of the PMMs to zero energy is faster. Since we cannot directly determine the exact value of h from experimental data, conclusions about the chain length needed to observe isolated MMs are challenging to draw. However, we can rule out the case of h = 0.00, since previous studies of the Mn/Nb(110) system found that Rashba SOC in this system is significant and can induce observable effects like splittings of the YSR states in antiferromagnetic dimers 4 or gaps in a YSR band of a Mn chain 1 . Therefore, isolated MMs are clearly expected to occur in long Mn chains. The observations of Supplementary Fig. 2 are in agreement with the expected inverse relationship between and the decay length of the Majorana wavefunction 6 : namely, for larger values of , the modes converge faster towards zero-energy. This is further validated by analyzing the long-range decay of the wavefunction's envelope. For this, we fit the envelope of the lowest-energy wavefunction in a chain with N = 800 sites to the expression with the decay parameter and a normalization factor . The result for the decay length −1 is shown in Supplementary Fig. 3b and clearly anti-correlates with the magnitude of the topological gap, i.e. −1 diverges when the gap closes and is small for large gaps. This further demonstrates that the topological nature of the end states manifests itself in the long-range decay of the wavefunction instead of the observed apparent localization on the terminal sites only.
Another interesting observation is the presence of avoided crossings between the lowest energy state and the higher energy excitations ( Supplementary Fig. 2). The degree of separation from the trivial states increases with h and could thus potentially serve as an estimator for the pairing magnitude if the topological gap itself is not resolved in short chains with PMMs. In order to reduce the MM interactions in future experiments, it would thus be beneficial to enhance the topological gap . Surprisingly, an increase of the Rashba SOC parameter h does not necessarily lead to an increase in . This is shown in Supplementary Fig. 3a, where the system's gap is plotted against h and F,0 . Supplementary Fig. 3b shows a line-cut through the gap diagram in Supplementary Fig. 3a at F,0 = 0.53 ⁄ , which was found to match the experiment. As already seen in Fig. 1e of the main text, there is a gap closing between two topologically non-trivial regions at F,0 ≈ 0.5 / . The reason lies in the k-dependence of the pwave pairing strength ( ) in our model, as it will be elaborated in the following. Interestingly, one finds that ( ) features zero crossings at = F,0 in the case of vanishing non-magnetic scattering = 0 7 . Thus, no topological gap can be opened by the pairing term if the Fermi wavevector F defined by the crossing of YSR bands through F equals the Fermi wavevector of the host, i.e. F = F,0 . In Supplementary Fig. 3c, we show the k-dependent superconducting pairing ( ), and the k-dependent hopping terms ℎ are shown in Supplementary  Fig. 3d for different magnitudes of the SOC parameter h . Notably, the bands ℎ indeed cross zero at F ≈ F,0 , i.e. where the pairing is zero. This explains why the topological gap in Supplementary Fig. 3b is small, albeit the pairing amplitude in Supplementary Fig. 3c reaches values of up to 1 meV. One can derive analytically 7 that F = F,0 for = 1, = 0 (i.e. for hybridizing YSR states, where an isolated state would be at exactly zero energy) and F,0 = 0.5 ⁄ . Thus, no band gap will open in this case (see also Fig. 1e of the main text). Since we have strong evidence that F,0 ≈ 0.5 ⁄ , it is reasonable that the system presented in the main text is close to the gapless condition F = F,0 . This could explain the small magnitude of the p-wave pairing in our chains although there has been evidence that the Rashba SOC strength in Nb is significant 4 and an additional, effective Rashba term can be introduced by the inevitable inversion-symmetry breaking at the surface where the chain is located (see Supplementary Note 5). By this, an additional constraint is introduced for realizing strongly protected, i.e. isolated MMs in future experiments.

Supplementary Note 4 | Overlap of the MMs with increasing chain length
As shown in Figs. 2b,c of the main text, localized zero energy states are found at the chain's ends for N = 32 whose phenomenology can be reproduced by the model for a single YSR band in Fig. 4b. These PMMs at approximately zero energy presented in the main text feature several similarities with isolated MMs, such that they could be misinterpreted as strongly protected MMs: They have a good localization on the chain's ends, near-zero energy and are separated from higher energy excitations by a gap FS . The important difference is, that these states cannot be decomposed perfectly into spatially separated MMs localized on either chain ends 8 Supplementary Fig. 4b is not a result solely from the topological nature of the end state. This is the reason for their weak protection against changes in the chain length. In contrast, the mode for N = 200 in Supplementary Fig. 4c can be separated into two well isolated MMs on the respective chain ends (Supplementary Fig. 4f). Interestingly, the separation ( Supplementary Fig. 4d) for the finite-energy state in Supplementary Fig. 4a is better than that for the near-zero-energy state in Supplementary Fig. 4b. Although counterintuitive, the oscillating Majorana overlap and the energy of the Dirac state have been shown to anticorrelate in one-dimensional topological superconductors in the short wire limit 9,10 .

Supplementary Note 5 | Majorana localization by inversion symmetry breaking
The good localization of the PMMs despite their strong, chain-length-dependent energy splitting is counterintuitive. In order to understand this behavior, we use a tight-binding toy-model approximating the lowenergy band structure of a one-dimensional ferromagnetic chain on two layers of a superconducting substrate, in addition to the effective one-band model mentioned above and in the main text. The calculations were performed using the Kwant code 11 . This model thus reproduces the correct geometry of the experimental setup, where the chain is located on the surface of the superconductor and not in the bulk. The Bogoliubov-de-Gennes Hamiltonian for the chain reads: using the basis = ( ↑ , ↓ , ↓ † , − ↑ † ) T . Here, c is the magnetic interaction strength, c is the hopping parameter in the chain, c is the proximity-induced superconducting order parameter in the chain, c is the strength of the Rashba SOC and c is the chemical potential in the chain. and are the Pauli matrices in particle-hole-and spin-space, respectively, and ⊗ represents the Kronecker product. < , > indicates the summation over nearest neighbor sites and , † and , are the creation and annihilation operators of (Dirac-) fermionic states with spin at site i of the chain.
Two layers of a cubic substrate lattice with 70 x 15 sites are additionally considered, and hopping between the chain and the substrate is taken into account, Similar to above, s is the hopping parameter in the substrate, s is the superconducting order parameter in the substrate, s is the chemical potential in the substrate and cs is the hopping parameter between chain and substrate.
We compute the energy-dependent LDOS at site x by diagonalizing the Hamiltonian ℋ = ℋ chain + ℋ substrate + ℋ hop and summing over all pairs of eigenvalues and eigenvectors : with the Fermi-Dirac function ( , ) simulating the experimental thermal broadening. The summation over both the particle-and hole-components of the solutions in this section is justified by the value of ( , ) ≈ 0.5 for the -YSR states studied in the main text.
We start by reproducing the even-odd modulation of the sub-gap states using the following parameters, which result in a band structure with F ≈ 0.5 ⁄ for c = 0: Calculating the topological invariant of the system 12,13 for these parameters, we find that it is non-trivial for all non-zero c . Adding an explicit Rashba SOC term c = 0.055 s = 0.42 s in the chain, the band structure becomes gapped (see Supplementary Fig. 5a). This corresponds to a Rashba SOC parameter of ℏ = 2 c = 0.0059 • Å, in good agreement with the model used in the main text. As shown in Supplementary Fig. 5b, we find localized PMMs in a chain of N = 32. However, similar to the model for YSR bands used in the main text and Supplementary Note 3, these PMMs split strongly in energy for chains with ≤ 45, as shown in Supplementary  Figs. 5c,d. This splitting is modulated with a period of ∆ ≈ 2, yielding the same even-odd phenomenology as in the experiment.
In addition, we compare the localization of the edge modes for different chain-substrate couplings. If we set cs to zero, the chain is effectively one-dimensional and the substrate does not contribute to the LDOS on the chain. We quantify the localization of the low-energy modes by the ratio of the LDOS on the first site and the averaged LDOS in the chain's bulk, of a chain with N = 32. The LDOS is symmetric, hence it is equivalent to use the LDOS at the last site of the chain instead of the first site in the numerator of equation (S6). In Supplementary Figs. 6a,b, we compare the wavefunctions of the low-energy modes for the effective 1D model ( cs = 0) and the extended model using a two-layer substrate ( cs = 0.7 s ). We find that the localization for chains coupled to a substrate is consistently better than the localization in the purely one-dimensional model (see Supplementary Fig. 6c). Similar findings were reported in a model for Fe chains on Pb 13 . Surprisingly, when including the substrate, the lowest energy modes are even localized for a vanishing explicit Rashba SOC term c = 0. This is explained by the broken inversion-symmetry at the surface 14 .
In the model used in the main text and Supplementary Note 3, the localization of the observed features can be explained by effective long-range hopping and pairing terms in the Hamiltonian, which are a result of an integration over the three-dimensional bulk of the superconducting host. Hence, the short-range decay of the MM wavefunction cannot be described as an exponential decay but is more complex 7 (see Supplementary Notes 3,4). Interestingly, in the limit of small coherence lengths and weak Rashba SOC, all interactions are local and the LDOS corresponding to the low-energy states only slightly deviate from plane waves. Here, the YSR-band model converges to the one-dimensional case that is shown in Supplementary Fig. 6a.

Supplementary Note 6 | Estimating the Fermi wavevector kF,0
One input parameter for the model described in the Methods section of the main text and Supplementary Note 3 is the Fermi wavevector F,0 of the substrate hosting the YSR atoms. We can try to estimate the value of F,0 by probing the wavefunction of the single atom YSR state. The wavefunction has been calculated to oscillate far away from the impurity with ( ) ∝ sin( F,0 ) 15,16 . Since the measured LDOS is proportional to the absolute squares of the wavefunction, it is expected to be modulated with = 2 F,0 . In Extended Data Fig. 3a, we show an example of a dI/dV map around a single Mn atom acquired at the energy of the + -YSR state 4 . We computed the Fourier transform (FFT) of the map in Extended Data Fig. 3a with the central part removed to avoid effects from the orbital structure of the YSR state. The result is shown in Extended Data Fig. 3b and exhibits peaks at ≈ (1.2 ± 0.2) ⁄ and thus F,0 ≈ (0.6 ± 0.1) ⁄ , compatible with ± 2 ⁄ as proposed by our results on Mn chains in the main text.

Supplementary Note 7 | Calculation of the YSR band structure within a single-band model based on first-principles simulations
We also determine the band structure of the YSR states using the model described in Ref. 4, which previously has been applied to Mn adatoms and dimers on Nb(110). The Hamiltonian of the system reads where ∥ denotes the wave number inside the one-dimensional Brillouin zone of the chain oriented along the [11 ̅ 0] direction. The three-dimensional wave vector in the reciprocal space of bulk Nb is given by = ( ⊥ , ∥ + ∥ ), where ⊥ is the perpendicular component of the wave vector and ∥ is a reciprocal lattice vector of the chain. The individual terms in the Hamiltonian represent the spectrum of bulk Nb nearest-neighbor Rashba-type SOC on the bcc(110) surface of Nb s-wave superconducting pairing and the impurity potential Note that ∥ is fixed in all terms and the summations are restricted to ⊥ and ∥ . In Eqs. (S7)-(S11), Φ = 1/√2( ↑ , ↓ , − ↓ † , − ↑ † ) is the Nambu spinor, with the electron annihilation operators . The matrices in Nambu space read where are the usual 2 × 2 Pauli matrices for = x, y, z and 0 is the unit matrix. The bulk Nb spectrum , the Rashba parameter R and the scattering potential shapes Ψ ( ) were determined based on ab initio simulations in Ref. 4, where details about the procedure can be found. A finer mesh is used in reciprocal space compared to Ref. 4 due to the restriction to fixed ∥ values because of the infinite chain. The order parameter Δ = 1.5 meV is based on the experimental value Nb . In the impurity potential, is the spin quantum number, ( , , ) is the direction of the localized magnetic moment, and indexes the different scattering channels. The parameters and /2 are determined by an approximative procedure for the adatom discussed in Ref. 4, and the same values are used here for modelling the chain with and without SOC. The chain was assumed to be ferromagnetic along the direction, and the calculations were restricted to the -YSR state.
The Green's function at complex energy is expressed as ( , , ′ ) = 0 ( , ) , ′ + 0 ( , ) ( , , ′ ) 0 ( , ). ( 13) Note that ( , , ′ ) and ( , , ′ ) are diagonal in the wave number of the chain ( ∥ = ∥ ′ ) due to the periodicity in real space, while the Green's function of the substrate is diagonal in ⊥ as well. The operator is defined as where 4 is the unit matrix in Nambu space.
The Shiba bands were determined by finding the poles of the operator, i.e. the real energy values | | < Δ for which the determinant of the matrix in Eq. (S15) vanishes. The results are summarized in Supplementary Fig. 7.
Note that due to the particle-hole symmetry of the Hamiltonian, the YSR states always appear in pairs at positive and negative energies, and solid lines connect only the positive-energy or the negative-energy solutions as a guide to the eye. The minimal distance in energy between positive-and negative-energy states 2 gap may be interpreted as an upper bound for the induced gap in the system, limited by the resolution in ∥ . We found gap = 3 µeV in the absence (Supplementary Fig. 7a) and gap = 53 µeV in the presence of SOC ( Supplementary  Fig. 7b), respectively. The inner gap is more than a factor of ten higher when SOC is taken into account, and the spectrum is also flatter close to the minima than without SOC. This indicates that in the absence of SOC, the gap will close as the resolution in ∥ is increased, but SOC may open an induced gap in the system. The magnitude of the gap here, which is based on first-principles calculations of the Rashba parameter, is comparable to the one obtained from the effective model in the main text (see also Supplementary Note 3).
were calculated on the positions of the atoms in the chain, for each state indexed by ∥ at the YSR energies determined above. In Eqs. (S16) and (S17), δ denotes a small imaginary part, which was set to B ⋅ 300 mK in order to approximate the experimentally observed energy resolution. The spin polarization of the states is defined as the ratio of the MDOS and the LDOS. In the absence of SOC, the YSR states are expected to be fully spin-polarized as → 0, with opposite spin polarizations of the particle and hole states. Therefore, a sign change of the spin polarization when looking only at the positive energy states indicates a YSR band crossing the Fermi level. In the presence of SOC, the particle and hole states should hybridize, but far from the avoided crossings the sign of the spin polarization may still be used as an indication for whether a band crossing has occurred. This is illustrated by the color-coding in Supplementary Fig. 7. With SOC, the spin polarization changes sign at around ∥ = 0.35 / , exactly at the position of the avoided crossing in Supplementary Fig. 7b. This value is not very far from the value ∥,F = 0.53 / used in the low-energy model calculations in Fig. 4 of the main text, taking into account that no fitting parameters were used for the chain calculations described here. Although the positiveand negative-energy bands are the closest for ∥ = 0.97 / in Supplementary Fig. 7b, there is no sign of an inversion of the spin polarization in this regime. In summary, the calculations indicate that the YSR bands cross the Fermi level only once for ∥ > 0, and this crossing is gapped out by the SOC. This supports the interpretation of the chain being in the topologically non-trivial regime.

Supplementary Note 8 | Interacting chains
Since we are constructing the chains by single atom manipulation, it is possible to create more complex nanostructures. As an example, we construct a junction of two MnN chains and control their interaction by the number of unoccupied sites ∅ between them (Fig. 5 of the main text). The full results for 0 ≤ ∅ ≤ 5 are shown in Extended Data Fig. 4. At large distances ( ∅ > 5), the chains do not interact strongly. Notice that the dI/dV line-profiles for the two respective chains with ∅ > 5 are extremely similar, underlining the excellent reproducibility of our results. Moving the chains closer together results in small energy changes of the end state. Most notably, the energies of the left and right end states of the same respective chain change equivalently and symmetrically. This proves that it is a collective mode of the chain instead of a local electronic state. For ∅ = 1, we observe additional in-gap features between the chains which stem from hybridization between the states on the individual chains. This happens in particular at energies around ≈ 800 µeV and demonstrates that the interaction between the two chains is significant for these distances. Such additional in-gap states could hamper the identification of MMs in networks of chains if the experimental energy resolution is not sufficient. Finally, for ∅ = 0, the two chains merge into one Mn24 chain. Here, there is apparently no excessive dI/dV signal at the chain junction but only at the remaining two chain ends. However, the length of N = 24 accommodates PMMs with a finite energy splitting, as already discussed in detail in the main manuscript. Interestingly, while the higherenergy states appear to be spatially inhomogeneous (see ≈ 600 − 800 µeV), which we attribute to local variations in the substrate's chemical potential, the lowest energy mode remains symmetric in space, further indicating that it is a single quantum state. All these observed effects are in agreement with the expectations for PMMs.

Supplementary Note 9 | Deconvolution of dI/dV spectra taken with superconducting tips
In all experiments presented in the main text and the Supplementary Information, W tips covered with clusters of superconducting Nb were used (see Methods). While this is known to increase the experimental energy resolution, the dI/dV spectra measured with superconducting tips are not, as usually valid for metallic tips, We compute the inverse A -1 of A and multiply the result with the vector [dI/dV] SIS , which approximately yields s ( ). In the deconvolution process, the value T = 320 mK is used, as read out by a temperature sensor on the 3 He pot, which is thermally well coupled to the STM unit 17 (see Methods). In addition, the deconvolution can only be performed if t ( ) is known.
Here, a broadened Dynes density of states with t = 1.42 meV is assumed for the superconducting Nb cluster on the tip: .
( 20) 0 is the normal conducting DOS, Re is the real part, and is the Dynes broadening parameter. = 0.01 meV was chosen for all tips used throughout this work, as validated by the observation of negative dI/dV features (see Supplementary Fig. 8), which are known to appear in the differential conductance when tunneling between very sharp features in s ( ) and t ( ). The values for t and do not change drastically when the tip apex is modified (e.g. by indenting the tip several nm into the Nb surface), thereby suggesting that the Nb coating of the tip is stable and thick.
The value of t can be extracted from measurements of multiple Andreev reflections (MARs), which appear as very weak features in dI/dV spectra already at comparably large tip-sample distances. This can be seen in Supplementary Fig. 8: At = ±( s + t ) = ±2.87 meV , the convolution of the sample's and tip's coherence peaks yields a large peak in dI/dV. Zooming into the central gap region, we find additional peaks at ±1.51 mV and ±1.36 mV, representing Andreev tunneling occurring at = ± s and = ± t for this particular tip 18 . Since the peak at 1.36 mV is observed to slightly change in energy for different microtips, we identify it with t and the peak at 1.51 mV with s , which is in good agreement with the bulk gap of Nb. For the tip used in the main text and the Supplementary Information, ±( s + t ) = ±2.93 meV and thus t = 1.42 meV is found. This value was used for equation (S20) and the numerical deconvolution process. Note that the MAR peaks are about three orders of magnitude smaller than the peaks at ±( s + t ), indicating that single-particle tunneling is the dominant transport channel. Therefore, MARs do not need to be considered in the deconvolution process.