Surface lattice vibration and electron–phonon interaction in Pb(111) ultrathin superconducting films from first principles: both with and without Si substrate

The lattice vibration and electron–phonon (EP) coupling of three types of Pb ultrathin films were calculated from first principles. For free-standing films, our calculations showed that softening and stiffening of the surface phonons coexist. The obtained EP coupling λ and the superconducting transition temperature TC are smaller than their bulk values. For (1×1)-Pb/Si(111) films, the obtained λ and TC values in hexagonally close-packed (hcp) stacking are larger than those in face-centered-cubic (fcc) stacking. The largest TC value is obtained for N=4 ML (monolayers) in hcp stacking, which is very much consistent with the experimental study. Hence, our calculations imply that epitaxial ultrathin Pb films on the surface of Si(111) may possess an hcp rather than an fcc structure. Finally, the surface reconstruction and EP coupling of the striped incommensurate phase with 4/3 ML Pb on Si were studied, and the estimated TC in this phase is very close to its experimental value. Our calculations show clearly that the dynamical properties and EP coupling of ultrathin films depend strongly on the interlayer relaxation and the bonding environment between overlayer and substrate atoms at the interface.


Introduction
With the development of advanced technology in material synthesis, ultrathin metal films even down to a thickness of a few atomic layers [1][2][3] have been grown epitaxially on a substrate. For ultrathin Pb films on Si, Guo et al [4] observed quantum oscillations of the superconducting order parameter as a function of film thickness. Eom et al [5] discovered that superconductivity remained surprisingly robust for films as thin as 5 ML (monolayers). More recently, superconductivity for even 1 and 2 ML Pb films grown epitaxially on Si(111) substrate was reported. Qin et al [6] found that there were two different types of 2 ML Pb films, which differed in subtle atomic reconstructions with superconducting transition temperatures of 4.8 and 3.8 K, respectively. Zhang et al [7] reported that the one-atomic-layer film showed a superconducting transition temperature of 1.83 K for Pb coverage of 4/3 ML and 1.52 K for Pb coverage of 6/5 ML. Observations of superconductivity at this ultimate limit will stimulate the investigation of two-dimensional (2D) superconductivity because of potential applications.
Yndurain and Jigato [8] performed a theoretical first-principles calculation for Pb(111) slabs. They found that the frequency of localized surface phonons and the deformation potential for the lowest unoccupied and highest occupied quantum well states had oscillatory behavior with a bilayer periodicity. However, their calculations were limited to the phonons at K = 0. More recently, Noffsinger and Cohen [9] reported the results of their first-principles calculations on the electron-phonon (EP) coupling for free-standing 2-6 ML Pb films and found that the transition temperature of Pb was suppressed by the presence of stiffened surface phonons. However, it is believed that the substrate is also a key factor in influencing the superconductivity of ultrathin films. To elucidate the role of the substrate, in this paper we perform slab calculations from density-functional perturbation theory to study surface phonons and EP interactions for both free-standing and Si-substrate-supporting Pb(111) films with thickness 1-5 ML. The (1 × 1) and √ 3× √ 3 structures of Pb on Si are studied.

Computational method
First-principles calculations based on density-functional theory in the local-density approximation (LDA) were performed using the PWSCF code of the Quantum-ESPRESSO distribution [10]. The all-electron potentials are replaced by the norm-conserving pseudopotentials as generated in the scheme of Troullier-Martins [11] with Pb (5d 10 6s 2 6p 2 ) and Si (3s 2 3p 2 ) electrons treated as valence states. The cut-off energy of 40 Ryd is used for the expansion of the electronic wave function in plane waves. The films are modeled using the slabs, which are separated by about seven atomic layers of vacuum. For the electronic structure calculations, the Brillouin zone integrations of the (1 × 1) surface unit cell are performed with a (12 × 12) grid by using the first-order Hermite-Gaussian smearing technique [12]. For the surface unit cell other than (1 × 1), equivalent k points have been used. Within the framework of the linear response theory, the dynamical matrices and the EP interaction coefficients are calculated for a (6 × 6) grid of special q points in the 2D irreducible Brillouin zone. The dense (24 × 24) grid is used in the Brillouin zone integrations in order to produce accurate EP interaction matrices.

