Phase Crystals

Superconductivity owes its properties to the phase of the electron pair condensate that breaks the $U(1)$ symmetry. In the most traditional ground state, the phase is uniform and rigid. The normal state can be unstable towards special inhomogeneous superconducting states: the Abrikosov vortex state, and the Fulde-Ferrell-Larkin-Ovchinnikov state. Here we show that the phase-uniform superconducting state can go into a fundamentally different and more ordered non-uniform ground state, that we denote as a phase crystal. The new state breaks translational invariance through formation of a spatially periodic modulation of the phase, manifested by unusual superflow patterns and circulating currents, that also break time-reversal symmetry. We list the general conditions needed for realization of phase crystals. Using microscopic theory we then derive an analytic expression for the superfluid density tensor for the case of a non-uniform environment in a semi-infinite superconductor. We demonstrate how the surface quasiparticle states enter the superfluid density and identify phase crystallization as the main player in several previous numerical observations in unconventional superconductors, and predict existence of a similar phenomenon in superconductor-ferromagnetic structures. This analytic approach provides a new unifying aspect for the exploration of boundary-induced quasiparticles and collective excitations in superconductors. More generally, we trace the origin of phase crystallization to non-local properties of the gradient energy, which implies existence of similar pattern-forming instabilities in many other contexts.


I. INTRODUCTION
The defining characteristic of superfluidity and superconductivity is spontaneous symmetry breaking of the global U(1) phase χ , associated with the order parameter = | | exp(iχ ). The phase and its spatial variations give rise to phenomena of importance for technological applications, such as type-II superconductivity where Abrikosov vortices are formed in an external magnetic field and in Josephson junctions [1]. Within the BCS paradigm [2], a uniform fixed value of the phase is directly tied to the finite amplitude | | of the macroscopic Cooper-pair wave function. If the phase is nonuniform, by Galilean invariance it results in superflow with superfluid velocity and momentum mv s = p s (R) = (h/2)∇χ (R), where m is the electron mass andh is the reduced Planck constant. Such phase variations and the associated condensate currents cost gradient energy where the gradient energy coefficient k > 0 should be computed from microscopic theory. A physical picture emerges where the phase is rigid and coherent over macroscopic distances and the superconducting state is stable. Thus, it would be surprising if there existed a more ordered state with a softer phase and spontaneous superflow with energy gain F sf < 0.
Here we propose that under certain conditions there exists a low-temperature superconducting state where the rigid phase acquires structure by breaking translational invariance. In this state, which we refer to as a phase crystalline state, a periodic pattern with wave vector q is formed where A q (R ⊥ ) is a function of coordinates orthogonal to q. The additional order parameter in the phase crystal is the finite Fourier amplitude C q . The superconducting ground state with spatially oscillating phase also breaks time-reversal symmetry and sustains a nontrivial periodic superflow pattern and circulating currents j(R), as illustrated in Fig. 1(a). Similar current patterns have been found in numerical work on mesoscopic grains of d-wave superconductors [3], and the unusual superflow field p s (R) was recently analyzed [4]. Here we establish that the physical origin of this surface state is phase crystallization. Breaking of continuous translational symmetry is particularly striking. Its reduction to discrete translations gives a multitude of crystals [5] and ultimately quasicrystals where translational symmetry is absent [6][7][8]. Crystal analogs in the  1. (a) The phase crystal has a periodic modulation of the superconducting phase χ (R) and a superflow p s (R) that forms a special vector field with a lattice of sources and sinks (filled circle), while the particle-conserving current j(R) forms a checkerboard pattern with opposite circulation flow. (b) This phase modulation is a result of four degenerate instability vectors {±q 0 , ±q 0 } with nonzero currents orthogonal to them [see Eq. (6)]. time dimension [9,10] have been observed recently [11,12]. Emergent multiparticle crystalline structures are predicted to appear in frustrated magnetic materials [13] and have been engineered in ultracold atoms interacting with light [14]. Superconducting states with periodically modulated amplitude (R) ∝ q cos(q · R) were first proposed to exist in ferromagnetic metals [15] and are currently investigated in a variety of systems ranging from cold Fermi gases with spin imbalance [16,17] to color superconductivity [18].
Several features make the phase crystal a distinctly different ground state from other nonuniform superconducting states. The amplitude-modulated state and its single-mode [19] counterpart (R) ∝ q e iq·R are both amplitude instabilities of the normal metal occurring at finite q and they do not carry currents. The phase crystal, on the other hand, is associated with a modification of the symmetry variable χ describing the degeneracy manifold of the superconducting state and can occur even when the order parameter amplitude | | is large, i.e., deep inside the superconducting state far from the normal to superconductor transition; the phase crystal does maintain nontrivial particle currents. Moreover, it is also different from the textures appearing in systems with multicomponent order parameters and a more complex degeneracy space, such as 3 He and liquid crystals [20][21][22]. In those systems the long-wavelength textures are a result of a competition between condensation and gradient terms involving different combinations of the order parameter components. The phase crystal is a result of a highly nonlocal superfluid response when sample surfaces, geometry, or other external influences impose a certain structure on the superfluid kernel itself. The patterns are formed on the much shorter coherence length scale ξ 0 =hv F /2π k B T c , where v F is the Fermi velocity, T c is the superconducting transition temperature, and k B is the Boltzmann constant (h = k B = 1 in the following). To describe this physics we ignore the amplitude gradient terms in the free energy and generalize the kinetic superflow energy in the limit of small p s as where we introduce a nonlocal superfluid density kernel K i j (R, R ) = K ji (R , R). Summation over repeating spatial indices is assumed. Higher-order gradient terms in F sf would determine the magnitude of spontaneous currents at temperatures below the transition temperature. Here we neglect those and focus on the instability analysis. 1 The energy change due to a small Galilean boost u, The physical χ and j are obtained by variational minimization of the free energy with respect to the phase. It gives the continuity equation −δF sf [∇χ ]/δχ (R) = ∇ · j(R) = 0.

