Full orbital decomposition of Yu-Shiba-Rusinov states based on first principles

Magnetism and conventional superconductivity are phenomena related to competing types of order. Theoretically, it was first shown by Abrikosov and Gor’kov Abrikosov and Gor’ Kov (1959) that at a paramagnetic impurity concentration of $\sim 1\%$ in a superconductor the spectral energy gap would no longer correspond to the ordering parameter $\Delta$. For higher concentrations of impurities the spectral gap can even vanish totally, leading to gapless superconductivity characterised by a finite T_C and $\Delta$ but with no spectral gap.Phillips (1963) Local, real space models were subsequently constructed by Yu Yu (1965), Shiba Shiba (1968) and Rusinov Rusinov (1969) using local, one band models around a classical impurity spin. They predicted the existence of a pair of localised, in-gap (YSR) states either side of the Fermi energy associated with the exchange splitting $J$ of the spin.

Subsequently, experiments investigating YSR states have predicted multiple pairs of in-gap states Ji et al. (2008); Hatter et al. (2015); Ruby et al. (2016); Choi et al. (2017); Liebhaber et al. (2020), where it was argued that the origin of multiple resonances may arise from magnetic anisotropy Hatter et al. (2015), the orbital character Ruby et al. (2016); Choi et al. (2017) or modulations in the charge density Senkpiel et al. (2019); Liebhaber et al. (2020). Ruby et al Ruby et al. (2016) investigated the (001) surface of Pb with a Mn impurity adsorbed onto the surface. They argue that the multiple YSR resonances originate from the crystal field splitting of the Mn d-orbitals. Using energy considerations and real space $dI/dV$ maps they were able to assign the relevant orbitals to the YSR resonances.

When investigating more complex superconductors like NbSe₂ Senkpiel et al. (2019); Liebhaber et al. (2020) or $\beta-$Bi₂Pd Choi et al. (2018) it becomes immediately obvious that disentangling hybridised YSR peaks, or YSR peaks entangled with coherence peaks will become increasing challenging. Symmetry arguments and energy consideration will not be sufficient to uniquely assign the large number of in-gap resonances. In this report we build upon previous work Saunderson et al. (2020a, b) to address this problem from first principles, considering all electrons and their magnetic orderings fully while capturing superconductivity within a one-parameter model of local BCS type pairing.

We use this formalism to investigate the 3d series of elements as impurities in fcc Pb. After a brief introduction of the methodological implementations we show our results for all magnetic 3d impurities discussing the distinct YSR resonance pairs arising from the anticipated t_2g and e_g orbitals. However, beyond that we will highlight the existence of $l=0$ (s-electron) YSR resonances enforcing the necessity of an all-electron description. Finally, we will analyse the spatial decay of the magnetism as well as the in-gap states within the superconducting Pb, presenting the orbital-resolved local densities.

The implementation of the Bogoliubov-de Gennes (BdG) equation in the Korringa-Kohn-Rostoker Greens function method Saunderson et al. (2020a) has been described earlier as well as the extension to real space impurity systems Saunderson et al. (2020b). Here, we present a further step namely the incorporation of collinear magnetism and in the following we will restrict the technical discussion to the equations relevant for this development.

In order to incorporate magnetism and superconductivity Oliveira et al. (1988) into density functional theory Dreizler and Gross (1990), three effective potentials are required, the electron potential $V_{eff}(\mathbf{r})$, the magnetic field $B_{eff}(\mathbf{r})$ and the effective pairing potential $\Delta_{eff}(\mathbf{r})$,

Here, $E_{xc}[\rho,m]$ is the exchange correlation functional for the normal state, $\rho(\mathbf{r})$, $m(\mathbf{r})$ are the usual charge and spin densities and $\chi(\mathbf{r})$ is the anomalous density. Finally, $\Lambda$ is the interaction parameter Saunderson et al. (2020a), which is the one free parameter in our description typically fixed to recover experimentally observed gap sizes. This framework, using a simplified phenomenological parameter, was introduced in Ref.  Suvasani et al. (1993) and subsequently implementation were presented in Refs. Csire et al. (2015, 2018a). It was already been shown to effectively describe gap anisotropy Saunderson et al. (2020a) and impurity scattering Saunderson et al. (2020b), along with complex superconducting order parameters in LaNiC2 Csire et al. (2018b) and LaNiGa2 Ghosh et al. (2020) Within the non-relativistic theory the densities are given by

Hence, the resulting spin BdG Hamiltonian is defined as,

where

Equation (7) can be brought into a block-diagonal form such that

