Robust and tunable coreless vortices and fractional vortices in chiral d -wave superconductors

,


I. INTRODUCTION
Two of the most outstanding issues in condensed matter physics are the direct identification of the superconducting pairing symmetry in unconventional superconductors and of the topological invariant in topologically non-trivial materials.These difficulties severely limit the ability to correctly interpret experiments and the applicability of newly discovered superconducting and topological materials.In few other systems is this as problematic as in multi-component superconductors, especially chiral superconductors, where both topology and superconducting symmetry need to be identified.Theoretically, chiral superconductors, and more generally chiral superfluids, are characterized by a non-trivial topology [1][2][3][4][5] and a discretely degenerate ground state that spontaneously breaks time-reversal symmetry [6].They belonging to the class of integer quantum Hall systems [7][8][9] with a finite Chern number generated by the winding of the superfluid order parameter [10][11][12][13][14][15][16][17], and with topologically protected chiral edge modes generating spontaneous surface currents and orbital angular momentum (OAM) [18][19][20][21][22][23].
In this complementary work, we demonstrate a strong tunability of the CV size and shape, also directly reflected in e.g. the LDOS.Furthermore, we provide extensive data that demonstrate a strong robustness of the results for a range of realistic models, over extensive parameter ranges, and in the presence of additional vortices or disorder.Overall, we relate the robustness of the experimental signatures of chirality, pairing symmetry, and Chern number, to the fact that they fundamentally stem from the parallel versus antiparallel alignment of vorticity and chirality, which are both topologically protected.In contrast, a non-chiral superconductor lacks this alignment possibility, since it lacks chirality.Our work therefore establishes CVs as a robust signature of spin-singlet chiral d-wave superconductivity, and furthermore the realization of fractional vortices in these materials.
This work is organized as follows.In Sec.II we summarize our model and methods, and describe basic properties of chiral d-wave superconductors.In Sec.III we introduce the basic properties of CVs, also discussing their overall stability and formation.In Sec.IV we demonstrate the large tunability of the CV size due to thermodynamics and electrodynamic interactions.Similarly, we study the interaction between CVs and other vortices in Sec.V and the behavior of CVs in confinement in Sec.VI, again demonstrating a tunability of both the CV size and shape as well as establishing strong robustness of CVs.We further demonstrate robustness against more general and anisotropic Fermi surfaces (FSs) in Sec.VII, non-degenerate order parameter components in Sec.VIII, and non-magnetic impurities in Sec.IX.Finally in Sec.X we briefly summarize our results.

