Breaking time-reversal and translational symmetry at edges of $d$-wave superconductors: microscopic theory and comparison with quasiclassical theory

We report results of a microscopic calculation of a second-order phase transition into a state breaking time-reversal and translational invariance along pair-breaking edges of $d$-wave superconductors. By solving a tight-binding model through exact diagonalization with the Bogoliubov-de~Gennes method, we find that such a state with current loops having a diameter of a few coherence lengths is energetically favorable below $T^*$ between 10%-20% of $T_{\mathrm{c}}$ of bulk superconductivity, depending on model parameters. This extends our previous studies of such a phase crystal within the quasiclassical theory of superconductivity, and shows that the instability is not qualitatively different when including a more realistic band structure and the fast oscillations on the scale of the Fermi wavelength. Effects of size quantization and Friedel oscillations are not detrimental. We also report on a comparison with quasiclassical theory with the Fermi surfaces extracted from the tight-binding models used in the microscopic calculation. There are quantitative differences in for instance the value of $T^*$ between the different models, but we can explain the predicted transition temperature within each model as due to the different spectral weights of zero-energy Andreev bound states and the resulting gain in free energy by breaking time-reversal and translational invariance below $T^*$.


I. INTRODUCTION
High-temperature superconductors have been shown to have a superconducting order parameter of d-wave symmetry. 1,2 Although this is well established, a lot of research is still ongoing to explain the full phase diagram, and the mechanism of superconductivity is still controversial. 3,4 In addition, there are remaining interesting, and somewhat controversial, questions on the properties of the d-wave superconducting phase, in particular in devices where surfaces and interfaces play an important role. 5,6 Surfaces and interfaces with a normal not exactly aligned with a crystallographic axis are pair breaking, with associated formation of zero-energy Andreev bound states. 7 These states play an important role in determining device physics. They show up in tunneling experiments as zero-bias conductance peaks, 8 and influence the current-phase relation of Josephson junctions. 9,10 From a fundamental physics point of view, it is interesting that the flat band of zero-energy Andreev bound states can be related to bulk topology, 11 but at the same time may lead to instabilities where additional symmetries are broken. 12 Early on it was suggested that timereversal symmetry may be broken by forming a subdominant component of the order parameter, [13][14][15][16] thereby shifting the Andreev states and splitting the zero-bias peak as sometimes seen experimentally. 17,18 Another suggestion is magnetic ordering at the surface causing a spin split. 19,20 Experimentally, breaking of time-reversal symmetry remains controversial. 2,6,[21][22][23][24] The zero-bias conductance peak does not always split at low temperature, and the associated currents and magnetic fields are so far unobserved or small. These works may also become relevant to the ongoing research into Sr 2 RuO 4 , where recent experiments 25,26 point towards a singlet superconductor with unconventional orbital symmetry. 27,28 Recently, we have reported on a different scenario of shifting the Andreev bound states and lowering the free energy in a more complicated manner, where both timereversal and translational invariance along the surface are broken. [29][30][31][32][33] In a second order phase transition, spontaneous current loops become energetically favorable at a temperature T * up to 20% of T c of bulk superconductivity. The inhomogeneous broadening along the surface may explain that the zero-bias conductance peak is not necessarily split in a tunneling experiment, but instead acquires a temperature-independent width for T < T * . In addition, since neighboring current loops have opposite circulation, the magnetic fields tend to cancel when averaged over one period of the current pattern, which is about 10ξ 0 , where ξ 0 is the superconducting coherence length. 29 This unexpected symmetry-broken state does not involve any subdominant pairing channel or any other mean-field order than the superconducting order parameter with pure d-wave symmetry. The free energy is lowered through Doppler shifts of Andreev states. The superflow causing the Doppler shift is associated with the order-parameter phase having an oscillatory behavior along the edge, resembling a crystal pattern that may be called a phase crystal. 33 The phase crystal has a higher T * than the state with translational invariant superflow. [34][35][36][37][38] The previous work 29- 33 and predictions were carried out using the quasiclassical theory of superconductivity. 39 The quasiclassical approximation involves integrating out effects relevant at the atomic scale. This requires a good separation between the atomic scale and the relevant superconducting scale, i.e. the superconducting coherence length. Typical high-temperature superconductors have very short coherence lengths, and the validity of the quasiclassical approximation can be questioned. A first comparison of results from a tight-binding Bogoliubov-deGennes theory with those of quasiclassical theory 40 were made analysing the unconventional charging effects in small superconducting YBCO single-electron tunneling devices 6 . In this paper we explore the phase crystal at edges of d-wave superconductors within a tight-binding model. We are interested in the effects of including the atomic scale oscillations, as well as the effects of the more realistic band structure and Fermi surface taken into account in the tight-binding model. So far 29-33 only a circular Fermi surface was used in the quasiclassical calculations. We will in this paper also extract the relevant Fermi surfaces predicted by the tight-binding models and compare with quasiclassical theory.
The paper is structured as follows. In Section II we introduce the real-space tight-binding model and how it is solved using the Bogoliubov-de Gennes exact diagonalization method. This section also includes a reciprocal space calculation to characterize the parameter spaces, set key model parameters, and extract the relevant Fermi surfaces. In Section III we present results obtained within the tight-binding model for the transition into the phase crystal. In Section IV we introduce a comparison between the tight-binding model and the quasiclassical formulation with the extracted Fermi surfaces, thereby highlighting the qualitative similarities but quantitative differences between the two approaches. Finally in Section V we summarize. A few calculations have been collected in the Appendix.

