Effects of strong spin-orbit coupling on Shiba states from magnetic adatoms using first-principles theory

Magnetic impurities at surfaces of superconductors can induce bound states referred to as Yu–Shiba–Rusinov states (i.e. Shiba states) within superconducting (SC) gaps. For superconductors with strong spin–orbit coupling (SOC), Shiba states arising from even single magnetic adatoms are too complex to be fully understood using effective models alone because SOC cannot be treated perturbatively and multiple orbitals are strongly mixed with spin projections. Here we investigate Shiba states of single magnetic adatoms at the surface of strongly spin-orbit coupled SC Pb, by solving the fully relativistic Dirac–Bogoliubov–de Gennes equations using multiple scattering Green’s function methods. For Fe and Co adatoms on Pb(110), we show that the Shiba states are better characterized by total angular momentum, J, and its projections on the z axis, m J . As a hallmark of the SOC effect, the Shiba states show a strong dependence of the orientation of the adatom moment. As the orientation of the Fe/Co moment changes, the deepest Shiba states merge at zero energy. This zero-energy state disappears with an additional non-magnetic adatom next to the magnetic adatom, although the other Shiba states unchange. For a Mn adatom on Pb, our Shiba states overall agree with experiments. The characteristics of our Shiba states are also observed with the similar energies and characters in the experiments. The deepest Shiba states that we compute, however, do not appear as close to the Fermi level as the experimental data. It would be interesting to compute the Shiba states with continuously varying vertical distances of the Mn adatom from the surface or with varying the charge state of the adatom, and to calculate the spatial dependence of the spectral density. Our findings will be also useful for understanding of Shiba states for dimers and longer spin chains on the Pb surface considering noncollinear magnetic structures in them.