II. PHASE INSTABILITY IN THE BULK
By using the nonlocal Ginzburg-Landau expression in Eq. (3) one can specify the general criteria when a nontrivial pattern of currents can emerge from the state with homogeneous phase χ 0 = 0. In a translationally invariant infinite system the superfluid free energy with kernelK (R − R ) has the following form in Fourier space: For the two-dimensional case, the kernel is a 2 × 2 Hermitian matrixK (q) =K † (q) with real eigenvalues κ 1,2 and corresponding eigenvectors e 1,2 . Their values depend on temperature and q. The instability at a particular wave vector q 0 can happen when q T 0K (q 0 )q 0 = κ 1 [e 1 · q 0 ] 2 + κ 2 [e 2 · q 0 ] 2 = 0. This equality can be satisfied if the eigenvalues have opposite signs and are tunable by temperature, or more generally by 1 We also drop corrections to the superflow due to the vector potential A(R) of the self-induced field ∇χ → ∇χ − 2π 0 A. These corrections result in energy terms that are smaller than the phase-gradient terms by a factor (ξ 0 /λ) 2 , which is small in type-II superconductors (see, e.g., Refs. [4,23]). some other parameter. To linear order in χ (q), the Fourier component of the current is j = j 0 iχ (q 0 ), where i = √ −1 and j 0 =K (q 0 )q 0 = e 1 κ 1 [e 1 · q 0 ] + e 2 κ 2 [e 2 · q 0 ]. For a nonzero current to appear at the q 0 = 0 transition, it must also satisfy the conservation law ∇ · j ∝ q 0 · j = 0. This implies an orthogonality constraint q 0 ⊥ j 0 , which is possible to fulfill if the eigenvectors e 1,2 are not collinear with q 0 [see Fig. 1(b)]. In this case we can write j 0 =x j 0x +ŷ j 0y , with j 0x / j 0y = −q 0y /q 0x . Since the phase χ (R) is real, the same conditions must be satisfied for −q 0 , which requires inversion symmetry. With two instability vectors q 0 and −q 0 we get an emerging phase χ (R) = C cos(q 0 · R) with stripes of current j(R) = Cj 0 sin(q 0 · R) running perpendicular to q 0 . Additional symmetries allow for other instability vectors. For example, reflection symmetry x → −x guarantees another pair of instability vectors q 0 and −q 0 , with q 0x = −q 0x . Diagonalization of the kernel at q 0 gives the same eigenvalues κ 1,2 as those at q 0 , while the eigenvectors e 1,2 are obtained from e 1,2 by flipping the x components, and the current amplitude is j 0 = e 1 κ 1 [e 1 · q 0 ] + e 2 κ 2 [e 2 · q 0 ]. In the four-harmonic state the phase and current are given by as plotted in Fig. 1(a). Higher-order terms O((∇χ ) 4 ) must be included to determine the energetics between two-and four-harmonic states. One notices that the loop currents in the phase crystal appear without phase winding and are not associated with topological defects. We conclude that realization of spontaneous periodic loop currents requires a superfluid density tensor with (i) spatial anisotropy, (ii) positive and negative eigenvalues that can be tuned by some parameter, and (iii) eigenvectors e 1,2 ∦ q 0 . Conditions (i) and (ii) can be satisfied simultaneously, for example, in an anisotropic-gap superconductor with an applied Zeeman field. Condition (iii) requires a mismatch between the momentum (related to q 0 ) and the velocity (determines the current response tensor) of the quasiparticles on the Fermi surface. To satisfy this last geometric condition, one would generally require a system with a spatial symmetry as low as possible. To formalize the analysis we can write a general Ginzburg-Landau expansion of the tensorK (q) in the superconducting state with orthorhombic symmetry C 2v (or D 2h in three dimensions). This symmetry is also required by condition (i) to have two eigenvectors of the kernel of different sign. The general form of the tensor is where the finite components are a 0 = K (0) xx = K (0) yy = b 0 , K (2) xxxx = a 2 , K (2) yyyy = b 2 , and K (2) xxyy = c 2 and all permutation of indices allowed. The configuration space of these five coefficients is large enough to allow for a set of instability wave vectors (q x , q y ) that do not lie along the high-symmetry directions and thus do not coincide with the direction of the current ( j x , j y ). Such a configuration would not be possible in a state with square symmetry that has only three independent coefficients a 0 = b 0 , a 2 = b 2 , and c 2 . The superfluid tensor will possess the C 2v symmetry in orthorhombic crystals, in nematically ordered systems, or in superconducting states with gap structure different along two principal axes, such as polar or planar states. The complete analysis of a crystallization transition with short-wavelength modulations is quite complex and has to include higher-order q terms. We leave this for future studies. We note that in typical weak crystallization theories the instability vectors are only given at the phenomenological level [7,8]. In the following we present the microscopic theory forK near pair-breaking surfaces and show how all these conditions are naturally satisfied and why a preferred ordering vector emerges.

