Effect of layer thickness on the superconducting properties in ultrathin Pb films

The effect of film thickness on the properties of superconducting states has been widely investigated in two-dimensional (2D) systems for several decades [1–3]. Those conducting research are motivated by the unique properties [4–7] and prospects for the potential technological application of these materials, for example in modern electronic and opto-electronic devices [8]. It has been found that a few intriguing aspects and novel physical phenomena exist in 2D superconductors. One of them is a transition from a superconducting to an insulating state with increasing disorder or magnetic field [9–13]. Additionally, quantum size effects such as oscillations of the critical temperature (T_C) and other superconducting parameters with the thickness have stimulated considerable theoretical and experimental interest [14–21].

Recently, using scanning tunnelling microscopy measurements, superconductivity with a T_C up to 1.83 K has been observed in one atomic layer of Pb film grown on a bulk-terminated Si(111) substrate [22]. Moreover, the transition temperature as a function of the film thickness shows oscillations; in particular with the film thickness decrease, an increase in the quantum oscillations of T_C was observed. Systematic studies also found that T_C decreases rapidly to very low values at the limit of 2 monolayers (ML) [23]. It is worth emphasizing that lead attracts more attention in superconductivity studies, because bulk Pb is a prototypical strong coupling electron–phonon elemental superconductor with a relatively high critical temperature of 7.19 K, energy gap ($\Delta (0)$) of 1.4 meV [24] and a ratio $2\Delta (0){/k}_{{\rm{B}}}{T}_{{\rm{C}}}$ of 4.52 [25], well above the canonical BCS value of 3.53 [26].

Motivated by the recent experimental progress in this area, we have carried out calculations to explore the thermodynamic properties of the superconducting state of Pb films in the ultrathin regime (5–10 ML). It is interesting to compare the superconducting properties of lead thin films with a bulk lead sample. In this paper we present the results of our calculations performed within the framework of the strong coupling approach.

The paper is organized as follows. Section 2 contains a short outline of the Eliashberg formalism and computational methods. In section 3, we discuss the thermodynamic properties of superconducting Pb(111) films consisting of five to ten monolayers. section 4 summarizes the obtained results.

The Eliashberg theory of superconductivity [27] allows prediction of the thermodynamic properties of conventional superconductors on the quantitative level. This is possible due to the fact that strong electron–phonon coupling and retardation effects implied by the time dependence of the pairing interaction are not omitted unlike in the case of BCS theory [26, 28].

The Eliashberg theory is based on the electron–phonon interaction Hamiltonian predicted by Fröhlich [29]:

$$\begin{eqnarray}&&H\equiv \displaystyle \sum _{{\bf{k}}\sigma }{\bar{\varepsilon }}_{{\bf{k}}}{c}_{{\bf{k}}\sigma }^{\dagger }{c}_{{\bf{k}}\sigma }+\displaystyle \sum _{{\bf{q}}}{\omega }_{{\bf{q}}}{b}_{{\bf{q}}}^{\dagger }{b}_{{\bf{q}}}+\displaystyle \sum _{{\bf{k}}{\bf{q}}\sigma }{g}_{{\bf{k}}}({\bf{q}}){c}_{{\bf{k}}+{\bf{q}}\sigma }^{\dagger }{c}_{{\bf{k}}\sigma }{\phi }_{{\bf{q}}}.\end{eqnarray} \tag{ 1 }$$

Here ${\bar{\varepsilon }}_{{\bf{k}}}\equiv {\varepsilon }_{{\bf{k}}}-\mu$, where ${\varepsilon }_{{\bf{k}}}$ denotes the electronic band dispersion and μ is the chemical potential. Symbols ${c}_{{\bf{k}}\sigma }$ and ${c}_{{\bf{k}}\sigma }^{\dagger }$ denote the electron annihilation and creation operators of momentum ${\bf{k}}$ and spin σ = $\uparrow$ or $\downarrow$. Symbol ${\omega }_{{\bf{q}}}$ stands for the phonon frequencies, ${b}_{{\bf{q}}}$ and ${b}_{{\bf{q}}}^{\dagger }$ are the phonon annihilation and creation operators of modes ${\bf{q}}$. The electron–phonon matrix element is given by ${g}_{{\bf{k}}}({\bf{q}})$ and ${\phi }_{{\bf{q}}}\equiv {b}_{-{\bf{q}}}^{\dagger }+{b}_{{\bf{q}}}$ [29]. Nonlinear phonon–phonon and quadratic electron–phonon coupling can be included in this model by using anharmonic oscillators [30–37].