where $\sigma=\{\uparrow,\downarrow\}$ and $\sigma^{*}$ represents the opposing spin to $\sigma$. The corresponding Green’s function is defined as

where $\hat{H}_{BdG}(\mathbf{r})=\langle\mathbf{r}|\hat{H}_{BdG}|\mathbf{r}\rangle$, and can be simplified into a block diagonal form accordingly

The relevant densities are expressed by the Green’s function

Implementing this extension in the corresponding bulk Saunderson et al. (2020a) and real-space impurity code Saunderson et al. (2020b) will enable us to address the coupling between the superconducting state and magnetism. Here, we will focus on the effect of magnetic impurities with the corresponding in-gap YSR states. As test scenario we consider superconducting bulk fcc Pb Saunderson et al. (2020a) with the interaction parameter $\Lambda$ tuned such to the experimental gap size from Ruby et al. Ruby et al. (2015). For the real-space impurity cluster embedded in the infinite periodic crystals we used 87 atoms as shown in panel Fig. 1 (a).

Starting with self-consistent scalar relativistic (SRA) solutions for the normal state, the resulting magnetic moments of the localised 3d impurity atoms are summarized in Fig. 2. For the Pb host we identify only V, Cr, Mn and Fe to be magnetic. For all other 3d elements only the non-magnetic solutions were established self-consistently. The resulting local density of states (LDOS) for all magnetic impurities is shown in Fig. 3. Notably, while for V it is the majority spin channel which is close to the Fermi energy, this impurity level shifts further away from the Fermi energy as we go through the 3d series. For the Fe impurity it is the minority spin channel which is situated right at the Fermi energy. For the elements in between we see a gradual transition between the two cases.

In the cubic fcc lattice the crystal field lifts the degeneracy of the d-orbitals and the impurity levels are split into the e_g and t_2g orbitals. Dresselhaus et al. (2007) As an example this is shown for the Fe minority level as inset in Fig. 3.

Within simplified models Yu (1965); Shiba (1968); Rusinov (1969); Balatsky et al. (2006) it has been shown that it is the energy splitting between the minority and majority levels which will determine the energy positioning of the in-gap states. For these models the impurity spin is considered as classical magnetic moment $S$ coupled to the electrons of the crystal via the exchange interaction $J$. The resulting energies are than expressed in terms of a T-matrix scattering approach. As the e_g and t_2g experience slightly different splitting this will result in a clear lifting of the degeneracy between those orbitals for the in-gap states. We have summarized the corresponding exchange energies

in Fig. 4. In the above definition $p^{\uparrow}$ and $p^{\downarrow}$ are the spin-up and spin-down peak positions, respectively. As expected, the trend follows the total magnetic moment being largest for Cr and Mn. Importantly, for V we find the strongest difference between the splitting in the e_g and t_2g states which would suggest that the corresponding in-gap states will be well separated in energy. In contrast, this is weakest for Fe as indicated by the smallest separation of the in gap states.

Extending the analysis including superconductivity gives access to the YSR in-gap states. For each step we perform the corresponding normal state calculations first and extend it to the superconducting case in the following. In all cases we perform fully self-consistent SRA calculations based on the BdG Hamiltonian Eq. 7 to finally generate the electronic and magnetic structure for the magnetic atoms embedded in superconducting Pb. In all cases we set the interaction parameter $\Lambda$ at the impurity site to zero and keep it at the bulk value for all other sites. Saunderson et al. (2020a, b) The results are summarized in Fig. 5 for all magnetic impurities. In all four cases we find pronounced in gap states, well separated into e_g and t_2g states. As expected from the previous analysis the splitting between the e_g and t_2g is largest for V and hard to resolve for Fe. Table 1 summarizes the associated energies of the bound states. A clear correlation between Table 1 and the normal state exchange energies, summarized in Fig. 4, is evident.

Furthermore, the superconducting symmetry implies that for each in-gap state there is a minority and a majority state symmetrically placed relative to the Fermi energy, which we indicated via the $\pm$ in Table 1. While the energetic positions are forced to be symmetric the height of the corresponding states is determined by the normal state LDOS at the Fermi energy. As the minority state of the Fe impurity is perfectly placed at the Fermi energy (see Fig. 3) this leads to the largest in-gap peak for the minority in-gap state (see Fig. 6). Similarly for V the majority level is closest to the Fermi energy leading to the corresponding in-gap state to be significantly larger than the minority peak. In contrast, for Cr both levels have similar distance to the Fermi energy in the normal state resulting in a very similar height for the Cr in-gap states.

