Supplementary Materials forBand-resolved Caroli–de Gennes–Matricon states of multiple-flux-quantavortices in a multiband superconductor
Thomas Gozlinski et al.
Corresponding author: Thomas Gozlinski, thomas.gozlinski@kit.edu
Sci. Adv. 9, eadh9163 (2023)
DOI: 10.1126/sciadv.adh9163

The PDF file includes:
Extended Methods
Supplementary Notes 1 to 5
Supporting Data 1 and 2
Table S1
Figs. S1 to S10
Legends for data S1 to S6
References
Other Supplementary Material for this manuscript includes the following:
Data S1 to S6

Contents
Extended Methods
Supplementary Note 1: Anomalous vortices in varying magnetic fields
Supplementary Note 2: Giant vortex
Supplementary Note 3: Particle-hole and time-reversal symmetry
Supplementary Note 4: Vortex at 4.3 K
Supplementary Note 5: Self-consistent calculation of ∆Supporting Data 1: Absence of ordered vortex pattern
Supporting Data 2: dI/dU data for Fig. 3 and 4

Extended methods
STM measurement conditions: Tab. S1 summarizes the most important STM/STS mea-

surement parameters used to obtain the figures in the main text and supplementary material.

3D Eilenberger calculations : We follow Ref. (57) in notation and describe how we solved
the quasiclassical Eilenberger equations (65, 66)−ℏvF ∇ĝ(r; pF , iϵn ) =
iϵn + vF eA(r)−∆(r, pF ), ĝ(r; pF , iϵn ) ,∆† (r, pF )−iϵn − vF eA(r)(1)

that hold for kF ξ ≫ 1. The quasiclassical Green’s function propagator

ĝ =


g(r, pF , iϵn )f (r, pF , iϵn ),−f ∗ (r, −pF , iϵn ) g ∗ (r, −pF , iϵn )(2)depends on spatial coordinate r, crystal momentum pF = ℏkF and energy ϵn , and must satisfy
the normalization condition ĝ 2 = 1̂. g and f are normal and anomalous quasiclassical Green’s
function propagators and ϵn are fermionic Matsubara frequencies. Incorporated in this formalism is the self-consistent calculation of pair potential ∆(r) and A(r), which is essential for1the description of a superconductor hosting vortices as these two fields are dependent on eachother. Although this approach is numerically much more feasible than direct diagonalizationof a BdG Hamiltonian when treating inhomegeneties, it is not necessary to solve the problem
completely self-consistently in this work. The reason for that is, that our experimental results
clearly show, that even though the sub-gap states in the vortex display highly anisotropic behaviour, the recovery of ∆(r) is isotropic in the plane around the vortex core. The local pairpotential is assumed to have s-wave symmetry and was therefore modeled by∆(r, pF ) = ∆(r) =rz∆0 tanh + Θ(z)W tanhξa

x + iy
r

m
,

(3)

