Anisotropic non-split zero-energy vortex bound states in a conventional superconductor

Majorana bound states (MBSs) in condensed matter systems have attracted considerable interest as a key element for fault-tolerant topological quantum computation.^1–3 They have been theoretically predicted to exist in topological superconductors, which may be realized by low-dimensional magnetic nanostructures, semiconductor nanowires, or topological insulators proximitized to s-wave superconductors.^4–8 However, due to experimental difficulties in the fabrication of these sophisticated nano-scale hybrid systems, significant efforts have been directed toward studies of intrinsic topological superconducting bulk materials, being a promising platform for hosting MBSs in the cores of vortices which form in the presence of an external magnetic field.⁹ Recently, several materials have emerged as potential topological superconductors possessing topological surface states, including iron-based superconductors^10–13 and transition-metal dichalcogenides.¹⁴ MBSs are expected to be observed at zero energy (i.e., the Fermi level) in the center of vortex cores. In particular, the observation of spatially non-split zero-energy states in the local density of states (LDOSs) emanating from the vortex core using scanning tunneling microscopy and spectroscopy (STM/STS) has been regarded as an essential experimental fingerprint for MBSs in these materials.^6,10,12,14

The identification of MBSs based exclusively on spectroscopic signatures is made challenging by the simultaneous presence of topologically trivial bound states. These include Yu–Shiba–Rusinov (YSR) states in magnet–superconductor hybrid systems,^15–17 Andreev bound states in semiconductor nanowire–superconductor hybrid systems,¹⁸ and Caroli–de Gennes–Matricon bound states (CBSs) in vortex cores;¹⁹ the latter were observed in topologically trivial and conventional superconductors such as NbSe₂^20,21 and on the (110) surface of bulk Nb²² as well as in unconventional superconductors.^23–26 The energy of CBSs for conventional s-wave superconductors is given by E_μ = μΔ/(k_Fξ),^19,27 where μ is a half-integer, Δ is the superconducting gap, k_F is the Fermi wave vector, and ξ is the superconducting coherence length. In principle, this makes it possible to distinguish the zero-energy MBS from CBSs which cannot be located at zero energy. However, the energy spacing between CBS is quite small in most superconductors due to the small value of Δ and the large value of k_Fξ. Recently, it has been demonstrated that for several Fe-based superconductors, having a topologically non-trivial band structure, the energy difference Δ/(2k_Fξ) separating the lowest-energy CBS from the MBS is sufficiently large to be resolved by low-temperature STM/STS, representing the so-called quantum-limit case.^10–12,28 Why zero-energy MBSs have only been observed for a fraction of vortices by STS in real space for a given sample remains an open question.^10–12,23

For materials not in the quantum limit, where Δ/(k_Fξ) is tiny, the spectral feature of a CBS is a pronounced peak at zero energy, which appears identical to that of a zero-energy MBS, making it difficult to distinguish them by simply observing the spectrum at the center of the vortex core. Such zero-energy CBSs have been observed for several topologically trivial superconductors.^20–22,24 CBSs may theoretically and experimentally be distinguished from MBSs even away from the quantum limit by their spatial distribution around the vortex:^6,29,30 the energy of CBSs increases away from the core, leading to a splitting in the spectra with increasing distance, while the MBS remains at zero energy. Based on these properties, the simultaneous observation of the anisotropic non-split and the split zero-energy states has been interpreted as a signature for the coexistence of MBSs and CBSs within a single vortex.^14,31 The difference between spatially splitting and non-splitting states as a criterion for identifying MBSs has not been challenged so far. However, the splitting of CBSs with distance is based on the assumption of an isotropic electronic structure for the superconductor,^32–34 while the materials studied in Refs. 14 and 31 have been demonstrated to possess an anisotropic Fermi surface. Calculations of quasiparticle LDOSs at vortices based on the quasi-classical approximation have revealed the crucial role of the anisotropy of the superconducting gap and the Fermi surface for its spatial distribution.^33,35

Here, we present an experimental proof for the presence of anisotropic non-split zero-energy bound states in the vortex cores of a conventional and topologically trivial superconductor lanthanum based on STM and STS investigations. The observed spatial distribution of the zero-energy states is reproduced by quasi-classical theory calculations taking into account the anisotropic Fermi surface formed by the quasi-two-dimensional surface states of La(0001),^36,37 thereby identifying them as topologically trivial CBSs. Accordingly, non-split zero-energy states can no longer be considered as evidence for the existence of MBSs within vortices, especially, in the non-quantum limit.

