Proximity Eliashberg theory of electrostatic field-effect-doping in superconducting films

We calculate the effect of a static electric field on the critical temperature of a s-wave one band superconductor in the framework of proximity effect Eliashberg theory. In the weak electrostatic field limit the theory has no free parameters while, in general, the only free parameter is the thickness of the surface layer where the electric field acts. We conclude that the best situation for increasing the critical temperature is to have a very thin film of a superconducting material with a strong increase of electron-phonon (boson) constant upon charging.


I. INTRODUCTION
In the last decade, electric double layer (EDL) gating has come to the forefront of solid state physics due to its capability to tune the surface carrier density of a wide range of different materials well beyond the limits imposed by solid-gate field-effect devices. The orderof-magnitude enhancement in the gate electric field allows this technique to reach doping levels comparable to those of standard chemical doping. This is possible due to the extremely large specific capacitance of the EDL that builds up at the interface between the electrolyte and the material under study 1-4 .
More recently however, several experimental studies have been devoted to the exploration of field effect in superconductors 25 with a large ( 1 · 10 22 cm −3 ) native carrier density. The interplay between two different ground states, namely superconductivity and charge density waves, was explored in titanium and niobium diselenides [26][27][28] . The thickness and gate voltage dependence of a high-temperature superconducting phase were studied in iron selenide, both in thin-film 29 and thin flake 30,31 forms. The effect of the ultrahigh interface electric fields achievable via EDL gating were also probed in standard BCS superconductors, namely niobium 32 and niobium nitride 33,34 .
With the exception of the work of Ref. 28 on niobium diselenide, all of these studies have been performed on relatively thick samples ( 10 nm) with a thickness larger than the electrostatic screening length. These systems are thus expected to develop a strong dependence of their electronic properties along the z direction (z being perpendicular to the sample surface). As a first approximation, this dependence can be conceptualized by schematizing the system as the parallel of a surface layer (where the carrier density is modified by the electric field) and an unperturbed bulk. The two electronic systems can be expected to couple via superconducting proximity effect, resulting in a complicated response to the applied electric field that goes well beyond a simple modification of the superconducting properties of the surface layer alone 34 and is strongly dependent on both the electrostatic screening length and the total thickness of the film.
So far, the only quantitative assessment of this phenomenon has been reported in the framework of the strong-coupling limit of the BCS theory of superconductivity 34 . A proper theoretical treatment for field effect on more complex materials, which can be described only by means of the more complete Eliashberg theory 35,36 , is lacking. Given the rising interest in the control of the properties of superconductors by means of surface electric fields, the development of such a theoretical approach would be very convenient not only to quantitatively describe the results of future experiments, but also to determine beforehand the experimental conditions (e.g., device thickness) most suitable for an optimal control of the superconducting order via electric fields.
In this work, we use the Eliashberg theory of proximity effect to describe a junction composed by the perturbed surface layer (T c = T c,s ), where the carrier density is modulated (with a doping level per unitary cell x), and the underlying unperturbed bulk (T c = T c,b ). Here s and b indicate "surface" and "bulk" respectively (see Fig. 1). Under the application of an electric field, T c,s = T c,b and the material behaves like a junction between a superconductor and a normal metal in the temperature range bounded by T c,s and T c,b . If the application of the electric field increases (decreases) T c,s , then the surface layer will be the superconductor (normal metal) and the bulk will be the normal metal (superconductor). We perform the calculation for lead since all the input parameters of the theory are well-known in the literature for this strong-coupling superconductor 35 .

