Superconducting size effect in thin films under electric field: mean-field self-consistent model

We consider effects of an externally applied electrostatic field on superconductivity, self-consistently within a BCS mean field model, for a clean 3D metal thin film. The electrostatic change in superconducting condensation energy scales as $\mu^{-1}$ close to subband edges as a function of the Fermi energy $\mu$, and follows 3D scaling $\mu^{-2}$ away from them. We discuss nonlinearities beyond gate effect, and contrast results to recent experiments.


I. INTRODUCTION
Quantum oscillations in superconducting properties due to size quantization in thin films were predicted early [1][2][3][4] , and they were later observed in metallic films [5][6][7][8] . Modification of superconducting properties by changing the electron density by electrostatic fields was also observed, [9][10][11][12][13][14] and is best studied in high-Tc superconductors where the charge density can be low enough to enable efficient gating. Generally, modifications of critical temperature T c and critical current I c have been reported. Modification of I c only was also recently reported in Refs. 14 and 15 in metallic thin-film samples, but the proper interpretation in the latter is still unclear.
Electrostatics of superconductors is an old problem (see e.g. Ref. 16 for a historical review), and the effect of electric fields on superconducting surfaces are theoretically discussed in several works. [17][18][19][20][21][22][23][24] In these, effects on the amplitude of superconductivity (T c ) are usually related to modulation of electronic density of states (DOS), which is also what contributes to the quantum size effects. A common approach is to consider "surface doping" and assume the DOS is modified within a Thomas-Fermi screening length from the surface. Self-consistently screened calculations in superconductors have been previously discussed in Refs. [25][26][27], in a different context. For the normal state, there is a large literature on microscopic calculations with surface screening, which are routine today e.g. using density functional theory. 24,28 Modification of I c on the other hand is often assumed to come from changes in the vortex surface pinning potential. 14,29 In a simple picture, a static electric field appears as a perturbation of the potential that confines electrons within the thin film. Static fields generally extend up to a screening length from the surface, and so their effect decreases towards high charge density. Although the effects increase with the applied electric field, achievable field magnitude is limited by electric breakdown (e.g. via field emission 30 ).
Electrostatic gating of superconductivity in the BCS mean field picture relies on electron-hole asymmetry within an energy window determined by the order parameter and Debye frequency centered at the Fermi level. [31][32][33][34] In a simple clean thin-film model, strong asymmetry naturally exists in the form of the steplike multiband 2D DOS, which also gives rise to the quantum size effect, and the picture also extends to weakly disordered samples. The only question is to what degree the DOS asymmetry is retained, even though sharp features in the DOS are smeared by disorder, 35 and when samples cannot be significantly gated (metallic regime 15), since the Fermi level is not necessarily fixed at a sensitive point. Regardless, sharp DOS features can increase the charge density range in which electrostatic effects are large enough to be observed. Motivated by the recent experimental results 15 where large effects were seen, we revisit the problem.
In this work, we write down and solve a simple meanfield model for superconductivity in thin films under electric fields, including self-consistent screening. We point out connections between the dependence of electrostatic energy on superconductivity and modulation of superconductivity by fields, and discuss applicability of "surface doping" models in this picture. We also discuss to what degree nonlinear effects beyond linear electrostatic gating could appear in strong fields. We conclude that effects such as observed in Ref. 15 likely are not present in the model considered.
The manuscript is structured as follows. In Sec. II we introduce the mean-field model considered and discuss results obtained for the electric fields and modulation of superconducting properties. Sec. III concludes with discussion.