Free-standing Pb(111) films
We first study the free-standing Pb(111) films with thickness 2-5 ML. The experimental lattice parameter a 0 = 4.95 Å is adopted and the in-plane lattice constant of the (111) film is √ 2a 0 /2. Starting from an ideal surface, the relaxation is performed by minimizing the total energy with respect to the atomic positions in the slab. When the force on each atom is smaller than 10 −4 (Ryd a.u. −1 ), we assume that the atoms are in their equilibrium positions. The results show that the contraction of the first to second layer spacing is unusually large. The calculated interlayer relaxations d 1,2 are listed in table 1. Here d 1,2 is given by d 1,2 = 100(d 1,2 − d 0 )/d 0 . d 1,2 is the relaxed interlayer distance between the top two layers and d 0 = √ 3a 0 /3 is the unrelaxed interlayer distance for (111) slabs. We can see that the calculated percentage of contraction shows a minimum for N = 4 ML. However, the extent of the contraction in our LDA calculations is larger than that for the previous general gradient approximation (GGA) calculations [13]. This may be due to the well-known underestimation of the volume of the unit cell in the LDA calculation.
After relaxation, we calculate the lattice dynamics and EP interaction for ultrathin Pb films. The EP coupling constant λ and the logarithmically averaged frequency ω ln are directly obtained by evaluating [14]  and where the EP spectral function α 2 F(ω) can be determined self-consistently by the linear response theory. The calculated phonon densities of states (DOS) F(ω) for Pb films with N = 2-5 MLs as well as for the bulk are shown in the left panel of figure 1. As a result of the symmetry breaking at the surface, the degeneracy of some modes is lifted, leading to rich structures in the curves. Furthermore, we note that the extending range of the phonon spectrum for ultrathin films is much larger than that for bulk Pb. Compared to the bulk phonon DOS, some DOS peaks for Pb films shift towards lower frequency on the low-frequency side; meanwhile, some others also shift towards higher frequency on the high-frequency side. Our results show clearly that softening and stiffening coexist in Pb ultrathin films. The softening of the surface phonon is due to the lack of bonding atoms in the vacuum side of the surface. In contrast, the stiffening of the surface phonon is the result of a stronger Pb-Pb bond arising from the contraction of the topmost interlayer distance. Calculations by Noffsinger and Cohen [9] showed that many soft phonon modes resulted in an unstable structure for 3 and 4 ML Pb films. However, the 5 softening of the phonon does not lead to dynamical unstability according to our calculations. This difference may stem from different lattice constants used in calculations. In this paper, we use the experimental lattice constant, while Noffsinger and Cohen used their optimized value [9]. The calculated EP spectral functions α 2 F(ω) are shown in the right panel of figure 1. In all cases, the shapes of F(ω) and α 2 F(ω) are very similar.
The superconducting temperature T C can be estimated from the Allen-Dynes modified McMillan equation [14]: where µ * is the Coulomb repulsion parameter. Our obtained λ = 1.27 for bulk Pb is slightly smaller than that in the previous calculation (1.41) [9]. When taking µ * = 0.1, the calculated superconducting transition temperature 6.04 K for bulk Pb is very close to the measured values for much thicker films [4,6]. Hence, we take µ * = 0.1 in the following calculations. It is well known that there are three main factors determining T C for Bardeen-Cooper-Schrieffer (BCS)type superconductors: the electronic DOS at the Fermi energy N (E F ), the logarithmic average frequency ω ln and the strength of the EP coupling λ. The calculated N (E F ), λ, ω ln and the estimated T C are summarized in table 1 for ultrathin Pb films. The calculated data for bulk face-centered cubic (fcc) Pb are also listed in table 1 for comparison. For Pb thin films, the calculated N (E F ) is the largest for N = 3 ML, which agrees well with [9]. However, we find that the variation of N (E F ) with thickness does not correlate with the trend in T C for Pb films.
Our obtained T C is the largest for N = 4 ML. This may be related to the behavior of the lattice dynamics. From figure 1, it is evident that the value of the maximum frequency for N = 4 ML is smaller than that for other films. The phonon modes are softer due to the smaller contraction of the topmost interlayer distance for N = 4 ML, as discussed above. Yndurain and Jigato [8] also found that the variation of surface phonon states with slab thickness correlated well with the interlayer relaxation in Pb(111) thin films. According to equation (1), the softening of phonons results in an increase of the EP coupling constant λ and is favorable to an enhancement of T C for N = 4 ML.
It is noted that the values of λ and T C obtained for free-standing ultrathin Pb films are smaller than their bulk values. However, in an early investigation [15], it was estimated that the regions of the film that exhibit the greatest increase in T C were within a few monolayers in both Al and Sn. Our recent first-principles calculations for free-standing Al films also confirmed this [16]. We think that this difference between Al and Pb ultrathin films is ascribed to their different behavior of lattice relaxation. For Al(100) and (111) films, the topmost interlayer distance expands [16,17]. Both this expanding and the low-coordinated atoms at the surface result in the softening of surface-layer phonons and an enhancement of superconductivity for Al films at small thickness [16]. However, for Pb(111) films, the topmost interlayer distance contracts. The above results clearly show that softening and stiffening of surface phonons coexist in Pb ultrathin films. If the effect of the stiffening from contraction of the interlayer distance dominates over the effect of the softening from low-coordinated atoms, EP interaction will be suppressed.

