Crystal electron binding energy and surface work function control of tin dioxide

Keith T. Butler,1 John Buckeridge,2 C. Richard A. Catlow,2 and Aron Walsh1,* 1Centre for Sustainable Chemical Technologies and Department of Chemistry, University of Bath, Claverton Down, Bath BA2 7AY, United Kingdom 2University College London, Kathleen Lonsdale Materials Chemistry, Department of Chemistry, 20 Gordon Street, London WC1H 0AJ, United Kingdom (Received 26 November 2013; revised manuscript received 13 March 2014; published 28 March 2014)


I. INTRODUCTION
The work function (φ) of a material is a critical parameter for determining the efficiency of charge transfer. As introduced by Bardeen, the work function depends on two independent quantities: (i) the binding energy of an electron in the bulk solid, sometimes termed the Galvani potential and (ii) the energy required to move the electron through an electrostatic double layer at the material surface [1]. The first is largely due to the electrostatic and bonding properties of the bulk material; the second is sensitive to surface structure, composition, and environment [2]. Understanding and controlling both contributions are important issues from technological and fundamental perspectives.
Transparent conducting oxides (TCOs) are a class of materials which are becoming ever more technologically relevant in a number of optoelectronic contexts [3][4][5], due to their combination of optical transparency and electrical conductivity. To date the choice of TCO for incorporation into device architectures has been dictated primarily by the bulk properties of the TCO. However, it is becoming increasingly apparent that rational design and optimization of novel devices must also consider the alignment of electronic energy levels at interfaces, e.g., in organic photovoltaics a high work function material is generally required for optimal performance.
Fermi level and band edge engineering in oxide materials, through doping and defect manipulation, is a well established process and computational modeling has been highly successful in the prediction of new doping strategies, in particular for SnO 2 [6][7][8][9][10][11][12][13][14][15][16]. The manipulation of the absolute electron energies (with respect to all other materials) is less well understood. The addition of dielectric layers and nanodots have been shown to improve performance and characteristics in several applications [17,18]; however, no consensus has emerged regarding the reasons for their success [19]. Several mechanisms have been put forward. It has been proposed that the layer can block metal induced gap states normally present at metal/semiconductor junctions [20,21]; alternatively, it has * a.walsh@bath.ac.uk been suggested that multipoles or fixed charges in the interface region result in a potential change across the interface, lowering the band offset [22,23].
Work functions and ionization potentials of oxide materials are extremely difficult to determine experimentally: Surface dipoles affect local vacuum levels, doping levels determine the Fermi energy, and the presence of defects alters both. The work function of SnO 2 is extremely surface sensitive [24] and has recently been shown to vary between 4.1 and 5.7 eV, depending on surface conditions and bulk doping [25]. Several theoretical schemes have been proposed for predicting the electronic energy level offsets of materials, from heuristic models, based on chemical electronegativities [26,27], to alignment of energy levels based on vacuum electrostatic potentials, determined from quantum mechanical calculations of two-dimensional (2D) slabs of the material [28][29][30][31], to explicit supercell simulation of materials interfaces [32,33].
In this work we investigate the fundamental factors which contribute to the binding energy of electrons in SnO 2 . Through the application of a recently developed multiscale modeling technique [34], as well as density functional theory (DFT), we are able to decouple the two quantities defined by Bardeen, i.e., the bulk electron binding energy and the surface contribution. This allows us to estimate the extent to which the surface controls and determines the work function. Furthermore, we resolve the surface effect into two additional categories: (i) a contribution from the crystal termination which affects energy levels in the bulk material and (ii) a contribution which is strictly confined to the surface region. By isolating these contributions we are able to place the SnO 2 electron energy levels on an absolute scale, allowing for their alignment with the energy levels of other technologically important materials, not including interface-specific effects. We explain the aforementioned reports of improved device performance through the inclusion of thin films and nanodots; moreover, we demonstrate the possibility of tuning energy levels through the inclusion of ultrathin films, similar to modifications using organic monolayers [35]. We consider a number of prototype situations of heteroepitaxial rutile capping layers on the SnO 2 (100) surface.