Bulk-like α-lanthanum films with thickness larger than the superconducting coherence length [ξ_lit = 36.3 nm (Ref. 38)] were epitaxially grown on a (0001)-oriented rhenium single crystal substrate (see Sec. IV for more details). Superconducting La-coated PtIr tips were prepared to improve the energy resolution in tunneling spectroscopic measurements beyond the Fermi–Dirac limit.^37,39,40 All STM/STS measurements were carried out at T = 1.65 K, which is well below the superconducting transition temperature of both the La films and the La-coated tip, T_c ∼ 5–6 K.⁴¹ 

We first characterized the atomic and electronic structure of the atomically flat surface of the La(0001) film in zero magnetic field. In Fig. 1(a), a local tunneling conductance (dI/dV) spectrum obtained on the clean La(0001) surface shows a clear peak at +110 meV reflecting the band edge of the surface states originating from the d_z²-like atomic orbitals of La,^36,37 together with a superconducting gap feature close to the Fermi energy. This peak feature is considerably different from the V-shaped dip in the vicinity of the Fermi level, which has been interpreted as a signature of Dirac cones in the dispersion of topological surface states.^10,14 The calculation of symmetry-based indicators (see Sec. IV and Refs. 42 and 43) of the band structure of La does not predict the presence of topologically protected surface states, either. According to our calculations within the framework of density-functional theory (DFT), the observed surface states of La(0001) crossing the Fermi level are almost spin degenerate, with a negligible Rashba splitting on the order of a few meVs. In contrast to this, the theoretical description of MBSs at vortex cores requires the formation of an effective p-wave pairing, facilitated by a single non-degenerate band crossing the Fermi level. The inset of Fig. 1(a) shows an atomic-resolution STM image with the hexagonal arrangement of surface La atoms. To resolve the surface electronic structure in momentum space, we measured energy-dependent differential tunneling conductance (dI/dV) maps on La(0001). In Fig. 1(b), we present a dI/dV map obtained at V = +2.0 mV showing quasiparticle interference (QPI) patterns resulting from surface electron scattering at residual atomic defects on the La(0001) surface. The corresponding Fourier transformation (FT) of the dI/dV map in Fig. 1(c) clearly shows a hexagonally shaped pattern which can be attributed to the anisotropic Fermi surface contour of the quasi-2D surface electronic band.³⁷ Q₁ and Q₂ denote the scattering vectors along the high-symmetry crystallographic directions where the scattering of surface quasiparticles occurs toward the vertex and the flat part of the hexagon, respectively. In Ref. 37, it was demonstrated using FT analysis at various energies in agreement with ab initio calculations that the hexagonal Fermi contour originates from intra-band scattering of the quasi-2D surface states, leading to the large spatial extension of YSR bound states around magnetic impurities on the La(0001) surface.

Next, we characterized the superconductivity of the La sample and the STM tip, and verified the presence of vortex bound states by applying a magnetic field B perpendicular to the La(0001) surface. The field-free tunneling spectrum (black, solid) at the bottom of Fig. 1(d) reflects the characteristics of a superconductor–superconductor (S–S) tunnel junction showing strong coherence peaks at E_coh = ±|Δ_S+Δ_T|= ±1.6 mV, where E_coh is the energy of the coherence peaks and Δ_{S (T)} is the superconducting gap of the La surface (tip). Note that Δ_T/e (e is the elementary charge) corresponds to the Fermi energy on the La(0001) surface, i.e., zero energy. By simulating the spectrum (orange, dotted) considering the superconducting state of the La-coated tip, we deduce values of Δ_S = 1.01 meV and Δ_T = 0.65 meV (see the supplementary material for more details on this procedure). They have been used for obtaining the deconvoluted surface DOS without the contribution of the tip DOS as shown in Fig. 1(e). By applying a magnetic field of B = 0.08 T perpendicular to the surface, we created vortices penetrating the La film. As superconductivity is quenched within the vortex cores, resulting in reduced intensity of the coherence peaks, we could visualize the spatial distribution of the vortices as depressions around the vortex centers by measuring the dI/dV map at E_coh as shown in Fig. 1(f). At the center of individual vortices, we observe pronounced peaks in the tunneling spectrum [red curve in Fig. 1(d)] at V_S = +Δ_T/e, i.e., equivalent to zero energy [red curve in Fig. 1(e)], which are absent from the spectrum obtained away from the vortices [blue curves in Figs. 1(d) and 1(e)]. The pronounced peaks in the vortex cores can be attributed to zero-energy vortex bound states or CBSs, which, in principle, are due to the coherent superposition of the Andreev reflected electron and hole from the inhomogeneous superconducting pair potential around the vortex. The spatial distribution of CBS can be visualized in the dI/dV map at V_S = +Δ_T/e showing an anisotropic star shape extending more than ∼100 nm along the six directions equivalent to the Q₂ direction [Fig. 1(g)]. A star-shaped LDOS modulation around vortices has been observed previously for the conventional superconductor NbSe₂ with an anisotropy of the superconducting gap.^20,21

