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

TL;DR

This paper develops a theoretical framework to understand how static electric fields influence the critical temperature of superconducting films, emphasizing the importance of film thickness and electron-phonon interactions.

Contribution

It introduces a proximity Eliashberg theory model for electrostatic field effects on superconductors, highlighting parameter-free predictions in weak fields and the impact of surface layer thickness.

Findings

Thin superconducting films with increased electron-phonon coupling can have higher critical temperatures.

The theory provides parameter-free predictions in the weak electrostatic field limit.

Optimal enhancement of Tc occurs with very thin, strongly charged superconducting layers.

Abstract

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.

Tables

Table 1. Table 1: Calculated input parameters with DFT and calculated critical temperature with Eliashberg theory without proximity effect.

$x\left(e^{-}/cell\right)$	$\lambda$	${\omega }_{ln}$ $\left(meV\right)$	$N\left(0\right)$ $states/\left(eVspin\right)$	$\Delta E_F$ $\left(meV\right)$	${\mu }^{\ast }$	$T_c$ $\left(K\right)$
0.000	1.5612	4.8431	0.25866	0.00	0.14164	7.2200
0.075	1.5582	4.8432	0.25754	108.42	0.14136	7.2197
0.150	1.6137	4.7176	0.25611	218.77	0.14116	7.3165
0.300	2.0237	4.2175	0.26770	435.07	0.14074	8.2862
0.400	2.5392	3.5668	0.27833	571.62	0.14048	8.9406

Equations

$${\omega }_n{Z}_b(i{\omega }_n)={\omega }_n+\pi T\sum _m{\Lambda }_b^Z(i{\omega }_n,i{\omega }_m)N_b^Z(i{\omega }_m)+$$

$$+{\Gamma }_{\mathrm{b}}{N}_s^Z(i{\omega }_n)$$

$$\displaystyle Z_{b}(i\omega_{n})\Delta_{b}(i\omega_{n})=\pi T\sum_{m}\big{[}\Lambda^{\Delta}_{b}(i\omega_{n},i\omega_{m})-\mu^{*}_{b}(\omega_{c})\big{]}\times$$

$$\times \Theta ({\omega }_c-|{\omega }_m|)N_b^{\Delta }(i{\omega }_m)+{\Gamma }_{\mathrm{b}}{N}_s^{\Delta }(i{\omega }_n)$$

$${\omega }_n{Z}_s(i{\omega }_n)={\omega }_n+\pi T\sum _m{\Lambda }_s^Z(i{\omega }_n,i{\omega }_m)N_s^Z(i{\omega }_m)+$$

$${\Gamma }_{\mathrm{s}}{N}_b^Z(i{\omega }_n)$$

$$\displaystyle Z_{s}(i\omega_{n})\Delta_{s}(i\omega_{n})=\pi T\sum_{m}\big{[}\Lambda^{\Delta}_{s}(i\omega_{n},i\omega_{m})-\mu^{*}_{s}(\omega_{c})\big{]}\times$$

$$\times \Theta ({\omega }_c-|{\omega }_m|)N_s^{\Delta }(i{\omega }_m)+{\Gamma }_{\mathrm{s}}{N}_b^{\Delta }(i{\omega }_n)$$

$${\Lambda }_{s(b)}(i{\omega }_n,i{\omega }_m)=2{\int }_0^{+\infty }d\Omega \Omega {\alpha }_{s(b)}^{2}F(\Omega )/[({\omega }_n-{\omega }_m{)}^2+{\Omega }^2]$$

$${\Gamma }_{s(b)}=\pi |t{|}^{2}A{d}_{b(s)}{N}_{b(s)}(0)$$

$$N_{s(b)}^{\Delta }(i{\omega }_m)={\Delta }_{s(b)}(i{\omega }_m)/\sqrt{{\omega }_m^2+{\Delta }_{s(b)}^2(i{\omega }_m)}$$

$$N_{s(b)}^Z(i{\omega }_m)={\omega }_m/\sqrt{{\omega }_m^2+{\Delta }_{s(b)}^2(i{\omega }_m)}$$

$${\lambda }_{s(b)}=2{\int }_0^{+\infty }d\Omega \frac{{\alpha }_{s(b)}^{2}F(\Omega )}{\Omega }$$