The Nambu formalism [38] and Greenʼs function theory [39] allow us to determine, directly from the Fröhlich Hamiltonian, the set of two coupled non-linear integral Eliashberg equations for the superconducting order parameter function ${\Delta }_{n}\equiv \Delta \left({\rm{i}}{\omega }_{n}\right)$ and the mass renormalization function ${Z}_{n}\equiv Z\left({\rm{i}}{\omega }_{n}\right)$ [40]. The one-band equations defined on the imaginary frequency axis take the following form:

$$\begin{eqnarray}&&{\Delta }_{n}{Z}_{n}=\displaystyle \frac{\pi }{\beta }\displaystyle \sum _{m=-M}^{M}\displaystyle \frac{K\left({\rm{i}}{\omega }_{n}-{\rm{i}}{\omega }_{m}\right)-{\mu }^{\star }\theta \left({\omega }_{{\rm{c}}}-| {\omega }_{m}| \right)}{\sqrt{{\omega }_{m}^{2}+{\Delta }_{m}^{2}}}{\Delta }_{m}\end{eqnarray} \tag{ 2 }$$

and

$$\begin{eqnarray}&&{Z}_{n}=1+\displaystyle \frac{1}{{\omega }_{n}}\displaystyle \frac{\pi }{\beta }\displaystyle \sum _{m=-M}^{M}\displaystyle \frac{K\left({\rm{i}}{\omega }_{n}-{\rm{i}}{\omega }_{m}\right)}{\sqrt{{\omega }_{m}^{2}+{\Delta }_{m}^{2}}}{\omega }_{m}.\end{eqnarray} \tag{ 3 }$$

Here ${\omega }_{n}\equiv \frac{\pi }{\beta }(2n-1)$, with $n=0,\pm 1,\pm 2,...,\pm M$ are Matsubara frequencies, where M = 1100. The inverse temperature is given by $\beta \equiv 1/\left({k}_{{\rm{B}}}T\right)$, where k_B represents the Boltzmann constant. Symbols θ and ${\omega }_{{\rm{c}}}$ denote the Heaviside function and the cut-off frequency which is set to ten times the maximum phonon frequency (${\Omega }_{\mathrm{max}}=10$ meV). The electron–phonon pairing kernel is defined as:

$$\begin{eqnarray}&&K\left({\rm{i}}{\omega }_{n}-{\rm{i}}{\omega }_{m}\right)\equiv \lambda \displaystyle \frac{{\Omega }_{{\rm{C}}}^{2}}{{\Omega }_{{\rm{C}}}^{2}-{({\rm{i}}{\omega }_{n}-{\rm{i}}{\omega }_{m})}^{2}}.\end{eqnarray} \tag{ 4 }$$

The symbol ${\Omega }_{{\rm{C}}}$ represents the characteristic phonon frequency, which is commonly used to reproduce, in the framework of the Eliashberg formalism, the critical temperature observed in experimental measurements [41]. At this point it should be noted, that the electron–phonon interaction in metals is generally considered to be isotropic. However, the necessity of an anisotropic theory has clearly been demonstrated in two important cases, magnesium diboride MgB₂, and the graphite intercalation compound CaC₆ [42, 43]. It has also been confirmed that bulk Pb is a two-gap superconductor [44–46]. The anisotropic properties of the superconducting state can be examined in the framework of the multi-band Eliashberg formalism [47, 48]. However, in the case of free-standing Pb(111) films the one-band Eliashberg formalism with very good accuracy reproduces the experimental results obtained using scanning tunneling spectroscopy [49]. A central quantity in the Eliashberg theory is the effective electron–phonon spectral function ${\alpha }^{2}F(\Omega )$ which delivers a description of the electron–phonon interaction by means of a spectral density. This function contains all the information about an investigated system needed to determine the thermodynamic properties of a superconducting state. The expression used to Scalculate the Eliashberg functions for free-standing Pb(111) films is given by

