format: text

text:
# Zero-Bias States and the Mechanism of the Surface 𝑑→𝑑+𝑖𝑠 Transition

## I Introduction

It has now been firmly established that the order parameter $`\mathrm{\Delta }`$ in homogeneous cuprate superconductors has a $`d`$-wave symmetry . It follows that inhomogeneities can scatter quasiparticles between directions that experience opposite signs of $`\mathrm{\Delta }`$. This effect is strongest at a specularly reflecting (110) surface, because $`\mathrm{\Delta }`$ changes sign along each quasiclassical trajectory upon reflection (see Fig. 1). As a consequence of the Atiyah-Patodi-Singer index theorem , the Andreev equation along each such trajectory has in its spectrum a bound state at exactly zero energy, irrespective of the detailed shape of the potential . Collectively, these states then make up a peak in the tunneling spectra, which has been observed experimentally .

It was predicted theoretically  and confirmed experimentally  that below a certain temperature of the order of 1K, this peak of the zero-bias states (ZBSs) splits. The theoretical interpretation of this splitting is that via a sub-dominant pairing interaction of a different symmetry, say $`s`$, a subdominant order parameter is induced close to the surface. It is phase-shifted by $`\pi /2`$ relative to the dominant $`d`$-wave, which gives the total order parameter “$`d+is`$” symmetry and indicates broken time-reversal symmetry.

The full self-consistent calculations in the Eilenberger formalism that are based on this interpretation  are in a good quantitative agreement with the experimental data. However, the mechanism of the transition is not manifest in the numerical solutions of the Eilenberger equations. We believe that understanding of the basic physics of the $`dd+is`$ transition is especially needed now in light of recent experiments that call this interpretation into question . That is the purpose of the work presented in this paper.

Below, we show by a simple argument based on first-order perturbation theory that the degrees of freedom driving the $`dd+is`$ transition are the ZBSs, and that we can neglect the effect of all the remaining states. Hence, to understand the mechanism of the transition, we have to deal only with the ZBSs, which is convenient since these states are least sensitive to the unknown surface details.

The paper is organized as follows: In Section II, we demonstrate our strategy on the familiar case of BCS instability. The main argument is presented in Section III after we have extended the BCS formalism to inhomogeneous systems and non-$`s`$-wave pairing. Based on this argument, we calculate $`\mathrm{\Delta }`$ at $`T=0`$ in Section IV and estimate the transition temperature to the $`d+is`$ state in Section V. In Section VI, we discuss the surface current. Finally, we discuss our results in Section VII.

## II BCS Instability

There are various ways to consider the energetic costs and benefits of the transition to the superfluid state. The one that has proven useful in our study of the $`dd+is`$ transition is to decouple the attractive four-fermion interaction by the Hubbard-Stratonovich (HS) transformation, and to make a saddle-point (mean-field) approximation. That way, we break up the total free energy of the system into free energy of single particle states, which is lowered by the gap $`\mathrm{\Delta }`$, and the extra term from the HS transformation, which grows (quadratically) with $`\mathrm{\Delta }`$. We then see that at small enough $`T`$, the system favors transition to the superfluid state.

We will demonstrate this on the familiar BCS case. The model Hamiltonian is ***Since there is no universal convention as to whether the attractive interaction term should have a plus sign with a negative coupling constant $`V`$ or a minus sign with positive $`V`$, we use $`|V|`$ which is unambiguously positive.

$$H=\underset{𝐤,\sigma }{}ϵ_𝐤c_{𝐤,\sigma }^{}c_{𝐤,\sigma }|V|\underset{𝐤,𝐤^{}}{}c_𝐤^{}c_𝐤^{}c_𝐤^{}c_𝐤^{},$$
(1)
which gives rise to the partition function

$$Z=𝒟\overline{c}_{𝐤\sigma }𝒟c_{𝐤\sigma }e^{\underset{0}{\overset{\beta }{}}𝑑\tau \left[\underset{𝐤\sigma }{}\overline{c}_{𝐤\sigma }(_\tau +ϵ_𝐤)c_{𝐤\sigma }|V|\underset{𝐤,𝐤^{}}{}\overline{c}_𝐤\overline{c}_𝐤c_𝐤^{}c_𝐤^{}\right]},$$
(2)
where the $`c`$’s are now $`\tau `$dependent Grassmann numbers. We perform the HS transformation by multiplying the partition function by the (infinite) constant

$$𝒟\overline{\varphi }_𝐤𝒟\varphi _𝐤e^{|V|\underset{0}{\overset{\beta }{}}𝑑\tau \underset{𝐤,𝐤^{}}{}(\overline{\varphi }_𝐤\overline{c}_𝐤\overline{c}_𝐤)(\varphi _𝐤^{}c_𝐤^{}c_𝐤^{})},$$

so

$$Z=𝒟\overline{\varphi }_𝐤𝒟\varphi _𝐤𝒟\overline{c}_{𝐤\sigma }𝒟c_{𝐤\sigma }e^S,$$
(3)
where

$$S=\underset{0}{\overset{\beta }{}}𝑑\tau \left[\underset{𝐤}{}\left(\begin{array}{cc}\overline{c}_𝐤& c_𝐤\end{array}\right)\left(\begin{array}{cc}_\tau +ϵ_𝐤& \mathrm{\Delta }\\ \overline{\mathrm{\Delta }}& _\tau ϵ_𝐤\end{array}\right)\left(\begin{array}{c}c_𝐤\\ \overline{c}_𝐤\end{array}\right)+\frac{|\mathrm{\Delta }|^2}{|V|}\right],$$
(4)
where we defined

$$\mathrm{\Delta }=|V|\underset{𝐤}{}\varphi _𝐤.$$

From the action (4), we can read off that in the mean-field approximation, the total free energy of the system is

$$F(|\mathrm{\Delta }|)=\underset{𝐤}{}((E_𝐤)+(E_𝐤))+\frac{|\mathrm{\Delta }|^2}{|V|}$$
(5)
upon minimization with respect to $`|\mathrm{\Delta }|`$. Here,

$$(E)=T\mathrm{ln}(1+e^{E/T}),$$
(6)
and

$$E_𝐤=\sqrt{ϵ_𝐤^2+|\mathrm{\Delta }|^2}.$$

We see the instability most clearly at $`T=0`$, where $`F(|\mathrm{\Delta }|)=E(|\mathrm{\Delta }|)`$. Then

$$(E)=\theta (E)E,$$

so

$$E(|\mathrm{\Delta }|)E(0)=N(0)\underset{\omega _D}{\overset{0}{}}𝑑ϵ(\sqrt{ϵ^2+|\mathrm{\Delta }|^2}ϵ)+\frac{|\mathrm{\Delta }|^2}{|V|},$$
(7)
where $`N(0)`$ is the density of states at the Fermi level, and $`\omega _D`$ is the Debye frequency. Direct calculation shows that the integral behaves as $`|\mathrm{\Delta }|^2\mathrm{ln}\frac{|\mathrm{\Delta }|}{\omega _D}`$ for $`|\mathrm{\Delta }|0`$, whose non-analytic decrease will win over the analytic increase of the second term for small enough $`\mathrm{\Delta }`$, no matter how weak the attractive interaction $`|V|`$ is. By the same calculation, we can also see that the integral becomes analytic if we do not integrate $`ϵ`$ all the way up to zero, but to a finite negative energy. This means that the states close to the Fermi energy drive the BCS transition—they benefit most from opening of the gap $`|\mathrm{\Delta }|`$. Similarly, we shall see that the states at zero energy, that is the ZBSs, will drive the $`dd+is`$ transition.

So far, the argument has shown the BCS instability only at $`T=0`$. At finite temperatures,

$$(E_𝐤)+(E_𝐤)=T\mathrm{ln}(2+2\mathrm{cosh}\frac{E_𝐤}{T}),$$

which, upon expansion in powers of $`|\mathrm{\Delta }|^2`$, gives

$$F(|\mathrm{\Delta }|)F(0)=|\mathrm{\Delta }|^2\left(\frac{1}{|V|}N(0)\underset{\omega _D}{\overset{0}{}}\frac{dϵ}{ϵ}\mathrm{tanh}\frac{ϵ}{2T}\right)+O(|\mathrm{\Delta }|^4).$$
(8)
This shows that the system is unstable to the BCS transition at temperatures below $`T_c`$ that satisfies

$$\frac{1}{|V|}N(0)\underset{\omega _D}{\overset{0}{}}\frac{dϵ}{ϵ}\mathrm{tanh}\frac{ϵ}{2T_c}=0.$$
(9)
In a similar way, we shall see below that $`T_s`$, the transition temperature into the $`d+is`$ state, is finite.

## III The $`d+is`$ Instability

### A Formalism

We now need to develop the formalism that will enable us to extend the strategy from Section II to the $`d+is`$ case. We shall consider a single (2-dimensional) CuO plane, and model it by the Hamiltonian