II. MEAN-FIELD MODEL
Self-consistent electrostatic screening and the size effect on superconductivity in a clean superconducting metal in a simple mean-field approximation is convenient to consider starting from a Hartree-Bogoliubov free energy. It can be obtained 34,36-38 by decoupling a longranged Coulomb and a (retarded) superconducting contact interaction via Hubbard-Stratonovich transformations, and considering only the classical saddle point in arXiv:1903.01155v1 [cond-mat.supr-con] 4 Mar 2019 the static limit: Here, G is the electron equilibrium Green function, U is a background potential, µ chemical potential, φ is equivalent to the static electric potential, ∆ the superconducting order parameter, and ρ ion and external charge density. The electron charge is −e and we use units with = k B = 1. The first term in the free energy is the electronic contribution, and the second part contains the electrostatic and superconducting mean-field contributions. In the absence of currents and magnetic field, at saddle point with suitable gauge ∆ can be chosen realvalued and the values of vector potential and phase are zero. Above, φ has to be taken as the saddle-point solution, which as typical for variational Poisson does not minimize F .
Variations vs. φ and ∆ give the Poisson and BCS selfconsistency equations: where G satisfies the Gor'kov equations G −1 G = 1 under the self-consistent potentials. We also here consider a BCS weak-coupling model, with ∆(ω) = ∆θ(ω c − |ω|), with the coupling λ taken as constant and the cutoff ω c similar to the Debye frequency. In bulk, the BCS gap equation is then directly ∆ = 2ω c e −1/(N0λ) with N 0 the DOS per spin at Fermi level.
For uniform system, expanding G in Eq. (3) to lowest order in φ results to RPA (q)q 2 φ(q) = δρ(q), where RPA (q) = 0 − e 2 q 2 Π(q; ∆) is the self-consistent static dielectric function of a clean superconductor. 39,40 In this model, the static fields are screened, and external charge density affects the electronic DOS, but not the Coulomb effect 41 to λ. The latter is due to considering mean-field with the decoupling assumed; corrections appear from fluctuations of φ (see e.g. Ref. 42 for explicit calculations), or on mean-field level with different decoupling 43 .
Aiming to describe effects on a qualitative level, we now consider a simplified model, similar to those used in several previous studies of the quantum size effect in superconducting thin films. 1,3 A confining potential U (r) is taken to be an infinite quantum well at |x| < L/2, which supports some number of populated 2D electronic subbands (see Fig. 1). In a static problem without currents, the electric field is perpendicular to the metal surface,  (6) in a superconductor is determined by the density of states and a occupation factor broadened by the superconducting interaction. and the problem is inhomogeneous only in x-direction. Moreover, we take as a variational Ansatz ∆ spatially constant 2 inside the well; the resulting energies will then be upper bounds to the exact solutions.
With these assumptions, the problem is elementary and mostly given by known results, 2,44 and can be solved without further approximations. First, where A N is the normal-state spectral function (per spin), u, v = [ 1 2 (1 ± ξ )] 1/2 and = ξ 2 + ∆ 2 . Due to the spatial symmetry, the problem reduces to one dimension. The normal-state DOS per volume is where ξ n are the 2D subband edges and V the film volume. The subbands and the potential φ are obtained from the Schrödinger-Poisson problem, Eq. (3) with where ψ n (x) are the transverse wave functions of the 2D subbands. Here, γ describes the contribution to the charge from each subband: where n(ξ) → u 2 f 0 ( ) + v 2 (1 − f 0 ( )) for ω c → ∞ and f 0 is the Fermi function. 45 The occupation factor n is broadened by the interactions in a window ∆ around the Fermi level, with the deviation from the Fermi function starting to decay more rapidly beyond the interaction range at |ξ| ω c . Variations in the DOS within this window contribute to the charge response of the amplitude of superconductivity (see Fig. 1). 31,32,46 To be specific, we assume an external charge density outside the sample (e.g. capacitor plates with constant charge density) such that the amplitudes of the electric fields at the surfaces are fixed, −∂ x φ(± L 2 ) = E ± . Numerically, the nonlinear Poisson problem can be solved iteratively, 47 for a fixed value of ∆.
The condensation energy f ns ])/V per volume for fixed ∆ now depends only on the density of states. 2,44 Via direct calculation, where we note that Here, g(y) → where η(y) = arsinh y + ( y 2 + 1 − |y|)y. The selfconsistent value ∆ * is attained at a solution of f ns (∆ * ) = 0.
It is also possible to express the electrostatic energy directly in terms of the self-consistent electric field, at small field strengths. Consider expansion of the electronic energy around a reference electric potential, considering small φ 1 = φ − φ 0 and ρ 1 = ρ − ρ 0 : where n e is the electron density and Π the density response function. Solving the resulting saddle-point equation for φ 1 and substituting the solution into F gives, after integration by parts: where C is independent of ∆. In this order of expansion in small φ 1 , the additional electrostatic field energy in (19) coincides with the standard expression. The linear term ∼ ρ 1 φ 0 describes a gate effect on superconductivity, which in this approach we see to be related to the ∆-dependence of the equilibrium potential φ 0 . Using Eq. (21) the quadratic part can be expressed as ∼ φ 1 RPA φ 1 . It corresponds to a (quantum) capacitance modulation 22,23 by superconductivity. The result (19) can be directly used for computing δf ns (∆) (if δφ ≡ φ[∆] − φ[0] is known) and is equivalent with (16) in the small-field limit. However, due the ∆-dependence of φ 0 it is not necessarily very practical to compute, as solving the nonlinear Poisson problem is still required. However, the above expressions can be used as a consistency check.
As noted above, we consider charge density ρ = ρ 1 +ρ 0 where ρ 1 outside the sample fixes the electric field at the surface. Finally, we need to specify the background ("ion") charge density ρ 0 . The electric potential due to ρ 0 together with U gives the pseudopotential for the electron system. 48 For simplicity, unless otherwise mentioned, below we assume ρ 0 = en e [∆ = 0, φ = 0, µ], which results to φ 0 = 0 as the solution in the normal state, and µ becoming the parameter that fixes the charge density in the normal state. This is of course a crude toy model of the surface electron behavior, even within Hartree-type models 28 , but likely modifies mainly the precise positions of the subbands but not the main qualitative features of the effect of the screening of external charges on superconductivity.