In order to get more insight into the energy-dependent spatial distribution of the anisotropic vortex bound states, we present the spatial evolution of the LDOS around an individual vortex based on tunneling spectroscopic measurements in Figs. 2(a)–2(d) (see also the supplementary material Movie 1). At zero energy [E = 0 or V_S = +Δ_tip/e, Fig. 2(a)], the vortex has an anisotropic star shape with six rays extending from the center of the vortex toward the directions of the six flat parts of the Fermi contour (Q₂) as shown in Fig. 1(c). As the energy is increased, each ray is split into two and they move away from the center while keeping the direction of the extension. Since the separation between two split rays becomes larger and larger at higher energy, this leads to the expansion of the vortex and simultaneously to the appearance of a depression of the LDOS at the vortex center as shown in Figs. 2(b)–2(d).

Figures 2(e)–2(g) show the spectral evolutions of the LDOS along the three radial lines with θ = 0°, 15°, and 30° passing through the center of the vortex as the origin, respectively. All spectroscopic maps show highly symmetric evolutions with respect to the vortex center. At θ = 0°, where the radial line is parallel to the Q₁ direction, the zero-energy peak of the vortex bound state at the center is split into two dispersive inner and outer branches, i.e., a linear and a curved X-shaped splitting is observed away from the center. Remarkably, the splitting of the linear branch becomes smaller by increasing θ [Fig. 2(f)]. Finally, as shown in Fig. 2(g), the inner lines merge into a single line, namely, a spatially non-split zero-energy state along the radial direction with θ = 30°, parallel to Q₂ pointing toward the flat parts of the Fermi contour. On the other hand, the outer branch exhibits nearly the same curved X-shaped splitting for all values of θ. In Fig. 2(h), the angle-dependent LDOS plot along a circular line (white dotted line in Fig. 2(a), 30 nm away from the vortex center) clearly shows that the spatially non-split zero-energy states appear at every radial line for θ = 30°+N·60° due to hexagonal symmetry (with N an integer), while it shows energy splitting elsewhere, θ ≠ 30°+N·60°.

Although a similar anisotropic extension of spatially non-split zero-energy bound states has been interpreted as an essential signature of MBSs in 2M-WS₂ and Bi₂Te₃/FeTe_0.55Se_0.45,^14,31 this explanation can be ruled out here due to the topologically trivial electronic structure of La. The difference between the energy levels of CBSs is Δ_S/(k_Fξ)$≈$ 10 μeV based on the measured Δ_S = 1.01 meV [see Fig. 1(d)], k_F ≈ 2π·0.5 nm⁻¹ [based on Fig. 1(c)], and ξ = 36.3 nm for bulk La. Since this is well below the experimental resolution, the considered system is not in the quantum limit and the pronounced zero-energy vortex bound state can be identified as a CBS. The star-shaped LDOS modulation around the vortex has been similarly observed in a conventional superconductor, NbSe₂,^20,21 and the modulation of the LDOS has been demonstrated to be anisotropic in vortices on the (110) surface of bulk Nb.²² Possible explanations of the anisotropic extension include vortex–vortex interactions,^44–46 the crystal field effect,^44,45,47 or the anisotropy of either the Fermi surface or of the superconducting gap,^33–35 regardless of topological properties. The model considering the crystal field effect has not been sufficient to understand the energy-dependent LDOS maps around a vortex.⁴⁶ The influence of vortex–vortex interactions in the present system appears to be negligible, since the extension of bound states is hardly influenced by inhomogeneities in the vortex lattice, as shown in Figs. 1(g) and 1(f). The field-free tunneling spectrum in Fig. 1(d) shows a nearly ideal single-component BCS-type superconducting gap, excluding the importance of gap anisotropy. This is in direct contrast to NbSe₂ for which previous theoretical studies explained the emergence of zero-energy peaks as well as the spatial distribution of LDOS around the vortex based on the anisotropy of the superconducting gap^33,34 (see the supplementary material).

