Enhanced chiral edge currents and orbital magnetic moment in chiral $d$-wave superconductors from mesoscopic finite-size effects

Chiral superconductors spontaneously break time-reversal symmetry and host topologically protected edge modes, supposedly generating chiral edge currents which are typically taken as a characteristic fingerprint of chiral superconductivity. However, recent studies have shown that the total edge current in two dimensions (2D) often vanishes for all chiral superconductors except for chiral $p$-wave, especially at low temperatures, thus severely impeding potential experimental verification and characterization of these superconductors. In this work, we use quasiclassical theory of superconductivity to study mesoscopic disc-schaped chiral $d$-wave superconductors. We find that mesoscopic finite-size effects cause a dramatic enhancement of the total charge current and orbital magnetic moment (OMM), even at low temperatures. We study how these quantities scale with temperature, spontaneous Meissner screening, and system radius $\mathcal{R} \in [5,200]\xi_0$ with superconducting coherence length $\xi_0$. We find a general $1/\mathcal{R}$ scaling in the total charge current and OMM for sufficiently large systems, but this breaks down in small systems, instead producing a local maximum at $\mathcal{R} \approx 10{\text{--}}20\xi_0$ due to mesoscopic finite-size effects. These effects also cause a spontaneous charge-current reversal opposite to the chirality below $\mathcal{R}<10\xi_0$. Our work highlights mesoscopic systems as a route to experimentally verify chiral $d$-wave superconductivity, measurable with magnetometry.

