Mapping the orbital structure of impurity bound states in a superconductor

A magnetic atom inside a superconductor locally distorts superconductivity. It scatters Cooper pairs as a potential with broken time-reversal symmetry, leading to localized bound states with subgap excitation energies, named Shiba states. Most conventional approaches regarding Shiba states treat magnetic impurities as point scatterers with isotropic exchange interaction. Here, we show that the number and the shape of Shiba states are correlated to the spin-polarized atomic orbitals of the impurity, hybridized with the superconductor. Using scanning tunnelling spectroscopy, we spatially map the five Shiba excitations found on subsurface chromium atoms in Pb(111), resolving their particle and hole components. While particle components resemble d orbitals of embedded Cr atoms, hole components differ strongly from them. Density functional theory simulations correlate the orbital shapes to the magnetic ground state of the atom, and identify scattering channels and interactions, all valuable tools for designing atomic-scale superconducting devices.

Y u-Shiba-Rusinov (Shiba) states [1][2][3] are identified in scanning tunnelling spectra as pairs of intra-gap resonances symmetrically positioned around zero-bias [4][5][6][7][8][9] . Each resonance corresponds to the injection of an electron or hole into the bound state 8 , thus representing the particle or hole components of the quasiparticle (QP) wavefunction. Since Shiba excitations lie inside the superconducting gap, their lifetime considerably exceed that of other QPs. This anticipates that Shiba peaks behave as a robust probe of scattering phenomena in superconductors, revealing intrinsic properties such as the distribution of the order parameter 5 , the QP band structure 10 , the effect of dimensionality 9 or Andreev tunnelling processes 8 . Owing to their long lifetime, Shiba excitations exhibit narrow lineshapes, which enables the study of magnetic phenomena in impurities with high energy resolution, such as magnetic anisotropy 11,12 and magnetic coupling 6,13 .
Theoretical models using classical spins with isotropic exchange fields have predicted multiple Shiba bound states, labelled by an angular momentum quantum number l (ref. 3). In fact, magnetic transition metal (TM) atoms in a superconductor show a varying number of Shiba bound states depending on the nature of both element and superconductor 5,6,8,9,14,15 . However, the exchange fields of atomic scatterers are not isotropic. Due to the orbital character of the scattering channels, the spin can only be exchange scattered by spin-polarized atomic states, which for 3d-TM atoms are those defined by the angular momentum l ¼ 2. The symmetry reduction by crystal fields, and hybridization with the host atomic lattice lift the degeneracy of the l ¼ 2 subshell. In this scenario, Shiba multiplets are predicted to appear, thus reflecting the occupation level of the atomic shell 16 . The resolution of their orbital character would render them as the ideal probe for identifying the magnetic ground state of a single impurity in a superconductor. Furthermore, the spatial extension of Shiba states also provides precise access to the coupling between impurities 6 and to basic properties of the superconducting host 5,9,15,17 .
Here we explore possible orbital components in Shiba states. We compare the QP local density of states (LDOS) of chromium atoms deposited on Pb(111), mapped with a low-temperature scanning tunnelling microscope (STM), with atomistic simulations based on density functional theory (DFT). We find that the atoms are embedded under the first Pb layer and show five intra-gap states due to QP scattering with the five d orbitals of Cr. The spectral maps of the Shiba states are very localized on the impurity and show a characteristic particle-hole asymmetry. Based on theoretical maps of particle and hole components we interpret their orbital origin and obtain hints on the interaction strength for each channel.

