Long-range focusing of magnetic bound states in superconducting lanthanum

Quantum mechanical systems with long-range interactions between quasiparticles provide a promising platform for coherent quantum information technology. Superconductors are a natural choice for solid-state based quantum devices, while magnetic impurities inside superconductors give rise to quasiparticle excitations of broken Cooper pairs that provide characteristic information about the host superconductor. Here, we reveal that magnetic impurities embedded below a superconducting La(0001) surface interact via quasiparticles extending to very large distances, up to several tens of nanometers. Using low-temperature scanning probe techniques, we observe the corresponding anisotropic and giant oscillations in the LDOS. Theoretical calculations indicate that the quasi-two-dimensional surface states with their strongly anisotropic Fermi surface play a crucial role for the focusing and long-range extension of the magnetic bound states. The quasiparticle focusing mechanism should facilitate the design of versatile magnetic structures with tunable and directed magnetic interactions over large distances, thereby paving the way toward the design of low-dimensional magnet–superconductor hybrid systems exhibiting topologically non-trivial quantum states as possible elements of quantum computation schemes based on Majorana quasiparticles.

The FS is defined in the reciprocal space of crystals, but its influence is also manifest in real space in the vicinity of impurities via screening of the impurity potential resulting in a modulation of the electron density, also known as Friedel oscillations 2 . Anisotropic FSs effectively focus and guide the scattered quasiparticles along specific directions, enabling the characterization of non-magnetic and magnetic impurities at large distances 3,4 .
In a conventional superconductor, magnetic impurities locally break Cooper pairs and generate quasiparticle excitations inside the superconducting gap, known as Yu-Shiba-Rusinov (YSR) bound states [5][6][7] . Arranging the magnetic impurities into low-dimensional arrays has recently attracted great interest in the pursuit of Majorana bound states (MBS), which are promising as a key element in topological quantum computation [8][9][10][11][12] . The long-range and anisotropic interaction of the impurities reflected in the spatial extent of the YSR states allows the tailoring of MBSs in low-dimensional magnetic systems on superconductors.
Scanning tunneling microscopy (STM) and spectroscopy (STS) measurements provide access to the YSR states at the atomic level with high spatial and energy resolution. Using these techniques, the spatial extension of YSR states for most magnetic impurities on surfaces of three-dimensional (3D) superconductors was found to be limited to the range of a few atomic lattice constants [12][13][14][15][16] . A significantly longer extension of a few nanometers was observed for Mn atoms on Pb(111), and has been attributed to the projection of the bulk FS onto this crystallographic direction 17 . Recently, Ménard et al. reported a slow decay of YSR states for subsurface Fe impurities inside the layered superconductor 2H-NbSe2, observable at distances up to several nanometers 18 . This was attributed to the essentially two-dimensional (2D) electronic structure caused by the weak hybridization between the layers, and it was argued that layered materials are preferable for studying long-range interactions between YSR states. 3 The sizeable interaction between YSR states around magnetic molecules on NbSe2 was indeed demonstrated by Kezilebieke et al 19 .
Here, we demonstrate that YSR bound states with anisotropic and extremely long-range extension, reaching up to several tens of nanometers, can be observed on the (0001) surface of the 3D elemental superconductor La using low-temperature STM/STS. This is remarkably longer than what has been reported based on previous STM/STS investigations on any 3D or even quasi-2D superconductor [12][13][14]16,18 . Furthermore, we confirm the presence of the interaction between impurities at distances of several nanometers. Corroborated by firstprinciples and model calculations, we show that this extraordinary long-range modulation of the local density of states (LDOS) due to YSR impurities provides direct information on the anisotropic FS formed by the quasi-2D surface states through the quasiparticle focusing effect.
We have studied bulk-like -lanthanum films epitaxially grown on a Re (0001)  The magnetic impurities are clearly identifiable in differential tunneling conductance (dI/dV) maps taken above the atomically flat La(0001) surface at a bias voltage corresponding to states inside the superconducting gap of La, as shown in Fig. 1 However, it is also possible that different magnetic elements being present as natural impurities even in high-purity La material may contribute to the different spatial and spectroscopic signatures of the observed YSR impurities. In the following, we mainly focus on YSR1-type impurities, characterized by the highest modulation intensity. and an oscillatory modulation characteristic of YSR states 6,7 . The oscillations are observable 4 up to ~30 nm away from the impurity along the beams, and up to ~10 nm distance along the directions halfway between the beams, both being considerably longer than in previously reported STM/STS studies of YSR bound states [12][13][14]16,18 .
In Fig. 1(d), the YSR quasiparticle excitations show up as a pair of resonances in the dI/dV spectra inside the superconducting gap of La (La,surface~1 meV). Taking into account the convolution of the Bardeen-Cooper-Schrieffer-type DOS of the sample and the La-coated superconducting tip in the measured tunneling spectrum 21 , the binding energies of the YSR states are found at EYSR=0.55 meV (red and blue arrows). Their symmetric position with respect to the Fermi level reflects the particle-hole symmetry of the superconductor.
The spatial extension of the YSR states shown in Figs. 1(b) and (c) has been quantitatively analyzed by taking differential tunneling conductance profiles ( Fig. 1(e)) and spectroscopic maps ( Fig. 1(f)) along lines crossing the center of the magnetic impurity. The period of the oscillations, corresponding to half of the Fermi wavelength (F/2) [Supplementary Note 5], is found to be ~0.99 nm along the beams (Q2 direction in Fig. 1(e)) and ~0.90 nm along the direction halfway between the beams (Q1 direction), respectively, for both positive and negative bias voltages. The phase shift between the spectroscopic maps at EYSR, determined by the binding energy 17,18 , is clearly conserved as far as 30 nm away from the impurity along the Q2 direction, confirming that even at this distant point the LDOS modulation can be attributed to the coherent particle-hole symmetric YSR states.
To address the origin of the anisotropic long-range extension of the YSR states, we turn our attention towards the surface electronic structure of La(0001). The dI/dV spectrum averaged over a flat terrace shows a sharp resonance peak at E=+110 meV, which reflects the band edge of the surface state of La(0001) as shown in the left panel of Fig Since both the quasiparticle scattering at the Fermi energy and the YSR states find their origin in the FS, it is reasonable to consider the analogies between them. They share the same period of oscillation along the two directions Q1 and Q2. Due to the quasiparticle focusing effect, both types of oscillations are extended to a longer range along the flat directions of the FS than along the vertices, as can be seen when comparing Fig. 2  intra-band scattering of quasiparticles originating from the surface band (inner hexagon in Fig.   2(d)). The QPI pattern and the YSR states appear to be insensitive to the bulk band (outer hexagon present in the BSF calculations in Fig. 2(d)), which is rotated by 30with respect to the surface band. This implies that the surface band of La(0001) plays a major role for the quasiparticle scattering at the surface while the contribution of the bulk band is weak or lacking.
Note that similarly to the present case, only one of the two FSs contributes to the observed YSR states in Pb(111) 17 as well as in NbSe2 18 .
The QPI patterns and the first-principles calculations thereby reveal two major contributions to the anisotropic long-range extension of the YSR states. The first is the quasi-2D electronic structure, because primarily electrons in the surface band contribute to the scattering processes.
The second is the quasiparticle focusing effect, causing beam-like extensions perpendicular to the flat parts of the anisotropic FS. In order to quantify the role of the shape of the FS on the spatial extension, we performed 2D atomistic model calculations parametrized by the firstprinciples results to determine the energy and the LDOS modulation of the YSR states [see Methods]. As shown in Fig. 3(a), an ideal circular FS leads to an isotropic decay of the YSR intensity, following a 1/r power law at intermediate distances where r is the distance from the impurity 17, 18 . For the hexagonal Fermi surface in Fig. 3(b), the LDOS along the beam directions perpendicular to the flat sides of the FS decays slower than 1/r, effectively reducing the dimensionality due to the strong focusing effect [see Supplementary Note 5 and Supplementary   6 Video] 3,24 . By considering the realistic FS of La(0001) obtained from the first-principles calculations, which is slightly rounded from the ideal hexagon as shown in Figs. 2(c) and (d), our numerical results (Fig. 3(c), left) are consistent with the experimentally obtained dI/dV map for the YSR1 impurity (Fig. 3(c), right) with regards to the shape and the range of the extension, as well as the oscillation period.
It is worth mentioning that the focusing due to the anisotropic Fermi surface and due to the dimensionality of the system cannot rigorously be distinguished, because they are essentially different manifestations of the same mechanism. For example, an isotropic two-dimensional  Supplementary Fig. 5]. However, we find that a circular two-dimensional or cylindrical three-dimensional Fermi surface could not fully account for the long-range extension, since they would lead to a characteristic 1/r decay that is faster than the spatial decay observed in the experiments (~1/r 0.92 ).
Finally, we demonstrate that the YSR impurities on La(0001) exhibit a significant longrange coupling between them. The interaction between YSR impurities results in a variation of EYSR compared to the case of a single impurity 19,[25][26][27][28] . Figure 4(a) presents four YSR1-type impurities: two of them (A, B) at about 3.6 nm distance roughly along the Q2 direction, and two isolated ones (I1, I2) being significantly further away. Interestingly, due to the interference between the YSR states of the impurities A and B, the LDOS shows a twin-star pattern with inhomogeneous modulations around the paired YSR1-type impurities. As shown in Fig. 4(b), the binding energy EYSR=0.56 meV for the isolated I1 and I2 impurities are consistent with the one for the YSR1 impurity in Fig. 1(d). However, on the A and B impurities a sizeable shift of around 90 eV is observed, being identical for the two atoms. This clearly indicates that the extended YSR bound states give rise to mutual long-range interactions between the magnetic impurities.
In conclusion, we have demonstrated how the anisotropic FS formed by the surface band on La(0001) causes magnetic bound states in the superconducting phase to extend coherently over tens of nanometers along preferential directions, the extension being about an order of 7 magnitude longer than what has previously been reported for other bulk elemental superconductors. The results confirm that the slow power-law decay and long-range extension of the YSR bound states, which may be interpreted as an effective reduction in dimensionality, can be observed not only in layered systems, but also in elemental bulk superconductors exhibiting distorted FSs with quasi-2D bands. Moreover, the large spatial extent of YSR states results in a coupling between the magnetic impurities at remarkably long distances. These