$$H=d^2r\underset{\sigma =,}{}\psi _\sigma ^{}(𝐫)ϵ(i)\psi _\sigma (𝐫)+d^2rd^2r^{}V(𝐫𝐫^{})\psi _{}^{}(𝐫)\psi _{}^{}(𝐫^{})\psi _{}(𝐫^{})\psi _{}(𝐫),$$
(10)
where $`ϵ`$ is the band energy and $`V`$ is the short-range interaction responsible for pairing. What makes this difficult problem tractable is the separation of energy scales (the Fermi energy $`E_F`$ is much bigger than the superconducting gap $`\mathrm{\Delta }`$), which gives rise to separation of length scales $`\lambda _F`$ (Fermi wave length) and $`\xi `$ (the coherence length). We may, therefore, expand in powers of the small parameter $`\lambda _F/\xi `$; keeping the lowest non-trivial order is called the quasiclassical approximation. This procedure is usually done at the level of Green’s function , which are thus transformed into Eilenberger functions that satisfy transport-like equations.

Since we want to understand the $`dd+is`$ transition in terms of quasiparticle eigenstates rather than Green’s functions, we will perform this separation of scales at the operator level instead. We denote as $`2\mathrm{\Lambda }`$ the width of the shell around the Fermi surface containing the states that take part in the pairing (see Fig. 2). We then factor out the fast Fermi-surface oscillations and define the slowly varying field operator $`\psi _{\sigma ,\theta }(𝐫)`$  by

$`\psi _\sigma (𝐫)`$ $`=`$ $`{\displaystyle \frac{d^2k}{(2\pi )^2}c_{𝐤\sigma }e^{i𝐤𝐫}}`$ (11)
$``$ $`{\displaystyle \underset{F.S.}{}}{\displaystyle \frac{dk_F(\theta )}{2\pi }}\left({\displaystyle \underset{\mathrm{\Lambda }}{\overset{\mathrm{\Lambda }}{}}}{\displaystyle \frac{dk_{}}{2\pi }}c_{𝐤\sigma }e^{ik_{}𝐧(\theta )𝐫}\right)e^{i𝐤_F(\theta )𝐫}`$ (12)
$``$ $`{\displaystyle \underset{F.S.}{}}{\displaystyle \frac{dk_F(\theta )}{2\pi }}\psi _{\sigma ,\theta }(𝐫)e^{i𝐤_F(\theta )𝐫}.`$ (13)
When we substitute this into (10), we obtain

$`H`$ $`=`$ $`{\displaystyle }d^2r[{\displaystyle \underset{\sigma }{}}{\displaystyle \underset{F.S.}{}}{\displaystyle \frac{dk_F(\theta )}{2\pi }}\psi _{\sigma ,\theta }^{}(𝐫)𝐯_F(\theta )(i)\psi _{\sigma ,\theta }(𝐫)+`$ (14)
$`+`$ $`{\displaystyle \underset{F.S.}{}}{\displaystyle \frac{dk_F(\theta )}{2\pi }}{\displaystyle \frac{dk_F(\theta ^{})}{2\pi }}V(\theta ,\theta ^{})\psi _\theta ^{}(𝐫)\psi _\theta ^{}(𝐫)\psi _\theta ^{}(𝐫)\psi _\theta ^{}(𝐫)],`$ (15)
where $`𝐯_F(\theta )`$ is the Fermi velocity at the point $`𝐤_F(\theta )`$, and

$$V(\theta ,\theta ^{})d^2re^{i(𝐤_F(\theta )𝐤_F(\theta ^{}))𝐫}V(𝐫).$$

The derivation of (14) is given in Appendix A. Note the linearized kinetic energy in (14), which will be crucial in the following.

The Hamiltonian (14) gives rise to a partition function, which we can write as a path integral over the fermion fields $`\psi _{\sigma ,\theta }(𝐫)`$. We can again decompose the interaction by the HS transformation, i.e., we can multiply the partition function by the constant

$`{\displaystyle }𝒟\overline{\varphi }_\theta (𝐫)𝒟\varphi _\theta (𝐫)\mathrm{exp}[{\displaystyle \underset{0}{\overset{\beta }{}}}d\tau {\displaystyle }d^2r{\displaystyle \underset{F.S.}{}}{\displaystyle \frac{dk_F(\theta )}{2\pi }}{\displaystyle \frac{dk_F(\theta ^{})}{2\pi }}V(\theta ,\theta ^{})\times `$ (16)
$`\times (\overline{\varphi }_\theta (𝐫)\overline{\psi }_\theta (𝐫)\overline{\psi }_\theta (𝐫))(\varphi _\theta ^{}(𝐫)\psi _\theta ^{}(𝐫)\psi _\theta ^{}(𝐫))].`$ (17)
In the mean-field approximation, the total free energy of the system equals the free energy given by the (single-particle) Hamiltonian

$`H`$ $`=`$ $`{\displaystyle }d^2r[{\displaystyle \underset{F.S.}{}}{\displaystyle \frac{dk_F(\theta )}{2\pi }}\left(\begin{array}{cc}\psi _\theta ^{}(𝐫)& \psi _\theta (𝐫)\end{array}\right)\left(\begin{array}{cc}𝐯_F(\theta )(i)& \mathrm{\Delta }_\theta (𝐫)\\ \mathrm{\Delta }_\theta ^{}(𝐫)& 𝐯_F(\theta )(i)\end{array}\right)\left(\begin{array}{c}\psi _\theta (𝐫)\\ \psi _\theta ^{}(𝐫)\end{array}\right)`$ (18)
$``$ $`{\displaystyle \underset{F.S.}{}}{\displaystyle \frac{dk_F(\theta )}{2\pi }}{\displaystyle \frac{dk_F(\theta ^{})}{2\pi }}V(\theta ,\theta ^{})\varphi _\theta ^{}(𝐫)\varphi _\theta ^{}(𝐫)],`$ (19)
upon minimization with respect to $`\varphi _\theta (𝐫)`$, where we defined

$$\mathrm{\Delta }_\theta (𝐫)=\underset{F.S.}{}\frac{dk_F(\theta )}{2\pi }V(\theta ,\theta ^{})\varphi _\theta ^{}(𝐫).$$
(20)
We shall write explicit formulae for the total energy and free energy in Sections IV and V (formulae (34) and (47)). Here we just note that to calculate the single-particle contribution to the free energy, we will have to find the spectra of the Andreev Hamiltonians labeled by $`\theta `$, i.e., we will need the energies $`E_{\theta ,n}`$ that satisfy

$$\left(\begin{array}{cc}𝐯_F(\theta )(i)& \mathrm{\Delta }_\theta (𝐫)\\ \mathrm{\Delta }_\theta ^{}(𝐫)& 𝐯_F(\theta )(i)\end{array}\right)\left(\begin{array}{c}f_{\theta ,n}(𝐫)\\ g_{\theta ,n}(𝐫)\end{array}\right)=E_{\theta ,n}\left(\begin{array}{c}f_{\theta ,n}(𝐫)\\ g_{\theta ,n}(𝐫)\end{array}\right).$$
(21)
We note that the linear kinetic energy in (14) makes this equation effectively one-dimensional, i.e., an independent equation for each line in the direction $`𝐯_F(\theta )`$. In the presence of the specularly reflecting boundary, we must find the Andreev spectra along reflected lines such as the one in Fig. 1a. Equivalently, we solve the equation on a straight line with the pairing potential $`\mathrm{\Delta }`$ shown in Fig. 1b. This is intuitively obvious; a derivation is given in Appendix B.

As we mentioned in the Introduction, the spectrum along each trajectory having opposite signs of $`\mathrm{\Delta }`$ at the two asymptotic ends will contain a zero-bias state. Its wave function is, up to a normalization constant,

$$\left(\begin{array}{c}f(\theta ,\rho )\\ g(\theta ,\rho )\end{array}\right)_{ZBS}=\left(\begin{array}{c}1\\ i\end{array}\right)\mathrm{exp}(\underset{0}{\overset{\rho }{}}𝑑\rho ^{}\mathrm{\Delta }(\theta ,\rho ^{})),$$
(22)
where the upper (lower) sign corresponds to $`\mathrm{\Delta }(\theta ,\rho =\mathrm{})<0`$, $`\mathrm{\Delta }(\theta ,\rho =+\mathrm{})>0`$ ($`\mathrm{\Delta }(\theta ,\rho =\mathrm{})>0`$, $`\mathrm{\Delta }(\theta ,\rho =+\mathrm{})<0`$), so that the wave function is normalizable. In our notation, we will freely interchange the dependence on $`𝐫`$ (actually only on $`x`$, since the system is translationally invariant in the $`y`$-direction) with the dependence on the angle $`\theta `$ and the coordinate $`\rho `$ along the trajectory. Their relationship is obvious from Fig. 1.

### B Argument

We now have all the tools needed to demonstrate the $`d+is`$ transition in a way that brings out its physical mechanism. We follow the same line of thought as in Section II: We go to the zero temperature, and look at the energy gains and losses when the $`s`$-wave component of $`\mathrm{\Delta }`$ appears.

For any $`s`$-wave pairing to appear, it is necessary that the part of the functional integral (16) over the $`s`$-component of $`\varphi `$ converge, i.e., that $`V`$ on top of the dominant $`d`$-wave attraction contain also an $`s`$-wave part We use again $`|V_s|`$ rather than $`V_s`$.,

$$V(\theta ,\theta ^{})=V_d(\theta ,\theta ^{})|V_s|.$$

