Spin-dependent tunneling between individual superconducting bound states

Magnetic impurities on superconductors induce discrete bound levels inside the superconducting gap, known as Yu-Shiba-Rusinov (YSR) states. YSR levels are fully spin-polarized such that the tunneling between YSR states depends on their relative spin orientation. Here, we use scanning tunneling spectroscopy to resolve the spin dynamics in the tunneling process between two YSR states by experimentally extracting the angle between the spins. To this end, we exploit the ratio of thermally activated and direct spectral features in the measurement to directly extract the relative spin orientation between the two YSR states. We find freely rotating spins down to 7mK, indicating a purely paramagnetic nature of the impurities. Such a non-collinear spin alignment is essential not only for producing Majorana bound states but also as an outlook manipulating and moving the Majorana state onto the tip.

Magnetic impurities on superconductors induce discrete bound levels inside the superconducting gap, known as Yu-Shiba-Rusinov (YSR) states.YSR levels are fully spin-polarized such that the tunneling between YSR states depends on their relative spin orientation.Here, we use scanning tunneling spectroscopy to resolve the spin dynamics in the tunneling process between two YSR states by experimentally extracting the angle between the spins.To this end, we exploit the ratio of thermally activated and direct spectral features in the measurement to directly extract the relative spin orientation between the two YSR states.We nd freely rotating spins down to 7 mK, indicating a purely paramagnetic nature of the impurities.Such a non-collinear spin alignment is essential not only for producing Majorana bound states but also as an outlook manipulating and moving the Majorana state onto the tip.
Magnetic impurities on superconductors give rise to a wealth of phenomena due to the subtle interplay between pairing processes of conduction electrons and their exchange interaction with the impurity [1][2][3][4][5][6].As a consequence, single spin-nondegenerate bound states inside the superconducting gap may emerge [7][8][9][10][11][12][13][14][15][16].A fundamental understanding of the role of spin in the coupling between such states is highly relevant [17], for example, in the context of creating Majorana bound states at the ends of a chain of Yu-Shiba-Rusinov (YSR) states [18,19].Aside from the conventional quasiparticle hopping from one site to the next, pair creation from the superconducting condensate can also couple two neighboring sites in a chain of YSR states.To reach the topological regime in such a chain, the coexistence of both types of coupling processes between neighboring sites is indispensable [20,21].
e spin provides a vital component selecting between the two coupling processes.Conventional hopping is favored when the spins of two sites are aligned, while pair creation is favored when neighboring spins are anti-aligned.
Past experiments observed end states on magnetic atom chains on superconductors using scanning tunneling microscopy [2,18,19], but the two tunneling processes between neighboring sites in the chain could not be individually accessed.To separate the two tunneling processes, control over the energy detuning between adjacent sites is needed.For this purpose, instead of using the tip as an external sensor [19,22], measuring a YSR state with a functionalized YSR tip turns the tunnel junction into an internal probe of the coupling between two YSR sites.
Here, we use a scanning tunneling microscope (STM) [23] to couple a YSR state on the vanadium tip and an intrinsic YSR state on the V(100) surface, where we independently resolve direct and thermally activated Shiba-Shiba tunneling (Fig. 1(a-c)).We rst show experimentally from the spin blockade of a characteristic Andreev re ection that an indi-vidual state is spin-nondegenerate.Measuring the two tunnel processes between two spin-nondegenerate superconducting bound states, we nd that both direct and thermal processes are present for all measured impurities.eir coexistence excludes the possibility that the impurity spins are collinear, ful lling the requirement for use in topological chains.e measured relative strength of the two tunnel processes enables us to extract the angle  of the relative spin orientation.Tracing the relative spin orientation from 1 K to 7 mK, we nd  = 90 • in all measurements indicating freely rotating spins of the YSR states in tip and sample with no detectable anisotropy.At the end, we draw an analogy between Shiba-Shiba tunneling and the coupling of neighboring states in a YSR chain.
When both tip and sample exhibit a YSR state, the tunneling between these discrete levels can be observed in the current-voltage  ( ) characteristic as isolated spectral features [24].At 10 mK, we only observe current peaks at bias voltages indicated by blue arrows labeled  ± in Fig. 1(d), which we refer to as direct Shiba-Shiba tunneling.is process is found at a bias voltage  of the sum of two YSR energies  = ±|  +  | ( is the electron charge,  , are YSR energies for the tip and sample, respectively).e corresponding process is plo ed schematically in Fig. 1(b).
At higher temperature (1 K), however, the YSR state may get thermally excited, such that thermally activated tunneling becomes possible.Figure 1(c) shows schematically how an excitation tunnels from one side to the other.Additional current peaks appear at the bias voltage  = ±|  −   | (see Fig. 1(e) indicated by orange arrows and labeled  ± ), which we call thermal Shiba-Shiba tunneling.
e intensity of the thermal Shiba-Shiba peaks is smaller compared to the direct Shiba-Shiba peak due to the small thermal Boltzmann factor at 1 K. Due to their appearance at distinct bias voltages, the direct and thermal Shiba-Shiba processes can be separately addressed in a single measured  ( ) spectrum.While the general dynamics of the tunneling processes of Shiba-Shiba tunneling has been discussed recently [24], the role of the spin in the transport process is more subtle, but plays a de ning part in the coupling between YSR states.
Conventional models describe YSR states as spinnondegenerate states, i.e., 100% spin polarized [8,26], but not necessarily with a preferred axis.By contrast, a dominant Kondo coupling with a Kondo temperature  K larger than the superconducting gap  ( K ) may lead to a spin-degenerate bound state.Experimentally, it should thus be veri ed that an alleged YSR resonance actually is spin polarized [27,28].Without resorting to magnetic elds, this can be checked by looking at the presence or absence of the single Andreev re ection process involving the same bound state [29] as shown in Fig. 2(a).e absence of such a spectral feature at a bias voltage of  = ± indicates that the bound state is a spin-nondegenerate YSR state [29].A high conductance spectrum of a single YSR state in the tunnel junction featuring multiple Andreev re ection peaks is shown in Fig. 2(b).At the expected bias voltage, the respective single Andreev process is absent in the measurement, i.e., there is no spectral feature at the voltage indicated by the symbol , which con rms the spin-nondegeneracy of the YSR state (for the assignment of all observed peaks and other details see the supplementary information [25]).
In order to understand the in uence of spin on the Shiba-Shiba peak intensities, we rst consider the limiting case of the YSR states being aligned either parallel ( = 0 • ) or antiparallel ( = 180 • ) as shown in Figs.3(a)-(d).If the spins of the corresponding parts of the YSR states are parallel, direct Shiba-Shiba tunneling is forbidden (Fig. 3(a)).is is due to spin conservation in the tunneling process and because the electron ( − ) and hole (ℎ + ) parts, which participate in the tunneling process, have opposite spins.e thermal Shiba-Shiba process, however, is allowed, because no spin blockade exists (Fig. 3(b)) for tunneling between two hole parts due to parallel spins.
Observing both direct and thermal Shiba-Shiba peaks in Fig. 1(e) directly leads to the conclusion that the spins in the measurement can neither be completely parallel nor com-FIG.3: eory showing thermal-direct Shiba-Shiba ratio as a measure of the relative spin orientation  .(a,b) e tip and sample YSR states have parallel spins ( = 0 • ).Note that the  − part and ℎ + part of the same YSR state always have opposite spins [26][27][28].Here, the direct Shiba-Shiba process is forbidden while thermal Shiba-Shiba tunneling is allowed due to spin conservation during tunneling.(c,d) e tip and sample YSR states have opposite spins ( = 180 • ).Here, the direct Shiba-Shiba process is allowed while thermal Shiba-Shiba tunneling is forbidden.(e) e continuous dependency of the direct and thermal Shiba-Shiba peak intensity on the relative spin angle.(f) e universal dependency of  / th on the relative spin angle.
pletely anti-parallel.Using the standard expression for the tunneling current explicitly including the spin [25] where  =↑, ↓ and σ =↓, ↑, the  ( ) spectra can be simulated expanding the above discussion to arbitrary angles between the spins of the two YSR states (0 • <  < 180 • ) (for details see supplementary information [25]).e result is shown in Fig. 3(e).e Shiba-Shiba peak intensities evolve smoothly between the extreme cases of parallel and anti-parallel orientations, con rming the limiting cases discussed before.In the measurement, we nd two peaks for each direct and thermal process (one at positive and one at negative bias voltages).To account for both peaks in the analysis, we average between them using the geometric mean which eliminates the electron-hole asymmetry and the absolute intensity of the peaks, so that the resulting ratio  only depends on the angle  between the spins of the two YSR states as well as the Boltzmann factors of the two YSR states where  is the temperature,  B is the Boltzmann constant,   with  =  ± ,  ± is the intensity of the corresponding Shiba-Shiba peak.Further, we de ne a parameter  th as the di erence of the two Boltzmann factors, which are experimentally accessible, such that we nd Note that here it is possible to use either current peak area or peak height, because the width of direct and thermal Shiba-Shiba is the same (for details see supplementary information [25]).e predicted evolution of  / th is plo ed in Fig. 3(f), which maps a unique angle for every ratio.If  / th → ∞, the spins are parallel; if  / th → 0, spins are anti-prallel.Equation (3) presents a universal curve that is independent of the details of each YSR state.Because the Shiba-Shiba peaks are typically well separated in the measured spectrum, both  and  th in Eq. ( 3) can be precisely derived from directly measurable quantities in a single measurement, such that the angle  can be experimentally determined.It is this ability to independently resolve the relative strength of two di erent Shiba-Shiba tunneling processes in the  ( ) spectrum, which enables an all-electrical measurement of the relative spin-orientation in our experiment -a striking di erence to conventional spin-valve experiments, where a polarizing magnetic eld is used to rotate (or ip) the relative orientation.
To extract the Boltzmann factors, we consider experimentally the thermal to direct ratio of the YSR state tunneling into the continuum [30].A statistics of such measurements at 1K is plo ed in Fig. 2(c) (for details see the supplementary information [25]).All Boltzmann ratios for both tip and sample YSR states lie on a line corresponding to 1 K ± 0.05 K. is agrees with our thermometer reading and also con rms that we have only one pair of YSR states inside the gap [30].
We have collected the experimental parameters  and  th from a number of di erent Shiba-Shiba systems in the linear tunneling regime using intrinsic defects on the V(100) surface and a V tip at 1 K.Note that in the tunneling regime the two-state system is minimally disturbed as higher order tunneling processes are suppressed.In Fig. 4(a), we plot the thermal-direct Shiba-Shiba ratio  against  th .Interestingly, for all Shiba-Shiba systems, the data points lie on the identity line ( =  th ), within a small error interval.According to Eq. (3), this means that the the relative angle between the two YSR spins is  = 90 • ± 3 • .Considering that every Shiba-Shiba system in Fig. 4(a) shows the same angle  , it is unlikely that these angles are time-independent.Instead, we are measuring the average angle of a time dependent spin orientation.e tunneling time is short compared to any motion of the spins, such that every tunneling event will observe a static snapshot of the system, so that the measured current is a time average over many independent tunneling events.erefore, we propose that the spins are freely rotating, giving an average e ect  =  th , where indicates averaging over an angle isotropically distributed over the unit sphere (for more details, see the supplementary information [25]).Note that although the spin is freely rotating, the electron and hole part of the same YSR state has opposite spin direction at any time.
A temperature of 1 K may be still too high for a magnetic system to develop an easy axis [31], which could explain the observed isotropy.We, therefore, trace the parameters  and  th as a function of temperature down to 0.6 K. Figure 4(b) shows the thermal-direct Shiba-Shiba peak ratio  for four di erent Shiba-Shiba systems between 0.6 K and 1 K. Upon cooling, the data points for the ratio  nicely follow the calculation from Eq. ( 2) (solid lines) assuming an angle  = 90 • . is indicates that  =  th with a still freely rotating spin scenario down to 0.6 K.
If we cool further down below 0.6 K, the ratio becomes more di cult to access because of the exponential suppression of thermal Shiba-Shiba peak.However, recalling Fig. 3(e), we note that the direct Shiba-Shiba peak intensity by itself critically depends on the relative spin angle  .To eliminate the conductance dependency as well as temperature e ects on the YSR lifetime and environmental broadening [24,[32][33][34][35], we use the peak area normalized by the conductance (see supplementary information for details [25]).Figure 4(c) shows a typical data set of the normalized direct Shiba-Shiba peak area from 1 K down to 7 mK, showing no signi cant change over two decades of temperature range.
ese measurements are consistent with a relative spin orientation rotating isotropically down to 7 mK without any anisotropy energy ( B  < 0.6 eV).We nd that for YSR states, the relative orientation of the spins determines whether direct or thermal tunneling is possible (cf.Fig. 3).
is is di erent from the spin valve e ect in normal conducting tunnel junctions, where an antialignment suppresses the tunneling current, but does not favor a di erent transport channel.
It is tempting to draw an analogy between the two coupled YSR states in our experiment and a link in a chain of YSR states [2].
e two types of coupling processes, direct and thermal Shiba-Shiba tunneling, are counterparts of the pair creation and hopping terms in the Kitaev Hamiltonian (see Fig. 4(d)) [20,24,25].In this analogy, Majorana end states require the coexistence of both coupling processes, which is guaranteed by the non-collinearity of the neighboring impurity spins.We surmise that the relative spin orientation between YSR states is equally important in the chain as in the tunnel junction.One di erence, however, is that in our experiment the two impurities are placed on di erent superconducting substrates relatively far away compared to the atomic chain.
is rules out magnetic interactions via the substrate, which have been suggested to induce helical order in Kitaev chains [36][37][38][39][40][41][42].Nevertheless, as an outlook, our results suggest an intriguing experiment to probe the edge of a magnetic impurity chain by a YSR tip.e tip YSR state will couple to the state at the edge of the chain by direct and thermal Shiba-Shiba tunneling and extend the chain by one site.
e Majorana end state will move from the substrate onto the tip.
In summary, we propose the ratio between the thermal and direct Shiba-Shiba tunneling peaks as a direct and unambiguous all-electrical probe for the relative spin orientation between YSR states.We conclude that for the Shiba-Shiba systems consisting of intrinsic YSR impurities on V(100) and YSR states on the vanadium tip, the YSR states derive from freely rotating quantum spin-1 ⁄ 2 impurities from 1 K down to 7 mK.e spin plays a de ning role in the transport between spin-nondegenerate superconducting bound states selecting between hopping and pairing mechanisms (i.e., direct and thermal Shiba-Shiba tunneling).Additionally, the analogy between the Shiba-Shiba junction and a chain link in a YSR chain a ributes a critical role to the spin for producing Ma-jorana end states.Further, a YSR tip in close proximity to a Majorana end state may induce a topological phase transition allowing the Majorana end state to be manipulated by the tip.Going beyond the immediate consequences, a Shiba-Shiba junction may also provide the basis for a minimal source to inject triplet Cooper pairs into a superconductor [

Tip and sample preparation
e sample was a V(100) single crystal with > 99.99% purity, and was prepared by standard UHV metal preparation procedure of multiple cycles of argon ion spu ering around 10 −6 mbar argon pressure with about 1 keV acceleration energy and annealing at around 700 • C. e sample was heated up and cooled down slowly, around 1 • C/s, to reduce strain in the crystal.
e tip was made from a polycrystalline vanadium wire of 99.8% purity, which was cut in air and prepared in ultrahigh vacuum by Argon spu ering.Subsequent eld emission on V(100) surface as well as standard tip shaping techniques were used to obtain a tip exhibiting clean bulk gap as well as good imaging capabilities.

YSR tip preparation
By controlled indentation, we can induce a YSR state to the apex of the vanadium tip, shown in detail in Ref. [1].We perform multiple cycles of dipping the tip into the sample to exchange material between tip and surface.In this way, YSR states with desired energy, intensity and asymmetry can be produced on the tip.e tip YSR state likely comes from a combination of simple elements such as oxygen and carbon that are picked up from the surface, but the exact mechanism requires further investigation experimentally and theoretically [1].
Reaching di erent temperature with the mK-STM e mK-STM is capable of reaching temperatures between 7 mK and 1 K.It has two continuous operation modes: one at 1 K where nearly all 3 He-4 He mixture is taken out from the dilution refrigerator, another at 10 mK where nearly all mixture is condensed in the mixing chamber of the dilution refrigerator.To reach around 0.6 K and 0.8 K, we can start from the 1 K mode and put a li le bit mixture inside the cycle to increase the cooling power.To reach a temperature from 10 mK to around 500 mK continuously, we can start from 10 mK and apply some resistive heating in the mixing chamber.
To reach below 10 mK, we can continue pumping the mixture but stop the mixture input into the mixing chamber to reduce heat load on the mixing chamber. is is a single shot mode because once the 3 He is pumped empty, the mixing chamber will warm up.Before that, the temperature can be lowered to between 6 mK and 7 mK.

EXPERIMENTAL DETAILS V(100) surface and intrinsic YSR defects
Vanadium is a type-II conventional Bardeen-Cooper-Schrie er (BCS) superconductor, with a transition temperature  C = 5.4 K and a gap parameter  760 μeV [2,3].Our preparation of tip and sample repeatably yields the measured gap parameter close to the bulk value for both tip and sample.In the temperature range between 7 mK and 1 K, both tip and sample are well superconducting.e V(100) surface typically shows a (5 × 1) oxygen reconstruction due to the di usion of bulk oxygen impurities to the surface through spu ering-annealing cycles during sample preparation [4][5][6].
e existence of such an oxygen layer does not have any e ect on the superconductivity measured on the surface [7].
ere are various types of impurities on the surface, with oxygen vacancies and carbon probably the most abundant [1], which show no magnetic signature for the most cases.
ere is a sparse distribution of intrinsic magnetic impurities (on the order of 0.04% of a monolayer) that show YSR states inside the superconducting gap on the surface, which is observed consistently throughout di erent samples or preparation cycles.e abundance of such impurities, however, does depend on the annealing temperature.e exact origin of such impurities requires further research, but they probably come from some uncommon arrangements of simple elements such as oxygen and carbon [1].e reverse triangles denote the expected energy position for the process depicted in (d), which is forbidden resulting in the absence of the peak.(f) e conductance dependence of the corresponding peaks from (e), showing di erent order   .e horizontal dahsed line is the noise limit.If there were peaks following the process in (d), the reverse triangles would follow a  2 dependence, which is not the observation here.

Spin-forbidden YSR-MARs
e spin polarization of YSR states has been shown previously by spin polarized STM [8], but here we use the signature in multiple Andreev Re ections (MARs) to demonstrate the full spin polarization of YSR states without needing a spin polarized probe (nor magnetic eld).MARs are higher order tunneling processes that arise when one or both electrodes are superconducting.MARs are responsible for many subgap features in the tunneling spectra between two superconductors.Depending on the number of quasiparticles being transferred across the tunnel junction, the intensity of the associated spectral feature scales with the conductance   according to a power law ∝    ,  = {2, 3, 4...}, where  is the order of MAR. e order of MAR is de ned here as the number of quasiparticles crossing the interface and, therefore, the lowest order MAR is second order (the simple Andreev re ection) with one re ection.erefore, from the peak position as well as the conductance dependence of a certain peak in the  / spectrum, we can unambiguously a ribute the relevant MAR and its order.
Due to the superlinear conductance dependence, MARs only become visible at relatively high conductance values.In the main part of this paper (Shiba-Shiba tunneling), we worked at very low conductance so that no MAR features are involved, simplifying the analysis to extract the relative spin angle.Here, however, we are going to show in the relatively high conductance regime that MARs can be utilized to demonstrate the full spin polarization of YSR states.e experimental data shown in Figs.S1(e,f) comes from the tunneling between an intrinsic sample YSR state on the V(100) surface and a clean superconducting vanadium tip, where the gap parameters are   = 700eV,   = 760eV for the tip and sample respectively and the YSR energy being  = 260eV extracted from the spectrum.Figure S1 (e) (i.e., main text Fig. 2(b), reproduced here for easier discussion and comparison) is measured at the normal state conductance of around 0.16 0 .All labels ( * , +, •, etc.) across Figs.
In the case of one YSR state on the sample, multiple MAR processes can occur (Fig. S1) which have been discussed in detail in theory [9].First, there are direct quasiparticle processes (Fig. S1(a)) responsible for the  / peaks at the gap edge (labeled with * ) as well as the conventional YSR-BCS tunneling peak (labeled with +).ese processes scale linearly with   at low conductance but enter a sublinear regime at high conductance (Fig. S1(f)), which has been discussed before in the context of multiple tunneling events and lifetime e ect [1,10].e second family is conventional MARs without involving the YSR states (Fig. S1(b)), and the peaks labeled with • in Fig. S1(e) is a ributed to be a lowest order MAR (second order) that scales quadratically with the conductance (Fig. S1(f)).
e third family comprises MARs connecting the YSR level with the clean superconductor coherence peaks (Fig. S1(c)), where the second order (labeled with ) as well as third order (labeled with ×) processes are observed in Figs.S1(e),(f) featuring the expected conductance dependence.ese MARs are usually quite strong because of the spectroscopic signi cance of YSR levels compare to the coherence peaks. is family of MARs does not feature any spin dependence, because it tunnels into the continuum, which is spin degenerate.
e fourth family is exotic because it starts from and ends in a YSR state (Fig. S1(d), labeled with ).Imagine a spin up electron tunnels to the right, and in order to generate a Cooper pair on the right hand side, it needs to take one spin down electron from the le to pair with. is process is not allowed because the YSR state is fully polarized and the electron part has only one spin species.Apart from this intuitive description, it has been shown in theory that this family of MARs is strictly forbidden due to the spin polarization of the YSR state [9].If the YSR state were a spin degenerate level or not fully spin polarized, we would expect to see the peaks originating from this family in the spectrum.ese features would be quite prominent due to YSR state being nearly a singularity in the density of states.As a side note, if the quasiparticle background of the superconductor is non-negligible, there could be a peak at the expected energy position, but will scale linearly with conductance, because it is direct quasiparticle tunneling from the YSR level to the continuous quasiparticle background of the other electrode rather than MARs [9,11].
In the experiment, we do observe the above spin ltering.We have seen all lowest order MARs (second order) as well as one third-order MAR in the other families, but there is no peak at the expected position for the lowest order process of the fourth MARs family (Figs.S1(d),(e) ).Plo ing the  / signal at the expected position against conductance in Fig. S1(f) ( ) does not show the expected ∝  2  dependence.Consequently, without applying a magnetic eld or using a spin-polarized probe, we demonstrate that the YSR states are fully spin polarized by showing that this "self Shiba-Shiba tunneling" (fourth family of MARs) is strictly forbidden.

Shiba-Shiba tunneling
Although Shiba-Shiba tunneling results in multiple peaks in the narrow superconducting gap even in the case of single YSR state on each electrode, the spectral features are generally well separated allowing for unambiguous a ribution.Since 0 <  , <  , , we have |  −   | < |  +   | < | , +  , |, where  , are the superconducting gap parameter of the sample and the tip, respectively. = | , +  , | is the lowest bias voltage possible to drive a direct tunneling process into the quasiparticle continuum outside the superconducting gap, which corresponds to the conventional tunneling from a YSR state into the continuum outside the gap of the other electrode.Consequently, the thermal Shiba-Shiba peak ( = ±|  −   |) is closest to zero voltage, separated from the direct Shiba-Shiba peak( = ±|  +   |), while all other features are more outside.
is becomes more complicated if the temperature is higher or the conductance is higher, where more peaks are expected inside the gap due to thermal excitations or higher order processes.
e relative spin orientation in relation to the quantum phase transition e spins of the YSR levels depend not only on the impurity spin, but also on which side of the quantum phase transition (QPT) the system is.When a YSR state moves across a QPT, e.g., when the spin-dependent sca ering strength crosses a threshold, the  − −ℎ + asymmetry will be exchanged, accompanied with reversed polarization of the levels above and below the Fermi energy [8,12,13].Consequently, the spin-parallel case in the discussion in the main text corresponds to two possibilities, either the impurity spins are parallel and the YSR states are on the same side of the QPT, or impurities have opposite spins and the YSR states are on the opposing side of the QPT.e spin-anti-parallel case would correspond to the opposing combinations, i.e., parallel impurity spins but at opposing sides of the QPT, or anti-parallel impurity spins at the same side of the QPT.  +   − should be done under appropriate tunneling conditions, because the theoretical model shown in the main text and in more detail later here in the supplementary material is valid only in the low conductance limit (linear regime [1]).At higher conductance, resonant processes as well as multiple Andreev re ections make the interpretation more complicated [9].erefore, in order to extract the Shiba-Shiba peak intensities correctly, the system has to be in the linear regime, where the Shiba-Shiba peak intensities depend linearly on the tunneling conductance and the YSR states have enough time to relax into the ground state before the next tunneling event.
As the intrinsic lifetime of the YSR states varies from system to system, it is not a priori obvious below which conductance threshold the system behaves linearly, so a conductance dependency measurement is necessary [1].To illustrate this, we have plo ed the Shiba-Shiba peak currents for a typical YSR tip and YSR sample system in Fig. S2(a) as a function of the normal state conductance.
e normal state conductance is extracted from the conductance measurement at 4 mV, which is much larger than the superconducting gap.
e measurement extends several orders of magnitude of the conductance to resolve the di erent regimes.e peak current values have been extracted as indicated in some selected spectra in Fig. S2(b).Only in the blue shaded area in Fig. S2(a), the system behaves linearly.Beyond this conductance threshold, higher-order e ects start to emerge forcing the system into a sublinear conductance dependence.As a result, to measure  correctly, one needs a conductance dependence experiment and take  in the linear regime (shaded area in Fig. S2(c)).
e second parameter  th can be extracted from the experiment in one of two ways.One way is to directly obtain the temperature  from the thermometer reading and the YSR energies  , from the position of the peaks in the spectra.e other method to determine  th experimentally is independent of the rst method.It relies on the ratio between the peak intensities of the direct and thermal transitions in the single impurity YSR tunneling into the continuum.Here, the ratio between the thermal peak (at  = ±( , −  , )) and the direct peak (at  = ±( , +  , )) in the di erential conductance is simply the Boltzmann factor  − , / B  (see theory part below and [10]).Fig. 2(c) of the main text shows the result of the ratio for various tip and sample YSR states separately at 1 K as a function of the respective YSR state energy using the la er procedure, showing good consistency.Consequently, we can characterize the Boltzmann factor of the sample YSR state by measuring the spectrum using a clean superconducting tip, and measure that of the tip YSR state by moving the YSR tip to a clean spot on the sample without YSR states and measure the spectrum.A er these measurements,  th in the Shiba-Shiba tunneling is then simply the di erence between the two measured Boltzmann factors as in Eq. (S62).
Error estimation of the thermal-direct Shiba-Shiba peak ratio e uncertainty   of the ratio  = √︃   +   −   +   − shown as shading in Fig. S2(c) comes from three factors.First, the measured Shiba-Shiba peak current   (especially the thermal Shiba-Shiba peak   ) in the linear regime can be as low as 10 fA so that the noise from the I-V converter becomes signi cant in the determination of the current peak height.Second, the linear regime may locate at quite small conductance making its determination di cult, which is especially challenging at very low temperature.ird, such a set of measurements takes time due to long sweep average as well as multiple spectra for conductance dependency, and relocating the YSR defect a er dri may contribute to additional error.
e rst error contribution can be quantitatively characterized, and the error estimate of each single point in Fig. S2(c) (yellow shading) considers this factor only, with the formula , where   is the uncertainty of the current measurement extracted from the spectrum.We assume that the current noise is independent of the bias voltage because at such low current, the noise from the I-V converter dominates which is independent of the bias.
e second and third error contribution can only be estimated roughly, and the total error estimate of  (gray shading in Fig. S2(c) and subsequently error bars in Fig. 2(c) and Figs.4(a),(b) of the main text) also try to include these e ects.Notice that the ratio is not prone to the third error because the thermal YSR and direct YSR state usually have similar spatial extension, but the absolute intensity depends critically on the spatial position, which may explain the deviation from constant of the temperature dependency of direct Shiba-Shiba peak area in Fig. 4(c) in the main text.
Shiba-Shiba tunneling below 0.6 K e ratio becomes di cult to access below 0.6 K because of the exponential suppression of thermal Shiba-Shiba peak.In addition, the YSR lifetime is enhanced when cooling down due to the suppression of thermal relaxation channels, resulting in the reduction of the linear regime to even lower conductance (cf.Fig. S2(a)), making both Shiba-Shiba processes harder to resolve.erefore, the direct Shiba-Shiba peak intensity is used to trace the relative spin angle  .Notice, however, unlike the ratio where both peak area and peak height can be used, only current peak area can be used here to eliminate the temperature e ect on the YSR lifetime and environmental broadening [1,[14][15][16][17]]. e reason is that environmental broadening as well as external noise preserves the peak area.

GREEN'S FUNCTION THEORY OF SHIBA-SHIBA TUNNELING
To simulate the tunneling spectra, we use the non-equilibrium Green's function theory formalism.First, we work in the 2 × 2 Nambu space to demonstrate some basic properties, then we introduce the 4 × 4 spin-Nambu space which is essential for arbitrary spin orientations.

General formalism in the 2×2 Nambu space
We describe a general superconductor-superconductor tunnel contact with the Hamiltonian [18] where Ĥ, are the unperturbed Hamiltonians of the two electrodes,  is the coupling parameter,  () is the time-dependent superconducting phase with  () =  0 + 2 0  ℏ ,  0 =  , and  is the bias voltage applied across the junction.e general expression of current between two electrodes can be wri en in the time domain in which  3 is the Pauli z matrix in Nambu space, Tr is the trace over Nambu degrees of freedom, and t () is the time-dependent Green's function of a YSR state e Hamiltonian of a YSR state in the classical Shiba model can be wri en as [12,13] where  is the exchange coupling,  is the spin and  is potential sca ering.e YSR Green's function in the 2 × 2 Nambu space can be derived from the above Hamiltonian with dimensionless parameters  =  0  and  =  0  , and  0 contains the normal-state density of state at the Fermi energy.
When  = 0,  = 0, it reduces to the usual BCS Green's function A substitution of  →  +  results in the retarded Green's function   and a substitution of  →  −  results in the advanced Green's function   , where  is a small real number.e unperturbed lesser and greater Green's functions can be obtained from Eqs. (S11) and (S12).e 2 × 2 Green's function above only gives one pole, whose energy is the Shiba energy e quantum phase transition (QPT, which is below the QPT (weak coupling), then  0 > 0. If  > √︁ 1 +  2 , which is above the QPT (strong coupling), then  0 < 0. Changing the sign of , which means ipping the spin of the impurity, also changes the sign of  0 .
We are interested in the Green's function near  =  0 , in which the Green's function of the YSR can be simpli ed to with Here, changing the sign of  exchanges  2 and  2 .Now we write the retarded and advanced Green's function in which the phenomenological parameter  stands for the intrinsic relaxation of the Shiba state.With Eqs.(S11) and (S12) in mind, we can write the lesser and greater Green's functions as well as the densities of states: ermal-direct ratio in YSR-BCS tunneling We can simluate the conventional YSR-BCS tunneling spectrum by plugging the BCS Green's function as one electrode (  ) and the YSR Green's function as another electrode (  ) into the expression for the tunneling current Eq.(S14).We know that e function in the integral is only non-zero around   0 , so we substitute  =  −  0 and limit the integration to  ∈ [−, ], where the parameter  is chosen such that    B  so that the Fermi function can be approximated as constant in the integration range.e existence of such ℎ is ensured by the assumption of   B   which is generally true for our experiment.Under this assumption,  () 0,  (−) 1.We assume also  0 > 0 for convenience.Given that  BCS () sharply peaks at  = ±, the rst term in Eq. (S25) corresponds to the peaks at  0 = − 0 ±  ( + and  − , where ,  stand for thermal and direct YSR-BCS tunneling, and +, − stand for the process at positive and negative voltage) and the second term corresponds to the peaks at  0 =  0 ∓  ( − and  + ).For  0 not too close to  nor to zero (not within ±, which also holds generally in the experiment) so that at each resonance only one of the two terms is dominant, the expressions at the resonances become where we utilize the even property of  BCS () such that It can be easily seen that the asymmetry of the thermal YSR-BCS peaks is reversed compared to the asymmetry of the direct peaks,  t+  t− =  d−  d+ = − 2  2 , which was also shown in Ref. [10].e thermal-direct ratio  YSR−BCS is which is the Boltzmann factor.Here we use the current peak to calculate the ratio, but the  / peak height can also be shown to yield the same ratio.Consequently, we can measure the thermal-direct ratio of a conventional YSR-BCS tunenling process ( / or current peak), and the ratio should be the Boltzmann factor of the YSR state.Experimentally this is done by moving the YSR tip to a clean position on the surface without sample YSR states and measure the spectrum at low conductance, and use a clean tip to measure the intrinsic YSR state on the surface.In this way, the Boltzmann factor of the tip and sample YSR state can be separately characterized experimentally.is can also be used as a temperature sensor because  0 can be easily read out in the spectrum, leaving  as the only parameter le .We actually used this property to calibrate the temperature sensor at this temperature range of the machine.Notice that below approximately 0.6 K, the thermal process is usually di cult to detect experimentally.
Shiba-Shiba tunneling in the 2 × 2 Nambu Space Here, we rst demonstrate the Shiba-Shiba tunneling in the 2 × 2 formalism established in the previous part, to understand the lineshape as well as some basic properties of Shiba-Shiba tunneling.is formalism is relatively simple, but it is limited in the way that the two YSR states can only have parallel or anti-parallel spin.
ere is only one pole in the YSR Green's function, which means that only one among the two processes (thermal or direct Shiba-Shiba) can exist in one spectrum.A 4 × 4 formalism is needed to overcome this limitation and re ect more general experimental situations which will be discussed later.
To derive a current expression in the 2 × 2 formalism, we need to insert two Shiba Green's functions as the le and right electrode.It is su cient to use the simpli ed formula in Eq. (S19), because the Shiba-Shiba peak is always separated from other spectral features and thus only the Green's function near the resonance participates: Plugging Eq. (S24) into Eq.(S14), we have which is a convolution of two Lorentzians and the Fermi function.is integral can be carried out analytically under the assumption which is generally true for Shiba-Shiba tunneling even at 10 mK [1].e rst term in the integral of Eq. (S33) is only nonzero near  =   ,  0 =   −   within a range of the order of  , .e second term is only nonzero near  =   ,  0 =   −   within a range on the order of  , .Due to the assumption in Eq. (S34), the Fermi function can be approximated as constant in this region and taken out of the integral.Using the de nite integral of two Lorentzian functions we arrive at the following expressions 1 0 and  2 0 have di erent signs of current and resonant voltage, meaning that they are responsible for the positive and negative Shiba-Shiba peaks.e current peak heights are e current peak areas are, however, independent of the intrinsic relaxation  , : In the temperature-dependent measurement of direct Shiba-Shiba peak, the peak area is used rather than the peak height, because in this way the temperature dependence of the intrinsic relaxation can be ruled out, leaving only the possible variation of spin dynamics.

Properties of Shiba-Shiba tunneling
Continuing from the discussion before about the sign of  and the relation between  and √︁ 1 +  2 , the situation of two YSR states with parallel spins (same side of QPT and parallel impurity spin or di erent side of QPT and opposite impurity spin) is equivalent to     > 0, and the situation of two YSR states with opposite spins (di erent side of QPT and parallel impurity spin or same side of QPT and opposite impurity spin) is equivalent to     < 0.
If     > 0,  0 is exponentially suppressed at low temperature by the prefactor  (  ) −  (  ), which refers to a thermal Shiba-Shiba process.
e bias voltage needed is e direct Shiba-Shiba process is missing in this situation, corresponding to  = 0 • in the main text.
Note that due to the limitation of the 2 × 2 formalism, the spins are either parallel or anti-parallel, so that the thermal Shiba-Shiba and direct Shiba-Shiba processes cannot coexist.To account for arbitrary angles, we need to extend to the 4 × 4 formalism as discussed in the following.

Extension to the 4×4 Nambu space
To account for arbitrary angles between tip and sample YSR spins, we extend the Green's function theory to the 4 × 4 spin-Nambu space.With the convention used in Ref. [12], the basis for the 4 × 4 spin-Nambu space is e general expression of the current can be extended from Eq. (S2) to the 4 × 4 spin-Nambu space [9]: is model consists of a chain of  fermionic sites, each characterized by energy   and corresponding creation (annihilation) operators given by  †  (  ), for  = {1, . . .,  }. e chain has two phases, a trivial phase where the spectrum has a gap and no edge modes are present, and a topological phase where the spectrum contains a single fermionic state at zero energy formed by combining two Majorana modes localized at the two edges of the chain.
e topological phase of the chain is sensitive to the two amplitudes   and   describing coupling between nearest neighbor sites.e amplitude   describes the pair creation process, while   is the amplitude of the tunnel hopping process.When both amplitudes exceed the energies of the sites,   ,     for all sites, the chain is in the topological phase.In contrast, when either   <   , or   <   , the topological phase is lost and the chain does not host Majorana modes at its edges.
e above Hamiltonian can be realized by a chain of tunnel coupled spin-polarized YSR states on a common superconducting substrate [20][21][22][23][24][25][26][27][28], each with energy   .If the dynamics of the magnetic impurity spins are decoupled from the electronic degrees of freedom, the above Hamiltonian describes the low energy excitations.In the analogy between the tunnel junction and the YSR chain,   is the amplitude of the direct Shiba-Shiba tunneling process and   the amplitude of the thermal Shiba-Shiba tunneling process.e YSR chain will reach the topological phase as long as both amplitudes exceed a certain threshold set by the disorder of the individual Shiba energies.It is, therefore, important to show that the relative spin-polarization of neighboring YSR states allows the two processes to coexist.
e mechanism leading to this coexistence in the case of the experiment Ref. 29 remains under theoretical debate.A number of works have proposed that the RKKY interaction in such a chain leads to a helical ordering [21,22,25,27,[30][31][32], while a counter view holds that anisotropy due to the superconducting surface may lead to ferromagnetic ordering, where the spin blockade of the pair creation process is li ed only if the chain exhibits strong spin-orbit coupling [23,26,28].To shed light on these discussions, separately detecting the two processes between neighboring sites experimentally and clarifying the role of spin there is desired.
Our results demonstrate the coexistence of the two processes for a pair of YSR impurities, one placed on the substrate and one on the tip of the STM, which paves the road for further research regarding the conditions the relative spin orientation in a chain of atomic impurities, such as realized in Ref. 29,33, exhibits the same lack of anisotropy.Our results allow us to conclude that a functionalized STM tip with an YSR state will couple to an impurity chain in such a way that the topological phase of the chain can be extended to the YSR state on the tip.

FIG. 1 :
FIG. 1: Direct Shiba-Shiba tunneling and thermal Shiba-Shiba tunneling.(a) Schematics of the experimental setup realizing tunneling between Yu-Shiba-Rusinov (YSR) states using a scanning tunneling microscope, with  denoting the angle between the two spins.b,c) Tunneling processes of direct (b) and thermal (c) Shiba-Shiba tunneling.In the direct process, the  − part and the ℎ + part of the two YSR states are aligned by the bias voltage.In the thermal process, the two  − parts or two ℎ + parts are aligned. , are the superconducting gap parameters of the sample and tip respectively, with both around 760 eV for vanadium.(d),(e) Typical  ( ) spectra of Shiba-Shiba tunneling at 10 mK (d) and 1 K (e).Blue arrows ( ± ) indicate direct Shiba-Shiba peaks, and orange arrows ( ± ) indicate thermal Shiba-Shiba peaks.

FIG. 2 :
FIG.2: Properties of the YSR states.(a) A family of multiple Andreev Re ection processes (MARs) connecting the levels of the same YSR state (the schematic here show only the lowest order process).Due to full spin polarization of the YSR state, this family of MARs is spin forbidden.(b) A typical dI/dV spectrum involving one YSR state tunneling into a clean superconductor at high conductance (here ∼ 0.16 0 , where  0 = 2 2 /ℎ is the conductance quantum) to reveal MARs.Detailed peak assignment can be found in the supplementary information[25].eexpected position for the lowest order MARs in (a) is labeled with symbol , where no peak is seen conrming the spin blockade.(c) Statistics of the thermal-direct ratio of the conventional YSR-BCS peak (in the low conductance regime) of various YSR impurities (tip or sample) with di erent YSR energy measured at 1K.All data points fall on top of the prediction of the Boltzmann factor at 1K within a narrow window (±50 mK).

FIG. 4 :
FIG. 4: Temperature dependency measurements.(a) Statistics of the thermal-direct Shiba-Shiba ratio  (of peak current) measured at 1K, plo ed against the di erence of the Boltzmann factor  th of tip and sample Shiba state.All data points lie on the identity line within a small error bar, indicating a spin averaging e ect.(b) e temperature dependence of thermal-direct Shiba-Shiba ratio  of selected examples with di erent YSR energies (data points) agreeing with the prediction from the di erence of the Boltzmann factor  th (curves of corresponding color).(c) e temperature dependence of the direct Shiba-Shiba current peak area normalized by the transmission  =   / 0 , which remains nearly constant, with no indication of spin freezing down to 7 mK.(d) e direct and thermal tunneling between two YSR states constitute the fundamental interaction mechanism of the formation of a chain of magnetic impurities on a superconductor.e co-existence of both processes is crucial for the chain to host topologically non-trivial end states.

FIG. S1 :
FIG. S1: Multiple Andreev Re ections (MARs) with YSR states.(a) Direct quasiparticle tunneling processes (no Andreev re ection).Le ( * ): tunneling between continua.Right (+): tunneling from YSR state to the continuum.(b) MARs without involving the YSR state.Here in the diagram (•), the second order MAR is shown (one re ection).(c) MARs connecting YSR levels to coherence peaks.Le ( ): second order process (one Andreev re ection).Right (×): third order process (two Andreev re ections).(d) MARs connecting levels of the same YSR state.Here in the diagram ( ): second order process (one reection).is family of MARs is strictly forbidden due to the full spin polarization of the YSR state.(e) (main text Fig. 2(b))  / spectrum of a YSR state on the sample surface at the normal state conductance of 0.16 0 , with all peaks labeled with the same marks in (a)-(d) indicating the origin of each peak.e reverse triangles denote the expected energy position for the process depicted in (d), which is forbidden resulting in the absence of the peak.(f) e conductance dependence of the corresponding peaks from (e), showing di erent order   .e horizontal dahsed line is the noise limit.If there were peaks following the process in (d), the reverse triangles would follow a  2 dependence, which is not the observation here.

FIG. S2:
FIG. S2: Conductance dependence measurements of the Shiba-Shiba tunneling.(a,b) Conductance dependency of thermal and direct Shiba-Shiba tunneling.(a) e evolution of the Shiba-Shiba current peak height with respect to normal state conductance, showing the linear regime in the low conductance limit as blue shaded area.(b) Selected  ( ) spectra at di erent conductance, with thermal and direct Shiba-Shiba peaks labeled respectively.(c) Conductance dependency of the thermal-direct Shiba-Shiba ratio, with the low conductance value indicated by the gray dashed line being the value  used in this paper.e vertical shadings are the estimated uncertainty.