A time-dependent molecular orbital approach to ion-solid surface collisions


Material Information

A time-dependent molecular orbital approach to ion-solid surface collisions
Physical Description:
viii, 150 leaves : ill. ; 29 cm.
Feng, Eric Quinn, 1951-
Publication Date:


bibliography   ( marcgt )
theses   ( marcgt )
non-fiction   ( marcgt )


Thesis (Ph. D.)--University of Florida, 1991.
Includes bibliographical references (leaves 144-149).
Statement of Responsibility:
by Eric Quinn Feng.
General Note:
General Note:

Record Information

Source Institution:
University of Florida
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
aleph - 001689210
notis - AJA1246
oclc - 25138709
System ID:

Full Text







To Susan and Michael and to My Mother's Memory


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.


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 Many-Electron 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
2.1 The Atom-Surface System ...................... 20
2.2 TDHF Equations of Density Matrix ................. 21
2.3 TDMOs as Linear Combinations of Travelling Atomic Orbitals 25
CHAPTER 3 ........ ................ ........... 28
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
4.1 The Linearization Procedure .................... 42
4.2 The Case without Electron-Electron Interaction .......... 48
CHAPTER 5 ................................. 51
5.1 The Partition Procedure ........................ 52
5.2 The Approximation in the Secondary Region ........... 55
CHAPTER 6 .. ................................ 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
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 Na-W(110) Model System ...........
8.1.1 Hamiltonian and Basis Functions .
8.1.2 Electronic Couplings ...........
8.2 Atom-Surface Interaction Potentials and the
8.2.1 Atom-Surface 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 .......................
APPENDIX A ......................
APPENDIX B ......................

APPENDIX C ..................


............ 101
............ 101
............ 102
Trajectory ...... 103
............ 103
............ 105
.......... 108
........... 111
. ........... 113
............ 125
. . 138
............ 138
. . 138
ND .......... 140
............ 142

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



Eric Quinn Feng

May 1991

Chairman: David A. Micha
Major Department: Physics

A time-dependent 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 time-dependent 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 near-resonance charge transfer mechanism. The time-dependence

of the density matrix provides insight on the dynamics of electron transfer along

the atomic trajectory.


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.


Charge transfer in ion-surface collisions may lead to neutralization or ioniza-

tion and can be expressed as

S+ A+ -- E+ A (1.1)


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 beam-foil 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 near-resonance

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 near-resonance 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 near-resonance processes. The Auger process involves

two or more electrons. In two-electron 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

where an ion affinity level is in resonance with the conduction band of the metal,

such as in Na++W(110), a near-resonance 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 ultra-high 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 ion-surface combinations, and is ascribed to near-resonance 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


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 ion-surface collisions.

1.2 A Survey of Theoretical Methods

Advances in spectroscopy techniques in the last two decades have allowed ex-

tensive experimental studies of ion-surface 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 many-electron processes with

a strong and time-dependent coupling between the electronic degrees of freedom

and nuclear motions. The traditional methods which deal with systems of periodic

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.

However, it cannot explain the oscillatory charge transfer behavior observed in

other systems which is believed to be associated with nonadiabatic charge transfer


1.2.2 Many-Electron 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 many-body processes accompanying charge exchange in ion-surface

collisions. The most widely used one is the approach based on the Anderson-

Newns Hamiltonian [Anderson, 61] to describe the ion-surface 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

degrees of freedom can be obtained from the Heisenberg equations (with h = 1)

i = -[H, Ca] = ca(t)Ca + Vak(t)Ck (1.7)

i = --[H, Ck] = Ck Ck + Vak(t)Ca (1.8)

Eliminating Ck, it yields a differential-integral 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 time-dependent Anderson-Newns Hamiltonian have

neglected the intra-atomic Coulomb repulsion to simplify the calculations. For

most alkali-metal collisions, which produce few negative ions, this approximation

seems to be justified. But for scatterings like H-W(110), there is a significant

