UFDC Home  Search all Groups  UF Institutional Repository  UF Institutional Repository  UF Theses & Dissertations  Vendor Digitized Files   Help 
Material Information
Subjects
Notes
Record Information

Full Text 
A TIMEDEPENDENT MOLECULAR ORBITAL AP PROACH TO IONSOLID SURFACE COLLISIONS By ERIC QUINN FENG A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 1991 To Susan and Michael and to My Mother's Memory ACKNOWLEDGMENTS I am greatly indebted to Professor David A. Micha for his guidance during my research and Ph.D. thesis work. His clear understanding of physical problems, his enthusiasm in facing challenges and his spirit of devotion to the work always influenced me. His financial support enabled me to complete this work. I thank the Quantum Theory Project and the Department of Physics of the University of Florida for providing me excellent facilities in the course of this work. I also thank all QTP members who have constantly helped me in the past few years. I thank my colleagues and friends, Keith Runge and Robert Asher who have offered me many valuable suggestions. I also thank my friend Bob Meier for his help on the English language. Finally, my special appreciation goes to my wife, Susan; her constant encouragement and support helped me finish this work. TABLE OF CONTENTS ACKNOWLEDGMENTS ............................. iii ABSTRACT..................................... vii CHAPTER 1 .................................... 1 INTRODUCTION ................................ 1 1.1 Collisional Charge Transfer at Surfaces ................ 1 1.2 A Survey of Theoretical Methods ... ................ 4 1.2.1 Binary Collision Theory of Charge Transfer ........... 5 1.2.2 ManyElectron Treatments .................... 6 1.3 Existing Problems ............................ 8 1.3.1 Electronic States and Electronic Couplings ........... 8 1.3.2 Treatments of Extended Systems ......... ...... 11 1.3.3 Effect of Nuclear Motion .................... 12 1.3.4 Atomic Orbital Polarization ................... 13 1.4 Outline of the Chapters ........................ 14 CHAPTER 2 .................................... 19 THE TIMEDEPENDENT HARTREEFOCK (TDHF) APPROACH 19 2.1 The AtomSurface System ...................... 20 2.2 TDHF Equations of Density Matrix ................. 21 2.3 TDMOs as Linear Combinations of Travelling Atomic Orbitals 25 CHAPTER 3 ........ ................ ........... 28 AVERAGE ELECTRONIC POPULATIONS, ELECTRIC MULTIPOLES AND ORBITAL POLARIZATION ................... 28 3.1 Coordinate Frames .......................... 31 3.2 Electric Multipoles ........................... 33 3.3 Alignment and Orientation Parameters ............... 35 3.4 Multipoles and Alignment and Orientation Parameters in a Subsystem ............................... 38 CHAPTER 4 ................................... 41 LINEARIZATION OF TDHF EQUATIONS ................ 41 4.1 The Linearization Procedure .................... 42 4.2 The Case without ElectronElectron Interaction .......... 48 CHAPTER 5 ................................. 51 PARTITION OF EXTENDED SYSTEMS ................. 51 5.1 The Partition Procedure ........................ 52 5.2 The Approximation in the Secondary Region ........... 55 CHAPTER 6 .. ................................ 58 ELECTRONIC BASIS FUNCTIONS AND MATRIX ELEMENTS 58 6.1 Generalized Wannier Functions ................... 60 6.1.1 Definition of Generalized Wannier Functions (GWFs) .... 60 6.1.2 Generalized Wannier Functions as Linear Combinations of Gaussians .............................. 63 6.1.3 Determination of Generalized Wannier Functions ....... 65 6.2 Generalized Wannier Functions for a Jellium Surface ....... 67 6.2.1 Results .......................... ... 72 6.3 Atomic Basis Functions ........................ 86 6.4 Overlap Matrix Elements ....................... 87 6.5 Hamiltonian and its Matrix Elements ................ 88 6.5.1 Pseudopotential for the Atomic Core .............. 88 6.5.2 A Corretction Term to the Hamiltonian ............ 90 6.5.3 Hamiltonian Matrix Elements ................. 92 CHAPTER 7 ................................... 94 INTEGRATION OF LINEARIZED TDHF EQUATIONS ........ 94 7.1 The Algorithm for Numerical Integration of TDHF Equations 94 7.2 Computation Program ......................... 95 7.3 Stability and Convergence of the Numerical Integration ..... 97 7.3.1 Tolerances ............................. 97 7.3.2 Initial and Final Distances .................... 98 CHAPTER 8 ................... APPLICATIONS ......... ....... .............. 100 ...... ......... 100 8.1 NaW(110) Model System ........... 8.1.1 Hamiltonian and Basis Functions . 8.1.2 Electronic Couplings ........... 8.2 AtomSurface Interaction Potentials and the 8.2.1 AtomSurface Interaction Potentials .. 8.2.2 Trajectory ................. 8.3 Charge Transfer ................ 8.4 Results ...................... 8.4.1 Evolution of Electronic Populations . 8.4.2 Electronic Populations after Collisions . CHAPTER 9 ....................... DISCUSSION AND CONCLUSIONS ....... APPENDIX A ...................... CALCULATION OF COEFFICIENT MATRIX B APPENDIX B ...................... CALCULATION OF TOTAL ENERGY OF APPENDIX C .................. ABAB .... ............ 101 ............ 101 ............ 102 Trajectory ...... 103 ............ 103 ............ 105 .......... 108 ............111 ........... 111 . ........... 113 ............125 ............ 125 . . 138 ............ 138 . . 138 ............140 ND .......... 140 ............ 142 CALCULATION OF ELECTRONIC PROPERTIES OF A FINITE SLAB ............................ ........ BIBLIOGRAPHY .............................. BIOGRAPHICAL SKETCH ........................ Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy A TIMEDEPENDENT MOLECULAR ORBITAL APPROACH TO IONSOLID SURFACE COLLISIONS By Eric Quinn Feng May 1991 Chairman: David A. Micha Major Department: Physics A timedependent molecular orbital method has been developed to study charge transfer in collisions of ions with metal surfaces at energies between I au and 170 au. A set of localized basis functions, consisting of generalized Wannier functions for the surface and s and p atomic functions for the ion, is used to separate the system into primary and secondary regions. An effective Hamiltonian and timedependent equations for the electron density matrix are obtained in the primary region, where most charge transfer occurs. The equations for the electron density matrix are solved with a linearization scheme. The method is suitable to study atomic orbital polarization for collisions of ions and surfaces. A model calculation for Na'+W(110) collisions with a prescribed trajectory is presented. The interaction potentials between the W(110) surface and Na 3s and 3p orbitals are calculated from a Na pseudopotential and a step potential for W(110). Results show that the yield of neutralized atoms oscillates with the collision energy as the result of the nearresonance charge transfer mechanism. The timedependence of the density matrix provides insight on the dynamics of electron transfer along the atomic trajectory. CHAPTER 1 INTRODUCTION 1.1 Collisional Charge Transfer at Surfaces Surfaces and various physical processes occurring near surfaces are receiving increased attention in present scientific research, due mainly to rapidly growing applications in modem analytical techniques. Researchers in fields of surface sci ence, chemical physics, solid state physics, atomic physics, electrical engineering, and even medical science are carrying on extensive, and sometimes overlapping, studies and analyses involving surfaces. One of the most important ways to un derstand the nature of surfaces is to study various electronic processes occurring in collisions of ions with them. For example, one can learn from scattering about the topology and electronic structures of surfaces, amount of impurities and ad sorbates present, surface electron and lattice relaxations, and interactions between surfaces and atoms and molecules. Particular attention has been paid to nona diabatic processes accompanying charge transfer between scattered particles and surfaces at low collision energies ranging from a few electron volts to a few thousand electron volts, which often involve strong coupling and energy transfer between electronic, translational, vibrational, and rotational degrees of freedom of projectiles and surfaces. 2 Charge transfer in ionsurface collisions may lead to neutralization or ioniza tion and can be expressed as S+ A+  E+ A (1.1) and E + A + A'F (1.2) respectively, where E is a surface such as W(110) or Ni(ll11), and A is an alkali atom such as Li, Na or K. Similar processes occur in beamfoil experiments in which charge exchange takes place when ions or atoms pass through thin foils, instead of being reflected from surfaces. The two main mechanisms through which charges are exchanged are the near resonance process and the Auger process [Hagstrum, 78]. In a nearresonance neutralization, an electron in the conduction band of the solid tunnels to an unoccupied level of the ion which lies near or below the Fermi level of the solid. In a nearresonance ionization, an electron from a filled excited state of a neutral atom transits to an empty level in the conduction band or above the Fermi level of the solid. The two energy levels involved in these processes are close, which characterizes the nearresonance processes. The Auger process involves two or more electrons. In twoelectron Auger neutralization, one electron jumps from the metal conduction band to the ion ground state below the band and, to conserve energy, the second electron is excited from another level in the band and may leave the surface and be detected by an Auger spectroscopy. For systems 3 where an ion affinity level is in resonance with the conduction band of the metal, such as in Na++W(110), a nearresonance process is more likely; for systems where an ion ground state is below the conduction band, such as H++Ni(l11), the Auger process is usually the dominant charge transfer mechanism. Complicated combinations of the two processes are also possible in certain circumstances. The rise of interest in these problems has promoted an extensive effort to develop ultrahigh vacuum (UHV) and modern spectroscopy techniques to observe and detect various physical properties in much greater detail in the last two decades than before [Hagstrum, 77; Hagstrum, 78]. In these spectroscopy experiments, one measures the intensity of scattering products (atoms, photons or Auger electrons) and their energy and angular momentum distributions which usually depend on the energy and angle of incident beams, the surface properties and, sometimes, on temperature and external fields. The ion scattering spectroscopy (ISS) [Smith, 71; Hulpke, 75; Taglauer and Hailand, 76; Overbosch et al., 80] measures the scattered ion intensity as a function of incident energy and reflected angle, and can be used to determine the topology of surfaces. Oscillation in the intensity of neutralization of ions in the low energy range has been observed for a number of ionsurface combinations, and is ascribed to nearresonance charge transfer [Erickson and Smith, 75; Rusch and Erickson, 77]. In ion neutralization spectroscopy (INS) experiments, one can gain knowledge of the band structure by measuring the kinetic energy of Auger electrons [Hagstrum, 75]. Time of flight (TOF) experiments detect survival probabilities of states of scattered atoms 4 and their energy distributions [Kasi et al., 88]. So far, most experiments have been done on noble and transition metals, simply because their surfaces are relatively easy to clean and to prepare. However, due to the strong interest in electronics, studies on semiconductor and other material surfaces have been carried out [Richard and Eschenbacher, 84]. In atomic collision studies, people have developed new instruments which have a time resolution of picoseconds to monitor the time evolutions of inter mediate atomic states. It has been possible to measure the angular distribution of scattering product orbitals [Hale et al., 84; Hertel et al, 85; Andersen, 87; Campbell and Hertel, 87]. These new techniques provide valuable information about the collision dynamics, and it is expected that they will be soon used in studies of ionsurface collisions. 1.2 A Survey of Theoretical Methods Advances in spectroscopy techniques in the last two decades have allowed ex tensive experimental studies of ionsurface collisions. However, the development of the theory lags behind that of experiment in this field, and there are far and away more experimental data than can be understood. This is because the systems we are dealing with are very challenging ones in the sense that they lack transla tion symmetry. Moreover, they involve dynamical manyelectron processes with a strong and timedependent coupling between the electronic degrees of freedom and nuclear motions. The traditional methods which deal with systems of periodic 5 lattices or of a few electrons are either not suitable or have to be modified. It seems, however, that theoretical studies are accelerating and are on the verge of reproducing some of the experimental data. 1.2.1 Binary Collision Theory of Charge Transfer In this theory [Rusch and Erickson, 77; Kasi et al., 89] the scattering of an ion from a surface is described by a single binary collision or a sequence of elastic collisions between the ion and surface atoms. In this theory, the scattering angle 0 satisfies the relationship El [cosO + (' sin20)1/2]2 (1.3) Eo0 (1 +p)2/ (3) where mi and m2 are the masses of the incidental ion and surface atom, =m2/ml, and Eo and Ei are the kinetic energies of the ion before and after the scattering. The yield of the scattered ions for a given scattering angle, Y(Eo, 0), is Y(Eo, 0) oc a(Eo, 0) P(Eo, 0) (1.4) where a(Eo, 0) is the elastic differential cross section and P(Eo, 0) is the probability for an ion to remain ionized after the scattering. In classical calculations a is a monotonically decreasing function of the incident energy and P is given by P oc exp(a/v.) (1.5) where a is a constant and v1 is the ion velocity perpendicular to the surface. This theory, in which the yield is a smooth function of the incident ion energy, fits the energy dependence of 4He+ scattered from Cu and certain other systems. 6 However, it cannot explain the oscillatory charge transfer behavior observed in other systems which is believed to be associated with nonadiabatic charge transfer mechanisms. 1.2.2 ManyElectron Treatments A variety of theoretical methods [Tully, 77; Brako and Newns, 81; Holloway and Gadzuk, 85; Hood et al., 85; Lee and George, 85] have been proposed to ac count for the manybody processes accompanying charge exchange in ionsurface collisions. The most widely used one is the approach based on the Anderson Newns Hamiltonian [Anderson, 61] to describe the ionsurface scattering [Blandin et al., 76; Norskov and Lundqvist, 79; Brako and Newns, 81; Lang, 83]. In this approach, neglecting spin and using second quantization notations, the Hamilton ian of the system is written as H = EaC +Ca Ci kC+k + [VakC gCk + VakC+Ca] (1.6) k k where C+a, Ca are the creation and destruction operators corresponding to atomic state Oa, Cj!, Ck are the creation and destruction operators corresponding to solid state ?k, ca and ek are the energies of ba and kk respectively, and Vak = (alVIk), where V is the perturbation due to coupling of the atom to the metal. The last two terms describe the electron hopping between 1ia and bk. It should be noted that the energy of the atomic level is time dependent through the ion's trajectory RA(t), that is, fa(t) = Haa [RA(t)]. A set of equations of motion of electronic 7 degrees of freedom can be obtained from the Heisenberg equations (with h = 1) i = [H, Ca] = ca(t)Ca + Vak(t)Ck (1.7) k i = [H, Ck] = Ck Ck + Vak(t)Ca (1.8) Eliminating Ck, it yields a differentialintegral equation for Ca which can be solved by numerical procedures. The number of electrons on the atom after the scattering is given by na(oo) = ('a IC+(oo)Ca(oo)Ila) (1.9) Most studies using the timedependent AndersonNewns Hamiltonian have neglected the intraatomic Coulomb repulsion to simplify the calculations. For most alkalimetal collisions, which produce few negative ions, this approximation seems to be justified. But for scatterings like HW(110), there is a significant fraction of H products after scattering, and the intraatomic Coulomb repulsion plays an important role. Instead of using the above Hamiltonian one should use H =1 EC aCa. + 1 fCkoCo or ko (1.10) + E [VakCoCko + VaC+Ca.] + U(t)naTa.n ko where a is a spin index, naT = CaCa, nal = C1CaJ, and U(t) is the effective intraatomic Coulomb repulsion which depends on the distance between the ion and the surface. Using the HartreeFock approximation, one obtains a set of 8 differentialintegral equations [Yoshimori et al., 84; Sulston et al., 88] which can be solved numerically. However, the qualitative behavior of the results are found to be not much different. There have been several attempts to go beyond the HartreeFock approximation [Sebastian, 83; 85; Kasai and Okiji, 87], by which some preliminary results were obtained. 1.3 Existing Problems Several problems inherent in the theoretical treatment of ionsurface scattering will be addressed here. We view them as some key points in furthering the theoretical study of ionsurface collisions. 1.3.1 Electronic States and Electronic Couplings Ionsurface systems are many electron systems which lack translation symme try. This brings both difficulty and challenge to theoretical studies of the problem. One way to deal with surfaces is to assume that the electronic wave functions in surfaces are the same as those of an infinite periodic system and to treat the colliding ion as a moving defect and the atomic states as localized states. In this way, one can relatively easily calculate surface properties and the scattering yield [Brako and Newns, 81; Lang, 83; Shindo and Kawai, 86]. However, this approximation is not very accurate, and usually can only be used as a starting point. At surfaces, as the results of broken symmetry and the strong perturbation from the colliding ion, the electronic bands are distorted and the electron behavior is different from that of bulks. One of the important features of the metal and 9 semiconductor electronic structure is the existence of the localized states. While s and pelectrons are generally described by continuous states, delectrons are more localized. At surfaces, localized states can also be created by impurities, adsorbates or chemisorbed layers. The presence of localized states greatly affects the local electronic environment; for example, an absorbed layer of alkali on metal surfaces lowers the work function [Gomer, 75; Medvedev et al., 70]. Correctly calculating or estimating the electronic couplings between the ion and surface is obviously important in theoretical calculation. But current theoreti cal methods do not provide a simple and yet accurate way to estimate the electron hopping matrix elements. People have to rely on semiempirical calculation. For example, in the timedependent AndersonNewns Hamiltonian, Vk measures the coupling between the atomic states and surface states and is related to the lifetime broadening of an atomic state, A(t). Assuming that k and t dependence of Va are separated, Vak [RA( ] = Vku(t) (1.11) where Vk depends only on k and u(t) only on t, one has A(t) = A(c)lu()2 (1.12) where A(e) is defined by A(e) = jVkl2b( k) (1.13) k which can be approximately evaluated by methods developed in chemisorption theory. But a common practice is to use parameterization of A(t). For instance, 10 one can use a simple exponential form which is independent of energy A(t) = Aoexp(ZA) (1.14) and determine parameters Ao and 7 from experimental data [Brako and Newns, 81; Lang, 83; Hood et al., 85] or from density functional calculation [Lang and Norskov, 83]. Because of the electronic coupling from the surface electrons, the atomic energy ea is distance dependent. Again, a simple function form is often used (with the electron charge e=l au), ea(za) = Ca(oo) (1.15) 4ZA where ZA is the distance between the ion core and the surface. This is based on the approximation that the change of ea is due to a classical image charge. Although the parameterization approach is widely used, it ignores the fact that the electronic coupling is strongly dependent on the rearrangement of charges during collisions, and its validity at short distances is questionable. Ideally, the treatment of the ionsurface systems should contain a selfconsistent electronic structure calculation which accounts for localized states. A selfconsistent lo calized orbital method has been developed to deal with transition metal surfaces and chemisorption on metal surfaces [Smith and Gay, 75; Smith et al., 80]. The density functional method [Hohenberg and Kohn, 64] has also been applied to ionsurface interactions [Lang and Williams, 76]. 11 1.3.2 The Treatment of Extended Systems From the theoretical point of view, it is easier to deal with either a system of a few particles, for example, small molecules, or a system with a periodic structure, such as single crystals. Ionsurface systems, on the other hand, are manybody systems which lack translation symmetry. To combine a full self consistent calculation of electronic interaction with a dynamic description of ion surface collisions is a difficult, if not impractical, task. In developing a better and more physical description of ionsurface systems a very important fact should be noticed, namely that during the collisions, ions interact mainly with local regions of the surface, while the remainder of the solid is relatively unperturbed. This fact provides a hint that it is possible to properly handle electronic motions by concentrating on these local regions. Olson and Garrison [Olson and Garrison, 85] have used a cluster in place of a surface; this enables them to employ the molecular orbital method developed in molecular scattering. Another method to deal with surfaces is the embedding technique [Grimley and Mola, 74; Kirtman and de Melo, 81; Feibelman, 85] used in chemisorption studies, which uses a higher level treatment of electronic interactions within a molecular complex (a small region of the surface plus adsorbates) and a simple description outside of the complex. de Melo et al. [de Melo et al., 87] have developed a selfconsistent method using a density matrix and applied it to the AndersonNewns Hamiltonian for a one dimensional system. McDowell [McDowell, 85] uses an embedding technique to obtain generalized 12 Langevin type equations for the spin orbitals in a primary zone which couples with the secondary zone through driving terms. These treatments show that, by concentrating on a small number of orbitals, it is possible to achieve self consistency. Another common feature in these treatments is the use of localized spin orbitals, which allows one to naturally deal with localized phenomena [Feng et al., 91]. 1.3.3 Effect of Nuclear Motion While the electronic motions are treated quantum mechanically, it is reason able to assume that the nuclear motion evolves according to the classical mechan ics for the collision energy in the hyperthermal range. In this classicalquantal approach the nuclear motion, RA(t), can be treated at several levels. In the simple classical treatment, the nuclear trajectory is fully determined by a classical ion surface potential. In an improved semiclassical treatment, the nuclear potential is coupled with the electronic motion. The classical trajectory is valid only in the high energy range [Tully, 76] as demonstrated in gasphase collisions. The reason is that at lower energies, trajectories are sensitive to the detail of chemical inter actions among the atoms and depend on the electronic states. It has been found that in Hp collisions, the classical trajectory would give rise to a significant error when the collisional energy is lower than 100 ev [Runge et al., 90]. Looking from another angle, the nuclear motion can be treated either by the socalled trajectory approximation or by a multichannel procedure. In the trajectory treatment, the nuclear trajectory is uniquely defined by a single fixed 13 potential, which can either be a pure classical one or an effective potential which contains the coupling with the electronic degrees of motion. One example is the straight line trajectory resulted from a hard wall potential. Runge et al. [Runge et al., 90] determine the ion trajectory by an effective force which is dependent on the gradient of the electronic density matrix. It has instead been suggested using a trajectory approximation, constructing multidimensional potential energy hypersurfaces describing various atomic configurations and allowing trajectories to hop back and forth between hypersurfaces [Tully, 77; Holloway and Gadzuk, 85; Newns, 85]. A useful procedure in dealing with molecule collisions is the eikonal method developed by Micha [Micha, 83]. In this method the nuclear variables are coupled with the electronic ones and both must be found selfconsistently. An application of the eikonal method to ionsurface collision is by Olson and Garrison [Olson and Garrison, 85] who model the surface by a small cluster. 1.3.4 Atomic Orbital Polarization In atomic scattering, the scattered atoms can be in different electronic states and their distribution of electronic angular momentum can be anisotropic. This causes orbital alignment and orientation [Hippler, 85]. Experimentally, this requires a partial or full determination of the atomic states after the collision in contrast to the conventional study which measures the differential cross section. Orbital polarization has been the subject of much theoretical work [Fano and Macek, 73; Andersen and Nielsen, 87; Nielsen and Andersen, 87], in which the 14 density matrix [Blum, 81] is extensively used. Although the use of the density matrix in scattering theory is not new, the study of orbital polarization provides a much more severe test of a theory than does the study of cross section, this is because not only the diagonal elements, but also the offdiagonal elements of the density matrix are required to determine the alignment and orientation parameters. Of course, the study of orbital polarization will lead to a deeper insight into the interaction and mechanisms involved in the collisions. Until now, no detailed work on orbital polarization in ionsurface collision had been reported either experimentally or theoretically. It is just a matter of time before the experimental techniques and theory are developed in this field. 1.4 Outline of the Chapters This dissertation presents a theoretical study of charge transfer in ionsurface collisions based on the timedependent molecular orbital approach. A partition ing technique is used to divide a ionsurface system into a primary region which contains a few centers on the surface and the scattering ion, and a secondary region containing the remainder of the surface. The charge transfer is calculated by solving the timedependent HartreeFock (TDHF) equation for the electron density matrix in the primary region which couples with the secondary region electronically. This approach permits determination of the atomic states during and after the scattering and calculation of the alignment and orientation parame ters. In the following chapters atomic units are used; therefore, electron charge e=l, electron mass me=l and F=1. 15 In Chapter 2, we present the basic framework of our approach to the ion surface collisions. TDHF equations are used to describe the time evolution of the molecular orbitals. The formalism is constructed in terms of an electronic density matrix because it can simplify computation and is convenient in the calculation of alignment and orientation parameters. In the later part of the chapter, the TDHF equation is rewritten in the basis of the travelling atomic orbitals (TOA). In Chapter 3, we define the parameters for average electronic population, electronic multipoles and polarization of atomic orbitals by introducing tensor operators and irreducible operators. The orientation and alignment parameters can be calculated from the density matrix. The orbital polarization parameters for a subsystem are also expressed by the density matrix. In Chapter 4, we present a local time linearization procedure for solving TDHF equations. The method is based on the assumption that, in a very short time interval, the solution of the TDHF equation is a linear perturbation caused by the motion of the atomic core, added to correct the evolution of the system first calculated as if the nuclei were fixed. Using this procedure, the linearized equations for the perturbations, which contain a timedependent driving term due to the nuclear motions, are constructed and solved in conjunction with the equations for fixed nuclear positions. The procedure is presented for a general case. In the later part of the Chapter, we apply this procedure to a simple case in which the electronelectron interaction is neglected and obtain analytical solutions for the linearized TDHF equations. 16 In Chapter 5, a partitioning procedure for ionsurface collisions is presented. The system is divided into a primary region which contains the ion and a few centers on the surface and a secondary region which is the remainder of the surface. The density matrix is approximated in the secondary region by an unperturbed one, which is justified when the collisional energy of the ion is not very low. The effective equations in the primary region are coupled with the secondary region through a driving term. To fit the asymptotic behavior of the system after the partitioning, a correction term is added to the partitioned Hamiltonian. Combining the linearization procedure described in Chapter 4, we obtain the effective equations in the primary region and their formal solutions. Again, we further illustrate this partition method by applying it to a special case when the electron interaction is neglected. In Chapter 6, we construct a set of localized basis functions to be used in the timedependent molecular orbital calculation. For the surface part, we introduce a set of generalized Wannier functions (GWFs) which are orthonormal and localized at centers of the surface. A variational procedure is employed to generate these functions. As an example, we calculate the generalized Wannier functions for a three dimensional jellium solid. The basis functions associated with the ion core are assumed to be atomic orbitals and are obtained from a pseudopotential for the ion core. Our basis functions are timedependent and in the form of linear combinations of Gaussians. The overlap matrix and Hamiltonian matrix are constructed in this basis set. 17 In Chapter 7, we discuss the computational aspects. At first, we present our algorithm for numerical integration of the timedependent HartreeFock equations as well as the flowchart of the computational program. Then the stability and the convergence of the numerical integration are discussed at length. Effects of some computational parameters are investigated, including step size, initial and final positions of the ion, and tolerances, among others. In Chapter 8, we present an application of this approach and its results to the neutralization in a normal collision between a Sodium ion and a Tungsten (110) surface. An atomsurface interaction potential is constructed to determine a prescribed trajectory. By using methods described in previous chapters, we obtain the time evolution of the density matrix from which the electronic population and orbital polarization parameters are calculated. We also study the effects of the collision energy on the final yield of the neutral atoms and their polarization after the collision, and compare the results with experimental data. We also briefly discuss the application of the approach to other systems, and to collisions at other scattering angles. In Chapter 9, we discuss the features of our timedependent molecular method and its applications to charge transfer in atomsurface collisions in hyperthermal energy range. Analysis and conclusions are related to the physical features and comparison with other theoretical methods in the field. We also discuss the approximations used in the present study. Finally, we offer suggestions for future research. 18 In Appendix A, we calculate matrix Bx which is the expansion coefficient matrix of the Wannier functions for free electrons in terms of Gaussians. In Appendix B, we calculate the total band energy of a jellium slab using the generalized Wannier functions. In Appendix C, we show how the Fermi energy is calculated for a finite jellium slab. CHAPTER 2 THE TIMEDEPENDENT HARTREEFOCK (TDHF) APPROACH The timedependent HartreeFock approach is suitable for studying the many body collision dynamics of our problem. Using a oneelectron effective field it can describe the evolution of collisional systems and calculate dynamic parameters to compare with experimental data. One of the advantages of the TDHF method is that the selfconsistency between the effective field and the orbital coefficients or density matrix is achieved automatically by solving timedependent differen tial equations without having to perform the selfconsistent iteration procedure required in the timeindependent case. TDHF has been extensively applied to collision problems in nuclear physics [Ring and Schuck, 80; Negele, 82; Davies et al., 82], and recently applications in atomic molecular collisions have been reported [Kulander et al., 82; Devi and Garcia, 83; Gazdy and Micha, 86; 87]. Micha and Gazdy [Micha and Gazdy, 87] have proposed a variational procedure which improves the accuracy of transition amplitude calculations using TDHF trial functions. Yoshimori et al. [Yoshimori et al., 84] have discussed using the TDHF method in charge transfer in ionsurface collisions. In this Chapter, we apply the TDHF method to charge transfer in ionsurface collisions. In Sect. 2.1, by introducing timedependent molecular orbitals, 20 we derive the TDHF equations for the electron density matrix. The reason we use the density matrix [McWeeny, 60; Cohen and Frishberg, 76] is that it simplifies calculation for a large system. Another advantage is that many interested properties of collisions can be defined and readily calculated from it. One of the key problems in applying the TDHF approach to ionsurface studies is choosing the basis set. With the great number of electrons in the systems, choosing a complete basis set would make any practical calculation out of the question. One has to truncate the complete basis set in a proper way so that the introduced error is minimal. However, we will leave this problem to Chapter 5. In this Chapter, we first describe the ionsurface system. In Sect. 2.2 we derive the basic formalism without specifying the particular form of the basis set except for the condition that each basis function is associated with a certain center of the system. In Sect. 2.3 the TDHF equation is rewritten in the basis of travelling atomic orbitals (TAO). 2.1 The AtomSurface System Let's consider a system of a scattering ion and a semiinfinite solid. We will assume the solid has a conduction band of continuous oneelectron states and some localized states (This is in contrast to some approaches which are limited to continuous states only). For the colliding atom we assume it has a core and several valence levels. Because the deeper levels of the solid and the core states of the ion are tightly bound and usually do not participate in charge transfer in the energy range we are interested in, they are not considered. For the same 21 reason, highly excited atomic states and those solid states highly above the Fermi level are not considered either. We also assume that the atomic levels are in the range of or close to the conduction band and that the nearresonance tunnelling of electrons is the dominant charge transfer process. The thermal motion of the nuclei in the solid can be neglected at the collision energies of interest, so that the only moving nucleus is that of the scattering ion. For simplicity, we assume that the collision energy is high enough that the trajectory approximation can be used. In this approximation, the motion of the ion nucleus can be described by a known function RA(t) determined by an effective interaction potential between the ion core and the surface. In Chapter 9 we will discuss the validity of this approximation and the possible improvements which are necessary at low collision energies. 2.2 TDHF Equations of Density Matrix Let T (t) be the electronic wavefunction which satisfies the timedependent Shrodinger equation i T(t) = H(t) W(t) (2.1) where H(t) is the total Hamiltonian of the system. Using the timedependent HartreeFock approximation, the TDHF wavefucntions have the form of I(t) = N"l1/2det[Oi(xj, t)] (2.2) 22 where Net is the number of electrons and Oi(x, t) is the ith spin orbital for electron variables z = (r', ). We write bi(x, t) as =i(x, t)(= (r, t)7i(() (2.3) where is the corresponding spin function. The electron density operator A is given by E= Z )(01 (2.4) iEocc where the summation is over the occupied orbitals, which is assumed to satisfy the TDHF equation [Dirac, 30] AnJp = p'rr rFr (2.5) where P is the Fock operator which can be written as [Pople and Beveridge, 70] P(z1, t) = HI(F, t) + Gr(x, t) (2.6) where 1 ftI(i, t) = vi +VA(l, t) + VM(i, t) (2.7) is the core Hamiltonian operator with VA the potential from the atomic core and VM the potential from atomic cores in the solid; and 1(2.8) G6(xI, t) = fdx20L7(x2,t) r[1 P i2]M(x2,t) (2.8) is the electron Coulomb potential operator, where P12 is the permutation operator exchanging electrons only between spinorbitals of spin 7. Introducing a set of 23 timedependent basis functions {((r, t)}, p=1, 2, .. N, the molecular orbital is written as a linear combination of the basis functions so that 7(F,,t) = (t)(rf,t) i= 1, 2, (2.9) where cU (t) is a timedependent coefficient. The density matrix P is defined as P = I)p (lI (2.10) where I) and (I[ are row and column matrices of C1 orbitals, respectively. The matrix element of P is [Pople and Beveridge, 70; Szabo and Ostlund, 82] ,(t)= ( C(t)cfv(t) (2.11) iEocc The Fock matrix F7 is defined as F" = (\FI ) (2.12) Inserting equation (2.6) into the above definition we have F' = H + G' P, P') (2.13) where H= K +VA+ VM (2.14) is the core Hamiltonian matrix, G7 the HartreeFock electronelectron interaction matrix, VA the atomic potential matrix and VM the surface potential matrix. The matrix elements of F' are F,(t) = H,,(t) + P ~ (t)(pl a) + P (t)(pvIXa} (2.15) Ao 24 where 7' is the spin opposite to 7 and with (pVII\A) = (uV\jA) (p\IjAv) (yIvAa)= d d'r2( t)(rt)=) 1 (Fi, t)o(rl, t) r12 (2.16) (2.17) To derive the TDHF equation for the density matrix, let's define a matrix a Substitute (2.10) into (2.5), the left hand side is .0 .a t FP a= ie g)hi P'( + ai + iil and the right hand side is (2.18) (2.19) (2.20) Pp / p = PF)P  ( Io)P  (vI Multiply this equation by (5 from the left and by I) from the right, we have inlP'S + iSP'S + iSP'tIf = FTPfS SPTrF (2.21) where S = (J) is the overlap matrix. Multiply both sides of the above equation by S1 we obtain the timedependent HartreeFock equation for the density matrix i P = S1FP'p PF'rS' iSlP'7 ipt'rtS1 (2.22) The last two terms on the right hand side of Eq. (2.22) arise from the time dependence of the basis functions. 25 2.3 TDMOs as Linear Combinations of Travelling Atomic Orbitals When the molecular orbitals are expanded in a basis set of atomic functions, it sometimes introduces artificial couplings at large distance which originates in the dependence of the atomic orbitals on the position of the moving nuclei [Bates and McCarroll, 1958]. To avoid this one can choose the atomic basis functions in the form of travelling atomic orbitals (TAO). We choose (4 to be a travelling atomic orbital (TAO) associated with the mth center at R(t), (F, t) = x, [F m(t)] Tm(f, t) (2.23) where r is the position vector with respect to the origin of the reference frame, X, is the atomic orbital centered at Rm(t) and the translation factor Tm(r, t) is defined by T(r, t) = exp im, (M i Iv2 dt' (2.24) tin where me is the electron mass and m = dAm/dt is the local velocity of the mth center in a spacefixed system. Thus, a "v = H i&') + n Vn vl) (2.25) and (2.26) Using (2.22) we find S= ime S' + (xITT. n V.)l,) (2.27) where (2.28) = d rX(rT.'(r')Tn(r) [( a. Vn)x(rv ] Noticing that Vx,) = VnIXv), the kinetic energy matrix in this basis is ((IKI v) = ((4 V 2 = m,2 s + i(plT* T Vn)Ix.) + (x.IKIx.) (2.29) = illv + KIj where Eq. (2.26) has been used and k = (xIx,) = d3r T, (F, t)T,( t) x(r, t)Kx,(r, t) (2.30) Substituting Eq. (2.30) into Eq. (2.22) it becomes iP7 = S1f'rp P'FtS1 (2.31) with a modified Focklike matrix F7 = I + G (2.32) where = + VA + VM (2.33) Equation (2.31) is our basic equation. It appears in a simpler form in the TAO basis, but the price paid is that the matrix F7 is no longer Hermitian. This equation can be solved numerically if the nuclear trajectory is known, which may or may not include the effect of the coupling between the electron 27 and nuclear degrees of freedom. In the eikonal treatment, this equation should be solved with the equations of motion for nuclei simultaneously. Runge et al. [Runge et al., 90] have applied this method to charge transfer in the collision between a hydrogen atom and a proton. So far, we have put no restrictions on the choice of basis functions except that the functions are localized at certain centers of the system. In Chapter 6 we will discuss in detail the choice of basis functions, which is a key part to our approach for handling extended systems. CHAPTER 3 AVERAGE ELECTRONIC POPULATIONS, ELECTRIC MULTIPOLES AND ORBITAL POLARIZATION One of the parameters which characterizes scattering phenomena is the cross section, which has the dimension of area and, roughly speaking, measures the size of electron clouds. The total cross section and the differential cross section for state to state transitions can be measured in experiments and provide fundamental information about excitation mechanisms. Thus extensive efforts have been made to develop theories and methods to calculate the cross section. On the other hand, during collisions with targets, absorption of energy will excite atoms to various states which can be anisotropic. Electron clouds not only change their sizes but also change their shapes and rotations, which is termed orbital polarization. The shape change and rotation of the electron cloud are characterized by alignment and orientation. Recent advances in experimental techniques have made it possible to completely determine the states of scattered atoms. For example, in the HtH experiment [Hippler et al., 86], H experiences a 2s to 2p excitation and the distribution of electrons on m=l and m=l states varies with the collision energy. The orbital polarization reflects some of the specific and subtle aspects of the collision processes which can not be fully unveiled in the cross section study. Obviously, knowledge of orbital polarization properties can offer a deeper insight 29 into the collision mechanism and a more detailed understanding of the dynamics of collisions. Such a study will also provide a sensitive test to the scattering theories and models. The density matrix method has been proved to be a powerful tool in inves tigations of orbital polarization in scatterings [Fano and Macek, 73; Andersen, 87; Hippler et al., 86]. Using the language of the density matrix, the diagonal elements are related to the cross section, and the offdiagonal elements contain the information about the shape and rotation of the electron cloud. Usually, the density matrix is constructed to calculate the orbital polarization parameters for the electron cloud associated with the scattered atom. In this way one assumes that the electron cloud contains only the electrons in the orbitals of the scattered atoms and that the contribution from the target electrons is excluded. Such an approach is suitable for most experimental situations where the scattered atoms are detected far away from the target but before or in the course of decaying from their excited states. However, at short or medium distances, the atomic orbitals are mixed with the surface orbitals and the electron cloud associated with the atom contains also the contribution from the surface electrons. The above approach does not seem satisfactory from our theoretical point of view, since we are interested in following the time evolution of the electron cloud which requires to describe the orbital polarization of the atom at all distances. The definitions of the orbital polarization parameters should be consistent with that for an isolated atom as the distance between the atom and the surface becomes infinitely large. 30 In this Chapter, we use the electronic density matrix to define parameters characterizing the orbital polarization. Unlike other studies [Fano and Macek, 73; Nielsen and Andersen, 87], we start with the electronic density matrix of the full system, i.e., the atom plus the surface, and then define the orbital parameters for a subsystem which either is the scattered atom alone or is a molecular complex containing the atom and a part of the surface. Because the mixing between the atomic orbitals and the surface orbitals are taken into account, the electron cloud associated with the scattered atom contains also the contribution from the surface electrons. This permits us to examine the evolution of the electron cloud and gives a dynamic picture of the orbital polarization at all the distances during the collision. As the distance between the atom and the surface becomes very large, the contribution from the surface electrons becomes insignificant and our polarization parameters for the atom orbitals are asymptotically equal to those defined in other approaches. In Section 3.1, we describe two coordinate systems: the collision frame and the natural frame. The collision frame seems suitable to describe the scattering; the natural frame is more convenient to picture orbital polarization. In Section 3.2, using coordinate tensor operators, we define electric multipoles which provide a picture of the spacial distribution of the electron cloud. In Section 3.3, the alignment and orientation parameters are defined by use of the irreducible tensor operators and expressed them in terms of the density matrix. In Section 3.4, the polarization parameters are defined for a subsystem which contains either the 31 scattered atom alone or the atom plus a small part of the surface. This is useful for comparison with experimental data from spectroscopies. The connection between the orbital polarization parameters for a subsystem and that for an isolated atom is discussed. 3.1 Coordinate Frames Since symmetry plays an essential role in the study of scatterings and orbital polarization, it is the first thing one should look at when choosing a coordinate frame. Let's consider the symmetry of the porbitals. The three eigenfunctions of the angular momentum operator for J=1, Ip+i), Ipo) and Ipl) correspond to the magnetic quantum number m=1, 0 and 1 respectively. We can also use three real porbital functions Ipx), IPy) and Ip,) which are symmetric along x, y and zaxes respectively. The orbital jpz) = Ipo) has a negative reflection symmetry with respect to the XY plan and others have positive reflection symmetry. The collision plane is determined by the incoming velocity vector in and outgoing velocity vector bout of the projectile. For an atomsurface collision, the collision plane is perpendicular to the solid surface. The collision frame (Xc, Yc, Z7) is defined such that the XcZc plane coincides with the collision plane with the Z. axis perpendicular to the surface and that the Xc axis parallel to the surface, and the XcYe plane coincides with the solid surface with the Yc axis perpendicular to the collision plane, see Fig. 1. Because many collision systems have certain type of symmetry about the Zc axis or with respect to the collision plane, this frame is convenient for collision problems [Andersen, 86]. Zc, XN ,,ut Figure 1.1. The coordinate frames used for the description of orbital polarization of the electron cloud. The collision frame (Xe, Yc, Zc) and the natural frame (Xn, Yn, Zn) are shown. The incoming velocity vector U'i and the outgoing velocity vector o ut are in the collision plane which is perpendicular to the solid surface and coincides with the XcZ7 plane of the collision frame and the XnYn plane of the natural frame. Xc, YN CZ N 33 In the natural frame (Xn, Yn, Zn), as shown in Fig. 1.1, the XnYn plane is in the collision plane with the Xn axis perpendicular to the solid surface, the YnZn plane coincides with the solid surface with the Zn axis perpendicular to the collision plane. This frame seems to be "natural" to the analysis of the angular momentum transfer and orbital polarization, since the expectation value of the zcomponent of angular momentum is related to the orientation and alignment [Andersen, 86]. 3.2 Electric Multipoles One way to describe the shape of the electron cloud is to consider its electric multipoles. For this purpose we introduce a set of tensor operators in the Cartesian coordinates, J(0) = 1 1) = ri (3.1) I()= 28. 3rirj where i, j=x, y, z, and the super indices in the parentheses indicate the ranks of the tensor operators. The electric multipoles are defined as P(k) = (I(k)) (3.2) where I() represents the tensor operator defined in Equation (3.1), and the symbol < > indicates the ensemble average which can be expressed in terms of the density 34 operator A and the electronic density matrix PT in a basis set {~} TrI(k)p EE lP (Mk) (I )= Tr(i) ~~' (3.3) Tr(A) EEC P ,^ S; where the trace is over both orbital indices p, v and spin index 7, and Spy is the overlap matrix element in this basis. The denominator Tr(/) = PS, = = n (3.4) 7 J' 7 is the total electronic population, and n = P ,S (3.5) ;tV is the average electronic population for spin 7. The electric multipoles contain the information about the spacial distribution of the electron cloud. For k=0, p(O) = 1. For k > 0, P() =Tr(/ri) 1 P,(ri). (3.6) STr(~) =n is the component of the dipole moment P (since the electron mass is 1 au) and p(2) = Tr [A(r2i 3rr)] 1 (3.7) 1i Tr[A3] P2,(r2 3riri) (3.7) is the component of the quadrupole tensor of the electron cloud. The electric multipoles can be defined either in the collision frame or in the natural frame. It is also useful to define them in a bodyfixed "atomic frame" centered at the scattered atom, in which the axes are along the major axes of the electron cloud and the quadrupole matrix appears diagonal. 35 3.3 Alignment and Orientation Parameters Another way to analyze the orbital polarization is to examine the anisotropy of the distribution of the orbital angular momentum of the electron cloud. For the orbitals of certain total angular momentum J, if the distributions of the orbitals for all IJM) states are the same the cloud is isotropic; if the orbital distributions on these states are different, the cloud is said to be oriented; if the orbital distribution for IJM) and IJ M) are equal but differ for different magnetic number M the cloud is said to be aligned. To precisely describe these polarization properties, we define the orientation and alignment parameters through rotational moments [Fano and Macek, 73; Blum, 81]. Introduce a set of spherical irreducible tensor operators j() = 1 j(1) z (2) j= J 2 (3.8) 2 J1 = :FJ+(2Jz + 1) j(2) (3J_ 2 J2) where J is the angular momentum operator and Jz is its zcomponent, J+ and J_ are the raising and lowing operators. The irreducible tensor operator Jqk) satisfies the following relations [Zare, 88] [Jz, qk)] = qjk) (3.9) [J, J)] = [k(k + 1) q(q + 1)]J 1 (3.10) J [= J (3.11) We define the orientation vector O01) and the alignment tensor A) as qOfj) 1 Re(Jq) (3.12) X J(J + 1) Aq)(J) = Re(J+2)) (3.13) /J(J + 1) where < > indicates the ensemble average for a given J. Using the density operator and density matrix, the average of the irreducible operator can be written as j ,) Tr(3J')k) 1 (Jq) Tr( ) n E E ) (J (3.14) The vector Oql)(J) describes a preferred direction of angular momentum. The tensor A )(J) describes the preferred spatial distribution of angular momentum corresponding to each IMI. O0l)(J) is also called the orientation parameter. In the natural frame, it is the measure of the average component of angular momentum perpendicular to the collision plane, JL. 37 The meaning of these parameters is more transparent and the calculation is easier in the natural frame. But some collision problems are more conveniently dealt with in the collision frame. One can first construct the electron density matrix in the collision frame and then calculate orientation and alignment parameters in the natural frame after a frame rotation. The density matrix contains all the information on the electron states and can be used to define other parameter to describe the orbital polarization. For example, we can define the alignment angle a(J) which measures the angle between the major axis of the J component of the electron cloud and the Xn axis. Choosing the basis set such that p = JM, we define, for J=1 [Nielsen and Andersen, 87] a(1) = [r + arg(P11,i)] (3.15) and for J=2 a(2) = [7r + arg(P22,2o + P20,22)] (3.16) where M = M. We can also define, for J>2, the octupolar angle shift r7(J). For instance, for J=2 7(2) = arg(P22,2) a(2) (3.17) The parameters Lj, a(J) and q(J) have practical meanings and can be compared with experimental date to test the theoretical model and methods used in obtaining the density matrix [Panev et. al., 87; Andersen et. al., 86]. 38 3.4 Multipoles and Alignment and Orientation Parameters in a Subsystem The parameters defined above describe the orbital polarization of the whole system. Sometimes we are more interested in the orbital polarization of a subsystem. For example, in atomsurface collision experiments, spectroscopies detect scattered atoms at the distances where the influence of the surface is negligible and only the polarization of the atomic orbitals is of interest. Even in the theoretical description of the collision states at short distance, only a small part of the surface should be taken into account, since the atom usually interacts mainly with a small region of the surface. Let's divide a system into a primary region and a secondary region which are labeled by p and s respectively. In the case of atomsurface collisions, the primary region can be chosen to contain the scattered atom and a small part of the surface and the interaction between the atomic states and those in the secondary region is assumed to be small. In the natural frame the ensemble average of Jq for the whole system can be written as (+ (Trj(k)) + ) pEa vEp pE E Ep vE s (3.18) The last term comes from the secondary region and is not of interest. The second and third terms contains the mixing between the primary region and the secondary region. 39 Because we are interested in the orbital polarization of the scattered atom we define a set of irreducible operators L(1)O) 1 L(1) = L 0 L L2 1L2 (3.19) 2 f L( = FL(2Lz 1) r(2) 1 2 2 S(3L2 L2) where L is the orbital angular momentum of the scattered atom, L its z component, and L. and L the raising and lowing operators. The average of L) is (L)) TrA (PLk))  (k) TrA()) = AE E E P (L()) (3.20) pEA vEp where trace is over the subspace of atomic states, and nA = TrA(/) = E E PSJ (3.21) 7 pEA vEP is the average electron population on the atom. The contribution to the electron cloud associated with the secondary region is not included since the mixing between the second region orbitals and the atomic states is small. In this way we can define the orientation and alignment parameters for the scattered atom as OI)(L)A = 1 Re(L1)) (3.22) V/L(L + 1) 40 A (L)A = i Re(L ) (3.23) /L(L +1) ( The alignment angle and octupolar angle shift can also be defined in the similar way. At small distance, the surface electrons contribute to the atomic orbital polarization. At large distance, mixing between the atomic orbitals and the surface orbitals becomes small and the summation of v vanishes unless for v E A, the polarization parameters are equal to those defined for an isolated atom. CHAPTER 4 LINEARIZATION OF TDHF EQUATIONS The TDHF equation developed in Chapter 2 can be solved straightforwardly only for simple systems requiring small basis sets. If the system is complicated, the TDHF equations can instead be solved by some linearization procedures. However, one has to be very careful when developing a linearization method for cases where the density matrix oscillates rapidly with time, such as in collisions accompanied by a charge transfer through a nearresonance process, since the solutions of the linearized TDHF may converge slowly or not converge at all. In this chapter, we describe a local time linearization procedure for solving the TDHF equations for the nearresonance charge transfer process. It tackles the rapid variation of the density matrix by defining a timedependent reference density matrix PY(t). In section 4.1, we develop the local time linearization procedure for a general case and obtain a linearized equation for P7 (t), which is a solution for P(t) when the ion position is fixed in space, and a linearized equation for matrix QT(t), which is the first order change of the density matrix due to the timedependent perturbation caused by the nuclear motion. The formal solutions for PoT(t) and QY(t) are obtained with use of an exponential transformation. In section 4.2, we apply this procedure to a system in which the electronelectron 42 Coulomb interaction is ignored. For this special case, the coefficient matrices in the equations of P1(t) and QT(t) are time independent and the analytical solutions can be obtained. 4.1 The Linearization Procedure In a collision which involves charge transfer through a nearresonance process, electrons jump back and forth among several close energy levels and this leads to a rapid oscillation of the density matrix. In the contrast, the change of the density matrix under the timedependent perturbation caused by the motion of the nuclei is relatively slow. We separate the two time scales in what follows. For a small time interval to 5 t < tl, we define a matrix P07(t) to satisfy the equation iP&(t) = So'gP0of'() PO'F(t)FytSo1 (4.1) and the initial condition P'(to) = P'(to), with FP = F"(to) and So = S(to). Hence Pt(t) can be interpreted as a solution of the TDHF equation when the ion nuclei are fixed at RA(to). The evolution POT(t) oscillates with time reflecting the charge transfer among the energy levels with frequencies proportional to the inverses of the energy differences of energy levels. On the other hand, if the time interval is small, the change of the density matrix due to the nuclear motion is small. To get a linearized equation for the this change, we define QT(t) = PT(t)Po'(t) (4.2) AS1(t) = S(i) So1 (4.3) and AF (t) = f(t) F (4.4) From (2.32) we have AF(t) = Ai + AeGf(t) + AnG(t) (4.5) The matrices AeG and AnG are related to the electronelectron interaction and are given by AGfeG(t) = C Q7(t)(VjIIIA)io + Q'A(t)( iIA)to (4.6) and A.G1,(t) = _P.(t)A(pvl JA), + P,'(t)A(pvl A)t (4.7) IA where (pVIcA)t and (yjllcA)t are the Coulomb and Coulomb symmetrized 2electron integrals respectively when the scattering atomic core is at RA(t). Inserting Eqs. (4.1)(4.5) into Eq. (2.31), we obtain to first order in Q7(t), iQ'(t) = Sol' Q'(t) Q'(t)FgtSo1 +So1AHi(t)Po(t) pO'(t)Af(t)tSo1 +Po'(t)FatAS1(t) ASz(t)qFPo(t) (4.8) +So1AeGr(t)PO'(t) po'(t)AeG7(t)tSo1 + So AnG'(t)POr(t) PO'(t)AnG'(t)tSo1 44 To introduce a shorthand notation, we use a double index notation to replace the single index notation in the following. Let upper case Roman letters refer to double indices, then, K = (pv) and L = (cA). With this notation, matrices P07, Q7, AeGY and A,G7 become column matrices P07, Q7, Aeg and And. Defining two square double index matrices X(t) and y(t) such that XKL(t) = Xv,s(t) = (pvjjIA)t (4.9) YKL(t) = Yp,,,(t) = (pvlKA) (4.10) we can rewrite Eqs. (4.6) and (4.7) as G(t) = E[XKLQ'(t) + YoKLQ'(t)] (4.11) L 9(t) = [AXKLPL'(t) + AYKLP''(t)] (4.12) L or in the matrix form, AeG7= XoQ + YoQ' (4.13) Angr = AX Pp0 + Ay p07' (4.14) where Xo = X(to), Yo = Y(to) and AX = X(t)  Xo (4.15) Ay = y(t) yo Taking advantage of using the double index notation, the linearized TDHF equations (4.1) and (4.8) can be written in a neat form ijorf =A"op' + .A'p0o' _ po'Art por'A7't iQ =A4 "Q+A"''Qf'' QA Q'A7'7t + V7 and (4.17) (4.18) where A" and A"' are square matrices in double index notation or four index matrices in single index notation, AKL= A",, = SoI s.,6A, + XOC9,APS) C9 (4.19) KAL ItvKAX = SoCYCOrCAP'. C17 and the column matrix Dr is a driving term with elements S= [So1(Af+ AG)PoY Pof(AHt+ AnGt)So1 +AS'1j~p'o PO'ftAS1], If we collect Po7 and Po0' into a column matrix PO, Q7 and Q7'' matrix Q, V7 and P7' into a column matrix V and let (4.20) (4.21) into a column ( A"' Y 77' A = Ay' A^' 4 (4.16) (4.22) then Eqs. (4.17) and (4.18) become ilp = Apo poAt (4.23) and iQ = AQ QAt + E (4.24) These are the linearized TDHF equations for the density matrix in to < t < ti. It should be noted that the matrix A is generally timedependent. To solve the above equations, we let U(t, to)=T[exp i A(t')dt' (4.25) U(t, to)=T[exp i A(t') dt' I (4.26) where T [exp{ ... }] denotes a timeordered exponential expansion. We consider transformations Po(t) = U(t, to)lP(t)U(t, to)t (4.27) and Q(t) = U(t, to)Qv(t)U(t, to)t (4.28) Replacing Eq. (4.27) into the left hand side of Eq. (4.23) and Eq. (4.28) into the left hand side of Eq. (4.24) we have iP0(t) = APO PoAt + iU(t, to)P,(t)U(t, to)t (4.29) and iQ(t) = AQ QAt + iU(t, to)Qp(t)U(t, o)t (4.30) Comparing them with the right hand side of Eqs. (4.23) and (4.24), we obtain iP(t) = 0 (4.31) and iQr(t) = iU(t, to)DV(t) [U(t, to)t ] (4.32) Their solutions are P(t) = Po(to) = P(to) (4.33) and t Qp(t) = Qp(to) i fu(t', to)'D (t')[U(t', to)t]dt' (4.34) to Using Eqs. (4.27) and (4.29) and noticing that Q(to) = 0, we obtain the formal solutions for Po(t) and Q(t) Po(t) = U(t,to)1o(t0) [l(t, ti)t] (4.35) and t Q(t) = iU(t, to) U(t', to)_ (t') [U(t', to) t dtU(t, to) i to (4.36) = i f u(t, t')D(t')U(t,t') dt' to 48 4.2 The Case without ElectronElectron Interaction In this section, we will show the use of our linearization procedure by applying it to a special case where the electronelectron is ignored. In this case G = 0 and F = Hf. Since the different spins are then uncorrelated, we return to the notation involving K, A, p and v, and drop the spin index 7 to write Eqs. (4.23) and (4.24) as iPo(t) = WP(t) p(t)Wt (4.37) iQ(t) = WQ(t) Q(t)Wt + D'(t) (4.38) w = SfIo Ho = fH(t0) (4.39) (4.40) Unlike in the independent. D'(t) = S1AH(t)P0(t) P t (t)S (4.41) +Po(t)IAS1(t) AS1(t)HoPo(t) general situation, the coefficient matrices W and Wt are time To solve this equation let us consider transformations P(t) = exp[iW(t to)]Po(t)ezp[iWt(t to)] where (4.42) and Q(t) = exp[iW(t to)]QD(t)exp[iWt( to)] (4.43) Inserting them into Eqs. (4.37) and (4.38) respectively, we get equations for PD and QD iPO(t) = 0 (4.44) and iQD(t) = exp[iW(t to)]D(t)exp[iWt(t to)] (4.45) The latter has a solution of t QD(t) = i exp[iW(t' to)]D'(t')ep[iW(t' to)]dt' (4.46) to since QD(to) = 0. From the above equations we can write the formal solutions for PO and Q P(t) = exp[iW(t to)]Po(to)expiWt(t to)] (4.47) and Q(t) = i exp i(t t') D' (t')exp[iW ( t')]dt' (4.48) to We assume that the matrix W can be diagonalized by a linear transformation, that is, there exists a matrix L such that, W = LwL1 (4.49) where w is a diagonal matrix. Then Q t(t) = i j Lk [ (L)kAkApl(t to )(L) L (4.50) where to The formalism developed in this section can be directly applied to systems with a few electrons. We have used it to calculated charge transfer in H+H' collisions [Runge et al., 90]. However, it must be modified when applied to extended systems. In Chapter 6 we describe a partitioning procedure which simplifies the treatment of the scattering of atoms by surfaces. CHAPTER 5 PARTITION OF EXTENDED SYSTEMS For an ionsurface system which involve a great number of electrons, the TDHF equations obtained in Ch. 2 usually require a huge basis set. Some technique for truncating the system as well as the basis set must be developed to solve the equations. Partitioning techniques have been used to treat various kinds of extended systems [Lowdin, 70; Ying 77; Williams et al, 82; Kirtman and de Melo 81; de Melo et al., 87; McDowell, 82, 85]. The main idea is to divide the system into two regions and to employ some relatively accurate treatments to deal with electronic interactions in a primary region, which is of main interest, while using certain approximations in a secondary region. The difference among a variety of methods is in the treatments of the secondary region and of the interaction between the two regions. For some systems the interaction between the two regions is small and the secondary region can be treated as a small perturbation. For other systems the secondary regions are often approximated by some simpler systems, for example, in studying Kondo effect [Kondo, 69; Heeger, 69; Anderson, 61], an approximation is to treat the conduction electrons surrounding isolated magnetic ions as a free electron sea. 52 For ionsurface systems, it has been found that, during collisions, electron rearrangement occurs mainly within local regions of surfaces [Grimley et al., 83; McDowell, 85]. Based on this observation we design a partitioning technique to truncate the ionsurface system into two parts and treat electronic interactions and electronic state evolution in these parts differently. This partitioning technique is suitable to deal with timedependent processes in extended systems with a strong local coupling. It requires a set of localized basis functions, but does not depends on the concrete form of the functions. All the formal derivations in this chapter are done with localized basis functions which can be of any form. In Section 5.1 we describe our partitioning method and derive the effective linearized TDHF equation for the density matrix in the primary region, which includes a driving term containing the coupling with the secondary region. In Section 5.2, we introduce an approximation to the secondary region and the couplings between the two regions, which enables us to solve the effective TDHF equations in the primary region. We apply this method to a simple case, where the electronelectron interaction is ignored, and obtain analytical solutions for the density matrix in the primary region. 5.1 Partition Procedure We divide an ionsurface system into two regions: (i) the primary region that consists of the scattering ion and a small impact area of the surface, where the perturbation by the ion is strongly felt during the collision; (ii) the secondary region which consists of the remainder of the surface, where the influence of the 53 ion is assumed to be felt indirectly through the coupling between the electrons in the two parts of the surface. Considering a set of localized basis functions for the solid and the ion, and letting p be the index of the primary region and s the index of the secondary region, the matrix PO is written as Po0. P pO) (5.1) Psp s and other matrices are partitioned in a similar way. In the case where the electron electron interaction is ignored, Eqs. (4.37) and (4.38) are then split into eight equa tions for Pop(t), Pos(t), POp(t), Ps (t), Qpp(t), Qps(t), Qp(t) and Q.(t) iP (t) = W P (t) Po(t)W (5.2) +WpP(t) Po (t)w iP (t) = WppPo(t) Pp(t)Wt (5.3) +WpP(t) PO(t)W8 iP (t) = WP p() Pp(t)W (5.4) +WP p(t) PO)(t)Wsp iP(t) = W P(t P (t)W p (5.5) +W PO (t) P (t)W iQpp(t)=WppQpp(t) Qpp(t)Wp (5.6) +WpQsp(t) Qps(t)Wtp + D' (t) iQpS(t)=WppQps(t) Qpp(t)Wp, (5.7) +WpQss(t) Qps(t)Wt + D',(t) iQsp(t)=WspQpp(t) Qsp(t)Wpp (5.8) +WsQsp(t) Qss(t)Wtp + D'p(t) iQs(t)=WspQps(t) Qsp(t)Wp, (5.9) +W ssQss(t) Qss(t)W, + D (t) where Dp(t), Dps(t), D p(t) and D,(t) are the submatrices of D'(t) defined in Eq. (4.21). For example, the matrix Dpp(t) is given by D',(t) = (So)ppAHfpp(t)Pp(t) Pp(t)AfIp(t)(Sol)p +(Sol) ppAHps(t)Pp( t)Ps,()A p()(Sol) +(Sol)p AItp(t)P p(t) Pgp(t)Ap' (t)(Sp1)p +(So l)p AHIs(t)P p(t) P,(t)AHL(t) (So1) (5.10) +Pp(t)ii (ASl)pp(t) (AS) pp(t)topp p (t) +Ppp(t)HI (AS1)s(t) (AS1)ps(t)IospPPp(t) +Pps(t) p(AS1)pp(t) (ASl)pp(t)opsPsp(t +P0s(t)HltIt(ASl )p(t) (AS')ps t)o P (t) The primary region where charge transfer, energy transfer and other collision phenomena take place is our main interest. But, solving the effective equations in the primary region, Eqs. (5.2) and (5.6), is not any easier than solving the original linearized TDHF equations, since they contain the density matrix elements in the secondary region and have to be solved simultaneously with the equations for density matrix elements involved the secondary region. Obviously, some 55 approximation has to be made in the secondary region to make this partition procedure practical. 5.2 The Approximation in the Secondary Region When considering the secondary region, we notice that during the collision, the charge exchange is mainly restricted to the primary region. To understand it, let us consider the interaction and the evolution of the electronic states in the two parts of the surfaces, i.e., a small impact area on the surface which is close to the scattering ion during the collision and the remainder of the surface. During the collision, the impact area experiences a strong and timedependent perturbation from the ion and the evolution of the electronic states in the area is significantly altered from that of the unperturbed surface. As the result, charge in the impact area may move to or from the ion. In the remainder of the surface, the direct interaction with the ion is much weaker, the electrons in this area feel the effect of the ion indirectly through their coupling with electrons in the impact area. Because the time scale of charge rearrangement in the solid, 7ear, is much longer than the duration of the collision at the surface, reol, at the collision energies of present interest, the electrons in the remainder of the surface do not have time to adjust to the charge rearrangement in the impact area, and the electronic states follow approximately the evolution pattern of an unperturbed system. Thus, it is reasonable to assume that the secondary region is basically unperturbed by the ion during the collision. In this approximation, the submatrices of Q associated with the secondary region, Qps(t), Qsp(t) and Qss(t), are found to be negligible 56 compared to Qpp, and the submatrices Pgp(t), P0p(t) and P((t) can be replaced by those of the uncoupled system, Pps(t), Psp(t) and Ps,(t). Other submatrices are treated in the similar way. For an unperturbed surface, the coefficient of the molecular orbital has the form of 4i,(t) = Ei,(tin)exp[ie(t ti,)] (5.11) where ei is the energy of the ith orbital. From Eq. (2.11) Pij(t) = 1 c41(t)A(t) = I c> (tin)ci4(tin) i=occ i=occ (5.12) = Pi,(ti,) As the result of these approximations, we need to solve only the effective equations in the primary region iPgp(t) = WppPpp(t) POp(t)W p (5.13) and pp(t)=WppQpp(t) Qpp(t)W p + pp(t) (5.14) The equation for Qpp(t) is coupled with the secondary region through terms within Dpp(t) which is Dpp(t) = (S)ppA pp(t)Pp(t) POp(t)A t pp(t)(So1)pp + Pp(t) tpp(AS1)p(t) (AS1)pp(t)ftopppp)(t) (5.15) +Pps( tin)0sp(AS1) P(t) (AS1)p(t)iops tin ) Given the matrix elements S' and Ho, the driving term Ib'p(t) can be calculated and the density matrix in the primary region can be obtained either analytically or numerically by using Eqs. (4.50) and (4.51). It should be pointed out that the validity of our approximation in the secondary region depends on two factors, the ratio of Tcol/Tre. and the size of the primary region. For metal surfaces and colliding atoms with their kinetic energy above leV, Tcol/Trer is small and our partitioning procedure should apply. As for the size of the primary region, the larger the primary region, the more accurate the result. But increasing its size will significantly enlarge the basis set in the primary region and increase the computing time. In Ch. 9 we will discuss the effect of the primary region size in more detail. CHAPTER 6 ELECTRONIC BASIS FUNCTIONS AND MATRIX ELEMENTS The use of the localized basis functions in the present study appears nec essary for our partition procedure introduced in Chapter 4, however, the real reason underlying it is the local nature of the phenomena associated with sur faces. For example, localized surface states, surface adsorptions, chemisorption and atomsurface collisions have their effects basically confined to small regions. Experimental results have also provided evidence that charge densities and local densities of states, while remarkably different from those of bulks on the top lay ers of surfaces, quickly recover to the bulk values on deeper layers. For example, the electron energy distribution results obtained from the ion neutralization spec troscopy, which probes only the top layer of atoms, show a qualitative difference from those of the bulk [Hagstrum and Becker, 73], but the results from the ul traviolet photoemission spectroscopy, which probes about four atom layers, are essentially dominated by the bulk density of states [Eastman and Grobman, 72]. Localized basis functions have been used in calculating electronic structure of bulk materials; the examples are Wannier functions and atomiclike basis functions used in the tightbinding method. The local nature of the surface phenomena suggests that some of the electronic behaviors of surfaces could 59 be more advantageously and conveniently described in terms of localized basis functions which return to the form of the bulk functions on deeper layers. Recently, some works have reported using atomic orbitals to calculate surface electronic properties. Smith and his colleagues use a set of atomiclike basis functions to calculate surface band structure for transition metals [Smith and Gay, 75; Smith et al., 80; Alinghaus, et al., 80]. Their basis functions are constructed by fitting them to the solutions found by solving the KohnSham equations [Kohn and Sham, 65]. Kohn and Onffory introduce the concept of the generalized Wannier functions which can be used in the theoretical study and calculation of electronic structures for extended systems lacking translational symmetry [Kohn and Onffroy, 73]. These functions are localized on the lattice sites and, as they result form a unitary transformation of the eigenfunctions of the system, are orthonormal and complete. Gay and Smith [Gay and Smith, 74] proposed a variation method for determining the generalized Wannier functions without having to first calculate the eigenfunctions of systems. This chapter is devoted to constructing localized functions as our basis function set for atomsurface systems, which consists of a set of localized functions for the surface and a set of atomic functions for the atom, and to calculating the relevant matrix elements in this basis. Although the localized basis functions used by Smith and his colleagues are successful in giving good results for transition metals, their calculation and application are complicated. As a test of our time dependent molecular orbital method and the partition procedure, we feel that a set 60 of simple localized functions would be adequate and hence choose the generalized Wannier functions as a part of our basis set. In Section 6.1, the definition of the generalized Wannier functions is intro duced and a variation procedure for determination of these functions is described. In Section 6.2, we apply the variation procedure to a jellium surface to obtain the generalized Wannier functions. For convenience in the calculation of the matrix elements, the generalized Wannier functions are written as linear combinations of Gaussians. In Section 6.3, the atomic functions for the atom are constructed by using a pseudopotential for the Sodium atom core. In Section 6.4, we calculate the matrix elements of the overlap and Hamiltonian. Some of the details of the calculation are given in Appendices A and B. 6.1 Generalized Wannier Functions 6.1.1 Definition of Generalized Wannier Functions (GWFs) We first consider an infinite periodic system with a Hamiltonian H0. Its eigenfunction (r) is a Bloch function labeled by the wave vector kI and band index i and satisfies the Schrodinger equation !ftoo(r) = Ceq() (6.1) where ie is the eigenvalue associated with i and k. The functions 0(rf) are then collected in a row matrix V = (00 ,0,.. ) (6.2) 61 and the above Schrodinger equation has a matrix form joi = iPEi (6.3) where E is a diagonal matrix with diagonal elements . For each energy band of a periodic lattice, Wannier functions are defined by [Wannier, 37] w9() = ( ) eiR'o (6.4) kEBZ where N is the number of the lattice points, Ri is the lattice vector labeled by a vector index in, the summation is over the Brillouin zone. The Wannier function w%(rF is localized about the lattice point RA. In matrix form, the above transformation can be written as w? = Uot (6.5) where wi= (wi w, ,. ) (6.6) and 1  U = expilk. ) (6.7) Since matrix UO is unitary, the Wannier functions are orthonormal and are an alternative set of basis functions in electronic structure calculation [See, for example, Slater and Koster, 54]. 62 For a semiinfinite solid with its surface parallel to the xy plane, its Hamil tonian H and eigenfunction ~i(F), which is characterized by the band index i and the wave vector k, satisfy the Schrodinger equation Hiq(rfj = e.i.ik(r) (6.8) where eE is the eigenvalue. Collecting the functions qig(f in a row matrix Ti = (i, ) (6.9) the above Schrodinger equation can be written in matrix form Ii = ,iEi (6.10) where Ei is a diagonal matrix with diagonal elements cg. In the following, we drop the band index i with an understanding that all the functions and matrices involved are referred to a particular band. As a results of breaking translation symmetry in the direction perpendicular to the surface, the eigenfunctions &(r) are no longer of Blochtype and the functions defined in Eq. (6.4) are no longer meaningful. Kohn and Onffroy point out that it is possible to construct a set of orthonormal functions for systems with defects, called generalized Wannier functions [Kohn and Onffroy, 73]. Following their idea, we define the generalized Wannier functions of an electronic band as a unitary transformation of the eigenfunctions of the system, that is, (6.11) where w = (wit, ., ) (6.12) is the row matrix of generalized Wannier functions and U is a unitary matrix to be determined. Obviously, w is orthonormal, i.e., (wlw) = I (6.13) It has been proved that Generalized Wannier functions are localized about the lattice sites R/ and are exponentially decaying away from their center [Kohn and Onffroy, 73]. The wellbehaved localization property of the generalized Wannier functions combined with their orthonormality make them particular useful to investigate the electronic structures or localized phenomena for nonperiodic systems. 6.1.2 Generalized Wannier Functions as Linear Combinations of Gaussians For the purpose of practical applications we write the generalized Wannier functions as linear combinations of Gaussian functions centered at different sites. For a certain energy band only Gaussian functions with proper symmetry should be used to construct the generalized Wannier function. Thus, W.h() = i Bni a g ar( ari, r) (6.14) where BnA is a coefficient, a=(nlm) is a compound index, and 9gt(t A, A f ) = brYim(0, p)rIep,(Ai,)' (6.15) 64 is the a type Gaussian primitive function centered at Rw,, where Y.m(0, p) is the spherical harmonic functions and be a normalization factor. Letting (6.16) (here we have dropped index a in g since a is the same for all the Gaussians for a given band), Eq. (6.14) can be written in a matrix form w=Bg (6.17) where (B)jj = Bjr (6.18) Substituting Eq. (6.17) into the normalization relation Eq. (6.13) yields BtGB = I where G = (glg) is the overlap matrix of Gaussians, we find that a possible choice for B is B= G where the matrix G satisfies GG = I G2GG = I (6.19) (6.20) (6.21) (6.22) g = (gSaii, ga2ii, ) 65 6.1.3 Determination of Generalized Wannier Functions The generalized Wannier functions can be constructed from the eigenfunctions ~(r) according to their definition (6.11). However, one of the important features of the generalized Wannier functions is that they can be directly determined from the system Hamiltonian without having to know the original eigenfunctions [Kohn and Onffroy, 73; Gay and Smith, 74], which allows freedom to construct them as desired. In this study we use a variational procedure [Gay and Smith, 74] to find the generalized Wannier functions w. Since E = (lfftPA ) (6.23) the total energy of a band EB = Tr(E)= Tr(({lIHI)) ( \ (6.24) = Tr((wlHfw)) = E[w] is a functional of the w's, and we have used in the second line the orthonormality Eq. (6.13). The variational principle proposed by Kohn and Onffory [Kohn and Onffroy, 73] says that the total energy of a band attains its minimum if a correct set of w's is used. Thus the generalized Wannier functions can be determined by minimizing EF[w], that is, min {E[w]} = (6.25) {w) while being subject to the constraint Eq. (6.13). 66 Once w is known we can insert Eq. (6.11) into the eigenequation of the system to have fHwU = wUE (6.26) or, using Eq. (6.13), (wlfllw)U = UE (6.27) This equation can be used to find the matrix U as well as the eigenfunction E(rF). By writing the generalized Wannier functions as linear combinations of Gaussians, the functional EB[w] becomes a function of /3,'s, 1(/3 2, .), and the problem of searching for a set of w becomes one of finding a set of optimized #,a's, i.e., min [(/3,1, /2, )], g/2, (6.28) One may recall that the above procedure for writing the generalized Wannier functions as linear combinations of Gaussians is similar to those used in molecular orbital calculation. However the difference is that the Gaussian function or the exponent /#, depend on the site position R4. For a surface, due to the periodicity in x and y directions, all the generalized Wannier functions on the same layer are the same and the exponent /3, changes only with layers. From now on we will denote /3, by /m,, where mz is the z component of the vector index ri, to indicate that it is the exponent for the mth layer. 67 Before proceeding to find Im, 's we notice an asymptotic behavior of the generalized Wannier functions; they approach the Wannier functions of a periodic system in an exponential manner, as proved by Kohn and Onffory [Kohn and Onffroy, 73]. This indicates that only a small number of the generalized Wannier functions need to be determined for nearperiodic lattices which loose translation symmetry only in local areas or in certain directions, such as imperfect crystals and surfaces. For a perfect bulk, the Wannier functions w (r)'s on all the sites R are the same and have the same exponent parameter po. For a semiinfinite lattice with its surface at me=l, the generalized Wannier function w,j(r) approaches w,(r's and /3,n approaches 30 when /R is deep inside the surface. There exists a cutoff M such that when mz > M it is found approximately that 'm. = /30 (6.29) The variation procedure (6.28) now becomes min [(/3i, A32, *)] = min [(/1, 32 3M)] 31, P2, ,PM (6.30) 01,02, 01,92,,8m Thus only M variational parameters, 31, #2, *, 3M, must be determined. The choice of the number M depends on the system and the accuracy requirement of calculations. 6.2 Generalized Wannier Functions for a Jellium Surface The method developed in the previous section can be applied to construct generalized Wannier functions for any system for which the effective one electron 68 Hamiltonian is known, although the calculation of the three dimensional matrix L G is a little timeconsuming [Gay and smith, 74]. In this section we apply the method to a simple model, a jellium slab with a step potential. A jellium is a system with the atom core lattice replaced by an uniform positive charge background and, as a result, it completely ignores the effect of the atomic lattice. Because of its simplicity in the electronic structure calculation, jellium is often used in condensed matter physics as a testing model for theory or methods. We also assume that there is only one band and use Is type Gaussians to construct generalized Wannier functions. We fit the band parameters such as Fermi level, and the energy of the bottom of the band to those of W(110). For a finite jellium slab with a thickness D in zdirection and a step potential {0 z>0) V() = Vo D < z < 0 (6.31) 0 z < D the Hamiltonian is H(F) = V' + V(r) 2 (6.32) = HZ(z) + Hy(y) + HZ(z) where HZ(x) and HY(y) equal to the kinetic energies in x and y directions and 1 a2 HZ(z) = + VZ(z) (6.33) with VZ(z) = V() (6.34) Labeling x, y and z components by superscripts x, y and z, the eigenfunction b&(r), which is assumed to satisfy cyclic boundary conditions in x and ydirections, and eigenvalue Er have the forms of O() = q4(z)x (y)qS,(z) (6.35) and EE = Ek, + Ek, + Ek, (6.36) respectively. In Appendix C we describe the evaluation of the Fermi energy. We now consider an artificial cubic mesh in the slab, which contains N,=NxNy~N sites with a distance d between sites, where Nx, Ny and Nz are the numbers of sites in x, y and zdirections, respectively. The generalized Wannier functions for the slab are associated with sites and are written as wg(r) = w,(x)wn,(y)w,(z) (6.37) Recalling Eq. (6.11), the matrix U can then be written as a product of x, y and z component matrices U = Ux'YUZ (6.38) with Ux, Uy and Uz lead to the transformations k() = UX ( t), (x) (6.39) k, w,(y) = yt) q o(y) (6.40) ,(z) = UZtk, (z) (6.41) k, 70 respectively. For this system, therefore, a three dimensional calculation is sim plified to a one dimensional one. For this system, the x and ycomponents of the eigenfunctions are 1 ,= exp(ikzxz) 1 (6.42) k,Y = 1 exp(ikxz) 4 (Ny 1)d the generalized Wannier functions in the x and ydirections are just the Wannier function for periodic lattice 1 sin 1 W =,() = (6.43) (Nn 1) d sin [r(sn,d) 1 sin *' dd) W (y) s (6.44) S(Nz 1) sin r(Ynyd) L (N1)d and the component transformation matrices Ux and UY are 1 =Un N ,= exp(ik.nxd) (6.45) Uk, = 1N lexp(ikynyd) (6.46) In the z direction, the generalized Wannier function for nzth layer, wn (z), is written as linear combination of Gaussians ,. (z) = B~.B,,,mg (,, z) (6.47) m, 71 where the normalized one dimensional Gaussian function in the z direction g,(3m,,,z) is centered at Zm. on mzth layer with an exponent /3m, and Bz = (GZ)7 with Gz = (gzlgZ) the overlap matrix of Gaussians in the z di rection. As discussed in the previous section, only the exponent parameters on the first M layers are assumed to be different; that is, I1 = 3N,, 32 = #N,1, *,/M1 = #N/M+2, but /M = /3Mi = = fN.M+I. In the x and y directions, the Wannier functions on the mzth layer are written as linear combinations of one dimensional Gaussians with an exponent 3m,. For example, the Wannier function in the x direction is wn',(x) = > Bmm,(m,) 1,(fm,, x) (6.48) m, where the superscript mz for the function is used to indicate the layer the function is on. The coefficient matrix element Bnzm,(am,) is obtained by the least square principle [Shavitt, 63] by fitting Eq. (6.48) with Eq. (6.43). The details of the calculation are given in Appendix A. From Eq. (6.24) the total energy is EB = (wlHlwM) = Z E Z ji (P(fm( )[HIgi( (6m( )) 649) = n(?, #2, * *, ; ) with the details of the calculation given in Appendix B. Minimizing 1 respect to the #'s gives a set of optimized O's. 72 To determine the matrix Uz, we substitute Eq. (6.41) into the eigenequation of k,(z) HEZ((z) = Ek, 0(z) (6.50) to have SUk.m.HzW (Z) = Ek, Uk.m.,<(z) (6.51) m, m, Multiplying it by w((z)* and integrating over z give the equation for transfor mation matrix Uz SUk,,m (wu,. IHZw ) = Ek,Uk,,, (6.52) 6.2.1 Results We have calculated the generalized Wannier functions for a finite jellium slab modeling W(110), with the depth of the potential well Vo=0.846 au. The distance between sites, d, is not equal to the lattice constant of W, and is chosen to be d=2.61 au so that the average number of electrons in the volume associated with a site is one. The optimum exponents of the Gaussians, fm,, which are determined by minimizing the total energy Eq. (6.49), are listed in Table 6.1. It has been found that the optimum exponent decreases for the first couple of layers but changes very little after the third layer. This shows that the effect of the surface on the generalized Wannier functions is basically limited to the several top layers. Gay and Smith also find that the exponents are different from that 73 Table 6.1 The optimum Gaussian exponents for the first three layers. Layer First Second Third Pm,, 0.53 0.52 0.30 of the bulk only on the top two or three layers even though they use different potentials [Gay and Smith, 75]. The generalized Wannier functions for the first four layers are plotted in Figs. 6.1, 6.2, 6.3 and 6.4, assuming #m,=0.30 for mz > 3. These figures show that the generalized Wannier functions are well localized about their centers with long decaying tails. For the first two layers, the effect of the surface is significant, the generalized Wannier functions are asymmetric and their tails oscillate little. On moving to the deeper layers, the generalized Wannier functions are more and more characteristic of an infinite bulk. For the fourth layer the tail of the generalized Wannier function towards the center of the slab is very similar to that of the bulk Wannier function. Although the surface is a severe perturbation, our results show that the generalized Wannier functions rapidly become the bulk Wannier functions when moving away from the surface, which is in agreement with the theoretical prediction for the generalized Wannier functions [Kohn and Onffroy, 73]. To exam the validity of using generalized Wannier functions as basis functions, we compare the eigenvalues and wavefunctions calculated from these functions with the exact ones. Table 6.2 lists both exact and calculated eigenvalues, Figs. 6.5, 6.6 and 6.7 depict both exact and calculated wavefunctions for Ek, =0.8444, 74 0.6906 and 0.3990 for states jz=1, 10, 17 respectively. The agreements are good for all the energies, and especially for lower energies. Since the pointbypoint comparison of wavefuntions is probably the most stringent test of the accuracy, the good agreement in our calculation will assure the accuracy in the calculation of observables such as charge density. In the x and y directions, the Wannier functions w ,(x) and wl,(y) on the mzth layer are approximated by linear combinations of Gaussians with the mzth layer exponent fm,. The accuracy of the approximation depends on the number of Gaussians used. It is found that for the larger fm,, (close to the surface in our case) more Gaussians are needed to achieve the same accuracy. Figs. 6.8, 6.9 and 6.10 show the exact and calculated Wannier functions for the first three layers when 21 Gaussians are used. The calculation and application of generalized Wannier functions for a jellium slab reveal some nice features of these functions. Firstly, they are well localized about their centers with their tails characterizing the translation symmetry of the systems. The functions rapidly recover to the bulk Wannier functions away from defects or surfaces, thus only a few generalized Wannier functions need to be determined. Secondly, the generalized Wannier functions can reproduce the wavefunctions and eigenvalues to a satisfactory accuracy. These attractive features make generalized Wannier function suitable for electronic structure calculations for surfaces or other systems with a broken translation symmetry, and makes them useful in studies of local phenomena. 0.75 0.50 0.25 0.00 0.25 12.0 9.0 6.0 3.0 0.0 3.0 z(D) Figure 6.1 Generalized Wannier function wz (z) on the first layer of a jellium slab with a square well potential and a potential depth Vo=0.846. The surfaces of the slab are located at z=0 and z=21d, where d=2.61 au is the distance between layers. 0D 0C~ 0.75 0.50 0.25 0.00 0.25  I, I I 12.0 9.0 6.0 3.0 0.0 z(D) Figure 6.2 Generalized Wannier function wa. (z) on the second layer of a jellium slab with a square well potential. d) Td *i^ F4 o3 3.0 0.75 ' 0.50 0.25 0.00 0.25 I I I 12.0 9.0 6.0 3.0 0.0 3.0 z(D) Figure 6.3 Generalized Wannier function w ,(z) on the third layer of a jellium slab with a square well potential. 0) r0 0a 0 0.75 0.50 0.25 0.00 0.25 AI 12.0 9.0 6.0 3.0 0.0 z(D) Figure 6.4 Generalized Wannier function w' (z) on the fourth layer of a jellium slab with a square well potential. T0 44 3.0 Thble 6.2 Exact eigenvalues and approximate eigenvalues calculated by using generalized Wannier functions for a slab with a square well potential given by Eq. (6.31). Exact eigenvalue Approximate eigenvalue 0.8444 0.8444 0.8398 0.8397 0.8320 0.8320 0.8211 0.8209 0.8071 0.8068 0.7900 0.7895 0.7698 0.7690 0.7465 0.7453 0.7201 0.7185 0.6906 0.6883 0.6580 0.6549 0.6224 0.6181 0.5837 0.5778 0.5420 0.5340 0.4973 0.4864 0.4496 0.4349 0.3990 0.3792 0.3454 0.3185 0.2890 0.2517 0.2299 0.1769 0.1682 0.1299 0.20 0.10 0.00 0.10 0.20 ., I I I I , 25.0 20.0 15.0 10.0 5.0 0.0 5.0 z(D) Figure 6.5 Exact and approximate wavefunctions associated with the exact eigenvalue 0.8444 au for an onedimensional square well. O C3 0.20 . 0.10 0.00 0.10 0.20 I 25.0 20.0 15.0 10.0 5.0 0.0 z(D) Figure 6.6 Exact and approximate wavefunctions associated with the exact eigenvalue 0.6906 au for an onedimensional square well. 0 ' 4 U 2, 5.0 0.20 . i 0.10  0.00  0.10  1 I 0.20 25.0 20.0 15.0 10.0 5.0 0.0 z (D) Figure 6.7 Exact and approximate wavefunctions associated with the exact eigenvalue 0.3990 au for an onedimensional square well. CI 6O U 5.0 0.75 0.50 0.25 0.00 0.25 9.0 9.0 6.0 3.0 0.0 3.0 6.0 9.0 x(D) Figure 6.8 Exact and approximate Wannier functions for a slab in the directions parallel to the surface. The latter is a linear combination of Gaussians with an exponent /?=0.53. 0 3 +.a 0.75 0.50 1 0.25 0.00 0.25 II I I I 9.0 6.0 3.0 0.0 3.0 6.0 9.0 x (D) Figure 6.9 Exact and approximate Wannier functions for a slab in the directions parallel to the surface. The latter is a linear combination of Gaussians with an exponent #=0.52. 0.75 0.50 0.25 0.00 0.25 L 9.0 6.0 3.0 0.0 3.0 6.0 9.0 x(D) Figure 6.10 Exact and approximate Wannier functions for a slab in the directions parallel to the surface. The latter is a linear combination of Gaussians with an exponent 6=0.30. d) 'u P4 E i 86 6.3 Atomic Basis Functions Since in low energy atomsurface collisions only the valence electrons of atoms are actively involved in charge exchange, it is a good approximation to ignore the innershell electrons and use a pseudopotential for the atomic core. In this section we construct the atomic basis functions for 3s and 3p orbitals of Sodium and a pseudopotential for its atomic core. We choose Slatertype orbitals (STO) 0a (r; C, RA) as basis functions, where a=3s, 3px, 3py, and 3pz, C is the exponent of the orbitals, and RA is the position of the atomic core. For the convenience of calculation, we use STOKG functions fa (r; C, RA) which are linear combinations of K Gaussians and satisfy S(T C, RA) = Ca (; 1,RA) (6.53) and are in the form of K 3, (F; 1, A) = E d3s,k 9 (; s k, RA) (6.54) k=l K *3pj (r 1, RA) = L d3p,k g2.p (f; a, RA) (6.55) k=1 where j=x, y and z, and gl, (r ak, RA) and g2pj (r k, RA) are is and 2p type Gaussian primitive functions respectively, g1 (;ak ,RA = ( exp ak RA2) (6.56) 92p(,; ak, A) = (128 ) rep(ak RaI2) (6.57) (rt; a, RA 73 ) all 87 We use K=3 and C=1.27. The data for ak, d3s,k and d3p,k [Hehre et al., 70] are listed in Table 6.3. Table 6.3 Coefficients and exponents for Gaussians in STO3G functions. k ak d3s,k d3p,t_ 1 5.27266(2) 9.00398(1) 4.62001(1) 2 1.34715(1) 2.25595(1) 5.95167(1) 3 4.82854(1) 2.19620(1) 1.05876(2) 6.4 Overlap Matrix Elements Our localized basis set consists of generalized Wannier functions localized at the sites of the slab and atomic functions centered at the atomic core of the colliding ion as described in the previous sections, that is {f } = {{U }, {1a}1 (6.58) Since both {qS} and {4.} are orthonormal sets, the diagonal elements of the overlap matrix S are 1 and all the offdiagonal elements are null except Sai (A) = ( 0j) = Sa (RA) which are functions of the atomic core position RA. Using Eqs. (6.23), (6.53) and (6.54) (6.59) K Sai (RA) = E da,k Bif aG(ak, m, RA) k=1 rnh where Ga (ak Pm, RA) = (9a (ak, RA) 9(m,)) (6.60) is the overlap matrix of Gaussians and is also a functions of RA 6.5 Hamiltonian and Its Matrix Elements In this section we construct the one electron Hamiltonian for the system of a Na atom and a jellium slab and calculate its matrix elements in the basis described in the previous sections. The oneelectron Hamiltonian of the system is i =Iv 2 i + VA (r,RA) + VM( (6.61) where VA (r, RA) is the atomic potential and VM(i) is the potential of the slab which we choose to as the one given by Eq. (6.31). However, due to the use of the localized basis functions and the partition method, the effective Hamiltonian in the primary region should be reconstructed, which includes the construction of the pseudopotential for the atomic core and adding a correction term to get correct electronic couplings. 6.5.1 Pseudopotential for the atomic core In the present study, since we are not intending a full ab initio calculation, a pseudopotential for the atomic core is used for VA(, RA to eliminate the innershell electrons from the calculation [Kahn et al., 76]. The pseudopotential usually has a form of VA (f, RA) = V' (r, A) 9 (6.62) 1=0 where S= lm)(lm (6.63) m=I is the projection operator on Ilm) and the functions V1 can be chosen in one of many possible forms [Szasz, 85] and usually contain some parameters whose values are chosen to give agreement with either the results of ab initio calculation or with that of experiments. In this study we choose the form suggested by Schwartz and Switalski [Schwartz and Switalski, 72], ,Z N, exp(KIi RA2) VI RA = e + A RA (6.64) r RAI RA where Z is the atomic number, Ne the number of electrons of the core, and Al and Ka are the pseudopotential parameters. The first term is the Coulomb potential of a charge of ZNe and the second term is a correction due to the electrons in the inner shells. Using the pseudopotential the eigenequation for an isolated Na atom is [V + VAF, A a F, RA = Ea RA) (6.65) We determine the parameter A, and al by fitting the eigenvalues ca with the experimental values 3s,=0.18884 au and e3p=0.11156 au [Callaway, 69]. An other constrain imposed when choosing the parameters is that the pseudopotentials should become Coulombic at the distance where the core dies off, for example, at 3r7, where 3ri is the radius of orbital 1. However, the variation of parameter At is not too critical [Schwartz and Switalski, 72] and can be chosen in a cer 90 tain range. The parameters for our pseudopotential are listed in Table 6.4. The pseudopotentials for 3s and 3p orbitals are plotted in Fig. 6.11. Table 6.4. Pseudopotential parameters for 3s and 3p orbitals of Na atom. AI Kl 3s 1.0 2.415 3p 3.0 0.515 6.5.2 A Correction Term to the Hamiltonian The partition procedure described in Ch. 5 can be used to describe the evolution of the electronic states and the collisional charge transfer at short distances because of the local nature of the problem. But it complicates the description of the interaction between the surface and the atom when the atom is far away from the surface and the electrons relax, due to the use of the localized basis functions and the partition of the surface. Ignoring this relaxation would introduce an error. Although the error is small and does not effect the dynamics of electronic states to a great extent at small distance, the error accumulated at large distances could cause problems. This error induced by the partition could be significant in the case of collisional neutralization of ions. To represent the relaxation of the electrons in the localized basis functions we add to the full Hamiltonian a correction term loc = (F Ea)lXA)(X\al (6.66) riEp 6.0 ' 1v3  3 ' el"V 4.0 V3  6d S2.0 0 O 0 0.0 GO ............... 2.0 0.0 1.0 2.0 3.0 4.0 5.0 Z (au) Figure 6.11 Pseudopotentials for 3s and 3p orbitals of the Na atom. 92 where the summation is over the primary region p and e, = (xM HM Ix) with HM = V2 + VM the Hamiltonian of the slab. This term is constructed such that the asymptotic values of the electronic Hamiltonian are correct With this correction term, the effective Hamiltonian in the primary region is H' = H + Vioc (6.67) which will be used in the TDHF calculation in place of H . 6.5.3 Hamiltonian Matrix Elements In the primary region the matrix elements of the Hamiltonian, Eq. (6.67), are H'V = HO, + E (F eAi)SIASAv (6.68) vifEp or explicitly, H'a, = H.a, + E (eF ,Sa)SacaSAa (6.69) 'IEp H'j = H,, + (EF Ea)San (6.70) Hi'i = Hai + E ei (6.71) H',j, = Hiil, n n (6.72) The matrix elements of H are calculated as follows Ha, = (a V2 + VA + VMIa') 2 (6.73) = EaSaa, + (alVMla') 