II. MODEL: PROXIMITY ELIASHBERG EQUATIONS
In general, a proximity effect at a superconductor/normal metal junction is observed as the opening of a finite superconducting gap in the normal metal together with its reduction in a thin region of the superconductor close to the junction. In our model we use the one band s-wave Eliashberg equations 35,36 with proximity effect to calculate the critical temperature of the system. In this case we have to solve four coupled equations for the gaps ∆ s,b (iω n ) and the renormalization functions Z s,b (iω n ), where ω n are the Matsubara frequencies. The imaginaryaxis equations with proximity effect 37-41 are: where µ * s(b) are the Coulomb pseudopotentials in the surface and in the bulk respectively, Θ is the Heaviside function and ω c is a cutoff energy at least three times larger than the maximum phonon energy. Thus we have and thus Γs and where α 2 s(b) F (Ω) are the electron-phonon spectral functions, A is the junction cross-sectional area, |t| 2 the transmission matrix equal to one in our case because the material is the same, d s and d b are the surface and bulk layer thicknesses respectively, such that (d s + d b = d where d is the total film thickness) and N s(b) (0) are the densities of states at the Fermi level E F,s(b) for the surface and bulk material. The electron-phonon coupling constants are defined as and the representative energies as The solution of these equations requires eleven input parameters: the two electron-phonon spectral fuctions α 2 s(b) F (Ω), the two Coulomb pseudopotentials µ * s(b) , the values of the normal density of states at the Fermi level N s(b) (0), the shift of the Fermi energy ∆E F = E F,s −E F,b that enters in the calculation of the surface Coulomb pseudopotential (as shown later), the value of the surface layer d s , the film thickness d and the junction crosssectional area A. The values of d and A are experimental data. The exact value of d s , in particular in the case of very strong electric fields at the surface of a thin film, is in general difficult to determine a priori 34 . Thus, we leave it as a free parameter of the model, and we perform our calculations for different reasonable estimations of its value.
Typically, the bulk values of α 2 b F (Ω), µ * b , N b (0) and E F,b are known and can be found in the literature. Thus, we assume that we need to determine only their surface values. In the next Section, we will use density functional theory (DFT) to calculate α 2 s F (Ω), ∆E F and N s (0). The value of the Coulomb pseudopotential in the surface layer µ * s can be obtained in the following way: in the Thomas-Fermi theory where the dielectric function is 42 ε(q) = 1 + k 2 T F q 2 and the bare Coulomb pseudopotential µ is the angular average of the screened electrostatic potential Hence we write Since a b can be calculated by numerically solving Eq. 13, and by remembering that the square of Thomas-Fermi wave number and thus The new Coulomb pseudopotential 43 in the surface layer is thus We note that, usually, the effect of electrostatic doping on µ * is very small and can be neglected. We can quantify the effect on T c of this small modulation of µ * by computing it in the case of maximum doping x = 0.40 and very thin film (d = 5 nm), i.e. when the effect is largest. As discussed in the next Section, the unperturbed Coulomb pseudopotential is µ * (x = 0) = 0.14164, while for the maximum doping Eqs. 12-16 give µ * (x = 0.4) = 0.14048. If we use d s = d T F we find respectively T c = 7.3770 K for the bulk (unperturbed) value of the Coulomb pseudopotential and T c = 7.3768 K for the surface value of the Coulomb pseudopotential. Thus, if we consider the Coulomb pseudopotential to be doping-independent we underestimate the critical temperature of a ∆T c | ∆µ * = −0.0002 K (∆T c | ∆µ * /T c = 0.0027 percent). However, a possible critical situation can appear when the applied electric field is very strong and the Thomas-Fermi approximation does not hold anymore. In such a case, µ * becomes ill-defined as the Thomas-Fermi dielectric function is no longer strictly valid for very large electric fields. Nevertheless, the true dielectric function ε(q) should still be a function of the ratio k T F /k F 44 , which in the free-electron model is independent on the normal density of states at the Fermi level. Thus, Eq. 11 should still be able to describe the behavior of the system as a first approximation. DFT calculations are performed within the mixedbasis pseudopotential method (MBPP) 47 . For lead a norm-conserving relativistic pseudopotential including 5d semicore states and partial core corrections is constructed following the prescription of Vanderbilt 48 . This provides both scalar-relativistic and spin-orbit components of the pseudopotential. Spin-orbit coupling (SOC) is then taken into account within each DFT selfconsistency cycle (for more details on the SOC implementation see 51 ). The MBPP approach utilizes a combination of local functions and plane waves for the basis set expansion of the valence states, which reduces the size of the basis set significantly. One d type local function is added at each lead atomic site to efficiently treat the deep 5d potential. Sufficient convergence is then achieved with a cutoff energy of 20 Ry for the plane waves. The exchange-correlation it treated in the local density approximation (LDA) 49 . Brillouin zone (BZ) integrations are performed on regular k-point meshes in conjuction with a Gaussian broadening of 0.2 eV. For phonons, 16 × 16 × 16 meshes are used, while for the calculations of density of states and electron-phonon coupling (EPC) even denser 32 × 32 × 32 meshes are employed.
Phonon properties are calculated with the densityfunctional perturbation theory 52,53 as implemented in the MBPP approach 50 , which also provides direct access to the electron-phonon coupling matrix elements. The procedure to extract the Eliashberg function is outlined in Ref. 51 . Dynamical matrices and corresponding EPC matrix elements are obtained on a 16 × 16 × 16 phonon mesh. These quantities are then interpolated using standard Fourier techniques to the whole BZ, and the Eliashberg functions are calculated by integration over the BZ using the tetrahedron method on a 40 × 40 × 40 mesh. SOC is consistently taken into account in all calculations including lattice dynamical and EPC properties. It is well known from previous work, that SOC is necessary for a correct quantitative description of both the phonon anomalies and EPC of undoped bulk lead 51 .
Charge induction is simulated by adding an appropriate number of electrons during the DFT selfconsistency cycle, compensated by a homogeneous background charge to retain overall charge neutrality. Electronic structure properties and the Eliashberg function are calculated for face centred cubic (fcc) lead with the lattice constant a = 4.89Å as obtained by optimization for the undoped case. For doping levels considered here, we found that to a good approximation charge induction does not change the band structure but merely results in a shift of E F . In a previous study, the variation of the EPC was studied as function of the averaging energy 54 . The present approach goes beyond this analysis by taking into account explicitly the effect of charge induction on the screening properties, which modifies both the phonon spectrum and the EPC matrix elements.
Finally we point out that, since this DFT approach simulates the effect of the electric field by adding extra charge carriers to the system together with a uniform compensating counter-charge (Jellium model 55 ) is unable to describe inhomogeneous distributions caused by the screening of the electric field itself. A more complete approach has been developed in Ref. 58, and requires the self-consistent solution of the Poisson equation; however, this method is currently unable to compute the phonon spectrum of the gated material, making it unsuitable for the application of the proximity Eliashberg formalism.

