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

TL;DR

This paper reveals that in standard BCS superconductors, the superconducting gap persists at the boundaries at higher temperatures than in the bulk, indicating multiple critical points and boundary-specific phase transitions.

Contribution

The study demonstrates, both analytically and numerically, that boundary regions in BCS superconductors can remain superconducting at higher temperatures than the bulk, challenging the traditional single critical point paradigm.

Findings

Boundary superconductivity persists at higher temperatures than bulk.

Multiple critical points exist for boundary and bulk transitions.

Analytical and numerical solutions confirm boundary states in BCS models.

Abstract

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.

Figures

Equations

$$\vec{n}\cdot \nabla \Delta =0$$

$$\frac{1}{V}\Delta (\text{r})=\int d{\text{r}}^{\prime }K(\text{r},{\text{r}}^{\prime })\Delta ({\text{r}}^{\prime })$$

$$K(\text{r},{\text{r}}^{\prime })=\int d\text{k}d{\text{k}}^{\prime }{F}_{\text{k},{\text{k}}^{\prime }}{w}_{\text{k}}^{\ast }(\text{r})w_{{\text{k}}^{\prime }}^{\ast }(\text{r})w_{\text{k}}({\text{r}}^{\prime })w_{{\text{k}}^{\prime }}({\text{r}}^{\prime })$$

$$F_{\text{k},{\text{k}}^{\prime }}=T\sum _{\omega }\frac{1}{{\xi }_{\text{k}}-\mathrm{i}\omega }\frac{1}{{\xi }_{{\text{k}}^{\prime }}+\mathrm{i}\omega }=\frac{1-n_F({\xi }_{\text{k}})-n_F({\xi }_{{\text{k}}^{\prime }})}{{\xi }_{\text{k}}+{\xi }_{{\text{k}}^{\prime }}}$$

$$w_{\text{k}}=\frac{1}{\sqrt{2\pi }}{e}^{\mathrm{i}\text{k}\cdot \text{r}},E_{\text{k}}=\frac{{\text{k}}^2}{2m}$$

$$\frac{1}{V}{\Delta }_k=D_k{\Delta }_k$$

$$\hat{k}=\frac{k}{2k_F},\hat{x}=2k_{F}x,\hat{T}=\frac{T}{\mu },\hat{V}=\frac{V}{4\pi \mu }.$$

$${\hat{T}}_{c1}=c{e}^{-\frac{1}{\hat{V}}}\mathrm{for}\mathrm{the}\mathrm{bulk}\mathrm{of}\mathrm{the}system.$$

$$H(\text{r})\equiv \Delta (\text{r})-\int K^0(\text{r},{\text{r}}^{\prime })\Delta ({\text{r}}^{\prime })d{\text{r}}^{\prime }$$