II. ELECTRON BINDING ENERGY IN THE BULK
To determine the absolute binding energy or ionization potential (IP) of an electron in bulk SnO 2 , excluding surface-specific effects, we employ a hybrid quantum mechanical/molecular mechanical (QM/MM) embedded cluster approach (see Fig. 1) [34]. The central cluster is treated at a QM level of theory using the PBE0 hybrid functional [36,37] and a correlation-consistent polarized valence-only double-zeta Gaussian basis set [38,39]. The value was also calculated with the metahybrid BBK1 functional, yielding a value within 0.1 eV of the PBE0 value, demonstrating the robustness of the methodology with respect to functional choice. The QM cluster is embedded within an external potential, provided by a larger cluster treated at an MM level of theory and a surrounding layer of point charges, fitted to reproduce the Madelung potential of the infinite crystal, which represents the system remainder [40,41]. The MM model is designed to reproduce accurately the structural, elastic, and dielectric properties of bulk SnO 2 (see the appendix for details on the force field). At the interface between the QM and MM regions, specially tailored effective-core pseudopotentials (ECPs) are placed on cationic sites to prevent spillage of electronic density into the MM region, and eliminate surface or interface effects [42] (see the appendix for further details).
The IP of the bulk material is determined from the total energy difference between the system in the neutral and positive charge states, allowing all electronic degrees of freedom to relax within a specified cut-off radius, beyond which long range polarization effects are accounted for [40]. Using different QM region cluster sizes (from 17 to 89 ions), we determine the IP to be 8.04 eV. In order to equate this quantity to the first contribution to Bardeen's definition of the work function [1], the bulk binding energy of an electron, it is necessary to define a reference average electrostatic potential in the material. We define this value as zero; it is equivalent to the potential in the region denoted "frozen potentials" in Fig. 1.