IV. RESULTS AND DISCUSSION
We start our calculations by fixing the input parameters for bulk lead according to the established literature. We set T c,b to its experimental value 35 T c,b = 7.22 K. The undoped α 2 b F (Ω) gives a corresponding electron-phonon coupling λ b = 1.5596. Assuming a cutoff energy ω c = 60 meV and a maximum energy ω max = 70 meV in the Eliashberg equations, we are thus able to determine the bulk Coulomb pseudopotential to be µ * b = 0.14164 to obtain the exact experimental critical temperature T c,b .
In Fig. 2a we show the calculated electron-phonon spectral functions α 2 F (Ω) resulting at the increase of the doping level x. Specifically, we plot the curves corresponding to x = 0.000, 0.075, 0.150, 0.300, 0.400 e − /unit cell. We calculate the spectral functions until x = 0.4 e − /cell because for larger values of doping an instability emerges in the calculation processes. We can see the phonon softening evidenced by a reduction of ω ln with increasing doping level. The increase of the carrier density gives rise to two competing effects: the value of ω ln (i.e. the representative phonon energy) decreases while the value of electron-phonon coupling costant λ increases (see Fig. 2b). Since the critical temperature is an increasing function of both ω ln and λ, in general this could result in either a net enhancement or suppression of T c , depending x(e − /cell) λ ω ln (meV ) N (0) states/(eV spin) ∆EF (meV ) µ * on which of the two contributions is dominant. Consequently the ideal situation for obtaining largest critical temperature in an electric field doped material is to have a strong increase of λ and ω ln concurrently. In the case of lead the contribution from the increase of λ is dominant over that from the reduction of ω ln , giving rise to a net increase of the superconducting critical temperature (as we report in Fig. 2c). In addition, in Table I we summarize all the input parameters of the proximity Eliashberg equations as obtained from the DFT calculations.
Having determined the response of the superconducting properties of a homogeneous lead film to a modulation of its carrier density, we can now consider the behavior of the more realistic junction between the perturbed surface layer and the unperturbed bulk. In order to do so, however, it is now mandatory to select a value for the thickness of the perturbed surface layer. Close to T c , the superfluid density is small 59 and the screening is dominated by unpaired electrons. Thus, a very rough approximation would be to set d s to the Thomas-Fermi screening length d T F , which for lead can be estimated to be 0.15 nm 56 . However, we have recently shown 34 that this assumption might not be satisfactory in the presence of the very large electric fields that build up in the electric double layer. Indeed, our experimental findings on niobium nitride indicated that the screening length increases for very large doping values 34 . However, it is reasonable to assume the exact entity of this increase to be specific to each material. Thus, while the qualitative behavior can be expected to be general, the exact values of d s determined for niobium nitride cannot be directly applied to lead.
In order not to lose the generality of our approach, we calculate the behavior of our system for three different choices of the behavior of d s . We start by expressing where m is a dimensionless parameter indicating how much d s expands beyond the Thomas-Fermi value for large doping levels, and x 0 is the specific doping value upon which this increase in d s takes place. By selecting x 0 = 0.2, we allow the upper half of our doping values to go beyond the Thomas-Fermi approximation. We then perform proximity-coupled Eliashberg calculations for m = 0, 1, 4 and five different film thicknesses d = 5, 10, 20, 30, 40 nm, always assuming the junction area to be A = 10 −7 m 2 . Note that the case m = 0 of course corresponds to the case where the ma-terial satisfies the Thomas-Fermi model for any value of doping: in this case, the model has no free parameters.
In Fig. 3 we plot the evolution of T c upon increasing electron doping and assuming that the Thomas-Fermi model always holds (m = 0 and d s = d T F ), for different values of film thickness. The calculations show that the qualitative increase in T c with increasing doping level that we observed in the homogeneous case is retained also in proximized films of any thinckess (see Fig. 3a). However, the presence of a coupling between surface and bulk induced by the proximity effect gives rise to a key difference with respect to the homogeneous case, namely, a strong dependence of T c on film thickness in the doped films. Indeed, the magnitude of the T c shift with respect to the homogeneous case is heavily suppressed already in films as thin as 5 nm. This behavior is best seen in Fig.  3b, where we plot the same data as a function of the total film thickness for all doping levels. As we can see the increase of critical temperature drops dramatically with increasing film thickness. We have not calculated the critical temperature for monolayer films since the approximations of the model would no longer apply in this case: in particular the unperturbed electron-phonon spectral function would have been different from the bulk-like one we employed in our calculations 57 .
We now consider the effect of the different degrees of confinement for the induced charge carriers at the surface of the films. We do so by allowing the perturbed surface layer to spread further in the depth of the film for large electron doping, i.e. by increasing the m parameter in the definition of d s . In Fig. 4 we plot the evolution of T c with increasing electron doping and for different film thicknesses, in the two cases m = 1 (d s is allowed to expand up to 2d T F = 0.3 nm) and m = 4 (d s is allowed to expand up to 5d T F = 0.75 nm). We can first observe how a different value of d s does not change the qualitative behavior of the films. The evolution of T c with increasing electron doping is still comparable to both the homogeneous case and the proximized films in the Thomas-Fermi limit. The suppression of the T c increase with increasing film thickness is also similar to the latter case. However, the magnitude of the T c shift for the same values of film thickness and doping level per unit cell is clearly the more enhanced the larger the value of d s . This is to be expected, as larger values of d s increase the fraction of the film that is perturbed by the application of the electric field and reduce the T c shift dampening operated by the proximity effect. In principle, for values of m large enough (or film thickness d small enough) one could reach the limit value d s ≃ d and recover the homogeneous case where the T c shift is maximum.
All the calculations we performed so far assumed that one could directly control the induced carrier density per unit volume, x, in the surface layer, without an explicit upper limit. However, this is not an experimentally achievable goal in a field-effect device architecture. In this class of devices, the polarization of the gate electrode allows one to tune the electric field at the interface and thus the induced carrier density per unit surface, ∆n 2D , required to screen it, i.e. ∆n 2D = ds 0 ∆n 3D dz within our model is distributed within a layer of thickness d s . In general, the determination of the exact depth profile of this distribution requires the self-consistent solution of the Poisson equation 58 ; however, as a first approximation we can consider this distribution to be constant, obtaining an effective doping level per unit volume simply as x = ∆n 2D /d s . This procedure allows one to employ the same DFT-Eliashberg formalism we developed before in order to simulate a field effect experiment on a superconducting thin film.
In addition, in the previous calculations we supposed that d s can only take on two values as a function of x, depending on the threshold value x 0 . When we consider the field-effect architecture, however, the parameters m and x 0 in the expression d s = d T F [1 + mΘ(x − x 0 )] are no longer independent as in the previous case. Moreover, according to our recent experimental findings on niobium nitride 34 , d s is a monotonically increasing function of ∆n 2D . We include this behavior in our calculations in the following way: Once the maximum doping level x 0 is selected, m = m(∆n 2D ) is automatically determined by the requirement d s (∆n 2D ) = ∆n 2D · x 0 for any x > x 0 . Fig. 5a shows the resulting dependence of the doping per unit volume x and surface layer thick- ness d s on the induced carrier density per unit surface ∆n 2D , for two different values of the maximum doping level x 0 = 0.3 and x 0 = 0.4 e − /unit cell. When ∆n 2D is small enough so that x < x 0 , the Thomas-Fermi screening holds, d s = d T F is constant and x linearly increases with ∆n 2D . As soon as ∆n 2D becomes large enough that x = x 0 is constant (∆n 2D (x 0 )), the Thomas-Fermi screening is no longer valid and d s increases linearly with ∆n 2D .
In Fig. 5b and 5c we plot the resulting modulation of T c for five different film thicknesses in the cases x 0 = 0.4 and 0.3 e − /unit cell respectively. In both cases we can readily distinguish between two regimes of ∆n 2D . When ∆n 2D ∆n 2D (x 0 ), Thomas-Fermi screening holds and we reproduce the behavior we observed in Fig. 3a. In this regime, the induced carrier density directly modulates x and thus the electron-phonon spectral function α 2 F (Ω). The T c modulation is thus a result of a direct modification of the material properties at the surface, with proximity effect simply operating a "smoothing" the larger the value of the film thickness. On the other hand, when ∆n 2D > ∆n 2D (x 0 ), the surface properties (α 2 F (Ω)) are no longer modified by the extra charge carriers, and the further modulation of T c originates entirely from the proximity effect as determined by the increase in d s .
We can also compare the T c shifts for different maximum doping levels x 0 . Fig. 5d shows the difference between the T c corresponding to x 0 = 0.4 and 0.3 e − /unit cell as a function of the total film thickness, for different values of ∆n 2D . We can clearly see how T c is always larger for the films with larger x 0 , for any value of film thickness, even if the associated values of d s are always smaller. This indicates that the maximum achievable value of x 0 is dominant with respect to the increase of d s to determine the final value of T c , also in the doping regime ∆n 2D > ∆n 2D (x 0 ).
Of course, in a real sample we don't expect the transition between the two regimes to be so clear-cut, as the saturation of x to x 0 would occur over a finite range of ∆n 2D . In this intermediate region, the modulation of α 2 F (Ω) and d s would both contribute in a comparable way to the final value of T c in the film. We stress, however, that in both regimes the proximity effect is funda-mental in determining the T c of the gated film. We also note that the proximized Eliashberg equations are able to account for a non-uniform scaling of the T c shift for different values of film thickness, unlike the models that use approximated analytical equations for T c .

