A Determination Method of the Work function using the Slab Model with a First-Principles Electronic Structure Calculation

We have reinvestigated the method to determine the work function of metals by the first-principles electronic structure calculation using the plane-wave expansion method with pseudo-potentials. Changing width of a vacuum layer of a slab model, we obtained the Fermi level measured from the effective Kohn-Sham potential at the center of the vacuum layer. By extrapolating the values of the Fermi level to the infinite vacuum-layer limit, the Fermi energy is determined. We propose a fitting function for the extrapolation. This function works well for determination of the highest occupied level and the lowest unoccupied level of materials with a surface relative to the vacuum level in the Kohn-Sham scheme. The obtained work functions for [100], [110] and [111] surfaces of three noble metals, Rh, Pd and Pt are comparable to the known theoretical estimation as well as the experimental values. [DOI: 10.1380/ejssnt.2008.103]


I. INTRODUCTION
Determination of the work function of metals using the slab model is possible, if a careful treatment of the vacuum layer as well as the metal layer is made.When we prepare a series of slab models by systematically changing the width of the vacuum layer, the difference between the Fermi level and the vacuum level has a well converged value in the limit of infinite vacuum layer.There were several researches on the quantum size effect coming from finite thickness of the metal layer [1][2][3][4][5][6].However, to our best knowledge, the finite-size behavior with respect to the thickness of the vacuum layer using the Kohn-Sham effective potential given by the density functional theory (DFT) [7,8] has been rarely discussed in the literature.This is due to slow convergence of the exchangecorrelation potential [5].
In this paper, we propose a fitting function to extrapolate the pseudo work function given in a finite slab model.The function is given by, where d is the thickness of the vacuum region.Three constants, a, b and c are fitting parameters.The Fermi level and the vacuum level are fitted by f (d).This function turns out to be applicable for estimation of the highest occupied level and the lowest unoccupied level of a material relative to the vacuum level.Thus this method should be applied for evaluation of the electron affinity of a semiconductor within the Kohn-Sham scheme except for difficulty in estimation of the band-gap within the local density approximation.As a test calculation, we performed estimation of the work function of several noble metals, Rh, Pd and Pt, which have the high catalytic activity.We adopted a first-principles calculation with a generalized gradient approximation [9].A slab model with a thin plate of a noble metal is considered.The plane-waveexpansion method with the ultrasoft pseudo-potentials * Corresponding author: ikuno@aquarius.mp.es.osaka-u.ac.jp is adopted.The actual calculation was done with the Quantum-espresso ver 3.0 package [10].The effective potential profile and the Fermi level are directly given by this computation software.
For a finite slab model, even at the center of the vacuum region, the full effective potential is slightly fluctuating on a plane parallel to the surface of a metal.We calculated the Fermi level relative to the averaged value of the effective potential on this plane.The length of the vacuum layer was changed from 12.5 Å to 100 Å.We fitted the obtained Fermi level to f (d) by optimizing three parameters a, b and c.Then the true Fermi energy and the vacuum level are estimated by taking the limit of d → ∞ to have f (∞) = c.The parameter c gives an estimation of the converged Fermi energy of the infinite vacuum layer.Using the Fermi energy, we evaluated the work function of noble metals, Pd, Rh and Pt.Obtained values were compared to the experimental values for several surfaces of these noble metals [12,13].
On the anisotropy of the work function with respect to the surface orientation, the most famous explanation may be the Smoluchowski model, which attributes the anisotropy to the smoothing of electron density at the surface [17].In order to discuss the Smoluchowski effect, the work functions were investigated for three different plane directions, [100], [110] and [111].We compared our results with other results obtained by the full-potential linear-muffin-tin orbital (FP-LMTO) method [14].
In the section II, we summarize the computational methods and conditions of calculations.In the section III, details of the present result are discussed.In the section IV, we summarize our conclusion.

II. CALCULATION METHOD
In the present calculation, the generalized gradient approximation (GGA) [9] was adopted.This is because the lattice constants of bulk noble metals were well reproduced in the self-consistent solutions of DFT.In the Quantum-espresso package, the Kohn-Sham orbital and the charge density were expanded in plane-wave basis sets.