III. SURFACE PHASE CRYSTAL
Using microscopic quasiclassical theory, we derive the general expression for the superfluid density kernel. The technical details of the calculation are given in Appendix A. We apply it first to the d-wave case and consider the s-wave case at the end of this section. The d-wave superconductor has an order parameter (R, p F ) = 0 (R)[2p xpy ] ≡ p , oriented as shown in Fig. 2(a). Thep = p F /|p F | is the unit vector pointing in the direction of momentum p F on the Fermi surface. The kernel between two points R and R in a semi-infinite system has two contributionsK (R, R ) = K 1 (R, R ) +K 2 (R, R ) that correspond to propagation of quasiparticles along the direct path or with a reflection at the surface. We set a uniform amplitude 0 (R) = 0 , which allows for analytic expressions (Appendix B). This assumption also demonstrates that the phase crystal is not caused by the suppression of the order parameter per se, but rather by the contribution from the symmetry-related surface Andreev bound states. The coordinate along a quasiparticle trajectory is denoted by s, with s = 0 at the reflection point. The kernel components are calculated in Appendix C, and for the direct path (p =p) they are where ε m = π T (2m + 1) are the Matsubara energies, κ u = 2 /v F , and = √ ε 2 m + 2p ; also s = s R − s R is the trajectory distance between the two points and s < = min(y, y )/|p y | is the trajectory coordinate of the point R or R closest to the surface. For the reflection path [p =p =p − 2ŷ(ŷ ·p)] where the overall minus sign is due to the fact that at the integration and observation points the order parameter has opposite signs p = − p . This reflection involving the sign (a) Microscopic model of the superfluid density tensor near a pair-breaking surface of a d xy superconductor. (b) Averaged local components (11) as a function of distance to the surface y and the modulation vector q x . The thinner dashed lines show direct path's contribution and the dotted lines the reflected path. The superfluid density far from the surface is determined by correlations between two points R and R through the direct path. This leads to positive superflow energy from diagonal components, favoring a uniform phase p s ∝ ∇χ 0 = 0. Near the surface the superflow energy is lowered by negative contributions of K xx and K yy coming from Andreev bound states, favoring the nonuniform phase crystal ∇χ = 0. change of the order parameter also leads to the zero-energy Andreev surface states [24]. The characteristic bound-state term, proportional to 2 p /ε 2 m , gives an overall 1/T temperature dependence of the kernel. The direct kernel in Eq. (8) may also show this 1/T dependence near the surface when the second term inside the square brackets dominates.
Pattern-forming instabilities are notorious for being technically challenging to analyze even at the level of linearized equations [25]. In what follows we work directly with the integral representation of the nonlocal physics. Since the unperturbed superconducting state is translationally invariant along the surface, we haveK (R, R ) =K (x 1 − x 2 , y 1 , 0, y 2 ), and we may write the superflow free energy in terms of Fourier components of the phase, χ (x, y) = C q x χ (y)e +iq x x , assuming the χ (y) profile to be real. We get where the prime denotes a derivative with respect to the y coordinate. The kernel is a complicated function of several variables K i j = K i j (q x , y 1 , y 2 ; T ). To describe its most important features we use a center coordinate representation y = (y 1 + y 2 )/2 and integrate over the relative coordinateȳ = y 1 − y 2 , This averaged response is shown in Fig. 2(b) as a function of distance from the surface y, where we also include the q x multiplication factors to directly relate the kernel to the free energy. For y L y ≈ 3ξ 0 -5ξ 0 , the response is dominated by the direct path. The off-diagonal components are zero and K xx and K yy are positive. Near the surface the diagonal components become negative, causing the instability, and large off-diagonal components appear. All components have the 1/T low-temperature dependence near the surface. The sign-changing nature of K i j and its T dependence lead to fulfillment of conditions (i) and (ii) for the phase crystal near the surface. Moreover, exponential decay of the bound states into the bulk creates an asymmetric environment at the surface with multiple q 0y components contributing to the instability. Condition (iii) is thereby also satisfied. We perform a variational analysis of Eq. (10) with an ansatz for the y dependence of the phase decaying into the bulk on the scale of y 0 , This choice is guided by considerations that there should be no currents deep in the sample, and we look for a state with no superflow in the y direction at the surface. The latter condition is not a strict requirement, since the physical condition of no current across the boundary j y (y = 0) = 0 is fulfilled automatically by the form of the total kernelK (R, R ). This guess gives a good semiquantitative result, but we note that to get the exact profile of χ (y) one has to perform a more sophisticated eigenvector analysis of the free energy (10). For each wave vector q x and temperature T we scan the variational parameter y 0 and find the minimum free energy. This minimum corresponds to the physical solution with currents satisfying ∇ · j = 0. The instability into the modulated-phase state with a nonzero C q x occurs at a temperature where the minimum of F sf crosses into negative values. The transition temperature T * (q x ) and the corresponding y 0 (q x ) are shown in Fig. 3(a) for the d-wave case. The highest transition temperature T * ∼ 0.3T c occurs at finite modulation q * x ≈ ξ −1 0 . By x → −x reflection symmetry there is degeneracy (q x , −q x ) that in the emerging state gives a real-value phase and superflow  In the vicinity of the optimal transition, the instability temperature behaves as Such dependence is a characteristic ansatz in theories of weak crystallization [7], where all the parameters are taken as phenomenological. We find T * ≈ 0.3T c , q * x ≈ 1.0/ξ 0 , and β ≈ 0.15T c ξ 2 0 . Here the appearance of a preferred finite phase modulation vector q * x is the result of an interplay between terms in the free energy (10) that in general have a different dependence on the y coordinates, T and q x . This physics can be crudely visualized by considering the superfluid freeenergy density, as shown in Figs. 3(b)-3(d). 2 The key element is the dependence of the phase decay length y 0 on q x [see Fig. 3(a), where we plot the inverse y −1 0 (q x )]. The superfluid response amplitudes grow with increasing q x . At the same time, the peaks in q 2 x K xx and q x K xy,yx move to smaller y [see Fig. 2(b)]. This requires a smaller y 0 to control the current components to satisfy ∇ · j = 0. Deviation of q x from its optimal value to smaller q x [compare Fig. 3(b) with Fig. 3(c)] leads to a longer extent away from the surface of the phase 2 The superfluid free-energy density cannot be uniquely defined in nonuniform, and especially nonlocal, systems. However, the two following definitions gave similar pictures: oscillations which increases the bulk energy cost from K xx and K yy . On the other hand, a deviation to larger q x gives a small y 0 , which results in a large cost due to off-diagonal K xy,yx components [compare Fig. 3(d) with Fig. 3(c)]. The instability for nonoptimal q x occurs at a lower temperature, where the K xx component becomes more negative near the surface by virtue of its 1/T dependence, which compensates for the energy increase in the other terms.
From this analysis we may conclude that the nonlocal multicomponent kernel leads to an intricate energy balance of the phase-gradient terms in the free energy. Because of the kernel structure, which fulfills the criteria (i)-(iii), a nontrivial phase crystallization occurs at a particular q * x ∼ 1/ξ 0 . To this broad class of phase instabilities belong several previously described surface states with paramagnetic surface currents caused by spectral displacement of Andreev states [26,27]. That work assumed translational invariance of the superflow and currents along the surface, which guaranteed particle conservation ∇ · j(R) = 0, but as a result required additional mechanisms of reducing superflow in the bulk. In semiinfinite systems one relies on the Meissner effect to screen the bulk superflow on the penetration depth length scale λ, which leads to T * ∼ (ξ 0 /λ)T c [23,28]. In slabs of width D < λ the bulk contribution is obviously limited, resulting in spontaneous superflow below T * ∼ (ξ 0 /D)T c [29]. In a similar fashion, we can interpret the phase crystal as self-screening of the loop currents over the surface region L y leading to T * ∼ (ξ 0 /L y )T c . A similar transition can appear in other anisotropic superconductors with reduced point group symmetry of the order parameter, such as the polar p wave [30], which hosts a flat band of zero-energy surface fermions, or some of the noncentrosymmetric three-dimensional superconductors that have finite areas of zero-energy states and prominent zero-energy peaks [31]. Interestingly, phase crystallization can happen in conventional s-wave superconductors, where orbital pair-breaking scattering is absent. In this case, magnetically active interfaces can provide the proper environment for the phase instability, for example, in superconductorferromagnetic structures. Such systems are being considered as important building blocks for spintronics applications, where nonlocality and quantum coherence will play important roles [32]. As described in Appendix C, a similar form of the superfluid density tensor appears for the ϑ = π spin-mixing angle. The phase diagram and the result of a self-consistent calculation are shown in Fig. 4.
The observable consequence of the spontaneous charge currents are magnetic fluxes near the surface. The associated reconstruction of the edge ground state is important from another perspective, since it can prevent realization of topological surface channels, as happens in topological insulators [34,35]. Moreover, softening of the surface superfluid density at some finite wave vector can result in special features of surface transport, even without a fully developed instability. This may be particularly relevant to transport in confined geometries.
Universal features of the pattern-formation phenomena in very different systems are manifested in the similarity of the phase diagram and the current patterns in Fig. 3 to those of the Rayleigh-Bénard convection instability, which is also a result of geometrical constraints and conservation laws. There the control parameter, instead of T , is the inverse Rayleigh ratio of buoyancy force to dissipative forces [36]. We note that the convection roll currents in that case is due to an instability in a nonequilibrium driven system, while the phase crystal is a second-order phase transition into a new ground state.

