- Email: [email protected]
A combined ab initio and photoelectron study of galena (PbS)
Surface Science 491 (2001) 226±238
www.elsevier.com/locate/susc
A combined ab initio and photoelectron study of galena (PbS) Joseph Muscat a,*, Craig Klauber b a
CSIRO Minerals, Box 312, Clayton South, Vic., 3169, Australia CSIRO Minerals, P.O. Box 90, Bentley, WA, 6982, Australia
b
Received 21 February 2001; accepted for publication 29 June 2001
Abstract The results of a study of the physical, elastic and electronic properties of PbS using both ®rst principles theoretical and experimental X-ray and ultraviolet photoelectron techniques are presented. We have investigated the in¯uence of the computational approximations such as the basis set, pseudopotential and the treatment of exchange and correlation on the accuracy of the theoretical results. The current study work con®rms that the binding energy of the S 3s state lies at 13 eV, in agreement with previous experimental work and contrary to earlier ab initio Hartree±Fock (HF) calculations. We establish that the reason for the discrepancy between previous experimental and theoretical determinations of the electronic structure of PbS is due to the treatment of the exchange interaction. We demonstrate that calculations using the HF approximation tend to overestimate the binding energies of the valence and core states yielding densities of states in poor agreement with experiment whereas calculations performed with density functional theory within either the local density or the generalised gradient approximation give a signi®cantly better description of the electronic structure. Finally, we con®rm that hybrid functional techniques such as B3LYP provide an excellent description of the physical and electronic structure of PbS. Ó 2001 Elsevier Science B.V. All rights reserved. Keywords: Ab initio quantum mechanical methods and calculations; Low index single crystal surfaces; Density functional calculations; Surface electronic phenomena (work function, surface potential, surface states, etc.); Sulphides; X-ray photoelectron spectroscopy; Visible and ultraviolet photoelectron spectroscopy
1. Introduction Lead remains an important industrial commodity, sourced primarily from galena (PbS) with minor contributions from cerussite (PbCO3 ) and anglesite (PbSO4 ). Prior to smelting, the galena component in natural ores is generally concentrated or bene®ciated via froth ¯otation. Galena,
* Corresponding author. Tel.: +61-39545-8500; fax: +6139562-8919. E-mail address: [email protected] (J. Muscat).
when freshly crushed, exposes naturally hydrophobic surfaces that readily attach to ¯otation bubbles, hence galena represents the classic example of bene®ciation via ¯otation. Weathering and other process variables can reduce this hydrophobicity and so the surface chemistry of galena, particularly its oxidation and the role played by collectors such as xanthates (to increase hydrophobicity) are key factors in a fundamental understanding of the ¯otation process. Considerable empirical eort over several decades has been directed toward an understanding of the surface chemistry of galena, in particular the predominant
0039-6028/01/$ - see front matter Ó 2001 Elsevier Science B.V. All rights reserved. PII: S 0 0 3 9 - 6 0 2 8 ( 0 1 ) 0 1 4 0 8 - X
J. Muscat, C. Klauber / Surface Science 491 (2001) 226±238
(1 0 0) cleavage plane. However, to eciently progress the understanding of such a mineral system by means of molecular modelling the ®rst essential step is to better understand the electronic structure. Once this baseline understanding has been achieved it should be possible to progress to modelling aspects such as the galena±xanthate adsorption or the oxidative reaction of the galena surface to a range of dierent valence states such as elemental sulphur, lead sulphate and other surface speci®c species. Galena is one of the most studied sulphide minerals due in part to the abundant availability of good quality single mineral samples and its simple NaCl-halite structure. There have been a number of experimental studies of the structural, elastic and electronic properties of PbS but many of these basic properties are not known accurately. For example, considerable uncertainty remains in the bulk modulus of PbS despite a number of measurements, with two of the most recent studies giving values of 48 GPa [1] and 73 GPa [2] respectively. Doubt has also been expressed as to the electronic structure PbS. This has recently been probed by photoemission experiments [3] and also calculated theoretically in subsequent ab initio simulations using Hartree±Fock (HF) theory [4]. The two studies yielded results that were generally in broad agreement with the exception of the determination of the binding energy (BE) of the S 3s band which was measured to lie at around 12 eV but predicted by the simulations to be about 17 eV. Mian et al. [4] note that the experimental data from Santoni et al. [3] did not extend to include the 17 eV BE region and suggest that the peak at about 12 eV could be due to HS-type species formed due to contamination of the surface by hydrogen. The photoelectron spectroscopies of X-ray photoelectron (XPS) and ultraviolet photoelectron (UPS) are excellent methods for obtaining core and valence level BE information for occupied levels. Both are ®nal state spectroscopies, hence the actual observed BEs are altered by ®nal state relaxation eects shielding the core hole. Such extra-atomic relaxation is almost independent of the core hole (e.g. spectra from noble gas atoms implanted in noble metals [5]) so that the relative
227
BE values are useful parameters to check the validity of electronic calculations. Earlier angle resolved HeI UPS work on PbS shows sharp structure in the 11±13 eV range which becomes less apparent at the higher photon energies probed in Santoni et al. [3]. In the later work the strongly dispersing feature at 12 eV is assigned by the authors to S 3s, at considerable variance to the HF computed value of Mian et al. [4] of 17 eV. Unfortunately, in the region 12±17 eV and below, artefacts from the steeply rising secondary electron emission (SEE) background, close to the vacuum zero, restricts the usefulness of the low UV energies in determining a reliable S 3s BE. This is well recognised in the earlier HeI work on galena [6,7] in which discussion of the electronic structure is limited to above 6 eV. In these studies, simple angle integrated HeII and MgKa photoemission spectra at normal emission are used to reassess the structure from the Pb 5d5=2 line to the Fermi level, in particular the S 3s. Since the work of Mian et al. [4], there have been signi®cant improvements in ab initio techniques, in particular the development and implementation of hybrid functional techniques [8,9] that combine elements of HF and density functional theory (DFT). These hybrid techniques have been shown to generally predict the electronic structure [10] and the structural and elastic properties [11] of a range of materials signi®cantly more reliably than the previous HF or DFT methods in common use. In the current study, we have performed ab initio calculations together with photoemission experiments of the valence band in order to clarify the situation as regards to the electronic structure of PbS. The calculations of electronic structure and the physical and elastic properties of PbS have been performed using the new three parameter B3LYP hybrid functional method [8,9] and the more traditional HF and DFT approximations in order to assess the performance of the dierent theoretical techniques and to enable direct comparisons with previous theoretical calculations. The in¯uence of other approximations inherent in ab initio calculations such as the basis set and the pseudopotential on the computed properties have also been examined as a prerequisite to future
228
J. Muscat, C. Klauber / Surface Science 491 (2001) 226±238
studies of the surface structure and chemistry of PbS surfaces.
2. Method 2.1. Computational details All the calculations have been performed using the C R Y S T A L 9 8 ab initio software package [12] which is based on the periodic, linear combination of atomic orbitals (LCAO) formalism where the Bloch orbitals of the crystal are expanded using atom centred Gaussian orbitals with s, p or d symmetry. A variety of self-consistent treatments of exchange and correlation have been used including HF, the LDA to DFT, the GGA to DFT and the B3LYP functional based on Becke's three parameter hybrid functional (B3LYP) employing a combination of HF and DFT exchange [8]. The DFT±LDA calculations were performed using LDA exchange [13] and the correlation functional parametrised from Monte-Carlo calculations [14] by Perdew and Zunger [15] whereas the GGA calculations employed the exchange and correlation functionals by Perdew et al. [16]. The B3LYP hybrid functional used in this study is a modi®ed version of Becke's original three-parameter exchange-correlation functional [8] as suggested by Stephens et al. [9] where the exchange and correlation is given by: EXC 1
a0 EXLDA a0 EXHF ax DEXB88
ac ECLYP 1
ac ECVWN
1
where ax DEXB88 is Becke's gradient correction to the exchange functional [17], and the correlation is given by a combination of the Lee, Yang and Parr (LYP) [18] and the Vosko et al. [19] functionals. Becke's original three empirical parameters (a0 0:2, ax 0:72 and ac 0:81) are used which were found to optimise the atomisation energies, ionisation potentials and proton anities of a number of small molecules [8]. Solid state calculations using this functional have been shown to yield band gaps in excellent agreement with experiment [10].
Due to the computational expense required to describe the large number of electrons and the problems associated with treating relativistic effects in core shells for heavy elements such as lead, the core electrons up to and including the 5p shell have been replaced with a pseudopotential. Recent studies have however shown that the pseudopotential approximation can be a source of signi®cant error [20]. We have examined the in¯uence of the pseudopotential on the computed results by performing calculations using the Hay and Wadt [21], the Durand [22] and the recently released Stoll and Preuss type pseudopotentials [23,24]. The sulphur ions have been treated using an all-electron basis set. In the current work, we have used basis sets that have been developed and optimised for use in recent studies of PbS [4]. We have performed tests using HF, DFT and B3LYP to determine if augmentation of the Pb or S basis sets by the addition of d-symmetry functions (denoted by the ``'' superscript) to the Pb or S basis set has a signi®cant in¯uence on the computed cell parameters or total energy. We ®nd that addition of d-functions to the S ions has a small, but signi®cant in¯uence on the predicted cell parameters and total energy but similar improvements to the Pb basis set make little dierence to the computed results (see Fig. 1). Calculations performed using the S basis set yield lattice parameters that are smaller than systematically about 0.03±0.06 A those due to the S basis set whereas structures computed using the Pb basis set dier from those . Furdue to the Pb basis set by less than 0.02 A thermore, the total energy computed using the S basis set is signi®cantly lower than that due to the S basis set (by 0.8±1.9 eV) whereas the dierence between the Pb and Pb basis sets is much smaller (around 0.03 eV). Since the use of the Pb basis set causes an increase in the computational cost of the calculations but only very small dierences in the computed results compared to the Pb basis set, the discussion in the current article will focus on the results calculated using the PbS basis set. Sampling of K-space has been performed using Pack±Monkhorst grids [12,25] of shrinking parameter eight, yielding 29 symmetry inequivalent K-points for bulk PbS. Tests reveal that increasing the number of K-points produces no signi®cant
J. Muscat, C. Klauber / Surface Science 491 (2001) 226±238
229
Fig. 1. The eect of the basis set and pseudopotential on the computed cell parameters.
dierence in the computed structure or energy of the crystals (the total energy of the bulk crystals is converged to within about 1 10 5 eV). C R Y S T A L 9 8 computes the matrix elements of the Coulomb and exchange terms by direct summation of the in®nite periodic lattice. The truncation of these summations is controlled by ®ve Gaussian overlap criteria, details of the control of these parameters is available elsewhere [12,26]. The values of the overlap criteria chosen in the current study were high (ITOLS parameters set to 10 8 , 10 8 , 10 8 , 10 8 and 10 16 ) in order to converge numerical errors to about 3 meV in the total energy [27]. The structural optimisations were converged to in cell parameters and 10 4 a tolerance of 0.01 A eV in the total energy using a modi®ed Broyden± Fletcher±Goldfarb±Shanno minimisation algorithm [28]. We have computed the bulk modulus of
PbS by ®tting a computed energy volume curve consisting of at least 30 points to the Murnaghan equation of state [29]. 2.2. XPS and UPS methodology Valence band spectra, excited by MgKa X-rays and the HeII UV line were obtained from two galena samples, one from Zacatecas (Mexico) and the other from Broken Hill (Australia). Samples were attached to stainless steel sample stubs by silver-dag and cleaved under a high purity nitrogen atmosphere in a glove-box attached directly to the spectrometer entry lock. The Broken Hill sample exhibited far better cleavage along the (1 0 0) face and only that data is reported here, although the Zacatecas spectra were identical. Both XPS and UPS data were obtained using a VG ESCALAB Mk II (nominal base pressure of about
230
J. Muscat, C. Klauber / Surface Science 491 (2001) 226±238
1 10 10 mbar) retro®tted with multiple channeltron (5) electron detection and an SSI AlKa monochromated source in addition to the original MgKa and AlKa twin anode source. The MgKa (1253.6 eV) spectra were obtained at a source power of 300 W irradiated over the whole sample and the monochromated AlKa (1486.6 eV) spectra (Pb 5d only) at 100 W from a nominal incident spot of 250 1000 lm. In both cases 6 mm analyser slits and a constant analyser pass energy of 20 eV were employed. The UV source used was a VG UVL-HI dierentially pumped discharge lamp which produced a stable HeII output with a main chamber helium pressure of 1 10 9 mbar. UPS spectra were obtained at a constant analyser retard ratio of 10. The spectrometer was calibrated using the Ag 3d5=2 line at 368.22 eV, the Au 4f7=2 line at 83.95 eV and a wide range linearity check using the Cu 2p line [30]. In all cases argon ion etched polycrystalline foils were used for calibration. Final adjustment for the XPS valence region of the galena was against the Au line which yielded a Pb 5d5=2 BE of 18:71 0:05 eV with respect to the Fermi level. In the case of the UPS valence region adjustment was with respect to the Fermi edge of a clean nickel foil which gave a Pb 5d5=2 BE of 18:74 0:05 eV. As the UPS established nickel Fermi edge zero was considered to be more precise the XPS valence region was shifted to align the Pb 5d5=2 peaks. 3. Results 3.1. Structural properties In Table 1, we present the optimal cell parameter for PbS computed with HF, LDA, GGA and B3LYP treatments of exchange and correla-
tion using the PbS basis set. The in¯uence of the pseudopotential on the computed cell parameter is signi®cant. Calculations performed using the Hay±Wadt pseudopotential predict lattice para larger meters that are systematically about 0.1 A than those computed using the Durand pseudopotential, whereas the Stoll±Preuss pseudopotential yields cell parameters which are typically in between those computed by the Hay±Wadt and the Durand pseudopotentials. These dierences are shown in Fig. 1. For comparison, results computed using the PbS, Pb S, PbS and Pb S basis sets are also included which show that the dierences in the computed cell parameter due to the pseudopotential are consistent across a range of basis sets. The treatment of exchange and correlation also has a systematic eect on the computed cell parameter, as can be seen from Table 1. The HF and B3LYP approximations overestimate the experimentally measured cell volume by around 1.9± 3.7% (HF) and 0.4±2.4% (B3LYP) depending on the pseudopotential used. The LDA to DFT underestimates the cell parameter and the GGA corrects for this eect leading to an increase in the cell parameter of about 1% compared to the LDA results. Similar trends arising due to the treatment of exchange and correlation have been reported elsewhere [20]. There is large uncertainty in the literature as to the magnitude of the bulk modulus of PbS, with the most recent measurements giving values ranging from 48 to 73 GPa. The computed bulk moduli, presented in Table 2 are all within this experimental range. The computed bulk modulus is very sensitive to the pseudopotential used; calculations performed using the Hay±Wadt pseudopotential systematically predict the bulk modulus to be approximately 10 GPa lower than
Table 1 The computed cell parameters of PbS. Figures in parentheses indicate the percentage deviation from the experimental cell parameter Durand Hay±Wadt Stoll±Preuss Experiment
HF
LDA
GGA
B3LYP
6.045 (1.85) 6.152 (3.66) 6.076 (2.37)
5.824 ( 1:87) 5.907 ( 0.47) 5.829 ( 1:79)
5.898 ( 0:62) 5.979 (0.74) 5.843 ( 1:55)
5.958 (0.39) 6.080 (2.44) 5.990 (0.93)
5.935
J. Muscat, C. Klauber / Surface Science 491 (2001) 226±238 Table 2 The computed bulk moduli for PbS compared to experimental measurements PS
HF
LDA
GGA
B3LYP
Durand Hay±Wadt Stoll±Preuss Experiment
61 51 55
73 62 72
63 58 60
60 50 56
48±73
that computed using the Durand pseudopotential. The Stoll±Preuss pseudopotential predicts the bulk modulus to be between these two extremes, similar to the trend reported for the computed cell parameter. HF and B3LYP predict very similar values for the bulk modulus of about 50 (Hay±Wadt), 55 (Stoll±Preuss) or 60 GPa (Durand). The GGA predicts bulk moduli that are slighly higher than the HF or B3LYP results whereas the LDA yields bulk moduli that are signi®cantly higher, over 10 GPa larger than the HF or B3LYP values. We have also performed tests to determine if the basis set has an in¯uence on the predicted bulk modulus and have found that the presence of (or lack of) d-functions makes little dierence to the computed bulk modulus.
231
3.2. Electronic structure As noted in Section 1, there is some uncertainty in the literature regarding the BE of the sulphur 3s peak. Recent UPS studies [3] observed a band of states at a BE of around 13 eV which is attributed to the S 3s states in bulk galena. Subsequent HF calculations [4] of galena (PbS) predicted the position of this peak to be 17 eV. The authors of the HF study suggested that the dierence between the experimental and calculated peak positions may be explained by considering the eect that contamination of the PbS surface with hydrogen (which is almost unavoidable in ultra-high vacuum experiments) could have on the photoemission spectrum. Adsorption of H onto the surface of Galena would lead to the formation of HS or H2 S-like species with a lower BE. Mian et al. postulated that the peaks observed by the UPS experiments were actually due to these surface HS-type species and that the states due to the S 3s electrons in galena actually lie closer to the computed BE of 17 eV. As the Santoni experiments did not probe BEs higher than 14 eV, it was not possible to validate this hypothesis. However, as the cross-section for H 1s
Fig. 2. MgKa XPS of the galena valence band region.
232
J. Muscat, C. Klauber / Surface Science 491 (2001) 226±238
Fig. 3. HeII UPS of the galena valence band region. Pb 5d ghost peaks due to the 48.37 eV radiation overlap part of the valence region. Dotted lines show the results from alternative treatments for removal of the Pb ghost peaks.
photoemission is very small, it is unlikely for it to be mistaken for the S 3s state, as supported by calculations of the adsorption of H2 on the PbS (0 0 1) surface [31]. The MgKa XPS spectrum is shown in Fig. 2 and the HeII UPS spectrum in Fig. 3. Both processed curves are compared in Fig. 4. Removal of MgKa X-ray line and satellite structure is straightforward [32], but correction for the He ghost line is not. As noted, a reliable determination in the S 3s region requires a photon energy above that of HeI radiation (21.22 eV). Maximising discharge lamp intensity for the more desirable HeII at 40.81 eV has the disadvantage of strengthening other minor resonance lines. Based on the valence band edges for a clean nickel foil sample the minor lines at 23.09, 23.74 and 48.37 eV [33] were detectable under lamp conditions for optimising HeII intensity. The proximity of the Pb 5d doublet to the valence band means that ghost peaks from the 48.37 eV line will overlap the bottom of the valence band structure. This is evident in the raw data of Fig. 3. Attempted removal of the ghost lines by deconvolution of a source function (as with the MgKa spectra) represented by the 40.81 and 48.37 eV lines proved unsatisfactory due to the high emis-
sion intensity induced at low kinetic energies by the dominant HeI. As the extent of valence band contamination from the 48.37 eV radiation is small (evidenced by the low intensity above the Fermi level in Fig. 3) the best method for removal of the ghost Pb 5d doublet was found to be replication of the HeII doublet and direct subtraction. The HeII Pb 5d5=2 peak was matched to a Voigt pro®le with FWHM of 0.537 eV. The position and intensity of the Pb 5d3=2 ghost peak can be ascertained and so removal of the Pb 5d5=2 ghost peak, assuming the same peak shape and width, simply requires a splitting and intensity ratio. The apparent splitting of the main peaks (D) with HeII is 2.4 eV (background distorted), for the ghost peaks 2.62 eV, for MgKa 2.54 eV and monochromatic AlKa photoemission 2.57 eV. The expected statistical intensity ratio (r) of 1.5 compares to 1.45 for MgKa , 1.38 from the monochromatic AlKa photo emission and 1.49 and 1.42 calculated for MgKa and AlKa [34]. With spectrometer transmission (under CRR conditions) also having an eect the optimum description for the Pb 5d5=2 ghost peak is unclear. In Fig. 3 alternative corrections for the Pb lines are shown (based on D 2:62 eV, r 1:50 and D 2:57 eV, r 1:38) to highlight
J. Muscat, C. Klauber / Surface Science 491 (2001) 226±238
233
Fig. 4. Comparison of the valence band photoemission spectra from MgKa and HeII sources for cleaved galena.
that signi®cant intensity remains at about 11 eV, irrespective of how the removal is performed. The ®rst correction based on observed ghost peak splitting and the statistical ratio is preferred. The processed valence band photoemission spectra from both MgKa and HeII sources for cleaved galena are compared in Fig. 4. The XPS is very similar to earlier AlKa work [6,7] although slightly more structure is displayed. In contrast, the normal emission HeII UPS diers both from the early HeI work [6,7] and substantially from the more recent synchrotron work [3]. The peaks near the Fermi energy (EF ) at 1.4 and 2.9 eV match the two peaks at normal emission for HeI [6,7], but, as would be expected, in all other respects the angleresolved work is far more structured. Dierences with the synchrotron work are more interesting. The Santoni et al. [3] spectra all display a consistently sharp feature about 0.5 eV below their valence band maximum (VBM). In their experimental band structure over CK this appears as a very ¯at band which does not correlate with their calculated structure. The authors assign it to photoelectrons scattered into normal emission from a band maximum along CK. Although sharp
features at about 1 eV below EF are evident in the early HeI angle-resolved work [6,7] they undergo shifts both in position and intensity. If the 41 eV spectrum for normal emission (Santoni et al. [3], their Fig. 2) is overlaid with the HeII spectrum in Fig. 4, aligning the peak at 1.4 eV with the VBM zero level, then there is a precise BE correlation with four of the main features at 5.8, 7.3, 8.5 and 12.9 eV. The intense forward scattered peak though is not reproduced. There would appear to be a discrepancy of 1.4 eV between our spectrometer Fermi level referenced values and the VBM of Santoni et al. [3]. The authors note that they found a subsequent charging error of 0.8 eV, but did not indicate whether it was corrected for in their Fig. 2. The core Pb 5d5=2 BE for galena is known to be invariant, independent of the carrier type for the PbS(1 0 0) surface [35], indicating pinning of the Fermi level. Neither the XPS nor UPS show any emission intensity at 17 eV. The densities of states (DOS) of bulk galena computed in the present study using HF, DFT and B3LYP treatments of exchange and correlation are presented in Fig. 5 and Table 3. We ®nd that computed properties such as the band gap, BE of
234
J. Muscat, C. Klauber / Surface Science 491 (2001) 226±238
Fig. 5. The computed densities of states of PbS using HF, LDA, GGA and B3LYP. The Fermi level has been set to 0 eV.
peaks and bandwidths are not signi®cantly in¯uenced by the presence or lack of d-functions on the S or Pb sites or the choice of the pseudopotential. However, DOS computed using DFT or B3LYP exhibit some signi®cant dierences to the HF results. The DOS computed within the HF approximation in the current study are identical to those reported by Mian et al. with the S 3s states at a BE
of about 17.4 eV (relative to the Fermi energy which has been set to 0 eV). This is as expected since the current HF DOS calculations have been performed with very similar computational conditions to the study by Mian et al. The DOS computed with DFT or B3LYP predict the 3s states at a BE of around 12.8 or 13.4 eV respectively (see Fig. 5). We have also compared the performance of the dierent theoretical approximations for the calculation of core states in S (see Table 3). The BE of the S 2p3=2 and S 2s have been measured (with respect to the Pd 5d5=2 at 18.74 eV) to be 160:9 0:1 and 225:3 0:1 eV respectively. The position of the S 2p states is most accurately reproduced by the B3LYP calculations which predict a value within 5 eV of experiment compared to HF theory which yields a BE that is about 10 eV too high and DFT about 10 eV too low. For the S 2s states, LDA or GGA underestimate the BE by over 20 eV; HF theory gives the closest agreement, overestimating by 10 eV with B3LYP underestimating by 15 eV. We note that in recent HF calculations of PbS [36], the BE of the S core states reported were corrected by applying an arbitrary energy shift to the BE setting the BE of the bulk S 2s and 2p states to the experimental values [37].
4. Discussion We have checked the in¯ence of the pseudopotential, basis set and the treatment of exchange and correlation on the computed properties of PbS. All the above approximations have a systematic eect on the computed structure as reported in Section 3.
Table 3 A summary of the main features of the computed and experimental DOS Technique
Band gap
S3p valence bandwidth
S3s peak
Pb6s/S3p
S2p
S2s
HF LDA/GGA conductor B3LYP HF (Mian) Experiment (Santoni et al. [3]) Experiment (current)
5.7
7.3 6 5.7 7.3 5
17.4 12.8 13.4 17.4 12 13
10.6 8.5 8.8 10.6 7.5 8.5
173 150 156
235 204 210
N/A 160.9
N/A 225.3
1.6 5.7
J. Muscat, C. Klauber / Surface Science 491 (2001) 226±238
For the computed bulk modulus, all three pseudopotentials result in values within the experimental range of 48±73 GPa. Calculations performed using the Hay±Wadt pseudopotentials consistently predict the lowest bulk moduli, around 10 GPa lower than with the Durand due to the tendency of the Hay±Wadt pseudopotential to yield larger cell volumes than with the Durand pseudopotential. The inclusion of d-functions in the basis set leads to very small changes in the computed bulk modulus. Finally, GGA and B3LYP give very similar bulk moduli whereas HF tends to underestimate and LDA overestimate, probably because the LDA yields the lowest cell volumes whereas HF predicts the highest. The comparison of calculated and experimental electronic properties is subject to a number of possible pitfalls. Not only does photoelectron spectroscopy probe the near surface bulk of the solid, usually to a greater extent than the surface (photon energy and detection angle dependent), but it is also a ®nal state spectroscopy. The latter means the relationship between calculated initial state BE values and experimentally observed values is clouded by the in¯uence of relaxation eects, both intra-atomic and inter-atomic, during the photoemission process. The question of core-level BE shifts between bulk and surface atoms or the surface core level shift (SCLS) is well known [38]. These are experimentally dicult to detect and of small enough magnitude to not have a signi®cant eect on the comparisons between the experimental and computed valence band in this work. For example, in the case of PbS, a SCLS of around 0.3 eV has been measured for the S 2p [39] with recent calculations of PbS surfaces supporting the existence of a small SCLS [36,31]. Calculations of the DOS of the clean (0 0 1) surface of PbS have been performed which show that the major features of the valence band, such as the valence band width and the positions of the S3s and Pb6s/S3p peaks are not signi®cantly changed from those of the bulk although the density at a given energy may vary somewhat [40]. The question of ®nal state relaxation eects is more important, but in terms of core level chemical shifts the initial state BE nonetheless dominates. As the photoemission process occurs in a
235
signi®cantly shorter time frame than the period of atomic vibrations we need only consider ®nal state electron relaxation. This in turn is a sum of three contributions, that from core and valence electrons of the emitting atom and from the inter- or extraatomic contribution [41]. For atoms in dierent environments only the latter two contributions should vary and the inter-atomic component should depend mainly on nearest neighbour dipole polarisabilities. An illustration of the size of such contributions can be obtained by comparing BE values for the same levels for a given element in a variety of compounds. Mansour [42] has acquired a series of standard photoelectron spectra for nickel and in a variety of forms; Ni, NiO, a-Ni(OH)2 , b-Ni(OH)2 , c-NiOOH, Ni2 O3 6H2 O and LiNiO2 . Across that series the formal oxidation state of the Ni atom is variously 0, 2 and 3 with the atom's environment also changing. The experimental BE is reported for the 2p1=2 , 2p3=2 , 3s and 3p levels for the whole series. As these are core levels the photoelectron from any level for one of the given compounds should experience approximately the same relaxation, i.e. measured DBE values between the levels would be zero in that case. Corrected for static charge the measured BE varies over the compound series by about 3 eV for the 2p levels and about 2 eV for the 3s and 3p levels. Empirically the DBE values between levels for the series ranges from 0.25 eV to about 1 eV. This implies that the relaxation eect can vary with the level involved, even for core electrons, and can probably reduce the initial state determined BE shift from 10% to 50% of that expected from the change in chemical state. This is a signi®cant eect when BE changes are considered, but only a minor eect terms of absolute initial state BE as determined by ab initio calculations. In terms of the approximations made in the ab initio calculations, the computed DOS is relatively insensitive to the basis set and the choice of pseudopotential. Tests reveal that the dominant eect on the DOS is due to the treatment of electronic exchange. We ®nd that calculations performed within the HF approximation which treats the non-local exchange interaction exactly, but neglects correlation, result in DOS that are in poor
236
J. Muscat, C. Klauber / Surface Science 491 (2001) 226±238
agreement with experiment with respect to the position of the S 3s peak which is overbound by about 4 eV compared to experiment. However, calculations using DFT within the LDA or GGA which only treat the local part of the exchange interaction (and its gradient in the case of the GGA) or B3LYP (which is a mixture of 20% nonlocal HF exchange and 80% LDA and GGA exchange) predict the position of the S 3s peak to within 1 eV of the experimental value. In Fig. 6, we present the DOS computed with a variety of treatments of exchange starting o with the B3LYP hybrid functional and varying the proportion of non-local exchange from 0% (pure DFT exchange) to 100% (pure HF). It is clear that the position of the S 2p and S 3s peaks is governed by the proportion of non-local HF to DFT exchange used. A ratio of 20% HF and 80% DFT exchange, as implemented in the B3LYP functional yields a description of the valence band in excellent agreement with experiment. Core states such as the
S 2p are also predicted to within about 5 eV with respect to experiment using B3LYP; by comparison, HF overstimates and LDA or GGA understimate the BE of the S 2p peak by over 10 eV. The treatment of electronic correlation, which has a signi®cant in¯uence on the computed structure [4], has an insigni®cant eect on the predicted DOS. DFT calculations performed with complete neglect of electronic correlation yield DOS that are indistinguishable from DOS computed with correlation. 5. Conclusions In the current study, we have performed a systematic investigation of the physical, elastic and electronic properties of PbS and have isolated the in¯uence of the basis set, pseudopotential and the treatment of exchange and correlation on these properties.
Fig. 6. The in¯uence of exchange on the computed binding energy of the sulphur 3s (top panel) and 2p (bottom panel) states of PbS.
J. Muscat, C. Klauber / Surface Science 491 (2001) 226±238
In the current study, we have used basis sets which have already been well tested on PbS. Augmenting the Pb basis set does not make a signi®cant dierence to the results whereas improving the S basis set leads to a small change in the total energy and cell volume. Other properties such as the bulk modulus or the computed DOS are not aected by the basis set. We ®nd that the choice of pseudopotential has a systematic eect on the computed cell volume and bulk modulus but an insigni®cant eect on the computed DOS. The cell volume increases for each type of exchange and correlation used in the order; Durand, Stoll±Preuss and Hay±Wadt. Similarly, bulk moduli computed with the Durand pseudopotential are the highest, those due to the the Hay± Wadt the lowest with the Stoll±Preuss lying in between. The cell volumes and bulk moduli predicted with HF, LDA and GGA treatments of exchange and correlation follow the well recognised trends, notably that HF and GGA tend to overestimate cell volumes with respect to LDA. Hybrid calculations using the B3LYP functional produce resuts that lie between the HF and GGA values. In the current study, we ®nd that the computed DOS for PbS is dramatically aected by the treatment of exchange employed. The current experimental work con®rms that the BE of the S 3s state lies at 13 eV, in agreement with previous experimental work. We reveal that the reason for the discrepancy between previous experimental and theoretical determinations of the electronic structure of PbS is due to the treatment of the exchange interaction. Our calculations show that exact HF non-local exchange tends to overestimate the BEs of the valence and core states yielding DOS in poor agreement with experiment whereas approximate DFT exchange, based on either the LDA or the GGA tends to underestimate BE. We ®nd that the hybrid B3LYP functional, which contains a mixture of exact non-local exchange and DFT exchange gives estimates of BE in excellent agreement with experimental observations. Our results show that DFT or hybrid functional methods give a good description of PbS and form the starting point for further studies of the surface structure and chemistry of this mineral.
237
Acknowledgements The calculations were performed on a Compaq XP-1000 workstation and parallel PC cluster [43]. J.M. would like to thank J. Gale and N. Harrison for useful discussions regarding this work.
References [1] V.C. Padaki, S.T. Lakhshmikumar, S.V. Subrahmanyam, E.S.R. Gopal, Pramana 71 (1981) 25. [2] P.B. Littlewood, J. Phys. C 13 (1980) 4855. [3] A. Santoni, G. Paolucci, G. Santoro, K.C. Prince, N.E. Christensen, J. Phys.: Condens. Matter 4 (1992) 6759. [4] M. Mian, N.M. Harrison, V.R. Saunders, W.R. Flavell, Chem. Phys. Lett. 257 (1996) 627. [5] P.H. Citrin, D.R. Haman, Chem. Phys. Lett. 22 (1973) 301. [6] T. Grandke, L. Ley, M. Cardona, Phys. Rev. Lett. 38 (1977) 1033. [7] T. Grandke, L. Ley, M. Cardona, Phys. Rev. B 18 (1978) 3847. [8] A.D. Becke, J. Chem. Phys. 98 (1993) 5648. [9] P.J. Stephens, F.J. Devlin, C.F. Chabalowski, M.J. Frisch, J. Phys. Chem. 98 (1994) 11623. [10] J. Muscat, A. Wander, N.M. Harrison, Chem. Phys. Lett. 342 (2001) 397. [11] N. Wilson, J. Muscat, Mol. Sim., submitted for publication. [12] R. Dovesi, V.R. Saunders, C. Roetti, M. Causa, N.M. Harrison, R. Orlando, E. Apra, C R Y S T A L 9 5 User's Manual, University of Turin, Turin, 1996. [13] W. Kohn, L.J. Sham, Phys. Rev. 140 (1965) A1133. [14] D.M. Ceperley, B.J. Alder, Phys. Rev. Lett. 45 (1980) 566. [15] J.P. Perdew, A. Zunger, Phys. Rev. B 23 (1981) 5048. [16] J.P. Perdew, K. Burke, M. Ernzerhof, ACS Symposium series 629 (1996) 453. [17] A.D. Becke, Phys. Rev. A 38 (1988) 3098. [18] C. Lee, W. Yang, R.G. Parr, Phys. Rev. B 37 (1988) 785. [19] S.H. Vosko, L. Wilk, M. Nusair, Can. J. Phys. 58 (1980) 1200. [20] J. Muscat, Ph.D. Thesis, University of Manchester, Manchester, 1999. [21] P.J. Hay, W.R. Wadt, J. Chem. Phys. 82 (1985) 299. [22] J.C. Barthelat, P. Durand, Gazz. Chim. Ital. 108 (1978) 225. [23] The Stoll-Preuss pseudopotential library, http://www. theochem.uni-stuttgart.de, 2000. [24] B. Metz, H. Stoll, M. Dolg, J. Chem. Phys. 113 (2000) 2563. [25] J.D. Pack, H.J. Monkhorst, Phys. Rev. B 16 (1977) 1748. [26] C. Pisani, R. Dovesi, C. Roetti, Hartree±Fock Ab Initio Treatment of Crystaline Systems, vol. 48, Springer, Berlin, 1988. [27] R. Dovesi, C. Roetti, C. Freyria-Fava, E. Apra, V.R. Saunders, N.M. Harrison, Phil. Trans. R. Soc. Lond. A 341 (1992) 203.
238
J. Muscat, C. Klauber / Surface Science 491 (2001) 226±238
[28] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, L-BFGS-B ± Fortran Subroutines for Large Scale Bound Constrained Optimisation, Department of Electrical Engineering and Computer Science, Northwestern University, Illinois, 1994. [29] F.D. Murnaghan, Proc. Natl. Acad. Sci. USA 30 (1944) 244. [30] M.P. Seah, I.S. Gilmore, G. Beamson, Surf. Interf. Anal. 26 (1998) 642. [31] J. Muscat, J. Gale, 2001, in preparation. [32] C. Klauber, Surf. Interf. Anal. 20 (1993) 703. [33] R.T. Poole, J. Liesagang, R.C.G. Leckey, J.G. Jenkin, J. Electron Spectrosc. Relat. Phenom. 5 (1974) 773. [34] J.H. Sco®eld, J. Electron Spectrosc. Relat. Phenom. 8 (1976) 129. [35] T. Grandke, M. Cardona, Surf. Sci. 92 (1980) 385.
[36] U. Becker, M.F. Hochella Jr., Geochim. Cosmochim. Acta 60 (1996) 2413. [37] U. Becker, Private communication. [38] P.S. Bagus, C.R. Brunble, G. Pacchioni, F. Parmigiani, Surf. Sci. Reports 19 (1993) 265. [39] J.A. Leiro, K. Laajalehto, I. Kartio, M.H. Heinonen, Surf. Sci. 412±413 (1998) L918. [40] U. Becker, K.M. Rosso, Am. Miner. 86 (2001) 862. [41] G. Moretti, J. Electron Spectrosc. Relat. Phenom. 95 (1998) 95. [42] A.N. Mansour, Surf. Sci. Spectra 3 (1996) 211±286. [43] J. Muscat, W. Burrage, Low Cost Unix Type Multiprocessing System, CSIRO Minerals, DMR no. 1684, Melbourne, 2001.
www.elsevier.com/locate/susc
A combined ab initio and photoelectron study of galena (PbS) Joseph Muscat a,*, Craig Klauber b a
CSIRO Minerals, Box 312, Clayton South, Vic., 3169, Australia CSIRO Minerals, P.O. Box 90, Bentley, WA, 6982, Australia
b
Received 21 February 2001; accepted for publication 29 June 2001
Abstract The results of a study of the physical, elastic and electronic properties of PbS using both ®rst principles theoretical and experimental X-ray and ultraviolet photoelectron techniques are presented. We have investigated the in¯uence of the computational approximations such as the basis set, pseudopotential and the treatment of exchange and correlation on the accuracy of the theoretical results. The current study work con®rms that the binding energy of the S 3s state lies at 13 eV, in agreement with previous experimental work and contrary to earlier ab initio Hartree±Fock (HF) calculations. We establish that the reason for the discrepancy between previous experimental and theoretical determinations of the electronic structure of PbS is due to the treatment of the exchange interaction. We demonstrate that calculations using the HF approximation tend to overestimate the binding energies of the valence and core states yielding densities of states in poor agreement with experiment whereas calculations performed with density functional theory within either the local density or the generalised gradient approximation give a signi®cantly better description of the electronic structure. Finally, we con®rm that hybrid functional techniques such as B3LYP provide an excellent description of the physical and electronic structure of PbS. Ó 2001 Elsevier Science B.V. All rights reserved. Keywords: Ab initio quantum mechanical methods and calculations; Low index single crystal surfaces; Density functional calculations; Surface electronic phenomena (work function, surface potential, surface states, etc.); Sulphides; X-ray photoelectron spectroscopy; Visible and ultraviolet photoelectron spectroscopy
1. Introduction Lead remains an important industrial commodity, sourced primarily from galena (PbS) with minor contributions from cerussite (PbCO3 ) and anglesite (PbSO4 ). Prior to smelting, the galena component in natural ores is generally concentrated or bene®ciated via froth ¯otation. Galena,
* Corresponding author. Tel.: +61-39545-8500; fax: +6139562-8919. E-mail address: [email protected] (J. Muscat).
when freshly crushed, exposes naturally hydrophobic surfaces that readily attach to ¯otation bubbles, hence galena represents the classic example of bene®ciation via ¯otation. Weathering and other process variables can reduce this hydrophobicity and so the surface chemistry of galena, particularly its oxidation and the role played by collectors such as xanthates (to increase hydrophobicity) are key factors in a fundamental understanding of the ¯otation process. Considerable empirical eort over several decades has been directed toward an understanding of the surface chemistry of galena, in particular the predominant
0039-6028/01/$ - see front matter Ó 2001 Elsevier Science B.V. All rights reserved. PII: S 0 0 3 9 - 6 0 2 8 ( 0 1 ) 0 1 4 0 8 - X
J. Muscat, C. Klauber / Surface Science 491 (2001) 226±238
(1 0 0) cleavage plane. However, to eciently progress the understanding of such a mineral system by means of molecular modelling the ®rst essential step is to better understand the electronic structure. Once this baseline understanding has been achieved it should be possible to progress to modelling aspects such as the galena±xanthate adsorption or the oxidative reaction of the galena surface to a range of dierent valence states such as elemental sulphur, lead sulphate and other surface speci®c species. Galena is one of the most studied sulphide minerals due in part to the abundant availability of good quality single mineral samples and its simple NaCl-halite structure. There have been a number of experimental studies of the structural, elastic and electronic properties of PbS but many of these basic properties are not known accurately. For example, considerable uncertainty remains in the bulk modulus of PbS despite a number of measurements, with two of the most recent studies giving values of 48 GPa [1] and 73 GPa [2] respectively. Doubt has also been expressed as to the electronic structure PbS. This has recently been probed by photoemission experiments [3] and also calculated theoretically in subsequent ab initio simulations using Hartree±Fock (HF) theory [4]. The two studies yielded results that were generally in broad agreement with the exception of the determination of the binding energy (BE) of the S 3s band which was measured to lie at around 12 eV but predicted by the simulations to be about 17 eV. Mian et al. [4] note that the experimental data from Santoni et al. [3] did not extend to include the 17 eV BE region and suggest that the peak at about 12 eV could be due to HS-type species formed due to contamination of the surface by hydrogen. The photoelectron spectroscopies of X-ray photoelectron (XPS) and ultraviolet photoelectron (UPS) are excellent methods for obtaining core and valence level BE information for occupied levels. Both are ®nal state spectroscopies, hence the actual observed BEs are altered by ®nal state relaxation eects shielding the core hole. Such extra-atomic relaxation is almost independent of the core hole (e.g. spectra from noble gas atoms implanted in noble metals [5]) so that the relative
227
BE values are useful parameters to check the validity of electronic calculations. Earlier angle resolved HeI UPS work on PbS shows sharp structure in the 11±13 eV range which becomes less apparent at the higher photon energies probed in Santoni et al. [3]. In the later work the strongly dispersing feature at 12 eV is assigned by the authors to S 3s, at considerable variance to the HF computed value of Mian et al. [4] of 17 eV. Unfortunately, in the region 12±17 eV and below, artefacts from the steeply rising secondary electron emission (SEE) background, close to the vacuum zero, restricts the usefulness of the low UV energies in determining a reliable S 3s BE. This is well recognised in the earlier HeI work on galena [6,7] in which discussion of the electronic structure is limited to above 6 eV. In these studies, simple angle integrated HeII and MgKa photoemission spectra at normal emission are used to reassess the structure from the Pb 5d5=2 line to the Fermi level, in particular the S 3s. Since the work of Mian et al. [4], there have been signi®cant improvements in ab initio techniques, in particular the development and implementation of hybrid functional techniques [8,9] that combine elements of HF and density functional theory (DFT). These hybrid techniques have been shown to generally predict the electronic structure [10] and the structural and elastic properties [11] of a range of materials signi®cantly more reliably than the previous HF or DFT methods in common use. In the current study, we have performed ab initio calculations together with photoemission experiments of the valence band in order to clarify the situation as regards to the electronic structure of PbS. The calculations of electronic structure and the physical and elastic properties of PbS have been performed using the new three parameter B3LYP hybrid functional method [8,9] and the more traditional HF and DFT approximations in order to assess the performance of the dierent theoretical techniques and to enable direct comparisons with previous theoretical calculations. The in¯uence of other approximations inherent in ab initio calculations such as the basis set and the pseudopotential on the computed properties have also been examined as a prerequisite to future
228
J. Muscat, C. Klauber / Surface Science 491 (2001) 226±238
studies of the surface structure and chemistry of PbS surfaces.
2. Method 2.1. Computational details All the calculations have been performed using the C R Y S T A L 9 8 ab initio software package [12] which is based on the periodic, linear combination of atomic orbitals (LCAO) formalism where the Bloch orbitals of the crystal are expanded using atom centred Gaussian orbitals with s, p or d symmetry. A variety of self-consistent treatments of exchange and correlation have been used including HF, the LDA to DFT, the GGA to DFT and the B3LYP functional based on Becke's three parameter hybrid functional (B3LYP) employing a combination of HF and DFT exchange [8]. The DFT±LDA calculations were performed using LDA exchange [13] and the correlation functional parametrised from Monte-Carlo calculations [14] by Perdew and Zunger [15] whereas the GGA calculations employed the exchange and correlation functionals by Perdew et al. [16]. The B3LYP hybrid functional used in this study is a modi®ed version of Becke's original three-parameter exchange-correlation functional [8] as suggested by Stephens et al. [9] where the exchange and correlation is given by: EXC 1
a0 EXLDA a0 EXHF ax DEXB88
ac ECLYP 1
ac ECVWN
1
where ax DEXB88 is Becke's gradient correction to the exchange functional [17], and the correlation is given by a combination of the Lee, Yang and Parr (LYP) [18] and the Vosko et al. [19] functionals. Becke's original three empirical parameters (a0 0:2, ax 0:72 and ac 0:81) are used which were found to optimise the atomisation energies, ionisation potentials and proton anities of a number of small molecules [8]. Solid state calculations using this functional have been shown to yield band gaps in excellent agreement with experiment [10].
Due to the computational expense required to describe the large number of electrons and the problems associated with treating relativistic effects in core shells for heavy elements such as lead, the core electrons up to and including the 5p shell have been replaced with a pseudopotential. Recent studies have however shown that the pseudopotential approximation can be a source of signi®cant error [20]. We have examined the in¯uence of the pseudopotential on the computed results by performing calculations using the Hay and Wadt [21], the Durand [22] and the recently released Stoll and Preuss type pseudopotentials [23,24]. The sulphur ions have been treated using an all-electron basis set. In the current work, we have used basis sets that have been developed and optimised for use in recent studies of PbS [4]. We have performed tests using HF, DFT and B3LYP to determine if augmentation of the Pb or S basis sets by the addition of d-symmetry functions (denoted by the ``'' superscript) to the Pb or S basis set has a signi®cant in¯uence on the computed cell parameters or total energy. We ®nd that addition of d-functions to the S ions has a small, but signi®cant in¯uence on the predicted cell parameters and total energy but similar improvements to the Pb basis set make little dierence to the computed results (see Fig. 1). Calculations performed using the S basis set yield lattice parameters that are smaller than systematically about 0.03±0.06 A those due to the S basis set whereas structures computed using the Pb basis set dier from those . Furdue to the Pb basis set by less than 0.02 A thermore, the total energy computed using the S basis set is signi®cantly lower than that due to the S basis set (by 0.8±1.9 eV) whereas the dierence between the Pb and Pb basis sets is much smaller (around 0.03 eV). Since the use of the Pb basis set causes an increase in the computational cost of the calculations but only very small dierences in the computed results compared to the Pb basis set, the discussion in the current article will focus on the results calculated using the PbS basis set. Sampling of K-space has been performed using Pack±Monkhorst grids [12,25] of shrinking parameter eight, yielding 29 symmetry inequivalent K-points for bulk PbS. Tests reveal that increasing the number of K-points produces no signi®cant
J. Muscat, C. Klauber / Surface Science 491 (2001) 226±238
229
Fig. 1. The eect of the basis set and pseudopotential on the computed cell parameters.
dierence in the computed structure or energy of the crystals (the total energy of the bulk crystals is converged to within about 1 10 5 eV). C R Y S T A L 9 8 computes the matrix elements of the Coulomb and exchange terms by direct summation of the in®nite periodic lattice. The truncation of these summations is controlled by ®ve Gaussian overlap criteria, details of the control of these parameters is available elsewhere [12,26]. The values of the overlap criteria chosen in the current study were high (ITOLS parameters set to 10 8 , 10 8 , 10 8 , 10 8 and 10 16 ) in order to converge numerical errors to about 3 meV in the total energy [27]. The structural optimisations were converged to in cell parameters and 10 4 a tolerance of 0.01 A eV in the total energy using a modi®ed Broyden± Fletcher±Goldfarb±Shanno minimisation algorithm [28]. We have computed the bulk modulus of
PbS by ®tting a computed energy volume curve consisting of at least 30 points to the Murnaghan equation of state [29]. 2.2. XPS and UPS methodology Valence band spectra, excited by MgKa X-rays and the HeII UV line were obtained from two galena samples, one from Zacatecas (Mexico) and the other from Broken Hill (Australia). Samples were attached to stainless steel sample stubs by silver-dag and cleaved under a high purity nitrogen atmosphere in a glove-box attached directly to the spectrometer entry lock. The Broken Hill sample exhibited far better cleavage along the (1 0 0) face and only that data is reported here, although the Zacatecas spectra were identical. Both XPS and UPS data were obtained using a VG ESCALAB Mk II (nominal base pressure of about
230
J. Muscat, C. Klauber / Surface Science 491 (2001) 226±238
1 10 10 mbar) retro®tted with multiple channeltron (5) electron detection and an SSI AlKa monochromated source in addition to the original MgKa and AlKa twin anode source. The MgKa (1253.6 eV) spectra were obtained at a source power of 300 W irradiated over the whole sample and the monochromated AlKa (1486.6 eV) spectra (Pb 5d only) at 100 W from a nominal incident spot of 250 1000 lm. In both cases 6 mm analyser slits and a constant analyser pass energy of 20 eV were employed. The UV source used was a VG UVL-HI dierentially pumped discharge lamp which produced a stable HeII output with a main chamber helium pressure of 1 10 9 mbar. UPS spectra were obtained at a constant analyser retard ratio of 10. The spectrometer was calibrated using the Ag 3d5=2 line at 368.22 eV, the Au 4f7=2 line at 83.95 eV and a wide range linearity check using the Cu 2p line [30]. In all cases argon ion etched polycrystalline foils were used for calibration. Final adjustment for the XPS valence region of the galena was against the Au line which yielded a Pb 5d5=2 BE of 18:71 0:05 eV with respect to the Fermi level. In the case of the UPS valence region adjustment was with respect to the Fermi edge of a clean nickel foil which gave a Pb 5d5=2 BE of 18:74 0:05 eV. As the UPS established nickel Fermi edge zero was considered to be more precise the XPS valence region was shifted to align the Pb 5d5=2 peaks. 3. Results 3.1. Structural properties In Table 1, we present the optimal cell parameter for PbS computed with HF, LDA, GGA and B3LYP treatments of exchange and correla-
tion using the PbS basis set. The in¯uence of the pseudopotential on the computed cell parameter is signi®cant. Calculations performed using the Hay±Wadt pseudopotential predict lattice para larger meters that are systematically about 0.1 A than those computed using the Durand pseudopotential, whereas the Stoll±Preuss pseudopotential yields cell parameters which are typically in between those computed by the Hay±Wadt and the Durand pseudopotentials. These dierences are shown in Fig. 1. For comparison, results computed using the PbS, Pb S, PbS and Pb S basis sets are also included which show that the dierences in the computed cell parameter due to the pseudopotential are consistent across a range of basis sets. The treatment of exchange and correlation also has a systematic eect on the computed cell parameter, as can be seen from Table 1. The HF and B3LYP approximations overestimate the experimentally measured cell volume by around 1.9± 3.7% (HF) and 0.4±2.4% (B3LYP) depending on the pseudopotential used. The LDA to DFT underestimates the cell parameter and the GGA corrects for this eect leading to an increase in the cell parameter of about 1% compared to the LDA results. Similar trends arising due to the treatment of exchange and correlation have been reported elsewhere [20]. There is large uncertainty in the literature as to the magnitude of the bulk modulus of PbS, with the most recent measurements giving values ranging from 48 to 73 GPa. The computed bulk moduli, presented in Table 2 are all within this experimental range. The computed bulk modulus is very sensitive to the pseudopotential used; calculations performed using the Hay±Wadt pseudopotential systematically predict the bulk modulus to be approximately 10 GPa lower than
Table 1 The computed cell parameters of PbS. Figures in parentheses indicate the percentage deviation from the experimental cell parameter Durand Hay±Wadt Stoll±Preuss Experiment
HF
LDA
GGA
B3LYP
6.045 (1.85) 6.152 (3.66) 6.076 (2.37)
5.824 ( 1:87) 5.907 ( 0.47) 5.829 ( 1:79)
5.898 ( 0:62) 5.979 (0.74) 5.843 ( 1:55)
5.958 (0.39) 6.080 (2.44) 5.990 (0.93)
5.935
J. Muscat, C. Klauber / Surface Science 491 (2001) 226±238 Table 2 The computed bulk moduli for PbS compared to experimental measurements PS
HF
LDA
GGA
B3LYP
Durand Hay±Wadt Stoll±Preuss Experiment
61 51 55
73 62 72
63 58 60
60 50 56
48±73
that computed using the Durand pseudopotential. The Stoll±Preuss pseudopotential predicts the bulk modulus to be between these two extremes, similar to the trend reported for the computed cell parameter. HF and B3LYP predict very similar values for the bulk modulus of about 50 (Hay±Wadt), 55 (Stoll±Preuss) or 60 GPa (Durand). The GGA predicts bulk moduli that are slighly higher than the HF or B3LYP results whereas the LDA yields bulk moduli that are signi®cantly higher, over 10 GPa larger than the HF or B3LYP values. We have also performed tests to determine if the basis set has an in¯uence on the predicted bulk modulus and have found that the presence of (or lack of) d-functions makes little dierence to the computed bulk modulus.
231
3.2. Electronic structure As noted in Section 1, there is some uncertainty in the literature regarding the BE of the sulphur 3s peak. Recent UPS studies [3] observed a band of states at a BE of around 13 eV which is attributed to the S 3s states in bulk galena. Subsequent HF calculations [4] of galena (PbS) predicted the position of this peak to be 17 eV. The authors of the HF study suggested that the dierence between the experimental and calculated peak positions may be explained by considering the eect that contamination of the PbS surface with hydrogen (which is almost unavoidable in ultra-high vacuum experiments) could have on the photoemission spectrum. Adsorption of H onto the surface of Galena would lead to the formation of HS or H2 S-like species with a lower BE. Mian et al. postulated that the peaks observed by the UPS experiments were actually due to these surface HS-type species and that the states due to the S 3s electrons in galena actually lie closer to the computed BE of 17 eV. As the Santoni experiments did not probe BEs higher than 14 eV, it was not possible to validate this hypothesis. However, as the cross-section for H 1s
Fig. 2. MgKa XPS of the galena valence band region.
232
J. Muscat, C. Klauber / Surface Science 491 (2001) 226±238
Fig. 3. HeII UPS of the galena valence band region. Pb 5d ghost peaks due to the 48.37 eV radiation overlap part of the valence region. Dotted lines show the results from alternative treatments for removal of the Pb ghost peaks.
photoemission is very small, it is unlikely for it to be mistaken for the S 3s state, as supported by calculations of the adsorption of H2 on the PbS (0 0 1) surface [31]. The MgKa XPS spectrum is shown in Fig. 2 and the HeII UPS spectrum in Fig. 3. Both processed curves are compared in Fig. 4. Removal of MgKa X-ray line and satellite structure is straightforward [32], but correction for the He ghost line is not. As noted, a reliable determination in the S 3s region requires a photon energy above that of HeI radiation (21.22 eV). Maximising discharge lamp intensity for the more desirable HeII at 40.81 eV has the disadvantage of strengthening other minor resonance lines. Based on the valence band edges for a clean nickel foil sample the minor lines at 23.09, 23.74 and 48.37 eV [33] were detectable under lamp conditions for optimising HeII intensity. The proximity of the Pb 5d doublet to the valence band means that ghost peaks from the 48.37 eV line will overlap the bottom of the valence band structure. This is evident in the raw data of Fig. 3. Attempted removal of the ghost lines by deconvolution of a source function (as with the MgKa spectra) represented by the 40.81 and 48.37 eV lines proved unsatisfactory due to the high emis-
sion intensity induced at low kinetic energies by the dominant HeI. As the extent of valence band contamination from the 48.37 eV radiation is small (evidenced by the low intensity above the Fermi level in Fig. 3) the best method for removal of the ghost Pb 5d doublet was found to be replication of the HeII doublet and direct subtraction. The HeII Pb 5d5=2 peak was matched to a Voigt pro®le with FWHM of 0.537 eV. The position and intensity of the Pb 5d3=2 ghost peak can be ascertained and so removal of the Pb 5d5=2 ghost peak, assuming the same peak shape and width, simply requires a splitting and intensity ratio. The apparent splitting of the main peaks (D) with HeII is 2.4 eV (background distorted), for the ghost peaks 2.62 eV, for MgKa 2.54 eV and monochromatic AlKa photoemission 2.57 eV. The expected statistical intensity ratio (r) of 1.5 compares to 1.45 for MgKa , 1.38 from the monochromatic AlKa photo emission and 1.49 and 1.42 calculated for MgKa and AlKa [34]. With spectrometer transmission (under CRR conditions) also having an eect the optimum description for the Pb 5d5=2 ghost peak is unclear. In Fig. 3 alternative corrections for the Pb lines are shown (based on D 2:62 eV, r 1:50 and D 2:57 eV, r 1:38) to highlight
J. Muscat, C. Klauber / Surface Science 491 (2001) 226±238
233
Fig. 4. Comparison of the valence band photoemission spectra from MgKa and HeII sources for cleaved galena.
that signi®cant intensity remains at about 11 eV, irrespective of how the removal is performed. The ®rst correction based on observed ghost peak splitting and the statistical ratio is preferred. The processed valence band photoemission spectra from both MgKa and HeII sources for cleaved galena are compared in Fig. 4. The XPS is very similar to earlier AlKa work [6,7] although slightly more structure is displayed. In contrast, the normal emission HeII UPS diers both from the early HeI work [6,7] and substantially from the more recent synchrotron work [3]. The peaks near the Fermi energy (EF ) at 1.4 and 2.9 eV match the two peaks at normal emission for HeI [6,7], but, as would be expected, in all other respects the angleresolved work is far more structured. Dierences with the synchrotron work are more interesting. The Santoni et al. [3] spectra all display a consistently sharp feature about 0.5 eV below their valence band maximum (VBM). In their experimental band structure over CK this appears as a very ¯at band which does not correlate with their calculated structure. The authors assign it to photoelectrons scattered into normal emission from a band maximum along CK. Although sharp
features at about 1 eV below EF are evident in the early HeI angle-resolved work [6,7] they undergo shifts both in position and intensity. If the 41 eV spectrum for normal emission (Santoni et al. [3], their Fig. 2) is overlaid with the HeII spectrum in Fig. 4, aligning the peak at 1.4 eV with the VBM zero level, then there is a precise BE correlation with four of the main features at 5.8, 7.3, 8.5 and 12.9 eV. The intense forward scattered peak though is not reproduced. There would appear to be a discrepancy of 1.4 eV between our spectrometer Fermi level referenced values and the VBM of Santoni et al. [3]. The authors note that they found a subsequent charging error of 0.8 eV, but did not indicate whether it was corrected for in their Fig. 2. The core Pb 5d5=2 BE for galena is known to be invariant, independent of the carrier type for the PbS(1 0 0) surface [35], indicating pinning of the Fermi level. Neither the XPS nor UPS show any emission intensity at 17 eV. The densities of states (DOS) of bulk galena computed in the present study using HF, DFT and B3LYP treatments of exchange and correlation are presented in Fig. 5 and Table 3. We ®nd that computed properties such as the band gap, BE of
234
J. Muscat, C. Klauber / Surface Science 491 (2001) 226±238
Fig. 5. The computed densities of states of PbS using HF, LDA, GGA and B3LYP. The Fermi level has been set to 0 eV.
peaks and bandwidths are not signi®cantly in¯uenced by the presence or lack of d-functions on the S or Pb sites or the choice of the pseudopotential. However, DOS computed using DFT or B3LYP exhibit some signi®cant dierences to the HF results. The DOS computed within the HF approximation in the current study are identical to those reported by Mian et al. with the S 3s states at a BE
of about 17.4 eV (relative to the Fermi energy which has been set to 0 eV). This is as expected since the current HF DOS calculations have been performed with very similar computational conditions to the study by Mian et al. The DOS computed with DFT or B3LYP predict the 3s states at a BE of around 12.8 or 13.4 eV respectively (see Fig. 5). We have also compared the performance of the dierent theoretical approximations for the calculation of core states in S (see Table 3). The BE of the S 2p3=2 and S 2s have been measured (with respect to the Pd 5d5=2 at 18.74 eV) to be 160:9 0:1 and 225:3 0:1 eV respectively. The position of the S 2p states is most accurately reproduced by the B3LYP calculations which predict a value within 5 eV of experiment compared to HF theory which yields a BE that is about 10 eV too high and DFT about 10 eV too low. For the S 2s states, LDA or GGA underestimate the BE by over 20 eV; HF theory gives the closest agreement, overestimating by 10 eV with B3LYP underestimating by 15 eV. We note that in recent HF calculations of PbS [36], the BE of the S core states reported were corrected by applying an arbitrary energy shift to the BE setting the BE of the bulk S 2s and 2p states to the experimental values [37].
4. Discussion We have checked the in¯ence of the pseudopotential, basis set and the treatment of exchange and correlation on the computed properties of PbS. All the above approximations have a systematic eect on the computed structure as reported in Section 3.
Table 3 A summary of the main features of the computed and experimental DOS Technique
Band gap
S3p valence bandwidth
S3s peak
Pb6s/S3p
S2p
S2s
HF LDA/GGA conductor B3LYP HF (Mian) Experiment (Santoni et al. [3]) Experiment (current)
5.7
7.3 6 5.7 7.3 5
17.4 12.8 13.4 17.4 12 13
10.6 8.5 8.8 10.6 7.5 8.5
173 150 156
235 204 210
N/A 160.9
N/A 225.3
1.6 5.7
J. Muscat, C. Klauber / Surface Science 491 (2001) 226±238
For the computed bulk modulus, all three pseudopotentials result in values within the experimental range of 48±73 GPa. Calculations performed using the Hay±Wadt pseudopotentials consistently predict the lowest bulk moduli, around 10 GPa lower than with the Durand due to the tendency of the Hay±Wadt pseudopotential to yield larger cell volumes than with the Durand pseudopotential. The inclusion of d-functions in the basis set leads to very small changes in the computed bulk modulus. Finally, GGA and B3LYP give very similar bulk moduli whereas HF tends to underestimate and LDA overestimate, probably because the LDA yields the lowest cell volumes whereas HF predicts the highest. The comparison of calculated and experimental electronic properties is subject to a number of possible pitfalls. Not only does photoelectron spectroscopy probe the near surface bulk of the solid, usually to a greater extent than the surface (photon energy and detection angle dependent), but it is also a ®nal state spectroscopy. The latter means the relationship between calculated initial state BE values and experimentally observed values is clouded by the in¯uence of relaxation eects, both intra-atomic and inter-atomic, during the photoemission process. The question of core-level BE shifts between bulk and surface atoms or the surface core level shift (SCLS) is well known [38]. These are experimentally dicult to detect and of small enough magnitude to not have a signi®cant eect on the comparisons between the experimental and computed valence band in this work. For example, in the case of PbS, a SCLS of around 0.3 eV has been measured for the S 2p [39] with recent calculations of PbS surfaces supporting the existence of a small SCLS [36,31]. Calculations of the DOS of the clean (0 0 1) surface of PbS have been performed which show that the major features of the valence band, such as the valence band width and the positions of the S3s and Pb6s/S3p peaks are not signi®cantly changed from those of the bulk although the density at a given energy may vary somewhat [40]. The question of ®nal state relaxation eects is more important, but in terms of core level chemical shifts the initial state BE nonetheless dominates. As the photoemission process occurs in a
235
signi®cantly shorter time frame than the period of atomic vibrations we need only consider ®nal state electron relaxation. This in turn is a sum of three contributions, that from core and valence electrons of the emitting atom and from the inter- or extraatomic contribution [41]. For atoms in dierent environments only the latter two contributions should vary and the inter-atomic component should depend mainly on nearest neighbour dipole polarisabilities. An illustration of the size of such contributions can be obtained by comparing BE values for the same levels for a given element in a variety of compounds. Mansour [42] has acquired a series of standard photoelectron spectra for nickel and in a variety of forms; Ni, NiO, a-Ni(OH)2 , b-Ni(OH)2 , c-NiOOH, Ni2 O3 6H2 O and LiNiO2 . Across that series the formal oxidation state of the Ni atom is variously 0, 2 and 3 with the atom's environment also changing. The experimental BE is reported for the 2p1=2 , 2p3=2 , 3s and 3p levels for the whole series. As these are core levels the photoelectron from any level for one of the given compounds should experience approximately the same relaxation, i.e. measured DBE values between the levels would be zero in that case. Corrected for static charge the measured BE varies over the compound series by about 3 eV for the 2p levels and about 2 eV for the 3s and 3p levels. Empirically the DBE values between levels for the series ranges from 0.25 eV to about 1 eV. This implies that the relaxation eect can vary with the level involved, even for core electrons, and can probably reduce the initial state determined BE shift from 10% to 50% of that expected from the change in chemical state. This is a signi®cant eect when BE changes are considered, but only a minor eect terms of absolute initial state BE as determined by ab initio calculations. In terms of the approximations made in the ab initio calculations, the computed DOS is relatively insensitive to the basis set and the choice of pseudopotential. Tests reveal that the dominant eect on the DOS is due to the treatment of electronic exchange. We ®nd that calculations performed within the HF approximation which treats the non-local exchange interaction exactly, but neglects correlation, result in DOS that are in poor
236
J. Muscat, C. Klauber / Surface Science 491 (2001) 226±238
agreement with experiment with respect to the position of the S 3s peak which is overbound by about 4 eV compared to experiment. However, calculations using DFT within the LDA or GGA which only treat the local part of the exchange interaction (and its gradient in the case of the GGA) or B3LYP (which is a mixture of 20% nonlocal HF exchange and 80% LDA and GGA exchange) predict the position of the S 3s peak to within 1 eV of the experimental value. In Fig. 6, we present the DOS computed with a variety of treatments of exchange starting o with the B3LYP hybrid functional and varying the proportion of non-local exchange from 0% (pure DFT exchange) to 100% (pure HF). It is clear that the position of the S 2p and S 3s peaks is governed by the proportion of non-local HF to DFT exchange used. A ratio of 20% HF and 80% DFT exchange, as implemented in the B3LYP functional yields a description of the valence band in excellent agreement with experiment. Core states such as the
S 2p are also predicted to within about 5 eV with respect to experiment using B3LYP; by comparison, HF overstimates and LDA or GGA understimate the BE of the S 2p peak by over 10 eV. The treatment of electronic correlation, which has a signi®cant in¯uence on the computed structure [4], has an insigni®cant eect on the predicted DOS. DFT calculations performed with complete neglect of electronic correlation yield DOS that are indistinguishable from DOS computed with correlation. 5. Conclusions In the current study, we have performed a systematic investigation of the physical, elastic and electronic properties of PbS and have isolated the in¯uence of the basis set, pseudopotential and the treatment of exchange and correlation on these properties.
Fig. 6. The in¯uence of exchange on the computed binding energy of the sulphur 3s (top panel) and 2p (bottom panel) states of PbS.
J. Muscat, C. Klauber / Surface Science 491 (2001) 226±238
In the current study, we have used basis sets which have already been well tested on PbS. Augmenting the Pb basis set does not make a signi®cant dierence to the results whereas improving the S basis set leads to a small change in the total energy and cell volume. Other properties such as the bulk modulus or the computed DOS are not aected by the basis set. We ®nd that the choice of pseudopotential has a systematic eect on the computed cell volume and bulk modulus but an insigni®cant eect on the computed DOS. The cell volume increases for each type of exchange and correlation used in the order; Durand, Stoll±Preuss and Hay±Wadt. Similarly, bulk moduli computed with the Durand pseudopotential are the highest, those due to the the Hay± Wadt the lowest with the Stoll±Preuss lying in between. The cell volumes and bulk moduli predicted with HF, LDA and GGA treatments of exchange and correlation follow the well recognised trends, notably that HF and GGA tend to overestimate cell volumes with respect to LDA. Hybrid calculations using the B3LYP functional produce resuts that lie between the HF and GGA values. In the current study, we ®nd that the computed DOS for PbS is dramatically aected by the treatment of exchange employed. The current experimental work con®rms that the BE of the S 3s state lies at 13 eV, in agreement with previous experimental work. We reveal that the reason for the discrepancy between previous experimental and theoretical determinations of the electronic structure of PbS is due to the treatment of the exchange interaction. Our calculations show that exact HF non-local exchange tends to overestimate the BEs of the valence and core states yielding DOS in poor agreement with experiment whereas approximate DFT exchange, based on either the LDA or the GGA tends to underestimate BE. We ®nd that the hybrid B3LYP functional, which contains a mixture of exact non-local exchange and DFT exchange gives estimates of BE in excellent agreement with experimental observations. Our results show that DFT or hybrid functional methods give a good description of PbS and form the starting point for further studies of the surface structure and chemistry of this mineral.
237
Acknowledgements The calculations were performed on a Compaq XP-1000 workstation and parallel PC cluster [43]. J.M. would like to thank J. Gale and N. Harrison for useful discussions regarding this work.
References [1] V.C. Padaki, S.T. Lakhshmikumar, S.V. Subrahmanyam, E.S.R. Gopal, Pramana 71 (1981) 25. [2] P.B. Littlewood, J. Phys. C 13 (1980) 4855. [3] A. Santoni, G. Paolucci, G. Santoro, K.C. Prince, N.E. Christensen, J. Phys.: Condens. Matter 4 (1992) 6759. [4] M. Mian, N.M. Harrison, V.R. Saunders, W.R. Flavell, Chem. Phys. Lett. 257 (1996) 627. [5] P.H. Citrin, D.R. Haman, Chem. Phys. Lett. 22 (1973) 301. [6] T. Grandke, L. Ley, M. Cardona, Phys. Rev. Lett. 38 (1977) 1033. [7] T. Grandke, L. Ley, M. Cardona, Phys. Rev. B 18 (1978) 3847. [8] A.D. Becke, J. Chem. Phys. 98 (1993) 5648. [9] P.J. Stephens, F.J. Devlin, C.F. Chabalowski, M.J. Frisch, J. Phys. Chem. 98 (1994) 11623. [10] J. Muscat, A. Wander, N.M. Harrison, Chem. Phys. Lett. 342 (2001) 397. [11] N. Wilson, J. Muscat, Mol. Sim., submitted for publication. [12] R. Dovesi, V.R. Saunders, C. Roetti, M. Causa, N.M. Harrison, R. Orlando, E. Apra, C R Y S T A L 9 5 User's Manual, University of Turin, Turin, 1996. [13] W. Kohn, L.J. Sham, Phys. Rev. 140 (1965) A1133. [14] D.M. Ceperley, B.J. Alder, Phys. Rev. Lett. 45 (1980) 566. [15] J.P. Perdew, A. Zunger, Phys. Rev. B 23 (1981) 5048. [16] J.P. Perdew, K. Burke, M. Ernzerhof, ACS Symposium series 629 (1996) 453. [17] A.D. Becke, Phys. Rev. A 38 (1988) 3098. [18] C. Lee, W. Yang, R.G. Parr, Phys. Rev. B 37 (1988) 785. [19] S.H. Vosko, L. Wilk, M. Nusair, Can. J. Phys. 58 (1980) 1200. [20] J. Muscat, Ph.D. Thesis, University of Manchester, Manchester, 1999. [21] P.J. Hay, W.R. Wadt, J. Chem. Phys. 82 (1985) 299. [22] J.C. Barthelat, P. Durand, Gazz. Chim. Ital. 108 (1978) 225. [23] The Stoll-Preuss pseudopotential library, http://www. theochem.uni-stuttgart.de, 2000. [24] B. Metz, H. Stoll, M. Dolg, J. Chem. Phys. 113 (2000) 2563. [25] J.D. Pack, H.J. Monkhorst, Phys. Rev. B 16 (1977) 1748. [26] C. Pisani, R. Dovesi, C. Roetti, Hartree±Fock Ab Initio Treatment of Crystaline Systems, vol. 48, Springer, Berlin, 1988. [27] R. Dovesi, C. Roetti, C. Freyria-Fava, E. Apra, V.R. Saunders, N.M. Harrison, Phil. Trans. R. Soc. Lond. A 341 (1992) 203.
238
J. Muscat, C. Klauber / Surface Science 491 (2001) 226±238
[28] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, L-BFGS-B ± Fortran Subroutines for Large Scale Bound Constrained Optimisation, Department of Electrical Engineering and Computer Science, Northwestern University, Illinois, 1994. [29] F.D. Murnaghan, Proc. Natl. Acad. Sci. USA 30 (1944) 244. [30] M.P. Seah, I.S. Gilmore, G. Beamson, Surf. Interf. Anal. 26 (1998) 642. [31] J. Muscat, J. Gale, 2001, in preparation. [32] C. Klauber, Surf. Interf. Anal. 20 (1993) 703. [33] R.T. Poole, J. Liesagang, R.C.G. Leckey, J.G. Jenkin, J. Electron Spectrosc. Relat. Phenom. 5 (1974) 773. [34] J.H. Sco®eld, J. Electron Spectrosc. Relat. Phenom. 8 (1976) 129. [35] T. Grandke, M. Cardona, Surf. Sci. 92 (1980) 385.
[36] U. Becker, M.F. Hochella Jr., Geochim. Cosmochim. Acta 60 (1996) 2413. [37] U. Becker, Private communication. [38] P.S. Bagus, C.R. Brunble, G. Pacchioni, F. Parmigiani, Surf. Sci. Reports 19 (1993) 265. [39] J.A. Leiro, K. Laajalehto, I. Kartio, M.H. Heinonen, Surf. Sci. 412±413 (1998) L918. [40] U. Becker, K.M. Rosso, Am. Miner. 86 (2001) 862. [41] G. Moretti, J. Electron Spectrosc. Relat. Phenom. 95 (1998) 95. [42] A.N. Mansour, Surf. Sci. Spectra 3 (1996) 211±286. [43] J. Muscat, W. Burrage, Low Cost Unix Type Multiprocessing System, CSIRO Minerals, DMR no. 1684, Melbourne, 2001.