Phase Crystals

TL;DR

This paper introduces the concept of a phase crystal, a novel non-uniform superconducting ground state characterized by a spatially periodic phase modulation, which breaks translational and time-reversal symmetries, expanding understanding of superconducting phases.

Contribution

It proposes and analyzes the phase crystal state, deriving conditions for its realization and demonstrating its relevance to surface states and previous numerical observations in unconventional superconductors.

Findings

Identification of phase crystallization as a new ordered state in superconductors.

Analytic expression for superfluid density tensor in non-uniform environments.

Prediction of phase crystal phenomena in superconductor-ferromagnetic structures.

Abstract

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.

Equations

$$F_{\mbox{\tiny sf}}=\frac{1}{2}\int d{\bf R}\;k|\Delta|^{2}\,|\mbox{\boldmath$\nabla$}\chi({\bf R})|^{2}\,,$$

$$\chi (\mathbf{R})=C_{\mathbf{q}}{A}_{\mathbf{q}}({\mathbf{R}}_{\perp })\cos (\mathbf{q}\cdot \mathbf{R}),$$

$$\displaystyle\begin{split}F_{\mbox{\tiny sf}}[\mbox{\boldmath$\nabla$}\chi]=\frac{1}{2}\iint d{\bf R}d{\bf R}^{\prime}\;\nabla_{i}\chi({\bf R})K_{ij}({\bf R},{\bf R}^{\prime})\nabla_{j}\chi({\bf R}^{\prime})\,,\end{split}$$

$$\begin{array}{r}{j}_i(\mathbf{R})=\frac{\delta F_{\text{mbox}sf}[{\mathbf{v}}_s]}{\delta p_{s,i}(\mathbf{R})}=\int d{\mathbf{R}}^{\prime }{K}_{ij}(\mathbf{R},{\mathbf{R}}^{\prime }){\nabla }_j\chi ({\mathbf{R}}^{\prime }).\end{array}$$

$$\begin{array}{r}{F}_{\text{mbox}sf}=\frac{1}{2}\int \frac{d^{2}q}{(2\pi {)}^2}\chi (-\mathbf{q})\left[{\mathbf{q}}^T\hat{K}(\mathbf{q})\mathbf{q}\right]\chi (\mathbf{q}).\end{array}$$

$$\begin{array}{rl}\chi (\mathbf{R})& =\cos ({\mathbf{q}}_0\cdot \mathbf{R})+\cos ({\underline{\mathbf{q}}}_0\cdot \mathbf{R})\propto \cos (q_{0x}x)\cos (q_{0y}y),\\ \mathbf{j}(\mathbf{R})& ={\mathbf{j}}_0\sin ({\mathbf{q}}_0\cdot \mathbf{R})+{\underline{\mathbf{j}}}_0\sin ({\underline{\mathbf{q}}}_0\cdot \mathbf{R})\\ & \propto \left(\begin{array}{c}{j}_{x0}\sin (q_{0x}x)\cos (q_{0y}y)\\ j_{0y}\cos (q_{0x}x)\sin (q_{0y}y)\end{array}\right),\end{array}$$

$$\begin{array}{rl}{K}_{ij}(q_x,q_y)& =K_{ij}^{(0)}+K_{ijlm}^{(2)}{q}_l{q}_m+\cdots =\\ & =\left(\begin{array}{cc}{a}_0+a_2{q}_x^2+c_2{q}_y^2& 2c_2{q}_x{q}_y\\ 2c_2{q}_x{q}_y& b_0+b_2{q}_y^2+c_2{q}_x^2\end{array}\right)\end{array}$$

$$\displaystyle\begin{split}K^{\mbox{\tiny1}⃝}_{ij}({\bf R},{\bf R}^{\prime})=&\left[{\hat{p}}_{i}{\hat{p}}_{j}\right]\,v_{\mbox{\tiny F}}^{2}N_{\mathrm{F}}\;4\pi T\sum_{\varepsilon_{m}>0}\frac{\Delta^{2}_{\hat{p}}}{\Omega^{2}}\frac{2}{v_{\mbox{\tiny F}}}\,\frac{e^{\displaystyle-\kappa_{u}|{\mbox{\footnotesize$\Delta$}}s|}}{2\pi|{\mbox{\footnotesize$\Delta$}}s|}\\ &\times\left[\left(1-e^{-\kappa_{u}|s_{<}|}\right)^{2}-\frac{\Omega^{2}}{\varepsilon_{m}^{2}}e^{-2\kappa_{u}|s_{<}|}\right],\end{split}$$

$$\displaystyle\begin{split}K^{\mbox{\tiny2}⃝}_{ij}({\bf R},{\bf R}^{\prime})=-\left[{\hat{p}}_{i}\underline{{\hat{p}}}_{j}\right]\,v_{\mbox{\tiny F}}^{2}N_{\mathrm{F}}\;4\pi T\sum_{\varepsilon_{m}>0}\frac{\Delta^{2}_{\hat{p}}}{\varepsilon_{m}^{2}}\frac{2}{v_{\mbox{\tiny F}}}\,\frac{e^{\displaystyle-\kappa_{u}|{\mbox{\footnotesize$\Delta$}}s|}}{2\pi|{\mbox{\footnotesize$\Delta$}}s|},\end{split}$$

$$\displaystyle\begin{split}F_{\mbox{\tiny sf}}=&\frac{1}{2}\int\frac{dq_{x}}{(2\pi)}|C_{q_{x}}|^{2}\int\limits_{0}^{\infty}dy_{1}\int\limits_{0}^{\infty}dy_{2}\;\times\\ &\Big{[}q_{x}^{2}K_{xx}\;\chi(y_{1})\chi(y_{2})+K_{yy}\;\chi^{\prime}(y_{1})\,\chi^{\prime}(y_{2})\\ &-iq_{x}K_{xy}\;\chi(y_{1})\,\chi^{\prime}(y_{2})+iq_{x}K_{yx}\;\chi^{\prime}(y_{1})\,\chi(y_{2})\Big{]},\end{split}$$