IV. CONCLUSION
We have described a superconducting state where the global U(1) phase spontaneously forms a modulation in space, breaking continuous translational invariance. The phase modulation results in a pattern of loop currents and breaking of time-reversal symmetry. We have identified the general criteria (i)-(iii) that have to be met in order to get a nonlocal superfluid density tensor that favors phase crystallization. Using microscopic theory, we showed that the circulating currents can appear at pair breaking surfaces of d-wave superconductors. In that case, quasiparticle reflections off the surface play a double role: (a) They lead to a flat band of zero-energy Andreev bound states controlling signs of the superfluid components and (b) they connect the y and x degrees of freedom at the level of the superfluid response resulting in preferred finite-q x modulation of the superflow. From previous numerical studies we know that this state remains stable in external magnetic fields [4] and survives significant reduction of spectral weight of bound states [37]. Thus, one should expect that similar phenomena will arise in other condensates with zero-energy surface states. To demonstrate this, we have stabilized the phase crystal in a conventional s-wave superconductor in contact with a magnetically active material, as can happen in hybrid superconductorferromagnet devices. One particularly interesting scenario, for the future, would be to generate this phase in a bulk system. The phase crystal presents an alternative vision of "supersolids" where phase-coherent states also spontaneously break translational symmetry, only in the amplitude of the order parameter [38][39][40][41]. More generally, our results indicate that nonlocal effects in broken-symmetry states, especially with multicomponent order parameters or competing orders, can lead to new states of matter. Such prospects are supported by early [42] and more recent [43] investigations of nonlocal physics in superconductors, as well as research into pattern formation due to long-range nonlocality in biological systems [44][45][46].