A. Model
We consider weak-coupling superconductivity in equilibrium and in two dimensions (2D), assuming spindegeneracy, all appropriate for a quantitative description of a spin-singlet d-wave superconductor.We start by studying clean superconductors shaped like discs, with an electron-doped and circular Fermi surface (FS).We then relax all these assumptions by studying systems with either different discrete rotational symmetries or completely irregular shapes, as well as hole-doped and anisotropic FSs.We also consider non-degenerate order parameters, as well as dirty superconductors with non-magnetic impurities.For the specific setup, we align the superconducting plane with the xy-axes and use a perpendicular (orbital) external magnetic flux density B ext = (Φ ext /A) ẑ with homogeneous flux Φ ext across the system area A to induce vortex states.We assume type-II superconductivity appropriate for most non-elemental or unconventional superconductors, but consider different penetration depths λ 0 ∈ [2, ∞), via the Ginzburg-Landau coefficient κ ≡ λ 0 /ξ 0 .The penetration depth sets the length scale and strength of flux screening, defined by We perform our numerical simulations using the opensource framework SuperConga [174], which is a stateof-the-art implementation of the quasiclassical theory of superconductivity [161][162][163][164][165][166][167][168][169][170][171][172][173], running on graphics processing units (GPUs), and with extensive documentation and unit testing [175,176].SuperConga solves self-consistently [177] for both the superconducting or-der parameter ∆(p F , R) and vector potential A(R) via the gap equation and Maxwell's equations, respectively.Here, p F = p F (cos θ F , sin θ F ) is the Fermi momentum with angle θ F on the FS, while R = R(cos ϕ, sin ϕ) is the in-plane center-of-mass coordinate with polar angle ϕ.SuperConga also solves for impurity self-energies selfconsistently using the well-established t-matrix approach [178].
The low-energy expansion results in quasiclassical propagators, which we express in Nambu (particle-hole) space as with quasiparticle propagator g(p F , R; z) and anomalous pair propagator f (p F , R; z), where "tilde" denotes particle-hole conjugation α(p Here, z is the quasiparticle energy associated with the corresponding propagator, and is generally complex valued.Specifically, the retarded propagators g R (p F , R; ε) are used for spectral quantities, evaluated at z R ≡ ε + iδ with real energy ε and small positive broadening δ.For all other quantities, we use the Matsubara propagators g [185][186][187][188][189][190].The propagators in Eq. ( 1) are obtained via the Eilenberger equation [161] together with the normalization condition ĝ2 = −π 2 1, where ĥ is the self energy and τi the Pauli matrices in Nambu space.The self energy in Nambu space is with mean-field superconducting order parameter ∆(p F , R), while the diagonal part in the present work is capturing electrodynamic interactions via Σflux (R) (described further below) and impurity scattering via Σimp (p F , R; z) (described in Sec.IX).We parametrize the even-parity spin-singlet order parameter ∆(p F , R) via where Γ labels the irreducible representations of the crystallographic point group and the basis function η Γ (p F ) encodes the pairing symmetry on the FS [191], also related to the attractive pairing interaction V via Here, V Γ is the pairing strength of the respective symmetry channel.We self-consistently compute ∆(p F , R) via the superconducting gap equation with cutoff energy Ω c [172], and FS average [192,193] ⟨. . .
The electrodynamics is modelled via where and to the induced (ind) magnetic-flux density B ind (R) (i.e.screening) from the total charge-current density j(R) via Ampère's law We compute j(R) via We further compute the LDOS via Finally, we note that the quasiparticle energies are effectively Doppler shifted by the vector potential and any phase gradients [194][195][196][197], seen by applying a unitary gauge transformation to the Eilenberger equation as in e.g.Refs.[174,198,199], modifying Eq. ( 9), with the gauge-invariant superfluid momentum (superflow) This allows phase gradients and vector potentials to be treated on an equal footing, and leads to the Doppler shifted quasiparticle energy z p = z − v F (p F ) • p s (R) in the Eilenberger equation Eq. ( 2) [200][201][202], thus also influencing the LDOS in Eq. (12).

C. Chiral superconductivity
We consider spin-singlet chiral d-wave superconductivity, modelled using an attractive pair potential for the two irreducible d-wave representations Γ ∈ {d Following the notation in Eq. ( 5), the resulting order-parameter components ∆ d x 2 −y 2 (p F , R) and ∆ dxy (p F , R) are referred to as the nodal components.We initially assume that these channels are degenerate, since such a degeneracy is guaranteed by symmetry in any material with a three-or six-fold rotationally symmetric lattice [17], relevant for many of the recently proposed chiral d-wave superconductors [105][106][107][108][109][110][111][112][113][114][115][116][117][118].Still, for sake of full completeness, we later relax this assumption.Furthermore, we note that our theoretical framework includes other pair correlations allowed by symmetry, e.g.s-wave [174], while the possibility of additional attractive interactions in other pair channels is left as an outlook [203].
In order to better quantify chiral superconductivity, we transform the nodal order parameters to the eigenbasis of the OAM operator with the chiral order parameter components which are the two degenerate ground states in a bulk chiral superconductor.Below T c the system spontaneously chooses one of these as the dominant bulk chirality, e.g ∆(p F , R) = ∆ + (p F , R), while the opposite chirality ∆ − (p F , R) becomes subdominant and vanishes asymptotically in the translationally invariant bulk [204].Thus, the ground state of a chiral superconductor is described by a complex-valued order parameter that spontaneously breaks time-reversal symmetry [10,205,206], with a fully gapped bulk spectrum and Cooper pairs with an OAM l orb z = ±|M |ℏ [21].Even (odd) |M | correspond to spin-singlet (spin-triplet), and |M | = 1, 2 generate chiral p, d-wave order parameters, respectively.In this work we focus on spin-singlet chiral d-wave superconductivity with √ 2, which when equating Eq. ( 5) with Eq. ( 15) yields the relation between the two parametrizations In a chiral superconductor, the topological invariant is the Chern number M corresponding to the winding of the superconducting order parameter on the FS and giving rise to |M | chiral edge modes traversing the bulk gap whenever the topoloigcal invariant changes, in particular at vacuum interfaces but also domain walls [10][11][12][13][14][15][16][17].While these edge modes are topologically protected, they generate chiral edge currents and OAM which are not [7,126,127].Furthermore, close to the edges, the opposite (subdominant) chirality is often also locally induced, such that the order parameter takes the more general form in Eq. (15).This extends more generally to other forms of spatial inhomogeneities such as domain walls and vortices, and we therefore always use the most general form in Eq. (15) in our calculations, allowing for a completely general spatial dependence of both amplitudes and phases.We note that this in principle allows the system to go into a different state, e.g. a nodal dwave or nematic d-wave state [207], but we always find the chiral state to be robust.
Chiral superconductors also hosts domain walls, which are topological defects separating regions of opposite dominant chirality [1,28,29].Domain walls thus have |M | chiral edge modes on each side with opposite winding [30], also generating chiral currents on either side.These currents, together with the exchange of chirality across the domain wall, lead to a slight increase in free energy and an effective line tension [70].This usually makes domain walls metastable, but they are often trapped and further stabilized by pinning, geometric effects, and vortices [208].
Just like any superconductor, a chiral superconductor can also host vortex defects.A chiral superconductor with total vorticity m is associated with an m × 2π quantized phase winding in the dominant chiral component [56], i.e. χ + (R) ≈ mϕ along any path sufficiently far from and encircling all vortex defects.Abrikosov vortices (antivortices [209][210][211][212][213][214]) correspond to m = −1 (m = +1) in positive external flux Φ ext > 0, and vice versa for negative flux, also with a corresponding 2π phase winding in each nodal component χ d x 2 −y 2 (R) and χ dxy (R ′ ) if the vortex cores are overlapping, R = R ′ .Spatially separating the nodal winding centres R ̸ = R ′ leads to a disassociation of the Abrikosov vortex into two fractional vortices, one for each winding center, and to a Josephson-like term in the free energy that usually grows with the separation distance [75,215,216], thus making the fractional vortices unstable.However, inside a domain wall such a separation typically becomes favorable instead [70].Furthermore, the slight suppression of the total order parameter in the domain wall acts as an attractive pinning center for Abrikosov vortices, providing a mutual stabilization of the domain wall and fractional vortices [217], and thereby a mechanism for forming a CV as demonstrated in the next section III A.
Finally, changing magnetic flux direction allows for the vorticity to either be aligned antiparallel or parallel with the chirality, which leads to inequivalent vortices and also to inequivalent CVs.We illustrate this by first considering the total angular momentum, Lz = Lorb z + Lc.m.  = (M + m)ℏ, and is therefore a superposition between the OAM generated by chirality (i.e.Chern number) and the c.m. angular momentum generated by vorticity (i.e.winding quantization).Thus, antiparallel (parallel) alignment of vorticity and chirality leads to a negative (positive) superposition of the total angular momentum.Similarly, the phase winding of the subdominant chirality also shows such a behavior.Close to a vortex defect, the subdominant chirality is generally induced with finite amplitude and phase χ − (R) ≈ pϕ.The quantized phase winding p is constrained according to the relation [56,136] here with integer n capturing higher-order harmonics generated by e.g. a non-circular system or anistropic FS [191].Despite such terms often being unimportant [56], we in this work include them for full completeness.Equation (19) shows that the phase winding of the subdominant component also is a superposition of the vorticity and Chern number and can therefore be minimized (maximized) for an antiparallel (parallel) alignment.

III. CORELESS VORTICES
We begin this section by briefly summarizing the basic structure and properties of CVs in spin-singlet chiral d-wave superconductors in Sec.III A. In Sec.III B we discuss the stability and formation of CVs, and that the most stable CV is typically quadruply quantized in chiral d-wave superconductors.

A. Coreless vortex structure
In this subsection we summarize the basic properties of antiparallel and parallel CVs.In comparison to our earlier work [136], we here choose to study a somewhat smaller system with slightly different parameters, to illustrate that the important qualitative features do not depend on such parameters.Note that the spatial inhomogeneities induced by the CVs therefore are significant compared to the system size.Still, when we use the term 'dominant bulk chirality', we refer to the spontaneously chosen ground state chirality in the absence of vorticity, or equivalently, the dominant chirality in a much larger but otherwise analogous system.For reference, see Appendix A showing that the important qualitative features discussed here remain in systems with radius R ≥ 150ξ 0 , i.e. an order of magnitude larger.
Figure 1 (Fig. 2) shows an antiparallel (parallel) CV in a disc-shaped system with radius R = 15ξ 0 , dominant bulk chirality ∆ + , temperature T = 0.1T c , penetration depth λ 0 = 10ξ 0 , and external flux Φ ext = 7.5Φ 0 (Φ ext = −7.5Φ0 ).The first (second) row shows the amplitudes and phases of the nodal (chiral) components, while the third row shows the charge-current density and induced magnetic-flux density.Each nodal component has four 2π-phase windings that suppress the corresponding nodal amplitude but not the other nodal component.Since they lie at different coordinates for the two nodal components, the total order parameter is everywhere non-singular with no normal-state regions or core.The vortices are therefore fractional, in contrast to singular Abrikosov vortices which have spatially overlapping phase windings.A total of eight fractional vortices lie on a circularly (octagonally) formed domain wall, the latter seen in the second row of Fig. 1 (Fig. 2), where it separates the outer and inner regions of dominant chiralities ∆ + and ∆ − , respectively.There is a total vorticity of m = ±4 in the disc, seen by the m × 2π winding of the dominant chirality χ + (R) = mϕ in the outer region, which means this is a quadruply quantized CV.There is no phase-winding in the inner region, since there the dominant phase is constant χ − = 0, indicating that the vorticity is distributed along the domain wall.We next turn to the subdominant phase, which shows a π-shift across the domain wall in Fig. 1, which further stabilizes the structure but is otherwise unimportant [136].Thus, apart from this phase shift, the phase χ − is completely trivial in Fig. 1.In contrast, Fig. 2 shows a total of eight winding centers in the subdominant component χ − in the outside region.These results are in full agreement with the phase winding constraint in Eq. ( 19), with p = −4 + 4 = 0 (p = 4 + 4 = 8) in Fig. 1 (Fig. 2) corresponding to a CV with antiparallel (parallel) alignment of vorticity and chirality.Importantly, Fig. 2 shows that the winding centers lie outside the CV, spontaneously breaking axial symmetry, as defined by the winding not being generated by rotation around a single central axis.This occurs in order to lower the free energy, since the hypo-   thetical axisymmetric state with p = 8 would correspond to a giant vortex with a large normal core [56], consequently suppressing superconductivity and increasing the free energy.We have verified that such a giant vortex is indeed unstable.In contrast, the axially symmetrybreaking CV is stable since it avoids the energy penalty, while still lowering the kinetic energy caused by the external flux.Thus, the antiparallel (parallel) CV is axisymmetric (non-axisymmetric) with a continuous (discrete eight-fold) rotational symmetry.Our earlier work showed that this leads to a smoking-gun signature in the LDOS of both the topologically protected and quantized Chern number and vorticity [136].The third row shows a corresponding rotation symmetry of the charge-current density and induced magnetic flux, with multiple sign changes due to the chiral edge modes, domain wall, and overlapping Meissner screening currents.The paramagnetism is maximal along the domain wall, leading to a characteristic ring-like magnetic structure, in contrast to a point-like structure of an Abrikosov vortex [136].
Overall these results demonstrate the structure and basic properties of quadruply quantized CVs in chiral dwave superconductors, i.e. quadruple-quantum vortices.This is the chiral d-wave extension of the double-quantum vortex in chiral p-wave superfluids [56][57][58][59][60][61][62][63][64][65][66].Beyond this comparison, we note that an extension between the two different systems can in general be very non-trivial due to the different spin-symmetries and angular momentum quantization, and therefore it is not a priori certain that the same kind of vortex defects are even stable in both systems, let alone have the same qualitative properties.For example, the parallel CV in Fig. 2 shows multiple sign changes, compared to no sign changes for the parallel CV in a chiral p-wave double-quantum vortex reported in Ref. [56].More generally, "p-wave is special" [128] in many regards compared to all systems with higher Chern number, e.g. when it comes to the chiral edge currents and OAM [127][128][129][130][131][132][133][134][135].

B. Stability and coreless vortex formation
We next discuss the stability and formation of CVs.The peculiar combination of a domain wall and vorticity in a CV allows the system to carry finite vorticity, which reduces the kinetic energy caused by external flux but without paying the price of a normal core.This is significant, since the superconducting state is per definition the most energetically favorable state below the second critical field B c,2 (T ).The CV will thus be energetically more favorable than Abrikosov vortices if this gain outweighs the cost of the domain wall.However, the CV is very robust even when there are other vortex configurations with a lower free energy, i.e. even when technically metastable.This is a general feature of both Abrikosov vortices and CVs, related to the fact that they are topological defects that cannot be trivially removed from the system.They typically have to enter and exit the system via the edges, but such entrance and expulsion is hampered by large energy barriers, e.g.geometric and Bean-Livingston barriers [174,[218][219][220][221][222].Moreover, vortex motion is hampered by pinning and dissipation associated with normal-state resistance.Thus, once a particular arrangement of vortex defects has entered the system, it can become extremely robust even far into the flux-temperature regime where other vortex arrangements technically have even a significantly lower energy.Summarized briefly, experiments to a large degree observe metastable states [223], and such behavior is also typical in self-consistency simulations.Thus, the most relevant question is not necessarily whether a particular vortex configuration has the lowest energy, but if the necessary conditions for its formation can be prepared [220].
We see the vortex stability repeatedly in our simulations, both for CVs and Abrikosov vortices.In particular, we find that CVs spontaneously enter the system instead of Abrikosov vortices in certain parameter regimes, or can easily form when both domain walls and vortices are present.In the latter scenario, we find that the domain wall attracts and pins the Abrikosov vortices and, upon entering the domain wall, they disassociate into fractional vortices that lowers the free energy [57,70,75,215,216].Conversely, to break the CV, the vortices have to exit the domain wall or the domain wall has to disappear.However, such vortex expulsion is prevented by the pinning, and more importantly, the instability of the fractional vortices outside the domain wall.Thus, fractional vortices typically first have to recombine to a regular Abrikosov vortex before expulsion, but such re-combination increases the free energy.For the domain wall to disappear, it either has to shrink to zero size or expand to the system edges.Such shrinking is however prevented by the strong repulsive interaction between vortices, while expansion of the CV is counter-balanced by the attractive interaction (line tension) caused by the domain-wall currents, as well as the repulsive interaction between the fractional vortices and system boundaries.Among the very rare instances where we find the CV becoming unstable, it is this latter scenario that seems the most plausible; the line tension is significantly modified by non-degenerate order-parameter components (e.g.competing nodal superconductivity discussed in Sec.VIII) or the CV expanding to the system edge combined with a flux-temperature combination very far from the energy minimum (e.g. in Sec.IV).Apart from these scenarios, we find the CV to be extremely robust in all our calculations and often spontaneously appearing, even in the presence of strong perturbations, disorder, and when there are other vortex configurations with considerably lower free energy.
Finally, we discuss the most stable CV, which we gener-ally find to be the quadruply quantized CV with |m| = 4, shown in Figs. 1 and 2, and discussed in our previous work [136].This is easy to understand for the antiparallel CV, since it corresponds to the special commensurate scenario, such that the phase winding of the subdominant chirality vanishes p = m+2M = 0 [Eq.( 19)].In contrast, a finite phase winding p > 0 would either suppresses superconductivity if axisymmetric (thus costing energy), or increase the phase winding generating a modified superfluid momentum and line tension if non-axisymmetric (also costing energy).Furthermore, beyond commensurability, there is also the matter of balancing the repulsive versus attractive interactions, which overall stabilize the CV and its finite size (as discussed in Sec.IV), which is then important for the parallel CV, since there cancellation in p is impossible by definition.Considering for example higher vorticity |m| > 2|M |, this leads to increased repulsion, but also modified line tension due to the additional phase windings in p.As a consequence, the CV becomes less energetically favorable and less robust.The latter is also true for lower vorticity |m| < 2|M |, as there here might no longer be enough vortices to stabilize the domain wall.We verify these arguments during the extensive self-consistency calculations of the present work, including the large parameter ranges and model comparisons.Although we have found some parameter regimes where CVs with higher or lower vorticity become metastable rather than completely unstable, these states were generally less favorable and significantly more difficult to get to appear in the system.In summary, the commensurate scenario m = −2M allows the antiparallel CV to be coreless with maximized order parameter, leading to quadruple-quantum vortices in chiral d-wave superconductors, and more generally 2|M |-quantum vortices for other chiral superfluids.

IV. TUNABLE CORELESS VORTEX SIZE
This section demonstrates the large tunability of the CV size, via easily accessible parameters in experiment such as external flux Φ ext and temperature T , but also via the penetration depth λ 0 and system size R.For all the parameter ranges considered in this section, we note that the overall qualitative features presented in Figs. 1  and 2 remain the same.
The CV has a finite radius, R CV , balanced by attractive and repulsive interactions, acting to contract and expand the CV, respectively [220].The attractive interaction is exerted by the effective line tension from the domain wall and its chiral currents [70], while there is a mutual repulsive interaction between the fractional vortices in the domain wall [75,215].Hence, a closed domain wall will typically collapse and disappear in the absence of vorticity (we have verified this in our self-consistent calculations) [70].Furthermore, anything influencing the currents or vortices will change the balance, and therefore also R CV .This is also further demonstrated by studying  the interaction between CVs and Abrikosov vortices in Sec.V or with the system edges in Sec.VI.We start by describing how to unambiguously define and calculate R CV for antiparallel and parallel CVs.The midpoint of the CV is always well-defined, and a straight line across this point will generally intersect the domain wall of the CV twice, i.e. in two different points with degeneracy |∆ + | = |∆ − | as indicated in Fig. 3(a).We note that these points generally coincide with the maximum of the zero-energy LDOS [136].For the antiparallel CV, the CV diameter is the distance between the intersection points and is independent of the angle, and R CV is therefore unambiguously defined as half this distance as in Fig. 3(a).For the parallel CV, we instead define R CV from the average half distance for all angles, as displayed in Figs. ).However, we find that R CV is practically unambiguously defined even for the parallel CV, since ∆R CV ≡ R max CV − R min CV ≲ 1ξ 0 in all our simulations across all parameter ranges.
In Fig. 3(c) we illustrate that R CV can be effectively tuned by an externally applied magnetic field.Specifically, R CV decreases as |Φ ext | increases, since the currents grow in magnitude, while the distance between fractional vortices reduce (hence an overall stronger contraction).This is in a sense analogous with how larger flux causes smaller vortex separation and denser vortex lattices in regular type-II superconductors [224].Similarly, a shorter penetration depth λ 0 also leads to a smaller vortex-vortex separation, implying a smaller effective repulsive interaction as seen in Fig. 3(d).The overall dependence on λ 0 can be divided into two regimes: λ 0 < R, where screening becomes considerable and strongly modifies R CV , and λ 0 > R, where the system is poorly screened and the effect is minimal.In the limit of small λ 0 such that ξ 0 ≲ λ 0 ≪ R, the CV radius is almost completely determined by the screening regardless of system size, while in the opposite limit of large λ 0 , R CV eventually reaches the asymptotic limit λ 0 → ∞ (zero screening).We here note that the penetration depth is a materials property, which can be modified by the inclusion of impurities, as non-magnetic and magnetic impurities typically increase and decrease the penetration depth, respectively [225].
Figure 3(e) shows that R CV decreases at lower temperatures.We interpret the overall temperature dependence to be directly proportional to the effective coherence length ξ eff ≡ ℏv F /|∆(T )|, which reduces but saturates at small temperatures (due to saturating |∆(T )|) and increases dramatically at large temperature (due to vanishing |∆(T )|).Note that this is consistent with stronger but saturating chiral currents at lower temperatures, hence increasing the contraction.Of course, R CV is strictly limited by the size of the system (relative to the coherence length), consistent with the observed small (large) temperature dependence in small (large) systems in Fig. 3(e).Specifically, in the small systems there is a strong overlap between the boundary and CV at all temperatures leading to saturation, while the much weaker overlap in larger systems leave considerably more room for variation in R CV with T .We note an overall trend that R CV → 0.4R for large temperature, for all system sizes R considered in our simulations.In other words, the system size, and more generally surrounding environment, can strongly influence the maximum R CV and its temperature dependence.
Figure 3(f) shows directly how R CV increases with system size R.This is a mesoscopic finite-size effect, which can be divided into two regimes, corresponding to small and large R.For small R, the CV-induced currents strongly overlap with the chiral edge currents of the system.More importantly, the system edges impose an energy barrier [174,[218][219][220][221][222] and an effective repulsive interaction (at least at sufficiently high flux), which contracts the CV.This effect is also seen in Sec.VI, and is well-known for vortex lattices, leading to a number of interesting mesoscopic finite-size and shape effects [148,154,210,211,[226][227][228][229][230][231][232][233][234], see also Ref. [174] and references therein.For large R, the repulsion from the edges eventually becomes negligible, but there is still a slow asymptotic behavior of R CV which we interpret to be due to a slow saturation also present in properties related to the spectrum and chiral currents surrounding the CV.
In summary, these results demonstrate a strong tunability of the CV size, traced back to the effective attractive and repulsive interactions balancing the finite size [220], but also to the effective coherence length and its dependence on the superconducting gap.We note that while there are significant differences in the LDOS for the antiparallel and parallel CVs, due to symmetry breaking for the latter, the overall CV size is roughly simlar for both CVs.In the next Sec.V and Sec.VI we further demonstrate a tunability of the CV shape in the presence of other vortices or anisotropic effects.

V. INTERACTION WITH ABRIKOSOV VORTICES
In this section, we address how CVs coexisting with Abrikosov vortices changes the CV shape.In Figs 4 and  5 we show results for an antiparallel CV, and in Fig. 6 for a parallel CV.Before describing these figures in detail, we note that unlike the mostly point-like Abrikosov vortex, the CV has an intrinsic structure whose shape is to a large degree set by the repulsive interaction between its fractional vortices and their interaction with the environment, as established in the previous section IV.For example, the CV interacts repulsively with other vortices, whether it is Abrikosov vortices or other CVs, or attractively with antivortices.This section also establishes robustness of both the CV and its distinctive LDOS signature in the presence of such strong perturbations as additional vortices, and furthermore also shows the distinctly different LDOS signatures of CVs versus Abrikosov vortices.We note that the section also essentially studies the interaction between fractional vortices and regular Abrikosov vortices.Apart from illustrating all these aspects, combinations of CVs and Abrikosov vortices are reasonable to expect in a chiral d-wave superconductor, as discussed in the end of the section.
In Fig. 4 we present in the different columns the order parameter amplitudes and phases, while each row represents a different configuration of one antiparallel CV with one or more Abrikosov vortices or an antivortex.We start by analyzing the amplitudes and the overall CV shape, and then analyze the phases.The figure explicitly shows how the Abrikosov vortices completely suppress both the nodal components |∆ d x 2 −y 2 | and |∆ dxy | at its core, as well as both the chiral components |∆ + | and |∆ − |.By contrast, the fractional vortices in the domain wall of the CV only suppress the corresponding nodal component.Importantly, the figure shows signifi-cant modification of the overall CV size and shape.To explain these results, we note that at a sufficiently high external flux, Abrikosov vortices and CVs are both repelled from the system edges, related to the geometric barriers and the Bean-Livingston barrier [174,[218][219][220][221][222].This confinement leads to relatively small distances between Abrikosov vortices and the CV, which deforms the CV due to the mutual repulsive interaction.The resulting shape depends on the exact number and spatial arrangement of the Abrikosov vortices.Hence, we find that different deformation modes appear, as clearly seen in rows three to six.However, we note that the CV still keeps an overall elliptical form, which is clearly traced back to its original unperturbed circular form.If the Abrikosov vortex is instead situated at the center of the CV (first row), it is trapped and the CV expands due to the mutual repulsive interaction.If instead an antivortex is situated inside the CV (second row), it attracts the CV, which then shrinks substantially.However, we find that this configuration is always unstable unless pinning centers are artificially added (thus stabilizing the configuration), since the slightest deviation will otherwise fully attract the antivortex into the CV domain wall where it will be annihilated against two of the fractional vortices.We note that all other results and scenarios considered here are very robust even without such pinning, and we only choose to plot the antivortex scenario as it clearly illustrates how competing attractive and repulsive interactions set the overall shape and size of the coreless vortex.Specifically, all other results show a fully converged self-consistent solution, stabilized and trapped in the system by large energy barriers, and corresponding to a minimum of free energy.
Next, we study the phases and note in particular that the dominant chiral phase (i.e.χ ± outside and inside the CV, respectively) always winds according to the vorticity m, both locally around each vortex defect, and globally around the perimeter of the disc.For example, consider positive external flux and a vortex defect located at (x, y) = (x 0 , y 0 ) with winding m, where m = ∓1 for Abrikosov vortices and antivortices, respectively, while m = −4 for the CV considered here.Close to (x 0 , y 0 ), the dominant phase is described by χ + (R) ≈ mϕ with polar angle ϕ.Far from all the vortices at the disk perimeter, the dominant chiral phase globally winds χ + (R) ≈ m tot ϕ, with total vorticity m tot = −(N V +4N CV )+N AV , where N V counts the number of vortices, N CV the number of CVs, and N AV the number of antivortices.Thus, from top to bottom row, m tot = −(1 + 4), m tot = −4 + 1, m tot = −(4 + 1), m tot = −(4 + 2), m tot = −(4 + 3), m tot = −(4 + 4).We also find that the winding constraint p = m + 2M from Eq. ( 19) for the subdominant chiral phase is always fulfilled.
In Fig. 5 we display the spatially-resolved LDOS for the exact same systems and solutions as in Fig. 4, where each column is taken for a different fixed subgap energy ε (i.e.bias voltage).Importantly, at low energies, each Abrikosov vortex appears as a point-like peak represent- ing the Caroli-de-Gennes-Matricon states [147,232,[235][236][237], which expands to a size of roughly ∼ 1ξ 0 at higher energies.By contrast, the CV appears like a ring-like peak that is an order of magnitude larger already at zero energy, R CV ∼ 10ξ 0 .The CV expands into two concentric rings at higher energies, corresponding to the combined superflow generated by vorticity and the edge modes on either side of the domain wall.The intensity of the subgap states in the CV and Abrikosov vortex are also separated by an order of magnitude, but the LDOS peak of the CV should still be observable as it can be significantly larger than the coherence peak and is tunable by both temperature and flux, as shown in our earlier work [136].Notably, as the CV is deformed by the Abrikosov vortices, we also see how the LDOS is correspondingly deformed in rows 2-6.Thus the LDOS is explicitly tracking the CV shape.
Figure 6 shows the LDOS for similar combinations of a CV with Abrikosov vortices, but now for a parallel CV (Φ ext < 0 such that m = +4 and p = 8 instead of m = −4 and p = 0), and without the antivortex scenario.For completeness, Appendix B contains a plot of the corresponding order parameter amplitudes and phases for these scenarios (i.e. the analogue of Fig. 4).The overall trend in Fig. 6 is similar to that of the antiparallel CV, but importantly, we note that the distinct LDOS signature of the axial symmetry breaking is clearly present, including the eightfold symmetry related to p = m + 2M = 8 for the CV (instead of p = 0 for the antiparallel CV).Hence, despite the strong local perturbation caused by the presence of Abrikosov vortices, the overall Doppler shift caused by finite superflow from the p = 8 winding centers remains clearly distinguishable, as is thus then the direct signature of the Chern number M .Interestingly, the last row of Fig. 6 was initialized with four vortices outside the CV (i.e. the same arrangement as the last row Fig. 5), but during the selfconsistency loop, one of the vortices was spontaneously absorbed into the center of the CV, thereby lowering the free energy.Notably, this is a self-consistent and robust solution, illustrating that it is realistic to study and expect the appearance of configurations with an Abrikosov vortex trapped inside the CV.Importantly, during the trapping of the Abrikosov vortex in the self-consistency loop, the vortex was at some point located at the domain wall of the CV where it could have been disassociated into fractional vortices, thus leading to a higher quantized CV with m = +5 and p = +9.However, the displayed solution with m = 4 and p = 8 was still preferred.Hence, this is another strong indication that the quadruply quantized CV is the most robust CV in chiral d-wave superconductors.
Finally, on more general grounds, we point out that studying a combination of CVs and Abrikosov vortices is relevant, since both are robust topological defects and can thus appear simultaneously in a sample.This is further supported by the high energy barriers associated with vortex dynamics, meaning that a particular vortex solution can be trapped in the system far into its metastable regime [174,[218][219][220][221][222], where another vortex solution technically has a lower energy but can still not enter the system.Generally, both Abrikosov vortices and domain walls can be "kicked" into the system by e.g.annealing and rapid quenches in temperature and flux, and they can be further stabilized and trapped by pinning centers and certain geometry [208,238].Indeed, we find that combinations of CVs and Abrikosov vortices spontaneously enter and stabilize in our self-consistency calculations for different flux-temperature combiniations.
In summary, these results show a robustness of the CV in the presence of Abrikosov vortices.At the same time, a tunability of the shape is demonstrated, although the CV shape can still be traced back to its original circular (octagonal) shape for the antiparallel (parallel) CV.Specifically, the LDOS at different energies appear as concentric and convex (concave) line segments, corresponding to the Doppler shifts caused by the axisymmetric (nonaxisymmetric) superflow, which in turn is generated by the internal (external and internal) phase windings for the antiparallel (parallel) CV.We also note that these results give rise to an even stronger experimental signature in the LDOS, as the point-like Abrikosov vortex is distinctly different from the line-like CV.Finally, we propose that similar deformations might be caused by other strong local electromagnetic perturbations, e.g. an appropriately prepared STM tip with strong magnetization.

VI. NON-CIRCULAR GEOMETRY AND STRONG CONFINEMENT
The previous two sections IV and V illustrated that the overall size and shape of the CV is balanced by effective attractive versus repulsive electrodynamical interactions, traced back to the domain wall currents and fractional vortices respectively.In this section, we further illustrate this via the interaction with the system edges, and show how confinement alone can induce asymmetric deformation modes in the CV.In addition, the results show that the LDOS signature remains robust and is not relying on the symmetry (or lack thereof) of the system itself.
Figure 7 shows an antiparallel (parallel) CV in odd (even) rows in systems with different shapes, where the columns show from left to right: magnitude of the chiral order parameter components, charge-current density, and LDOS at different fixed subgap energies.The first two rows show a sample shaped like a pentagon, importantly illustrating that the overall circular versus octagon rotation symmetries of the CVs remain, even when incommensurate with the rotation symmetry of the system.Furthermore, this is an example of a system with higher-order harmonics discussed in Sec.II C, where Eq. ( 19) is modified with an additional term, such that p = m + 2M + n, here with integer n = −5 due to the five-fold rotational symmetry of the superconducting grain.This leads to additional phase gradients and therefore superflow, which in turn generates additional current components.This effect is responsible for helping the current turn the sharp corners of the system, which is a well-known effect in chiral superfluids [14].Furthermore, the additional phase gradients and superflow also leads to a locally enhanced LDOS at the corners at finite energies, again via the Doppler shift discussed in relation to Eq. ( 13).As a result, the subdominant phase χ − of the antiparallel CV has integer winding p = −4+4+n = −5, while p = 4 + 4 + n = 3 for the parallel CV.Hence, the higher-order harmonics is superimposed with the antiparallel versus parallel vorticity and chirality, especially seen by the additional signatures with five-fold rotational symmetry in the LDOS at high energies in the last column Fig. 7(h).Importantly, the higher-order harmonics still does not modify the overall strong signature of vorticity and chirality in the LDOS.In other words, the strong LDOS distinction between parallel and antiparallel CVs remains robust.Next, the third and forth rows show a completely irregular system without any rotation symmetry.Again, the LDOS signature is robust, but the sharp wedges together with the overall asymmetry between x and y directions cause a slight deformation of the CVs.The last two rows show a rectangular system, with an even stronger asymmetry between x and y directions.Due to the effectively repulsive interaction between the system edges and the fractional vortices in the CV, the resulting CV shape is strongly deformed, with a clear x and y asymmetry.The effective repulsive interaction with the system edges is related to the energy barriers for vortex entrance and expulsion at sufficiently high external flux [174,[218][219][220][221][222].In summary, this section illustrates both a tunability of the CV shape due to mesoscopic confinement, and most importantly that despite these CV shape changes, the experimental signatures in the LDOS are robust at all subgap energies and do not rely on the overall rotation symmetry of the system.

VII. NON-CIRCULAR FERMI SURFACES
So far, we have assumed a circular and electron-doped FS as in our previous work [136].Here we show that our main results and conclusions do not depend on the shape of the FS or particular doping level.In particular, we consider FSs formed in a hole doped material and with weak to strong deviation from a circular shape, and also with anisotropy between k x -and k y -momentum directions, to further mimic possible broken symmetries in the normal state.In particular, we parametrize a noncircular FS via the momentum k = k x kx + k y ky through the normal-state dispersion ϵ k on a square lattice in terms of the lattice constant a 0 , nearest-neighbor hopping t > 0 (which we use as a natural unit for all tightbinding energies), next-nearest neighbor hopping t ′ , nextnext nearest neighbor hopping t ′′ , and with anisotropy α xy between k x and k y [239].We here consider four different tight-binding models taken from the literature [147,240,241], labeled as FS #1 to #4, defined in Table I and illustrated in Fig. 8. Here, the non-circular FS leads to a modified v F (p F ) entering the Eilenberger equation ( 2), thus modifying the propagators and all other quantities defined in Sec.II correspondingly.See Ref. [193] for further details on parametrizing such a microscopic tightbinding Fermi surface within quasiclassical theory of superconductivity.We further note that all of these FSs are hole-doped, corresponding to being centered around (k x , k y ) = (π/a 0 , π/a 0 ), and show either a four-fold (FSs #1, #2, #4) or two-fold (FS #3) discrete rotational symmetry.There is a weak to strong deviation from circular Table I.Parametrization of tight-binding Fermi surfaces (FSs) for the normal-state dispersion in Eq. ( 20) with nearest neighbor hopping t, next-nearest neighbor hopping t ′ , next-next nearest neighbor hopping t ′′ , chemical potential µ, and hopping anisotropy αxy.Resulting FSs are illustrated in Fig. 8.  [193].
shape in changing between FS #1 to FS #4.In contrast, an electron-doped FS is centered around (k x , k y ) = (0, 0).We show the antiparallel (parallel) CV computed with these FSs in odd (even) rows in Fig. 9 for an octagonal sample and in Fig. 10 for a square sample.In all figures, the columns show from left to right the chiral order parameter amplitudes, charge-current density, and LDOS at different subgap energies ε.Overall, we find that both types of CVs show traces of the underlying symmetry of the FS, which can be explained in terms of  I, FS #1 to FS #4, with antiparallel (parallel) CV in odd (even) rows.Columns, from left to right, show magnitude of the chiral order parameter components, charge-current density, and LDOS at different fixed subgap energies.Here, ∆+ is the dominant bulk chirality with T = 0.1Tc, λ0 = 80ξ0, and with Φext = ±8Φ0 for antiparallel and parallel CVs respectively.
higher-order harmonics superimposed on the CV as discussed in Sec.II C and studied for non-circular samples in Sec.VI, but here they are instead originating from the FS.For example, a four-fold rotational symmetry of the FS (or sample) leads to a corresponding four-fold rotational symmetry with nodes and kinks developing in the CV.Similarly, an anisotropic FS with two-fold rotational symmetry (FS #3) deforms the otherwise circular CV into an ellipse, due to suppression (enhancement) of v F along k y (k x ), as illustrated in Fig. 8(c).Inter-estingly, the elliptical deformation occurs along opposite directions for the antiparallel and parallel CVs, which we interpret to be due to opposite signs of v F (p F )•p s (R) for the two CVs, which can be traced back to opposite signs of A (i.e.opposite external field directions) entering p s in Eq. (13).Furthermore, we note here that FS #4 is very distorted compared to a circular FS, leading to also very strong distortions in the CV.Despite these symmetry-breaking terms in the FS causing distortion of the CV, we find that both CV solu- tions are always robust, and the important asymmetry remains clear in the LDOS.Specifically, the CV with antiparallel vorticity and chirality (odd rows) generates convex and concentric lines in the LDOS, from the axisymmetric angular momentum and superflow.In contrast, the CV with parallel vorticity and chirality (even rows) always generates characteristic concave LDOS patterns due to the multiple phase winding centers which are non-overlapping (i.e.axial symmetry-breaking).Moreover, the emergent discrete rotational symmetry and interweaving resonances at higher energies is a direct experimental signature of the quantized vorticity and Chern number, due to it tracing back to winding superposition p = m + 2M , which is robust due to quantization and topology [136].These results establish that the signatures of CVs are robust even for a highly anisotropic FS, reflecting broken symmetries of the normal state.

VIII. NON-DEGENERATE NODAL COMPONENTS
In all other sections and our previous work [136] we assumed degeneracy between the two nodal d-wave pairing symmetries, d x 2 −y 2 and d xy , such that their transition temperatures are the same.Such an exact degeneracy is experimentally relevant: it is enforced by the symmetry and group theory in any material with a three-or six-fold rotationally symmetric lattice [17].This includes triangular, hexagonal, and honeycomb materials and is as such guaranteed in many of the materials currently proposed as chiral d-wave superconductors [105][106][107][108][109][110][111][112][113][114][115][116][117][118].Still, for sake of full completeness, we in this section show that our main results and conclusions in addition hold for systems where this degeneracy is somehow broken.Specifically, we consider non-degenerate pairing interactions modeled by different coupling constants resulting in different transition temperatures, quantified by the ratio Hence, we set the d xy -component to be subdominant for all α < 1, resulting also in different bulk amplitudes |∆ dxy | < |∆ d x 2 −y 2 |.However, apart from inserting these different coupling strengths, we do not constrain the order parameter components in any way, and solve for both of them completely self-consistently.For example, performing self-consistent calculations without any vorticity, we still find that the chiral d-wave state is the ground state even for highly non-degenerate systems with α < 0.8, thus surviving a strong suppression of the d xycomponent.Notably, such a state is still fully gapped in the bulk, with a Chern number M = ±2.Hence, the possibility of antiparallel versus parallel superposition of vorticity and chirality in a CV is still possible.
In order to investigate the effects on the CV from nondegenerate d-wave nodal components, we begin by summarizing the scenario of full degeneracy (α = 1) studied so far.Here, both the axially symmetric CV with antiparallel vorticity and chiralty, and the axial symmetrybreaking CV with parallel vorticity and chirality, are extremely robust solutions over a large range of temper-atures and flux.Notably, for degenerate nodal d-wave components in a disc-shaped system and FS, the total superconducting order parameter has full rotation symmetry for the antiparallel CV, and thus physical properties such as currents and magnetic fields generally do not reflect the four-fold symmetry of the individual nodal components.However, as the degeneracy between the nodal components is broken, it is reasonable to expect that the nodal four-fold rotational symmetry will be imprinted also on the antiparallel CV.
In Fig. 11 we study antiparallel and parallel CVs from weak non-degeneracy α = 0.99 (top two rows) and continue to strong non-degeneracy α = 0.8 (lowest two rows).For this full range of asymmetry, we find that both CVs are still very robust, but over a slightly narrower range of flux.By decreasing α we find that the broken degeneracy and four-fold nodal symmetry become more apparent in the CV, as expected.For example, along the domain wall, the suppression of ∆ d x 2 −y 2 (∆ dxy ) now occurs over a smaller (larger) region, as compared to α = 1.As α is further reduced, the dominant ∆ d x 2 −y 2 covers nearly the whole domain wall, except at four isolated points.Consequently, these points become the only locations in the domain wall where ∆ dxy is finite.The fractional vortices are however still well-separated, and the CV structure is notably still intact.This spatial structure of the individual nodal components leads to signatures also in all other quantities, including the chiral order parameters components, and also the currents, induced flux (not shown here), and LDOS.Still, the results in Fig. 11 show that the overall conclusions and experimental signatures established in the rest of the work for the degenerate case remain robust and clear also with a strong asymmetry between the two nodal d-wave components.In particular, the LDOS for the antiparallel CV keeps its overall concentric and convex circular lines due to the order parameters and currents also exhibiting such a profile, with non-degeneracy only turning the circles more square-like.Meanwhile, the parallel CV still shows concave octagonal structures, with eight-fold interweaving resonances at higher energy due to non-trival additional phase winding in the subdominant chirality (but now overlapped with strong four-fold structure).Thus the antiparallel CV keeps the axisymmetry, while the parallel CV does not, just as established in the rest of the results.This robustness in the different LDOS patterns between the two CVs is expected: after all, the LDOS patterns stem directly from a positive versus negative superposition of the quantized and topologically protected Chern number (OAM from chirality) and vorticity (c.m. angular momentum), as introduced in Sec.II C. Thus, as long as there is a chiral state, the positive versus negative superposition generates completely different scenarios, but possibly superimposed with higher-order contributions, in this case due to the additional broken nodal degeneracy.
Finally, for completeness, let us address the extreme limit of non-degeneracy, although this is not expected for stable chiral d-wave superconductors and thus not of central importance here.Eventually, the local variation of the current leads to a modified line tension, modifying the stability.Indeed, as α → 0.6, the parameter-space region of stable CV shrinks rapidly.At some point, the CV also becomes unstable and multiple regular Abrikosov vortices are instead stabilized.We define the critical ratio where this occurs as α * (T, Φ ext , λ 0 , R), hence possibly depending on all parameters such as temperature, flux, penetration depth, and system size.The full parameter space is of course far beyond the scope of the present work, but we consider a subset of the parameter space for illustrative purposes.For example, at fixed temperature T = 0.1T c , external flux |Φ ext | = 8Φ 0 , and λ 0 = 80ξ 0 , we find that the antiparallel CV is unstable below α * ≈ 0.70 at R = 25ξ 0 , α * ≈ 0.60 at R = 50ξ 0 , and α * ≈ 0.55 at R = 75ξ 0 , while the parallel CV is unstable below α * ≈ 0.73 at R = 25ξ 0 , α * ≈ 0.64 at R = 50ξ 0 , and α * ≈ 0.61 at R = 75ξ 0 .The CV is thus less stable in smaller systems, especially for increased non-degeneracy.We interpret this to stem from that smaller systems exhibit significant overlap between opposite system edges where both of the nodal components are suppressed, as well as between the CV and the system edge.This suppression is naturally enhanced by the non-degeneracy.As a result, the chiral state competes with both the normal state and a nodal d-wave state, which effectively hampers the formation of the chiral state, and consequently therefore also the formation of domain walls and CVs.

IX. NON-MAGNETIC IMPURITIES
Our earlier work [136] demonstrated that the LDOS signatures of the CV are robust under the inclusion of a phenomenological energy broadening δ of the spectrum.Such an energy broadening can be caused by e.g.disorder, impurity scattering, fluctuations, or interfaces with finite transparency [242][243][244][245][246][247].Here we exemplify this by studying dirty systems with non-magnetic impurities, and show that the CV as well as the LDOS signatures, are robust.
We model the non-magnetic impurities using the wellestablished t-matrix approach within the quasiclassical theory of superconductivity [178,248].The diagonal impurity self-energy from Eq. ( 4) then takes the form with dilute impurity concentration n i and impurityscattering matrix t fulfilling the additional selfconsistency equation  with scattering potential û0 (p F , p ′ F ; z).Equation ( 23) results from a diagrammatic expansion describing multiple scattering of quasiparticles and pairs by an impurity, connecting different scattering channels with momenta p F and p ′ F on the FS (here integrated over the momentum p ′′ F ).Here we assume equilibrium, a non-crossing approximation, and an isotropic scattering potential, such that û0 with scattering energy Γ u = n i /(πN F ) and scattering phase shift δ 0 = arctan(πu 0 N F ).We solve Eq. ( 24) self-consistently, together with the gap equation and Maxwell's equation.We define the "pair-breaking energy" as Γ = Γ u sin 2 δ 0 , related to the normal-state mean-free-path l = ℏv F /(2Γ).We consider two extreme limits, namely the weak-scattering Born limit (δ 0 → 0 and Γ u → ∞, such that Γ is constant) and the strong-  scattering unitary limit (δ 0 → π/2 and u 0 → ∞, such that Γ = Γ u ).In these limits, the equilibrium solutions simplify to respectively.We vary the scattering energy over orders of magnitude, γ u ≡ Γ/(2πk B T c ) ∈ [0.002, 0.1].By comparison, the zero-temperature bulk gap is roughly |∆ 0 |/(2πk B T c ) ≈ 0.280.We still use a phenomenological broadening δ/(2πk B T c ) = 0.0005 to avoid divergent LDOS for small γ u , but this value is an order of magnitude smaller than used in the rest of this work, δ/(2πk B T c ) = 0.005.Figures 12 and 13 show the resulting LDOS in the presence of non-magnetic impurities for an antiparallel and parallel CV, respectively, with left and right columns showing the Born and unitary limits, respectively, with the figure panels (c-d) showing line-cuts across the CV at zero energy and panels (e-f) showing the LDOS at the domain wall.As expected, the LDOS peaks are broadened when increasing the scattering energy, eventually becoming almost completely broadened for γ u → 0.1 as indicated by red lines in (c-f).This result is expected because such strong γ u is comparable with the bulk gap.
The broadening is also naturally not as strong in the Born limit (left) as in the unitary limit (right).Consequently, panels (a,c) illustrate that both CVs are strongly distinguishable in the spatially resolved LDOS even for γ u = 0.1 in the Born limit, while panels (b,d) show that both CVs are distinguishable at γ u = 0.05 in the unitary limit.We note that the peak at the disc center in Fig. 12 is a resonance related to the perfect rotation symmetry [136].
Finally, we note that the antiparallel CV radius R CV slightly increases by 1ξ 0 (2ξ 0 ) as γ u changes from 0.002 to 0.1 for the Born (unitary) limit, as indicated in Figs.12(c,d).For the parallel CV, the increase in R CV is even smaller, about 0.5ξ 0 , see Figs. 13(c,d).We expect R CV to increase more significantly with γ u in systems where R ≥ λ 0 .This is because non-magnetic impurities generally increase the penetration depth λ 0 [225], which in turn increases R CV as shown in Sec.IV.Here, in contrast, λ 0 = 80ξ 0 is much larger than R, explaining the small variation in R CV .In summary, we find that the LDOS signatures of the Chern number, superconducting pairing symmetry, and chirality is robust against strong (moderate) scattering energy in the Born (unitary) limit, despite a corresponding broadening of the LDOS peaks.Furthermore, we find that the CV itself is very robust in all cases considered.This demonstrates the viability of CVs and its signatures to identify chiral superconductivity also in dirty systems.

X. CONCLUSIONS
In this work, we show a strong tunability of CVs in spin-singlet chiral d-wave superconductors, as well as a robustness of their experimental signature for a large range of material models, parameter regimes, perturbations, anisotropy, and disorder.
In terms of tunability, we find that the finite size of the CV is balanced by the attractive and repulsive interactions exerted by its domain-wall currents and fractional vortices, respectively.Thus, we show that the overall size is easily tuned directly by changing an externally applied magnetic flux and the temperature, but also depend on system size and penetration depth, the latter generally tunable by artificially adding impurities.We also find that the overall shape is tunable and deforms in an anisotropic environment, e.g.due to other vortices, an irregular system shape, or an anisotropic FS.
For the experimental signatures, our earlier work established that the LDOS host distinct signatures for the two inequivalent CVs, with antiparallel or parallel chirality and vorticity, and that this can be used to clearly identify chirality, Chern number, and even the superconducting pairing symmetry [136].More specifically, the antiparallel CV is axisymmetric with a continuous rotation symmetry, associated with LDOS peaks appearing as concentric and convex circular lines.The parallel CV spontaneously breaks axial symmetry, generat-ing additional winding centers outside the CV, deforming its shape into a concave structure with discrete rotation symmetry, directly related to the Chern winding number.At zero energy (bias voltage), the LDOS directly probes the domain wall structure of the CV and its overall rotation symmetry and thereby the Chern number.At higher energies, there are additional interweaving resonances between the additional winding centers, even more clearly exhibiting the rotation symmetry and Chern number.This forms strong experimental signatures, directly measurable with STS and STM.In this work we establish that all of these signatures are robust for a large range of possible perturbations, system, and model changes.In particular, we demonstrate how the results hold in systems with incommensurate or no rotation symmetry, at strong confinement, for both electron doped and hole doped FSs or anisotropic FSs, non-degenerate nodal d-wave components, as well as when non-magnetic impurities are present.We also find robustness for large ranges of temperatures, external flux strength, penetration depths, and system sizes, as well as when additional Abrikosov vortices are present.
In conclusion, our work establishes CVs as a tunable and robust experimental signature of spin-singlet chiral d-wave superconductivity, which furthermore provide a platform to study fractional vortices.and in systems with opposite dominant bulk chirality ∆ − .In particular, Fig. 14 shows various quantities for a parallel CV in a disc with radius R = 150ξ 0 , dominant bulk chirality ∆ − , external flux Φ ext = 8Φ 0 , temperature T = 0.1T c , and penetration depth λ 0 = 80ξ 0 , to be compared to Fig. 2.There is still four well-separated fractional vortices in each nodal component, an octagonalshaped domain wall in the chiral amplitudes, and integer phase windings m = −4 in the dominant phase χ − and p = −8 in the subdominant phase χ + (reversed signs due to reversed bulk chirality).Furthermore, the currents still show the same overall spatial profile and number of sign changes.Interestingly, the additional p = −8 phase windings in the subdominant phase χ + remains at a small distance outside the CV.We note that for an antiparallel CV in such a large system, the π-shift in the dominant chirality instead remains closer to the edge (not shown here).This appendix contains additional numeric results for the interaction between CVs and Abrikosov vortices.In particular, Fig. 15 shows the order parameter amplitudes and phase windings for the same system as in  Fig. 6, i.e. this is the analogue of Fig. 4 but for a parallel (symmetry-broken) CV caused by negative external flux.Importantly, we note that despite all the additional Abrikosov vortices generating phase windings that overlap with the p = 8 winding centers of the CV, the latter still generates the distinct shape of the parallel CV studied in the rest of the work.
with elementary charge e = −|e|, Fermi velocity v F on the FS, normal-state density of state N F on the FS (per spin), and speed of light c.Here, our natural length unit is ξ 0 ≡ ℏv F /(2πk B T c ), sometimes referred to as an effective superconducting coherence length over which superconductivity spatially varies, with Planck constant ℏ, Boltzmann constant k B , and superconducting transition temperature T c .We study superconducting systems with a diameter or side length D ∈ [20, 300]ξ 0 , for different temperatures T ∈ [0.01, 0.99]T c , and external fluxes Φ ext ∈ [−15, 15]Φ 0 with flux quantum Φ 0 ≡ hc/2|e|.We keep all parameters fixed during the self-consistency simulations.
winding quantization.Here, Lorb z is the OAM generated by chirality as explained earlier in this subsection, while Lc.m. z = (ℏ/i)∂ ϕ is the generator of c.m. angular momentum with eigenvalue l c.m. z = mℏ for a state with vorticity m.Thus, the total angular momentum of the Cooper pair is l z = l orb z + l c.m. z
3(b) where the thick red line shows the numerically extracted point |∆ + | = |∆ − | and arrows show the minimum (maximum) radius R min CV (R max CV

Figure 4 .
Figure 4. Systems containing both an antiparallel CV and Abrikosov vortices or antivortices, in a disc-shaped system with radius R = 25ξ0, dominant bulk chirality ∆+, at T = 0.1Tc and λ0 = 80ξ0.Columns show, from left to right, the chiral and nodal amplitudes, then the chiral and nodal phases.In the first (second) row, an Abrikosov (anti)vortex trapped inside the CV, remaining rows an increasing number of Abrikosov vortices (one to four) outside the CV.External flux is from top to bottom row: 8Φ0, 8Φ0, 12Φ0, 12Φ0, 14Φ0, and 14Φ0 in order to stabilize the various configurations.

Figure 5 .
Figure 5. Same as Fig. 4, but with each column showing the LDOS at a fixed subgap energy ε, with gap roughly 1.76kBTc.

Figure 6 .
Figure 6.Same as Fig. 5, but for a parallel CV, without any antivortex scenario, and where one of the four vortices in the last row has spontaneously been trapped at the CV center.The system has dominant bulk chirality ∆+, at T = 0.1Tc and λ0 = 80ξ0, and is exposed to negative external fluxes, from top to bottom: −8Φ0, −12Φ0, −12Φ0, −14Φ0, and −14Φ0.

Figure 7 .
Figure 7. CVs in various non-circular samples, with antiparallel (parallel) CVs in odd (even) rows.Columns, from left to right, show magnitude of the chiral order parameter components, charge-current density, and LDOS at different fixed subgap energies.Here, dominant bulk chirality is ∆+, with T = 0.1Tc, λ0 = 80ξ0, and with Φext = ±8Φ0 for antiparallel and parallel CVs respectively.

Figure 8 .
Figure 8. Normal-state band structures for the non-circular tight-binding hole-doped FSs defined in Table.I. Colors indicate band energy E, solid lines denote FS (E = 0), with arrows indicating the Fermi velocity vF(pF) used as input in the quasiclassical parametrization[193].

Figure 9 .
Figure 9. CVs in an octagon-shaped sample.Different rows show the different FSs defined in Table I, FS #1 to FS #4, with antiparallel (parallel) CV in odd (even) rows.Columns, from left to right, show magnitude of the chiral order parameter components, charge-current density, and LDOS at different fixed subgap energies.Here, ∆+ is the dominant bulk chirality with T = 0.1Tc, λ0 = 80ξ0, and with Φext = ±8Φ0 for antiparallel and parallel CVs respectively.

Figure 10 .
Figure 10.Same as Fig. 9 but for a square-shaped sample.

Figure 15 .
Figure 15.Same as in Fig.4, but for a parallel CV and without the antivortex scenario.