In order to make the connection between the normal state impurity levels and the superconducting in-gap states more quantitative we follow previous models established for the YSR states Yu (1965); Shiba (1968); Rusinov (1969); Balatsky et al. (2006). For isotropic (l=0) scattering the energy can be approximated to

where

Here, $\Delta_{0}$ is the bulk quasiparticle gap, $N_{0}$ is the density of states of Pb at the Fermi level in the normal state, and V is the non-magnetic scattering potential. Ignoring the correction from the scalar potential, $\beta=0$, but generalizing to the case where $N_{0}$ is different for e_g and t_2g states, the equation reduces to

where $a$ stands for either e_g and t_2g states and $b$ is the index for the impurity V, Cr, Mn, Fe. Correspondingly, $\alpha^{b}_{a}$ generalises to

As the energies depend crucially on $N^{a}_{0}$ we decided to determine this parameter by fitting Eq. 22 to the full ab initio calculations in the case of Mn. Here, we fixed $(JS)^{Mn}_{a}$ to the results shown in Fig. 4. The resulting values are $N^{e_{g}}_{0}=0.966~{}({\mathrm{e}V})^{-1}$ and $N^{t_{2g}}_{0}=0.593~{}({\mathrm{e}V})^{-1}$. In the following we kept $N^{a}_{0}$ the same for all impurities as it should capture the properties of the nonmagnetic host Pb and calculated the in-gap states using $(JS)^{b}_{a}$ for all the other impurities. The results in comparison to the directly extracted energies are shown in Fig. 6. While the agreement is far from perfect the trends are correctly reproduced. However, especially for V and Fe, the elements with the largest and smallest splitting between the t_2g and e_g states, respectively, the model fails to reproduce the quantitative results of the full calculations. This highlights the importance of the full ab initio description as the model fails to capture the quantitative details making direct comparisons to experimental findings difficult.

A further conventionally made approximation is the restriction to $d$-orbitals only. Ruby et al. (2016); Schrieffer (1967) Given that the magnetism is dominated by the $d$ electrons and it is the magnetism which induces the in-gap states it appears natural to follow this approximation. However, in any real material all orbitals will hybridize and magnetism arising in the $l=2$ orbitals will ultimately induce spin-polarization in all other orbitals as well. Within our calculations this is particularly visible for Cr and Mn for which additional YSR peaks are visible near the Pb coherence peak (see Fig. 5) at energies just above $(\epsilon-\epsilon_{F})\sim 1.3meV$. In this particular case it turns out they are $l=0$ orbital contributions. This finding highlights the fact that while in a first approximation it appears natural to reduce the discussion to the d orbitals, in details especially near the edge of the superconducting gap other orbitals might play a dominant role. Furthermore, determining which impurity has a large response from other orbitals will not be easy to predict without all-electron calculations.

Finally, we would like to analyse the spatial dependence of the in gap states. In Fig. 7 the radial decay of the magnetic moment is shown for the case of a V impurity. It drops quickly to almost zero in the first shell already and oscillates weakly up to the 6th shell. For the total magnetic moment, there is no visible difference between the normal and the superconducting state in the spatial decay. In order to visualize the behaviour of the in-gap states we present the atom-resolved charge densities at the energy associated with the in-gap state in Fig. 8. In the case of the V impurities these energies are for the e_g states $\epsilon=0.95$meV and for the t_2g states $\epsilon=0.35$meV. As discussed earlier, the cubic lattice leads to a splitting into e_g and t_2g states but does not lift the degeneracy of the $d_{z^{2}}$ and $d_{x^{2}-y^{2}}$ orbitals within the e_g-level nor for the $d_{yz}$, $d_{xz}$ and $d_{xy}$ orbitals within the t_2g level. For this reason any visualisation for the two levels would preserve the cubic symmetry. Resolving all orbitals separately (see Fig. 8) highlights the power of full ab initio calculations. While the spatial resolution is limited to atomic sites the orbital characters are nevertheless clearly visible. The $d_{yz}$ and $d_{xz}$ reduce to a two-fold rotational symmetry around the z-axis rotated by $90$ degrees relative to each other. The other three orbitals show the corresponding four-fold rotational symmetry around the z-axis. Furthermore, the larger spatial extension of the $d_{z^{2}}$ orbital in the $z$ direction is clearly visible. Finally, as the the $d_{xy}$ and $d_{x^{2}-y^{2}}$ orbitals are rotated by $45$ degrees relative to each other this results in the $d_{x^{2}-y^{2}}$ pointing along the nearest neighbour bonds in the $z$ plane of the impurity atom. In contrast, the $d_{xy}$ orbitals show the largest contribution out of plane as there the nearest neighbours are in direction of the orbital lobes. At a surface the degeneracy of these orbitals would be lifted and the orbital induced angular dependence of the density can be resolved with STM experiments. Ruby et al. (2016); Choi et al. (2017); Ruby et al. (2018); Küster et al. (2021)