fraction of H- products after scattering, and the intra-atomic 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
where a is a spin index, naT = CaCa, nal = C1CaJ, and U(t) is the effective

intra-atomic Coulomb repulsion which depends on the distance between the ion

and the surface. Using the Hartree-Fock approximation, one obtains a set of

differential-integral 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

Hartree-Fock 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 ion-surface scattering

will be addressed here. We view them as some key points in furthering the

theoretical study of ion-surface collisions.

1.3.1 Electronic States and Electronic Couplings

Ion-surface 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

semiconductor electronic structure is the existence of the localized states. While

s- and p-electrons are generally described by continuous states, d-electrons 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 semi-empirical calculation. For

example, in the time-dependent Anderson-Newns 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)
which can be approximately evaluated by methods developed in chemisorption

theory. But a common practice is to use parameterization of A(t). For instance,

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)

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 ion-surface systems should contain a self-consistent electronic

structure calculation which accounts for localized states. A self-consistent 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

ion-surface interactions [Lang and Williams, 76].

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. Ion-surface systems, on the other hand, are

many-body 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 ion-surface 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 self-consistent method using a density matrix

and applied it to the Anderson-Newns Hamiltonian for a one dimensional system.

McDowell [McDowell, 85] uses an embedding technique to obtain generalized

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 classical-quantal

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 gas-phase 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 H-p 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 so-called trajectory approximation or by a multichannel procedure. In the

trajectory treatment, the nuclear trajectory is uniquely defined by a single fixed

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 self-consistently. An application

of the eikonal method to ion-surface 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

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 off-diagonal 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 ion-surface 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 ion-surface

collisions based on the time-dependent molecular orbital approach. A partition-

ing technique is used to divide a ion-surface 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 time-dependent Hartree-Fock (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.

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 time-dependent 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 electron-electron interaction is neglected and obtain analytical solutions

for the linearized TDHF equations.


In Chapter 5, a partitioning procedure for ion-surface 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

time-dependent 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 time-dependent and in the form of

linear combinations of Gaussians. The overlap matrix and Hamiltonian matrix

are constructed in this basis set.


In Chapter 7, we discuss the computational aspects. At first, we present our

algorithm for numerical integration of the time-dependent Hartree-Fock 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 atom-surface 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 time-dependent molecular method

and its applications to charge transfer in atom-surface 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



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.


The time-dependent Hartree-Fock approach is suitable for studying the many-

body collision dynamics of our problem. Using a one-electron 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 self-consistency between the effective field and the orbital coefficients

or density matrix is achieved automatically by solving time-dependent differen-

tial equations without having to perform the self-consistent iteration procedure

required in the time-independent 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 ion-surface collisions.

In this Chapter, we apply the TDHF method to charge transfer in ion-surface

collisions. In Sect. 2.1, by introducing time-dependent molecular orbitals,


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 ion-surface 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 ion-surface 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 Atom-Surface System

Let's consider a system of a scattering ion and a semi-infinite solid. We will

assume the solid has a conduction band of continuous one-electron 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

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 near-resonance 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


2.2 TDHF Equations of Density Matrix

Let T (t) be the electronic wavefunction which satisfies the time-dependent

Shrodinger equation

i T(t) = H(t) W(t) (2.1)

where H(t) is the total Hamiltonian of the system. Using the time-dependent

Hartree-Fock approximation, the TDHF wavefucntions have the form of

I(t) = N"l1/2det[Oi(xj, t)]


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)

is the corresponding spin function. The electron density operator A is given by

E= Z )(01 (2.4)
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)

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

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 spin-orbitals of spin 7. Introducing a set of

time-dependent 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 time-dependent 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)
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)


H= K +VA+ VM (2.14)

is the core Hamiltonian matrix, G7 the Hartree-Fock electron-electron 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)

where 7' is the spin opposite to 7 and


(pVII\A) = (uV\jA) (p\IjAv)

(yIvAa)= d d'r2( t)(rt)=) 1 (Fi, t)o(rl, t)



