Coupled Yu-Shiba-Rusinov states in molecular dimers on NbSe$_2$

Magnetic impurities have a dramatic effect on superconductivity by breaking the time-reversal symmetry and inducing so-called Yu-Shiba-Rusinov (YSR) low energy bound states within the superconducting gap. The spatial extent of YSR states is greatly enhanced in 2D systems, which should facilitate the formation of coupled states. Here, we observe YSR states on single cobalt phthalocyanine (CoPC) molecules on a 2D superconductor NbSe$_2$ using low-temperature scanning tunneling microscopy (STM) and spectroscopy (STS). We use STM lateral manipulation to create controlled CoPc dimers and demonstrate the formation of coupled YSR states. The experimental results are corroborated by theoretical analysis of the coupled states in lattice and continuum models. Our work forms an important step towards the realization of exotic topological states in designer magnetic lattices.

Ferromagnetic coupling between the two impurities results in YSR states that hybridize and split into bonding and anti-bonding states [6,[15][16][17].Atomic chains with a suitable spintexture have been suggested to support 1D topological superconductivity [18] and Majorana modes at each end [19], which have been recently realized in iron chains on Pb(110) [12,13].
In general, controlled coupling of YSR states should enable realizing designer quantum materials with novel topological phases [20].However, the experimentally observed YSR states on s-wave superconductors (Pb and Nb) are localized within a few atomic distances from the impurity centre and their energies vary strongly depending on the adsorption site of the impurity.This greatly hampers forming assemblies of controllably coupled YSR states.
Recently, Ménard et al. demonstrated that the spatial extent of the YSR states depends on the dimensionality of the superconductor and is greatly enhanced in 2D systems, where the impurity bound state can extend over several nanometres away from the impurity [21].
This was demonstrated on iron impurities embedded into a niobium diselenide (NbSe 2 ) substrate, which is layered material with a superconducting transition at 7.2 K. Due to the weak van der Waals interaction between adjacent layers, it behaves essentially as a 2D system.However, coupling between YSR states on 2D superconductors has not yet been demonstrated.Here, we observe YSR states on single cobalt phthalocyanine (CoPC) molecules on a 2D superconductor NbSe 2 using low-temperature scanning tunneling microscopy (STM) and spectroscopy (STS).We use STM lateral manipulation to create controlled CoPc dimers and demonstrate the formation of coupled YSR states.The experimental results are corroborated by theoretical analysis of the coupled states in lattice and continuum models.Our work forms an important step towards the realization of exotic topological states in designer magnetic lattices [20,22]. of the outermost Se atoms, and the well-known 3 × 3 charge-density wave (CDW) superstructure [21,23].Figure 1b displays a topographic STM image of two isolated cobalt phthalocyanine (CoPc) molecules on NbSe 2 , where the principal directions of the substrate have been marked with arrows.We have determined the adsorption site of CoPc from atomically resolved images (Supplementary Information), which is relevant as the molecular properties (e.g.spin state and more importantly, the magnetic coupling to the superconductor) might depend on this.Our STM measurements indicate that independent of the molecular orientation, the central Co atom always sits on top of an underlying Se atom.This experimentally observed adsorption configuration is confirmed by density-functional theory (DFT) calculations, with Fig. 1c showing the most stable adsorption structure.dI/dV spectroscopy was performed with NbSe 2 -coated SC tip to increase the energy resolution beyond the thermal limit [6,7,23].Figure 1d (upper panel) shows a typical dI/dV spectrum acquired on bare NbSe 2 substrate.It exhibits a gap that is twice the size of the SC gap of NbSe 2 , with sharp coherence peaks at ±2∆ ∼ ±2 meV [24,25].The small peak close to V = 0 is due to the tunneling of thermally excited quasiparticles.dI/dV curve taken on a CoPC molecule (Fig. 1d, lower panel) shows that the coherence peaks at ±2∆ are replaced by pronounced peaks of asymmetric heights within the SC gap and two replicas of those peaks close to V = 0, similarly to an earlier report on MnPC on Pb surface [7].These asymmetric peaks in the SC gap arise from the interaction between an isolated spin on CoPC and Cooper pairs in NbSe 2 , i.e. the formation of YSR states.Despite CoPc adsorption being driven by the weak van der Waals forces, the magnetic coupling is still sufficiently strong to result in the formation of YSR states.DFT calculations predict that CoPC on NbSe 2 has the same spin as in the gas phase, S = 1/2.The inset of Fig. 1c shows the spin density of CoPC on NbSe 2 surface, which is derived from the cobalt d z 2 orbital, as expected.
Because of the 2D nature of NbSe 2 , the YSR states have been shown to have a large spatial extent and to decay more slowly than on a 3D superconducting substrate [21].We have probed the spatial extent of the YSR states in our system by recording spectra along a line over a single CoPc as shown in Fig. 1e (zero corresponds to the centre of the molecule).It is seen from the figure that the YSR states persist > 2 nm from the centre of the molecule.
In addition to the very slight energy variations of the YSR resonance over the molecule (distances < 1 nm), the smooth energy variation of the resonances at larger distances is likely to result from tip-molecule interactions (see Supplementary Information for more details).
The YSR wavefunctions on NbSe 2 are expected to have six-fold symmetry with larger spatial extent along the principal crystal directions of the surface [21].We have tested this effect by recording another set of line spectra at a 45 • angle with respect to the data shown in Fig. 1e, and carried out grid spectroscopy experiments (Supplementary Information).These experiments suggest that YSR wavefunctions are anisotropic also in our system.The coupling of subgap states is illustrated in Fig. 2a, which shows schematically how two nearby YSR states hybridize and form bonding and antibonding combinations.In this case, phenomenological theory predicts a splitting that depends on the dimensionless coupling strength α (i.e the energy of an individual YSR state), the relative orientations of the two (classical) spins and the separation between the spins [4,15,20,26] (Supplementary Information).α can be determined from experiments on a single impurity and we estimate a value of α ≈ 0.5.The splitting oscillates with a period determined by k F (See Fig. 2b) and obtains its maximum value for a ferromagnetic alignment of the spins.Coupling of antiparallel spins (dotted line in Fig. 2b) does not result in energy splitting but the degenerate energy level is slightly shifted from the individual YSR energy.
Using STM lateral manipulation [27], we successfully constructed molecular dimers with different separations (see Supplementary Information for details).In all cases, we have manipulated one molecule of the dimer and have recorded the spectrum on the one which has not been moved.In this way, we have made sure that the target molecule is always at the same adsorption side.Figure 2c shows point spectra on the molecular dimers with different separations (STM images of the dimers are shown in the insets).Unlike in the case of a single CoPc molecule where we observed only one pair of YSR resonances in the SC gap, now, the main YSR peaks are split to two peaks (left panel).In order to remove the influence of the superconducting density of states of the tip, numerical deconvolution was performed to directly extract the local density of states (LDOS) of the YSR states due to the CoPc molecules.Figure 2d shows the measured LDOS, where the split YRS states (left panel) can be compared to the single YSR peaks at ±0.63 meV.It is evident that the energy positions of split peaks change depending on distance between the dimers.
We observe a maximum splitting of ∼ 0.5 meV, which is very similar to the predictions from a simple continuum model for experimentally relevant distances (Fig. 2b).Magnetic coupling between the spins (through e.g.RKKY interaction or superexchange) is expected to be weak compared to kT at T = 4.2K.Therefore the spins in a dimer are randomly oriented (thermal average).Since majority of random orientations give rise to energy splittings comparable to the maximum splitting, we would observe experimentally essentially the same splitting as the maximum value (see Supplementary Information for details).
Other possible reasons for splitting of the YSR states include different angular momentum scattering contributions, individual d orbitals acting as separate scattering potentials or lowenergy excitations due to magnetic anisotropy or vibrations [8][9][10][28][29][30][31].All of the above scenarios are predicted on a single magnetic impurity, where we always observe only a single pair of YSR states.This leaves the magnetic coupling between the impurities through the SC medium and the formation of coupled YSR states as the natural explanation of the observed spectra.
The right panel of Figure 2c shows a different set of dI/dV spectra on the molecular dimers.While not identical, the intra-dimer distances are over a similar range as in left panel of Fig. 2c.Surprisingly, small changes in the distance result in drastic changes in the spectra and we observe alternating single and split YSR state behaviour.Once the impurity separations exceed ∼ 2.5 nm, we only observe response consistent with single impurity YSR states.As indicated above, we do not expect the spins to be strictly antiparallel, and the expected distance dependence (period of k F ) is not consistent with these rapid changes between split and non-split dimer states.
In order to understand this behaviour, we have to go beyond the continuum description of the SC substrate.The Fermi surface of the NbSe 2 is anisotropic [32], which gives rise to star-shaped structures in the LDOS around a magnetic impurities and vortices [21].In addition, the YSR wavefunctions have atomic scale oscillations arising from the Bloch part of the wavefunction.As the coupling of the YSR states is determined by the wavefunction overlap, these oscillations are a potential source for the atomic scale variations in the observed YSR splitting (see below).We model these effects using a next-nearest neighbour tight-binding model [15] (see Supplementary Information for details).Calculated LDOS of a single impurity is shown in Fig. 3a.The wavefunction strongly reflects the crystal symmetry and it is easy to see that the coupling might strongly depend on the relative orientation of the dimer.We have calculated coupled dimers for different positions of the magnetic impurities and the splitting of the YSR states is shown in Fig. 3b.In addition to a very clear six-fold symmetry stemming from the crystal lattice, there are strong atomic scale variations.This variations result from strong changes in the wavefunction overlap (due to the Bloch part of the wavefunction), when one of the impurities is moved by a single lattice site.This is highlighted in Fig. 3c, which shows the YSR splitting as a function of the dimer separation over the experimentally relevant range.It is seen to oscillate wildly, in agreement with the experimental data (shown with solid red squares).
In Fig. 3d,e we have illustrated the bonding and antibonding wavefunctions of two dimers, where one of the impurities has been moved by a single lattice site from (5,0) to (4,1) (the other impurity is at (0,0)).The calculated YSR splitting in these two dimers is 0.35 meV (5,0) and 0.02 meV (4,1).The latter is far below our experimental energy resolution and would not result in an observable splitting of the YSR resonances.The calculated wavefunctions reflect this: while they are delocalized on both impurities in the strongly coupled dimer (Fig. 3d, analogous to H 2 molecule), they are mostly localized on a single impurity in the weakly coupled dimer (Fig. 3e).This highlights that the atomic scale details are important for a detailed understanding of the coupled YSR states.
In conclusion, we have demonstrated the formation and coupling of YSR states on CoPc molecules on NbSe 2 .Using STM lateral manipulation, we have constructed well-defined CoPc dimers and observed coupled YSR states.Experimentally, we find strong variations of the coupling strength depending on the detailed geometry of the CoPc dimer, which can be understood based on the details of the wavefunction overlap of the two impurity states.
The demonstration of coupled YSR states is promising for realization of novel topological states predicted in YSR lattices.[36].Before computing the electronic structure, all atomic forces were relaxed to < 0.01 eV/ Å.The adsorption site of CoPc can be determined from atomically resolved STM images (Fig. S2).The STM feedback current was increased at the top and bottom parts of the image to allow resolving the NbSe 2 lattice.While scanning over the molecules, the setpoint current was reduced in order not to accidentally manipulate the molecules.Extrapolating the atomically resolved NbSe 2 lattice (Se is visible) allows estimating the CoPc adsorption site.Both molecules are adsorbed with the cobalt centre directly on top of a selenium atom in agreement with the DFT calculations.

SPECTRA OVER A SINGLE COPC MOLECULE
We have carried out experiments to probe the nature of the individual YSR states.In addition to the spectra along line profile 1 shown in the main manuscript, we have additional data along profile 2 at a 45 • angle w.r.t.profile 1 (Fig. S2).This data suggests that the YSR state has shorter decay along this direction.It can also been seen that the YSR resonances shift towards the gap edge, which is difficult to consolidate with the picture of the YSR states being eigenstates of the impurity-surface complex.This effect is also seen (to a lesser extent) on the profile 1 in Fig. 1e of the main text.To check if this effect could be caused by the interaction with the STM tip, we have measured spectra at different tip-sample distances in the middle of a CoPc molecule (Fig. S3).We start all the experiments at a distance determined by the set-point conditions and approach the tip by a distance of z offset before recording the spectrum.dI/dV spectra acquired at different tip-sample distances show a shift of the YSR resonance to higher bias starting at roughly z offset = 30 − 40 pm.At around z offset = 80 pm, the YSR resonance has merged with the superconducting coherence peak at the gap edge.The spectrum recorded at z offset = 100 pm has significantly higher noise and broader resonances compared to the other spectra, and the molecule becomes unstable under the tip at larger values of z offset .On the bare substrate, varying the tip-sample distance has no effect on the shape or position of the superconducting coherence peaks.This measurement clearly indicates that there are considerable interactions between the tip and the CoPc molecules at reduced tip-sample distances.This interaction modifies the coupling with the underlying NbSe 2 substrate, as evidenced by the continuous shift of the YSR resonances as a function of the tip-sample distance.Specifically, as the YSR resonance shifts towards the SC gap edge, the tip-molecule interaction reduces the coupling of the magnetic moment with the superconducting substrate.Speculating, we are likely to be in the attractive force regime of the tip-sample interactions and as the molecule is weakly (van der Waals) bonded on the surface, the tip-sample interaction could have an effect on its adsorption height.Alternatively, the screening from the metallic tip could have an effect on the scalar potential at the impurity site, which would also have an effect on the YSR energy.These experiments allow us to conclude that under our normal conditions (spectra in the manuscript are recorded in conditions similar to z offset = 0), the tip-sample interactions do not play a significant role when we carry out the dI/dV spectroscopy in the middle of the molecule.However, as the tip moves towards the sample close to the edges of the molecule, some shifts of the YSR resonances may occur due to the tip-sample interaction.
In order to shed further light into the YSR states on single CoPc molecules, we have mapped them out by performing grid spectroscopy (recording a complete dI/dV spectrum at each scan point).These experiments are quite demanding due to the mobility of the molecules and we carried them out in STM feedback to have enough signal on the substrate and not to have too much current on the molecule.In addition, we used as low currents as possible to minimize tip-molecule interactions and to reduce the changes of accidental lateral  The YSR peaks are mostly localized in the centre of the molecule (Fig. S4e), with faint tails in different directions that seem to coincide with the principal lattice directions of the underlying NbSe 2 substrate and not with the molecular symmetry.There is a corresponding dip in the SC coherence peak intensity (Fig. S4f).While the normalization removes most of the effect of the molecular backbone (which will has an effect on the tunneling barrier between the tip and sample), this is still faintly visible in (Fig. S4e and f CoPc molecules are weakly adsorbed on the NbSe 2 surface, making it easy to laterally manipulate them by the STM tip (Fig. S5.We successfully constructed molecular dimers with different separations by STM manipulation.Manipulation of the CoPc was carried out by placing the tip above the centre of the molecule at 0.1 V bias voltage and increasing the current to 1 nA with the feedback engaged.The tip was then dragged towards the desired location.In order to avoid the potential variability caused by the adsorption orientation, we always manipulated one molecule of the dimer and recorded the spectrum on the one which had not been moved.In this way, we have made sure that the target molecule is always in the same adsorption configuration.Figure S5 shows a series of manipulation steps demonstrating that we only move the target molecule on the surface.The structural models of the CoPc molecules overlaid on the last panel show that the molecules are still not in "contact".

DECONVOLUTION OF THE dI/dV SPECTRA
The tunneling current I at bias V can be calculated from [37] where ρ t and ρ s are the tip and substrate densities of states and f is the Fermi function.
NbSe 2 has an anisotropic gap structure [38], which we approximated by sum of gapped DOS with some broadening and an additional gaussian component, similar to the expressions used before for modelling STM experiments with SC tips [7,39] ρ where A 1,2 , ∆ 1,2 and γ 1,2 are fitted from the spectra measured on a clean NbSe 2 substrate with a superconducting tip.We can fit the experimental spectra extremely well as shown in Fig. S6 using the same parameters for the bulk NbSe 2 substrate and for the superconducting tips prepared by controlled contacts with the clean substrate.After this, using the same tip with a known ρ t , we extract the substrate DOS on the CoPc molecules by a direct numerical deconvolution of Eq. ( 1).This is an example of a Fredholm integral equation of the first kind and can be solved numerically through a matrix equation [40].

Lattice model
In this supplement, we outline in detail the theoretical modelling of coupled Yu-Shiba-Rusinov states on the surface of NbSe 2 .The starting point of the analysis is a lattice description of an s-wave superconductor with Hamiltonian where t ii is an on-site potential and t ij (i = j) are hopping elements between the lattice sites.The second term describes superconducting pairing of electrons with the pairing gap ∆.The operators c † iσ , c jσ create and destroy electrons at lattice site i with spin σ and obey the usual fermionic anticommutation relations.We assume a triangular lattice where the band structure of NbSe 2 can be reproduced by an appropriate choice of parameters t ij .The magnetic impurities are assumed to be local and described by the Hamiltonian where J i , S i are a magnetic coupling and a classical spin vector of an impurity at lattice position n i .Here, we have introduced Pauli matrices σ = (σ x , σ y , σ z ) and the second-quantized spinor operators C n i = (c i↑ , c i↓ ) T .In addition to the magnetic coupling, the impurities may also perturb the superconductor with additional scalar potential parametrized by V 1 and V 2 .
The standard method of solving the eigenstates of the Hamiltonian H 0 + H imp is to generalize the problem to particle-hole space with the basis Ψ i = (c i↑ , c i↓ , c † i↓ , −c † i↑ ) T and to diagonalize the Bogoliubov-de Gennes Hamiltonian which is a 4N 1 ×4N 2 matrix where N 1 and N 2 are the number of lattice sites in the direction of the primitive lattice vectors of a triangular lattice.We have solved the Bogoliubov-de Gennes problem H BdG Ψ = EΨ for a tight-binding model of NbSe 2 with a finite on-site, nearestneighbour and next-nearest neighbour hoppings.The energy spectrum is symmetric with respect to zero with a gap 2∆.Inside the gap we recover two pairs of states, corresponding to the bonding and antibonding combinations of individual YSR states.In the calculations presented in the main text, we have used values of -100 meV for the on-site and -125 meV for the nearest-and next-nearest neighbour hopping parameters [15] and ∆ = 1 meV, where ε k = k 2 2m − µ is the kinetic energy measured from the Fermi surface while the other terms are straightforward counterparts of those present in the lattice model.In the case of a single impurity, the standard calculation leads to a pair of subgap states with energies E = 1−α 2 1+α 2 , with the dimensionless coupling α = πνJS which contains the density of states at the Fermi level ν.Since we are not including a spin-orbit coupling, the single impurity results does not depend on the orientation of the impurity moment.While the two-spin problem does not admit a simple closed form analytical expression, the eigenstates can be solved, for example, by the methods of Ref. [41].This leads to four subgap states

Figure
Figure 1a shows an atomic resolution STM image of a NbSe 2 surface at 4.2 K (See Supporting Information for details on the experiments).It reveals the hexagonal arrangement

FIG. 2 .
FIG. 2. Formation of coupled YSR states on CoPc dimers.(a) Schematic of the formation of coupled YSR states.(b) Calculated evolution of the YSR state energies as a function of the impurity dimer separation for a ferromagnetic (solid lines) and antiferromagnetic (dotted lines) dimer.(c) Set of dI/dV spectra showing split (left panel) and non-split YSR states (right panel).A spectrum measured on an isolated CoPc (black line) is shown in both panels for comparison (feedback setpoint: V = 100 mV, I = 50 pA, z offset = 100 pm).The dotted line shows the position of the SC gap edge at ±2∆/e.(d) The LDOS extracted by numerical convolution from the experimental dI/dV curves.

FIG. 3 .
FIG. 3. Atomic scale details of the coupling between two magnetic impurities.(a) Calculated LDOS of single impurity in the next-nearest neighbour tight-binding model.(b) Calculated splitting as a function of the impurity position (the other impurity is at (0,0)).(c) Calculated (open symbols) and measured (closed symbols) splitting of the YSR states as a function of the distance between the impurities.(d,e) Calculated bonding (left panels) and antibonding (right panels) state LDOS for a strongly coupled (5,0) (panel d) and weakly coupled (4,1) (panel e) dimer.
FIG. S1.STM image of two CoPc molecules (V = 50 mV, I = 1 nA on NbSe 2 and V = 0.7 V, I = 3 pA on CoPc).The STM feedback current was reduced of the CoPc to avoid their accidental manipulation.
FIG. S3. dI/dV spectra (normalized) recorded over the middle of a CoPc molecule at different tip-molecule distances.z offset = 0 pm corresponds to the set-point conditions at V = 20 mV, I = 300 pA.
manipulation of the molecule.The results are shown in Fig. S4, which displays raw data dI/dV slices at the energies corresponding to the YSR (panel c) and the superconducting coherence peaks (panel d) and the same results normalized to take into account the varying tip height (panels e and f).
FIG. S4.Grid spectroscopy on an individual CoPc molecule.(a) Topography image acquired simultaneously with the dI/dV spectra.(b), Topography at enhanced contrast to show the atomically resolved underlying NbSe 2 lattice.(c) dI/dV slice at the bias corresponding to the YSR peak (positive bias).Gray scale is between 0 − 0.03 µS.(d) dI/dV slice at the bias corresponding to the superconducting coherence peak (positive bias).Gray scale is between 0 − 0.02 µS.(e) Maps corresponding to panels c and d, where the dI/dV signal has been normalized by the current at the beginning of the spectrum.Set-point V = 500 mV, I = 5 pA, z offset = 100 pm.
FIG. S5.Lateral molecular manipulation of CoPc on NbSe 2 .Subsequent STM images of CoPc molecules after one of the molecules was laterally manipulated by the STM tip.

±E 1 and
FIG. S7.Angular averaging of the YSR splitting.(a) YSR state energies as a function of the angle of the magnetic moments of the two impurities according to the continuum model (k F a = 3.8, α = 0.5).(b) Same data presented as a histogram, which shows that majority of the angles result in a YSR splitting very close to the maximum value.This would result in clearly split resonances in the experimental spectra.