$$w_k(x)=\{\begin{array}{ll}\sqrt{\frac{2}{\pi }}\sin kx,& \text{if}x\ge 0\\ 0& \text{if}x<0\end{array},k>0,E_k=\frac{k^2}{2m}.$$

$$\frac{1}{V}{\Delta }_k=D_k{\Delta }_k-{\int }_0^{\infty }{A}_{k,k^{\prime }}{\Delta }_{k^{\prime }}d{k}^{\prime }$$

$${\Delta }_k\propto {\delta }_{\xi }(k)\equiv \frac{2}{\pi }\frac{\xi }{1+{\xi }^2{k}^2}\text{for}k\ll \delta k_F$$

$$\frac{1}{\hat{V}}={\hat{D}}_{\mathrm{i}/\hat{\xi }}\simeq \ln \frac{c}{\hat{T}}+\frac{\alpha }{{\hat{T}}^2}\frac{1}{{\hat{\xi }}^2}$$

$$\frac{\alpha }{\hat{T}}={\int }_0^{\infty }{\hat{A}}^2(0,\hat{k})d\hat{k},\text{with}\alpha =\frac{7\zeta (3)}{{\pi }^2}\simeq \mathrm{0.85.}$$

$${\Delta }_{\hat{k}}\simeq {\hat{D}}_0{\delta }_{\hat{\xi }}(\hat{k})-{\hat{A}}_{\hat{k},0}.$$

$$\frac{1}{\hat{V}}-\ln \frac{c}{\hat{T}}=-{\int }_0^{\infty }{\hat{A}}_{0,\hat{k}}\frac{{\Delta }_{\hat{k}}}{{\Delta }_0}d\hat{k}\simeq \frac{\pi }{2\hat{\xi }}\left(\frac{\alpha }{\hat{T}\ln \frac{c}{\hat{T}}}-1\right).$$

$$\frac{1}{\hat{V}}-\ln \frac{c}{\hat{T}}=\frac{{\pi }^2}{4\alpha }{\left(\frac{\alpha }{\ln \frac{c}{\hat{T}}}-\hat{T}\right)}^2.$$

$$\delta \hat{T}\equiv \frac{{\hat{T}}_{c2}-{\hat{T}}_{c1}}{{\hat{T}}_{c1}}\simeq \exp \left(\frac{{\pi }^2}{4\alpha }{\left(\alpha \hat{V}-c{e}^{-\frac{1}{\hat{V}}}\right)}^2\right)-1.$$

$$\frac{1}{\hat{V}}=\frac{{\int }_0^{\infty }{\hat{D}}_{\hat{k}}|{\Delta }_{\hat{k}}{|}^{2}d\hat{k}-{\int }_0^{\infty }{\Delta }_{\hat{k}}{\hat{A}}_{\hat{k},{\hat{k}}^{\prime }}{\Delta }_{{\hat{k}}^{\prime }}d\hat{k}d{\hat{k}}^{\prime }}{{\int }_0^{\infty }|{\Delta }_{\hat{k}}{|}^{2}d\hat{k}}\le \frac{1}{{\hat{V}}_{\text{exact}}}.$$

$$H=\sum _{i,j=1,\sigma =\pm }^{i,j=N}{c}_{i\sigma }^{†}{h}_{ij}{c}_{j\sigma }+\sum _{i=1}^N\left({\Delta }_i{c}_{i\uparrow }^{†}{c}_{i\downarrow }^{†}+{\Delta }_i^{\ast }{c}_{i\downarrow }{c}_{i\uparrow }\right)$$

$$c_{i\sigma }=\sum _{n=1}^{2N}{u}_{i\sigma }^n{\gamma }_n-\sigma v_{i\sigma }^{n\ast }{\gamma }_n^{†},$$

$$c_{i\sigma }^{†}=\sum _{n=1}^{2N}{u}_{i\sigma }^{n\ast }{\gamma }_n^{†}-\sigma v_{i\sigma }^n{\gamma }_n$$

$${\Delta }_i=\frac{V}{2}\sum _{n=1}^{2N}{u}_{i\uparrow }^n{v}_{i\downarrow }^{n\ast }\tanh \left(\frac{E_n}{2k_{B}T}\right)$$

$$c_{i\sigma }=\sqrt{\frac{2}{N+1}}\sum _{n=1}^N\sin \left(k_{n}i\right)d_{n\sigma }$$

$$H\simeq -\sum _{n\sigma }{\xi }_n{d}_{n\sigma }^{†}{d}_{n\sigma }+\sum _n\left({\tilde{\Delta }}_n^{\prime }{d}_{n+}^{†}{d}_{n-}^{†}+({\tilde{\Delta }}_n^{\prime }{)}^{\ast }{d}_{n-}{d}_{n+}\right)$$

$$\frac{2}{N+1}\sum _{i=1}^N\sin (k_{n_1}i)\sin (k_{n_2}i){\Delta }_i\simeq {\delta }_{n_1,n_2}{\tilde{\Delta }}_{n_1}^{\prime }.$$

$${\Delta }_i=\frac{1}{\sqrt{N+1}}\sum _{n=0}^N\cos (2k_{n}i){\tilde{\Delta }}_n.$$

$${\tilde{\Delta }}_n=-\frac{V}{2}\frac{1}{\sqrt{N+1}}\frac{⟨\Delta ⟩}{\sqrt{{\xi }_n^2+⟨\Delta {⟩}^2}}\tanh \frac{\sqrt{{\xi }_n^2+⟨\Delta {⟩}^2}}{2T}$$

$$1=\frac{V}{2}\frac{1}{N+1}\sum _{n=1}^N\frac{\tanh \frac{\sqrt{{\xi }_n^2+⟨\Delta {⟩}^2}}{2T}}{\sqrt{{\xi }_n^2+⟨\Delta {⟩}^2}}$$

$${\rho }_i(E)=-\sum _n\left[|u_{i\uparrow }^n{|}^2{n}_F^{\prime }(E_n-E)+|v_{i\downarrow }^n{|}^2{n}_F^{\prime }(E_n+E)\right].$$

$$\text{non periodic:}{w}_n(i)=\sqrt{\frac{2}{N+1}}\sin (k_{n}i),E_n=2t\cos k_n$$

Full text

I Introduction

The behavior of the gap function near a boundary of a BCS superconductor Bardeen et al. (1957) is one of the most basic questions in superconductivity. The problem was widely believed to be solved after the works of Caroli, De Gennes and Matricon (CdGM) de Gennes (1964); Caroli et al. (1962, 1963); De Gennes and Guyon (1963). These works derived the famous boundary conditions for the order parameter de Gennes (1964); Caroli et al. (1962) using linearized gap equations. The conclusion was reached that the gap at the boundary should vanish at the same temperature as the gap in the bulk, implying that the presence of a boundary in a large system does not alter the conclusions of BCS theory regarding the critical temperature associated with the disappearance of resistivity. For the simplest vacuum-superconductor surface the CdGM boundary condition reads

[TABLE]

where $\vec{n}$ is the unit vector normal to the boundary and $\Delta$ is the order parameter. This boundary condition is a cornerstone concept for the interpretation of standard surface probes of a superconductor such as Scanning Tunneling Microscopy.

While the standard textbook picture may leave the impression that this is extremely well theoretically and experimentally understood physics, we show below that the actual structure of the superconductor-vacuum boundary is principally different. First, let us emphasize that, as CdGM acknowledge themselves, de Gennes (1964); Caroli et al. (1962) the condition (1) breaks down in a thin region near the boundary. Moreover, the studies in Refs. de Gennes (1964); Caroli et al. (1962) were carried out only under the assumption of a superconducting bulk. However, equality of the critical temperatures for the bulk and the boundary is not a priori guaranteed. Below, we carry out a more general analysis than CdGM and show that the conclusions in the works de Gennes (1964); Caroli et al. (1962) are not correct. Instead, even in the simplest BCS model, near a boundary there is a strong increase of the order parameter which in general takes place on a macroscopic length scale. Hence a superconductor forms a superconducting boundary state. The tail of this state extends inside the bulk over a large length scale which makes the CdGM boundary condition incorrect. The consequence of the gap behavior near the boundary that we find is that the superconducting state gets multiple critical temperatures: above the temperature at which the bulk of the system becomes normal, the boundary retains superconductivity.

This effect is consistent with several experimentally observed facts. First, it is a commonly known fact that resistivity, diamagnetic response and thermodynamic measurements do not in general give the same value for critical temperature $T_{c}$. By contrast, a commonly observed situation is that a weak superconducting response occurs at slightly higher temperature than a specific-heat jump. Usually this is attributed to the effect of inhomogeneities giving some areas of the sample higher $T_{c}$. However, the interpretation of these measurements in terms of possible boundary superconductivity in systems with normal bulk has also been voiced Lortz et al. (2006); Janod et al. (1993); Butera (1988). Secondly, possible superconductivity of boundaries in systems with normal bulk was previously invoked in connection with the increase of critical temperatures observed in granular elemental superconductors Deutscher et al. (1973); Cohen and Abeles (1968); Cohen et al. (1967). The basis of this interpretation is that the effect is observed even for relatively large and relatively well-connected granules. The theoretical mechanisms suggested for both cases are (i) the conjecture by V.L. Ginzburg Ginzburg (1964); Naugle et al. (1973) that a surface may modify phonons, and (ii) the conjecture of different chemical composition of the surface. In these scenarios the appearance of boundary states is attributed to a modification of the surface, and hence to the appearance of a different effective Hamiltonian with larger coupling constants near the boundaries, which to the best of our knowledge was not demonstrated at the ab-initio level. Another example of a situation in which superconductivity survives at elevated temperatures only on the boundary while the bulk is normal was found recently for a pair-density-wave state Barkman et al. (2019). However, the mechanism discussed in Ref. Barkman et al. (2019) doesn’t hold for conventional superconductors. We demonstrate below that more robust superconductivity appears on a clean surface of a BCS superconductor quite generically, and does not require the conjecture of increased coupling constants.

II Boundary superconductivity in the model of almost free fermions with attractive interaction

Let us first recap the standard calculation of the critical temperature for a superconductor. This can be done for example in Bogoluibov-de Gennes (BdG) formalism de Gennes (1964); De Gennes and Guyon (1963); Caroli et al. (1962, 1963), in Gorkov formalism Gorkov (1960) or directly from a path-integral formulation of the model. Since the transition is of second order, the order parameter should be small near $T_{c}$. Hence the self-consistency gap equation can be expanded in powers of $\Delta$. To first order this yields de Gennes (1964):

[TABLE]

[TABLE]

[TABLE]

where $\omega=2\pi T(\nu+\frac{1}{2})$, $\xi_{\textbf{k}}=E_{\textbf{k}}-\mu$, $\mu=\frac{k_{F}^{2}}{2m}$, $n_{F}$ is the Fermi-Dirac distribution and the $w_{\textbf{k}}$ are single-electron wave functions with eigenvalues $E_{\textbf{k}}$.

In the bulk of a superconductor, single-electron wave functions can be approximated as plane waves:

[TABLE]

with $\int_{-\infty}^{\infty}w_{\textbf{k}}(\textbf{r})w_{\textbf{k}}^{*}(\textbf{r}^{\prime})d\textbf{k}=\delta(\textbf{r}-\textbf{r}^{\prime})$. Then in one dimension, in momentum space the linear gap equation (2) for $\Delta_{k}=\int_{-\infty}^{\infty}w_{k}(x)\Delta(x)dx$ becomes:

[TABLE]

where $D_{k}=\frac{1}{2\pi}\int_{-\infty}^{\infty}F_{k-k^{\prime},k^{\prime}}dk^{\prime}$ is maximal at $k=0$. For this one-dimensional model we can introduce dimensionless quantities as

[TABLE]

Hence the only parameters will be $\hat{T}$ and $\hat{V}$. Then $\hat{D}_{0}\simeq\ln\frac{c}{\hat{T}}$ with $c=\frac{8e^{\gamma}}{\pi}$ for $\hat{T}\to 0$ (namely, $\hat{T}\lesssim 0.1$ or $\hat{V}\lesssim 0.2$). This gives the usual bulk critical temperature:

[TABLE]

Note that this means that $\Delta_{k}=\delta(k)$ and hence $\Delta(x)=\text{const}$.

The linear gap equation (2) corresponds to system of linear equations of the form $M_{ij}\Delta_{j}=0$, where we for simplicity discretize $\Delta(\textbf{r}),\ K(\textbf{r},\textbf{r}^{\prime})\to\Delta_{i},\ K_{i,i^{\prime}}$ and matrix $M=\delta_{i,i^{\prime}}-K_{i,i^{\prime}}$. Hence for $T>T_{c}$ we obtain $\det M\neq 0$, which gives the solution $\Delta=0$, i.e. the normal state. Exactly at $T_{c}$ $\det M=0$, which means that there exist a zero eigenvalue of $M$. The corresponding eigenvector gives the configuration of the order parameter up to a multiplicative constant. However, the calculation of $T_{c}$ is in general nontrivial. While the standard assumption of a homogeneous order-parameter configuration at first glance appears to be reasonable, we show below that it is not accurate. It is important to keep in mind that the conclusion of CdGM about the boundary condition (1) was derived solely on the basis of the linear gap equation, de Gennes (1964); Caroli et al. (1962) which does not allow one to assess the configuration of the order parameter at $T<T_{c}$.

Let us consider the CdGM estimate for the derivative of the order parameter. It is derived in Ref. Caroli et al. (1962), equations III.22 – III.27. To calculate III.26 one needs to estimate $\int d\textbf{r}H(\textbf{r})$, defined by:

[TABLE]

where $K^{0}(\textbf{r},\textbf{r}^{\prime})$ is defined so that $\int K^{0}(\textbf{r},\textbf{r}^{\prime})d\textbf{r}^{\prime}=1$ and $K^{0}(\textbf{r},\textbf{r}^{\prime})=K^{0}(\textbf{r}^{\prime},\textbf{r})$. These properties are used by CdGM to go from equation III.24 to equation III.26. By using the same properties, it is possible to show that the derivation by CdGM of the equation III.26 is not correct. Namely, let us compute the integral: $\int d\textbf{r}H(\textbf{r})=\int d\textbf{r}\Delta(\textbf{r})-\int d\textbf{r}d\textbf{r}^{\prime}K^{0}(\textbf{r},\textbf{r}^{\prime})\Delta(\textbf{r}^{\prime})=\int d\textbf{r}\Delta(\textbf{r})-\int d\textbf{r}^{\prime}\Delta(\textbf{r}^{\prime})=0\nsim\frac{d\Delta}{dz}$. This means that in fact it is not possible to obtain the derivative of the order parameter near the boundary, contrary to the conclusions of Ref. Caroli et al. (1962).

Now let us demonstrate that there is a second critical temperature $T_{c2}$ in a semi-infinite superconductor. This effect appears due to the existence of energetically preferred inhomogeneous solutions of the linear gap equation near the boundary. For simplicity let’s consider a one-dimensional semi-infinite superconductor positioned at $x>0$. The sample is usually modeled to have a periodic potential inside and a finite potential barrier outside. Then single-electron wave functions are given by standing Bloch waves inside and are exponentially suppressed outside. We begin by considering an infinite potential barrier and the almost-free-electron approximation. Thus the single-electron wave function is nullified outside the superconductor. Hence, near the boundary, a single-electron wave function is formed by incoming and reflected electrons, which give standing waves:

[TABLE]

We transform to momentum space with $\Delta_{k}=\int_{0}^{\infty}\phi_{k}(x)\Delta(x)dx$, where $\phi_{k}(x)=\sqrt{\frac{2}{\pi}}\cos kx$, which satisfies $\int_{0}^{\infty}\phi_{k}(x)\phi_{k}^{\prime}(x)=\delta(k-k^{\prime})+\delta(k+k^{\prime})$. Then the one-dimensional gap equation (2) becomes:

[TABLE]

where $A_{k,k^{\prime}}=\frac{1}{2\pi}F_{\frac{k+k^{\prime}}{2},\frac{k-k^{\prime}}{2}}$ and $D_{k}=\int_{0}^{\infty}A_{k,k^{\prime}}dk^{\prime}$ is the same as in the bulk equation (6). Equation (11) has to be solved, for a given interaction strength $\hat{V}$, in order to determine whether the boundary has different critical temperature $\hat{T}_{c2}$ than the bulk. Equivalently, one can fix $\hat{T}$ and look for the biggest eigenvalue, which gives $\frac{1}{\hat{V}}$.

We obtained solutions numerically by discretizing the integral in Eq. (11). The solutions show that there are superconducting boundary states. A typical example of the configuration of the order parameter in such a state is shown in Fig. 1.

Note that the type of configuration shown in Fig. 1b is superficially similar to non-superconducting Tamm-Shockley single-electron surface states Tamm (1932); Shockley (1939). These have the form of rapidly oscillating standing waves exponentially decaying into the bulk. In contrast to the 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.

Let us now consider how the boundary state in question is formed. We have three different scales at play here: (i) The shortest one is the Fermi scale $\frac{1}{2k_{F}}$, which sets the scale of the rapid oscillations of $\Delta(x)$ near the boundary, see Fig. 1b. (ii) There are oscillations with slightly different values of $k$. They dephase and disappear over the dephasing length scale $\lambda_{d}\propto\frac{1}{\delta k_{F}}$, where $\delta k_{F}$ is the width of the $k=2k_{F}$ peak in Fig. 1a. (iii) There is an increase of the order parameter near the boundary. It serves as a source of superconducting correlations in the bulk. It should asymptotically decay to zero over the *bulk coherence length * $\xi$.

Next, we consider the limit of $T_{c2}\to 0$, or equivalently $V\to 0$. For $x\gg\lambda_{d}$, the order parameter decays as $\Delta(x)\propto e^{-x/\xi}$. In momentum space this corresponds to

[TABLE]

where $\lim\limits_{\xi\to\infty}\int_{0}^{\infty}g(k)\delta_{\xi}(k)dk=g(0)$. Let us estimate $\xi$. Exponential decay in real space can be obtained by the substitution $k\to\mathrm{i}/\xi$. Then, neglecting the second term on the right-hand side of (11), we have:

[TABLE]

where

[TABLE]

We are looking for the biggest eigenvalue $\frac{1}{V}$ of the integral operator on the right hand side of (11). It is possible to find it by the power-iteration method. In other words, we can iteratively apply this operator, until convergence is obtained. In practice, however, starting from Eq. (12) even the first iteration gives a configuration of the order parameter which is in good agreement with the numerical solution:

[TABLE]

Now the boundary critical temperature $T_{c2}$ can be estimated by using the gap equation (11) for $k=0$:

[TABLE]

By substituting the coherence length from Eq. (13) we obtain:

[TABLE]

Then using the expression for the bulk critical temperature (8) we obtain the following expression for the relative boundary-bulk critical temperature difference for $\hat{T}_{ci}\to 0$:

[TABLE]

Note that the effect studied here is very easy to miss in an approximate calculation. For example, $\delta\hat{T}=0$, if one follows the common procedure and averages the oscillations over the Fermi scale $\frac{1}{k_{F}}$. This is because these oscillations are encoded in $\hat{A}_{\hat{k},0}$ in (15), which is essential for correct calculation of the difference in the critical temperatures.

To ascertain that we are not overestimating the critical temperature of the boundary states in (18) we can use the following inequality for the approximated value of $\frac{1}{\hat{V}}$:

[TABLE]

Then from Eqs. (19) and (6) we can obtain a lower bound on the critical temperature difference $\delta T$. Combining this result with numerical solutions we obtain the dependence of $\delta T$ on interaction strength $\hat{V}$ shown in Fig. 2.

The physical interpretation of the above derivation is the following: the Cooper pairs near a boundary are made of electrons in a superposition of incoming and reflected waves (i.e. a standing wave). This type of pairing wins over the conventional Cooper pairing of electrons with opposite momenta. This means that in order to observe an increase of the critical temperature one should study surfaces or other kinds of interfaces that efficiently reflect electrons, resulting in the formation of standing waves. An example is a superconductor-to-insulator boundary.

III Boundary superconductivity in the tight-binding model of fermions

In the previous section we demonstrated the existence of boundary states in the simplest model of almost free electrons with weak attractive interaction. In order to make a more general assessment of the presence or absence of this effect, in this section we consider the case of lattice fermions in a tight-binding model. In this section we utilize the Bogoliubov-de-Gennes approach. In the tight-binding limit the standard Hamiltonian for an s-wave superconductor reads Zhu (2016):

[TABLE]

where $h_{ij}=-t_{ij}-\mu\delta_{ij}$ and the hopping coefficient is nonzero only for nearest neighbors $t_{\langle ij\rangle}=t$.

First we will consider one-dimensional systems. We have to work with a sufficiently large system size (with number of sites up to $N=100$) to avoid finite-size effects. Coefficients were set to $t=1$, $\mu=\frac{1}{2}$ and $V=2$.

We iteratively diagonalize the Hamiltonian with Bogoliubov transformations:

[TABLE]

And recompute $\Delta_{i}$ from the self-consistency equation:

[TABLE]

After $\sim 10^{3}$ iterations we obtain converged states depicted in Fig. 3. These states are global minima and have energy lower than the normal state. Observe that for $T<T_{c1}$ the superconducting gap is uniform in the bulk but has an increase near the boundary. We find that this increase gradually vanishes in the bulk on the length scale of the bulk coherence length $\xi$.

The numerical calculations confirm that the largest extent of the boundary state takes places near the critical temperature associated with the bulk transition $T\simeq T_{c1}$. Indeed, at $T_{c1}$ the bulk coherence length $\xi$ diverges and the boundary state has a power-law decay into the bulk. For $T>T_{c1}$ a superconductor is normal in the bulk. However, in the tight-biding regime we similarly find that the gap remains non-zero near the boundaries above the temperature of the bulk transition. Note that in Fig. 3 the order-parameter field varies relatively smoothly since it’s defined only on sites, which corresponds to averaging or evaluating at the peaks of the oscillations of Fig. 1b. Note that the situation is similar to that for almost free electrons, considered in the previous section. Namely, the free boundary conditions here correspond to perfect reflection of the particles. This allows the formation of boundary states due to the increase of the density of states near the boundaries.

To understand the oscillation of the gap near the boundary, let us now make an analytic approximation for $T<T_{c1}$. It is possible to diagonalize the hopping term in the tight-binding Hamiltonian (20) using standing waves by performing the following transformation Jullien (1981); Busch and Penson (1987):

[TABLE]

where $k_{n}=\frac{\pi}{N+1}n$. In the limit $T\to 0$, the terms that depend on $\Delta$ are diagonalized and we obtain the Hamiltonian:

[TABLE]

where $\xi_{n}=2t\cos(k_{n})+\mu$ and

[TABLE]

Next the solution for $\Delta_{i}$ can be composed of standing waves and a constant as follows:

[TABLE]

Since a superconductor has a uniform gap in its bulk, for $T\to 0$ the gap is almost constant except a small region near the boundary and hence $\tilde{\Delta}_{0}\gg\tilde{\Delta}_{n\neq 0}$ and $\tilde{\Delta}^{\prime}_{n}\simeq\langle\Delta\rangle\simeq\frac{\tilde{\Delta}_{0}}{\sqrt{N+1}}$, where $\langle\Delta\rangle$ is the average order parameter in the bulk. The gap equation (22) gives for $n\neq 0$:

[TABLE]

where $\langle\Delta\rangle$ is defined through the usual gap equation:

[TABLE]

Then we use Eq. (27) to transform back to real space. By employing Eq. (26) we obtain a real-space $\Delta_{i}$ configuration which is quite close to the numerically obtained solution shown in Fig. 3 for $T=0$. The crucial point in this analysis is that when performing the summation (26), essentially, one represents an almost constant $\Delta_{i}\simeq\langle\Delta\rangle$ order parameter by some set of standing waves. Thus, clearly, oscillations inevitably should appear on the boundaries. In a general context this type of effect is known as the Wilbraham-Gibbs phenomenon Wilbraham (1848); Gibbs (1899). Similarly, the electron density exhibits oscillations and increases near a boundary. Another way to put it is that the Friedel oscillations of electron density are strongest close to the boundary. Most importantly, there are oscillations and a sharp increase in the local density of states (LDOS). These oscillations in the LDOS are presented in Fig. 4, calculated using the equation Zhu (2016):

[TABLE]

The spikes in the LDOS produce the increase in the gap, which then relaxes into the bulk at the length scale of the bulk coherence length $\xi$.

Although the gap function given by (27) exhibits characteristic oscillations at $2k_{F}$ as in the continuous case (Fig. 1a), it doesn’t give a good approximation for the critical temperature. Therefore, in order to calculate the critical temperature we use the linear gap equation for the tight-binding case. This is just a discrete version of (2, 3, 4) with the single-electron wave function and energies given by:

[TABLE]

[TABLE]

Here $N$ is the number of sites in the corresponding direction, $p_{n}=\frac{2\pi n}{N}$ and $k_{n}$ was defined in (23).

Firstly, let’s consider the periodic case. We transform to momentum space by $\tilde{\Delta}_{n}=\sum_{i=1}^{N}(w_{n}^{(p)}(i))^{*}\Delta_{x}$ with $\sum_{i=1}^{N}(w_{n}^{(p)}(i))^{*}w_{n^{\prime}}^{(p)}(i)=\delta_{n,n^{\prime}}$ for $n,n^{\prime}\in[1,N]$. Then the one-dimensional linear gap equation becomes:

[TABLE]

where

[TABLE]

Here $(p)$ means that the Eq. (31) is used. Equations (32) and (33) are trivially generalized to higher dimensions. The interaction strength $V$ is numerically calculated for given bulk critical temperature $T_{c1}$ and chemical potential $\mu$ from $\frac{1}{V}=D_{N}^{(p)}=D_{0}^{(p)}$.

Note that $t$ can be scaled away. Hence, we set $t=1$.

Now let us consider the case where in one of the directions the system is not periodic, i.e. it has open boundary conditions. Using the transformation (26) to go to momentum space, we obtain the following one-dimensional linear gap equation:

[TABLE]

where we used Eq. (30). $D_{n}$ and $A_{n,n^{\prime}}$ are defined as:

[TABLE]

[TABLE]

This equation is easy to generalize to higher dimensions by adding periodic directions. For example, for a two-dimensional system that is non-periodic in $x$ and periodic in $y$, this amounts to replacing $\xi\to E_{n_{x}}+E^{(p)}_{n_{y}}-\mu$ and performing the sum $\frac{1}{N_{y}}\sum_{n_{y}=1}^{N_{y}}$ in (35).

One can find $T_{c2}$ by solving the equations with fixed interaction. We followed an equivalent procedure that gives the same solution in a different way. We fixed the boundary critical temperature $T_{c2}$ and chemical potential $\mu$ to various values and numerically solved the right hand side of (34) for the biggest eigenvalue, which equals $\frac{1}{V}$. Combining this with data from Eq. (32), we obtain the relative difference in critical temperatures for a range of $V$ and $\mu$, see Fig. 5. For this we numerically simulated systems of $N_{x}=1000$ and $10000$ in one dimension, $N_{x}=1000,\ N_{y}=1000$ in two dimensions and $N_{x}=1000,\ N_{y}=100,\ N_{z}=100$ in three dimensions. The results are shown in Fig. 5. The boundary states we find are most prominent in one dimensional systems and the relative critical temperature difference is an order of magnitude smaller in two- and three- dimensional systems. The reason why the effect is smaller in higher dimensions is the following: when introducing another periodic direction, $A_{n,n^{\prime}}$ becomes effectively averaged over $[-2t-\mu,2t-\mu]$ in (34). This, as can be seen from Fig. 5a, will decrease the critical temperature difference. This leads to the conclusion that superconductivity can survive to substantially higher temperatures in low-dimensional systems. This includes the edges and corners of three-dimensional cuboidal samples and boundaries with large curvature in general. That is, these results show that if one increases the temperature in a cuboid sample, then at mean-field level, the following sequence of the transitions will take place: first, the gap will close in bulk, but will remain on the surfaces, increasing temperature further will make the gap vanish on the surfaces, but the gap will remain on the edges of the sample. Increasing temperature further will result in the gap closure at the edges, but the gap will survive to a higher temperature in the corners of the sample.

IV Transparency of interfaces

In this section we briefly comment on the implication of our findings for the other types of interfaces. The Fig. 6a shows an example of a numerical solution for the simplest interface between a superconductor and a normal material obtained numerically in the BdG model. There we change the onsite potential $\mu$ for the normal half of the system. When $\mu$ is small the interface is transparent and the standard proximity effect behavior de Gennes (1964) is observed. When increasing the onsite potential, which makes the interface more reflective, the behavior changes gradually and eventually the gap function develops a spike near the interface. We note that the width of the solution, both in this case and for the boundary state is macroscopic: i.e. it is of the scale of bulk coherence length. Therefore we believe that the effect may in principle be described in the Ginzburg-Landau or quasiclassical Zaitsev (1984); Kuprianov and Lukichev (1988) approaches. However, a naive quasiclassical approach will not capture the effect: width of the solution is given by coherence length but, the source of enhanced superconducting correlations is at microscopic length scale. Thus to describe this effect in these approaches a careful treatment of microscopic degrees of freedom is required. This problem of compatibility of the results with quasiclassical theory will be addressed elsewhere.

The problem of superconductor-normal interface is related to the the problem of different type of superconducting surfaces. Namely, in a real material a surface can be oxidized, have different chemical composition or different phonons, resulting in weaker coupling constant. On Fig. 6b we show how depending on the characteristic of the altered boundary layer, the gap can be suppressed or enhanced at a macroscopic length scale away from the boundary.

V Conclusions

In conclusion, we showed that in the BCS model there are several critical temperatures at the mean-field level. The superconducting gap survives at the faces of a sample up to a higher temperature than in the bulk. Furthermore, at mean field level it survives to a higher temperature at the edges and even higher at vertices. This suggests an explanation of the frequently observed disparity between critical temperatures determined by specific-heat and diamagnetic-response probes. A resolution of the phenomena, that we find, in specific-heat measurements would require a homogeneous cuboid sample with a perfect surface. A number of experiments have explicitly claimed observation of boundary superconductivity above the temperatures where the samples lost superconductivity in the bulk, including hints of such behavior in the specific heat alone Lortz et al. (2006); Janod et al. (1993); Butera (1988). The most advanced experimental analysis of this effect was performed on YBCO. There the conclusion about superconducting surfaces was reached by cutting samples and observing the same splitting of the phase transitions. In these experiments the difference between $T_{c1}$ and $T_{c2}$ is consistent with our calculations. We calculated, using several approaches, $T_{c2}$ for a clean, absolutely reflecting boundary. 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 but still would contribute to the diamagnetic response. Apart from that, the found boundary states are important for interface superconductors and interpretation of Scanning Tunneling Microscopy probes. The solution for the boundary gap and the sequence of phase transitions can be directly studied in fermionic ultracold atoms in a box potential Gaunt et al. (2013); Shin et al. (2007); Schirotzek et al. (2008). Beyond mean-field approximation, the bulk transition will be in the 3D XY universality class and the surface transition will be of Berezinskii-Kosterlitz-Thouless type, which can also be experimentally probed. On the other hand we show that going beyond the homogeneous BCS model and adding even a small layer with weaker coupling, i.e. due to oxidation, different chemical composition, or different electron-phonon coupling can lead to a suppression rather than enhancement of the gap.

Note added: During the review process we were pointed to the article Giamarchi et al. (1990) that analyzed a related model and arrived at the opposite conclusion that there is no boundary superconductivity in a system with an attractive BCS fermion-fermion interaction. Another opposite conclusion in Giamarchi et al. (1990) is the behavior of the order parameter near a boundary when the bulk of the system is superconducting. According to Giamarchi et al. (1990), for a BCS system with attractive interaction the gap exhibits only Friedel oscillations. In other words, it is enhanced only at a microscopic length scale, while averaging over microscopic length scales the order parameter stays constant at the length scale of the coherence length. Next, only by going to a different, specifically modified model was it possible to achieve a relative $T_{c}$ difference of the order $10^{-3}$. We investigated the discrepancy by performing momentum space calculations, similar to Giamarchi et al. (1990), and obtained a numerical solution for the unmodified model in one dimension. Based on that calculation, our conclusions are opposite: we do observe a higher critical temperature near the boundaries and the enhancement of the average gap at a macroscopic length scale. Our interpretation of the reason for the discrepancy is the following: The conclusion in Giamarchi et al. (1990) was obtained numerically. We note that discretizing the integral in the linear gap equation is a subtle problem – the lower $\hat{V}$ is the harder it is to obtain precise results. In our one-dimensional case we had to go up to $6\times 10^{4}$ points to obtain somewhat exact results for $\hat{V}\simeq 0.2$.

Acknowledgements

We thank Mats Barkman, Emil Blomquist, Andrea Benfenati, Filipp N. Rybakov, Attila Geresdi, Wolfgang Ketterle and Vladimir Krasnov for discussions. We especially thank Vadim Grinenko for pointing out the experiments claiming the observation of two critical temperatures after this work was published as an eprint. The work was supported by the Swedish Research Council Grants No. 642-2013-7837 and No. VR2016-06122, 2018-0365 and the Goran Gustafsson Foundation for Research in Natural Sciences and Medicine.