$$ln({\omega }_{ln,s(b)})=\frac{2}{{\lambda }_{s(b)}}{\int }_0^{+\infty }d\Omega ln\Omega \frac{{\alpha }_{s(b)}^{2}F(\Omega )}{\Omega }$$

$$\mu =\frac{1}{4{\pi }^2\hslash v_F}{\int }_0^{2k_F}\frac{V(q)}{\epsilon (q)}qdq$$

$$\mu =\frac{k_{TF}^2}{8k_F^2}ln(1+\frac{4k_F^2}{k_{TF}^2}).$$

$${\mu }_b=\frac{a_b^2}{2}ln(1+\frac{1}{a_b^2}).$$

$$a_s^2=a_b^2(\frac{N_s(0)}{N_b(0)})/(\frac{1+\Delta E_F}{E_{F,b}})$$

$${\mu }_s=\frac{a_s^2}{2}ln(1+\frac{1}{a_s^2}).$$

$${\mu }_s^{\ast }({\omega }_c)=\frac{{\mu }_s}{1+{\mu }_{s}ln((E_{F,b}+\Delta E_F)/{\omega }_c)}$$

Full text

Field effect, Proximity effect, Eliashberg equations

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 order-of-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 FujimotoReview2013 ; UenoReview2014 ; GoldmanReview2014 ; SaitoReview2016 .

EDL gating was first exploited to control the surface electronic properties of relatively low-carrier density systems, where the electric field effect is more readily observable. Field-induced superconductivity was first demonstrated in strontium titanate UenoNatureMater2008 and zirconium nitride chloride YeNatureMater2010 ; SaitoScience2015 , and subsequently on other insulating systems such as perovskites UenoNatureNano2011 and layered transition-metal dichalcogenides YeScience2012 ; JoNanoLett2015 ; LuScience2015 ; ShiSciRep2015 ; YuNatNano2015 ; CostanzoNatNano2016 ; SaitoNatPhys2016 . Significant effort was also invested in the control of the superconducting properties of cuprates bollinger11 ; LengPRL2011 ; LengPRL2012 ; MaruyamaAPL2015 ; JinSciRep2016 ; FeteAPL2016 ; BurlachkovPRB1993 ; GhinovkerPRB1995 ; walter16 , although in this class of materials the mechanism behind the carrier density modulation is still debated walter16 .