A. Size effect in electric field
In the same way as the variation in thickness, 1,2 gating by a surface electric field can in principle make a single subband edge ξ n to cross the Fermi level, which results to a jump in superconducting properties. Such response can be larger than in bulk material, and is not captured by "surface doping" models often used for the electric field effect, 17,18 where the LDOS ν(x, ξ) is assumed to be modified in a surface layer of thickness of a screening length λ T F according to bulk relations. In addition, the field screening is not exactly Thomas-Fermi type, but this causes less relevant changes than the difference in the DOS.
The order of magnitude of δf ns can be estimated in a Thomas-Fermi approximation. Taking where k n = πn/L and q(z) = z 2 /(1+z 2 ). From Eq. (17), and keeping only the smallest |ξ n | < ω c , where f ns,3D = − 1 4 mk F π 2 ∆ 2 is the bulk 3D condensation energy and a 0 = 4π 0 /(me 2 ) the Bohr radius. The above result is valid in the leading order in ∆, as φ is assumed to be independent of it. The factor q(2λ T F k n ) in reality depends on details of the screening, and below we consider it as a constant of order of magitude 1.
Including the next-order eigenvalue perturbation δξ (2) n in (16) and considering terms of order E 2 gives the second-order correction, where n = k means the limit ξ k → ξ n . This energy contribution is associated with the change Π  24) is not the only contribution to the E 2 term, and solving the selfconsistent electrostatic problem is in general required. 49 Conversely, calculation of the effect of superconductivity on the dielectric function requires taking the selfconsistency of ∆ = ∆ * [φ] into account. 40 The overlap factor q above depends on how accurate the Thomas-Fermi screening assumption is close to the surface. For the simple problem here, we can solve the Poisson equation numerically. Such a solution is illustrated in Fig. 2a  screening is not fully exponential, but the electric potential exhibits 1/k F oscillations. The correction δφ = φ(∆) − φ(∆ = 0) to the equilibrium electrostatic potential from superconductivity is small in the high charge density regime considered here. The chemical potential is chosen to be close to a subband edge in the figure.
The corresponding dependence of δf ns on the electric field magnitude is shown in Fig. 2b, together with the corresponding result from Eq. (23). The electrostatic energy expression (20) is also shown, and coincides with the exact result in the small-field regime. Generally, the electric field effect is appreciable only for E + L µ. In the estimate from Eq. (23), we here set q(z) = 1, to account for the expectation that likely for the true screening length λ T F k F 1. The second-order correction (24) is neglibile for these parameters, being higher order in λ 2 T F eE/(Lω c ), and the nonlinearity visible in the result originates from g(ξ).
The linear gating effect can be suppressed with a charge-neutral field configuration E + = E − = E, which corresponds to an experiment using floating gate electrodes (i.e. placing the system inside a plate capacitor). The result from Poisson equation for this case is shown in Fig. 3, together with the result for δf ns . Imposing the field on both sides produces a larger δφ. However, as the linear contribution to the free energy cancels, the modulation δf ns arises from the next-order effect and is an order of magnitude smaller than with the gate effect. Although the energy can still be expressed also via Eq. (20), the eigenvalue perturbation result (24) does not agree as well, as expected.
Whether electrostatic effects are significant depends on how large the modulation δf ns is compared to f ns . The dependence of their ratio on the chemical poten-  tial, and hence charge density, is shown in Fig. 4, at a relatively large external field. The magnitude of δf ns depends strongly on whether the chemical potential is located near a band bottom, where the effect is amplified (see Fig. 1), which produces the oscillations visible in Fig. 4. When µ is close to a subband bottom, the magnitude appears to be captured well by Eq. 23 (dashed line). When the chemical potential is not close to a band bottom, depending on the ratio between the subband spacing and the cutoff ω c , the electric field effect can vary by order of magnitude. Note that as long as |ξ n | ω c for some n, the result is dominated by the smallest ξ and the cutoff ω c < ∞ is of limited importance. The sum (16) is convergent also for ω c → ∞. However, these results are based on the simple weak-coupling model for superconductivity, and the precise shape of the modulation may be sensitive to details of the interaction. Regardless, from the above results one can see that the relative magnitude at resonance scales as ∝ ∆/µ, and not as (∆/µ) 2 as one would expect for the amplitude response in 3D bulk 39,40 . Away from the subband edge resonances, δf ns ∝ (∆/µ) 2 .
The self-consistent value of ∆ * , f ns (∆ * ) = 0, is shown in Fig. 5a as a function of the film thickness, showing the well-known quantum size effect. 1,2 The corresponding dependence on the chemical potential is shown in Fig. 5b, for several values of the external electric field. In this figure, it is obvious that the electrostatic field simply gates the system: the size effect physics is dominated by the ξ n closest to the chemical potential, so that the gate-induced shift δξ n is identical to a chemical potential shift −δµ.
The above discussion can be compared to a surface doping model, where the DOS is assumed to change in a surface layer of thickness λ T F , and the system is considered a superconducting bilayer in the Cooper limit. In such a model,