To derive the TDHF equation for the density matrix, let's define a matrix


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




P-p -/ 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


where S = (J) is the overlap matrix. Multiply both sides of the above equation

by S-1 we obtain the time-dependent Hartree-Fock equation for the density matrix

i P = S-1FP'p P-F'rS-' iS-lP'7 ipt'rtS-1


The last two terms on the right hand side of Eq. (2.22) arise from the time-

dependence of the basis functions.

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)
where me is the electron mass and -m = dAm/dt is the local velocity of the mth

center in a space-fixed system. Thus,

"v = H i&') + n Vn vl) (2.25)



Using (2.22) we find

S= -ime S' + (xITT. n V.)l,) (2.27)


= 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 = S-1f'rp P'FtS-1 (2.31)

with a modified Fock-like matrix

F7 = I + G (2.32)


= + 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


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.


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 Ht-H

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

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 off-diagonal 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.

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

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 p-orbitals. The three eigenfunctions

of the angular momentum operator for J=1, Ip+i), Ipo) and Ip-l) correspond to

the magnetic quantum number m=1, 0 and -1 respectively. We can also use three

real p-orbital functions Ipx), IPy) and Ip,) which are symmetric along x-, y- and

z-axes respectively. The orbital jpz) = Ipo) has a negative reflection symmetry

with respect to the X-Y 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 atom-surface collision, the

collision plane is perpendicular to the solid surface. The collision frame (Xc,

Yc, Z7) is defined such that the Xc-Zc 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 Xc-Ye 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


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 Xc-Z7 plane of
the collision frame and the Xn-Yn plane of the natural frame.

Xc, YN


In the natural frame (Xn, Yn, Zn), as shown in Fig. 1.1, the Xn-Yn plane

is in the collision plane with the Xn axis perpendicular to the solid surface, the

Yn-Zn 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

z-component 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

J(0) = 1

1) = ri
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

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)
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)
ST-r(~) =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 body-fixed "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.

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)

J1 = :FJ+(2Jz + 1)

j(2) (3J_ 2 J2)

where J is the angular momentum operator and Jz is its z-component, 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.

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].

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 atom-surface 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 atom-surface 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
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


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)


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.


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 near-resonance 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 near-resonance charge transfer process. It tackles

the rapid variation of the density matrix by defining a time-dependent 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

time-dependent 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 electron-electron

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 near-resonance 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 time-dependent 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)


AS-1(t) = S-(i) So1 (4.3)


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 electron-electron interaction and
are given by

AGfeG(t) = C Q7(t)(VjIIIA)io + Q'A(t)( -iIA)to (4.6)


A.G1,(t) = _P.(t)A(pvl JA), + P,'(t)A(pvl A)t (4.7)
where (pVI|cA)t and (yjllcA)t are the Coulomb and Coulomb symmetrized
2-electron 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)FatAS-1(t) AS-z(t)qFPo(t) (4.8)

+So1AeGr(t)PO'(t) po'(t)AeG7(t)tSo1

+ So AnG'(t)POr(t) PO'(t)AnG'(t)tSo1

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)

9(t) = [AXKLPL'(t) + AYKLP''(t)] (4.12)
or in the matrix form,

AeG7= XoQ + YoQ' (4.13)

Angr = AX P-p0 + Ay p07' (4.14)

where Xo = X(to), Yo = Y(to) and

AX = X(t) -- Xo


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'-A-rt por'A7't

iQ- =A-4 "Q+A"''Qf'' QA Q'A7'7t + V7




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)


and the column matrix D-r is a driving term with elements

S= [So1(Af+ AG)PoY Pof(AHt+ AnGt)So1

+AS-'1j~p'o PO'ftAS-1],

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



into a column

( A"' Y 77'
A = Ay' A^-'- 4



then Eqs. (4.17) and (4.18) become

ilp = Apo poAt (4.23)


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 time-dependent.

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 time-ordered exponential expansion. We consider


Po(t) = U(t, to)lP(t)U(t, to)t (4.27)


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



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)


iQr(t) = iU(t, to)-DV(t) [U(t, to)t ]- (4.32)

Their solutions are

P(t) = Po(to) = P(to) (4.33)

Qp(t) = Qp(to) -i fu(t', to)'D (t')[U(t', to)t]-dt' (4.34)
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)

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'