In this Subsection, we will show that this condition is also sufficient: At zero temperature, the system will favor the $`d+is`$ state for an arbitrarily weak attraction $`V_s`$.

With both $`d`$\- and $`s`$-wave pairing present,

$$\varphi _\theta (𝐫)=\varphi _{d\theta }(𝐫)+\varphi _s(𝐫).$$

(The $`s`$-components of both $`V`$ and $`\varphi `$ are angle-independent.) We begin with the second term in (18), which then is

$``$ $`{\displaystyle \underset{F.S.}{}}{\displaystyle \frac{dk_F(\theta )}{2\pi }}{\displaystyle \frac{dk_F(\theta ^{})}{2\pi }}(V_d(\theta ,\theta ^{})|V_s|)(\varphi _{d\theta }^{}(𝐫)+\varphi _s^{}(𝐫))(\varphi _{d\theta ^{}}(𝐫)+\varphi _s(𝐫))=`$ (23)
$`=`$ $`{\displaystyle \underset{F.S.}{}}{\displaystyle \frac{dk_F(\theta )}{2\pi }}{\displaystyle \frac{dk_F(\theta ^{})}{2\pi }}V_d(\theta ,\theta ^{})\varphi _{d\theta }^{}(𝐫)\varphi _{d\theta ^{}}(𝐫)+`$ (24)
$`+`$ $`|V_s|{\displaystyle \underset{F.S.}{}}{\displaystyle \frac{dk_F(\theta )}{2\pi }}{\displaystyle \frac{dk_F(\theta ^{})}{2\pi }}\varphi _s^{}(𝐫)\varphi _s(𝐫),`$ (25)
where we used the orthogonality of the $`s`$ and $`d`$ components:

$$\underset{F.S.}{}\frac{dk_F(\theta )}{2\pi }V_d(\theta ,\theta ^{})=\underset{F.S.}{}\frac{dk_F(\theta ^{})}{2\pi }V_d(\theta ,\theta ^{})=\underset{F.S.}{}\frac{dk_F(\theta )}{2\pi }\varphi _{d\theta }(𝐫)=0.$$

We can also split up (20) into components and define

$`\mathrm{\Delta }_{d\theta }(𝐫)`$ $`=`$ $`{\displaystyle \underset{F.S.}{}}{\displaystyle \frac{dk_F(\theta ^{})}{2\pi }}V_d(\theta ,\theta ^{})\varphi _{d\theta ^{}}(𝐫)`$ (26)
$`\mathrm{\Delta }_s(𝐫)`$ $`=`$ $`{\displaystyle \underset{F.S.}{}}{\displaystyle \frac{dk_F(\theta )}{2\pi }}(|V_s|)\varphi _s(𝐫).`$ (27)
The $`d`$-component of $`\mathrm{\Delta }`$ was established well above $`T_s`$, so the change of the second term in (18) due to the opening of a (small) $`s`$-wave gap will be

$$\frac{|\mathrm{\Delta }_s(𝐫)|^2}{|V_s|}$$
(28)
just as in the BCS case. Due to the translational invariance in the $`y`$-direction, we will from now on write $`\mathrm{\Delta }_s(𝐫)\mathrm{\Delta }_s(x)`$. Along the quasiclassical trajectory, $`x`$ depends on both $`\rho `$ and $`\theta `$ (see Fig.1), so we will then write $`\mathrm{\Delta }_s(\theta ,\rho )`$.

To examine the effect of the small $`s`$ wave component on the quasiparticle energies, we need to look at the change of the spectra of the $`1D`$ Andreev problems

$$\left(\begin{array}{cc}iv_F(\theta )_\rho & \mathrm{\Delta }_d(\theta ,\rho )\\ \mathrm{\Delta }_d(\theta ,\rho )& iv_F(\theta )_\rho \end{array}\right)\left(\begin{array}{c}f_n(\theta ,\rho )\\ g_n(\theta ,\rho )\end{array}\right)=E_{\theta ,n}\left(\begin{array}{c}f_n(\theta ,\rho )\\ g_n(\theta ,\rho )\end{array}\right).$$
(29)
upon $`\mathrm{\Delta }_d(\theta ,\rho )\mathrm{\Delta }_d(\theta ,\rho )+\mathrm{\Delta }_s(\theta ,\rho )`$. As $`\mathrm{\Delta }_s`$ is small, it can be treated as a perturbation; then the change of the quasi-particle energies to the lowest order is

$`E_{\theta ,n}^{(1)}[\mathrm{\Delta }_s]`$ $`=`$ $`{\displaystyle \underset{\mathrm{}}{\overset{+\mathrm{}}{}}}𝑑\rho \left(\begin{array}{cc}f_n^{}(\theta ,\rho )& g_n^{}(\theta ,\rho )\end{array}\right)\left(\begin{array}{cc}0& \mathrm{\Delta }_s(\theta ,\rho )\\ \mathrm{\Delta }_s^{}(\theta ,\rho )& 0\end{array}\right)\left(\begin{array}{c}f_n(\theta ,\rho )\\ g_n(\theta ,\rho )\end{array}\right)=`$ (30)
$`=`$ $`{\displaystyle \underset{\mathrm{}}{\overset{+\mathrm{}}{}}}𝑑\rho [f_n^{}(\theta ,\rho )g_n(\theta ,\rho )\mathrm{\Delta }_s(\theta ,\rho )+g_n^{}(\theta ,\rho )f_n(\theta ,\rho )\mathrm{\Delta }_s^{}(\theta ,\rho )].`$ (31)
Let us first look at the change of energy of the zero-energy bound states. Then from (22)

$$g_{ZBS}(\theta ,\rho )=if_{ZBS}(\theta ,\rho ),$$
(32)
so

$$E_{\theta ,ZBS}^{(1)}[\mathrm{\Delta }_s]=\pm \underset{\mathrm{}}{\overset{+\mathrm{}}{}}𝑑\rho |f(\theta ,\rho )|^22Im\mathrm{\Delta }_s(\theta ,\rho ),$$
(33)
where the upper (lower) sign corresponds to the $`+y`$\- ($`y`$-)moving trajectory. We notice several things by looking at (33):

$``$ It depends only on $`Im\mathrm{\Delta }_s`$, since $`Re\mathrm{\Delta }_s`$ just changes the position of the node in the total $`\mathrm{\Delta }(\theta ,\rho )`$, in which case the bound state remains at zero energy. Hence, we will assume $`Re\mathrm{\Delta }_s=0`$, and write $`\mathrm{\Delta }_s(\theta ,\rho )=is(\theta ,\rho )`$.

$``$ $`E_{\theta ,ZBS}^{(1)}[\mathrm{\Delta }_s]`$ is non-zero due to the form of the bound-state wave function (32) and due to the fact that $`s(\theta ,\rho )`$ does not change sign along the quasiclassical trajectory (by virtue of the $`s`$-symmetry). Out of the two possibilities for the sign of $`s`$, we will choose $`s(\theta ,\rho )>0`$ in the following, which means all the $`+y`$-moving states are shifted up in energy, whereas the $`y`$-moving states are shifted down.

Since we are at zero temperature, only the states that move down from zero energy will be occupied. We can then argue similarly as in the BCS case: opening of the additional $`s`$-wave gap costs the system energy $`s^2/|V_s|`$ (from (28)) but the quasiparticles save energy $`s`$. The lowering of the quasiparticle energy is only linear in $`s`$, $`i.e.`$, not as dramatic as the non-analytic decrease in the BCS case, but nevertheless it beats the quadratic increase for small enough $`s`$. Thus, for an arbitrarily small but non-zero interaction $`|V_s|`$, $`s=0`$ cannot be a minimum of the total energy, and the additional $`s`$-wave gap phase shifted by $`\pi /2`$ from the $`d`$-wave gap will appear. From the formula (33), we see that the superconductor will benefit from opening up the gap only close to the surface where $`|f|^2`$ is effectively non-zero, so the transition into the $`d+is`$ state is a surface effect. The decay into the bulk will be discussed more quantitatively in the next section.

We should also note that the remaining states on the quasiclassical trajectories do not change this situation, that is, they do not contribute linearly to the change of the total quasiparticle energy. Due to the time-reversal symmetry in the pure $`d`$-wave state, every state on a given quasiclassical trajectory corresponds to a state of the same energy on a reversed trajectory. Indeed, if we label the coordinate along the trajectory reversed to the one in (29) as $`\stackrel{~}{\rho }=\rho `$, then the Hamiltonian on the reversed trajectory is

$$\left(\begin{array}{cc}iv_F(\theta )_{\stackrel{~}{\rho }}& \mathrm{\Delta }_d(\theta ,\stackrel{~}{\rho })\\ \mathrm{\Delta }_d(\theta ,\stackrel{~}{\rho })& iv_F(\theta )_{\stackrel{~}{\rho }}\end{array}\right)=\left(\begin{array}{cc}iv_F(\theta )_\rho & \mathrm{\Delta }_d(\theta ,\rho )\\ \mathrm{\Delta }_d(\theta ,\rho )& iv_F(\theta )_\rho \end{array}\right)$$

since $`v_F(\theta )=v_F(\theta )`$, so