II. TIGHT-BINDING DESCRIPTION OF A D-WAVE SUPERCONDUCTOR
A. Real space formulation The normal state of the material is described by a single band, where the bandstructure is determined by a few hopping integrals t ij . The single-particle Hamiltonian readŝ where t ij includes nearest (t), next-nearest (t ), and next-next-nearest neighbor (t ) hopping parameters, see Fig. 1. The indices i and j enumerate the lattice sites, σ labels spin, and δ ij is the Kronecker delta-function. The operatorĉ jσ annihilates an electron with spin σ on site j. The chemical potential µ is set by the doping level, but is here treated as a parameter of the model. We take the hopping parameters from literature 41,42 where they have been extracted either from relevant experimental data or from density functional theory. The Hamiltonian describing weak coupling d-wave superconductivity on a lattice iŝ where the mean field order parameter lives on nearest neighbor links where V is the coupling constant. In the bulk, ∆ ij are positive for links along the crystallographicâ-axis, R i − R j = ±aâ, and negative for links along theb-axis, where a is the lattice constant. We note that depending on the parameters of the model (see next subsection), an extended s-wave order may be stable instead of the d-wave, in which case the link order parameter ∆ ij is positive in both crystallographic directions. In this paper, we will focus on the part of parameter space where d-wave order is stable in the bulk. After performing a Bogoliubov rotation, the Bogoliubov-de Gennes (BdG) equation for a superconductor is obtained where E n is the eigenvalue and u T is the corresponding eigenfunction. Since we are considering a singlet superconductor, we are here suppressing the spin indices but consider both positive and negative energy eigenvalues, thereby avoiding double counting. Note that if the u-amplitudes are spin up, then the v-amplitudes are spin down, see e.g. Ref. 43. The most straightforward strategy is to resort to direct diagonalization and obtain all eigenvalues and eigenvectors. The only complication is that the link order parameter must be computed selfconsistently through the gap equation where T is the temperature. Note that in this paper we set the reduced Planck constant , the Boltzmann constant k B , and the electron charge e, all to unity. Once self-consistency of ∆ ij has been achieved we may compute observables, such as currents and local density of states. The local current density is in the tight-binding model defined on links coupling nodes through the hopping integrals t ij in Eq. (1). The current from site j to site i is computed through the formula We use as an extra convergence test that the currents flowing into and out of all nodes in the grain are conserved.
The local density of states at position R j is computed as where η > 0 is a small imaginary part of the energy. The total density of states of the grain is obtained by summation over sites.
Once we have the full eigenvalue spectrum we can also study the thermodynamic properties. The free energy is given as 44 while the entropy has the well-known form (9) We also use the convergence of the free energy as a check for our solutions. In principle, Eq. (8) contains an additional term, that stems from the internal energy E g . This term gives, when non-zero, the same contribution in both the normal and in the superconducting states and is omitted, since we will present results for the free energy difference between the superconducting and normal states.