Introduction
When magnetic impurities are present at surfaces of superconductors, quasiparticle excitations referred to as Yu-Shiba-Rusinov states [1][2][3] (i.e. Shiba states) can appear within superconducting (SC) gaps. Spin-up and spin-down electrons in forming the Cooper pairs experience different exchange interaction with the magnetic moment of the impurity. This interaction weakens or even breaks the bonding between the Cooper pair elements, which results in states within the SC gap.
These excitations were studied in various pairing symmetries of SC states [4][5][6][7]. Being able to calculate the properties of Shiba states and to describe their formation conditions has been shown to be pivotal in understanding of topological end modes [8][9][10][11][12][13][14][15][16][17][18][19] in the presence of spin-orbit coupling (SOC) or spiral magnetic order. Very recently, a scanning tunneling microscopy/spectroscopy experiment [20] showed zero-energy end modes in dilute spin chains in proximity to an s-wave SC substrate. These end modes turned out to be topologically trivial Shiba states. Furthermore, Shiba states were proven to be very useful to probe magnetic nanostructures at surfaces with high resolution by using SC STM tips [21].
In many theoretical reports, these excitations (i.e. Shiba states) have been studied using single-orbital spin S = 1/2 effective models in the absence of SOC [1-4, 9, 22-26]. This approach does not capture major properties of Shiba states arising from transition-metal adatoms and chains, or Shiba states when the SC substrates have strong SOC. Recently, characteristics of Shiba states have been studied using a multi-orbital effective model consisting of five l = 2 orbitals in the presence of crystal field splitting without SOC [27].
In an effective-model approach for Shiba states, SOC was considered at different levels in different models [8,11,[28][29][30][31][32]. In [11], within a single-orbital model, Rashba SOC was included at the perturbative limit in ferromagnetic chains but not in the SC substrate, and the properties of Shiba states and Majorana zero-energy modes are discussed. In [28], Rashba SOC was taken into account in the SC substrate using a model comprising three l = 1 orbitals. It was found that the positions of Shiba states within the SC gap change as a function of SOC strength and direction of the magnetic moment [28]. In [29], SOC was inherently considered in terms of magnetic anisotropy for transition-metal magnetic impurities with high spin S > 1/2. This study showed phase diagrams of the number of Shiba states and the projection of the total spin as a function of pairing strength and uniaxial magnetic anisotropy parameter. All of the aforementioned model Hamiltonian studies are neither self-consistent calculations nor based on realistic parameters for the magnetic impurities and the SC substrates.
In order to consider the effects of SOC and multi-orbitals on Shiba states within self-consistent calculations with realistic band structures, we previously developed relativistic first-principles methods [33], and applied them to Shiba states arising from single Mn impurities and Mn dimers at a semi-infinite Nb substrate [33]. Within this methodology, we were able to characterize multiple Shiba states and showed that the splitting of Shiba states from ferromagnetic Mn dimers quantitatively agrees with experimental data from [34]. The first-principles study in [33] differs from the recent study using a similar first-principles approach [35], where quasiparticle excitations were studied for single magnetic impurities embedded deeply inside three-dimensional bulk Pb without SOC.
In this work, we investigate characteristics of Shiba states from single transition-metal magnetic adatoms (Fe, Co, Mn) on a Pb substrate as well as resonances from the neighboring Pb atoms, by performing first-principles calculations including SOC in the SC state. In our simulations, s, p, and d orbitals of the magnetic adatoms and the band structure of Pb are all included. We find that each Shiba state of d-orbital character exhibits strong coupling among different 3d orbitals and spin-up and spin-down projections due to strong SOC of the Pb substrate, in contrast to a SOC-free or weak-SOC case where each Shiba state corresponds to an individual 3d orbital. Therefore, we characterize Shiba states using the total angular momentum J and its projection along the magnetic moment direction m J , and show unique features of Shiba states in which SOC-induced spin mixing is strong. We study how the Shiba states are affected when the orientation of the adatom moment changes. We also discuss why the Shiba states from the Fe and Co adatoms appear deep inside the SC gap, and compare our results with available experimental data.

Methods
Regarding our first-principles calculations of the normal and the SC states, we use the methods developed in [33]. Nyári et al [33] summarizes all the essential formulas and necessary procedures to solve the Dirac-Bogoliubov-de Gennes (DBdG) equation based on first-principles multiple scattering theory for embedded magnetic impurities on surfaces of superconductors. Therefore, we here briefly describe only the relevant details of the implementation for systems of interest.
First we perform first-principles calculations of bulk Pb and a semi-infinite Pb slab in the normal state with the experimental lattice constant (4.95 Å [36]) using the fully relativistic screened Korringa-Kohn-Rostoker (SKKR) Green's function method [33,37] within Kohn-Sham-Dirac density functional theory in the atomic sphere approximation. For the Pb(110) and Pb(001) slab calculations, we consider four Pb atomic layers and four vacuum layers in the interface region treated self-consistently. For the Pb(111) slab calculation, we consider six Pb atomic layers and three vacuum layers in the interface region. The local spin density approximation [38] is used for exchange-correlation functional. Relativistic effects are taken into account in terms of the Kohn-Sham-Dirac equation. Scaling of the SOC is employed as proposed by Ebert et al [39]. In all calculations unless specified otherwise, 2450 k points are used with an orbital angular momentum cutoff of l max = 2 in the irreducible wedge of the two-dimensional Brillouin zone. The global x axis is along the [110] direction, and the global y axis is along the [001] direction of Pb(110) surface (see figure 1(a)). The global z axis is normal to the surface. The Fermi level of bulk Pb and vacuum potential of the Pb(110) slab are determined from self-consistent calculations.
We employ the embedded cluster method [40] self-consistently within the SKKR formalism to a single Fe, Co, or Mn magnetic adatom at the surface of the semi-infinite Pb(110) slab, as well as a Mn adatom on Pb(111) and Pb(001) slabs. For the Pb(110) slab, the Fe, Co, or Mn atom is located at the hollow site of the surface with the vertical distance of 1.75 Å from the Pb surface, unless specified otherwise. This distance is an interlayer distance of a pristine Pb(110) slab. Similarly, for the Pb(111) and P(001) slabs, the Mn adatom is located at the hollow site with the vertical distances of 2.86 and 2.48 Å, respectively. In the embedded cluster method, the system is divided into an impurity cluster (IC) and a host. The former consists of the magnetic impurity atom, neighboring Pb atoms, and surrounding empty spheres which represent vacuum layers, while the latter corresponds to the rest of the Pb atoms. Considering only nearest neighboring Pb atoms and vacuum sites to the magnetic adatom, we use a 13-site IC. Following the similar logic, for a Fe-S impurity, we use a 22-site IC, where sulfur (S) is located right next to the Fe site along the x axis in the same atomic layer. When the Fe impurity site is located at the Pb surface layer, we consider a 13-site IC. We confirm that our results do not change with a further increase of the IC size. Then using the SKKR electrostatic potential of the bare semi-infinite Pb slab, we compute the potential of the IC self-consistently in the normal state. Normal-state impurity density of states (DOS) plots are computed from Bloch spectral functions with 1250 k points, using a broadening of 1.0 mRy and an energy resolution of 5.0 mRy.
The converged normal-state potentials are used to solve the DBdG equations for the Green's function matrix [33,37,41]. We use the pairing potential ∆ of bulk and interface Pb layers as the experimental SC gap of bulk Pb, 0.1 mRy (=1.36 meV) [36], with s-wave pairing symmetry, while the pairing potential of the Fe, Co, or Mn atom is set to zero. The DOS ρ I (ϵ) of atom I, i.e. local DOS (LDOS), at energy ϵ in the IC is computed non-self consistently from the electron-and hole-components of the Green's function matrix G ee L,σ;L,σ and G hh L,σ;L,σ with a broadening of 0.5 µRy and an energy resolution of 1 µRy: where L and σ are orbital and spin indices, respectively. Here the integral runs over an atomic sphere around atom I with its Wigner-Seitz radius. Spin-decomposed DOS is calculated from projection of the total LDOS onto spin-up and down components along the rotating magnetic moment direction. Orbital characters are obtained from the global coordinates which do not change with the direction of the magnetic moment. We also decompose the SC-state LDOS into |J, m J ⟩ states using the {κ, µ} representation [41], where κ is the spin-orbit quantum number and µ = m J , when the magnetic moment is oriented along the z axis.

Properties of the normal state
We investigate properties of the normal state of the Fe adatom at the hollow site above the Pb surface (figure 1(a)) since they provide important insights into properties of Shiba states in the SC state. Let us first discuss the case that the Fe magnetic moment vector points normal to the surface (i.e. along the z axis). In this case, we obtain the Fe spin and orbital moments of 3.518 and 0.906 µ B , respectively, which significantly differ from those of an isolated Fe atom, i.e. 4.0 and 2.0 µ B , respectively. This is due to substantial charge transfer from the semi-infinite Pb bulk to the Fe adatom as well as orbital hybridization between them. Furthermore, small magnetic moments are induced in Pb atoms close to the Fe adatom. Each of the nearest-neighboring Pb atoms has a spin moment of 0.026 µ B with a direction opposite to the Fe spin moment. The rest of the Pb atoms in the IC have spin moments of an order of 0.001 µ B or less. Figure 1(b) shows spin-decomposed LDOS of the Fe adatom in the normal state. The Fe spin-up and spin-down LDOS spectra are broadened due to orbital hybridization between the Fe adatom and the Pb substrate. The energy difference between the spin-up DOS peak and spin-down DOS peak is about 0.212 Ry. The position of the spin-down DOS peak is slightly above the Fermi level E F . Considering the large peak-to-peak separation between the spin-up and spin-down DOS compared to their linewidths, we observe that the Fe spin-up LDOS is almost fully occupied and that the LDOS at E F predominantly has spin-down character. Figure 1(d) shows the orbital-decomposed spin-down LDOS of the Fe adatom. All d orbitals of the Fe adatom substantially contribute to the LDOS at E F , although d x 2 −y 2 orbital has the largest contribution. The comparable d orbitals are expected to strongly couple to the spin via hybridization with the strongly spin-orbit coupled Pb substrate. In addition to the d orbitals, there exists small LDOS of s-and p-orbital character with unequal spin-up and spin-down weights at the Fermi level. These features along with the charge-transfer effect may be difficult to include in an effective model to study Shiba states.
One important consequence of SOC is that electronic and magnetic properties can be affected by the direction of the Fe magnetic moment. We first discuss rotations of the magnetic moment by angle θ from the  Table 1. Normal-state spin and orbital moments of the transition-metal (TM) adatoms and a nearest-neighboring Pb atom, magnetic exchange energy, and existence of zero-energy state (ZES) for the rotations of the moment vectors. In the Fe and Co cases, two different vertical positions of the Fe or Co adatom are considered such as above the Pb surface (TM/Pb(110)) and at the Pb surface (TM-surface/Pb). The TM spin and orbital moments are computed when the TM magnetic moment points along the direction indicated in the parenthesis such as x, y, and z. The magnetic exchange energy is estimated from the spin-down LDOS peak to the spin-up LDOS peak in the normal state. The number next to ZES indicates the number of ZESs with the rotation of the TM moment vectors. z axis in the xz plane. Due to twofold rotation symmetry around the z axis and time-reversal symmetry, a quarter of the xz plane suffices to study the effect of moment rotation. The Fe spin-down LDOS at the Fermi level noticeably increases as θ increases (figure 1(g)) with a small change in the spin-up LDOS peak value (figures 1(b) and (c)). The orbital moment decreases to 0.780 µ B when the moment vector points along the x axis, although the spin moment is not affected (see table 1). Since the amount of charge transfer from the Pb substrate to the Fe adatom does not change with the moment orientation for a given vertical distance of the adatom, both the LDOS changes at the Fermi level and the sizeable decrease in the orbital moment imply remarkable changes in the hybridization with the Pb substrate as a function of the moment orientation. This scenario seems to be supported by a significant influence of the moment orientation on orbital-decomposed spin-down LDOS (figure 1(e) vs (d)). Now we discuss the moment rotation in the xy plane by angle ϕ from the x axis. The C 2v point-group symmetry dictates that a quarter of the xy plane suffices to examine the effect of moment rotation. The spin-down LDOS at the Fermi level is subtly modified as a function of rotation angle. The orbital moment slightly further decreases to 0.743 µ B when the moment vector points along the y axis, without any changes in the spin moment. Orbital-decomposed spin-down LDOS is more noticeably influenced by the Fe moment direction (figure 1(f) vs (d)). When SOC is turned off in both the Fe adatom and the Pb substrate, the Fermi level slightly increases by 8.12 mRy and the Fe spin magnetic moment is slightly lowered by 0.021 µ B compared to the SOC-included case. The spin-up Fe LDOS has a narrower linewidth than in the SOC-included case, and it is again well separated from the spin-down Fe LDOS (figure 1(h)). Although the Fermi level is shifted upwards, the spin-down LDOS at the Fermi level is almost the same as in the SOC-included case. Comparing figure 1(i) to figures 1(d)-(f), we find that there are enough differences between the orbital-decomposed LDOS values with and without SOC at the Fermi level, which suggests quite different Shiba states as discussed in sections 3.2 and 3.3.

Shiba states without SOC
Let us first discuss our calculated LDOS of the Fe impurity in the SC state without SOC. The total Fe LDOS, ρ e+h , reveal five symmetric pairs of peaks, i.e. Shiba states, at energy E/∆ = ±0.22, ±0.41, ±0.69, ±0.77, and ±0.94 relative to E F (zero energy), as shown in figure 2(a). The Shiba states at ±0.22∆ and ±0.41∆ originate from d xy and d z 2 orbitals, respectively, while the Shiba states at ±0.69∆ and ±0.77∆ arise from {d xz , d yz } and d x 2 −y 2 orbitals, respectively. For the Shiba states at ±0.69∆, the contributions of d xz and d yz orbitals are too close to be resolved within the numerical accuracy. The Shiba states at ±0.94∆ have very small amplitudes and they originate mainly from s orbital. Therefore, without SOC, each Shiba state is well characterized by a single orbital. The fact that the Shiba states of s-orbital character appear near the SC gap edges is consistent with a very small s-orbital LDOS at the Fermi level in the normal state leading to a very weak exchange coupling to the conduction electron spins. The normal-state LDOS values at the Fermi level are somewhat reflected in the ordering of the Shiba states, which agrees with the prediction by Shiba [2]. For example, in the normal state, at the Fermi level, the LDOS of d xy orbital is larger than that of d z 2 orbital which is larger than the LDOS of d xz orbital. In the SC state, the Shiba states of d xy orbital appear closest to the zero energy and the next two Shiba states closer to the zero energy are d z 2  77∆ comprise large spin-down hole components and very small spin-up electron components. The characteristics of the Shiba states at the corresponding positive energies are dictated by particle-hole symmetry. Without SOC, a given Shiba state at a negative or positive energy has predominantly either a spin-down electron or a spin-down hole component (figures 2(d)-(g)) rather than having comparable spin-down electron and spin-down hole components. As shown later, that is not the case when SOC is turned on.

Shiba states with SOC for moment direction along the z axis
Now we turn on SOC in the Fe adatom and the SC substrate when the moment vector is oriented along the z axis. Figure 2(h) shows a total Fe LDOS in the SC state. For convenience, five symmetric pairs of Shiba states are denoted as 1-5. The electron components at negative energies are labeled as 1-5, while those at positive energies are labeled as 1 ′ -5 ′ (figure 2(i)). Similarly, the hole components are denoted following particle-hole symmetry (figure 2(j)). We find that the energies of the Shiba states are distinct from those without SOC, except for pair 1. Comparing figure 2(a) to figure (h), we observe that with SOC, the four pairs 2-5 are shifted toward the SC gap edges and multiple orbitals comparably contribute to each of these Shiba states, in contrast to the SOC-free case. Interestingly, two pairs 3 and 5 have comparable spin-down electron and spin-down hole components at negative or positive energies, which is not the case for the Shiba states without SOC.
The above features are due to large mixing between orbital and spin degrees of freedom in the presence of strong SOC. In this case, it may be more convenient to analyze characters of the Shiba states using the total angular momentum, J, and its projection onto the z axis, m J (or more specifically |l, s = 1 2 ; J, m J ⟩, where l and s are the orbital and spin angular momenta). Table 2 summarizes dominant |l, 1 2 ; J, m J ⟩ contributions to the electron and hole components of each Shiba state. Pair 1 arises from |0, 1 2 ; 1 2 , ± 1 2 ⟩ and so it is not affected by SOC. Pair 2 originates mainly from |2, where m l and m s are projections of orbital and spin angular momenta along the z axis. Since only one spin channel is dominant, for this pair SOC-induced spin mixing is highly suppressed. The fact that only one of the states like 2 ′ has a large amplitude rather than both states (2, 2 ′ ) having substantial amplitudes, is similar to the SOC-free case. For pair 3, both |2, 1 2 ; 3 2 , + 3 2 ⟩ and |2, 1 2 ; 5 2 , + 3 2 ⟩ significantly contribute to state 3, while |2, 1 2 ; 5 2 , + 1 2 ⟩ and |2, 1 2 ; 3 2 , + 1 2 ⟩ predominantly contribute to state 3 ′ . Using the Clebsch-Gordon coefficients, state 3 can be rewritten as Similarly, state 3 ′ can be also rewritten as √ Therefore, both states 3 and 3 ′ have two spin channels with large amplitudes, and so strong SOC-induced spin mixing is plausible. As a result, the amplitudes of both states (3, 3 ′ ) are large. Similarly, for pair 5, strong spin-mixing induced by SOC is expected. For pair 4, somewhat more complex |l, 1 2 ; J, m J ⟩ contributions are found (see table 2). The ratio between the amplitude of state 4 and that of state 4 ′ is small, which indicates weak spin mixing induced by SOC.

Resonance states from neighboring Pb atoms
When the Fe adatom is located above the Pb surface as discussed earlier, there are five neighboring Pb atoms closest to the Fe adatom and equally separated from the adatom, i.e. four Pb atoms at the topmost Pb layer and one Pb atom below the Fe adatom in the subsurface Pb layer. Considering the orientation of the Fe moment along the z axis, we compute the LDOS of these five Pb atoms in the SC state. Some amplitudes at the SC band gap edges must be continued to the continuum DOS outside the SC gap. Firstly, the energies of these resonance states (figure 3(a)) are the same as those of the Shiba states discussed in section 3.3, although their amplitudes are much more reduced than those of the Shiba states. A similar result was obtained for resonances of neighboring boron atoms when a Mn impurity was placed above the surface of MgB 2 [42]. Secondly, at either positive or negative energies of all resonance states, both electron and hole components

Shiba states as a function of Fe moment direction
In the presence of SOC, since the orientation of the Fe moment significantly affects the orbital moment and the normal-state LDOS at the Fermi level, apparent changes are expected in energies and characters of the Shiba states. Figure 4(a) shows how the energies of the Shiba states are varied as a function of θ, where states 1 and 1 ′ are not included since they have mainly s-orbital character. As the Fe moment is rotated from the z axis, three interesting features are observed: (i) the deepest Shiba states are noticeably shifted toward the Fermi level (i.e. the zero energy); (ii) there is a zero-energy peak in the LDOS around θ = 85 • ; (iii) pair 2 at ±0.86∆ for θ = 0 splits into two pairs and one of them is merged into another pair. The first feature indicates an increase of effective exchange coupling between the Fe magnetic moment and the conduction electron spins as θ increases. This may be attributed to both the noticeable increase in the normal-state LDOS at the Fermi level with θ ( figure 1(g)) and the fact that the Pb substrate has strong Rashba SOC which gives in-plane spin texture. In order to have a better understanding of the zero-energy peak, we compute the electron component of the LDOS as a function of θ (see figure 4(b)). At θ = 0 the pair of the deepest Shiba states includes a negative-energy electron component with a large amplitude and a positive-energy electron component with a small amplitude. As θ increases, the latter becomes to pick up a larger amplitude than the former. Around θ = 85 • the negative-energy electron component of the pair merges with its positive-energy electron component, and then the two components are exchanged for larger θ values (85 • < θ ⩽ 90 • ). This process is schematically shown in figures 4(c)-(e). An analysis of the orbitals contributed to the deepest Shiba states confirms the exchange of the two components. Figure 5(a) shows how the energies of the Shiba states are varied as a function of ϕ. Overall, the energies of the Shiba states do not change as dramatically as those for the xz-plane rotation, although merging and splitting of Shiba states appear. This feature is consistent with the subtle changes of the normal-state LDOS at the Fermi level as a function of ϕ ( figure 1(g)). Again we find a zero-energy peak in the LDOS around ϕ = 8 • .
Similarly to the case of the xz rotation, we compute the electron components of the Shiba states as a function of ϕ ( figure 5(b)). From an analysis of them, we find that when the magnetic moment points along the y axis, the electron component of the deepest Shiba states at the negative energy has a negligible amplitude, while the corresponding component at the positive energy has a large amplitude. As discussed earlier in section 3.3, this is a signature of strong suppression of spin-mixing. Note that when the moment points along the z or x axis, such suppression occurs for the Shiba states closer to the SC gap edges rather than near the zero energy.

Effect of additional non-magnetic adatom on Shiba states
In section 3.5 we have shown that the zero-energy peaks appear in the LDOS at some critical directions of the Fe moment vector. It would be interesting to study if the zero-energy peak remains robust over a small non-magnetic impurity close to the Fe adatom. To this end, we consider a case that a foreign non-magnetic adatom is present right next to the Fe adatom on the Pb substrate. A type of the non-magnetic adatom may not be that important for this study. We simulate a sulfur (S) atom located at the nearest neighboring site of the Fe adatom along the x axis on Pb(110). Both the Fe and S adatoms are at the same vertical distance from the Pb surface. When the Fe moment vector is rotated, the S moment vector is also rotated almost antiparallel to the Fe moment. In the normal state, the spin and orbital magnetic moments of the S adatom as well as the changes of the Fe moments are only an order of 0.01% of the spin and orbital moments of the Fe adatom in the S-free case, independently of the moment direction. Therefore, we expect a minimum effect of the S adatom on the Shiba states.
We calculate the Fe LDOS of this system in the SC state as a function of rotation angle θ ( figure 6). This LDOS can be compared to the Fe LDOS in the S-free case (section 3.5), since the contribution of the S LDOS is negligible. For an analysis, we compare the electron components of the Shiba states in this case to those in section 3.5 ( figure 6 vs 4(b)). As expected, overall, the energies of the Shiba states barely change for most θ values. However, a zero-energy peak is now absent in the Fe LDOS at all considered angles. In order to show that this is not due to numerical accuracy, we focus on the electron components of the deepest Shiba states sufficiently close to the zero energy, i.e. for 80 • ⩽ θ ⩽ 90 • . For a zero-energy peak to appear in this range of θ, the amplitudes and the orbital character of the electron components at the negative and positive energies must be reversed relative to a critical angle, as illustrated in figures 4(c)-(e). Such a feature is not observed at any θ values in the range ( figure 6). This suggests that the tiny change made by a non-magnetic atom can be enough for avoiding the Fermi level crossing, and consequently for the zero-energy state to disappear.

Effect of Fe vertical distance on Shiba states
Due to different atomic sizes of Fe and Pb, there is a possibility that the Fe adatom may be located somewhat below the hollow site which we have so far considered. In order to study an effect of Fe vertical distance from the surface, we consider a Fe adatom located at the Pb surface layer, i.e. replacing one of the Pb surface atoms by Fe, which is, henceforth, referred to as Fe-surface site. The position of the Fe adatom discussed in sections 3.1-3.6 is referred to as Fe-above site. In the normal state, when the Fe moment is oriented along the z axis, the spin-up and spin-down Fe LDOS in the Fe-surface case are as well-separated as those in the Fe-above case (see section 3.1). The spin-down LDOS at the Fermi level in the Fe-surface case is about 88% of that in the Fe-above case, which is attributed to slightly more broadened LDOS in the former. The spin and orbital moments of the Fe adatom are now 3.464 and 0.665 µ B . The orbital moment decreases by 0.341 µ B , while the spin moment decreases only by 0.054 µ B , compared to the moments for the Fe-above case. The more broadened LDOS and moment changes can be understood from an increase in the charge transfer from the Pb substrate to the Fe adatom as well as larger hybridization with the Pb substrate than in the Fe-above case. When the Fe moment is oriented along the x or y axis, the orbital moment slightly changes but the spin moment remains almost the same as that for the z axis (see table 1). There are also small induced spin and orbital moments on neighboring Pb atoms. Now let us discuss the Fe LDOS in the SC state when the Fe moment is oriented along the z axis. We identify six pairs of Shiba states at ±0.95, ±0.78, ±0.69, ±0.54, ±0.47, and ±0.26∆, which are labeled as 1-6, respectively, as shown in figure 7(a). Similarly to figures 2(h)-(n) in section 3.3, each pair consists of electron and hole components at positive and negative energies. For example, pair 5 comprises electron and hole components 5 and 5 ′ at ±0.47∆. The Shiba states in the Fe-surface case have some similarities to those in the Fe-above case. Pair 1 at ±0.95∆ mainly arises from the s orbital, which is similar to pair 1 in figure 2(h) in the Fe-above case. In addition, pairs 1, 3, 5, and 6 in the Fe-surface case have similar orbitals and |l, s = 1 2 ; J, m J ⟩ characteristics to pairs 1, 2, 4, and 5 in the Fe-above case, respectively (figure 7(a) vs 2(h)). The Shiba states in the Fe-surface case, however, have several major differences from those in the Fe-above case. Some differences are that there is one more pair of Shiba states and that the deepest pair of Shiba states is much closer to the zero energy, compared to the Fe-above case. Another noticeable difference is that in the Fe-surface case, electron (or hole) components of all the Shiba states do not have significant amplitudes at both positive and negative energies (figures 7(b)-(e)). For instance, in the Fe-surface case, for pair 4, the electron component of the Shiba states has a large amplitude at the positive energy (4 ′ , +0.54∆) and a tiny amplitude at the corresponding negative energy (4, −0.54∆). This implies that SOC-induced spin mixing is highly suppressed in the Fe-surface case. Note that for pairs 3 and 5 in the Fe-above case, the electron (or hole) components of the Shiba states have large or substantial amplitudes at both positive and negative energies (figures 2(h)-(j)), which indicates strong spin mixing induced by SOC . When the orientation of the Fe moment changes from the z axis in the xz plane, the evolution of the energies of the Shiba states as a function of θ is similar to that in section 3.5 ( figure 7(f) vs 4(a)). There is also a zero-energy peak near θ = 86 • below and above which the electron components of the deepest Shiba states at the negative and positive energies are exchanged, as shown in figure 7(h). Regarding the rotation of the Fe moment in the xy plane, the energies of the Shiba states varies as a function of ϕ similarly to the Fe-above case (figure 7(g) vs 5(a)) for small ϕ values, whereas the angular dependence is modified to great extent for intermediate and large θ values. Interestingly, there is a large window of ϕ where a zero-energy state or nearly-zero-energy state appears (see figure 7(i)).
To sum up, for the Fe-surface site, there are larger charge transfer and stronger hybridization with the Pb substrate than the Fe-above site, which entails the noticeable reduction in the Fe orbital moment and the deepest Shiba states closer to the zero energy. When the Fe moment vector points along the z axis, the characteristics of some Shiba states (i.e. the states closer to the zero energy and states near the SC gap edges) remain very similar for the two Fe vertical distances. The energies of the Shiba states in the whole range of θ and at small ϕ values are not very sensitive to the vertical distance of the Fe adatom. A zero-energy state also appears as the orientation of the moment is varied for the Fe-surface site.

Properties of the normal state
When a Co adatom is at the hollow site above the Pb surface, some features of the Co LDOS in the normal state are similar to those of the Fe adatom. Therefore, we discuss mainly differences between the two cases. When the Co moment is oriented along the z axis, the spin and orbital moments of the Co adatom are computed to be 2.183 and 1.215 µ B , respectively. They are much smaller than the nominal values of an isolated Co atom which are 3.0 µ B for both spin and orbital moments. This deviation is much more significant than that in the Fe case, despite the same vertical distance of the Fe and Co adatoms. Therefore, the Co adatom seems to be much more hybridized with the Pb substrate than the Fe adatom.
Since the Fermi level lies almost at the peak of the spin-down LDOS for the Co adatom, the spin-down LDOS at the Fermi level in the Co case is higher than in the Fe case ( figure 8(a) vs 1(b)). Thus the exchange coupling to the conduction electron spins is expected to be stronger than in the Fe case. The peak-to-peak energy difference between the spin-up and spin-down Co LDOS is about 0.078 Ry smaller than that of the Fe LDOS. The Co LDOS shows a somewhat weaker dependence on the moment orientation than the Fe LDOS ( figure 8(b) vs 1(g)). The decomposition of the Co d orbitals at the Fermi level substantially differs from that of the Fe d orbitals with and without SOC (figures 8(c)-(f) vs 1(d)-(f),(i)).

Shiba states without SOC
In the SC state, without SOC, the total Co LDOS shows six pairs of Shiba states at E/∆ = ±0. 15, ±0.18, ±0.22, ±0.51, ±0.72, and ±0.98, as shown in figure 9(a). They originate from {d z 2 ,d x 2 −y 2 }, d yz , d xz , {d x 2 −y 2 ,d z 2 }, d xy , and s orbitals, respectively. The small peaks at ±0.98∆ are marked by arrows in figure 9(a). The energies and characters of the Shiba states qualitatively differ from those of the Fe case, except for the Shiba state of s orbital at ±0.98∆. Comparing the Shiba states with the orbital-decomposed normal-state Co LDOS ( figure 8(f)), we find that the d orbital with the largest LDOS at the Fermi level, d x 2 −y 2 , contributes to the deepest pair of Shiba states and that the d orbital with the smallest LDOS at the Fermi level, d xy , contributes to the farthest pair of Shiba states except for the s-orbital, within the SC gap. Overall, for each d orbital, the order of the Shiba states within the gap is consistent with the normal-state LDOS of each d orbital at the Fermi level. This observation suggests that higher normal-state LDOS at the Fermi level can be considered as an indicator of stronger effective exchange coupling to the conduction electron spin. This is consistent with the prediction by Shiba [2]. Similarly to the Fe case, without SOC, for a given negative or positive energy, either an electron or a hole component of the Shiba state is dominant rather than both components (see figures 9(b)-(e)).

Shiba states with SOC for Co moment along the z axis
When the Co moment vector is oriented along the z axis, five pairs of Shiba states appear at E/∆ = ±0.23, ±0.34, ±0.48, ±0.51, and ±0.95, which are respectively labeled as 5, 4, 3, 2, and 1, as shown in figure 9(f). Table 3    θ > 39 • , the electron component of the deepest Shiba state with a smaller amplitude is exchanged with that with a larger amplitude and the two electron components merge again around 56 • ( figure 10(c)). However, when the Co moment vector rotates in the xy plane, a zero-energy state is not realized ( figure 10(d)). The evolution of the Shiba-state energies as functions of θ and ϕ is distinct from that in the Fe case, which may be due to much stronger hybridization of the Co adatom to the Pb substrate.

Effect of different Co vertical distance on Shiba states
Now let us discuss a case that the Co adatom sits in the Pb surface layer (referred to as Co-surface site), when the moment is oriented along the z axis. Shiba states in this case overall appear deeper in the SC gap than in the case of the Co adatom above the Pb surface (referred to as Co-above site), as shown in figure 10(e) vs (a). This result is consistent with the fact that the lower Co vertical distance, the stronger exchange coupling to the conduction electron spins. Figures 10(e) and (f) shows how the energies of the Shiba states evolve as functions of rotation angles θ and ϕ for the Co-surface case. For the rotation in the xz plane, a zero-energy state now appears around 76 • , while for the rotation in the xy plane, it occurs around 54 • . Note that a zero-energy state is not formed for the rotation in the xy plane in the Co-above case. The angular dependence of the Shiba states in the Co-surface case qualitatively differs from that for the Co-above case (section 4.4) except for the outermost two pairs of the Shiba states. This angular dependence in the Co-surface case is also distinct from that in the Fe-surface case (section 3.7) and that in the Fe-above case (section 3.5).

Properties of the normal state
When the Mn adatom is at the hollow site above the Pb surface, the spin and orbital moments of the Mn adatom are computed to be 4.812 and 0.027 µ B , respectively, which are much closer to those of an isolated Mn atom (5.0 and 0.0 µ B ) than those for the Fe or0 Co case (see table 1). Therefore, we expect that hybridization with the Pb substrate is much weaker than that of the Fe or Co adatom. The Mn adatom has  also well-separated spin-up and spin-down LDOS in the normal state. The peak-to-peak energy difference between the spin-up and spin-down LDOS is about 0.275 Ry, which is significantly larger than that in the Fe or Co case. Since the Fermi level lies at the tail of the Mn spin-down LDOS ( figure 11(a)), the spin-down LDOS at the Fermi level is about one order of magnitude lower than the corresponding value in the Fe or Co case. The LDOS at the Fermi level does not depend much on the orientation of the Mn moment ( figure 11(c)).

Shiba states with SOC
In the SC state, for the Mn adatom above the Pb surface, we observe four pairs of Shiba states at ±0.76∆, ±0.91∆, ±0.96∆, and ±0.98∆ ( figure 12(a)), when the Mn moment is oriented along the z axis. Compared to the Fe or Co case, all of the Shiba states appear very close to the SC gap edges. This is attributed to a very small exchange coupling to the conduction electron spins, which is consistent with the very small normalstate LDOS at the Fermi level, despite the largest spin moment of the Mn adatom. The Shiba states at ±0.98∆ arise from large and small contributions of d x 2 −y 2 and d z 2 orbitals, respectively. This is in contrast to the case of Fe and Co adatoms, where the pair of Shiba states closest to the SC gap edges originates from s orbital. The Shiba states at ±0.96∆ and ±0.90∆ have characters of d z 2 and d xy orbitals, respectively. The fact that the Shiba states of d xy character appear a bit farther from the SC gap edges than those of d z 2 character is consistent with the orbital decomposition of the spin-down LDOS at the Fermi level ( figure 11(b)). Interestingly, in the Mn case, the s-orbital contributions to the Shiba states occur not too close to the SC gap edges such as at ±0.76∆. All of the Shiba states do not noticeably depend on the orientation of the Mn moment in the xz plane ( figure 12(b)) and in xy plane (not shown). Even when the Mn adatom sits in the Pb surface layer, the energies of the Shiba states barely change with the orientation of the Mn moment ( figure 12(c)).
In order to show an importance of large normal-state impurity DOS, we calculate the Mn LDOS in the SC state for the Mn adatom above the Pb surface when the chemical potential in the normal state is raised to the peak of the spin-down LDOS marked in figure 11(a). Figures 12(e) and (f) show Shiba states as functions of θ and ϕ in this case. First of all, now for θ = 0, the deepest Shiba states appear very close to the zero energy, i.e. strong exchange coupling to the conduction electrons. Secondly, the calculated Shiba states strongly depend on the orientation of the Mn moment. Thirdly, for the rotation in the xy plane, there is a zero-energy state around ϕ = 75 • ( figure 12(d)), similarly to the Fe or Co case.

Comparison to experiments
There are no experimental data on Shiba states arising from single Fe, Co, and Mn adatoms on Pb(110). It was reported in [10] [26,[43][44][45] using SC STM tips. (In our work, we studied Fe, Co, and Mn adatoms on Pb(110) with an intention of examining dimers and longer chains on Pb(110) in the future which were experimentally studied.) Different surface orientations have different point group symmetries which play an important role in Shiba states. Therefore, we calculate Shiba states of a Mn adatom on Pb(111) and Pb(001) in the presence of SOC, when the magnetic moment of the Mn adatom points normal to each surface, using the same methods discussed in section 2. Then we compare our calculated Shiba states with the corresponding experimental results, assuming that experimental differential conductance, dI/dV, is proportional to calculated LDOS of the adatoms. In our theory, effects of the SC STM tips are not included. Note that the SC gap of Pb, ∆, is 1.36 meV in our calculations. Therefore, differences between the experimental energies and the calculated energies of the Shiba states are of an order of 0.1 meV.

Mn adatom on Pb(111)
We compute the Shiba states of the Mn adatom on a Pb(111) surface. Our calculation shows three pairs of Shiba states within the gap ( figure 13(a)). The Shiba peaks near the gap have d z 2 -orbital character, while the deepest Shiba peaks have mainly s-orbital character with a small contribution from d z 2 orbital. The other Shiba peaks appearing at ±0.94∆ have d xz and d yz orbital characters. The orbital characterization is consistent with the point group symmetry of the (111) surface, C 3v .
Experimental STM data [45] showed two pairs of Shiba states labeled as l = 0 and l = 1 peaks (red bars in figure 13(b)). Another experimental STM data [26,43] showed three pairs of Shiba states identified as d xz,yz , d x 2 −y 2 ,xy , and d z 2 orbital contributions (blue bars in figure 13(b)). The Shiba peaks from the two experimental data have some similarities and some differences. The Shiba peaks in [45] are shown at ±0.38(±0.02) and ±1.02(±0.02) meV, while the peaks in [26] appear at ±0.22, ±0.77, and ±1.18 meV. For easier comparison to our results, the experimental energies of the Shiba states in figure 13(b) are normalized by the SC gap, 1.52 (1.35) meV, as reported in [45] ( [26,43]).
First of all, some calculated Shiba peaks with d xz,yz and {s, d z 2 } orbital characters ( figure 13(a)) are also found in the experimental data [26,43,45]. In our calculation, however, the Shiba states of d x 2 −y 2 and d xy character found in one of the experiments [43] do not appear within the SC gap. The energies of the deepest Shiba pair are quite different between our theory and the two experimental data. The most noticeable discrepancy between our theory and the experiments is the asymmetry in the experimental Shiba peak heights relative to the Fermi level (i.e. zero). It was proposed that the asymmetry in the spectral density occurs when on-site Coulomb interaction of the adatom (or non-magnetic potential scattering by the adatom) is greater than the SC gap [4]. For a normal STM tip, it was theoretically shown [46] that the dI/dV peak heights become asymmetric when the magnetic impurity interacts with other impurities or a normal reservoir. For a SC STM tip, a theory for the asymmetry was proposed in [26].  [44]. For easier comparison between our theory and the experimental data, in (b) and (d), only dI/dV peak heights vs sample bias deconvoluted from the SC STM tip are shown, where ∆exp is the reported experimental SC gap. It is 1.52 (1.35) meV for the red (blue) bars in (b) and is 1.25 meV for (d). In [45], normalized dI/dV was reported, while in references [26,43,44], dI/dV was in units of µS. The orbital characters are given for the Shiba peaks in (a)-(d).

Mn adatom on Pb(001)
Now we discuss the case of a Mn adatom on Pb(001) in the SC state. Our calculation shows four pairs of Shiba states within the gap ( figure 13(c)). The Shiba peaks near the gap have contributions from three orbitals such as d x 2 −y 2 , d xy , and d z 2 . The (001) surface has C 4v symmetry. The orbitals belonging to different irreducible representations are hybridized due to strong SOC. Two split pairs of Shiba states of d xz and d yz hybridized orbitals are found. The deepest pair of Shiba states have mainly s orbital character with a small d z 2 -orbital contribution.
The experimental energies of the Shiba states are ±0.31, ±0.82, and ±1.21 meV that are normalized by the reported SC gap, 1.25 meV [43,44]. Two of the four Shiba pairs such as d x 2 −y 2 -orbital and s-orbital characters are also found in the experimental data (figure 13(d) [43,44]) even with almost the same energies as the calculated values, if we allow a possibility of hybridization between d z 2 and s orbitals. Interestingly, for the Pb(001) surface, the experimental Shiba peak heights (dI/dV values) are symmetric relative to the Fermi level just like our calculated LDOS values. Again, the experimental deepest Shiba states appear much closer to the Fermi level than the calculated states.

Conclusions
We studied Shiba states of Fe, Co, and Mn adatoms on Pb(110) and a Mn adatom on Pb(111) and Pb(001) by using fully relativistic first-principles calculations including the band structure of the Pb substrate and all orbitals of the adatoms in the normal and SC states. In the normal state, the LDOS at the Fermi level is the highest (lowest) for the Co (Mn) adatom, while the spin moment of the adatom has an opposite trend.
For the Fe or Co adatom, five or six pairs of the Shiba states were spread out within the SC gap. Due to strong SOC of the Pb surface, the Shiba states cannot be characterized by separate orbitals categorized by irreducible representations of the point group but by J and m J . Another hallmark of the SOC effect is that the Shiba states greatly change when the orientation of the adatom magnetic moment varies. By using a spin-polarized STM tip, the orientation of the adatom moment may be rotated within a surface plane.
Additionally, the deepest pair of the Shiba states can be merged to form a zero-energy state with the rotation of the adatom moment, whether the Fe or Co adatom is located above the surface or at the surface. Interestingly, this zero-energy state is found to be so fragile that it can be easily destroyed by an additional non-magnetic adatom right next to the magnetic adatom, although the other Shiba states remain unchanged. This kind of fragility test may be performed in STM experiments. This zero-energy state seems to be similar to the zero-bias peaks observed in Fe impurities on SC Ta(100)-(3×3)O [17], in magnetic molecules on a SC Pb substrate [47,48], in impurities on SC vanadium [49], and in Gd dimers/trimers in Bi-coated SC Nb [50], although their origins differ from our case.
For the Mn adatom on Pb(110), we show that all Shiba states are observed quite far from the zero energy, independently of the orientation of the moment. This is attributed to the tiny LDOS at the Fermi level, despite the largest Mn spin moment among the three adatoms. Therefore, we do not expect formation of a zero-energy state from the rotation of the Mn moment unless the LDOS at the Fermi level becomes significant.
Compared to the experimental data on Mn adatoms [26,[43][44][45], overall our Shiba states are in good agreement with experimental Shiba states. Our Shiba states show a significant dependence on the surface orientation, similarly to the experiments. Some calculated Shiba states are also found in the experimental data with the similar orbital characteristics and/or similar energies. There are, however, two small discrepancies between our theory and the experimental data for both Pb(111) and Pb(001). Firstly, the experimental deepest Shiba states [26,[43][44][45] appear much closer to the Fermi level than the calculated Shiba states. We cautiously speculate that this may be due to a large charge transfer between the adatom and the substrate in the experimental samples, which could increase the normal-state LDOS value at the Fermi level, similarly to figures 12(d)-(f). Secondly, the Shiba states with strong s-orbital character have not been identified in references [26,43,44], although it was reported for [45].
In order to further compare with experiments, we propose the following first-principles calculations/developments. Firstly, it would be interesting to compute a spatial distribution of the Shiba states along different k directions and compare it to experiments, similarly to references [51][52][53][54]. The experimental dI/dV intensity was shown to oscillate with a periodicity of 5.8 Å, for the Shiba states from the Mn adatom on Pb(111) [43]. Secondly, it would be useful to compute the Shiba states with either varying the vertical distances continuously, or varying the charge state of the adatom. Experimental data from different groups show somewhat different Shiba states even for the same Mn adatoms on the same Pb surface orientation [26,43,45]. As [43] pointed out, this is likely due to different vertical heights of the adatoms from the surface. The vertical height from the surface determines an amount of hybridization between the adatom and the substrate states. In addition, it can affect the charge state of the adatom via charge transfer between the adatom and the substrate [55,56]. Thirdly, it would be interesting to compute the Shiba states for dimers (with different distances between them) and longer chains on Pb surfaces with noncollinear magnetism and study effects of canting and noncollinear magnetism on the Shiba states in long magnetic atomic chains.

Data availability statement
The data cannot be made publicly available upon publication because they are not available in a format that is sufficiently accessible or reusable by other researchers. The data that support the findings of this study are available upon reasonable request from the authors.