Ikuno and Kusakabe
For the calculation of bulk metals, the kinetic energy cutoff for the orbitals and the energy cutoff for the charge density were set to be 25 Ry and 100 Ry, respectively.The Brillouin-zone (BZ) integration was performed with a k-mesh of 18 × 18 × 18.Typical values of relative errors in the optimized lattice constants from the experimental values were less than 2%.
For the slab calculations, the kinetic energy cutoff for the orbitals and the energy cutoff for the charge density expansion were set to 25 Ry and 220 Ry, respectively.The values were confirmed to be enough by checking convergence in the total energy as well as the structure parameters of lattices.The BZ integrations were performed with the Marzari-Vanderbilt cold smearing [15] with a smearing parameter of 0.03 Ry.The sampling k-mesh for the slab model was 8 × 8 × 1.
In the slab calculations, we have studied three noble metals, Rh, Pd and Pt, with three plane directions, [100], [110], and [111].These surfaces were modeled by repeated slabs consisting of five atomic layers and a vacuum layer, whose thickness d was chosen to be 12.5, 25, 50, 75 and 100 Å.Here, d denotes the spacing between the top atom of the metal layer to the bottom atom of the next repeated metal layer.As for the in-plane lattice constant, we used the bulk lattice constant determined by the selfconsistent-field calculation with GGA.For most of our slab calculations, the top two layers of slab were then optimized using the Broyden-Fletcher-Goldfarb-Shanno algorithm [16] and the other three layers were fixed.The optimization was done until the force on each atom was below 10 −3 [Ry]/a B [a.u.].This method to optimize surface layers only is adopted to reduce the total computation time.However, we have to be careful about induced dipole moments caused by asymmetric atomic structure of the slab.We will give discussion on errors coming from this set up for the slab model.Actually, we can show that the effect from this asymmetric structure is negligible by comparing the result to that for the corresponding symmetric slab model.We didn't consider the surface reconstruction of Pt found in low temperatures.
The work function φ is defined as the energy required to extract one electron from a metal surface and to let it go infinitely away from the surface.If we denote the total energy of the N electron system in a metal as E(N ), and if V vac is the electrostatic potential far from the metal surface, φ is given by E(N − 1) + V vac − E(N ).The value is known to be equivalent to [11], where E ∞ and E F are the true vacuum level and the Fermi level of the metal.In the Kohn-Sham scheme, the total effective potential is given as a summation of three terms as, Here v Hartree (r), v xc (r) and v ext (r) are the Hartree potential, the exchange-correlation potential and the external ionic potential, respectively.Using a finite slab, the pseudo work function has been estimated by the formula where E vac (d) is estimated by averaging the full effective potential v eff (r) on a plane parallel to the surface at

III. RESULTS AND DISCUSSION
Our calculated lattice constants of bulk Rd, Pd and Pt are listed in TABLE I. Our lattice constants overestimated the true lattice constants by 1-2 %.In order to show validity of our calculation, we compared the present result to a former simulation result and to experimental observation with respect to the work function φ, ∆d 12 , ∆d 23 , of the Rh[110] surface [18].Here, ∆d 12 or ∆d 23 are deviation of inter-layer distance between the 1st and http://www.sssj.org/ejssnt(J-Stage: http://www.jstage.jst.go.jp/browse/ejssnt/) e-Journal of Surface Science and Nanotechnology  2nd layers or the 2nd and 3rd layers from inter-layer distance in the bulk, and σ is the surface energy defined by σ = (E slab (5) − 5E bulk ) /2 for the present five layer system.Here E slab (5) and E bulk are the total energy of the slab model and the bulk, respectively.The results given in Tables I and II certify accuracy of the present results.We estimated the work function of three noble metals, Rh, Pd and Pt with three different plane directions, [100], [110] and [111] (see Table III).Our values are about 10% smaller than the corresponding values by Methfessel et al. [14].If the present result follows the Smoluchowski rule, the work function should align as φ([110]) < φ([100]) < φ ([111]).However, we see that the order is not reproduced for Pt.
To show that the essential feature does not depend strongly on the structure of the slab, we compared the converged Fermi energy in two cases.(See Table IV.)The first case is given by a slab with top two layers relaxed and the second is another slab with all layers relaxed.We have almost the same results irrespective of the way to optimize the structure, where the plane direction is one of [100], [110] and [111].
In the present work, we defined the vacuum level as the mean value of the total effective potential at the center of the vacuum region.There is another possibility to count E F (d) by letting the value of the electrostatic potential,  The same fitting function explains also the ddependence of the conduction band top and the valence band minimum of semiconductors, as tested in Hydrogenterminated hexagonal Boron Nitride structures.These examples showing that specified energy levels such as the highest occupied level and the lowest unoccupied level of a material always change following f (d) suggest that the effective potential v eff (r) itself behaves as f (d), when the inter slab distance d is increased.This behavior is found in self-consistent solutions of v eff (r), which is determined by the fully relaxed electron charge density.The density forming the double layer at each surface is affected by the next slab.We need to note that the total repeated slabs create the full v eff (r).

IV. CONCLUSION
By determining the extrapolated value of the Fermi energy in the limit of the infinite vacuum layer, we have redetermined the work function of noble metals and their dependence on the surface orientation.We found rough agreement between experimental values and the calculated ones.The obtained results for Pt, however, do not follow the Smoluchowski rule, althouth we did not considered the surface reconstruction.The fitting function of Eq. ( 1) captures change in the Fermi energy against the thickness of the vacuum layer d.

FIG. 1 : 4 )
FIG. 1: The pseudo Fermi energy of Rh [110] surface.The horizontal axis represents inverse of the vacuum layer width, 1/(d − b), while the vertical axis represents EF (d) − c.The parameters a = 136.25 [eV/ Å] and c = −4.6147[eV] as well as b = −5.4698[ Å] are determined by fitting the raw data shown by the circles to f (d).The linear behavior of the data certifies validity of the fitting function, f (d) = a/(d − b) + c shown by the dashed line.

TABLE I :
Calculated results of lattice constant for Rh, Pd, Pt.The difference is shown in percentage.

TABLE III :
Calculated results of work funciton

TABLE IV :
Structure dependence of the work function of Rh.Two slab models, one with two surface layers relaxed and the other with full relaxed layers, are compared.Hartree (r)+v ext (r), at the center of the vacuum region the origin of the energy.When the vacuum region becomes thick enough, the values given by these two definitions should be the same.The present result shows that E F (d) actually converges to the same value, even if we consider v Hartree (r) + v ext (r) in place of v eff (r).