$${\overline{K}}_{ij}(q_x,y;T)=\underset{-2y}{\overset{2y}{\int }}d\bar{y}{K}_{ij}\left(q_x,y+\frac{1}{2}\bar{y},y-\frac{1}{2}\bar{y};T\right).$$

$$\chi (y)=\left(1+\frac{y}{y_0}\right)e^{-\frac{y}{y_0}},{\chi }^{\prime }(y)=-\frac{y}{y_0^2}{e}^{-\frac{y}{y_0}}.$$

$$\begin{array}{rl}& \chi (x,y)\propto -\left(1+\frac{y}{y_0}\right)e^{-y/y_0}\cos q_{x}x,\\ & {\mathbf{p}}_s(x,y)\propto \left[q_x\left(1+\frac{y}{y_0}\right)\sin q_{x}x,\frac{y}{y_0^2}\cos q_{x}x\right]e^{-y/y_0},\end{array}$$

$$T^{\ast }(q_x)=T^{\ast }-\beta (q_x-q_x^{\ast }{)}^2.$$

$$[(i\varepsilon_{m}-{\bf v}_{f}\cdot{\bf p}_{s})\hat{\tau}_{3}-\hat{\Delta}({\bf R},{\bf p}_{\mbox{\tiny F}})\,,\,\hat{g}]+i{\bf v}_{f}\cdot\mbox{\boldmath$\nabla$}\hat{g}=0$$

$$\displaystyle\begin{split}&i{\bf v}_{\mbox{\tiny F}}\cdot\mbox{\boldmath$\nabla$}\gamma+2[i\varepsilon_{m}-{\bf v}_{\mbox{\tiny F}}\cdot{\bf p}_{s}]\gamma+\gamma\tilde{\Delta}\gamma+\Delta=0,\\ &i{\bf v}_{\mbox{\tiny F}}\cdot\mbox{\boldmath$\nabla$}\tilde{\gamma}-2[i\varepsilon_{m}-{\bf v}_{\mbox{\tiny F}}\cdot{\bf p}_{s}]\tilde{\gamma}+\tilde{\gamma}\Delta\tilde{\gamma}+\tilde{\Delta}=0.\end{split}$$

$$\tilde{\gamma }(\mathbf{R},{\mathbf{p}}_{\text{mbox}F};{\epsilon }_m)=\gamma (\mathbf{R},-{\mathbf{p}}_{\text{mbox}F};{\epsilon }_m{)}^{\ast },$$

$$\mathbf{j}(\mathbf{R})=2T\sum _{{\epsilon }_m>0}2N_{\mathrm{F}}\text{mbox}Re{⟨v_{\text{mbox}F}\hat{p}\delta g(\mathbf{R},{\mathbf{p}}_{\text{mbox}F};{\epsilon }_m)⟩}_{\hat{p}},$$

$$\delta g(\mathbf{R},{\mathbf{p}}_{\text{mbox}F};{\epsilon }_m)=2i\pi \text{mbox}sgn({\epsilon }_m)\frac{{\gamma }_1{\tilde{\gamma }}_0+{\gamma }_0{\tilde{\gamma }}_1}{(1+{\gamma }_0{\tilde{\gamma }}_0{)}^2}.$$

$$\displaystyle\begin{split}{\hat{p}}\cdot\mbox{\boldmath$\nabla$}\frac{\gamma_{1}}{1+\gamma_{0}\tilde{\gamma}_{0}}+\kappa\frac{\gamma_{1}}{1+\gamma_{0}\tilde{\gamma}_{0}}=-2i\;{\hat{p}}\cdot{\bf p}_{s}\;\frac{\gamma_{0}}{1+\gamma_{0}\tilde{\gamma}_{0}}.\end{split}$$

$$\kappa (\mathbf{R},\hat{p};{\epsilon }_m)\equiv \frac{2}{v_{\text{mbox}F}}\left[{\epsilon }_m+\frac{{\gamma }_0{\tilde{\Delta }}_0-{\tilde{\gamma }}_0{\Delta }_0}{2i}\right]=\tilde{\kappa },$$

$$\begin{array}{r}\frac{{\gamma }_1}{1+{\gamma }_0{\tilde{\gamma }}_0}(\mathbf{R},\hat{p};{\epsilon }_m)=-2i\underset{-\infty }{\overset{s_{\mathbf{R}}}{\int }}ds\exp \left(-{\int }_s^{s_{\mathbf{R}}}\kappa (\rho )d\rho \right)\\ \times \hat{p}(s)\cdot {\mathbf{p}}_s({\mathbf{R}}^{\prime }(s))\frac{{\gamma }_0}{1+{\gamma }_0{\tilde{\gamma }}_0}(s).\end{array}$$

$$C({\mathbf{R}}_2,{\mathbf{R}}_1)=\frac{1}{2\pi |{\mathbf{R}}_2-{\mathbf{R}}_1|}\frac{2{\epsilon }_m}{v_{\text{mbox}F}}\exp \left(-{\int }_{{\mathbf{R}}_1}^{{\mathbf{R}}_2}\kappa (\rho ,\hat{\rho })d\rho \right),$$

$$j_i(\mathbf{R})=\int d^2{R}^{\prime }{K}_{ij}(\mathbf{R},{\mathbf{R}}^{\prime })p_{s,j}({\mathbf{R}}^{\prime }).$$

$$\begin{array}{r}{K}_{ij}(\mathbf{R},{\mathbf{R}}^{\prime })=v_{\text{mbox}F}^2{N}_{\mathrm{F}}8\pi T\sum _{{\epsilon }_m>0}\sum _{\text{mbox}1⃝,\text{mbox}2⃝}\frac{1}{4{\pi }^2{\epsilon }_m}\times \\ \times \text{mbox}Re[{\hat{p}}_i{\tilde{f}}_0(\mathbf{R},\hat{p})C(\mathbf{R},{\mathbf{R}}^{\prime })f_0({\mathbf{R}}^{\prime },{\hat{p}}^{\prime }){\hat{p}}_j^{\prime }+\\ +{\hat{p}}_j^{\prime }{\tilde{f}}_0({\mathbf{R}}^{\prime },-{\hat{p}}^{\prime })C({\mathbf{R}}^{\prime },\mathbf{R})f_0(\mathbf{R},-\hat{p}){\hat{p}}_i],\end{array}$$