Results
Multiple Shiba states of Cr atoms. We deposited Cr atoms on a Pb(111) film grown on SiC(0001) (thickness \100 nm) at 15 K. The atoms appear in STM images (measured at a base temperature of 1.2 K) as protrusions with a small apparent height of B50 pm (Fig. 1a). DFT calculations reveal that the most stable configuration corresponds to the Cr atom underneath the uppermost Pb layer, as sketched in Fig. 1b,c (Methods section and Supplementary Note 2). The barrier for reaching this state is as low as 21 meV, which is smaller than, for example, the surface lateral diffusion barriers (Supplementary Table 1). Since the Cr atoms arrive hot to the Pb(111) surface, they can easily access the embedded site before they thermalize to the substrate temperature. The simulated STM topography agrees with the small experimental height of the STM images.
The QP LDOS associated to subsurface Cr impurities was obtained from differential tunnelling conductance (dI/dV) spectra using superconducting Pb-terminated tips to increase the energy resolution beyond the thermal limit 7 . Figure 2a shows a typical dI/dV spectrum acquired on bare Pb(111) that exhibits a doubled superconducting (SC) gap of Pb with sharp coherence peaks at ±2D ¼ ±2.7 meV (Methods section). The dI/dV curve taken on a Cr atom in Fig. 2b reveals a rich spectral structure inside the superconducting gap, with six intra-gap peaks at each polarity (that is, six peak pairs), in addition to the original Pb coherence peaks. To obtain the QP LDOS, we deconvoluted the experimental dI/dV spectrum on chromium ( peaks, labelled as E h;p n (n ¼ 1, 2, 3, 4, 5) in Fig. 2d. We attribute these pairs of peaks to five Shiba states with excitation energies The additional pair of peaks, denoted V h 6 and V p 6 in the dI/dV spectrum of Fig. 2d, appears at bias voltages ± (D À E 5 )/e and, thus, corresponds to a thermal replica of the Shiba pair E h 5 /E p 5 (ref. 8). Shiba multiplets arise naturally as different angular momentum components in isotropic exchange fields. However, the resolution of a fivefold Shiba multiplet achieved here enables us to reinterpret their origin in terms of atomic properties of the impurity 11,16 .
Spatial shape of Shiba states. To probe the orbital origin of the Shiba bound states, we map their intensity in a small region around the Cr impurity. Figure 2d shows dI/dV maps obtained for the same Cr atom at the voltages of all twelve peaks found inside the SC gap (Methods section). In contrast to the featureless protrusion in topography (inset of Fig. 2d), the conductance maps reveal various features around the Cr impurity. Each Shiba state appears with a different shape, which is smoother for the first two, and shows specially marked intra-atomic features for states E 3 and E 4 , and E 5 . Interestingly, the maps show a clear dependence with polarity, signalling a particle-hole asymmetry in the wavefunction of the corresponding Shiba states. This is particularly clear for states E p 3 and E p 4 , and E p 5 , which are markedly different than their corresponding hole states. Such spatial asymmetry is also observed in the thermal replicas at V h 6 ¼ À ðD À E p 5 Þ=e and V p 6 ¼ðD À E h 5 Þ=e, which show the shapes of their mirror states E p 5 and E h 5 , respectively. The different shape of particle and hole components of every state explains their different peak amplitude in point spectra like in Fig. 2b. In fact, dI/dV spectra averaged over the whole extension around the Cr impurity levels their intensity at both polarities ( Supplementary Fig. 4).
To find out the origin of the number and shape of the Shiba states, we simulated with DFT (Methods section) the scattering channels of the embedded Cr impurity and modelled their effect in producing Shiba bound states. The channels are very sensitive to atomic-scale details of the Cr atom and its environment. We employed the most stable Cr site in the subsurface plane described in Fig. 1b, and found that this site produced results compatible with the experiment.
The identification of Cr-derived scattering channels is difficult; the embedded atom does not lie in a fully symmetric position inside the Pb FCC crystal, while it induces a significant displacement on its first Pb neighbours. Consequently, the Cr d To facilitate their identification in the band's continuum, we calculated first the single-electron energy levels of the cluster formed by the Cr atom and the 7 first neighbours with apparent interaction (shown in Fig. 1c). We found that, at least in this reduced ensemble, there are five non-degenerate singly occupied states with clear d character centred in the Cr atom. Figure 3a plots isosurfaces of wavefunction amplitude for the five eigenstates inside the atomic cluster. Their shape partly resembles some of the original symmetries of d orbitals. This implies that these states will mostly exchange-scatter l ¼ 2 host electrons 18 . However, the degeneracy of the d multiplet is lifted and the m l orbitals are mixed by the anisotropic Pb crystal field. In the full Pb system, these states appear as broad resonances around À 2.6 eV (Fig. 3c shows the DOS of the Pb slab projected on the five cluster states of Fig. 3a) due to the mixture with sp bands of the lead film, but their spin polarization remains very strong. We deem that the intra-gap bound states are created by exchange scattering of host electrons into these five one-electron states.
Simulation of Shiba state images. The spatial shapes of the five Cr-states of the cluster, projected on the surface (Fig. 3b), show already some level of agreement with the experimental Shiba maps. For example, the maps of peaks E h 3 , E p 4 or E p 5 resemble the shapes of states Cr 3 , Cr 4 or Cr 5 , respectively. This validates the subsurface atomic configuration obtained by DFT because the mixture of Cr-d orbitals is very sensitive to fine details of the atomic crystal field environment. To correlate orbital states with subgap impurity bound states, we computed the corresponding Shiba states by considering that each of them scatters with conduction electrons via a one-electron potential. The effect of an impurity on the total Hamiltonian can be divided in a potential scattering term (K(r), no spin degrees of freedom) and a spin contribution 3 , such that where J(r) is the exchange-coupling between the impurity spin S and the 1/2 spin of conduction electrons r. In a first approximation, we consider that the five cluster eigenstates of Fig. 3a also diagonalize these contributions, giving rise to five spin-polarized scattering channels. Since their degeneracy is lifted, they have a different exchange interaction J m with the host and produce a different bound state. They are then responsible for the appearance of five positive-energy (particle) and five negative-energy (hole) Shiba states in the QP LDOS, symmetrically aligned with respect to the Fermi energy 16 . The DFT results indicate that for this system the potential scattering terms K m in equation (1) are very small. In the absence of a potential scattering term, a large degree of particle-hole symmetry in the weight and amplitude of the Shiba states is expected 5,19 . This is indeed the case after adding up all differential conductance spectra measured over the surface region around the Cr impurity, as commented above and shown in Supplementary Fig. 4. The spatial distribution of Shiba amplitudes can be depicted by the squared modulus of the Bogoliubov QP coefficients, ju m ðrÞj 2 and jv m ðrÞj 2 , representing the LDOS of their particle and hole components, respectively 20 . Figure 3d plots the quantities ju m ðrÞj 2 and jv m ðrÞj 2 over the embedded Cr impurity, obtained as described in Supplementary Note 3. The results reproduce the different shape of particle and hole states from the experiments, especially for those resonances with marked intra-atomic features. We find that the particle components u m 2 resemble the shape of Cr d states in Fig. 3c. This is a result of the extended Fermi surface of lead, which provides a continuous set of scattering wave-vectors that partly reconstruct the shape of the scattering state. However, the hole components v m 2 strongly deviate for the orbital shapes. This is due to the dominating spin contribution to the potential: QP scattering with hole components produces a phase reversal for wavevectors of the Fermi surface, resulting in a clear distortion of the shape of scattering channels.

Discussion
We can recognize several features from the simulated particle and hole components of Bogoliubov QPs (Fig. 3d) in the experimental conductance maps of Shiba states (Fig. 2d). In particular, the particle components ju m j 2 of orbitals Cr 3 and Cr 5 resemble the positive bias maps of Shiba states E p 3 and E p 5 , respectively, whereas their jv m j 2 fits well with the corresponding negative-bias maps, E h 3 and E h 5 . This situation is indicative of an exchange constant J m smaller than the superconducting paring energy, so that tunnelling into the particle component of the bound superconducting pair requires extra positive energy. An unexpected particle-hole reversal is found for the Shiba bound state E 4 , where the calculated particle component (and orbital shape) clearly matches the dI/dV map at negative bias. When J m is larger than the pairing energy, the bound superconducting pair breaks and the new ground state is a Bogoliubov QP, while the scattering spin channel becomes screened (s ¼ 0). A fingerprint of this new ground state is an energy reversal of its particle and hole excitations 21 , while their crossing through zero energy signals the transition to this new ground state 7,12,21 . Thus, the resemblance of E h 4 to ju 4 j 2 denotes that the corresponding orbital experiences a larger interaction with the Pb bands (that is, a larger hybridization, compatible with our DFT results, Supplementary Note 4), becoming fully screened and not contributing to the total atomic spin.
The sharp maps of Shiba states change shape rapidly between peaks, in barely o200 meV, a much greater accuracy than any other method of orbital imaging one-electron states, and in spite of the large degree of hybridization of the buried atom. This is ensured by the superconducting gap in the QP spectrum, which keeps Shiba states with long excitation lifetimes and decoupled from other QPs (except indirectly via thermal effects and/or via Andreev processes 8 ). Moreover, Shiba states are fairly unperturbed by the presence of nearby impurities. Our measurements, shown in Supplementary Fig. 5, reveal that two impurities separated by only 4 surface lattice parameters (1.4 nm) display largely undisturbed Shiba conductance maps, what is attributed to the large localization of Shiba states in three dimensions 3,9 . This suggests that coupling impurity-induced bound states together into extended magnetic structures 6,22,23 requires impurity separations well in the sub-nm range. In these structures, identification of the characteristic orbital symmetries of impurity bound states is thus a unique fingerprint to unveil their (channel-specific) magnetic ground state, their magnetic coupling to other nearby impurities, and to follow the formation of extended Shiba bands.
Methods STM measurements. The experiments were conducted in a commercial SPECS GmbH Low-Temperature (1.2 K) STM, under Ultra-High Vacuum conditions. Pb(111) films (thickness \100 nm) on SiC(0001) substrates appear as crystalline grains with diameters 4300 nm, which show bulk-like superconductivity. To increase the energy resolution, the STM tip was repeatedly embedded in the Pb film until a superconducting tip was obtained. The full superconducting state of the tip was proved by performing STS spectra on bare Pb regions and showing that the two sharp coherence peaks appear at ± 2D, where D ¼ 1.35 meV, as in Fig. 2a. Contrary to single-crystal measurements 24 the spectra on bare Pb films never showed the characteristic double gap structure, probably due to the small film thickness compared to the Pb coherence length scale 25 , which washes out the anisotropy of the Fermi Surface 26 . Cr atoms were deposited on the Pb(111) surface at a surface temperature of B15 K. The shown spectra were obtained using a lockin amplifier, with modulation V rms ¼ 10 mV at 938.5 Hz. The Shiba conductance maps were obtained from a matrix of 52 Â 52 dI/dV(V) spectra measured in a region of 4.8 nm 2 over the Cr impurity. The corresponding topography images do not show any of the features observed in the Shiba maps (see inset of Fig. 2d). Analysis of STM and STS data were performed with the WS Â M 27 and SpectraFox 28 (http://www.spectrafox.com) software packages.
Theoretical details. Standard calculations using DFT were performed to reproduce the total energy and the electronic structure of different configurations of a Cr atom on Pb(111) (Supplementary Note 1). The calculations unequivocally predict a subsurface configuration under the surface bridge site, with considerable distortion of the surface Pb layer. A minimal Pb-Cr cluster is used to determine the local electronic structure near the Cr atom. Using the wavefunctions of this cluster, the Bogoliubov-de Gennes equations are approximated, in the presence of an exchange field, obtaining the equivalent of Rusinov's equations using generalized scattering channels and a non-free-electron normal metal. The equations are further approximated to yield the spatial distribution of the Shiba states, as detailed in the Supplementary Information. Data availability. The data that support the findings of this study are available from the corresponding authors upon request.