$$\left(\begin{array}{c}f_n(\theta ,\stackrel{~}{\rho })\\ g_n(\theta ,\stackrel{~}{\rho })\end{array}\right)=\left(\begin{array}{c}g_n(\theta ,\rho )\\ f_n(\theta ,\rho )\end{array}\right)$$

will also have energy $`E_{\theta ,n}`$. Now (30) implies that to the first order, a small $`Im\mathrm{\Delta }_s`$ will shift the energies of the two corresponding states by an equal amount with opposite signs. Hence, the only way they can linearly contribute to the total energy at $`T=0`$ is when one of them crosses zero and thus changes its occupancy, which happens only when their original energy (in absolute value) is smaller than the $`s`$-wave gap. But as $`s0`$, there will be fewer and fewer such states in smaller and smaller neighborhoods of the $`d`$-wave nodes. It is only the ZBSs that change their occupancy for arbitrarily small $`s`$. We thus conclude that the onset of the transition into the $`d+is`$ state is driven by these states.

## IV $`s`$-wave gap at $`T=0`$

In the study of the instability of the pure $`d`$ state in the last section, we used first-order perturbation theory since $`s0`$ at the onset of $`d+is`$. Now we will argue that this theory holds up to the the actual value of $`s`$, $`i.e.`$ $`s<<|\mathrm{\Delta }_d|`$. As we will show in Section V, $`s(T=0)T_s`$, the transition temperature into the $`d+is`$ state. Because $`T_s`$ $`7`$ K, it is much smaller than $`T_d`$, the superconducting transition temperature of the order of 100K, which sets the scale for $`\mathrm{\Delta }_d`$. Figure 3 shows the magnitude of the two gaps as a function of angle around a quarter of the Fermi surface. We see that the required inequality $`s<<|\mathrm{\Delta }_d|`$ holds for most of the Fermi surface except for small neighborhoods of the nodes. First-order perturbation theory certainly breaks down there, but upon averaging over the Fermi surface, the nodes will only introduce an error of the order $`T_s/T_d`$. Thus, we will use that theory to obtain $`s(x)`$ at $`T=0`$. As discussed at the end of Subsection III B, first-order perturbation theory implies that we have to look only at the zero-energy states. Also, since $`s`$ is a small perturbation, we shall neglect its effect on $`\mathrm{\Delta }_d`$.

Now we can write down the energy due to $`s`$ per unit length of the surface (the $`y`$-direction) as a functional of $`s(x)`$:

$$E[s(x)]=\underset{0}{\overset{+\mathrm{}}{}}𝑑x\frac{s^2(x)}{|V_s|}+\underset{\theta ϵ(\pi /2,0)}{}\frac{dk_F(\theta )}{2\pi }E_\theta [s(\theta ,\rho )]\mathrm{cos}\theta ,$$
(34)
where $`E_\theta [s]`$ is given by (33); for the rest of this section, we shall drop the superscript “(1)”, since we shall be using only the first-order formula. We freely interchange $`s(x)`$ for $`s(\theta ,\rho )`$; the relation between the two is discussed below (28). Note the correct dimensions: the $`x`$-integration makes $`s^2/|V_s|`$ from energy per unit area into energy per unit length. In the second term, the integrand is energy and the dimension of the measure is $`k_F`$, $`i.e.`$, inverse length. The extra factor of $`\mathrm{cos}\theta `$ in the second integral accounts for the difference of the density of trajectories along the $`y`$-direction compared to their angle-independent intrinsic density (measured perpendicularly to their direction), as shown in Fig. 4. In the second term in (34), we sum up only the occupied $`y`$-moving states for which $`E_\theta <0`$ according to (33).

We obtain $`s(x)`$ by minimizing (34). Let us first make an order-of-magnitude estimate

$$E[s]\xi \frac{s^2}{|V_s|}k_Fs,$$
(35)
since $`s`$ will extend into the bulk only as far as the coherence length $`\xi =\mathrm{}v_F/\mathrm{\Delta }_o`$ ($`\mathrm{\Delta }_o`$ is the amplitude of the $`d`$ wave), and from (33), we see $`E_\theta [s]s`$. The angular averaging will, up to numerical factors of order unity, multiply $`E_\theta [s]`$ by $`k_F`$. Minimization of (35) will give

$$s\frac{k_F|V_s|}{\xi }.$$
(36)
By taking $`s1`$meV from the experiment, $`k_F1`$Å-1, and $`\xi 10`$Å, we get an estimate for the strength of the $`s`$-wave pairing

$$|V_s|10\text{meVÅ}^2.$$

We minimize (34) exactly by solving

$$\frac{\delta E[s]}{\delta s(x)}=0,$$

$`i.e.`$,

$$2\frac{s(x)}{|V_s|}\underset{\theta ϵ(\pi /2,0)}{}\frac{dk_F(\theta )}{2\pi }\mathrm{cos}\theta \underset{\mathrm{}}{\overset{+\mathrm{}}{}}𝑑\rho 2|f(\theta ,\rho )|^2\frac{\delta s(\theta ,\rho )}{\delta s(x)}=0.$$

Now

$`{\displaystyle \frac{\delta s(\theta ,\rho )}{\delta s(x)}}`$ $`=`$ $`\delta (x\rho \mathrm{cos}\theta )+\delta (x+\rho \mathrm{cos}\theta )`$ (37)
$`=`$ $`{\displaystyle \frac{1}{\mathrm{cos}\theta }}\left(\delta (\rho {\displaystyle \frac{x}{\mathrm{cos}\theta }})+\delta (\rho +{\displaystyle \frac{x}{\mathrm{cos}\theta }})\right),`$ (38)
since $`\mathrm{cos}\theta >0`$, and $`\rho `$, unlike $`x`$, can be both positive and negative. The factors of $`\mathrm{cos}\theta `$ cancel, and we obtain

$$s(x)=|V_s|\underset{\theta ϵ(\pi /2,0)}{}\frac{dk_F(\theta )}{2\pi }\left(|f(\theta ,\frac{x}{\mathrm{cos}\theta })|^2+|f(\theta ,\frac{x}{\mathrm{cos}\theta })|^2\right).$$
(39)
Physically, we get two terms on the right-hand side because for each angle $`\theta `$, there are two trajectories contributing to $`s`$ at a given point as shown in Fig. 5.

We should remark here that we also obtain the formula (39) when we calculate the contribution from the occupied ($`y`$-moving) bound states to the pairing potential in the gap equation. This is done in Appendix C. The result is

$$\mathrm{\Delta }_s(x)_{ZBS}=i|V_s|\underset{\theta ϵ(\frac{\pi }{2},0)}{}\frac{dk_F(\theta )}{2\pi }\left[|f(\theta ,\frac{x}{\mathrm{cos}\theta })|^2+|f(\theta ,\frac{x}{\mathrm{cos}\theta })|^2\right]$$
(40)
in agreement with (39). This formula, however, shows more clearly the internal consistency of the picture: For $`\mathrm{\Delta }_d`$ in Fig. 1, the additional $`is`$ potential pushes down the $`y`$-moving states if $`s>0`$. As (40) shows, these states, in turn, give rise to $`\mathrm{\Delta }_s=i\times positive`$.

We should note here that $`\mathrm{\Delta }_s`$ is absent on the right-hand side of (40), so the gap equation in this case (unlike in the BCS theory) is an explicit formula for the gap. The physical reason for this is that $`s(x)`$ is considered small, so we neglect the change of the bound-state wave functions due to its presence. The only effect of $`s(x)`$ we are taking into account is the change of the occupancy of the zero-energy states, which, by (33), depends only on the sign of $`s`$, not on its detailed shape. This is why $`s(x)`$ does not feed back into the right-hand side of (40).

To estimate the decay of $`\mathrm{\Delta }_s`$ into the bulk, we shall assume $`\mathrm{\Delta }_d`$ to be constant in space and with the angular dependence

$$\mathrm{\Delta }_{d,\theta }(𝐫)=\mathrm{\Delta }_o\mathrm{sin}2\theta ,$$
(41)
which should hold for

$$x>\xi \frac{\mathrm{}v_F}{\mathrm{\Delta }_o}.$$

Also, we shall assume a spherical (circular) Fermi surface,

$$dk_F(\theta )=k_Fd\theta .$$

Then the wave function of the $`y`$-moving bound states, including the normalization, will be

$$\left(\begin{array}{c}f(\theta ,\rho )\\ g(\theta ,\rho )\end{array}\right)=\sqrt{\frac{|\mathrm{sin}2\theta |}{2\xi }}\left(\begin{array}{c}1\\ i\end{array}\right)e^{\frac{|\rho \mathrm{sin}2\theta |}{\xi }},$$
(42)
so

$`s(x)`$ $`=`$ $`{\displaystyle \frac{k_F|V_s|}{\xi }}{\displaystyle \underset{\pi /2}{\overset{0}{}}}{\displaystyle \frac{d\theta }{2\pi }}2\times {\displaystyle \frac{|\mathrm{sin}2\theta |}{2}}e^{\frac{2}{\xi }\left|\frac{x}{\mathrm{cos}\theta }\mathrm{sin}2\theta \right|}=`$ (43)
$`=`$ $`{\displaystyle \frac{k_F|V_s|}{\pi \xi }}{\displaystyle \underset{0}{\overset{\frac{\pi }{2}}{}}}𝑑\theta \mathrm{sin}\theta \mathrm{cos}\theta e^{4\mathrm{sin}\theta \frac{x}{\xi }}.`$ (44)
We can do the integral by substitution $`\mathrm{sin}\theta =t`$, which gives

