Band-resolved Caroli–de Gennes–Matricon states of multiple-flux-quanta vortices in a multiband superconductor

Superconductors are of type I or II depending on whether they form an Abrikosov vortex lattice. Although bulk lead (Pb) is classified as a prototypical type-I superconductor, we show that its two-band superconductivity allows for single-flux-quantum and multiple-flux-quanta vortices in the intermediate state at millikelvin temperature. Using scanning tunneling microscopy, the winding number of individual vortices is determined from the real space wave function of its Caroli–de Gennes–Matricon bound states. This generalizes the topological index theorem put forward by Volovik for isotropic electronic states to realistic electronic structures. In addition, the bound states due to the two superconducting bands of Pb can be separately detected and the two gaps close independently inside vortices. This yields strong evidence for a low interband coupling.


I. INTRODUCTION
The classical solutions in Ginzburg-Landau theory (GL) allow a thermodynamic classification of superconductors into type-I and type-II.Decisive for their behaviour in magnetic field is the interface energy between the superconducting and normal phase driven by the ratio of the London magnetic penetration depth λ L and the superconducting coherence length ξ.For Ginzburg-Landau parameters κ = λ L /ξ < 1/ √ 2, type-I behaviour and for > 1/ √ 2 type-II behaviour was predicted [1].A type-I superconductor is characterized by a positive interface energy and an attractive vortex-vortex interaction favouring an intermediate state with large normal domains [2].A type-II superconductor is characterized by a negative interface energy and a repulsive vortex-vortex interaction that favours an Abrikosov lattice of single-flux-quantum vortices in the mixed phase [3].
Flux quantization Φ = nΦ 0 in units of the flux quantum Φ 0 = h/2e is one of the most fundamental traits of the underlying off-diagonal long-range order of the superconducting condensate [4].n is the integer winding number of the vortex.This number of confined flux quanta is expected to affect the size and shape of the vortex.More quantitatively, Volovik demonstrated for superconductors with isotropic Fermi surfaces that one can determine n from an index theorem, provided one can measure the wave function of the quasiparticle bound states that form in the vortex core [5].The vorticity is particularly ambiguous, and hence interesting near the Bogomol'nyi point, κ = 1/ √ 2.Then, vortices near T c behave like non-interacting particles and the vortex configuration is infinitely degenerate [6][7][8][9] such that multiple-flux-quanta (or giant) vortices may be stabilized [10].This degeneracy is lifted below T c and a transitional region, in which the su-perconductor cannot be categorized into either of the types described above, emerges [11].When leaving the Ginzburg-Landau limit towards lower temperatures and especially, when considering multiple band superconductors, microscopic interactions are predicted to become increasingly important: the vortex-vortex interaction energy can become non-monotonic in distance through the existence of multiple, distinct superconducting coherence lengths ξ i [12,13] and topological hysteresis due to transitions between flux tubes and laminar pattern influences the flux patterns of the intermediate state [14].
In the past, the crucial parameter κ has been tuned toward the type-II regime by either an increase of λ L or reduction of the effective ξ, i.e. by using thin films below a critical thickness [15,16], by incorporation of impurities (for relevant experiments see Ref. [11] and references within) or by interface scattering, e.g.Pb/Si(111) [17].This approach, however, has the drawback that the quasiparticle bound states in the vortex core are considerably smeared out such that the index theorem cannot be applied [17][18][19][20][21].
Here, we take the alternative approach and study the twoband superconductor Pb in the clean limit in form of a bulk single crystal Pb(111) at 45 mK.In this respect, the multiple band superconductor Pb, which is closest to the Bogomol'nyi point of all elemental superconductors, is a good candidate to study the transitional phase at temperatures well below T c .We map the quasiparticle bound states for the two superconducting gaps of Pb with high energy and spatial resolution and use the index theorem including realistic band structures to determine the winding number of single-and multiple-flux-quanta vortices.Ultimately, we show that our observations allow to investigate the inter-band coupling of the two gaps.