ACKNOWLEDGMENTS
The computations were performed on resources at Chalmers Centre for Computational Science and Engineering provided by the Swedish National Infrastructure for Computing. We thank the Swedish Research Council (VR) for financial support. P.H. acknowledges Chalmersska Forskningsfonden for travel support.

APPENDIX A: SUPERFLUID DENSITY NEAR A SURFACE
To find the superfluid response tensor we use a microscopic approach based on quasiclassical theory [47]. Our starting point is the Eilenberger equation for the quasiclassical propagatorĝ, In this equation a spatially varying phase χ of the order parameter = | |e iχ (R) was eliminated in favor of the superflow field p s = 1 2 ∇χ . This can always be done, if needed, by a gauge transformationĝ →ÛĝÛ † withÛ = e iτ 3 χ/2 . The superflow is a function of position p s = p s (R), and we consider a singlet mean-field order parameter = (R, p F ). The commutator-based Eilenberger equation is transformed into the Riccati-type equations for the coherence amplitudes [48] iv These amplitudes conveniently parametrize the quasiclassical propagator [49] and are functions of position, momentum, and energy, γ = γ (R, p F ; ε m ). The two coherence amplitudes are related by symmetry, which also applies to other tilde-related functions. For the singlet real order parameter˜ (R, p F ) ≡ * (R, −p F ) = (R, p F ). We look at the current response due to a small but arbitrary superflow field p s = p s (R), starting from a currentless background state 0 (R, p F ) and the corresponding coherence amplitudes γ 0 (R, p F ; ε m ). The following linear-response calculation is valid for any spatial profile of γ 0 (R, p F ; ε m ), and we specify in the end its particular form. The current at a point R near the surface is calculated from the correction to the diagonal propagator δg, with g = −iπ sgn(ε m ) 1−γγ 1+γγ as where N F is density of states at the Fermi level per spin projection and · · · p = dp/2π · · · denotes a cylindrical Fermi-surface average (Fig. 5). In terms of linearized coherence amplitudes γ = γ 0 + γ 1 , the propagator change due to small superflow is We first neglect the effect of the superflow on the amplitude of the order parameter, assuming that (R) = 0 (R) even in the current-carrying state, and linearize Eqs. (A2) to find transport equations for the function γ 1 /(1 + γ 0γ0 ), We get a similar equation for the tilde analog. The parameter determines the correlation length of the response. In a uniform state it reduces to κ = 2v −1 The solution of Eq. (A6) along a quasiclassical trajectory s is found, for positive ε m , by integration forward along the trajectory starting from zero value in the bulk γ 1 (s = −∞) = 0, where there is no superflow. We get To write the current at the observation point R we need to integrate over all trajectories arriving at point R. By introducing a correlation function connecting two points R 1 and R 2 via a 013104-7 dp F R dA F.S.