$$s(x)=\frac{k_F|V_s|}{\pi \xi }(t\frac{\xi }{4x}\left(\frac{\xi }{4x}\right)^2)e^{4t\frac{x}{\xi }}]_{t=0}^1.$$
(45)
We can neglect the contribution from the upper limit because it is effectively non-zero only for $`x<\xi /4`$, where our assumption of constant $`\mathrm{\Delta }_d`$ does not hold. The lower limit should have been at $`T_s/T_d`$, rather than at 0, to exclude the trajectories close to the nodes where the first-order perturbation theory breaks down. That cuts off the lower-bound contribution at $`x\frac{T_d}{4T_s}\xi 100`$Å, beyond which we would need a more refined theory for the behavior of the quasiparticles around the nodes. For $`x`$ much smaller than this distance, we can neglect the first term on the right hand side of (45), and replace the exponential by 1. We conclude, therefore, that

$$s(x)\frac{k_F|V_s|\xi }{16\pi x^2}$$
(46)
for

$$\xi <x<<\frac{T_d}{T_s}\xi .$$

We see that $`s=k_F|V_s|/\xi `$ times a function that is of order unity for $`x<\xi `$, and decays fast for $`x>\xi `$, as expected.

## V Transition Temperature

So far, we have shown the instability $`dd+is`$ only at $`T=0`$. Just as in the BCS case, it remains to be demonstrated that the transition temperature $`T_s`$ is finite. We therefore must study the free energy of the system, which we obtain from (34) when we replace $`E_\theta [s]`$ by $`(E_\theta [s])`$, the free energy of a single fermion level (see (6)), that is,

$$F[s]=\underset{0}{\overset{\mathrm{}}{}}𝑑x\frac{s^2(x)}{|V_s|}+k_F\underset{\pi /2}{\overset{\pi /2}{}}\frac{d\theta }{2\pi }\mathrm{cos}\theta (T)\mathrm{ln}(1+e^{E_\theta [s]/T}).$$
(47)
Minimization of this functional will give an equation for $`s(x)`$ that again agrees with the contribution to the gap equation from the ZBSs. As we see from (47), the variational equation for $`s`$ will now be very non-linear; it will no longer be an explicit formula for $`s`$. The reason is that at finite temperatures, the occupancy of a given state depends on the value of its energy. Even in first-order perturbation theory, this value depends on the shape of $`s(x)`$, not just its sign, so $`s(x)`$ enters through the Fermi function into the right-hand side of the gap equation, making it non-linear and therefore difficult to solve.

We still can make an order-of-magnitude estimate of $`F`$ as follows

$`\mathrm{ln}(1+e^{E_\theta [s]/T})+\mathrm{ln}(1+e^{E_\theta [s]/T})`$ $`=`$ $`\mathrm{ln}[(1+e^{E_\theta [s]/T})(1+e^{E_\theta [s]/T})]`$
$`=`$ $`\mathrm{ln}(2+2\mathrm{cosh}{\displaystyle \frac{E_\theta [s]}{T}})`$
$``$ $`\mathrm{ln}4+{\displaystyle \frac{1}{4}}\left({\displaystyle \frac{E_\theta [s]}{T}}\right)^2+O(E_\theta [s]^4)`$
$``$ $`\mathrm{ln}4+{\displaystyle \frac{1}{4}}{\displaystyle \frac{s^2}{T^2}}+O(s^4),`$
so

$$F[s]F[0]s^2\left(\frac{\xi }{|V_s|}\frac{k_F}{T}\right)+O(s^4).$$
(48)
From (48) we see that the system is unstable to the transition to the $`d+is`$ state below the temperature $`T_sk_F|V_s|/\xi `$, which is therefore of the same order of magnitude as $`|\mathrm{\Delta }_s|_{T=0}`$. (see (36)).

Following , we can trade the coupling constant $`V_s`$ for the transition temperature, $`T_{cs}`$ of a BCS superconductor with this coupling, $`T_{cs}e^{1/|V_s|}.`$ Then

$$T_s\frac{1}{\mathrm{ln}T_{cs}}.$$
(49)
Hence, $`T_s`$ increases sharply close to $`T_{cs}=0`$, which is consistent with the numerical results .

## VI Current

To study the surface current in the $`d+is`$ state, we shall go back to $`T=0`$ for simplicity. We observe that the states on the $`+y`$-moving quasiclassical trajectory from Fig. 1 will, upon the transition into the $`d+is`$ state with $`s>0`$, feel the pairing potential shown in Fig. 6. As $`\rho `$ goes from $`\mathrm{}`$ to $`+\mathrm{}`$, the twist of the phase $`\phi `$ of the order parameter is clockwise (from $`\pi `$ to 0) for an $`+y`$-moving trajectory and counterclockwise (from 0 to $`\pi `$) for a $`y`$-moving one. In both cases, this implies current flowing in the $`y`$ direction. This agrees with our previous calculations that showed that the $`y`$\- ($`+y`$-)moving bound states will be (un)occupied if $`s>0`$.

The agreement is quantitative, as we can easily check. In the linearized Andreev formalism, the contribution to the current density from a given state is

$$j_n^{(1d)}(\theta ,\rho )=ev_F(|f_n(\theta ,\rho )|^2+|g_n(\theta ,\rho )|^2),$$
(50)
that is, the charge of the state times its (Fermi) velocity times the occupation of that state. In our case, all of the current is carried by the occupied bound states because the contributions from the remaining pairs of corresponding countermoving states cancel each other out (see the end of Subsection III B). To calculate the total current density in the $`y`$-direction, we again have to include the contribution from both the incoming and the outgoing part of each $`y`$-moving trajectory (see Fig. 5), and we have to project onto the $`y`$-direction

$$\left(j_{ZBS}(x)\right)_y=k_F\underset{\pi /2}{\overset{0}{}}\frac{d\theta }{2\pi }\mathrm{sin}\theta \left[j_{ZBS}^{(1d)}(\theta ,\frac{x}{\mathrm{cos}\theta })+j_{ZBS}^{(1d)}(\theta ,\frac{x}{\mathrm{cos}\theta })\right].$$
(51)
On the other hand, in terms of the one-dimensional density $`n^{(1d)}=k_F/\pi `$ and the phase of the order parameter $`\phi (\theta ,\rho )`$ along the trajectory,

$`j_{o.p.}^{(1d)}(\theta ,\rho )`$ $`=`$ $`e{\displaystyle \frac{1}{2m}}n^{(1d)}_\rho \phi (\theta ,\rho )=`$ (52)
$`=`$ $`{\displaystyle \frac{ev_F}{2\pi }}_\rho \phi (\theta ,\rho ),`$ (53)
since for a spherical Fermi surface $`v_F=k_F/m`$. The formula for the total surface-current density will be the same as (51), except that we now have to integrate over both $`+y`$\- and $`y`$-moving trajectories

$$\left(j_{o.p.}(x)\right)_y=k_F\underset{\pi /2}{\overset{\pi /2}{}}\frac{d\theta }{2\pi }\mathrm{sin}\theta \left[j_{o.p.}^{(1d)}(\theta ,\frac{x}{\mathrm{cos}\theta })+j_{o.p.}^{(1d)}(\theta ,\frac{x}{\mathrm{cos}\theta })\right].$$
(54)
We do not expect the current densities (51) and (54) to be the same at a given point because the formula (52) has corrections, which are higher-order derivatives of $`\phi `$. Those corrections will not, however, contribute to the total surface current

$$I_y=_0^{\mathrm{}}𝑑xj_y(x),$$
(55)
which should then come out the same in the two calculations. Indeed, the bound states give us

$`(I_{ZBS})_y`$ $`=`$ $`k_F{\displaystyle \underset{\pi /2}{\overset{0}{}}}{\displaystyle \frac{d\theta }{2\pi }}\mathrm{sin}\theta \mathrm{cos}\theta {\displaystyle \underset{\mathrm{}}{\overset{+\mathrm{}}{}}}𝑑\rho j_{ZBS}^{(1d)}(\theta ,\rho ),`$ (56)
$`=`$ $`{\displaystyle \frac{ev_Fk_F}{4\pi }},`$ (57)
since

$$\underset{\mathrm{}}{\overset{+\mathrm{}}{}}𝑑\rho j_{ZBS}^{(1d)}(\theta ,\rho )=ev_F$$

due to the normalization of the wave functions; the minus sign indicates that the current is flowing in the $`y`$ direction. The formula for the total current in terms of the order-parameter phase is

$$(I_{o.p.})_y=k_F\underset{\pi /2}{\overset{\pi /2}{}}\frac{d\theta }{2\pi }\mathrm{sin}\theta \mathrm{cos}\theta \underset{\mathrm{}}{\overset{+\mathrm{}}{}}𝑑\rho j_{o.p.}^{(1d)}(\theta ,\rho ).$$
(58)
Now