4.2 The Case without Electron-Electron Interaction

In this section, we will show the use of our linearization procedure by

applying it to a special case where the electron-electron 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


iQ(t) = WQ(t) Q(t)Wt + D'(t)


w = S-fIo

Ho = fH(t0)



Unlike in the


D'(t) = S1AH(t)P0(t) P t (t)S
+Po(t)IAS-1(t) AS-1(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)]




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)


iQD(t) = exp[iW(t to)]D(t)exp[-iWt(t to)] (4.45)

The latter has a solution of
QD(t) = i exp[iW(t' to)]D'(t')ep[-iW(t' to)]dt' (4.46)
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)


Q(t) = -i exp -i(t t') D' (t')exp[iW ( t')]dt' (4.48)

We assume that the matrix W can be diagonalized by a linear transformation,

that is, there exists a matrix L such that,

W = LwL-1 (4.49)

where w is a diagonal matrix. Then

Q t(t) = -i j Lk [ (L-)kAkApl(t to )(L-) L (4.50)



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.


For an ion-surface 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.

For ion-surface 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 ion-surface system into two parts and treat electronic interactions and

electronic state evolution in these parts differently. This partitioning technique is

suitable to deal with time-dependent 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 electron-electron interaction is ignored, and obtain analytical solutions for the

density matrix in the primary region.

5.1 Partition Procedure

We divide an ion-surface 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

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
+WpP(t) Po (t)w

iP (t) = WppPo(t) Pp(t)Wt
+WpP(t) PO(t)W8

iP (t) = WP p()- Pp(t)W
+WP p(t) PO)(t)Wsp

iP(t) = W P(t- P (t)W p
+W PO (t) P (t)W

iQpp(t)=WppQpp(t) Qpp(t)Wp
+WpQsp(t) Qps(t)Wtp + D' (t)

iQpS(t)=WppQps(t) Qpp(t)Wp,
+WpQss(t) Qps(t)Wt + D',(t)

iQsp(t)=WspQpp(t) Qsp(t)Wpp
+WsQsp(t) Qss(t)Wtp + D'p(t)

iQs(t)=WspQps(t) Qsp(t)Wp,
+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)
+Pp(t)ii (AS-l)pp(t) (AS-) pp(t)topp p (t)

+Ppp(t)HI (AS-1)s(t) (AS-1)ps(t)IospPPp(t)

+Pps(t) p(AS-1)pp(t) (AS-l)pp(t)opsPsp(t
+P0s(t)HltIt(AS-l )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

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 time-dependent 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, 7-ear, 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

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)


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(AS-1)p(t) (AS-1)pp(t)ftopppp)(t) (5.15)

+Pps( tin)0sp(AS-1) P(t) (AS-1)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.


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 atom-surface 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 atomic-like basis

functions used in the tight-binding method. The local nature of the surface

phenomena suggests that some of the electronic behaviors of surfaces could

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 atomic-like 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 Kohn-Sham

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 atom-surface 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

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)

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(-) = ( ) e-iR'o (6.4)

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)


wi= (wi w, ,. ) (6.6)

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].

For a semi-infinite solid with its surface parallel to the x-y 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 Bloch-type 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,



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 well-behaved localization property of the generalized

Wannier functions combined with their orthonormality make them particular useful

to investigate the electronic structures or localized phenomena for non-periodic


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)rIe-p,(-Ai,)'


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