The topologically non-trivial bulk order parameter of chiral superconductors leads to topologically protected chiral edge modes through the bulk-boundary correspondence [43][44][45][46][47][48][49].The chiral edge modes in turn should presumably generate chiral edge currents and the condensate pairs also carry a relative orbital angular momentum (OAM) l z = M ℏ with reduced Planck constant ℏ [50][51][52][53][54][55][56][57][58][59][60][61][62][63][64][65][66].However, while chiral superconductors more generally belong to the class of integer quantum Hall systems [67][68][69], the chiral currents (and OAM) are actually not topologically protected [67,70], unlike in e.g.Chern insulators.Specifically, charge is not a conserved quantity in superconductors, making all chargerelated quantities unprotected, and chiral superconductors also lack a proper Chern-Simons action related to a current-current correlation but instead correspond to a current-density correlation [62,63].As a result, the overall contribution of the edge modes to the OAM and chiral current, i.e. both charge and mass currents, have been shown to be highly variable and non-intrinsic, depending on e.g.boundary conditions, impurities, gap anisotropy, pairing symmetry, and band effects [65,[71][72][73][74][75][76][77][78][79][80][81][82][83][84][85][86].Several of these studies have shown that both the OAM and the total current per edge (i.e.integrated current density) are typically much smaller in chiral d-wave superconductors compared to in chiral p-wave superconductors.In fact, it was shown that with a sharp confining potential in 2D, the OAM and total current even vanish at low temperature in the BCS limit (|∆| ≪ E F with Fermi energy E F ) for all chiral states, except the isotropic chiral p-wave state [73,75].A vanishing current has also been found in certain chiral p-wave systems [76] and in the microscopic limit on e.g. a square lattice [73], although chiral d-wave superconductivity on the honeycomb lattice still generate finite edge charge currents [70].For more details, see Ref. [84] and references therein.This severe or even complete suppression of the edge currents, hampers experimental verification and characterization of proposed chiral d-wave superconductors, since the chiral currents are typically taken as a fingerprint of chiral superconductivity.It therefore becomes important to identify situations where the currents can be generally enhanced, or to find alternative fingerprints [87,88].In this work we focus on the former scenario, showing that the chiral charge currents and orbital magnetic moment arXiv:2308.15258v1[cond-mat.supr-con]29 Aug 2023 (OMM) (i.e. the charged analogue of the OAM) can in fact be enhanced and even have a maximum at low temperatures in the BCS limit, due to mesoscopic finite-size effects.Moreover, with these being quantities related to electric charge and magnetism, they offer straightforward possibilities for experimental detection.
In particular, in this work we study how the chiral charge currents and OMM depend on the system radius R, temperature T , and spontaneous Meissner screening (self-screening) in mesoscopic chiral d-wave superconductors.We find two size regimes dominated by completely different physical behavior.For large system sizes, R > 20ξ 0 , we find a 1/R scaling for the total charge current and OMM, thus reproducing previous results of vanishing current in semi-infinite systems (i.e. as R → ∞) [84].Here, our natural length unit is ξ 0 ≡ ℏv F /2πk B T c , which is the effective superconducting coherence length [89] over which superconducting phenomena typically vary, with Fermi velocity v F (p F ) = |v F (p F )| on the FS and Boltzmann's constant k B .However, in smaller systems R ≤ 20ξ 0 , we find that the spectrum is highly modified by hybridization between opposite edges and current conservation enforcing a faster suppression of the current density as function of distance from the edge.The combination of these effects even produce a sign reversal of the total charge current below R < 10ξ 0 .The existence of these two distinct finite-size regimes lead to a strongly enhanced charge current and OMM with local maxima at R ≈ 20ξ 0 and R ≈ 10ξ 0 , respectively, and also notably with a maximum at low temperatures, which is in stark contrast to the vanishing current previously reported at the lowest temperatures for larger systems [81].We further find that extremely strong self-screening causes a suppression of the charge current and OMM, but also a strong enhancement of the local induced magnetic-flux density which can be measured with magnetometry techniques e.g.via a superconducting quantum interference device (SQUID) [90][91][92][93][94][95][96][97][98][99][100][101][102][103][104][105].Our results thus demonstrate generically finite charge currents and OMM in chiral dwave superconductors enhanced by mesoscopic finite-size effects, specifically edge-edge hybridization and effects of current conservation.As a consequence, mesoscopic superconductors may be a viable route to experimental verification of chiral superconductivity.
The rest of this work is outlined as follows.We present our theoretical model and methods in Sec.II.Section III gives a background with general properties of chiral dwave superconductivity and its order parameter.Section IV gives a general background of chiral edge modes and compares them with flat bands of Andreev bound states (ABS), followed by a study of edge-edge hybridization in finite systems.Section V presents our main results, consisting of the system-size dependence and temperature dependence of the charge-current density, total charge current per edge, and OMM.Section VI studies the effects of self-screening on these quantities for different system sizes.The work is concluded in Sec.VII.Appendices contain further details on our theoretical for-malism, analytic calculations in bulk and at interfaces, and additional numerical results.

II. MODEL AND METHODS
Chiral d-wave superconductivity has been proposed in a number of materials, with widely different atomic structures and properties [20-42, 49, 106-110].In this work, we aim to study properties that are naturally emergent from the chiral d-wave pairing symmetry.This section describes our model and theoretical framework.

A. Model
We are primarily interested in the combined influence of mesoscopic finite-size effects, temperature, and spontaneous Meissner screening on the steady-state charge current and its signature in spin-singlet chiral d-wave superconductors.We therefore appropriately assume equilibrium and spin degeneracy without additional spinorbit interactions, instead leaving such effects as an interesting outlook.Since many proposed chiral d-wave superconductors are layered materials [20][21][22][23][24][25][26][27][28][29][30][31][32][33][34][35][36][37] and some of these are weakly coupled, we consider weak-coupling superconductivity in 2D aligned with the xy-plane in our coordinate system, further modeling a cylindrically symmetric FS such that there is translational invariance along ẑ [111].For clarity, we further assume clean systems with specular surfaces (i.e.perfect specular reflection boundary conditions that conserve the surfaceparallel momentum [112]), but effects of diffuse scattering are also briefly discussed.We do not consider the possibility of boundary-enhanced superconductivity as studied in Refs.[113][114][115][116][117] which is primarily relevant for temperatures T close to the critical temperature T c , while the most interesting effects we find happen in the range T ∈ [0, T c /2].We consider effects of spontaneous Meissner screening of the chiral currents, due to a finite London-penetration depth λ 0 ≡ c 2 / (4πe 2 v 2 F N F ) setting the effective length-scale of magnetic screening and induction, where c is the speed of light, e = −|e| the elementary charge, and N F the normal-state density of states (per spin) at the Fermi level.We assume type-II superconductivity, appropriate for layered materials, thin films, non-elemental materials or most unconventional superconductors.Type-II superconductivity is typically quantified via the Ginzburg-Landau parameter κ ≡ λ 0 /ξ 0 > 1/ √ 2 (or more appropriately via the critical fields) [118][119][120][121].

B. Quasiclassical theory: separation of scales
We focus our attention on the low-energy physics close to the FS, where superconductivity is most important, and note that there is typically a separation in the relevant length-and energy-scales in many superconductors.In particular, the superconducting gap |∆| is usually much smaller than the Fermi energy E F , while the superconducting coherence length ξ 0 is often much larger than both the Fermi wavelength ℏ/p F and the atomic length scale a 0 .In such superconductors, the low-energy (long-wavelength) physics can be separated from the high-energy (short-wavelength) contributions.The quasiclassical theory of superconductivity [112,[122][123][124][125][126][127][128][129][130][131] is a controlled expansion in the resulting small parameters, e.g.|∆|/E F and ℏ/(p F ξ 0 ), and where the leading-order terms are obtained from the low-energy bands close to the FS.Higher-energy contributions can still be inserted into the theory, e.g. through microscopic boundary conditions [132][133][134][135][136][137][138][139].
In quasiclassical theory, spatial variations typically occur on the mesoscopic length scales set by ξ 0 and λ 0 .In particular, the propagation of quasiparticles and superconducting pairs are captured by the quasiclassical propagators g(p F , R; z) and f (p F , R; z), respectively.Here, z is the complex-valued energy of the corresponding propagator.In this work, the retarded propagator g R (p F , R; ε) with z R ≡ ε + iδ (real energy ε and small positive broadening δ) is used to calculate the local density of states (LDOS) with total density of states (DOS) N (ε) = A d 2 RN (R, ε) integrated over the system area A, and where the angle brackets denote the FS average [111] ⟨. . .
which is a line integral in 2D averaging over the Fermimomentum direction.Similarly, the normal-state density of states (per spin) can further be expressed as and the Fermi velocity is parametrized via the FS ac- are used for all other quantities as described in the following, evaluated on the imaginary axis via the Matsubara energies z M ≡ iε n = iπk B T (2n + 1) with temperature T and integer n.In Nambu (particle-hole) space, we construct where "tilde" denotes particle-hole conjugation We briefly drop the arguments (p F , R; z) for a more compact notation and follow the Eilenberger formulation where ĝ is solved from the Eilenberger equation [122] iℏv F • ∇ĝ + zτ 3 − ĥ, ĝ = 0, (5) with the normalization condition ĝ2 = −π 2 τ0 .Here, τi are the 4 × 4 Pauli-spin matrices in Nambu space with identity matrix τ0 , while ĥ are the self-energies, which we divide into diagonal ( Σ) and off-diagonal ( ∆) parts, Here, ∆ is the mean-field superconducting order parameter as described in Sec.II C, while Σ generally includes different interactions and other diagonal self-energies.In this work, Σ captures the electrodynamics, see Sec.II D.

C. Superconductivity
Assuming an even-parity spin-singlet superconductor, the effective pairing interaction V (p F , p ′ F ) is decomposed into the symmetry channels of the corresponding crystallographic point group [140], Here, Γ labels the even-parity spin-singlet irreducible representations, V Γ the pairing strength of the respective symmetry channel, and η Γ (p F ) is the basis function encoding the pairing symmetry on the FS.Similarly, the total superconducting order parameter ∆(p F , R) is written as where each symmetry channel is associated with an order parameter component ∆ Γ (p F , R) with amplitude |∆ Γ (R)| and phase χ Γ (R).We use the efficient "Ozaki summation" [141] based on the Matsubara technique [142][143][144][145][146] to solve the order parameter self-consistently from the superconducting gap equation where the cutoff Ω c is effectively the bandwidth of the pairing interaction [131].
In this work, we study chiral d-wave superconductivity, requiring finite pairing channels Γ ∈ {d where additional pair correlations of other symmetries are automatically included in the theory (e.g.s-wave), while the influence of other attractive pairing channels [147] is left as an outlook.We assume a degenerate pairing with equal V Γ in both d-wave channels which is enforced by symmetry in any material with a three-or sixfold rotationally symmetric lattice [49], thus relevant for many of the recently proposed chiral d-wave superconductors [22][23][24][25][26][27][28][29][30][31][32][33][34][35].We thus consider a total order parameter of the form To better highlight the chiral properties, it is useful to rewrite the superconducting order parameter from Eq. (10) in the eigenbasis η ± (p F ) of the orbital (orb) angular-momentum operator, Lorb z ≡ (ℏ/i)∂ θF , via the linear combination We see that ∆ ± corresponds to a relative phase shift ±π/2 between the components ∆ d x 2 −y 2 and ∆ dxy .There is no loss in generality in this basis change (it is valid even in the non-chiral state, see Sec.III for definition of chiral state).Equating Eq. ( 10) with Eq. ( 11) yields the transformation between the two parametrizations For comparison, we can generalize Eq. ( 12) for a superfluid of angular order M (often referred to as the chiral order), i.e. the Chern number [43][44][45][46][47][48][49], where in 2D even and odd M correspond to spin-singlet and spin-triplet (e.g.M = ±1, ±2, ±3 for p, d, f -wave respectively), ) with eigenvalue l orb = ℏM This corresponds to a condensate which carries an OAM along the ±ẑ-axis [58].In contrast, a non-chiral state, like a nodal or nematic d-wave state [148], is not an eigenstate of Lorb z .We from here on refer to ∆ d x 2 −y 2 and ∆ dxy as the nodal components, and ∆ ± as the chiral components.
Finally, note that apart from applying a start guess to the order parameter, we do not constrain the order parameter components ∆ Γ or their phases χ Γ in any way.Instead, we let them evolve completely independently in the self-consistency.This in principle allows the system to find a non-chiral state, but we always find the chiral state to be stable.

D. Self-consistent electrodynamics
Electrodynamical interactions coupling to the orbital degrees of freedom and caused by external magnetic fields (ext) or magnetic induction (ind) enter the diagonal selfenergies Σ via where A(R) = A ext (R) + A ind (R) is the electromagnetic gauge field.The magnetic flux-density, B(R) = B ext (R) + B ind (R), is related to the vector potential via B(R) = ∇ × A(R).We do not consider external flux in this work, such that B ext = 0 and A ext = 0. Still, A ind is an induced vector potential due to spontaneous charge currents in the system, related via Ampère's law The charge-current density is in turn given by which we express in units of Here j(R) is the total charge-current density with contributions from both quasiparticles and superconducting pairs (hence j is conserved with ∇ • j = 0).The induction gives rise to an additional selfconsistency equation, Eq. ( 16), to be solved together with the self-consistent superconducting gap equation, Eq. ( 9).We note that the induction is effectively weighted by a factor (λ 0 /ξ 0 ) −2 ≪ 1 [149], which for extreme type-II superconductors (λ 0 ≫ ξ 0 ) often makes the induction negligible when λ 0 ≫ R. Still, we always solve both ∆ and A fully self-consistently when λ 0 is finite, using Eqs.(9) and (16).
The orbital magnetic moment (OMM) m = m z ẑ (per 2D layer) is computed from the charge-current density via [112] with natural units m 0 = N ℏ |e| m * = 2µ B N with Bohr magneton µ B = ℏ|e|/2m * , particle number N , and effective quasiparticle mass m * defined via p F = m * v F .We note the analogue between the OMM in Eq. ( 18) and the orbital angular momentum (OAM) L orb = A d 2 R R × j m (R), typically expressed in units L 0 = N ℏ/2, with mass-current density j m [65].
Finally, we introduce the area-averaged induced flux density (or equivalently total induced flux) with flux quantum Φ 0 ≡ hc/2|e| and Planck constant h.

E. Calculations details
We use the Riccati formalism [125,126,128,129,131,135,139] described in Appendix A to solve the Eilenberger equation, Eq. ( 5), numerically, specifically using the open-source framework SuperConga [112] which is free to download from its Gitlab repository [150] with extensive documentation [151].See Ref. [112] for extensive implementation details.The Eqs. ( 5), ( 9) and ( 16) are solved self-consistently in an iterative process by starting from an appropriate start guess, and proceeding until the global error ϵ G of quantity O i at itera- R), A(R), j(R), Ω} and tolerance ϵ tol .In the present work, we set the tolerance to ϵ tol = 10 −7 , use a spatial resolution of 20 discrete points per coherence length, and use an energy cutoff Ω c ≳ 100k B T c .We parametrize the FS with 256 discrete points, except at low temperature or when computing the LDOS, in which cases we use 512 and 720 discrete points, respectively.We find these values sufficient, as we do not notice any difference with finer resolution.
In order to investigate the influence of finite geometry, we simulate discs with radii R ∈ [5,200]ξ 0 , temperatures T ∈ [0.01, 0.99]T c , and penetration depths λ 0 via κ ≡ λ 0 /ξ 0 ∈ [1, ∞).We keep the same spatial resolution when changing R, and keep all parameters fixed during the process of converging towards the self-consistent solution.In order to avoid divergences in the LDOS calculation, we use a small, phenomenological, energy broadening δ ∼ 0.03k B T c in all the LDOS calculations.
In order to improve the understanding and interpretation of our numerical results, we provide supplementary analytic calculations in Appendices B and C for bulk and surface scenarios, respectively.