$`{\displaystyle \underset{\mathrm{}}{\overset{+\mathrm{}}{}}}𝑑\rho j_{o.p.}^{(1d)}(\theta ,\rho )`$ $`=`$ $`{\displaystyle \frac{ev_F}{\pi }}(\phi (\theta ,+\mathrm{})\phi (\theta ,\mathrm{}))=`$ (59)
$`=`$ $`{\displaystyle \frac{ev_F}{2}}\text{sgn}(\theta ),`$ (60)
so

$$(I_{o.p.})_y=\frac{ev_Fk_F}{4\pi }=(I_{ZBS})_y,$$
(61)
since $`\text{sgn}(\theta )\mathrm{sin}\theta `$ is an even function, so the factor 1/2 in (59) compensates for the doubling of the integration domain of $`\theta `$ in (54) compared to (51). To get an order-of-magnitude estimate, we put

$`e`$ $``$ $`10^{19}\text{C}`$
$`v_F`$ $``$ $`10^5\text{m/s}`$
$`k_F`$ $``$ $`10^{10}\text{m}^1`$
and get $`|I_y|10^5`$A per CuO plane. From the approximate form of the bound state wave functions introduced in the previous section, we can also estimate the spatial distribution of the current density

$`j_y(x)`$ $`=`$ $`4ev_Fk_F{\displaystyle \underset{\pi /2}{\overset{0}{}}}{\displaystyle \frac{d\theta }{2\pi }}\mathrm{sin}\theta \left|f(\theta ,{\displaystyle \frac{x}{\mathrm{cos}\theta }})\right|^2=`$ (62)
$`=`$ $`{\displaystyle \frac{4ev_Fk_F}{\xi }}{\displaystyle \underset{0}{\overset{\pi /2}{}}}{\displaystyle \frac{d\theta }{2\pi }}\mathrm{sin}^2\theta \mathrm{cos}\theta e^{4\frac{x}{\xi }\mathrm{sin}\theta }=`$ (63)
$`=`$ $`{\displaystyle \frac{2ev_Fk_F}{\pi \xi }}{\displaystyle \underset{0}{\overset{1}{}}}𝑑tt^2e^{\frac{4x}{\xi }t}.`$ (64)
By the same argument as presented in the last section, we find that for $`\xi <x<<\frac{T_d}{T_s}\xi ,`$

$$j_y(x)\frac{ev_Fk_F}{16\pi \xi }\left(\frac{\xi }{x}\right)^3.$$
(65)
The extra power of $`x`$ in the denominator in (65), compared to (46), comes from the directional sine in (51).

The surface current induces magnetic field, which will be screened by the diamagnetic current in the superconductor. The total current density therefore is

$$(j_{tot}(x))_y=(j_{ZBS}(x))_y+(j_{dm}(x))_y.$$
(66)
According to (65), the current is localized within the distance $`\xi `$ from the surface, which is much smaller than $`\lambda `$, the in-plane penetration depth, because the cuprates are strongly type-2 superconductors. Hence, the diamagnetic response does not resolve the internal structure of $`(j_{ZBS}(x))_y`$, and we can estimate

$$(j_{dm}(x))_y=(j_{dm}(0))_ye^{x/\lambda }.$$
(67)
The surface current will be screened completely, because the magnetic field it induces is smaller than $`B_{c1}^c`$, the lower critical field in the $`c`$-direction. Indeed, even if the current flows in the same direction along all the CuO planes (as shown in Fig. 7), the magnetic field at distances $`x>d_c`$ (the interplane spacing $`10`$Å) from the surface will be

$$B=\frac{2\pi }{c}\frac{I}{d_c}10^2\text{G},$$
(68)
which is smaller by an order of magnitude than $`B_{c1}^c`$ for YBCO (see ). The complete screening implies

$$(I_{tot})_y\underset{0}{\overset{\mathrm{}}{}}𝑑x(j_{tot}(x))_y=0,$$
(69)
which together with (61) gives

$$(j_{dm}(0))_y=\frac{ev_Fk_F}{4\pi \lambda }.$$
(70)
For $`\xi <x<<\frac{T_d}{T_s}\xi ,`$ we can approximate the exponential in (67) by 1, so

$$(j_{tot}(x))_y\frac{ev_Fk_F}{16\pi }\left[\frac{4}{\lambda }\frac{1}{\xi }\left(\frac{\xi }{x}\right)^3\right].$$
(71)
This changes sign at distance

$$x_o\xi \kappa ^{1/3},$$
(72)
where $`\kappa \lambda /\xi `$ is the Ginzburg-Landau parameter. Due to the one-third power, $`x_o\xi `$ for reasonable values of $`\kappa `$ (say, between 50 and 500). This is consistent with the numerical results .

Note that the two calculations of the surface current agree (see (61)) because the ZBSs moving in the direction of the current are shifted down in energy and thus occupied, whereas those moving against the current are shifted up and unoccupied. We wish to stress that this is exactly opposite to the sign of the Doppler shift: the states moving along the current would be Doppler-shifted up, whereas those moving against the current would be Doppler-shifted down.

This point is further supported by the analogy between the ZBSs and low-lying excitations in a core of an $`s`$-wave vortex. We consider an idealized case: $`\mathrm{\Delta }=0`$ inside the vortex (at distances from the center smaller than $`R`$), and $`|\mathrm{\Delta }|=const`$ outside with the phase winding counterclockwise once around. We look at a quasiclassical trajectory passing close to the center of the vortex. We denote the coordinate along the trajectory as $`\rho `$ again and the phase at the intersection point with the vortex edge as $`\phi _\pm `$, see Fig. 8. Then the energy of a low-lying excitation moving from $`\rho =\mathrm{}`$ to $`\rho =+\mathrm{}`$ on that trajectory is

$$E=\frac{v_F}{4R}[(\phi _+\phi _{})\pi ]_{mod2\pi }.$$
(73)
For a trajectory passing through the center, $`\phi _+\phi _{}=\pi `$, so $`E=0`$, and the low-lying excitation is a ZBS. If we now shift the trajectory slightly to the left as shown in Fig. 8, then $`\phi _+\phi _{}>\pi `$ and $`E>0`$. In a real vortex, $`\mathrm{\Delta }`$ would be non-zero even inside, and the phase of $`\mathrm{\Delta }`$ would wind clockwise as we go from $`\rho =\mathrm{}`$ to $`\rho =+\mathrm{}`$, so we are going against the current. Moreover, $`\phi (\rho =+\mathrm{})\phi (\rho =\mathrm{})=\pi \text{ mod }2\pi `$, so $`\mathrm{\Delta }`$ behaves the same way as in Fig. 6. Hence, the $`dd+is`$ transition is analogous to shifting the quasiclassical trajectory away from the vortex center. In both cases, the ZBS will have a positive energy if it is moving against the current and negative energy if it is moving in the direction of the current.

## VII Discussion

In order to understand the basic physical mechanism of the $`dd+is`$ transition, we considered the change of the free energy of a CuO half-plane when an $`s`$-wave component of $`\mathrm{\Delta }`$ appears close to the 110 surface. By the Hubbard-Stratonovich transformation and the mean-field approximation, we decomposed the total free energy into the contribution from the single-particle states, which is decreased by Im $`\mathrm{\Delta }_s`$, and the HS term, which increases quadratically with $`\mathrm{\Delta }_s`$. Using first-order perturbation theory, we saw that the system favors the $`d+is`$ state at $`T=0`$, and that the transition is driven by the zero-bias states. Based on this argument and on the separation of energy scales associated with the $`d`$\- and $`s`$-wave components of $`\mathrm{\Delta }`$, we then calculated $`\mathrm{\Delta }_s`$ at zero temperature and estimated the transition temperature. Finally, we discussed the surface current in the $`d+is`$ state. We saw that it is carried by the occupied ZBSs; these states are not Doppler-shifted by the current.

## Acknowledgements

I want to give special thanks to M. Stone for suggesting the problem and discussing it with me. I wish to thank A.J. Leggett for helping me figure out the diamagnetic response in Sec.VI. I have also benefited from discussions with I. Adagideli, H. Aubin, L.H. Greene, D.E. Pugel, R. Ramazashvili, S. Sachdev, M. Turlakov, and H. Westfahl. I am grateful to C. Elliott for proofreading the manuscript. The project was supported by the grant NSF-DMR-98-17941.

## A Derivation of the Quasiclassical Hamiltonian

We derive (14) from (10) by expressing the original field $`\psi _\sigma (𝐫)`$ in terms of the slowly varying field $`\psi _{\sigma \theta }(𝐫)`$, see (11). We first substitute into the kinetic-energy operator:

$`ϵ(i)\psi _\sigma (𝐫)`$ $`=`$ $`{\displaystyle \underset{F.S.}{}}{\displaystyle \frac{dk_F(\theta )}{2\pi }}e^{i𝐤_F(\theta )𝐫}ϵ(𝐤_F(\theta )i)\psi _{\sigma ,\theta }(𝐫)=`$ (A1)
$`=`$ $`{\displaystyle \underset{F.S.}{}}{\displaystyle \frac{dk_F(\theta )}{2\pi }}e^{i𝐤_F(\theta )𝐫}\left[ϵ(𝐤_F(\theta ))+_𝐤ϵ(𝐤)|_{𝐤_F}(\theta )(i)+\mathrm{}\right]\psi _{\sigma ,\theta }(𝐫)`$ (A2)
$``$ $`{\displaystyle \underset{F.S.}{}}{\displaystyle \frac{dk_F(\theta )}{2\pi }}e^{i𝐤_F(\theta )𝐫}𝐯_F(\theta )(i)\psi _{\sigma ,\theta }(𝐫),`$ (A3)
since