$$\begin{array}{r}\gamma =i\frac{\Delta }{|{\epsilon }_m|+\sqrt{{\Delta }^2+{\epsilon }_m^2}}\text{mbox}sgn({\epsilon }_m),\\ {\epsilon }_m-i\gamma \Delta =\text{mbox}sgn({\epsilon }_m)\sqrt{{\epsilon }_m^2+{\Delta }^2}\end{array}$$

$$\begin{array}{r}{\Omega }_i=\sqrt{{\Delta }_i^2+{\epsilon }_m^2},{\kappa }_{u,i}=\frac{2}{v_f}\sqrt{{\Delta }_i^2+{\epsilon }_m^2},\\ {\Omega }_f=\sqrt{{\Delta }_f^2+{\epsilon }_m^2},{\kappa }_{u,f}=\frac{2}{v_f}\sqrt{{\Delta }_f^2+{\epsilon }_m^2}.\end{array}$$

$$\begin{array}{r}{\gamma }_i=i\frac{{\Delta }_i}{{\epsilon }_m+{\Omega }_i},{\gamma }_f=i\frac{{\Delta }_f}{{\epsilon }_m+{\Omega }_f},\\ {\tilde{\gamma }}_i=-i\frac{{\tilde{\Delta }}_i}{{\epsilon }_m+{\Omega }_i},{\tilde{\gamma }}_f=-i\frac{{\tilde{\Delta }}_f}{{\epsilon }_m+{\Omega }_f}\end{array}$$

$$\begin{array}{r}{\gamma }_0(s<0)={\gamma }_i\to {\gamma }_0(s=+0)={\Gamma }_f⟶\\ {\gamma }_0(s>0)={\gamma }_f+\frac{(1+{\gamma }_f{\tilde{\gamma }}_f)({\Gamma }_f-{\gamma }_f)e^{-{\kappa }_{u,f}s}}{1+{\gamma }_f{\tilde{\gamma }}_f+({\Gamma }_f-{\gamma }_f){\tilde{\gamma }}_f(1-e^{-{\kappa }_{u,f}s})}\end{array}$$

$$\begin{array}{r}{\tilde{\gamma }}_0(s>0)={\tilde{\gamma }}_f\to {\tilde{\gamma }}_0(s=-0)={\tilde{\Gamma }}_i⟶\\ {\tilde{\gamma }}_0(s<0)={\tilde{\gamma }}_i+\frac{(1+{\gamma }_i{\tilde{\gamma }}_i)({\tilde{\Gamma }}_i-{\tilde{\gamma }}_i)e^{{\kappa }_{u,i}s}}{1+{\gamma }_i{\tilde{\gamma }}_i+{\gamma }_i({\tilde{\Gamma }}_i-{\tilde{\gamma }}_i)(1-e^{{\kappa }_{u,i}s})}\end{array}$$

Full text

I Introduction

The defining characteristic of superfluidity and superconductivity is spontaneous symmetry breaking of the global $U(1)$ phase $\chi$, associated with the order parameter $\Delta=|\Delta|\exp(i\chi)$. 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 Tinkham (1985). Within the BCS paradigm Bardeen et al. (1957), a uniform fixed value of the phase is directly tied to the finite amplitude $|\Delta|$ of the macroscopic Cooper-pair wavefunction. If the phase is non-uniform, by Galilean invariance it results in superflow with superfluid velocity and momentum m{\bf v}_{s}={\bf p}_{s}({\bf R})=(\hbar/2)\mbox{\boldmath\nabla}\chi({\bf R}), where $m$ is the electron mass and $\hbar$ is the reduced Planck constant. Such phase variations and the associated condensate currents cost gradient energy

[TABLE]