III. ELECTRON BINDING ENERGY AT THE SURFACE
We now calculate the ionization potentials in the presence of (110) and (100) surfaces, relative to a reference vacuum level (plateau in the Hartree potential), with a slab representation of the material, repeating periodically in two dimensions and terminating to a vacuum in the third. Here the IP is equivalent to the work function (the Fermi level is located at the top of the valence band); although, it should be noted that undoped SnO 2 is usually n type due to oxygen substoichiometry and hence the Fermi level will be found close to the conduction band.
Slab structures were created from bulk SnO 2 with cell parameters and ion positions relaxed (energy difference < 0.001 eV) using the PBEsol functional [43], projector augmented pseudopotentials [44], and a cutoff energy of 500 eV, with k-point sampling defined as an evenly spaced grid in reciprocal space with a density scaled to the unit cell size to achieve uniform sampling with a target length cutoff of 10Å, as described by Moreno and Soler [45]. All slab calculations were performed using the VASP code [46]. The surfaces chosen have been studied previously [47,48] and are known to be the two most stable surfaces in rutile SnO 2 .
The IPs were determined by hybrid functional calculations using 25% screened exact exchange [49]. The Hartree potential profile was plotted using a freely available code [50] based on the MATPLOTLIB package [51].
Capped surfaces were generated by replacing the Sn atoms in the uppermost layer of the (100) surface with a series of 115320-2 isovalent metal atoms, which also form rutile oxides (M = Si, Ti, Pb). The capping layer was generated on both surfaces of the 2D slab, to ensure electrostatic symmetry and to avoid the formation of a macroscopic dipole. The slab and vacuum layer widths were increased for each system until the Hartree potential in vacuum was fully converged. The use of a monolayer capping oxide means that the surface layer is below the critical thickness for reconstruction or formation of dislocations.
The effect of the surface on the band energies can be separated into two contributions (Fig. 1): (i) surface electrostatics (D s ) and (ii) intrinsic band bending due to the presence of evanescent surface states or changes in ion coordination. The first contribution arises because the electron density at a surface penetrates into the vacuum, resulting in a reduction in electron density immediately below the surface. The excess of electrons in the vacuum and the deficit of electrons immediately below the surface results in a multipolar layer, causing a potential step across the interface, penetrating the bulk of the material. The second contribution arises because the coordination of atoms at the surface is different from the bulk, resulting in electronic states characteristic of the surface (often within the band gap of the material) and from shifts in the energy levels of atoms close to the surface; the intrinsic band bending effect is strictly a surface effect.
To estimate contribution of the surface to the ionization potential we apply the following procedure: (1) Calculate the energy gap between the O 1s eigenvalues ( b s ) and the valence band maximum (E b VBM ) in bulk SnO 2 , with no surface effects.
(2) Calculate the bulk IP, from a total energy difference in QM/MM, as described previously.
(3) Calculate the O 1s eigenvalues at the center of the slab ( s s ) and the vacuum Hartree potential (V ) for the slab configuration.
(4) Evaluate the valence band maximum of the slab, without the influence of surface states, by comparison to the bulk calculation in (i): where s is the core-level shift, which is the difference between core s electrons in the bulk and the slab ( b s − s s ). (5) The slab IP excluding the influence of intrinsic band bending is evaluated from (6) The slab IP including the influence of intrinsic band bending (IP surf ) is evaluated as the difference between V and the highest occupied eigenstate of the slab ( h ).
(7) Finally the surface multipolar shift is evaluated as The values of IP for both surfaces are given in Table I. The value of 8.76 eV for IP surf of the most stable (110) surface, which contains band bending and surface electrostatic effects, is within the experimental range of ∼ 7.9-8.9 eV [25].
The values presented in Table I demonstrate the extent to which the surface determines the overall ionization potential,  (100) surface can re-arrange to a greater degree than those bonded to undercoordinated Sn at the (110) surface. Therefore, at the (100) surface, O electron density can stabilize the lower coordination environment more than at the (110) surface, resulting in evanescent surface states with lower energy and reduced intrinsic band bending at the (100) surface. The electrostatic effects at both surfaces differ very little and the value of the potential in the crystal bulk tends towards similar values in the presence of both surfaces.

IV. SURFACE MODIFICATION WITH HETEROLAYERS
In light of the results for the bulk material and the pristine surfaces, we now investigate how the surface contribution may be harnessed to control the ionization potential of a material. Ionization potentials and electrostatic potential shifts in the presence of hetero-oxide capping layers are reported in Table I. There is a significant change of the electron energies relative to the clean slab. This finding explains the aforementioned effect of ultra-thin dielectrics in improving device performance. By realigning the contact energy levels, the band offset at the interface can be tuned and reduced. The results demonstrate how, despite the capping layers consisting of isovalent isostructural metal oxides, the effect on the electron energies of the slab can vary by almost 1 eV. The effects of a capping layer depend on both the ionic and electronic structures, affecting both the local and long-range band edge positions in the substrate.
Charge density profiles for the (100) and capped surfaces are plotted in Fig. 3. The extension of charge density beyond the surface into the vacuum results in a pronounced decrease in the density at the surface O sites. The charge density close to the crystal surface reconstructs in an attempt to smear out the net positive charge remaining in the slab. In the density profile of the (100) slab, the charge density below the surface shows a reconstructed shape. The different capping layers result in different arrangements of charge at the surface and, consequently, different electrostatic fluctuations, as reported in Table I. SiO 2 has the largest effect on the surface dipole shift (∼ 1.1 eV), owing to both ionic and electronic rearrangements. is very similar to the clean (100) slab; however, the greater electronegativity of Pb, compared to Sn, results in a pronounced trough in the electron density just below the surface. Ti(IV) has the same oxidation state as Si(IV); however, it is significantly less electronegative and also has a smaller size mismatch with Sn(IV); therefore, the charge density, and surface dipole, are less affected by the Ti layer; the changes for Ti and Pb are ∼ 0.3 eV.
The relative bulk band edges of the different systems considered are shown in Fig. 4, demonstrating the prospect of tuning the bulk energy levels of a material by the inclusion of a thin capping layer. The systems calculated here already suggest the application of such capping layers in organoelectronic applications, where high IPs are required for contacting to deep molecular levels [5]. Currently Sn-doped In 2 O 3 (ITO) is used as an electrical contact, due to its high work function. The application of an SiO 2 capping layer could be used to engineer the band energies of fluorine doped tin oxide, making it a sustainable replacement for ITO.