Notes on site versus link quantities
In the tight-binding model it is important to keep in mind that some quantities are defined on lattice sites, while others are defined on the links. For instance, when plotting the current it is sometimes convenient to define a vector field with arrows residing on the lattice sites. Let us introduce the following notation for the link current from site R j to site R i = R j + δ, which is illustrated in Fig. 1 for a nearest neighbor link δ = aâ. When we in this way single out site R j , we see that current will flow along links to neighbors R i for which the hopping integrals t ij in Eq.
(1) are non-zero. Positive and negative values means that current flows out of or into site j. Current conservation requires Consider now the following vector field whereδ = δ/|δ| are unit vectors. This vector field gives a nice overview of current flow patterns. However, it is not fulfilling current conservation 45 and is in a strict sense an unphysical quantity. Nevertheless, we will for convenience present also this vector field.
In the same way, in many works 46-50 the d-wave order parameter is presented as a site quantity where the signature function s d (δ) is equal to +1 for links along theâ-axis and −1 for links along theb-axis, and ∆ j (δ) is defined in analogy with Eq. (11). A caveat is that since the coupling constant V for superconductivity resides on the links, it may favor an extended s-wave order instead. That would correspond to summing with plus signs in both directions In the bulk with translational invariance, these two signatures reflect the two representations B 1 and A 1 of the crystal lattice. In a finite size system, however, both necessarily coexist which simply reflects that there is no translational invariance and the system is inhomogeneous across the grain, for instance at the corners. Although the nodal quantity ∆ s (R j ) necessarily is non-zero at all temperatures below T c , it is always small for the parameter space studied here. We do not consider this a subdominant order parameter component in the sense used in continuum models, 15 where there is a second coupling constant.
In Section III we will find spontaneous currents and a complex valued link order parameter, ∆ ij = |∆ ij | exp [iχ ij ], with the superconducting phase χ ij oscillating along the edge. It would then be natural to apply a gauge transform to make the order parameter real and extract the superfluid momentum, given by the gradient of the phase 31 . However, the gauge transform is not well defined on a finite size lattice because the number of link phases χ ij is not enough to uniquely compute the node phases needed [51][52][53] to define the gauge transform. This problem does not appear for an infinite lattice [51][52][53] and is peculiar to lattices with edges. See Ref. 54 for a similar problem appearing when attempting to convert between the link order parameter and the quantities in Eqs. (14)- (15). In conclusion, we will in this paper focus on the self-consistently calculated link order parameter only. When showing results, the order parameter ∆ ij is assigned to the midpoint R = (R i + R j )/2 between nodes, as indicated in Fig. 1.

B. Reciprocal space formulation: analysis of parameter spaces
Before presenting results of the real space calculation it is useful to extract material parameters for an extended system both in the normal and in the superconducting states. In addition, we extract information about the Fermi surfaces that we will use in Section IV to compare the tight-binding and quasiclassical results.

Normal state
The normal state dispersion for a single band is where t, t , and t are the hopping integrals from Eq. (1). In the following we use t > 0 as the natural unit of energy in the tight-binding model. The chemical potential is set by µ, and we introduce We can then find the Fermi surface, i.e. the set of k F that satisfies ξ kF = 0. This has to be done numerically.
The velocity is defined as The expression can straightforwardly be calculated analytically. The Fermi velocity is then obtained as v kF using the calculated k F . We will concentrate on two sets of parameters, 41,42 as summarised in Table I. The normal state characteristics of the two bandstructures are shown in Fig. 2. In Fig. 2(a) and Fig. 2(b) we show the bandstructures as density plots, with the Fermi surfaces indicated by black lines. The Fermi velocities vary around the Fermi surfaces as indicated by the arrows. Clearly, parameter set #1 leads to an almost circular Fermi surface, while parameter set #2 gives a more realistic model of a typical high-temperature superconductor.
The normal state density of states is obtained through where N k is the number of k-points we include in the first Brillouin zone (1.BZ). The resulting bulk density of states for the two models are shown in Fig. 2(c). We outline in Appendix A how these Fermi surface parameters are fed into the quasiclassical calculations.