where the gradient energy coefficient $k>0$ should be computed from microscopic theory. A physical picture emerges where the phase is rigid, 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_{\mbox{\tiny 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, that we denote a phase crystalline state, a periodic pattern with wavevector ${\bf q}$ is formed

[TABLE]

where $A_{\bf q}({\bf R}_{\perp})$ is a function of coordinates orthogonal to ${\bf q}$. The additional order parameter in the phase crystal is the finite Fourier amplitude $C_{\bf q}$. The superconducting ground state with spatially oscillating phase also breaks time-reversal symmetry and sustains a non-trivial periodic superflow pattern and circulating currents ${\bf j}({\bf R})$, as illustrated in Fig. 1a. Similar current patterns have been found in numerical work on mesoscopic grains of $d$-wave superconductors Håkansson et al. (2015), and the unusual superflow field ${\bf p}_{s}({\bf R})$ was recently analyzedHolmvall et al. (2018). 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 Powell (2010) and ultimately quasicrystals where translational symmetry is absent Senechal (1995); *symm_crystal; Kats et al. (1993); Martin et al. (2016). Crystal analogues in the time dimension Wilczek (2012); Yao et al. (2017) have been recently observed Zhang et al. (2017); Choi et al. (2017). Emergent multi-particle crystalline structures are predicted to appear in frustrated magnetic materials,Kamiya and Batista (2014) and have been engineered in ultracold atoms interacting with light.Ostermann et al. (2016) Superconducting states with periodically modulated amplitude $\Delta({\bf R})\propto\Delta_{\bf q}\cos({\bf q}\cdot{\bf R})$ were first proposed to exist in ferromagnetic metals Larkin and Ovchinnikov (1964); *Larkin:1965wj, and are currently investigated in a variety of systems ranging from cold Fermi-gases with spin imbalance Kinnunen et al. (2018); Dutta and Mueller (2017) to color superconductivity Casalbuoni and Nardulli (2004).

Several features make the phase crystal a distinctly different ground state from other non-uniform superconducting states. The amplitude-modulated state and its single-mode Fulde and Ferrell (1964) counterpart $\Delta({\bf R})\propto\Delta_{\bf q}e^{i{\bf q}\cdot{\bf R}}$, are both amplitude instabilities of the normal metal occurring at finite ${\bf q}$, and they do not carry currents. The phase crystal, on the other hand, is associated with a modification of the symmetry variable $\chi$ describing the degeneracy manifold of the superconducting state, and can occur even when the order parameter amplitude $|\Delta|$ is large, i.e. deep inside the superconducting state far from the normal to superconductor transition; the phase crystal does maintain non-trivial particle currents. Moreover, it is also different from the textures appearing in systems with multi-component order parameters and a more complex degeneracy space, such as 3He and liquid crystalsD.Vollhardt and P.Wölfle (1990); Chaikin and Lubensky (1995); de Gennes and Prost (1995). 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 non-local 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 $\xi_{0}=\hbar v_{\mbox{\tiny F}}/2\pi k_{\mathrm{B}}T_{\mathrm{c}}$, where $v_{\mbox{\tiny F}}$ is the Fermi velocity, $T_{\mathrm{c}}$ is the superconducting transition temperature and $k_{\mathrm{B}}$ is the Boltzmann constant ($\hbar=k_{\mathrm{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 ${\bf p}_{s}$ as

[TABLE]

where we introduce a non-local superfluid density kernel $K_{ij}({\bf R},{\bf R}^{\prime})=K_{ji}({\bf R}^{\prime},{\bf R})$. Summation over repeating spatial indices is assumed. Higher order gradient terms in $F_{\mbox{\tiny sf}}$ would determine the magnitude of spontaneous currents at temperatures below the transition temperature. Here we neglect those and focus on the instability analysis. 111 We also drop corrections to the superflow due to the vector potential ${\bf A}({\bf R})$ of the self-induced field \mbox{\boldmath\nabla}\chi\to\mbox{\boldmath\nabla}\chi-\frac{2\pi}{\Phi_{0}}{\bf A}. These corrections result in energy terms that are smaller than the phase-gradient terms by factor $(\xi_{0}/\lambda)^{2}$, which is small in type-II superconductors. See e.g. Refs. Holmvall et al., 2018; Barash et al., 2000

The energy change due to a small Galilean boost ${\bf u}$, $F_{\mbox{\tiny sf}}[{\bf v}_{s}-{\bf u}]=F_{\mbox{\tiny sf}}[{\bf v}_{s}]-m{\bf j}\cdot{\bf u}$, defines the particle current

[TABLE]

The physical $\chi$ and ${\bf j}$ are obtained by variational minimization of the free energy with respect to the phase. It gives the continuity equation, -{\delta F_{\mbox{\tiny sf}}[\mbox{\boldmath\nabla}\chi]}/{\delta\chi({\bf R})}=\mbox{\boldmath\nabla}\cdot{\bf j}({\bf R})=0.

II phase instability in the bulk

By using the non-local Ginzburg-Landau expression in Eq. (3) one can specify the general criteria when a non-trivial pattern of currents can emerge from the state with homogeneous phase $\chi_{0}=0$. In a translationally-invariant infinite system the superfluid free energy with kernel $\hat{K}({\bf R}-{\bf R}^{\prime})$ has the following form in Fourier space

[TABLE]

For the two-dimensional case, the kernel is a two-by-two Hermitian matrix $\hat{K}({\bf q})=\hat{K}^{\dagger}({\bf q})$ with real eigenvalues $\kappa_{1,2}$ and corresponding eigenvectors ${\bf e}_{1,2}$. Their values depend on temperature and ${\bf q}$. The instability at a particular wavevector ${\bf q}_{0}$ can happen when ${\bf q}_{0}^{T}\hat{K}({\bf q}_{0}){\bf q}_{0}=\kappa_{1}[{\bf e}_{1}\cdot{\bf q}_{0}]^{2}+\kappa_{2}[{\bf e}_{2}\cdot{\bf q}_{0}]^{2}=0$. This equality can be satisfied if the eigenvalues have opposite signs and are tunable by temperature, or more generally by some other parameter. To linear order in $\chi({\bf q})$, the Fourier component of the current is ${\bf j}={\bf j}_{0}\,i\,\chi({\bf q}_{0})$, where $i=\sqrt{-1}$ and ${\bf j}_{0}=\hat{K}({\bf q}_{0}){\bf q}_{0}={\bf e}_{1}\kappa_{1}[{\bf e}_{1}\cdot{\bf q}_{0}]+{\bf e}_{2}\kappa_{2}[{\bf e}_{2}\cdot{\bf q}_{0}]$. For a non-zero current to appear at the ${\bf q}_{0}\neq 0$ transition, it must also satisfy the conservation law \mbox{\boldmath\nabla}\cdot{\bf j}\propto{\bf q}_{0}\cdot{\bf j}=0. This implies an orthogonality constraint ${\bf q}_{0}\perp{\bf j}_{0}$, which is possible to fulfill if the eigenvectors ${\bf e}_{1,2}$ are not collinear with ${\bf q}_{0}$, see Fig. 1b. In this case we can write ${\bf j}_{0}=\hat{\bf x}j_{0x}+\hat{\bf y}j_{0y}$ with $j_{0x}/j_{0y}=-q_{0y}/q_{0x}$. Since the phase $\chi({\bf R})$ is real, the same conditions must be satisfied for $-{\bf q}_{0}$, which requires inversion symmetry. With two instability vectors ${\bf q}_{0}$ and $-{\bf q}_{0}$ we get an emerging phase $\chi({\bf R})=C\cos({\bf q}_{0}\cdot{\bf R})$ with stripes of current ${\bf j}({\bf R})=C{\bf j}_{0}\sin({\bf q}_{0}\cdot{\bf R})$ running perpendicular to ${\bf q}_{0}$. Additional symmetries allow for other instability vectors. For example, reflection symmetry $x\to-x$ guarantees another pair of instability vectors, $\underline{{\bf q}}_{0}$ and $-\underline{{\bf q}}_{0}$, with $\underline{q}_{0x}=-q_{0x}$. Diagonalization of the kernel at $\underline{{\bf q}}_{0}$ gives the same eigenvalues $\kappa_{1,2}$ as those at ${\bf q}_{0}$, while the eigenvectors $\underline{{\bf e}}_{1,2}$ are obtained from ${\bf e}_{1,2}$ by flipping the $x$-components, and the current amplitude is $\underline{{\bf j}}_{0}=\underline{{\bf e}}_{1}\kappa_{1}[\underline{{\bf e}}_{1}\cdot\underline{{\bf q}}_{0}]+\underline{{\bf e}}_{2}\kappa_{2}[\underline{{\bf e}}_{2}\cdot\underline{{\bf q}}_{0}]$. In the four-harmonics state the phase and current are given by

[TABLE]

as plotted in Fig.1a. Higher order terms {\cal O}[(\mbox{\boldmath\nabla}\chi)^{4}] must be included to determine the energetics between two- and four-harmonics 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

1. (i)

spatial anisotropy, 2. (ii)

positive and negative eigenvalues that can be tuned by some parameter, 3. (iii)

eigenvectors ${\bf e}_{1,2}\nparallel{\bf 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 symmetry of the Fermi surface and the quasiparticle excitations in momentum space, and the symmetry of the current response tensor. To satisfy this last geometric condition, one would generally require a system with as lower spatial symmetry as possible. To formalize the analysis we can write a general Ginzburg-Landau expansion of the tensor $\hat{K}({\bf q})$ in the superconducting state with orthorhombic symmetry $C_{2v}$. 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

[TABLE]

where finite components are $a_{0}=K_{xx}^{(0)}\neq K_{yy}^{(0)}=b_{0}$, $K_{xxxx}^{(2)}=a_{2}$, $K_{yyyy}^{(2)}=b_{2}$, $K_{xxyy}^{(2)}=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 wavevector $(q_{x},q_{y})$ that do not lie along the high symmetry directions, and thus do not coincide with direction of the current $(j_{x},j_{y})$. Such 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 a short-wavelength modulations is quite complex, and has to include higher order ${\bf q}$-terms. We leave this for future studies. We note that in typical weak crystallization theories the instability vectors are only given at phenomenological level.Kats et al. (1993); Martin et al. (2016) In the following we write down the microscopic theory for $\hat{K}$ near pairbreaking 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 moved to 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 $\Delta({\bf R},{\bf p}_{\mbox{\tiny F}})=\Delta_{0}({\bf R})\,[2{\hat{p}}_{x}{\hat{p}}_{y}]\equiv\Delta_{\hat{p}}$, oriented as shown in Fig. 2a. The $\hat{p}={\bf p}_{\mbox{\tiny F}}/|{\bf p}_{\mbox{\tiny F}}|$ is the unit vector pointing in the direction of momentum ${\bf p}_{\mbox{\tiny F}}$ on the Fermi surface. The kernel between two points ${\bf R}$ and ${\bf R}^{\prime}$ in a semi-infinite system has two contributions, $\hat{K}({\bf R},{\bf R}^{\prime})=\hat{K}^{\mbox{\tiny1}⃝}({\bf R},{\bf R}^{\prime})+\hat{K}^{\mbox{\tiny2}⃝}({\bf R},{\bf R}^{\prime})$, that correspond to propagation of quasiparticles along the direct path or with a reflection at the surface. We set a uniform amplitude $\Delta_{0}({\bf R})=\Delta_{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 (${\hat{p}}^{\prime}={\hat{p}}$) they are

[TABLE]

where $\varepsilon_{m}=\pi T(2m+1)$ are the Matsubara energies, $\kappa_{u}=2\Omega/v_{\mbox{\tiny F}}$ and $\Omega=\sqrt{\varepsilon_{m}^{2}+\Delta_{\hat{p}}^{2}}$; also {\mbox{\footnotesize\Delta}}s=s_{\bf R}-s_{{\bf R}^{\prime}} is the trajectory distance between the two points, and $s_{<}=\mbox{min}(y,y^{\prime})/|{\hat{p}}_{y}|$ is the trajectory coordinate of the point, ${\bf R}$ or ${\bf R}^{\prime}$, closest to the surface. For the reflection path (${\hat{p}}^{\prime}=\underline{{\hat{p}}}={\hat{p}}-2\hat{y}(\hat{y}\cdot{\hat{p}})$)

[TABLE]

where the overall minus sign is due to the fact that at the integration and observation points the order parameter has opposite signs $\Delta_{\underline{{\hat{p}}}}=-\Delta_{\hat{p}}$. This reflection involving the sign-change of the order parameter also leads to the zero-energy Andreev surface states.Hu (1994) The characteristic bound states term, proportional to $\Delta_{\hat{p}}^{2}/\varepsilon_{m}^{2}$, 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 linearised equations Pesch and Kramer (1996). In what follows we work directly with the integral representation of the non-local physics. Since the unperturbed superconducting state is translationally invariant along the surface, we have $\hat{K}({\bf R},{\bf R}^{\prime})=\hat{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, $\chi(x,y)=C_{q_{x}}\chi(y)e^{+iq_{x}x}$, assuming the $\chi(y)$-profile to be real. We get

[TABLE]

where the prime denotes a derivative with respect to the $y$-coordinate. The kernel is a complicated function of several variables $K_{ij}=K_{ij}(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 $\bar{y}=y_{1}-y_{2}$,

[TABLE]

This averaged response is shown in Fig. 2b as 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\gtrsim L_{y}\approx 3\div 5\xi_{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_{ij}$, and its $T$-dependence, lead to fulfilment 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}$,

[TABLE]

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 kernel $\hat{K}({\bf R},{\bf R}^{\prime})$. This guess gives a good semi-quantitative result, but we note that to get the exact profile of $\chi(y)$ one has to perform a more sophisticated eigenvector analysis of the free energy Eq. (10). For each wave vector $q_{x}$ and temperature $T$ we scan the variational parameter $y_{0}$ and find the minimum of the free energy. This minimum corresponds to the physical solution with currents satisfying \mbox{\boldmath\nabla}\cdot{\bf j}=0. The instability into the modulated-phase state with a non-zero $C_{q_{x}}$ occurs at a temperature where the minimum of $F_{\mbox{\tiny sf}}$ crosses into negative values. The transition temperature $T^{*}(q_{x})$ and the corresponding $y_{0}(q_{x})$ are shown in Fig. 3a, for the $d$-wave case. The highest transition temperature $T^{*}\sim 0.3T_{c}$ occurs at finite modulation $q_{x}^{*}\approx\xi_{0}^{-1}$. By $x\to-x$ reflection symmetry there is degeneracy $(q_{x},-q_{x})$ that in the emerging state gives a real-valued phase and superflow

[TABLE]

with the superflow exhibiting critical points ${\bf p}_{s}=0$ at the surface, as marked in Figs. 3b-d by filled orange circles.

In the vicinity of the optimal transition, the instability temperature behaves as

[TABLE]

Such dependence is a characteristic ansatz in theories of weak crystallizationKats et al. (1993), where all the parameters are taken as phenomenological. We find $T^{*}\approx 0.3T_{c}$, $q_{x}^{*}\approx 1.0/\xi_{0}$ and $\beta\approx 0.15T_{c}\xi_{0}^{2}$. Here the appearance of a preferred finite phase modulation vector $q_{x}^{*}$ is the result of an interplay between terms in the free energy Eq. (10) that in general have different dependence on the $y$-coordinates, $T$ and $q_{x}$. This physics can be crudely visualized by considering the superfluid free energy density, as shown in Fig. 3b-d. 222 The superfluid free energy density cannot be uniquely defined in non-uniform, and especially non-local, systems. However, the two following definitions gave similar pictures: $f_{1}({\bf R})=\int d{\bf r}\;{\bf p}_{s}({\bf R}_{+})^{T}\,\hat{K}({\bf R}_{+},{\bf R}_{-})\,{\bf p}_{s}({\bf R}_{-})$ with ${\bf R}_{\pm}={\bf R}\pm{\bf r}/2$, and $f_{2}({\bf R})={\bf p}_{s}({\bf R})\cdot{\bf j}({\bf R})={\bf p}_{s}({\bf R})^{T}\int d{\bf R}^{\prime}\,\hat{K}({\bf R},{\bf R}^{\prime})\,{\bf p}_{s}({\bf R}^{\prime})$.

The key element is the dependence of the phase decay length $y_{0}$ on $q_{x}$, see Fig. 3a where we plot the inverse $y_{0}^{-1}(q_{x})$. The superfluid response amplitudes grow with increasing $q_{x}$. At the same time, the peaks in $q_{x}^{2}K_{xx}$ and $q_{x}K_{xy,yx}$ move to smaller $y$, see Fig. 2b. This requires a smaller $y_{0}$ to control the current components to satisfy \mbox{\boldmath\nabla}\cdot{\bf j}=0. Deviation of $q_{x}$ from its optimal value to smaller $q_{x}$, compare Fig. 3b with Fig. 3c, leads to a longer extent away from the surface of the phase 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. 3d with Fig. 3c. The instability for non-optimal $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 non-local multi-component kernel leads to an intricate energy balance of the phase gradient terms in the free energy. Because of the kernel structure, that fulfills the criteria (i)-(iii), a non-trivial phase crystallization occurs at a particular $q_{x}^{*}\sim 1/\xi_{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.Fogelström et al. (1997); Higashitani (1997) That work assumed translational invariance of the superflow and currents along the surface, which guaranteed particle conservation \mbox{\boldmath\nabla}\cdot{\bf j}({\bf R})=0, but as a result required additional mechanisms of reducing superflow in the bulk. In semi-infinite systems one relies on the Meissner effect to screen the bulk superflow on the penetration depth length scale $\lambda$, which leads to $T^{*}\sim(\xi_{0}/\lambda)T_{c}$.Barash et al. (2000); Löfwander et al. (2000) In slabs of width $D<\lambda$ the bulk contribution is obviously limited, resulting in spontaneous superflow below $T^{*}\sim(\xi_{0}/D)T_{c}$. Vorontsov (2009) 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^{*}\sim(\xi_{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 polar $p$-wave which may also host a flat band of zero-energy surface fermions. Interestingly, phase crystallization can happen in conventional $s$-wave superconductors, where orbital pairbreaking scattering is absent. In this case, magnetically active interfaces can provide the proper environment for the phase instability, for example in superconductor-ferromagnetic structures. Such systems are being considered as important building blocks for spintronics applications, where non-locality and quantum coherence will play important roles.Eschrig (2011) As described in Appendix C, a similar form of the superfluid density tensor appear for $\vartheta=\pi$ 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 Novelli et al. (2019); Wang et al. (2017). Moreover, softening of the surface superfluid density at some finite wavevector 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 with 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.Cross and Hohenberg (1993) We note that the convection roll currents in that case is due to an instability in a non-equilibrium driven system, while the phase crystal is a second-order phase transition into a new ground state.

IV Conclusions

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 non-local 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 Holmvall et al. (2018) and survives significant reduction of spectral weight of bound states Holmvall et al. (2019). 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 superconductor-ferromagnet 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.Boninsegni and Prokof’ev (2012); Léonard et al. (2017); Böttcher et al. (2019); Chomaz et al. (2019) More generally, our results indicate that non-local effects in broken-symmetry states, especially with multi-component order parameters or competing orders, can lead to new states of matter. Such prospects are supported by earlyPippard and Bragg (1953) and more recentKoyama and Machida (2013) investigations of non-local physics in superconductors, as well as research into pattern formation due to long-range non-locality in biological systems.Tanaka and Kuramoto (2003); Bressloff and Kilpatrick (2008); García-Morales and Krischer (2008)

V Acknowledgements

The computations were performed on resources at Chalmers Centre for Computational Science and Engineering (C3SE) provided by the Swedish National Infrastructure for Computing (SNIC). We thank the Swedish Research Council 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.Serene and Rainer (1983) Our starting point is the Eilenberger equation for the quasiclassical propagator $\hat{g}$

[TABLE]

In this equation a spatially varying phase $\chi$ of the order parameter $\Delta=|\Delta|e^{i\chi({\bf R})}$, was eliminated in favor of the superflow field {\bf p}_{s}=\frac{1}{2}\mbox{\boldmath\nabla}\chi. This can always be done, if needed, by a gauge transformation $\hat{g}\to\hat{U}\hat{g}\hat{U}^{\dagger}$ with $\hat{U}=e^{i\hat{\tau}_{3}\chi/2}$. The superflow is a function of position ${\bf p}_{s}={\bf p}_{s}({\bf R})$, and we consider a singlet mean-field order parameter $\Delta=\Delta({\bf R},{\bf p}_{\mbox{\tiny F}})$. The commutator-based Eilenberger equation is transformed into the Riccati-type equations for the coherence amplitudesEschrig (2000)

[TABLE]

These amplitudes conveniently parametrize the quasiclassical propagator,Schopohl and Maki (1995); *Shelankov2000 and are functions of position, momentum, and energy, $\gamma=\gamma({\bf R},{\bf p}_{\mbox{\tiny F}};\varepsilon_{m})$. The two coherence amplitudes are related by symmetry,

[TABLE]

that also applies to other tilde-related functions. For the singlet real order parameter $\tilde{\Delta}({\bf R},{\bf p}_{\mbox{\tiny F}})\equiv\Delta^{*}({\bf R},-{\bf p}_{\mbox{\tiny F}})=\Delta({\bf R},{\bf p}_{\mbox{\tiny F}})$. We look at the current response due to a small but arbitrary superflow field ${\bf p}_{s}={\bf p}_{s}({\bf R})$, starting from a current-less background state $\Delta_{0}({\bf R},{\bf p}_{\mbox{\tiny F}})$ and the corresponding coherence amplitudes $\gamma_{0}({\bf R},{\bf p}_{\mbox{\tiny F}};\varepsilon_{m})$. The following linear response calculation is valid for any spatial profile of $\gamma_{0}({\bf R},{\bf p}_{\mbox{\tiny F}};\varepsilon_{m})$, and we specify in the end its particular form. The current at a point ${\bf R}$ near the surface is calculated from the correction to the diagonal propagator $\delta g$, with $g=-i\pi\mbox{sgn}(\varepsilon_{m})\frac{1-\gamma\tilde{\gamma}}{1+\gamma\tilde{\gamma}}$, as

[TABLE]

where $N_{\mathrm{F}}$ is density of states at the Fermi level per spin projection, and $\left\langle\dots\right\rangle_{{\hat{p}}}=\int d{\hat{p}}/2\pi\,\dots$ denotes a cylindrical Fermi surface average, Fig. 5. In terms of linearised coherence amplitudes $\gamma=\gamma_{0}+\gamma_{1}$ the propagator change due to small superflow is

[TABLE]

We first neglect the effect of the superflow on the amplitude of the order parameter, assuming that $\Delta({\bf R})=\Delta_{0}({\bf R})$ even in the current-carrying state, and linearise Eqs. (16) to find transport equations for the function $\gamma_{1}/(1+\gamma_{0}\tilde{\gamma}_{0})$,

[TABLE]

We get a similar equation for the tilde-analogue. The parameter

[TABLE]

determines the correlation length of the response. In a uniform state it reduces to $\kappa=2v_{\mbox{\tiny F}}^{-1}\sqrt{\Delta_{\hat{p}}^{2}+\varepsilon_{m}^{2}}\sim 1/\xi_{0}$.

The solution of Eq. (20) along a quasiclassical trajectory $s$ is found, for positive $\varepsilon_{m}$, by integration forward along the trajectory starting from zero value in the bulk $\gamma_{1}(s=-\infty)=0$, where there is no superflow. We get

[TABLE]

To write the current at the observation point ${\bf R}$ we need to integrate over all trajectories coming into point ${\bf R}$. By introducing a correlation function connecting two points, ${\bf R}_{1}$ and ${\bf R}_{2}$, by a quasiclassical trajectory $\hat{\rho}=({\bf R}_{2}-{\bf R}_{1})/|{\bf R}_{2}-{\bf R}_{1}|$,

[TABLE]

one can combine the Fermi surface average at the observation point and integration along trajectories into integration over all space ${\bf R}^{\prime}$, see Fig. 6, and write the current response as

[TABLE]

Inserting (22) into (19) and using definition (23), the superfluid kernel is then given by

[TABLE]

where $f_{0}$ and $\tilde{f}_{0}$ are off-diagonal propagators in the unperturbed state. In terms of coherence amplitudes $f_{0}=-2i\pi\,\mbox{sgn}(\varepsilon_{m})\frac{\gamma_{0}}{1+\gamma_{0}\tilde{\gamma}_{0}}$. This kernel connects the observation point ${\bf R}$ to the integration point ${\bf R}^{\prime}$. For each pair of points there are two paths, one direct \mbox{\tiny1}⃝ and one involving reflection at the surface \mbox{\tiny2}⃝, where we assumed mirror-like reflection, see Fig. 5. The momentum direction ${\hat{p}}$ at the observation point is given by the trajectory direction ${\bf R}^{\prime}\to{\bf R}$, and similarly for momentum at the integration point ${\hat{p}}^{\prime}$ (Fig. 5). These directions are different for the direct and reflected paths.

Appendix B Coherence amplitudes and propagators with a step-like order parameter

Neglecting the suppression of the order parameter at the surface allows us to proceed further analytically. The bulk uniform coherence amplitude is

[TABLE]

Now consider, Fig. 7, a (straightened) trajectory that for $s<0$ is in a region with the order parameter $\Delta_{\underline{k}}=\Delta_{i}$, and for $s>0$ is in the region with $\Delta_{k}=\Delta_{f}$ (e.g. for the most pairbreaking surface $\Delta_{i}=-\Delta_{f}$). Denote

[TABLE]

Far away from the interface, the coherence amplitudes have their uniform bulk values (we assume $\varepsilon_{m}>0$, otherwise understand $\varepsilon_{m}=|\varepsilon_{m}|$ and add $\mbox{sgn}(\varepsilon_{m})$ in front)

[TABLE]

For a sudden-step order parameter the amplitudes $\gamma_{0},\tilde{\gamma}_{0}(s)$ can be found analytically, integrating Riccati equations (16) in forward or backward direction, correspondingly. Including the sudden jump of the amplitudes at the surface according to the boundary condition, we get

[TABLE]

and for tilde-function integrating backward:

[TABLE]

The propagators on the trajectory are (e.g. for $s>0$)

[TABLE]

and the off-diagonal component that enters the expression for the current response is

[TABLE]

where we wrote 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 $\tilde{\gamma}_{f}$ in both terms since it is coming from $s=+\infty$, but the $\gamma$-amplitude can be either $\gamma_{f}$ far from the reflection point or $\Gamma_{f}\leftarrow\gamma_{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, that mixes $\Gamma_{f}$ and $\tilde{\gamma}_{f}$ in denominator, is the one that mainly determines bound states effects. In both diagonal and off-diagonal items the continuum and the bound states 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 $\kappa$ that determines the correlations extent in the current response:

[TABLE]

Here we consider an order parameter orientation such that the amplitudes on the incoming and reflected parts of the trajectory are the same, so $\kappa_{u,i}=\kappa_{u,f}=\kappa_{u}$. The generalization for different amplitudes can be easily carried out retaining indices $\Omega_{i,f}$, $\kappa_{u;i,f}$ etc. This expression for $\kappa(s)$ is quite general and easy to integrate along trajectories, as required for correlation functions $C({\bf R},{\bf R}^{\prime})$ and $C({\bf R}^{\prime},{\bf R})$. In both these functions integration goes from initial to final point as determined by the momentum direction, and is shown in Fig. 8.

For the case (a) both $s_{1}$ and $s_{2}$ are on the same side of the interface and $s_{2}$ is further away from the interface than $s_{1}$, we have

[TABLE]

If we reverse the trajectory the signs of $s$ change (so that $s_{1}$ and $s_{2}$ determine absolute distance to the surface)

[TABLE]

For the (c) case we break the integral into two parts for in and out $s-in-out:$

[TABLE]

The denominators in (34-36) will cancel numerators in some of the $f$-functions (32) when combined in the kernel expression (25). The numerators in (34-35) can be written as

[TABLE]

For any given points ${\bf R}$ and ${\bf R}^{\prime}$ we define two paths, direct and reflected, and each will have ${\bf R}\to{\bf R}^{\prime}$ and ${\bf R}^{\prime}\to{\bf R}$ contributions, $\tilde{f}({\hat{p}},{\bf R})C({\bf R},{\bf R}^{\prime})f({\hat{p}},{\bf R}^{\prime})+\tilde{f}(-{\hat{p}},{\bf R})C({\bf R}^{\prime},{\bf R})f(-{\hat{p}},{\bf R})$. Let’s denote by ${\hat{k}}$ momentum away from the surface, and in this case we identify indices $f={\hat{k}}$, $i=\underline{{\hat{k}}}$. The trajectory we are integrating $\gamma$-function goes from $s_{1}=s_{<}$ (point closest to the interface) to $s_{2}=s_{>}$ (point farthest from interface). For reverse trajectory we have $f=-\underline{{\hat{k}}}$, $i=-{\hat{k}}$ and integration happens from $-s_{2}$ to $-s_{1}$.

The two terms give, after mentioned cancellations, for direct path

[TABLE]

For the reflected path this sum has a more compact form that directly reflects the bound states factors

[TABLE]

Note, that to generalize for inequivalent gap size on in-out trajectories we need to use appropriate $\kappa_{u}$ along given directions, e.g. $\kappa_{u}|s^{\prime}_{>}+s^{\prime}_{<}|\to\kappa_{u,{\hat{k}}^{\prime}}s^{\prime}_{>}+\kappa_{u,\underline{{\hat{k}}}^{\prime}}s^{\prime}_{<}$ for trajectory $\underline{{\hat{k}}}^{\prime}\to{\hat{k}}^{\prime}$ with reflection. These are completely general expressions for the one-component order parameters, where we neglect suppression of OP amplitude near the surface, and assume specular scattering.

We apply the developed formalism and approximations to a $d$-wave superconductor with maximally pairbreaking surface. In this case we have $\Delta_{\hat{k}}=-\Delta_{\underline{{\hat{k}}}}$ for all incident trajectories, and $\gamma_{-{\hat{k}}}=\gamma_{{\hat{k}}}=-\gamma_{\underline{{\hat{k}}}}$, $\tilde{\gamma}_{-{\hat{k}}}=\tilde{\gamma}_{{\hat{k}}}=\tilde{\gamma}_{\underline{{\hat{k}}}}$, $\Gamma_{\hat{k}}=\gamma_{\underline{{\hat{k}}}}=-\gamma_{\hat{k}}$, and two important combinations of the coherence amplitudes are

[TABLE]

The correlation coefficient Eq. (33) along a trajectory $s$ is

[TABLE]

where $\kappa_{u}=2\Omega/v_{\mbox{\tiny F}}$ and $\Omega=\sqrt{\varepsilon_{m}^{2}+\Delta_{\hat{k}}^{2}}$. The distance along a trajectory, measured from the surface, is $s=y/{\hat{k}}_{y}$. One uses these relations for coherence amplitudes in combinations (38) and (39) to find the kernel (25) components, as given in the main text, for the direct path, Eq. (8), and the reflection path, Eq. (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 amplitudesEschrig (2009)

[TABLE]

with ${\cal M}=e^{i\vartheta\hat{m}\cdot{\bf\sigma}/2}$ and $\tilde{\cal M}={\cal M}^{*}$. Magnetic spin mixing leads to the bound states $\varepsilon_{b}=\pm\Delta\cos(\vartheta/2)$, that result in zero energy states for $\vartheta=\pi$ and the boundary condition for coherence amplitudes $\Gamma_{{\hat{k}}}=-\gamma_{\underline{{\hat{k}}}}$.