$$ϵ(𝐤_F)=0\text{, and}_𝐤ϵ(𝐤)|_{𝐤_F}=𝐯_F.$$

The kinetic energy therefore is

$`{\displaystyle d^2r\underset{\sigma =,}{}\psi _\sigma ^{}(𝐫)ϵ(i)\psi _\sigma (𝐫)}`$ $`=`$ $`{\displaystyle }d^2r{\displaystyle \underset{\sigma }{}}{\displaystyle \underset{F.S.}{}}{\displaystyle \frac{dk_F(\theta )}{2\pi }}{\displaystyle \underset{F.S.}{}}{\displaystyle \frac{dk_F(\theta ^{})}{2\pi }}e^{i(𝐤_F(\theta )𝐤_F(\theta ^{}))𝐫}\times `$ (A4)
$`\times `$ $`\psi _{\sigma ,\theta ^{}}^{}(𝐫)𝐯_F(\theta )(i)\psi _{\sigma ,\theta }(𝐫).`$ (A5)
Now $`e^{i(𝐤_F(\theta )𝐤_F(\theta ^{}))𝐫}`$ oscillates with a wavelength much shorter than the lengthscale on which $`\psi _{\sigma ,\theta }(𝐫)`$ changes. Thus, the integral will be zero unless $`𝐤_F(\theta )𝐤_F(\theta ^{})=0`$, so we can effectively drop one integration over the Fermi surface, and obtain the kinetic energy of the form

$$d^2r\underset{\sigma }{}\underset{F.S.}{}\frac{dk_F(\theta )}{2\pi }\psi _{\sigma ,\theta }^{}(𝐫)𝐯_F(\theta )(i)\psi _{\sigma ,\theta }(𝐫).$$
(A6)
The potential energy is now given by

$`{\displaystyle }d^2rd^2r^{}{\displaystyle \frac{dk_F(\theta _1)}{2\pi }}{\displaystyle \frac{dk_F(\theta _2)}{2\pi }}{\displaystyle \frac{dk_F(\theta _3)}{2\pi }}{\displaystyle \frac{dk_F(\theta _4)}{2\pi }}V(𝐫𝐫^{})\times `$ (A7)
$`\times \mathrm{exp}[i(𝐤_F(\theta _1)𝐫𝐤_F(\theta _2)𝐫^{}+𝐤_F(\theta _3)𝐫^{}+𝐤_F(\theta _4)𝐫)]\times \psi _{\theta _1}^{}(𝐫)\psi _{\theta _2}^{}(𝐫^{})\psi _{\theta _3}(𝐫^{})\psi _{\theta _4}(𝐫)`$ (A8)
$`{\displaystyle }d^2rd^2r^{}{\displaystyle \frac{dk_F(\theta _1)}{2\pi }}{\displaystyle \frac{dk_F(\theta _2)}{2\pi }}{\displaystyle \frac{dk_F(\theta _3)}{2\pi }}{\displaystyle \frac{dk_F(\theta _4)}{2\pi }}V(𝐫𝐫^{})e^{i(𝐤_F(\theta _2)𝐤_F(\theta _3))(𝐫𝐫^{})}\times `$ (A9)
$`\times \psi _{\theta _1}^{}(𝐫)\psi _{\theta _2}^{}(𝐫)\psi _{\theta _3}(𝐫)\psi _{\theta _4}(𝐫)\times \mathrm{exp}[i(𝐤_F(\theta _1)𝐤_F(\theta _2)+𝐤_F(\theta _3)+𝐤_F(\theta _4))𝐫]`$ (A10)
$`={\displaystyle }d^2rd^2r^{}{\displaystyle \frac{dk_F(\theta _1)}{2\pi }}{\displaystyle \frac{dk_F(\theta _2)}{2\pi }}{\displaystyle \frac{dk_F(\theta _3)}{2\pi }}{\displaystyle \frac{dk_F(\theta _4)}{2\pi }}V(\theta _2,\theta _3)\times `$ (A11)
$`\times \psi _{\theta _1}^{}(𝐫)\psi _{\theta _2}^{}(𝐫)\psi _{\theta _3}(𝐫)\psi _{\theta _4}(𝐫)\times \mathrm{exp}[i(𝐤_F(\theta _1)𝐤_F(\theta _2)+𝐤_F(\theta _3)+𝐤_F(\theta _4))𝐫],`$ (A12)
where we used the assumption that $`V`$ changes on a much shorter length scale than $`\psi _{\sigma ,\theta }(𝐫)`$, performed the $`r^{}`$-integration, and introduced the Fourier transform

$$V(\theta ,\theta ^{})d^2re^{i(𝐤_F(\theta )𝐤_F(\theta ^{}))𝐫}V(𝐫).$$

Since $`V(𝐫)=V(𝐫)`$, we see that $`V(\theta ,\theta ^{})=V(\theta ^{},\theta )`$. Again, the integral vanishes unless the sum of the four momenta is zero. Out of the various ways that this may happen, we pick only the one that contributes to the singlet pairing, namely

$`𝐤_F(\theta _1)+𝐤_F(\theta _2)`$ $`=`$ $`𝐤_F(\theta _3)+𝐤_F(\theta _4)=0,i.e.`$ (A13)
$`\theta _1+\theta _2`$ $`=`$ $`\theta _3+\theta _4=0\text{mod}2\pi .`$ (A14)
Since we now have two constraints, we can drop two Fermi-surface integrations, and obtain the potential energy of the form

$$d^2r\frac{dk_F(\theta )}{2\pi }\frac{dk_F(\theta ^{})}{2\pi }V(\theta ,\theta ^{})\psi _\theta ^{}(𝐫)\psi _\theta ^{}(𝐫)\psi _\theta ^{}(𝐫)\psi _\theta ^{}(𝐫).$$

The total Hamiltonian in the quasiclassical approximation then is

$`H`$ $`=`$ $`{\displaystyle }d^2r[{\displaystyle \underset{\sigma }{}}{\displaystyle \underset{F.S.}{}}{\displaystyle \frac{dk_F(\theta )}{2\pi }}\psi _{\sigma ,\theta }^{}(𝐫)𝐯_F(\theta )(i)\psi _{\sigma ,\theta }(𝐫)+`$ (A15)
$`+`$ $`{\displaystyle \underset{F.S.}{}}{\displaystyle \frac{dk_F(\theta )}{2\pi }}{\displaystyle \frac{dk_F(\theta ^{})}{2\pi }}V(\theta ,\theta ^{})\psi _\theta ^{}(𝐫)\psi _\theta ^{}(𝐫)\psi _\theta ^{}(𝐫)\psi _\theta ^{}(𝐫)].`$ (A16)
The derivation of (A15) from (10) is far from rigorous. We could give a somewhat better, although much longer, argument. We believe the approximations used here are equivalent to the approximation in the Eilenberger formalism, because the Eilenberger equations can now be rigorously (apart from the mean-field approximation) derived from (A15).

## B Boundary Conditions at the Surface

We show the effect of the boundary on the Andreev spectrum. In general, the boundary will cause mixing of different $`\theta `$’s. For each $`\theta `$, though, we have a different Andreev equation, so adding the solutions of (21) for different $`\theta `$’s does not make sense. However, the Andreev wave functions describe only the slow variation of our excitations (changes on the length-scale $`\xi `$). The full wave functions containing the rapid oscillations as well are

$$\left(\begin{array}{c}f_{\theta ,n}(𝐫)\\ g_{\theta ,n}(𝐫)\end{array}\right)e^{i𝐤_F(\theta )𝐫},$$
(B1)
and these describe the single-particle excitations of the same Hamiltonian (10) (in the mean-field approximation), so those can be added. If we assume a specularly reflecting boundary, then the wave function will contain only two terms:

$$\left(\begin{array}{c}f_{\theta _{in},n}(𝐫)\\ g_{\theta _{in},n}(𝐫)\end{array}\right)e^{i𝐤_F(\theta _{in})𝐫}+\left(\begin{array}{c}f_{\theta _{out},n}(𝐫)\\ g_{\theta _{out},n}(𝐫)\end{array}\right)e^{i𝐤_F(\theta _{out})𝐫},$$

such that

$$\theta _{in}+\theta _{out}=\pi \text{mod}2\pi ,$$
(B2)
since the angles are measured from the positive-$`x`$ semiaxis, see Fig. 1a. The Dirichlet boundary condition gives

$$\left(\begin{array}{c}f_{\theta _{out},n}(𝐫)\\ g_{\theta _{out},n}(𝐫)\end{array}\right)=\left(\begin{array}{c}f_{\theta _{in},n}(𝐫)\\ g_{\theta _{in},n}(𝐫)\end{array}\right)\left(e^{i(𝐤_F(\theta _{in})𝐤_F(\theta _{out}))𝐫}\right)_{𝐫ϵsurface},$$