Superconducting state
Turning to the bulk superconducting state, we focus on the d-wave, i.e. the B 1 link order parameter. The pairing interaction and order parameter are where the d-wave orbital basis function is The normalization is chosen as The temperature-dependent gap equation takes the form where E k = ξ 2 k + ∆ 2 k . The superconducting coherence length can be defined in different ways, but we will use where The parameters describing the superconducting state for the different tight-binding realisations are listed in Table II. For the model parameters we study in this paper, the bulk order parameter symmetry is d-wave, see the variation of T c with µ for the two tight-binding models in Fig. 2(d).
The A 1 extended s-wave channel has a basis function Y s (k) = cos(k x a) + cos(k y a) and is stable for negative values of µ in Fig. 2(d). If we would consider this order parameter symmetry as well in the comparison with a quasiclassical calculation, we should feed in the strength of that subdominant order through its bulk T c . We see in Fig. 2(d) that it is zero in model #1 and ∼ 10 −3 in model #2 and can therefore be neglected. Finally, we utilize the relevant Fermi surface basis function Y d (k F ) in the quasiclassical calculations in Section IV.

III. RESULTS: TIGHT-BINDING MODEL
We have numerically solved the Bogoliubov-de Gennes equation, Eq. (4), by direct diagonalization to extract eigenvalues and eigenvectors for each guess of the linkorder parameter in Eq. (5). The procedure is iterative and continues until self-consistency between the order parameter and the eigenvalues and eigenvectors has been achieved. We have studied a range of grain sizes and shapes. For larger systems we have also used the corresponding Green's function formalism with the same tight-binding Hamiltonians. In this case, we use both the recursive method and the Chebyshev polynomial expansion method 55 . For the recursive method, we use our own implementation 56 of the knitting algorithm 57 . All these methods give the same final result.
A subtle detail in the numerics is the initial guess for the order parameter that is needed to find the minimum of the free energy. By guessing a purely real order parameter, the metastable phase without spontaneous currents is always found. We have utilized a few different strategies in order to stabilize the regular pattern of circulating currents. One is to use the ansatz for the phase in Eq. (26) below. Another one is to include for the first few iterations an on-site s-wave order parameter with a phase winding along the edge and then throw it away. The link order parameter then picks up the phase oscillations. Yet another possiblity is to have a region near the edge where the link order parameter is guessed to have a Larkin-Ovchinnikov 58 type of amplitude oscillation, but with equal magnitudes of the real and imaginary parts. These guesses give different paths to the free energy minimum. In fact, depending on the period of the guessed oscillations, one may stabilize different numbers of current loops in the grain. But there is only one solution that has minimal free energy. We note that a random seed of a complex part of the order parameter results in a lot of noise in the system and it is much harder to find the correct free energy minimum. To summarize, it is important to carefully check for good convergence of the order parameter such that the free energy is minimized and the current is conserved. For the results we show here, current conservation as in Eq. (12) is upheld to a relative accuracy ∼ 10 −8 at all individual sites.
We show the link order parameter ∆ ij (R) for a triangular grain with 12621 sites in Fig. 3. Results for tightbinding models #1 and #2 are displayed in the left and right columns, respectively. The grain has a single pair breaking [110] edge at y = 0. At the [110] edge the wellknown zero-energy states are formed 7-9 . In the tightbinding model they show up as a large number of eigenvalues at zero energy, see red squares in Fig. 4(a)-(b). The corresponding eigenvectors have large amplitudes at the pair breaking edge, and the order parameter is suppressed there. This flat band of zero-energy states is energetically unfavorable, and as discussed in the introduction a phase transition can be induced at T * where spontaneous current loops appear at the edge. Such currents create phase gradients and superflow that Doppler shifts the zero-energy states away from zero energy thereby lowering the free energy, see also Section III A below. The shifts in the energy eigenvalues are also clearly seen in Fig. 4(a)-(b), black diamonds and indigo triangles. The resulting total density of states of the grain are shown for both models in Fig. 4(c)-(d). The large zero energy peak for T > T * due to the flat band of zero-energy eigenvalues is clearly broadened for T < T * . At the same time, the high energy spectra, including the oscillations due to the finite size of the grain, remain unchanged when lowering the temperature from above to below T * .
For T < T * the order parameter phase χ ij acquires oscillations along the [110] edge with a period ∼ 50a ∼ 10ξ i.e. a few superconducting coherence lengths to fit a pair of counter-flowing loop currents, see Fig. 3(e)-(f). We introduce the wavevector 1/q x ∼ ξ of this oscillation to adapt to the notation in Ref. 33, where a variational ansatz was introduced for the phase oscillations near T * within quasiclassical theory. It is and fits to the tight-binding results are plotted as black dash-dotted lines in Fig. 3(e)-(f). Near the edge, the amplitude of the phase oscillations is rather large, of the order π, but there is no phase winding and no topological defects in the order parameter. The phase decays to zero away from the edge on the scale y 0 ∼ ξ. The amplitude |∆ ij (R)| also acquires a small oscillation parallel to the edge on the same scale 1/q x , although this effect is more pronounced in tight-binding model #1 than in tight-binding model #2. The additional fast oscillations on the scale of the lattice constant a are inherent to the tight-binding model. There are also oscillations of the amplitude on the scale of the inverse Fermi wave vector 1/k F ∼ 3a to 4a when looking at the amplitude as a function of distance from the edges. These also appear at the non-pairbreaking edges, 59,60 see indigo dash-dotted lines in Fig. 3(c)-(d). These oscillations are a result of Friedel oscillations at the edge and are tem- perature independent.