$$\begin{eqnarray}\begin{array}{rcl}{\alpha }^{2}F(\omega ) & = & \displaystyle \frac{1}{{\rm{\hslash }}N({E}_{{\rm{F}}})}\displaystyle \sum _{{\bf{q}},\nu }\delta (\omega -{\omega }_{{\bf{q}},\nu })\displaystyle \sum _{{\bf{k}},{\bf{k}},n\prime }\delta ({\epsilon }_{{\bf{k}},n}-{E}_{{\rm{F}}})\\ & & \times {| {g}^{n,n\prime }({\bf{k}},{\bf{q}},\nu )| }^{2}\delta ({\epsilon }_{{\bf{k}}+{\bf{q}},n\prime }-{E}_{{\rm{F}}}),\end{array}\end{eqnarray} \tag{ 5 }$$

where, $N({E}_{{\rm{F}}})$ is the electronic density of states at the Fermi level, $| {g}^{n,n\prime }({\bf{k}},{\bf{q}},\nu )|$ is an electron–phonon matrix element, and ${\epsilon }_{{\bf{k}},n}$ and ${\epsilon }_{{\bf{k}}+{\bf{q}},n\prime }$ are electronic energies [50]. The calculation for each film thickness was carried out using first-principles calculations in the density functional perturbation formalism, where the free-standing films are modelled by a supercell with N-layer Pb(111) films separated by a vacuum gap of $22 {\mathrm \AA}$ [50].

The information about the layer thickness is contained in the shape of the Eliashberg spectral function, therefore, the number of layers directly influences the value of the electron–phonon coupling constant which is given by the integration of the ${\alpha }^{2}F(\Omega )$ function:

$$\begin{eqnarray}&&\lambda =2{\displaystyle \int }_{0}^{{\Omega }_{\mathrm{max}}}\displaystyle \frac{{\alpha }^{2}F(\Omega )}{\Omega }{\rm{d}}\Omega ,\end{eqnarray} \tag{ 6 }$$

If $\lambda \gt 1$, the compound is classified in the group of strong coupling superconductors and should be described using the strong coupling Eliashberg approach. In the case of free-standing Pb(111) films consisting of five to ten atomic layers we have $\lambda \gt 1.5$ [50]. For conventional superconductors the electron–phonon coupling is counteracted by Coulomb repulsion. The depairing Coulomb interactions are modelled by the semiempirical screened Coulomb pseudopotential parameter ${\mu }^{\star }$. In the presented paper the values of λ are taken from paper [50] and ${\mu }^{\star }=0.13$ (the value obtained for bulk Pb) [24].

The Eliashberg theory allows us to calculate the properties of all conventional superconductors once λ and ${\mu }^{\star }$ are known [51]. To be able to calculate the real physical quantities, the imaginary frequency axis solutions, computed from equations (2) and (3), should be implemented as input parameters to the Eliashberg equations in the mixed representation (defined both on the real and imaginary frequency axes) [52, 53]. In order to determine all of the significant thermodynamic properties of ultrathin Pb films, we conducted numerical calculations based on a self-consistent iteration method [54], which has been used successfully in our previous papers [55–58].

The critical temperatures of ultrathin Pb films have been determined from temperature-dependent gap measurements in paper [49]. On this basis, we were able to calculate the values of ${\Omega }_{{\rm{C}}}$ for each investigated thin film layer thickness (5-10 ML). In particular, in the Eliashberg equations we have assumed that $T={T}_{{\rm{C}}}$ and then we have decreased the value of the parameter ${\Omega }_{{\rm{C}}}$ until we have reached the equality ${\Delta }_{m=1}({\Omega }_{{\rm{C}}})=0$. The full dependence of ${\Delta }_{m=1}$ on ${\Omega }_{{\rm{C}}}$ is presented in the main part of figure 1. The inset presents ${\Omega }_{{\rm{C}}}$ as a function of the layer thickness. It is easy to observe that the critical value of ${\Omega }_{{\rm{C}}}$ oscillations is in the range of 3.06 to 3.94 meV.

Zoom In Zoom Out Reset image size

In figure 2, the dependence of the order parameter on the temperature is presented for film layer thicknesses in the range of 5 to 10 ML. Moreover, the above results are supplemented by the results obtained for bulk Pb [59]. It can be observed that the value of ${\Delta }_{m=1}$ strongly decreases with the growth of the temperature and varies depending on the thickness. This fact can be exactly parameterized using the following formula:

$$\begin{eqnarray}&&{\Delta }_{m=1}={\Delta }_{m=1}(0)\sqrt{1-{\left(\displaystyle \frac{T}{{T}_{{\rm{C}}}}\right)}^{\Gamma }},\end{eqnarray} \tag{ 7 }$$