B. Superfluid weight
The effect of the electrostatic field on the phase fluctuations can be studied via the superfluid weight D s ij , 50 which describes the free energy cost of superflow ∆(r) ∝ e 2iA·r : The "vector potential" A describing the superflow can be introduced in Eq. (2) by replacementk →k + A. The calculation of D s is standard for multiband BCS superconductor. Since the current operators in y/z directions do not here couple different bands, the result for i, j ∈ {y, z} is D s ij = δ ij n n s (ξ n )/(mL) where n s (ξ n ) is the BCS superfluid density: 50 where n(ξ) is given by Eq. (8) and = (ξ ) 2 + ∆ 2 . As well known, n s (ξ) → n e (ξ) = 2mγ(ξ) at T = 0. The electrostatic modulation of the superfluid stiffness is then similar to that of the charge density, i.e., small in the metallic regime. Similar conclusion then applies to the phase stiffness, and quite likely also to the phase-slip energy barrier 51 . These results however apply in the clean limit.

III. DISCUSSION AND CONCLUSIONS
We discussed an elementary BCS/Hartree-Bogoliubov mean field model for the size effect under self-consistent electrostatic fields in superconducting thin films, and studied it based on numerically exact solutions. As the size modulation in superconducting properties decays relatively slowly with increasing charge density, it increases the response to applied electric fields, effectively changing the small parameter from (∆/µ) 2 to ∆/µ for fine-tuned values of µ, also in films thick compared to the screening length.
The mean-field approach likely is not useful in describing atomically thin, or strongly disordered and resistive samples, where fluctuation effects matter. Phaseplasmon fluctuation effects in principle can be included in the approach above in a standard way by expanding in Re ∆, Im ∆ and V = −iφ around the mean-field solution. A priori, in view of some existing results, 38,42 however, it's not clear why such corrections would depend strongly on the external electric field.
Large electrostatic size effects in thin-film systems are expected to be visible mainly in relatively low charge densities, e.g. semiconducting materials. As noted in previous works, 29 it appears likely this is a main effect in high-Tc superconductors. The modulation of screening by superconductivity will also appear in proximity systems, e.g. in semiconductor/superconductor hybrids 27,52 recently considered as Majorana fermion platforms.
With regard to the large modification of superconducting critical current by electric fields reported in Refs. 15, it then appears somewhat less likely these results can be understood with electrostatic effects of the type discussed above. At metallic densities ∆/µ ∼ 10 −4 , electrostatic effects in the model here, even at a sharp DOS feature, likely can only give |δf ns /f ns,3D | 10 −2 , which is too small to cause large measurable effects. It appears unlikely this is easily rectified by lifting some of the approximations we made. This is simply a manifestation of the "Anderson theorem": 44 the amplitude of conventional superconductivity is insensitive to time-reversal symmetric perturbations, and suppressing it requires perturbations large compared to µ, which are usually not achievable in the metallic regime below the electrical breakdown field. Also, as the linear gate effect generally should dominate nonlinearities, whether superconductivity is suppressed or enhanced depends on the sign of the electric field, quite unlike in Refs. 15. Previously, reduction in the critical current by an applied field was attributed to modification of vortex pinning. 14,29 In Ref. 15 effects appear also in aluminum strips with lateral size ξ, making this explanation less favorable. In the clean-limit model here, it also appears unlikely the phase slip rates would be significantly affected.
In summary, we considered effects of electrostatic fields on superconductivity self-consistently within a BCS model, connected them to questions about electrostatic energy, and commented on their relation to recent experiments. We obtain results for the size and external electric field modulation of superconductivity, and contrast results to a surface doping model. Expanding about this mean field solution, considering electric field effects on phase slips and phase fluctuations is possible. Experimentally, the effects are best visible in low charge density systems, e.g. semiconductor hybrid structures. state spectral function for a potential well is where k p = πp L , E p = k 2 p 2m . Then we have, which can be evaluated. Here, the terms p = q imply the limit E p → E q .