A. Thermodynamics
We present in Fig. 5 the low-temperature thermodynamics of a square d-wave superconducting grain with all four edges at a 45 • angle from the ab-axes. This grain has 10255 sites and is slightly smaller than the one in Fig. 3. In particular, each edge is much shorter than in the triangular grain. The free energy of the metastable phase, which has no currents and a purely real order parameter [red open squares in Fig. 5(a)-(b)], shows a pronounced upturn at low temperature, signalling the energy cost of the zero-energy Andreev bound states. For larger grains, the upturn is less visible since the weight of the edge is diminished compared with the grain interior. It is interesting that we can stabilize different number of current loops along each edge: i) four current loops, purple triangles, and ii) two current loops, black diamonds. The two configurations have the same T * . For lower temperature, T < T * , the configuration with four current loops (purple triangles) has lower free energy. From the entropy, Fig. 5(c)-(d), one can see from the abrupt change of slope at T * that the phase transition is of second order and has the same value within our numerical accuracy irrespective of number of current loops. To enhance the visibility of the knee in the entropy as function of temperature, we present the difference in entropy between the phase with currents and the metastable phase which has no currents, S S − S ms .
The difference in T * between the two tight-binding models can be understood as due to the different number of zero-energy eigenvalues. In Fig. 4 we see that there are more eigenvalues and a higher spectral weight in tightbinding model #1. The Doppler shifts then become more energetically favorable, and consequently, T * is higher, see Table II. The higher energy cost of the zero-energy states in model #1 is also visible in the larger upturn in the free energy of the metastable state, see Fig. 5(a)-(b).
Thus, we can conclude that the phase crystallization that was studied in Refs. 29-33 within the quasiclassical approximation can be found also in the more general mean field tight-binding Bogliubov-de Gennes theory. Early work 49,50 on this type of tight-binding model also found structures in the order parameter and spontaneous currents. At that time only small grains could be studied, probably because of limited computational resources, and the regular amplitude and phase oscillations that we find here were not seen. In those works, the appearance of edge currents was associated with the development of the A 1 site quantity in Eq. (15). Here, instead, we associate the spontaneous currents with phase crystallization and identify the resulting Doppler shifts as the microscopic mechanism behind the lowering of the free energy below the phase transition temperature T * . We also note that, at least from quasiclassical theory, a subdominant s-wave order parameter at the edge would favor a translational invariant current flow along each edge, see Ref. 29. In addition, for a subdominant s-wave order parameter we would expect a spectral gap around zero energy, which is not present in Fig. 4. To further support these conclusions, we compare in the next section with results that we have obtained with the quasiclassical theory.

