Observation of Coexistence of Yu-Shiba-Rusinov States and Spin-Flip Excitations

We investigate the spectral evolution in different metal phthalocyanine molecules on NbSe2 surface using scanning tunnelling microscopy (STM) as a function of the coupling with the substrate. For manganese phthalocyanine (MnPc), we demonstrate a smooth spectral crossover from Yu-Shiba-Rusinov (YSR) bound states to spin-flip excitations. This has not been observed previously and it is in contrast to simple theoretical expectations. We corroborate the experimental findings using numerical renormalization group calculations. Our results provide fundamental new insight on the behavior of atomic scale magnetic/SC hybrid systems, which is important, for example, for engineered topological superconductors and spin logic devices.


Geometries and spin densities from DFT calculations
Density functional theory calculations were performed with the FHI-AIMS computational package 1,2 and the PBE generalized gradient approximation for the exchange-correlation functional. 3 We used the standard "light" numerical settings and basis sets of numeric atomic-centered orbitals tested and recommended by FHI-AIMS. Periodic NbSe 2 supercells were sampled with a 2 × 2 k-point grid centred on the Γ point. Van der Waals interactions were included by the post-SCF Tkatchenko-Scheffler correction. 4 Before computing the electronic structure, all atomic forces were relaxed to < 0.01 eV/Å. a b c Figure S1: (a,b,c) Results for CoPc (a), FePc (b), and MnPc (c) computed with density functional theory (DFT). 1,2 The total computed spins for the systems are S = 0.41, S = 1.22, and S = 1.72, respectively. The localized metal atom-projected (Co, Fe, or Mn) spins are S = 0.49, S = 1.10, and S = 1.70, respectively. The localized atom-projected spins differ by < 0.1 compared to their gas phase counterparts with the same DFT parameters and numerical settings.

Experiments on FePc in external magnetic field
We can verify that the observed transitions in FePc do indeed correspond to inelastic spin excitations by acquiring dI/dV spectra at different magnetic (B) field strengths aligned perpendicular to the sample surface. The energies of the first and second feature change with the external field, demonstrating that these are spin excitations (see Fig. S2). The zero-field splittings (ZFS) are well described by the following spin Hamiltonian: where D is the axial ZFS constant to determine the magnetic anisotropy, E is the transverse anisotropy, and S z , S x , S y are the spin operators. The selection rule for the spin excitations implies that they can only occur between the states differing by m S =±1. For S = 1 spinstate, we would expect only one step in the absence of transverse anisotropy (E = 0) and at B = 0T. Transverse anisotropy (E > 0) mixes the states and therefore a second step can occur. Hence, the dI/dV spectra on FePc on NbSe 2 suggest spin-state S = 1 with transverse anisotropy. We fitted our data using a phenomenological spin Hamiltonian 5-8 where g is the Landgé g factor, µ B the Bohr magneton, and S γ is the component of spin

Experiments on MnPc in external magnetic field
It is conceivable that there would be inelastic excitations in the energy corresponding to the features seen outside the superconducting gap on MnPc that correspond to molecular vibrations (phonons). To support our assignment of this feature to a magnetic excitation,

Numerical renormalization group calculations
The impurity model is solved using the numerical renormalization group (NRG) method. 9,10 The NRG is a numerical procedure for solving quantum impurity models that is based on a logarithmic discretization of the continuum of electrons (here itinerant electrons from the substrate hybridized with the molecular states). The discretization is controlled by the parameter Λ > 1 that determines the coarseness of the frequency grid, ω n ∼ Λ −n . The The energies of the sub-gap states can be directly extracted from the renormalizationgroup flow diagrams in the NRG. This approach is significantly more accurate than determining the peak positions in the spectra, and does not suffer from any broadening artifacts.
Example results are shown in Fig. S5 The  Figure S5: Energies of the many-body states, referenced to the ground-state energy of the system. This diagram shows the energies of all states below the edge of the continuum excitations. The lowest state is hence by definition the ground state, while the excited states give rise to the YSR peaks in the spectra. Black lines correspond to the states resulting from the S = 3/2 multiplet after anisotropy splitting (i.e., the S z = ±1/2 pair). Red lines correspond to the states resulting from the impurity state "screened" by binding one Bogoliubov quasiparticle from the continuum. These states originate from the S = 1 multiplet and are separated into a S z = 0 (lower energy) and a S z = ±1 (higher energy) subsets. The S z = 0 YSR state is present for all J, the S z = ±1 pair starts emerging from the continuum at ρJ ≈ 0.08. At ρJ ≈ 0.12, a quantum phase transition between the S z = ±1/2 and S z = 0 many-body states occurs. equilibrium spectral function, this procedure is somewhat ill-defined, and furthermore suffers from possible broadening artifacts. In theoretical calculation using the NRG one can, however, use a more robust and direct approach for very accurately extracting the energies of spin excitations. It consists of computing the transverse part of the dynamical spin susceptibility function The spin excitation energies can then be directly read off from the peak position, which is unique and very well defined. This is illustrated in Fig. S6 S9  Figure S7: Spectral functions for a range of exchange coupling strengths. These plots correspond to horizontal line-cuts of the density plot in Fig. 4a of the main text. They show more clearly how the YSR state detaches from the band edge of the continuum of Bogoliubov excitations, and how the spin-flip line-shape evolves (in particular for ρJ ≥ 0.06, i.e., after the emergence of YSR as a well-defined sub-gap state). At ρJ ∼ 0.1, a new subgap feature starts to detach from the band edge: this is an additional YSR state (of type |S z = 1 ± |S z = −1 ) that results from the magnetic anisotropy splitting of the S = 1 YSR multiplet. This parameter range is, however, not experimentally relevant.
The insensitivity of the YSR state energy on D 0 is illustrated in Fig. S8, which enables us to calibrate the exchange coupling based on the experimental YSR energies.
Figure S10: Calculated impurity spectral function for normal-state and superconducting substrate. We compare the low-energy parts of the spectral functions computed for an anisotropic Kondo impurity model (S = 3/2, longitudinal anisotropy D 0 ) in the limit of very small exchange coupling J, so that the Kondo temperature is negligibly small. For normal-state substrate, the spectrum shows steps due to inelastic excitations for bias voltage beyond the spin-flip excitation threshold ω sf = 2D 0 . The significant width of the step is a broadening artifact of the numerical method (lower broadening parameter would lead to stronger artifacts, which are already visible in these results in the form of weak oscillatory features). For superconducting substrate, a gap is formed, while the spin-flip features are shifted to ω sf = ∆+2D 0 . Furthermore, the spectral shape of the inelastic excitations inherits the form of the density of states of the superconductor near the threshold of the band of Bogoliubov excitations.

Extraction of inelastic-excitation step heights
Extracting the amplitude of the inelastic excitations in the case of a normal-state substrate is simple: extract the step amplitude and reference it to the asymptotic value. For a superconducting substrate, the inelastic excitation corresponds to a spectral feature with a complicated shape that also changes with the increasing value of the exchange coupling.
In order to compare these results with the experiments, we have used two different ap-  Properties of the anisotropic S = 1 Kondo model Figure S12: Theoretical spectral functions for a S = 1 impurity with varying exchange coupling J. Spectral function for the anisotropic S = 1 Kondo impurity with magnetic anisotropy parameters D = 5.5 meV and E = 1.5 meV (spin-flip excitation energies D − E = 4 meV and D + E = 7 meV). At low J, the superconducting band edges are hardly visible, but their amplitude increases with increasing J. At still higher J (not shown here), a S = 1/2 sub-gap YSR state detaches from the gap edge.