To elucidate the major role of the Fermi surface for the spatial distribution of vortex bound states, we calculated the LDOS around a single vortex based on a quasi-classical Green's function method within the Kramer–Pesch approximation (KPA), which is suitable for low-energy quasiparticle excitations within the vortex core.⁴⁸ The quasi-classical approximation is justified since the considered system is not in the quantum limit as mentioned above. Energy-resolved QPI measurements together with DFT calculations of the surface states indicated that the magnitude of Fermi velocity v_F, and consequently the coherence length ξ, is similar along the Q₁ and Q₂ directions,³⁷ but the quasiparticle trajectories are focused along the latter due to the shape of the Fermi contour.

Figures 3(a)–3(d) show the calculated results of quasiparticle LDOS maps around an isolated superconducting vortex by considering the realistic Fermi contour of the La(0001) surface [the inset of Fig. 3(a)]. The nearly-hexagonal Fermi contour of the La(0001) surface has been extracted from first principles and can be compared with the quasiparticle interference pattern at the Fermi energy [see Fig. 1(c) and Ref. 37]. Remarkably, the sixfold star-shaped extension observed at zero energy [Fig. 3(a)] splits into twelve pairwise parallel rays at higher energies [Figs. 3(b)–3(d) and see also the supplementary material Movie 2], reflecting the trajectory of scattered quasiparticles with a higher impact parameter measured from the vortex center. Since the strongly anisotropic Fermi surface selects preferential directions for the quasiparticle velocities, the trajectories primarily run parallel to Q₂ and symmetrically equivalent directions at all energies. We found an enhancement of the quasiparticle LDOS at the positions where the extension lines are interconnected with each other. Furthermore, Figs. 3(e)–3(g) show the calculated LDOS maps as a function of energy and distance from the vortex center along radial lines for θ = 0°, 15°, and 30° angles measured from the Q₁ direction, respectively (see also the supplementary material Movie 3). As the crossing point of quasiparticle trajectories with an angle of 60° between them is moved away from the vortex center with increasing energy, a cross-shaped pattern is formed in the spectral evolution for θ = 0° with increasing CBS energies away from the center, as reflected by high-intensity hexagons in Figs. 3(b)–3(d). The branches of the cross approach each other for θ = 15° and merge for θ = 30°, resulting in a non-split zero-energy bound state. All the above features of the simulated spectra are in excellent agreement with the experimental observations in Fig. 2. Within the quasi-classical description, the non-split zero-energy state is a reflection of the delocalized quasiparticle trajectories at zero energy due to the heavily weighted Fermi velocity along the direction where the Fermi contour is flat.^37,49 On the other hand, for an isotropic Fermi contour, the spatial distribution of vortex bound states is isotropic (or circular), see the supplementary material. In addition to the anisotropic distribution of the CBS, the shape of the Fermi surface plays an important role in the long-ranged extension of the CBS due to the effective reduction of dimensionality compared to the case of an isotropic Fermi surface (see Ref. 37 and the supplementary material). Since we found nearly identical spectroscopic signatures for the topologically trivial elemental superconductor La as reported before for other types of superconductors with an anisotropic Fermi surface, this clearly indicates that in systems far from the quantum limit a non-split zero-energy bound state in a superconducting vortex is not a characteristic signature of a MBS, but rather an effect of the anisotropic Fermi surface on the spatial distribution of the CBS.

In summary, we have studied the spatially and energy-resolved LDOS of vortex bound states in a conventional and topologically trivial superconducting La(0001) sample using high-resolution STM/STS. We have shown that the spatial distribution of the vortex bound states reflects the quasiparticle trajectory around the vortex via the anisotropic Fermi surface of La(0001). Anisotropic non-split zero-energy vortex bound states coexisting with split branches at finite energy within the vortex core have been observed experimentally, consistent with the calculated LDOS based on quasi-classical Green's function theory. Our findings stimulate a revision of the interpretation of previously reported non-split zero-energy bound states within vortex cores of different classes of superconducting materials having an anisotropic electronic band structure.