IV. RESULTS: COMPARISON WITH QUASICLASSICAL THEORY
The quasiclassical theory of superconductivity 39,61 is used to study triangular grains, following the methods described in Refs. 29-33, but with modified Fermi-surface averages 62 as explained in App. A. It is noted that Fermi-surface effects were studied 63,64 for a similar spontaneous-TRSB phase 65,66 . In addition to the Fermi surfaces defined by tight-binding models #1 and #2, a circular Fermi surface is also investigated. It shows very similar results to model #1 due to the similar Fermi surface shape, and is therefore omitted from some of the figures. We begin by comparing the self-consistent orderparameter and current obtained with quasiclassics versus with BdG. This is followed by an analysis of the thermodynamics and the phase transition, where numerical and analytical calculations are combined to explain the results in terms of the Fermi-surface features.
A. Self-consistent observables As the system enters the symmetry-breaking phase, all three models follow Eq. (26), with a phase that varies sinusoidally with wavenumber q x ∼ 1/ξ 0 parallel to the interface, and that decays exponentially over distance y 0 ∼ ξ 0 away from the interface. Here, we use the definition ξ 0 ≡ |v F | FS /(2πT c ). Such a superflow gives rise to smooth and round currents 29,31,33 , very similar to those found with the BdG approach. As the temperature is lowered far below the transition temperature, however, the shape of the phase is significantly modified, from sinusoidal to triangle-wave-like. This ensures long and constant phase gradients, and hence constant superflow p s , which leads to the most efficient Doppler shift of midgap states. This is seen in the deviation from the fit to Eq. (26) in Fig. 6, where the order-parameter amplitude and phase are plotted at T = 0.1T c < T * , with results for models #1 and #2 displayed in the left and right columns, respectively. The same deviation from Eq. (26) seems to be somewhat present in the BdG results as well, but is highly obscured by the atomic-scale oscillations parallel to the interface (see Fig. 3). The resulting currents are qualitatively very similar in quasiclassics and BdG, see Fig. 7.
It is noted that the size of the unit cell illustrated in the BdG results is roughly 50a × 30a, which with a coherence length of ξ ∼ 5a corresponds to the unit cell of roughly 10ξ 0 × 6ξ 0 illustrated in the quasiclassical results. Furthermore, the same peculiar pattern of edge-defects (i.e. sources, sinks and saddle-points) in the superflow vector field p s as reported in Ref. 31 are still present in all Fermi surface models (not shown here).

B. Thermodynamics of the phase transition
The pair-breaking interface leads to zero-energy midgap states (MGS) with an associated energy cost, in the form lost condensation energy, that scales as 1/T . In the presence of a superflow p s , the zero-energy states are Doppler shifted to finite energies δ ∝ v F · p s , consequently reducing the energy cost. However, the super- flow also shifts the continuum states in the bulk, where it instead costs kinetic energy. The optimal form of the superflow is therefore an exponential decay away from the pair-breaking interface 67 , as in Eq. (26). Due to the 1/T -dependence of the MGS energy cost, the energy gain of spontaneous superflow at the interface eventually outweighs the cost of its tail in the bulk, below a transition temperature T * . This transition temperature is extracted from the self-consistent numerics via the quasiclassical free energy 68 , by comparing the free energy of systems with and without spontaneous superflow. The latter is referred to as meta-stable (ms) as it exhibits a higher free energy below T * . Table III presents the numerically obtained T * and average spectral weight (A) for all three Fermi-surface models, together with evaluated low-temperature analytic expressions of the gap and free-energy terms (explained further in App. B). Model #1 (#2) shows a lower (higher) T * than that of a circular Fermi surface, which correlates with a lower (higher) cost of the MGS, and a lower (higher) energy gain of Doppler-shifting them. The density of states also follows this trend, with a lower (higher) spectral weight for model #1 (#2) above T * . This is because there are less (more) zero-energy states at the interface, that ex- Here, the depairing current is defined as j d ≡ 4πTc|e|NFvF.
tend over a shorter (longer) distance into the bulk due to a smaller (larger) suppression of the order-parameter amplitude |∆|, see Fig. 8 (a). The spectral weight is obtained by integrating the local density of states, N (R; ), in a small window around zero energy, and spatially averaging over one unit cell (u.c.) of the periodic phase (i.e. the region shown in Fig. 6), according to where A u.c. = u.c. dR. Figure 8 shows that in the presence of superflow below T * , the Doppler shift is lower (higher) for tight-binding model #1 (#2) compared to the circular Fermi surface, leading to a higher (lower) spectral weight close to zero energy. To summarize, model #2 has the highest spectral weight of zero-energy states, and the most efficient Doppler shift, yielding the highest T * . Note that this is in contrast to the BdG results, where model #1 has the highest spectral weight and T * .
The quasiclassical results are further understood by analyzing the dependence on the Fermi velocity, v F (k F ), in the angular averages entering the observables (explained in App. A). The Fermi velocity is constant across the circular Fermi surface, while it varies for models #1 and #2, with maxima and minima which either align with the nodes or the antinodes (lobes) of the d x 2 −y 2wave basis function, as is seen in Fig. 2(a)-(b). Model #1 has its Fermi velocity maxima along the nodes, and its minima along the antinodes, while it is the opposite 69 for model #2. The full analysis is non-trivial, but in short, this behaviour ensures that model #2 and #1 weighs the MGS contribution higher and lower, respectively (see App. B for further details). This trend was verified by also investigating the cost of the midgap states for a d xywave basis function with a pair-breaking [100] interface (not shown here). Here, the Fermi-surface shows the opposite features, such that the opposite trend is expected for T * . Indeed, in this case, model #1 and #2 have instead a larger and smaller transition temperature than the circular Fermi surface, respectively.
To conclude, the quasiclassical results thus highlight that the zero-energy states and the transition temperature can change quite significantly with the Fermi surface TABLE III. Quasiclassical results obtained with full selfconsistent numerics (T * and A), and obtained by evaluating analytic low-temperature expressions (∆ and δΩ, see App. B), for the three different Fermi surface models. Here, A and Ams are the spectral-weights with superflow and in a meta-stable (ms) system without superflow, respectively, integrated over one superflow period (∼ 10ξ0 × 6ξ0) at the interface following Eq. (27), and evaluated at temperature T = 0.1Tc < T * . The spectral energy range is /Tc ∈ [−5η, 5η], where η = 0.05 is the smearing (i.e. magnitude of the small imaginary part added to the energy when computing the LDOS), and A0 = TcNF. For the analytic results, ∆ bulk 0 denotes the bulk gap, δΩ BCS the BCS free energy (i.e. without any pair-breaking interface), δΩ MGS the energy-cost of the midgap states, and δΩ DS the energy-gain of Doppler shifting the midgap states to finite energies. Finally, the energies and superflow are given in units ΩA ≡ ANFT 2 c , ΩL ≡ Lξ0NFT 2 c and p0 ≡ Tc/vF, with A ≡ dxdy, L = dx andvF ≡ |vF| FS. parameters, and that it is possible to pinpoint the dependence on these parameters for various observables and quantities. In the tight-binding results, there are additional microscopic effects complicating the situation, but the end result is the same in the sense that the transition is governed by the spectral weight of zero-energy states.