V. CONCLUSIONS
In this work, we have developed a general method for the theoretical simulation of field-effect-doping in superconducting thin films of arbitrary thickness, and we have benchmarked it on lead as a standard strong-coupling superconductor. Our method relies on ab-initio DFT calculations to compute how the increasing doping level x per unit volume modifies the structural and electronic properties of the material (shift of Fermi level ∆E F , density of states N (0), and electron-phonon spectral function α 2 F (Ω)). The Coulomb pseudopotential µ * is determined by simple calculations from some of these parameters. The properties of the pristine thin film (critical temperature T c , device area A and total film thickness d) can be obtained either from the literature or experimentally from standard transport measurements. For doping values where the Thomas-Fermi theory of screening is satisfied, the perturbed surface layer thickness is constant (d s = d T F ) and the theory has no free parameters.
Once all the input parameters are known, our method allows to compute the transition temperature T c for arbitrary values of film thickness d and electron doping in the surface layer x by solving the proximity-coupled Eliashberg equations in the surface layer and unperturbed bulk. On the other hand, if no reliable estima-tions of the surface layer thickness d s are available, our method allows one to determine d s (x) by reproducing the experimentally-measured T c (x). This allows to probe deviations from the standard Thomas-Fermi theory of screening in the presence of very large interface electric fields.
We also show how, even in the case where the Thomas-Fermi approximation breaks down and the doping level x can no longer be increased, the transition temperature T c of a thin film can still be indirectly modulated by the electric field by changing the surface layer thickness d s . For what concerns artificial enhancements of T c in superconducting thin films, we conclude that very thin films (d d s , in order to minimize the smoothing operated by the proximity effect) of a superconductive material characterized by a strong increase of the electron-phonon (boson) coupling upon changing its carrier density are required to optimize the effectiveness of the field-effectdevice architecture.
Finally, our calculations indicate that sizable T c enhancements of the order of ∼ 0.5 K should be achievable in thin films of a standard strong-coupling superconductor such as lead, for easily realizable thicknesses of ∼ 10 nm and doping levels routinely induced via EDL gating in metallic systems. These features may open the possibility for superconducting switchable devices and electrostatically reconfigurable superconducting circuits above liquid helium temperature.