III. BACKGROUND: CHIRAL d-WAVE SUPERCONDUCTIVITY
To set the stage for our studies of the chiral edge modes and charge currents in the subsequent sections, we here review the basic properties of a chiral d-wave order parameter in a finite system.
A chiral d-wave superconductor is characterized by the existence of two degenerate ground states of opposite chiralities, here denoted ∆ + (p F , R) and ∆ − (p F , R) defined in Eqs. ( 11)- (13).At the onset of chiral superconductivity, the bulk superconductor spontaneously chooses one of the two chiral states (∆ ± ) as the dominant order, while the other state (∆ ∓ ) becomes subdominant and vanishes asymptotically [152].Hence, the order parameter takes the approximate form ∆(p F , R) ≈ ∆ ± (p F , R) in these domains, which leads to a fully gapped state [1][2][3][4].Since the two states ∆ + and ∆ − are related via time reversal symmetry, the appearance of a dominant chirality leads to spontaneous time-reversal symmetry-breaking [43,153].The overall global features of the superconductor is determined by this dominant component.The subdominant component can still carry additional information about local features as it is induced within a few coherence lengths close to spatial inhomogeneities, such as interfaces, defects and vortices [154,155].Hence, in  the presence of such inhomogeneities, the order parameter is not in a pure chiral state, and takes the more general form (11).In Fig. 1, we demonstrate the behavior of the chiral d-wave order parameter solved self-consistently as described in Sec.II, in a mesoscopic system shaped like a disc with dominant component ∆ − , but where the subdominant component ∆ + arises close to the edge.For comparison, we show results with both parametrizations (a,b) ∆ = ∆ d x 2 −y 2 + ∆ dxy and (c,d) ∆ = ∆ + + ∆ − , while (e,f) show the spatial dependence of the order parameter amplitudes at different temperatures, starting at the edge (y = 0) where we find |∆ + | = |∆ − |.This degeneracy is related to the surface suppression of one of the nodal components (∆ dxy at a [100] interface) due to its sign change, leading to surface bound states that locally enhances the other nodal component (∆ d x 2 −y 2 ) [65].We study edge states in detail in Sec.IV.
We note that bulk chiral superconductivity recovers over a characteristic and effective length scale r 0 ≈ 6-10ξ 0 (depending on e.g.temperature), typical for many inhomogeneous superconducting phenomena in clean and mesoscopic systems [149,156,157].We trace the length scale r 0 as naturally emerging from the typical spatial dependence for the ballistic propagators in non-uniform environments.For example, at a [100] interface with translational invariance along x (or equivalently large radius of curvature compared to ξ 0 ), this spatial dependence is captured by the appearance of factors exp(−y/ξ y ) in the quasiparticle propagators and pair propagators.Here, ξ y denotes the effective "non-constant coherence length" [157], which we derive analytically in Appendix C1 as where v F,y (p F ) is the y-component (i.e.perpendicular to the surface) of the Fermi velocity v F (p F ).To compute most quantities, the propagators are averaged over the FS (as described in Sec.II), which means that the factor exp(−y/ξ y ) essentially capture multiple confinement scales ξ y (p F ; ε n ; T ) [65], and illustrates how the boundaries couple the real-space ballistic trajectories to the momentum-space dependence on the FS, as well as the temperature-and energy-dependence.In averaging the propagators over the FS and summing over Matsubara energies, the temperature-dependent decay length r 0 over which bulk superconductivity recovers emerges through the exponential factors.However, due to the non-trivial averages and sums, a direct relation between r 0 and ξ y seems intractable, and we therefore treat r 0 as an emergent length scale which ultimately depends on the superconducting coherence length via ξ y .Finally, we note that as the system radius R becomes comparable to r 0 , there is a strong wave function overlap between opposite pairbreaking surfaces, consequently causing significant reduction of both nodal d-wave components in the entire system, see Appendix D 1 where we include additional results explicitly showing how the order parameter changes with R. In Sec.IV we show that this overlap is also associated with an edge-edge hybridization of the spectrum and the chiral edge modes.Thus, these are mesoscopic finite-size effects that strongly influence the physics in mesoscopically sized superconductors, as we illustrate throughout this work.

IV. CHIRAL EDGE MODES
Chiral superfluids are known to host chiral edge modes (Weyl fermions [13,68]) which are related to the bulk Chern number via the bulk-boundary correspondence [3,4,11].To provide further background for the mesoscopic finite-size effects on the chiral currents, we begin this section by revisiting the well-studied problem of the chiral spectrum and explore the edge modes in a chiral dwave superconductor.Specifically, we focus on the temperature dependence of the LDOS and point out clear differences against surface-bound states in a regular dwave superconductor.Finally, we compute the LDOS for different system sizes and illustrate how the edge modes are influenced by mesoscopic finite-size effects, specifically edge-edge hybridization.

A. Edge modes and bulk LDOS
A chiral superfluid with Chern number M has |M | topologically protected edge modes [43][44][45][46][47][48][49].In Fig. 2(a) we show the subgap spectrum at a semi-infinite surface of a chiral d-wave superconductor (derived analytically without self-consistency in Appendix C 2, while all subsequent results in the main text are numeric and fully self-consistent).There are |M | = 2 chiral branches with dispersion ε(p ∥ ) as a function of the surface-parallel momentum p ∥ [70].The broken time-reversal symmetry of the chiral state is explicitly seen, since for each occupied state ε(p ∥ ), there is no time-reversed partner at ε(−p ∥ ) [43,153].Figure 2(b) shows the fully selfconsistent LDOS at the center (x, y) = (R, R) and edge (x, y) = (R, 0) in a disc with radius R = 30ξ 0 , at different temperatures.In the center, a bulk-like state The bulk state is fully gapped, see also Appendix B for analytic derivation of the bulk gap.At the edge, the chiral edge modes generate a nearly constant subgap LDOS, see also the zoom of the edge LDOS in Fig. 2(c).Finally in Fig. 2(d) we show that at zero energy, the edge modes decay into the bulk over the same characteristic distance r 0 ∼ 6-10ξ 0 as the order parameter variations in Sec.III.Note that at higher temperatures, the superconducting gap is sup-  pressed due to thermal excitations, both at the edge and in the bulk.This leads to a larger effective coherence length, Eq. ( 20), and the edge modes consequently reach further into the bulk at higher temperatures.