where the parameter Γ, in all cases is equal to 3.7. The values of T_C and ${\Delta }_{m=1}(0)$ are shown in the insets (A) and (B) of figure 2. We notice that the order parameter function allows us to determine the physical value of the energy gap at the Fermi level $\Delta (0)$. For this purpose, we have numerically solved the Eliashberg equations in the mixed representation by using the solutions from the imaginary axis as input parameters. As a result we obtained 1.32, 1.26, 1.26, 1.21, 1.23 and 1.20 meV for 5, 6, 7, 8, 9, and 10 ML. It should be emphasized that a very good consistency of the experimental results [49] with the presented theoretical ones for $\Delta (0)$ is observed. This proves the correctness of the theoretical model and numerical methods used in this paper. Let us note that in the case of bulk Pb, the above value is much higher: $[\Delta (0){]}_{\mathrm{bulk}}^{\mathrm{Pb}}=1.40$ [24, 25].

Zoom In Zoom Out Reset image size

The energy gap function also allows us to calculate the total normalized density of states [40]:

$$\begin{eqnarray}&&\mathrm{NDOS}\equiv \displaystyle \frac{{\mathrm{DOS}}_{S}(\omega )}{{\mathrm{DOS}}_{N}(\omega )}=\mathrm{Re}\left[\displaystyle \frac{\omega }{\sqrt{{\omega }^{2}-{\Delta }^{2}(\omega )}}\right],\end{eqnarray} \tag{ 8 }$$

where the symbols ${\mathrm{DOS}}_{S}$ and ${\mathrm{DOS}}_{N}$ denote the density of states of the superconducting and normal state, respectively. The plot of the total normalized density of states is presented in the lower panel of figure 2.

In contrast to the order parameter, the value of the mass renormalization function (${Z}_{m=1}$) grows with increasing temperature and the maximum of this function is at T_C, as presented in figure 3. The inset shows the dependence of ${Z}_{m=1}({T}_{{\rm{C}}})$ as a function of the layer thickness. We can observe the saturation effect of oscillations for higher values of thickness. From the physical point of view, the ${Z}_{m=1}$ function with a very good approximation reproduces the value of the electron effective mass: ${m}_{{\rm{e}}}^{\star }\simeq {Z}_{m=1}{m}_{{\rm{e}}}$, where the symbol m_e denotes the electron band mass.

Zoom In Zoom Out Reset image size

To explore the influence of film thickness on the thermodynamic properties of the investigated systems, in the first step we calculated the free energy difference between the superconducting and the normal state ($\Delta F\equiv {F}^{{\rm{S}}}-{F}^{{\rm{N}}}$) [40]:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\Delta F}{N({E}_{{\rm{F}}})} & = & -\displaystyle \frac{2\pi }{\beta }\displaystyle \sum _{n=1}^{M}\left(\sqrt{{\omega }_{n}^{2}+{\Delta }_{n}^{2}}-\left|{\omega }_{n}\right|\right)\\ & & \times \left({Z}_{n}^{{\rm{S}}}-{Z}_{n}^{{\rm{N}}}\displaystyle \frac{\left|{\omega }_{n}\right|}{\sqrt{{\omega }_{n}^{2}+{\Delta }_{n}^{2}}}\right),\end{array}\end{eqnarray} \tag{ 9 }$$

where symbols ${Z}_{n}^{{\rm{S}}}$ and ${Z}_{n}^{{\rm{N}}}$ denote the mass renormalization functions for the superconducting and for the normal state, respectively. Then, on the basis of equation (9), the thermodynamic critical field H_C was determined:

$$\begin{eqnarray}&&\displaystyle \frac{{H}_{{\rm{C}}}}{\sqrt{N({E}_{{\rm{F}}})}}=\sqrt{-8\pi \left[\Delta F/N({E}_{{\rm{F}}})\right]}\end{eqnarray} \tag{ 10 }$$

In figure 4, the detailed dependence of the thermodynamic critical field and free energy difference on the temperature for selected values of film layer thickness is plotted.

Zoom In Zoom Out Reset image size

From the second derivative of the $\Delta F$ function the specific heat difference between the superconducting and the normal state ($\Delta C\equiv {C}^{{\rm{S}}}-{C}^{{\rm{N}}}$) was calculated:

$$\begin{eqnarray}&&\displaystyle \frac{\Delta C}{{k}_{{\rm{B}}}N({E}_{{\rm{F}}})}=-\displaystyle \frac{1}{\beta }\displaystyle \frac{{{\rm{d}}}^{2}\left[\Delta F/N({E}_{{\rm{F}}})\right]}{{\rm{d}}{\left({k}_{{\rm{B}}}T\right)}^{2}},\end{eqnarray} \tag{ 11 }$$

where the specific heat in the normal state is defined as: ${C}^{{\rm{N}}}=\gamma /\beta$. Symbol γ denotes the Sommerfeld constant: $\gamma \equiv (2/3){\pi }^{2}(1+\lambda ){k}_{{\rm{B}}}N({E}_{{\rm{F}}})$. In figure 5 we can trace the course of the specific heat function in the superconducting and normal states.

Zoom In Zoom Out Reset image size

The thermodynamic functions determined in this paper allow us to estimate the dimensionless ratios: ${R}_{\Delta }\equiv 2\Delta (0){/k}_{{\rm{B}}}{T}_{{\rm{C}}}$, ${R}_{{\rm{C}}}\equiv \Delta C\left({T}_{{\rm{C}}}\right)/{C}^{{\rm{N}}}\left({T}_{{\rm{C}}}\right)$ and ${R}_{{\rm{H}}}\equiv {T}_{{\rm{C}}}{C}^{{\rm{N}}}\left({T}_{{\rm{C}}}\right)/{H}_{{\rm{C}}}^{2}(0)$ which play a very important role in the theory of superconductivity because they can be determined in experimental measurements and compared with theoretical predictions. In the framework of the BCS theory, ${R}_{\Delta }$, R_C and R_H are universal and equal to 3.35, 1.43 and 0.168, respectively [26, 28, 40]. The results obtained for ultrathin Pb(111) films are presented in table 1. We can see that the values computed in the framework of the Eliashberg theory strongly differ from the predictions of the BCS theory. This situation is caused by the strong coupling and retardation effects which can be written approximately as the ratio ${k}_{{\rm{B}}}{T}_{{\rm{C}}}{/\omega }_{\mathrm{ln}}$, where ${\omega }_{\mathrm{ln}}$ denotes the average phonon frequency which relates to the dynamics of the superconducting state [60]. For lead ${\omega }_{\mathrm{ln}}$ is no the order of 5.5 meV [61]. In the weak coupling BCS limit ${k}_{{\rm{B}}}{T}_{{\rm{C}}}{/\omega }_{\mathrm{ln}}\to 0$, while in the case of the Pb superconductor this parameter exceeds the value of 0.1. It should be noted that the results obtained for thin films are generally higher than the corresponding bulk values; in particular for bulk Pb we have ${[{R}_{\Delta }]}_{\mathrm{bulk}}^{\mathrm{Pb}}=4.53$, ${[{R}_{{\rm{C}}}]}_{\mathrm{bulk}}^{\mathrm{Pb}}=2.63$ and ${[{R}_{{\rm{H}}}]}_{\mathrm{bulk}}^{\mathrm{Pb}}=0.129$ [24, 59].

Table 1.  The values of the dimensionless parameters: ${R}_{\Delta }$, R_C and R_H for selected values of film thickness in ML.

Thickness	5	6	7	8	9	10
${R}_{\Delta }$	5.02	4.62	4.76	4.63	4.66	4.62
RC	3.76	3.21	3.46	3.27	3.32	3.24
RH	0.126	0.134	0.131	0.134	0.133	0.134

The Eliashberg theory has been used with great success to calculate the thermodynamic properties of ultrathin Pb films in the superconducting state. In particular, the measured superconducting critical temperature was reconstructed by using the matching parameter ${\Omega }_{{\rm{C}}}$. Then the energy gap, free energy difference between the superconducting and normal state, thermodynamic critical field and the specific heat were determined for a wide range of film thicknesses (5–10 ML). We observed that for thin films, the thermodynamic parameter exhibits an oscillatory behaviour connected with a quantum size effect. Next, on the basis of the above thermodynamic functions determined for Pb films, it was proven that the values of the dimensionless ratios ${R}_{\Delta }$, R_C and R_H differ significantly from the predictions of the BCS theory and are generally higher than the corresponding bulk values. This is connected with the existence of strong coupling and retardation effects in the investigated systems.

The results, presented in this paper are expected to stimulate further experimental and theoretical exploration or discovery of new 2D superconductors.

The author acknowledges financial support from the Czȩstochowa University of Technology under Grant No. BS/MN-203-303/2015.