whereas the Neumann boundary condition gives

$$\left(\begin{array}{c}f_{\theta _{out},n}(𝐫)\\ g_{\theta _{out},n}(𝐫)\end{array}\right)=\left(\begin{array}{c}f_{\theta _{in},n}(𝐫)\\ g_{\theta _{in},n}(𝐫)\end{array}\right)\left(e^{i(𝐤_F(\theta _{in})𝐤_F(\theta _{out}))𝐫}\right)_{𝐫ϵsurface}$$

since

$$𝐧(𝐤_F(\theta _{in})𝐤_F(\theta _{out}))=0$$

for $`𝐧`$ perpendicular to the surface, and we used

$$k_F\left(\begin{array}{c}f\\ g\end{array}\right)>>𝐧_r\left(\begin{array}{c}f\\ g\end{array}\right)$$

to neglect the gradient of the Andreev wave function. Since (21) is linear, we can drop multiplicative constants and simply assume

$$\left(\begin{array}{c}f_{\theta _{out},n}(𝐫)\\ g_{\theta _{out},n}(𝐫)\end{array}\right)=\left(\begin{array}{c}f_{\theta _{in},n}(𝐫)\\ g_{\theta _{in},n}(𝐫)\end{array}\right)$$
(B3)
at the surface for either choice of the boundary condition. As $`\theta _{in}`$ is uniquely determined by $`\theta _{out}`$ through the relation (B2), we shall label the potential $`\mathrm{\Delta }`$ along the trajectory as well as the solutions of the corresponding Andreev equation by $`\theta _{out}`$. We shall drop the subscript “$`out`$” everywhere except in Appendix C, where we will need to distinguish $`\theta _{out}`$, the label for a trajectory as in Fig. 1, from $`\theta `$, the label for a position on the Fermi surface as in Fig. 2.

## C Gap Equation

We obtain the gap equation by substituting for $`\varphi _\theta (𝐫)`$ in (20) its mean-field value, that is, the pairing amplitude $`\varphi _\theta (𝐫)\psi _\theta (𝐫)\psi _\theta (𝐫)`$. To calculate this amplitude, we expand the field operators into energy eigenstates

$$\left(\begin{array}{c}\psi _\theta (𝐫)\\ \psi _\theta ^{}(𝐫)\end{array}\right)=\underset{n}{}\gamma _{\theta ,n}\left(\begin{array}{c}f_n(\theta ,\rho )\\ g_n(\theta ,\rho )\end{array}\right).$$
(C1)
Equation (C1) gives at $`T=0`$

$`\psi _\theta (𝐫)\psi _\theta (𝐫)`$ $`=`$ $`{\displaystyle \underset{n,n^{}}{}}f_n(\theta ,\rho )g_n^{}^{}(\theta ,\rho )\gamma _{\theta ,n^{}}^{}\gamma _{\theta ,n}=`$ (C2)
$`=`$ $`{\displaystyle \underset{n}{}}\mathrm{\Theta }(E_{\theta ,n})f_n(\theta ,\rho )g_n^{}(\theta ,\rho ).`$ (C3)
Note that in (C1), we explicitly sum over both positive and negative energies, unlike the Bogoliubov-de Gennes (BdG) formalism where we can sum over positive energies only, using the fact that

$$\left(\begin{array}{c}u_n(𝐫)\\ v_n(𝐫)\end{array}\right)\text{and}\left(\begin{array}{c}v_n^{}(𝐫)\\ u_n^{}(𝐫)\end{array}\right)$$
(C4)
are both solutions of the BdG equations with energies equal in absolute value and opposite in sign. However, this symmetry is lost here because the Andreev wave functions corresponding to the BdG wave functions (C4) live on different quasiclassical trajectories.

Close to the surface, we have to remember again that each line contributes to the pairing amplitude at a given point for two directions $`\theta `$; see Fig. 5. We will, therefore, have to distinguish between the label of the trajectory $`\theta _{out}`$ and the label for the pairing amplitude $`\theta `$. Specifically,

$`\theta _{out}`$ $`=`$ $`\theta \text{for}\theta ϵ({\displaystyle \frac{\pi }{2}},0)\text{and}`$ (C5)
$`\theta _{out}`$ $`=`$ $`\pi \theta \text{for}\theta ϵ(\pi ,{\displaystyle \frac{\pi }{2}}),`$ (C6)
so the contribution to the pairing amplitude from the $`y`$-moving bound states will be

$`\psi _\theta (𝐫)\psi _\theta (𝐫)_{ZBS}`$ $`=`$ $`f(\theta _{out},{\displaystyle \frac{x}{\mathrm{cos}\theta _{out}}})g^{}(\theta _{out},{\displaystyle \frac{x}{\mathrm{cos}\theta _{out}}})=`$ (C7)
$`=`$ $`f(\theta ,{\displaystyle \frac{x}{\mathrm{cos}\theta }})g^{}(\theta ,{\displaystyle \frac{x}{\mathrm{cos}\theta }})\text{for}\theta ϵ({\displaystyle \frac{\pi }{2}},0)`$ (C8)
and

$`\psi _\theta (𝐫)\psi _\theta (𝐫)_{ZBS}`$ $`=`$ $`f(\theta _{out},{\displaystyle \frac{x}{\mathrm{cos}\theta _{out}}})g^{}(\theta _{out},{\displaystyle \frac{x}{\mathrm{cos}\theta _{out}}})=`$ (C9)
$`=`$ $`f(\pi \theta ,{\displaystyle \frac{x}{\mathrm{cos}\theta }})g^{}(\pi \theta ,{\displaystyle \frac{x}{\mathrm{cos}\theta }})\text{for}\theta ϵ(\pi ,{\displaystyle \frac{\pi }{2}}).`$ (C10)
Substituting (C7) and (C9) into (20) gives

$`\mathrm{\Delta }_\theta (x)_{ZBS}`$ $`=`$ $`{\displaystyle \underset{\theta ^{}ϵ(\pi ,0)}{}}{\displaystyle \frac{dk_F(\theta ^{})}{2\pi }}V(\theta ,\theta ^{})\psi _\theta ^{}(𝐫)\psi _\theta ^{}(𝐫)_{ZBS}=`$ (C11)
$`=`$ $`{\displaystyle \underset{\theta _{out}^{}ϵ(\frac{\pi }{2},0)}{}}{\displaystyle \frac{dk_F(\theta _{out}^{})}{2\pi }}[V(\theta ,\theta _{out}^{})(i)|f(\theta _{out}^{},{\displaystyle \frac{x}{\mathrm{cos}\theta _{out}^{}}})|^2+`$ (C12)
$`+`$ $`V(\theta ,\pi \theta _{out}^{})(i)|f(\theta _{out}^{},{\displaystyle \frac{x}{\mathrm{cos}\theta _{out}^{}}})|^2],`$ (C13)
where

$$V(\theta ,\theta _{out}^{})=|V_s|+V_d(\theta ,\theta _{out}^{}),$$

and we used (32). The contribution to the $`s`$-wave component of the pairing potential from the occupied bound states therefore is

$$\mathrm{\Delta }_s(x)_{ZBS}=i|V_s|\underset{\theta _{out}^{}ϵ(\frac{\pi }{2},0)}{}\frac{dk_F(\theta _{out}^{})}{2\pi }\left[|f(\theta _{out}^{},\frac{x}{\mathrm{cos}\theta _{out}^{}})|^2+|f(\theta _{out}^{},\frac{x}{\mathrm{cos}\theta _{out}^{}})|^2\right]$$
(C14)
For the calculation of the $`d`$-wave component of $`\mathrm{\Delta }_{ZBS}`$, we will assume that the unperturbed $`\mathrm{\Delta }_d`$ is antisymmetric around its vertical node,

$$\mathrm{\Delta }_{d,\theta }(𝐫)=\mathrm{\Delta }_{d,\pi \theta }(𝐫).$$
(C15)
Presumably, $`\mathrm{\Delta }_d`$ arises from an antisymmetric interaction

$$V_d(\theta ,\theta ^{})=V_d(\theta ,\pi \theta ^{}).$$
(C16)
Along the quasiclassical trajectory, (C15) means

$$\mathrm{\Delta }_d(\theta ,\rho )=\mathrm{\Delta }_d(\theta ,\rho ),$$
(C17)
which, by (22), implies

$$|f(\theta ,\rho )|^2=|f(\theta ,\rho )|^2.$$
(C18)
Thus, under these assumptions,

$`\mathrm{\Delta }_{d,\theta }(x)_{ZBS}`$ $`=`$ $`(i){\displaystyle \underset{\theta _{out}^{}ϵ(\frac{\pi }{2},0)}{}}{\displaystyle \frac{dk_F(\theta _{out}^{})}{2\pi }}\left[V_d(\theta ,\theta _{out}^{})+V_d(\theta ,\pi \theta _{out}^{})\right]|f(\theta _{out}^{},{\displaystyle \frac{x}{\mathrm{cos}\theta _{out}^{}}})|^2=`$ (C19)
$`=`$ $`0.`$ (C20)