Data availability
The main datasets generated and analyzed during the current study are available from the authors upon reasonable request.

Code availability
The computational codes for the tight-binding model calculations of the band dispersion and the spatially-resolved YSR bound states are available from the authors upon reasonable request.

Methods (within 3000 words)
Preparation of the sample and tip. The rhenium single crystal was prepared by repeated cycles of O2 annealing at 1400 K followed by flashing at 1800 K to obtain an atomically flat Re(0001) surface 33 . The over 50 nm thick lanthanum films were prepared in situ by electron beam evaporation of pure La pieces (99.9+%, MaTeck, Germany) in a molybdenum crucible at room temperature onto a clean Re(0001) surface, followed by annealing at 900 K for 10 minutes 20 .
The films were thicker than the coherence length of bulk La (=36.3 nm) 34

STM/STS measurements.
All the experiments were performed in a low-temperature STM system (USM-1300S, Unisoku, Japan) operating at T=1.65 K under ultra-high vacuum conditions. Tunneling spectra were obtained by measuring the differential tunneling conductance (dI/dV) using a standard lock-in technique with a modulation voltage of 30 Vrms and a frequency of 1075 Hz with opened feedback loop. The bias voltage was applied to the sample and the tunneling current was measured through the tip using a commercially available controller (Nanonis, SPECS).
Electronic structure calculations. The ab initio calculations were performed in the dhcp structure of La, using the lattice constants of a0=3.772 Å and c0=12.144 Å. As a first step, slab calculations were performed using the Vienna Ab initio Simulation Package (VASP) [36][37][38] to determine the relaxation of the top atomic layer on the La(0001) surface. The slab consisted of 8 atomic layers and an about 24 Å thick vacuum layer to avoid interactions between the periodic images of the supercell. The potential was parametrized based on the Perdew-Burke-Ernzerhof method 39 , and a 15×15×1 Monkhorst-Pack mesh was used to sample the Brillouin zone. It was found that the outermost layer relaxes inwards by about 5.8% of the interlayer distance c0/4.
where is the energy, ∥ is the in-plane wave vector, is the layer index, + is the retarded Green's function, and integration is performed over the atomic sphere. As illustrated in Fig where , ∈ {↑, ↓} is a fermionic annihilation operator. is the quasiparticle spectrum in the normal state, which was determined from the ab initio calculations; see Supplementary   Table 1