Microscopic derivation of superconductor-insulator boundary conditions for Ginzburg-Landau theory revisited. Enhanced superconductivity at boundaries with and without magnetic field

Using the standard Bardeen-Cooper-Schrieffer (BCS) theory, we revise microscopic derivation of the superconductor-insulator boundary conditions for the Ginzburg-Landau (GL) model. We obtain a negative contribution to free energy in the form of surface integral. Boundary conditions for the conventional superconductor have the form $\textbf{n} \cdot \nabla \psi = \text{const} \psi$. These are shown to follow from considering the order parameter reflected in the boundary. The boundary conditions are also derived for more general GL models with higher-order derivatives and pair-density-wave states. It shows that the boundary states with higher critical temperature and the boundary gap enhancement, found recently in BCS theory, are also present in microscopically-derived GL theory. In the case of an applied external field, we show that the third critical magnetic-field value $H_{c3}$ is higher than what follows from the de Gennes boundary conditions and is also significant in type-I regime.


I. INTRODUCTION
Superconductivity in the Ginzburg-Landau (GL) model [1] is described by a complex-valued field ψ(r), which is called an order parameter or gap. In the bulk of the sample, ψ is found as a minimum of the free-energy functional F bulk [ψ]. The form of this functional was microscopically derived first by Gor'kov [2]. To solve for ψ near a boundary of superconductor one has to take into account the influence of the material outside the sample. This is done by a microscopically derived boundary condition for ψ or, equivalently, by an additional surface term F surf [ψ] in free-energy functional. Namely, dropping vector potential (it is restored in Section VII), the free energy F of a superconductor placed in Ω is given by: The order parameter is found as a minimum of F . Hence the boundary condition at ∂Ω is: where n is unit vector pointing outside of the sample. We begin by reviewing how the value of γ changes the gap behavior near the boundary in different systems. A calculation by de Gennes [3,4] gave that γ > 0 at a boundary between a superconductor and a normal metal. Hence superconductivity is suppressed near such a boundary. One also defines the extension length Λ = K/γ > 0, which controls the range of induced superconducting correlation in the metal. On the other hand when γ < 0 superconductivity is enhanced. Such a situation is realized in several cases. Namely, for a contact with a superconductor with a higher critical temperature [5], for a local increase of electron-phonon coupling constants near the surface [6] and for superconductivity on twinning planes [7,8]. A phenomenological model with γ < 0 was analyzed in a number of works, see e.g. [5,7,[9][10][11][12][13][14]. Boundary conditions for the interface between superconductor and insulator were studied microscopically in [3,4,15,16], yielding the conclusion that it is a good approximation to set γ = 0. We will call the corresponding boundary conditions the de Gennes boundary conditions. However, the situation for the superconductorinsulator interface is not trivial. Namely, it was recently shown microscopically that boundaries of superconductors can have (i) higher critical temperature and (ii) the gap can be enhanced at the scale of the bulk coherence length. We call this enhancement of superconductivity 1 the boundary state. It was found in one, two, and three-dimensional superconductors in the tightbinding BCS model and for the one-dimensional continuous BCS model [19]. An earlier study of the standard three-dimensional continuous BCS model concluded that the boundary state is absent: the averaged gap near the boundaries is neither suppressed nor enhanced [17] or weakly suppressed [18]. Coulomb repulsion-induced boundary states were reported to form in certain cases in a three-dimensional continuous model where the interaction was attractive for small energies and repulsive for higher energies [17]. Superconductivity enhancement in the form of pair-density-wave (PDW) boundary states 1 Note, that, near a superconductor-insulator boundary, the gap oscillates with period 1/k F , where k F is the Fermi momentum. Some works term the first peak of these Friedel oscillations "enhancement" [17,18], even in the case where the averaged gap is not enhanced, which corresponds to the situation where de Gennes boundary conditions apply in a GL theory. Here we mean by enhancement an increase in the gap averaged over a length scale that is much larger than 1/k F . was found in spin-imbalanced superconductors [20,21] as well.
The situation is also controversial experimentally. Evidence for a substantially enhanced superconductivity near the boundary was reported in some elemental and high-temperature superconductors, see, e.g. [5,[22][23][24][25][26][27][28][29][30]. The enhanced surface superconductivity was also mentioned in the context of enhancement of critical temperature observed in granular elemental superconductors [31][32][33]. The effect was interpreted as the surface being described by a different Hamiltonian, based on a conjecture of different chemical composition, which may indeed be the case, especially in complex compounds or with enhanced phonon interaction [6,34]. However, higher critical temperature of the surface was reported also for clean elemental superconductors [28,29]. Recently the claim of direct evidence for surface superconductivity was reported using the newly developed direct probe [30]. Also, results in [35] hint for possible interpretation in terms of surface critical temperature, which indeed should, in general, depend on the orientation of the boundary relative to the crystal axes.
The results from microscopic calculations and experiments show that the superconductor-insulator interface is nontrivial. This is important for various applications. For example, boundaries play a big role in quantum devices such as superconductors-based qubits and singlephoton detectors, see, e.g. [36,37]. Moreover, the GL model remains the only nonlinear model amenable to the numerical solution at a substantially large length scale, required for modeling such devices. The latter provides additional motivation for this work to revise the derivation of boundary conditions in the GL model. Additionally, we resolve the ambiguity in boundary conditions when terms of higher-order in derivatives are added, see the discussion in [21]. This paper is organized in the following way: In Section II we set up a microscopic BCS model, which is used to derive the GL model. In Section III we derive that boundary conditions can be found from the mirror reflecting the order parameter in the boundary. The surface term is neglected. In Section IV we microscopically obtain surface term F surf for several models. The result is used in subsequent sections. In Section V we solve for a phase diagram that includes boundary states in the GL model for the spin imbalanced system. In Section VI we obtain the difference of bulk and boundary critical temperatures in the GL model Eq. (1). In Section VII we introduce the magnetic field and give a microscopic assessment of how γ < 0 enhances the third critical magnetic field H c3 .

II. THE MICROSCOPIC MODEL
Consider continuous-space fermionic theory with the BCS type local attractive interaction given by strength V > 0. We regularize the interaction by the Debye fre-quency ω D such that only electrons with Matsubara frequency < ω D interact. The path integral formulation of this model is given by the action S and the partition function Z (see Chapter 6.4 in [38]): (3) where a σ (τ, r), a † σ (τ, r) are Grassmann fields that correspond to fermionic creation and annihilation operators and depend on imaginary time τ , d dimensional space coordinates r and spin σ. Next, ε σ ≡ E − µ σ , where µ σ is the chemical potential, and T is the temperature. The single-electron energy is E ≡ E(i∇) with E(0) = 0, which is E(k) = k 2 2m for free electrons. It is assumed that E depends only on the modulus of k so that E(k) ≡ E(|k|). We consider a superconductor positioned in the Ω domain and an ideal insulator positioned everywhere else. To model the insulator we assume that µ σ is finite in Ω and µ σ → −∞ elsewhere.
We perform a Hubbard-Stratonovich transformation in the Cooper channel by introducing an auxiliary bosonic field ∆(τ, r): By introducing Nambu spinors A † = a † ↑ , a ↓ and A = a ↑ a † ↓ we rewrite the partition function as Then, by performing the Berezin integral, we integrate out the fermionic degrees of freedom: where F is the free energy. Next, we make mean-field assumptions: that ∆ is a classical field (does not depend on τ ) and that it does not fluctuate thermally and is found as a minimum of F . In this approximation the problem simplifies: where Matsubara frequencies are ω n = 2πT (n + 1/2) and we used Green's functions for spin σ electrons are determined from: Since µ σ → −∞ in the insulator, the single-electron wave functions will be zero there. This results in boundary conditions for Green's functions in the following way: for coordinate r b lying on the boundary ∂Ω we get Now let us consider Ω : x > 0. In that case the Green's function can be obtained in the form: where g σ (r) is bulk Green's function, r ≡ r − 2xx and k ≡ |k|.
Below we assume the following: where T c is the critical temperature, v F is the Fermi velocity, and we are interested in a slow-varying order parameter with momentum |q| Q (this justifies the expansion in q that we do later). Close to the transition to the normal state, ∆ → 0. Hence we can truncate the expansion in ∆ in Eq. (7). For usual superconductors, µ ↑ = µ ↓ . We also consider a superconductor with a spin imbalance. There, GL expansion is done near the tricritical point associated with the bulk Fulde-Ferrel-Larkin-Ovchinnikov (FFLO) state [39,40] (in which case |µ ↑ − µ ↓ | is of the order of T c ).
Combining Eq. (6) and Eq. (7) we obtain the quadratic-in-∆ part of the free energy: Figure 1. Illustrations of free-energy terms which are secondorder in the order parameter Eq. (16). Orange (green) lines denote bulk Green's functions of spin-up (spin-down) free electrons (reflection from the boundary is marked by a vertical line). Squares denote the order parameter ∆. In all cases, the particle-particle bubble, consisting of two Green's functions, is explicitly given by f (k, k − q), see Eq. (17). The top two diagrams with zero and two boundary reflections lead to usual bulk-like terms, where order parameters with the same momentum are coupled. They give terms proportional to D(q) in Eq. (16). The bottom two diagrams have one reflected from the boundary electron. The latter results in coupling between order parameters of different momenta (the last term in free energy Eq. (16)). As shown below, this gives rise to boundary states.
To simplify F 2 we need to perform Fourier transform. However, integrals over x in F 2 are over half space. To extend them to full space we note that Together with defining ∆(r) = ∆(r), without loss of generality we obtain: Performing the Fourier transform using Eq. (11) and ∆(r) = 1 (2π) d +∞ −∞ ∆(q)e iqr dq we obtain, see Fig. 1: , and hence ε(k) = E(k) − µ. Then the Fermi momentum k F : ε(k F ) = 0 and the Debye momentum is nonzero when k and k − q are on a Fermi sphere of radius k F and thickness 2k D .
Usually it is assumed that ∆ varies slowly with momentum q Q k D Eq. (12). However, this is true only for directions parallel to the boundary. By contrast, in the x direction, fast oscillations with q 2k F are present, the gap exhibits Friedel oscillations [17,19]. The existence of oscillations, means that f (k, k − q) has contributions from three different types of points on Fermi sphere illustrated in Fig. 2.
In this work, we are interested in the description of the boundary of a superconductor at the level of the GL model. In the GL approximation, the order-parameter field is coarse-grained and thus varies slowly in real space, which corresponds to small momentum |q| Q, below it is denoted by ψ(q) ≡ ∆(q). Whereas the order parameter that changes fast in the x direction (large momenta

III. SMALL-Q CONTRIBUTION TO BOUNDARY CONDITIONS: DE GENNES APPROXIMATION
We begin by reproducing the de Gennes microscopic boundary conditions under an approximation similar to that used in [3,4,15,16] considering the small-q contribution to Eq. (16). That is, in this section, we consider only the contribution from a slowly varying order parameter ψ. In this approximation we restrict |q| < Q in ∆(q). Then we obtain Since q is small in Eq. (12) we expand D(q) as: where c i are the usual bulk coefficients [41]. Here we carry out the calculation in a more general form that is also applicable for the case where one has to keep higher-order derivative terms, such as the case of spin imbalanced superconductors, including those in the FFLO Note that, equivalently, we could have defined ω D cutoff for energies |ε| < ω D and let the sum over Matsubara frequencies be unrestricted.
Fermi sphere and three types of configurations of momenta (a)-(c), that give rise to boundary states. In all these cases, particle-particle bubble f (k, k − q) (defined by Eq. (17)) is large. This is due to the fact that q (blue) and k (red) lie in the Fermi sphere (black circles). In all cases |k| − kF ∈ [−kD, kD]. By kx we denote momentum orthogonal to the boundary and by k || momentum parallel to it. Configurations can be described by the following values of kx and q: (a) |kx| kD and q is large. Namely, |q − 2kxx| Q. (b) |kx| kD and |q| Q. (c) |kx| kD and |q| Q. Note, that, for two-and three-dimensional systems all three appear, whereas for one-dimensional systems, only cases (a) and (b) appear with k || = 0. state [42,43]. Transforming Eq. (19) into real space gives: where terms with derivatives can be written in any integrated-by-parts form like −c 2 ψ∇ 2 ψ * + c 2 ψ∇ 4 ψ * + ... This is possible since integral over x in Eq. (20) is over (−∞, +∞) with ψ reflected in the boundary as ψ(r) = ψ(r). Hence when going back to the actual system Ω : x > 0 one obtains boundary conditions at x = 0, which depend on the chosen form of c i terms.
For example, if the c 4 term is not included, the c 2 term gives the usual [3,4,15,16] n · ∇ψ = 0 (21) This is the de Gennes boundary condition. Applying it dictates that the critical temperature of bulk and boundaries are identical. To derive it, note that ψ is continuous across x = 0 but not necessarily smooth. Hence, writing c 2 |∇ψ| 2 has no additional terms at the surface and boundary conditions Eq. (21) are obtained by variation as usual. Alternatively we could have picked −c 2 ψ * ∇ 2 ψ.
In that case Integrating by parts we obtain again c 2 |∇ψ| 2 , since the delta function in Eq. (22) compensates for integration by parts. Consider now the case where the c 4 term is included, such as, for example, in the GL model of a spin imbalanced uniform and FFLO systems [42,43]. Then c 4 |∇ 2 ψ| 2 with Eq. (22) shows that in order for energy to be finite, we need first of the two boundary conditions: where the second boundary condition is obtained from variation of the energy. This provides microscopically derived boundary conditions in de Gennes approximation for moderately spinimbalanced systems and FFLO systems. Boundary conditions Eq. (23) were used before for these systems, see e.g. [44,45]. Note that these boundary conditions eliminate the PDW boundary states discussed in [20,21] and contradict the phenomenological GL boundary conditions used there. At the same time, the PDW boundary states are unambiguously demonstrated in the full microscopic model [21]. We resolve that question below.
Note, that this boundary condition of the order parameter being reflected, ∆(r) = ∆(r), follows from the general property of the Green's function Eq. (14). Hence it is easy to generalize it to different systems. For example, in noncentrosymmetric superconductors [46] with local interaction Eq. (3) the analog of Eq. (13) just has a matrix G and ∆ [47], while the property Eq. (14) follows from the boundary condition for the Green's function Eq. (10). Hence the counterpart of de Gennes boundary conditions, in that case, is obtained from reflecting the fields at the boundary as well. For the simplest model they are equivalent to boundary conditions obtained from variation of the free-energy function, see [47].
To summarize this section: we reported generalized to GL models with higher-derivatives derivation of de Gennes boundary conditions. However, these conditions do not reproduce the superconducting boundary states. On the other hand, it was shown microscopically that these states exist on the macroscopic length scale [21] and hence they should be reproducible in microscopically derived GL models. This problem is resolved in the next section.

IV. THE BOUNDARY CONDITIONS BEYOND THE DE GENNES APPROXIMATION
In this section, we consider whether there is a nonvanishing contribution from the terms in Eq. (16) coming from an averaging of the fast-oscillating order parameter. We denote them F 1 2 such that F 2 = F 0 2 + F 1 2 . These terms have a large-q counterpart of F 0 2 Eq. (19) and (a), (b), (c) parts of F 2 , see Fig. 2: where Q Q and q || are q components parallel to the boundary if there are any. By dp we denote integral +∞ −∞ dp excluding |p| k D . We defined: Eq. (24) is simplified to: where we used +Q −Q dqx 2π ψ(q x , r || ) = ψ(0, r || ). Note, that the full field ∆(x) is zero at the boundary since electrons are perfectly reflected from it and hence Green's functions are zero there. However, the GL order parameter ψ represents only a slowly varying part of the pairing field, which can be nonzero at the boundary.
By varying Eq. (26) with respect to ψ f we obtain the solution for fast-oscillating part of the order parameter: Inserting it back into Eq. (26) we get the surface term: Now let us analyze that contribution in various dimensions.

A. Superconducting wire
For one-dimensional systems, we have no contribution associated with configuration (c), see Fig. 2, in energy Eq. (24) and hence there should be no f 1 (0, 0) term in the expression for the boundary term Eq. (28). However we can use the same formula Eq. (28) since for d = 1 we have f 1 (p, q) = f (p, q) and hence f 1 (0, 0) = 0. We simplify Eq. (28) as: and estimate D(2p) where the density of states at the Fermi level is N = 1 πv F . By performing the integral over p we obtain and Ψ (n) are polygamma functions of order n. For h = 0 Eq. (30) reduces to: where ζ is the Riemann zeta function. Therefore, in one-dimensional GL theory, there is a boundary term for the interface between a superconductor and a vacuum Eq. (28). The term has a negative microscopically derived prefactor γ. This implies that the gap is increased near the boundary and there are superconducting boundary states [19]. The conclusion applies both to quasifree and band fermions.

B. Planar superconductor
Consider the one-dimensional boundary of a twodimensional sample. We estimate f 1 (p, p) for k F − k D |p| as: with where α is defined in Eq. (42), is the bulk critical temperature, γ E is Euler gamma, Ψ is digamma function and N is density of states at Fermi level. In this work, we consider the case of d dimensional system in the BCS limit where N V is small. For the two dimensional case, N = k F 2πv f . The function D(2p) for k F − k D |p| k D is then given by: Eq. (32) and Eq. (34) allow us to compute γ in Eq. (28) up to logarithmic accuracy: Therefore, similarly to the one-dimensional case of almost free fermions, we recover the boundary states at the level of the GL theory.

C. Three dimensional isotropic sample
Next we consider the two-dimensional boundary of a three-dimensional sample. In this case f 1 can be estimated as: The density of states of a three dimensional supercon- and D(2p) 0 for |p| k F . This allows us to calculate the integral in Eq. (28), which gives: Here c is the cut off parameter of order one defined by dp = 2 ∞ ck D dp. Eq. (38) shows that we cannot justifiably calculate γ using this analytical approach since it depends on c. Namely, as seen from Eq. (28), in our approximation we can obtain γ only of order f 1 (0, 0) or larger, whereas the leading-order contribution to the integral gives: dp 2π f 1 (0, 0). Hence, this level of approximation indicates that γ is very small for this model: of order Eq. (38) or smaller, or could be zero. This is the property of the special case: a twodimensional boundary of an isotropic continuous BCS model in three dimensions. Note, that the situation is different in the three dimensional tight-binding BCS model [19].

D. The boundary states and anisotropy in three dimensions
Since many of the superconducting materials of current interest are strongly anisotropic, let us explicitly consider the effects of anisotropy. Consider the three-dimensional model that has single electron energy E (k x /a) 2 + k 2 || , which is anisotropic for a = 1. Then γ is given by the same formula Eq. (28) with the replacement p → p/a inside the single electron energies E. Hence, using Eq. (36) we obtain the leading-order estimate: where k F is the Fermi momentum parallel to the boundary. Hence, for a > 1, the superconductivity is enhanced near the boundary and there are boundary states. In contrast, the gap is suppressed for a < 1. In other words, boundary states are present if the Fermi sphere is stretched in the direction orthogonal to the boundary. Note that [35] reported a difference in T c for samples with different orientations of the surfaces relative to crystal axes using the stiffnessometer experiment. The stifnessometer setup was proposed to resolve surface superconductivity [30].

V. BOUNDARY STATES IN SUPERCONDUCTOR WITH IMBALANCED FERMIONS
Now we consider the case where superconducting pairing takes place in a model with spin imbalance, i.e., unequal densities of spin components. In an infinite system, when there is a critical disparity of Fermi momenta of spin-up and spin-down fermionic components, the system undergoes a phase transition into an inhomogeneous FFLO state [39,40]. In such a state, the system has a modulation in the phase or modulus of the orderparameter field. At the level of the Ginzburg-Landau theory, a phase transition into such a state manifests itself through the coefficient in front of the quadratic gradient term becoming negative. Therefore, for the energy to be bounded from below, one needs to retain higher-order gradient terms with positive prefactors.
Combining the results of the microscopic derivation for bulk [42] and boundary Section II we obtain the GL model with spin imbalance: 40) Let us now study the problem of the boundary conditions in the presence of the higher-order derivative terms. These conditions are obtained from considering the mirror reflected model. The condition of finiteness of energy and variation of F with respect to ψ * gives us the boundary conditions: where n is a unit vector pointing outside of the sample. Note, that boundary conditions Eq. (41) differ from the phenomenological boundary conditions used in [20,21]. Whereas the de Gennes boundary conditions, previously used for FFLO systems [44,45] correspond to setting γ = 0 in our boundary conditions Eq. (41). We derive that γ is not zero. In Eq. (40) and Eq. (41) γ is given by Eq. (28) while the other parameters are [42]: Let us consider now the boundary physics of a spinimbalanced superconductor. First, we determine when the superconductor transitions to a normal state. We assume that this transition is of second order and hence ψ → 0 at the transition. Hence it is sufficient to solve linearized GL equations that follow from the variation of Eq. (40): This gives us that bulk transitions to the normal state at, see The system however has a boundary superconducting state with critical temperature which is higher than the bulk critical temperature. Let us now consider the phase transition from a superconducting boundary state to a normal state. For that matter consider sample positioned at x > 0. Then Eq. (43) should be solved together with the boundary conditions Eq. (41) and the requirement that the order parameter goes to zero at infinity. We obtain the order parameter configuration of the boundary state, see Fig. 4: which satisfies the second condition in Eq. (41) at the transition from the superconducting boundary to the normal state, see Fig. 3: (for example, a system in an in-plane magnetic field). Lines denote superconducting phase transitions. The bulk phase transition between normal and superconducting states is denoted by solid lines according to Eq. (44). The boundary remains superconducting for higher temperatures. The dashed lines denote a phase transition from surface superconducting to normal state according to Eq. (46) and the dot-dashed line according to Eq. (48). When the spin imbalance h is large enough, bulk and boundary states turn from sign definite (red) to periodically modulated in space states (blue). Gray hatching denotes the region where the usual GL model Eq. (1) is bounded from below (it works best for T → Tc1 and h → 0). Brown hashing shows the region where the GL model with higher-order derivatives Eq. (40) is convergent (the best at the tricritical point). Note, that this phase diagram has a boundary PDW state extending beyond the tricritical point and there is a smooth transition to a non-PDW boundary state, which agrees with the phase diagram obtained in a Bogoliubov-de Gennes formalism [21].

VI. BOUNDARY STATES IN A CONVENTIONAL SUPERCONDUCTOR
Next, we consider solutions for the conventional GL model Eq. (1) with that derived in the above boundary conditions Eq. (2). The model can be obtained from the more general expression Eq. (40) by setting K = K 1 = ν = 0. Hence bulk transition from superconducting to normal state takes place at α = 0, see Eq. (44). When the derived in the above microscopic boundary conditions are used the model exhibits superconducting boundary states. The transition from it to a normal state is obtained from the boundary conditions Eq. (2) and Eq. (43) at K = 0. The superconducting boundary-state solution that follows from that equation is: Hence transition to a normal state is given by the following condition, which corresponds to the boundary condition Eq. (2): This is the GL approximation of the superconducting boundary state obtained earlier as a solution of the full microscopic theory [19]. Note, that the coarse-grained GL field is smoothly varying and the enhanced pairing correlations are modeled as a source in the form of a boundary integral, yielding the boundary conditions Eq.
From Eq. (48) we obtain that for zero population imbalance, the boundary state transitions to normal at T c2 that is larger than bulk critical temperature T c1 : Hence, using Eq. (31) for a one dimensional system, τ is equal: For N V → 0 this is equal to τ = 7ζ (3) 4 (N V ) 2 , which we previously obtained in a one dimensional model without the Debye frequency [19] (note that, in [19] there is a typo and the rescaled interaction is actuallyV = N V ).
For two a dimensional system from Eq. (35) we obtain: For a three-dimensional anisotropic system with a > 1, we obtain from Eq. (39): A. Analogy to wetting As a side note, one can imagine that the plot of ψ in Fig. 4 for the usual GL model is a vertical cross-section of a tank filled with water, with ψ being the height of the surface of the water. This is not a coincidence. The 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 the energy per unit surface, dl is the length of the surface, ρ is the density of water, and g is the gravitational constant. Note, that we implicitly redefined ψ → λψ, where λ is some dimensional constant, so that for new [ψ] = [x]. Then the total energy is which is similar to the GL model Eq. (1), if we substitute σ → 2K and ρg → 2α.
Then the problem of the boundary states in a superconductor can be related to the problem of adhesion of water to a wall. In the latter case, the boundary condition is set by fixed contact angle, which is equivalent to fixing n · ∇ψ = const. This boundary condition is similar to Eq. (2) for a superconductor. Note that in this analogy, the superconductor-insulator interface behaves like a hydrophilic surface. For other interfaces, like superconductor-normal metal, or an interface with certain types of different boundary layers [19], the gap can be suppressed near the boundary, which corresponds to a hydrophobic surface.

VII. Hc3 FROM MICROSCOPICALLY DERIVED GL THEORY REVISITED
For conventional boundary conditions γ = 0, superconductivity in type-II materials survives near surfaces at the magnetic fields up to H c3 , which is higher than the critical field associated with the disappearance of superconductivity in the bulk H c2 . The boundary conditions that we derived have direct implications for the third critical magnetic field H c3 .
In the conventional picture [48] it is described by solving the linearized GL equation by using the standard de Gennes boundary conditions n · ∇ψ = 0. Note that it has been observed earlier that numerical solution in the fully microscopic Bogoliubov-de Gennes theory (i.e., obtained beyond quasiclassical approach) is not consistent with this picture [49], but should be more robust, however, that work did not determine H c3 . We are in a position now to calculate H c3 . Note that the problem that we will study below, namely H c3 in a GL model with an included surface term have been studied on phenomenological grounds in the past in the context of superconductors with modified surface layers, and superconductors with enhanced superconductivity on twinning planes [5,7,8,12]. Remarkably, enhanced H c3 was observed for ordinary surfaces of the elemental superconductors and experimental papers have been explicitly raising the question of whether that originates in enhanced superconducting properties of surfaces of unknown origin [28,29]. Our goal here is to calculate H c3 in the microscopically derived GL theory corresponding to a regular boundary of a standard BCS superconductor.
To include magnetic field one simply replaces the derivative with a covariant derivative in Eq. (3) as ∇ → ∇ − ieA, where A is the magnetic vector potential and e is the electron charge. Next, Green's function with magnetic field G A σ is obtained as where G σ is the Green's function for zero magnetic field defined in Eq. (11) and φ(r, r ) eA · (r − r ). Note that, in this approximation, it is assumed that A is very slowly changing and hence it can be A A(r) or A A(r ) or anywhere close to r, r . For details see [15,41,50]. In our calculation, however, it is convenient to choose φ(r, r ) eA(r) · r − eA(r ) · r (55) and to extend A(r) = A(r) 3 . Then property Eq. (14) is satisfied for G A σ as well. It means that the derivation for the model with the magnetic field is similar to the one outlined in Section II. The only difference is that now ∆ is replaced by: Hence, for the resulting GL model the magnetic field amounts to replacing In an external magnetic field, the boundary condition with γ < 0 leads to an increased critical magnetic field H c3 . It has a different dependence on temperature compared with the standard textbook derivation [48], see,for example, the phenomenological discussion in [5]. Similar observations were made in phenomenological studies of superconductors with twinning planes.
Here we compute H c3 with the microscopic boundary conditions derived above. Consider a two-or threedimensional system and assume that the external magnetic field is directed along the z direction and equals H. Hence we can set a gauge for vector potential so that only nonzero component is A y . Then the transition to the normal state is obtained by solving the linearized GL equation in terms of ψ(x): together with the boundary conditions Eq. (2). We set A y = (x + x 0 )H, where x 0 is to be optimized to get the highest H. Then the solution to Eq. (58) is where D ν (x) is parabolic cylinder function. To find x 0 , it is convenient to calculate the first dψ integral of Eq. (58). By using boundary conditions Eq. (2), we obtain (60) This equation sets the relation between x 0 and H, when the solution for ψ Eq. (59) is inserted. Since we search for the largest H, a derivative of Eq. (60) with respect to x 0 should be zero. It gives: Here we picked the minus sign so that Eq. (60) can be satisfied. Next, in order to find H we can either solve Eq. (60) or the boundary condition Eq. (2). We solve the latter, which using the expression for x 0 Eq. (61) and rescaling, amounts to solving numerically for η for a given a in: where H ν (x) is a Hermite polynomial, a = 1− αK γ 2 and η = . The H obtained then is equal to H c3 . See Fig.  5 for a plot of H c3 . The H c3 is, therefore, higher than in the original Saint-James de Gennes derivation [48]. Also, note that H c3 should exist for type-I superconductors in significant temperature range.

VIII. CONCLUSIONS
We considered the generic BCS model for spinbalanced and spin-imbalanced fermions. From that N H as a function of temperature T . Bulk transitions to the normal state for fields higher than Hc2 = − α 2|e|K . In the derivation [48] with zero boundary term (γ = 0), the surface transitions to the normal state at higher field Hc3 1.69Hc2. From Eq. (62) we obtain that, if in Eq. (1) the boundary term is present with γ < 0, then surface superconductivity is enhanced and boundary transitions to normal at field Hc3 higher than that for γ = 0. Here we have chosen γ = −0.05 N v F T c1 .
model, we derived the boundary conditions for the GL theory for the interface between a superconductor and an insulator. We showed that the free energy of a superconductor acquires an additional term given by the surface integral of γ|ψ| 2 . The physical origin of this term is the fact that near a well-reflecting boundary, the total gap oscillates with k F momentum. This oscillatory part is coupled to the averaged gap ψ in GL theory, leading to an additional surface term. We obtained that γ < 0 for one-and two-dimensional continuous BCS models. For the three-dimensional isotropic BCS model, the surface term is beyond the resolution of our analytical approach. Whereas for an anisotropic three-dimensional model we have shown that γ can be positive or negative. The negative γ leads to enhanced superconductivity near boundaries. Note that, for the tight-binding BCS model the surface superconductivity exists in all dimensions [19], which also implies boundary conditions with γ < 0.
To obtain boundary conditions we showed that the following procedure can be applied. Consider the GL model that in general has the highest in derivatives term of order k. Then one should reflect the order parameter in the boundary and write the free energy as an integral over the whole space. This reflection automatically applies proper boundary conditions. Namely, one can write kinetic terms in any integrated-by-parts form and search for a minimum of the total free energy. Boundary conditions are obtained from the condition on energy to be finite and from the variation of this functional with respect to the order parameter. As a result, we obtained that usually all normal to boundary odd derivatives of ψ of the order less than k − 1 should be zero, whereas k − 1 derivative will be proportional to γψ.
For the standard GL model with second-order derivatives, k = 2 and hence we obtain boundary conditions n · ∇ψ = − γ K ψ, see Eq.
(2). GL model for spinimbalanced systems requires taking into account fourthorder gradient terms and hence has boundary conditions n · ∇ψ = 0, n · ∇ 3 ψ = γ K ψ, see Eq. (41). Using these new microscopically derived boundary conditions we obtained superconducting boundary states in the conventional GL model and revised the calculation of PDW boundary states in the spin-imbalanced GL model.
The obtained boundary conditions allow GL theory to account, in a microscopically accurate way for boundary states that were found earlier in fully microscopic solutions of BCS theory [19]. Namely in the GL model for γ < 0, in zero external magnetic field, the superconducting gap is larger near the surface than in the bulk, and superconductivity survives for higher temperatures. Since microscopic calculations show that superconductiv-ity is more enhanced at the edges and corners of a threedimensional sample, to model these effects one should add analogous extra contributions for corners and edges.
By adding an external magnetic field we revised the theory of the third critical magnetic field H c3 for a BCS superconductor. The surface effects make this field larger and extending in a type-I regime, compared with results obtained using de Gennes boundary conditions. We note, that these surface effects can be described in a quasiclassical approach if one augments the theory by taking into account higher-momentum contributions to the boundary conditions.