V. CONCLUSIONS
We have presented a methodology for estimating the surface contribution to the crystal binding energies of electrons in SnO 2 . Furthermore, the surface contribution is separated into effects that are localized in the surface region and effects which penetrate into the bulk of the material. The ability to control the various components of the IP allows for the engineering of band energies through surface modification, a possibility demonstrated by the effects of ultrathin oxide films on the energy levels of SnO 2 . The same procedure can be extended, e.g., to include quasiparticle corrections to the band structure. The design principles of applying ultrathin films for ionization potential tuning described can be extended to any semiconductor, facilitating rational design of materials for optoelectronic applications.  The authors acknowledge A. A. Sokol and A. L. Shluger for constructive discussions. We also acknowledge support from the EPSRC (Grants No. EP/J017361/1 and No. EP/I01330X/1) and the Royal Society. The work benefited from the University of Bath's High Performance Computing Facility, and access to the HECToR supercomputer through membership of the UK's HPC Materials Chemistry Consortium, which is funded by EP-SRC (Grant No. EP/F067496). The authors acknowledge the use of the IRIDIS High Performance Computing Facility, and associated support services at the University of Southampton, in the completion of this work.

APPENDIX A: EMBEDDED CLUSTER-QM REGION
The GAMESS-UK [52] code was used to treat the QM region. A correlation-consistent polarized valence-only double-zeta Gaussian basis set was used for Sn and O ions, with 28 core electron ECPs on Sn [38,39]. The PBE0 hybrid functional [36,37] was used to model electron exchange and correlation.

APPENDIX B: EMBEDDED CLUSTER-MM REGION
We have fitted an interatomic potential model to treat the MM region in the embedded crystal, based on the Born model of ionic solids [53]. We simulate ion-ion interactions using a sum of four two-body terms and a three-body term. The first two-body term is a Coulomb sum: where U ij is the energy of interaction and r ij is the separation between ions i and j , and q i is the charge on ion i; the second is a Buckingham potential, of the form where the parameters A and ρ depend on the species of i and j . The third is a Lennard-Jones potential where B and C depend on the species of i and j . The fourth is a Morse potential of the form where D e , a, and r 0 depend on the species of i and j . The three-body term is a Bcoscross-type potential of the form where the parameters k, r 0 , m, and n depend on the species of ions i, j , and k. The polarizability of the ions is taken into account using the shell model of Dick and Overhauser [54], where each ion is separated into a core and shell, with the massless shell (charge Y ) connected to the core by a spring. The total charge of the core shell equals the formal charge of the ion. The energy is given by where K is the spring constant and r c−s is the distance between the core and shell. The parameters used are given in Table II.

APPENDIX C: EMBEDDED CLUSTER-INTERFACE REGION
To treat the interface between the QM and MM regions, a specially designed local ECP was placed on Sn sites located within a range of 5Å from the edge of the QM region. The ECP U p (r) has the form r 2 U p (r) = A 1 r exp(−Z 1 r 2 ) + A 2 r 2 exp(−Z 2 r 2 ) where the parameters A i and Z i were fitted in order to minimize the gradients on the ions in the QM and interface region, and the spread of deep core levels in the energy spectrum. The parameters are given in Table III.