A. Superconducting gaps and intermediate state
After several cycles of sputtering and annealing, we obtained a clean Pb(111) surface with wide terraces and monoatomic steps, as shown in Fig. 1(a).Upon zero-field cooling the Pb(111) sample to 45 mK it enters its superconducting state below T c ≈ 7.2 K [23].Due to our low electronic temperature of less than 100 mK [24] we are able to resolve the two gaps [25] in the density of states by scanning tunneling spectroscopy, even with a normal conducting tip.Fig. 1(c) shows the differential conductance in the superconducting state as black dots, including a temperature broadened two-gap BCS fit in orange.We determine the superconducting gaps to be ∆ 1 = (1.26± 0.02) meV (smaller gap) and ∆ 2 = (1.40 ± 0.02) meV (larger gap) in good agreement with previous measurement of the difference of the two gaps [25].The intensity difference of the two coherence peaks has previously been attributed to the k-dependent tunneling ma-trix elements and the larger gap has been assigned to the tubular Fermi surface sheet [25].This is in contrast to Bogoliubovde Gennes (BdG) based Korringa-Kohn-Rostoker (KKR) calculations [26], which deduced an opposite band-to-gap assignment.As will be discussed below, our study of the quasiparticle bound states in the vortices confirms the band-to-gap assignment of Saunderson et al. [26].We will from here on index the bands and Fermi surfaces according to their superconducting gap, i.e. the tubular Fermi surface responsible for ∆ 1 as Fermi surface 1 (FS 1) and the compact Fermi surface responsible for ∆ 2 as Fermi surface 2 (FS 2) (Fig. 1(c)).
After applying a perpendicular magnetic field of B = 85 mT, which is above the critical field µ 0 H c ∼ 80 mT [27], magnetic flux enters the sample from the sides and completely destroys superconductivity.Upon decreasing the field again below H c , the Landau intermediate state is reached.It is detected by recording the differential conductance at the coherence peak of ∆ 2 while ramping the field down.Once a jump to the superconducting state below the tip is detected, the ramp is stopped.This ensures that one typically finds both, super-conducting and normal conducting, areas in the scan range of the STM of 1.4 × 1.4 µm 2 .The intermediate state is characterized by large normal and superconducting domains.The shapes and sizes of these domains in the intermediate state of lead have been extensively studied by magneto-optical methods revealing the strong dependencies on temperature, sample shape and magnetic protocol [14,[28][29][30].
A typical domain wall in the intermediate state is shown in the dI/dU map in Fig. 1(b).At a tunneling bias of U t = 1.3 mV the normal conducting domains show up as areas of low conductance (purple) and the superconducting domains as areas of high conductance (green/yellow).Note that atomic step edges of the surface cause a contrast in dI/dU as typically found in STM experiments and is illustrated by the inset on the right.A cross-sectional line scan across the domain wall, as in Fig. 1(d), shows how both gaps change from zero to their maximum on the length scale of the coherence length.The local recovery of superconductivity agrees well with reported coherence lengths of ξ = 87 nm [23].For a detailed analysis of the coherence length ξ 1,2 of the two bands measured inside vortices, we refer to the next section.

B. Single-flux-quantum vortices
Inside the normal-core region of vortices, electronic bound states that lie within the superconducting gap are localized.These in-gap states were first theoretically studied by Caroli, de Gennes and Matricon (CdGM) in 1964 [31].The CdGM states for isotropic bands can be characterized by their orbital angular momentum number µ and their energy.The energy spacing of the discrete CdGM states is of the order of ∆ 2 /E F and the discrete states thus form a quasicontinuum or branch of CdGM states for most superconducting materials.The CdGM states with low µ are confined closer to the vortex core than the ones with high µ thus leading to a one-to-one correspondence between the angular momentum and the real-space behavior of the bound-state wave function.The higher the angular momentum µ or energy E(µ), the further away from the vortex center are the states localized [32].In case of a m-flux quanta vortex, m individual CdGM branches exist [5,[33][34][35][36] leading to a topological index theorem: The number of zeroenergy crossings of CdGM state branches with varying angular momentum is directly related to the vorticity [5,33,34].This theorem translates to real space and the number of zeroenergy crossings of CdGM state branches with varying radius from the vortex core is related to the vorticity.
STM allows to measure the variation of the local density of states (LDOS) inside the vortex and thus to determine the winding number of the vortex using the index theorem.In 1989, Hess et al. experimentally confirmed that CdGM states exist in single-flux-quantum vortices by scanning tunneling microscopy (STM) [37], but for vortices with multiple flux quanta, although studied in thin films with electron holography [38] and scanning Hall probe microscopy [39], their predicted bound states still lack experimental verification [35,36,40,41].
Using the detection method described in the previous sec-tion we are also able to find isolated, round normal conducting domains (appearance at eU = ∆ 2 ) of ≈ 100 nm in diameter.An example is shown in Fig. 2(h).As will be shown later, these are vortices in the superconductor with integer number of flux quanta.The finding of such small normal conducting domains is surprising considering that the domain wall energy in type-I superconductors is positive and the system thus tries to maximize its domain size.Even more surprising is the fact, that we found this shape after ramping the field down from the critical field.Magneto-optical measurements of cylindrical shaped intermediate state lead samples at 4.5 K reveal that normal domains are only tubular upon increasing magnetic field; after ramping down from the critical field the preferred structure is laminar [28].A deciding factor for the intermediate state domain structure on a microscopic scale could be the effect of flux branching [42][43][44].Since the overall domain structure in our experiment, however, consists of domains of various shapes and sizes, we argue that this finding is only consistent with circumstances under which vortices are weakly-interacting and have a non-monotonous interaction energy causing hysteresis, i.e. a superconductor in the transitional phase (close to the Bogomol'nyi point) at temperatures well below T c .The role of pinning of these vortices at bulk defects below the surface remains unclear.Repeating our magnetic protocol several times, the vortices may appear at similar positions suggestive of some kind of pinning.However, we also find the vortices to be mobile when varying the magnetic field (see next section).
To determine the amount of flux carried by the small normal domains, we record differential conductance (dI/dU ) maps at sub-gap energies, which essentially show the LDOS of the CdGM states.At zero bias voltage, we find a threefold symmetric state in form of a star with a maximum in the star's centre (Fig. 2(a)) inside the normal domain similar to starshaped CdGM states seen by Hess et al. in 2H-NbSe 2 [45].The quasiparticle density of states stretches over 100 nm in the 2 11 directions (blue/green).Additionally, three weak arms (dark blue) along the 1 21 directions are visible.With increasing energy (independent of sign of bias voltage) the star's arms split into two with increasing splitting distance, while the central peak splits nearly isotropically to a ring shape (Fig. 2(b-f)).For E ∆ 1 the strong arms are still visible and the ring reaches its maximal size (Fig. 2(g)).For E ∼ ∆ 2 the vortex shows as a relatively round area of low conductance with ≈ 100 nm in diameter (Fig. 2(h)).
For bands with anisotropic Fermi velocity, like in our case, the index theorem is not straight forward applicable since the radial symmetry is removed and a radial dependent measurement does not ensure crossing all diabolical points [5,33].Diabolical points here mean the points where the semiclassical particle and hole spectrum meet.The degenerate gapless fermionic excitations or "zero-modes" at these points carry the topological charge (singularity in the phase) and owe their name to the double-conical (diabolo) dispersion in parameter space [46].Here, a realistic band structure needs to be considered.We, thus, simulated the quasiparticle trajectories inside a vortex carrying one flux quantum within the quasiclassical Eilenberger theory including the Fermi velocity of Star shaped CdGM states have first been found by Hess et al. [45] in 2H-NbSe 2 , which due to the crystal symmetry have sixfold rotational symmetry.Pb crystallizes in the fcc structure and belongs to the point group F m3m.The electronic structure therefore only carries a threefold rotational symmetry about the [111] axis.Due to the discrete rotational symmetry, the angular momentum is not a good quantum number and the states mix.In general, the mixing results in states with different lateral confinement for different angles.At zero bias, the low angular momentum states carry the largest weight and lead to a maximum in the center of the vortex.At higher bias, states of larger angular momentum become increasingly im-portant leading to a movement of the maxima away from the center, i.e. the star splits.The large variation in quasiparticle localization depending on the angle from ∼ 10 nm to over 100 nm agrees with the strong anisotropy of the 3D Fermi velocity in the compact band.
In fact, star shaped vortices with long arms have been recently observed in the Abrikosov lattice of La(0001) [47], owing to the large anisotropy of the responsible band's in-plane Fermi velocity.In Pb with its two bands, one may expect that the individual gaps will close independently when approaching the vortex center.Fig. 3 displays the evolution of the gaps, i.e. the energy of the coherence peaks which appear as bright lines in the second derivative of dI/dU with respect to U , as function of the distance from the vortex center and direction.For large distances from the center, the gaps ∆ 1 and ∆ 2 decrease towards the center in parallel with roughly the same length scale ξ 1 ∼ ξ 2 ∼ 45 nm.In the whole range, ∆ 2 follows a simple tanh function.At 50 nm, however, ∆ 1 crosses ∆ 2 and stays larger than ∆ 2 .Thus, for the band 1, the gap size deviate from the tanh function near the center and decreases on a shorter length scale.This observation for a single band superconductor is known as the Kramer-Pesch effect [48].It was quantified by numerical calculations by Gygi and Schlüter [32] for the vortex core size of type-II superconductors.The slope of ∆ 1 near the core corresponds to a core size ξ dr −1 of only ∼ 10 nm.In our self-consistent calculation of the pair-potential ∆(r) for an isotropic vortex in the quasiclassical theory (see Supplementary Material), this Kramer-Pesch shrinking effect is also present and leads to substantial deviation from a tanh function with one universal ξ.Note that theory in the clean limit, however, predicts a core shrinking proportional to T /T c when lowering the temperature.At very low temperatures T T c the slope of the order parameter d∆/dr at the vortex centre is even predicted to become infinite, which would show as a jump of ∆(r) at r = 0 that is smoothed out over the distance ξT /T c [33].Experimentally, we do not find this extreme shrinking.The "squeezing" of low angular momentum states in the vortex core and thus the Kramer-Pesch effect is absent for the star, i.e.CdGM states of band 2. Consequently, ξ which indicates that intra-band coupling dominates over inter-band coupling, i.e. the two bands are sufficiently decoupled from each other, despite the gap sizes being not too different [49].For an isotropic two-band superconductor, the shrinking of the core region for only one of the bands has also been predicted in the case of rather weakly coupled bands [50].Again, this hints towards a low inter-band coupling in Pb.

C. Anomalous single-flux-quantum vortices
Besides the regular vortex situation, we find vortices, in which the two sets of CdGM states are laterally displaced.Figure 4 shows such an anomalous vortex.Both the ring centre and the star centre of the CdGM states are independently movable by a change in magnetic field and their relative dis-placement can be manipulated into different configurations, even back to the normal configuration from Figure 2 (see Supplementary Material).We explain this by two effects.First, a change in magnetic field laterally moves the vortex and with it, the two sets of CdGM states.Second, a change in magnetic field can lead to a tilting or bending of the flux lines away from a normal direction to the surface.This leads to a breaking of the cylinder symmetry and can displace the two sets of CdGM states relative to each other.The individual sets of CdGM states behave as those of the regular vortices, except for their relative displacement (see Figure 4 (a-h)).
Interestingly, the lateral displacement allows an independent probing of the state sets and thus, an independent identification of their bands.By looking at single bias spectra at the star centre (black) and the ring centre (red) in Figure 4(i) we find that the amplitude of the peak at zero bias is three times larger for the ring than for the star, which is supported by our separate band simulations from earlier and can be explained by the larger lateral confinement of low angular momentum states in band 1 compared to band 2. Figures 4(j,k) reveal maxima in the LDOS along the cross-sections marked in the right panel of (i) in the form of heat maps of the differential conductance's second derivative with respect to bias voltage.It becomes apparent that ∆ 2 (cyan line) closes entirely whereas ∆ 1 (red line) does not completely close in the star's centre (black circle).Instead, ∆ 1 (red line) closes about 20 nm away from the star centre (red circle).Consequently, the superconducting gap ∆ 1 is linked to the ring and ∆ 2 to the star.The fact that the quasiparticles of ∆ 1 and ∆ 2 can be independently displaced with respect to each other rules out rigid inter-band coupling.

D. Multiple-flux-quanta vortices
We also observed larger vortices in the experiments.Figure 5 displays two examples.(a-c) shows a vortex with two flux quanta.Our semiclassical calculations indicate that for each flux quantum and each band, a branch of CdGM states is present.As a result, for the vortex with two flux quanta, the band 2 causes a structure with two arms per direction at zero bias (see Fig. 5(a)) that individually split into two arms with bias voltage (see Fig. 5(b)) just as in the single-flux-quantum vortex.An analogous behaviour is observed in the vortex with three flux quanta shown in Fig. 5(d-f).Both vortices are larger than the single flux quantum vortex (compare Fig. 2(h) with Fig. 5(c,f)).Further, they deviate from a round shape and laterally grow with the number of flux quanta.The number of flux quanta in the vortices does not seem to be limited.We observed giant vortices with over 10 flux quanta (see Supplementary Material).For band 1, the problem of multiple-fluxquanta vortices is similar to that of a single band superconductor with a spherical Fermi surface and has been studied in detail by Volovik et al. [5,33].In essence, an m-flux quanta vortex results in m branches of circular CdGM states of different radii.Due to symmetry, at zero bias and for odd m, a central spot is formed by one branch of the CdGM states and the other states form pairs of increasing ring diameters.The CdGM states at energies away from the Fermi level evolve by changing the radii.Thus, the central CdGM state turns from a spot of zero radius to a ring of finite radius, and the pairs of CdGM states split in their radius.For an even m, no central spot is present at the Fermi energy but only pairs of CdGM states with distinct radii exist, that again split when moving away from the Fermi energy.This is the essence of the topological index number theory of Volovik.In our case, we find a central spot of high differential conductance at zero bias voltage for the vortices with m = 1 and m = 3, while there is no central spot for the vortex with m = 2 in agreement with the topological index number theory.When going away from the Fermi energy, the ring-shaped and spot-like CdGM states overlap with the star shaped states such that their distinction becomes impractical.Tunneling spectra at selected locations inside the vortex showing the very same effects are shown in Figure 5(g,h).
While the original formulation of Volovik's index number theorem is not directly applicable in real space (that is counting the number of zero modes in the vortex core as function of radius) in our case of an anisotropic Fermi surface, we still believe that the number of parallel arms of our star-shaped CdGM state pattern at zero energy exactly reflects the winding number m of the vortex.The reason for our confidence lies in the fact that the diabolical points in k-space, that need to be counted, still lie at zero energy, it is just not obvious which path in real space crosses each of them exactly once.Comparing the axisymmetric problem with the problem at hand, it becomes obvious that a parametrization in terms of the set of quantum numbers (k r , k φ , k z ) has to be replaced by an irreducible representation of the crystallographic point group.The flat parts of Fermi surface 2 focus the quasiparticles into high-symmetry directions and instead of a localization of different CdGM state branches at different radii, they now localize at different impact parameters.
At last, we tested whether vortices would also be present at temperatures much closer to T c .Our measurements at 4.3 K revealed that it is substantially harder to trap a vortex in our scan frame, but we managed to in a single case (see Supplementary Material).

III. CONCLUSIONS
In conclusion, we report the observation of single-fluxquantum and multiple-flux-quanta vortices in a traditional type-I bulk superconductor, i.e. single crystal Pb(111), by low-temperature scanning tunneling microscopy.We demonstrate a robust determination method for the winding number of the vortex by usage of a topological index theorem which relates the number of flux quanta and CdGM state branches.The spatial anisotropy of the quasiparticle states inside the vortex reflects the crystal symmetry, and its shape is governed by the anisotropic Fermi velocity in the superconducting bands.An influence of neighboring flux lines can be ruled out due to the absence of an ordered flux pattern.In addition, we could show how CdGM states from two weakly coupled bands interact with each other in a single flux line.While the generically unavoidable inter-band coupling leads to the formation of a single coherence length and order parameter right at T c [51], the physics below T c is considerably richer, in particular in the regime of weak inter-band coupling [49,50].Of particular interest is the emergence of mixed gap modes [50] or the possibility of the observation of the rather elusive Leggett mode in the limit of weak inter-band coupling [52], i.e. a collective fluctuation of the relative phase of the two bands, for superconducting Pb.We experimentally demonstrate that for systems with multi-band pairing, the transitional region of the hysteretic behaviour expands due to microscopic interactions and also allows for interesting vortex configu- rations, like vortex clusters or multiple-flux-quanta vortices [53].The emergence of such vortices in the intermediate state of a prototypical type-I superconductor is very surprising and demonstrates that the low temperature flux patterns of twoband superconductors rank on a spectrum between type I and II.Not only does the existence of vortices in Pb allow for an alternative test ground of vortex physics in multi-component superconducting systems apart from MgB 2 , the variability of the vortex winding number also suggests a combination with topological crystal defects that was recently predicted to result in topological quasiparticle states, like Majorana zero modes, under certain circumstances [54,55].

METHODS
a. Experimental details : The experiments were performed with a home-built scanning tunneling microscope with dilution refrigeration, which can reach a base temperature of 25 mK in a magnetic field of up to 7.5 T [24].In our setup, the bias voltage U t is applied between sample and common machine ground so that a positive bias voltage probes the unoccupied states of the sample.The STM chamber is kept at a base pressure of 1 × 10 −10 mbar.The single crystal Pb(111) (miscut angle: ±0.1 • , purity: 99.999%) has been purchased from MaTecK GmbH.It has cylindrical (hat) shape with a diameter of 8 mm and a thickness of 2 mm.At a base pressure of 1 × 10 −10 mbar the Pb crystal was prepared in cycles of sputtering with Ar + ions of 3 keV and subsequent anneal- were all performed below 45 mK.After zero-field cooling the Pb crystal, the vortices were formed by ramping the perpendicular magnetic field from 0 mT to 85 mT and back down to a constant value below the critical field.Note that 80 mT is the critical magnetic field for Pb [27].The differential conductance was measured using a Lock-in amplifier at a frequency of 3.4 − 3.6 kHz and AC peak amplitude U PK ac between 10 and 100 µV (for details, see Supplementary Material).dI/dU maps at sub-gap energies were recorded in a multi-pass configuration: In the "record" phase the tip records the z-profile at the feedback condition (constant tunneling current of I t at eU t > ∆) and in the "play" phase the z-profile is repeated at a different bias voltage.In order to increase the signal for subgap energies an offset is often added to the z-profile bringing the tip closer to the surface.This record and play phase alternation is performed line for line until the entire area has been scanned.
b. Calculation details : In order to obtain the simulated LDOS for vortices containing an arbitrary number of flux quanta, we used the Riccati parametrisation of the quasiclassical 3D Eilenberger equations as proposed in Ref. [56] and numerically solved the one-dimensional differential equations under appropriate boundary conditions (for details, see Sup-plementary Material) for each band separately.Motivated by our experimental findings, we assumed a radial symmetric local pair potential ∆(r) in the plane perpendicular to a vortex line with s-wave symmetry that vanishes in the vortex centre.We used the ratio ∆ 2 (∞)/∆ 1 (∞) obtained from the experiment.We respected the broken translation symmetry at the crystal surface by a work function term.The magnetic vector potential was set to zero for all calculations shown in the main text.An inclusion of a magnetic vector potential of appropriate form only yielded small quantitative deviations from the zero-field case (see Supplementary Material), yet drastically increased the required computation time, which is why we refrained from it for the LDOS maps.
Density functional calculations of the electronic structure of Pb were carried out in the framework of the mixed-basis pseudopotential method [57].The electron-ion interaction was represented by norm-conserving relativistic pseudopotentials [58].Spin-orbit coupling was incorporated within the pseudopotential scheme via Kleinman's formulation [59] and was consistently taken into account in the charge selfconsistency cycle using a spinor representation of the wave functions.Further details of the spin-orbit coupling implementation within the mixed-basis pseudopotential method can be found in a previous publication [60].For higher accuracy, 5d semicore states were included in the valence space.The deep d potential is efficiently treated by the mixed-basis approach, where valence states are expanded in a combination of plane waves and local functions.Here, local functions of d type at the Pb sites were combined with plane waves up to a kinetic energy of 20 Ry.Brillouin-zone integration was performed by sampling a 32×32×32 k-point mesh (corresponding to 2992 k points in the irreducible part of the Brillouin zone) in conjunction with a Gaussian broadening of 0.2 eV.The exchange-correlation functional was represented by the local-density approximation in the parameterization of Perdew-Wang [61].
This DFT technique was applied to obtain Fermi surface properties entering the Eilenberger equations.Band energies were calculated on fine radial grids for a cylindrical coordinate system taking the [111] direction as the z-axis, to determine Fermi momenta k F for each of the two relevant bands.At each k F , 3-dimensional Fermi velocities v F were then calculated taking numerical derivatives of band energies around this point.The optimized lattice parameter a = 4.89 Å was used throughout.TABLE S1.Measurement parameters.I t denotes the feedback condition of the tunnelling current at the bias voltage U t .U PK ac is the peak AC voltage amplitude added via the Lock-in amplifier.In multi-pass dI/dU maps z off is the distance the tip is brought closer to the surface compared to the feedback condition.
B is the magnetic flux density in vacuum and T is the temperature.I ensures, that quasiparticles travelling to the surface are decaying into the vacuum with the right damping factor.Since an isotropic ∆ implies an isotropic in-plane current density, the magnetic field profile around the vortex was described by a vector potential of the form [4] that is cylinder symmetric in the bulk and deviates from the bulk value near and above the surface.
Here, κ = √ k 2 + λ −2 , λ is the magnetic penetration depth which was chosen to be λ = ξ/ √ 2 and J 1 (x) is a Bessel function of first order.
Nils Schopohl and Kazumi Maki [5] could show that the Eilenberger equations can always be solved along a characteristic line and that the solution is universal for each point along this line.
The line simply has to be parallel to the Fermi velocity vector v F .Points on this line are then characterized by the variable X and two impact parameters Y and Z.
= x x + y ŷ + z ẑ Transformation from the mobile frame of reference ( û, v, ŵ), where û v F , to the fixed coordinate system ( x, ŷ, ẑ) is done by chaining rotation matrices using Euler angles: In order to align the axis û of the mobile frame with the axis x of the static coordinate system, no rotation about ŵ is needed and thus ψ = 0.The other two rotational angles are defined by the components of the Fermi velocity vector as follows: Using this transformation, the vector potential and gap parameter can be expressed in the variables of the mobile frame: A(r(X, Y, Z)) and ∆(r(X, Y, Z)).On a characteristic line defined by the impact parameters Y p and Z p , the functions from Eq. ( 1) become [1] ∆ The parametrisation of the Eilenberger equations on this 1D line is called the Riccati parametrisation.Eq. ( 1) formally reduces to the solution of two initial value problems, where the differential equations are scalar and of first order.The Eilenberger propagator is parametrised in terms of two scalar complex functions a(X) and b(X): The differential equations that need to be numerically solved for a(X) and b(X) are [5] v with boundary conditions In order to solve them at a certain energy E < ∆ 0 , the analytical continuation i n → E + i0 + was used.Finally, the local density of states at r was obtained from the real part of g(r(X)), i.e.
In order to obtain the total local density of states one still needs an integration over the Fermi surface to calculate all trajectories that are present due to the various Fermi velocity vectors that exist, plus an integration over the impact factors Y in order to account for trajectories that do not traverse the vortex centre.An integration over Z is redundant since the solution in the bulk is the same for every Z.Near the surface this is strictly not the case anymore but the effect is small and with a large enough variety of Fermi velocity vectors (which we have) these trajectories are not missed.
c. Influence of vector potential The inclusion of a non-zero vector potential A in the calculations increases the average time needed to solve the Eilenberger equations because the Matsubara frequencies in Eq. ( 14) gain a position dependent term that requires a coordinate transformation between the two reference frames.Therefore, the simulations shown in the main text are performed without vector potential.It is already visible from Eq. ( 14) that in a calculation where A and ∆ are not solved self-consistently, A only enters the equation as an effective energy term.
With a vector potential in the azimuthal direction like in Eq. ( 5) the scalar product with the Fermi velocity is only expected to yield significant contribution for large impact parameters (for Y p = 0, v F is perpendicular to A).That means trajectories with increasing impact parameter have large LDOS already for smaller distances than in the field free case.The splitting star arms in the LDOS maps should be squeezed to smaller distances from the core.In fact, this is what is seen in the calculations with vector potential at higher energies, as shown in Fig. S1.This proves that even though there are quantitative differences to the case without vector potential, in the general characteristics, the LDOS patterns remain unchanged.

FIG. 1 .
FIG. 1. Superconducting properties and intermediate state.(a) Topographic scan image of the Pb(111) surface after the cleaning procedure.(b) dI/dU map at Ut = 1.3 meV showing a typical domain wall in the intermediate state at B = 23 mT.Right overlay: Corresponding topographic scan image.(c) Differential conductance in the superconducting state (or superconducting domain) (black points) including a two-gap BCS fit (orange) with ∆1 being the smaller and ∆2 being the larger gap.Inset: 3D models of Fermi surface sheets 1 and 2 (from Ref. [22]) responsible for superconducting gaps ∆1 and ∆2 respectively.(d) dI/dU spectra along a cross-section from normal conducting to superconducting domain.The area is marked in white in (a) and the profile direction is indicated by a red arrow.Individual spectra are offset with respect to each other by 0.7 µS The spectra were locally averaged over a straight part of the domain boundary and recorded in distance increments of ∆d = 9.19 nm.

FIG. 2 .
FIG. 2. Single-flux-quantum vortex signature.(a-h) dI/dU maps (B = 0 mT) of a single-flux-quantum vortex at different bias voltages displaying the quasiparticle density at the surface.Right panels in (a) and (e) show the simulated LDOS of the bands 1 and 2 for the respective energy.The inset in (a) shows a 2D-FFT filtered topographic image of the Pb(111) lattice.Bulk crystal directions (red/yellow arrows) have been determined from glide planes.

FIG. 3 .
FIG. 3. Coherence lengths in the normal single-flux-quantum vortex.Angle dependent radial measurement of the vortex core states from the vortex in Fig. 2. Displayed is the second derivative of the differential conductance −∂ 2 σ/∂U 2 in order to highlight maxima in the LDOS.The point distance of single spectra is ∆d = 2.57 nm.Marked are the opening of ∆1 (red line), the opening of ∆2 (cyan line) and the CdGM states of band 2 (cyan dashed line).The inset shows the direction of the y-axis.

FIG. 4 .
FIG. 4. Anomalous vortex signature.(a-h) dI/dU maps (B = 19 mT) at different bias voltage for an anomalous vortex (compared to the normal vortex configuration from Fig. 2).(i) Right panel: Enlarged zero bias dI/dU map of the anomalous vortex including the position of line spectra (white arrows) and single spectra locations (red/black circle).Left Panel: Single dI/dU spectra in the star centre (black) and ring centre (red) revealing ZBPs of different amplitude.(j,k) Heat maps of the second derivative of the differential conductance −∂ 2 σ/∂U 2 along the cross-sections marked in (i).Red and cyan lines follow ∆1 and ∆2 respectively.

FIG. 5 .
FIG. 5. Multiple-flux-quanta vortex signature.(a-f) dI/dU maps of an m = 2 vortex at B = 0 mT (a-c) and an m = 3 vortex at B = 33 mT (d-f) at selected bias voltages.Right panels in (a) and (d) show the simulated LDOS of the bands 1 and 2 for the respective m quanta vortex and energy.(g,h) Bias spectroscopies (bottom panel) recorded with the tip at specific locations marked in the zero bias maps (top panel) of the m = 2 (g) and m = 3 (h) vortex.Individual spectra are offset by 0.4 µS for clarity.Black dashed lines indicate their respective zero conductance baselines.

FIG. S1 .
FIG. S1.Influence of Vector Potential: The LDOS obtained from solutions of the 3D Eilenberger equations for a single-flux-quantum vortex at energy eU = 0.8 meV without (a) and with vector potential (b),as formulated in Eq. (5), show only quantitative differences.With non-zero vector field, the star arms are still split, yet the CdGM states at this energy are squeezed into a smaller area around the vortex core.