Boundary states with elevated critical temperatures in Bardeen-Cooper-Schrieffer superconductors

Bardeen-Cooper-Schrieffer (BCS) theory describes a superconducting transition as a single critical point where the gap function or, equivalently, the order parameter vanishes uniformly in the entire system. We demonstrate that in superconductors described by standard BCS models, the superconducting gap survives near the sample boundaries at higher temperatures than superconductivity in the bulk. Therefore, conventional superconductors have multiple critical points associated with separate phase transitions at the boundary and in the bulk. We show this by revising the Caroli-De Gennes-Matricon theory of a superconductor-vacuum boundary and finding inhomogeneous solutions of the BCS gap equation near the boundary, which asymptotically decay in the bulk. This is demonstrated for a BCS model of almost free fermions and for lattice fermions in a tight-binding approximation. The analytical results are confirmed by numerical solutions of the microscopic model. The existence of these boundary states can manifest itself as discrepancies between the critical temperatures observed in calorimetry and transport probes.

Bardeen-Cooper-Schrieffer theory [1] describes superconducting transition as a single critical point where the gap function or, equivalently the order parameter vanishes uniformly in the entire system. Behavior of the gap function near a boundary of superconductor and vacuum is one of the most basic questions in superconductivity which is widely believed to be fully understood since the works of C. Caroli, P.G. De Gennes, J. Matricon (CdGM) [2][3][4][5]. Using linear gap equation they derived famous boundary conditions for order parameter [2,3], which indicate that the gap at the boundary should vanish at the same temperature as the gap in the bulk. For the simplest vacuum-superconductor surface CdGM boundary condition reads n · ∇∆ = 0 (1) where n is normal to boundary and ∆ is order parameter. This boundary condition is a cornerstone concept for interpretation of standard surface probes of a superconductor such as Scanning Tunnelling Microscopy.
While standard textbook picture may leave an impression that it is a extremely well theoretically and experimentally understood physics, we show below that the picture is not so simple. First, as CdGM acknowledge themselves, [2,3] this condition (1) breaks down in a thin region near the boundary. Moreover the study was carried out only under the assumption of superconducting bulk. On the other hand, the same critical temperature both for bulk and for the surfaces are not a priori guaranteed.
Below, we carry out more general analysis than CdGM and show that near surface there is an upshoot for order parameter. Then it forms surface state with a tail going inside the bulk. In this case CdGM boundary condition are violated and superconducting state gets multiple critical temperatures: while the bulk goes normal a surface of a sample retains superconductivity. We show that this is true even for the basic BCS superconducting states with superconductor-vacuum boundary, where upshoot is a consequence of Friedel oscillations of particle density near surface.
The effect is consistent with several experimentally observed facts. Firstly, it is a commonly occurring experimental situation where weak superconducting response occurs slightly before specific heat jump. Usually it is attributed to the effect of inhomogeneities giving some areas of a sample higher T c , but is also consistent with our picture of two critical temperatures. Secondly, possible superconductivity of surfaces with normal bulk was previously invoked in the interpretation of increase of critical temperatures of granular elemental superconductors [6][7][8]. Since the effect is observed even for relatively large and relatively well connected granules, the V.L. Ginzburg proposal of the enhancement of superconductivity near surfaces was often invoked, [9,10]. However, the Ginzburg proposal was based on a conjecture that a surface may modify phonons, thus near a surface one may get another effective Hamiltionan with larger coupling constants. By contrast we do not make this modification of the Hamiltonian. We demonstrate below that more robust superconductivity appears on a surface of a BCS superconductor, quite generically without new pairing glue. Another example of a situation where superconductivity survives at elevated temperatures only on the surface while bulk is normal was found recently for pair-density-wave state [11]. However the mechanism discussed in [11] doesn't hold for conventional superconductors.
Lets us first recap canonical way of obtaining critical temperature for a bulk of superconductor. One can start, for example, from Bogoluibov-de Gennes (BdG) [2][3][4][5], Gorkov [12] or directly from path integral formulation of the model. Then assuming that transition is of second order near T c , order parameter should be small. Hence self consistency gap equation can be expanded in powers of ∆, to first order yielding [2]: where ω = 2πT (ν + 1 2 ), ξ n = E n − µ, n F is Fermi-Dirac distribution and w n are orthonormal one electron wave functions with eigenvalues E n . arXiv:1904.10942v1 [cond-mat.supr-con] 24 Apr 2019 In this simple model in the bulk of a superconductor w n = ae iknx is a plane wave with E n = k 2 n 2m . Then canonical assumption is that the highest critical temperature corresponds to ∆ = const, which leads to the usual equation for T c1 : 1 = V n F n,n 1 L d , with an implicit assumption that spatial gradients in ∆ cost energy and thus cannot increase T c . Then the standard BCS procedure gives: for the bulk of the system (5) where ω D is Debye frequency, N F is density of states at Fermi level and constant c = 2e γ π 1. Note, that linear gap equation (2) corresponds to system of linear equations of the form M ij ∆ j = 0, where we for simplicity discretiseze ∆(r), K(r, r ) → ∆ i , K i,i and matrix M = δ i,i − K i,i . Hence for T > T c we obtain det M = 0, which gives solution ∆ = 0, i.e. normal state. Exactly at T c det M = 0, which allows for non zero order parameter. Hence configuration of the ∆ could be found very approximately for T T c , and the only reliable result from linear gap equation (2) is the value of T c , which strongly depends on the guess for order parameter. It is important to keep in mind that CdGM conclusion about boundary condition (1) was derived based solely on the linear gap equation [2,3].
Never the less lets consider CdGM estimate for derivative of order parameter presented in [3], see formulas III.22 -III.27. To calculate III.26 one needs to estimate drH(r) which is defined by: where K 0 (r, r ) is defined so that K 0 (r, r )dr = 1 and K 0 (r, r ) = K 0 (r , r). These properties are used by CdGM to go from III.24 to III.26. Also using the same properties it is possible to show that CdGM derivation of III.26 is not correct. Namely, lets compute the integral: drH(r) = dr∆(r) − drdr K 0 (r, r )∆(r ) = dr∆(r) − dr ∆(r ) = 0 d∆ dz . Which means that in fact it's not possible to make any conclusions about the derivative of order parameter. Now we demonstrate that there is the second critical tempeature T c2 in a semi-infinite superconductor due to existence of energetically preferred inhomogeneous solutions of the gap equation near a boundary. For simplicity let's consider one dimensional superconductor positioned at x ∈ [0, L], below we implicitly assume the limit L → ∞. Superconductor is usually modeled to have periodic potential inside and finite potential barrier outside. Then single electron wave functions are given by standing Bloch waves inside and exponentially suppressed outside. To demonstrate that surface states exist even in the most unfavorable for them approximation, we begin by considering an infinite potential barrier and free electrons approximation inside. I.e. single electron wave function is nullified at the boundary of superconductor. Hence single electron eigenstates (standing waves) and eigenvalues will be: If we relax the assumption ∆ = const, the equation (2) becomes challenging to tackle. One can assess the situation by integrating it over space. Hence we obtain: where we have used orthogonality of w n . Note that usual assumption of constant gap gives the following value of the factor I n = I (1) 1. Hence to get critical temperature T > T c1 one should search for an inhomogeneous solution near a surface for which I n > I (1) . To do so we make the following ansatz for order parameter: which is motivated by the fact that 1) ∆ → 0 in the bulk on the length scale of coherence length ξ ∼ 1/q, i.e. T > T c1 and bulk is normal; 2) according to (2), ∆ should be composed of rapidly oscillating standing waves so that ∆ ∝ w n (x)w m (x), which in the simplest case corresponds to ∆ ∝ sin 2 (k F x).
Note, that this type of configuration (9) is superficially similar to Tamm-Shockley one-electron surface states [13,14], which are given by rapidly oscillating standing wave exponentially decaying into the bulk. In contrast to Tamm-Shockley states our states are obtained for a superconducting order parameter. Another principal difference is the absence of an assumption of periodic potential in the bulk.
Gap exponentially upshoots near surface in strong contrast to the nearly constant gap assumption in the standard derivation. Whether this solution wins over the constant gap solution can be determined by calculating its critical temperature. To calculate I n factor we estimate k k F , which gives: where Q = q 2 8m and we used the fact that Q, ξ n µ. Hence second critical temperature is: F V for the surface of the system (11) Lets generalize this result to three dimensional superconductors. For T > T c1 we expect superconducting surfaces with normal bulk. Hence close to x = 0 boundary we can assume that ∆(x, y, z) ∼ w 2 n (y)w 2 n (z)∆(x), so that in y, z directions it is, as usually, composed of single electron wave functions, but in x direction has the form of (9). Then we obtain the same formula (11) for second critical temperature T c2 , where surface states disappear. For T > T c2 superconductor is normal in the bulk and surfaces but on the edges it will be still superconducting up to third critical temperature T c3 . To see that consider edge defined by superconductor present at x > 0 and y > 0. Then we can estimate order parameter as ∆(x, y, z) ∼ w 2 n (y)∆(x)∆(y). Hence the I factor is simply I (3) = I (2) 2 = 9 4 > I (2) . Similarly after that is expected transition to corner states with I (4) = I (2) 3 > I (3) . So in total for three dimensional rectangular superconductor, at the mean-field-level, one would expect sequence of phase transitions: superconducting bulk → superconducting surf aces → superconducting edges → superconducting vertexes given by critical temperatures: with Let us now demonstrate the existence of multiple critical temperatures numerically in a one dimensional tightbinding model. We will do it within Bogoliubov-de-Gennes approach. Using Wannier functions we obtain the following mean field tight binding Hamiltonian for s-wave superconductor [15]: where h ij = −t ij − µδ ij and hopping coefficient is non zero only for nearest neighbors t ij = t. For numerical simulations we chose sufficiently big system size to get rid of all finite size effects and fixed number of sites to N = 100. Coefficients were set to t = 1, µ = 1 2 and V = 2.
We iteratively diagonalize Hamiltonian with Bogoliubov transformations: And recompute ∆ i from self-consistency equation: Numerically after ∼ 10 3 iterations we obtain converged states depicted on Fig. 1. These states are global minima and have energy lower than normal state. Observe that for T < T c1 superconducting gap is uniform in the bulk but has an up-shoot on the surface. This upshoot is relaxed on the length scale of coherent length ξ into bulk. Note, that ξ is maximal at T T c1 in agreement with usual behavior of ξ in the bulk. For T > T c1 superconductor is normal in the bulk but still order parameter is non zero and exhibits superconducting surface states. Order parameter varies smoothly since it's defined only on sites, which corresponds to averaging or evaluating at picks of oscillations of (9). To highlight nature of the upshoot on the boundary lets make analytic approximation for T < T c1 . It is possible to diagonalize hopping term in the tight binding Hamiltonian (13) using standing waves by performing the following transformation [16,17]: with k n = π N +1 n. For T → 0 terms depending on ∆ are simultaneously diagonalised and we obtain Hamiltonian: where ξ n = 2t cos(k n ) + µ and Next ∆ i can be composed of standing waves and constant in the following expansion: Since superconductor is uniform in bulk for T → 0 gap is almost constant and hence∆ 0 ∆ n =0 and∆ n ∆ ∆ 0 √ N +1 , where ∆ is average order parameter (∆ in the bulk). Gap equation (15) gives for n = 0: and ∆ is defined through usual gap equation: Then using (20) to transform back to real space by means of (19) one obtains ∆ i configuration quite close to numerically obtained Fig. 1 for T = 0. Crucial point is that when performing summation (19) one basically represents almost constant ∆ i ∆ order parameter by some set of standing waves. When doing so inevitably oscillations appear on the boundaries, known as Wilbraham-Gibbs phenomenon [18,19]. Similarly electron density exhibits oscillations and upshoot near surface. Other way to put it is that Friedel oscillations of electron density are strongest closer to boundary and hence lead to the upshoot in the density of states and hence upshoot for gap, which then relaxes into bulk on coherence length ξ.
We demonstrated above the effect in the microscopic model. Next we address the question of what are the implications for the Ginzburg-Landau boundary conditions. The standard Ginzburg-Landau formalism does not have the effect of the upshoot of the order parameter near the boundary if one uses the boundary conditions according to C. Caroli, P.G. De Gennes, J. Matricon [2,3]. Namely for T > T c1 global minimum of canonical free energy would be ∆ = 0, i.e. normal state which is inconsistent with microscopic calculations presented above. Lets revise usual microscopic derivation of GL free energy and demonstrate how it fails close to the surfaces.
One can start from path integral formulation of weakly interacting electron gas. Then perform Hubbard-Stratonovich transformation to introduce bosonic field ∆ and integrate out fermions. Obtained model is treated classically. Next, near T c Hamiltonian is expanded in powers of small order parameter ∆. From now on we will use usual Ginzburg-Landau formalism notation for the Ψ ∝ ∆. At this stage Euler-Lagrange equations in lowest order would be of the form (2). Final step is the assumption that Ψ varies slowly in space and hence the following expansion is performed [20]: Hence, for example, Ψ(r 1 )Ψ * (r 2 ) term corresponds to −Ψ(r 1 )∇ 2 Ψ * (r 1 ). Generally for the simplest possible GL Hamiltonian we get: (23) where K 1 = 2a(1 − a)K, K 2 = a 2 K and K 3 = (1 − a) 2 K, with a being arbitrary parameter which sets the point of expansion in derivatives (a = 1 corresponds to (22), i.e. r 1 point). Microscopically it's possible to retain values only for K = K 1 + K 2 + K 3 , α and β.
Where β and K are positive, but α changes its sign at the bulk phase transition. Note, that it's not possible to set K 2 = K 3 = 0. The standard calculations do not explicitly discuss that K 2 and K 3 terms are integrated by parts to obtain K 1 term. The integration by parts is a truly minor aspect for a bulk, where we indeed obtain the usual term K|∇Ψ| 2 . However if a system has a boundary, the integration by parts of the terms K 2 , K 3 gives additional surface integral: where d S = ndS with n unit vector orthogonal to surface.
Moreover free energy (23) now becomes unbounded from below near surface. To see this for simplicity lets set a = 1 and hence K 1 = K 3 = 0, K 2 = K. Then ansatz Ψ = Ψ 0 e −kx results in total energy: (25) which is divergent for any finite non zero Ψ 0 since F → −∞ for k → ∞. Thefore this more careful derivation agrees with the microscopic theory that gives an upshoot that we see on the surface scale 1/k F . Except, in the Ginzburg-Landau approximation it naturally corresponds to a divergence, because the theory has no small length scale cutoff of the microscopic theory 1/k → 0. However regularizing the model by a finite cutoff such as the one present in a lattice Ginzburg-Landau model, minimization of the energy (23) numerically yields the same type of solutions as in microscopic model Fig. 1.
In the bulk usual Euler-Lagrange second order differential equations are valid: Hence it is still possible to use Ginzburg-Landau formalism to estimate size of the surface states derived above microscopically. First of all lets assume that there is some upshoot on the surface of length 1/k F → 0 and magnitude Ψ 0 . Then since canonical GL is still valid in the bulk, equation (26) should be satisfied. Hence to the first order in Ψ we get solution for half infinite superconductor positioned at x > 0: Which agrees with microscopic solution averaged over rapid oscillations of order k F q (9). Although we do recover the upshoot effect if the integration by parts in the microscopic Ginzburg-Landau expansion is treated slightly more carefully, the above arguments highlight the fact that the usual microscopic derivation of GL formalism based on "bulk" Green functions that disregards the integrated by parts terms is inconsistent and should be revisited for superconductors near surface. We therefore expect it to give additional surface integral or lower symmetry terms in bulk energy.
There is one analogy to an ostensibly unrelated phenomenon one can drive here. Consider, for example, configuration of order parameter for T /T c1 = 0.995, Fig. 1. One can imagine that plot of ∆ i is actually a vertical cross section of a tank filled with water with ∆ being surface of the water. Lets show that this is not a coincidence. Energy of a thin column of water of width dx and height ∆ is composed of surface tension energy σdl and gravitational energy ρg 2 ∆ 2 dx. Where σ is energy of unit surface, dl is length of the surface, ρ is density of the water and g is gravitational constant. Note, that we implicitly redefined ∆ → λ∆, where λ is some dimensional constant, so that for new [∆] = [x]. Next if ∆ varies slowly. Then total energy: which is Ginzburg-Landau model if we substitute σ → K and ρg → α.
Then the problem of the surface states in superconductor can be related to the problem of adhesion of water to the wall. Note, that contact angle θ < π/2 (defined by tan θ = − 1 ∆ ) reflects the fact that molecules of the water are attracted stronger to the molecules of the tank than to themselves. Similarly we can interpret surface states in superconductor as effective increase N V → IN V with I > 1 near boundaries.
According to Fig. 1 contact angle θ is rather small and we can even try to assume that boundary conditions are in fact opposite to CdGM (1): Which implies infinite derivative at boundary or θ = 0. Even though it is not possible to apply condition (29) in Ginzburg-Landau theory directly since it is expanded in powers of derivatives, interestingly in the case of water adhesion there is an exact solution with infinite derivative (zero contact angle) but finite height ∆ 0 at the boundary given by [21]: In conclusion, we showed that in the BCS model, superconductivity survives at the surfaces, edges and vertexes of a sample to a significantly higher temperature than in the bulk. In the introduction we mention some examples where experimental results led to conjectures of surface superconductivity of samples which are normal in the bulk. Similarly this physics should be relevant for interface superconductors. However there are several reasons that can obscure the observation of the effect. First we considered a clean surface only. The derivation of T c2 directly applies only to a perfect surface. For a rough surface the gaps will be higher on convex parts, while concave parts will represent weak links, so that a rough surface will not necessarily be able to have high critical current. Secondly, the difference between the mean-field values of T c1 and T c2 can be diminished by fluctuations. However if one heats the system through the bulk phase transition, phase fluctuations in the bulk can induce remanent vorticity on the surfaces, yielding a resistive state immediately above the bulk transition. While resolving spatial profile of the gap in the surface state might be experimentally difficult in supercondictors, the effect can be tested indirectly by a fabrication of structures with large surface to volume ratio. The solution for the surface gap and the sequence of the phase transitions can be directly studied in fermionic ultracold atoms in a box potential [22].