V. SUMMARY
In summary, we have studied two tight-binding models of d-wave superconducting grains with pair-breaking [110] edges. At a phase-transition temperature T * , the order parameter develops phase gradients, and circulating currents appear near the edges. This extends our earlier results 29-33 within the quasiclassical approximation to a more general theory with a realistic Fermi surface and where fast oscillations on the scale of the Fermi wavelength are taken into account. We find that the phase transition and the qualitative characteristics of the symmetry-broken phase are universal, independent of including the microscopic details or not. On the other hand, some quantitative details, such as the predicted T * , depend on the model. For instance, T * of tight-binding model #2 is lower than #1 within BdG, while it is higher within the quasiclassical calculation. This deviation can be due to several things: e.g. electron-hole asymmetry not included in quasiclassics, or details of the eigenfunctions and Friedel oscillations included in BdG. But as in the symmetry-broken phase. The average is taken within one superflow period (∼ 10ξ0 × 6ξ0) at the interface, as shown in Fig. 6 (a)-(b). The gray lines mark the integration interval of the spectral weight, Eq. (27). The inset of (a) shows the zero-energy peak within this interval.
we have shown we can still understand the variation of T * in all cases (two models within both BdG and quasiclassics) in terms of the spectral weight of zero-energy Andreev bound states. Higher spectral weight leads to a higher T * , since the gain in free energy due to the Doppler shifts then increases. Within quasiclassics, the variation of the Fermi velocity around the Fermi surface is also important. The universality can be understood since the scale of variations of the current patterns and the phase variations is a few coherence lengths, which is much larger than the fast 1/k F oscillations. When fast oscillations, and the high energy parts, are integrated out when making the quasiclassical approximation, the physics behind the second-order phase transition at T * is kept. Since the tight-binding model is more general than quasiclassical theory, we may assume that the predictions would be more reliable. On the other hand, the tight-binding model we have studied also has its approximations, including the weak-coupling mean field approximation. As an outlook, it would therefore be of interest to go beyond mean-field and include effects of electron correlations. 70 Within the quasiclassical theory of superconductivity, 39,61 the calculation of an observable or a self-energy involves a Fermi surface average. In our case of a two-dimensional system the average is over the Fermi lines in Fig. 2. Each of them defines a contour C. The average of an arbitrary function f (k F ) then takes the form of a line integral The total density of states at the Fermi energy is Let us extend the bandstructures in Fig. 2 to a repeated zone scheme, choose origo at k = (π, π), and use polar coordinates. For the parameters of our tight-binding models, the contour then becomes a closed non-circular loop that can simply be parameterized by the azimuth ϕ ∈ (−π, π]. The elementary arc length is and the line integral takes the form Thus, we extract k F (ϕ) numerically from the bandstructure of the chosen tight-binding model. The Fermi velocity v F (ϕ) and the derivative k F (ϕ) can be computed from analytic formulas with k F as input. 62 We may then compute the total density of states at the Fermi energy N F . Note that when viewing the Fermi surface in the first Brillouin zone, as in Fig. 2, we must translate the extracted Fermi momenta k F back from higher Brillouin zones.
To compute observables and self-energies, we need to compute the quasiclassical Green's functionĝ, see the Methods section in Ref. 32. That involves solving firstorder differential equations along trajectories defined by the extracted Fermi velocities v F (ϕ), see Eqs. (13)- (14) in Ref. 32. It is then important to note that the Fermi momentum k F (ϕ) and the Fermi velocity v F (ϕ) are not parallel. At the starting points and end points of the trajectories, i.e. at the edges, trajectories couple through a boundary condition. We use a specular boundary condition, and limit ourselves to grains with high-symmetry edges, either along crystal axes or 45 • rotated, like the [110]-edge in Fig. 1. This means that a trajectory parameterized by ϕ couples to the same trajectories ϕ as for a circular Fermi surface.
Appendix B: Analytic quasiclassical expressions at low temperatures To get an analytic handle on the energetics of the midgap Andreev states, analytic expressions are presented for a bulk d-wave superconductor at low temperatures (derived in Ref. 72). These expressions are used to obtain the gap and the energies in Tab. III. We use a gague transformation that eliminates the phase of the order parameter in favor of explicit dependence on superflow field p s 31,33 , hence rendering the order parameter real. In the absence of superflow, the bulk order parameter at zero temperature can be expressed as where γ E is the Euler-Mascheroni constant. For a d x 2 −y 2 order parameter with circular Fermi surface, the FS integral is taken analytically and reduces to the familiar result ∆ 0 /T c = √ 2π exp(−γ E − 1/2) ≈ 1.51. For general tight-binding parameters, the FS average takes the form of Eq. (A1), and is evaluated numerically for e.g. tightbinding models #1 and #2, again, see Tab. III. Having calculated the gap, the BCS free energy is where A = dxdy. Assuming an infinite and maximally pair-breaking interface with negligible order-parameter suppression, the cost of the midgap states is given by where L = dx,v F ≡ |v F | FS , andn =ŷ is the surface normal. For the circular Fermi surface, this reduces to δΩ MGS /(Lξ 0 N F T c ) = π( √ 8/3)∆ 0 . Introducing a small and homogeneous superflow p s = p sx parallel to the interface, the midgap states are Doppler shifted by an energy δ ∝ v F · p s , leading to an over-all energy gain Again, this can be taken analytically for the circular Fermi surface, yielding the dimensionless expression δΩ DS /(Lξ 0 N F T 2 c ) = −π |p s | /p 0 , where p 0 ≡ T c /v F is a natural scale for the superflow. Note that the energy gain in Eq. (B4) is linear in |p s |, while the energy cost of bulk superflow scales as |p s | 2 , with the pair-breaking (pb) superflow of the order p pb s = ∆ 0 /v F , see e.g. Refs. 72 and 73. Finally, the angular dependence of the integrands in Eqs. (B2)-(B4) are shown together with the Fermi velocity, v F , in Fig. 9. The figure shows that having a large Fermi velocity along the antinodes is conducive to large midgap-state energy costs, as well as large energy gains of Doppler-shifting those states, within the quasiclassical theory of superconductivity. The integrands are defined such that Ω X = 2π 0 I X dϕ, and the quantities ΩA and ΩL are the same as in Table III. In this choice of coordinate system the basis-function antinodes (lobes) are situated along ϕ = (n + 1/2)π/2, and the nodes along ϕ = nπ/2, with n ∈ Z. The interface normal is parallel to ϕ = 0.