(1×1)-Pb/Si(111) films
Next, we go beyond the hypothetical free-standing Pb approximation, and perform calculations for (1 × 1)-Pb/Si(111) thin films. The experimental lattice parameter of bulk Pb is also used  figure 2(a). The results of the total energy after structural relaxation show that the T1 site is the most stable structure among the other configurations. This result is in agreement with experiment [18] and other calculations [19].
Hence, we will focus on this phase alone in the following calculations. We first investigate the fcc ABC ABC . . . stacking sequence for the Pb films. There is compressive lateral strain in the Si substrate because of using the experimental lattice parameter of bulk Pb. In order to release this strain, the results of geometrical relaxation show increased interlayer spacings for the Si substrate with respect to its bulk value. However, the interlayer distance between the first deposited Pb layer and the Si substrate decreases substantially, showing the existence of a strong covalent bond at the interface. The experimental range of the covalent bond length between Si and Pb d Pb−Si is 2.65-2.93 Å [20]. Our calculated covalent bond length d Pb-Si is 2.69 Å for N = 1 ML Pb and about 2.74 Å for N 2 ML Pb. The strong covalent bond between Si and Pb pulls the first deposited Pb layer from the second Pb layer, leading to a substantially expanded interlayer spacing. On the other hand, the interlayer distance between the first two Pb layers on the side of vacuum contracts, as does that for the free-standing films. The calculated interlayer relaxations d 1,2 and d 1,2 are listed in table 2. Here d 1,2 and d 1,2 correspond to the relaxed results near the vacuum and substrate, respectively. The calculated phonon DOS F(ω) and the EP spectral functions α 2 F(ω) for the (1 × 1)-Pb/Si(111) with N = 4 ML Pb are shown in figure 3, as an example. The vibrational modes of Pb are below 90 cm −1 , while the vibrations of substrate Si atoms occupy the high-frequency region extending up to 600 cm −1 . From the curve of the EP spectral functions α 2 F(ω), we can see that the low-frequency modes from Pb make the main contribution to EP coupling, whereas the high-frequency modes from the interface Si make a small contribution to EP coupling. When the different numbers of layers of Pb are deposited on Si substrate, the overall structures of the phonon DOS are similar. But the subtle features of vibrational structures from Pb atoms are different due to the different extent of interlayer relaxation.
For the (1 × 1)-Pb/Si(111) films with N = 1-5 MLs, the calculated N (E F ), λ, ω ln and the estimated T C are summarized in table 2. We can see that these properties oscillate in a perfect even-odd fashion from 1 to 5 ML. However, only the calculated N (E F ) shows the bilayer oscillating pattern in the case of free-standing Pb films. The calculated N (E F ), λ and T C show a local maximum for even monolayers, while the calculated ω ln shows the inverse trend. When only one ML is deposited on Si, the strong covalent bond between Pb and Si at the interface leads to a strong stiffening of Pb phonon modes and the smallest EP coupling constant λ. As in the case of free-standing Pb films, the obtained λ and T C for N = 4 ML are larger than those for other values of the number of layers, but are still smaller than their bulk values.
When the number of Pb layers deposited exceeds two, we also investigate the hexagonally close-packed (hcp) stacking for Pb layers on the Si substrate. Some calculated results are listed in table 2. The calculated ω ln for each layer in hcp stacking is smaller than that in fcc stacking, showing that the phonons probably soften due to the different interlayer relaxations in the hcp case. The softening modes lead to a smaller ω ln , but make a large contribution to λ. So the obtained T C in hcp stacking is larger than that in fcc stacking when N 3 ML. For N = 4 ML, the softening phonon modes and the high N (E F ) are both favorable to superconductivity, resulting in the maximum T C , which is larger than the calculated bulk value. This result is very much consistent with the experimental study by Qin et al [6]. The transition temperature T C as a function of film thickness showed a maximum at N = 4 (see figure 1 of [6]). So our calculations imply that the epitaxial ultrathin Pb films on Si(111) surface may possess an hcp rather than an fcc structure.  Finally, we study the one-atomic-layer film of Pb on Si(111) called the striped incommensurate (SIC) phase, which has a Pb coverage of 4/3 ML and is modeled by a √ 3× √ 3 unit cell, as sketched in figure 2(b). There are four Pb atoms per three surface Si atoms. The experimental lattice parameter of Si is used, corresponding to the compression of Pb film by 5% in the SIC phase. The initial positions of four Pb atoms we put are (0, 0), (0.5, 0), (0, 0.5) and (0.5, 0.5) in relative coordinates of the primitive lattice vectors. After structural relaxation via the total energy minimization, we find that one Pb atom is still located at the H3 (0, 0) site and the other three are located at the bridge sites close to the T1 site (i.e. the position of a neighbor surface Si atom). The Pb overlayer forms a slightly corrugated surface. The H3 Pb adatom is about 0.21 Å lower than the other three atoms. The distance between the H3 Pb adatom and the Si substrate atom is 3.31 Å, whereas it is 2.81 Å for the distance between the other three Pb atoms and Si substrate atoms. According to the experimental range of the covalent bond length d Pb−Si (2.65-2.93 Å) [20], our structure data show that three of the four Pb atoms each form a covalent bond with an underlying Si atom, leaving one Pb atom without bonding to the Si substrate.
The calculated phonon DOS F(ω) and the EP spectral functions α 2 F(ω) for the SIC phase are shown in figure 4. The results for (1 × 1)-Pb/Si(111) with N = 1 ML Pb are also shown in figure 4 for comparison. Obviously, for the SIC phase, there are enhancements in the low-frequency range for both the F(ω) and α 2 F(ω) curves. Furthermore, the extending range of frequency is smaller in the SIC phase. Although the in-plane lattice constant is compressed by 5% in the SIC phase, the vibrations are softer than those in the (1 × 1) phase resulting from larger distances between Pb and Si atoms. The softening modes make a large contribution to the EP coupling. The calculated λ (0.737) for the SIC phase is much larger than that for the (1 × 1) phase with N = 1 ML (0.431), whereas the calculated ω ln (65.12 K) for the former is much smaller than that for the latter (129.10 K). The estimated T C is 2.56 K for this SIC phase, which is close to the experimental value (1.83 K) [7]. Our comparative study clearly shows that the dynamical properties and EP coupling at the interface depend strongly on the bonding environment between the overlayer and substrate atoms.

