Strained epitaxial interfaces of metal (Pd, Pt, Au) overlayers on nonpolar CdS ( 1 0 1 0 ) surfaces from first-principles

The depositions of (1 1 1) and (1 0 0) overlayers of Pd, Pt and Au on the CdS (1 0 1 0) surface are studied within epitaxial mismatches of 6% – 7%, using spin-polarized density functional theory. For both compressively strained and tensile-strained interfaces, the (1 0 0) overlayers were found to be thermodynamically more stable owing to better interfacial matching, and higher surface uncoordination resulting in higher reactivity. Pt(1 1 1) exhibits slip dislocations even for five-atomic-layer thick Pt slabs. Along with the leading metal-S interaction, the interfacial charge transfers indicate a weak metal-Cd interaction which decreases in strength in the order Pd > Pt ∼ Au. For the same substrate area, the accumulation of electronic charge for Pt overlayers is ∼ 1.5 – 2 times larger than that of Pd and Au. The n-type Schottky barriers of Au overlayers with the minimum mismatch are within 0.1 eV of the predictions of Schottky – Mott rule, indicating a relatively ideal, scantily reactive interface structure. This is in clear contrast to the Pt epitaxial overlayers which deviate by 0.6 – 0.8 eV. of (Pd, Pt, overlayers on nonpolar CdS ( 1010 ) from first-principles in the 10.1088/1361-648X/ab3919

S Supplementary material for this article is available online (Some figures may appear in colour only in the online journal) Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. theory predictions of the bandgap values are severely underestimated; the leading cause being the self-interaction error or the delocalization error of the local and semi-local exchange-correlational functionals [6,7]. To overcome these deficiencies, self-interaction correction schemes [8], manybody perturbation theory [9] or orbital-dependent Hubbardcorrection schemes [10] are employed. Each of these alternatives either throw challenges in computational expense or in slow conv ergence [11], or are limited to relying on an empirical or fixed U-value over the chemical process [12]. Typically, for bulk semiconductors, the prediction of valence and conduction band edges from local density approximation (LDA), semi-local Perdew-Burke-Ernzerhof (PBE) or hybrid (HSE) exchange-correlational functionals are benchmarked against many-body perturbation theory calcul ations in the GW approximation (G for Greens function and W for screened Coulomb interaction) or higher approximations. For the case of bulk wurtzite CdS, the band gap predictions by semi-local functionals like PBE underestimate the bandgap by 40%-50% depending on the chosen lattice parameters, while HSE and GW 0 predictions underestimate by ∼15% [13,14]. On the other hand, in attempting to obtain the experimental band gaps with an orbital-dependent DFT+U scheme, our tests show that one needs to go to unreasonably high U-values. In this work, we have relied on standard DFT to obtain accurate conduction band minima for the CdS, which is further elaborated on in section 3. This enables to make a reliable prediction of n-type Schottky barrier heights at the metal-CdS interfaces.

Computational settings
All calculations were performed with the plane-wave implementation of density functional theory, within the VASP package [15,16]. The projector augmented-wave potential sets [17,18] were used to model the core electrons. A planewave cutoff of 400 eV was used for the one-electron wavefunctions, and 560 eV for the augmentation charges. We have used the PBEsol exchange-correlational developed by Perdew and co-workers [19], for predicting properties of solids and the slowly varying electron density near the surfaces. Dipole corrections [20] were added to cancel the spurious electrostatic interactions between neighboring periodic cells. The electron density was self-consistently converged within 10 −5 eV, while the ions were relaxed by the Hellman-Feynman forces to less than 0.02 eV Å −1 . The Brillouin zone sampling was done using the k-mesh of 6 × 4 × 1, 6 × 2 × 1 and 2 × 2 × 1 in the increasing order of three supercell sizes shown in figure S1 of the supporting information (available online at stacks.iop.org/JPhysCM/31/505001/mmedia) (SI).
In modeling epitaxial interfaces it is essential to have lattice parameters, and their resulting mismatches, close to their experimental values. In our work, this is naturally met with the use of PBEsol exchange-correlational functional [19]. It is the most suited GGA-level functional to predict the equation of states, workfunctions and surface energies for 4d and 5d metals [21][22][23]. One downside which comes with this functional is its poor performance in describing the localized 3d metallic systems [24]; which also limits this work to examine interfaces of only 4d (Pd) and 5d (Pt and Au) metals.
The CdS (1 0 1 0) surface unit slab is six bilayers thick with at least 15 Å of vacuum, and is passivated by pseudohydrogens [25,26] for bulk-like termination. The (1 1 1) and (1 0 0) metal-overlayers, which are deposited over the CdS substrate, are five atomic-layer thick. All the atoms in the simulation cell are relaxed.
Among the metals considered here, it is theoretically [27][28][29][30][31][32] and experimentally [33,34] well known that strained Pd crystals and nanostructures exhibit ferromagnetism, as being isoelectronic to Ni, it has a high density of states at the Fermi level. Owing to this, our geometry relaxation calculations of Pd-overlayers and their interfaces have been carried out within spin-polarized theory. It should be cautioned that the spin-polarization in Pd slabs (both, isolated and deposited) make the potential energy surfaces very complex, where several isomer states exist with varied magnetic moments, within few tens of meV. The general trend is that the tensile strain (or expansion) of Pd slabs and crystals from their equilibrium geometry marks the onset of ferromagnetism, while in the compressive regime, the magnetization is severely quenched. The Au and Pt slabs are non-magnetic and are treated within standard DFT. The CdS slabs in themselves are non-magnetic, both before and after forming the interfaces with metal overlayers.

Results and discussion
The lattice constants for both components of the interface are reported in table 1. The agreement of calculated values (a DFT ) is within 1% difference of the experimental data a exp , where values for the metals are corrected for zero-point vibration effects [21]. Our lattice constants of Pd, Pt and Au are also in close agreement with previous theoretical reports using the PBEsol exchange-correlational functional [22,23]. From the bulk unit cells of these face-centered cubic (FCC) crystals, orthogonal surface unit cells of (1 1 1) and (1 0 0) orientations were constructed using the ASE package [35]. Further below in table 1, the calculated workfunctions (φ DFT ), E vac − E F , are shown for clean metal slabs. Here, E vac is derived using the level at which electrostatic potential energy is constant (just outside the slab) in the direction perpendicular to the surface, and E F is the Fermi energy of the composite system. The surface unit slabs were calculated with a k-point sampling of 30 × 16 × 1 and 30 × 30 × 1 for orthogonal (1 1 1) and (1 0 0) unit cells. Our tests for an eight-layer unit slab resulted in workfunction values within 0.05 eV of the five-layer slabs used in this study. Comparing these with the workfunctions reported by Perdew et al [23], calculated with the PBEsol functional and eight-layer slab thickness, the differences are within 0.18 eV. In comparing the theoretically calculated workfunctions with the experiment, one should be careful that the error bars in the experimental measurements can be as large as 0.3 eV [36,37], which can be caused by the experimental technique, sample anisotropy, intrinsic surface defects or adsorbed species. The experimental values (φ exp ) are reported from the theoretical literature [36,37], wherein the original experimental works or relevant reviews are discussed in detail. It is generally known that the workfunctions of densely packed (1 1 1) surfaces are higher than those of the (1 0 0) surfaces [31,37,38]. This is observed in the present study for five-layer Pd and Au, while for the Pt(1 1 1) and Pt(1 0 0) the difference in φ DFT (5 meV) is within the convergence error bars caused due to limited slab thickness.

Clean CdS(1 0 1 0) and strained metal surfaces
A schematic of the surface unit cell of a six-bilayered thick CdS (1 0 1 0) substrate, which has [1 2 1 0] and [0 0 0 1] as its lateral directions, is shown in figure 1(a). The suitability of our modeled substrate can be seen from figure 1(b), where the band-edges of CdS (1 0 1 0) are shown as a function of substrate model thickness from 4-20 bilayered slabs. The bandedges are referenced from the respective E vac for each slab, where the vacuum thickness is at least 30 Å . The increase in the substrate thickness is marked by the increase in bulklike states at the conduction band edge, which leads to the lowering of the conduction band minima (CBM), while the valence band maxima are almost constant. The dotted line marks the experimentally reported CBM for the (1120) surface at E CBM = 4.79 eV [40], while the CBM of a six bilayered thick CdS slab of the same orientation has E CBM = 4.71 eV , in excellent agreement with the experiments. In our study, we have used the (1 0 1 0) slab of the same thickness, passivated with pseudohydrogens, which do not affect the band-edge positions. At this thickness, the (1 0 1 0) surface bandgap is 1.25 eV and the (1 1 2 0) bandgap is 1.35 eV, which is a ∼45% underestimation in comparison with experimental reports [39,40]. Figure 1 shows that much of the error in the CdS bandgap prediction occurs mainly at the valence band edge, in agreement with recent theoretical reports [14,41]. This finding also facilitates estimating the n-type Schottky barrier height from the CBM edge, within the standard DFT treatment, rather than using experimental bandgaps, as in the early ab initio literature [42].
the substrate supercells were limited to ∼2.5 nm, and chosen such that the metal-overlayers make minimal misfits with one of the supercell directions and less than ∼7% in the other direction. The interfacial mismatches with CdS (1 0 1 0) are listed in the second column of table 2. The lattice mismatches in one of the directions are < 1%, making the overall strain almost uniaxial. In the following discussion, we will focus on the larger of the two lattice mismatches of the interface, with |m i | > 1, as the one characterizing the interface. For the cases studied here, the variation in φ is about ∼0.2 eV, which is known to increase if multiple strain modes are existing at the interface [45]. The experimental studies on strained metals are difficult to correspond exactly with the theoretical studies. This is because the length scales of experimental samples are usually in the millimeter range, and witness the local variation of stresses, which are difficult to include in simulations of nanosized structural models. This poses difficulties in correlating the mechanical and electronic responses of the studied system. Second, in the plastic regime dislocations and surface roughness can increase, and cause significant changes in the surface dipoles. However, these defects cannot be compared with theoretical models with certainty, and thus, impede clear correlations with the predictions made by theory. In the elastically strained surfaces, the theory and the experiments still show qualitatively similar trends [46,47]. Figure 1(c) shows that within the strained regimes examined here, φ (1 1 1) > φ (1 0 0) holds for Pd and Au surfaces, and is very similar for Pt surfaces. The workfunction values of unstrained five-layer metal slabs which were shown in table 1 are at the zero of the abscissa. The variation in φ with respect to lattice mismatch is almost linear for both Pt, Pd surfaces and Au(1 0 0), with a slight deviation for the Au (1 1 1) surface. For the strained Pt(1 1 1) surfaces, the two data points not connected by the line, represent the cases where dislocations were formed upon relaxation. These two cases, which are energetically more stable than their elastically deformed counterparts, will be later discussed when discussing the interfaces. Here, the point to note is that the dislocated surfaces have lower φ in comparison to smooth surfaces (plotted along the line). This is in agreement with experiments in the plastic regime [48], and a simple electrostatic model [49]. The workfunctions for the dislocated surfaces are lower owing to the emerging surface dipoles, which lead to a shallower surface potential, and a lower φ.

Metal-CdS(1 0 1 0) interfaces
For the interface calculations, the metal overlayers were directly introduced to the substrate supercell, without relaxing the metal surface separately. In this way, from the onset of relaxation we expose the metal-overlayer to the asymmetric stress between the bottom and top metal layer, as well as to the overall stress of supercell misfits. As an initialized condition, the bottom layer of the metal chemically interacts with the substrate, while the top layers witnesses the shear stress due to the lattice misfit and also due to the chemical anchoring of the interacting bottom layer. This is the natural way to relax an epitaxial interface since all the asymmetries of the growth process are initialized in the model calculation. Another way of relaxing an epitaxial interface is by first relaxing the overlayer in the cell (having the substrate dimensions), and subsequently introducing the substrate. In this alternative approach, the shear stresses applied on the metal are symmetric, in the sense, that both sides of the metal overlayers witness no effect of the substrate except for the lattice misfit of the supercell. The geometries from this symmetric-stress calculation result in elastically sheared metal, chemically interacting with the substrate below it, without any deformation. Indeed, in both relaxation approaches the metal lattice constants in the directions perpendicular to misfit also change due to the Poisson effect. We have also used this symmetricstrain method to examine the energy differences of dislocated and elastically deformed surfaces. Since (1 1 1) planes in FCC systems usually undergo slipping, we have tested this alternative symmetric-stress calcul ation for (1 1 1) overlayers. For the thickness of five-layer slabs, the geometries derived for the (1 1 1) metal-overlayers using the two methods give  (1 1 1). These interfaces with dislocations are shown in figure 2, where the surface dislocations are seen to emerge owing to the initialized asymmetric stresses and periodic boundary conditions imposed on the metal-overlayer. The slip planes of the dislocations are the [1 1 0] planes. The energetics of these cases will be discussed below. The interfaces are listed in table 2 with their adhesion energies ( E adh ) in units of energy per interfacial area (J m −2 ) and energy per metal atom in the interfacial layer. Here, the E adh is defined as: where E int is the total energy of the interface, and E M and E CdS are the energies of the individual comp onents. The interfacial energy can be defined either per atom at the top surface layer or per unit interfacial area. As the area is determined by the substrate cell dimensions, it is the same for both types of overlayers in each case. However, differences in the planar density of the metal overlayers warrant a thermodynamic assessment also in terms of eV/atom. Comparing the cases with the same substrate area (same sign of m i ), the (1 0 0) metal overlayers have a lower lattice mismatch than those of the (1 1 1) overlayers for Pd and Pt. The lower interfacial mismatches result in better exposure to the substrate, and consequently, enhances chemical bonding. This fact, along with the undercoordination of the (1 0 0) metal atoms, results in higher interface adhesion energies (E adh , both, in terms of J m −2 and eV/atom) for all the metals. In general, table 2 shows that, for epitaxies with periodicity less than the length regimes considered here, the (1 0 0) surfaces are expected to form thermodynamically more stable interfaces with CdS(1 0 1 0). Despite the lower planar density of atoms in (1 0 0) overlayers, and having m (1 0 0) ∼ m (1 1 1) for the Au interfaces with the same misfit sign, the adhesion energies (J m −2 ) of (1 0 0) overlayers are significantly higher than their (1 1 1) counterparts. For Pt(1 1 1) the values in the parentheses are the adhesion energies for the elastically sheardeformed interfaces, while the ones not in the parentheses are for the dislocated Pt surfaces (shown in figure 2). Figure 2(a) shows the Pt(1 1 1) overlayer in compressive strain and the panel (b) in tensile strain. The total interface adhesion energies of the two dislocated surfaces are higher than or equal to the elastically deformed ones. Intermetallic comparison of the interfaces indicates that Pd and Pt have comparable adhesion energies, while Au overlayers have reduced adhesion due to deeper and less reactive d − bands, resulting in weaker chemical interaction at the Au-CdS interfaces.
For assessing the strength of chemical interaction and assigning charge transfers between the metal and substrate we have used the Bader charge partitioning and integration scheme [50]. In figure 3 we show the Bader charge differences for the (1 0 0) interfaces of three metals as a function of interfacial area. First, all the metals show electronic charge accumulation and it increases as expected with the interfacial area. Further, by the basic sum of electronic charges, for every interface, we calculated the Bader charge differences of three interacting species (metal, Cd and S) adding up to zero. The values for each of the species represent the Bader charge differences between the relaxed interfacial geometries shown in figure S1 and their single-point geometry without the overlayer/substrate. For all the interfaces of Pd, Pt and Au, as |Q S | > |Q Cd |, it is clear that metal-S interaction is the  leading one; whereas, the metal-Cd interaction is a weak one. In our previous work, we had demonstrated the trend of the Ni-group metals to form a complex with Cd, where the interaction decreases with the decrease in the localization of the d-bands (i.e. down from 3d to 4d to 5d) [3]. In figure 3 the data-points of the three metals with the same cross-sectional area of 167 Å 2 , show slightly larger charge accumulation at Cd (i.e. stronger interaction) for the Pd case, in comparison to Pt and Au. This result is also in agreement with the aforementioned trend observed for the Ni-group metals in our previous work [3].
Further, in figure 4, we show an intermetallic comparison of sum of Bader charge differences (Q M ) between the deposited and isolated metal-overlayers, for the case of the same interfacial area of 8.27 × 20.19 Å 2 . Pt overlayers show the highest charge accumulation, while Pd and Au show similar charge accumulation at almost half of the Pt value. Although the (1 0 0) and (1 1 1) overlayers have 21 and 24 metal atoms per layer for this substrate dimensions, it is the (1 0 0) overlayer which results in higher charge accumulation for Pt and Pd. On average, this suggests a stronger interaction per atom of the (1 0 0) overlayers owing to higher undercoordination than in the case of (1 1 1) overlayers, resulting in higher reactivity. The inset of figure 4 shows Q M as a function of the interfacial area of all the interfaces considered on nonpolar CdS substrate. The dashed line indicates the interface which was discussed in the bar plots of the main panel. Other cases also indicate that charge accumulation of the (1 0 0) overlayers is, in general, slightly greater than or equal to that of deposited (1 1 1) metal overlayers. Increasing the interfacial area for an epitaxial interface of given stress would result in a linear increase in charge transfer. Since there are only two data-points for (1 0 0) and (1 1 1) overlayers of each metal, in the inset we have plotted a linear fit (orange colored line) using the data-points for Pd and Au interfaces. The sum of squared residuals is only 0.05e, indicating that the Au and Pd interfaces show very similar charge transfers. It is important to note that the uncertainty in determining the accurate ground state charge density of Pd slabs, due to several magnetic isomers, naturally bears errors in estimating the charge transfers at its interface. Our tests show that the variation in the Q M is ∼0.05e.
Finally, owing to the excellent prediction of CBM of the CdS substrate within the level of theory used here, we calculated the n-type Schottky barrier heights (SBH) of these interfaces. The n-type SBH is theoretically derived using the difference in CBM and the Fermi level of the relaxed interface. These energy levels were obtained for each of the interfaces with a single-point calculation with a dense mesh of k-points: 8 × 5 × 1, 8 × 4 × 1 and 4 × 4 × 1 in the order of increasing supercell area. The CBM of the substrate is not trivial to estimate as the metal-induced gap states populate the bandgap and have a contribution from the metal overlayer and a few top substrate layers. Here, we obtain the substrate CBM by using the electrostatic bulk average potential as a reference, which is derived from the deeper layers of the substrate slab of the relaxed interface geometry. Figure 5 plots the SBH of the interfaces as a function of the workfunction of the strained and isolated metal overlayers (shown in figure 2). This allows a comparison with the simplest, chemically unreactive interface prediction of the Schottky-Mott rule, which equates the SBH to the difference of the two energy levels (with respect to vacuum) in the metal and semiconductor's isolated forms. This is indicated by the dotted black line, y = x -EA, where EA = 4.71 eV. The only interface with SBH within 0.1 eV of the Schottky-Mott values are the well-matched interfaces of Au, with m ∼ 1.9%. The calculated SBH values are severely underestimated in comparison to the experimental reports on bulk-like clustered interfaces [51]. The Pd and Pt interfaces, as well as the strained Au interfaces, indicate a relatively reactive interfacial structure, as the differences with the ideal Schottky-Mott rule increase substantially. The lowest SBH is for the Pd(1 0 0) interface with 0.1 eV. The variation in the SBH of (1 1 1) interfaces due to the change in the strain from compressive to tensile, is within 85 meV for all the metals. The Au(1 0 0) interfaces also show a variation of only 34 meV, while the Pd(1 0 0) and Pt(1 0 0) interfaces show a 0.2 eV difference.

Conclusions
We show by means of first-principles calculations that (1 0 0) overlayers of Pt, Pd and Au are preferred over (1 1 1) overlayers on CdS(1 0 1 0) surfaces, at least within periodicity dimensions of 2 nm. Both a better match of surface lattice constants and the undercoordination of surface atoms at (1 0 0) overlayers contribute to the interface stability. We also report that even for five-atomic-layer thick Pt(1 1 1) slabs, dislocations are formed in the overlayers both for compressive and tensile strains. The metal-surface interaction has a leading contribution from the S-bonds and a weaker interaction with Cd, where the latter's strength decreases from 4d to 5d-metals. The minimally strained Au interface with misfits of only 1.9%, shows agreement with the simple Schottky-Mott rule. Also, unlike the (1 0 0) interfaces which are more reactive, the Schottky barrier heights for (1 1 1) interfaces does not substantially change with different epitaxial strains.

Supplementary material
The supplementary material file shows all the studied metal-CdS interfaces for different mismatches.