Investigation of the Yu-Shiba-Rusinov states of a multi-impurity Kondo system
A. Kamlapure, L. Cornils, R. Žitko, M. Valentyuk, R. Mozara, S. Pradhan, J. Fransson, A. I. Lichtenstein, J. Wiebe, R. Wiesendanger
** [email protected] † [email protected] Investigation of the Yu-Shiba-Rusinov states of a multi-impurity Kondo system
A. Kamlapure *, L. Cornils , R. Žitko , M. Valentyuk , R. Mozara , S. Pradhan , J. Fransson , A. I. Lichtenstein , J. Wiebe , and R. Wiesendanger Department of Physics, University of Hamburg, Jungiusstrasse 11, D-20355 Hamburg, Germany Jožef Stefan Institute, Jamova 39, SI-1000 Ljubljana, Slovenia Faculty of Mathematics and Physics, University of Ljubljana, Jadranska 19, SI-1000 Ljubljana, Slovenia Department of Theoretical Physics and Applied Mathematics, Ural Federal University, 19 Mira Street, Yekaterinburg, 620002, Russia Department of Physics and Astronomy, Uppsala University, P.O. Box 516, Uppsala SE-751 21, Sweden
Recent studies of mutually interacting magnetic atoms coupled to a superconductor have gained enormous interest due to the potential realization of topological superconductivity . The Kondo exchange coupling (cid:1836) (cid:3012) of such atoms with the electrons in the superconductor has a pair-breaking effect which produces so-called Yu-Shiba-Rusinov (YSR) states within the superconducting energy gap, whose energetic positions are intimately connected with the requirements for topological superconductivity. Here, using the tip of a scanning tunneling microscope, we artificially craft a multi-impurity Kondo system coupled to a superconducting host consisting of an Fe adatom interacting with an assembly of interstitial Fe atoms on an oxygen-reconstructed Ta(100) surface and we experimentally investigate the signatures of Kondo screening and the YSR states. With the help of numerical renormalization group (NRG) calculations, we show that the observed behavior can be qualitatively reproduced by a two-impurity Kondo system whose inter-impurity antiferromagnetic interaction (cid:1836) is adjusted by the number of interstitial Fe atoms in the assembly. When driving the system from the regime of two decoupled Kondo singlets (small (cid:1836) ) to that of an antiferromagnetic dimer (large (cid:1836) ), the YSR state shows a characteristic cross-over in its energetic position and particle-hole asymmetry.
Introduction
Mutually interacting magnetic impurities which are Kondo exchange coupled to the electrons forming the Cooper pairs in a superconducting host constitute building blocks for the realization of a topological superconductor that can potentially host Majorana states on its edges . While the Kondo exchange coupling (cid:1836) (cid:3012) affects the position of the Yu-Shiba-Rusinov (YSR) states inside the gap of the superconductor (2∆) , the mutual magnetic interactions (cid:1836) between the impurities dictate how strongly these YSR states split or hybridize into bands , and govern the collective magnetic state of the magnetic impurity cluster. The investigation of the effects of varying (cid:1836) (cid:3012) and (cid:1836) on the YSR states in such systems is, thus, of fundamental interest towards the end of realizing topological superconductivity. If (cid:1836) (cid:3012) is sufficiently large, the impurity enters the strong-coupling Kondo regime where the impurity spin forms a many-body singlet ground state with the conduction electrons for temperatures below the characteristic Kondo temperature , (cid:1846) (cid:3012) . For (cid:1863) (cid:3003) (cid:1846) (cid:3012) ≳ ∆ , the spectroscopic signature of the Kondo screened spin, the so-called Kondo resonance centered at the Fermi energy and having a width of order (cid:1863) (cid:3003) (cid:1846) (cid:3012) , can coexist with the YSR states . For normal metallic substrates, the behavior of two Kondo impurities interacting via the same conduction electron system that produces the Kondo screening, the so-called two-impurity Kondo problem , has been studied in great detail . For antiferromagnetic (AFM) inter-impurity interaction (cid:1836) of increasing strength, the system passes over from two decoupled Kondo singlets (small (cid:1836) ) to an AFM dimer (large (cid:1836) ), which causes the Kondo resonance to split. However, the analog of this two-impurity Kondo problem for superconducting host systems has not been studied experimentally to a large extent so far , and it is an interesting question, how the YSR states behave as a function of (cid:1836) . In this work, we experimentally study artificial Fe clusters consisting of a single Fe adatom magnetically interacting with an assembly of interstitial Fe atoms (IFA) of different shapes and sizes on an oxygen reconstructed superconducting Ta(100) surface. Through comparison with NRG calculations, our spectroscopic investigation shows that the system qualitatively behaves as a two-impurity spin-1/2 Kondo system coupled to a superconductor with (cid:1863) (cid:3003) (cid:1846) (cid:3012) ~(cid:1836) ≫ ∆ , where the first Kondo impurity is associated to the Fe adatom and the second to the IFA assembly. As a function of the size of the IFA assembly, the system shows a smooth transition between the weak and strong exchange coupling regimes (i.e., independent Kondo singlets vs. magnetically aligned spins) observed through the formation of a spectral depletion in the Kondo peak, i.e., an exchange-splitting gap. Such simple behavior is surprising, since the impurities have multiple orbitals, potentially producing high-spin states. Furthermore, it is not a-priori expected that multiple IFAs can be conflated into a single effective spin degree of freedom. Interestingly, simultaneously with the transition observed in the shape of the Kondo peak on the large energy scale of (cid:1863) (cid:3003) (cid:1846) (cid:3012) ~(cid:1836) , the YSR state energy whose scale is small compared to ∆, as well as the particle-hole asymmetry of this state, reveal a characteristic crossover that can likewise be assigned to the variation of (cid:1836) with the help of the NRG calculations. Results
Construction of the Fe clusters
The details of the preparation of a (3 x 3) oxygen-reconstructed Ta(100) surface, in the following called TaO, and of the cold-deposition of Fe can be found in Refs. [4,51,52]. We use vertical atom manipulation to build various clusters consisting of a single Fe adatom and an assembly with different numbers of IFAs in close proximity (see Methods). Figure 1a shows a representative surface of TaO where an Fe adatom was placed at the center of the (3 x 3) plaquette and a single interstitial Fe atom was positioned between two (3 x 3) plaquettes . Figure 1c shows the result of a geometry optimization of the surface through DFT calculations revealing the stable positions of the Fe adatom and the interstitial Fe atom. Within the DFT calculations, we find a finite magnetic moment for both the Fe adatom as well as the IFA (see methods). The schematic top views and STM images of the other investigated clusters with 0 – 4 IFAs are shown in Figs. 1d,e (see Supplementary Fig. la for other examples of similar clusters). To build such clusters, initially, we moved Fe atoms to the interstitial locations to form the assembly, followed by placing the Fe adatom in close proximity (case 1 – 4). We also studied the single Fe adatom without an assembly of IFAs for comparison (case 0). Spectroscopy of the Fe clusters
Spectra measured on these clusters with a superconducting tip in a wide bias range (Fig. 2a) show a resonance around the Fermi energy. Here we have used a large amplitude of the modulation voltage such that the superconducting gap and features within it are washed out. For the bare adatom, the resonance is very broad (half-width at half maximum (HWHM) (cid:1985)~ (cid:1985)~12 (±1) mV]. We therefore assume that the resonance is of many-body origin (i.e., due to the Kondo effect), and interpret the differences of the spectral width of the resonance as an effect of the neighboring atoms on the Fe adatom’s Kondo coupling (cid:1836) (cid:3012) to the substrate electrons. Moreover, it is remarkable to see that, when the adatom is coupled to an assembly with two or more IFAs, the resonance peak appears to split, resulting in a zero bias anomaly (ZBA). A careful look at the spectra for adatoms showing ZBA reveals that, during the formation of the gap-like structure at the Fermi level, the overall width of the Kondo resonance remains unchanged. More precisely, the asymptotic tails of the peaks in the sequence of 1,2,3,4 IFAs overlap to a good approximation, indicating that J K remains constant, while an increasingly wide gap is forming at the peak center. In order to quantify the width of the zero bias resonance (cid:1985) and the gap width (cid:1831) (cid:3046) , we fit each dI/dV spectrum in Fig. 2a to a Fano curve with an additional multiplicative factor featuring a Lorentzian dip for the spectra with a ZBA (see Methods and dashed lines in Fig. 2a). It is evident from the fits that the details of the spectral shape are excellently captured using this functional form. Similar fits performed for a number of other clusters are shown in Supplementary Fig. 1b. The Kondo temperatures (cid:1846) (cid:3012) extracted from the fitted values of (cid:1985) (see Methods) are on the order of 25 K - 40 K. It is interesting to note that the Kondo energy scale in our system is very large compared to the superconducting gap ( (cid:1863) (cid:3003) (cid:1846) (cid:3012) ≫ ∆ ), and that the width of the ZBA, (cid:1831) (cid:3020) , is of similar size as (cid:1863) (cid:3003) (cid:1846) (cid:3012) . We will show below that the ZBA can be interpreted as an exchange-splitting gap due to a strong interaction between the Fe adatom and the assembly of IFAs. Next, dI/dV spectra taken on the same clusters in the low-bias range with small ac modulation will be discussed (Fig. 2b). They reveal sharp peaks inside the energy gap of the superconducting substrate due to YSR states . Since the measurement is performed using a superconducting tip, we extract the corresponding local electron density of states (LDOS) using a standard deconvolution method . The corresponding curves (Fig. 2c) show two sharp peaks, one with larger intensity located at E YSR and one with smaller intensity located at - E YSR . We observe a large variation of E YSR as the Fe cluster size and its configuration changes. In particular, E YSR tends to decrease and reverse its sign when the number of Fe atoms in the assembly is increased (see, in particular, configurations 1, 2, and 3). This is also apparent from the change in the sign of the asymmetry in the intensities of the two peaks (see Supplementary Fig. 2 for the definition of the asymmetry and data of all clusters). In order to establish a possible correlation between the YSR states and the shape of the Kondo resonance, we plot E YSR normalized to ∆ as a function of the reduced Kondo temperature, (cid:1863) (cid:2886) (cid:1846) (cid:2895) ∆⁄ (Fig. 3a) and, as a function of E S (Fig. 3b) extracted from the fitting. Correlation between Kondo resonance width and YSR energy
From Fig. 3a, it is evident that the Kondo temperature mostly changes by a factor of 1/2 (from ((cid:1863) (cid:3003) (cid:1846) (cid:3012) ∆⁄ ) ~ 10 to 5) when a single IFA is added to the Fe adatom, but does not further reduce with the number or configuration of IFAs in the assembly. At the same time, E YSR strongly changes from ~∆ to −∆ when the number of IFAs in the assembly is increased. Consequently, there is almost no correlation between E YSR and (cid:1846) (cid:3012) . At first glance, this result seems surprising, as similar trends of experimentally observed YSR energies were usually assigned to the quantum phase transition of the system induced by an increasing Kondo coupling . In order to exclude this effect as the major reason for the observed changes to E YSR , we simulated spectra within Slave-boson mean-field theory (SBMFT), which is ideally suited for the strong-coupling Kondo regime [ ((cid:1863) (cid:3003) (cid:1846) (cid:3012) ∆⁄ ) ≫ 1 ] and thus valid for the experimentally investigated system . We consider a simple model consisting of a single spin-1/2 Anderson impurity coupled to a superconducting substrate (see Methods). The SBMFT simulation of E YSR as a function of the Kondo temperature is plotted as a violet line in Fig. 3a. The deviation of the experimental data points from the simulation infers that the trend in E YSR in our case is not governed by changes to the Kondo coupling alone. This is also evident from the comparison of our data to experimental data and NRG calculations for another system (Supplementary Fig. 3), where we see that the system of Fe clusters on TaO investigated here always stays on one side of the Kondo-coupling induced quantum phase transition. Only the initial decrease of (cid:1846) (cid:3012) by a factor of 2 when a single IFA is added to the Fe adatom is the result of a change of its Kondo coupling, presumably through a small structural relaxation. On the other hand, from Fig. 3b, we observe that E YSR is strongly correlated with the ZBA in the large-bias spectra. In particular, while the YSR state energy is positive for most of the Fe clusters with a small number of IFAs in the assembly which do not show a ZBA in the Kondo resonance ( E S = 0), it is negative for those with a large number of IFAs which have a non-zero E S . The ZBA of the Kondo resonance can be naturally interpreted as an exchange splitting induced by the competition of Kondo screening ( (cid:1846) (cid:3012) ) and exchange coupling of the Fe adatom with the assembly of IFAs ( (cid:1836) ) which probably changes by the number of IFAs in the assembly . In order to qualitatively test this hypothesis, we consider the following minimal model. We regard our system as a realization of the two-impurity Kondo problem where the assembly of the IFAs is represented by a single effective spin-1/2 impurity interacting via exchange interaction (cid:1836) with the adatom as a second spin-1/2 impurity. Within this model, for small AFM interaction, the two spins are independently screened by the conduction electrons resulting in a single peak in the spectral function. For very strong interaction ( (cid:1836) ≥ (cid:1836) (cid:3030) ), the two spins undergo a crossover to a local AFM singlet state and are no longer screened independently . The spectral manifestation of this phenomenon is the observed splitting of the Kondo peak which for small and moderately large values of (cid:1836) takes the appearance of a “spin-gap” of width E S inside the resonance. The reason behind this particular line-shape is that the Kondo screening proceeds unperturbed at high energy scales, until it reaches the scale of (cid:1836) where the residual magnetic moments rapidly bind antiferromagnetically. This matches the evolution of the Kondo peak shape observed in the experiment. The width of the gap E S seen in the spectra with ZBA is, therefore, a direct measure of the effective strength (cid:1836) of the exchange interaction between the adatom and the assembly of the IFAs taken as a whole. Moreover, the order of magnitude of E S of 2 meV - 8 meV, implies a large AFM interaction between the adatom and the assembly of IFAs. The hierarchy of parameters is thus (cid:1836)~ (cid:1846) (cid:3012) ≫ ∆ . NRG calculations
In order to check whether the competition of Kondo screening ( T K ) and exchange coupling (cid:1836) within this minimal model can also explain the experimentally observed correlation between the width of the exchange-splitting gap E S and the YSR state binding energy E YSR (Fig. 3b), we have performed NRG calculations that produce reliable spectral functions on all energy scales, including the region inside the superconducting gap . We have used a two-impurity Anderson model to describe the system (see Methods section), with one impurity orbital representing the IFA assembly, and one representing the Fe adatom chosen to be particle-hole asymmetric. This model can also account for the effects of potential scattering on the adatom. The results show that while the ground state of the system is a spin singlet for all values of exchange coupling ( (cid:1836) ), its nature smoothly changes from a Kondo singlet to an AFM singlet as indicated e.g. by the increasing magnitude of the spin-spin correlation function. This is reflected by the emergence of the Kondo peak splitting for (cid:1836)~(cid:1836) (cid:3030) = 5∆ ( (cid:1836) (cid:3030) being the critical exchange interaction) (Fig. 4a) and the characteristic evolution of the YSR peaks (Fig. 4b). To understand the latter, we remind that for two decoupled impurities, the sub-gap spectrum consists of a singlet (Kondo) ground state S , two doublet excitations D (unscreened impurity states of either impurity), as well as a degenerate singlet S and triplet T , which correspond to the excited states where both impurities are unscreened. With increasing (cid:1836) , the singlet S decreases in energy, until it mixes with the Kondo state S leading to an anti-crossing shape in the energy diagram. For large (cid:1836) , the Kondo state S is pushed to higher energies together with the two doublet excitations D , while S is the new ground state. This leads to the observed evolution of YSR peaks in Fig. 4b, which for small (cid:1836) slightly decrease in energy, then for large (cid:1836) rapidly increase in energy and eventually merge with the continuum (Fig. 4c). The change of direction of YSR shifting occurs at (cid:1836) ~ (cid:1836) (cid:3030) which is very close to (cid:1836) where the exchange splitting becomes observable in the Kondo peak (compare Fig. 4c with 4d). This behavior is in qualitative agreement with the trend observed in the experiment, as shown by the comparison of the plot of E YSR versus E S extracted from the NRG-calculated spectra with the same parameters extracted from the experimental data (see Fig. 3b). We also observe that the sign change in the asymmetry of the YSR peak intensities occurs across (cid:1836) ~ (cid:1836) (cid:3030) , in line with the experiment (see Supplementary Fig. 2). Conclusions
In conclusion, the experimentally observed correlation between the small energy scale position of the YSR state and the large energy scale exchange-splitting gap of the Fe adatom coupled to the assembly of IFAs can be assigned to changes in the effective exchange coupling (cid:1836) between the adatom and the assembly in the regime where (cid:1836) ≶ (cid:1846) (cid:3012) ≫ ∆ . We have shown that the Fe adatom coupled to a single IFA has a reduced Kondo coupling with respect to the isolated Fe adatom, while further addition of IFAs to the assembly only increases (cid:1836) . We compared our experimental results to NRG calculations employing a minimal model that shows how with increasing exchange interaction the system of Fe clusters undergoes a crossover from the regime of an independent Kondo cloud to a local magnetically aligned state. This crossover results in the opening of a gap in the Kondo peak and in the anti-crossing of the sub-gap YSR states. Our results also indicate that introducing the potential scattering results in the change of the particle-hole asymmetry of the YSR peaks across this crossover, thereby facilitating its observation. Our work establishes the important role of exchange interaction on the experimentally observed shifts in YSR peaks, which were previously supposed to be prevalently determined by the Kondo coupling of the impurity with the substrate. Our results, therefore, challenge the current understanding of the intimate correlation between Kondo effect and YSR state and motivate further studies. In particular, it is becoming clear that while a phenomenological description in terms of hybridizing classical YSR states is possible and valid, the calculation of its parameters should be based on a microscopic quantum model that needs to properly incorporate exchange coupling (including its transverse parts that correspond to spin-flip terms). This calls for further development of efficient theoretical tools for solving sizable assemblies of quantum spins in interaction with a superconducting environment. 0
Methods
Experimental methods:
All the experiments were carried out in a custom-built commercial cryostat (SPECS) at a temperature of T = 1.2 K in ultra-high vacuum. Details of the sample cleaning and cold Fe deposition can be found in Refs. . To make different clusters of Fe atoms, we use vertical atom manipulation where Fe atoms were first transferred from the sample surface to the tip and then dropped back to the desired location. Details of the vertical atom manipulation and also the interstitial atom manipulation are described in Ref. . For spectroscopy, we used a superconducting tip which was prepared by coating a Cr tip with Ta via a controlled dipping of the tip into the Ta substrate. The d I /d V spectra were measured with the standard lock-in technique with a modulation frequency of f = 827 Hz and a typical modulation voltage of V mod = 20 μV (2 - 3 mV) for low (high) bias range spectra. To extract the local density of states from the measured spectra in the low-bias range, we use a standard deconvolution technique . Fitting of the Kondo resonance
In order to determine the Kondo temperature, we fit d I /d V spectra taken on the adatom to a Fano function given by : (cid:2026)((cid:1848)) ∝ (cid:2026) (cid:2868) + (cid:2026) (cid:2869) (cid:4672)(cid:1869) + (cid:1848) − (cid:1848) (cid:2868) (cid:1985) (cid:4673) (cid:2870) (cid:2868) (cid:1985) (cid:4673) (cid:2870) where (cid:1869) is the form factor which determines the line shape of the spectra, (cid:1985) is the half width at half maximum (HWHM), and (cid:1848) (cid:2868) is the bias value at the maximum of the spectra. (cid:2026) (cid:2868) and (cid:2026) (cid:2869) are the offset and the amplitude of the spectra, respectively. For those spectra showing the resonance and a ZBA, we fit the dI/dV curves to the product of a Fano function and an inverted Lorentzian, given by: 1 (cid:2026)((cid:1848)) ∝ (cid:2026) (cid:2868) + (cid:2026) (cid:2869) (cid:3430)(cid:4672)(cid:1869) + (cid:1848) − (cid:1848) (cid:2868) (cid:1985) (cid:4673) (cid:2870) (cid:2868) (cid:1985) (cid:4673) (cid:2870) (cid:3434) × ⎣⎢⎢⎡1 − (cid:2026) (cid:2870) (cid:3046) (cid:1985) (cid:3046) (cid:4673) (cid:2870) ⎦⎥⎥⎤ Here, (cid:1831) (cid:3046) = (cid:1857)(cid:1985) (cid:3046) is the half width in energy of the ZBA gap, (cid:1848) (cid:3046) is the position of the gap, and (cid:2026) (cid:2870) is a constant. Finally, we use the HWHM, (cid:1985) , to extract the Kondo temperature of each assembly and use Wilson’s definition to define the Kondo temperature as (cid:1846) (cid:3012) = 0.27(cid:1857)(cid:1985) (cid:1863) (cid:3003) ⁄ , where k B is the Boltzmann constant. Theoretical methods: DFT method:
To determine the low-temperature composition of the Fe cluster with IFA, we performed a set of DFT simulations using the VASP package . A minimal supercell of the bcc(001) crystal was built from five layers of Ta atoms with a lattice constant of 3.30 Å. The supercell of the surface was chosen to reproduce the area of one (3 x 3) plaquette, as seen in STM images. All simulations were done using spin-polarized functionals using the generalized gradient approximation (GGA+U) for the exchange-correlation part . The standard calculation was performed with U = 4.3 eV, J = 0.9 eV for the Fe d -states and U = 6 eV, J = 0.8 eV for the O p -states. An enlarged cut-off energy was set up to 500 eV and the number of bands was increased by 50%. The structural optimization was performed until forces were less than 0.01 eV/Å on the k -points grid. We searched for the possible positions of IFAs in the first interspace of the TaO substrate, i.e. between the first and the second layer of the bcc(001)-arranged Ta atoms. These configurations of the surface with a single IFA were then compared by total energies and the induced type of reconstructions (see Supplementary Fig. 4). 2 The adsorption of the Fe adatom at the center of the (3 x 3) plaquette was studied separately to account for van-der-Waals (vdW) contributions at the surface. The slab was transformed to the symmetric setup of five layers with two identical surfaces to avoid an interaction with charge images. The midpoint of the (3 x 3) plaquette coincides with the central Ta atom, which is marked as a particularly polarized atom in the surface electronic structure . To account for the polarizability of the local structures on the surface, we performed a set of structural optimizations including dispersion correction energies, with an accuracy of atomic forces of less than 0.001 eV/Å. First, the slab with an adatom in the central position was relaxed using the Tkatchenko-Scheffler (TS) method , where all surface atoms, including the adatom, were free to relax. To account for possible corrections from the ionic states of the surface, we also included the Hirschfield iterative (HI) partitioning scheme (TS/HI) with a self-consistent cycle . After relaxation, the Fe adatom has a tiny in-plane offset from the symmetric central position, with a magnitude of the vector of only 0.001 Å. The extracted position of the adatom was then used to compare between different methods in simulations with both the Fe adatom and the IFA. Numerical renormalization group method:
We have used the “NRG Ljubljana” implementation of the numerical renormalization group. The calculations have been performed for a two-impurity Anderson model given by: (cid:1834) (cid:3080) = (cid:3533) (cid:2035) (cid:3038) (cid:1855) (cid:3038)(cid:3097)(cid:3080)(cid:2993) (cid:1855) (cid:3038)(cid:3097)(cid:3080)(cid:3038)(cid:3097) − ∆ (cid:3533)(cid:3435)(cid:1855) (cid:3038)↑(cid:3080)(cid:2993) (cid:1855) (cid:2879)(cid:3038)↓(cid:3080)(cid:2993) + (cid:1834). (cid:1855). (cid:3439) (cid:3038) + (cid:3533) (cid:2035) (cid:3080) (cid:1856) (cid:3097)(cid:3080)(cid:2993) (cid:1856) (cid:3097)(cid:3080)(cid:3097) + (cid:1847)(cid:1866) ↑(cid:3080) (cid:1866) ↓(cid:3080) + (cid:3533) (cid:1848) (cid:3038)(cid:3080) (cid:3435)(cid:1856) (cid:3080)(cid:2993) (cid:1855) (cid:3038)(cid:3080) + (cid:1834). (cid:1855). (cid:3439) (cid:3038)(cid:3097) (cid:1834) (cid:2869)(cid:2870) = (cid:1836)(cid:1871) (cid:2869) ⋅ (cid:1871) (cid:2870) − (cid:1872) (cid:3533)(cid:3435)(cid:1856) (cid:2869)(cid:3097)(cid:2993) (cid:1856) (cid:2870)(cid:3097) + (cid:1834). (cid:1855). (cid:3439) (cid:3097) (cid:1834) = (cid:1834) (cid:2869) + (cid:1834) (cid:2870) + (cid:1834) (cid:2869)(cid:2870) . 3 Here α indexes the two impurities ( α = 1, 2), σ is the spin (↑, ↓), k is the electron momentum, (cid:2035) (cid:3038) the dispersion of band electrons, and (cid:1848) (cid:3038)(cid:3080) is the hybridization between the impurity α and the superconducting band ( Γ = (cid:2024)(cid:2025)|(cid:1848)| (cid:2870) , where (cid:2025) is the density of states in the band in the absence of superconductivity). Finally, (cid:1866) (cid:3097)(cid:3080) = (cid:1856) (cid:3097)(cid:3080)(cid:2993) (cid:1856) (cid:3097)(cid:3080) is the occupancy operator for impurity α , while (cid:1871) (cid:3080) is the impurity spin operator defined as (cid:1871) (cid:3080) = (1 2⁄ )(cid:1856) (cid:3036)(cid:3080)(cid:2993) (cid:2252) (cid:3036)(cid:3037) (cid:1856) (cid:3037)(cid:3080) , where (cid:2252) is the vector of Pauli matrices, while i,j are the internal spin indexes. We use the parameters (cid:1847) (cid:2869) = (cid:1847) (cid:2870) = Γ (cid:2869) = Γ (cid:2870) = (cid:2035) (cid:2869) = −4 , (cid:2035) (cid:2869) = −5 , with the gap fixed at Δ = (cid:1872) = J . The calculations were performed for the discretization parameter Λ = (cid:1840) (cid:3053) = Slave-boson mean-field theory calculations:
The total Hamiltonian for an Anderson impurity coupled to a superconducting substrate is (cid:1834) =(cid:1834) (cid:3020) + (cid:1834) (cid:3005) + (cid:1834) (cid:3021) . Here the substrate is described by the Hamiltonian (cid:1834) (cid:3020) = ∑ (cid:2035) (cid:3038)(cid:3097) (cid:1855) (cid:3038)(cid:3097)(cid:2993) (cid:1855) (cid:3038)(cid:3097)(cid:3038)(cid:3097) +△∑ (cid:1855) (cid:3038)↑(cid:2993) (cid:1855) (cid:2879)(cid:3038)↓(cid:2993)(cid:3038)(cid:3097) , where we assume that the substrate density of states in the normal state is constant within a bandwidth of 2D and △ is the energy gap of the superconductor. The tunneling Hamiltonian is given by (cid:1834) (cid:3021) = ∑ (cid:2021) (cid:3038) (cid:3435)(cid:1855) (cid:3038)(cid:3097)(cid:2993) (cid:1856) (cid:3036)(cid:3097) + ℎ. (cid:1855). (cid:3439) (cid:3038)(cid:3097) and the tunneling amplitude is given by Γ (cid:3036) ((cid:2033)) = (cid:2024) ∑ |(cid:2029) (cid:3036)(cid:3038) | (cid:2870) (cid:2012)((cid:2033) − (cid:2035) (cid:3038) ) (cid:3038) . Here we drop the energy dependency of the tunneling amplitude 4 and assume that it is constant over the energy window Γ (cid:3036) ((cid:2033)) = Γ (cid:3036) = 0.016(cid:1830) . The Hamiltonian for the impurity is (cid:1834) (cid:3005) = ∑ (cid:2035) (cid:3031) (cid:1856) (cid:3097)(cid:2993) (cid:1856) (cid:3097)(cid:3097) , where (cid:2035) (cid:3031) is the energy level of the individual impurity. We solve this model using Slave-boson mean-field theory (SBMFT) . We restrict the Hilbert space such that there is no double occupancy in the impurity due to a large Coulomb repulsion. This is done by introducing a bosonic and fermionic degree of freedom for each electronic operator, namely (cid:1854) (cid:3036)(cid:3097)(cid:2993) = (cid:1858) (cid:3036)(cid:3097)(cid:2993) (cid:1854) (cid:3036) , where (cid:1858) (cid:3036)(cid:3097) is fermion and (cid:1854) (cid:3036) is boson creation operator. No double occupancy is achieved by the constraint (cid:1854) (cid:3036)(cid:2993) (cid:1854) (cid:3036) + ∑ (cid:1858) (cid:3036)(cid:3097)(cid:2993) (cid:1858) (cid:3036)(cid:3097)(cid:3036)(cid:3097) = 1 . Within the mean-field approximation, the bosonic operator is replaced by a number and the constraint is fulfilled on an average by introducing a chemical potential. In the large T K limit it is favorable to break the Cooper pair such that the impurity spin can be screened. SBMFT is supposed to be valid in this regime , and thus, in particular, for our experimental setup k B T K /∆ >> 1. The spectrum of a single magnetic impurity in a superconductor consists of a continuum at higher energy and a bound state within the gap. This bound state appears close to the Fermi energy for small T K and shifts towards higher energy for larger T K . Figure 3a shows the evolution of bound states for different values of T K . 5 Acknowledgements
The authors thank Markus Ternes for fruitful discussions. This work was supported by the European Research Council Advanced Grant ADMIRE (project no. 786020). R.W. and J.W. acknowledge funding by the Cluster of Excellence ‘Advanced Imaging of Matter’ (EXC 2056 - project ID 390715994) of the Deutsche Forschungsgemeinschaft (DFG). R.Ž. acknowledges support by the Slovenian Research Agency (ARRS) under Program No. P1-0044. M.V. and R.M. acknowledge financial support from the Deutsche Forschungsgemeinschaft (DFG) through Project No. SFB668-A3 and A1, and from the European Union’s Horizon2020 research and innovation program under grant agreement No. 696656 - GrapheneCore1. The work of M.V. was additionally funded through 2019 Equal Opportunity Fund of University of Hamburg. The DFT computations were performed with resources provided by the North-German Supercomputing Alliance (HLRN). S.P. and J.F thank Stiftelsen Olle Engqvist Byggmästare and Vetenskapsrådet for financial support.
Author contributions
A.K., L.C. and J.W. conceived and designed the experiments. A.K. and L.C. performed the experiments. A.K., J.W. and R.Ž. analyzed the experimental data. NRG calculations were performed by R.Ž.. M.V. and R.M. carried out the DFT calculations. S.P. and J.F. performed SBMFT calculations. A.K. wrote the manuscript. All authors discussed the results and commented on the manuscript. 6
Figures
Figure 1 | Fe adatom and interstitial Fe assemblies at the TaO surface. a , Constant-current STM image of the TaO surface with an Fe adatom and an interstitial Fe atom (IFA) indicated by the red arrow. b , Height profile across the IFA and the Fe adatom in a . c, Perspective view of the DFT-optimized fully relaxed positions of an Fe adatom, an IFA, and the atoms in the TaO surface. d , e , Schematic top-view ( d ) and the corresponding STM images ( e ) of an Fe adatom with assemblies of the given different numbers of IFAs ( V = 50 mV, I = 100 pA). In d , the green lines represent the rows of oxygen atoms, and the larger (smaller) spheres represent the adatom (the IFAs). 7 Figure 2 | dI/dV spectroscopy on Fe adatoms for various interstitial Fe assemblies. a , Large-bias dI/dV spectra (open dots) showing the Kondo resonance. The dashed lines are fits to a Fano function, which is multiplied by a Lorentzian dip for the spectra exhibiting the zero-bias anomaly (ZBA), i.e. the exchange splitting of the Kondo resonance (see Methods). The pictograms and numbers on the left side of the spectra show the corresponding structure and numbers of the IFA assemblies. Spectra are shifted vertically for clarity. Stabilization parameters for STS: V stab = 100 mV, I stab = 100 pA, V mod = 2 - 4 mV. b , Small-bias dI/dV spectra (open dots) showing peaks inside the superconducting gap due to YSR states. Dashed lines over each spectrum are fits obtained using a numerical deconvolution method (see Methods). Spectra are shifted vertically for clarity. Stabilization parameters for STS: V stab = 2.5 mV, I stab = 100 pA, V mod = 20 μV. The corresponding deconvoluted LDOS is plotted in c . The green curve in c corresponds to the reference substrate density of states. 8 Figure 3 | Correlation of Kondo temperature, exchange splitting, and YSR state energy. a , Scatter plot of the energy of the larger YSR peak versus the Kondo temperature, both extracted from the experimentally investigated Fe adatom with different IFA assemblies. Various symbols in the plot represent various types of IFA assemblies as indicated in c . Horizontal error bars indicate the estimated errors in the fitting of the Kondo-resonance width. The violet curve in a represents the slave-boson mean-field theory calculation. b , Scatter plot of the energy of the larger YSR peak versus the exchange splitting of the Kondo resonance, both extracted from the experimentally investigated Fe adatom with different IFA assemblies (assignment of the symbols to the different assemblies, see c ). Results of the NRG calculations, which are plotted as open black dots, showing qualitative agreement with the experiments. All energy scales are normalized to the energy gap (∆) of the superconducting substrate. 9 Figure 4 | Numerical renormalization group calculations. a , Spectral function of the impurity representing the Fe adatom, which is magnetically coupled to a second impurity representing the IFA assembly, on the energy scale of the Kondo peak for various magnetic interaction strengths J . The superconducting gap is washed out here using a convolution which mimics the experimental lock-in-modulation broadening. Red and blue colors indicate the nature of the ground state of the system (Kondo vs. AFM) controlled by the strength of J . b , Discrete (sub-gap) part of the spectral function; a narrow Gaussian kernel is used to broaden the delta peaks. c , Energy of the YSR state as a function of J . d , Width of the exchange-splitting gap as a function of J . The insets illustrate the two regimes of decoupled Kondo singlets (small J ) and AFM dimer (large J ). 0 References
1. Nadj-Perge, S. et al.
Observation of Majorana fermions in ferromagnetic atomic chains on a superconductor.
Science , 602–607 (2014). 2. Ruby, M., Heinrich, B. W., Peng, Y., von Oppen, F. & Franke, K. J. Exploring a Proximity-Coupled Co Chain on Pb(110) as a Possible Majorana Platform.
Nano Lett. , 4473–4477 (2017). 3. Kim, H. et al. Toward tailoring Majorana bound states in artificially constructed magnetic atom chains on elemental superconductors.
Sci. Adv. , eaar5251 (2018). 4. Kamlapure, A., Cornils, L., Wiebe, J. & Wiesendanger, R. Engineering the spin couplings in atomically crafted spin chains on an elemental superconductor. Nat. Commun. , 3253 (2018). 5. Choy, T.-P., Edge, J. M., Akhmerov, A. R. & Beenakker, C. W. J. Majorana fermions emerging from magnetic nanoparticles on a superconductor without spin-orbit coupling. Phys. Rev. B , 195442 (2011). 6. Martin, I. & Morpurgo, A. F. Majorana fermions in superconducting helical magnets. Phys. Rev. B , 144505 (2012). 7. Braunecker, B. & Simon, P. Interplay between Classical Magnetic Moments and Superconductivity in Quantum One-Dimensional Conductors: Toward a Self-Sustained Topological Majorana Phase. Phys. Rev. Lett. , 147202 (2013). 8. Klinovaja, J., Stano, P., Yazdani, A. & Loss, D. Topological Superconductivity and Majorana Fermions in RKKY Systems.
Phys. Rev. Lett. , 186805 (2013). 9. Nadj-Perge, S., Drozdov, I. K., Bernevig, B. A. & Yazdani, A. Proposal for realizing Majorana fermions in chains of magnetic atoms on a superconductor.
Phys. Rev. B , 1 020407 (2013). 10. Pientka, F., Glazman, L. I. & von Oppen, F. Topological superconducting phase in helical Shiba chains. Phys. Rev. B , 155420 (2013). 11. Vazifeh, M. M. & Franz, M. Self-Organized Topological State with Majorana Fermions. Phys. Rev. Lett. , 206802 (2013). 12. Ruby, M. et al.
End States and Subgap Structure in Proximity-Coupled Chains of Magnetic Adatoms.
Phys. Rev. Lett. , 197204 (2015). 13. Pawlak, R. et al.
Probing atomic structure and Majorana wavefunctions in mono-atomic Fe chains on superconducting Pb surface. npj Quantum Inf. , 16035 (2016). 14. Schecter, M., Flensberg, K., Christensen, M. H., Andersen, B. M. & Paaske, J. Self-organized topological superconductivity in a Yu-Shiba-Rusinov chain. Phys. Rev. B , 140503 (2016). 15. Feldman, B. E. et al. High-resolution studies of the Majorana atomic chain platform.
Nat. Phys. , 286–291 (2017). 16. Jeon, S. et al. Distinguishing a Majorana zero mode using spin-resolved measurements.
Science , 772–776 (2017). 17. Röntynen, J. & Ojanen, T. Topological Superconductivity and High Chern Numbers in 2D Ferromagnetic Shiba Lattices.
Phys. Rev. Lett. , 236803 (2015). 18. Li, J. et al.
Two-dimensional chiral topological superconductivity in Shiba lattices.
Nat. Commun. , 12297 (2016). 19. Rachel, S., Mascot, E., Cocklin, S., Vojta, M. & Morr, D. K. Quantized charge transport in chiral Majorana edge modes. Phys. Rev. B , 205131 (2017). 2 20. Ménard, G. C. et al. Two-dimensional topological superconductivity in Pb/Co/Si(111).
Nat. Commun. , 2040 (2017). 21. Palacio-Morales, A. et al. Atomic-scale interface engineering of Majorana edge modes in a 2D magnet-superconductor hybrid system.
Sci. Adv. , eaav6600 (2019). 22. Yu, L. Bound State in Superconductors With Paramagnetic Impurities. Acta Phys. Sin. , 75 (1965). 23. Shiba, H. Classical Spins in Superconductors. Prog. Theor. Phys. , 435–451 (1968). 24. Rusinov, A. I. On the Theory of Gapless Superconductivity in Alloys Containing Paramagnetic Impurities. Sov. J. Exp. Theor. Phys. , 1101 (1969). 25. Malavolti, L. et al. Tunable Spin–Superconductor Coupling of Spin 1/2 Vanadyl Phthalocyanine Molecules.
Nano Lett. , 7955–7961 (2018). 26. Farinacci, L. et al. Tuning the Coupling of an Individual Magnetic Impurity to a Superconductor: Quantum Phase Transition and Transport.
Phys. Rev. Lett. , 196803 (2018). 27. Liebhaber, E. et al.
Yu-Shiba-Rusinov states in the charge-density modulated superconductor NbSe2. arXiv
Nano Lett. , 2311–2315 (2018). 29. Flatté, M. E. & Reynolds, D. E. Local spectrum of a superconductor as a probe of interactions between magnetic impurities. Phys. Rev. B , 14810–14814 (2000). 30. Morr, D. K. & Stavropoulos, N. A. Quantum interference between impurities: Creating novel many-body states in s -wave superconductors. Phys. Rev. B , 020502 (2003). 3 31. Yao, N. Y. et al. Phase diagram and excitations of a Shiba molecule.
Phys. Rev. B , 241108 (2014). 32. Ruby, M., Heinrich, B. W., Peng, Y., von Oppen, F. & Franke, K. J. Wave-Function Hybridization in Yu-Shiba-Rusinov Dimers. Phys. Rev. Lett. , 156803 (2018). 33. Ji, S.-H. et al.
High-Resolution Scanning Tunneling Spectroscopy of Magnetic Impurity Induced Bound States in the Superconducting Gap of Pb Thin Films.
Phys. Rev. Lett. , 226801 (2008). 34. Choi, D.-J. et al.
Influence of Magnetic Ordering between Cr Adatoms on the Yu-Shiba-Rusinov States of the β-Bi2Pd Superconductor.
Phys. Rev. Lett. , 167001 (2018). 35. A. Hewson.
The Kondo Problem to Heavy Fermions . (Cambridge University Press, 1993). 36. Franke, K. J., Schulze, G. & Pascual, J. I. Competition of Superconducting Phenomena and Kondo Screening at the Nanoscale.
Science , 940–944 (2011). 37. Bauer, J., Pascual, J. I. & Franke, K. J. Microscopic resolution of the interplay of Kondo screening and superconducting pairing: Mn-phthalocyanine molecules adsorbed on superconducting Pb(111).
Phys. Rev. B , 075125 (2013). 38. Hatter, N., Heinrich, B. W., Ruby, M., Pascual, J. I. & Franke, K. J. Magnetic anisotropy in Shiba bound states across a quantum phase transition. Nat. Commun. , 8988 (2015). 39. Jones, B. A. & Varma, C. M. Study of two magnetic impurities in a Fermi gas. Phys. Rev. Lett. , 843–846 (1987). 40. Jones, B. A., Varma, C. M. & Wilkins, J. W. Low-Temperature Properties of the Two-Impurity Kondo Hamiltonian. Phys. Rev. Lett. , 125–128 (1988). 4 41. Wahl, P. et al. Exchange Interaction between Single Magnetic Adatoms.
Phys. Rev. Lett. , 056601 (2007). 42. Spinelli, A. et al. Exploring the phase diagram of the two-impurity Kondo problem.
Nat. Commun. , 10046 (2015). 43. Jayaprakash, C., Krishna-murthy, H. R. & Wilkins, J. W. Two-Impurity Kondo Problem. Phys. Rev. Lett. , 737–740 (1981). 44. Prüser, H. et al. Interplay between the Kondo effect and the Ruderman–Kittel–Kasuya–Yosida interaction.
Nat. Commun. , 5417 (2014). 45. Bork, J. et al. A tunable two-impurity Kondo system in an atomic point contact.
Nat. Phys. , 901–906 (2011). 46. Žitko, R. Numerical subgap spectroscopy of double quantum dots coupled to superconductors. Phys. Rev. B , 165116 (2015). 47. Žitko, R., Bodensiek, O. & Pruschke, T. Effects of magnetic anisotropy on the subgap excitations induced by quantum impurities in a superconducting host. Phys. Rev. B , 054512 (2011). 48. Yao, N. Y., Glazman, L. I., Demler, E. A., Lukin, M. D. & Sau, J. D. Enhanced Antiferromagnetic Exchange between Magnetic Impurities in a Superconducting Host. Phys. Rev. Lett. , 087202 (2014). 49. Saldaña, J. C. E. et al.
Two-Impurity Yu-Shiba-Rusinov States in Coupled Quantum Dots. arXiv et al.
Andreev molecules in semiconductor nanowire double quantum dots.
Nat. Commun. , 585 (2017). 5 51. Cornils, L. et al. Spin-Resolved Spectroscopy of the Yu-Shiba-Rusinov States of Individual Atoms.
Phys. Rev. Lett. , 197002 (2017). 52. Mozara, R. et al.
Atomically thin oxide layer on the elemental superconductor Ta(001) surface.
Phys. Rev. Mater. , 094801 (2019). 53. Fano, U. Effects of Configuration Interaction on Intensities and Phase Shifts. Phys. Rev. , 1866–1878 (1961). 54. Coleman, P. New approach to the mixed-valence problem.
Phys. Rev. B , 3035–3044 (1984). 55. Avishai, Y., Golub, A. & Zaikin, A. D. Superconductor-quantum dot-superconductor junction in the Kondo regime. Phys. Rev. B , 041301 (2003). 56. Borkowski, L. S. & Hirschfeld, P. J. Low-temperature properties of anisotropic superconductors with kondo impurities. J. Low Temp. Phys. , 185–206 (1994). 57. Simon, P., López, R. & Oreg, Y. Ruderman-Kittel-Kasuya-Yosida and Magnetic-Field Interactions in Coupled Kondo Quantum Dots. Phys. Rev. Lett. , 086602 (2005). 58. Heersche, H. B. et al. Kondo Effect in the Presence of Magnetic Impurities.
Phys. Rev. Lett. , 017205 (2006). 59. Craig, N. J. et al. Tunable nonlocal spin control in a coupled-quantum dot system.
Science , 565–7 (2004). 60. Žitko, R. & Pruschke, T. Energy resolution and discretization artifacts in the numerical renormalization group.
Phys. Rev. B , 085106 (2009). 61. Bulla, R., Costi, T. A. & Pruschke, T. Numerical renormalization group method for quantum impurity systems. Rev. Mod. Phys. , 395–450 (2008). 6 62. Khajetoorians, A. A. et al. Tuning emergent magnetism in a Hund’s impurity.
Nat. Nanotechnol. , 958–964 (2015). 63. Kresse, G. & Furthmüller, J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Phys. Rev. B , 11169–11186 (1996). 64. Blöchl, P. E. Projector augmented-wave method. Phys. Rev. B , 17953–17979 (1994). 65. Perdew, J. P., Burke, K. & Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. , 3865–3868 (1996). 66. Dudarev, S. L., Botton, G. A., Savrasov, S. Y., Humphreys, C. J. & Sutton, A. P. Electron-energy-loss spectra and the structural stability of nickel oxide: An LSDA+U study. Phys. Rev. B , 1505–1509 (1998). 67. Cococcioni, M. & de Gironcoli, S. Linear response approach to the calculation of the effective interaction parameters in the LDA+U method. Phys. Rev. B , 035105 (2005). 68. Tkatchenko, A. & Scheffler, M. Accurate Molecular Van Der Waals Interactions from Ground-State Electron Density and Free-Atom Reference Data. Phys. Rev. Lett. , 073005 (2009). 69. Bučko, T., Lebègue, S., Hafner, J. & Ángyán, J. G. Improved Density Dependent Correction for the Description of London Dispersion Forces.
J. Chem. Theory Comput. , 4293–4299 (2013). 70. Bučko, T., Lebègue, S., Ángyán, J. G. & Hafner, J. Extending the applicability of the Tkatchenko-Scheffler dispersion correction via iterative Hirshfeld partitioning. J. Chem. Phys. , 034114 (2014). 71. Tkatchenko, A., DiStasio, R. A., Car, R. & Scheffler, M. Accurate and Efficient Method 7 for Many-Body van der Waals Interactions.
Phys. Rev. Lett. , 236402 (2012). 72. Bučko, T., Lebègue, S., Hafner, J. & Ángyán, J. G. Tkatchenko-Scheffler van der Waals correction method with and without self-consistent screening applied to solids.
Phys. Rev. B , 064110 (2013). 73. Coleman, P. Mixed valence as an almost broken symmetry. Phys. Rev. B , 5072–5116 (1987). 74. Schwab, P. & Raimondi, R. Andreev tunneling in quantum dots: A slave-boson approach. Phys. Rev. B , 1637–1640 (1999). 75. Aono, T., Eto, M. & Kawamura, K. Conductance through Quantum Dot Dimer Below the Kondo Temperature. J. Phys. Soc. Japan , 1860–1863 (1998). Supplementary Figures
Supplementary Figure 1 | dI/dV spectroscopy on Fe adatoms of additional Fe clusters. a,
Schematic top view diagram and corresponding STM images of different clusters of Fe atoms. b , Large bias range dI/dV spectra (open dots) showing the Kondo resonance. The dashed lines are fits to a Fano function, which is multiplied by a Lorentzian dip (see Methods). Stabilization parameters for STS: V stab = 100 mV, I stab = 100 pA, V mod = 2 - 4 mV. c , Small bias range dI/dV spectra measured with a superconducting tip (open dots), showing the peaks due to the YSR states. Stabilization parameters for STS: V stab = 2.5 mV, I stab = 100 pA, V mod = 20 μV. 9 Supplementary Figure 2 | Asymmetry of the YSR peaks as a function of bound state energy.
We observe that the intensities of the YSR peaks are asymmetric and as a general trend reverse sign with increasing number of IFAs in the assembly. To characterize this effect, we define the asymmetry as:
𝐴𝑠𝑦𝑚𝑚𝑒𝑡𝑟𝑦 = ℎ 𝑅 −ℎ 𝐿 ℎ 𝑅 +ℎ 𝐿 , where ℎ 𝑅 (ℎ 𝐿 ) is the amplitude of the YSR peak on the positive (negative) bias side. The corresponding experimental values are shown as filled symbols, where different symbols represent the various Fe clusters defined in Figure 3c of the main manuscript. The plot also shows the results of the NRG calculations as open dots for comparison, showing qualitative agreement with the experimental trend. 0 Supplementary Figure 3 | Comparison of YSR energy versus Kondo temperature with literature.
Blue data points and black lines are adapted from the experimental data and numerical renormalization group (NRG) calculations, respectively, for the system of MnPc molecules on Pb(111) . Brown filled circles and violet line are our experimental data points and Slave boson mean-field theory (SBMFT) calculations, respectively, for the system of Fe clusters on TaO. Supplementary Figure 4 | Interstitial Fe atom (IFA) within DFT approach. a,
Suggested positions with one IFA per (3 x 3) plaquette (initial sites) that are marked by letters (i-iv). b, Total energy comparison per unit cell as obtained within GGA+U (U
Fe-d | U
O-p ) simulations and compared between different compositions (i-iv). We have also considered an additional position of an IFA in a hollow cite of the Ta island. The latter configuration leads to significant reconstruction of the surface and a loss of the (3 x 3) plaquette structure, so we dismissed it. Among the other positions, one can see that, after relaxation, the IFA that was initially located on (iii) converged to location (i). In this setup, the IFA occupies approximately the center of the Ta octahedra with irregular base and two apexes, connected by the bridge oxygen atom. Among all considered configurations, only case (i) satisfies the criteria of lowest energy and smallest distortion of the surface (3 x 3) plaquette structure. To control the charge localization in the system, we used an on-site Coulomb parameter for Fe and O states, which was found to improve the convergence of the electronic and ionic minimizations significantly but did not change the structure qualitatively. The main surface change, induced by the presence of such an IFA in the first interlayer, is an elevated oxygen atom in the bridge position together with a slight upwards shift of the nearest Ta atoms. This change would be reflected in an STM image as an enhancement of the cross-like shape of the (3 x 3) plaquette near the IFA , which is indeed seen in the experiment (see Fig. 1a of the main manuscript). 2 Supplementary Figure 5 | VdW approximations for the Fe adatom at the center of the (3 x 3) plaquette.
Comparisons of the heights and total energies (in the slab with two surfaces) for the system of one Fe adatom at the center of the (3 x 3) plaquette with respect to different vdW frameworks, implemented in the VASP package: no vdW (no van-der-Waals interaction), TS (Tkatchenko-Scheffler method), TS/HI (TS scheme with iterative Hirschfield partitioning) and TS/HI-SCS (self-consistent scheme of TS/HI). The Fe adatom was allowed to move only along the z -direction, while the two in-plane adatom coordinates were fixed to the center position resulting from the full free relaxation of the surface within the TS/HI scheme. The reduction of the repulsive interaction between the adatom due to the dispersion forces, acting between the adatom and the polarizable substrate, is clearly visible here. This effect indicates the presence of a strong dynamical dipole-dipole interaction, that should be accounted for in the adsorption process on such surfaces, as it was predicted earlier . Within a full surface relaxation, the screening TS/HI-SCS method or switching the vdW part off (no vdW) lead to the relaxation of the adatom to the nearest hollow position, which is off-center. We attribute this error to a rather effective way of the inclusion of dynamics of the interaction between fluctuating charges as part of a pseudopotential and its possible reduction in the DFT-vdW-SCS method. The obtained height of the Fe adatom over the surface is approx. 2.5 Å (TS/HI), that fits to the experimental results . 3 Supplementary Figure 6 | Magnetization of the Fe adatom in the presence of an IFA.
Density of states for the d -states of the Fe adatom ( a ) and the IFA ( b ) in the slab with both Fe atoms. Here, continuous lines refer to the slab with 2 atoms (Fe adatom and IFA) and dashed lines refer to the slab with one Fe atom in absence of the other (either adatom or IFA). Red and blue colors denote spin-up and spin-down densities, respectively. The position of the adatom was extracted from the vdW calculations , while the IFA occupies the space beneath the bridge oxygen atom (case (i) in Supplementary Figure 4). The system is relaxed within the spin-polarized GGA+U scheme with the same parameters used for the simulations with the IFA. The distance between the IFA and the adatom is 6.5 Å after relaxation. One can see that the IFA has a small but non-zero magnetic moment (0.9 𝜇 𝐵 ) which does not change upon removal of the Fe adatom. In contrast, the Fe adatom has a large magnetic moment of 3.1 𝜇 𝐵 , and a splitting of the peaks near the Fermi energy appears due to the magnetic interaction with the IFA. Supplementary references
1. Bauer, J., Pascual, J. I. & Franke, K. J. Microscopic resolution of the interplay of Kondo screening and superconducting pairing: Mn-phthalocyanine molecules adsorbed on superconducting Pb(111).
Phys. Rev. B , 075125 (2013). 2. Mozara, R. et al. Atomically thin oxide layer on the elemental superconductor Ta(001) surface.
Phys. Rev. Mater. , 094801 (2019). 3. Cornils, L. et al. Spin-Resolved Spectroscopy of the Yu-Shiba-Rusinov States of Individual Atoms.
Phys. Rev. Lett. , 197002 (2017). 4. Kamlapure, A., Cornils, L., Wiebe, J. & Wiesendanger, R. Engineering the spin couplings in atomically crafted spin chains on an elemental superconductor.
Nat. Commun.9