B. Comparison with nodal d-wave state
In contrast to the full gap of the bulk chiral d-wave state, a regular d-wave superconductor has a gapless nodal spectrum and can host surface ABS [156,[158][159][160], the latter we also include in Fig. 2(b) for comparison (gray dash-dotted line).There exist several important differences between the ABS in a nodal d-wave superconductor, and the chiral edge modes in a chiral d-wave superconductor.First, the ABS correspond to a flat band (in k-space) [161], while the chiral edge modes are dispersive [70].Second, the ABS peak typically appears as a Lorentzian function centered around the Fermi level and with a width influenced by e.g.boundary conditions and typical energy-broadening effects [156,[162][163][164], while the chiral edge mode has a nearly constant subgap LDOS.Third, we find that the zero-energy LDOS of the ABS is roughly two orders of magnitude larger than the chiral edge mode LDOS in Fig. 2(b).Fourth, the chiral LDOS is largely independent of the surface orientation, while the ABS in a d x 2 −y 2 -wave superconductor appear mainly at the pairbreaking [110]-interfaces but are absent at [100] interfaces [156].Fifth, the ABS are highly degenerate and thermodynamically unstable towards spontaneous symmetry breaking [165][166][167][168][169][170][171][172][173][174], possibly triggering a phase transition into a state known as a "phase crystal" [175] at a relatively high temperature T * ∼ 0.2-0.5Tc [175][176][177][178][179][180][181][182].We find that the low spectral weight in the chiral d-wave system in Fig. 2(b) does not facilitate such a thermodynamic instability.Finally, we point out that the ABS have a significant paramagnetic response [86,[183][184][185][186][187], which we find can be very different from the response of the chiral edge modes (not shown).

C. Mesoscopic finite-size effects: hybridization
Next, we investigate how the finite size influences the subgap LDOS, and in particular the chiral edge modes.Figure 3 shows the subgap LDOS versus energy and coordinate along the diameter (top row) and line cuts through the same data at fixed coordinates (middle row) and at fixed energies (bottom row).Different columns correspond to different disc radii R. Based on these results we find two different regimes with significantly different behavior: R > 20ξ 0 , and R < 20ξ 0 .
For all radii R > 20ξ 0 , the edge modes behave the same as in Fig. 3(a-c): the edge modes decay into the bulk over the same characteristic length scale r 0 ≈ 6-10ξ 0 as the order parameter in Figs.1(i-j).We note that this is similar to the spatial dependence of ABS at pairbreaking [110] interfaces of a nodal d x 2 −y 2 -wave superconductor [156,179].However, Fig. 3(c) shows that the spatial dependence is largely independent of subgap energy, which stands in contrast to the ABS.The overall spatial profile of the edge mode LDOS is smooth except for some very small oscillation, but we find that these oscillations reduce monotonically with higher temperature and larger system size, see additional plots in Appendix D 2.
Figures 3(d-l) show that as the disc radius shrinks below R < 20ξ 0 , significant finite-size effects develop, causing a qualitatively different and highly non-trivial spatial dependence.At R = 15ξ 0 , small kinks start to develop in the LDOS, see Fig. 3(d-f).These kinks develop into full peak structures as R < 10ξ 0 , see Figs. 3(g-l).We also find that the typical decay length of the subgap states decreases with R, reaching ∼ 2-3ξ 0 at R = 5ξ 0 (depending on energy), see Figs. 3(j,l).Hence, the subgap states are "compressed" to a smaller spatial region, and in the process obtain a higher magnitude.Furthermore, as R decreases, we find both that the largest subgap LDOS peak moves slightly away from the edge and that the gap close to the edge reduces, see in particular Fig. 3(j-l).
We interpret the finite-size results in Fig. 3 to be due to hybridization between edge states on opposite edges, and due to resonances between the condensate and the edge arising when the system size becomes comparable with the effective coherence length [188].These hybridization effects are naturally enhanced at smaller system sizes, but also with higher temperatures due to a larger effective coherence length.

V. SPONTANEOUS CHARGE CURRENTS
In the previous section IV, we showed how mesoscopic finite-size effects cause a strong local enhancement of the edge mode LDOS (see Fig. 3).Here, we study the consequent effects on the chiral charge currents.Specifically, the occupation of the dispersive edge modes leads to a finite charge-current density.It might be tempting to conjecture that since the two chiral branches have roughly the same group velocity they will have a constructive contribution to the total charge current.However, such a simple analysis based on the group velocity can have little bearing on the real situation, and the branches can in fact have the opposite contributions [81,84].Consequently, the charge-current density changes direction close to the interface such that the total integrated charge current (and OAM) can vanish.This occurs in particular for all higher-momentum chiral superfluids beyond p-wave in the BCS limit and at low temperatures [73,81], in a microscopic square lattice [81], but also in some chiral pwave systems [76].The cancellation can be understood in terms of the multiple confinement scales and the current response of surface subgap states versus the condensate [65], and the contribution from paired versus unpaired fermions [75], see also Ref. [84].In contrast, we here show that the mesoscopic finite-size effects studied in the previous sections completely modify the charge-current density, causing instead a strong enhancement of both the total charge current (per edge) and the OMM.

A. Charge-current density
We start by analyzing the charge-current density in a chiral d-wave superconductor.In Fig. 4(a) we plot the azimuthal (surface-parallel) component of the chargecurrent density j x (y) versus perpendicular distance y from the surface, at different temperatures in a system with radius R = 20ξ 0 .We find an overall monotonic decrease of the magnitude with higher temperatures, which we attribute to thermal suppression of superconductivity.We find two different signs of the charge current as a function of distance from the edge (barely visible at T = 0.9T c ), with an overall spatial decay into the bulk over the same characteristic length scale r 0 ≈ 6-10ξ 0 as the order parameter and the edge modes in previous sections.In Figs.4(b,c) we analyze the dependence on system size R by showing the spatial dependence of the charge-current density for R ≥ 20ξ 0 and R ≤ 20ξ 0 , respectively.In both regimes, there is always a sign change.The current-density magnitude also increases monotonically with decreasing disc radius in agreement with the larger LDOS of the dispersive low-energy edge modes (see Fig. 3).Apart from these similarities between large and small systems, the exact spatial dependence varies between the two regimes.
In larger systems R ≥ 20ξ 0 shown in Fig. 4(b), the charge-current density varies smoothly with the distance from the edge, with the sign change occuring at roughly the same distance y ≈ 2.7ξ 0 for all R. We find that the overall decay length from the edge remains roughly the same as that of the chiral edge modes (r 0 ).As the disc radius increases, the charge-current density reduces monotonically, but seems to saturate around R ≈ 150ξ 0 .This is similar to the saturation of the small oscillations in the LDOS which also occurs at R ≈ 150ξ 0 , see the discussion in Sec.IV C. We note that at R = 20ξ 0 (R = 150ξ 0 ), the portion of positive and negative charge-current density is very dissimilar (similar) implying a large (small) integrated charge current, discussed further in Sec.V B.
In smaller systems R ≤ 20ξ 0 shown in Fig. 4(c), a reduction of the radius leads to a notable increase in the magnitude of the charge-current density, as well as the sign change occurring at smaller distances, approaching y ≈ 1ξ 0 at R = 5ξ 0 .We point out two important factors contributing to this behavior.The first is the spatial decay of the edge modes themselves.Specifically, as we established in Sec.IV C, a smaller disc radius effectively leads to a "compression" of the dispersive edge modes to a smaller region resulting in a higher LDOS magnitude, due to edge-edge hybridization.The occupation of this locally increased LDOS of dispersive states naturally leads to a locally increased charge-current density.The second important factor is that current conservation enforces the charge-current density to sum to zero across the whole system diameter, such that j x (y) is antisymmetric across the disc center (x, y) = (R, R) where it necessarily goes to zero.For sufficiently small discs this forces j x (y) to decay to zero below the characteristic decay length r 0 .As a consequence of these finite-size effects, the portion of positive current becoming smaller but with a larger maximum magnitude, while the portion of negative current increases significantly.
To further supplement the above data, we provide in Appendix D 3 the analogue plots of Fig. 4(a) for different R, corroborating that the above analysis holds qualitatively at different temperatures.
In summary, the charge-current density j(R) in large systems R > 20ξ 0 shows the same decay length r 0 as the edge modes.As R increases, the spatial profile of j(R) slowly approaches an asymptotic form.In small systems R < 20ξ 0 , j(R) is strongly modified with a shorter decay length but with a much larger magnitude, which we relate to a similar behavior of shorter decay length and larger LDOS magnitude of the dispersive edge modes.

B. Total charge current and magnetic moment
We showed in Fig. 4 that the charge-current density j(R) is finite on the length scale of the coherence length.Experimentally resolving signatures of charge currents on such a small scale is already extremely challenging.More feasible techniques are based on measuring e.g.magnetic signatures of a total charge current averaged over some area.However, the sign-change in the charge-current density reduces such averages.Consequently, the exact balancing of the positive versus negative portions of the charge-current density and its modification by the mesoscopic finite-size effects become important for experimental signatures of the total charge current, and thus for the experimental verification of chiral superconductivity.In this sub-section, we therefore consider the net charge current contribution per edge [189], i.e. the total charge current (sheet charge current) as a function of disc radius R and compare against the semi-infinite scenario (R → ∞) [81].We also consider the OMM m z [Eq.(18)] which is averaged over the entire superconducting sample.Figures 5 and 6 show the total charge current I x and OMM m z , respectively, at λ 0 = ∞.In both figures, (a) and (b) show the temperature dependence at fixed system size R for the two regimes R ≥ 20ξ 0 and R ≤ 20ξ 0 , respectively, while (c) shows the R dependence at fixed temperature T .
In large systems R > 20ξ 0 shown in Figs.5(a) and 6(a), increasing R monotonically reduces both I x and m z at all temperatures T , and the reduction shows a slow asymptotic behavior above R ≈ 150ξ 0 .The reduction is explained by a similar reduction in the magnitude of the charge-current density j(R) with larger R, and more importantly, that the positive and negative portions of j(R) become similar in area thus reducing the spatial integral, see Fig. 4(b).The slow asymptotic behavior in I x and m z is in turn understood from the similar slow asymptotic behavior in the exact spatial form of j(R) above R ≈ 150ξ 0 .Next, we note that I x has a maximum at T ≈ 0.4T c and a minimum at T → 0 which are both in agreement with calculations in semi-infinite systems [81].Similarly, m z also seems to develop a slight minimum at low temperature, but only for the largest R. The minima in I x and m z tends slowly to zero as R → ∞, see analysis further below.Most importantly though, both quantities show a dramatic enhancement close to R = 20ξ 0 which is in complete contrast to the semi-infinite scenario, which we explain by the strongly enhanced positive portion of j(R) [Fig.4(b)], in turn due to the strong local enhancement of the LDOS of the dispersive edge modes.
In small systems R < 20ξ 0 shown in Figs.5(b) and 6(b), I x has a maximum at R ≈ 15-20ξ 0 while m z has a maximum at R ≈ 10ξ 0 .Below these maxima, both quantities decrease monotonically with smaller R.However, this decrease is not described by a decrease in the charge-current density j(R) since it in fact grows monotonically with smaller R, but is instead completely described by the area becoming similar for its positive and negative portions, see Fig. 4(c).Interestingly, the negative portion of j(R) eventually obtains a larger area than the positive portion and starts dominating at low temperatures below R < 7ξ 0 , causing a complete reversal of the total charge current, thus suddenly flowing in the opposite direction to the chirality.In contrast, m z shows no such reversal.These different behaviors in I x 0.0 0.2 0.4 0.6 0.8 1.0 T /T c 0 0.02 0.04  and m z come from their different integrands in Eqs. ( 21) and (18), specifically m z always weighs the positive portion of the charge-current density at the edge stronger due to the integrand R × j(R) (note that R is the c.m. coordinate with origin in the disc center).
Figures 5(c) and 6(c) show I x and m z , respectively, as functions of R to directly illustrate the scaling behavior.We find that a fit to f fit (R) = c 0 + c 1 /R describes the asymptotic behavior for large R rather well, where c 0 and c 1 are fit parameters.We choose such a simple model to reduce over-fitting and we verify that fitting only a few data points 30ξ 0 ≤ R ≤ 80ξ 0 (gray region) gives a good extrapolation to all other data points in the interval R ≥ 20ξ 0 .See Appendix E for further details about the fitting procedure.Importantly, we find that as T → 0 and R → ∞, both I x and m z tend to very small values, albeit very slowly.We naturally do not find an exact zero, since our smallest temperature is finite T = 0.02T c , we have a finite maximum size R = 200ξ 0 , our fit model is simple, and due to finite numerical accuracy.Furthermore, we point out that our disc-shaped system has a curved surface in contrast to the straight semi-infinite interface.Thus, even at R ≈ 150ξ 0 , the surface curves significantly over the distance of a few coherence lengths.Next, in the opposite limit as R reduces towards the small regime R ≤ 20ξ 0 , the 1/R dependence leads to a significantly increased I x and m z due to an increase in the positive portion of j(R), but then a reduction in I x and m z as the negative portion of j(R) starts increasing.As explained previously, this can be directly traced back to the edge-edge hybridization of the dispersive edge modes.Eventually these mesoscopic finite-size effects also cause a breakdown of the 1/R scaling in I x and m z .
In summary, we find a 1/R scaling with increased system size in both the total charge current I x and OMM m z .The scaling enhances both quantities in smaller systems, but eventually the scaling breaks down due to the mesoscopic finite-size effects, specifically the edge-edge hybridization and effects of current conservation, causing a local maximum around R ≈ 20ξ 0 .As a consequence and in contrast to the semi-infinte systems where the current is minimal at low temperatures, we find that both the total charge current and OMM can be large at low temperatures in small systems, pointing to mesoscopic systems as a potential platform to experimentally measure signatures of the chiral currents.

VI. SPONTANEOUS MAGNETIC INDUCTION AND MEISSNER SCREENING
Any finite charge-current density j(R), i.e. like the one we showed in Figs. 4, might induce a magneticflux density (induction) according to Eq. ( 16) which couples back to the superconductor and generates additional screening currents (self-screening).This spontaneously induced flux serves as a potential experimental observable, measurable with e.g.nano-SQUIDS and other magnetometry setups [90][91][92][93][94][95][96][97][98][99][100][101][102][103][104][105].Generally, the induction and screening scale inversely with the penetration depth λ 0 [112,149], tending to zero in the limit κ ≡ λ 0 /ξ 0 → ∞ (the limit studied so far), but may become important when λ 0 < R. It is reasonable to expect that most proposed chiral d-wave superconductors fall in the former limit, since most non-elemental or unconventional superconductors are strongly type-II, implying a relatively large penetration depth λ 0 , as is also the case for most thin films and layered materials [118][119][120][121].However, for sake of completeness, we in this section study the influence of finite penetration depth in the whole type-II interval λ 0 ∈ [1, ∞)ξ 0 , since the penetration depth also depends on materials properties, e.g.typically increasing (decreasing) with the concentration of non-magnetic (magnetic) impurities in the system [121].
We begin by summarizing that we find only a small variation in the order parameter profile and in the LDOS for different values of λ 0 (see Appendices D 1 and D 2, respectively for analysis).We find that the charge-current density j(R) is unmodified for λ 0 > R but changes significantly as λ 0 ≪ R consequently obtaining an additional sign change (see Appendix D 3).The total charge current I x and OMM m z show a stronger dependence with λ 0 as they depend more crucially on the exact balancing between positive and negative portions of j(R).This we illustrate in Figs.7 (a) and (b), respectively, while (c) shows the total induced flux Φ ind defined in Eq. (19).Quite generally, the ratio between the penetration depth λ 0 and system size R essentially determines the importance of the screening [190].For λ 0 ≫ R, we find that the system is poorly screened, with both the total charge current I x and OMM m z reaching the asymptotic value of zero screening (λ 0 → ∞).For λ 0 < R, the screening becomes stronger, and the charge-current density re-organizes to minimize the total charge current in the system.In particular, the chiral charge currents lead to a finite superfluid momentum, which in turn generates a finite kinetic energy.The Meissner screening acts to reduce this kinetic energy.Thus, as λ 0 becomes small, both I x and m z become vanishingly small.This means that in strongly screened systems, the total charge current and OMM is difficult to measure.On the other hand, the charge-current density might still be finite locally for such small λ 0 (see Appendix D 3 for further analysis), leading to an induced magneticflux density B ind via Ampère's law in Eq. ( 16).Rewriting this equation on a dimensionless form in our nat-   19)] (c), versus finite penetration depth λ0 at T = 0.1Tc.Note the semi-log scale in (c) in units of flux quantum Φ0 ≡ hc/2|e|.Line colors denote system radius R, while dashed lines denote asymptotes of infinite penetration depth (zero screening).
To summarize, we find that mesoscopic finite-size effects dramatically increase different experimental observables that can be used for verifying chiral superconductivity, specifically the total charge current I x and OMM m z for large λ 0 .This scenario is reasonable to expect, since λ 0 is expected to be large in most proposed chiral d-wave superconductors.Even if λ 0 is for some reason small, our results show that the induced flux Φ ind instead becomes large.

VII. CONCLUDING REMARKS
Chiral superconductors are characterized by dispersive chiral edge modes that are presumed to generate spontaneous currents close to interfaces and surfaces.However, with charge not being a conserved quantity in superconductors, the charge current is not topologically protected.
The charge-current density j(R) may be finite on the superconducting coherence length scale, ξ 0 , but previous studies have shown that it in many systems also has a sign change close to the surface in all chiral d-wave superconductors [73-76, 81, 84].This sign change often causes the total charge current per edge, I x (i.e. the integrated charge-current density), to vanish, thereby inhibiting experimental verification of chiral superconductivity.
In this work, we revisit the issue of edge charge currents in chiral d-wave superconductors, focusing on mesoscopic sizes relevant for quantum devices, and also computing the OMM m z .We study the dependence of I x and m z on the system radius R, temperature T , and spontaneous Meissner screening due to a finite penetration depth λ 0 .We show that there are essentially two regimes with qualitatively different behavior, namely large systems R > 20ξ 0 and small systems R < 20ξ 0 , with the border inbetween comparable to r 0 ≈ 10ξ 0 , the characteristic decay (recovery) length of edge (bulk) quantities, e.g. the order parameter, edge modes, and charge-current density.
For large systems R > 20ξ 0 , we show that the total charge current and OMM scale as 1/R, recovering previous results of vanishingly small total current as R → ∞ [73,81,84].In smaller systems we find that the 1/R scaling causes a dramatic increase in I x and m z , but this scaling eventually breaks down below R < 20ξ 0 due to mesoscopic finite-size effects becoming dominant.These effects correspond to current conservation forcing the current to decay faster than r 0 , as well as edge-edge hybridization causing the dispersive edge modes to compress to a smaller region but with a much larger LDOS locally.As a result, the charge-current density j(R) shows a similar shorter decay length and also with a much larger magnitude locally as R reduces.Notably, the portions of positive and negative signs in j(R) re-balance, such that the positive portion closest to the surface shrinks, while the negative portion further away grows.This causes a local maximum to appear in I x at R ≈ 15-20ξ 0 and in m z at R ≈ 10ξ 0 , since they are both described by a spatial integral of j(R).Importantly, these local maxima survive in the limit of zero temperature, which is in contrast to semi-infinite systems where instead I x vanishes in the same setup [81].For even smaller systems R ≲ r 0 , the re-balancing of the portions of opposite signs in j(R) eventually causes a sign change in the total charge current, i.e. the net current spontaneously changes direction due to the mesoscopic finite-size effects, flowing in the opposite direction to the chirality.A reversal of the edge charge current in a chiral superconductor has previously been shown to occur e.g.due to FS nesting [72], surface disorder [79,86] or contributions from non-edge states [191], while we here show that it also occurs due to mesoscopic finite-size effects.Furthermore, we show that screening due to finite penetration depth λ 0 leads to suppression of the charge current and OMM if ξ 0 ∼ λ 0 ≪ R, but we still find an increase in the spontaneously induced magnetic flux Φ ind , measurable with magnetometry [90][91][92][93][94][95][96][97][98][99][100][101][102][103][104][105].Thus, our results show that mesoscopic finite-size effects significantly enhance the total charge current I x and OMM m z in systems with large λ 0 (relative to the system size), while the induced magnetic flux Φ ind is still significantly large even in systems with small λ 0 .As a consequence, our results highlight finite mesoscopic systems as a promising platform for enhancing experimental signatures of chiral superconductivity.This Appendix introduces analytic solutions to a bulk chiral d-wave superconductor in the same model as described in Sec.II.These results are used for comparison in the main text against the self-consistent numeric results, as well as in the subsequent Appendix C against analytic results in a semi-infinite system.
In a uniform environment, ∆(p F , R) = ∆(p F ) and v F • ∇γ = 0.The scalar Riccati equations (A7)-(A8) then take the polynomial form From these equations, the homogeneous (h) solutions to the Riccati equations are obtained as where For a chiral d-wave order parameter, we note that In contrast to the numeric calculations where both amplitudes are allowed to evolve independently, we here assume equal amplitudes appropriate for a degenerate bulk environment Inserting the above solutions and definitions into Eqs.(A9)-(A10) yields the bulk propagators noting in particular that g 0 (p F ; z) = g 0 (z) is independent of p F (i.e.isotropic).We next use the bulk propagators from Eq. (B12) and (B13) to derive the DOS, the current, as well as the lowtemperature gap and free energy.These quantities are used for comparison and as scales throughout the main text.The LDOS is given by Eq. ( 1) in terms of the causal (retarded in time) propagator with z R = ε + iδ, where δ is a small positive broadening.Inserting the homogeneous propagator from Eq. (B12) and taking the trivial FS average, we get the DOS i.e. same as a fully gapped s-wave superconductor.
Next for the charge-current density Eq. ( 17), we note that the homogeneous Matsubara propagator with z = iε n is given by which is purely imaginary.The current effectively depends on the real part of this propagator [157] and is therefore zero.Hence, there are no currents in a uniform chiral d-wave superconductor The breaking of additional symmetries, e.g.translational symmetry by an interface, can induce chiral edge modes and a finite chiral charge-current density j ̸ = 0.
We next turn our attention to the pair propagator, which can be used to solve the low-temperature gap equation in a uniform environment.At low temperatures, we denote the amplitude |∆| → ∆ 0 , and get the analytic solution which is the same as for an s-wave superconductor [190] and where γ E is the Euler-Mascheroni constant.The lowtemperature order parameter is consequently with |∆(T = 0, p F )| = ∆ 0 constant across the whole FS.Finally, we consider the bulk free energy, i.e. the BCS condensation energy.At low temperature, it can be derived from the Luttinger-Ward functional [124,[192][193][194] and we find it to be which again is the same as for an s-wave superconductor [190].Here, V = d 3 R is the volume of the sample.
In summary, we have seen that in a uniform environment, the properties of the chiral d-wave superconductor are very similar to that of a conventional s-wave superconductor.

Appendix C: Analytic solutions at vacuum interfaces and domain walls
In this Appendix, we present analytic solutions to the quasiclassical propagators close to interfaces, taking different boundary conditions into consideration.These solutions are used to interpret the spatial dependence of the numeric results in the main text, in particular in terms of the effective coherence length.Furthermore, we show that perfect specular reflection leads to chiral edge modes and that such a surface is equivalent to a transparent domain wall.Surface retro-reflection and back scattering, on the other hand, lead to a reduced weight of the chiral edge modes.Additionally, the solutions provide a starting point for further analytic studies of the properties of chiral d-wave superconductors, such as the thermodynamics, and the spectrum and its full contribution to the angular momentum and equilibrium currents.

Inhomogeneous coherence functions and effective coherence length
We start by deriving the form of the coherence functions in an inhomogeneous environment and show the natural appearance of the effective coherence length discussed in the main text.
Consider a semi-infinite superconductor, as depicted in Fig. 8(a).For simplicity, we let the interface normal y be aligned with the high-symmetry axes of the crystal lattice, but note that the orientation does not influence the physics for a chiral order parameter.This is in contrast to e.g. a nodal d-wave superconductor, where [110] and [100] interfaces are distinctly different [156].Assuming that any spatial variations in the order parameter can be approximated by a piecewise constant function, the inhomogeneous Riccati equations (A7)-(A8) can be solved analytically.In the absence of vorticity and other inhomogeneities along the interface coordinate x, there is translational invariance ∂ x γ = ∂ x γ = 0.The only inhomogeneity is then introduced by the boundary, and the inhomogeneous solutions take the explicit form with integration constants C and C, and where the homogeneous solutions γ h and γh , as well as Ω, are given by Eqs.(B3)-(B5).We note that even in the absence of translational invariance along the interface, the solutions take this form, but with modified exponential factors.
Here, ξ y gives the effective superconducting coherence along y (i.e.perpendicular to the interface) which can diverge e.g. at the coherence peaks.We determine the integration constants C and C by imposing continuity wherever there is an inhomogeneitiy, in this case at the boundary of each piecewise constant part in ∆.This yields the relation Γ(p ′ F , y = 0; z) = Γ(p F , y = 0; z) and similar for Γ.Here, ) denotes the trajectory impinging on the boundary, while p F = p F (cos θ F , sin θ F ) with θ F ∈ [0, π) denotes the outgoing trajectory.From here on, we introduce the subscripts "in" and "out" to denote the impinging and scattered trajectories at the interface, with arguments (p ′ F , y = 0; z) and (p F , y = 0; z), respectively.This determines the integration constants where our notation specifies Γ in and Γin via Eqs.(C1) and (C2), and γ h,out and γ h,out from Eqs. (B3) and (B4).
From here on, we assume that the important contributions lie in the boundary condition and quasiparticle scattering, rather than in the exact spatial form of the order parameter, and therefore approximate a uniform order parameter ∆(p F , R) ≈ ∆(p F ).This leads to simplification Γ in = γ h,in and Γin = γh,in .Note that we of course relax this assumptions (and several others) in our numeric simulations.Effects of the spatial dependence of the order parameter on the analytics is discussed e.g. in Ref. [83].We apply the same unitary gauge transformation as in Appendix B, such that again and ∆(p F ) ∆(p F ) = ∆∆ * = |∆| 2 (i.e.spin-singlet).
In the following derivations, we drop the arguments: (p F ) from ∆(p F ), (p F ; z) from Ω(p F ; z), (p F , R; z) from γ(p F , R; z), and similar for tilde quantities.

Specular boundary conditions: reflective vacuum interfaces and transparent domain walls
We next proceed to derive the surface propagators and their spatial dependence at interfaces with perfect specular boundary conditions, followed by expressions for the chiral edge modes and chiral edge currents.These expressions form a base line to compare against fully selfconsistent results in the main text.
Perfect specular scattering at an interface is characterized by the momentum change (p ∥ , p ⊥ ) → (p ∥ , −p ⊥ ), connecting two trajectories with angles θ ′ F + θ F = 2π, as shown in Fig. 8(a).In both the chiral p-wave and d-wave systems, this leads to a sign-change in one of the two order parameter components, while the other is uninfluenced.In our notation, this means that from an incoming to an outgoing trajectory, (∆ 1 , ∆ 2 ) → (∆ 1 , −∆ 2 ), and consequently ∆ → ∆ * .In other words, a quasiparticle undergoing such a scattering sees a chiral inversion, in exactly the same way as if it transmits through a transparent domain wall.The two scenarios are treated analogously.This leads to Inserting these solutions into Eqs.(C4)-(C5) yields the integration constants such that after a bit of algebra, we find the inhomogeneous solutions Using these solutions, we can solve for the surface propagators via Eqs.(A9) and (A10), but note that we have to treat the direction of v F (p F ) carefully.For θ F ∈ [0, π) we get γ s = Γ and γs = γh , while we for θ F ∈ [π, 2π) we get γ s = γ h and γs = Γ.We find that the solutions can be generalized for all θ F by introducing the terms The surface propagators then become which are valid at a specular reflective interface and at a transparent domain wall.We see that far from the interface, where e −y/ξy → 0, the bulk propagators in Eqs.(B12)-(B13) are recovered.Close to the surface (e −y/ξy → 1), the propagators take the form We next take a closer look at the edge modes and currents related to these surface propagators.(C17) This corresponds to bound states at energy ε bs = s∆ 1 which is a propagating chiral mode with momentum p F,x = p F cos(θ F ).The spectrum of the chiral edge modes can essentially be derived via the momentum-resolved LDOS, i.e. omitting the angular average in Eq. ( 1).This yields the spectrum plotted in Fig. 2(a), thus with two edge modes for s = ±1.

b. Chiral surface currents
In equilibrium, the Fermi-Dirac distribution together with the chiral edge modes lead to the occupation of a state with a net momentum, breaking time-reversal symmetry and generating a charge-current density.We next study this charge-current density, defined by Eq. ( 17), in terms of the Matsubara propagators g M (p F , R; ε n ) = g 0 (p F , R; z = iε n ) with Matsubara energies ε n = πk B T (2n + 1).We find that the real part of the Matsubara quasiparticle propagator is such that the chiral edge current is given by Here, the overall sign is set by that of ∆ 2 , i.e. the chirality directly sets the direction of the current.The above current shows a sign change at ∼ 1ξ 0 , similar to the fully self-consistent numerical calculations in Fig. 4.However, the overall magnitude differs significantly, which we interpret to be due to the approximation in the spatial dependence of the order parameter [83], and other nontrivial effects omitted in the analytical but not in the numerical calculations.

Retro-reflection and back scattering
We next study interfaces with retro-reflection and back scattering (defined below), and show that it is an interesting scenario that can be used to further distinguish between spin-singlet and spin-triplet superconductors, e.g. between chiral d-wave and chiral p-wave states.We note that retro-reflection becomes experimentally relevant even in superfluid 3 He [195].
Surface roughness and diffusivity causes a finite probability of back scattering.In the extreme limit of perfect retro-reflection, e.g. at a surface with a mesoscopic saw-tooth profile of perfect reflectors with length L (where a 0 ≪ L ≪ ξ 0 ) as depicted in Fig. 8(b), then (p ∥ , p ⊥ ) → (−p ∥ , −p ⊥ ).A spin-singlet order parameter is symmetric with respect to such a momentum reversal, ∆(−p F ) = ∆(p F ), which leads to Γ in = γ h,out , and consequently C = C = 0, such that Thus, despite the presence of the interface, this yields the same bulk propagators and solutions presented in Appendix B. In contrast, a spin-triplet system is instead antisymmetric ∆(−p F ) → −∆(p F ), and the retro-reflective boundary condition therefore leads a Jackiw-Rebbi zero mode [196] associated with surface pairbreaking and the appearance of flat bands of ABS.Hence, perfect retroreflection technically causes a removal of the chiral edge modes since quasiparticles essentially scatter between the same Chern numbers, but this result is highly unstable as any infinitesimal modification from perfect retroreflection leads to a recovery of the chiral edge modes, see also Ref. [65] for further discussion.Thus, for more realistic back scattering, there is instead a finite suppression of the LDOS of the chiral edge modes, in favor of a bulk (flat-band ABS) LDOS in spin-singlet (spin-triplet) systems.
The above results point to an interesting scenario to distinguish between spin-singlet and spin-triplet order parameters, since the spin-singlet system effectively behaves like a bulk system close to a retro-reflecting interface, while the spin-triplet system leads to order parameter suppression and a large density of dispersionless surface ABS with energies at the Fermi level.

Appendix D: Additional numerical results
This appendix presents additional numerical results for a chiral d-wave superconductor to supplement the analysis in the main text.The first subsection studies the order parameter, the second subsection studies the LDOS of the zero-energy edge modes, while the third subsection studies the charge-current density.

Order parameter
We next study how the spatial dependence of the order parameter is influenced by mesoscopic finite-size effects in Fig. 9 and by finite Meissner screening in Fig. 10, to be compared with the analogous Fig. 1 in the main text.
Figure 9 shows the amplitude of the nodal (chiral) order parameter components in the top (bottom) row, as a function of perpendicular distance y from a [100] interface, in a disc-shaped superconductor with radius R, dominant chirality ∆ − , and penetration depth λ 0 = ∞.The first (second) column shows large (small) discs at low temperature T = 0.1T c , while the third (fourth) column shows large (small) discs at high temperature T = 0.9T c .The nodal component ∆ d x 2 −y 2 (∆ dxy ) is enhanced (suppressed) due to fermionic bound states at the surface [65], but recover to their bulk values over a characteristic length r 0 , i.e. the same as in the main text, which is also the decay length of the bound states.Similarly, the dominant (subdominant) chiral component ∆ − (∆ + ) is suppressed (enhanced) at the edge, but also recover to their bulk value over the length scale r 0 .In large systems R > 20ξ 0 , the order parameter profiles remain the same for all disc radii, with the same r 0 ≈ 6-10ξ 0 (depending on e.g.temperature).As R becomes comparable with r 0 , the influence of opposite edges becomes important, leading to hybridization as discussed in Sec.IV C. Therefore, in small systems R < 20ξ 0 , the order parameter profiles vary strongly with disc radius, with the typical decay length reducing monotonically to ∼ 2ξ 0 at R = 5ξ 0 .Finally, we note that increased temperatures lead to an overall suppression of the order parameter due to thermal excitations, which leads to a longer effective coherence length as discussed in Sec.III.Hence, finite-size effects are significantly enhanced at higher temperatures.This leads to a complete destruction of superconductivity at T = 0.9T c in the smallest disc with R = 5ξ 0 .Thus, the above results show a strong dependence of the exact order parameter profile with system size.
Figure 10 shows the spatial dependence of the order parameter amplitude in a disc with radius R = 10ξ 0 , for different λ 0 , illustrating negligible influence for all λ 0 ∈ [1, ∞)ξ 0 at both high and low temperatures.The influence is even smaller for larger R (not shown).Hence, the self-screening effects due to small λ 0 discussed in Sec.VI does not correspond to a modified order parameter amplitude, but rather the subgap states and the exact shape of the charge-current density as discussed next.

Local density of states
We next study how the spatial dependence of the zeroenergy LDOS approaches an asymptotic form for large disc radius R in Fig. 11, and how the spatial dependence is influenced by finite Meissner screening in Fig. 12, to be compared with e.g. the bottom row in Fig. 3 in the main text.
Figure 11 plots the LDOS at zero energy as a function of distance from the edge, for different disc radii R, at fixed T = 0.5T c and λ 0 = ∞.Importantly, the spatial decay length r 0 ≈ 6-10ξ 0 is roughly the same in all large systems R > 20ξ 0 and there are no strong kinks or additional peak structures like in the small systems R ≤ 20ξ 0 .Hence, the overall spatial dependence is quite smooth, apart from a small oscillations.These oscillations reduce significantly with larger R and subsides almost completely at R ≈ 150ξ 0 .We note that this is the size where the charge-current density reaches an asymptotic value, see Fig. 4(b).However, it is not clear at the present if this similar asymptotic behavior with large R in the edge modes and charge-current density are related, or if this is more of a coincidence.
Figure 12 shows the spatial dependence of the zeroenergy edge modes at T = 0.1T c but for different penetration depths, in a large disc R = 30ξ 0 (a), and a small disc R = 10ξ 0 (b).In the large disc, different λ 0 mainly influence the small oscillations in the LDOS, but the edge mode decay is otherwise not influenced.In the small disc, a reduction of λ 0 causes a smaller LDOS magnitude and slight increase in the decay length from the edge.Since the LDOS corresponds to dispersive edge modes, these results suggest that a stronger screening reduces the currents, which is consistent with the results found in the main text.

Currents
We next study how the spatial dependence of the charge-current density depends on temperature for different disc radii R in Fig. 13, i.e. the direct analogue of Fig. 4(a).We also study how the charge-current density is influenced by finite screening in Fig. 14, supplementing the analysis in Sec.VI.
Figure 13 shows the spatial dependence of the chargecurrent density j(R) at λ 0 = ∞ and at different temperatures, where each panel corresponds to a different system radius R. Below R ≤ 10ξ 0 the system radius is smaller than the typical decay length r 0 of the currents.This leads to a significantly increased magnitude due to larger LDOS of dispersive edge modes as discussed in Sec.IV C, and a more complicated spatial dependence with kink-like structures (higher-order harmonics) due to the edge-edge hybridization and effects of current conservation as discussed in Sec.V A. We also note that for large R, the size of the positive and negative portions of j(R) become more similar in size implying a smaller to- Order parameter amplitudes as a of distance from a edge in a disc with radius R and dominant chirality ∆− at λ0 = ∞, for the nodal components (top row) and chiral components (bottom row).First (second) column: low temperature in large (small) discs.Third (fourth) column: high temperature in large (small) discs.tal integrated current I x , while they are very dissimilar due to a larger (smaller) positive portion at R = 20ξ 0 (R = 5ξ 0 ) implying a large positive (negative) current.These effects are most pronounced at low temperatures, which is in agreement with the main text in Sec.V B.
Figure 14 shows the charge-current density versus distance from the edge in the large regime R > 20ξ 0 (a), and the small regime R < 20ξ 0 (b).Different lines correspond to different penetration depths λ 0 .Above λ 0 > R, there is negligible effect of screening (small difference between finite and infinite λ 0 ).As λ 0 reduces below R, the magnitude of the positive portion decreases, while the negative portion increases, which is seen most clearly in Fig. 14(a).Hence, the total current reduces, consistent with a stronger diamagnetic Meissner screening from the bulk condensate.However, as λ 0 becomes comparable with ξ 0 , a qualitatively different behavior develops.The current obtains a more complicated spatial depen- dence with kinks, a smaller negative region, and an additional sign change far from the edge in agreement with semi-infinite systems studied in Ref. [81].This leads to a modified balance for the integrated current, such that it changes sign as shown in the main text in Sec.VI.
Appendix E: Fit of asymptotic behavior for large system size In this appendix, we provide additional details on the fit procedure used in Sec.V B for the total charge current I x in Fig. 5(c) and the OMM m z in Fig. 6(c) as a function of R.
We initially compare different fit functions, specifically exponential decay and various polynomial decays.We however do not find a good agreement with exponential decay and focus in the following on the polynomial function.We limit ourselves to only considering polynomials with a few terms to reduce over-fitting, using the fit function with fit parameters c 0 , c 1 , and c 2 .We find that the fit function which most accurately describe the asymptotic behavior of I x and m z for large R to be the polynomial 1/R (best fit with c 2 ≈ 1 ± 0.1), while the higher-order polynomials (e.g.1/R 2 ) result in very poor fits.Moreover, we verify that fitting only a few data points, e.g.30ξ 0 ≤ R ≤ 80ξ 0 [gray region in Figs.5(c) and 6(c)] , gives a good extrapolation to all other data points in the interval R ≥ 20ξ 0 .Table I summarizes the best fit parameters for I x and m z for fixed c 2 = 1.From these fit parameters, we find that as T → 0 and R → ∞, both I x and m z tends to a very small, negligible, value.An exact zero is of course not found due to finite temperature T /T c = 0.02 and finite numerical accuracy.The numerical accuracy causes some uncertainty in c 0 for each scenario depending on fit region, but c 1 remains more or less the same with c 1 ≈ 0.7 ± 0.1.The fact that both I x and m z have the same limiting behavior for large R is understood from their expressions in Eqs. ( 21) and ( 18): the total current is the integral of j(R) over the radius, while the magnetic moment is the integral over R × j(R) over the area.

Figure 1 .
Figure 1.Ground-state order parameter of a chiral d-wave superconductor with disc radius R = 25ξ0, dominant chirality ∆−, T = 0.1Tc, λ0 = ∞.(a,b) Amplitudes of nodal components, and (c,d) chiral components.(e,f) Temperature dependence of the order parameter amplitudes, as a function of the distance from a [100] edge, along the dashed line (x = 25ξ0) in (a).Temperatures vary from T = 0.1Tc (blue) to T = 0.9Tc (red), in steps of ∆T = 0.1Tc.

Figure 3 .
Figure3.First row: subgap LDOS versus energy ε and coordinate y across the disc diameter at T = 0.1Tc, λ0 = ∞, and different disc radii R (columns), with bulk gap ∆0 ≈ 1.76kBTc.Second row: same but line cuts at fixed y from the disc edge y = 0 (blue) to the disc center y = R (red).Third row: same but line cuts at fixed ε from ε = 0 (purple) to ε = kBTc (orange).

Figure 4 .
Figure 4. (a) Spatial dependence of the surface-parallel component of the charge-current density jx(y) in a disc with radius R = 20ξ0 at λ0 = ∞ with temperatures from T = 0.1Tc (blue) to T = 0.9Tc (red).Horizontal gray line is a guide to the eye marking jx = 0. Inset: heatmap of the current-density magnitude at T = 0.1Tc.Charge-current density jx(y) at T = 0.1Tc in systems with radii R ≥ 20ξ0 (b) and R ≤ 20ξ0 (c).Note the different ranges of the axes.

Figure 5 .
Figure 5.Total charge current Ix [Eq.(21)] versus temperature in large (b) and small (b) systems, and versus system size (c).Markers are calculated data points, solid lines are a guide to the eye, while dashed line is a fit.Here, λ0 = ∞ and I0 ≡ j0ξ0 = ℏ|e|NFv 2 F .

Figure 8 .
Figure 8. Superconductor-vacuum interface aligned with the crystal ab-axes of a semi-infinte superconductor (SC).(a) Translationally invariant interface with perfect specular reflection, connecting incoming (outgoing) quasiparticle momentum p ′ F (pF) according to p ′ F,x = pF,x and p ′ F,y = −pF,y.(b) Interface consisting of mesoscopic perfect retro-reflectors of size L with a0 ≪ L ≪ ξ0, such that surface scattering connects incoming and outgoing quasiparticle momenta according to p ′ F = −pF.

Figure 10 .
Figure 10.Order parameter amplitudes as a function of tance from a[100] edge in a disc with dominant chirality ∆−, for the nodal components (top row), and chiral components (bottom row).First (second) column: low (high) temperature.

Figure 11 .
Figure 11.Spatial dependence of the edge modes at zero energy, T = 0.5Tc and λ0 = ∞, for different disc radii R.

Figure 13 .Figure 14 .
Figure 13.Azimuthal (surface-parallel) component of the charge-current density as a function of distance from the edge in a disc with dominant chirality ∆− at λ0 = ∞.Different panels correspond to different R, with temperatures T = 0.1Tc (blue) to T = 0.9Tc (red).Horizontal line (gray) marks jx = 0. Note different x-and y-axis scales in top versus bottom rows.