R'
ds=dR' FIG. 6. Connection between spatial integral and the trajectory-Fermi-surface integral. A volume element d 2 R in space can be written in cylindrical coordinates as d 2 R = dAds = |s R − s R |dp F ds, where |s R − s R | is the distance between points R and R along a trajectory and dp F is the angular integration over the Fermi surface.
one can combine the Fermi-surface average at the observation point and integration along trajectories into integration over all space R (see Fig. 6) and write the current response as Inserting (A8) into (A5) and using the definition (A9), the superfluid kernel is then given by where f 0 andf 0 are off-diagonal propagators in the unperturbed state. In terms of coherence amplitudes f 0 = −2iπ sgn(ε m ) γ 0 1+γ 0γ0 . This kernel connects the observation point R to the integration point R . For each pair of points there are two paths, one direct 1 and one involving reflection at the surface 2 , where we assumed mirrorlike reflection (see Fig. 5). The momentum directionp at the observation point is given by the trajectory direction R → R, and similarly for momentum at the integration pointp (Fig. 5). These directions are different for the direct and reflected paths.

APPENDIX B: COHERENCE AMPLITUDES AND PROPAGATORS WITH A STEPLIKE ORDER PARAMETER
Neglecting the suppression of the order parameter at the surface allows us to proceed further analytically. The bulk The coherence amplitudes can be found analytically if we ignore suppression of the order parameter at the interface. For each trajectory the order parameter sharply changes between i and f at s = 0. In this case, γ i on the incoming trajectory is a constant; then a boundary condition γ i → f gives an initial value that evolves to γ f on the outgoing part of the trajectory. For typical nonmagnetic specular scattering f = γ i . uniform coherence amplitude is Now consider (Fig. 7) a (straightened) trajectory that for s < 0 is in a region with the order parameter k = i and for s > 0 is in the region with k = f , e.g., for the most pair-breaking Far away from the interface, the coherence amplitudes have their uniform bulk values [we assume ε m > 0 and otherwise understand ε m = |ε m | and add sgn(ε m ) in front] For a sudden-step order parameter the amplitudes γ 0 andγ 0 (s) can be found analytically, integrating the Riccati equations (A2) in the forward or backward direction, correspondingly. Including the sudden jump of the amplitudes at the surface according to the boundary condition, we get and for tilde function integrating backward, The propagators on the trajectory are (e.g., for s > 0) and the off-diagonal component that enters the expression for the current response is where we write the functions in several different ways to cancel some terms later on. Notice the physical interpretation of the propagator form. For example, for f 0 (s > 0) we have the sameγ f in both terms since it is coming from s = +∞, but the γ amplitude can be either γ f far from the reflection point or f ← γ i close to reflection points and they give rise to the two different terms in f 0 . All other expressions for f functions follow the same pattern. The second term, which mixes f andγ f in the denominator, is the one that mainly determines bound-state effects. In both diagonal and off-diagonal items the continuum and the bound-state contribution are nicely separated.