More recently however, several experimental studies have been devoted to the exploration of field effect in superconductors libro with a large ($\gtrsim 1\cdot 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 LiNature2016 ; YoshidaAPL2016 ; XiPRL2016 . The thickness and gate voltage dependence of a high-temperature superconducting phase were studied in iron selenide, both in thin-film ShiogaiNaturePhys2015 and thin flake LeiPRL2016 ; HanzawaPNAS2016 forms. The effect of the ultrahigh interface electric fields achievable via EDL gating were also probed in standard BCS superconductors, namely niobium ChoiAPL2014 and niobium nitride PiattiJSNM2016 ; PiattiNbN .

With the exception of the work of Ref.XiPRL2016, on niobium diselenide, all of these studies have been performed on relatively thick samples ($\gtrsim 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 PiattiNbN 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 PiattiNbN . A proper theoretical treatment for field effect on more complex materials, which can be described only by means of the more complete Eliashberg theory carbibastardo ; ummarinorev , 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}\neq 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 carbibastardo .

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 carbibastardo ; ummarinorev with proximity effect to calculate the critical temperature of the system. In this case we have to solve four coupled equations for the gaps $\Delta_{s,b}(i\omega_{n})$ and the renormalization functions $Z_{s,b}(i\omega_{n})$, where $\omega_{n}$ are the Matsubara frequencies. The imaginary-axis equations with proximity effect Mc ; Carbi1 ; Carbi2 ; Carbi3 ; kresin are:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\mu^{*}_{s(b)}$ are the Coulomb pseudopotentials in the surface and in the bulk respectively, $\Theta$ is the Heaviside function and $\omega_{c}$ is a cutoff energy at least three times larger than the maximum phonon energy. Thus we have

[TABLE]

[TABLE]

and thus $\frac{\Gamma_{s}}{\Gamma_{b}}=\frac{d_{b}N_{b}(0)}{d_{s}N_{s}(0)}$,

[TABLE]

and

[TABLE]

where $\alpha^{2}_{s(b)}F(\Omega)$ 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

[TABLE]

and the representative energies as

[TABLE]

The solution of these equations requires eleven input parameters: the two electron-phonon spectral fuctions $\alpha^{2}_{s(b)}F(\Omega)$, the two Coulomb pseudopotentials $\mu^{*}_{s(b)}$, the values of the normal density of states at the Fermi level $N_{s(b)}(0)$, the shift of the Fermi energy $\Delta 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 cross-sectional 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 PiattiNbN . 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 $\alpha^{2}_{b}F(\Omega)$, $\mu^{*}_{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 $\alpha^{2}_{s}F(\Omega)$, $\Delta E_{F}$ and $N_{s}(0)$.

The value of the Coulomb pseudopotential in the surface layer $\mu^{*}_{s}$ can be obtained in the following way: in the Thomas-Fermi theory where the dielectric function is UmmaC60 $\varepsilon(q)=1+\frac{k^{2}_{TF}}{q^{2}}$ and the bare Coulomb pseudopotential $\mu$ is the angular average of the screened electrostatic potential

[TABLE]

Since $V(q)=\frac{4\pi e^{2}}{q^{2}}$ Morel it turns out that

[TABLE]

Hence we write

[TABLE]

with $a_{b}=2k_{TF,b}/k_{F,b}$. Since $a_{b}$ can be calculated by numerically solving Eq. 13, and by remembering that the square of Thomas-Fermi wave number $k_{TF,b(s)}$ is proportional to $N_{b(s)}(0)$, we have

[TABLE]

and thus

[TABLE]

The new Coulomb pseudopotential Morel in the surface layer is thus

[TABLE]

We note that, usually, the effect of electrostatic doping on $\mu^{*}$ is very small and can be neglected. We can quantify the effect on $T_{c}$ of this small modulation of $\mu^{*}$ 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 $\mu^{*}(x=0)=0.14164$, while for the maximum doping Eqs. 12-16 give $\mu^{*}(x=0.4)=0.14048$. If we use $d_{s}=d_{TF}$ 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 $\Delta T_{c}|_{\Delta\mu^{*}}=-0.0002$ K ($\Delta T_{c}|_{\Delta\mu^{*}}/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, $\mu^{*}$ becomes ill-defined as the Thomas-Fermi dielectric function is no longer strictly valid for very large electric fields. Nevertheless, the true dielectric function $\varepsilon(q)$ should still be a function of the ratio $k_{TF}/k_{F}$ pastore , which in the free-electron model is independent on the normal density of states at the Fermi level. Thus, Eq. 9 should still be able to describe the behavior of the system as a first approximation.

III CALCULATION OF $\alpha^{2}_{s}F(\Omega)$, $\Delta E_{F}$ and $N_{s}(0)$

DFT calculations are performed within the mixed-basis pseudopotential method (MBPP) meyer . For lead a norm-conserving relativistic pseudopotential including $5d$ semicore states and partial core corrections is constructed following the prescription of Vanderbilt vanderbilt . This provides both scalar-relativistic and spin-orbit components of the pseudopotential. Spin-orbit coupling (SOC) is then taken into account within each DFT self-consistency cycle (for more details on the SOC implementation see heid10 ). 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) hedin . Brillouin zone (BZ) integrations are performed on regular k-point meshes in conjuction with a Gaussian broadening of $0.2$ eV. For phonons, $16\times 16\times 16$ meshes are used, while for the calculations of density of states and electron-phonon coupling (EPC) even denser $32\times 32\times 32$ meshes are employed.

Phonon properties are calculated with the density-functional perturbation theory baroni ; zein as implemented in the MBPP approach heid99 , which also provides direct access to the electron-phonon coupling matrix elements. The procedure to extract the Eliashberg function is outlined in Ref. heid10 . Dynamical matrices and corresponding EPC matrix elements are obtained on a $16\times 16\times 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\times 40\times 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 heid10 .

Charge induction is simulated by adding an appropriate number of electrons during the DFT self-consistency 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 sklyadneva . 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 Pines ) is unable to describe inhomogeneous distributions caused by the screening of the electric field itself. A more complete approach has been developed in Ref. Brumme, , 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 carbibastardo $T_{c,b}=7.22$ K. The undoped $\alpha^{2}_{b}F(\Omega)$ gives a corresponding electron-phonon coupling $\lambda_{b}=1.5596$. Assuming a cutoff energy $\omega_{c}=60$ meV and a maximum energy $\omega_{max}=70$ meV in the Eliashberg equations, we are thus able to determine the bulk Coulomb pseudopotential to be $\mu^{*}_{b}=0.14164$ to obtain the exact experimental critical temperature $T_{c,b}$.

In Fig. 2a we show the calculated electron-phonon spectral functions $\alpha^{2}F(\Omega)$ 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 $\omega_{ln}$ with increasing doping level. The increase of the carrier density gives rise to two competing effects: the value of $\omega_{ln}$ (i.e. the representative phonon energy) decreases while the value of electron-phonon coupling costant $\lambda$ increases (see Fig. 2b). Since the critical temperature is an increasing function of both $\omega_{ln}$ and $\lambda$, in general this could result in either a net enhancement or suppression of $T_{c}$, depending 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 $\lambda$ and $\omega_{ln}$ concurrently. In the case of lead the contribution from the increase of $\lambda$ is dominant over that from the reduction of $\omega_{ln}$, giving rise to a net increase of the superconducting critical temperature (as we report in Fig. 2c). In addition, in Table 1 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 HirschPRB2004 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_{TF}$, which for lead can be estimated to be 0.15 nm Basavaiah . However, we have recently shown PiattiNbN 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 PiattiNbN . 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 $d_{s}=d_{TF}[1+m\Theta(x-x_{0})]$, 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}$ m2. Note that the case $m=0$ of course corresponds to the case where the material 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_{TF}$), 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 Pratappone .

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_{TF}=0.3$ nm) and $m=4$ ($d_{s}$ is allowed to expand up to $5d_{TF}=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}\simeq 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, $\Delta n_{2D}$, required to screen it, i.e. $\Delta n_{2D}=\int_{0}^{d_{s}}\Delta 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 Brumme ; however, as a first approximation we can consider this distribution to be constant, obtaining an effective doping level per unit volume simply as $x=\Delta 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_{TF}[1+m\Theta(x-x_{0})]$ are no longer independent as in the previous case. Moreover, according to our recent experimental findings on niobium nitride PiattiNbN , $d_{s}$ is a monotonically increasing function of $\Delta n_{2D}$. We include this behavior in our calculations in the following way: Once the maximum doping level $x_{0}$ is selected, $m=m(\Delta n_{2D})$ is automatically determined by the requirement $d_{s}(\Delta n_{2D})=\Delta n_{2D}\cdot x_{0}$ for any $x>x_{0}$. Fig. 5a shows the resulting dependence of the doping per unit volume $x$ and surface layer thickness $d_{s}$ on the induced carrier density per unit surface $\Delta n_{2D}$, for two different values of the maximum doping level $x_{0}=0.3$ and $x_{0}=0.4$ e-/unit cell. When $\Delta n_{2D}$ is small enough so that $x<x_{0}$, the Thomas-Fermi screening holds, $d_{s}=d_{TF}$ is constant and $x$ linearly increases with $\Delta n_{2D}$. As soon as $\Delta n_{2D}$ becomes large enough that $x=x_{0}$ is constant ($\Delta n_{2D}(x_{0})$), the Thomas-Fermi screening is no longer valid and $d_{s}$ increases linearly with $\Delta 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 $\Delta n_{2D}$. When $\Delta n_{2D}\lesssim\Delta 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 $\alpha^{2}F(\Omega)$. 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 $\Delta n_{2D}>\Delta n_{2D}(x_{0})$, the surface properties ($\alpha^{2}F(\Omega)$) 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 $\Delta 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 $\Delta n_{2D}>\Delta 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 $\Delta n_{2D}$. In this intermediate region, the modulation of $\alpha^{2}F(\Omega)$ 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 fundamental 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 $\Delta E_{F}$, density of states $N(0)$, and electron-phonon spectral function $\alpha^{2}F(\Omega)$). The Coulomb pseudopotential $\mu_{*}$ 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_{TF}$) 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 estimations 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\lesssim 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-effect-device architecture.

Finally, our calculations indicate that sizable $T_{c}$ enhancements of the order of $\sim 0.5$ K should be achievable in thin films of a standard strong-coupling superconductor such as lead, for easily realizable thicknesses of $\sim 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.

Acknowledgements.

The work of G. A. U. was supported by the Competitiveness Program of NRNU MEPhI.