The rhenium single crystal was prepared by repeated cycles of O₂ annealing at 1400 K followed by flashing at 1800 K to obtain an atomically flat Re(0001) surface.⁵⁰ More than 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 min. The films were thicker than the coherence length of bulk La (ξ = 36.3 nm) in order to suppress the inverse proximity effect from the Re. The surface cleanliness and the thickness of the La films were checked by STM after transferring the sample into the cryostat. A La-coated PtIr tip was used for STM/STS measurements to improve energy resolution in the spectra at a given experimental temperature via the superconducting tip. To fabricate the La-coated tip in situ, a mechanically polished PtIr tip was intentionally indented into the clean La film. The superconducting La–La junctions were confirmed by observing the Josephson tunneling current at V = 0.0 mV as well as the spectral signature of multiple Andreev reflections (see supplementary Fig. S1).

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 μV_rms 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).

The surface electronic structure of La(0001) was investigated using the screened Korringa–Kohn–Rostoker method.⁵¹ The Fermi surface formed by the quasi-2D surface bands was determined via the calculation of the Bloch spectral function in the top La layer. Details of the calculations are reported in Ref. 37.

The symmetry-based indicators were calculated using the Check Topological Mat. program available on the Bilbao Crystallographic server.⁴² The topological classification of bulk La is available in existing online databases.^42,43 However, the nominal valence bands of La used in these databases are above the Fermi level in the vicinity of the Γ point, where the quasi-2D surface bands on the La(0001) surface are found. To deduce the topological character of these surface states, the nominal number of valence bands has to be reduced by four compared to the database value. The electronic structure was recalculated using the Vienna ab initio simulation package (VASP);^52–54 see the supplementary Fig. 6 for the band structure along high-symmetry lines. La was found to be an enforced semimetal also for lower nominal band filling; see the supplementary material Data 1 for the input file of Check Topological Mat. and supplementary material Data 2 for the reduction of bands according to the irreducible representations, using the naming conventions of the Bilbao Crystallographic server. These results do not predict a topological nature of the observed surface states.

In the quasi-classical approximation, the Andreev scattered quasiparticles confined inside the vortex can be described by considering the quasiparticle trajectory which passes through the position characterized by the impact parameter.⁵⁵ The binding energy of the quasiparticles is a monotonically increasing and continuous function of the impact parameter, which is coupled to the angular momentum of the quasiparticle. At each spatial point r, the LDOS of the quasiparticles can be calculated by integrating over the number of quasiparticle trajectories passing through r at a given impact parameter. Since the trajectory follows the direction of Fermi velocity, which reflects the electronic band structure, this facilitates the correlation of the contribution of the Fermi surface with the spatial distribution of the quasiparticle LDOS around the vortex. The local density of states is calculated on the basis of the Kramer–Pesch approximation in the quasi-classical Eilenberger framework, which is appropriate in the low-energy regime.⁴⁸ The Fermi surface for these calculations was chosen to approximate the results of the ab initio calculations by interpolating between the perfectly circular and the ideal hexagonal Fermi surfaces using a deformation parameter. Details of this interpolation procedure are given in Ref. 37.

See the supplementary material for details on the numerical analysis for the experimental results, and the theoretical simulations including corresponding notes, figures, additional data and videos.

H.K. and R.W. conceived and designed the experiments. H.K. and D.S. carried out the STM/S experiments. H.K. analyzed the data. L.R. performed the first-principles calculations and their analysis. Y.N. performed the quasi-classical Green's function calculations. R.W. supervised the project. H.K., L.R., and R.W. wrote the paper. All authors discussed the results and contributed to the paper.

We thank J. Wiebe, D. Wang, and D. Morr for fruitful and useful discussions, as well as E. Simon, M. V. Valentyuk, and A. I. Lichtenstein for discussions and support with the first-principles calculations. This work was supported by the European Research Council via Project No. 786020 (ERC Advanced Grant ADMIRE) at the University of Hamburg. L.R. is supported by the Alexander von Humboldt Foundation. This work was partially supported by the “Topological Materials Science” (No. 18H04228) JSPS-KAKENHI on Innovative Areas.

The data that support the findings of this study are available within the article and its supplementary material.