with Θ being a Heaviside step function, ∆0 the maximum gap size, ξ the coherence length, athe lattice constant, m the winding number of the vortex and W the work function. The second
term ensures, that quasiparticles travelling to the surface are decaying into the vacuum with theright damping factor. Since an isotropic ∆ implies an isotropic in-plane current density, themagnetic field profile around the vortex was described by a vector potential of the form (67)A = A(r, z)êφ ,Z∞
mΦ0
J1 (kr)A(r, z) =dk 2
S(k, z) ,22πλ
k + λ−2
0
(
κ
e−kz
z>0
S(k, z) = k+κ k κz
1 − k+κ ez ≤ 0,

(4)
(5)

(6)

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 J1 (x) is a Bessel function of first order.
Nils Schopohl and Kazumi Maki (68) could show that the Eilenberger equations can always
be solved along a characteristic line and that the solution is universal for each point along thisline. The line simply has to be parallel to the Fermi velocity vector vF . Points on this line are2then characterized by the variable X and two impact parameters Y and Z.
r(X) = X û + Y v̂ + Z ŵ,(7)= xx̂ + y ŷ + z ẑ(8)Transformation from the mobile frame of reference (û, v̂, ŵ), where û ∥ vF , to the fixed
coordinate system (x̂, ŷ, ẑ) is done by chaining rotation matrices using Euler angles:
 
 
vx
vX
vy  = Rzyz  vY 
vz
vZ
Rzyz = Rz (η)Ry (χ − π/2)Rz (ψ)




cos η − sin η 0
sin χ 0 cos χ
cos ψ − sin ψ 01
0   sin ψ cos ψ 0 .
=  sin η cos η 0  0
001− cos χ 0 sin χ
0
0
1

(9)

(10)

In order to align the axis û of the mobile frame with the axis x̂ of the static coordinate system,
no rotation about ŵ is needed and thus ψ = 0. The other two rotational angles are defined by
the components of the Fermi velocity vector as follows:vzη = arccos
,
|v|

vy

vx > 0arctan vx ,vy
χ = arctan vx + π , vx < 0 .
sgn(vy ) π2 ,vx = 0(11)(12)

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 definedby the impact parameters Yp and Zp, the functions from Eq. (1) become (57)∆(X) = ∆(r(X, Yp , Zp )) ,(13)iϵ̃n (X) = iϵn + vF eA(r(X, Yp , Zp )) ,(14)

ĝ(X) = ĝ(r(X, Yp , Zp ), pF , iϵn ).

(15)

3

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 interms of two scalar complex functions a(X) and b(X):1
1 − a(X)b(X)
2ia(X)ĝ(X) =.
−2ib(X)−1 + a(X)b(X)
1 + a(X)b(X)(16)

The differential equations that need to be numerically solved for a(X) and b(X) are (68)ℏvF a′ (X) + [2ϵ̃n + ∆† (X)a(X)]a(X) − ∆(X) = 0 ,(17)ℏvF b′ (X) − [2ϵ̃n + ∆(X)b(X)]b(X) + ∆† (X) = 0(18)

with boundary conditions
∆(−∞)
p
,
ϵn + ϵ2n + |∆(−∞)|2∆† (+∞)
p.
b(+∞) =
ϵn + ϵ2n + |∆(+∞)|2a(−∞) =

(19)
(20)

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.
N (pF ) = N0 (pF )Re

1 − a(X)b(X)
1 + a(X)b(X)


.

(21)

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 thatexist, plus an integration over the impact factors Y in order to account for trajectories that do nottraverse 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.
4Influence of vector potential: The inclusion of a non-zero vector potential A in the cal-culations increases the average time needed to solve the Eilenberger equations because theMatsubara frequencies in Eq. (14) gain a position dependent term that requires a coordinatetransformation between the two reference frames. Therefore, the simulations shown in themain 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 aneffective 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 a substantial contribution for large impact parameters (for Yp = 0, vF 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 shownin 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.

Supplementary Note 1: Anomalous vortices in varying magnetic fields
The direction and magnitude of spatial displacement between the sets of CdGM state branchesfrom the two superconducting bands found in anomalous vortices could not be tied to any
crystal direction and appears to be almost random. By slowly varying the static magnetic field,
the displacement can be modified, as shown in the dI/dU maps of the same anomalous vortex in
Fig. S2 at different fields. We do, however, see an additional blue stripe in the zero-bias dI/dU
maps of anomalous vortices that appears on the same side of the star pattern as the ring state.
For clarity, the ring state is enclosed by a white-dashed circle and the borders of the mentioned
stipe are highlighted by white-dashed contour lines in Fig. S2. These stripes might hint at the5direction in which the magnetic field lines are tilted beneath the surface.
As Fig. S3 demonstrates, a vortex core can be moved by a slow change in magnetic field
and the sets of CdGM states from the two superconducting bands are displaced independently.
Consequently, we are able to transform an anomalous vortex back to a normal vortex provingthat the two types are not inherently different and that the anomalous type may well be explainedby flux line tilting.

Supplementary Note 2: Giant vortex
In a single case, a giant vortex containing m > 10 flux quanta was found near a large sputtering
defect. The zero-bias dI/dU map of this vortex, shown in Fig. S4A, features more than 10 armsof CdGM states stretching in each ⟨21̄1̄⟩ direction, which again split into pairs at larger bias
(Fig. S4B). Unfortunately this particular giant vortex was partly outside the scan frame of thefine motion piezo tube, which prevented us from a definite determination of m.

Supplementary Note 3: Particle-hole and time-reversal symmetry
We find that the differential conductance maps show no qualitative change (compared to thefigures in the main text) when the electric field in the tunnelling junction is reversed (inverse bias
voltage sign U ↔ −U ) or the magnetic field direction is reversed (B êz ↔ −B êz ). As a point of
proof, Fig. S5A and B show a normal vortex stabilized at −14 mT (after saturation at −85 mT)
for U ≥ 0 mV and Fig. S5C and D show an anomalous vortex stabilized at 18 mT for U ≤0 mV. The invariance of the LDOS maps under reversal of the electric field demonstrates theparticle-hole symmetry of the CdGM states and justifies their treatment as excitations of a BCS
ground state, i.e. describing their dynamics by a mean-field BdG Hamiltonian or Eilenberger’s
quasiclassical Green’s function propagators. The invariance of the LDOS maps under reversal
6

of the magnetic field is the logical consequence of time-reversal symmetry of the CdGM states.

Supplementary Note 4: Vortex at 4.3 K
At 4.3 K the stabilization of single vortices within the STM scan frame was substantially harder.
Upon repeating the magnetic protocol in small steps seven times, we only found a vortex within
our scan frame in a single case. The vortex signature is shown in Fig. S6. At U = 1.8 mV ≈∆/e (A), it appears as a round depression in differential conductance that is slightly smaller
than at base temperature which is to be expected due to the temperature dependence of thepcoherence length ξ(T ) = ξ0 1 − T /Tc . At zero bias voltage (B), the CdGM states are smearedout. However, a conductance maximum is still found in the centre of the vortex (C) indicatingthat it is a vortex with odd winding number and, considering its size and shape, most likely
m = 1. A distinction between ∆1 and ∆2 is not possible anymore at this temperature due to thetemperature broadening of the dI/dU spectra (C).

Supplementary Note 5: Self-consistent calculation of ∆We followed Ref. (19) and solved the Eilenberger equations in 2D self-consistently for T =0.1/7.2 Tc in the clean limit, i.e. we refined the pair potential using the self-consistency equation(19, 65):∆(r) = 2πT Λ

Nc
X

SF−1

I

f (iϵn , r, vF (k))d2 k

(22)

FS

n=0

with
"
Λ = log(T /Tc ) +NcX1
n + 1/2
n=0

#−1
.

(23)

SF is the Fermi surface area, ϵn = 2(n+1)πT the fermionic Matsubara frequencies and Λ is thecoupling constant. We choose Nc such that ϵNc = 5 Tc . Setting the vector potential A = 0 andassuming an isotropic vF we reach convergence for ∆ with an accuracy of 10−4 Tc . The result
7

for ∆(r) is displayed in Fig. S7 along with the function ∆(r) = ∆0 tanh r/ξ0 . For a distancefrom the vortex centre of r < 2ξ0 the recovery of ∆(r) deviates from the function ∆0 tanh r/ξ0 .
The slope of ∆(r) close to r = 0 is roughly four times as steep as expected from a simple tanh
h
i−1
d∆(r)
(c)
behaviour with universal ξ0 . This leads to a core size ξ = ∆(∞) limr→0 drthat isroughly ξ0 /4.

Supporting Data 1: Absence of ordered vortex pattern
We provide supporting data that clearly shows that our bulk Pb single crystal does not host an
Abrikosov lattice (or any other kind of ordered pattern) of vortices, but that vortices coexist
with large superconducting and normal domains in the intermediate state. In addition, we provide topographic images that show the atomically flat Pb surface below these vortices with noindication for flux pinning by defects.
Fig. S8 shows three connected areas for which we show the topographic scan image in A andthe corresponding dI/dU map recorded at Ut = 1.4 mV in B. The topographic image shows atypical Pb surface after our cleaning procedure. Here, we see six atomically flat and clean
terraces (a seventh terrace begins just at the left edge of the leftmost scan frame) without larger
defects at the surface. In the corresponding dI/dU map we see two vortices (round blue/purple
objects marked by arrows), one larger normal conducting domain (purple area marked as “NC”)
and a large superconducting domain (green/red area). The vortices obviously do not form anordered pattern but coexist with large normal and superconducting domains in the intermediate
state.

Supporting Data 2: dI/dU data for Fig. 3 and 4
We provide supporting data that shows the differential conductance (dI/dU ) spectra from which
Fig. 3 and Fig. 4J-K of the main text were derived. In Fig. 3 and Fig. 4J-K of the main text the8second derivative of dI/dU with respect to the bias voltage is shown. This method enhances
visibility of peaks in the spectra and is commonly used for photoemission data. We now also
show dI/dU as a function of bias voltage and the radial distance to the vortex centre as a 3D
surface plot in Fig. S9 for the four different angles (0◦ , 20◦ , 40◦ , 60◦ ) which were defined in themain text and Fig. 3. The same we do for the two plots of Fig. 4J-K of the main text and showthese 3D surface plots in Fig. S10. Since it is hard to follow the CdGM states as well as thegaps in one view, we also provide these interactive 3D surface plots in .html format in 3D-Data
S1 to S6.

9

Table S1: Measurement parameters. It denotes the feedback condition of the tunnellingPK
is the peak AC voltage amplitude added via the Lock-in
current at the bias voltage Ut . Uac
amplifier. In multi-pass dI/dU maps zoff 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 thetemperature.
Fig.
1

2
3
4

5

S2
S3
S4
S5
S6

A
B
C
D
A-G
H
A-G
H
I
J
K
A-B
C
D-E
F
G
H
A-D
A-B
C-D
A-B
A-B
C-D
A
B
C

PK
It (nA) Ut (mV) Uac(µV) zoff (pm) B(mT) T (mK)
0.2
1
< 45
0.05
1.3
50
23
< 45
1
3
20
0
39
2
3
20
23
48
0.5
1.4
50
20
0
< 45
0.5
1.4
50
0
< 45
0.5
3
20
0
< 45
1
1.4
50
20
19
< 45
1
1.4
50
19
< 45
1
3
20
19
40
1
1.5
10
19
40
1
2
20
19
40
0.5
1.4
50
20
0
< 45
0.5
1.4
50
0
< 45
1
1.4
50
20
33
< 45
1
1.4
50
33
< 45
1
3
20
0
< 45
1
3
20
33
< 45
1
1.4
50
INDa
< 45
0.2
1.4
50
IND
< 45
1
1.4
20
IND
< 45
1
3
50
20
0
< 45
0.5
1.4
50
IND
36
1
1.4
20
IND
42
0.1
1.8
50
0
4300
0.1
1.8
50
20
0
4300
0.1
6
100
0
4300
a

This parameter is indicated in the sub-figure itself.

10

A

B

Figure 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) andwith 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.

11

A

B

C

D

Figure S2: Manipulation of anomalous vortex pattern. Zero bias dI/dU maps of an anomalous vortex at different magnetic fields. The relative positional shift between the star centre andring centre (red) is manifold and not tied to crystal directions. The position of the ring centre ismarked by a white-dashed circle. The bright stripe mentioned in the main text is marked at its
borders by white dashed lines.

12

A

C

B

D

Figure S3: Manipulation from anomalous to normal vortex. (A,C) Topographic imagesof the corresponding quasiparticle patterns in (B,D) showing the same location on the surface.
(B,D) Zero bias dI/dU map of the anomalous vortex before (B) and after the magnetic field ramp
(D). After the decrease of magnetic field, the anomalous vortex in (B) moved and transformed
into a normal type (D). In the process, both sets of CdGM states moved.

13

A

B

Figure S4: Giant vortex. dI/dU maps of a giant vortex containing m > 10 flux quanta. It
shows more than 10 arms at zero bias (A) that individually split in two at sub-gap energies
away from the Fermi level (B). This vortex is located at the edge of a large sputtering defect.

14

A

B

C

D

Figure S5: Reversal of magnetic and electric field. Applying the magnetic field or electricfield in the opposite direction has no effect on the LDOS pattern inside the vortex. (A,B) dI/dU
maps of a normal single-flux-quantum vortex in reversed magnetic field. (C,D) dI/dU maps ofan anomalous vortex showing that reversing the sample bias yields indentical LDOS patterns.

15

A

B
C

Figure S6: Vortex at higher temperature. (A,B) dI/dU maps of a vortex at T = 4.3 K and
B = 0 mT. (C) Bias spectroscopies far away from the vortex (black squares) and in the centreof the vortex (blue triangles) at 4.3 K. The red line shows a fit of Dynes form to the differential
conductivity in the superconducting phase with a gap size of ∆ = 1.24 meV.

16

Figure S7: Kramer-Pesch effect: The self-consistent calculation of the pair potential ∆ fora vortex with isotropic Fermi velocity at T = 0.1/7.2 Tc exhibits a shrinking of its core size
according to the Kramer-Pesch effect, i.e. a steeper recovery of ∆ close to the vortex centre that
does not follow tanh r/ξ0 .

17

Figure S8: Absence of ordered vortex pattern: (A) Topographic image of three connectedareas showing six atomically flat terraces. (B) dI/dU maps at Ut = 1.4 mV and B = 22 mT
corresponding to the areas in (A) which show 2 vortices (round blue/purple objects marked byarrows) next to a larger normal conducting domain (purple area marked as “NC”). The green/red
areas in the map belong to the superconducting parts (marked as “SC”).

18

Figure S9: dI/dU spectra to Fig. 3: 3D surface plots of dI/dU (r, U ) for the four different
angles as defined in Fig. 3 (0◦ , 20◦ , 40◦ , 60◦ ). Bias voltage U is shown as x-coordinate, radial
distance to the vortex centre r as y-coordinate and dI/dU as color coded z-coordinate. The
color encoding changes exponentially over the equally spaced values of the z-axis in order toemphasize small in-gap peaks corresponding to CdGM states and coherence peaks of larger
intensity at the same time.

19

Figure S10: dI/dU spectra to Fig. 4: 3D surface plots of dI/dU (r, U ) for the two line profiles
defined in Fig. 4I (differential form of spectra shown in Fig. 4J-K). Bias voltage U is shown as
x-coordinate, radial distance to the vortex centre r as y-coordinate and dI/dU as color coded
z-coordinate. For large range spectra the color encoding changes exponentially over the equally
spaced values of the z-axis in order to emphasize small in-gap peaks corresponding to CdGM
states and coherence peaks of larger intensity at the same time. (A) dI/dU (r, U ) spectra to
Fig. 4J. (B) dI/dU (r, U ) spectra to Fig. 4K.

20

3D-Data S1: 3D data for Fig. 3. dI/dU as a function of bias voltage and the radial distance tothe vortex centre as a 3D surface plot for 0◦ angle as defined in the main text in Fig. 3.
3D-Data S2: 3D data for Fig. 3. dI/dU as a function of bias voltage and the radial distance tothe vortex centre as a 3D surface plot for 20◦ angle as defined in the main text in Fig. 3.
3D-Data S3: 3D data for Fig. 3. dI/dU as a function of bias voltage and the radial distance tothe vortex centre as a 3D surface plot for 40◦ angle as defined in the main text in Fig. 3.
3D-Data S4: 3D data for Fig. 3. dI/dU as a function of bias voltage and the radial distance tothe vortex centre as a 3D surface plot for 60◦ angle as defined in the main text in Fig. 3.
3D-Data S5: 3D data for Fig. 4J. dI/dU as a function of bias voltage and the radial distanceto the vortex centre as a 3D surface plot for the profile in Fig. 4J of the main text.
3D-Data S4: 3D data for Fig. 4K. dI/dU as a function of bias voltage and the radial distanceto the vortex centre as a 3D surface plot for the profile in Fig. 4K of the main text.

21

REFERENCES AND NOTES
1. V. L. Ginzburg, L. D. Landau, On the theory of superconductivity. Zh. Eksp. Teor. Fiz. 20, 1064–1082(1950).
2. L. D. Landau, Theory of supraconductivity. Zh. Eksp. Teor. Fiz. 7, 371 (1937).
3. A. A. Abrikosov, On the magnetic properties of superconductors of the second group. Sov. Phys. JETP
5, 1174 (1957).
4. C. N. Yang, Concept of off-diagonal long-range order and the quantum phases of liquid He and of
superconductors. Rev. Mod. Phys. 34, 694–704 (1962).
5. G. E. Volovik, Localized fermions on quantized vortices in superfluid 3He-B. J. Phys. Condens. Matter
3, 357–368 (1991).
6. E. B. Bogomol’nyi, A. I. Vainshtein, Stability of strings in gauge Abelian theory. Sov. J. Nucl. Phys.
23:5, 7327784 (1976).
7. L. Jacobs, C. Rebbi, Interaction energy of superconducting vortices. Phys. Rev. B 19, 4486–4494 (1979).
8. E. J. Weinberg, Multivortex solutions of the Ginzburg-Landau equations. Phys. Rev. D 19, 3008–3012(1979).
9. C. H. Taubes, ArbitraryN-vortex solutions to the first order Ginzburg-Landau equations. Commun. Math.
Phys. 72, 277–292 (1980).
10. A. V. Efanov, Exact N-vortex solutions to the Ginzburg-Landau equations for κ=1/√2. Phys. Rev. B 56,
7839–7842 (1997).
11. I. Luk’yanchuk, Theory of superconductors with κ close to 1/√2. Phys. Rev. B 63, 174504 (2001).
12. V. Moshchalkov, M. Menghini, T. Nishio, Q. H. Chen, A. V. Silhanek, V. H. Dao, L. F. Chibotaru, N.
D. Zhigadlo, J. Karpinski, Type-1.5 superconductivity. Phys. Rev. Lett. 102, 117001 (2009).

13. E. Babaev, J. Carlström, J. Garaud, M. Silaev, J. M. Speight, Type-1.5 superconductivity in multiband
systems: Magnetic response, broken symmetries and microscopic theory – A brief overview. Phys. C
Supercond. 479, 2–14 (2012).
14. R. Prozorov, R. W. Giannetta, A. A. Polyanskii, G. K. Perkins, Topological hysteresis in theintermediate state of type-I superconductors. Phys. Rev. B 72, 212508 (2005).
15. G. D. Cody, R. E. Miller, Magnetic transitions of superconducting thin films and foils. I. Lead. Phys.
Rev. 173, 481–493 (1968).
16. G. J. Dolan, J. Silcox, Critical thicknesses in superconducting thin films. Phys. Rev. Lett. 30, 603–606(1973).
17. T. Cren, D. Fokin, F. Debontridder, V. Dubost, D. Roditchev, Ultimate vortex confinement studied byscanning tunneling spectroscopy. Phys. Rev. Lett. 102, 127005 (2009).
18. A. Kamlapure, M. Simonato, E. Sierda, M. Steinbrecher, U. Kamber, E. J. Knol, P. Krogstrup, M. I.
Katsnelson, M. Rösner, A. A. Khajetoorians, Tuning lower dimensional superconductivity withhybridization at a superconducting-semiconducting interface. Nat. Commun. 13, 4452 (2022).
19. M. A. Silaev, V. A. Silaeva, Self-consistent electronic structure of multiquantum vortices insuperconductors at T ≪ Tc. J. Phys. Condens. Matter 25, 225702 (2013).
20. T. Yamane, Y. Nagai, K. Tanaka, N. Hayashi, Impurity scattering effect on the zero-energy peak of thelocal density of states in a multi-quantum vortex core. Phys. C Supercond. 494, 128–130 (2013).
21. Y. Masaki, Y. Kato, Impurity effects on Caroli–de Gennes–Matricon mode in vortex core insuperconductors. J. Physical Soc. Japan 85, 014705 (2016).
22. T. Cren, L. Serrier-Garcia, F. Debontridder, D. Roditchev, Vortex fusion and giant vortex states in
confined superconducting condensates. Phys. Rev. Lett. 107, 097202 (2011).
23. C. K. Poole, H. A. Farach, R. J. Creswick, Handbook of Superconductivity (Academic Press, 2000).

24. T. Balashov, M. Meyer, W. Wulfhekel, A compact ultrahigh vacuum scanning tunneling microscopewith dilution refrigeration. Rev. Sci. Instrum. 89, 113707 (2018).
25. M. Ruby, B. W. Heinrich, J. I. Pascual, K. J. Franke, Experimental demonstration of a two-band
superconducting state for lead using scanning tunneling spectroscopy. Phys. Rev. Lett. 114, 157001(2015).
26. T. G. Saunderson, J. F. Annett, B. Újfalussy, G. Csire, M. Gradhand, Gap anisotropy in multiband
superconductors based on multiple scattering theory. Phys. Rev. B 101, 064510 (2020).
27. G. Chanin, J. P. Torre, Critical-field curve of superconducting lead. Phys. Rev. B 5, 4357–4364 (1972).
28. R. Prozorov, Equilibrium topology of the intermediate state in type-I superconductors of different
shapes. Phys. Rev. Lett. 98, 257001 (2007).
29. R. Prozorov, A. F. Fidler, J. R. Hoberg, P. C. Canfield, Suprafroth in type-I superconductors. Nat. Phys.
4, 327–332 (2008).
30. R. Prozorov, J. R. Hoberg, Dynamic formation of metastable intermediate state patterns in type-I
superconductors. J. Phys. Conf. Ser. 150, 052217 (2009).
31. C. Caroli, P. G. De Gennes, J. Matricon, Bound Fermion states on a vortex line in a type IIsuperconductor. Phys. Lett. 9, 307–309 (1964).
32. F. Gygi, M. Schlüter, Self-consistent electronic structure of a vortex line in a type-II superconductor.
Phys. Rev. B 43, 7609–7621 (1991).
33. G. E. Volovik, Vortex core anomaly from the gapless fermions in the core. JETP Lett. 58, 455 (1993).
34. G. E. Volovik, Topological superfluids. J. Exp. Theor. Phys. 129, 618–641 (2019).
35. S. M. M. Virtanen, M. M. Salomaa, Multiquantum vortices in superconductors: Electronic andscanning tunneling microscopy spectra. Phys. Rev. B 60, 14581–14584 (1999).

36. S. M. M. Virtanen, M. M. Salomaa, Microscopic structure of multiquantum vortices insuperconductors. Phys. B Condens. Matter 284-288, 741–742 (2000).
37. H. F. Hess, R. B. Robinson, R. C. Dynes, J. M. Valles, J. V. Waszczak, Scanning-tunneling-microscopeobservation of the abrikosov flux lattice and the density of states near and inside a fluxoid. Phys. Rev.
Lett. 62, 214–216 (1989).
38. S. Hasegawa, T. Matsuda, J. Endo, N. Osakabe, M. Igarashi, T. Kobayashi, M. Naito, A. Tonomura, R.
Aoki, Magnetic-flux quanta in superconducting thin films observed by electron holography and digitalphase analysis. Phys. Rev. B 43, 7631–7650 (1991).
39. J. Ge, J. Gutierrez, J. Cuppens, V. V. Moshchalkov, Observation of single flux quantum vortices in theintermediate state of a type-I superconducting film. Phys. Rev. B 88, 174503 (2013).
40. Y. Tanaka, A. Hasegawa, H. Takayanagi, Energy spectrum of the quasiparticle in a quantum dot
formed by a superconducting pair potential under a magnetic field. Solid State Commun. 85, 321–326(1993).
41. K. Tanaka, I. Robel, B. Jankó, Electronic structure of multiquantum giant vortex states in mesoscopic
superconducting disks. Proc. Natl. Acad. Sci. U.S.A. 99, 5233–5236 (2002).
42. L. D. Landau, The intermediate state of supraconductors. Nature 141, 688 (1938).
43. L. D. Landau, On the theory of the intermediate state of superconductors. J. Phys. USSR 7, 99 (1943).
44. S. Conti, F. Otto, S. Serfaty, Branched microstructures in the Ginzburg-Landau model of type-I
superconductors. arXiv:1507.00836 (2015).
45. H. F. Hess, R. B. Robinson, J. V. Waszczak, Vortex-core structure observed with a scanning tunnelingmicroscope. Phys. Rev. Lett. 64, 2711–2714 (1990).
46. M. V. Berry, M. Wilkinson, Diabolical points in the spectra of triangles. Proc. R. Soc. Lond. A Math.
Phys. Sci. 392, 15–43 (1997).

47. A. Odobesko, F. Friedrich, S.-B. Zhang, S. Haldar, S. Heinze, B. Trauzettel, M. Bode, Anisotropic
vortices on superconducting Nb(110). Phys. Rev. B 102, 174502 (2020).
48. H. Kim, Y. Nagai, L. Rózsa, D. Schreyer, R. Wiesendanger, Anisotropic non-split zero-energy vortex
bound states in a conventional superconductor. Appl. Phys. Rev. 8, 031417 (2021).
49. L. Kramer, W. Pesch, Core structure and low-energy spectrum of isolated vortex lines in clean
superconductors at T ≪Tc. Z. Phys. 269, 59–64 (1974).
50. M. Ichioka, V. G. Kogan, J. Schmalian, Locking of length scales in two-band superconductors. Phys.
Rev. B 95, 064512 (2017).
51. M. Silaev, E. Babaev, Microscopic theory of type-1.5 superconductivity in multiband systems. Phys.
Rev. B 84, 094515 (2011).
52. V. G. Kogan, J. Schmalian, Ginzburg-Landau theory of two-band superconductors: Absence of type1.5 superconductivity. Phys. Rev. B 83, 054515 (2011).
53. A. J. Leggett, Number-phase fluctuations in two-band superconductors. Prog. Theor. Phys. 36, 901–
930 (1966).
54. V. H. Dao, L. F. Chibotaru, T. Nishio, V. V. Moshchalkov, Giant vortices, rings of vortices, andreentrant behavior in type-1.5 superconductors. Phys. Rev. B 83, 020503 (2011).
55. A. V. Samokhvalov, V. D. Plastovets, A. S. Mel’nikov, Topological transitions in electronic spectra:
Crossover between Abrikosov and Josephson vortices. Phys. Rev. B 102, 174501 (2020).
56. S. Rex, R. Willa, A topological flux trap: Majorana bound states at screw dislocations. New J. Phys. 24,053057 (2022).
57. N. Schopohl, Transformation of the Eilenberger equations of superconductivity to a scalar Riccati
equation. arXiv:cond-mat/9804064 (1998).
58. C. Elsässer, N. Takeuchi, K. M. Ho, C. T. Chan, P. Braun, M. Fähnle, Relativistic effects on groundstate properties of 4D and 5D transition metals. J. Phys. Condens. Matter 2, 4371–4394 (1990).

59. B. Meyer, M. Fähnle, Ab initio calculation of the formation energy and the formation volume ofmonovacancies in Mo. Phys. Rev. B 56, 13595–13598 (1997).
60. D. Vanderbilt, Optimally smooth norm-conserving pseudopotentials. Phys. Rev. B 32, 8412–8415(1985).
61. L. Kleinman, Relativistic norm-conserving pseudopotential. Phys. Rev. B 21, 2630–2631 (1980).
62. R. Heid, K.-P. Bohnen, I. Y. Sklyadneva, E. V. Chulkov, Effect of spin-orbit coupling on the electronphonon interaction of the superconductors Pb and Tl. Phys. Rev. B 81, 174527 (2010).
63. J. P. Perdew, Y. Wang, Accurate and simple analytic representation of the electron-gas correlationenergy. Phys. Rev. B 45, 13244–13249 (1992).
64. T.-S. Choy, J. Naset, J. Chen, S. Hershfield, C. Stanton, A database of Fermi surface in virtual reality
modeling language (VRML). Bull. Am. Phys. Soc. 45, L36–L42 (2000).
65. G. Eilenberger, Transformation of Gorkov’s equation for type II superconductors into transport-likeequations. Z. Phys. 214, 195–213 (1968).
66. A. I. Larkin, Y. N. Ovchinnikov, Quasiclassical method in the theory of superconductivity. Sov. Phys.
JETP 28, 1200 (1969).
67. G. Carneiro, E. H. Brandt, Vortex lines in films: Fields and interactions. Phys. Rev. B 61, 6370–6376(2000).
68. N. Schopohl, K. Maki, Quasiparticle spectrum around a vortex line in a d-wave superconductor. Phys.
Rev. B 52, 490–493 (1995).