(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




(B)jj = Bjr


Substituting Eq. (6.17) into the normalization relation Eq. (6.13) yields

BtGB = I


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
G-2GG- = I





g = (gSaii, ga2ii, )

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.


E = (lfftPA ) (6.23)

the total energy of a band

EB = Tr(E)= Tr(({lIH|I))
( \ (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)

while being subject to the constraint Eq. (6.13).

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 R-4. 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.

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 near-periodic 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 semi-infinite 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

Hamiltonian is known, although the calculation of the three dimensional matrix
G- is a little time-consuming [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 z-direction 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)


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 y-directions, and

eigenvalue Er have the forms of

O(-) = q4(z)x (y)qS,(z) (6.35)


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 z-directions, 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)

w,(y) = yt) q o(y) (6.40)

,(z) = UZtk, (z) (6.41)

respectively. For this system, therefore, a three dimensional calculation is sim-

plified to a one dimensional one.

For this system, the x- and y-components of the eigenfunctions are

,= exp(ikzxz)
1 (6.42)
k,Y = 1 -exp(ikxz)
4 (Ny 1)d
the generalized Wannier functions in the x- and y-directions are just the Wannier

function for periodic lattice

1 sin 1-
W =,() = (6.43)
(Nn 1) d sin [r(s-n,d)
1 sin -*' dd)

W (y) s (6.44)
S(Nz 1) sin r(Y-nyd)
L (N-1)d
and the component transformation matrices Ux and UY are

=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)

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, *,/M-1 = #N/-M+2, but /M = /3M-i = = 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)

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( (6-m( )) 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.

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

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,


-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 point-by-point

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.






-9.0 -6.0 -3.0 0.0



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.







-0.25 -- I, I I
-12.0 -9.0 -6.0 -3.0 0.0


Figure 6.2 Generalized Wannier function wa. (z) on the
second layer of a jellium slab with a square well potential.





0.75 '




-0.25 I I I
-12.0 -9.0 -6.0 -3.0 0.0 3.0


Figure 6.3 Generalized Wannier function w ,(z) on the
third layer of a jellium slab with a square well potential.







-0.25 AI
-12.0 -9.0 -6.0 -3.0 0.0


Figure 6.4 Generalized Wannier function w' (z) on the
fourth layer of a jellium slab with a square well potential.



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
-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 ., I I I I ,
-25.0 -20.0 -15.0 -10.0 -5.0 0.0 5.0


Figure 6.5 Exact and approximate wavefunctions associated with the
exact eigenvalue -0.8444 au for an one-dimensional square well.



0.20 .




-0.20 I
-25.0 -20.0 -15.0 -10.0 -5.0 0.0


Figure 6.6 Exact and approximate wavefunctions associated with the
exact eigenvalue -0.6906 au for an one-dimensional square well.

' -4



0.20 .


0.10 -

0.00 -

-0.10 -

1 I
-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 one-dimensional square well.








-6.0 -3.0 0.0 3.0 6.0



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.





1 0.25-


-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.25 L-

-6.0 -3.0 0.0 3.0 6.0



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.




6.3 Atomic Basis Functions

Since in low energy atom-surface collisions only the valence electrons of
atoms are actively involved in charge exchange, it is a good approximation to
ignore the inner-shell 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 Slater-type 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 STO-KG 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
3, (F; 1, A) = E d3s,k 9 (; s k, RA) (6.54)

*3pj (r 1, RA) = L d3p,k g2.p (f; a, RA) (6.55)
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

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 STO-3G 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


Since both {qS} and {4.} are orthonormal sets, the diagonal elements of

the overlap matrix S are 1 and all the off-diagonal 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)


Sai (RA) = E da,k Bif aG(ak, m, RA)
k=1 rnh


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 one-electron 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

inner-shell electrons from the calculation [Kahn et al., 76]. The pseudopotential

usually has a form of

VA (f, RA) = V' (r, A) 9 (6.62)


S= |lm)(lm| (6.63)
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)

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 Z-Ne 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-

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.

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)

6.0 '

1v3 -
3 '
4.0- V3 -



0 0.0
GO ...............

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.

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)
or explicitly,

H'a, = H.a, + E (eF ,Sa)SacaSAa (6.69)

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')