We have extended previous work Saunderson et al. (2020a, b) of implementing the BdG equations into the KKR method with substitutional impurities and collinear magnetism. As a model system we consider 3d impurities in a Pb three dimensional crystal inspired by experiments performed by Ruby et al Ruby et al. (2016) investigating magnetism induced in-gap states at a superconducting surface. As predicted by simple models, the position and height of the induced in-gap states is strongly related to the normal state exchange splitting as observed for the magnetic impurities. However, while the over all correlation is clearly visible in the details the simplified model fails to make quantitatively correct predictions for the varying 3d impurity atoms. The most significant limitations of the model are the restriction to isotropic scattering and ignoring any hybridization of the d-electrons to other orbitals. Naturally, these models can be extended to more complex scenarios however rendering them more intractable as the parameter space increases. Within our real space superconducting DFT descriptions we are able to capture all orbital induced variations in their full complexity.

Beyond the consideration of the d-orbitals we established clear signatures of $l=0$ in-gap states induced via s-d hybridization. While these states are in close vicinity to the superconducting coherence peaks and as such difficult to observe in experiment, their existence highlights the complexity of any quantitative interpretation of experimental observation. Within our calculations these $l=0$ in-gap states were particular pronounced for Cr and Mn but hardly visible for V and Fe. This again points to the importance of a correct all-electron description of the underlying band structure. This is especially relevant as similar YSR resonances have already been investigated Senkpiel et al. (2019).

Finally, we investigate the radial decay of the magnetism and in-gap states inside the superconducting crystal. Interestingly, the decay of the magnetic moment is largely unaffected by the superconducting state. This is related to the fact that even in the normal state the magnetic moment decays quickly for the nearest neighbour atoms already. Nevertheless, by investigating the decay of the in-gap states for each orbital separately it was clearly possible to resolve the distinct orbital symmetries. This final step will enable us to make direct contact to STM experiments Ruby et al. (2016) where the breaking of further spatial symmetries at the surface of the material will lift the degeneracies among the e_g and t_2g states.

In summary, we have performed fully self consistent calculations for magnetic impurities within the superconducting state using the BdG equations implemented within the KKR formalism. We have discussed qualitatively and quantitatively the formation of YSR resonances associated with the magnetic moment from the e_g and t_2g orbitals according to the cubic symmetry. Furthermore, we have established the existence of $l=0$ YSR resonances, highlighting the need for an all-electron description for even a qualitatively correct description of real materials. This is further strengthened by the fact that the position of the in-gap states is at best described qualitatively with simple models even for the dominant e_g and t_2g resonances. With the introduction of more impurities we could investigate the hybridisation of the YSR resonances resulting in the formation of YSR bands. Zhang et al. (2020); Ruby et al. (2017) Similarly, the use of the KKR framework will enable us to extend the method to concentrated alloys exploiting the coherent potential approximation (CPA) as outlined previously. Martin et al. (1999); Moradian et al. (2000a, b) Finally, the incorporation of the fully relativistic BdG equations Csire et al. (2018a) including spin-orbit coupling will, in a next step, enable us to investigate the existence of Majorana zero modes Nadj-Perge et al. (2013); Röntynen and Ojanen (2015); Kim et al. (2018); Schneider et al. (2020) at the surface of conventional superconductors.

This work was carried out using the computational facilities of the Advanced Computing Research Centre, University of Bristol - http://www.bris.ac.uk/acrc/. The above work was supported by the Centre for Doctoral Training in Condensed Matter Physics, funded by EPSRC Grant No. EP/L015544/1. G.C. gratefully acknowledges support from the European Union’s Horizon 2020 research and innovation program under the Marie Sklodowska-Curie Grant Agreement No. 754510. This work was supported by Spanish MINECO (the Severo Ochoa Centers of Excellence Program under Grant No. SEV- 2017-0706), Spanish MICIU, AEI and EU FEDER (Grant No. PGC2018-096955-B-C43), and Generalitat de Catalunya (Grant No. 2017SGR1506 and the CERCA Program). The work was also supported by the European Union MaX Center of Excellence (EU-H2020 Grant No. 824143). M.G. thanks the visiting professorship program of the Centre for Dynamics and Topology ad Johannes Gutenberg-University Mainz. The authors would like to thank M.-H. Wu and R. Gupta for many helpful discussions.