Summary
In summary, we have studied with first-principles methods the lattice dynamical properties and EP coupling of three types of Pb(111) films with a thickness of several MLs. For free-standing ultrathin Pb(111) films, our calculations have shown that softening and stiffening of the surface phonons coexist. But the effect of the stiffening from contraction of the topmost interlayer distance dominates over the effect of the softening from low-coordinated atoms, resulting in the decrease of λ and T C , which are smaller than their bulk values. For the (1 × 1)-Pb/Si(111) films with N = 1-5 MLs, the calculated N (E F ), λ, ω ln and T C all oscillate in a perfect even-odd fashion. When N 3 ML, the obtained λ and T C in the hcp stacking are larger than the values in the fcc stacking. Furthermore, the estimated T C is the largest and larger than the bulk value for N = 4 ML in hcp stacking, which is very much consistent with the experimental study. Our calculations therefore imply that the epitaxial ultrathin Pb films on Si(111) surface may possess an hcp rather than an fcc structure. For the SIC phase, our structure data show that three of the four Pb atoms form a covalent bond with an underlying Si atom, leaving one Pb atom without bonding to the Si substrate. The vibrations in the SIC phase are softer than those in the (1 × 1) structure of 1 ML Pb on Si, resulting from larger distances between Pb and Si atoms in the former. The estimated T C in the SIC phase is close to its experimental value. Our calculations clearly show that the dynamical properties and EP coupling of ultrathin films depend strongly on the interlayer relaxation and the bonding environment between the overlayer and substrate atoms at the interface.