APPENDIX C: CURRENT KERNEL WITHOUT THE ORDER PARAMETER SUPPRESSION
We use the results of Appendix B to calculate the current response kernel. First, we find κ, which determines the correlations extent in the current response, Here we consider an order parameter orientation such that the amplitudes on the incoming and reflected parts of the trajectory are the same, so κ u,i = κ u, f = κ u . The generalization for different amplitudes can be easily carried out retaining indices i, f , κ u;i, f , etc. This expression for κ (s) is quite general and easy to integrate along trajectories, as required for correlation functions C(R, R ) and C(R , R). In both these functions integration goes from the initial to the final point as determined by the momentum direction and is shown in Fig. 8.
For the case in Fig. 8(a) both s 1 and s 2 are on the same side of the interface and s 2 is farther away from the interface than s 1 ; thus we have, for s out, If we reverse the trajectory the signs of s change (so that s 1 and s 2 determine absolute distance to the surface); thus we have, for s in, For the Fig. 8(c) case we break the integral into two parts for s in and out: For any given points R and R we define two paths, direct and reflected, and each will have R → R and R → R contributions f (p, R)C(R, R ) f (p, R ) +f (−p, R)C(R , R) f (−p, R). Let us denote byk momentum away from the surface, and in this case we identify indices f =k and i =k. The trajectory integrating the γ function goes from s 1 = s < (point closest to the interface) to s 2 = s > (point farthest from interface). For the reverse trajectory we have f = −k and i = −k and integration happens from −s 2 to −s 1 .
The two terms give, after the mentioned cancellations, for the direct path, (1 − e −κ u s < ) +γk 1 + kγk e −κ u s < e −κ u |s > −s < | γk 1 + γkγk (1 − e −κ u s < ) + k 1 + kγk e −κ u s < For the reflected path this sum has a more compact form that directly reflects the bound states factors γ 0 1 + γ 0γ0 C(R ← R ) γ 0 1 + γ 0γ0 +γ Note that to generalize for the inequivalent gap size on in-out trajectories we need to use the appropriate κ u along given the directions, e.g., κ u |s > + s < | → κ u,k s > + κ u,k s < for the trajectoryk →k with reflection. These are completely general expressions for the one-component order parameters, where we neglect suppression of the order parameter amplitude near the surface and assume specular scattering. We apply the developed formalism and approximations to a d-wave superconductor with a maximally pair-breaking surface. In this case we have k = − k for all incident trajectories, and γ −k = γk = −γk, γ −k =γk =γk, and k = γk = −γk; two important combinations of the coherence amplitudes are The correlation coefficient (C1) along a trajectory s is where κ u = 2 /v F and = √ ε 2 m + 2 k . The distance along a trajectory, measured from the surface, is s = y/k y . One uses these relations for coherence amplitudes in combination with (C6) and (C7) to find the kernel (A11) components, as given in the main text, for the direct path (8) and the reflection path (9), correspondingly.
Similar expressions for the superfluid density are valid for an s-wave superconductor with scattering at a specular magnetically active surface. We use the boundary conditions for coherence amplitudes [50] k iσ 2 = Mγ k iσ 2M , with M = e iϑm·σ /2 andM = M * . Magnetic spin mixing leads to the bound states ε b = ± cos(ϑ/2), which result in zeroenergy states for ϑ = π and the boundary condition for coherence amplitudes k = −γk.